地面与井中联合观测的微地震定位和速度结构同时反演

黄国娇　巴　晶　钱　卫

(河海大学地球科学与工程学院地球探测研究所，江苏南京 211100)

1　引言

水力压裂技术是非常规油气资源(如致密油气、页岩油气)开发过程中重要的增产措施，该技术通过将流体高压注入地层，使岩石破裂以提高地层渗透率[1]。微地震监测技术通过观测、定位及分析微地震事件监测生产活动及地下状态，可用于水力压裂效果分析[2-6]，在储层开发监测中起到重要作用。目前，微震监测技术应用的主要困难在于： ①水力压裂诱导的微地震事件通常信噪比较低，且观测数据量大，初至旅行时难拾取； ②即使能够较准确地拾取初至旅行时(旅行时信息包含了地震事件的发震时刻和旅行时信息)，微地震事件的定位仍然存在误差，原因是缺乏精确的地下速度模型。微地震反演中采用的速度模型一般根据地震、声波测井、VSP等资料建立[6]，并采用压裂前的射孔资料校正速度模型[7,8]。但是，射孔反演速度模型需已知射孔位置和发震时刻，而射孔事件的准确发震时刻往往是未知的。再者，射孔数据一般在压裂前采集，而在水力压裂过程中，高压流体会产生裂缝，新增裂缝和孔隙压力的改变都会引起地层速度的变化[9]。因此如果采用射孔资料校正后的速度模型进行微地震定位，可能存在较大误差[10]。

针对微地震初至拾取困难的问题，Song等[11]建立了一个微地震强信号的波形数据库，通过分析强弱微地震信号的波形相关性提取弱微震事件的初至； 宋维琪等[12]提出了独立分量分析与压缩感知联合理论的微地震弱信号的提取方法； 王晨龙等[13]和李振春等[14]采用地震干涉逆时定位算法，不需要拾取初至信息； 王鹏等[15]提出基于熵值的微震事件识别和恢复方法，以解决微地震监测数据信噪比低的问题； 李稳等[16]将稀疏码收缩法和小波变换相结合处理井下微地震数据，以保证有效地震信号的识别和提取； 曹俊海等[17]采用多道匹配跟踪算法提高微地震信号的信噪比，进而提高微地震信号识别的准确性。

为得到较准确的速度模型，减小速度模型对定位结果的影响，Grechka等[18]同时反演了介质的各向异性参数和微地震位置，但假设地下为均匀介质； Li等[19]将双差定位法推广至同时确定水平层状各向异性参数及微地震定位； 谭玉阳等[20]提出基于初至旅行时差的微地震速度模型反演方法，可在射孔起始时间未知的情况下反演地层速度模型。

在井下微地震监测中，微地震事件与检波器之间的距离通常在几十至数百米的范围内，因此将地下介质假设为水平层状模型[7，21]，提出了地面与井中观测旅行时联合的微地震定位和速度结构同时反演方法： 首先推导了旅行时关于速度、界面深度以及震源参数的一阶偏导数，给出了Jacobi矩阵计算方法； 然后用最短路径算法[22]追踪射线路径和计算旅行时；再采用共轭梯度算法求解带约束的阻尼最小二乘问题，更新地下层状介质的结构(层速度和层界面深度)以及微地震定位(包括震源位置和发震时刻)； 最后讨论了方法的抗噪声能力。

2　方法原理

2.1　反演问题

一般而言，基于旅行时数据的非线性同时反演问题，可归结为下列带约束的阻尼最小二乘最优化问题，其反演问题的L2范数目标函数为

(1)

式中： δt=(δt1，δt2，…，δtM)是微地震事件的实际观测旅行时与理论数据的旅行时残差，M则为拾取的初至数量；δm=(λ1δV，λ2δD，λ3δX)，其中δV和δD分别是地下介质各层速度和地层界面深度的相对改变量，δX为微地震事件的震源坐标和发震时刻的相对改变量，λ1、λ2、λ3为各类参数的更新权重系数；A是Jacobi矩阵；μ是阻尼因子；a和b分别是反演解中更新各类参数更新量的上、下边界值。

式(1)解的一般形式[23]为

[μCm+ATCdA]Zmδm=ATCdδt

(2)

式中：Cm、Cd为模型和数据空间的协方差矩阵；Zm为分段约束算子，即

(3)

式(2)是非线性问题局部线性化的基本反演公式，其解具有局域解的特征。为了得到具有实际物理意义的解，可采用Cm、Cd及Zm的先验信息进行约束。因此，选择不同的组合(Cm、Cd、Zm)可得到不同形式的反演解。

本文采用Tarantola等[24]提出的广义带约束的阻尼最小二乘反演解，即Cm、Cd分别取模型和数据空间的协方差矩阵的逆，则解为

(4)

上述矩阵方程可基于共轭梯度法采用迭代算法求解[25]，其关键问题是如何求解具有偏导性质的Jacobi矩阵。

2.2　Jacobi矩阵元素计算

当地下介质结构(层速度和层深度)和微地震事件同时反演时，Jacobi偏导矩阵包含了三部分：一是旅行时关于速度变化的导数；二是旅行时关于层深度变化的导数； 三是旅行时关于微地震震源参数的导数。

(5)

式中： Δti,j代表第i个震源第j条射线的旅行时扰动；vk和dk分别是第k层内的速度分布和第k层界面的深度； Δvk和Δdk则分别是相应速度和层界面深度的扰动量；N为射线穿过的层数；K是与射线相交的界面数； Δxi,k(k=1，2，3，4) 是第i个微地震事件分别在x、y、z坐标及发震时刻的扰动量。

式(5)中第一项，即旅行时关于速度的一阶偏导数为

(6)

式中L为射线在第k层内的长度。

式(5)中的第二项，即旅行时关于层界面深度变化的偏导数参照Rawlison[26]的方法根据图1分析并推导。当层界面的深度变化Δd时(从图1中实线到虚线位置)，引起的旅行时变化为

(7)

则旅行时关于层界面深度变化的偏导数为

(8)

式中：向量wj、wj+1和wz如图1所示；vj和vj+1则分别是层界面两侧的速度；e和f分别为界面扰动后界面上、下两侧射线长度的变化量。

式(5)中的第三项，即旅行时关于微地震事件参数(位置和发震时刻)变化的偏导数，可由以下步骤计算。首先，旅行时(t=ta-t0)关于发震时刻的导数为

(9)

其次，旅行时关于微地震事件位置变化的偏导数的计算，可参见图2所示。若微地震震源的深度扰动值为Δzh，则旅行时关于震源位置xh、yh、zh的偏导数[26]为

(10)

式中：θh如图2所示；vh为震源处的速度值；φh、ψh分别为射线水平面投影后与x和y轴的夹角。

图2　旅行时关于震源位置变化偏导数计算示意图

3　数值模拟实例

为验证该算法的有效性及实用性，本文选用一个三维模型同时进行地下介质的层析成像及微地震事件的定位试验。如图3所示，模型尺寸为200m×200m×50m，网格尺寸为4m×4m×5m。模型是四层的水平层状介质，速度从地表向下依次为2.0、2.4、2.8和3.2km/s，界面的深度分别为16、30和40m。为模拟地面和井中微地震，采用实际微地震中常用的地表“米”字型和井中联合观测方式(如图3中三角形所示)，地表采集空间间隔为20m，井中采集空间间隔为10m，共46个检波器。设定了4个虚拟微地震事件，其空间坐标分别为S1(20m，40m，42m)、S2(100m，100m，45m)、S3(150m，180m，48m)、S4(170m，30m，26m)，激发时刻分别为10、15、20和5ms。

图3　三维模型及微地震观测系统三角形代表观测点，圆点代表微地震事件

在理想条件下(正确的速度模型、无噪声)采用该反演算法进行微地震事件S1、S2和S3定位。定位结果如图4所示，即使初始震源参数与真实值相差较大，反演的震源位置(图4a～图4c)和发震时刻(图4d)也都与理论参数完全一致。

然而实际数据处理中不可能满足理想条件，即地下介质的速度结构(层速度和层界面位置)未知、数据含噪声等。下面给出本文反演算法同时进行速度模型反演和微地震定位的结果，并讨论其对噪声的敏感性。首先选用两种不同的微地震事件组合测试该反演算法同时进行速度模型反演和微地震定位。第一种采用三个微地震事件(S1、S2和S3)的初至P波旅行时，第二种采用四个微地震事件(S1、S2、S3和S4)的初至P波旅行时，同时反演的结果如图5(速度结构反演结果)和图6(微地震定位结果)所示。

图4　理想条件下的微地震定位结果

从图5 可以看出，当仅采用三个微地震事件的旅行时数据时，第三层的速度结构(层界面和层速度)层析的结果(图5中蓝线)与真实模型(图5中红线)差别较大；而当采用四个微地震事件的旅行时数据时，层析结果明显改善，无论是层速度值还是层界面位置(图5中绿线)都接近真实值。出现这种情况是因为增加的微地震事件S4位于第二层介质中，与微地震事件S1、S2和S3位置有较大区别，射线覆盖及射线交错概率有所增加，提高了层析成像的分辨率。所以当速度结构反演结果得到改善时，相应的微地震定位结果精度也有所提高，特别是震源深度(图6b和6c)的定位结果。因此，当射线的角度覆盖较好时，同时进行速度层析成像和微地震定位是可行的，反演精度较高。

图5　不同微地震事件组合同时反演的速度结构

图6　不同微地震事件组合同时反演的微地震定位结果

考虑拾取初至旅行时可能存在误差，在理论的初至P波旅行时数据中加入幅度为0.4ms的随机误差，四个事件同时反演的结果如图7和图8所示。结果显示：加入随机噪声后，对微地震定位精度影响不大，且微地震震源水平位置的定位结果(图8a)明显好于震源深度的定位(图8b和图8c)。表1给出了具体的定位误差，可见未加随机噪声之前水平定位误差在几厘米的范围内，加入随机噪声后水平定位误差在十厘米左右，而深度的定位误差为1～2m，但仍在可接受的误差范围内。速度结构层析的结果变化较大(图7)，层速度和层界面深度都存在一定程度的偏差。这是由于当速度结构层析和微地震定位同时进行时，两者的反演结果实际上是一定程度的折中，并且层速度和层界面的深度在反演时也存在平衡的问题。总体而言，微地震定位对旅行时中的随机噪声不敏感。

图7　噪声敏感性试验速度层析结果

x误差/my误差/mz误差/mt0误差/ms无噪加噪无噪加噪无噪加噪无噪加噪S10．002-0．013-0．026-0．1301．3991．975-0．077-0．109S2-0．0010．1680．0010．0500．8240．341-0．053-0．002S30．0140．118-0．0640．0240．8541．014-0．055-0．040S4-0．092-0．0620．1430．145-0．1690．233-0．103-0．227

图8　噪声敏感性试验微地震定位结果

上述的地面和井中联合反演结果均迭代了20次左右。

4　结束语

为降低速度模型估计不准对微地震定位精度的影响，本文提出了对水平层状介质同时进行速度层析成像和微地震定位的反演算法。首先采用最短路径算法计算初至波的射线路径和旅行时，导出旅行时关于层速度、层界面深度及震源参数的偏导数，然后结合共轭梯度算法求解带约束的阻尼最小二乘问题，发展了一种同时反演速度结构(层速度和层界面)和震源参数(震源位置和发震时刻)的方法。地面和井中联合观测条件下的数值模拟结果表明：该反演算法能有效反演地下介质的速度结构并进行微地震定位，而且微地震定位结果对随机噪声不敏感。

本文基于模型正演无噪和加入随机噪声的数据，验证了所提出方法的有效性。在实际生产中，还存在地层起伏不平、介质非均质性强、数据信噪比低等诸多因素，因此还有必要结合实际数据与现场地质情况，应用、测试方法的可靠性和稳健性，提高实际定位精度。

[1]张光明.水平井水力压裂数值模拟研究[学位论文].安徽合肥：中国科学技术大学，2010.

Zhang Guangming.A Numerical Simulation Study on Hydraulic Fracturing of Horizontal Well[D].University of Science and Technology of China,Hefei,Anhui,2010.

[2]Rutledge J T and Phillips W S.Hydraulic stimulation of natural fractures as revealed by induced micro-earthquakes,Carthage Cotton Valley gas field,east Texas.Geophysics,2003,68(2)：441-452.

[3]Maxwell S C,Rutledge J,Jones R et al.Petroleum reservoir characterization using downhole microseismic monitoring.Geophysics,2010,75(5)：A129-A137.

[4]容娇君,李彦鹏,徐刚等.微地震裂缝检测技术应用实例.石油地球物理勘探,2015,50(5)：919-924.

Rong Jiaojun,Li Yanpeng,Xu Gang et al.Fracture detection with microseismic.OGP,2015,50(5)：919-924.

[5]Li J L,Kuleli S,Zhang H J et al.Focal mechanism determination of induced microearthquakes in an oil field using full waveforms from shallow and deep seismic networks.Geophysics,2011,76(6)：WC87-WC101.

[6]张山,刘清林,赵群等.微地震监测技术在油田开发中的应用.石油物探,2002,41(2)：226-231.

Zhang Shan,Liu Qinglin,Zhao Qun et al.Application of microseismic monitoring technology in development of oil field.GPP,2002,41(2)：226-231.

[7]Warpinski N R,Sullivan R B,Uhi J et al.Improved microseismic fracture mapping using perforation timing measurements for velocity calibration.Proceedings of the SPE Annual Technical Conference and Exhibition,Denver,Colorado,2003,84488.

[8]Van Dok R,Fuller B,Engelbrechet L et al.Seismic anisotropy in microseismic event location analysis.The Leading Edge,2011,30(7)：766-770.

[9]Block L V,Cheng C H,Fehler M C et al.Seismic imaging using microearthquakes induced by hydraulic fracturing.Geophysics,1994,59(1)：102-112.

[10]张晓林,张峰,李向阳等.水力压裂对速度场及微地震定位的影响.地球物理学报,2013,56(10)：3552-3560.

Zhang Xiaolin,Zhang Feng,Li Xiangyang et al.The influence of hydraulic fracturing on velocity and microseismic location.Chinese Journal of Geophysics,2013,56(10)：3552-3560.

[11]Song F X,Kuleli H S,Toksoz M N et al.An improved method for hydrofracture-induced microseismic event detection and phase picking.Geophysics,2010,75(6)：A47-A52.

[12]宋维琪,李艳清,刘磊.独立分量分析与压缩感知微地震弱信号提取方法.石油地球物理勘探,2017,52(5)： 984-989,1041.

Song Weiqi,Li Yanqing,Liu Lei.Microseismic weak signal extraction based on the independent component analysis and compressive sensing.OGP,2017,52(5)：984-989,1041.

[13]王晨龙,程玖兵,尹陈等.地面与井中观测条件下的微地震干涉逆时定位算法.地球物理学报,2013,56(9)：3184-3196.

Wang Chenlong,Cheng Jiubing,Yin Chen et al.Microseismic events location of surface and borehole observation with reverse-time focusing using interfero-metry technique.Chinese Journal of Geophysics,2013,56(9)：3184-3196.

[14]李振春,盛冠群,王维波等.井地联合观测多分量微地震逆时干涉定位.石油地球物理勘探,2014,49(4)：661-666,671.

Li Zhenchun,Sheng Guanqun,Wang Weibo et al.Time-reverse microseismic hypocenter location with interferometric imaging condition based on surface and downhole multi-components.OGP,2014,49(4)：661-666,671.

[15]王鹏,常旭,王一博.基于时频稀疏性分析法的低信噪比微地震事件识别与恢复.地球物理学报,2014,57(8)：2656-2665.

Wang Peng,Chang Xu,Wang Yibo.Automatic event detection and event recovery in low SNR microseismic signals based on time-frequency sparseness.Chinese Journal of Geophysics,2014,57(8)：2656-2665.

[16]李稳,刘伊克,刘保金.基于稀疏分布特征下的井下微地震信号识别与提取方法.地球物理学报,2016,59(10)：3869-3882.

Li Wen,Liu Yike,Liu Baojin.Downhole microseismic signal recognition and extraction based on sparse distribution features.Chinese Journal of Geophysics,2016,59(10)：3869-3882.

[17]曹俊海,顾汉明,尚新民.基于局部相关谱约束的多道匹配追踪算法识别微地震信号.石油地球物理勘探,2017,52(4)：704-714.

Cao Junhai,Gu Hanming,Shang Xinmin.Microseismic signal identification with multichannel matching pursuit based on local coherence spectrum constraint.OGP,2017,52(4)：704-714.

[18]Grechka V,Singh P and Das I.Estimation of effective anisotropy simultaneously with locations of microseismic events.Geophysics,2011,76(6)：WC141-WC155.

[19]Li J L,Zhang H J,Rodi W L et al.Joint microseismic location and anisotropic tomography using differential arrival times and differential back azimuths.Geophy-sical Journal International,2013,195(3)：1917-1931.

[20]谭玉阳,何川,张洪亮.基于初至旅行时差的微地震速度模型反演.石油地球物理勘探,2015,50(1)：54-60.

Tan Yuyang,He Chuan,Zhang Hongliang.Microseismic velocity model inversion based on moveouts of first arrivals.OGP,2015,50(1)：54-60.

[21]Rutledge J T,Phillips W S and Schuessler B K.Reservoir characterization using oil-production-induced microseismicity,Clinton County,Kentucky.Tectonophysics,1998,289(1)：129-152.

[22]Bai C Y,Tang X P and Zhao R.2D/3D multiply transmitted,converted and reflected arrivals in complex layered media with the modified shortest path me-thod.Geophysical Journal International,2009,179(1)：201-214.

[23]Menke W.Geophysical Data Analysis：Discrete In-verse Theory.Academic Press Inc.,San Diego,California,1984,40-126.

[24]Tarantola A,Valette B.Generalized non-linear inverse problems solved using the least squares criterion.Reviews of Geophysics Space Physics,1982,20(2)：219-232.

[25]Zhou B,Sinadinovski C,Greenhalgh S A.Iterative algorithm for the damped minimum norm,least-squares and constrained problem in seismic tomography.Exploration Geophysics,1992,23(3)：497-505.

[26]Rawlison N and Sambridge M S.Seismic traveltime tomography of the crust and lithosphere.Advances in Geophysics,2003,46：81-197.