一种基于新差分模板和无穷范数的震源波场重建方法

包乾宗 戴 雪 梁 雪

(①长安大学地质工程与测绘学院地球物理系，陕西西安 710054；②海洋油气勘探国家工程研究中心，北京 100028；③东方地球物理公司，河北涿州 072750)

0 引言

伴随状态法[1]通过两次正演运算获取模型参数的梯度，可避免对所有震源和接收点格林函数的直接计算和存储，已广泛应用于地震波动方程偏移和反演，如逆时偏移(RTM)[2-5]、全波形反演(FWI)[6-10]、最小二乘逆时偏移(LSRTM)[11-13]、波动方程旅行时反演(WETI)[14-16]等。

伴随状态法通过震源波场和伴随波场间的相互运算(如相关)计算参数梯度/敏感核。但震源波场和伴随波场分别沿时间正向和逆向传播。为同时访问震源波场和伴随波场，震源波场需先保存在计算机内存中或在波场反向传播中进行重建。存储所有时刻、所有空间网格点的震源波场是不现实的，特别是三维情况。最优检测点方法(Optimal CheckPointing)在波场正向传播中仅需保存几个时刻(检测点)的震源波场，在波场反向传播中将存储的检测点波场作为初值重新计算相邻检测点之间的震源波场[17-19]。最优检测点方法可大幅减少计算机内存消耗，但需额外的正演运算。边界值方法采用位于边界区域几层(与差分精度有关)网格点的波场和最后两个时刻的快照重建震源波场[20-22]。常规边界值方法需要存储每个时刻M层网格点的边界波场(M为有限差分算子长度)，存储量仍然较大。Feng等[23]采用变阶数有限差分算子重建震源波场，只需存储一层空间网格点的边界波场。但由于降阶处理，该方法的重建精度不高。Tan等[24]提出基于拉格朗日多项式或高阶波动方程外推法，该方法需保存一层或两层网格点的边界波场，内存需求约为常规方法的 37.5%。Liu等[25]采用边界波场的线性组合重建震源波场，所需存储量仅为常规方法的1/M。段沛然等[26]将该方法推广到交错网格有限差分波场重建。Liu等[25]的方法存在如下问题：基于L2范数准则求取差分系数，有效波数范围/带宽较窄；重建系数的自由度较小(最小为2)，很难精确逼近重建频散关系。为保证重建误差小于 0.0025，每个波长需要的采样点数至少为 5。Ren等[27]设计了新的差分模板，采用边界区域的N层波场和M-N层波场的线性组合来重建内部区域的震源波场(0 ≤N≤M)。该方法通过调整N的大小实现重建精度和内存需求的折衷。当N=0时，退化为Liu等[25]的方法。Ren等[27]的方法每个波长只需3个采样点就可保证重建误差小于0.001。

Ren等[27]的方法是基于规则网格设计的，无法直接用于变密度声波或弹性波速度—应力方程RTM或FWI。本文发展了基于交错网格有限差分的震源波场重建方法，设计了新的交错网格差分模板，仍通过存储边界区域N层波场和M-N层波场的线性组合来重建震源波场；构建了无穷范数极小化目标函数，采用Remez交换算法[28-29]优化了重建系数。基于改进的差分模板和优化的重建系数，新方法可大幅提高震源波场的重建精度，且内存需求仅为常规方法的(N+1)/M。本文通过精度和稳定性分析、声波RTM和弹性波FWI算例验证了方法的有效性。

1 原理方法

1.1 交错网格有限差分震源波场重建方法

交错网格有限差分离散格式为[30]

(1)

式中：ai为差分系数；h为空间网格间隔；p为地震波场(压强或速度)。常规边界值方法采用边界区域的边界波场p-i+0.5(i=1,2,…,M)重建内部区域的震源波场[20-22]。为减小内存需求，设计如图1所示的差分模板。内部区域各层(x=jh，j=0,1,…,M-1)的空间导数通过下式计算

图1 交错网络差分模板示意图

(2)

式中：ai(i=1,2,…,M)和bj,i(j=1,2,…,M-1,i=1,2,…,N+j+1)统称为重建系数。式(2)可改写为

(3)

波场正向延拓中存储边界波场p-i+0.5(i=1,2,…,N)和波场的线性组合A，波场反向延拓中基于式(3)重建内部区域不同层(x=jh，j=0,1,…,M-1)的空间偏导数。所提方法只需保存N层边界波场和波场线性组合A，内存需求为常规边界值方法的(N+1)/M。该方法可通过调节N控制重建精度和存储需求。随着N增大，重建系数的个数增加(自由度增大)，精度提高，但存储量随之增加。当N=M或0时，本文方法退化为常规边界值方法[20-22]或Liu等[25]的方法。

基于平面波理论推导频散关系，将波场表示为

pi=p0eikxh

(4)

式中kx为x方向的波数。将式(4)代入式(3)，有

(5)

式中j=1,2,…,M-1。式(5)为新重建方法的频散关系，定义如下的相对误差用于后续重建精度分析

(6)

(7)

1.2 重建系数优化

(8)

式(8)为约束方程，可保证低波数成分的重建精度。

将式(8)代入式(6)，可得

δ0(kxh)=φ1(kxh)+

(9)

δj(kxh)=φ1(kxh)+

ψ1(kxh)]-1

(10)

建立基于无穷范数的极小化目标函数

(11)

式中βj为最大波数(j=0,1,…,M-1)。利用式(11)无法直接求取导数，需采用Remez交换算法求解，详细步骤如下。

(1)求解线性方程组

δ0(kxlh)=(-1)lλ0l=1,2,…,M

(12)

得到ai(i=2,3,…,M)和λ0。式中：λ0为待求常数；kxl为kx的第l个采样点。M个初始采样点通过对波数范围[0,β0]均匀采样得到。

(2)采用式(11)计算E0。当E0>η(η为最大允许误差)，执行步骤(3)；当E0≤η，执行步骤(6)。

(3)式(12)中重建误差在采样点处正负交替，M个采样点之间存在M-1个根，M-1个根将[0,β0]分成M个区间，搜索得到M个区间的极值点。

(4)将M个极值点作为新的采样点，重复步骤(1)～步骤(3)，直到满足收敛条件E0≤η，输出优化的ai(i=2,3,…,M)和β0。

(5)若达到最大迭代次数，减小β0，重复步骤(1)～步骤(4)。

(6)基于约束方程式(8)得到优化的a1。

(7)求解下式得到bj,i(i=2,3,…,N+j+1)和λj(j=1,2,…,M-1)。

δj(kxlh)=(-1)lλjj=1,2,…,M-1；

l=1,2,…,N+j+1

(13)

式中λj为待求常数。初始采样点(N+j+1个)可通过对波数范围[0,βj]均匀采样得到。

(8)采用式(11)计算Ej(j=1,2,…,M-1)。当Ej>η，执行步骤(9)；当Ej≤η，执行步骤(12)。

(9)式(13)中重建误差在采样点处正负交替，N+j+1个采样点之间存在N+j个根，N+j个根将[0,βj]分成N+j+1个区间，搜索得到N+j+1个区间的极值点。

(10)将N+j+1个极值点作为新的采样点重复步骤(7)～步骤(9)，直到满足收敛条件Ej≤η，输出优化的bj,i和βj。

(11)若达到最大迭代次数，减小βj，重复步骤(7)～步骤(10)。

(12)基于约束方程式(8)得到优化的bj,1。

对于不同的层x=jh(j=0,1,…,M-1)，优化的最大波数βj是不同的。用于评价所提方法的重建精度的有效波数范围/带宽为

(14)

上述优化步骤可保证所提方法的重建误差在[0,βf]范围内均小于最大允许误差η。表1和表2分别给出N=0和1时的重建系数。

表1 本文方法的重建系数(N=0，M=6和η=0.001)

表2 本文方法的重建系数(N=1，M=6和η=0.001)

2 精度分析

采用式(6)分析重建方法的精度。图2和图3分别为N=0和1时的频散曲线，可见本文方法的频散曲线呈波纹状分布。与δ0和δj(j=2,3,…,M-1)相比，δ1(x=h层的精度)的幅值最先超出最大允许误差(η= 0.001)。因此，δ1直接决定重建方法的精度。当M增大时，x=0和x=jh(j=2,3,…,M-1)层的精度增加，但x=h层的精度几乎不变；当N增大时，x=h层的精度大幅改善。为了更好地比较，表3给出不同方法的有效带宽βf。常规方法中采用优化差分系数[29]。随着N的增大，有效波数范围变大；当N=1时，本文方法的带宽与常规方法(存储M层)非常接近(约为2)。当一个波长采两个样点时，β=kxh达到尼奎斯特极限(等于π)。通过G=2π/βf将带宽βf转换为一个波长需要的采样点数G，表4给出不同方法的G值。由表可知：本文方法的G值略大于常规方法。随着N增大，G逐渐减小，当N=1时，一个波长只需要3个采样点就可使重建误差小于0.001。

图2 本文方法不同有限差分算子的频散曲线(N=0，η=0.001)

图3 本文方法不同有限差分算子的频散曲线(N=1和η=0.001)

表3 不同方法的有效带宽βf(η=0.001)

表4 不同方法一个波长需要的采样点数G(η=0.001)

3 稳定性分析

基于冯·诺依曼分析[27]可推导本文方法的稳定性条件。二维稳定性因子为

(15)

式中：j=1,2,…,M—1，s为允许的最大库朗数，稳定性随着s增大逐渐改善。表3给出不同方法的稳定性因子。与常规方法相比，本文方法的稳定性条件更加严格。随着N的增加，稳定性略微变差。

表5 不同方法的稳定性因子s

4 模型算例

将本文方法应用于声波逆时偏移和弹性波全波形反演进行检验。时间二阶、空间十二阶优化交错网格有限差分[29]用于波场正、反向延拓。常规边界值重建方法采用优化差分系数[29]。正向震源波场、常规方法得到的像和反演结果作为参考解，计算不同方法结果的最大绝对误差εMAV和均方根误差εRMS。

4.1 声波逆时偏移

声波Marmousi模型如图4,时间步长为1.0ms，记录时间为4.0s，空间网格大小为10m×10m，网格点数为501×353，震源为20Hz主频的Ricker子波。51炮和501个检波点均匀分布于地表。图5为不同方法重建的波场，炮点位于(9.5km,0)。表6给出不同方法重建波场的εMAV和εRMS。由表可知，常规方法的重建误差非常小，本文方法的εMAV和εRMS分别在10-2和 10-4数量级。重建波场的误差随着N的增大有所减小。图6为逆时偏移中采用不同重建方法得到的像，表7为本文方法的εMAV和εRMS，可见，N=0或1都能保证成像结果的均方根误差小于10-3。当N=0或1时，本文方法的存储量为常规方法的16.6%或33.3%。为兼顾重建精度和内存需求，声波RTM建议采用N=0。

图6 声波Marmousi模型不同方法的逆时偏移结果

表6 声波Marmousi模型不同方法重建波场的ε MAV和ε RMS

表7 声波Marmousi模型本文方法不同N的ε MAV和ε RMS

图4 声波Marmousi模型

图5 声波Marmousi模型不同方法重建的波场

4.2 弹性波全波形反演

弹性波Sigsbee2A模型如图7所示。空间网格尺寸为16m×16m，网格点数为351×186，时间步长为1.5ms，记录长度为4.5s。爆炸震源激发主频为15Hz的Ricker子波。36个震源和351个检波点均匀分布于地表。初始纵波和横波速度模型为随深度线性变化的一维模型(图8)。采用多尺度反演(0～5Hz、0～10Hz和0～15Hz)策略，每个尺度上迭代10次。图9给出不同方法得到的波场，炮点位于(2720m，0)，表8给出不同方法重建波场的εMAV和εRMS。常规方法(存储M层)的误差在10-6数量级，本文方法的误差随N的增大而减小。图10为采用不同重建方法进行全波形反演得到的纵、横波速度模型。表9为本文方法反演结果的εMAV和εRMS。当N=0时，反演结果的精度较差；当N=1时，反演精度得到大幅提高(εRMS<0.5 %)。本文方法的存储需求仅为常规方法的(N+1)/M。波场重建误差对模型参数梯度的精度有影响，为保证反演精度，全波形反演建议采用N=1。

图7 弹性Sigsbee2A模型

图8 弹性Sigsbee2A模型的初始模型

图9 弹性Sigsbee2A模型不同方法重建的1.5s波场

图10 弹性Sigsbee2A模型不同方法重建波场的全波形反演结果

表 8 弹性Sigsbee2A模型不同方法重建波场的ε MAV和ε RMS

表9 弹性Sigsbee2A模型本文方法不同N时反演结果的ε MAV和ε RMS

与常规方法相比，本文方法可降低存储量，但重建波场的误差有所增加(详见精度分析和模型算例)。实际应用中，可通过调整N的大小进行重建精度和内存需求的折衷。由式(1)和式(3)可知，常规方法的计算复杂度为3M2，本文方法计算复杂度为(7M2+2MN+3M-4)/2。当N为0或1时，本文方法的计算量约为常规方法的7/6倍。与波场正、反向延拓相比，本文震源波场重建方法增加的计算量微不足道。

5 结论

本文提出了一种交错网格有限差分震源波场重建方法。该方法通过存储N层边界波场和M—N层波场的线性组合重建震源波场。为提高重建精度，基于无穷范数准则和Remez交换算法优化了重建系数。详细分析了提出方法的精度和稳定性，并将其应用于声波逆时偏移和弹性波全波形反演。新方法可通过调整N的大小实现重建精度和存储需求的平衡。新方法能够获得较好的重建波场、像和反演结果，且内存需求仅为常规方法的(N+1)/M。

本文方法可直接应用于三维逆时偏移和全波形反演或其他基于伴随状态法的地球物理问题中，如最小二乘逆时偏移、波动方程旅行时反演或全球尺度波形层析等。但新方法无法用于基于黏声或黏弹波动方程的成像和反演中。如何避免衰减介质波场重建的不稳定性需要进一步考虑。