不依赖于子波的全波形反演方法研究【推荐论文】 .doc

精品论文不依赖于子波的全波形反演方法研究刘宝龙，董良国，迟本鑫（同济大学海洋地质国家重点实验室，上海 200092）5摘要：地震全波形反演（FWI）在理论上具有最高的反演分辨率，但目前在实际应用中遇到 诸多挑战，挑战之一是传统的 FWI 需要已知每炮的震源子波。但是，实际震源子波的准确反演是比较困难的，目前常规 FWI 通常是对震源子波做某种估计，这样 FWI 反演结果易受 所估计的不准确的震源子波的影响。本文通过将观测数据和正演模拟数据分别与各自的参考10波场相除或交叉相乘后再构建目标函数，使得反演方法不再依赖于地震子波（这两种方法分 别称为褶积法和反褶积法），避免了实际子波估计不准确而对全波形反演结果的影响。关键词：全波形反演；不依赖于子波；褶积法；反褶积法中图分类号：P315.815Method of Source-independent Full Waveform InversionLIU Baolong, DONG Liangguo, CHI Benxin(State Key Laboratory of Marine Geology,Tongji University, Shanghai 200092)Abstract: Full-waveform inversion of seismic data possesses the highest resolution. However, it has met several challenges in practice. One of them is the lack of precise knowledge of the source20wavelet. Since the field source wavelet is hard to invert and currently available approaches involvesome form of approximations to the source, the FWI inversion results are subject to the validity of the estimated source. In this paper, both the observed data and the forward-modeled data are multiplied by the respective reference wavefields on the cross or divided by their respective reference wavefields. Then the two terms of corresponding normalized wavefields are subtracted25to construct an objective function, which is called convolution and deconvolution method. These two modifications make the FWI method to be independent with the source wavelet. So the influence of the inexact estimated wavelet to FWI has been avoided.Key words: full waveform inversion; source independence; convolution;deconvolution300引言自从地震波形反演方法在 1983 年第一次被 Lines 和 Kelly 提出以来，因其能给出精细的 地下结构，FWI 方法开始被大量研究（Pratt 等，1998）1。特别是 Lailly 2和 Tarantola 3提 出可以用伴随状态法计算梯度以后，FWI 就被广泛的应用于时域和频域旅行时层析成像反 演和波形反演中4-11。35在这样的全波形反演中，在每步迭代中除了需要更新速度模型之外，都需要迭代更新子 波。然而，在野外实际数据中，通常没有足够的有效震源子波信息，因此很多情况下不能反 演出有效的实际震源子波。为此，Lee 和 Kim12与 Zhou 和 Greenhalgh13提出不依赖子波的 FWI 算法，其中他们用每炮第一道或所有道波场的平均值对波场做归一化处理，缺点在于 他们直接计算 Jacobian 矩阵。40在本文中，我们发展了不依赖子波的归一化的 FWI 算法（褶积法和反褶积法）。通过 对波场进行归一化来构建一个与子波无关的目标函数。基于在频域模拟中所构建的伴随状态基金项目：博士点基金项目（20090072110030）；自然科学基金项目（41274116） 作者简介：刘宝龙（1988），男，同济大学硕士研究生，主要研究方向：地震波全波形反演 通信联系人：董良国（1966-），男，同济大学教授，主要研究方向：地震波传播与模拟以及地震波反演. E-mail:dlg- 12 -的矩阵形式，通过反传归一化的野外数据和模型响应间的残差，并跟虚拟源（在正演模拟的处理中计算出来的）相乘，就能隐式地计算出梯度，这样在该梯度的计算中就削弱和消除了45常规因估计出的可能不正确的震源子波对 FWI 造成的影响，从而有效的减少了 FWI 对子波 的依赖性。1反褶积法实际地震勘探记录到的波场可以表示为如下形式:d (t) = sd (t) * I d (t) ， i = 1, 2,., ns ， j = 1, 2,., nr（1）ij i ij50其中 sd (t ) 是 t 时刻的震源子波，I d (t) 是波动方程的脉冲响应，上标 d 表示记录到的数i ij据，下标 i 和 j 分别表示第 i 炮和第 j 个检波器，* 表示褶积。 设w 是角频率，那么（1）式在频域的形式为：d (w) = sd (w)I d (w)（2）ij i ij类似的，正演模拟的波场可以表达为:55u(w) = su (w)I u (w)（3）ij i ij其中上标 u 表示理论计算得到的正演模拟数据。设下标 R 表示同一炮的某一参考道，利用该参考道构建不依赖于子波的归一化的“观 测波场”：d (w) d(w) = S d (w)I d (w)S d (w)I d (w) = I d (w)I d (w)（4）ij iR i ij i iR ij iR60同样的，模拟数据的归一化的波场为：u (w) u(w) = Su (w)I u (w)Su (w)I u (w) = I u (w)I u (w)（5）ij iR i ij i iR ij iR用归一化的观测数据和模拟数据构造不依赖于地震子波的残差：DDij = dij (w)diR (w) - uij (w)uiR (w)（6）跟常规波形反演一样，利用（6）式构建 L2 范数目标函数：wmaxns nrE = 1 (DD) × (DD * )652 å ååij ij（7）-wmaxi =1j =1其中* 表示复共轭，wmax 表示最大角频率， ns 表示炮点数， nr 表示检波点数。 用伴随状态法计算目标函数（7）对速度参数 m 的梯度，过程如下：¶Ewmaxns nr¶DD *ij= Re å å å() × (DDij )¶m-wmaxi =1j =1 ¶m= Rewmaxns nruij× ¶uiR × ( diju- ij )* -1 ¶u d× ij × ( iju- ij )*í å å å u 2¶m d u u¶m d uý （8）-wmaxi =1j =1iR iR iR iR iR iR70为获得等式（8）中的偏导数项，引入频域有限差分法波动方程的矩阵形式：Au = s其中 A 是复阻抗矩阵， u 是波场矢量， s 是源矢量。 对等式（9）求其对速度参数 m 的导数，得（9）¶u = A-1 (- ¶A u) = A-1 f¶m ¶m（10）75其中 f 是虚拟源矢量（Pratt 等，1998）1。用这种关系，把等式（8）重新写为如下形式：ì é 1 ( di1 - ui1 )* ù üï ê ú ïïé 0 ùï ê M úê uiRêdiRMuiR ú ïú ïêï ê úê ú ï¶E = Re( A-1 f )u d u( -ú)* -( A-1 f )1 ( diR - uiR )*wmaxï nsTnrnsij ij ijT ê ú ïå íåêå 2ú å ê ú ý¶m -wmax ïi =1ê j =1uiR diR uiR úi =1uiR diR uiRê úïï ê M ú M ïï ê úê ú ïï 0 ê 1( dinr - uinr )* ú ïï ê ú ïî ë uiRdiRuiR û þ其中上标T 表示矩阵的转置。利用 ( A-1 )T = A-1 的对称性，得反褶积法 FWI 梯度的简洁形 式：wmax ns ¶A¶E / ¶m = Re å åuT ( )T A-1 (DI - DII )*（11）80其中*-wmax i =1 ¶mé 0 ùé (1 / uiR )(di1 / diR ) - (ui1 / uiR ) ù ê M úê ú ê úê M ú*ê 0 úê nrúDI = ê (1 / uiR )(diR / diR ) - (uiR / uiR ) ú ，êå ú2 *ê ú DII =ê (uij / uiR ) (dij / diR ) - (uij / uiR )Rth Row úê M ú*j =1ê ú0iR inr iR inriRê(1 / u )(d / d ) - (u / u )ú ê úë û ê M úê ú0从（11）式可以发现，将波场残差 DI - DII 反传，计算反传波场和正演波场的卷积就 可以求得不依赖于震源子波的 FWI 的梯度。2褶积法85在第二种方法中，用频域相乘构建残差。在第 i 个炮道集中，将模拟的第 j 道数据和观测的第 R 个参考道数据在频域相乘，将观测的第 j 道的数据和模拟的第 R 个参考道数据在 频域相乘，然后两者再相减得到新的波场残差：DDij = dij uiR - diRuij将（2）（3）等式带入到等式（12）的右端项得DD = s I s I - s I s I90d d u u d d u u iji ij i iRi iR i ij（12）（13）i i从等式（13）看出，在等式右端的两项中，公共项 s d su 相当于作为第 i 炮新的震源子波， 这就形成了对估计出的震源子波的依赖性有所降低的波场残差。利用第二类归一化数据残差（13），构建 L2 范数目标函数：E = 1wmaxns nrDD ) × (DD * )2 å åå(ij ij（14）-wmaxi =1j =195可以用伴随状态法计算目标函数（14）对速度参数的梯度：¶Ewmaxns nr *¶ = Re å å ådij (¶uiR¶m) - diR (¶uij¶m)(dij uiR - diRuij )mwmax- wmax i =1ns nrj =1*wmaxiR ij iR iR ijns nr *ij ij iR iR ij= Re å å å (¶uiR¶m) êéëd (d u - d u )úùû - Re å å å (¶uij¶m) êéëd (d u - d u ) úùû- wmax i =1j =1-wmax i =1j =1利用等式（10）把上式写为如下形式：ì é 0 ù é* ù üï ê Mú ê diR ( di1uiR - diRui1 ) ú ï¶Ewmaxï nsT ê nr* ú ns T ï= Re å íå ( A-1 f )êå d(d u- d u) ú - å ( A-1 f ) ê d( d u- d u)* ú ý¶m-wmaxï i =1ê j =1ij ij iR iR ijúi =1ê iR iR iR iR iR úï*þê M úï ê M úê ú ïï ê úîïêë0 úûdiR ( dinr uiR - diRuinr )ï100类似于反褶积法 FWI 梯度的推导，可以得到褶积法 FWI 的梯度：wmax ns ¶A¶E / ¶m = Reå åuT ( )T A-1 (DI - DII )*（15）-wmax其中i =1 ¶mé 0 ù*é diR( di1uiR - diRui1 ) ùê M úêúê úêMú*ê 0 úê nrúDI = ê diR( diRuiR - diRuiR ) ú ，ê * (- )úêúDII = êå dijdij uiRdiRuijRth Row úêMú*j =1ê úê 0 údê iR ( dinr uiR - diRuinr )úëûê M úê ú0ûú105类似于逆时偏和常规波形反演，将波场残差 DI - DII 反传，计算反传波场和正演波场 的卷积就可以求得不依赖于震源子波的 FWI 的梯度。3反演方案在波形反演中，通过用如下的迭代更新速度模型：mn+1 = mn + d ma md m = -H -1Ñ E（16）（17）110其中 mn+1 和 mn 是在第 n +1次和第 n 次迭代中的物性参数，d m 是第 n 次迭代时的参数 更新量。H a 是 Hessian 矩阵，Ñm E 是利用等式（11）或（15）中计算出的梯度。近似 Hesssian 能够表达为：aH = ReJ T J *（18）115其中， J 是 Jacobian 矩阵，上标T 表示转置，* 表示复共轭。由于近似 Hessian 是大规模稀 疏矩阵14，用由 Shin 等（2001）15描述的拟 Hessian 矩阵将梯度进行正则化：wmax ns-1 é ( % T-1)ùT-1 ( )*H a Ñm E =å å ëdiag real F F + l I-wmax i =1û f ADI - DII（19）ë û-1其中，édiag real ( F% T F ) + l I ù所有虚拟源矢量 f 的矩阵。4数值试验是具有衰减因子 l 的拟 Hessian 矩阵对角元素的逆，F 是1201251301351404.1One inclusion 模型先用简单的 One inclusion 模型（图 2）验证上述两种不依赖于子波的 FWI 算法。模型 宽 8km 深 2km，上层背景速度 3000m/s，中间异常球体的速度是 3500m/s，下层强反射层的 速度是 4000m/s，网格间距 20m，x 方向和 z 方向网格点数为 401 和 101。在地表每隔 80m 均匀放炮和布设检波点，每炮 100 道，共 99 炮。在频域用声波方程的有限差分法进行正演 模拟，所用子波为脉冲子波。反演采用从低频到高频连续反演的策略，共 7 个频率，分别为4Hz，5Hz，7 Hz，10 Hz，13 Hz，16 Hz，20 Hz，每个频率反演 12 次，共反演 84 次。 选择图 1 所示模型为反演初始模型。常规的 FWI（图 3）采用已知的精确脉冲震源子波，而反褶积法和褶积法 FWI 中都采用与实际子波类型不同的 Ricker 子波，结果如图 3，图 4， 图 5 所示。可以发现，对于两种不依赖于震源子波的全波形反演方法，尽管在反演迭代过程 中每步正演模拟中用到的 Ricker 子波与实际数据震源子波完全不同，但得到的反演结果（图4 和图 5）与常规的已知震源子波的 FWI 反演结果（图 3）基本一致。图 1 初始速度模型图 2 真实速度模型Fig. 1 Initial velocity modelFig. 2 True velocity model图 3 常规 FWI 反演结果Fig. 3 Inverted velocity model by conventional FWI图 4 反褶积法 FWI 反演结果Fig. 4 Inverted velocity model by source-independent FWI using deconvolution145150155160图 5 褶积法 FWI 反演结果Fig. 5 Inverted velocity model by source-independent FWI using convolution4.2Overthrust 模型再用复杂的 Overthrust 模型验证两种步依赖于子波 FWI 算法。如图 7 所示，Overthrust 模型长 20km 深 4.65km，网格间距 25m，网格大小为 801*187。在地表每隔 100m 均匀放炮 和布设检波点，每炮 200 道，共 199 炮。利用 Ricker 子波，通过频域用声波方程有限差分 法合成数据。反演采用从低频到高频连续反演的策略，共 7 个频率，分别为 3.5Hz，4.7Hz，7.2Hz，9.6 Hz，13.3 Hz，16.9 Hz，7.2Hz，20.6 Hz，每个频率反演 12 次，共反演 84 次。选 择图 6 所示的模型为初始模型。常规的 FWI（图 8）中采用的子波与合成地震数据时完全一 样的 Ricker 子波，而反褶积法和褶积法 FWI 都采用与合成地震数据时完全不同的脉冲子波。可以发现，两种不依赖于震源子波的全波形反演方法得到的反演结果（图 9 和图 10） 与常规的已知震源子波的 FWI 反演结果（图 8）基本一致。从抽取的 X=7.5km 和 X=6.25km 的速度剖面上（图 11 和图 12），也可以得到同样的结论。图 6 初始 overthrust 模型图 7 真实 overthrust 模型Fig. 6 Initial overthrust modelFig. 7 True overthrust model图 8 常规 FWI 反演结果Fig. 8 Inverted overthrust model by conventional FWI图 9 反褶积法 FWI 反演结果结果Fig. 9 Inverted overthrust model by source-independent FWI using deconvolution165170175图 10 褶积法 FWI 反演结果Fig. 10 Inverted overthrust model by source-independent FWI using convolution图 11 三种 FWI 方法反演的 X=7.5km 处速度与真实速度对比Fig. 11 Comparison of true and inverted velocity profile at X=7.5km.图 12 不同 FWI 方法反演的 X=6.25km 处速度与真实速度对比Fig. 12 Comparison of true and inverted velocity profile at X=6.25km.4.3Marmousi 模型最后用复杂的 Marmousi 模型验证两种算法。如图 14 所示，模型长 12.22km 深 2.98km， 网格间距 20m，网格大小为 612*150。在地表每隔 100m 均匀放炮和布设检波点，每炮 200180185190195道，共 199 炮。使用脉冲子波，通过频域声波方程的有限差分法合成地震数据用于地震反演。反演采用从低频到高频连续反演的策略，共 7 个频率，分别为 3.5Hz，4.7Hz，7.2Hz，9.6 Hz，13.3 Hz，16.9 Hz，7.2Hz，20.6 Hz，每个频率反演 12 次，共反演 84 次。选择图 13 的模型 作为初始模型。常规的 FWI 中采用与合成地震数据完全相同的脉冲子波，而反褶积法和褶 积法 FWI 都采用与合成地震数据完全不同的 Ricker 子波。可以发现，两种不依赖于震源子波的全波形反演方法得到的反演结果（图 16 和图 17） 与常规的已知震源子波的 FWI 反演结果（图 15）基本一致。从抽取的 X=7km 和 X=8km 的 速度剖面上（图 18 和图 19），也可以得到同样的结论。图 13 初始 Marmousi 模型图 14 真实 Marmousi 模型Fig. 13 Initial Marmousi modelFig. 14 True Marmous model图 15 常规 FWI 反演结果Fig. 15 Inverted Marmousi model by conventional FWI图 16 反褶积法 FWI 反演结果Fig. 16 Inverted Marmousi model by source-independent FWI using deconvolution图 17 褶积法 FWI 反演结果Fig. 17 Inverted Marmousi model by source-independent FWI using convolution图 18 不同 FWI 方法反演的 X=7km 处速度与真实速度对比Fig. 18 Comparison of true and inverted velocity profile at X=7km2005结论图 19 不同 FWI 方法反演的 X=8km 处速度与真实速度对比Fig. 19 Comparison of true and inverted velocity profile at X=8km205210通过构建归一化的目标函数，使得全波形反演的梯度不依赖于震源子波，从而尽量削弱 和降低因估计震源子波不准对反演造成的不利影响。反褶积法中的归一化，能在梯度的表达式中完全消除子波的影响，但频谱相除有可能会 造成计算上的不稳定。而卷积法中的归一化是广义的归一化，梯度中涉及到目标函数构建的 公共子波，只是在梯度计算中削弱了因对震源子波可能不正确的估计带来的影响，由于只涉 及到波场频域相乘，计算上比起反褶积法会更加稳定。文中所述方法在反演过程中能有效减少对子波的依赖，但参考道的选取非常重要。通过 本文的试验研究发现，选择每炮的最近偏移距对应的地震道为参考道，反演效果最好。数值 试验的结果比起常规的 FWI 差不多，证明了算法的可靠性。参考文献 (References)2152202252302352401 Pratt, R.G., Shin, C.S., and Hicks, G.J. Gauss-Newton and full Newton methods in frequency-space seismic waveform inversionJ. Geophysics. J. Int., 1998, 133: 341-362.2 Lailly, P., The seismic inverse problem as a sequence of before stack migrationA. Bednar, J. B. Conference on inverse scattering: Theory and applicationC. Oklahoma: Soc. Ind. Appl. Math., 1983. 206-220.3 Tarantola, A. Inversion of seismic reflection data in the acoustic approximationJ. Geophysics, 1984, 49:1259-1266.4 Kolb, P., Collino, F., and Lailly, P. Prestack inversion of a 1-D mediumJ. Proceedings of the IEEE, 1986, 74:498-508.5 Gauthier, Odile, Viriux, Jean and Tarantola, A. Two dimensional nonlinear inversion of seismic waveform: Numerical resultsJ. Geophysics, 1986, 51: 1387-1403.6 Mora, P.Nonlinear two-dimensional elastic inversion of multi-offset seismic dataJ. Geophysics, 1987, 52:1211-1228.7 Mora, P. Inversion=migration + tomographyJ. Geophysics, 1989, 54: 1575-1586.8 Sun, R. and McMechan, G. A. Nonlinear reverse-time inversion of elastic offset vertical seismic profile dataJ. Geophysics, 1988, 53: 1295-1302.9 Sun, R. and McMechan, G. A. Full-wavefield inversion of wide-aperture SH and Love wave dataJ. Geophysical Journal International, 1992, 111: 1-10.10 Xu, T., McMechan, G. A. and Sun, R. 3-D prestack full-wavefield inversionJ. Geophysics, 1995, 60:1805-1818.11 Pratt, R. G., and Worthington, M. H. Inverse theory applied to multi-source cross-hole tomography I: Acoustic wave-equation methodJ. Geophys. Prosp., 1990, 38: 287-310.12 Lee, K. H., and H. J. Kim. Source-independent full-wave form inversion of seismic dataJ. Geophysics, 2003,68: 2010-2015.13 Zhou, B., and S. A. Greenhalgh. Crosshole seismic inversion with normalized full-waveform amplitude dataJ. Geophysics, 2003, 68: 1320-1330.14 Chavent, G. and Plessix, R. An optimal true-amplitude least-squares prestack depth-migration operatorJ.Geophysics, 1999, 64: 508-515.15 Shin, C.S., Yoon, K.J., Marfurt, K.J., Park, K., Yang, D., Lim, H.Y., Chung. S., and Shin, S. Efficient calculation of a partial-derivative wavefield using reciprocity for seismic imaging and inversionJ. Geophysics,2001, 66: 1856-1863.