基于各向异性分析的微地震震源矢量场重建和裂缝解释

宋维琪，徐奔奔，喻志超，秦晅，张宇

中国石油大学（华东）地球科学与技术学院，青岛 266580

1 引言

微地震监测技术在国内外（Prugger and Gendzwill，1988；Rutledge et al.，1998；Jupe et al.，2000；姜福兴等，2006；叶根喜等，2008；王晨龙等，2013）已经研究应用多年，但是到目前为止，微地震裂缝解释还没有系统的理论基础和有效的解释方法.尽管天然地震研究领域提供了一些震源的解释方法（姚振兴等，1994；王卫民等，2013；Stanton，2013；刘培洵等，2014），但由于研究尺度、震源受力的方式、微地震事件的多样性等诸多方面的不同，照搬天然地震领域的理论和方法，往往会导致微地震裂缝解释出现较大的偏差，甚至会得出错误的解释结果.

当然，天然地震的震源机制理论（靳平等，1998）、声发射效应理论、岩石破裂理论是微地震研究的重要的基础理论.本研究在前人（Stoffa and Sen，1991；Furumura et al.，1998）研究基础上，结合实际生产应用中的具体问题，考虑到微地震压裂的地下环境，从分析震源辐射（刘恩儒等，2006）特征入手，进行岩石破裂和震源响应特征关系研究，这对后续的裂缝优化解释有着理论指导意义，进而基于速度各向异性构建格林函数场（张钋等，2000），从实际观测记录中还原震源矢量场信息，进而提出一种全新的裂缝解释方法.

2 岩石破裂、能量释放及地震波特征研究

影响岩石破裂的主要因素有：岩石的受力方式、非均匀性、结构、岩性与物性、岩石厚度与埋深、应力场分布等.一般情况下，岩石破裂程度与施加力的一般规律是：剪切力大于拉张力，拉张力大于压缩力；一般岩石越均匀，越不容易破裂；如果岩石发育了纹理、微裂缝结构，增大了岩石的破裂点，则更容易破裂；不同岩石类型和不同物性的岩石，破裂的难易程度不同，硬度偏大与硬度偏小的岩石难于破裂；岩石厚度大埋深大更不容易破裂；对于微地震监测的水力压裂，是由于地层岩石在液体压力、不同方向地应力作用下，原有力平衡被打破，岩石发生破裂.

2.1 岩石破裂能量释放特征分析

（1）释放能量大小和裂缝破裂大小有关，具体与裂缝前缘面积和推进长度有关.V＝S×L.其中S是裂缝的前缘面积，L是瞬时推进长度.

（2）释放能量大小和岩石的类型有关，一般情况，脆性岩石破裂释放的能量大，塑性岩石破裂的能量小.

（3）岩石能量释放最大方向和岩石裂缝破裂方向有关，岩石地层裂缝开时的朝向不同，能量释放的主方向不同，一般岩石破裂裂缝能量释放的主方向和裂缝的开启方向一致.

2.2 岩石破裂过程及微地震事件类型、特点分析

通过大量的研究实践发现，塑性岩石在达到破裂点之前，能够观测到其能量体现在振幅变化上.岩石达到破裂点后，由于能量的释放、围压的减小，使得裂缝发生部分闭合，闭合过程中又释放一些能量，又形成新的微地震事件记录.

因此，在同一个位置点处，可以产生三种不同类型和特点的微地震事件，即开启事件、闭合对偶事件及压裂液瞬时变化产生的微地震事件，它们的位置是不变的，因此走时是不变的，但振幅是变化的.

3 震源矢量场的构建

3.1 各向异性条件VTI介质下走时计算

3.1.1 各向异性条件下地层速度

Thomsen给出了VTI介质中的相速度（朱光明等，2008）的近似表达式为

其中，θ为相角，Vp0，Vsv0，Vsh0分别为 qP、qSV 和qSH波垂直方向传播速度，δ，ε，γ为各向异性参数，Vp，Vsv，Vsh分别表示 VTI介质中qP、qSV、qSH 波的相速度.

4.中英文的缩写。中英文的缩写是指用英文单词的首字母表示该单词，用每个汉字汉语拼音的第一个字母表示该字。例如：

群角和相角关系为

群速度与相速度的关系为

3.1.2 各向异性条件下走时计算

在各向异性条件下，求取地震波的路径走时，关键问题是如何求取相速度、群速度以及相角、群角.在不考虑各向异性条件下，近似设定相速度和群速度相等.通过以往的射线追踪方法（Moser，1991；Leidenfrost et al.，1999），可以求得各个地层分界面的入射角和出射角.在考虑近似情况下，约定相角等于群角.得到不同地震射线的入射角和出射角以后，如果再求得各向异性参数，则代入公式（1）、（2）、（3）可计算得到各向异性条件下的相速度.群速度取以上各向同性条件下的速度.在求取了群速度、相速度及相角以后，代入公式（4）可计算出群角，求得相速度后，利用公式（5）可以计算群速度.各向异性条件下走时的计算公式为

式中Ri为各地层射线路径.

各向同性和各向异性条件下走时的模拟记录结果如图1和图2所示.检波器垂直排列接收，在地层中从上到下，对应的记录道序号是从小到大.震源在最后一级检波器下面.分析图中结果，在两图的左侧对应的小的道号即在上边的检波器道，离震源较远，对应较小的入射角，各向同性和各向异性两种情况的走时几乎相同，但是在两图的右侧，离震源较近，入射角较大的时候，各向同性和各向异性两种情况的走时，不再相同.

3.2 各向异性条件透射系数求取

Ruger A通过研究给出了VTI介质考虑各向异性透射系数的近似计算公式，包括各向同性项和各向异性项，其中各向同性项为（李幼铭，1988）

图1 各向同性走时模拟结果Fig.1 Travel－time simulation result of isotropy

图2 各向异性走时模拟结果Fig.2 Travel－time simulation result of anisotropy

为了得到更准确的透射系数，这里把左普里兹方程组解出的透射系数作为各向同性系数，各向异性项透射系数的计算根据上述近似公式计算.如图3，图中实线为各向同性条件结果，虚线为各向异性条件结果.横坐标为入射角度，纵坐标为透射系数大小.从结果图中可以看出，各向异性对透射系数的影响随入射角的不同而不同，在小入射角附近，各项异性和各向同性透射系数差别不大，随着入射角的增大二者差别明显变大，特别是在超临界角情况下.

3.3 格林函数微地震记录正演模拟

在考虑各向异性条件下，求得了地层模型的透射系数、群速度、相速度及事件位置等各参数后，利用高斯束射线追踪方法进行地震记录的正演合成（符力耘等，1994；王山山等；1997），求得均匀辐射单位点震源的VTI介质速度模型各向同性及各向异性条件下的微地震波场模拟结果.如图4和图5所示.比较分析两模拟结果，其波场特征不同.

3.4 震源矢量场重建

图3 透射系数随入射角的变化Fig.3 Variation of transmission coefficient with incidence angle

图4 各向同性模拟结果Fig.4 Simulation result of isotropy

图5 各向异性模拟结果Fig.5 Simulation result of anisotropy

为了得到震源极化的全矢量场，充分利用随机分布的多级检波器水平分量信息.把格林函数模拟场结果投影到多级检波器水平分量和垂直分量上，把多级检波器的观测结果减去格林函数多级检波器分量投影结果，得到震源场在多级检波器的投影分量结果，称为剩余场分量.然后对剩余场分量利用偏振分析方法，求得特征矢量场，得到相对较全的震源矢量场.

4 裂缝解释方法

4.1 根据震源矢量场估计裂缝开启方向

把经过去掉地层效应的实际观测的记录结果和上述模拟波场记录进行对比分析，通过偏振分析，研究其特征向量，根据三个特征向量的大小，判断分析源能量的辐射大小和方向.

利用正演模拟结果，建立不均匀震源相对均匀震源方位角偏离大小的定量关系；建立多级检波器Z分量与震源均匀与非均匀辐射的定量变化关系；综合利用上述研究得到的震源矢量场的特征矢量场、Z分量的空间变化关系及方位角的变化关系，确定裂缝的开启方向.具体方法是，首先通过震源特征矢量场找出最大特征特征矢量，把该特征矢量的方向确定为裂缝的开启方向，然后再利用以上建立的两种关系，进行进一步分析，最终确定裂缝的开启方向.

4.2 裂缝宽度、延伸长度估计方法

对重建的震源矢量场，根据事件的不同位置和地层岩性，进行空间能量扩散校正、吸收衰减校正，然后计算总能量.对于开启方向相同的裂缝，能量越大则裂缝宽度越大.

如果裂缝开启方向都有微地震事件，则所有事件的空间位置分布，就指示了裂缝带的延伸范围大小.

4.3 裂缝解释与优化

压裂裂缝的生长具有一定的不确定性，如果发育多条裂缝，则整个压裂区域形成一个不确定的裂缝网络.针对此问题，研究提出了微地震定位结果的简约优化方法.其方法的基本思路是：根据震源矢量场提取其特征属性，如最大主能量大小和方向、总能量、事件的方位、事件的空间位置，利用提取的震源属性，采用聚类分析方法，把相同类的多个事件合并约简，得到最终优化后的裂缝网络.

4.4 实际结果分析

以下为3个不同地区的压裂微地震监测结果，A、C地区是直井压裂施工监测，A地区压裂目的层段深度是2135m，岩石岩性为粉砂岩.C地区压裂目的层段深度2300m，岩石岩性为灰岩.B地区是水平井压裂施工监测，目的层深度2567m，水平井段长度是800m，实施分段压裂，间隔100m.

图6 A地区监测结果平面图（a）初始定位结果；（b）优化定位结果.Fig.6 Planar monitoring result of A area（a）Initial positioning result；（b）Optimal positioning result.

图7 B地区监测结果平面图（水平井压裂的其中两段）（a）初始定位结果；（b）优化定位结果.Fig.7 Planar monitoring result of B area（a）Initial positioning result；（b）Optimal positioning result.

图8 C地区监测结果平面图（a）初始定位结果；（b）优化定位结果.Fig.8 Planar monitoring result of C area（a）Initial positioning result；（b）Optimal positioning result.

首先对微地震监测结果进行了一系列的处理、事件拾取及反演定位工作，得到了微地震监测定位结果.为了更好地进行裂缝发育结果解释，采用研究的方法对定位结果进行了优化.优化过程采用定位结果和实际构造分布及压裂工艺综合裂缝网络优化控制方案.利用定位结果的空间（这里分析平面情况）点的分布密度，提取裂缝带的轮廓特征.在网络优化时，除了采用裂缝网络和局部裂缝带优势方向一致控制外，还实施构造发育方向控制，如果有定向加压压裂，还需加入定向方向控制.这里研究的3个地区都是常规压裂，因此，在裂缝网络优化过程中主要实施了构造发育方向控制方案.

图6a为压裂监测定位结果，图6b为对其优化结果.优化过程中，在实施局部裂缝带优势方向控制下，还加入了该地区东北方向的构造发育方向控制.对比分析图中结果看到，裂缝网络的优势发育方向为东北方向，并且优化结果更加聚团、收敛.图7为水平井压裂定位结果，原来定位结果中两段界线模糊，经过网络优化处理后，两段的界线清楚.分析图8a和b，通过对比分析，优化前结果定位点分布凌乱，裂缝带发育的优势方向不清楚；优化后网络的优势方向为北西方向，在该图的上部，又有两条东北方向的局部裂缝带发育.

5 结论

（1）地层的各向异性引起地震走时和振幅异常，随着地震射线入射角的不同，这种异常变化呈现明显的非线性特性，特别是在超临界角情况下这种异常变化尤为明显.

（2）微地震裂缝解释的正确与否和震源辐射场特征具有密切的关系，要得到完整的震源矢量场，必须首先去掉地层影响的传播效应.利用多级检波器分量构建完整水平分量场，利用重建后的场再定位，定位结果无论精度还是稳定性方面都具有较大程度的提高.

（3）利用重建的震源矢量场，根据裂缝破裂产生的矢量场特征分析，形成了裂缝破裂方向、裂缝破裂面积大小震源矢量场定量解释方法.在分析岩石受力破裂过程的基础上，提出了同位置点微地震事件的概念及其同位置点事件特征和识别方法，结合微地震定位结果的精度限制分析，建立了裂缝及裂缝网络优化新方法.

实际资料的测试分析，证明了各项研究方法的应用效果.研究的技术方法为微地震解释提供了重要的解释理论指导.

Fu L Y，Mou Y G.1994.Boundary element method for elastic wave forward modeling.Acta Geophysica Sinica （in Chinese），37（4）：521－529.

Furumura T，Kennett B，Takenaka H.1998.Parallel 3－D pseudospectral simulation of seismic wave propagation.Geophysics，63（1）：279－288，doi：10.1190／1.1444322.

Jiang F X，Yang S H，Cheng Y H，et al.2006.A study on microseismic monitoring of rock burst in coal mine.Chinese Journal of Geophysics （in Chinese），49（5）：1511－1516.

Jin P，Xu G M，Lou W T.1998.Elastic waves from a point force in transversely isotropic media.Acta Geophysica Sinica（in Chinese），41（4）：525－536.

Jupe A，Cowles J，Jones R，et al.2000.Micro－seismic monitoring：Listen and see the reservoir.World Oil，12：171－173.

Leidenfrost A，Ettrich N，Gajewski D，et al.1999.Comparision of six different methods for calculating travel times.Geophysical Prospecting，47（3）：269－297.

Li Y M.1988.Transparent wave field reconstruction and the thickness inversion of the layered media.Acta Geophysica Sinica （in Chinese），31（6）：649－656.

Liu E R，Yue J H，Pan D M.2006.Frequency－dependent anisotropy：Effects of multiple fracture sets on shear wave polarisations.Chinese Journal of Geophysics（in Chinese），49（5）：1401－1409.

Liu P X，Chen S Y，Guo Y S，et al.2014.Moment tensor inversion of acoustic emission.Chinese Journal of Geophysics（in Chinese），57（3）：858－866，doi：10.6038／cjg20140315.

Moser T J.1991.Shortest path calculation of seismic rays.Geophysics，56（1）：59－67.

Prugger A F，Gendzwill D J.1988.Microearthquake location：A nonlinear approach that makes use of a simplex stepping procedure.Bull.Seism.Soc.Am，78（2）：799－815.

Rutledge J T，Phillips W S，House L S，et al.1998.Microseismic mapping of a Cotton Valley Hydraulic Fracture using decimated downhole arrays.International Exposition and Sixty Eighth Annual Meeting.New Orleans，Louisiana.

Stanton A，Mauricio D S.2013.Vector reconstruction of multicomponent seismic data.Geophysics，78（4）：V131－V145，doi：10.1190／geo2012－0448.1.

Stoffa P L，Sen M K.1991.Nonlinear multi－parameter optimization using genetic algorithms：Inversion of plane－wave seismograms.Geophysics，56（11）：1794－1810.

Wang C L，Cheng J B，Yin C，et al.2013.Microseismic events location of surface and borehole observation with reverse－time focusing using interferometry technique.Chinese Journal Geophysics（in Chinese），56（9）：3184－3196，doi：10.6038／cjg20130931

Wang S S，Deng Y Q.1997.Hybrid algorithm for forward modeling of elastic waves.Oil Geophysical Prospecting （in Chinese），32（6）：804－808.

Wang W M，Hao J L，Yao Z X.2013.Preliminary result for rupture process of Apr.20，2013，Lushan Earthquake，Sichuan，China.Chinese Journal of Geophysics（in Chinese），56（4）：1412－1417，doi：10.6038／cjg20130436.

Yao Z X，Zheng T Y，Wen L X.1994.Moment tensor inversion method for determining the earthquake process by the use of P－wave form date.Chinese Journal of Geophysics（in Chinese），37（1）：34－44.

Ye G X，Jiang F X，Yang S H.2008.Possibility of automatically picking first arrival of microseismic wave by energy eigenvalue method.Chinese Journal of Geophysics（in Chinese），51（5）：1574－1581，doi：10.3321／j.issn：0001－5733.2008.05.033.

Zhang P，Yang C C，Li Y M.2000.3－D multi－valued seismic wave field reconstruction and Green′s function computation.Chinese Journal of Geophysics （in Chinese），43（2）：241－250，doi：10.3321／j.issn：0001－5733.2000.02.012.

Zhu G M，Li G H，Zhang W B.2008.Numeric simulation of 3－Component cross－hole seismic surveying wave field in VTI media.Oil Geophysical Prospecting （in Chinese），43（2）：206.

附中文参考文献

符力耘，牟永光.1994.弹性波边界元法正演模拟.地球物理学报，37（4）：521－529.

姜福兴，杨淑华，成云海等.2006.煤矿冲击地压的微地震监测研究.地球物理学报，49（5）：1511－1516.

靳平，徐果明，楼沩涛.1998.点力源在横向各向同性介质中激发的弹性波.地球物理学报，41（4）：525－536.

李幼铭.1988.透射波波场重建和介质层厚反演研究.地球物理学报，31（6）：649－656.

刘恩儒，岳建华，潘冬明.2006.地震各向异性——多组裂隙对横波偏振的影响.地球物理学报，49（5）：1401－1409.

刘培洵，陈顺云，郭彦双等.2014.声发射矩张量反演.地球物理学报，57（3）：858－866，doi：10.6038／cjg20140315.

王晨龙，程玖兵，尹陈等.2013.地面与井中观测条件下的微地震干涉逆时定位算法.地球物理学报，56（9）：3184－3196，doi：10.6038／cjg20130931.

王山山，邓玉琼.1997.弹性波正演模拟的混合算法.石油地球物理勘探，32（6）：804－808.

王卫民，郝金来，姚振兴.2013.2013年4月20日四川芦山地震震源破裂过程反演初步结果.地球物理学报，56（4）：1412－1417，doi：10.6038／cjg20130436.

姚振兴，郑天愉，温联星.1994.用P波波形资料反演中强地震地震矩张量的方法.地球物理学报，37（1）：34－44.

叶根喜，姜福兴，杨淑华.2008.时窗能量特征法拾取微地震波初始到时的可行性研究.地球物理学报，51（5）：1574－1581，doi：10.3321／j.issn：0001－5733.2008.05.033.

张钋，杨长春，李幼铭.2000.三维多值走时地震波场重建及格林函数计算.地球物理学报，43（2）：241－250，doi：10.3321／j.issn：0001－5733.2000.02.012.

朱光明，李桂花，张文波.2008.VTI介质三分量井间地震观测波场数值模拟.石油地球物理勘探，43（2）：201－207.