石油地球物理勘探  2020, Vol. 55 Issue (6): 1292-1304  DOI: 10.13810/j.cnki.issn.1000-7210.2020.06.015
0
文章快速检索     高级检索

引用本文 

吴如悦, 李晗, 常旭, 严红勇, 王一博. 页岩各向异性对被动源微地震波场的影响. 石油地球物理勘探, 2020, 55(6): 1292-1304. DOI: 10.13810/j.cnki.issn.1000-7210.2020.06.015.
WU Ruyue, LI Han, CHANG Xu, YAN Hongyong, WANG Yibo. Impact of shale anisotropy on the microseismic wavefield of passive sources. Oil Geophysical Prospecting, 2020, 55(6): 1292-1304. DOI: 10.13810/j.cnki.issn.1000-7210.2020.06.015.

本项研究受国家自然科学基金青年科学基金项目“微地震震源机制与页岩储层衰减系数(Q值)联合反演)”(41804050)和国家科技重大专项课题“陆相页岩油甜点地球物理识别与预测方法”(2017ZX05049-002)联合资助

作者简介

吴如悦  硕士研究生, 1996年生; 2018年毕业于吉林大学地球探测科学与技术学院, 获地球物理学专业学士学位; 现于中国科学院地质与地球物理研究所攻读硕士学位。主要从事被动源微地震各向异性数值模拟方面的学习和研究

常旭, 北京市朝阳区北土城西路19号中国科学院地质与地球物理研究所, 100029。Email:changxu@mail.iggcas.ac.cn

文章历史

本文于2020年1月22日收到,最终修改稿于同年9月30日收到
页岩各向异性对被动源微地震波场的影响
吴如悦①②③ , 李晗①② , 常旭①② , 严红勇①② , 王一博①②     
1 中国科学院地质与地球物理研究所中国科学院油气资源研究重点实验室, 北京 100029;
2 中国科学院地球科学研究院, 北京 100029;
3 中国科学院大学, 北京 100049
摘要:页岩储层具有较强的各向异性特征,忽略页岩各向异性对被动源微地震波(水力压裂产生)传播的影响,可能会降低储层表征结果的可信度。为此,提出一种定量评估地震波场旅行时和振幅偏差的方法,并基于理论数值算例,研究页岩各向异性对被动源微地震波场的影响。首先,利用矩张量模型描述爆炸型和剪切型微地震震源,并基于三维各向异性介质弹性波动方程和交错网格有限差分算法,正演模拟被动源的微地震波场;然后,对比不同震源、不同程度各向异性介质的理论波场,计算其P波和S波相对于各向同性介质波场的旅行时及振幅偏差,定性分析和定量评估不同震源机制下页岩各向异性对微地震波场特征的具体影响。数值算例的结果表明:对于所采用的地面微地震监测系统和平层VTI介质模型,页岩各向异性对微地震波场旅行时和振幅的影响均较为显著,旅行时偏差可达子波周期的30%,振幅相对偏差则高达60%,且旅行时偏差随炮检距增加而增大。此外,不同震源机制下,微地震波场特征受页岩各向异性影响的规律相似。
关键词微地震    页岩各向异性    震源机制    直达波    振幅    旅行时    数值模拟    
Impact of shale anisotropy on the microseismic wavefield of passive sources
WU Ruyue①②③ , LI Han①② , CHANG Xu①② , YAN Hongyong①② , WANG Yibo①②     
1 Key Laboratory of Petroleum Resource Research, Institute of Geology and Geophysics, Chinese Academy of Sciences, Beijing 100029;
2 Institution of Earth Science, Chinese Academy of Sciences, Beijing 100029;
3 University of Chinese Academy of Sciences, Beijing 100049
Abstract: Shale reservoir is featured by strong anisotropy. Ignoring the impact of shale anisotropy on the microseismic wavefield of passive sources (induced by fracturing stimulation) could affect the credibility of reservoir characterization. We proposes a method for quantitatively evaluating the deviation in travel time and amplitude of seismic wavefields, and study the influence of shale anisotropy on the microseismic wavefield of passive sources based on synthetic numerical test. Firstly, we use moment tensor to define explosive and shear microseismic sources, and caculate the microseismic wavefield of passive sources based on the elastic wave equation for three-dimensional anisotropic media and the staggered-grid finite-difference algorithm. Secondly, by comparing theoretical wavefields induced by different sources in different anisotropy media, we calculate the deviation in travel time and amplitude of P wave and S wave relative to isotropic models, and quantitatively and qualitatively analyze the impact of shale anisotropy on the microseismic wavefields under different source mechanisms. Numerical test show that for a surface microseismic monitoring system and a horizontally layered VTI medium model, shale anisotropy has significant influence on the travel time and amplitude. The travel time can deviate by 30%, and the deviation of amplitude can be up to 60%. The deviation of the travel time increases with increasing offset. In addition, under different source mechanisms, the impact of shale anisotropy on the microseismic wavefield are similar to each other.
Keywords: microseismic    shale anisotropy    source mechanism    direct wave    amplitude    travel time    numerical test    
0 引言

非常规油气,特别是页岩油气的勘探开发在美国获得巨大成功,成为改变化石能源短缺现状的重要后备资源,受到社会广泛重视[1-2]。页岩是具有大量片状或片状层理分布的沉积岩,岩石物理学研究表明,页岩具有不可忽视的各向异性特征(不同方向的速度差异可达35%)[3-4]。目前开发的大多数页岩油气储层具有明显的各向异性特征[5]

页岩油气的开采需要实施水力压裂改造储层。储层压裂产生微地震,受储层强各向异性的影响,这种被动源的微地震波在页岩中传播时,信号的相应特征(包括传播时间、振幅、相位和频率等)可能会发生改变。波场特征信息的变化,会直接影响微地震监测(定位成像、震源机制反演等)的结果[6-8],进而影响对储层的评价。因此,定量评估介质各向异性对被动源微地震波场的影响,对非常规油气的勘探开发有着重要意义。Carcione等[9]研究了介质各向异性对反射震相(PP和PS)振幅的影响,认为在AVO分析中考虑各向异性是非常必要的。被动源微地震监测的研究表明,如果忽略页岩各向异性,会导致微地震事件的定位出现偏差,因此需对各向异性造成的影响以及各向异性参数的反演开展进一步的研究[7]。Tsvankin等[10]总结了各向异性介质中地震波数值模拟、数据处理及反演的最新进展,认为各向异性的影响对微地震波(特别是剪切波)的影响是不可忽视的。通过基于VTI声波方程的微地震波场正演研究,孙伟家等[11]发现,与各向同性介质相比,页岩的强各向异性会导致地震波的旅行时和振幅出现显著偏差。Li等[12]研究了页岩各向异性引起的反射地震波振幅及相位偏差,定量分析了不同程度各向异性对反射波波场特征的影响。

上述关于页岩各向异性对地震波影响的研究,更多是利用声波方程,或者是利用弹性波方程研究反射波。基于三维弹性波动方程、对被动源微地震波场的研究相对较少。相比声波,弹性波更接近真实的地震波传播特性。本文基于三维介质弹性波动方程,引入Thomsen参数[13]表征页岩各向异性,并采用交错网格有限差分算法[14-19],正演模拟不同各向异性模型、不同震源机制[20]的微地震波场记录;通过计算不同波场记录的旅行时偏差和振幅相对偏差,定性分析和定量评估页岩各向异性对被动源微地震波场特征的影响。

1 方法原理 1.1 三维各向异性介质微地震波正演理论

一般认为,弹性波模拟更符合实际地层中地震波的传播。三维介质弹性波波动方程的表达式为

$ \rho \frac{{{\partial ^2}{u_i}}}{{\partial {t^2}}} = \frac{\partial }{{\partial {x_j}}}\left( {{c_{ijkl}}\frac{{\partial {u_k}}}{{\partial {x_l}}}} \right) + \rho {F_i} $ (1)

式中:u代表质点的位移场;ρF分别表示介质的密度和所受外力;xt分别表示空间和时间变量;cijkl表示4阶Hooke定律中的刚度矩阵元素,描述了应力与应变的对应关系,下标i, j, k, l=1、2、3,分别表示空间的xyz方向。由于弹性系数张量的对称性,刚度矩阵中一共只有21个独立变量。如果是VTI介质,弹性刚度矩阵可简化为

$ \left( {\begin{array}{*{20}{c}} {{c_{11}}}&{{c_{12}}}&{{c_{13}}}&0&0&0\\ {{c_{12}}}&{{c_{11}}}&{{c_{13}}}&0&0&0\\ {{c_{13}}}&{{c_{13}}}&{{c_{33}}}&0&0&0\\ 0&0&0&{{c_{44}}}&0&0\\ 0&0&0&0&{{c_{44}}}&0\\ 0&0&0&0&0&{{c_{66}}} \end{array}} \right) $ (2)

式(2)中的弹性参数无明确的物理意义。因此引入Thomsen参数[13]

$ \left\{ {\begin{array}{*{20}{l}} {\varepsilon = \frac{{{c_{11}} - {c_{33}}}}{{2{c_{33}}}}}\\ {\gamma = \frac{{{c_{66}} - {c_{44}}}}{{2{c_{44}}}}}\\ {\delta = \frac{{{{({c_{13}} + {c_{44}})}^2} - {{({c_{33}} - {c_{44}})}^2}}}{{2{c_{33}}({c_{33}} - {c_{44}})}}} \end{array}} \right. $ (3a)
$ \left\{ {\begin{array}{*{20}{l}} {{V_{{\rm{Pz}}}} = \sqrt {\frac{{{c_{33}}}}{\rho }} }\\ {{V_{{\rm{Sz}}}} = \sqrt {\frac{{{c_{44}}}}{\rho }} }\\ {{V_{{\rm{Pn}}}} = {V_{{\rm{Pz}}}}\sqrt {1 + 2\delta } } \end{array}} \right. $ (3b)

式中:ε表征纵波各向异性强度;γ表征横波各向异性的强度,同时也表征横波分裂的强度;δ为变异系数,表征地层纵波速度非对称轴方向的变化;VPzVSz分别代表纵波与横波的相速度;VPn表示地层内纵波的NMO速度。由式(3)可得[18]

$ \left\{ {\begin{array}{*{20}{l}} {{c_{11}} = \rho (1 + 2\varepsilon )V_{{\rm{Pz}}}^2}\\ {{c_{33}} = \rho V_{{\rm{Pz}}}^2}\\ {{c_{44}} = \rho V_{{\rm{Sz}}}^2}\\ {{c_{12}} = \rho V_{{\rm{Pz}}}^2[1 + 2\varepsilon - 2(1 - {f^\prime })(1 + 2\gamma )]}\\ {{c_{13}} = \rho V_{{\rm{Pz}}}^2\sqrt {{f^\prime }({f^\prime } + 2\delta )} - \rho V_{{\rm{Sz}}}^2}\\ {{c_{66}} = \rho (1 + 2\gamma )V_{{\rm{Sz}}}^2 = \frac{1}{2}({c_{11}} - {c_{12}})} \end{array}} \right. $ (4)

其中f′=$1-\frac{V_{\mathrm{Sz}}^{2}}{V_{\mathrm{Pz}}^{2}}=\frac{c_{33}-c_{44}}{c_{33}}$表示横纵波速比。

由式(4)可看出,VTI介质的弹性参数矩阵一共有6个参数,其中有一个(c66)是非独立的,所以只需要5个Thomsen参数就可以完全描述VTI介质。

1.2 三维微地震波交错网格有限差分算法

由式(1)~式(4)建立了VTI介质的弹性波正演理论。具体模拟时,需将波场分量沿时间和空间方向离散化到相应数值网格上。本文采用经典的交错网格有限差分法[5, 21-24],该方法具有高效率、低内存需求的特点,已被广泛用于勘探地球物理的地震波场模拟研究中。

交错网格有限差分算法原理为,计算一组空间点的压力及另一组空间点的速度[21],用网格点上的差商代替偏微分方程的微商建立差分方程,通过解方程组的形式得到波动方程的数值解。交错网格方法的特点是速度独立于应力而更新,因而数值模拟效率更高效[22]

弹性波波动方程(式(1))用速度—应力表示为

$ \rho (\mathit{\boldsymbol{x}})\frac{{\partial {v_i}}}{{\partial t}}(x,t) = \sum\limits_{j = 1}^3 {\frac{{\partial {T_{ij}}}}{{\partial {x_j}}}} (x,t) + \rho (x){F_i}(x,t) $ (5)

式中vi(x, t)、Tij(x, t)分别表示质点速度和应力张量。VTI介质中应力—应变关系式为

$ \left[ {\begin{array}{*{20}{c}} {\frac{{\partial {T_{11}}}}{{\partial t}}}\\ {\frac{{\partial {T_{22}}}}{{\partial t}}}\\ {\frac{{\partial {T_{33}}}}{{\partial t}}}\\ {\frac{{\partial {T_{23}}}}{{\partial t}}}\\ {\frac{{\partial {T_{13}}}}{{\partial t}}}\\ {\frac{{\partial {T_{12}}}}{{\partial t}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{c_{11}}}&{{c_{12}}}&{{c_{13}}}&0&0&0\\ {{c_{12}}}&{{c_{11}}}&{{c_{13}}}&0&0&0\\ {{c_{13}}}&{{c_{13}}}&{{c_{33}}}&0&0&0\\ 0&0&0&{{c_{44}}}&0&0\\ 0&0&0&0&{{c_{44}}}&0\\ 0&0&0&0&0&{{c_{66}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {\frac{{\partial {v_1}}}{{\partial {x_1}}}}\\ {\frac{{\partial {v_2}}}{{\partial {x_2}}}}\\ {\frac{{\partial {v_3}}}{{\partial {x_3}}}}\\ {\frac{{\partial {v_2}}}{{\partial {x_3}}} + \frac{{\partial {v_3}}}{{\partial {x_2}}}}\\ {\frac{{\partial {v_1}}}{{\partial {x_3}}} + \frac{{\partial {v_3}}}{{\partial {x_1}}}}\\ {\frac{{\partial {v_1}}}{{\partial {x_2}}} + \frac{{\partial {v_2}}}{{\partial {x_1}}}} \end{array}} \right] $ (6)

其中刚度矩阵可由式(4)给出。

交错网格有限差分算法中偏导数可基于以下原理展开

$ \begin{array}{l} \frac{{\partial f}}{{\partial a}}(a) \approx \frac{1}{{\Delta a}}\sum\limits_{m = 1}^N {{k_m}} \times \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left[ {f\left( {a + \frac{{2m - 1}}{2}\Delta a} \right) - f\left( {a - \frac{{2m - 1}}{2}\Delta a} \right)} \right] \end{array} $ (7)

式中:km是有限差分系数(通过求解范德蒙矩阵来计算[23]);m=1, 2, …, N, 差分阶数为2N(在本文的正演模拟中,时间采用二阶,空间使用十阶)。

根据式(7),式(5)可以离散化为(这里示例给出v1)

$ \begin{array}{*{20}{l}} \begin{array}{l} {v_1}\left( {{x_1},{x_2},{x_3},t + \Delta t} \right)\\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} = {v_1}\left( {{x_1},{x_2},{x_3},t} \right) + {F_1}\Delta t + \frac{{\Delta t}}{{\Delta x}}\frac{1}{\rho }\sum\limits_{m = 1}^N {{k_m}} \times \\ \qquad \begin{array}{*{20}{l}} {\left[ {{T_{11}}\left( {{x_1} + \frac{{2m - 1}}{2}\Delta x,{x_2},{x_3},t + \frac{1}{2}\Delta t} \right) - } \right.}\\ {{T_{11}}\left( {{x_1} - \frac{{2m - 1}}{2}\Delta x,{x_2},{x_3},t + \frac{1}{2}\Delta t} \right) + }\\ {{T_{12}}\left( {{x_1},{x_2} + \frac{{2m - 1}}{2}\Delta x,{x_3},t + \frac{1}{2}\Delta t} \right) - }\\ {{T_{12}}\left( {{x_1},{x_2} - \frac{{2m - 1}}{2}\Delta x,{x_3},t + \frac{1}{2}\Delta t} \right) + }\\ {{T_{13}}\left( {{x_1},{x_2},{x_3} + \frac{{2m - 1}}{2}\Delta x,t + \frac{1}{2}\Delta t} \right) - }\\ {\left. {{T_{13}}\left( {{x_1},{x_2},{x_3} - \frac{{2m - 1}}{2}\Delta x,t + \frac{1}{2}\Delta t} \right)} \right]} \end{array} \end{array} \end{array} $ (8)

根据式(7),式(6)可以离散化为(这里示例给出T11)

$ \begin{array}{l} \begin{array}{*{20}{l}} {{T_{11}}({x_1},{x_2},{x_3},t + \Delta t)}\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} = {T_{11}}({x_1},{x_2},{x_3},t) + \frac{{\Delta t}}{{\Delta x}}\sum\limits_{m = 1}^N {{k_m}} \times } \end{array}\\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \begin{array}{*{20}{c}} {\left\{ {{c_{11}}\left[ {{v_1}\left( {{x_1} + \frac{{2m - 1}}{2}\Delta x,{x_2},{x_3},t + \frac{1}{2}\Delta t} \right) - } \right.} \right.}\\ {\left. {{v_1}\left( {{x_1} - \frac{{2m - 1}}{2}\Delta x,{x_2},{x_3},t + \frac{1}{2}\Delta t} \right)} \right] + }\\ {{c_{12}}\left[ {{v_2}\left( {{x_1},{x_2} + \frac{{2m - 1}}{2}\Delta x,{x_3},t + \frac{1}{2}\Delta t} \right) - } \right.}\\ {\left. {{v_2}\left( {{x_1},{x_2} - \frac{{2m - 1}}{2}\Delta x,{x_3},t + \frac{1}{2}\Delta t} \right)} \right] + }\\ {{c_{13}}\left[ {{v_3}\left( {{x_1},{x_2},{x_3} + \frac{{2m - 1}}{2}\Delta x,t + \frac{1}{2}\Delta t} \right) - } \right.}\\ {\left. {\left. {{v_3}\left( {{x_1},{x_2},{x_3} - \frac{{2m - 1}}{2}\Delta x,t + \frac{1}{2}\Delta t} \right)} \right]} \right\}} \end{array} \end{array} $ (9)

式中的Δx和Δt分别是空间和时间步长。交错网格有限差分算法中,速度场和应力场在时间和空间上是交错的[5, 22, 24]

本文正演模拟的稳定性条件为

$ \left\{ {\begin{array}{*{20}{l}} {\Delta t \le \frac{1}{{\sqrt 3 }} \cdot \frac{{\Delta x}}{{{v_{\max }}}}}\\ {\Delta x \le \frac{{{\lambda _{\min }}}}{5}} \end{array}} \right. $ (10)

式中:vmax是介质中的最大速度;λmin是地震波场的最小波长。此外,在处理边界问题时,使用了吸收边界条件,即对边界区域应用以下衰减函数

$ A(i) = {\rm{exp}}\left[ { - d{{\left( {\frac{i}{M}} \right)}^2}} \right] $ (11)

式中:Mi分别为是衰减区域中网格点总数和序号; d是衰减系数(本文中,M=100,d=0.1)。

1.3 被动源微地震波场受各向异性影响的评估

微地震信号易受噪声干扰,实际一般只能清晰识别出直达波震相。据此,本文主要研究各向异性对直达纵波(P波)、直达横波(S波)的影响。考虑到Thomsen参数中εδ是纵波各向异性强度的度量,γ表征横波各向异性强度,本文的研究思路为:首先基于均匀及层状介质模型,建立符合页岩特征的不同程度各向异性介质模型,分别分析εδ对直达P波、γ对直达S波旅行时和振幅的影响。

计算各向异性模型与各向同性模型地震记录的互相关函数,表征各向异性对微地震波旅行时的影响,即

$ {r_{i,j}}(k) = \frac{{\sum\limits_{n = 1}^N {{D_i}} (n){D_j}(n + k)}}{{\sqrt {\sum\limits_{n = 1}^N {D_i^2} (n)\sum\limits_{n = 1}^N {D_j^2} (n)} }} $ (12)

式中:Di表示选定的参考道,即各向同性模型的地震波记录;Dj为加入各向异性后的地震记录。假设互相关函数ri, j(k)的最大值对应的时刻为Δti, j,则表明当Dj平移Δti, j时,DiDj波形相似度最大,故Δti, j即为求得的各向异性模型与各向同性模型地震记录之间的旅行时偏差。

表征各向异性对微地震波振幅的影响时,使用每个检波器地震信号的最大振幅(含正负,符号表示地震波的极性),以及各向异性与各向同性模型微地震信号的最大振幅偏差百分比(Maximum Amplitude Difference Percentage,MADP)

$ {\rm{MADP}} = \left| {\frac{{{A_{{\rm{ani}}}}(n) - {A_{{\rm{iso}}}}(n)}}{{{A_{{\rm{iso}}}}(n)}}} \right| \times 100\% $ (13)

式中:Aani(n)为各向异性介质模型第n个检波器中的最大振幅(含正负);Aiso(n)为各向同性介质模型第n个检波器中的最大振幅。MADP可定量描述各向异性对微地震信号振幅的影响程度。

页岩各向异性对被动源微地震波场的影响评估按如下步骤进行:

(1) 根据页岩储层的地质特征,建立均匀各向异性介质模型M0、ME、MD和MG(具体参数见表 1),各向异性参数的选择借鉴了其他页岩各向异性研究成果[25-26],基于控制变量法分析不同各向异性参数的影响;

表 1 均匀模型的各向异性参数

(2) 使用交错网格有限差分正演方法,分别模拟模型M0、ME、MD和MG的地震波场;

(3) 计算模型ME、MD和MG相对于模型M0的直达P波、直达S波旅行时偏差和最大振幅偏差。

2 数值模拟分析

分别针对均匀和层状VTI介质,建立不同各向异性参数模型,分析不同各向异性对微地震波场特征的影响。

2.1 均匀介质模型

设计一个网格节点数为201×201×201的均匀模型:网格大小为10m×10m×10m,VP0为3500m/s,VS0为2020m/s,密度为2.3g/cm3。震源采用主频为20Hz的雷克子波,置于模型内部(1000m, 1000m, 1250m)。为了模拟不同震源机制的影响,设立一组爆炸源(简写为ISO源)及两组剪切源(分别简写为DC1、DC2震源),矩张量及机制球表现形式见表 2。模型介质各向异性参数同表 1。在地表平面内以震中为中心,沿xy方向分别布设201个检波器,检波器间距均为10m,沿x方向检波器序号为1~201,沿y方向检波器序号为202~402(图 1)。时间采样间隔为0.5ms,采样时长为1.25s。

表 2 三种不同机制的微地震震源

图 1 均匀模型、微地震震源及观测系统示意图
2.1.1 各向异性对旅行时的影响

考虑到地震勘探通常使用单分量(一般为z分量)检波器,本文主要分析z分量微地震波场。图 2t=0.3s时三种震源机制x-z切面的波场快照;图 3为各向异性模型ME相对于各向同性模型M0z分量初至P波旅行时的偏差图形;图 4为不同各向异性模型地面接收波形图。分析图 2~图 4,可以看出:

图 2 不同震源、不同各向异性均匀介质模型x-z切面波场快照(t=0.3s) (a)模型M0; (b)模型ME;(c)模型MD;(d)模型MG。图中从上到下的震源分别为ISO、DC1、DC2

图 3 不同震源机制均匀介质模型ME与M0初至P波z分量的旅行时偏差 (a)ISO震源;(b)DC1震源;(c)DC2震源

图 4 不同震源、不同各向异性均匀介质模型的微地震波形记录 图中从上到下分别对应震源ISO、DC1、DC2

(1) 各向异性对微地震波旅行时的影响非常显著。在图 3中,旅行时偏差最大为0.015s(30个采样点),相当于子波周期(0.05s)的30%。

(2) 各向异性影响直达P波在不同方向上的传播。同一震源机制下,不同各向异性介质模型微地震波z分量初至P波相位不同(图 4箭头处)。此外,ME模型的波场快照(图 2b)中P波在垂直各向异性对称轴方向(x方向)上传播最快,在平行各向异性对称轴方向(z方向)上传播最慢。Thomsen各向异性假设下,纵、横波速度可表示为

$ \left\{ {\begin{array}{*{20}{l}} {{V_{\rm{P}}}(\theta ) = {V_{{\rm{Pz}}}}(1 + \delta {{\sin }^2}\theta {{\cos }^2}\theta + \varepsilon {{\sin }^4}\theta )}\\ {{V_{{\rm{SV}}}}(\theta ) = {V_{{\rm{Sz}}}}\left[ {1 + \frac{{V_{{\rm{Pz}}}^2}}{{V_{{\rm{Sz}}}^2}}(\varepsilon - \delta ){{\sin }^2}\theta {{\cos }^2}\theta } \right]}\\ {{V_{{\rm{SH}}}}(\theta ) = {V_{{\rm{Sz}}}}(1 + \gamma {{\sin }^2}\theta )} \end{array}} \right. $ (14)

式中θ表示波前法向与垂直方向的夹角。沿z方向传播时,θ=0°, VP(0°)=VPz; 沿x方向传播,θ=90°, VP(90°)=VPz(1+ε)。本文基于页岩储层实际情况,ME模型ε取值为0.25,所以P波在x方向上传播更快。

(3) 各向异性对直达S波的影响与P波不同。从图 2中可以看到,震源机制相同、不同各向异性介质模型(每一行)得到的z分量直达S波,在平行和垂直各向异性对称轴上旅行时相同。通过分析检波器三分量记录的直达S波属性可知,z分量对应SV震相(表 3所示)。对于M0、ME、MD和MG四种模型,沿z方向传播时,θ=0°, VSV(0°)=VSz;沿x方向传播时,θ=90°, VSV(90°)=VSz,因而在这两个方向上的SV波传播速率相同。

表 3 检波器三分量记录直达S波属性

(4) 各向异性对微地震波旅行时的影响与炮检距有关。从图 3可看出,同一震源机制下,随着炮检距的增加,ME与M0的地震波旅行时偏差随之增大。由于炮检距增大,地震波的传播路径变长,各向异性的影响累积增加;其次,随着炮检距增大,检波器接收到地震波场的水平分量能量增强,导致总波场的能量分配从垂直分量占绝对优势变为垂直、水平分量均等,考虑到各向异性导致P波水平方向速度大于垂直方向,故由各向异性造成的旅行时偏差增大。

(5) 各向异性对微地震波旅行时的影响与震源机制无关。图 2中,模型相同、震源机制不同时(每一列),波场旅行时无明显差异。从图 3也可看出震源机制对微地震波旅行时偏差无影响。地震波旅行时仅由传播距离及速度决定,传播距离一定,速度由式(14)确定,因此震源机制对微地震波旅行时无关。

2.1.2 各向异性对振幅的影响

图 2中同一震源机制(每一行)均采用相同增益,通过横向对比可以看出,各向异性模型的微地震波振幅明显不同于各向同性模型。

首先分析各向异性对直达P波振幅的影响。图 5为不同震源(ISO、DC1、DC2)模型ME、MD和M0z分量初至P波最大振幅及相应偏差(MADP)分布图(沿y方向测线,序列号202~402),从中可以看出:

图 5 不同震源机制不同各向异性均匀介质模型z分量初至P波振幅对比 (a)最大绝对振幅;(b)MADP值。图中从左到右震源分别对应ISO、DC1、DC2

(1) εδ对微地震P波振幅的影响较为显著。图 5b中,DC1震源激发,ME和MD的MADP最高可达60%,说明考虑各向异性影响后,微地震波振幅最高可增大0.6倍。

(2) 从爆炸源的最大绝对振幅(图 5a左)发现,M0、ME、MD三组模型的最大绝对振幅均沿震源左右对称;剪切源的最大绝对振幅(图 5a中、图 5a右)相对于震源中心对称,且极性出现反转。剪切源与爆炸源最大绝对振幅分布的明显不同是由震源机制造成的。此外,图 2中不同方向微地震波的振幅不同,沿各向异性对称轴方向(z方向)振幅最强,垂直方向(x方向)振幅最弱,这是因为沿x方向微地震波传播速度最快,传播路径最远,因而能量衰减最多,振幅最弱。这与波的空间传播效应、各向异性导致的垂直和水平能量分配不一致有关。

(3) 不同各向异性参数对微地震波场的影响程度不同。从图 5b可看出,在近炮检距处(检波器序号250~350之间,对应炮检距约小于1350m),δ对微地震波振幅影响明显强于ε;远炮检距处(检波器序号202~250和350~402之间,对应炮检距大于1350m),ε对微地震波振幅的影响相对较强。这表明εδ对振幅的影响程度随炮检距发生变化。

(4) 不同震源机制各向异性对微地震波振幅的影响规律可能相似。从图 5b可看出,不同震源机制ME、MD相较于M0的MADP数值不同,表明各向异性对振幅能量的影响存在差异。但是,图 5b中MADP随炮检距的变化曲线形态相似,这是因为MADP反映的是相对振幅偏差,在一定程度上减弱了震源机制的影响。

爆炸源在各向异性介质中不产生S波,故本文仅分析剪切震源各向异性对微地震直达S波振幅的影响。图 6为剪切源MG和M0模型的z分量直达S波最大振幅分布图(沿y方向测线)。从图 6可以看出,MG与M0的最大绝对振幅图形重合,偏差为零。这是因为直达S波z分量为SV震相(表 2),VSV不受γ影响(式(14)可以证明), 即γ对微地震直达S波z分量振幅无影响。

图 6 均匀介质模型不同震源初至S波SV分量最大绝对振幅 (a)DC1震源; (b)DC2震源。MG(蓝色)与M0(红色)最大绝对振幅重合

根据表 2可知,y方向检波器的直达S波SH震相分布在水平分量,考虑地震勘探通常使用单分量(一般为z分量)检波器,不常使用水平分量, 故本文暂不分析SH震相。关于γ对直达S波SH震相振幅的影响可作为后续研究方向。

2.2 层状VTI介质模型

设计的层状VTI介质模型如图 7所示。模型网格节点数为201×201×201,网格大小为10m×10m×10m。震源采用主频为20Hz的雷克子波,位于(1000m, 1000m, 1250m)处。模型中间层为各向异性页岩层,各向异性参数同表 1,其他参数见表 4。观测系统同图 1:在地表过震中布设两条垂直交叉测线(分别沿xy方向,图 7),每条测线布有201个检波器,间隔10m。基于三组震源(ISO、DC1、DC2震源,表 1),对四组介质模型分别正演模拟了402个检波器的地震波形记录;时间采样间隔为0.5ms,采样时长为1.25s。

图 7 层状VTI介质模型及震源、观测系统示意图

表 4 层状VTI页岩模型的介质参数
2.2.1 各向异性对地震波场旅行时影响

图 8t=0.3s时不同震源机制x-z切面的微地震波场快照;图 9为不同震源机制ME与M0z分量初至P波、S波旅行时偏差图;图 10为不同震源、不同各向异性模型的波形记录图。分析图 8~图 10,可以看出:

图 8 层状VTI介质模型不同震源、不同各向异性x-z切面的波场快照(t=0.3s) (a)模型M0;(b)模型ME;(c)模型MD;(d)模型MG。从上到下震源分别为ISO、DC1、DC2

图 9 层状VTI介质模型ME与M0不同震源机制z分量初至波旅行时偏差图 (a)ISO震源;(b)DC1震源;(c)DC2震源。上图为初至P波旅行时偏差,下图为初至S波旅行时偏差

图 10 不同震源、不同各向异性层状VTI介质模型的微地震波形记录 从上到下分别对应震源ISO、DC1、DC2

(1) 图 9中P波旅行时偏差最大为0.012s(相当于24个采样点,达子波周期(0.05s)的24%,S波旅行时偏差最大为0.015s(相当于30个采样点),达子波周期(0.05s)的30%,表明各向异性对微地震波旅行时的影响较为显著。

(2) 同一震源机制,不同各向异性介质z分量初至P波(图 10红色箭头)、初至S波(图 10蓝色箭头)相位均不同。此外, 从图 8b可以看出,ε>0时,P波在x方向上传播最快,在z方向上传播最慢。

(3) 同一震源机制,不同各向异性介质模型微地震波z分量的初至SV波,在xz方向上的旅行时相同。

(4) 随着炮检距的增加,地震波的传播路径变长,各向异性对微地震波旅行时的影响累积增加。

(5) 不同的震源机制,微地震波初至旅行时无明显区别(图 8图 9),各向异性对微地震波旅行时的影响与震源机制无关。

与均匀介质各向异性相比,层状VTI介质模型各向异性相与各向同性模型P波的旅行时偏差值变小。考虑到震源位置相同,检波器布设在地表,均匀介质(图 1)和层状VTI模型(图 7)两者的差异主要体现在建立各向异性模型的方式(均匀介质为全空间各向异性,VTI介质仅中间页岩层存在各向异性)不同。震源位于模型中心,因此两种模型微地震波的差异主要由第一层(图 7)介质参数(速度、密度、各向异性)的不同而引起。相较于均匀介质,层状VTI介质模型的总体各向异性程度减弱,因此VTI模型P波的旅行时偏差(子波周期的24%)小于均匀介质(子波周期的30%)。

2.2.2 各向异性对地震波场振幅影响

首先考虑εδ对微地震直达P波振幅的影响。图 11为不同震源ME、MD和M0模型z分量初至P波的最大绝对振幅值及相应偏差图(沿y方向测线,序列号202~402)。通过分析可得到以下认识:

图 11 不同震源、层状VTI介质模型初至P波z分量最大绝对振幅及偏差图 (a)最大绝对振幅;(b)MADP值。从左到右震源分别对应ISO、DC1、DC2

(1) 各向异性对微地震波振幅有明显影响。从图 11b可看到,DC1震源的MADP最高可达30%。

(2) 爆炸源的最大绝对振幅(图 11a)沿震源左右对称,随炮检距的增加而逐渐减小,在震源处最大;剪切源的最大绝对振幅(图 11a中、11a左)出现极性反转,沿震源中心对称。剪切源与爆炸源最大绝对振幅图形的不同,本文认为这与震源机制有关。

(3) 不同震源,ME和MD的MADP变化规律相似(图 11b),表明各向异性对振幅能量的影响规律相似,这也是由于由式(13)计算的MADP反映相对偏差在一定程度上减弱了震源机制的影响。

(4) 观察同一震源模型ME和MD的MADP(例如图 11b左)变化可知,εδ对振幅的影响程度随炮检距增加而发生变化。在近炮检距处(检波器序号260~340之间,对应炮检距小于1310m),δ对微地震波振幅的影响明显强于ε;远炮检距处(检波器序号202~260和340~402之间,对应炮检距大于1300m),ε对微地震波振幅的影响较强。

(5) 与均匀模型相对比,层状VTI介质模型MADP数值普遍较小,不超过30%。这是因为层状VTI介质模型各向异性仅存在于中间页岩层处,各向异性程度总体减弱,故MADP数值降低。

(6) 相比均匀模型,层状VTI模型εδ对振幅影响程度发生反转的位置不同,对应的炮检移由约1350m处减小到约1310m。这是因为检波器位于地表,接收到的两个模型微地震波信号特征的差异主要由第一层介质不同导致。

同前面均匀介质模型相似,DC1、DC2震源直达S波z分量(图 12)MG与M0得到的最大绝对振幅图形重合,偏差为零,说明微地震直达SV波(表 3)振幅不受γ影响。

图 12 层状VTI介质模型不同震源z分量初至S波最大绝对振幅 (a)DC1震源; (b)DC2震源。MG(蓝色)与M0(红色)最大绝对振幅重合
3 结论与建议

本文针对页岩储层各向异性对被动源微地震波场特征的影响展开研究。首先基于三维各向异性介质弹性波动方程和交错网格有限差分算法,正演模拟了被动源的微地震波场。通过对比不同震源、不同程度各向异性介质模型的理论波形记录,计算其旅行时、振幅相对偏差,定量分析了均匀及层状VTI模型中各向异性对微地震波场特征的影响。

根据实际微地震事件的记录特点,本文主要评估了垂直z分量微地震波受各向异性的影响。数值算例的结果表明,对于所采用的地面微地震监测系统及均匀(或平层状VTI)介质模型,得到如下结论:(1)各向异性的微地震记录与各向同性介质的波场记录相比,旅行时和振幅相对偏差均较为显著(旅行时偏差可达子波周期的30%,振幅相对偏差则高达60%);(2)各向异性对微地震波传播速度具有方向性的影响。加入各向异性后(ε>0),直达P波在垂直各向异性对称轴的方向上传播最快,平行对称轴方向上传播最慢,直达SV波在垂直、平行各向异性对称轴方向上传播速度相同;(3)炮检距越大,各向异性对波形的影响越显著;(4)不同各向异性参数对波场影响的程度不同,近炮检距处(炮检距约小于1310m),δ对微地震波振幅影响明显强于ε,远炮检距处(炮检距约大于1310m),ε对微地震波振幅的影响较强;γ对直达S波z分量(SV震相)的振幅无影响;(5)微地震波场特征随各向异性的变化规律与震源机制无关。

各向异性对地震记录的影响,由地层产状、各向异性对称轴方向及地震波传播路径等因素决定。进行页岩气开采和利用微地震信息进行储层评价时,需要充分考虑各向异性的影响。忽略各向异性,将影响微地震震源机制的正确判断,这直接影响储层评价工作的进行。关于各向异性对被动源微地震波场影响的研究仍存在很多问题,本文建立的模型较为简单,具有倾斜对称轴的TTI模型可能更符合实际情况。此外,各向异性对直达S波SH震相的影响也需进一步研究。本文为定量评估页岩各向异性的影响提供了一种研究方法和思路,在实际微地震传播特征研究中应具有良好应用前景。

参考文献
[1]
张金川, 徐波, 聂海宽, 等. 中国页岩气资源勘探潜力[J]. 天然气工业, 2008, 28(6): 136-141.
ZHANG Jinchuan, XU Bo, NIE Haikuan, et al. Exploration potential of shale gas resources in China[J]. Natural Gas Industry, 2008, 28(6): 136-140.
[2]
徐刚, 李德旗, 王适择, 等. 微地震监测在地震地质工程一体化中的应用[J]. 石油地球物理勘探, 2018, 53(增刊2): 115-123.
XU Gang, LI Deqi, WANG Shize, et al. Microseismic monitoring in seismic and geology integration[J]. Oil Geophysical Prospecting, 2018, 53(S2): 115-123.
[3]
Vernik L, Nur A. Ultrasonic velocity and anisotropy of hydrocarbon source rocks[J]. Geophysics, 1992, 57(5): 727-735. DOI:10.1190/1.1443286
[4]
Johnston J E, Christensen N I. Seismic anisotropy of shales[J]. Journal Geophysical Research, 1995, 100(B4): 5991-6003. DOI:10.1029/95JB00031
[5]
刘财, 邓馨卉, 郭智奇, 等. 基于岩石物理的页岩储层各向异性表征[J]. 石油地球物理勘探, 2018, 53(2): 339-346.
LIU Cai, DENG Xinhui, GUO Zhiqi, et al. Shale reservoir anisotropic characterization based on rock physics[J]. Oil Geophysical Prospecting, 2018, 53(2): 339-346.
[6]
Virieux J. SH-wave propagation in heterogeneous media:Velocity-stress finite-difference method[J]. Geophysics, 1984, 49(11): 1933-1942. DOI:10.1190/1.1441605
[7]
陈颙, 黄庭芳, 刘恩儒. 岩石物理学[M]. 安徽合肥: 中国科学技术大学出版社, 2009.
CHEN Yong, HUANG Tingfang, LIU Enru. Rock Physics[M]. Hefei, Anhui: Press of University of Science and Technology of China, 2009.
[8]
常旭, 王一博. 微地震反演研究[M]. 北京: 科学出版社, 2019.
CHANG Xu, WANG Yibo. Microseismic Inversion[M]. Beijing: Science Press, 2019.
[9]
Carcione, JoséM, Helle H B, et al. Effects of attenuation and anisotropy on reflection amplitude versus offset[J]. Geophysics, 1998, 63(5): 1652-1658. DOI:10.1190/1.1444461
[10]
Tsvankin I, Gaiser J, Grechka V, et al. Seismic aniso-tropy in exploration and reservoir characterization:An overview[J]. Geophysics, 2010, 75: A15-A29.
[11]
孙伟家, 符力耘, 管西竹, 等. 页岩气地震勘探中页岩各向异性的地震模拟研究[J]. 地球物理学报, 2013, 56(3): 961-970.
SUN Weijia, FU Liyun, GUAN Xizhu, et al. A study on anisotropy of shale using seismic forward modeling in shale gas exploration[J]. Chinese Journal of Geophysics, 2013, 56(3): 961-970.
[12]
Li H, Liu X, Chang X, et al. Impact of shale anisotropy on seismic wavefield[J]. Energies, 2019, 12(23): 4412. DOI:10.3390/en12234412
[13]
Thomsen L. Weak elastic anisotropy[J]. Geophysics, 1986, 51(10): 1954-1966. DOI:10.1190/1.1442051
[14]
Levander A R. Fourth-order finite-difference P-SV seismograms[J]. Geophysics, 1988, 53(11): 1425-1436. DOI:10.1190/1.1442422
[15]
Saenger E H, Gold N, Shapiro S A. Modeling the pro-pagation of elastic waves using a modified finite-difference grid[J]. Wave Motion, 2000, 31(1): 77-92. DOI:10.1016/S0165-2125(99)00023-2
[16]
Saenger E H, Bohlen T. Finite-difference modeling of viscoelastic and anisotropic wave propagation using the rotated staggered grid[J]. Geophysics, 2004, 69(2): 583-591. DOI:10.1190/1.1707078
[17]
吴国忱, 梁锴, 印兴耀. TTI介质弹性波相速度与偏振特征分析[J]. 地球物理学报, 2010, 53(8): 1914-1923.
WU Guochen, LIANG Kai, YIN Xingyao. The analysis of phase velocity and polarization feature for elastic wave in TTI media[J]. Chinese Journal of Geophy-sics, 2010, 53(8): 1914-1923.
[18]
王璐琛, 常旭, 王一博. TTI介质的交错网格伪P波正演方法[J]. 地球物理学报, 2016, 59(3): 1046-1058.
WANG Luchen, CHANG Xu, WANG Yibo. Forward modeling of pseudo P waves in TTI medium using staggered grid[J]. Chinese Journal of Geophysics, 2016, 59(3): 1046-1058.
[19]
姚振岸, 孙成禹, 唐杰, 等. 基于不同震源机制的黏弹各向异性微地震波场模拟[J]. 石油地球物理勘探, 2017, 52(1): 63-70.
YAO Zhen'an, SUN Chengyu, TANG Jie, et al. Micro-seismic forward modeling in viscoelastic anisotropic media based on different focal mechanisms[J]. Oil Geophysical Prospecting, 2017, 52(1): 63-70.
[20]
唐杰, 温雷, 王浩, 等. 正交各向异性介质中的剪张源震源机制与矩张量特征[J]. 石油地球物理勘探, 2018, 53(6): 1247-1255.
TANG Jie, WEN Lei, WANG Hao, et al. Focal mechanisms and moment tensor in orthorhombic anisotropic media[J]. Oil Geophysical Prospecting, 2018, 53(6): 1247-1255.
[21]
Danggo M Y, Mungkasi S. A staggered grid finitedifference method for solving the elastic wave equations[J]. Journal of Physics:Conference Series, 2017, 909: 012047. DOI:10.1088/1742-6596/909/1/012047
[22]
Graves R W. Simulating seismic wave propagation in 3D elastic media using staggered-grid finite diffe-rences[J]. Bulletin Seismological Socitey of America, 1996, 86(4): 1091-1106.
[23]
Liu Y, Sen M K. An implicit staggered-grid finite-difference method for seismic modelling[J]. Geophy-sical Journal International, 2009, 179(1): 459-474. DOI:10.1111/j.1365-246X.2009.04305.x
[24]
Bossy E, Talmant M, Laugier P. Three-dimensional simulations of ultrasonic axial transmission velocity measurement on cortical bone models[J]. J Acoust Soc Am, 2004, 115(5): 2314-2324. DOI:10.1121/1.1689960
[25]
Liu X, Guo Z, Liu C, et al. Anisotropy rock physics model for the Longmaxi shale gas reservoir, Sichuan Basin, China[J]. Applied Geophysics, 2017, 14(1): 21-30.
[26]
Qin Y, Zhang Z, Xu S. CDP mapping in tilted transversely isotropic (TTI) media.Part Ⅱ:Velocity analysis by combining CDP mapping with a genetic algorithm[J]. Geophysical Prospecting, 2003, 51(4): 325-332. DOI:10.1046/j.1365-2478.2003.00372.x