2. 中科院地质与地球物理研究所, 北京 100029
2. Institute of Geology and Geophysics, Chinese Academy of Sciences, Beijing 100029, China
地震波在地层中传播时由于受到地层的吸收衰减作用而产生能量耗散, 因此用黏弹性介质模型描述地震波的传播过程更加符合实际情形,黏弹性介质中传播的地震波携带了介质的相关物性参数信息.地球从大尺度的地质构造到小尺度的孔隙结构都存在物性参数的不均匀性, 实际地层的岩石通常是由岩石骨架和各种黏性的孔隙填充物组成的多相黏弹性介质(Kuster and Toksöz, 1974),由于岩石的内摩擦作用以及填充物的黏滞作用,地震波在介质中传播时将会发生速度频散和衰减现象(王海洋等,2012).相关研究表明地下介质对地震波的吸收衰减作用与介质的孔隙度、孔隙填充物有密切联系(Winkler and Nur, 1982; Pride et al., 2004; ller et al., 2010).地层对地震波的吸收衰减效应用品质因子Q来表征, 它是介质非完全弹性的度量, 是地球介质的本质属性之一, Varela研究了这种衰减规律(Varela et al., 1993),Jacobsen研究了沉积岩中地震波衰减和速度弥散以及频率之间的相互关系(Jacobsen, 1987).
地球物理反演算法大致可以分为线性反演和非线性反演两类:在线性反演中,各种不同的线性近似被引入到反演算法中,例如Born近似(Oristaglio, 1985)、Rytov近似(Devaney, 1981)等,这些算法只在某些特定情况下有效.非线性反演算法没有改变反演问题的非线性本质,例如全波形反演方法(Symes, 2008;张文生等,2015)、遗传算法(Stoffa and Sen, 1991)、神经网络方法(Zhang and Paulson, 1997)等,这些算法通过不断迭代的方式对地下构造进行反演成像.地球物理反演可以在时间域进行(Moghaddam and Chew, 1992),也可以在频率域进行(Pratt, 1999),Tarantola对考虑衰减效应的波形反演理论进行了研究(Tarantola, 1998),地震波在地层介质中的传播速度是地球探测活动中非常重要的参数之一, 在黏弹性介质速度反演研究方面,龙桂华等利用失配函数二范数最小准则,用预条件梯度类方法对黏弹性声波介质的速度结构进行了逐频反演(龙桂华等,2009),Causse等在黏弹性介质全波形反演过程中引入了两种预处理方法(Causse et al., 1999),取得了较好的效果.
随着地震波勘探活动范围的不断扩大,勘探目标越趋复杂,在地表接收不到某些构造的反射波,但可以接收到散射波,基于散射波场的反演方法是地球物理领域中一种非常有效的反演方法之一(黄联捷和杨文采,1991;杨晓春等,2005).Van den Berg等进行了基于电磁散射波场的反演研究(Van den Berg and Kleinman, 1997; Van den Berg et al., 1999; Van den Berg and Abubakar, 2001),该方法是源积分类方法的拓展(Habashy et al., 1994),该反演算法以背景模型是均匀介质为基本前提,Abubakar将该方法进一步扩展为适应非均匀背景介质的情形(Abubakar et al., 2008).在最初的电磁波反演问题中,反演尺度范围较小,没有考虑电磁波在介质中的吸收衰减效应,反演过程基于完全弹性声波方程.在地震波勘探活动中,地球介质尺度非常大,吸收衰减效应将不能忽略,本文反演过程基于黏弹性声波方程, 更加符合地震波在地球介质中的真实传播过程.
2 正演模拟 2.1 控制方程逆散射反演问题示意图如图 1所示,我们定义反演域D和数据域S,其中炮点和检波器位于数据域S中,两者的组合称为总计算区域T.地层对地震波的吸收衰减效应用品质因子Q来表征,地层不同位置对应不同的Q值,c表示地震波传播速度.相对于真实速度模型,假设有一个背景模型,引入反演变量χ,该变量表示已知背景模型和未知真实模型物性参数之间的相对扰动量.反演域中的每一个网格点都可以看作是相对于背景模型的一个散射点,通过求解一个最优化问题来不断更新χ,从而得到模型的速度分布.
本文假设反演过程中介质密度ρ为常数,考虑地层吸收衰减效应时地震波总波场pj(r)在频率-空间域满足Helmholtz方程为
(1) |
其中∇2表示Laplace算子,pj(r)表示总波场,j为地表不同位置炮点的标号, Source(r, rjs)表示震源项,rjs表示震源位置,T表示地下半空间,k(r)为复波数,公式为
(2) |
式中ω为地震波频率,
(3) |
其中c为地震波在介质中传播的实速度,当不考虑吸收衰减效应时用c代替式(2)中的
(4) |
利用复速度与复波数间的关系,我们可以正演得到受地层吸收衰减影响的波场数据,本文参考频率-空间域9点差分格式模拟地震波在吸收衰减介质中的传播过程(Jo et al., 1996).
2.2 吸收边界条件在对地震波传播过程进行有限差分数值模拟过程中,为了有效吸收入射到计算区域边界的波场,采用吸收边界条件是必要的.PML吸收边界条件是非常有效的吸收边界反射的方法(Berenger, 1994),本文使用复数延伸坐标PML吸收边界条件来消除计算边界区域产生的人工反射(Chew and Weedon, 1994).在吸收边界区域式可以写为
(5) |
其中
将考虑地层吸收衰减的地震波传播总波场写成背景波场和散射波场之和为
(6) |
背景波场pjbac满足方程:
(7) |
其中kb(r)为背景模型对应的复波数,此处应当注意,当反演过程中考虑地层吸收衰减效应时,kb(r)的表达式中地震波速度用复速度
将式从式中减去,得到考虑吸收衰减效应时散射波场满足的方程为
(8) |
此处引入一个变量wj(r),公式为
(9) |
其中χ(r)定义为
(10) |
式中
传统反演算法的优化目标函数为检波器接收数据与正演模拟数据的残差值,本文方法的优化目标函数由数据域误差和反演域误差以及基于反演参数总变差值的正则项组合而成,反演域误差项和正则项的引入在一定程度上提高了算法的稳定性和抗干扰性.算法利用非线性共轭梯度法(Fletcher and Reeves, 1964; Gilbert and Noedal, 1992)最小化目标优化函数从而不断更新w,总波场p以及χ的值,最终得到地下介质的速度参数分布.
本文算法的反演目标优化函数为
(11) |
其中fj表示第j炮地表检波器接收到的散射波场数据,Lb表示频率-空间域声波方程正演算子.ρj,n表示数据域误差,rj,n表示反演域误差,
(12) |
利用非线性共轭梯度法,w迭代更新公式为
(13) |
其中αj, nw表示迭代步长,vj, n表示Polka-Ribiere搜索方向,公式为
(14) |
由Frechet导数的定义,目标优化函数关于wj的梯度gj, nw为
(15) |
其中归一化系数
更新完w以后,利用式(6)、式(8)更新总波场p值,然后利用非线性共轭梯度法更新χ值,公式为
(16) |
正则化目标优化函数Fn(χ, wj)关于χ的梯度为
(17) |
其中
(18) |
重复上面的迭代过程直到满足预先设定的误差限,最后由χn得到模型的速度分布,本文利用反传播方法初始化wj(0), 总波场pj(0)和χ(0)的取值(Habashy et al., 1994).
4 数值算例本文利用二维Marmousi模型对反演算法进行测试,该速度模型非常复杂,真实速度模型见图 2a,模型速度范围为1500~5500 m·s-1, 在x和z方向被离散为384×122个等间距网格,网格间距为15 m,震源采用主频为7.5 Hz的雷克子波.我们采用四个频率(4 Hz, 12 Hz, 20 Hz, 24 Hz)依次进行反演,地震数据由均匀分布于地表的64炮产生,炮间距为90 m,检波点布置在地表所有网格点上.假设初始背景模型速度分布为由上到下线性增加,速度范围为1500~4200 m·s-1,初始背景模型见图 2b.
在实际地震波勘探中,地表检波器接收到的波场数据是受地层吸收衰减作用后的波场数据,首先我们在反演过程中使用的正演算子考虑了地层的吸收衰减效应,反演中将低频反演结果作为下一个高频反演的背景模型.单频反演结果见图 3所示,从图中可以看出反演频率较低时,反演结果给出了模型的大尺度构造,随着反演频率的逐渐增加,重建速度模型的分辨率也逐步提高,Marmousi模型中的一些小尺度细节构造也能够较好的重建.
图 4所示为考虑地层吸收衰减效应的各频率反演速度模型与真实速度模型、初始速度模型在地表 1485 m和3435 m位置处沿深度方向的速度比较图,从速度比较曲线看出,初始速度模型为线性递增的,4 Hz反演结果在初始模型的基础上给出了真实模型的大尺度轮廓,12 Hz、20 Hz、24 Hz反演结果不断逼近真实速度模型.在中浅层反演模型与真实模型基本完全一致,由于Marmousi模型深度1400 m以下的速度纵横向变化剧烈,包含大量薄地层结构,异常高速体下面又隐藏有低速体构造,从反演结果可以看出,模型深部反演结果并不理想,这是由于地震波在地层中传播时由于地层吸收衰减效应使得地震波能量不断减弱,地震波频谱中高频成分的衰减较低频成分更为严重,导致反演结果的分辨率降低.
作为对比,本文对在反演过程中不考虑地层吸收衰减效应的反演结果进行了分析,即在反演过程中正演算子不包含品质因子Q的影响,反演结果如图 5所示,可以看出,由于在反演过程中所采用的正演算子不符合地震波在地层中的真实传播过程,与前面反演结果相比,在4 Hz时反演得到了模型的大尺度速度结构,但随着反演频率的提高,反演所得速度模型的分辨率没有得到相应提高,12 Hz、20 Hz、24 Hz反演结果几乎一样,反演结果与真实模型有较大差距.
图 6所示为相同地表水平位置处,不考虑地层吸收衰减效应时沿深度方向的速度比较曲线,可以看出反演速度模型与真实模型有较大差距,随着反演频率的提高,模型细节尺度的构造无法得到重建,甚至在模型浅层也没能给出正确的重建速度模型结果.
图 7为地表 3435 m处沿深度方向真实模型与两者24 Hz反演结果的速度分布比较曲线,由于地表检波器接收到的波场数据为受地层吸收衰减作用以后的数据,可以看出考虑地层衰减效应的反演结果更加接近于真实速度模型,因为反演过程中的正演算子更符合地震波在地层中的真实传播过程,相反,反演中不考虑衰减效应则无法得到正确的速度重建结果.
为了进一步提高反演计算效率,本文采用基于MPI的并行算法对多炮地震波数据进行并行计算,反演过程中使用了10个独立进程,单频反演所耗计算时间约为3 h.图 8所示为反演迭代过程中目标函数残差值随迭代次数的变化曲线,可以看出计算过程收敛很快,计算效率较高.
背景模型的选取对反演精度会产生一定的影响,在前面线性背景模型的基础上,本文又采用将原模型进行平滑处理后作为初始背景模型进行了逐频反演,初始背景模型如图 9所示,此外为检验方法的抗噪能力,反演过程中在每个单频反演中将检波器散射波数据加入5%的随机噪声,公式为
(19) |
其中r表示地表检波器的位置,pj′(r)表示第j炮产生的含有随机噪声的位于地表r处的检波器散射波数据,pj(r)表示不含噪声的散射波数据,a和b为(-1, 1)之间的随机数,单频反演结果如图 10所示.
从图 10可以看出,选取平滑化初始背景模型得到的反演重建结果可以得到模型的细节尺度构造,与前面线性背景模型反演结果相比,重建模型深部的反演精度和分辨率有了一定的提升,图 11为模型地表 1485 m、2235 m、3435 m、4485 m处速度分布随地层深度的变化曲线,从曲线比较中也可以得到相同的结论.此外从反演结果可以看出,在地震波数据包含随机噪声的情况下,由于正则化处理使得反演方法具有一定的抗噪能力.
我们利用式(20)对不同频率反演重建模型相对于真实模型的反演误差进行统计分析,表 1为两种不同初始背景模型反演结果相对于真实模型的误差值,公式(20)为
(20) |
从表 1可以看出,背景模型的选取对最终反演结果的精度将产生较大的影响,选择平滑化背景模型反演精度比线性背景模型高,表明背景模型的选取越接近于真实模型,反演结果精度越高.
5 结论地层对地震波的吸收衰减效应是地球介质的固有属性之一,地表检波器接收到的地震波数据是经过地层吸收衰减作用之后的波场数据,地层的吸收衰减效应对于地球内部物性构造的准确反演会产生较大的影响,如果在反演过程中不考虑这一效应,就无法得到正确的地层速度模型,导致对地层结构做出错误的判断,本文基于黏弹性声波方程的逆散射反演方法更接近于真实情形.本文在考虑地层衰减效应的基础上,对地震波速度进行反演研究,反演结果在地层中浅层分辨率较高,在地层深部由于地震波能量的衰减导致反演重建模型分辨率不太理想.同时背景模型的选取对反演结果的精度会产生一定的影响,背景模型越接近真实模型,反演结果精度越高,在实际地震勘探过程中利用其他勘探手段对地质模型进行初探,选择较为接近真实地层结构的背景模型,可以得到更高精度的反演结果.
地球的地层结构非常复杂,同时影响地震波传播过程的吸收衰减效应也很复杂,对于衰减吸收的有关特性还有待进一步的深入研究,本文假设品质因子Q与地震波速度之间满足一定的经验关系,后续研究可以进一步推广到地震波速度和地层品质因子Q同时反演的情形以及探寻如何提高地层深部反演分辨率的有效措施.
Abubakar A, Hu W, Van den Berg P M, et al. 2008. A finite-difference contrast source inversion method. Inverse Problems , 24 (6) : 065004. DOI:10.1088/0266-5611/24/6/065004 | |
Berenger J P. 1994. A perfectly matched layer for the absorption of electromagnetic waves. Journal of Computational Physics , 114 (2) : 185-200. DOI:10.1006/jcph.1994.1159 | |
Causse E, Mittet R, Ursin B. 1999. Preconditioning of full-waveform inversion in viscoacoustic media. Geophysics , 64 (1) : 130-145. DOI:10.1190/1.1444510 | |
Chew W C, Weedon W H. 1994. A 3D perfectly matched medium from modified Maxwell's equations with stretched coordinates. Microwave and Optical Technology Letters , 7 (13) : 599-604. DOI:10.1002/(ISSN)1098-2760 | |
Devaney A J. 1981. Inverse-scattering theory within the Rytov approximation. Optics Letters , 6 (8) : 374-376. DOI:10.1364/OL.6.000374 | |
Fletcher R, Reeves C M. 1964. Function minimization by conjugate gradients. Computer Journal , 7 (2) : 149-154. DOI:10.1093/comjnl/7.2.149 | |
Gilbert J C, Noedal J. 1992. Global convergence properties of conjugate gradient methods for optimization. SIAM Journal on Optimization , 2 (1) : 21-42. DOI:10.1137/0802003 | |
Habashy T M, Oristaglio M L, de Hoop A T. 1994. Simultaneous nonlinear reconstruction of two-dimensional permittivity and conductivity. Radio Science , 29 (4) : 1101-1118. DOI:10.1029/93RS03448 | |
Huang L J, Yang W C. 1991. Inversion of the acoustic wave equation by inverse scattering in a medium of linear reference velocity. Acta Geophysica Sinica , 34 (1) : 89-98. | |
Jacobsen R S. 1987. An investigation into the fundamental relationships between attenuation, phase dispersion, and frequency using seismic refraction profiles over sedimentary structures. Geophysics , 52 (1) : 72-87. DOI:10.1190/1.1442242 | |
Jo C H, Shin C, Suh J H. 1996. An optimal 9-point, finite-difference, frequency-space, 2-D scalar wave extrapolator. Geophysics , 61 (2) : 529-537. DOI:10.1190/1.1443979 | |
Kuster G T, Toks z M N. 1974. Velocity and attenuation of seismic waves in two-phase media:part 1.Theoretical formulations. .Geophysics , 39 (5) : 587-606. DOI:10.1190/1.1440450 | |
Long G H, Li X F, Zhang M G, et al. 2009. Visco-acoustic transmission waveform inversion for velocity structure in space-frequency domain. Acta Seismologica Sinica , 31 (1) : 32-41. | |
Moghaddam M, Chew W C. 1992. Nonlinear two-dimensional velocity profile inversion using time domain data. IEEE Transactions on Geoscience and Remote Sensing , 30 (1) : 147-156. DOI:10.1109/36.124225 | |
Müller T M, Gurevich B, Lebedev M. 2010. Seismic wave attenuation and dispersion resulting from wave-induced flow in porous rocks:a review. Geophysics , 75 (5) : 75A. | |
Oristaglio M L. 1985. Accuracy of the Born and Rytov approximations for reflection and refraction at a plane interface. J. Opt. Soc. Am , 2 (11) : 1987-1993. DOI:10.1364/JOSAA.2.001987 | |
Pratt R G. 1999. Seismic waveform inversion in the frequency domain:Part I-Theory, and verification in a physical scale model. Geophysics , 64 (3) : 888-901. DOI:10.1190/1.1444597 | |
Pride S R, Berryman J G, Harris J M. 2004. Seismic attenuation due to wave-induced flow. Journal of Geophysical Research , 109 : B01201. | |
Ravaut C, Operto S, Improta L, et al. 2004. Multiscale imaging of complex structures from multifold wide-aperture seismic data by frequency-domain full-waveform tomography:Application to a thrust belt. Geophysical Journal International , 159 (3) : 1032-1056. DOI:10.1111/gji.2004.159.issue-3 | |
Stoffa P L, Sen M K. 1991. Nonlinear multiparameter optimization using genetic algorithms:inversion of plane-wave seismograms. Geophysics , 56 (11) : 1794-1810. DOI:10.1190/1.1442992 | |
Symes W W. 2008. Migration velocity analysis and waveform inversion. Geophysical Prospecting , 56 (6) : 765-790. DOI:10.1111/gpr.2008.56.issue-6 | |
Tarantola A. 1988. Theoretical background for the inversion of seismic waveforms including elasticity and attenuation. Pure and Applied Geophysics , 128 (1-2) : 365-399. DOI:10.1007/BF01772605 | |
Tian S R. 1990. Estimation of Q value in inverse Q filtering with Li formula. Oil Geophysical Prospecting , 25 (3) : 354-361. | |
Van den Berg P M, Abubakar A. 2001. Contrast source inversion method:state of art. Progress in Electromagnetics Research , 34 : 189-218. DOI:10.2528/PIER01061103 | |
Van den Berg P M, Kleinman R E. 1997. A contrast source inversion method. Inverse Problems , 13 (6) : 1607-1620. DOI:10.1088/0266-5611/13/6/013 | |
Van den Berg P M, Van Broekhoven A L, Abubakar A. 1999. Extended contrast source inversion. Inverse Problems , 15 (5) : 1325-1344. DOI:10.1088/0266-5611/15/5/315 | |
Varela C L, Rosa A L R, Ulrych T J. 1993. Modeling of attenuation and dispersion. Geophysics , 58 (8) : 1167-1173. DOI:10.1190/1.1443500 | |
Wang H Y, Sun Z D, Chapman M. 2012. Velocity dispersion and attenuation of seismic wave propagation in rocks. Acta Petrolei Sinica , 33 (2) : 332-342. | |
Winkler K W, Nur A. 1982. Seismic attenuation:Effects of pore fluids and frictional-sliding. Geophysics , 47 (1) : 1-15. DOI:10.1190/1.1441276 | |
Yang X C, Li X F, Zhang M G. 2005. The fixed point principles of nonlinear seismic scattering inversion. Progress in Geophysics , 20 (2) : 496-502. | |
Zhang W S, Luo J, Teng J W. 2015. Frequency multiscale full-waveform velocity inversion. Chinese J. Geophys , 58 (1) : 216-228. | |
Zhang Y C, Paulson K V. 1997. Magnetotelluric inversion using regularized Hopfield neural networks. Geophysical Prospecting , 45 (5) : 725-743. DOI:10.1046/j.1365-2478.1997.660299.x | |
黄联捷, 杨文采. 1991. 参考波速线性变化时的声波方程逆散射反演. 地球物理学报 , 34 (1) : 89–98. | |
龙桂华, 李小凡, 张美根, 等. 2009. 频率域粘弹性声波透射波形速度反演. 地震学报 , 31 (1) : 32–41. | |
田树人. 1990. 用李氏经验公式估算反Q滤波中的Q值. 石油地球物理勘探 , 25 (3) : 354–361. | |
王海洋, 孙赞东, ChapmanM. 2012. 岩石中波传播速度频散与衰减. 石油学报 , 33 (2) : 332–342. | |
杨晓春, 李小凡, 张美根. 2005. 地震波散射非线性反演的不动点理论研究. 地球物理学进展 , 20 (2) : 496–502. | |
张文生, 罗嘉, 滕吉文. 2015. 频率多尺度全波形速度反演. 地球物理学报 , 58 (1) : 216–228. | |