一种雷米兹与拉格朗日耦合的有限差分系数优化方法

彭炜颋,黄建平,2

(1.中国石油大学(华东),山东青岛266580;2.海洋国家试验室海洋矿产资源评价与探测技术功能实验室,山东青岛266071)

有限差分法[1-2]是一种经典的数值模拟方法,能够兼顾计算效率与模拟精度,目前已经被广泛应用于勘探地震波场的数值模拟中。有限差分方法是用差分算子逼近空间偏导和时间偏导的方法,该方法可以实现离散模型的数值模拟,但是如果在对模型离散的过程中使用粗网格,模拟结果就会产生频散误差[3-4]。在有限差分方法中,常规的有限差分系数由泰勒级数展开推导而得[5],该差分系数在低波数段能有效且准确地模拟地震波场,但是在高波数段,模拟结果会出现严重的频散误差[6]。随着地震技术的不断发展,对地震波场正演模拟精度的要求不断提高,优化有限差分系数可以在不增加计算量,不改变差分格式的前提下,提高模拟的精度和计算效率。所以如何通过优化差分系数来减小频散误差是有限差分方法研究领域一个重要的研究方向。

在优化差分系数的研究中,用优化算法最小化时空域频散关系或者波数域频散关系的误差来获得优化差分系数是一种常用策略。雍鹏等[7]和邹强等[8]在时空域中,通过最小化给定波数范围内的频散误差来计算优化差分算子。雍鹏等[9]在时空域中实现了时空差分算子的同步优化,并且采用共轭梯度法增加求解过程的稳定性。LIU[10-11]提出用最小二乘法来计算优化差分系数,该算法的目标函数是波数域频散关系误差或时空域频散关系误差。该方法通过求解一定波数范围内的非迭代解,实现了对频散误差的全局优化。最小二乘方法虽然在低波数段产生了频散误差,却有效拓宽了有限差分算子的有效带宽,使得数值模拟中的高波数段频散显著降低。MIAO等[12]用交替方向乘子法解决了基于一范数目标函数的差分系数优化问题,且通过数值实验和理论分析验证了相较于二范数和无穷范数函数,一范数目标函数会获得更小的低波数段频散误差。ZHANG等[13-14]采用模拟退火算法最小化波数域频散关系绝对误差来计算优化差分系数。用模拟退火算法得到的优化差分系数比用最小二乘方法得到的优化差分系数具有更宽的有效带宽,但是由于该方法在迭代过程中的不稳定性,很难用该方法计算高阶的优化差分系数。YANG等[15]应用抽样方法来近似频散关系,该方法的优化差分系数在数值模拟中有着良好的表现。该抽样方法类似无迭代的雷米兹交换算法。AN[16]用类雷米兹交换算法来计算真实频散关系的统一近似,但是没有根据误差限来调整有效带宽,导致该方法在数值模拟过程中产生较大的数值误差。ERIK等[17]提出了利用振幅等纹波特性和最小二乘法利用计算任意样本位置的任意阶导数的优化差分系数方法。该方法基于复值雷米兹交换算法,将复值雷米兹交换算法应用于3个目标函数(总误差、相对误差和群速度误差)来获得优化差分系数。YANG等[18]利用雷米兹交换算法约束交错网格中的频散关系,得到了一种全新的差分系数,该差分系数能显著拓宽模拟带宽。但是,由于YANG等[18]对零波数处的严格约束,导致频散曲线不具有振幅误差等波纹的特性。

HE等[19]通过修改近零波数条件,用雷米兹交换算法约束波数域频散关系的绝对误差。该优化差分系数的频散曲线具有振幅误差等波纹的特性,使得该优化差分系数在数值模拟过程中具有最宽的有效带宽和最大的低波数段误差。为了兼顾雷米兹交换算法较宽的带宽特性并且减少低波数段误差,本文将拉格朗日乘数法[20-21]引入到差分系数优化过程中,用拉格朗日乘数法求解优化问题,得到新的优化差分系数,其中优化问题由二范数目标函数和由雷米兹交换算法导出的约束条件组成;再通过频散测试分析该方法保留较宽的带宽特性和降低低波数段频散误差的效果;最后用均匀模型和改进的Marmousi模型的模型数据验证本文方法的有效性。

1 雷米兹交换法

首先考虑波场空间二阶导数的差分形式,即:

(1)

式中:u为标量波场;h代表空间网格大小;am是有限差分方法的差分系数;M代表差分阶数的一半。

再根据平面波理论,可将波场表示为:

(2)

式中:k为波数;τ为时间采样间隔;ω为角频率。

将公式(2)代入公式(1),并用欧拉公式化简,可得到空间二阶导数波数域频散关系的绝对误差近似式:

(3)

式中:β=kh。在奈奎斯特采样定理下,β∈[0,π]。公式(3)代表了频散误差与波数之间的关系,可以由此计算差分系数的波数域频散关系曲线。

应用雷米兹交换算法之前,需要确定一个准确的近零波数条件。该近零波数条件[19]为:

(4)

雷米兹交换法是一种寻找连续函数的最佳逼近多项式的方法,而最佳逼近多项式和连续函数之间的差值会在n+1个点上正负交替变化且绝对值相等[19],这种关系可以表示为:

(f-p)x1=-(f-p)x2=…=

(-1)n+1(f-p)xn+1

(5)

根据公式(3)可得到波数域频散关系的绝对误差公式,即:

(6)

结合公式(5)和公式(6),雷米兹交换算法计算优化差分系数的线性方程组[19]为:

E(βi)=(-1)iAi=1,2,…,M+1

(7)

式中:A是变量;βi(i=1,2,…,M+1)是在[0,khmax_R]中的采样点,khmax_R是有效波数的最大值。

整理公式(7),可得计算优化差分系数的线性方程组:

(8)

采用雷米兹交换算法求解方程(8)的过程分为以下两步。

1) 将采样点βi代入方程(8)中,并用高斯消元法求解差分系数。求解差分系数am时,同时求出A。当求得的A比误差限小时,khmax_R增加(误差限是人为设置的,即在有效带宽内频散误差曲线的最大值,khmax_R的值就是有效带宽的大小);当求得的A比误差限大时,khmax_R减小。

2) 将步骤1)得到的差分系数代入公式(6),计算频散关系的绝对误差曲线。在频散关系的绝对误差曲线上取每个区间内(由零点划分的区间)最大绝对值点的波数值作为下一次迭代的采样点。在频散关系的绝对误差函数(公式6)中,方程E(β)=0有M+1个零点(包括β=0)。这些零点将绝对误差曲线分成M个区间,每个区间内最大绝对值点的波数值作为下一次迭代的采样点,但是这样只有M个新采样点,而下一次迭代需要M+1个新采样点,为了满足下一次迭代的条件,将第M+1个采样点设置为β=khmax_R。在第一次迭代中,采样点通过在波数范围[0,khmax_R]上均匀采样得到。

重复步骤1)和步骤2)直到A和误差限之间的差值小于一个足够小的值。当循环终止时,通过对优化差分系数计算绝对误差曲线(E(β))得到在该曲线上的M+1个零点。

2 拉格朗日乘数法

空间二阶导数波数域频散关系的绝对误差可从公式(3)变换为如下形式:

(9)

将公式(9)改写为:

(10)

并将公式(10)在波数区间[0,khmax_L]上积分,得到目标函数如下:

(11)

其中,khmax_L与雷米兹交换算法迭代终止时的khmax_R相同。

如前文所述,当迭代循环终止时,曲线E(β)具有M+1个零点,用前M个零点构建约束条件。先将前M个零点代入公式(3)中并化简可得:

(12)

且将公式(12)累加组成约束条件。约束条件形式如下:

(13)

目标函数和约束条件都是基于波数域频散关系的绝对误差,两者必须基于相同的频散误差形式,这样可以保持约束条件和目标函数的一致性。

之后,用拉格朗日乘数法将目标函数和约束条件组合:

ψ(am)=φ(am)+λφ(am)

(14)

其中,λ是拉格朗日算子。将公式(14)对am(m=1,2,…,M)和λ求导,可得如下表达式:

(15)

将公式(15)展开:

(16)

在波数区间[0,khmax_L]上,频散关系绝对误差E的最大值为:

(17)

表1 基于本文方法的二阶导数的优化差分系数(误差限为0.0001)

3 理论误差分析

为对比本文方法的优缺点,定义波数域频散误差为:

(18)

有效带宽是优化的波数范围的大小,本文方法的有效带宽是公式(11)中的khmax_L。有效带宽也和频散曲线上最后一个零点的波数大小成正比,频散曲线上最后一个零点的波数值越大,有效带宽越宽。

如图1所示,在误差限(0.001)相同的条件下,由本文方法的优化差分系数有效带宽小于由雷米兹交换的优化差分系数有效带宽,但是明显宽于最小二乘法的优化差分系数有效带宽。图2是在相同带宽下本文方法与最小二乘方法在不同差分阶数时的对比结果,可以看出,在相同的有效带宽下,本文方法在低波数段的频散误差比最小二乘法更大。

图1 3种优化方法16阶(a)和12阶(b)差分系数的波数域频散曲线对比结果

图2 在相同带宽下本文方法和最小二乘法的频散曲线对比结果a 24阶优化差分系数频散曲线; b 20阶优化差分系数频散曲线; c 16阶优化差分系数频散曲线; d 12阶优化差分系数频散曲线

图3为在相同误差限下,模拟退火法、雷米兹交换算法和本文方法在12阶(M=6)和16阶(M=8)时的频散曲线对比结果。当M=6和M=8时,黑线(本文方法)和红线(模拟退火方法[13])最后一个零点的波数值大小几乎一致,代表这两种方法得到的差分系数有效带宽近乎一致,蓝线(雷米兹交换算法[19])的最后一个零点的波数值大于其它两种方法,这是因为雷米兹交换算法的误差振幅符合等波纹条件,根据切比雪夫准则,该算法具有最大的有效带宽。模拟退火法的优化差分系数具有较宽的有效带宽[13],在该试验中(如图3),两种方法有效带宽几乎相同的现象表明了本文方法具有宽带宽的特性。模拟退火算法由于其不稳定性导致其无法求解高阶差分系数,而本文方法对于高阶差分系数的稳定求解使其在发展和应用前景方面优于模拟退火法。

图3 相同误差限下模拟退火法、雷米兹交换算法和本文方法的频散曲线对比结果a 12阶优化差分系数频散曲线; b 16阶优化差分系数频散曲线

4 数值模拟

4.1 均匀模型

如图4所示,在均匀模型数值模拟试验中,子波的主要能量集中分布在某一波数段内,而子波能量的分布又受主频和网格间距的影响。一般而言,网格间距越小,子波能量越趋向于低波数段。另一方面,通过优化方法得到的优化差分系数在低波数段都会产生频散,而用泰勒展开方法得到的常规差分系数在低波数段的频散趋近于0。所以,在均匀介质模拟且网格间距很小时,常规差分系数的模拟结果优于优化差分系数的模拟结果。常规差分系数的计算公式为:

(19)

为了检验本文方法,本文采用不同网格间距下的子波能量分布图和频散曲线对比图对均匀模型模拟结果进行预测分析。图4a展示了3种优化方法在误差限为0.00001时的8阶优化差分系数频散曲线。图4b是不同网格间距下雷克子波能量归一化结果。由图4b可以看出,当网格间距为5m时,雷克子波的能量主要集中于kh=0.4波数段,并且在该波数段内,最小二乘法的频散误差比其它算法的频散误差更小,所以最小二乘法模拟结果的误差会小于其它两种方法的误差。对于雷米兹交换法来说,在kh=0.4波数段上,雷米兹交换法比本文方法的频散误差更大,因此,当网格间距为5m时,雷米兹交换算法模拟结果的误差大于本文方法的误差。当网格间距为10m时,雷克子波的能量主要集中于kh=0.9波数段,根据图4a所示,当kh=0.9时,3种方法模拟结果的误差比较接近。当网格间距为15m时,雷克子波的能量主要集中于kh=1.2波数段,根据图4a所示,在该波数段上,最小二乘法的频散误差也最大,最小二乘法模拟结果的误差也最大。在图4a中,雷米兹交换法的有效带宽略宽于本文方法的有效带宽,都在1.2左右,且两种方法在kh=1.2附近的频散曲线相似,所以这两种方法模拟结果的误差非常相近。

图4 3种方法的8阶优化差分系数频散曲线(a)以及不同网格间距下的雷克子波能量归一化结果(b)

为了验证前文的分析结果,本文在200×200的网格上定义一个二维均匀模型,网格间距分别为5,10,15m。P波速度为1500m/s。震源位于模型中心。雷克子波的主频为20Hz。图5a是均匀模型,网格间距为5m,传播时间为0.4s,时间步长为0.0005s的波场快照对比结果。图5b和图5c是均匀模型,网格间距分别为10m和15m,传播时间为2.0s,时间步长为0.0005s的波场快照对比结果。图5d是网格间距为10m时的波场快照残差对比结果(波场快照残差等于参考波场快照减去对应方法的波场快照)。均匀模型测试中将30阶常规差分算子的模拟结果作为参考波场。

图5 4种方法的波场快照和波场残差对比结果a 网格间距为5m时的波场快照对比结果; b 网格间距为10m时的波场快照对比结果; c 网格间距为15m时的波场快照对比结果; d 网格间距为10m时的波场残差对比结果

为了直观地对比不同优化方法的频散压制能力,本文计算了这些波场快照频散误差总和与波场快照总误差(图6)。频散误差总和是波场快照残差沿着深度方向的总和。总误差就是将频散误差总和相加。图6a、图6b、图6c分别是图5a、图5b、图5c中不同方法的波场快照频散误差总和对比结果。图6d是图6a、图6b、图6c的总误差对比结果。由图5和图6可知,均匀模型数值模拟结果基本符合本文图4的预测结果。

在网格间距为5m时(图6a),由于子波主要能量集中于低波数段,所以模拟结果误差与低波数段频散呈正相关,相较于雷米兹交换算法,本文方法的模拟结果误差较小,说明本文方法具有比雷米兹交换算法更小的低波数段误差。图6c中,由于子波能量主要集中于高波数段,所以模拟结果误差与带宽大小呈正相关,在此前提下,本文方法和雷米兹交换算法具有相近的模拟结果误差,说明本文方法与雷米兹交换算法具有相似的有效带宽。均匀模型的测试结果验证了本文方法相较于雷米兹交换算法具有较低的低波数段频散和相似的有效带宽。

图6 不同网格间距的频散误差总和和总误差对比结果a 网格间距为5m时的频散误差总和对比结果; b 网格间距为10m时的频散误差总和对比结果; c 网格间距为15m时的频散误差总和对比结果; d 不同网格间距下的总误差对比结果

4.2 改进的Marmousi模型

为了在炮记录中更加直接和准确地观察直达波的频散,本文在标准的Marmousi模型上加了一层速度为1500m/s的地层。整个模型被离散为1360×700个网格点。改进的Marmousi模型如图7所示。图8、图9和图10分别为网格间距为10,20,30m时,不同传播时间的波场快照残差以及炮记录快照残差对比结果。数值模拟试验采用的时间步长为0.0005s,主频为20Hz,炮点位置为模型表面的中点。数值模拟采用吸收边界条件处理边界反射。图8、图9和图10 展示的炮记录残差由参考炮记录减去优化差分算子的模拟炮记录得到,波场快照残差由参考波场快照减去优化差分算子的波场快照所得到。其中,优化差分算子均是在误差限为0.0001时计算得到的20阶优化差分算子。参考波场快照和参考炮记录是40阶常规差分系数的模拟结果。

图7 改进的Marmousi速度模型

图8 网格间距10m时波场快照残差(a)和炮记录快照残差(b)对比结果

图9 网格间距20m时波场快照残差(a)和炮记录快照残差(b)对比结果

图10 网格间距30m时波场快照残差(a)和炮记录快照残差(b)对比结果

当网格间距为10m时,本文方法的波场快照残差和炮记录残差小于雷米兹交换法和最小二乘法的波场快照残差和炮记录残差,说明当网格间距为10m时,本文方法的优化效果优于其它两种方法的优化结果(图8)。但当网格间距继续增大至20m时,本文方法的波场快照残差和炮记录残差与雷米兹交换法的波场快照残差和炮记录残差相似,但是这两种方法的波场快照残差和炮记录残差小于最小二乘法方法的结果,说明当网格间距为20m时,本文方法的优化效果与雷米兹交换算法的优化效果相似,且两种方法的优化效果优于最小二乘法的优化效果(图9)。当网格间距为30m时,3种方法的波场快照残差和炮记录残差相似,说明网格间距为30m时3种优化差分算子的优化效果趋近于一致(图10)。在改进的Marmousi模型的模拟过程中,本文方法在网格间距较小的情况下,对于复杂模型中频散误差的压制效果优于最小二乘法和雷米兹交换算法的压制效果。

5 结论

为了兼容雷米兹算法的优点并减少该算法在低波数段产生的数值频散,本文提出了一种拉格朗日乘数法与雷米兹交换算法耦合的差分系数优化方法,引入雷米兹交换算法和拉格朗日乘数法来优化显式有限差分系数。频散分析和数值模拟试验结果表明,本文方法具有较宽的有效带宽和更低的低波数段数值频散,这些特性拓宽了雷米兹交换算法的适用范围。