Lippmann-Schwinger积分方程广义超松弛迭代解法及其收敛特性

徐杨杨， 孙建国， 商耀达， 孟祥羽， 魏脯力

吉林大学地球探测科学与技术学院， 长春 130026

0 引言

近年来，随着陆地反射地震工区从简单区到复杂区的扩展以及海洋反射地震从浅水向深水的转变，如何利用散射波进行偏移成像就逐渐进入了研究者们的视野(沈鸿雁,2010；卢明辉等,2010).相应地，散射波数值模拟方法研究也得到了前所未有的重视(孙建国，2006；Jakobsen and Wu,2016).一方面，对散射波场数值模拟的研究为认识散射波的时空分布特征与复杂构造之间的关系提供了有效的途径.另一方面，散射波场的快速计算方法为实现全波场偏移和全波形反演提供了坚实的理论基础和强有力的技术支撑(王本锋等,2016；Huang et al.,2020).

与反射波携带了关于反射面的几何形态和反射面两侧的岩石物理信息类似，散射波也携带着与复杂构造和复杂岩性有关的几何和物理信息.因此，如何能在利用常规反射波实现成像的同时，利用散射波实现复杂构造成像就显得十分重要.

利用扰动理论，可将描述散射问题的Helmholtz方程转化成为第二类Fredholm型积分方程，即Lippmann-Schwinger方程；简记为LS方程.从而，散射波数值模拟问题就转化成为了LS方程的求解问题.与常规积分方程一样，对LS方程可以利用各种方法求解，例如Nyström方法、矩量法(MOM)以及迭代法(SOR)等.如果利用迭代法解LS方程，所得到的级数称为Neumann级数.在散射研究领域，也将Neumann级数称为Born级数，简记为BS.如果只取Born级数的前两项，则可得到著名的Born近似.物理上，Born近似用散射体内的入射场代替散射体内的总场.

Born级数仅在小散射体或弱散射的情况下收敛.因此，Born近似是一种弱散射理论.如何将Born近似从弱散射扩展到强散射一直是散射研究中的核心问题之一.为解决这个问题，Wu和Huang(1995)、Wu(2003)将De Wolf近似引入到地震波散射级数中，对散射级数进行重整化处理，建立了声波近似下的MFSB(Multiple Forescattering Single Backscattering)近似.在这种近似中，利用多屏单向波传播算子进行双域计算，有效地解决了中等强度的前向散射计算问题.然而，MFSB近似要求反向散射场远远小于正向散射场，从而并不能用于解决前向散射问题.Yu等(2010)、Yu和Fu(2013)通过将边界元法(BEM)生成的边界积分方程矩阵分解为两部分：自作用算子和外推算子，其中外推算子由Born级数表示，并利用BEM+BS法来求解新的方程组.该算法可以克服边界元法在复杂散射介质情况下离散矩阵无法求逆的困境，并通过层状非均匀介质的地震波场数值模拟验证了该算法的可行性.Jakobsen(2012)将量子物理中和电磁学中常用的的重整化理论应用于地震正演模拟中，采用T算子方法求解Lippmann-Schwinger方程以改善Born级数的收敛性.Jakobsen和Wu(2016)利用De Wolf级数和T传播算子重新推导了重整化后的De Wolf散射级数，并对其收敛特征进行了详细分析.为改善Born级数的收敛性以解决强光学散射问题，Osnabrugge等(2016)通过在背景波数中引入虚部分量，得到了带阻尼因子的背景格林函数.并通过选取合适的预解因子来保证Born级数迭代的收敛性.该方法在光学散射数值模拟中表现出了精度高、计算效率快的优点.Huang等(2020)从重整化角度重新阐释了Osnabrugge得到的收敛Born级数，并通过数值算例验证了该方法在强扰动散射介质的地震波场正演模拟中的适用性.

在20世纪90年代，Kleinman等(1990a，b)、Kleinman和Van den Berg(1992)将解泛函方程的超松弛理论用于解决光学中的强散射问题，得到了收敛的修正Born级数迭代解.在松弛因子等于1时，该修正的Born级数退化为常规的Born级数.根据Kleinman等给出的数值计算结果，超松弛修正Born级数能够有效地解决强光学散射中的正反演问题.Yu和Fu(2014)将基于Krylov子空间投影的广义最小余量法应用于LS方程的求解，并通过与Yu等(2010)、Yu和Fu(2013)的BEM+BS法对比，分析了广义最小余量法的收敛特征.

相较于源于量子力学的散射级数重整化理论，广义超松弛迭代法理论简洁且具有更坚实的数学基础.然而，大多数光学系统中折射率范围在1～2之间，使得光学散射问题中的散射体扰动相对较小，仅根据Kleinman等人的工作难以对LS方程超松弛迭代级数解在地震散射问题中强速度扰动介质波场数值模拟中的可应用性和收敛性进行分析和评价.

本文的目的就是要解决上述问题，对Kleinman提出的LS方程广义超松弛迭代法在反射地震正演问题中的适应性和收敛特点进行分析和评价.在下文中，首先介绍了广义超松弛法的基本思想和基本公式，然后通过数值分析及与有限差分法数值结果进行对比，分析了其在强散射复杂介质下的散射波场数值模拟中的应用效果及适应性，并对算法收敛特点及计算效率进行分析，最后讨论了松弛迭代算法进一步改进的几个方向.

1 LS方程的广义超松弛迭代

频率域标量波动方程(Helmholtz)的解可以形式地表示为(A8)，即：

U(r)=Ub(r)+Us(r)

(1)

式中，U(r)为位移场，Ub(r)为背景波场，Us(r)为散射波场，ω为角频率，vb(r)为背景速度，α(r)为速度扰动，G(r,r′)为背景介质格林函数.

在文献中，式(1)被称为Lippmann-Schwinger方程，简称为LS方程(见附录A).形式上，LS方程 属于第二类Fredholm型积分方程，从而在原则上可以利用迭代级数法或数值分析法求解.

(2)

则LS方程(1)式的算子形式为

LU(r)=Ub(r),

(3)

由于格林函数G(r,r′)中的Hankel函数在r=r′点包含In(k0|r-r′|)项而奇异.根据Fu等(1997)，(2)式右端的积分项在点r∈Ω上是一致连续的，即积分算子K为弱奇异Fredholm型积分算子，LS方程(1)式成为第二类Fredholm型弱奇异积分方程.由紧算子理论可知，包含In(k0|r-r′|)项的弱奇异积分算子K是映L2(Ω)到L2(Ω)的紧算子(详细证明见：吕涛和黄晋，2013，P37定理1.5.4).因此，积分算子L=I-K也是L2(Ω)上的紧算子.应用Fredholm择一定理，可知LS方程存在唯一解.

根据附录B推论1和定理1(Petryshyn,1963)，选取合适的算子B和C使得：

ρ(I-C-1BL)<1，

(4)

则LS方程存在如下收敛的迭代格式：

un=C-1[(C-BL)un-1+BUb]， ∀u0∈L2(Ω).

(5)

1.1 LS方程的超松弛迭代

下面将通过选取不同的算子B和C，介绍LS方程的几种超松弛迭代格式.首先将LS算子方程的迭代格式(5)式改写为如下递推形式：

(6)

且残差满足如下递推公式：

rn=rn-1-LC-1Brn-1.

(7)

若取C=I；B=I，则得到Born散射级数(BS) (Kleinman and van den Berg,1991)：

(8)

根据式(8)，在第n次迭代后的结果为

(9)

若在式(9)中取U0=Ub且仅迭代一次，即可得到 广泛使用的Born近似.当扰动强度较小或散射体较小时，Born级数收敛，从而Born近似具有与较高的精度.然而，当速度扰动强度较强或散射体尺度较大时，Born级数一般为发散级数，Born近似的精度较低.

若取C=I；B=αI(α为常数)即迭代过程中保持因子α不变，则得到松弛因子为常数的超松弛迭代(GSOR-α)格式(Kleinman and van den Berg,1991)：

(10)

其第n次迭代的级数表示为

(11)

若取线性算子C=I；B=αnI(αn为复数)，即每次迭代过程中按照某一收敛原则逐次修正因子αn，则得到逐次超松弛迭代(GSOR-αn)递推形式(Kleinman and van den Berg,1991)：

(12)

及其第n次迭代的级数表示：

(13)

其中，(10)—(11)式和(12)—(13)式中的α和αn分别为复常数超松弛因子和逐次修正复松弛因子.从(6)式可以看出，迭代近似解Un是初始解和残差的线性组合，从空间映射的角度来阐释，即迭代近似解Un存在于初始解和残差项的仿射空间里.

若广义迭代格式(6)式右端项中的修正项不是简单的残差项的线性组合，而是在迭代过程中递归地构造残差，即第n步的残差rn是通过算子L的某个多项式与r0相乘得到且使得所构造的残差在某种内积意义下相互正交时，迭代格式则可以演变为更为复杂Krylov子空间迭代法.若Krylov子空间表示为Km(L,r0)=span{r0,Lr0,L2r0,…,Lm-1r0}，当残差rn与子空间Km正交时，则称为共轭梯度法(Kleinman and van den Berg，1991)，当残差rn取最小范数时称为广义最小残差法(详见Yu and Fu，2014).

通常情况下，散射问题中的积分算子是非自伴紧随算子，而常规的共轭梯度法仅适用于自伴随算子，需要通过对积分算子L乘上他的伴随算子L*来实现对称化，但算子LL*的条件数是算子L条件数的平方倍，或者通过构造预条件算子来加速收敛.广义超松弛迭代(6)式中若取合适预条件算子C使得‖I-C-1BL‖<1，则可以达到与共轭梯度法相当的收敛性质.

1.2 超松弛因子选取

Petryshyn(1963)、Chertock(1968)、Kleinman和Roach(1988)等将超松弛迭代法应用于积分方程的求解，并给出了积分算子L自伴、非自伴随情况下松弛因子的选取范围公式.这些松弛因子的选取公式均需要计算线性积分算子L的谱信息，而在实际问题中，计算线性积分算子L的谱信息十分困难.

为了避免求解积分算子L的谱信息，根据Kleinman和van den Berg(1991)的超松弛迭代在光学散射中的应用经验，通过在迭代过程中最小化残差‖rn‖来得到相对较优的逐次复松弛因子αn.为求出满足收敛性的复松弛因子αn，考虑U0=Ub时的第n次迭代结果，取αn使得‖rn‖在Hilbert空间中取极小.

定义U(r)∈L2(Ω)上的范数和内积为

(14)

根据‖L‖2范数的定义，有：

‖rn‖=‖rn-1-αnLrn-1‖

=min，

(15)

满足式(15)的必要条件为

(16)

即：

(17)

从而得出：

(18)

2 L-S方程超松弛迭代算法

广义超松弛迭代法求解LS积分算子方程，首先需要对弱奇异积分算子L进行离散，常用的离散方法有基于插值投影算子的配置法、基于正交投影算子的Galerkin法以及基于求积公式的Nyström法等.配置法和Galerkin法生成离散矩阵元素的计算量较大，尤其是Galerkin法的每个矩阵元素都需要计算Ω×Ω上的积分，配置法的每个矩阵元素也要计算Ω上的积分.基于求积公式的Nyström法则可以避免生成离散矩阵的大量积分计算，具有较低的计算成本，并且也可以获得满意的效果.因此，在本文中采用Nyström法对积分方程(1)进行离散.下面具体讨论离散过程以及超松弛迭代算法离散化以及超松弛迭代算法实现流程.

2.1 L-S方程离散化

由于积分算子L具有弱奇异性，为降低计算复杂度，本文中的数值试验采用基于常元插值基函数的Nyström法(吕涛和黄晋，2013)来求解LS方程构成的第二类弱奇异Fredholm型积分方程.

(19)

积分号内的被积函数的连续部分用其插值函数代替，即：

(20)

其中权函数为

(21)

由此得出LS方程的Nyström近似Un(r)满足线性方程组：

(22)

采用超松弛迭代法解线性方程组(22)式，求得Un(r)后回代LS式即可得到观测面上各接收点上的散射波场.对于多维奇异积分方程而言，常元Nyström积分法对于区域剖分较灵活性且离散方程构造更简单.因此，本文采用于常元插值基函数的Nyström法来求解LS方程.

(23)

则常元插值算子为

(24)

其中，Mj为Vj的中心.此时：

(25)

超松弛迭代法求解声波散射积分方程，当速度模型复杂度较高、散射体不规则度较大时，散射体积分方程的积分区域Ω选为全空间计算来保证计算精度.而当散射体较为规则且速度模型复杂度较低时，背景介质上速度扰动为零，则积分区域Ω可由全空间可简化为在散射异常体Ωscatter上的积分，即只需要离散散射体即可，从而提高计算效率.

2.2 数值算法描述

超松弛迭代法求解声波散射积分方程并计算散射波场的算法流程图如图1所示，迭代算法的具体实现步骤如下：

(1)选定待计算速度模型，数值模拟参数初始化：离散网格点数Nz×Nx，空间采样间隔(dx,dz)，震源子波主频fpeak，时间采样间隔dt，时间采样点数nt，源检位置(rs,rg)以及最大迭代次数Nitmax.

(2)开始循环频率：计算震源rs到地下介质任意网格点上的入射波场Ub(rs,r)，地下介质中任意两点之间格林函数G(r,r′)，根据积分算子A的离散矩阵形式求出离散积分算子L.

(4)开始循环迭代计算：

①niter=niter+1；

②根据迭代递推格式计算地下所有网格点上的波场Un和残差rn；

④判断niter是否达到最大值：若是，则退出；否则继续步骤(4)；

(5)开始循环接收点：计算背景速度下接收点rg到地下任意网格点上的格林函数G(r,rg)，回代计算当前频率下所有接收点的波场值U(r,rg).

(6)继续循环频率，执行步骤(2)—(5)，直至求得所有频率下的波场值U(r,rg)，进行 FFT逆变换，输出时域散射波场Us.

图1 算法实现程序流程图Fig.1 Flow chart of algorithm implementation

3 数值试验

3.1 球体模型

对简单球状散射体模型进行正演模拟，速度模型为均匀背景介质中含有一半径200 m的高速球状散射体(如图2a).模型大小为3000 m×3000 m，球散射体速度为2500 m·s-1，背景速度vb=2000 m·s-1，震源子波为主频15 Hz的雷克子波，震源位置rs=(1500 m,10 m)，检波器置于z=10 m的界面上，迭代法空间采样间隔为dx=dz=5 m，时间采样间隔为4 ms，为避免空间假频，有限差分法(FD)空间采样间隔为dx=dz=1 m，时间采样间隔为0.2 ms.

分别采用有限差分(FD)、广义逐次超松弛迭代(GSOR-αn)、复常数广义超松弛迭代(GSOR-α)以及Born级数迭代(BS)四种数值方法计算了单炮记录(如图2e、f、g、h)，其中GSOR-α迭代法中松弛因子α=α1=(r0,Lr0)/‖Lr0‖2且在迭代过程中保持不变.从图中可以看出速度扰动强度较低时，三种迭代方法都收敛，得出的数值计算结果均与有限差分正演结果具有较好的一致性.抽出GSOR-αn法的第80道、151道、230道与FD法进行比较(如图3)，结果显示，两种算法得到的结果差别较小，均在允许误差之内.

为进一步分析广义逐次超松弛迭代法(GSOR-αn)的收敛性质，取计算频率f=30 Hz，复松弛因子αn-迭代次数、归一化残差Error-迭代次数(Error<10-6)的变化关系如图2c、d.如图2c，随着迭代次数增加，残差无限接近于0时，复松弛因子在αn=1.0-0.2i附近浮动，趋向于最优选择.从图2d中可以看出，当满足残差Error<10-6时，GSOR-αn法的收敛速度较GSOR-α和BS法更快,所需要的迭代次数更少.图2b的频率-迭代次数关系曲线表明，在0～60 Hz频率范围内，尤其是高频点处，GSOR-αn法迭代收敛速度要明显优于GSOR-α和BS法.

3.2 规则多体散射模型

考虑地下介质中沿水平方向、垂直方向分布大小不一的多个散射体的情况，速度模型如图4a.模型大小为3000 m×1000 m，散射体速度为3000 m·s-1，背景速度vb=2000 m·s-1，震源子波为主频15 Hz的雷克子波，震源位置rs=(1500 m,10 m)，检波器置于z=10 m的界面上，迭代法空间采样间隔为dx=dz=5 m，时间采样间隔为4 ms，有限差分法空间采样间隔为dx=dz=2 m，时间采样间隔为0.4 ms.

图4b为频率f=30 Hz时GSOR-αn、GSOR-α和BS三种迭代方法归一化残差Error随迭代次数的变化关系，结果显示，在迭代相同次数时，GSOR-αn迭代法的残差与BS迭代法的残差小两个数量级.图4c中，频率与迭代次数的次数的关系曲线表明(残差满足Error<10-6)，当f<30 Hz时，三种迭代方法收敛速度性差较小，当f>30 Hz时，BS迭代法收敛速度变慢，当频率f=60 Hz时，GSOR-αn和GSOR-α迭代法所需要的迭代次数约为BS迭代法的1/4.

图2 球状散射体模型的有限差分和迭代法数值结果 (a) 速度模型； (b) GSOR-αn(实线)、GSOR-α(点划线)和BS(虚线)法频率-迭代次数关系； (c) GSOR-αn法复松弛因子αn-迭代次数关系(实部：实线，虚部：虚线)； (d) GSOR-αn(实线)、GSOR-α(点划线)和BS(虚线)归一化残差Error-迭代次数关系； (e) FD数值结果； (f) GSOR-αn迭代结果； (g) GSOR-α迭代结果; (h) BS迭代结果.Fig.2 Numerical results using the FD method and SOR method on spherical model (a) Velocity model; (b) Frequncy-iteration curves; (c) αn-iteration curve; (d) Error-iteration curve; (e) Scattering wavefield using FD method; (f) Scattering wavefield using GSOR-αn method; (g) Scattering wavefield using GSOR-α method; (h) Scattering wavefield using BS method.

图3 FD(实线)、GSOR-αn(虚线)数值结果单道对比图Fig.3 Comparison of single-trace seismic record using FD and GSOR-αn method

图4 多散射体模型的有限差分和迭代法数值结果 (a) 速度模型； (b) GSOR-αn(实线)、GSOR-α(点划线)和BS(虚线)归一化残差Error-迭代次数关系； (c) GSOR-αn(实线)、GSOR-α(点划线)和BS(虚线)法频率-迭代次数关系； (d) FD数值结果； (e) GSOR-αn迭代结果； (f) GSOR-α迭代结果; (g) BS迭代结果.Fig.4 Numerical results using the FD method and SOR method On multi-scatter model (a) Velocity model; (b) Error-iterations curves; (c) Frequncy-iteration curves; (d) Scattering wavefield using FD method; (e) Scattering wavefield using GSOR-αn method; (f) Scattering wavefield using GSOR-α method; (g) Scattering wavefield using BS method.

图5 FD(实线)、GSOR-αn迭代法(虚线)数值结果单道对比图Fig.5 Comparison of single-trace seismic records using FD and GSOR-αn method

图6 盐丘模型的有限差分和迭代法数值结果 (a) 速度模型； (b) GSOR-αn(实线)、GSOR-α(点划线)和BS(虚线)法频率-迭代次数关系； (c) GSOR-αn法复松弛因子αn-迭代次数关系(实部：实线，虚部：虚线)； (d) GSOR-αn(实线)、GSOR-α(点划线)和BS(虚线)归一化残差Error-迭代次数关系； (e) FD数值结果； (f) GSOR-αn迭代结果； (g) GSOR-α迭代结果; (h) BS(1 order)迭代结果.Fig.6 Numerical results using the FD method and SOR method on salt model (a) Velocity model; (b) Frequncy-iterations curve; (c) αn-iterations curve; (d) Error-iterations curve; (e) Scattering wavefield using FD method; (f) Scattering wavefield using GSOR-αn method; (g) Scattering wavefield using GSOR-α method; (h) Scattering wavefield using BS method.

图7 FD(实线)、GSOR-αn(点虚线)、GSOR-α(短虚线)、一阶BS(点划线)数值结果单道对比图Fig.7 Comparison of single-track seismic records using FD, GSOR-αn,GSOR-α and BS(1 order) methods

图4d、e、f、g分别为FD、GSOR-αn、GSOR-α以及BS四种数值方法得到的单炮地震记录，数值计算结果表明三种迭代方法均可以得到与有限差分法相当的结果，且可以清晰分辨出各个散射体所对应的同相轴.随机抽取GSOR-αn和FD法的3道进行对比(如图5)，可以看出，对于散射体较多的地质模型，相对于精细网格下的有限差分数值结果，粗网格的GSOR-αn迭代法同相轴在局部出现同相轴拉伸现象，这一点可以通过加密网格或增加迭代次数进一步减小残差来改善计算结果.

3.3 盐丘模型

下面通过将广义超松弛迭代法应用于散射体尺度较大且速度扰动较强的地下介质中来验证超松弛迭代方法对强散射介质的适应性.模型参数如图6a,背景速度vb=2000 m·s-1，震源子波为主频15 Hz的雷克子波，震源位置rs=(800 m,10 m)，检波器置于z=10 m的界面上，迭代法空间采样间隔为dx=dz=10 m，时间采样间隔为4 ms，为提高计算精度，有限差分法空间采样间隔为dx=dz=1 m，时间采样间隔为0.1 ms.

图6c、d分别为f=30 Hz时的复松弛因子αn的实虚部-迭代次数关系曲线和三种迭代方法的残差Error-迭代次数关系曲线.如图6c，随着迭代次数增加，残差Error不断趋于0的过程中，复松弛因子在αn=0.6-0.5i附近浮动.图6d表明，Born散射级数发散的情况下，超松弛迭代法依旧收敛，且GSOR-αn迭代法收敛速度较快，所需要的迭代次数更少.从图6b的频率域迭代次数的关系曲线可以看出，对于强散射介质，当f>10 Hz时，Born散射级数已经不再收敛，而GSOR-α和GSOR-αn迭代法仍然收敛，并且GSOR-αn迭代法保持了较好的收敛性质.四种数值方法得到的单炮地震记录见图6e、f、g、h，其中图6h为Born散射级数法仅进行一次迭代后得到的散射波场(Born近似法)，数值结果表明有限差分法得到的波场在盐丘上界面散射能量较强，但在岩体上侧的尖点以及盐丘下界面的散射能量较弱，GSOR-α和GSOR-αn迭代法在岩体上下界面及尖点处散射能量较强且同相轴清晰，受高速盐丘散射体的屏蔽影响较小.随机抽取3道进行对比，如图7所示，从图中可以看出，广义超松弛迭代法得到的结果在盐丘上界面与有限差分法子波时间和相位一致，而在尖点散射以及岩体下界面处有限差分法得到的结果振幅较弱.此外，Born近似得到的散射波场在盐丘下界面处的子波出现相位延迟现象，误差较大.

4 计算时间对比

数值试验过程中，球状散射体模型和盐丘模型采用局部积分法，多散射体模型采用全空间积分法进行数值模拟计算.数值计算过程均在微机上实现，CPU：4核，内存：8GB.

四种数值方法计算时间对比如表1中，可以看出，当地下介质中散射体分布情况较为简单时(如，单体金属矿体、盐丘等模型)，由于散射体周围介质中扰动为0，积分法中的积分域退化为散射体本身，超松弛迭代法相对有限差分法计算量较小，具有较高的计算效率.数值计算过程中，虽然GSOR-αn法更新松弛因子会增加一定计算量，但其收敛速度较快所需要的迭代次数较少，相对于传统的BS法并未增加过多的计算时间.而在散射体较多且分布情况较为复杂的情况下，积分法中的积分域扩大，生成线性方程组系数矩阵会占用大量计算时间，相对于有限差分法，广义超松弛迭代法计算效率较低，有待进一步改进算法以提高计算效率.

表1 FD法与GSOR-αn迭代法计算时间对比Table 1 Comparison of calculation time of FD method and GSOR-αn method

5 分析与讨论

文中采用积分方程法求解频率域标量波动方程，通过广义超松弛迭代法来解散射问题的Lippmann-Schwinger积分方程，并将该方法应用于强扰动介质模型的地震散射波场正演模拟计算.首先引入复松弛因子αn来对Lippmann-Schwinger积分方程的积分算子进行修正，得到新的算子方程，并给出了收敛的修正Born散射级数迭代形式.松弛因子αn的选择是根据Kleinman和Van den Berg(1992)给出的通过在迭代过程中最小化剩余函数rn来逐步修正超松弛因子，以加快修正Born散射级数的收敛速度.地震散射波场正演模拟数值结果表明通过广义超松弛迭代法得到的修正Born散射级数在散射强度较小、复杂度较低的速度模型的上(如球状散射体模型)，其收敛速度比传统的Born散射级数更快；在遇到散射强度较大的速度模型时(如盐丘模型)，仍保持了较好的收敛性质，但需要更多的迭代次数来保证收敛性.文中给出的松弛因子是Kleinman等根据一般矩阵理论推出的满足收敛要求的计算公式，并未考虑实际特定物理问题的渐进解和近似解.在今后研究方向上，可以根据实际地震散射问题，寻找更加适合的预解条件以及松弛因子，以提高复杂问题数值模拟的计算精度.

另外，数值实现过程中需在在每次迭代计算之前更新松弛因子αn以加速收敛，故相对于传统的Born级数会增加计算量.并且波场值的精度依赖于迭代的次数，该算法存在一定程度上的计算精度与计算效率上的折中.为节约计算时间和内存、提高计算效率，可从以下几个方面进行改进：(1)考虑到积分算子A的形式为卷积积分形式，可对其进行空间傅里叶变换到波数域中来避免计算复杂卷积积分；(2)频率计算的并行化；(3)迭代计算中进行数据划分与迭代空间相结合，实现超松弛迭代法的模块并行化.

6 结论

Lippmann-Schwinger方程的广义超松弛迭代法是一种针对强散射的数值模拟方法，主要用于解决强光学散射中的正反演问题.传统的Born散射级数受弱散射假设限制，且只适用于短程传播，在解决速度扰动较强、散射体尺度较大的地震散射问题时并不能得到理想的效果.为了改善Born散射级数的收敛性，本文将广义超松弛迭代法应用于地震散射波场正演模拟中，给出了LS方程所满足的Banach空间下线性算子方程的广义超松弛迭代格式以及松弛因子的选取原则，并设计了超松弛迭代算法的具体实现流程.通过在理论模型上的应用效果，验证了超松弛迭代法具有良好的收敛特性，受数值频散影响较小，在网格间距、时间采样间隔较大的情况下可以得到与有限差分相当的数值模拟结果.与散射级数重整化相比，超松弛法具有更坚实的数学基础，可以有效地克服Born级数的弱点，较好的解决地震波强散射问题.超松弛法在散射地震波数值模拟、全波场偏移成像和全波形反演等方面具有较大的应用潜力，值得进行进一步的研究和改进.未来研究拟整合积分方程FFT快速算法和MPI+openmp并行框架，进一步降低存储量、提高计算效率，为后续进行高效的地震正反演方法研究奠定基础.

附录A Lippmann-Schwinger(L-S)方程

考虑频率域标量波动(Helmholtz)方程：

(A1)

式中；ω为角频率，v(r)为介质波速，U(r)为位移场，s(ω)δ(r-rs)为Dirac函数表示的震源项.对速度v(r)进行分解，使得：

(A2)

vb(r)为背景速度，α(r)为速度扰动.

将(A2)代入(A1)，得到：

-s(ω)δ(r-rs)，

(A3)

为得到方程(A3)的形式积分解，将总场U(r)写为U(r)=Ub(r)+Us(r)的形式，Ub(r)为背景波场，Us(r)为散射波场.则式(A3)可以写为

(A4)

令Ub(r)满足方程：

(A5)

则散射场Us(r)满足方程：

(A6)

均匀背景介质空间的格林函数G(r,r′)满足Helmhlotz方程:

(A7)

利用第二格林定理，全空间假设下总场的积分方程表达式为

U(r)=Ub(r)+Us(r)

(A8)

式中：G(r,r′)为二维全空间格林函数，表达式为

(A9)

(A8)式为散射问题的形式解，称为Lippmann-Schwinger(LS)方程.

附录B 解算子方程的超松弛迭代法

迭代法是求解线性算子方程Au=f近似解的一种重要方法，其中算子A可以是任意线性赋范空间X中的矩阵、积分算子或其他抽象算子，f为X中的已知项.迭代法的基本思想是：从任一初始解u0出发，按照某一规则，逐次构造一个迭代解un，当un收敛于u*时，使u*是所给算子方程Au=f的解.对于迭代级数解u*存在性和唯一性，有如下定理Petryshyn(1962)：

定理1对任意f∈X，方程Au=f存在唯一解u*，当且仅当存在连续可逆算子C和有界线性算子B，使得级数：

(B1)

定义1设线性算子T的特征值为λi(i=1,2,…)且(λI-T)不存在有界逆算子，则称λi为T的谱点，用σ(T)表示全体谱点集合；称ρ(T)=supλi∈σ(T)|λi|为T的谱半径.

推论1若线性算子T满足下列任一条件：

(1) ‖T‖<1;

则对任意f∈X，方程Au=f存在唯一解u*.

定理2如果对任意f∈X，方程Au=f存在唯一解u*，则迭代格式：

un=C-1(Dun-1+Bf)，∀u0∈X，

(B2)

收敛于u*，且有误差估计：

‖un-u*‖≤‖(BA)-1D‖‖un-un-1‖.

(B3)

迭代格式(2)为求解线性方程组的广义迭代格式，通过选取不同的算子B和C就可以得到不同的迭代方法.从公式(B2)可以看出，线性可逆算子C相当于预条件算子，构造合适的预条件算子使得T=I-C-1BA=I-Lε(‖Lε‖<1)，将会加快迭代格式(B2)收敛速度，如取B=αI(α为常数)，此时预条件算子C可以看作是A的逆算子.相较于Krylov子空间迭代法(如共轭梯度法、广义最小余量法等)，广义超松弛迭代法数学实现简单，计算量较小，通过选择合适的线性算子B和C，该方法可以达到与共轭梯度法相当的收敛速度.如果把Lippmann-Schwinger方程中的积分算子视为Banach空间中的映射算子，可以将广义超松弛迭代法用于地震散射问题的LS积分方程.