散射波积分方程的Adomian分解解法

汪燚林，董良国

同济大学海洋地质国家重点实验室, 上海 200092

0 引言

在目前的地震勘探中，主流反演框架还是基于合成地震数据与观测地震数据的最佳匹配，最经典的反演方法仍然是基于最优化理论的梯度导引类方法.这类地震反演框架经常要面对的一个关键问题是：已知背景模型和背景地震波场，如何更精确地计算一定模型参数扰动下的摄动波场.

地震波散射场与介质参数扰动之间的关系，通常可以通过Green函数的积分方程表达，最为常见的是Lippmann-Schwinger积分方程与Ricatti积分方程.从反演角度看，这些积分方程是非线性的，通过对积分项中的Green函数及总波场进行某些近似来获得其近似解，其中，Born和Rytov近似应用最为广泛.这两种近似在散射波模拟(Devaney, 1981; Snieder and Lomax, 1996; Spetzler and Snieder, 2001, 2004)、层析成像(Woodward, 1992; Liu et al., 2009; Feng et al., 2020b)以及走时和波形反演(Luo and Schuster, 1991; Tarantola, 1984; Pratt et al., 1998)中被大量应用.

然而，对散射场和参数摄动之间实施这样简单的线性近似，直接影响了散射波场的模拟精度，必然会降低反演的精度和分辨率，也直接导致了目前地震波形反演方法对初始模型的强烈依赖这个痼疾.在模拟散射场时，如果能够根据不同情况和需求，利用精度更高的模拟手段，突破传统的对于散射波场的一阶近似，有助于提高复杂介质中散射场的模拟精度，就有可能克服上述成像反演中存在的问题(Innanen, 2009; Albertin et al., 2013; Wu and Zheng, 2014; Jakobsen and Wu, 2016).

要将散射序列中的高阶项应用到后续的反演成像中，需要研究散射序列的求解方法，并从正问题的角度认识散射序列中高阶项的物理含义和作用.

在现有文献中，关于Born散射序列和Rytov散射序列的求解，主要是基于Helmholtz方程或者积分方程，通过相应物理量的级数展开而得到的.对非线性Lippmann-Schwinger积分方程进行迭代求解，就可以将非线性的积分方程表达为显式求和过程的Born散射序列解(Yura et al., 1983; Tsihrintzis and Devaney, 2000a; Osnabrugge et al., 2016; 符力耘, 2010; 董兴朋等, 2016), 这种Born散射序列的求解过程相对简单，但理论支撑不足.基于波场的Rytov变换表达式，结合非均匀介质Helmholtz方程和复相位摄动函数的级数展开，可以得到Rytov散射序列解(Schmeltzer and Robert, 1967; Yura et al., 1983; Tsihrintzis and Devaney, 2000b; Kim and Tinin, 2009),这一类求解方法思路简单，但推导求解过程比较繁琐.另外，Marks(2006)从Rytov变换的极限表达式出发并结合复相位函数的级数展开，推导了一种包括了传统的Born序列和Rytov序列在内的混合序列，通过控制混合参数可以得到不同的散射序列.然而，该方法的推导过程复杂，且并未直接得到Born散射序列和Rytov散射序列.最近，Feng等(2018, 2020a)基于Ricatti非线性积分方程发展了一种在时间域迭代求解Ricatti积分方程的新方法，且得到了一种形式的Rytov序列解.

由Adomian发展起来的分解算法是一种高效求解线性或者非线性积分方程的方法(Adomian, 1990)，它的基本思想就是将积分方程的解表达成一个无穷级数和的形式，并利用一个多项式迭代求解，从低阶至高阶独立的求得级数表达式中的各项，进而得到不同精度的逼近解.该方法很好地保留了原方程的物理性质，自该方法提出以来，很多学者对这一算法进行了应用并证明了其解的收敛性(Adomian, 1990; Cherruault et al., 1992; Cherruault and Adomian, 1993; Wazwaz, 1997).同时，针对该方法在计算量、收敛性以及对于复杂问题的处理等方面存在的一些问题，国内外学者也相继提出了一些改进的方法(Wazwaz, 1999; Wazwaz and El-Sayed, 2001; Duan and Rach, 2011; Zhang and Liang, 2018; 解烈军, 2013; 郭银凤, 2019).总的来说，在面临一些非线性问题时，Adomian分解方法被认为是一种适应范围广泛、计算过程简单且收敛速度较快的方法(陈一鸣等, 2013; 牛红玲等, 2013).

在勘探地震领域，研究如何更好地数值求解散射波积分方程、更精确地计算地震波散射场，不但对于认识和理解地震波传播过程中的散射机理，更好地研究散射波场和介质扰动之间的关系有重要的理论意义，对于提高反演和成像精度也有重要的现实意义.因此，本文将Adomian分解算法应用至散射波积分方程的求解，实现了在统一理论框架下高效求解Lippmann-Schwinger积分方程和Ricatti积分方程，得到了地震波Born散射序列和Rytov散射序列，为下一步发展非线性反演方法提供了精度更高的散射波场正演手段.

1 散射场的积分表达

在空间位置r处介质速度为c(r)的地下模型中，在rs处激发的地震波场u(r;rs)满足标量波方程

(1)

u(r;rs)=

式中的g(r;r′)为背景介质下的格林函数.

从反演介质参数的角度看，Lippmann-Schwinger积分方程(2)为非线性Fredholm积分方程.但从波场正演的角度看，由于(2)式中的积分函数与u(r′;rs)之间为线性关系，因此Lippmann-Schwinger方程实际上是第一类线性Fredholm积分方程.

设φ(r;rs)为复相位摄动函数，将总波场和背景波场的关系表达为

(3)

则(1)式的解可以表达为下列Ricatti积分方程的形式(Wu, 2003)：

(4)

上述Ricatti积分方程中的积分函数与φ(r;rs)之间为非线性关系.因此，无论是从波场正演还是从地震反演的角度看，Ricatti方程都是均匀第一类非线性Fredholm积分方程，这一点与Lippmann-Schwinger积分方程不同.

通过求解上述Lippmann-Schwinger积分方程，可以确定散射场为Δu=u-u0;通过求解上述Ricatti积分方程，可以确定散射场为Δu=u0(eφ-1).

2 散射波积分方程的Adomian分解解法

2.1 求解非线性积分方程的Adomian分解方法简介

Adomian分解方法(Adomian decomposition method)是计算数学中求解Fredholm积分方程的常用方法(Adomian, 1990)，下面对该方法进行简要介绍.

将非线性Fredholm积分方程

(5)

中的u(x)表达为下列分解序列：

(6)

而非线性项F(u(t))是u(t)的非线性函数，将其表达为下列形式的Adomian多项式

(7)

其中，n=0,1,2,….

将(6)和(7)式代入(5)式，有：

(8)

进而可以得到非线性Fredholm积分方程(5)的序列解：

(9)

2.2 利用Adomian分解方法求解Lippmann-Schwinger积分方程

由于线性积分方程是非线性积分方程的一个特例，下面我们把上述Adomian分解方法应用于求解(2)式的Lippmann-Schwinger方程的散射波序列解.

利用(7)式的Adomian多项式，可得：

(10)

将以上各式代入(9)式，即利用Adomian分解方法，可得：

u0(r;rs)=u0(r;rs),

⋮

(11)

上式中，n满足n≥0.这样，就可以得到Lippmann-Schwinger积分方程的解为

+u2(r;rs)+…,

(12)

上述(12)式实际上就是标量波动方程(1)的Born序列解.其中，u0(r;rs)为背景场，而u1,u2,u3,…分别称为背景场u0(r;rs)的一阶、二阶、三阶以及更高阶的校正项，即总的散射场可以表达为

+u3(r;rs)….

(13)

可以看到，一阶校正项u1(r;rs)即为常规的Born近似散射场，而u2(r;rs)、u3(r;rs)、…、un(r;rs)则称为Born散射序列解中的高阶校正项.

2.3 利用Adomian分解方法求解Ricatti积分方程

下面利用Adomian分解方法求解(4)式的Ricatti方程的散射波序列解.

(4)式的Ricatti积分方程又可写为

(14)

利用(7)式可得：

(15)

将以上各式代入(9)式，即利用Adomian分解方法，可以得到非线性Ricatti积分方程的Rytov序列解：

(16)

由上述方程(16)依次计算得到各阶校正项φm(r;rs)(m≥0)后，就可以得到包含Rytov散射序列中M阶项的散射场：

(17)

其中，M≥1.

当M=1时,(17)式即为常规的Rytov近似散射场Δu(r;rs)=u0(r;rs){exp[φ0(r;rs)]-1}.而当M>1时，(17)式就是对传统Rytov近似散射场进一步校正的高阶Rytov序列逼近的结果.显然，计算Rytov散射序列中第n阶项要用到Rytov序列中第0项至第n-1项，而计算Born散射序列中的第n阶项un(r;rs)只需要用到Born散射序列中的第n-1阶项un-1(r;rs)，因为un(r;rs)是当前背景模型下的散射源un-1(r′;rs)ε(r′)产生的散射场.

3 散射序列解的数值计算方法

利用波动方程(1)，对真实速度模型和背景速度模型分别进行两次地震波模拟，两次模拟结果之差即为该模型扰动下的真实散射场.

理论上，从Born序列解(11)式可以发现，在当前背景模型c0(r)和介质扰动ε(r)的基础上，分别以ui(r;rs)(i=0,1,2…)作为背景场进行一次地震波模拟，就可以得到散射场的Born序列解中的校正项ui+1(r;rs)(i=0,1,2…).从Rytov序列解(16)式可以发现，在当前背景模型c0(r)和介质扰动ε(r)的基础上，以φi(r;rs)(i=0,1,2,…)的不同组合作为复合震源进行一次地震波模拟，就可以得到散射场的Rytov序列解中的校正项φi+1(r;rs)(i=0,1,2,…).

实际计算中，结合Born散射序列表达式(11)和Rytov散射序列表达式(16)，可以建立起两个散射序列解中的各阶校正项之间的对应关系：

(18)

在(18)式中，各等式左侧为Born序列中的各阶校正项(u1,u2,u3,…)与背景波场(u0)，右侧为Rytov序列中的各阶校正项(φ0,φ1,φ2,…).在本文中，先计算出Born序列解中的各阶校正项ui，再根据(18)式，可以求得Rytov序列解对应的各阶校正项φi，再分别根据(13)和(17)式便可求得对应不同阶的Born序列逼近和Rytov序列逼近下的散射场.

4 数值算例

显然，通过Adomian分解解法得到的散射波场的Born序列解和Rytov序列解，其中的一阶校正项就是传统的散射波场的Born近似解和Rytov近似解.在满足一定的条件下，在一阶近似解的基础上，引入Born序列解和Rytov序列解的高阶项，可以更精确地描述地震波散射场.

下面利用三个模型试验，通过和传统一阶近似结果对比，分析散射序列解中的高阶项对一阶近似散射场的校正作用.

4.1 模型一

首先选取了一个一维模型(图1)进行测试，震源和接收点分别位于模型的顶端和底端，因此这里仅分析前向散射波场.模型深度1.5 km，真实速度在1900 m·s-1和2100 m·s-1之间变化(图1中蓝线)，背景模型速度为常数2000 m·s-1(图1中红虚线).

图1 一维速度扰动模型与观测方式(蓝线：真实速度；红虚线：背景速度)模型顶端星号表示炮点位置，模型底端三角形表示检波点位置.

各阶Born、Rytov序列逼近的前向散射波形和真实散射波形对比结果如图2、图3所示.由于多次前向散射能量相对较弱，因此将1～2 s的散射波记录单独显示在图2b和图3b中.可以看到，无论是对Born散射序列还是Rytov散射序列，不管是能量比较强的一次前向散射还是能量相对较弱的多次前向散射，在一阶近似(即Born近似和Rytov近似)的基础上，添加各级高阶项会使得序列解(包括走时和振幅)更加逼近真解，散射场的描述精度得到提高.

另外，从图2和图3还可以看出，一阶近似几乎没有刻画多次前向散射的能力.考虑高阶校正项后，除了可以对一次前向散射进行校正以外(见图2a和图3a)，还可以明显提高对多次前向散射的刻画能

图2 各阶Born散射序列逼近解与散射波真解对比

图3 各阶Rytov散射序列逼近解与散射波真解对比

力(见图2b和图3b).

4.2 模型二

第二个试验选取了一个二维高斯球速度模型(图4)进行测试.半径为150 m的高速异常体位于模型中心，背景为2000 m·s-1的常速模型，异常的速度扰动最大值为300 m·s-1.点震源位于模型上方(1500,100)m处，在距离高速异常体中心半径为1300 m的右半圆周上大约每10°设置一个检波器，共记录19个不同角度的散射波信号，以便观测不同散射角度上的散射波场.

图4 高斯球速度模型与观测系统

由于不同散射角度散射振幅相差很大，为更清晰地对比各个角度不同阶散射波形和真实散射波形的匹配情况，对不同角度的散射波进行f′i(t)=fi(t)/max(abs(fi(t)i=1,2,3))的归一化.其中，i=1,2,3分别表示真实散射场、一阶近似散射场和高阶散射序列逼近的散射场，fi(t)和f′i(t)分别是同一散射角度归一化前后的散射波形，max(abs(fi(t)))为该散射角不同散射波数据绝对值的最大值.

图5、图8分别显示了19道不同角度的不同阶Born和Rytov序列逼近的散射场和真实散射场的归一化后的波形对比结果，图6、图9分别对应图5、图8中高阶Born序列和高阶Rytov序列逼近解的各道归一化因子max(abs(fi(t)i=3))的对数随角度的变化，它体现了不同散射角度的散射场振幅的相对变化.可以发现，在常速背景模型情况下，这个高斯球速度异常在不同的散射角度上产生的散射波场能量差别很大，小角度的前向散射波振幅和大角度的背向散射波振幅相差近4个数量级.小角度的前向散射波能量强、频带宽、高频成分丰富；随着散射角的增大，散射波能量逐渐变弱、频带变窄并向低频移动.

图5 不同散射角归一化后的真实散射场、一阶Born近似和高阶Born序列逼近散射场

图6 图5中高阶Born序列逼近解的各道归一化因子的对数随角度的变化

图7 40°～90°真实散射场与不同阶Born逼近的散射场波形对比

图8 不同散射角归一化后真实散射场、一阶Rytov近似和高阶Rytov序列逼近的散射场

图9 图8中高阶Rytov序列逼近解的各道归一化因子的对数随角度的变化

为进一步显示地震散射序列中各高阶项对于一阶Born和Rytov近似的校正作用，图7、图10分别显示了散射角约40°～90°的不同阶Born和Rytov序列逼近的散射波形和真实散射波形的对比结果.可以看到，高斯球速度扰动模型下，对于大部分散射角而言，在一阶近似(即Born近似和Rytov近似)的基础上(红线)，添加散射序列的二阶项后(蓝线)的结果与真实散射场更加逼近，再添加散射序列的三阶项后(绿线)精度进一步提高.可见，对该模型而言，在许多散射角度上，逐步添加各级高阶项计算得到的散射波场(包括走时和振幅)会逐渐逼近真解，散射波场的模拟精度逐渐提高，说明了考虑散射序列中的高阶项是很必要的.特别是在一阶近似(即Born和Rytov近似)的基础上，二阶项和三阶项的作用尤其明显.

图10 40°～90°真实散射场与不同阶Rytov序列逼近的散射场波形对比

4.3 模型三

第三个试验选取了一个相对复杂的Sigsbee速度模型(图11b)进行测试.背景速度为真实模型的平滑模型(图11a).点震源位于模型上方(2,0.1)km位置处，这里，我们主要关注前向散射波，因此，在模型下方深度为1.4 km位置处间隔10 m共放置了400个检波器(图11b)，用于记录模型的前向散射序列波场.

图11 (a)Sigsbee平滑背景与(b)真实速度模型及观测方式

图12显示了Sigsbee模型在当前背景速度下的真实前向散射波场记录，图13为不同阶Born序列解的前向散射波形记录.当前背景速度条件下，可以发现，一阶Born近似所描述的波现象和真实散射波场记录基本一致，另外，Born散射序列的二阶项、三阶项对于散射场的刻画仍然有一定贡献，且能量上呈逐渐减弱的趋势.进一步从图14中的两个检波器接收的散射波记录可以看到，Born散射序列的高阶项是对于一阶Born近似的校正，不断添加Born序列的高阶项，可以得到更精确的散射波场模拟结果.

图12 真实前向散射波场记录

图13 (a)一阶Born近似前向散射波场(u1)；(b)Born序列二阶前向散射波场(u2)；(c)Born序列三阶前向散射波场(u3)(均在图12所示的同一色标范围内显示)

图14 在水平位置1.0 km(a)及3.5 km(b)处检波器接收到的真实散射场与Born序列逼近波场

图15a显示了一阶Rytov近似前向散射波记录,图15b显示了二阶Rytov逼近和一阶Rytov近似的前向散射波记录的差，图15c显示了三阶Rytov逼近和二阶Rytov逼近的前向散射波记录的差.当前背景速度条件下，和Born散射序列的前三项类似，Rytov散射序列的前三项体现在散射波记录上的能量逐渐减弱(图15a—c).另外，进一步从图16中两个检波器接收到散射波记录可以看到，当前背景速度下，Rytov散射序列的高阶项是对于一阶Rytov近似的一种校正，不断添加Rytov散射序列的高阶项，可以得到更精确的散射波场模拟结果.

图15 (a)一阶Rytov近似前向散射波场；(b)二阶Rytov逼近与一阶Rytov近似前向散射波场之差；(c)三阶Rytov逼近与二阶Rytov逼近的前向散射波场之差(均在如图12所示的同一色标范围内显示)

图16 在水平位置1.0 km(a)及3.5 km(b)处检波器接收到的真实散射场与Rytov序列逼近波场

在模型相对复杂，背景模型较好的情况下，逐步添加散射序列的高阶项计算得到的散射波场(包括走时和振幅)会逐渐逼近真解，散射波场模拟的精度逐渐提高，说明了在一定条件下，复杂介质情况下考虑散射序列中的高阶项是必要的.

5 讨论

本文通过Born及Rytov散射序列来表达散射波场，这必然会涉及序列是否收敛的问题，只有序列级数收敛才有实际物理意义.上面三个数值试验说明：在速度异常体尺度不太大、速度扰动量不太大、模型本身没有非常复杂以及背景速度模型不太差的条件下，散射场的这两种序列解还是稳定收敛的，也就是说，通过引入Born序列解和Rytov序列解的高阶项，可以更精确地描述真实的地震波散射场.

我们发现，这两个散射序列级数的收敛性与介质模型的复杂程度、背景介质模型的精度、介质扰动的几何尺度和扰动强度均有关，也与观测的空间位置(如不同散射角)有关.譬如，对散射角较大的背向散射和散射角较小的前向散射，这两个散射序列级数就表现出不同的收敛性质，可见，地震波场散射序列的收敛性是一个十分复杂的问题.Jakobsen和Wu(2016)用强扰动介质的数值实验揭示了Born散射序列在某些条件下不收敛的现象，但是，文章也没有定量的给出Born散射序列的收敛条件，另外，文章没有涉及本文讨论的Rytov散射序列.Born和Rytov散射序列在理论上本身存在一定的关联性，我们发现，两者在收敛条件方面也具有很大的一致性.然而，影响散射序列收敛性这一问题的因素众多，不同介质情况下，很难定量给出散射序列能够对散射场进行精确描述(散射序列收敛)的条件.总的来说，在介质扰动的几何尺度不太大、扰动量不太大(或背景模型不太差)，散射波强度不太强的条件下，散射序列还是趋于稳定收敛的.因此，在实际应用中，应当根据具体的物理问题，仔细研究散射序列的收敛性质，明确散射序列的实际物理含义，再加以合理运用.

结合前文中数值算例部分的内容，我们在该部分进一步设计实验，仅考察速度扰动量的大小这一因素对于散射序列收敛性的影响，对散射序列存在的收敛性问题进行举例论证，以明确散射序列的收敛具有严格的收敛条件.

选取类似模型测试二中的二维高斯球速度模型进行对比测试.其他参数与模型测试二相同，只是将图4中模型的速度最大扰动由300 m·s-1(15%)提高到1600 m·s-1(80%).在这个较差的初始速度模型条件下，图17显示了部分散射角的不同阶Born散射序列逼近和真实散射波的记录对比结果，图18显示了部分散射角的不同阶Rytov散射序列逼近和真实散射波的记录对比结果.在模型扰动量增大至80%的情况下，首先可以看到，从小角度散射到中角度散射，散射波的能量(振幅)逐渐降低，其次，在散射角相对较大(50°和60°)、散射波能量相对较弱时，此时的一阶近似(Born近似及Rytov近似)相对真解引入了较大的近似误差，但此时引入相应的高阶项仍然实现了对于一阶近似的校正.然而，在当前介质强扰动情况下，在散射角较小(10°～40°)、散射波能量较强时，引入散射序列中的高阶项使得散射序列解更加偏离了真解，说明这时的散射序列不再稳定收敛，散射序列逼近解也失去了相应的物理含义.

图17 10°～60°真实散射场与不同阶Born序列逼近的散射场波形对比

图18 10°～60°真实散射场与不同阶Rytov序列逼近的散射场波形对比

目前，关于散射序列问题的研究主要集中于对正问题的讨论，一些学者尝试通过对传统散射序列实施重整化(Wu and Zheng, 2014; Jakobsen and Wu, 2016)、数学变换(Osnabrugge et al., 2016; Eftekhar et al., 2018)或者发展新的散射序列(Feng et al., 2018, 2020a; Jakobsen et al., 2020)等思路，以改善散射序列的相关特性.其中，Jakobsen 和 Wu(2016)也是从正问题的角度，考虑介质参数存在强扰动时，设计数值实验系统的讨论了De Wolf序列在收敛特性(精度)方面相较于Born散射序列的优势，并指出，这种改进的散射序列表达不仅有助于我们理解强散射介质中的波现象，也是进一步发展非线性反演方法的基础.

反问题方面，目前的地震成像和反演通常基于的Born近似，相当于一阶线性近似，从而产生了线性Fréchet导数，这就是目前占统治地位的基于局部线性化的迭代反演，其中，全波形反演(FWI)就是一个典型代表.然而，目前地震波形反演中出现的几个致命问题，譬如对初始模型的强烈依赖问题，其根源就在于对介质参数扰动和波场摄动之间的关系实施了简单的线性近似(Born近似或Rytov近似).实际地震波场与介质参数之间是强烈的非线性关系，如果能够很好地利用散射序列中的高阶项，就有可能规避目前成像和反演中存在的诸多问题.

6 结论

地震反演的精度很大程度上取决于正问题的描述精度，建立更精确的模型扰动和波场扰动之间的关系十分关键.目前，基于Lippmann-Schwinger积分方程和Ricatti积分方程近似下建立的散射场的一阶近似方法直接影响了后续成像和反演的精度.本文将计算数学领域中的Adomian分解方法应用至求解散射场的Lippmann-Schwinger和Ricatti积分方程中，简洁、高效地得到了散射场积分方程的序列解.数值试验结果表明，如果采用经典的一阶Born或者Rytov近似，某些情况下难以保证近似散射场描述的精度.在一定的条件下，散射序列解收敛时，与传统的Born和Rytov近似解相比，引入散射序列中的高阶项可以更精确地描述地震波散射场.