基于最小范数优化交错网格有限差分系数的波动方程数值模拟

2021-10-23文晓涛王文化

石油地球物理勘探 2021年5期

唐超文晓涛*② 王文化

(①成都理工大学地球物理学院，四川成都 610059；②成都理工大学地球勘探与信息技术教育部重点实验室，四川成都 610059)

0 引言

波动方程正演模拟是研究全波形反演及逆时偏移技术的基础。目前常用的波动方程数值模拟方法主要有伪谱法[1-2]、有限元法[3-4]、有限差分法[5-9]等。其中有限差分法具有内存占用小、易于实现和计算效率高的优点，因此它在波动方程数值模拟领域中广受欢迎。但是由于差分网格离散引起的数值频散无法避免，直接影响数值模拟的精度。

时域有限差分方法求解波动方程在时间域和空间域同时进行，按照数值相速度和真实传播速度的相对大小，其固有的频散误差可分为时间频散和空间频散[10]。Koene等[11]提出的时间频散变换，在远远超过稳定性极限所允许的时间步长条件下，能很好地去除时间频散。因此，本文主要针对空间频散问题开展研究。通常用于压制空间数值频散的方法有以下几种：其一，采用较小的空间步长[12]，但这会大大增加计算量和所需计算内存；其二，空间方向采用高阶差分格式[13]；其三，利用FCT(Flux-Corrected Transport)技术[14-15]削弱数值频散影响；其四，利用优化算法求取新的差分系数[16-20]进行波动方程数值模拟。在相同阶数下，相对于传统泰勒级数展开(TE)方法，利用优化算法求取新的差分系数在较大波数域区间精度更高，因此能取得更好的数值频散压制效果。Holberg[16]首次提出通过最小化给定空间频带内群速度的峰值相对误差优化空间有限差分(FD)算子，但是对群速度的最大相对误差没有严格的限制，导致模拟结果与精确解有较大偏差。Liu[17]基于L2范数建立目标函数，采用最小二乘法(LS)优化近似求取差分系数。此外，Zhang等[18]、Yang等[19]、He等[20]基于L∞范数建立目标函数求解有限差分系数，用于波动方程数值模拟。Zhang等[18]用模拟退火算法求得了声波方程空间导数的有限差分系数，并且给出了0.0001的误差容限。Yang等[19]采用Remez迭代算法求解空间一阶偏导数交错网格差分系数。He等[20]在Yang等[19]的基础上引入新的约束条件，得到“等波纹”解，获得理论上最宽有效频带。

现有的L2范数和L∞范数方法只考虑峰值误差，缺乏对中低波数的误差约束，因此实际数值测试中的误差积累较大。虽然上述方法都能有效减少数值频散，但是应该更加关注有限差分方法在波场传播中的误差积累，而不是一味地追求给定误差容限下的最大覆盖带宽[20-21]。Miao等[21]基于L1范数建立新的目标函数，然后通过交替方向乘子法算法(ADMM)[22]获得了空间二阶偏导数的规则网格有限差分系数。数值实验表明，在不改变误差容限的情况下，此方法尽管略微降低了有效带宽的覆盖范围，但是能加强对中低波数的误差约束，减小有效带宽内的总误差。

本文首先将ADMM优化方法推广到空间一阶偏导数交错网格有限差分系数的求解，将其用于一阶应力—速度声波方程的交错网格有限差分数值模拟。然后对基于三种不同范数优化方法进行频散曲线误差分析。最后利用均匀介质模型和复杂模型的长时程数值实验，证明本文L1范数优化方法在减小误差积累方面的优越性。

1 方法原理

1.1 交错网格有限差分格式

非均匀介质的二维一阶应力—速度声波方程组可表示为

(1)

式中：P为应力(标量)；Vx和Vz分别表示x和z方向上质点振动速度分量；ρ为介质密度；v为地震波速度。

交错网格有限差分格式在同等条件下比规则网格有限差分格式精度更高。通常，使用高阶有限差分系数可以提高空间导数的精度[7]，空间一阶导数的2M阶有限差分交错网格定义为

(2)

式中：cm是有限差分系数；h是空间网格间距。根据平面波理论有

u(x)=u0ejkx

(3)

式中：u0是常数；j为虚数单位；k为波数。令β=kh/2，代入式(2)得到

(4)

本文的目的是通过在给定范围内最小化绝对误差的总和计算优化的差分系数。绝对误差的总和可以表示为

(5)

将上式中β离散为β(i)=ikmaxh/N,其中N为离散度，kmax为最大波数。在忽略常数项kmaxh/N后，目标函数E可写成

(6)

1.2 目标函数的重构和正则化

将式(6)改写成线性方程组的形式

(7)

(8)

Wang等[23]证明了直接优化目标函数F(c)的常系数矩阵c可能导致数值振荡。因此，使用正则化技术[24]克服不稳定解。建立正则化目标函数为

ψ(c)=F(c)+J(c)

(9)

(10)

式中：α=0.0001为正则化参数；D为单位矩阵。本文的目的是通过求解正则化的问题来找到最优解，即

(11)

1.3 求解正则化最小模型

式(11)为非线性问题，目前可以运用多种方法求解。但是为了高效和精确地求得常系数，采用交替方向乘子法[22]将原始问题转换为约束问题

(12)

式中d=Ac+b，对应的增广拉格朗日函数可以写成

(13)

式中：u为拉格朗日变量；η为软阈值。后面的模型测试中取η=40。

ADMM求解最小化模型主要分为三步：更新c、更新d、更新u。

第一步：更新c

(14)

(15)

第二步：更新d

(16)

上式为一个套索问题[26]，可以用次微分学的封闭形式解决。上式的解可以表达成

d(k+1)=S1/η(Ac(k+1)+b+u(k))

(17)

式中软阈值算子S定义为

(18)

第三步：更新u

u(k+1)=u(k)+(Ac(k+1)+b-d(k+1))

(19)

1.4 搜索有限频带范围的最大值kmax

为了得到带宽的上限kmax，本文计算不同kmax和M条件下的最大误差

(20)

(21)

为了寻找最优解，从一个小的kmax开始，它产生一个最大误差εmax。当εmax小于给定的容错允差时，稍微增大kmax并重复上述过程，直到εmax等于(或小于)给定的容错允差。表1为M=8时不同范数优化后求解的差分系数。L1范数优化后求解的差分系数如表2所示。

表1 不同优化方法的差分系数(M=8)

表2 一阶导数交错网格有限差分系数

从图1中不同差分阶数的频散曲线可以看出，采用增大差分阶数可以获得更广的频带覆盖范围，且随着差分阶数的增加，低波数区间的频散误差逐渐减小。为了对各种方法进行比较，分别基于L1范数、L2范数[12]和L∞范数[20]求解声波方程的交错网格有限差分系数，并对其进行频散分析。图2中，在M=8(16阶)时，在误差阀值为0.0001的条件下，相对于L2范数和L∞范数来说，尽管基于L1范数优化方法频带覆盖范围略微变窄，但能更好地约束在中低波数区间的频散误差。

图1 L1范数优化后的频散曲线

图2 三种范数约束下的频散曲线(M=8)

1.5 稳定性分析

二维一阶应力—速度声波方程交错网格有限差

分离散格式的稳定性条件[26]为

r≤s

(22)

式中：r为库郎数；s为稳定性因子[27]

(23)

泰勒展开与L1、L2、L∞约束下的稳定性因子曲线如图3所示。

图3 不同方法的稳定性因子随差分阶数(2M)的变化

总体来看，无论是传统TE方法还是各种范数的优化方法，稳定性因子s都会随着差分阶数的增大而减小。此外，在相同差分阶数时，TE方法稳定性因子的值大于其他优化方法。值得注意的是，本文利用L1范数约束优化求得差分系数的稳定性因子比另外两种范数大，证明其稳定性限制条件比另外两种方法宽松。

2 理论模型试算

在本次正演模拟实验中，分别对均匀介质模型和Marmousi Ⅱ模型进行正演数值模拟。

2.1 均匀介质模型

均匀介质模型中，纵波速度为2000m/s，网格大小为401×401，网格间距h=5m，时间步长为0.2ms，使用主频为30Hz的雷克子波，震源放在模型的中心位置(1005m，1005m)。为获得声波长时长的传播记录，在未加边界条件下，总采样时间设为2s。不同优化方法得到二维交错网格有限差分声波方程数值模拟的波场快照如图4所示；图5为基于三种不同范数优化方法与采用传统泰勒120阶展开算法的模拟结果(参考解)的差值。在接收点Ra(1500m，500m)和Rb(1000m，500m)的单点部分接收记录如图6所示；图7为基于三种不同范数优化方法结果与参考解的差值。

图4 均匀介质模型波场快照(左图：t=0.5s；右图t=2.0s)

图5 均匀介质模型波场快照与参考解的差值(左图：t=0.5s；右图t=2.0s)

图6 均匀介质模型Ra和Rb单点部分时段接收记录

图7 均匀介质模型Ra和Rb单点部分接收时段记录与参考解的差值

从图5中波场快照的差值对比可看出，三种不同范数优化算法的模拟结果中，对误差控制最好的为L1范数，其次是L2范数，最差的是L∞范数，这与频散关系(图2)分析结果是一致的。图5中t=2.0s时刻的差值比t=0.5s时刻更大，说明随着波场传播时间的增加，误差积累更严重。L1范数对误差积累的控制较好，从Ra点和Rb点接收记录分析(图6、图7)也可得出这一结论。

2.2 复杂模型

复杂模型Marmousi Ⅱ速度模型如图8所示，纵波速度范围为1500m/s～4800m/s。震源函数使用主频为25Hz的雷克子波，震源在(1800m，250m)处，接收线与炮点同深度且水平排列。空间步长5m，时间采样间隔0.2ms，总采样时长4s。用不同优化方法得到二维交错网格有限差分声波方程数值模拟的单炮记录，提取第200道和第500道的部分时段地震记录如图9所示，图10为基于三种不同范数优化方法模拟结果与参考解的差值。不同优化方法得到的单炮炮集如图11所示，与参考解的炮集差如图12所示。