频率域八阶NAD有限差分模拟及全波形反演

韩如冰 郎 超

(①中国地质科学院地质研究所自然资源部深地动力学重点实验室，北京 100037； ②中国地质大学(武汉)地球物理与空间信息学院，湖北武汉 430074)

0 引言

全波形反演(Full-waveform Inversion)以波动方程作为数学模型模拟地震波的传播规律，能够更充分地利用地震数据信息，具有成像更准确、分辨率更高的特点[1-2]。Tarantola[1]首先提出了基于广义最小二乘的时间域全波形反演理论和方法，随后Pratt[3-4]将其推广到频率域内。相比于时间域，频率域全波形反演具有计算过程稳定、占用内存小、不存在累计误差、缓和局部极小、降低反演的非线性、易于处理衰减频散效应、易于并行计算等特点[5-10]，近些年得到广泛发展与应用[11-12]。

在频率域计算中，选择合适的数值方法离散波动方程，处理区域边界、求解大型线性方程组是必要的。常见的数值离散方法包括有限差分法[13-14]、有限元法[15-16]、谱元法[17-18]等。当网格不够密集时，常规的有限差分(Ordinary Finite-difference，OFD)方法容易出现明显的数值频散现象，尽管通过加密离散网格可以在一定程度上压制频散，但是相应的计算量增大，所需时间和内存也会增加[19]。

针对以上问题，Yang等[20-21]引入了近似解析离散化(Nearly Analytic Discrete，NAD)方法，并应用于时间域地震波模拟，随后Lang等[8]将其引入到频率域正演模拟中，其基本思想是利用波场位移与其梯度值的联合近似得到位移的高阶偏导数，进而离散波动方程。相比于其他差分方法，NAD算子携带了更多的波场信息，特别是含有刻画波场变化特征的梯度信息。前人研究表明，针对较粗的离散网格，NAD方法在压制数值频散方面具有很好的效果[8]，此外，与同阶差分格式相比，NAD数值格式计算量小，易于并行计算。然而现有的频率域四阶NAD(NAD4)方法对于计算效率的提升还很有限，需要考虑在此基础上构造更加有效的数值方法。

为提高全波形反演的计算效率，本文构造了频率域八阶NAD(NAD8)格式离散二维声波方程，采用完全匹配层吸收边界(Perfectly Matched Layer，PML)[22]条件。首先详细介绍了NAD8格式的具体离散过程；然后使用一类不精确旋转分块三角预处理算子(IRBTP)加速广义极小残量方法(GMRES)[23]求解离散后的大型稀疏线性方程组；通过数值频散分析以及波场模拟与六阶NAD(NAD6)、八阶OFD(OFD8)方法进行对比，结果表明NAD8方法在压制频散和提高计算效率等方面具有优势；最后，基于共轭梯度法[24]，首次构造了基于NAD8方法的频率域全波形反演算法，对两种典型分层介质模型和经典Marmousi模型进行了反演试算，获得了高保真、高分辨率的成像结果，进而验证了所提方法的有效性和适用性。

1 方法原理

归结起来，频率域NAD方法的离散过程可概括为以下几步：对频率域波动方程关于各个方向求偏导，并在计算区域的边界处施加吸收边界条件，进而得到偏微分方程组，采用NAD网格差分模板离散其中的高阶偏导数项，再按一定节点排序规则形成大型稀疏线性方程组。

假定地球为完全弹性介质时，可以用声波方程研究地震波的传播规律。常密度介质下的二维频率域声波方程为

(1)

式中：u(x，z，ω)表示位移或者压力场；ω=2πf表示角频率；c为地震波波速；s表示震源项。

对式(1)两端分别关于x方向与z方向求偏导，有

(2)

与时间域不同，在频率域施加PML边界条件的过程中，不需要进行反褶积计算，实现更方便、准确。首先引入复坐标控制衰减，以x方向为例，有

(3)

将衰减函数取为

(4)

(5)

式中

(6)

与NAD4方法[19]不同，对于NAD8格式，任一点(i,j)的高阶偏导数需要由(i-2，j)、(i-1，j)、(i，j)、(i+1，j)、(i+2，j)五个点处的波场及其空间梯度的加权组合近似。求解微分方程中高阶偏导数的离散格式为

(7)

(8)

(9)

(10)

128ui-1,j+1-128ui+1,j-1+128ui-1,j-1-

7ui+2,j-2+7ui-2,j-2)+

(11)

1408ui+1,j+1-1408ui-1,j+1-62ui+2,j-

2816ui+1,j+2816ui-1,j+62ui-2,j+

31ui+2,j-2+1408ui+1,j-1-1408ui-1,j-1-

(12)

1408ui+1,j+1-2816ui,j+1+1408ui-1,j+1-

1408ui+1,j-1+2816ui,j-1-1408ui-1,j-1-

31ui+2,j-2+62ui,j-2-31ui-2,j-2)+

(13)

式中h表示空间离散步长。将离散结果依次代入式(5)，并写成统一线性方程形式

(14)

(15a)

(15b)

(16a)

(16b)

(17a)

(17b)

由式(14)的离散结果发现，NAD8差分格式共涉及周围17个节点，显然算子长度小于同阶精度的其他格式[25]，具体差分模板如图1所示。

为快速有效求解离散后的大型稀疏线性方程组，本文采用一类不精确旋转分块三角(IRBT)预处理算子加速广义极小残量方法(GMRES)。该方法是一种专门用于非对称线性方程组求解的Krylov子空间方法，充分利用了系数矩阵T的稀疏结构，结合不精确预处理的思想，在迭代次数并未显著增加的前提下，可节省较多的运算时间与存储量，从而能明显提高计算效率[23]。

图1 NAD8差分网格模板

2 正演计算

从数值频散分析和波场模拟两方面考查NAD8方法相对于NAD6和常规八阶有限差分方法的模拟效率。震源选取雷克(Ricker)子波，根据Fourier变换的性质，其频率域表达式[9]为

(18)

式中：A表示振幅，波场模拟时统一设置为1；f0是主频。根据频谱分析[8]，主要能量分布在0～3f0。定义网格频率(Gf)为每个波长包含的离散网格点数，以度量网格剖分的大小[26]，其计算公式为

(19)

式中w表示波长。从式(19)可以看出，对于同一计算区域，Gf越大，每个波长内网格点数目会越多，从而计算精度越高，但是相应地增大了计算量。

2.1 理论频散曲线分析

首先从理论上论证各种数值方法在压制数值频散方面的能力。基本思想为：通过计算数值速度与实际波速的偏差考察各种离散方法的准确性，偏差越小说明频散越小。

为计算NAD方法的数值速度，首先引入均匀介质中的波动方程及其关于x、z方向的偏导

(20)

与式(5)的离散过程类似，通过NAD方法可将式(20)离散为矩阵形式

(21)

(22a)

(22b)

式中：a=exp(ikxh)；b=exp(ikzh)。

图2 θ=0(a)和θ=π/6(b)时三种数值方法的频散曲线

2.2 简单模型的波场模拟

以均匀介质模型(M1)和双层介质模型(M2)，进一步讨论NAD8方法的计算精度与数值效率。

M1模型的速度v=4km/s； M2模型的速度v1=4km/s、v2=5km/s，界面位于z=4km处，其他参数设置如表1所示。

表1 M1和M2模型的基本参数

图3为频率f=10和30Hz时， 本文构造的NAD8方法针对以上两种模型得到的频率域单频波波场快照。图4为NAD8、NAD6和OFD8三种方法模拟的t=0.5s时刻时间域波场快照。从图中可以看出，在相同网格规模的情况下， NAD8方法有效地压制了数值频散，而后两种方法则出现了明显的频散现象，OFD8方法频散最为严重，表明在相同网格频率下， NAD8方法较其他方法精度较高。根据Shannon采样定理[27]，为了准确地刻画波场，网格频率至少大于2.0， 而NAD8方法突破了这个极限，Gf最小值可到达1.98，这是因为NAD方法除了利用位移信息，还利用了波场梯度信息，大大提高了波场模拟精度，降低了对网格频率的要求。

图3 由八阶NAD方法计算的M1(左)和 M2(右)模型频率域单频波波场快照(a)10Hz； (b)30Hz

图4 三种方法计算的M1(左)和M2(右) 模型t=0.5s时刻的时间域波场快照(a)NAD8； (b)NAD6； (c)OFD8

最后，针对以上两个模型，在不产生数字频散的情况下，统计三种方法运行时间(表2，10次计算结果的平均)，其中NAD8(Gf=1.98)用时最少，比NAD6方法(Gf=2.40)约省时12%，比OFD8方法(Gf=3.50)约省时25%。

表2 不同方法简单模型波场模拟运行时间对比

2.3 复杂模型的波场模拟

选取经典的Marmousi模型[28](图5)检验NAD8方法的适用性。模型尺寸为17.25km×5.63km；网格间距h=18.75m，网格频率Gf=3.0；四周均为10层PML边界，震源位于(8.63km，0.38km)，主频为20Hz。

图6为NAD8方法模拟的频率f=10Hz和15Hz时Marmousi模型的频率域单频波波场快照，图7为t=1.00、1.67s时刻的时间域波场快照。可以看出，对于复杂的Marmousi模型，没有明显的频散，可见NAD8方法适用于复杂介质模型的正演模拟。

图5 Marmousi速度模型

图6 NAD8方法模拟的Marmousi模型单频波波场快照(a)f=10Hz； (b)f=15Hz

图7 NAD8方法模拟的Marmousi模型时间域波场快照(a)t=1.00s； (b)t=1.67s

3 频率域全波形反演

3.1 基于NAD8方法的频率域全波形反演

所谓反演即从初始模型出发，进行正演模拟，匹配计算波场与真实波场，然后不断修正模型，最终达到特定的精度。本文基于NAD8正演算法对经典模型进行反演，反演采用非线性共轭梯度法，其优势在于：与最速下降法相比，不用引入显著的正演计算量便可以得到更加准确、稳健的计算结果，与牛顿类方法相比，不用计算反演目标函数的二阶偏导数矩阵(Hessian矩阵)，收敛性较强。将离散后的频率域声波方程(式(14))改写为

T(ω)u(ω)=s(ω)

(23)

式中：T表示系数矩阵；u表示波场项；s表示震源项，且三者都与角频率ω有关。合成波场与观测波场数据残差的L2范数可表示为

(24)

式中：m为模型参数； “*”表示复共轭；ns为炮数；nr为每炮接收点数。δdij表示第i炮的第j个检波点对应的模拟波场值u与观测波场d的残差。选择合适的方法求解E(m)的最小值是问题的关键，基于梯度算法，通过推导计算可得相邻步数之间的迭代基本格式为

(25)

通过对波动方程关于模型参数求偏导，进而可导出

(26)

式中v表示残差反传波场。可见计算反演目标函数梯度时候，不需要计算矩阵J的具体表达式，整个问题只需计算TTv=δd*。复杂的矩阵计算可以变成简单的求解一次线性方程组的问题，大大降低了计算量。

通过以上过程计算出梯度之后，进而可以确定搜索方向。本文选用计算效率较高的抛物线拟合方法选择步长，即近似认为目标函数是关于步长的二次函数，确定一个试探步长r1，使其满足在目标函数的下降方向上，不断重复搜索步骤，直到目标函数不降反增，确定另一步长r2； 结合初始步长r0，三点确定抛物线，求取极值点便为搜索步长。最后，反演频率的选择参照非等距频率选择策略[29]，这样频率间距随着频率的增加变得越来越大，从而整体上越来越松散，但不损失反演结果的分辨率，使用更少的频率参与反演，可节省计算时间。

此外，定义标准化残差的概念来衡量反演误差的大小，即

(27)

式中：mt表示真实模型； ‖·‖2表示L2范数。显然给定某一频率，κ是一个随迭代步数变化的曲线，曲线越趋近于0表示反演结果越好。

3.2 简单模型反演

基于NAD8方法的正演程序，依据上文所述的频率域反演原理，首先对两个简单的速度模型进行反演与分析。

模型一(图8)参数设置如下：模型尺寸为2.5km×2.5km，背景速度为4.0km/s，中央有一速度为4.5km/s、边长为0.5km的异常体，网格间距h=25m； 将四周设置为10层PML边界，主频设置为10Hz，震源振幅为1.0×105(为了避免波场值太小而引起的机器出错)，初始速度设置为4.0km/s。炮点与接收点布设两组，位于计算区域内。在频率域接收点的多少并不会增大反演的计算量，而震源数目的增多会使计算量成倍增大。为保证计算效率，需要选择适量的震源数目，而尽量多布设接收点，横向间距为一个网格距离。为了得到深部更高分辨率的成像结果，分别设置两排震源和接收器，设置震源数目为上、下各11个，共计22个，纵向上分别位于0.275km和2.225km处，横向上为0.25～2.25km，间距为0.2km。上、下接收点同时接收传播信号。

图8 速度模型一及震源与接收点布设示意图

首先选取5个频率(1、7、13、19、25Hz)，低频勾勒大致轮廓，高频细化内部结构，每个频率的最大迭代步数为30步，将反演过程记为S1，最终f=25Hz的反演结果如图9a 所示。可以看出反演的速度模型与理论模型整体轮廓基本一致，但是分辨率相对较低，特别在靠近速度异常体周边的数值与理论值差异较大。进一步加密震源，每行放置21个，即2×21个，增加反演频数到30个(1～30Hz，间隔1Hz)，每个频率最大的迭代步数为50，将反演过程记为S2，最终f=30Hz的反演结果如图9b所示。与图9a相比，图9b的分辨率明显提高，细节刻画更清晰，对异常体的刻画更准确。

图9 模型一反演结果(a)S1； (b)S2

图10中蓝色曲线是频率取25Hz时S1反演的误差曲线，可以看出，当步数到达最终30时，误差曲线还有下降趋势；红色曲线为频率取25Hz时S2反演的误差曲线，表现为前几步急剧下降，后期逐步平缓，停滞在小于0.1范围内且基本不再下降，1～30Hz中其他频率的误差曲线类似。

图10 模型一反演误差曲线

模型二(图11)为双层模型，尺寸为2.5km×2.5km，上层速度为4.00km/s，下层速度为4.50km/s在分界面处有一速度为4.25km/s、边长为0.5km的异常体，网格间距为h=25m。将四周设置10层的PML边界，震源主频设为10Hz，震源振幅设置为1.0×105。反演的初始速度设置为4.00km/s。接收器尽可能多的放置，间距为一个网格距离，震源数目为上、下各布设41个，共计82个，纵向上分别位于0.275km和2.225km处，横向上从0.25km到2.25km，间距为0.05km。上、下接收点同时接收传播信号。

设置反演频数为30个(1～30Hz，间隔为1Hz)，每个频点的最大迭代次数为100。图12a为f=30Hz的反演结果，可以看出，整体分辨率较高，细节刻画清晰，特别是在速度分界面与异常体的形态反映较为真实，整体效果较为理想；图12b为f=30Hz的反演误差曲线，开始时下降明显，40次后稳定在0.1以内，符合高精度反演的误差曲线形态。可见，本文的基于NAD8的频率域反演方法对于简单的模型一和模型二取得了较为理想的结果。

图11 速度模型二及震源与接收器布设示意图

3.3 Marmousi模型反演

针对Marmousi模型(图5)的反演，参数设置如下：网格间距h=15m，反演区域为3.47km×1.14km，PML层数为10，主频设置为10Hz，震源振幅设为1.0×105。初始模型(图13)是对真实模型进行Gauss平滑后的结果，只具有大致的分层结构。震源布设44个，纵向上位于0.17km处，横向上从0.15km到3.35km， 间距为0.075km。仅在地表放置一排。接收器排列与震源的间隔相同，纵向位置上位于0.15km处，横向上往外多扩了两个网格距离。

选择30个频点(1～30Hz，间隔为1Hz)、每个频点迭代200次进行反演，图14为10和30Hz的反演结果。由图可以看出，10Hz低频反演结果基本给出了模型的大体轮廓，但是整体分辨率较低，出现了多处异常亮点，特别是深部高速异常体的精度较低。30Hz反演结果分辨率明显提高、细节更加清晰，特别是深部高速异常体的刻画更加精确。

图12 模型二反演结果(a)及反演误差曲线(b)

图13 初始Marmousi速度模型及震源与接收点布设示意图

图14 Marmousi模型10(a)和30Hz(b)的速度反演结果

图15是x=2.0km处、频率为30Hz的速度反演曲线与真实、初始速度曲线的对比，可以看出，整体的反演结果较为理想，特别是浅部低速体保真度较高，但是深层精度略有下降，导致反演速度与真实模型速度之差的L2范数略高(约为0.8)，这是因为震源位于地表处，波传播到深层时振幅降低、误差累积，从而分辨率下降。

图16为频率为30Hz的反演误差曲线，可以看出，整体上一开始下降明显，并逐渐稳定在0.1以内，反演效果较好，其他频率反演误差曲线基本相似。可见，基于NAD8的频率域全波形反演方法针对复杂模型依然适用，且分辨率、保真度都较高，效果良好。

图15 x=2.0km处反演曲线与真实模型、初始模型的对比

图16 Marmousi模型反演误差曲线

4 结论

为进一步提高频率域全波形反演的计算效率，本文首次构造了频率域NAD8差分格式离散波动方程，给出了PML边界条件与频率域NAD方法相结合的正演算法，并详细推导总结了整个数值离散过程，得到了大型线性代数方程组，根据方程组系数矩阵的数学结构，采用一类不精确旋转分块三角预处理算子加速Krylov迭代方法，提高了正演计算效率。数值频散分析以及均匀介质和双层介质模型的波场模拟结果表明：①与OFD8方法、NAD6方法相比，NAD8方法计算精度最高，在压制频散方面效果最明显； ②在有效压制频散的情况下，NAD8方法所需网格频率最小(Gf=1.98)，突破了Shannon采样定理可达到的理论最小值(Gf=2.0)，所用的计算时间最短，大约比NAD6方法(Gf=2.4)省时12%，比OFD8方法(Gf=3.5)省时25%。

对于两种典型的分层介质模型及经典Mar-mousi模型，本文基于NAD8格式的频率域全波形反演算法获得了较高分辨率和保真度的结果，结合反演误差曲线验证了方法的正确性和适用性。