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

摘要：為了改善Born散射級數(shù)解決地震強散射問題時的收斂性，將帶有虛部分量的復(fù)波數(shù)格林函數(shù)引入到求解格林函數(shù)Lippmann–Schwinger（LS）積分方程數(shù)值解的廣義逐次超松弛迭代法中，弱化格林函數(shù)的奇異性。引入預(yù)條件算子降低系數(shù)矩陣的條件數(shù)，加速迭代級數(shù)的收斂速度，給出了復(fù)波數(shù)LS方程的預(yù)條件廣義逐次超松弛（preconditioned generalized successive over-relaxation， Pre-GSOR）迭代格式。通過數(shù)值分析和收斂性分析重新選取合適的衰減因子和預(yù)條件算子，得到了滿足地震強散射條件的收斂Born級數(shù)，并將其用于地震強散射問題中數(shù)值格林函數(shù)的計算。數(shù)值結(jié)果表明：復(fù)波數(shù)LS方程Pre-GSOR迭代法可以得到與實波數(shù)LS方程直接法相匹配的數(shù)值模擬結(jié)果；復(fù)波數(shù)LS方程Pre-GSOR迭代法系數(shù)矩陣條件數(shù)在高頻時僅為原系數(shù)矩陣條件數(shù)的10%，相同迭代次數(shù)下歸一化收斂殘差可降低3個數(shù)量級以上，且對高頻適應(yīng)性強，可有效改善實波數(shù)LS方程廣義超松弛迭代法在強散射介質(zhì)中的收斂停滯問題。

關(guān)鍵詞：地震散射波場；格林函數(shù)；廣義超松弛迭代；預(yù)條件算子

doi：10.13278/j.cnki.jjuese.20230249

中圖分類號：P631.4

文獻標(biāo)志碼：A

徐楊楊，商耀達，孫建國. 基于預(yù)條件廣義逐次超松弛迭代法的數(shù)值格林函數(shù)計算方法. 吉林大學(xué)學(xué)報（地球科學(xué)版），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—），女，博士研究生，主要從事地震散射波場數(shù)值模擬、成像技術(shù)與快速算法研究，E-mail：846208564@qq.com

通信作者：孫建國（1956—），男，教授，博士生導(dǎo)師，主要從事地下波動理論與成像技術(shù)、計算地球物理、科學(xué)計算方法與技術(shù)、海洋反射地震資料處理、地下天線理論以及巖石物理等方面的教學(xué)和研究工作，E-mail：sun_jg@jlu.edu.cn

基金項目：國家自然科學(xué)基金項目（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" 引言

格林函數(shù)的表示和計算通常是實現(xiàn)基于Kirchhoff型、Born型等積分方程地震正反演與偏移成像方法的核心，格林函數(shù)表示方法不同則可得到不同的正反演與成像方法［1］。常用的計算格林函數(shù)的方法主要有基于漸近射線理論的方法［29］、基于微分算子的純數(shù)值方法［1014］和基于散射理論的方法［1520］等。漸近射線理論以波動方程的高頻近似假設(shè)為前提，利用零階幾何射線理論或高斯波束疊加等波場高頻表示來近似格林函數(shù)，即高頻漸近格林函數(shù)［6］。對于與地震波長相比較大的光滑模型，漸近射線理論非常有效，具有運算效率高、可控性好的優(yōu)點，但漸近射線理論認為只有一次反射才是有效信號，無法利用多次反射、繞射及散射等波動現(xiàn)象所攜帶的有效信息［1］。而純數(shù)值方法和散射理論均可以用于計算包含完整波形的格林函數(shù)。相比于純數(shù)值方法，散射理論具有更高的計算效率［18］。

基于散射理論表示格林函數(shù)時，通常將描述散射問題的Helmholtz方程轉(zhuǎn)化為Lippmann–Schwinger（LS）積分方程，如果利用迭代法來求解LS方程則可以得到廣泛使用的Born散射級數(shù)。然而，Born級數(shù)僅在小散射體或弱散射的情況下收斂；在強散射介質(zhì)中，通常需要對散射級數(shù)進行某種重整化來得到收斂的Born級數(shù)［2123］。Wu等［16］利用De Wolf近似實現(xiàn)了地震波散射級數(shù)的重整化，同時利用單向波傳播算子的屏近似算子進行雙域計算，實現(xiàn)了中等速度擾動強度下的前向散射波場計算。Jakobsen［21］將量子物理學(xué)中的重整化理論引入到地震散射波場的正演模擬中，利用T矩陣（傳輸矩陣）傳播算子重新表示LS方程并得到了關(guān)于T矩陣的迭代級數(shù)，從而改善了散射級數(shù)的收斂性。Jakobsen等［18］將T矩陣傳播算子引入到格林函數(shù)表示的De Wolf散射級數(shù)中，給出了重整化后的De Wolf級數(shù)，并將利用重整化De Wolf級數(shù)計算得到的格林函數(shù)應(yīng)用于強散射介質(zhì)的正演模擬中。同時，Jakobsen等［24］將利用T矩陣重整化后的De Wolf級數(shù)應(yīng)用于全波形反演中，并通過基于散射路經(jīng)矩陣的域分解技術(shù)來加速反演計算，得到了較好的反演效果。為了進一步改善Born級數(shù)的收斂性以處理強散射問題，Jakobsen等［2526］基于重整化群理論推導(dǎo)出了重整化的散射級數(shù)解，該級數(shù)解在非均勻性較強的散射體介質(zhì)中顯示出了較好的收斂性；并利用重整化群和同倫延拓方法導(dǎo)出了LS方程的新散射級數(shù)解，該解在任意大的對比體積（擾動）下保證收斂。Huang等［19］基于多重散射理論、高斯波束和非線性Born近似，提出了用于地下速度重建的逆散射方法。該方法是一種新的Born迭代逆散射方法。Osnabrugge等［27］為改善強光學(xué)散射中Born級數(shù)的收斂性，引入了帶虛部分量的復(fù)背景格林函數(shù)，得到了復(fù)波數(shù)LS方程的Born級數(shù)，并選取合適的對角預(yù)條件算子進一步加快Born散射級數(shù)的收斂。在Osnabrugge等的工作基礎(chǔ)上，Huang等［28］直接利用復(fù)波數(shù)LS方程得到的收斂Born散射級數(shù)來處理地震散射問題，并從重整化角度重新闡釋了Osnabrugge的收斂Born級數(shù)，通過數(shù)值算例驗證了該方法在強擾動散射介質(zhì)地震波場正演模擬中的適用性。徐楊楊等［2930］將解決強光學(xué)散射問題的廣義超松弛迭代方法用于計算LS方程的數(shù)值解，通過分析該算法在反射地震條件下的適用性和收斂性得到了收斂的迭代級數(shù)解。但該方法在強速度擾動介質(zhì)的情況下，當(dāng)頻率較高時，迭代級數(shù)收斂速度較慢，難以在有限的迭代次數(shù)內(nèi)得到滿足收斂誤差的數(shù)值解。

為了改善廣義超松弛迭代法求解LS方程時迭代級數(shù)高頻收斂速度慢的問題，本文在背景波數(shù)中引入虛部分量并將帶虛部分量的復(fù)波數(shù)背景格林函數(shù)引入到格林函數(shù)LS方程的廣義逐次超松弛（generalized successive over-relaxation， GSOR）迭代算法中；同時引入預(yù)條件算子來降低離散系數(shù)矩陣的條件數(shù)，進一步改善數(shù)值迭代級數(shù)的收斂性。通過分析衰減因子和預(yù)條件算子對數(shù)值格林函數(shù)解的影響，選取合適的衰減因子和預(yù)條件算子，給出適應(yīng)地震任意散射強度的收斂迭代級數(shù)解，并將該方法用于地震強散射問題中高頻數(shù)值格林函數(shù)的計算。

1" 格林函數(shù)LS方程

1.1" 格林函數(shù)實波數(shù)LS方程

頻率域中標(biāo)量波動方程格林函數(shù)滿足如下Helmholtz方程：

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

式中：G（x，x′）為格林函數(shù)，表示單位點源x′到場點x的波場（地下空間點x′，x∈RD，RD為線性向量空間，D為所處理問題的維數(shù)）；ω為角頻率；v（x）為地震波傳播速度；δ（x－x′）為x′點處激發(fā)的點脈沖源。

將波速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）等號右端為等效源。背景格林函數(shù)G（0）（x，x′）滿足

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

式中，k0=ω/v0（x），為背景波數(shù)。當(dāng)背景介質(zhì)為均勻介質(zhì)時，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函數(shù)。應(yīng)用表示定理，可以得到由格林函數(shù)表示的LS方程：

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

式中，Ω為O（x″）≠0的散射區(qū)域（x″∈Ω）。

將線性積分算子與Dirac函數(shù)兼容，式（5）的連續(xù)矩陣乘積形式［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 "格林函數(shù)復(fù)波數(shù)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）

定義復(fù)背景波數(shù)k^0=ω2/v20（x）+iε=k20+iε，相應(yīng)地，復(fù)散射位O^（x）=ω2O（x）/v20（x）－iε?；诒硎径ɡ?，可以得到格林函數(shù)表示的復(fù)波數(shù)LS方程：

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

式中，G^（0）（x，x′）為復(fù)波數(shù)背景格林函數(shù)。復(fù)波數(shù)LS方程（式（9））與實波數(shù)LS方程（式（5））在數(shù)學(xué)上等價，但在數(shù)值分析上并不等效。G^（0）（x，x′）滿足

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

相應(yīng)地，在均勻背景介質(zhì)中，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）

當(dāng)k^0含有正的虛部分量時，Helmholtz方程的基本解格林函數(shù)G^（0）（x，x′）隨距離呈指數(shù)衰減，使得背景格林函數(shù)的總能量有限，而虛部分量iε帶來的能量損失由散射位O^（x）的虛部分量來補償，這使得Helmholtz方程的解保持不變。復(fù)背景波數(shù)虛部分量iε的物理意義為：波在背景介質(zhì)中傳播是有能量損耗的，而這種能量損耗由散射位通過等量增益來補償。

2" 復(fù)波數(shù)LS方程的預(yù)條件廣義逐次超松弛迭代

根據(jù)格林函數(shù)表示的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為單位矩陣算子），并引入連續(xù)可逆算子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為迭代次數(shù)；L2（Ω）為完備內(nèi)積空間（Hilbert空間）；殘差向量r0=G（0）－LG0，rn=G（0）－LGn。構(gòu)造合適的算子C和B，使得譜半徑ρ（I－C－1BL）lt;1，則LS方程的廣義迭代格式（式（13）（14））收斂。這里線性可逆算子C相當(dāng)于預(yù)條件算子［29］。

當(dāng)C－1=I，B=αnI（αn為松弛因子）時，則得到GSOR迭代格式［29］。其中αn通過在迭代過程中最小化殘差向量范數(shù)‖rn‖得到［29， 32］?？紤]G0=G（0）時的第n次迭代結(jié)果，取

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

使‖rn‖在Hilbert空間中極小。相應(yīng)地，在背景波數(shù)中引入虛部分量iε后，復(fù)波數(shù)LS方程的算子形式可以表示為

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

式中，G^（0）、V^分別為G^（0）（x，x′）、V^（x，x″）=O^（x）δ（x－x″）的矩陣算子表示。復(fù)波數(shù)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^。

對于復(fù)波數(shù)LS方程的廣義迭代格式（式（18）），若C－1=γo=iεV^（γo為預(yù)條件矩陣算子），B=I，則可以得到Osnabrugge等［27］給出的強光學(xué)散射條件下的散射級數(shù)：

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=γ（γ為任意預(yù)條件矩陣算子），B=α^nI，則可以得到預(yù)條件廣義逐次超松弛（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）

在迭代計算中，矩矢相乘的運算量用時比重最大。根據(jù)標(biāo)量格林函數(shù)具有平移不變性及非方向性，由背景格林函數(shù)離散生成的矩陣為二重Toeplitz矩陣。根據(jù)系數(shù)矩陣的分塊Toeplitz結(jié)構(gòu)，只需要存第一列（第一行）元素即可生成整個矩陣，并且可以利用二維快速傅里葉變換（fast fourier transform， FFT）來加速迭代級數(shù)計算中的矩矢乘積［30］。

由于復(fù)波數(shù)背景格林函數(shù)隨距離衰減帶來的能量損失需由復(fù)散射位的虛部分量來補償，才能使得復(fù)波數(shù)LS方程的解與原LS方程的解保持不變。因此，在地震勘探尺度和強速度擾動介質(zhì)條件下，采用Pre-GSOR方法求解復(fù)波數(shù)LS方程計算數(shù)值格林函數(shù)時，有必要對衰減因子和預(yù)條件算子對數(shù)值格林函數(shù)解的影響進行詳細分析，進而選取合適的衰減因子和預(yù)條件算子，以得到地震強散射條件下的數(shù)值格林函數(shù)解。

3" 衰減因子對數(shù)值格林函數(shù)解的影響

下面通過給定不同的衰減因子來計算2D復(fù)波數(shù)背景格林函數(shù)，并分析它對復(fù)波數(shù)LS方程數(shù)值解的影響。

根據(jù)式（11），設(shè)計v0=3000 m/s的均勻介質(zhì)模型，模型大小為3.0 km×3.0 km，頻率為20 Hz，源點位置x′=（1500 m，0 m），分別計算ε=0， 0.05時2D均勻介質(zhì)中的復(fù)波數(shù)背景格林函數(shù)。此時，背景格林函數(shù)振幅隨距離的衰減曲線如圖1所示，2D均勻介質(zhì)中的背景格林函數(shù)如圖2所示。從圖1和圖2中可以看出，在地震勘探尺度上，背景速度和頻率給定后（k0為常數(shù)），當(dāng)ε=0.05時背景格林函數(shù)的能量就已經(jīng)出現(xiàn)明顯的衰減現(xiàn)象。

分析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的強散射介質(zhì)，模型大小為3.0 km×3.0 km，頻率為20 Hz，

源點位置x′=（1500 m，0 m），分別計算ε=0.05k20Omax和ε=k20Omax時的復(fù)波數(shù)背景格林函數(shù)。

此時，背景格林函數(shù)振幅隨距離的衰減曲線如圖3所示，2D強散射介質(zhì)中的背景格林函數(shù)如圖4所示。

從圖3和圖4中可以看出，當(dāng)ε=k20Omax時，其能量在源點附近迅速衰減。

a. 實部；b. 虛部。

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

根據(jù)復(fù)波數(shù)LS方程（式（9））及其離散線性方程組（式（16））可以看出，背景格林函數(shù)中引入虛部分量iε后帶來的能量損失需由復(fù)散射位中的虛部分量以算子相乘的形式來補償。因此，復(fù)波數(shù)背景格林函數(shù)在模型空間有效計算區(qū)域內(nèi)必須滿足|G^（0）（x，x′）|mingt;Θ（Θ表示浮點數(shù)0，數(shù)據(jù)類型float型Θ=1.0×10－6，double型Θ=1.0×10－15）。如此，當(dāng)背景格林函數(shù)與復(fù)散射位生成的矩陣算子相乘時在相應(yīng)的空間計算點上矩陣元素才有效，由虛部分量iε帶來的能量損失才能夠得到補償。對于地震強散射介質(zhì)，在計算頻率較高或模型尺度較大時，選取Osnabrugge等［27］給出的衰減因子，背景格林函數(shù)能量衰減過快，導(dǎo)致背景格林函數(shù)的能量損失無法得到補償，復(fù)波數(shù)Helmholtz方程得到的格林函數(shù)數(shù)值解無法保證與原方程解相等。

本文在Osnabrugge等［27］研究的基礎(chǔ)上，保留衰減因子與速度擾動之間的關(guān)系，通過引入常數(shù)a來給出地震勘探尺度下衰減因子的選取原則：

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

下文通過數(shù)值算例測試不同參數(shù)a的選取對迭代級數(shù)收斂性和解的影響，給出合適的參數(shù)a，使得在所有計算頻率上，有效模型空間點上復(fù)波數(shù)背景格林函數(shù)G^（0）（x，x′）gt;Θ。

4" 對角預(yù)條件算子對系數(shù)矩陣的影響

對于強散射介質(zhì)，采用迭代方法求解LS方程離散后的線性方程組時，隨著頻率變大，往往收斂速度會停滯不前甚至出現(xiàn)發(fā)散。廣義超松弛迭代法雖然能保證迭代級數(shù)的收斂性，但在求解高頻格林函數(shù)數(shù)值解時，仍然難以在有限的迭代次數(shù)內(nèi)獲得滿足收斂誤差的解［29］。

對于任意線性方程組Au=f，其數(shù)值解u～與精確解u的相對誤差與殘差向量r和系數(shù)矩陣條件數(shù)滿足以下關(guān)系［33］：

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

當(dāng)利用迭代法求解線性方程組時，對于條件數(shù)比較大的系數(shù)矩陣，即便殘差已經(jīng)很小，也難以得到滿足精度的數(shù)值解，因此，有必要通過選取合適的預(yù)條件算子，使方程組的系數(shù)矩陣具有更好的頻譜（特征值）或較小的條件數(shù)，來改善廣義超松弛迭代法收斂性。

下面通過分析Osnabrugge等［27］給出的預(yù)條件算子C－1=γo=iεV^對系數(shù)矩陣的作用，進一步給出符合地震勘探尺度下強散射介質(zhì)的預(yù)條件算子選取原則。將復(fù)波數(shù)LS方程的算子形式（式（16））改寫為線性方程：

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

對式（27）采用左側(cè)預(yù)條件算子［33］，則等效方程組表示為

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

對Osnabrugge等［27］給出的預(yù)條件算子進行化簡：

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

當(dāng)ε=k20Omax時，預(yù)條件算子為

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

從式（28）可以看出，在系數(shù)矩陣L^中對背景格林函數(shù)離散矩陣G^（0）右乘復(fù)散射勢對角算子V^，相當(dāng)于對G^（0）的矩陣元素按列縮放。而預(yù)條件算子γo同樣為對角矩陣算子，對L^使用左側(cè)預(yù)條件算子γo，其作用是對L^進行按行縮放來平衡矩陣元素，達到改善系數(shù)矩陣條件數(shù)的目的［33］。根據(jù)Osnabrugge等［27］給出的預(yù)條件算子表達式，為了保證收斂性，本文在其基礎(chǔ)上引入常數(shù)b，進一步給出地震尺度下預(yù)條件算子的表達式：

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

下文通過數(shù)值算例測試不同參數(shù)b選取生成不同對角預(yù)條件算子對迭代級數(shù)收斂性的影響，針對不同的模型選取合適的參數(shù)b，并給出相應(yīng)的預(yù)條件算子。

5" 數(shù)值實驗

下面通過數(shù)值實驗測試復(fù)波數(shù)LS方程Pre-GSOR迭代法計算強速度擾動模型數(shù)值格林函數(shù)的有效性，并分析其收斂特性和計算效率。LS方程的數(shù)值離散采用弱奇異核的Nystrm法，迭代求解過程中采用FFT加速空間褶積計算，通過并行算法獲得所有頻率下的數(shù)值格林函數(shù)，將得到的格林函數(shù)應(yīng)用于地震散射波場的正演模擬中?；谝陨纤悸?，為了與精確數(shù)值解進行對比，并分析LS方程生成離散系數(shù)矩陣的數(shù)學(xué)性質(zhì)，比較迭代法與直接法的數(shù)值結(jié)果，設(shè)計模型尺度較小的強速度擾動單體鹽丘模型和重采樣的SEG/EAGE（society of exploration geophysicists/ European association of geoscientists amp; engineers）鹽丘模型，并主要測試Pre-GSOR迭代法計算中高頻率格林函數(shù)時的穩(wěn)定性和收斂性，最后分析快速算法在各種模型下存儲和計算效率的優(yōu)勢。數(shù)值實驗環(huán)境：Linux操作系統(tǒng)，四個物理CPU，每個CPU為16核，CPU型號為Intel（R） Xeon（R） Gold 6130 CPU @ 2.10 GHz。

5.1" 單體鹽丘模型

將復(fù)波數(shù)廣義超松弛迭代算法應(yīng)用于強速度擾動的單體鹽丘模型（圖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。

對實波數(shù)LS方程離散后的線性方程組采用直接法求解LG=G（0），獲得精確的數(shù)值格林函數(shù)解。為了保證數(shù)值解的精確性，采用受系數(shù)矩陣特征值（或條件數(shù)）影響較小的奇異值分解（singular value decomposition，SVD）法來求解。復(fù)波數(shù)LS方程生成的線性方程組采用Pre-GSOR迭代法求解γL^G=γG^（0）。

表1給出了單體鹽丘模型選取不同a和b、滿足歸一化殘差r=‖G^（0）－L^Gn‖/‖G^（0）‖lt;1.0×10－6時所需要的迭代次數(shù)（設(shè)定最大迭代次數(shù)為1 000）。

令a=0.3，b=1.0，即ε=0.3k20Omax，γ=1+i［O（x）/Omax］，當(dāng)頻率為30、50 Hz時，數(shù)值格林函數(shù)解的實部和虛部如圖6、圖7所示。從圖6、圖7中可以看出，對背景波數(shù)引入虛部分量并采用Pre-GSOR迭代法來求解復(fù)波數(shù)LS方程，得到的格林函數(shù)數(shù)值結(jié)果與原LS方程的精確數(shù)值解（SVD）具有較好的一致性。

繪制格林函數(shù)數(shù)值解歸一化殘差與迭代次數(shù)的關(guān)系，對Pre-GSOR迭代法在中高頻率的收斂性進行分析。圖8為頻率為30、50 Hz時，采用Pre-GSOR迭代法和GSOR迭代法得到的歸一化殘差與迭代次數(shù)的關(guān)系曲線。顯然，Pre-GSOR迭代法的初始迭代誤差較小，即便在高頻時也表現(xiàn)出了較好的收斂性。

將復(fù)波數(shù)LS方程的Pre-GSOR格林函數(shù)數(shù)值解用于鹽丘模型的散射波場計算，共71道，道間距10 m。圖9、圖10分別為30、50 Hz接收點處的單頻散射格林函數(shù)Gsg。從圖9、圖10中可以看出，對于強散射介質(zhì)，復(fù)波數(shù)LS方程Pre-GSOR迭代法不僅具有較好的穩(wěn)定性，并且與原LS方程的精確數(shù)值解（SVD）匹配程度較高。

計算單頻散射波場usg=s（ω）Gsg（s（ω）為震源子波）并將所有頻率的散射波場變換到時間域，即可得到單炮記錄。將數(shù)值結(jié)果與有限差分（finite difference， FD）法得到的數(shù)值結(jié)果進行對比，結(jié)果如圖11所示。從圖11中可以看出，復(fù)波數(shù)LS方程Pre-GSOR法可以得到與FD相近的數(shù)值結(jié)果。隨機抽取第10道和第30道進行單道對比，結(jié)果見圖12。從圖12中可以看出，兩種數(shù)值方法獲得的單道信號相位和幅值均具有較好的一致性。

為進一步考察預(yù)條件算子的作用，計算系數(shù)矩陣2范數(shù)條件數(shù)（cond（L）2=‖L‖2‖L－1‖2=σmax（L）/σmin（L），σmax和σmin分別為系數(shù)矩陣L的最大和最小奇異值），并分析系數(shù)矩陣條件數(shù)隨頻率

的變化規(guī)律（表2）。從表中2中可以看出：隨著頻率的增加，原LS方程的系數(shù)矩陣L的條件數(shù)不斷增大，對應(yīng)的線性方程組變成病態(tài)方程組，這使得計算高頻數(shù)值格林函數(shù)時，迭代法收斂性變差、收斂速度停滯不前，無法得到符合收斂誤差的數(shù)值解；而對復(fù)波數(shù)LS方程系數(shù)矩陣L^使用左預(yù)條件算子γ后，得到的新系數(shù)矩陣γL^條件數(shù)較小，并隨計算頻率升高基本保持穩(wěn)定，因此該方法可以有效改善迭代級數(shù)的收斂性。

5.2" SEG/EAGE鹽丘模型

考慮更復(fù)雜的SEG/EAGE鹽丘模型。為了與直接法（SVD）求解實波數(shù)LS方程的數(shù)值解進行對比，對SEG/EAGE鹽丘速度模型進行重采樣（圖13），背景速度v0=1500 m，震源子波為主頻15 Hz的雷克子波，震源位置x′=xs=（900 m，0 m），檢波器置于z=0 m的界面上，迭代法空間采樣間隔為dx=dz=10 m，時間采樣間隔為4 ms。

同樣，對實波數(shù)LS方程離散后的線性方程組采用SVD求解LG=G（0），對復(fù)波數(shù)LS方程生成的線性方程組采用Pre-GSOR迭代法求解γL^G=γG^（0），來計算數(shù)值格林函數(shù)。

表3給出了選取不同a和b、歸一化殘差滿足rlt;1.0×10－3時所需要的迭代次數(shù)（設(shè)定最大迭代次數(shù)為3 000）。從表3中可以看出，當(dāng)參數(shù)a越大，即衰減越強時，方程組的收斂性越好；但復(fù)波數(shù)系數(shù)矩陣L^和右端項G^（0）較原實波數(shù)的系數(shù)矩陣L和右端項G（0）產(chǎn)生的擾動也相應(yīng)變大，當(dāng)系數(shù)矩陣L的條件數(shù)過大時將會造成迭代數(shù)值解與精確數(shù)值解之間的誤差變大。

由于高頻實波數(shù)系數(shù)矩陣條件數(shù)過大（表4），等效方程組中系數(shù)矩陣和右端項的小擾動也會引起解的誤差增大。為了保證收斂性的同時得到有效的數(shù)值解，取ε=0.03k20Omax，γ=1+i［O（x）/（8.0Omax）］，當(dāng)頻率為30、60 Hz時，SEG/EAGE鹽丘模型數(shù)值格林函數(shù)解的實部和虛部如圖14、圖15所示。從圖14、圖15中可以看出，復(fù)波數(shù)LS方程Pre-GSOR迭代法得到的格林函數(shù)數(shù)值解與原LS方程的精確數(shù)值解（SVD）之間的誤差較小。

選取特定頻率，繪制格林函數(shù)數(shù)值解的歸一化殘差與迭代次數(shù)的關(guān)系圖。圖16為頻率為30、60 Hz時，復(fù)波數(shù)LS方程的Pre-GSOR迭代法和實波數(shù)LS方程的GSOR迭代法得到的歸一化殘差與迭代次數(shù)的關(guān)系曲線?？梢钥闯觯啾葘嵅〝?shù)LS方程的GSOR迭代法，復(fù)波數(shù)LS方程的Pre-GSOR迭代法具有更快的收斂速度。

將復(fù)波數(shù)LS方程Pre-GSOR得到的格林函數(shù)數(shù)值解用于鹽丘模型的散射波場計算，共181道，

道間距10 m。圖17、圖18分別為30、60 Hz時，實波數(shù)LS方程直接法（SVD）和復(fù)波數(shù)LS方程Pre-GSOR迭代法得到的檢波點處的單頻散射波場。從圖17、圖18中可以看出：當(dāng)頻率為30 Hz時，兩種方法得到的數(shù)值結(jié)果基本一致，誤差較?。欢?dāng)頻率為60 Hz時，復(fù)波數(shù)LS方程Pre-GSOR迭代法

得到的數(shù)值解可能需要更多的迭代次數(shù)才能與實波數(shù)LS方程的精確數(shù)值解完全匹配。對于強散射介質(zhì)，復(fù)波數(shù)LS方程Pre-GSOR迭代法不僅具有較好的穩(wěn)定性，并且在有限的迭代次數(shù)之后與原LS方程的精確數(shù)值解（SVD）匹配程度較高，可以在有限的迭代次數(shù)之內(nèi)獲得有效的數(shù)值格林函數(shù)。

將所有頻率的散射波場變換到時間域，得到單炮記錄，與FD法得到的數(shù)值模擬結(jié)果進行對比，結(jié)果如圖19所示；隨機抽取第10道和第30道進行單

道對比，結(jié)果如圖20所示。從圖19、圖20中可以看出，兩種數(shù)值方法獲得的單道信號相位和振幅仍然具有較好的一致性。

從表4可以看出：隨著頻率的增加，系數(shù)矩陣L的條件數(shù)不斷增大，在頻率60 Hz時，其條件數(shù)已經(jīng)超過了2.2×104，這使得計算高頻數(shù)值格林函數(shù)時，迭代法超過一定迭代次數(shù)之后，收斂殘差幾乎不再減小，收斂速度停滯不前，無法在有限的迭代次數(shù)內(nèi)得到有效的數(shù)值解；而對復(fù)波數(shù)系數(shù)矩陣L^使用左預(yù)條件矩陣γ后，得到的新系數(shù)矩陣γL^條件數(shù)較小，可以有效改善迭代級數(shù)的收斂性質(zhì)。

6" 計算效率與內(nèi)存分析

當(dāng)?shù)叵陆橘|(zhì)比較復(fù)雜，如包含較大尺度的復(fù)雜散射體且速度擾動更強時，積分方程法應(yīng)用于此類介質(zhì)模型時需要在全空間中進行數(shù)值離散。格林函數(shù)離散后的矩陣為N=NxNz階方陣，普通的微機難以獲取相應(yīng)的棧內(nèi)存來存儲系數(shù)矩陣，且計算效率不滿足實際應(yīng)用需求。通過分析背景格林函數(shù)離散矩陣的分塊Toeplitz性質(zhì)，可將其存儲由O（N2）

降為O（N）。另外，利用分塊Toeplitz矩陣與循環(huán)矩陣的關(guān)系，可以利用二維FFT加速迭代過程中的矩矢乘積運算，將矩矢乘積的計算復(fù)雜度由O（N2）降低為O（NlogN）。表5為單頻數(shù)值格林函數(shù)計算的平均CPU時間與單節(jié)點內(nèi)存消耗?？梢钥闯觯帽疚奶岢龅目焖偎惴▉碛嬎銛?shù)值格林函數(shù)，計算過程中占用的內(nèi)存空間較小且計算效率較高，為后續(xù)格林函數(shù)在散射波場的正演模擬和偏移成像中的應(yīng)用提供了有效的快速數(shù)值算法。

7" 結(jié)論

1）本文提出一種預(yù)條件廣義逐次超松弛迭代法求解復(fù)波數(shù)LS方程的格林函數(shù)數(shù)值模擬方法，有效改善了強散射介質(zhì)下廣義超松弛迭代級數(shù)的收斂性。

2）通過分析衰減因子和對角預(yù)條件算子格林函數(shù)數(shù)值解的影響，給出了反射地震條件下衰減因子和預(yù)條件算子的選取原則。

3）利用復(fù)背景格林函數(shù)離散矩陣的塊Toeplitz性質(zhì)，可將內(nèi)存由O（N2）降為O（N），迭代計算中矩陣向量乘積的計算復(fù)雜度由由O（N2）降低為O（NlogN）。

4）將復(fù)波數(shù)LS方程的預(yù)條件廣義逐次超松弛迭代法應(yīng)用于強散射介質(zhì)下高頻數(shù)值格林函數(shù)的計算及散射波場正演模擬中，得到與實波數(shù)LS方程直接法相一致的格林函數(shù)數(shù)值解。

綜上，本文實現(xiàn)了強散射介質(zhì)格林函數(shù)和散射波場的快速數(shù)值模擬，為后續(xù)進行3D全波場數(shù)值模擬和反演成像方法的研究堅定了基礎(chǔ)。

參考文獻（References）：

［1］" 孫建國.高頻漸近散射理論及其在地球物理場數(shù)值模擬與反演成像中的應(yīng)用：研究歷史與研究現(xiàn)狀概述以及若干新進展［J］.吉林大學(xué)學(xué)報（地球科學(xué)版），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］. 地球物理學(xué)報，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型偏移理論的研究歷史、研究現(xiàn)狀與發(fā)展趨勢展望：與光學(xué)繞射理論的類比、若干新結(jié)果、新認識以及若干有待于解決的問題［J］.吉林大學(xué)學(xué)報（地球科學(xué)版）， 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］" 張志佳，孫章慶，王瑞湖，等. 復(fù)雜地表起伏多重變加密網(wǎng)格有限差分法波動方程數(shù)值模擬［J］. 吉林大學(xué)學(xué)報（地球科學(xué)版），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］" 張盼，韓立國，鞏向博，等. 金屬礦地震勘探方法技術(shù)研究進展［J］. 吉林大學(xué)學(xué)報（地球科學(xué)版），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］. 地球物理學(xué)報， 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］. 地球物理學(xué)報， 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.