水力压裂改造技术已成为页岩气和致密砂岩气等非常规低渗透性储层有效的增产技术, 对压裂过程中产生的微地震活动进行监测, 是评价压裂改造效果的一种重要手段^[1-2]。微地震监测主要有井中监测与地面监测两种方式。井中监测方式采集的数据质量较高, 但受检波器数量与观测范围限制, 对震源参数(震源位置与机制)反演约束较差, 而且如果在压裂井附近缺少观测井的情况下则无法进行监测^[3]; 地面监测方式观测范围更大, 覆盖次数更高, 因而可以避免此类问题的出现, 但其主要问题是受地表环境因素影响, 资料品质较差、信噪比较低^[4-7]。

传统的地面微地震监测资料的处理通常分为两个步骤:先进行震源定位, 然后利用定位信息进行震源机制反演^[8]。因为地面微地震监测资料信噪比较低, 初至拾取较为困难^[9-10], 所以格点搜索类方法成为震源定位反演的一种常用方法^[11]。该方法对震源位置空间进行格点离散并搜索, 目标函数采用叠加函数的形式, 叠加能量最大的格点即为震源位置。实际资料显示微地震震源机制有很强的双力偶震源成分, 其振幅极性在地表呈四象限分布^[12]。在进行震源定位时, 目标函数若采用直接叠加的方式, 振幅极性的反转会造成微震叠加过程中能量的抵消。采用绝对值叠加函数可避免这个问题, 但噪声水平会同步放大, 并不能很好地突出微震事件的叠加能量。叠加函数的选择会直接影响震源定位的精度, 在定位偏差较大的情况下进行震源机制反演, 会将定位误差带入到震源机制反演中, 造成震源机制反演误差较大, 这对后续的微地震监测资料解释造成困扰。

目前主要有两种方法解决振幅极性反转带来的问题。一种是基于互相关的地震干涉偏移类方法^[13]。ZHEBEL等^[14-15]通过对相邻道做互相关处理来消除振幅极性反转的影响, 从而实现了震源定位与震源机制联合反演; LI等^[16]基于该方法提出了加权弹性波干涉成像方法, 并应用于多波多分量记录的微地震定位。另外一种方法是通过震源定位与震源机制联合反演, 利用震源机制来校正微地震记录的极性。LIANG等^[17]采用格点扫描叠加方法, 对离散震源位置与震源机制组成的全空间进行离散扫描, 并用震源机制校正记录极性。ANIKIEV等^[18]在每个搜索格点进行震源机制反演来获得震源机制, 并用此震源机制来计算极性校正后的叠加函数。

本文介绍了一种格点扫描的震源定位与震源机制联合反演的迭代搜索算法。首先根据常规震源定位方法获得震源初始位置, 采用P波初动震源机制反演获得震源机制, 用该震源机制校正微地震记录并进行定位, 迭代此过程以获得较准确的震源位置, 然后采用P波辐射花样约束的震源机制反演方法, 以获得较准确的震源机制。最后给出了模型测试结果, 验证该方法的有效性和优越性; 并将该方法应用于中国四川盆地某压裂井的地面微地震监测实际资料的处理。

1 方法原理

传统的地面微地震监测处理流程是先进行震源定位再进行震源机制反演, 这种方式会造成误差的传递; 而对由震源位置与震源机制组成的震源参数全空间进行格点搜索的联合反演方式, 由于搜索样本数目太大, 对计算能力有较高的要求^[17]。

因此本文提出了一种震源位置与震源机制迭代搜索方法, 技术流程如图 1所示。采用这种方法既可以避免震源参数全空间同时离散造成的搜索格点急剧增长, 也可以叠加目标函数振幅校正使其更具有实际物理意义, 并且根据迭代过程中对震源机制精度需求的不同, 本文采用了两种震源机制反演方法, 进一步提高反演速度。

我们首先采用常规的定位算法对微地震事件进行初始定位, 获得震源初始位置后, 进行震源定位和震源机制迭代反演, 具体如下:首先进行P波初动极性反演获得震源机制, 然后用此震源机制来校正微地震记录, 并采用极性校正叠加函数重新进行定位, 最后迭代计算直至震源位置收敛, 此时可以获得较为准确的震源位置与初步的震源机制结果。

由于地面观测系统覆盖区域有限, 可能无法对某一震源机制节面位置进行采样, 因而P波初动极性震源机制反演存在多解性问题^[19]。但即使得不到准确的震源机制解, 也可以提供正确的检波点初动极性, 故不会影响定位结果。在获得准确的定位结果后, 我们可进一步采用更精确震源机制反演方法来反演震源机制, 本文采用的是辐射花样约束的震源机制反演方法。

在采用扫描叠加搜索方法进行震源定位时, 叠加目标函数一般有以下3种方式:

$ {{\mathit{f}}_{\text{1}}}\left( \mathit{t} \right)\text{=}\frac{\left| \sum\nolimits_{\mathit{i}=1}^{\mathit{N}}{{{\mathit{u}}_{\mathit{i}}}}\left( \mathit{t}-{{\mathit{\tau }}_{\mathit{i}\text{, }\mathit{x}}} \right) \right|}{\mathit{N}} $

(1)

$ {{\mathit{f}}_{2}}\left( \mathit{t} \right)\text{=}\frac{\sum\nolimits_{\mathit{i}=1}^{\mathit{N}}{\left| {{\mathit{u}}_{\mathit{i}}}\left( \mathit{t}-{{\mathit{\tau }}_{\mathit{i}\text{, }\mathit{x}}} \right) \right|}}{\mathit{N}} $

(2)

$ {{\mathit{f}}_{\text{3}}}\left( \mathit{t} \right)\text{=}\frac{\sum\nolimits_{\mathit{i}=1}^{\mathit{N}}{{{\mathit{s}}_{\mathit{i}}}{{\mathit{u}}_{\mathit{i}}}}\left( \mathit{t}-{{\mathit{\tau }}_{\mathit{i}\text{, }\mathit{x}}} \right)}{\mathit{N}} $

(3)

式中:u_i为第i道的记录; τ_{i, x}为震源在x点到第i检波点的旅行时; s_i为第i道P波初动极性; f₁(t)为叠加函数的绝对值; f₂(t)为绝对值叠加函数; f₃(t)为极性校正叠加函数。

假设微地震震源机制主要成分为双力偶型震源, 则其振幅在地表呈四象限分布。若采用公式(1)直接叠加方式, 则导致正负区域振幅能量相互抵消。甚至当震源倾角为近90°, 观测系统在震源上方对称分布时, 即使事件震级很大, 但因f₁(t)值为零, 故无法对其进行准确定位。如若目标函数采用公式(2)绝对值叠加方式, 虽然可以避免出现振幅抵消的情况, 但该函数会同时提高噪声水平, 不利于目标函数极值拾取。本文在迭代反演过程中采用的是公式(3)极性校正叠加:利用震源机制得到的极性来校正微地震记录, 然后进行叠加。这种方法可以避免采用公式(1)与公式(2)带来的问题, 以获得较好的叠加效果^[17-18]。

在整个反演流程中用到了两种震源机制反演方法。首先在震源定位和震源机制迭代反演过程中采用了P波初动极性震源机制反演算法^[19], 这种算法只利用了P波初动极性信息, 故计算速度较快, 而且可以提供一个大致准确的震源机制, 以满足定位时振幅校正的需要。该算法对震源坐标系在地平坐标系中的微震进行全空间扫描, 计算矛盾极性比ψ:

$ \mathit{\psi }\text{=}\frac{{{\mathit{N}}_{\text{C}}}}{{{\mathit{N}}_{\text{T}}}} $

(4)

式中:N_C为矛盾极性数, N_T为极性总数。当矛盾极性比ψ最小时, P波初动极性与观测极性拟合最佳, 即为求取的一组震源机制节面。

获得准确震源定位后, 为了避免P波初动极性震源机制反演的多解性, 获得更为精确的震源机制, 本文进一步采用P波辐射花样约束下的震源机制反演算法^[20]。该算法通过对震源机制全空间扫描, 求取与实际P波辐射花样吻合度最好的一组震源机制节面, 采用的反演目标函数是实际观测和理论计算之间的相关系数R(吻合度)^[6]:

$ \mathit{R}\text{=}\frac{\sum\limits_{\mathit{i}\text{=1}}^{\mathit{N}}{\text{(}{{\mathit{U}}_{\mathit{i}}}-\mathit{\bar{U}}\text{ )(}{{\mathit{V}}_{\mathit{i}}}-\mathit{\bar{V}}\text{)}}}{\sqrt{\sum\limits_{\mathit{i}\text{=1}}^{\mathit{N}}{{{\text{(}{{\mathit{U}}_{\mathit{i}}}-\mathit{\bar{U}}\text{ )}}^{\text{2}}}}}\cdot \sqrt{\sum\limits_{\mathit{i}\text{=1}}^{\mathit{N}}{{{\text{(}{{\mathit{V}}_{\mathit{i}}}-\mathit{\bar{V}}\text{)}}^{\text{2}}}}}} $

(5)

式中:U_i为第i道检波器记录到的直达P波振幅; U为所有N道的平均振幅; V_i为在扫描格点上计算得到的第i道直达P波振幅; V为计算得到的平均振幅; R越接近于1, 说明U_i与V_i两者吻合度越高。

2 模型测试

为了验证本文方法的可靠性, 我们利用理论模型数据进行了测试。由于受上、下地层与围压因素影响, 目的层埋深超过1 000 m时, 水力压力裂缝呈高角度裂缝^[21], 所以设计微地震震源模型采用的裂缝参数:方位、倾角、滑动角分别为20°, 90°, 40°, 位于(0, 0, 1 500 m)处。速度模型和观测系统如图 2所示, 在深度为1 000 m处有一高速层; 检波阵列呈米字状分布, 共320道检波器, 检波距为50 m, 最大偏移距为2 000 m, 最大偏移距与震源深度比约为1.3, 与实际水力压裂地面微地震监测观测系统分布相似。本文采用CHEN等^[22]的广义反射-透射系数正演模拟算法生成理论地震记录, 如图 3所示。

从图 3中拾取初至P波同相轴并截取初至P波波形, 通过与理论子波互相关, 计算各道初至P波的初动极性和振幅, 利用本文提出的联合反演方法反演震源位置与震源机制参数。其定位结果为(0, 0, 1 500 m), 震源机制反演结果为(21°, 90°, 40°), 如表 1中所示的信噪比为∞时的结果。而采用传统反演方法, 得到的位置为(0, 0, 1 508 m), 震源机制为(21°, 90°, 40°), 如表 2中所示的信噪比为∞时的结果。可以看出, 在没有噪声的情况下, 两种方法都得到较准确的结果, 但联合反演方法与传统方法相比, 垂直方向上定位精度更高。图 4给出了模型理论P波振幅与联合反演得到的振幅对比结果, 图中圆形表示理论振幅, 方形表示反演结果, 图标大小表示振幅相对大小, 红色和蓝色分别表示正负极性。可以看出, 理论振幅与反演振幅无论极性还是振幅相对大小分布都非常吻合, 验证了本文提出方法的有效性。

为了验证本文方法的抗噪性, 我们从实际野外地面微地震监测记录中截取一段背景噪声(图 5), 并与图 3所示的理论地震记录合成不同噪声水平的地震记录, 并对此记录进行反演测试。表 1给出了本文提出的联合反演方法的反演结果, 表 2给出了传统反演方法的反演结果。图 6对比了两种方法定位结果的误差, 其中蓝线为联合反演的定位误差, 红线为传统方法的定位误差, 可以看出, 相比于传统反演方法, 本文的联合反演方法定位精度更高, 且抗噪性能更好。对比表 1与表 2可以发现, 联合反演方法在垂直方向上具有更高的精度, 说明该方法可以提高垂直方向上的反演约束。图 7对比了不同信噪比条件下利用两种反演方法进行震源机制反演的结果误差, 图 7a为传统反演方法震源机制反演误差随信噪比变化结果, 图 7b为联合反演方法震源机制反演误差随信噪比变化结果。可以看到, 无论反演精度还是稳定性, 联合反演方法都要比传统方法好。

传统反演方法将震源定位与震源机制反演分开进行, 并且在震源定位时采用绝对值叠加函数, 使得幅值正负分布的噪声变成幅值为正的噪声。随着噪声水平的增加目标函数信噪比逐步降低, 计算目标函数极值误差逐步增加, 所以其定位误差也逐步增加。在定位误差增大的情况下, 再进行震源机制反演, 必然会将定位误差引入, 导致震源机制反演误差增大。而在联合反演方法中, 定位时采用的是极性校正叠加函数, 而不是绝对值叠加, 这样叠加时噪声仍呈正负随机分布, 通过叠加可以有效压制噪声, 使得定位精度的提高, 并有助于震源机制反演精度的提高。

另外本文提出的联合反演方法迭代收敛非常迅速, 一般2到3个循环即可收敛。原因在于迭代过程中只要反演获得的P波初动极性与观测极性拟合较好, 其定位结果会迅速收敛。P波初动极性震源机制反演的主要目的是为震源定位提供一个准确的极性校正函数, 而该反演算法也可以快速完成。所以联合反演方法所用的运算时间大概是传统反演方法2~3倍, 而与全空间离散扫描方法相比, 反演速度提高了10倍以上。

3 实际资料应用

采用本文提出的联合反演方法对四川盆地某压裂井地面微地震实际监测资料进行震源定位与震源机制联合反演研究。利用声波测井曲线建立速度模型如图 8a所示, 压裂监测所采用的观测系统如图 8b所示, 射孔段深度约为3 294~3 369 m。共10条测线, 呈放射状分布, 相邻测线方位角相差约36°, 每条测线有125个检波点, 道间距为25 m, 最小偏移距为300 m, 最大偏移距为3 400 m。本次监测共识别了1 073个微地震事件, 从中选择109个较强微地震事件进行处理分析。

图 9给出了一个典型微地震事件记录, 表 3给出了对该微地震事件分别采用传统反演方法和本文提出的反演方法的反演结果。可以看出, 利用联合反演方法得到的震源机制的反演结果吻合度更高, 可以认为反演结果更为可靠。图 10对比了联合反演结果振幅与实际振幅, 可以看出极性、振幅大小以及分布吻合得很好。图 11显示了所处理的109个微地震事件的震源机制吻合度, 平均值为0.945 2, 与传统方法得到的震源机制吻合度平均值0.938 0^[6]相比有所提升, 可以认为本次反演结果精度整体有所提升。

图 12对比了联合反演方法与传统反演方法震源定位结果, 结合图 13显示的由联合反演方法得到的震源位置与震源机制结果, 可以看出联合反演方法得到的震源节面连续性更好, 水力压裂产生的缝网发育方向以北偏东方向为主, 可以认为其更能准确反映裂缝的延展走向。图 14给出的震源机制分布显示, 压裂裂缝面的延展方向也是北偏东方向(10°~20°), 两者较为一致。

4 结论

本文介绍了一种基于格点扫描方式的震源定位与震源机制联合反演的迭代搜索算法。与传统震源参数反演方法相比, 该方法采用定位结果与震源机制结果互相约束, 可以同时提高震源位置与震源机制的反演精度; 并且采用迭代搜索方式, 避免了同时搜索震源位置与震源机制空间造成的计算量过大问题, 提高了反演速度。

如果在进行格点搜索时考虑采用多尺度搜索策略, 可以进一步降低搜索格点数目, 提高计算效率。另外若将所有格点旅行时预先计算好并存储起来, 可避免迭代过程中重复计算, 搜索效率将进一步提升。

在获得震源位置与震源机制信息后, 可以进一步结合地质资料、地层应力背景以及岩心信息来选择具体裂缝破裂面, 用于裂缝空间几何展布和缝网演化过程等压裂解释工作。