基于预条件广义逐次超松弛迭代法的数值格林函数计算方法

摘要：为了改善Born散射级数解决地震强散射问题时的收敛性，将带有虚部分量的复波数格林函数引入到求解格林函数Lippmann–Schwinger（LS）积分方程数值解的广义逐次超松弛迭代法中，弱化格林函数的奇异性。引入预条件算子降低系数矩阵的条件数，加速迭代级数的收敛速度，给出了复波数LS方程的预条件广义逐次超松弛（preconditioned generalized successive over-relaxation， Pre-GSOR）迭代格式。通过数值分析和收敛性分析重新选取合适的衰减因子和预条件算子，得到了满足地震强散射条件的收敛Born级数，并将其用于地震强散射问题中数值格林函数的计算。数值结果表明：复波数LS方程Pre-GSOR迭代法可以得到与实波数LS方程直接法相匹配的数值模拟结果；复波数LS方程Pre-GSOR迭代法系数矩阵条件数在高频时仅为原系数矩阵条件数的10%，相同迭代次数下归一化收敛残差可降低3个数量级以上，且对高频适应性强，可有效改善实波数LS方程广义超松弛迭代法在强散射介质中的收敛停滞问题。

关键词：地震散射波场；格林函数；广义超松弛迭代；预条件算子

doi：10.13278/j.cnki.jjuese.20230249

中图分类号：P631.4

文献标志码：A

徐杨杨，商耀达，孙建国. 基于预条件广义逐次超松弛迭代法的数值格林函数计算方法. 吉林大学学报（地球科学版），2024，54（5）：16961710. doi：10.13278/j.cnki.jjuese.20230249.

Xu Yangyang， Shang Yaoda， Sun Jianguo. A Preconditioned Generalized Successive Over-Relaxation Iterative Method for the Numerical Green’s Function Method. Journal of Jilin University （Earth Science Edition）， 2024， 54 （5）： 16961710. doi：10.13278/j.cnki.jjuese.20230249.

收稿日期：20231008

作者简介：徐杨杨（1991—），女，博士研究生，主要从事地震散射波场数值模拟、成像技术与快速算法研究，E-mail：846208564@qq.com

通信作者：孙建国（1956—），男，教授，博士生导师，主要从事地下波动理论与成像技术、计算地球物理、科学计算方法与技术、海洋反射地震资料处理、地下天线理论以及岩石物理等方面的教学和研究工作，E-mail：sun_jg@jlu.edu.cn

基金项目：国家自然科学基金项目（41974135）

Supported by the National Natural Science Foundation of China （41974135）

A Preconditioned Generalized Successive Over-Relaxation Iterative Method for the Numerical Green’s Function Method

Xu Yangyang， Shang Yaoda， Sun Jianguo

College of GeoExploration Science and Technology， Jilin University， Changchun 130026， China

Abstract：

Born scattering series is often limited by weak scattering assumptions when solving strongly seismic scattering problems， resulting in slow convergence or divergence. A simple and effective way is to improve the iterative algorithm using numerical analysis. One such method is the generalized successive over-relaxation （GSOR） iterative method， which can be applied to solve the Lippmann-Schwinger （LS） equation and obtain the desired convergent Born scattering series. However， in strongly heterogeneous media， the GSOR iterative method may also face the challenge of slow convergence speed while calculating the high-frequency Green’s function. In this paper， the complex wavenumber Green’s function is utilized with the GSOR iterative method to numerically solve the LS equation of the Green’s function. The complex wavenumber has imaginary components that enable localizing the energy of the background Green’s function and exponential decay， reducing the singularity of the background Green’s function. To reduce the condition number of the coefficient matrix， we further introduce the preconditioning operator and provide a preconditioned generalized successive over-relaxation （Pre-GSOR） iteration format. The convergent iteration series is obtained by selecting an appropriate damping factor and preconditioning operator. Then it is used to calculate the numerical Green’s function in the seismic strongly scattering media. Numerical results indicate that the Pre-GSOR iteration method for the complex wavenumber LS equations can produce simulation results consistent with those obtained by direct methods for real wavenumber LS equations. The condition number of the coefficient matrix in the Pre-GSOR iterative method for the complex wavenumber LS equation is only 10% of the original condition number at high frequencies. Under the same number of iterations， the normalized convergence residual obtained by this method can be reduced by more than three orders of magnitude. The new method exhibits lower convergence error， better convergence， and strong adaptability to high frequencies， effectively mitigating the convergence stagnation problem encountered by the generalized over-relaxation iterative method for the real wavenumber LS equation in strongly scattering media.

Key words：

seismic scattering wave field; Green’s function; generalized over-relaxation iteration; preconditioning operator

0" 引言

格林函数的表示和计算通常是实现基于Kirchhoff型、Born型等积分方程地震正反演与偏移成像方法的核心，格林函数表示方法不同则可得到不同的正反演与成像方法［1］。常用的计算格林函数的方法主要有基于渐近射线理论的方法［29］、基于微分算子的纯数值方法［1014］和基于散射理论的方法［1520］等。渐近射线理论以波动方程的高频近似假设为前提，利用零阶几何射线理论或高斯波束叠加等波场高频表示来近似格林函数，即高频渐近格林函数［6］。对于与地震波长相比较大的光滑模型，渐近射线理论非常有效，具有运算效率高、可控性好的优点，但渐近射线理论认为只有一次反射才是有效信号，无法利用多次反射、绕射及散射等波动现象所携带的有效信息［1］。而纯数值方法和散射理论均可以用于计算包含完整波形的格林函数。相比于纯数值方法，散射理论具有更高的计算效率［18］。

基于散射理论表示格林函数时，通常将描述散射问题的Helmholtz方程转化为Lippmann–Schwinger（LS）积分方程，如果利用迭代法来求解LS方程则可以得到广泛使用的Born散射级数。然而，Born级数仅在小散射体或弱散射的情况下收敛；在强散射介质中，通常需要对散射级数进行某种重整化来得到收敛的Born级数［2123］。Wu等［16］利用De Wolf近似实现了地震波散射级数的重整化，同时利用单向波传播算子的屏近似算子进行双域计算，实现了中等速度扰动强度下的前向散射波场计算。Jakobsen［21］将量子物理学中的重整化理论引入到地震散射波场的正演模拟中，利用T矩阵（传输矩阵）传播算子重新表示LS方程并得到了关于T矩阵的迭代级数，从而改善了散射级数的收敛性。Jakobsen等［18］将T矩阵传播算子引入到格林函数表示的De Wolf散射级数中，给出了重整化后的De Wolf级数，并将利用重整化De Wolf级数计算得到的格林函数应用于强散射介质的正演模拟中。同时，Jakobsen等［24］将利用T矩阵重整化后的De Wolf级数应用于全波形反演中，并通过基于散射路经矩阵的域分解技术来加速反演计算，得到了较好的反演效果。为了进一步改善Born级数的收敛性以处理强散射问题，Jakobsen等［2526］基于重整化群理论推导出了重整化的散射级数解，该级数解在非均匀性较强的散射体介质中显示出了较好的收敛性；并利用重整化群和同伦延拓方法导出了LS方程的新散射级数解，该解在任意大的对比体积（扰动）下保证收敛。Huang等［19］基于多重散射理论、高斯波束和非线性Born近似，提出了用于地下速度重建的逆散射方法。该方法是一种新的Born迭代逆散射方法。Osnabrugge等［27］为改善强光学散射中Born级数的收敛性，引入了带虚部分量的复背景格林函数，得到了复波数LS方程的Born级数，并选取合适的对角预条件算子进一步加快Born散射级数的收敛。在Osnabrugge等的工作基础上，Huang等［28］直接利用复波数LS方程得到的收敛Born散射级数来处理地震散射问题，并从重整化角度重新阐释了Osnabrugge的收敛Born级数，通过数值算例验证了该方法在强扰动散射介质地震波场正演模拟中的适用性。徐杨杨等［2930］将解决强光学散射问题的广义超松弛迭代方法用于计算LS方程的数值解，通过分析该算法在反射地震条件下的适用性和收敛性得到了收敛的迭代级数解。但该方法在强速度扰动介质的情况下，当频率较高时，迭代级数收敛速度较慢，难以在有限的迭代次数内得到满足收敛误差的数值解。

为了改善广义超松弛迭代法求解LS方程时迭代级数高频收敛速度慢的问题，本文在背景波数中引入虚部分量并将带虚部分量的复波数背景格林函数引入到格林函数LS方程的广义逐次超松弛（generalized successive over-relaxation， GSOR）迭代算法中；同时引入预条件算子来降低离散系数矩阵的条件数，进一步改善数值迭代级数的收敛性。通过分析衰减因子和预条件算子对数值格林函数解的影响，选取合适的衰减因子和预条件算子，给出适应地震任意散射强度的收敛迭代级数解，并将该方法用于地震强散射问题中高频数值格林函数的计算。

1" 格林函数LS方程

1.1" 格林函数实波数LS方程

频率域中标量波动方程格林函数满足如下Helmholtz方程：

SymbolQC@2G（x，x′）+ω2v2（x）G（x，x′）=－δ（x－x′） 。（1）

式中：G（x，x′）为格林函数，表示单位点源x′到场点x的波场（地下空间点x′，x∈RD，RD为线性向量空间，D为所处理问题的维数）；ω为角频率；v（x）为地震波传播速度；δ（x－x′）为x′点处激发的点脉冲源。

将波速v（x）表示为背景速度v0（x）与一个扰动的叠加，并定义速度扰动O（x）=v20（x）/v2（x）－1，则式（1）可以写为

SymbolQC@2G（x，x′）+ω2v20（x）G（x，x′）=" －δ（x－x′）－ω2O（x）v20（x）G（x，x′） 。（2）

式（2）等号右端为等效源。背景格林函数G（0）（x，x′）满足

SymbolQC@2G（0）（x，x′）+k20G（0）（x，x′）=－δ（x－x′） 。（3）

式中，k0=ω/v0（x），为背景波数。当背景介质为均匀介质时，G（0）（x，x′）的解析表达式为

G（0）（x，x′）=exp（ik0x－x′）4πx－x′，" x′，x∈R3；i4H（1）0（k0x－x′），" x′，x∈R2；i2k0exp（ik0x－x′），" x′，x∈R1。 （4）

式中，H（1）0为第一类零阶Hankel函数。应用表示定理，可以得到由格林函数表示的LS方程：

G（x，x′）=G（0）（x，x′）+ ""k20∫ΩG（0）（x，x″）O（x″）G（x″，x′）dx″ 。（5）

式中，Ω为O（x″）≠0的散射区域（x″∈Ω）。

将线性积分算子与Dirac函数兼容，式（5）的连续矩阵乘积形式［18］为

G（x，x′）=G（0）（x，x′）+""" ∫Ω∫ΩG（0）（x，x″）V（x，x″）G（x″，x′）dxdx″ 。（6）

其中，

V（x，x″）=k20O（x）δ（x－x″） 。（7）

1.2 "格林函数复波数LS方程

对方程（2）左右两端同时加上iεG（x，x′）（衰减因子εgt;0），得到

SymbolQC@2G（x，x′）+ω2v20（x）G（x，x′）+iεG（x，x′）=－ω2O（x）v20（x）G（x，x′）－δ（x－x′）+iεG（x，x′） 。 （8）

定义复背景波数k^0=ω2/v20（x）+iε=k20+iε，相应地，复散射位O^（x）=ω2O（x）/v20（x）－iε。基于表示定理，可以得到格林函数表示的复波数LS方程：

G（x，x′）=G^（0）（x，x′）+" "∫ΩG^（0）（x，x″）O^（x″）G（x″，x′）dx″ 。" （9）

式中，G^（0）（x，x′）为复波数背景格林函数。复波数LS方程（式（9））与实波数LS方程（式（5））在数学上等价，但在数值分析上并不等效。G^（0）（x，x′）满足

SymbolQC@2G^（0）（x，x′）+k^20G^（0）（x，x′）=－δ（x－x′） 。（10）

相应地，在均匀背景介质中，G^（0）（x，x′）的解析表达式为

G^（0）（x，x′）=

exp（ik^0x－x′）4πx－x′，" x′，x∈R3；i4H（1）0k^0x－x′，" x′，x∈R2；i2k^0exp（ik^0x－x′），" x′，x∈R1。（11）

当k^0含有正的虚部分量时，Helmholtz方程的基本解格林函数G^（0）（x，x′）随距离呈指数衰减，使得背景格林函数的总能量有限，而虚部分量iε带来的能量损失由散射位O^（x）的虚部分量来补偿，这使得Helmholtz方程的解保持不变。复背景波数虚部分量iε的物理意义为：波在背景介质中传播是有能量损耗的，而这种能量损耗由散射位通过等量增益来补偿。

2" 复波数LS方程的预条件广义逐次超松弛迭代

根据格林函数表示的LS方程的矩阵连乘形式（式（6）），其算子形式［18］可以表示为

G=G（0）+G（0）VG。（12）

式中，G、G（0）、V分别为G（x，x′）、G（0）（x，x′）、V（x，x′）的矩阵算子表示。定义矩阵算子L=I－G（0）V（I为单位矩阵算子），并引入连续可逆算子C和有界线性算子B，式（12）的广义迭代格式［29， 31］为

Gn=C－1［（C－BL）Gn－1+BG（0）］，G0∈L2（Ω），（13）

或

Gn=（I－C－1BL）Gn－1+C－1BG（0）=

Gn－1+C－1B（G（0）－LGn－1）=

Gn－1+C－1Brn－1。（14）

式中：n为迭代次数；L2（Ω）为完备内积空间（Hilbert空间）；残差向量r0=G（0）－LG0，rn=G（0）－LGn。构造合适的算子C和B，使得谱半径ρ（I－C－1BL）lt;1，则LS方程的广义迭代格式（式（13）（14））收敛。这里线性可逆算子C相当于预条件算子［29］。

当C－1=I，B=αnI（αn为松弛因子）时，则得到GSOR迭代格式［29］。其中αn通过在迭代过程中最小化残差向量范数‖rn‖得到［29， 32］。考虑G0=G（0）时的第n次迭代结果，取

αn=∫rn－1Lrn－1dx∫Lrn－1Lrn－1dx=（rn－1，Lrn－1）‖Lrn－1‖2，（15）

使‖rn‖在Hilbert空间中极小。相应地，在背景波数中引入虚部分量iε后，复波数LS方程的算子形式可以表示为

G=G^（0）+G^（0）V^G。 （16）

式中，G^（0）、V^分别为G^（0）（x，x′）、V^（x，x″）=O^（x）δ（x－x″）的矩阵算子表示。复波数LS方程的广义迭代格式可以表示为

Gn=C－1［（C－BL^）Gn－1+BG^（0）］， G0∈L2（Ω），（17）

或

Gn=（I－C－1BL^）Gn－1+C－1BG^（0）=

Gn－1+C－1B（G^（0）－L^Gn－1）=

Gn－1+C－1Br^n－1。（18）

其中：

r^0=G^（0）－L^G0；

r^n=G^（0）－L^Gn；

L^=I－G^（0）V^。

对于复波数LS方程的广义迭代格式（式（18）），若C－1=γo=iεV^（γo为预条件矩阵算子），B=I，则可以得到Osnabrugge等［27］给出的强光学散射条件下的散射级数：

Gn=（I－γoL^）Gn－1+γoG^（0） =

Gn－1+γo（G^（0）－

L^Gn－1）=Gn－1+γor^n－1。 （19）

即

Gn=Gn－1－γo［Gn－1－（G^（0）V^Gn－1+G^（0））］。（20）

若C－1=γ（γ为任意预条件矩阵算子），B=α^nI，则可以得到预条件广义逐次超松弛（preconditioned generalized successive over-relaxation， Pre-GSOR）迭代格式：

Gn=（I－γα^nL^）Gn－1+γα^nG^（0）=

Gn－1+γα^n（G^（0）－L^Gn－1）=

Gn－1+γα^nr^n－1，（21）

或

Gn=Gn－1－γα^n［Gn－1－（G^（0）V^Gn－1+G^（0））］。（22）

此时，

α^n=∫r^n－1L^r^n－1dx∫L^r^n－1L^r^n－1dx=（r^n－1，L^r^n－1）‖L^r^n－1‖2。（23）

在迭代计算中，矩矢相乘的运算量用时比重最大。根据标量格林函数具有平移不变性及非方向性，由背景格林函数离散生成的矩阵为二重Toeplitz矩阵。根据系数矩阵的分块Toeplitz结构，只需要存第一列（第一行）元素即可生成整个矩阵，并且可以利用二维快速傅里叶变换（fast fourier transform， FFT）来加速迭代级数计算中的矩矢乘积［30］。

由于复波数背景格林函数随距离衰减带来的能量损失需由复散射位的虚部分量来补偿，才能使得复波数LS方程的解与原LS方程的解保持不变。因此，在地震勘探尺度和强速度扰动介质条件下，采用Pre-GSOR方法求解复波数LS方程计算数值格林函数时，有必要对衰减因子和预条件算子对数值格林函数解的影响进行详细分析，进而选取合适的衰减因子和预条件算子，以得到地震强散射条件下的数值格林函数解。

3" 衰减因子对数值格林函数解的影响

下面通过给定不同的衰减因子来计算2D复波数背景格林函数，并分析它对复波数LS方程数值解的影响。

根据式（11），设计v0=3000 m/s的均匀介质模型，模型大小为3.0 km×3.0 km，频率为20 Hz，源点位置x′=（1500 m，0 m），分别计算ε=0， 0.05时2D均匀介质中的复波数背景格林函数。此时，背景格林函数振幅随距离的衰减曲线如图1所示，2D均匀介质中的背景格林函数如图2所示。从图1和图2中可以看出，在地震勘探尺度上，背景速度和频率给定后（k0为常数），当ε=0.05时背景格林函数的能量就已经出现明显的衰减现象。

分析Osnabrugge等［27］给出的衰减因子选取原则：

a. 实部；b. 虚部。

a. 实部（ε=0）；b. 实部（ε=0.05）；c. 虚部（ε=0）；d. 虚部（ε=0.05）。

ε≥maxxk2（x）－k20=

maxxk20O（x）=k20Omax。（24）

式中，Omax为O（x）在所有空间点上的最大绝对值。取v0=1500 m/s、v=4482 m/s的强散射介质，模型大小为3.0 km×3.0 km，频率为20 Hz，

源点位置x′=（1500 m，0 m），分别计算ε=0.05k20Omax和ε=k20Omax时的复波数背景格林函数。

此时，背景格林函数振幅随距离的衰减曲线如图3所示，2D强散射介质中的背景格林函数如图4所示。

从图3和图4中可以看出，当ε=k20Omax时，其能量在源点附近迅速衰减。

a. 实部；b. 虚部。

a. 实部（ε=0.05k20|O|max）；b. 实部（ε=k20|O|max）；c. 虚部（ε=0.05k20|O|max）；d. 虚部（ε=k20|O|max）。

根据复波数LS方程（式（9））及其离散线性方程组（式（16））可以看出，背景格林函数中引入虚部分量iε后带来的能量损失需由复散射位中的虚部分量以算子相乘的形式来补偿。因此，复波数背景格林函数在模型空间有效计算区域内必须满足|G^（0）（x，x′）|mingt;Θ（Θ表示浮点数0，数据类型float型Θ=1.0×10－6，double型Θ=1.0×10－15）。如此，当背景格林函数与复散射位生成的矩阵算子相乘时在相应的空间计算点上矩阵元素才有效，由虚部分量iε带来的能量损失才能够得到补偿。对于地震强散射介质，在计算频率较高或模型尺度较大时，选取Osnabrugge等［27］给出的衰减因子，背景格林函数能量衰减过快，导致背景格林函数的能量损失无法得到补偿，复波数Helmholtz方程得到的格林函数数值解无法保证与原方程解相等。

本文在Osnabrugge等［27］研究的基础上，保留衰减因子与速度扰动之间的关系，通过引入常数a来给出地震勘探尺度下衰减因子的选取原则：

ε=ak20Omax（0≤alt;1.0） 。（25）

下文通过数值算例测试不同参数a的选取对迭代级数收敛性和解的影响，给出合适的参数a，使得在所有计算频率上，有效模型空间点上复波数背景格林函数G^（0）（x，x′）gt;Θ。

4" 对角预条件算子对系数矩阵的影响

对于强散射介质，采用迭代方法求解LS方程离散后的线性方程组时，随着频率变大，往往收敛速度会停滞不前甚至出现发散。广义超松弛迭代法虽然能保证迭代级数的收敛性，但在求解高频格林函数数值解时，仍然难以在有限的迭代次数内获得满足收敛误差的解［29］。

对于任意线性方程组Au=f，其数值解u～与精确解u的相对误差与残差向量r和系数矩阵条件数满足以下关系［33］：

‖u～－u‖‖u‖≤cond（A）‖r‖f；r=f－Au～。（26）

当利用迭代法求解线性方程组时，对于条件数比较大的系数矩阵，即便残差已经很小，也难以得到满足精度的数值解，因此，有必要通过选取合适的预条件算子，使方程组的系数矩阵具有更好的频谱（特征值）或较小的条件数，来改善广义超松弛迭代法收敛性。

下面通过分析Osnabrugge等［27］给出的预条件算子C－1=γo=iεV^对系数矩阵的作用，进一步给出符合地震勘探尺度下强散射介质的预条件算子选取原则。将复波数LS方程的算子形式（式（16））改写为线性方程：

L^G=G^（0）。 （27）

对式（27）采用左侧预条件算子［33］，则等效方程组表示为

C－1L^G=C－1（I－G^（0）V^）G=C－1G^（0）。 （28）

对Osnabrugge等［27］给出的预条件算子进行化简：

γo=iεV^=iε［k20O（x）－iε］=1+ik20O（x）ε。（29）

当ε=k20Omax时，预条件算子为

γo=1+iO（x）Omax。（30）

从式（28）可以看出，在系数矩阵L^中对背景格林函数离散矩阵G^（0）右乘复散射势对角算子V^，相当于对G^（0）的矩阵元素按列缩放。而预条件算子γo同样为对角矩阵算子，对L^使用左侧预条件算子γo，其作用是对L^进行按行缩放来平衡矩阵元素，达到改善系数矩阵条件数的目的［33］。根据Osnabrugge等［27］给出的预条件算子表达式，为了保证收敛性，本文在其基础上引入常数b，进一步给出地震尺度下预条件算子的表达式：

C－1=γ=1+iO（x）bOmax（b≥1.0）。（31）

下文通过数值算例测试不同参数b选取生成不同对角预条件算子对迭代级数收敛性的影响，针对不同的模型选取合适的参数b，并给出相应的预条件算子。

5" 数值实验

下面通过数值实验测试复波数LS方程Pre-GSOR迭代法计算强速度扰动模型数值格林函数的有效性，并分析其收敛特性和计算效率。LS方程的数值离散采用弱奇异核的Nystrm法，迭代求解过程中采用FFT加速空间褶积计算，通过并行算法获得所有频率下的数值格林函数，将得到的格林函数应用于地震散射波场的正演模拟中。基于以上思路，为了与精确数值解进行对比，并分析LS方程生成离散系数矩阵的数学性质，比较迭代法与直接法的数值结果，设计模型尺度较小的强速度扰动单体盐丘模型和重采样的SEG/EAGE（society of exploration geophysicists/ European association of geoscientists amp; engineers）盐丘模型，并主要测试Pre-GSOR迭代法计算中高频率格林函数时的稳定性和收敛性，最后分析快速算法在各种模型下存储和计算效率的优势。数值实验环境：Linux操作系统，四个物理CPU，每个CPU为16核，CPU型号为Intel（R） Xeon（R） Gold 6130 CPU @ 2.10 GHz。

5.1" 单体盐丘模型

将复波数广义超松弛迭代算法应用于强速度扰动的单体盐丘模型（图5）。背景速度v0=2000 m/s，盐丘速度v=4500 m/s，震源子波为主频15 Hz的雷克子波，震源位置x′=xs=（350 m，0 m），检波器置于z=10 m的界面上，空间采样间隔dx=dz=10 m，时间采样间隔为4 ms。

对实波数LS方程离散后的线性方程组采用直接法求解LG=G（0），获得精确的数值格林函数解。为了保证数值解的精确性，采用受系数矩阵特征值（或条件数）影响较小的奇异值分解（singular value decomposition，SVD）法来求解。复波数LS方程生成的线性方程组采用Pre-GSOR迭代法求解γL^G=γG^（0）。

表1给出了单体盐丘模型选取不同a和b、满足归一化残差r=‖G^（0）－L^Gn‖/‖G^（0）‖lt;1.0×10－6时所需要的迭代次数（设定最大迭代次数为1 000）。

令a=0.3，b=1.0，即ε=0.3k20Omax，γ=1+i［O（x）/Omax］，当频率为30、50 Hz时，数值格林函数解的实部和虚部如图6、图7所示。从图6、图7中可以看出，对背景波数引入虚部分量并采用Pre-GSOR迭代法来求解复波数LS方程，得到的格林函数数值结果与原LS方程的精确数值解（SVD）具有较好的一致性。

绘制格林函数数值解归一化残差与迭代次数的关系，对Pre-GSOR迭代法在中高频率的收敛性进行分析。图8为频率为30、50 Hz时，采用Pre-GSOR迭代法和GSOR迭代法得到的归一化残差与迭代次数的关系曲线。显然，Pre-GSOR迭代法的初始迭代误差较小，即便在高频时也表现出了较好的收敛性。

将复波数LS方程的Pre-GSOR格林函数数值解用于盐丘模型的散射波场计算，共71道，道间距10 m。图9、图10分别为30、50 Hz接收点处的单频散射格林函数Gsg。从图9、图10中可以看出，对于强散射介质，复波数LS方程Pre-GSOR迭代法不仅具有较好的稳定性，并且与原LS方程的精确数值解（SVD）匹配程度较高。

计算单频散射波场usg=s（ω）Gsg（s（ω）为震源子波）并将所有频率的散射波场变换到时间域，即可得到单炮记录。将数值结果与有限差分（finite difference， FD）法得到的数值结果进行对比，结果如图11所示。从图11中可以看出，复波数LS方程Pre-GSOR法可以得到与FD相近的数值结果。随机抽取第10道和第30道进行单道对比，结果见图12。从图12中可以看出，两种数值方法获得的单道信号相位和幅值均具有较好的一致性。

为进一步考察预条件算子的作用，计算系数矩阵2范数条件数（cond（L）2=‖L‖2‖L－1‖2=σmax（L）/σmin（L），σmax和σmin分别为系数矩阵L的最大和最小奇异值），并分析系数矩阵条件数随频率

的变化规律（表2）。从表中2中可以看出：随着频率的增加，原LS方程的系数矩阵L的条件数不断增大，对应的线性方程组变成病态方程组，这使得计算高频数值格林函数时，迭代法收敛性变差、收敛速度停滞不前，无法得到符合收敛误差的数值解；而对复波数LS方程系数矩阵L^使用左预条件算子γ后，得到的新系数矩阵γL^条件数较小，并随计算频率升高基本保持稳定，因此该方法可以有效改善迭代级数的收敛性。

5.2" SEG/EAGE盐丘模型

考虑更复杂的SEG/EAGE盐丘模型。为了与直接法（SVD）求解实波数LS方程的数值解进行对比，对SEG/EAGE盐丘速度模型进行重采样（图13），背景速度v0=1500 m，震源子波为主频15 Hz的雷克子波，震源位置x′=xs=（900 m，0 m），检波器置于z=0 m的界面上，迭代法空间采样间隔为dx=dz=10 m，时间采样间隔为4 ms。

同样，对实波数LS方程离散后的线性方程组采用SVD求解LG=G（0），对复波数LS方程生成的线性方程组采用Pre-GSOR迭代法求解γL^G=γG^（0），来计算数值格林函数。

表3给出了选取不同a和b、归一化残差满足rlt;1.0×10－3时所需要的迭代次数（设定最大迭代次数为3 000）。从表3中可以看出，当参数a越大，即衰减越强时，方程组的收敛性越好；但复波数系数矩阵L^和右端项G^（0）较原实波数的系数矩阵L和右端项G（0）产生的扰动也相应变大，当系数矩阵L的条件数过大时将会造成迭代数值解与精确数值解之间的误差变大。

由于高频实波数系数矩阵条件数过大（表4），等效方程组中系数矩阵和右端项的小扰动也会引起解的误差增大。为了保证收敛性的同时得到有效的数值解，取ε=0.03k20Omax，γ=1+i［O（x）/（8.0Omax）］，当频率为30、60 Hz时，SEG/EAGE盐丘模型数值格林函数解的实部和虚部如图14、图15所示。从图14、图15中可以看出，复波数LS方程Pre-GSOR迭代法得到的格林函数数值解与原LS方程的精确数值解（SVD）之间的误差较小。

选取特定频率，绘制格林函数数值解的归一化残差与迭代次数的关系图。图16为频率为30、60 Hz时，复波数LS方程的Pre-GSOR迭代法和实波数LS方程的GSOR迭代法得到的归一化残差与迭代次数的关系曲线。可以看出，相比实波数LS方程的GSOR迭代法，复波数LS方程的Pre-GSOR迭代法具有更快的收敛速度。

将复波数LS方程Pre-GSOR得到的格林函数数值解用于盐丘模型的散射波场计算，共181道，

道间距10 m。图17、图18分别为30、60 Hz时，实波数LS方程直接法（SVD）和复波数LS方程Pre-GSOR迭代法得到的检波点处的单频散射波场。从图17、图18中可以看出：当频率为30 Hz时，两种方法得到的数值结果基本一致，误差较小；而当频率为60 Hz时，复波数LS方程Pre-GSOR迭代法

得到的数值解可能需要更多的迭代次数才能与实波数LS方程的精确数值解完全匹配。对于强散射介质，复波数LS方程Pre-GSOR迭代法不仅具有较好的稳定性，并且在有限的迭代次数之后与原LS方程的精确数值解（SVD）匹配程度较高，可以在有限的迭代次数之内获得有效的数值格林函数。

将所有频率的散射波场变换到时间域，得到单炮记录，与FD法得到的数值模拟结果进行对比，结果如图19所示；随机抽取第10道和第30道进行单

道对比，结果如图20所示。从图19、图20中可以看出，两种数值方法获得的单道信号相位和振幅仍然具有较好的一致性。

从表4可以看出：随着频率的增加，系数矩阵L的条件数不断增大，在频率60 Hz时，其条件数已经超过了2.2×104，这使得计算高频数值格林函数时，迭代法超过一定迭代次数之后，收敛残差几乎不再减小，收敛速度停滞不前，无法在有限的迭代次数内得到有效的数值解；而对复波数系数矩阵L^使用左预条件矩阵γ后，得到的新系数矩阵γL^条件数较小，可以有效改善迭代级数的收敛性质。

6" 计算效率与内存分析

当地下介质比较复杂，如包含较大尺度的复杂散射体且速度扰动更强时，积分方程法应用于此类介质模型时需要在全空间中进行数值离散。格林函数离散后的矩阵为N=NxNz阶方阵，普通的微机难以获取相应的栈内存来存储系数矩阵，且计算效率不满足实际应用需求。通过分析背景格林函数离散矩阵的分块Toeplitz性质，可将其存储由O（N2）

降为O（N）。另外，利用分块Toeplitz矩阵与循环矩阵的关系，可以利用二维FFT加速迭代过程中的矩矢乘积运算，将矩矢乘积的计算复杂度由O（N2）降低为O（NlogN）。表5为单频数值格林函数计算的平均CPU时间与单节点内存消耗。可以看出，利用本文提出的快速算法来计算数值格林函数，计算过程中占用的内存空间较小且计算效率较高，为后续格林函数在散射波场的正演模拟和偏移成像中的应用提供了有效的快速数值算法。

7" 结论

1）本文提出一种预条件广义逐次超松弛迭代法求解复波数LS方程的格林函数数值模拟方法，有效改善了强散射介质下广义超松弛迭代级数的收敛性。

2）通过分析衰减因子和对角预条件算子格林函数数值解的影响，给出了反射地震条件下衰减因子和预条件算子的选取原则。

3）利用复背景格林函数离散矩阵的块Toeplitz性质，可将内存由O（N2）降为O（N），迭代计算中矩阵向量乘积的计算复杂度由由O（N2）降低为O（NlogN）。

4）将复波数LS方程的预条件广义逐次超松弛迭代法应用于强散射介质下高频数值格林函数的计算及散射波场正演模拟中，得到与实波数LS方程直接法相一致的格林函数数值解。

综上，本文实现了强散射介质格林函数和散射波场的快速数值模拟，为后续进行3D全波场数值模拟和反演成像方法的研究坚定了基础。

参考文献（References）：

［1］" 孙建国.高频渐近散射理论及其在地球物理场数值模拟与反演成像中的应用：研究历史与研究现状概述以及若干新进展［J］.吉林大学学报（地球科学版），2016，46（4）：12311259.

Sun Jianguo. High-Frequency Asymptotic Scattering Theories and Their Applications in Numerical Modeling and Imaging of Geophysical Fields： An Overview of the Research History and the State-of-the-Art， and Some New Developments［J］. Journal of Jilin University （Earth Science Edition）， 2016， 46（4）： 12311259.

［2］" Cerveny V. Seismic Ray Theory［M］. Cambridge： Cambridge University Press， 2001.

［3］" Bleistein N. Mathematical Methods for Wave Phenomena［M］. San Diego： Academic Press， 2012.

［4］" Popov M M， Semtchenok N M， Popov P M， et al. Depth Migration by the Gaussian Beam Summation Method［J］. Geophysics， 2010， 75（2）： S81S93.

［5］" 石秀林，孙建国，孙辉，等. 基于delta波包叠加的时间域深度偏移［J］. 地球物理学报，2016， 59（7）： 26412649.

Shi Xiulin， Sun Jianguo， Sun Hui， et al. The Time-Domain Depth Migration by the Summation of Delta Packets［J］. Chinese Journal of Geophysics， 2016， 59（7）： 26412649.

［6］" 孙建国.Kirchhoff型偏移理论的研究历史、研究现状与发展趋势展望：与光学绕射理论的类比、若干新结果、新认识以及若干有待于解决的问题［J］.吉林大学学报（地球科学版）， 2012， 42（5）：15211552.

Sun Jianguo. The History， the State of the Art and the Future Trend of the Research of Kirchhoff-Type Migration Theory： A Comparison with Optical Diffraction Theory， Some New Result and New Understanding， and Some Problems to Be Solved［J］. Journal of Jilin University （Earth Science Edition）， 2012， 42（5）： 15211552.

［7］" Huang X， Sun H， Sun J. Born Modeling for Heterogeneous Media Using the Gaussian Beam Summation Based Green’s Function［J］. Journal of Applied Geophysics， 2016， 131： 191201.

［8］" Huang X， Sun J， Sun Z. Local Algorithm for Computing Complex Travel Time Based on the Complex Eikonal Equation［J］. Physical Review E， 2016， 93（4）： 043307.

［9］" Huang X. Integral Equation Methods with Multiple Scattering and Gaussian Beams in Inhomogeneous Background Media for Solving Nonlinear Inverse Scattering Problems［J］. IEEE Transactions on Geoscience and Remote Sensing， 2020， 59（6）： 53455351.

［10］" Schuster G T. Reverse-Time Migration=Generalized Diffraction Stack Migration［C］// 72nd Annual International Meeting. Oklohoma： SEG， 2002： 32123216.

［11］" Schuster G T. Seismic Inversion［M］. Tulsa： Society of Exploration Geophysicists， 2017.

［12］" Zhang G. Generalized Diffraction-Stack Migration［D］. Utah： The University of Utah， 2012.

［13］" Zhang G， Dai W， Zhou M， et al. Generalized Diffraction-Stack Migration and Filtering of Coherent Noise［J］. Geophysical Prospecting， 2014， 62（3）： 427442.

［14］" 张志佳，孙章庆，王瑞湖，等. 复杂地表起伏多重变加密网格有限差分法波动方程数值模拟［J］. 吉林大学学报（地球科学版），2023，53（2）：598608.

Zhang Zhijia， Sun Zhangqing， Wang Ruihu， et al. Numerical Simulation of Wave Equation Using Finite Difference Method with Rugged Variable Refined Multigrid in Complex Surface［J］. Journal of Jilin University （Earth Science Edition）， 2023， 53 （2）： 598608.

［15］" Wu R S. Wave Propagation， Scattering and Imaging Using Dual-Domain One-Way and One-Return Propagators［J］. Pure and Applied Geophysics，2003， 160（3/4）： 509539.

［16］" Wu R S， Xie X B， Wu X Y. One-Way and One-Return Approximations （De Wolf Approximation） for Fast Elastic Wave Modeling in Complex Media［J］. Advances in Geophysics， 2007， 48： 265322.

［17］" Innanen K A. Born Series Forward Modelling of Seismic Primary and Multiple Reflections： An Inverse Scattering Shortcut［J］. Geophysical Journal International， 2009， 177（3）： 11971204.

［18］" Jakobsen M， Wu R S. Renormalized Scattering Series for Frequency-Domain Waveform Modelling of Strong Velocity Contrasts［J］. Geophysical Journal International， 2016， 206（2）： 880899.

［19］" Huang L， Fehler M， Zheng Y， et al. Seismic-Wave Scattering， Imaging， and Inversion［J］. Communications in Computational Physics， 2020， 28（1）： 140.

［20］" 张盼，韩立国，巩向博，等. 金属矿地震勘探方法技术研究进展［J］. 吉林大学学报（地球科学版），2023，53（6）：19691982.

Zhang Pan， Han Liguo， Gong Xiangbo， et al. Research Progress on Seismic Exploration Methods and Technologies for Metal Mines［J］. Journal of Jilin University （Earth Science Edition）， 2023， 53 （6）： 19691982.

［21］" Jakobsen M. T-Matrix Approach to Seismic Forward Modelling in the Acoustic Approximation［J］. Studia Geophysica et Geodaetica， 2012， 56（1）： 120.

［22］" Eftekhar R， Hu H， Zheng Y. Convergence Acceleration in Scattering Series and Seismic Waveform Inversion Using Nonlinear Shanks Transformation［J］. Geophysical Journal International， 2018， 214（3）： 17321743.

［23］" Huang K. A Critical History of Renormalization［J］. International Journal of Modern Physics A， 2013， 28（29）： 1330050.

［24］" Jakobsen M， Wu R S. Accelerating the T-Matrix Approach to Seismic Full-Waveform Inversion by Domain Decomposition［J］. Geophysical Prospecting， 2018， 66（6）： 10391059.

［25］" Jakobsen M， Wu R S， Huang X. Seismic Waveform Modeling in Strongly Scattering Media Using Renormalization Group Theory［C］//88nd Annual International Meeting. Oklohoma： SEG， 2018： 50075011.

［26］" Jakobsen M， Wu R S， Huang X. Convergent Scattering Series Solution of the Inhomogeneous Helmholtz Equation via Renormalization Group and Homotopy Continuation Approaches［J］. Journal of Computational Physics， 2020， 409： 109343.

［27］" Osnabrugge G， Leedumrongwatthanakun S， Vellekoop I M. A Convergent Born Series for Solving the Inhomogeneous Helmholtz Equation in Arbitrarily Large Media［J］. Journal of Computational Physics， 2016， 322： 113124.

［28］" Huang X， Jakobsen M， Wu R S. On the Applicability of a Renormalized Born Series for Seismic Wavefield Modelling in Strongly Scattering Media［J］. Journal of Geophysics and Engineering， 2020， 17（2）： 277299.

［29］" 徐杨杨，孙建国，商耀达，等. Lippmann-Schwinger积分方程广义超松弛迭代解法及其收敛特性［J］. 地球物理学报， 2021， 64（1）： 249262.

Xu Yangyang， Sun Jianguo， Shang Yaoda， et al. The Generalized Over-Relaxation Iterative Method for Lippmann-Schwinger Equation and Its Convergence［J］. Chinese Journal of Geophysics， 2021， 64（1）： 249262.

［30］" 徐杨杨，孙建国，商耀达.一种利用Nystrm离散与FFT快速褶积的散射地震波并行计算方法［J］. 地球物理学报， 2021， 64（8）： 28772887.

Xu Yangyang， Sun Jianguo， Shang Yaoda. A Parallel Computation Method for Scattered Seismic Waves Using Nystrm Discretization and FFT Fast Convolution［J］. Chinese Journal of Geophysics， 2021， 64（8）： 28772887.

［31］" Petryshyn W V. On A General Iterative Method for the Approximate Solution of Linear Operator Equations［J］. Mathematics of Computation， 1963， 17： 110.

［32］" Kleinman R E， van den Berg P M. Iterative Methods for Solving Integral Equations［J］. Radio Science， 1991， 26（1）： 175181.

［33］" Golub G H， van Loan C F. Matrix Computations ［M］. 4nd ed. Baltimore： Johns Hopkins University Press， 2013.