2. 海洋国家实验室海洋矿产资源评价与探测技术功能实验室, 青岛 266071
2. Laboratory of Marine Mineral Resources Evaluation and Exploration Technology, Qingdao 266071, China
全波形反演(FWI)作为一种非线性反演方法,通过优化迭代的方式使模拟数据和观测数据之差最小化来估计地下物性参数(Lailly, 1983; Tarantola, 1984; Pratt et al., 1998).由于该方法基于双程波动方程,可以模拟所有波形(包括层间多次波),因此相比于常规基于射线和单程波动方程的成像方法具有更高的成像精度.然而,目前FWI并未得到广泛应用,特别是在高密度地震采集中,巨大的计算量成为一个瓶颈.
目前,编码技术是解决多炮高效成像和建模的有效方法之一.Romero等(2000)在叠前深度偏移中使用了相位编码方法,并测试了多种编码方式的效果.Krebs等(2009)将随机相位编码策略应用于时间域FWI,在每次迭代中随机生成新的编码矩阵以压制非相干正、反传残差波场产生的串扰噪音.Ben-Hadj-Ali等(2011)在频率域FWI中测试了随机相位编码的反演效果,编码后对随机噪音更为敏感并且敏感程度取决于组合炮的数量.Anagaw和Sacchi(2014)研究发现在每次迭代中改变编码函数和震源位置可以高效地压制串扰噪音.此外,使用不同的射线参数将炮集记录转化为平面波记录也可以显著降低计算量,适用于海洋观测系统且对随机噪音不敏感(Liu et al., 2006; Tao and Sen, 2013).Vigh和Starr(2008)将时间域平面波FWI应用于合成数据和实际3D海洋数据反演,而Tao和Sen(2013)表示在频率域可以通过相移合成平面波记录,并在实现2D频率域平面波FWI的基础上分析了射线参数的选取对反演结果的影响.与上面不同的是,Godwin和Sava(2010)提出了一种确定性编码方法称为截断奇异值编码(TSV)并用矩阵的形式表示了常规和联合震源成像的过程,对加窗后的近似单位矩阵进行奇异值分解以获得编码矩阵.
除编码策略外,一些学者从其他方面对减少FWI计算量进行了有益的研究.Wang等(2013)将L-BFGS方法和不精确牛顿方法相结合更高效地估计Hessian算子的逆,从而加快收敛以减少迭代次数.曲英铭等(2016)在时间域将双变网格技术应用于FWI,在不同区域采用不同网格间距,并采取局部变时间采样来减少计算量.Van Leeuwen和Herrmann(2013)提出了一种分批次增加震源的反演方法,反演初期参与的震源较少,具有较高的迭代效率,随着反演的进行,与常规优化方法越来越相近,因而收敛速度加快.2D频率域FWI中,直接解法(如LU分解)由于具有求解多炮的优势因而得到广泛应用,其计算量主要集中在大型稀疏矩阵的分解上.Abubakar等(2009)将有限差分对比源反演(FDCSI)方法应用于频率域FWI用以减少矩阵分解的次数.该方法只需在每个频率开始时进行一次阻抗矩阵分解,其后的迭代中保持背景模型不变,矩阵分解的结果可以重复使用,因而具有很高的迭代效率.He等(2015)利用两步反演策略将FDCSI方法扩展到纵波速度、密度双参数反演.随后通过对多分量弹性波场进行线性变换,形成拟守恒弹性波动方程,将FDCSI方法应用于弹性波方程反演(He et al., 2016).由于FDCSI方法的第一步是计算每炮的对比源,当炮数过多时,其计算效率将受到挑战;其次,该反演方法形式比较特殊,与常规FWI方法相差较大,所以可移植性差,从而制约了其发展.
借鉴FDCSI方法的思想,考虑到方法的计算量与可移植性问题,在常规频率域FWI理论的基础上,本文发展了一种基于波场迭代解法的高效FWI方法.该方法同样只需要在每个频率开始时进行一次阻抗矩阵的分解运算,其后迭代中波场的正传和残差波场反传都是基于共轭梯度法的矩阵回带运算.配合高效的编码策略,理论分析与实际测试都证明了新方法的计算效率要显著高于常规编码FWI.
1 方法原理 1.1 非线性最小二乘反演在频率域,二维各向同性常密度声波方程可以表示为(Marfurt, 1984)
(1) |
其中,r =(x, z)是笛卡尔坐标系下的地下空间位置;rs为震源空间位置;ω是角频率;▽2表示拉普拉斯算子;f(ω)代表频率域的震源项;δ(r - rs)为狄拉克函数;u(r, rs, ω)表示频率域波场;k(r, ω)=ω/v(r)代表波数,v(r)是地下介质的纵波速度.
作为一个最小二乘反演问题,全波形反演正是通过优化的方法使模拟数据和观测数据的误差最小来对模型进行迭代更新(Lailly, 1983; Tarantola, 1984).其基于L2范数的目标函数定义如下:
(2) |
其中,δd(rg, rs, ω)表示检波器位置处观测数据uobs(rg, rs, ω)和模拟数据usyn(rg, rs, ω)之差组成的列向量,rg指检波器所在位置;p(r)是模型参数(如纵波速度);上标t和*分别表示转置和共轭.
目标函数的梯度表达式为
(3) |
式中,Re表示取实部;矩阵J是Fréchet导数矩阵.
当模型参数个数m较大时,直接计算Fréchet导数矩阵的计算量巨大.将Tarantola的残差波场反传理论引入到频率域波形反演,可以高效地得到目标函数对第i个模型参数的梯度(Pratt, 1999;曹书红和陈景波, 2014):
(4) |
其中,A是与模型参数、频率、离散格式以及边界条件有关的阻抗矩阵;
在获取梯度以后,可以使用最速下降法(Choi et al., 2005)或进一步利用共轭梯度法(Kamei et al., 2014)及截断牛顿法(王义和董良国, 2015)等构造优化方向q(r).模型参数的更新迭代公式如下:
(5) |
上式中,n表示迭代次数;步长α可以通过线搜索的方式获得(Gauthier et al., 1986; Pica et al., 1990).
1.2 迭代法生成波场本文将总波场u分解为背景波场ubac和散射场usct之和的形式(u=ubac+usct),并且背景波场满足方程:
(6) |
其中,kb表示背景介质的波数.从方程(1)中减去方程(6)可以得到散射场满足以下方程:
(7) |
上式中,η(r, ω)=kb2(r, ω)-k2(r, ω).为了简化表示,引入线性算子Lb,则散射场可以表示为
(8) |
由此可以定义用于波场迭代的L2型目标函数(为了方便,下文中不再刻意强调坐标位置r、角频率ω等):
(9) |
其中,线性算子Lb由背景模型的阻抗矩阵分解结果确定并且在同一个频率内保持不变.选定需要进行波场模拟的模型后,通过对(9)式最小化可以求得对应波场u,从而避免了对模型阻抗矩阵的重复分解计算(直接解法).
本文采用共轭梯度法对(9)式进行迭代更新,其梯度公式如(10)式所示:
(10) |
(11) |
式中,Lb⊥表示Lb的伴随算子.通过梯度信息构造Polak-Ribière共轭梯度方向w,并令
(12) |
其中,〈 〉表示内积运算.由于步长β可以利用公式(12)显式表达,省去了对步长的估算,从而进一步提高了方法的效率.
上面的推导过程中解决了正传波场的求解问题,而全波形反演中另一个重要部分就是残差波场的反传.用残差波场的反传结果δu替换u,得到下面的目标函数:
(13) |
上式中,
上文中,完成了波场正传以及残差波场反传的计算,线性算子Lb和Lb⊥只与背景模型有关(背景模型保持不变),因此,采用直接解法(LU分解)是一种理想的选择.在这过程中只需对背景模型进行一次阻抗矩阵的分解,其余计算主要是基于共轭梯度法的矩阵乘法与加法运算.
假设离散的未知量个数为N,震源数量为NS,模型迭代次数为Niter,用于生成正传波场和残差反传波场的迭代次数均为niter,LU分解计算量的数量阶为O(2N1.5),向前向后回带的计算量为O(NlogN) (Davis and Duff, 1997).本文采用Polak-Ribière共轭梯度法对模型参数进行更新,并利用抛物插值法估算步长α(韩雨桐等, 2016).由于每次模型更新中至少额外进行两次正演模拟,所以新方法与传统CG方法的计算量分别为
(14) |
(15) |
由于直接法求解波动方程的回带阶段所需要的时间不到LU分解阶段CPU计算时间的5%(Pratt and Worthington, 1990),所以本文将计算量主要集中于回带阶段可以获取较高的计算收益.另外,新方法只需几次迭代就可以合成出高精度的波场(niter较小),再结合高效的编码策略和硬件并行加速方法进一步压缩NS,可以获得更高的反演效率.图 2中绘制了理论计算量随模型参数个数N变化的关系曲线(Niter、niter和NS分别取20、8和1),本文方法的计算量显著低于常规FWI,并且在N较大时,优势更为明显.
为了验证本文提出的波场迭代解法的有效性,采用了简单洼陷模型进行波场正传和残差波场反传测试,并与GMRES方法进行对比以突出高效性.然后将新方法与截断奇异值编码策略(TSV)(Godwin and Sava, 2010)进行了结合,利用Marmousi模型测试了常规编码FWI和niter取不同值时新方法的反演效果.
2.1 迭代法生成波场真实洼陷模型及其初始模型(利用Madagascar的sfsmooth函数平滑所得,平滑窗口为20)如图 3所示,参数如下:水平网格数nx=150,垂直网格数nz=100,水平和垂直网格间距均为10 m.检波器设置为地表全接收,震源位于x=750 m、z=20 m处,震源子波采用主频为10 Hz的雷克子波.
利用图 3a中的真实模型进行LU分解计算得到图 4a中的正传波场(8 Hz实部数据,下同)作为标准波场,同样,图 4b是由图 3b中平滑初始模型进行LU分解计算所得.采用本文的波场迭代方法,以图 4b作为背景波场进行迭代计算,10次迭代后的结果如图 4c所示.图 4d是真实洼陷模型使用GMRES方法求解的波场结果(以图 4b作为初始值),其迭代时间与图 4c一致(约0.7 s).图 4中的第3行(e,f,g)是第2行(b,c,d)分别与图 4a作差的结果,从中可以看出GMRES方法迭代以后的波场残差(图 4g)相比于图 4e有所减小,而本文方法迭代后的残差(图 4f)趋近于零,与标准波场基本一致.
图 5a是检波器接收到的残差波场在真实洼陷模型中的反传结果(8 Hz实部数据,下同),而图 5b是在初始模型中的反传结果.此外,图 5中图形的排列方式与图 4完全一致,可以看出经过迭代以后,GMRES方法的波场残差(图 5g)明显大于本文方法(图 5f).由此说明,在提供一个背景波场的前提下,借助背景模型的阻抗矩阵分解结果,新方法可以高效地获得已知速度场的正传和残差反传波场,从而为下一步计算梯度打下基础.
由于本文提出的波场迭代方法需要配合震源编码技术来压缩NS,才能充分发挥新方法的计算效率.而当NS很小时,常用的迭代解法(如:GMRES)就成为了一种可以考虑的选择.为了比较新方法和GMRES方法的迭代效果,在图 4和图 5中对比了两种方法在相同时间下的迭代结果.图 6 (a—b) 分别是正传和反传两种情况下,两种迭代方法的归一化波场误差随迭代时间变化曲线,本文方法的收敛速度明显高于GMRES方法(新方法不包含背景模型矩阵分解时间,因为一个频率反演中只需分解一次).主要原因是新方法是一种基于CG方法的无约束优化方法,具有明确的优化方向和步长值,而GMRES方法需要进行合适的预处理,而这通常是难以满足的(Virieux and Operto, 2009).图 6 (c—d)分别是正、反传时新方法的波场误差随迭代次数变化曲线,大约进行6次迭代后残差基本收敛.
为了测试新方法对初始模型的适应性,分别采用图 7a中平滑程度更高的初始模型(平滑窗口为30)和图 8a中线性梯度初始模型进行波场合成.图 7(b—d)分别是初始模型的模拟波场以及2次和10次迭代后的波场结果,而图 7 (e—g)则分别是图 7 (b—d)与图 4a中标准波场的差值.图 8的排列方式与图 7一致,并且可以看出,新方法对初始模型具有较好的适应性,随着迭代的进行,残差稳定减小.值得注意的是,本次测试中是以初始模型的波场作为迭代的初始值,与真实波场差异较大,实际反演中相邻迭代间模型差异较小,并以上一次迭代的波场作为初始值,因此需要的迭代次数通常更少.
为了验证本文方法的反演效果,采用Marmousi模型(Anagaw and Sacchi, 2014)进行合成数据测试.所使用的真实和初始模型如图 9所示(初始模型由真实模型平滑所得),网格大小为126×384,横、纵向空间采样间隔均为15 m.震源和检波器位于地表,第一炮位于15 m处,炮间距为30 m,共191炮,检波器设置为地表全接收,震源子波为主频15 Hz的雷克子波.
反演中所使用的频率为4、6、8、10、12、15、18、21、24、28、32和36 Hz共12个频率,每个频率迭代20次,低频的反演结果作为相邻高频的输入.本文使用截断奇异值编码(TSV)方法,将191炮编码为15炮,并采取并行反演策略,图 10是相应的串扰噪音矩阵.常规编码FWI反演结果如图 11a所示,niter分别取10、8和6时,新方法的反演结果对应于图 11(b—d).由于串扰噪音的影响,反演结果的浅部低速层受影响比较明显,但总体上说,图 11反演结果较好并且比较接近.
图 12给出了图 11中4个反演结果的归一化收敛曲线,常规编码FWI和新方法niter等于10时的收敛效果相差不大,并且常规编码FWI的模型误差最小,而新方法随着niter增大收敛性有所改善.图 13记录了反演时间随迭代次数变化曲线,常规编码FWI的计算用时要显著高于本文方法,并且niter越小新方法计算时间越少.本次测试中模型四周设有20个网格的PML边界,实际计算用时与图 2中的理论用时比较吻合,常规编码FWI计算时间约为niter等于8时新方法的8倍(N≈7×104).
本文提出了一种高效的迭代解法生成频率域波场方法,将其应用于声波常密度方程FWI,并给出了详细的反演流程.在实现算法的基础上,通过洼陷模型和Marmousi模型测试,得到如下几点认识:
(1) 新方法将频率域波场计算转化为一种无约束优化问题,通过几次迭代便能生成高精度的正传和残差反传波场,其误差收敛速度显著高于未经预处理的GMRES方法.
(2) 基于新的波场迭代法的FWI,每个频率只需进行一次阻抗矩阵的分解,其计算量与波场迭代次数niter成正比.配合高效的震源编码技术,新方法的计算时间相比常规编码FWI有明显降低(niter=8、N≈7×104时,计算时间约为常规编码FWI的1/8).
(3) niter越大,新方法的反演效果与常规编码FWI越接近,但是当生成的波场精度达到一定条件后,niter的增大对反演结果的改善并不明显.因此,需要在计算量与反演效果之间做好权衡.
虽然本文提出的基于波场迭代法FWI可以获得与常规FWI基本一致的反演效果,并且在计算效率上有显著优势,但是该方法还不是十分完善,也存在一些需要改进的地方,比如:(1)当NS较大时,新方法的计算效率将受到严重影响,因此必须配合高效的震源编码技术,这样不可避免地会引入一些串扰噪音;(2)新的波场迭代法本质上是一种基于直接解法的迭代方法,当模型维度过大时,同样存在存储量的问题.
Abubakar A, Hu W Y, Habashy T M, et al. 2009. Application of the finite-difference contrast-source inversion algorithm to seismic full-waveform data. Geophysics, 74(6): WCC47-WCC58. DOI:10.1190/1.3250203 |
Anagaw A Y, Sacchi M D. 2014. Comparison of multifrequency selection strategies for simultaneous-source full-waveform inversion. Geophysics, 79(5): R165-R181. DOI:10.1190/geo2013-0263.1 |
Ben-Hadj-Ali H, Operto S, Virieux J. 2011. An efficient frequency domain full-waveform inversion method using simultaneous encoded sources. Geophysics, 76(4): R109-R124. |
Cao S H, Chen J B. 2014. Studies on complex frequencies in frequency domain full waveform inversion. Chinese Journal of Geophysics, 57(7): 2302-2313. DOI:10.6038/cjg20140724 |
Choi Y, Shin C, Min D J, et al. 2005. Efficient calculation of the steepest descent direction for source-independent seismic waveform inversion:An amplitude approach. Journal of Computational Physics, 208(2): 455-468. DOI:10.1016/j.jcp.2004.09.019 |
Davis T A, Duff I S. 1997. An unsymmetric-pattern multifrontal method for sparse LU factorization. SIAM Journal on Matrix Analysis and Applications, 18(1): 140-158. DOI:10.1137/S0895479894246905 |
Gauthier O, Virieux J, Tarantola A. 1986. Two-dimensional nonlinear inversion of seismic waveforms:Numerical results. Geophysics, 51(7): 1387-1403. DOI:10.1190/1.1442188 |
Godwin J, Sava P. 2010. Simultaneous source imaging by amplitude encoding, Technical Report CWP-645, Center for Wave Phenomena, Colorado School of Mines.
|
Han Y T, Huang J P, Qu Y M, et al. 2016. Encoded full waveform inversion method by an optimum inversion step length for aliasing data. Progress in Geophysics, 31(1): 198-204. DOI:10.6038/pg20160123 |
He Q L, Han B, Chen Y, et al. 2015. Application of the finite-difference contrast source inversion method to multiparameter reconstruction using seismic full-waveform data. Journal of Applied Geophysics, 124: 4-16. |
He Q L, Chen Y, Han B, et al. 2016. Elastic frequency-domain finite-difference contrast source inversion method. Inverse Problems, 32(3): 035009. DOI:10.1088/0266-5611/32/3/035009 |
Kamei R, Pratt R G, Tsuji T. 2014. Misfit functionals in laplace-fourier domain waveform inversion, with application to wide-angle ocean bottom seismograph data. Geophysical Prospecting, 62(5): 1054-1074. DOI:10.1111/gpr.2014.62.issue-5 |
Krebs J R, Anderson J E, Hinkley D, et al. 2009. Fast full-wavefield seismic inversion using encoded sources. Geophysics, 74(6): WCC177-WCC188. DOI:10.1190/1.3230502 |
Lailly P. 1983. The seismic inverse problem as a sequence of before stack migration. //SIAM Conference on Inverse Scattering: Theory and Applications. Philadelphia: Society for Industrial and Applied Mathematics, 206-220.
|
Liu F Q, Hanson D W, Whitmore N D, et al. 2006. Toward a unified analysis for source plane-wave migration. Geophysics, 71(4): S129-S139. DOI:10.1190/1.2213933 |
Marfurt K J. 1984. Accuracy of finite-difference and finite-element modeling of the scalar and elastic wave equation. Geophysics, 49(5): 533-549. |
Pica A, Diet J P, Tarantola A. 1990. Nonlinear inversion of seismic reflection data in a laterally invariant medium. Geophysics, 55(3): 284-292. DOI:10.1190/1.1442836 |
Pratt R G, Worthington M H. 1990. Inverse theory applied to multi-source cross-hole tomography. Geophysical Prospecting, 38(3): 287-310. DOI:10.1111/gpr.1990.38.issue-3 |
Pratt R G, Shin C, Hick G J. 1998. Gauss-Newton and full Newton methods in frequency-space seismic waveform inversion. Geophysical Journal International, 133(2): 341-362. DOI:10.1046/j.1365-246X.1998.00498.x |
Pratt R G. 1999. Seismic waveform inversion in the frequency domain, Part 1:Theory and verification in a physical scale model. Geophysics, 64(3): 888-901. DOI:10.1190/1.1444597 |
Qu Y M, Li Z C, Huang J P, et al. 2016. Full waveform inversion based on multi-scale dual-variable grid in time domain. Geophysical Prospecting for Petroleum, 55(2): 241-250. |
Romero L A, Ghiglia D C, et al. 2000. Phase encoding of shot records in prestack migration. Geophysics, 65(2): 426-436. DOI:10.1190/1.1444737 |
Tao Y, Sen M K. 2013. Frequency-domain full waveform inversion with plane-wave data. Geophysics, 78(1): R13-R23. DOI:10.1190/geo2012-0267.1 |
Tarantola A. 1984. Inversion of seismic reflection data in the acoustic approximation. Geophysics, 49(8): 1259-1266. DOI:10.1190/1.1441754 |
VanLeeuwen T, Herrmann F J. 2013. Fast waveform inversion without source-encoding. Geophysical Prospecting, 61(S1): 10-19. |
Vigh D, Starr E. 2008. 3D prestack plane-wave, full-waveform inversion. Geophysics, 73(5): VE135-VE144. DOI:10.1190/1.2952623 |
Virieux J, Operto S. 2009. An overview of full-waveform inversion in exploration geophysics. Geophysics, 74(6): WCC1-WCC26. DOI:10.1190/1.3238367 |
Wang Y, Dong L G, Liu Y Z. 2013. Improved hybrid iterative optimization method for seismic full waveform inversion. Applied Geophysics, 10(3): 265-277. DOI:10.1007/s11770-013-0389-x |
Wang Y, Dong L G. 2015. Multi-parameter full waveform inversion for acoustic VTI media using the truncated Newton method. Chinese J. Geophys. (in Chinese), 58(8): 2873-2885. DOI:10.6038/cjg20150821 |
曹书红, 陈景波. 2014. 频率域全波形反演中关于复频率的研究. 地球物理学报, 57(7): 2302-2313. DOI:10.6038/cjg20140724 |
韩雨桐, 黄建平, 曲英铭, 等. 2016. 基于优化步长的混叠数据编码全波形反演方法. 地球物理学进展, 31(1): 198-204. DOI:10.6038/pg20160123 |
曲英铭, 李振春, 黄建平, 等. 2016. 基于多尺度双变网格的时间域全波形反演. 石油物探, 55(2): 241-250. |
王义, 董良国. 2015. 基于截断牛顿法的VTI介质声波多参数全波形反演. 地球物理学报, 58(8): 2873-2885. DOI:10.6038/cjg20150821 |