Lippmann-Schwinger積分方程廣義超松弛迭代解法及其收斂特性

徐楊楊， 孫建國， 商耀達， 孟祥羽， 魏脯力

吉林大學(xué)地球探測科學(xué)與技術(shù)學(xué)院， 長春 130026

0 引言

近年來，隨著陸地反射地震工區(qū)從簡單區(qū)到復(fù)雜區(qū)的擴展以及海洋反射地震從淺水向深水的轉(zhuǎn)變，如何利用散射波進行偏移成像就逐漸進入了研究者們的視野(沈鴻雁,2010；盧明輝等,2010).相應(yīng)地，散射波數(shù)值模擬方法研究也得到了前所未有的重視(孫建國，2006；Jakobsen and Wu,2016).一方面，對散射波場數(shù)值模擬的研究為認識散射波的時空分布特征與復(fù)雜構(gòu)造之間的關(guān)系提供了有效的途徑.另一方面，散射波場的快速計算方法為實現(xiàn)全波場偏移和全波形反演提供了堅實的理論基礎(chǔ)和強有力的技術(shù)支撐(王本鋒等,2016；Huang et al.,2020).

與反射波攜帶了關(guān)于反射面的幾何形態(tài)和反射面兩側(cè)的巖石物理信息類似，散射波也攜帶著與復(fù)雜構(gòu)造和復(fù)雜巖性有關(guān)的幾何和物理信息.因此，如何能在利用常規(guī)反射波實現(xiàn)成像的同時，利用散射波實現(xiàn)復(fù)雜構(gòu)造成像就顯得十分重要.

利用擾動理論，可將描述散射問題的Helmholtz方程轉(zhuǎn)化成為第二類Fredholm型積分方程，即Lippmann-Schwinger方程；簡記為LS方程.從而，散射波數(shù)值模擬問題就轉(zhuǎn)化成為了LS方程的求解問題.與常規(guī)積分方程一樣，對LS方程可以利用各種方法求解，例如Nystr?m方法、矩量法(MOM)以及迭代法(SOR)等.如果利用迭代法解LS方程，所得到的級數(shù)稱為Neumann級數(shù).在散射研究領(lǐng)域，也將Neumann級數(shù)稱為Born級數(shù)，簡記為BS.如果只取Born級數(shù)的前兩項，則可得到著名的Born近似.物理上，Born近似用散射體內(nèi)的入射場代替散射體內(nèi)的總場.

Born級數(shù)僅在小散射體或弱散射的情況下收斂.因此，Born近似是一種弱散射理論.如何將Born近似從弱散射擴展到強散射一直是散射研究中的核心問題之一.為解決這個問題，Wu和Huang(1995)、Wu(2003)將De Wolf近似引入到地震波散射級數(shù)中，對散射級數(shù)進行重整化處理，建立了聲波近似下的MFSB(Multiple Forescattering Single Backscattering)近似.在這種近似中，利用多屏單向波傳播算子進行雙域計算，有效地解決了中等強度的前向散射計算問題.然而，MFSB近似要求反向散射場遠遠小于正向散射場，從而并不能用于解決前向散射問題.Yu等(2010)、Yu和Fu(2013)通過將邊界元法(BEM)生成的邊界積分方程矩陣分解為兩部分：自作用算子和外推算子，其中外推算子由Born級數(shù)表示，并利用BEM+BS法來求解新的方程組.該算法可以克服邊界元法在復(fù)雜散射介質(zhì)情況下離散矩陣無法求逆的困境，并通過層狀非均勻介質(zhì)的地震波場數(shù)值模擬驗證了該算法的可行性.Jakobsen(2012)將量子物理中和電磁學(xué)中常用的的重整化理論應(yīng)用于地震正演模擬中，采用T算子方法求解Lippmann-Schwinger方程以改善Born級數(shù)的收斂性.Jakobsen和Wu(2016)利用De Wolf級數(shù)和T傳播算子重新推導(dǎo)了重整化后的De Wolf散射級數(shù)，并對其收斂特征進行了詳細分析.為改善Born級數(shù)的收斂性以解決強光學(xué)散射問題，Osnabrugge等(2016)通過在背景波數(shù)中引入虛部分量，得到了帶阻尼因子的背景格林函數(shù).并通過選取合適的預(yù)解因子來保證Born級數(shù)迭代的收斂性.該方法在光學(xué)散射數(shù)值模擬中表現(xiàn)出了精度高、計算效率快的優(yōu)點.Huang等(2020)從重整化角度重新闡釋了Osnabrugge得到的收斂Born級數(shù)，并通過數(shù)值算例驗證了該方法在強擾動散射介質(zhì)的地震波場正演模擬中的適用性.

在20世紀90年代，Kleinman等(1990a，b)、Kleinman和Van den Berg(1992)將解泛函方程的超松弛理論用于解決光學(xué)中的強散射問題，得到了收斂的修正Born級數(shù)迭代解.在松弛因子等于1時，該修正的Born級數(shù)退化為常規(guī)的Born級數(shù).根據(jù)Kleinman等給出的數(shù)值計算結(jié)果，超松弛修正Born級數(shù)能夠有效地解決強光學(xué)散射中的正反演問題.Yu和Fu(2014)將基于Krylov子空間投影的廣義最小余量法應(yīng)用于LS方程的求解，并通過與Yu等(2010)、Yu和Fu(2013)的BEM+BS法對比，分析了廣義最小余量法的收斂特征.

相較于源于量子力學(xué)的散射級數(shù)重整化理論，廣義超松弛迭代法理論簡潔且具有更堅實的數(shù)學(xué)基礎(chǔ).然而，大多數(shù)光學(xué)系統(tǒng)中折射率范圍在1～2之間，使得光學(xué)散射問題中的散射體擾動相對較小，僅根據(jù)Kleinman等人的工作難以對LS方程超松弛迭代級數(shù)解在地震散射問題中強速度擾動介質(zhì)波場數(shù)值模擬中的可應(yīng)用性和收斂性進行分析和評價.

本文的目的就是要解決上述問題，對Kleinman提出的LS方程廣義超松弛迭代法在反射地震正演問題中的適應(yīng)性和收斂特點進行分析和評價.在下文中，首先介紹了廣義超松弛法的基本思想和基本公式，然后通過數(shù)值分析及與有限差分法數(shù)值結(jié)果進行對比，分析了其在強散射復(fù)雜介質(zhì)下的散射波場數(shù)值模擬中的應(yīng)用效果及適應(yīng)性，并對算法收斂特點及計算效率進行分析，最后討論了松弛迭代算法進一步改進的幾個方向.

1 LS方程的廣義超松弛迭代

頻率域標量波動方程(Helmholtz)的解可以形式地表示為(A8)，即：

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

(1)

式中，U(r)為位移場，Ub(r)為背景波場，Us(r)為散射波場，ω為角頻率，vb(r)為背景速度，α(r)為速度擾動，G(r,r′)為背景介質(zhì)格林函數(shù).

在文獻中，式(1)被稱為Lippmann-Schwinger方程，簡稱為LS方程(見附錄A).形式上，LS方程 屬于第二類Fredholm型積分方程，從而在原則上可以利用迭代級數(shù)法或數(shù)值分析法求解.

(2)

則LS方程(1)式的算子形式為

LU(r)=Ub(r),

(3)

由于格林函數(shù)G(r,r′)中的Hankel函數(shù)在r=r′點包含In(k0|r-r′|)項而奇異.根據(jù)Fu等(1997)，(2)式右端的積分項在點r∈Ω上是一致連續(xù)的，即積分算子K為弱奇異Fredholm型積分算子，LS方程(1)式成為第二類Fredholm型弱奇異積分方程.由緊算子理論可知，包含In(k0|r-r′|)項的弱奇異積分算子K是映L2(Ω)到L2(Ω)的緊算子(詳細證明見：呂濤和黃晉，2013，P37定理1.5.4).因此，積分算子L=I-K也是L2(Ω)上的緊算子.應(yīng)用Fredholm擇一定理，可知LS方程存在唯一解.

根據(jù)附錄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散射級數(shù)(BS) (Kleinman and van den Berg,1991)：

(8)

根據(jù)式(8)，在第n次迭代后的結(jié)果為

(9)

若在式(9)中取U0=Ub且僅迭代一次，即可得到 廣泛使用的Born近似.當擾動強度較小或散射體較小時，Born級數(shù)收斂，從而Born近似具有與較高的精度.然而，當速度擾動強度較強或散射體尺度較大時，Born級數(shù)一般為發(fā)散級數(shù)，Born近似的精度較低.

若取C=I；B=αI(α為常數(shù))即迭代過程中保持因子α不變，則得到松弛因子為常數(shù)的超松弛迭代(GSOR-α)格式(Kleinman and van den Berg,1991)：

(10)

其第n次迭代的級數(shù)表示為

(11)

若取線性算子C=I；B=αnI(αn為復(fù)數(shù))，即每次迭代過程中按照某一收斂原則逐次修正因子αn，則得到逐次超松弛迭代(GSOR-αn)遞推形式(Kleinman and van den Berg,1991)：

(12)

及其第n次迭代的級數(shù)表示：

(13)

其中，(10)—(11)式和(12)—(13)式中的α和αn分別為復(fù)常數(shù)超松弛因子和逐次修正復(fù)松弛因子.從(6)式可以看出，迭代近似解Un是初始解和殘差的線性組合，從空間映射的角度來闡釋，即迭代近似解Un存在于初始解和殘差項的仿射空間里.

若廣義迭代格式(6)式右端項中的修正項不是簡單的殘差項的線性組合，而是在迭代過程中遞歸地構(gòu)造殘差，即第n步的殘差rn是通過算子L的某個多項式與r0相乘得到且使得所構(gòu)造的殘差在某種內(nèi)積意義下相互正交時，迭代格式則可以演變?yōu)楦鼮閺?fù)雜Krylov子空間迭代法.若Krylov子空間表示為Km(L,r0)=span{r0,Lr0,L2r0,…,Lm-1r0}，當殘差rn與子空間Km正交時，則稱為共軛梯度法(Kleinman and van den Berg，1991)，當殘差rn取最小范數(shù)時稱為廣義最小殘差法(詳見Yu and Fu，2014).

通常情況下，散射問題中的積分算子是非自伴緊隨算子，而常規(guī)的共軛梯度法僅適用于自伴隨算子，需要通過對積分算子L乘上他的伴隨算子L*來實現(xiàn)對稱化，但算子LL*的條件數(shù)是算子L條件數(shù)的平方倍，或者通過構(gòu)造預(yù)條件算子來加速收斂.廣義超松弛迭代(6)式中若取合適預(yù)條件算子C使得‖I-C-1BL‖<1，則可以達到與共軛梯度法相當?shù)氖諗啃再|(zhì).

1.2 超松弛因子選取

Petryshyn(1963)、Chertock(1968)、Kleinman和Roach(1988)等將超松弛迭代法應(yīng)用于積分方程的求解，并給出了積分算子L自伴、非自伴隨情況下松弛因子的選取范圍公式.這些松弛因子的選取公式均需要計算線性積分算子L的譜信息，而在實際問題中，計算線性積分算子L的譜信息十分困難.

為了避免求解積分算子L的譜信息，根據(jù)Kleinman和van den Berg(1991)的超松弛迭代在光學(xué)散射中的應(yīng)用經(jīng)驗，通過在迭代過程中最小化殘差‖rn‖來得到相對較優(yōu)的逐次復(fù)松弛因子αn.為求出滿足收斂性的復(fù)松弛因子αn，考慮U0=Ub時的第n次迭代結(jié)果，取αn使得‖rn‖在Hilbert空間中取極小.

定義U(r)∈L2(Ω)上的范數(shù)和內(nèi)積為

(14)

根據(jù)‖L‖2范數(shù)的定義，有：

‖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)進行離散.下面具體討論離散過程以及超松弛迭代算法離散化以及超松弛迭代算法實現(xiàn)流程.

2.1 L-S方程離散化

由于積分算子L具有弱奇異性，為降低計算復(fù)雜度，本文中的數(shù)值試驗采用基于常元插值基函數(shù)的Nystr?m法(呂濤和黃晉，2013)來求解LS方程構(gòu)成的第二類弱奇異Fredholm型積分方程.

(19)

積分號內(nèi)的被積函數(shù)的連續(xù)部分用其插值函數(shù)代替，即：

(20)

其中權(quán)函數(shù)為

(21)

由此得出LS方程的Nystr?m近似Un(r)滿足線性方程組：

(22)

采用超松弛迭代法解線性方程組(22)式，求得Un(r)后回代LS式即可得到觀測面上各接收點上的散射波場.對于多維奇異積分方程而言，常元Nystr?m積分法對于區(qū)域剖分較靈活性且離散方程構(gòu)造更簡單.因此，本文采用于常元插值基函數(shù)的Nystr?m法來求解LS方程.

(23)

則常元插值算子為

(24)

其中，Mj為Vj的中心.此時：

(25)

超松弛迭代法求解聲波散射積分方程，當速度模型復(fù)雜度較高、散射體不規(guī)則度較大時，散射體積分方程的積分區(qū)域Ω選為全空間計算來保證計算精度.而當散射體較為規(guī)則且速度模型復(fù)雜度較低時，背景介質(zhì)上速度擾動為零，則積分區(qū)域Ω可由全空間可簡化為在散射異常體Ωscatter上的積分，即只需要離散散射體即可，從而提高計算效率.

2.2 數(shù)值算法描述

超松弛迭代法求解聲波散射積分方程并計算散射波場的算法流程圖如圖1所示，迭代算法的具體實現(xiàn)步驟如下：

(1)選定待計算速度模型，數(shù)值模擬參數(shù)初始化：離散網(wǎng)格點數(shù)Nz×Nx，空間采樣間隔(dx,dz)，震源子波主頻fpeak，時間采樣間隔dt，時間采樣點數(shù)nt，源檢位置(rs,rg)以及最大迭代次數(shù)Nitmax.

(2)開始循環(huán)頻率：計算震源rs到地下介質(zhì)任意網(wǎng)格點上的入射波場Ub(rs,r)，地下介質(zhì)中任意兩點之間格林函數(shù)G(r,r′)，根據(jù)積分算子A的離散矩陣形式求出離散積分算子L.

(4)開始循環(huán)迭代計算：

①niter=niter+1；

②根據(jù)迭代遞推格式計算地下所有網(wǎng)格點上的波場Un和殘差rn；

④判斷niter是否達到最大值：若是，則退出；否則繼續(xù)步驟(4)；

(5)開始循環(huán)接收點：計算背景速度下接收點rg到地下任意網(wǎng)格點上的格林函數(shù)G(r,rg)，回代計算當前頻率下所有接收點的波場值U(r,rg).

(6)繼續(xù)循環(huán)頻率，執(zhí)行步驟(2)—(5)，直至求得所有頻率下的波場值U(r,rg)，進行 FFT逆變換，輸出時域散射波場Us.

圖1 算法實現(xiàn)程序流程圖Fig.1 Flow chart of algorithm implementation

3 數(shù)值試驗

3.1 球體模型

對簡單球狀散射體模型進行正演模擬，速度模型為均勻背景介質(zhì)中含有一半徑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)、復(fù)常數(shù)廣義超松弛迭代(GSOR-α)以及Born級數(shù)迭代(BS)四種數(shù)值方法計算了單炮記錄(如圖2e、f、g、h)，其中GSOR-α迭代法中松弛因子α=α1=(r0,Lr0)/‖Lr0‖2且在迭代過程中保持不變.從圖中可以看出速度擾動強度較低時，三種迭代方法都收斂，得出的數(shù)值計算結(jié)果均與有限差分正演結(jié)果具有較好的一致性.抽出GSOR-αn法的第80道、151道、230道與FD法進行比較(如圖3)，結(jié)果顯示，兩種算法得到的結(jié)果差別較小，均在允許誤差之內(nèi).

為進一步分析廣義逐次超松弛迭代法(GSOR-αn)的收斂性質(zhì)，取計算頻率f=30 Hz，復(fù)松弛因子αn-迭代次數(shù)、歸一化殘差Error-迭代次數(shù)(Error<10-6)的變化關(guān)系如圖2c、d.如圖2c，隨著迭代次數(shù)增加，殘差無限接近于0時，復(fù)松弛因子在αn=1.0-0.2i附近浮動，趨向于最優(yōu)選擇.從圖2d中可以看出，當滿足殘差Error<10-6時，GSOR-αn法的收斂速度較GSOR-α和BS法更快,所需要的迭代次數(shù)更少.圖2b的頻率-迭代次數(shù)關(guān)系曲線表明，在0～60 Hz頻率范圍內(nèi)，尤其是高頻點處，GSOR-αn法迭代收斂速度要明顯優(yōu)于GSOR-α和BS法.

3.2 規(guī)則多體散射模型

考慮地下介質(zhì)中沿水平方向、垂直方向分布大小不一的多個散射體的情況，速度模型如圖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隨迭代次數(shù)的變化關(guān)系，結(jié)果顯示，在迭代相同次數(shù)時，GSOR-αn迭代法的殘差與BS迭代法的殘差小兩個數(shù)量級.圖4c中，頻率與迭代次數(shù)的次數(shù)的關(guān)系曲線表明(殘差滿足Error<10-6)，當f<30 Hz時，三種迭代方法收斂速度性差較小，當f>30 Hz時，BS迭代法收斂速度變慢，當頻率f=60 Hz時，GSOR-αn和GSOR-α迭代法所需要的迭代次數(shù)約為BS迭代法的1/4.

圖2 球狀散射體模型的有限差分和迭代法數(shù)值結(jié)果 (a) 速度模型； (b) GSOR-αn(實線)、GSOR-α(點劃線)和BS(虛線)法頻率-迭代次數(shù)關(guān)系； (c) GSOR-αn法復(fù)松弛因子αn-迭代次數(shù)關(guān)系(實部：實線，虛部：虛線)； (d) GSOR-αn(實線)、GSOR-α(點劃線)和BS(虛線)歸一化殘差Error-迭代次數(shù)關(guān)系； (e) FD數(shù)值結(jié)果； (f) GSOR-αn迭代結(jié)果； (g) GSOR-α迭代結(jié)果; (h) BS迭代結(jié)果.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(虛線)數(shù)值結(jié)果單道對比圖Fig.3 Comparison of single-trace seismic record using FD and GSOR-αn method

圖4 多散射體模型的有限差分和迭代法數(shù)值結(jié)果 (a) 速度模型； (b) GSOR-αn(實線)、GSOR-α(點劃線)和BS(虛線)歸一化殘差Error-迭代次數(shù)關(guān)系； (c) GSOR-αn(實線)、GSOR-α(點劃線)和BS(虛線)法頻率-迭代次數(shù)關(guān)系； (d) FD數(shù)值結(jié)果； (e) GSOR-αn迭代結(jié)果； (f) GSOR-α迭代結(jié)果; (g) BS迭代結(jié)果.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迭代法(虛線)數(shù)值結(jié)果單道對比圖Fig.5 Comparison of single-trace seismic records using FD and GSOR-αn method

圖6 鹽丘模型的有限差分和迭代法數(shù)值結(jié)果 (a) 速度模型； (b) GSOR-αn(實線)、GSOR-α(點劃線)和BS(虛線)法頻率-迭代次數(shù)關(guān)系； (c) GSOR-αn法復(fù)松弛因子αn-迭代次數(shù)關(guān)系(實部：實線，虛部：虛線)； (d) GSOR-αn(實線)、GSOR-α(點劃線)和BS(虛線)歸一化殘差Error-迭代次數(shù)關(guān)系； (e) FD數(shù)值結(jié)果； (f) GSOR-αn迭代結(jié)果； (g) GSOR-α迭代結(jié)果; (h) BS(1 order)迭代結(jié)果.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(點劃線)數(shù)值結(jié)果單道對比圖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四種數(shù)值方法得到的單炮地震記錄，數(shù)值計算結(jié)果表明三種迭代方法均可以得到與有限差分法相當?shù)慕Y(jié)果，且可以清晰分辨出各個散射體所對應(yīng)的同相軸.隨機抽取GSOR-αn和FD法的3道進行對比(如圖5)，可以看出，對于散射體較多的地質(zhì)模型，相對于精細網(wǎng)格下的有限差分數(shù)值結(jié)果，粗網(wǎng)格的GSOR-αn迭代法同相軸在局部出現(xiàn)同相軸拉伸現(xiàn)象，這一點可以通過加密網(wǎng)格或增加迭代次數(shù)進一步減小殘差來改善計算結(jié)果.

3.3 鹽丘模型

下面通過將廣義超松弛迭代法應(yīng)用于散射體尺度較大且速度擾動較強的地下介質(zhì)中來驗證超松弛迭代方法對強散射介質(zhì)的適應(yīng)性.模型參數(shù)如圖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時的復(fù)松弛因子αn的實虛部-迭代次數(shù)關(guān)系曲線和三種迭代方法的殘差Error-迭代次數(shù)關(guān)系曲線.如圖6c，隨著迭代次數(shù)增加，殘差Error不斷趨于0的過程中，復(fù)松弛因子在αn=0.6-0.5i附近浮動.圖6d表明，Born散射級數(shù)發(fā)散的情況下，超松弛迭代法依舊收斂，且GSOR-αn迭代法收斂速度較快，所需要的迭代次數(shù)更少.從圖6b的頻率域迭代次數(shù)的關(guān)系曲線可以看出，對于強散射介質(zhì)，當f>10 Hz時，Born散射級數(shù)已經(jīng)不再收斂，而GSOR-α和GSOR-αn迭代法仍然收斂，并且GSOR-αn迭代法保持了較好的收斂性質(zhì).四種數(shù)值方法得到的單炮地震記錄見圖6e、f、g、h，其中圖6h為Born散射級數(shù)法僅進行一次迭代后得到的散射波場(Born近似法)，數(shù)值結(jié)果表明有限差分法得到的波場在鹽丘上界面散射能量較強，但在巖體上側(cè)的尖點以及鹽丘下界面的散射能量較弱，GSOR-α和GSOR-αn迭代法在巖體上下界面及尖點處散射能量較強且同相軸清晰，受高速鹽丘散射體的屏蔽影響較小.隨機抽取3道進行對比，如圖7所示，從圖中可以看出，廣義超松弛迭代法得到的結(jié)果在鹽丘上界面與有限差分法子波時間和相位一致，而在尖點散射以及巖體下界面處有限差分法得到的結(jié)果振幅較弱.此外，Born近似得到的散射波場在鹽丘下界面處的子波出現(xiàn)相位延遲現(xiàn)象，誤差較大.

4 計算時間對比

數(shù)值試驗過程中，球狀散射體模型和鹽丘模型采用局部積分法，多散射體模型采用全空間積分法進行數(shù)值模擬計算.數(shù)值計算過程均在微機上實現(xiàn)，CPU：4核，內(nèi)存：8GB.

四種數(shù)值方法計算時間對比如表1中，可以看出，當?shù)叵陆橘|(zhì)中散射體分布情況較為簡單時(如，單體金屬礦體、鹽丘等模型)，由于散射體周圍介質(zhì)中擾動為0，積分法中的積分域退化為散射體本身，超松弛迭代法相對有限差分法計算量較小，具有較高的計算效率.數(shù)值計算過程中，雖然GSOR-αn法更新松弛因子會增加一定計算量，但其收斂速度較快所需要的迭代次數(shù)較少，相對于傳統(tǒng)的BS法并未增加過多的計算時間.而在散射體較多且分布情況較為復(fù)雜的情況下，積分法中的積分域擴大，生成線性方程組系數(shù)矩陣會占用大量計算時間，相對于有限差分法，廣義超松弛迭代法計算效率較低，有待進一步改進算法以提高計算效率.

表1 FD法與GSOR-αn迭代法計算時間對比Table 1 Comparison of calculation time of FD method and GSOR-αn method

5 分析與討論

文中采用積分方程法求解頻率域標量波動方程，通過廣義超松弛迭代法來解散射問題的Lippmann-Schwinger積分方程，并將該方法應(yīng)用于強擾動介質(zhì)模型的地震散射波場正演模擬計算.首先引入復(fù)松弛因子αn來對Lippmann-Schwinger積分方程的積分算子進行修正，得到新的算子方程，并給出了收斂的修正Born散射級數(shù)迭代形式.松弛因子αn的選擇是根據(jù)Kleinman和Van den Berg(1992)給出的通過在迭代過程中最小化剩余函數(shù)rn來逐步修正超松弛因子，以加快修正Born散射級數(shù)的收斂速度.地震散射波場正演模擬數(shù)值結(jié)果表明通過廣義超松弛迭代法得到的修正Born散射級數(shù)在散射強度較小、復(fù)雜度較低的速度模型的上(如球狀散射體模型)，其收斂速度比傳統(tǒng)的Born散射級數(shù)更快；在遇到散射強度較大的速度模型時(如鹽丘模型)，仍保持了較好的收斂性質(zhì)，但需要更多的迭代次數(shù)來保證收斂性.文中給出的松弛因子是Kleinman等根據(jù)一般矩陣理論推出的滿足收斂要求的計算公式，并未考慮實際特定物理問題的漸進解和近似解.在今后研究方向上，可以根據(jù)實際地震散射問題，尋找更加適合的預(yù)解條件以及松弛因子，以提高復(fù)雜問題數(shù)值模擬的計算精度.

另外，數(shù)值實現(xiàn)過程中需在在每次迭代計算之前更新松弛因子αn以加速收斂，故相對于傳統(tǒng)的Born級數(shù)會增加計算量.并且波場值的精度依賴于迭代的次數(shù)，該算法存在一定程度上的計算精度與計算效率上的折中.為節(jié)約計算時間和內(nèi)存、提高計算效率，可從以下幾個方面進行改進：(1)考慮到積分算子A的形式為卷積積分形式，可對其進行空間傅里葉變換到波數(shù)域中來避免計算復(fù)雜卷積積分；(2)頻率計算的并行化；(3)迭代計算中進行數(shù)據(jù)劃分與迭代空間相結(jié)合，實現(xiàn)超松弛迭代法的模塊并行化.

6 結(jié)論

Lippmann-Schwinger方程的廣義超松弛迭代法是一種針對強散射的數(shù)值模擬方法，主要用于解決強光學(xué)散射中的正反演問題.傳統(tǒng)的Born散射級數(shù)受弱散射假設(shè)限制，且只適用于短程傳播，在解決速度擾動較強、散射體尺度較大的地震散射問題時并不能得到理想的效果.為了改善Born散射級數(shù)的收斂性，本文將廣義超松弛迭代法應(yīng)用于地震散射波場正演模擬中，給出了LS方程所滿足的Banach空間下線性算子方程的廣義超松弛迭代格式以及松弛因子的選取原則，并設(shè)計了超松弛迭代算法的具體實現(xiàn)流程.通過在理論模型上的應(yīng)用效果，驗證了超松弛迭代法具有良好的收斂特性，受數(shù)值頻散影響較小，在網(wǎng)格間距、時間采樣間隔較大的情況下可以得到與有限差分相當?shù)臄?shù)值模擬結(jié)果.與散射級數(shù)重整化相比，超松弛法具有更堅實的數(shù)學(xué)基礎(chǔ)，可以有效地克服Born級數(shù)的弱點，較好的解決地震波強散射問題.超松弛法在散射地震波數(shù)值模擬、全波場偏移成像和全波形反演等方面具有較大的應(yīng)用潛力，值得進行進一步的研究和改進.未來研究擬整合積分方程FFT快速算法和MPI+openmp并行框架，進一步降低存儲量、提高計算效率，為后續(xù)進行高效的地震正反演方法研究奠定基礎(chǔ).

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

考慮頻率域標量波動(Helmholtz)方程：

(A1)

式中；ω為角頻率，v(r)為介質(zhì)波速，U(r)為位移場，s(ω)δ(r-rs)為Dirac函數(shù)表示的震源項.對速度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)

均勻背景介質(zhì)空間的格林函數(shù)G(r,r′)滿足Helmhlotz方程:

(A7)

利用第二格林定理，全空間假設(shè)下總場的積分方程表達式為

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

(A8)

式中：G(r,r′)為二維全空間格林函數(shù)，表達式為

(A9)

(A8)式為散射問題的形式解，稱為Lippmann-Schwinger(LS)方程.

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

迭代法是求解線性算子方程Au=f近似解的一種重要方法，其中算子A可以是任意線性賦范空間X中的矩陣、積分算子或其他抽象算子，f為X中的已知項.迭代法的基本思想是：從任一初始解u0出發(fā)，按照某一規(guī)則，逐次構(gòu)造一個迭代解un，當un收斂于u*時，使u*是所給算子方程Au=f的解.對于迭代級數(shù)解u*存在性和唯一性，有如下定理Petryshyn(1962)：

定理1對任意f∈X，方程Au=f存在唯一解u*，當且僅當存在連續(xù)可逆算子C和有界線性算子B，使得級數(shù)：

(B1)

定義1設(shè)線性算子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相當于預(yù)條件算子，構(gòu)造合適的預(yù)條件算子使得T=I-C-1BA=I-Lε(‖Lε‖<1)，將會加快迭代格式(B2)收斂速度，如取B=αI(α為常數(shù))，此時預(yù)條件算子C可以看作是A的逆算子.相較于Krylov子空間迭代法(如共軛梯度法、廣義最小余量法等)，廣義超松弛迭代法數(shù)學(xué)實現(xiàn)簡單，計算量較小，通過選擇合適的線性算子B和C，該方法可以達到與共軛梯度法相當?shù)氖諗克俣?如果把Lippmann-Schwinger方程中的積分算子視為Banach空間中的映射算子，可以將廣義超松弛迭代法用于地震散射問題的LS積分方程.