基于預條件廣義逐次超松弛迭代法的數值格林函數計算方法

摘要：為了改善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.