1 引言
地震波形反演利用地震波场反演介质的弹性性质,即在利用波动方程描述地震波传播过程的基础上,从地震波形包含的信息中获取介质的性质(Pica et al., 1990),其目的是利用观测数据与模拟数据的最佳匹配建立地下地质模型.由于充分利用了地震记录所包含的全部信息,因此该方法可以有效地 提高反演精度和分辨率(Tarantolla,1984; Zhou and Gerard, 1993). 20世纪80年代以来,很多学者对基于最小二乘目标泛函的波形反演方法进行了研究(Woodward,1992; Devanvey,1984; Wu and Toksoz, 1987; Pratt and Goulty, 1991; Schuster and Aksel, 1993).Tarantolla(1984)导出可以通过剩余波场反向传播与正传波场的互相关,计算反问题的负梯度方向,进而实现波形迭代反演,这一发现推动了波形反演研究的快速发展;Mora(1987)利用弹性波波形反演技术,理论上解决了更为复杂的地质构造反问题;Zhou和Gerard(1993)、Pratt等(1998)和吴永栓等(2006)将波形反演应用于井间成像,取得了较好的效果;Pratt等(1998)在频率域使用高斯牛顿方法进行全波形反演,推动了频率域波形反演的发展.近年来,伴随状态法被广泛用于FWI中梯度方向的求取(Plessix,2006).
准确的近地表速度结构对静校正、基准面延拓、起伏地表偏移成像等至关重要.传统的射线走时层析成像方法被广泛用于表层速度结构的反演,但受高频射线理论的制约,该方法所得到的反演结果缺乏高频地质信息.刘玉柱等(2009)基于有限频理论(Dahlen et al., 2000)提出了菲涅尔体地震走时层析成像方法反演表层速度结构.该方法基于波动理论导出层析核函数,利用的信息为初至波主能量走时,理论模型与实际资料处理结果表明该方法可以有效提高表层的反演精度与分辨率.Bae等(2012)对比了Laplace全波形反演与菲涅尔体走时层析的核函数,指出Laplace全波形反演虽然重点利用的仍是初至波的走时信息,但可以获得更深部的速度结构.Sheng等(2006)指出,虽然初至波波形基本不包含深部速度信息,但利用其却可以获得精细的近地表速度结构,同时初至波波形反演可以克服全波形反演对初始模型的强烈依赖性,缓解波形反演的多解性与不稳定性.但Sheng等(2006)利用带阻尼的波动方程模拟初至波场,与传统的全波形反演方法相比,正演计算量并没有降低.
波形反演的每次迭代过程中都需要计算理论波场,因此正演计算量非常大.刘玉柱(2011)指出,可以利用高斯束实现格林函数与理论波场的高效计算,当然也可以合成初至地震记录.因此,为了提高初至波波形反演的计算效率,本文利用高斯束作为正演计算的工具.但利用高斯束难以结合伴随状态法计算目标函数的方向与步长.Woodward(1992)分别基于Born近似和Rytov近似导出了两种核函数,称为Born波路径和Rytov波路径,该核函数后来被成功地应用于波动方程层析、有限频层析(Dahlen et al., 2000)与菲涅尔体层析(Spetzler and Snieder, 2004; Liu et al., 2009).但由于存储核函数需要占用大量内存,其一直没有被应用于波形反演.为此,本文在反演部分提出了一种基于Born波路径的矩阵分解算法以实现方向与步长的快速计算,同时避免了庞大核函数的存储.
将该方法应用于理论模型实验,并与声波方程全波形反演和传统的射线初至波走时层析成像方法相对比.实验结果表明,本文提出的基于Born波路径的高斯束初至波波形反演方法能够有效地获得近地表精细速度结构.
2 反演方法在波形反演中,通常定义一目标函数如下:
其中,v为速度,u o是观测数据,u c是满足正演算子的理论合成波场.上述反演问题的目标即找到使目标函数Φ(v)达到极小的解v*.根据最优化理论,通常采用迭代方法获得目标函数Φ(v)的极小值.每次迭代过程中速度的更新量Δv可以表示为一方向向量 p 与步长t的乘积.方向向量 p 可以直接或间接通过求解目标函数Φ(v)的梯度获得.常用的方向和步长如下所示:
其中,p 1、 p 2分别表示一阶方向和二阶方向,p 为任意方向,可以是 p 1或 p 2.t1、t2为对应的一阶步长和二阶步长.Δ u = u o- u c表示残差波场.Ha是目标函数Φ(v)关于模型参数v的近似二阶导数,称为近似Hessian矩阵.∂uc/∂v是Fréchet导数,它反映了速度扰动对波场扰动的影响.T代表复矩阵的共轭转置,-1表示矩阵的逆.利用(1)式直接求方向与步长是非常困难的,因 为Fréchet导数与Hessian矩阵往往难以计算.Woodward(1992)基于声波方程Born近似给出了炮检对(g,s)所对应的波场扰动与速度扰动之间的近似线性关系,
其中,g代表检波点的位置,s代表炮点位置,r表示空间任意一点.v0是扰动前的介质速度,Δv是扰动速度.G0(g,r)是无扰动介质中r点产生的g点处的格林函数,u0(r,s)是由s点产生的r点处的无扰动波场.需要指出的是,Born近似属于一阶近似,即只对波场进行一阶展开,这也就意味着Born近似无法描述散射体的高阶散射,因此将其应用于反演需要基于一个较好的初始模型,且采用迭代策略,以逼近介质与波场之间的非线性关系.根据(6)式可以容易得到描述波场扰动与速度扰动之间一阶线性关系的核函数表达式,
注意,标量k是Fréchet矩阵 K 中炮检对(g,s)对应的核函数在位置r处的分量.将(7)式代入(2)式得到
于是方程(2)—(5)可以重写为如果用m表示数据的个数,n表示模型参数的个数,那么 p 是n×1的列向量,Δ u 是m×1的列向量,K 是一个m×n的矩阵.可见,在实际情况下矩阵 K 是难以存储的,因此很难实现 p 1与 p 2的直接计算.Virieux和Operto(2009)利用伴随状态法(Plessix,2006)实现了一阶方向的间接求解,如(14)式所示
其中,u t表示正传波场,∂ A /∂v是辐射模式,A -1Δ u *为反传的残差波场,*代表共轭算子.这样就避免了 K 的计算与存储.目前,伴随状态法已成为波形反演中计算方向向量的常用方法. 3 高斯束合成初至波场利用初至波波形反演可以实现近地表精细速度结构反演.同时,与全波形反演方法相比,初至波波形反演对初始模型的依赖性与多解性较弱,反演稳定性强(Sheng et al., 2006).但波形反演的每一次迭代过程中都需要计算理论波场,为了减少初至波波形反演的正演计算量,本文拟利用高斯束实现初至波波场与格林函数的快速计算.高斯束基于高频射线理论,其计算精度只有在高频情况下才能得到保证,随着频率的降低其计算误差将会增加,但在反演分辨率可以接受的情况下采用射线理论实现理论波场与格林函数的快速计算是可行的(赵崇进等,2012).
既然高斯束是弹性动力学方程在中心射线附近的高频渐进解,它同样代表了中心射线附近的高频格林函数.为了在得到射线路径上格林函数的同时得到傍轴格林函数,本文采用赵崇进等(2012)给出的高频格林函数公式
其中,T(s)表示波沿中心射线的传播时间,向量 q 为中心射线坐标系下的傍轴坐标,s为沿中心射线距离初始点的弧长,M(s)为走时场在中心射线坐标系下相对于空间坐标的二阶偏导数.为了验证高斯束正演的效果,将一背景速度为3500 m·s-1的均匀模型离散为201×201个网格,网格大小为25 m×25 m.在该模型中进行高斯束正演,模拟结果如图 1、2所示.
由于本文采用高斯束方法计算理论波场与格林函数,因此难于利用伴随状态法(Plessix,2006)实现初至波波形反演中梯度的计算.为了提高计算效率,减少内存占用量,本文基于Born波路径提出了如下矩阵分解算法实现本文反演方法中方向与步长的直接计算.
利用高斯束模拟得到G0(r,g)和u0(r,s)后,一阶方向 p 1即可根据公式(7)、(10)利用矩阵列分解(16)式求取.该矩阵向量乘算法虽然形式简单,但分解后矩阵 K T中的每一列(即方框标注的部分)具有明确的物理含义,即为某个单一炮检对所对应的核函数的共轭.因此,方向 p 1的求取不需要事先计算并存储整个 K 矩阵,一次只需要计算并存储一个核函数,然后通过如下累加计算就可以实现,这样可以大大节省内存占用量,使得Woodward(1992)提出的FWI框架能够应用于实际.
同样地,t2中 Kp 的计算可以直接通过行分解求取,如(17)式所示, 式中,矩阵 K 的每一行(即方框标注的部分)为某个单一炮检对所对应的核函数.理论上,二阶方向 p 2的计算难度较大,目前广泛采用L-BFGS方法(Nocedal,1980)实现其近似计算,但误差较大.Métivier等(2012)提出了迭代求取二阶方向的反演算法,但每一次迭代过程中都需要额外进行二次正演计算,因此计算量非常大.采用本文矩阵分解算法同样可以实现二阶方向的迭代反演,且每一次迭代过程中不需要额外正演计算,因此计算量非常小,关于该方法的具体实现作者将另文撰写.由于Hessian矩阵的存储与精确求逆仍未解决,同时Hessian矩阵理论上是主对角线占优的,因此目前还广泛采用对角Hessian近似计算二阶方向 p 2,该方向称为对角拟牛顿方向(Choi et al., 2007).物理上,主对角线元素代表了观测系统对地下介质的照明.参考(3)式或(11)式不难理解,这样做本质上是对一阶方向做照明补偿,因此这种实现方式又称为一阶方向的预条件.本文采用矩阵分解也可以实现对角近似Hessian H a的计算,
式中,每个对角阵都可根据某一单一炮检对所对应的核函数(即方框标注的部分)算得.因此,对角近似Hessian H a的求取同样不需要事先计算并存储整个 K 矩阵,一次只需要计算并存储一个核函数,然后通过(18)式的累加计算即可以实现,然后将其代入式(3)或(11)即可实现对角拟牛顿方向的计算.利用矩阵分解得到一阶方向 p 1和二阶方向 p 2后,一阶步长t1和二阶步长t2便可以利用(12)、(13)式简单求得.
5 理论模型实验为了验证本文反演方法的有效性,基于Inclusion1二维理论模型(Operto等,2007)进行了初至波射线走时层析、高斯束初至波波形反演和伴随状态法频率域声波方程全波形反演的对比实验.理论模型如图 3所示,模型离散为201×201个网格 点,网格大小为25 m×25 m,背景速度为3500 m·s-1,异常体速度为3700 m·s-1.为了排除观测角度的影响,84个炮点和84个检波点均匀布设于模型内侧距模型边界500 m处,因此炮间距和道间距均为200 m.两种波形反演方法采用一阶方向和一阶步长,共用10个频率成分,初始频率为2 Hz,频率间隔为2 Hz,每个频率迭代10次.初至波射线走时层析与频率无关,为了便于对比,迭代次数设为20.三种方法都以背景速度场作为初始模型,最终反演结果如图 4—6所示.为了定量对比三种反演方法在该模型上的反演效果,将z=2500 m与x=2500 m处的速度剖面提取出来,如图 7所示.
从图 3—7可以看出,三种方法都得到了较好的反演结果,这是由于模型比较简单,且观测系统覆盖均匀.但仔细观察仍然可以发现,由于初至波波形反演方法只利用了初至波波形,而且用高斯束代替声波方程进行理论波场的计算,所以反演精度与分辨率略低于全波形反演结果,但优于射线走时层析反演结果.表 1为三种反演方法的计算效率与内存占用量对比,从中不难看出本文提出的高斯束初至波波形反演方法的优势.需要注意的是,与频率域声波方程全波形反演方法相同,频率域初至波波形反演方法的计算量正比于所利用的频率数和每个频率的迭代次数,即利用的频率越多、迭代次数越多,计算量越大,但反演精度也越高.因此,将射线走时层析与高斯束初至波波形反演的迭代次数归一化后,两者的计算效率基本相当.
为了验证该方法处理复杂模型的有效性,本文基于Overthrust模型(Operto et al., 2007)进行了对比实验.理论模型如图 8所示,该模型被离散为801×187个网格点,网格大小为25 m×25 m.199个炮点均匀分布于模型表面以下25 m处,炮间距为100 m.200个检波点均匀分布于地表,道间距为100 m.第一炮位于水平方向50 m处,第一个检波点位于水平方向25 m处.为了提高反演精度与收敛速率,两种波形反演方法采用二阶方向和二阶步长进行反演计算.共用10个频率成分,初始频率为2 Hz,频率间隔为2 Hz,每个频率迭代10次.初至波射线走时层析与频率无关,为了便于对比,迭代次数设为20.三种方法所采用的初始模型是理论模型的高斯平滑结果(如图 9所示),反演结果如图 10— 12所示.为了进行定量对比,提取的地表以下500 m、 750 m、1000 m处的速度剖面如图 13所示.
通过图 8、图 10—13的对比可知,基于地表观测数据,采用高斯束初至波波形反演方法,可较准确地恢复出模型表层精细速度结构,但由于只利用了初至波波形信息,其在浅部的反演效果略劣于全波形反演,深部的反演效果差的较多.但与传统的初至波走时层析成像方法相比,高斯束初至波波形反演方法 在没有增加太多计算量的情况下(见表 2),不仅获得了明显优于前者的浅部反演结果,而且揭示出更深处的速度结构,这也体现了波形反演的优势.从收敛曲线图 14可知,高斯束初至波波形反演收敛稳定,甚至在高频段比全波形反演具有更好的收敛性.
初至波波形反演由于只利用了初至波波形信息,因此相对于全波形反演它对初始模型的依赖性更弱(Sheng et al., 2006),这也是本文方法的另一个优势.为了验证这一点,本文基于Overthrust理论模型采用大尺度高斯平滑得到了另一个较差的初始模型(如图 15所示),可见该初始模型已非常接近一维梯度模型.基于该初始模型的初至波波形反演 与全波形反演结果如图 16、17所示,地表以下500 m 处的速度剖面定量对比如图 18所示.可见,在初始模型较差的情况下,全波形反演迅速收敛到了一个局部极值,离真解非常远,而初至波波形反演结果在浅部仍具有较好的纵向与横向分辨率,但与图 11相比精度也有所降低.
6 结论
本文提出了基于Born波路径的高斯束初至波 波形反演方法.该方法利用观测数据的初至波波形信息,正演部分基于高斯束实现格林函数和空间波场的快速计算,反演部分提出基于Born波路径的矩阵分解算法降低波形反演对内存的占用量.理论模型实验对初至波射线走时层析、高斯束初至波波形反演与声波全波形反演方法的对比表明,在没有增加太多计算量的情况下,本文提出的方法不仅可以获得较传统射线初至波走时层析更精确、更深的地下速度结构,而且表层速度反演质量与声波方程全波形反演方法相当.同时,该方法计算效率高,内存占用少,对初始模型依赖程度低,计算稳定.该方法可以方便地扩展到三维,并且上述优势表明其在三维实际资料处理中具有很大的应用潜力.
[1] | Bae H S, Pyun S, Shin C, et al. 2012. Laplace-domain waveform inversion versus refraction-traveltime tomography. Geophysical Journal International, 190(1): 595-606. |
[2] | Choi Y, Shin C, Min D J. 2007. Frequency-domain elastic full waveform inversion using the new pseudo-Hessian matrix: elastic Marmousi-2 synthetic test. SEG Expanded Abstract, 1908-1912. |
[3] | Dahlen F A, Hung S H, Nolet G. 2000. Fréchet kernels for finite frequency traveltimes-I. Theory. Geophysical Journal International, 141(1): 157-174. |
[4] | Devanvey A J. 1984. Geophysical diffraction tomography. IEEE Transactions on Geoscience and Remote Sensing, GE-22(1): 3-13. |
[5] | Liu Y Z. 2011. Theory and applications of Fresnel volume seismic tomography [Ph. D. thesis] (in Chinese). Shanghai: Tongji University. |
[6] | Liu Y Z, Dong L G, Wang Y W, et al. 2009. Sensitivity kernels for seismic Fresnel volume tomography. Geophysics, 74(5): U35-U46. |
[7] | Liu Y Z, Dong L G, Li P M, et al. 2009. Fresnel volume tomography based on the first arrival of the seismic wave. Chinese J. Geophys. (in Chinese), 52(9): 2310-2320. |
[8] | Métivier L, Brossier R, Virieux J, et al. 2012. Toward gauss-newton and exact newton optimization for full waveform inversion. EAGE Expanded Abstract, P016. |
[9] | Mora P. 1987. Nonlinear two-dimensional elastic inversion of multi-offset seismic data. Geophysics, 52(9): 1211-1228. |
[10] | Nocedal J. 1980. Updating quasi-Newton matrices with limited storage. Mathematics of Computation, 35(151): 773-782. |
[11] | Operto S, Virieux J, Sourbier F. 2007. Documentation of FWT2D program (version 4.8): Frequency-domain full-waveform modeling/inversion of wide-aperture seismic data for imaging 2D acoustic media. Technical report N0007-SEISCOPE project. |
[12] | Pica A, Diet J P, Tarantola A. 1990. Nonlinear inversion of seismic reflection data in a laterally invariant medium. Geophysics, 55(3): 284-292. |
[13] | Plessix R E. 2006. A review of the adjoint-state method for computing the gradient of a functional with geophysical applications. Geophysical Journal International, 167(2): 495-503. |
[14] | Pratt R G, Goulty N R. 1991. Combining wave-equation imaging with travel-time tomography to from high-reso1ution images from cross-hole data. Geophysics, 56(2): 208-224. |
[15] | Pratt R G, Shin C S, Hicks G J. 1998. Gauss-Newton and full Newton methods in frequency space seismic waveform inversion. Geophysical Journal International, 133(2): 341-362. |
[16] | Schuster G T, Aksel Q B. 1993. Wave-path Eikonal travel-time inversion: Theory. Geophysics, 58(9): 1314-1323. |
[17] | Sheng J M, Leeds A, Buddensiek M, et al. 2006. Early arrival waveform tomography on near-surface refraction data. Geophysics, 71(4): U47-U57. |
[18] | Spetzler G, Snieder R. 2004. The Fréchet volume and transmitted waves. Geophysics, 69(3): 653-663. |
[19] | Tarantola A. 1984. Inversion of seismic data in acoustic approximation. Geophysics, 49(8): 1259-1266. |
[20] | Virieux J, Operto S. 2009. An overview of full-waveform inversion in exploration geophysics. Geophysics, 74(6): WCC1-WCC26. |
[21] | Woodward M J. 1992. Wave-equation tomography. Geophysics, 57(1): 15-26. |
[22] | Wu R S, Toksoz M N. 1987. Diffraction tomography and multisource holography applied to seismic imaging. Geophysics, 52(1): 11-25. |
[23] | Wu Y S, Cao H, Chen J G, et al. 2006. Crosswell seismic waveform inversion and its application. Progress in Exploration Geophysics (in Chinese), 29(5): 337-340. |
[24] | Zhao C J, Liu Y Z, Yang J Z. 2012. Realization of Fresnel volume tomography by Gaussian beam. Oil Geophysical Prospecting (in Chinese), 47(3): 371-378. |
[25] | Zhou C X, Gerard T. 1993. Waveform inversion of sub well velocity structure. SEG Expanded Abstracts, 106-109. |
[26] | 刘玉柱, 董良国, 李培明等. 2009. 初至波菲涅尔体地震层析成像. 地球物理学报, 52(9): 2310-2320. |
[27] | 刘玉柱. 2011. 菲涅尔体地震层析成像理论与应用研究[博士论文]. 上海: 同济大学. |
[28] | 吴永栓, 曹辉, 陈国金等. 2006. 井间地震波形反演与应用. 勘探地球物理进展, 29(5): 337-340. |
[29] | 赵崇进, 刘玉柱, 杨积忠. 2012. 基于高斯束的初至波菲涅尔体地震层析成像. 石油地球物理勘探, 47(3): 371-378. |