2. 中国科学院计算数学与科学工程计算研究所, 北京 100190;
3. 同济大学海洋与地球科学学院波现象与反演成像研究组, 上海 200092
2. Institute of Computational Mathematics and Scientific/Engineering Computing, Beijing 100190, China;
3. WPI, School of Ocean and Earth Science, Tongji University, Shanghai 200092, China
勘探地震学的研究本质是一个反问题,研究的客观对象即地球物理介质.从正方面研究,即是研究地球介质中地震波场的传播现象,归纳不同介质中地震波场传播的规律.而从反方面研究,既是根据人工能够观测得到的地震数据(包括地表观测数据、海底观测数据、井中观测数据、井间观测数据等),反向推演波场传播,并估计介质地球物理参数的过程.因此,从地震数据反推地下地球物理参数的过程都可以称为地震反演.地震波反演成像要达到三个层次的目标:地下岩石的几何结构描述,保振幅的处理结果和岩石参数估计.这三个目标依次提高.传统的偏移成像解决的基本是几何结构描述的问题[1].地震反演成像的目标是解决岩石参数估计的问题[2-3].第二个问题:保振幅,需要把偏移成像纳入反演成像来寻找答案,即以地震波方程描绘传播过程的程度和地震反演的理论来统一思考.从偏移成像的角度来看,影响最终成像结果振幅的因素主要有:球面扩散效应、地表一致性问题、波场传播算子的精度和介质的吸收衰减等问题.其中,介质的吸收衰减问题不受声波-弹性波方程所控制,可以用吸收衰减的补偿策略加以矫正[4].其它几个问题的解决可以从改进偏移算子的角度处理[5].但地震成像到了保振幅的阶段必须放在反演的框架中论述,才能找到更为合理的物理意义.最小二乘偏移是反演框架下的偏移结果[6-8],它是对偏移成像过程和结果的校正.这种校正的桥梁就是Hessian算子.Hessian算子中包括了能被波动方程控制的、影响偏移振幅的所有信息.最小二乘偏移的难点就在于求解一个精确的Hessian并求逆.当然,在同样的地震反演框架下,基于牛顿/高斯-牛顿法的全波形反演问题[9]中,Hessian算子具有相同的地位和作用.作为预条件,Hessian算子作用在全波形反演的过程还能够提高反演迭代的效率.问题是,构造一个什么样的泛函格式加入反演框架中才能兼顾计算效率和精度.本文将在前期最小二乘偏移中角度域Hessian的研究基础[7]之上进一步提出两种分别适用于全波形反演和最小二乘偏移的Hessian算子的变换格式.
2 地震反演与Hessian算子地震波场的正传播过程可以记为
![]() |
(1) |
式中,mtrue为地震地球物理参数矢量,如速度、密度、弹性参数等,dobs为观测到的地震数据,L(·)描述了依赖于mtrue的地震波场正传播过程,它表达了地震波场传播的系统.在此系统下,研究特定模型mtrue中地震波的传播过程以得到数据dobs的过程为地震正演.反之,由数据dobs反推模型参数mtrue的过程为地震反演.反演的过程可以记为
![]() |
(2) |
式中,L-1(·)即描述了利用各种数学工具进行地震反演以得到估计的模型参数mest的过程.相对地,地震反演问题比正演问题要复杂,这是由于地震正传播得到的数据只把模型的部分信息带入了数据中.而且,模型和正问题的传播过程都只是实际情况的近似.地震反演理应遵从普遍的数学反演方法规律.法国数学家JacquesHadamard[10]将反演的问题归结为三点:
·存在性.对于给定数据d,存在一个模型m能够满足d=L(m).
·唯一性.只存在一个m.
·稳定性.数据中一个小的扰动只能在解中带来一个微小的扰动.
如果反演问题不能满足上述条件中的至少一个,该问题就是一个病态问题.而在地震反演中,上述条件常常是不能同时满足的.如,由于实际地震观测孔径、频带范围、采集步长等造成数据量不足,使得地震反演的过程是不唯一和不稳定的.再如,观测数据中的噪声,并不受波动方程所控制,这会给反问题的存在性带来挑战.
2.1 牛顿类反演方法地震反演问题本质上是一个非线性问题,即正传播算子随着模型参数非线性变化.但人们在研究非线性问题时,常常在局部甚至全局线性化,以简化问题的求解过程.为求解这个反问题,在二范数意义下定义误差泛函:
![]() |
(3) |
式中,上标T表示矩阵转置.m在声波方程中是速度相关量,在弹性波方程中是弹性波参数的集合.方程(3)定义了基础的最小二乘目标函数,定义Fréchet微商为
![]() |
(4) |
在非线性问题中,Fréchet微商是依赖于模型的.牛顿法认为误差函数在初始模型附近满足二次型,因此可以用泰勒公式将误差函数在模型附近展开:
![]() |
(5) |
式中,
![]() |
(6) |
H即为Hessian算子,它有两部分构成.其中,Δd∂FT/∂m为其非线性项.这是因为对于非线性问题,Fréchet微商是依赖于模型的.而对线性问题这一项退化为0,Hessian算子退化到线性项:
![]() |
(7) |
数学上,Hessian算子就是误差泛函对模型的二次导数.在物理上,我们基于波动方程求取Hessian算子的公式来加以分析.在全牛顿法反演中,精确的Hessian算子包括两项:线性近似项Ha和非线性项R.大部分情况下,由于R本身计算量的需要十分巨大,并且,在局部线性假设的条件下R趋向于0.因此,首先需要研究近似的Ha.但是,单就Ha项而言,Hessian矩阵仍然是一个巨大的矩阵,对它的求逆是高斯-牛顿反演和最小二乘偏移中的一个难点.
从声波方程来看,
![]() |
(8) |
式中,f(xs,ω)为炮点xs处激发的能量在圆频ω上的谱.σ(x)=1/v(x)为介质点x上的慢度,v(x)为介质点集合x上的速度.这里,波场在正演过程中是一个随着慢度变化的量.取v(x)为反演参数,则牛顿反演中,梯度项和Hessian项分别展开为
![]() |
(9) |
![]() |
(10) |
可以看到,Hessian矩阵的元素H(x,y)对应着两组成像点x和y.高斯牛顿法中省略的非线性项R即是方程(10)中右端加号后面的部分.
在局部线性化的Born近似的意义下[11],接收点xr处的合成数据可以写为
![]() |
(11) |
其中,r(x)=σ2(x),把合成波场公式代入误差函数,由式(8),式(9)和式(10)可得
![]() |
(12) |
![]() |
(13) |
偏移成像的过程可以看作是一个偏移算子作用到地震数据上得到模型参数(成像)的过程.相反地,地震数据的产生可以看作是一个正演算子作用到模型上得到地震数据的过程.从这个意义上,波动方程又可以表示为
![]() |
(14) |
式中,L即是正演算子,在地震参数和正演数据之间建立的一个联系.在Born近似的条件下,正演算子的每个元素为
![]() |
(15) |
因此,线性Hessian矩阵就可以表示为
![]() |
(16) |
梯度项则变成:
![]() |
(17) |
而基于牛顿法求解地震反演问题的一般性公式为
![]() |
(18) |
这样,基于线性Hessian矩阵的高斯-牛顿法方程就可以表示为
![]() |
(19) |
假设初始反演参数m(0)=0,则可以得到第一次迭代过程公式为
![]() |
(20) |
这样,就建立了反演参数与正演算子的关系.事实上,一次迭代实现最小二乘偏移成像的过程中高斯-牛顿反演的梯度项为LTd.将地震传播算子的共轭作用到地震数据上,即对地震数据的反传播.因此此梯度项就是直接偏移成像的结果[12].可以记为
![]() |
(21) |
可以把这个介于零次迭代和一次迭代的中间结果称为是1/2次迭代结果.而事实上,这个作用过程就是偏移成像.把LT称为是偏移算子,它是正演算子L的共轭转置.式(20)的迭代过程就是通常说的最小二乘偏移成像.最小二乘偏移成像为真振幅偏移成像打开了一扇大门.这里,线性Hessian矩阵的逆作为加权因子作用到成像结果m(1/2)上.可以看到,Hessian矩阵就包含了振幅修正的意义.
3 Hessian算子的性质与近似Hessian算子在上述反演框架中居于重要位置,但其规模过于庞大(模型参数个数的平方倍元素),无论是在高斯-牛顿法波形反演过程中还是在最小二乘偏移中,人们都需要面对其庞大的求取计算量和更庞大的求逆计算量.因此,许多Hessian矩阵的近似求解方法被提出.Yu等[13]提出了一个基于v(z)介质的Hessian求逆方法;Lecomte[14]提出了一种基于射线理论的Hessian求逆方法;Plessix和Mulder[5]在空间域推导了Hessian矩阵公式,并提出了对角Hessian矩阵的四种近似解.Tang [15]用相位编码(Phase-Encoding)方法快速的计算Hessian矩阵.对于叠前深度偏移,最为直接的就是利用Hessian矩阵的对角特性.从前期的分析可知,线性Hessian矩阵是一个主对角占优的带状矩阵[7],其特征值主要由对角元素贡献.因此,可以只利用线性Hessian矩阵的主对角元素代替整个矩阵以实现不甚精确的反演.当然,如前面所述,这一假设引入的误差可以通过增加迭代次数来弥补.而一次迭代的最小二乘偏移则不能完全达到逆Hessian的效果.线性Hessian矩阵的对角元素为
![]() |
(22) |
(22)式仍然是张在整个模型空间x上的矩阵,它是线性Hessian矩阵的主对角元素的集合,称之为对角Hessian矩阵.可以看到,对角Hessian矩阵是正演算子L(xs,xr,x,ω)的平方.因此可以设想,它能够部分表述波场传播过程中的能量衰减及炮检系统对地下介质的照明程度.
3.1 平面波Hessian算子与平面波最小二乘偏移最小二乘偏移成像,可以很大程度上消除观测孔径、采集脚印、地震波场衰减等因素的影响,改善地震成像的振幅效应.同时,基于反演框架构建的最小二乘偏移公式,能够适用于各种偏移算子.只要提取出相应的正演算子L和共轭的偏移算子L*,则最小二乘偏移中的线性Hessian矩阵和梯度值都可以通过矩阵运算得到求解.
偏移距外推(OffsetContinuation)是地震数据映射方法研究的一个热点问题.Bagaini[16]推导出了炮域的偏移距外推方程.Formel[17]对前人的偏移距外推方程进行了修改,给出了更精确的外推方程.基于偏移距外推的思想,可以对地震波场进行平面波偏移.首先,把炮集数据分解为平面波:
![]() |
(23) |
式中,xm为共中心点(CMP)坐标,h为偏移距,ph为偏移距平面波射线参数.上式为一个单频ω数据u(h,xm;ω)变换到平面波域数据
叠前时间偏移的双平方根算子,在频率-波数域可以表示为
![]() |
(24) |
其中,kx和kh分别为对应共中心点xm和偏移距h的波数.τ为叠前时间偏移中的时间深度.vτ和kτ分别为对应时间深度τ的速度和波数,在水平层状介质假设条件下它们与xm不相关.上式可以被写为
![]() |
(25) |
由文献[18]可知,tanθ=-vτkh/(2kτ).这里角度θ为入射波和反射波之间的夹角.因此有如下关系:
![]() |
(26) |
双平方根算子进一步写为
![]() |
(27) |
而偏移距射线参数ph和θ之间的关系为
![]() |
(28) |
因此,双平方根算子对应的频率-偏移距射线参数域的波动方程变为
![]() |
(29) |
方程(29)可以通过有限差分的方法来求解.这样,就建立了从平面波数据
![]() |
(30) |
最小二乘偏移首先需要构造一个能够矩阵表达的偏移算子.式(23)-(30)通过偏移距外推方法进行平面波偏移的过程能够提取出偏移算子.其过程可以写为
![]() |
(31) |
那么,偏移算子L*(τ,ph;ω)的元素由什么构成呢?从(31)式可以看出,如果构造一个虚拟的地震数据
![]() |
(32) |
可以看到,在此条件下的线性Hessian矩阵反应了特定炮检系统和速度背景条件下不同射线参数之间的相关性.当ph=ph′时,即线性Hessian矩阵的对角线元素Ha(ph,ph)反映了单一射线参数的照明能量.利用上式给定的平面波Hessian,就可以构造平面波最小二乘偏移.
3.2 偏移距Hessian算子与拟牛顿法全波形反演在牛顿法全波形反演框架中,应用Hessian算子作为预条件可以提高全波形反演的计算效率.然而,如何构建与梯度求取相匹配的Hessian成为一个难点.在前期的研究中,与小波域波场传播算子匹配,我们发展了局部角度域Hessian算子[7]:
![]() |
(33) |
这一Hessian近似形式能够成功应用在最小二乘偏移中.但这一基于频率域角度分解提出的Hessian格式较难以应用在全波形反演中.这里,我们提出一种与角度Hessian类似的地下偏移距域Hessian算子:
![]() |
(34) |
由于这一格式不需要在频率波数域提取角度,因此能够方便的与扩展成像条件得到的成像道集[20-21]联合起来,形成迭代格式.我们期望其可以应用在全波形反演中,以加快反演速度.当然,考虑计算效率的提升,还可以将公式(34)中的检波点所对应的两个格林函数略去,可以得到仅利用炮照明的Hessian格式,我们将在模型实验部分对其效果进行测试.
4 基于近似Hessian的偏移与反演 4.1 平面波最小二乘偏移首先在一个较为简单的模型上检验用平面波Hessian进行最小二乘偏移的情况.在图 1a给出的实验速度模型上,在红五角星标识的位置以主频为30Hz的Ricker子波激发生成5炮数据.利用第3.1节的方法对其进行最小二乘偏移实验.选取图中标识的两条虚黄线标识的位置比较其反射系数.图 1b和图 1c中分别给出了直接平面波深度偏移成像结果和平面波最小二乘偏移成像结果.图 1d为所求取的平面波Hessian矩阵(纵横坐标分别为射线参数).可以看到利用平面波Hessian进行的加权使得叠前偏移成像的结果无论是在分辨率和振幅均衡性上均取得改善.地下照明不足的区域得到了明显的补偿,而偏移剖面的分辨率也得到了提高.当然,平面Hessian的图样也显示了该矩阵的强对角性.在较为简单的模型中,完全可以利用对角Hessian取代全Hessian应用到最小二乘偏移成像中.图 1e和图 1f分别给出的是两个参考道反射系数和成像值的对比结果.两图中的第一道为按照速度场沿垂向入射波场的反射系数,第二道为直接偏移得到的成像结果,第三道为最小二乘偏移成像得到的结果.其中,成像结果根据反射系数进行了振幅的归一化以消除量纲差异.两个参考道的对比都显示了照明减少的下面一层反射界面的成像值在最小二乘偏移结果中得到补偿.而成像子波的宽度也显示了最小二乘偏移在成像分辨率上的提高.
![]() |
图 1 平面波最小二乘偏移实验 (a)凹陷模型时间深度域速度模型;(b)平面波PSTM成像结果;(c)平面波LSM成像结果;(d)平面波Hessian矩阵;(e)参考道1反射系数、直接成像结果和LSM成像结果;(f)参考道2反射系数、直接成像结果和LSM成像结果. Fig. 1 Planewave least-squares migration test (a) Time-depth velocity model of Aoxian model; (b) planewave PSTM imaging result; (c) planewave LSM imaging result; (d) planewave Hessian matrix; (e) reflectivity, direct imaging and LSM imaging result of reference trace 1;(f) reflectivity, direct imaging and LSM imaging result of reference trace 2. |
全波形反演广泛地被认为是对地下复杂介质地震成像的强有力的工具,理论上是最为精确的参数估计的方法.然而其巨大的计算量仍给其广泛的应用带来了挑战.采用式(34)中的偏移距Hessian对全波形反演的过程进行加权以测试其对计算效率的改善.检测标准就是误差泛函的收敛效率.利用广为应用的Marmousi模型进行测试.在一个光滑的初始模型基础上(图 2a)分别进行常规的共轭梯度法(CG)全波形反演(FWI)、经炮照明加权的CG全波形反演和经偏移距Hessian加权的全波形反演的计算效率进行测试(图 2b-2d).从误差泛函收敛的曲线(图 2e)可以看出,经偏移距Hessian加权后的全波形反演收敛速度明显快于常规的CG法全波形反演.而只计算偏移距Hessian中的炮波场照明成分进行加权的FWI结果也取得较快的收敛效率.当然,从图中可以看出,偏移距Hessian加权的全波形反演过程中存在一个突变点.这说明了这一方法本身的稳定性问题值得注意.类似问题出现在角度域Hessian进行最小二乘偏移中[7].但图中曲线反映出的总体收敛趋势是稳定可以接收的.
![]() |
图 2 偏移距Hessian用于全波形反演实验 (a)初始模型;(b)常规CG法FWI270次迭代结果;(c)利用炮照明Hessian加权的CG法FWI270次迭代结果;(d)利用偏移距Hessian加权的CG法FWI270次迭代结果;(e)误差收敛曲线.红色为常规CG方法;蓝色为利用炮照明预处理的CG方法;绿色为利用偏移距Hessian矩阵预处理CG方法. Fig. 2 Tests on the offset-Hessian used to full waveform inversion (a) Initial model; (b) Traditional CG-FWI result after 270 iterations; (c) Shot-illuminated CG-FWI result after 270 iterations; (d) Offset-Hessian preconditioned CG-FWI result after 270 iterations; (e) Restraintion curve of misfit function, red line-Traditional; Blue line-Shot illuminated; Green line-Offset hessian preconditioned. |
本文从牛顿类地震反演问题出发,详细分析和讨论了Hessian算子在地震反演成像中的数学和物理意义.在数学上,Hessian算子是误差泛函对模型参数的二阶导数.因此,误差泛函对模型的变化实际上是否为二次型又成为一个问题.前人的研究表明了地震反问题的强非线性,然而在局部,仍可将此问题简化成二次型甚至线性问题.这种简化体现在基于波动方程的Born近似中,也是许多地震反问题如走时层析的基础假设之一.正因为这种假设,本文所讨论的最小二乘偏移和全波形反演都是基于局部寻优的数学方法导出的,因而这些方法依赖于一个较好的初始模型参数.声波方程条件下导出的Hessian算子的格林函数表达形式反映出了炮检系统对模型的照明关系,因此这种关系的逆过程即Hessian算子的逆就是对观测系统照明不均的补偿.因此,一个精确的逆Hessian算子作用到梯度(成像结果)上,理应得出一个保振幅的解,而本文基于平面波偏移提出的平面波Hessian的格式证明了最小二乘偏移在保幅方面的能力.同样,应用Hessian进行二次型的牛顿反演也比基于一次线性搜索的梯度法更加接近地震反演的非线性性质,这才是牛顿类全波形反演具有更快的计算效率的基础.地下偏移距的Hessian算子在全波形反演中的应用证明了这一观点.总之,Hessian算子在地震反演成像中的具有重要的意义,继续深入研究这一算子对保幅偏移和地震反演都具有重要的价值.对计算效率和计算稳定性的进一步提高是本项研究的下一个目标.
致谢感谢匿名评审提出的宝贵修改意见.特别感谢中国石油勘探开发研究院西北分院对本项研究的支持.
[1] | Claerbout J F. Toward a unified theory of reflector mapping. Geophysics , 1971, 36(3): 467-481. DOI:10.1190/1.1440185 |
[2] | Tarantola A. Inversion of seismic reflection data in the acoustic approximation. Geophysics , 1984, 49(8): 1259-1266. DOI:10.1190/1.1441754 |
[3] | 任浩然, 王华忠, 黄光辉. 地震波反演成像方法的理论分析与对比. 岩性油气藏 , 2012, 24(5): 12–18. Ren H R, Wang H Z, Huang G H. Theoretical analysis and comparison of seismic wave inversion and imaging methods. Lithologic Reservoirs (in Chinese) , 2012, 24(5): 12-18. |
[4] | 任浩然, 王华忠, 张立彬. 沿射线路径的波动方程延拓吸收与衰减补偿方法. 石油物探 , 2007, 46(6): 557–561. Ren H R, Wang H Z, Zhang L B. Compensation for absorption and attenuation using wave equation continuation along ray path. Geophysical Prospecting for Petroleum (in Chinese) , 2007, 46(6): 557-561. |
[5] | Etgen J, Gray S H, Zhang Y. An overview of depth imaging in exploration geophysics. Geophysics , 2009, 74(6): WCA5-WCA17. DOI:10.1190/1.3223188 |
[6] | Plessix R E, Mulder W A. Frequency-domain finite-difference amplitude-preserving migration. Geophysics Journal International , 2004, 157(3): 975-987. DOI:10.1111/gji.2004.157.issue-3 |
[7] | Ren H R, Wu R S, Wang H Z. Wave equation least square imaging using the local angular Hessian for amplitude correction. Geophysical Prospecting , 2011, 59(4): 651-661. DOI:10.1111/gpr.2011.59.issue-4 |
[8] | Ribodetti A, Operto S, Agudelo W, et al. Joint ray+Born least-squares migration and simulated annealing optimization for target-oriented quantitative seismic imaging. Geophysics , 2011, 76(2): R23-R42. DOI:10.1190/1.3554330 |
[9] | Pratt G, Shin C, Hicks G J. Gauss-Newton and full Newton methods in frequency-space seismic waveform inversion. Geophysical Journal International , 1998, 133(2): 341-362. DOI:10.1046/j.1365-246X.1998.00498.x |
[10] | Hadamard J. An Essay on the Psychology of Invention in the Mathematical Field. Princeton: Princeton University Press, 1945 . |
[11] | 陈生昌, 曹景忠, 马在田. 稳定的Born近似叠前深度偏移方法. 石油地球物理勘探 , 2001, 36(3): 291–296. Chen S C, Cao J Z, Ma Z T. Stable pre-stack depth migration method with Born approximation. Oil Geophysical Prospecting (in Chinese) , 2001, 36(3): 291-296. |
[12] | Lailly P. The seismic inverse problem as a sequence of before stack migrations.//Bednar J B, ed. Conference on Inverse Scattering: Theory and Application. Philadelphia: Society for Industrial and Applied Mathematics, 1983: 206-220. http://www.oalib.com/references/18989235 |
[13] | Yu J, Hu J, Schuster G T, et al. Prestack migration deconvolution. Geophysics , 2006, 71(2): S53-S62. DOI:10.1190/1.2187783 |
[14] | Lecomte I. Resolution and illumination analyses in PSDM: A ray-based approach. The Leading Edge , 2008, 27(5): 650-663. DOI:10.1190/1.2919584 |
[15] | Tang Y. Target-oriented wave-equation least-squares migration/inversion with phase-encoded Hessian. Geophysics , 2009, 74(6): WCA95-WCA107. DOI:10.1190/1.3204768 |
[16] | Bagaini C, Spagnolini U. 2-D continuation operators and their applications. Geophysics , 1996, 61(6): 1846-1858. DOI:10.1190/1.1444100 |
[17] | Formel S. Theory of differential offset continuation. Geophysics , 2003, 68(2): 718-732. DOI:10.1190/1.1567242 |
[18] | Stolt R H, Benson A K. Seismic Migration Theory and Practice: Theory and Practice. London: Geophysical Press, 1986 . |
[19] | 冯波, 王华忠. 三维偏移距平面波有限差分偏移叠前时间偏移. 地球物理学报 , 2011, 54(11): 2916–2925. Feng B, Wang H Z. 3D offset plane-wave finite-difference pre-stack time migration. Chinese J. Geophys. (in Chinese) , 2011, 54(11): 2916-2925. |
[20] | Sava P C, Fomel S. Angle-domain common-image gathers by wavefield continuation methods. Geophysics , 2003, 68(3): 1065-1074. DOI:10.1190/1.1581078 |
[21] | 刘守伟, 王华忠, 程玖兵, 等. 时空移动成像条件及偏移速度分析. 地球物理学报 , 2008, 51(6): 1883–1891. Liu S W, Wang H Z, Cheng J B, et al. Space-time-shift imaging condition and migration velocity analysis. Chinese J. Geophys. (in Chinese) , 2008, 51(6): 1883-1891. |