基于LSQR算法的二维声波方程频率域正演模拟与数值实现

张鑫磊， 陈建宇

(中国石油大学(北京) 地球科学学院，北京 102200)

0 引言

近几十年人们对地震波数值模拟进行了许多研究，取得了不少成就。大多数的波场模拟都是在时间域进行的，时间域计算方法是按时间片递推计算，由前一个时间片的波场值推出下一个时间片的波场值。因此前一时刻的计算误差会累加到下一时刻，如果计算的时间切片数目较多，计算误差将会累计，导致计算结果信噪比下降、精度太低而无法为接下来的波形反演提供研究基础。近三十年来，学者做了大量工作，频率域正演模拟最早由Lysmer等[1]提出；Shin等[2]进一步发展了这种方法，将频率域有限元法应用于地震波形反演; Jo等[3]提出了最优化9点差分格式，极大地减小了数值频散效应;Stekl等[4]在Jo的基础上，将最优化9点差分格式推广到变密度声波方程和黏弹性介质的弹性波方程; Min等[5]提出了一种频率域标量波动方程25点差分格式，通过计算最优化系数来减小数值频散，减少了空间采样点数;Hustedt[6]研究了交错网格上导数算子的四阶精度有限差分格式，给求解方程组带来更大带宽。在国内，近些年来吴国忱教授[7-8]做了大量研究，将25点优化差分算子理论应用于VIT介质和TTI介质的正演模拟，取得了较好的模拟效果;王华忠教授[9]也在这方面做了大量研究，利用最优化9点差分算子对标量波方程做了数值模拟。

频率域正演模拟涉及到大型稀疏矩阵的求解问题，求解的方法主要分为直接法(高斯消元法和LU分解法等)与迭代法(LSCG最小二乘共轭梯度法和LSQR最小二乘正交分解法等)。直接法通常占用内存大，资源消耗严重，因此笔者采用LSQR(Least Square QR-factorization)最小二乘正交分解法，LSQR算法是Sanders和Paige[10-11]提出的，它是利用Lanczos法求解最小二乘问题的一种方法。刘伊克等[12-13]将其运用到地震层析成像处理中；随后杨文采[14]利用该算法完善了地球物理反演理论；刘劲松等[15]又对该算法进行了并行改进，使其计算效率大大增加。LSQR具有占用储存空间少、计算速度快等优点，可以利用矩阵稀疏性简化运算，适合用来求解大型稀疏矩阵。与其他迭代方法相比，LSQR具有更快的收敛速度及更精确的结果。

笔者采用的差分格式为优化17点差分算子格式[15]，适合高精度正演模拟。该差分格式是在四阶精度9点差分格式的基础上，结合0°系和45°坐标系构造而来，具有计算精度高，时间短，矩阵带宽小等优点。在17点差分格式基础上同时采用稀疏矩阵LSQR算法计算波场，并使用PML完全匹配层吸收边界反射，最终达到提高正演模拟效率的目的。

1 基本原理

1.1 频率域17点有限差分格式

17点差分格式与传统的最优9点差分格式类似，将各项同性介质的频率域波动方程表示为0°和45°坐标系下的加权，用最优化方法计算得到加权系数[16]。

首先构造四阶精度的9点差分格式，时间域二维声波方程对时间t做傅里叶变换可以得到频率域声波方程:

(1)

地震波场上任意位置拉普拉斯算子可以写成加权和的形式：

(2)

对式(1)采用4阶精度有限差分格式可得0°坐标系9点差分格式：

(3)

其中:Δ=Δx=Δz为网格间距，对四阶9点差分格式在45°坐标系下进行扩展。则构造出的17点差分格式可表示为：

图1 17点差分格式频散曲线Fig.1 Dispersion curves of 17-point scheme

F(ω)=β1Um-2,n-2+β2Um-2,n+β1Um-2,n+2+

β3Um-1,n-1+β4Um-1,n+β3Um-1,n+1+

β2Um,n-2+β4Um,n-1+γUm,n+β4Um,n+1+

β2Um,n+2+β3Um-1,n+1+β4Um+1,n+

β3Um+1,n+1+β1Um+2,n-2+β2Um+2,n+

β1Um+2,n+2

(4)

1.2 频散特征

加权系数可以通过频散方程采用最优化方法得到[17-18]，将平面简谐波表达式代入式(4)，可得到17点格式频散方程：

4c-4d-4e)D)]/{(3aΔ2[30-16A+

B]+3Δ2(1-a)[15-16C+D])}

(5)

相速度的定义为：vρh=ω/k，代入式(5)可得到相速度的表达式：

2(1-b-4c-4d-4e)D)}/{(3aΔ2[30-

16A+B]+3Δ2(1-a)[15-16C+D])}

(6)

传播角度的取值范围为[0,π/4]，1/G的最大取值为0.4，为求出系数最优解，令k*=1/G，并使式(6)的值接近“1”，定义相速度的误差为：

ε(θ,k*,a,b,c,d,e)= [1-(vρh(θ,k*,a,

b,c,d,e)/v)]2

(7)

则目标函数为：

c,d,e)/v]2dk*dθ

(8)

依据式(8)可以解出最终的优化系数：a=1.067 3、b=0.887 5、c=0.025 1、d=0.023 7、e=-0.020 4，图2给出了 17点格式的频散曲线，可以看出相速度误差基本可控制在1%以内。

图2 完全匹配层边界示意图Fig.2 PML border diagram

1.3 完全匹配层

为了消除边界反射的影响必须引入人工边界条件，笔者采用频率-空间域的完全匹配层PML边界条件，在最外层网格区域的四个方向上构建PML[19-20]。

设计算区域从x= -d开始，则x=0处为PML层边界，PML层的厚度为0至 -d，设为L，对平面波方程加入衰减系数，则在边界处入射波经过了两次PML层的衰减，可表示为：

U(kx,kz,ω)=U(kx,kz,ω)

(9)

在x=0处的反射波方程为：

Ur(kx,kz,ω)=RUr(kx,kz,ω)

(10)

则在x=0处的反射系数为：

(11)

我们可以通过设置计算区域衰减因子为零值，边界外匹配层衰减因子为非零值而达到消除边界反射的目的[21-22]。对式(11)两边做傅里叶反变换，并求x的二阶导数，可得：

(12)

(13)

(14)

其中：

(15)

(16)

同理可得纵向匹配层关系式为式(17)。

(17)

角点区域与边界区域的处理方式不同，如图2所示，在四个角点区域对波场的衰减是对x方向和z方向衰减的叠加。则PML条件下的声波方程可写为式(18)。

(18)

结合式(4)可得PML条件17点有限差分格式为式(19)

(19)

图3 正演模拟的计算流程Fig.3 Calculation flow of forward modeling(a)直接法正演流程；(b)LSQR算法正演流程

图4 PML效果图Fig.4 PML effect diagram(a)5点网格厚度；(b)10点网格厚度

1.4 LSQR算法求解波场

将波场空间划分为Nx×Nz的网格区域，可以令l=Nx×Nz，将式(17)改写为：

F(ω)=S(v,ω)U(ω)

(20)

其中：S为复阻抗矩阵，大小为l×l；U为频率域的离散波场值，是l×1的列向量；F为震源项，大小也是l×1。首先将复系数矩阵转化为实系数分块矩阵，则式(18)可改写为：

(SRe+iSIm)(URe+iUIm)=FRe+iFIm

(21)

式中：SRe和SIm分别为阻抗矩阵的实部和虚部；URe和UIm分别为波场矩阵的实部和虚部，i为虚数单位。将(19)式展开可以得到一个方程组为式(22)。

(22)

写成分块矩阵为式(23)。

(23)

将复阻抗矩阵S按照式(21)排列，只存储非零元，然后采用LSQR法求解方程，就可以快速高效的求解出各网格点的波场值。

将要求解的式(21)简写为式(24)。

Am×mxm×1=bm×1

(24)

设x0为式(21)要求解的波场x的初始值，β0为式(21)已知震源信息的二范数，则LSQR算法步骤总结如下：

1)初始化。

(25)

2)迭代循环。

图5 PML效果图(18点网格厚度)Fig.5 PML effect diagram (18 point mesh thickness)

图6 均匀层状模型Fig.6 Uniform layer model

(26)

图7 180 ms波场快照Fig.7 Snapshot of 180 ms wave field(a)频率域LSQR法波场；(b)常规时间域波有限差分波场

3)参数修改。

(27)

其中:式(25)中的xi+1为最终输出的波场虚实分解的列向量，将其虚实结合重构即可求出最终波场值，图3为正演模拟的计算流程，左为直接法正演流程，右为LSQR算法正演流程，主要区别在于求解频率域波场的方法上，LSQR算法采用阻抗矩阵虚实分解和稀疏矩阵存储方法，具有更快的计算速度和更少的资源占用。

2 模型试算

2.1 PML厚度取值

PML厚度不同取值效果见图4、图5，均为均质模型180 ms波场快照，由图4、图5可以看出当取值达到18点网格厚度时已经有较好的吸收效果，所以最终取值为18点网格厚度。

2.2 均匀层状模型

试算模型为100×100网格的均匀层状模型(图6)，上层速度为1 500 m/s，下层速度3 000 m/s，网格间距Δ=Δx=Δz=10 m，震源采用Rick子波，主频30 Hz，位于(21，51)点，采样间隔1 ms，记录长度1 024 ms，检波器接收位置位于x=20，每隔一个网格点放置一个，边界采用18网格厚度的PML边界，求解波场使用LSQR算法(使用不完全LU分解算法作计算效率对比)，同时使用二阶时间精度四阶空间精度常规时域有限差分法取同时刻波场做结果对比，参数一致[24-27]。图7、图8为傅里叶反变换后时间域的波场快照与时间域有限差分波场快照对比，图9为炮集记录。

图8 180 ms波场快照局部Fig.8 Snapshot of 180 ms local wave field(a)频率域LSQR法波场；(b)常规时间域波有限差分波场

图9 共炮点道集记录Fig.9 Common shot point gather

图10 Overthrust模型Fig.910 Overthrust model

从图7(a)与图7(b)对比以及图8(a)与图8(b)对比图中可以看到，频率域LSQR差分波场和时间域差分波场差别很小，PML对边界反射波的吸收效果很好，没有数值频散，共炮点道集也显示清晰的直达波同相轴和双曲线形态的反射波同相轴，表1为不完全LU分解、LSQR和常规时间域有限差分法的计算效率对比，可以看出LSQR占用内存最小，计算时间最短，节省了大量系统资源。

2.3 Overthrust模型

模型采用较为复杂的Overthrust模型中的部分区域(图10)，网格范围为120×120，最大速度为6 000 m/s，最小速度为2 674 m/s，平均速度为4 292 m/s，网格间距取Δx=Δz=25 m，震源位置网格点为(27，27)，检波器位置x=22，其他参数与上个模型相同，表2为计算效率对比，图11、图12为波场快照对比，图13为共炮点道集记录，在波场对比可以发现频率域LSQR有限差分波场与常规时间域波场只有细微差异，产生差异主要原因是在大网格间距下，四阶精度时间域差分波场会产生轻微数值频散，影响其模拟精度，而优化17点差分波场拥有更好的频散特性，因此模拟精度更高[16]。炮集记录可以看到清晰的直达波和杂乱的反射波曲线，旅行时也与计算模型相吻合，证明了正演模拟的正确性。

表1 均匀层状模型计算效率对比

图11 250 ms波场快照Fig.11 Snapshot of 250 ms wave field(a)频率域LSQR法波场；(b)常规时间域波有限差分波场

图12 250 ms波场快照局部Fig.12 Snapshot of 250 ms local wave field (a)频率域LSQR法波场；(b)常规时间域波有限差分波场

3 结论

本次正演模拟采用优化17点差分算子格式，运用稀疏矩阵LSQR算法，并使用PML完全匹配层吸收边界反射。优化17点格式在大网格间距的条件下表现出很小的频散特性，可以节省大量的系统资源。在简单层状模型和Overthrust模型上的数值实验也证明了LSQR在快速和准确的求解频率域正演波场的同时只消耗很小的内存，高于直接解法的计算效率，因此在全波形反演领域将具有很好的应用前景。

表2 Overthrust模型计算效率对比

图13 共炮点道集记录Fig.13 Common shot point gather