徐楊楊， 孫建國 ， 商耀達

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

0 引言

類似于光學散射和電磁散射問題的數(shù)值模擬正演方法，常規(guī)的聲波散射數(shù)值模擬方法也分為兩大類：基于微分方程理論的有限差分法(Finite Difference)、有限元法(Finite Element)等和基于積分方程理論的邊界元法、矩量法、Nystr?m法等(Fu et al.,1997;Fu and Wu,2000; Hsiao et al.,2000；劉立彬等，2020；邱長凱等，2020；徐楊楊等，2021).微分方程類的數(shù)值算法具有容易實現(xiàn)、計算速度快等優(yōu)點，但是網(wǎng)格剖分需要在全空間上進行，對于較小的異常體也會生成規(guī)模非常大的離散矩陣，并且容易受到數(shù)值頻散的影響.而積分方程類的數(shù)值算法具有半解析性質(zhì)、只在散射異常體上進行網(wǎng)格離散即可.因此，對小型散射異常體的散射波場數(shù)值模擬，積分方程法生成的離散矩陣規(guī)模小，降低了計算和存儲的要求，其計算速度和精度要明顯優(yōu)于基于微分方程的數(shù)值算法.對地震散射問題中的大尺度復雜地質(zhì)條件而言，由于積分方程法形成的線性方程組系數(shù)矩陣的密實性，并且離散矩陣元素的計算過程復雜，計算耗時、存儲量大，相比于計算電磁學，基于積分方程的數(shù)值方法在地球物理中的應用受到了非常大的限制(王德智，2015).

針對積分方程法處理散射問題時所存在的問題，國內(nèi)外學者在積分方程快速算法方面作了大量研究工作.由于求解積分方程數(shù)值法生成的大型線性方程組非常耗時，有學者提出采用近似解法來避免解方程組.Zhdanov和Fang(1996)、Zhdanov等(2000)提出了一種擬線性的近似解法，認為散射異常場是背景場與反射系數(shù)的乘積，反射系數(shù)是最小二乘問題的解，該方法在電磁散射正反演中取得較好的應用效果.Sun(2005)根據(jù)計算電磁學的擬解析線性近似技術(shù)，將電反射率張量的概念引入直流電場積分方程數(shù)值法中，給出了擬線性級數(shù)的表達形式，提供了一種新的直流電場的正反演方法.孫建國(2006)借鑒電磁散射問題中的擬線性近似技術(shù)，通過假設(shè)散射波場與背景波場之間的線性關(guān)系，建立了聲波散射的擬解析近似表達式.擬解析近似是半解析解，只需計算數(shù)值積分，該方法可以有效的避開了復雜介質(zhì)情況下大型線性方程組的求解.

殷長春等(2018)利用插值技術(shù)將準線性近似算法應用于三維情況下航空電磁數(shù)據(jù)的數(shù)值模擬計算，在稀疏網(wǎng)格上計算散射場和點反射系數(shù)，然后進行線性插值，最后根據(jù)擬線性近似技術(shù)來求解電磁異常體的散射波場.另一類方法利用格林函數(shù)的Toeplitz性質(zhì)，通過快速傅氏變換(Fast Fourier Transform，簡稱FFT)來加速線性方程組的迭代類求解方法中的矩陣向量乘積，稱為積分方程快速傅里葉變換法(Integral Equation Fast Fourier Transform，簡稱IE-FFT)(Peters and Volakis,1988; Volakis and Barkeshli,1991)，達到降低存儲和計算復雜的目的，如共軛梯度快速傅里葉變換(Conjugate Gradient Fast Fourier Transform，簡稱CG-FFT).Osnabrugge 等(2016)在背景波數(shù)中引入虛部分量，得到了帶阻尼因子的背景格林函數(shù).通過加入預條件因子進一步加速了迭代的收斂性，同時利用FFT實現(xiàn)了循環(huán)卷積的快速算法.該方法在光學散射數(shù)值模擬中表現(xiàn)出了精度高、計算效率快的優(yōu)點.

此外，傳統(tǒng)的求解積分方程的矩量法計算復雜，離散后的線性方程組系數(shù)矩陣每個元素都需要計算內(nèi)積積分，計算量大，效率低.雖然快速多極子方法(Fast Multiple Method，簡稱FMM)和多層快速多極子方法(Multievel Fast Multiple Algorithm，簡稱MLFMA)(Rokhlin,1990; ?akir,2006; ?akr,2009；崔曉娜和劉洪，2020)能夠有效地降低存儲和計算復雜度，但對積分核依賴性強，且實現(xiàn)過程較為復雜.

除了上述方法之外，在計算數(shù)學中還有一種解積分方程的方法，稱為Nystr?m法.該方法利用適當?shù)臄?shù)值求積公式來近似積分算子(Nyst?m離散)，可以避免生成系數(shù)矩陣元素時的大量積分計算，相較于傳統(tǒng)的矩量法計算成本較低.如果積分方程的核函數(shù)只由Green函數(shù)組成，則基于數(shù)值積分公式的Nystr?m法所生成系數(shù)矩陣只對積分核進行簡單計算并且具有Toeplitz性質(zhì)，計算復雜度低且易實現(xiàn)并行化.對于奇異性的處理，可以采用階數(shù)較高的插值基函數(shù)、奇異性減法等技術(shù)均可處理奇異性(Canino et al.，1998).

然而，在地震散射問題中出現(xiàn)的L-S方程的積分核除了Green函數(shù)之外還包含有速度擾動項.因此，為了使得積分方程的核函數(shù)只包含Green函數(shù)，首先將L-S方程改寫為一個與原積分方程等價的方程(等價L-S方程)，并對該方程進行逐點歸一化，進而得到在離散后其系數(shù)矩陣具有Toeplitz性質(zhì)的歸一化等價L-S方程.然后，對歸一化等價L-S方程進行Nystr?m離散，從而得到其系數(shù)矩陣具有Toeplitz性質(zhì)的線性方程組，極大地降低了存儲量.同時，為了提高空間褶積的計算效率，將歸一化等價L-S方程變換到波數(shù)域，在其傅里葉譜空間中實現(xiàn)矩陣和向量之間的乘積，進而實現(xiàn)對線性方程組的快速求解.通過簡單的分析可知，上述兩個措施可將內(nèi)存和計算復雜度由O(N×N)分別降為O(N)和O(NlogN)(N=NxNz，Nx、Nz分別為x、z方向上的網(wǎng)格剖分數(shù)量).此外，因數(shù)值模擬在頻率域進行，各頻率下的波場值計算是獨立的.因此，為進一步提高計算效率，設(shè)計了基于MPI+OpenMP的進程級和線程級混合并行的實現(xiàn)模式，在有效降低其內(nèi)存消耗的基礎(chǔ)上最大限度地提高積分方程法散射波場數(shù)值模擬的計算效率.

1 等價與歸一化等價L-S方程

利用擾動理論，可將描述散射問題的Helmholtz方程轉(zhuǎn)化成為Lippmann-Schwinger(L-S)積分方程：

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

(1)

由于格林函數(shù)G(r,r′)中的Hankel函數(shù)在r=r′點包含In(k0|r-r′|)項而奇異，考慮Cauchy主值意義下的收斂性，(1)式右端的積分項在點r∈Ω上一致連續(xù)(Fu et al.,1997)，即L-S方程具有弱奇異性.應用Fredholm擇一定理，可知L-S方程存在唯一解.

(2)

(3)

這是一個積分核為Green函數(shù)的積分方程.因此，如果首先對方程(2)進行某種方式的數(shù)值離散，然后再進行形如方程(3)的逐點歸一化處理，則可得到其系數(shù)矩陣具有Toeplitz性質(zhì)的線性方程組.

2 等價L-S方程的Nystr?m離散

令Pn是插值算子，{rj,j=1,2,…,n}?Ω是插值基點，{lj(r′),j=j=1,2,…,n}是插值基函數(shù)，滿足lj(r′i)=δij.于是：

(4)

定義弱奇異積分算子：

(5)

則其近似算子滿足：

(6)

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

(7)

(8)

(9)

定義常元插值算子：

(10)

其中，Mj為剖分單元Vj的中心.此時，需要計算弱奇異積分：

(11)

從而得到權(quán)系數(shù).由于當r∈Vj且r=r′時G(r,r′)奇異，所以必須要對(11)進行奇異性消除處理，具體見劉寧和孫建國(2007)或張金會和孫建國(2009) .

3 線性方程組的逐次超松弛迭代

Lippmann-Schwinger積分方程的常元Nystr?m近似生成的系數(shù)矩陣為稠密矩陣，并且每一個矩陣元都必須計算格林函數(shù)在剖分元上的弱奇異積分.當散射體規(guī)模較大時，網(wǎng)格剖分數(shù)量也很大，在二維情況下，系數(shù)矩陣由N2個元素組成(N=NxNz，Nx、Nz分別為x、z方向上的網(wǎng)格剖分數(shù)量).采用高斯消元法或LU分解等直接法求解時需要存儲整個矩陣及其逆矩陣，內(nèi)存需求為O(N2)并且矩陣求逆計算復雜度為O(N3).而迭代法無需存儲逆矩陣，其復雜度下降為O(N2).常用的迭代算法有Gauss-Seisel迭代、超松弛迭代(SOR)、共軛梯度法、廣義最小殘差法以及其他Krylov子空間法(Kleinman et al., 1990a,b; Kleinman and Van Den Berg,1991; Yu and Fu, 2014).本文采用逐次超松弛迭代法(Kleinman and Van Den Berg, 1991)來求解L-S方程離散生成的線性系統(tǒng).

(12)

(13)

式中，βn為復逐次超松弛因子.根據(jù)Kleinman和Van Den Berg(1991)的超松弛迭代在光學散射中的應用經(jīng)驗，通過在迭代過程中最小化殘差‖rn‖ 來得到逐次復松弛因子βn，即取βn使得 ‖rn‖ 在Hilbert空間中取極?。?/p>

(14)

4 矩陣向量乘積的FFT實現(xiàn)

(15a)

(15b)

其中：

(15c)

對于原L-S積分方程，核函數(shù)為α(r′)G(r,r′).由于速度擾動函數(shù)在每個面元上的值不同，導致帶速度擾動權(quán)的積分核經(jīng)數(shù)值離散后的系數(shù)矩陣不具有Toeplitz性質(zhì).因此，有必要對L-S方程進行改寫，使得積分核僅包含Green函數(shù)，保證離散后的系數(shù)矩陣具有Toeplitz性質(zhì).

任給一n階循環(huán)矩陣C=circ(c0,c1,…,cn)可以被傅里葉矩陣F對角化(Davis，2013)，即:

(16)

(17)

對于一個n階對稱Toeolitz矩陣T，可以通過將它嵌到2n階循環(huán)矩陣C中，來計算矩陣向量乘積，即:

(18)

從而有:

(19)

根據(jù)式(19)可以得到矩陣向量乘積Tx.

(20)

其中:

0≤i≤Nx-1，

(21)

5 頻率域散射波場正演的并行化

MPI并行標準庫通過消息傳遞來實現(xiàn)程序的并行化，是一種進程級的并行工具.它具有效率高、可移植性強等優(yōu)勢，是目前實際并行計算中普遍采用的并行編程環(huán)境.OpenMP作為一種在內(nèi)存共享基礎(chǔ)上實現(xiàn)多個線程間的并行計算編程接口，通過簡單的高級語言指令即可實現(xiàn)多線程并行化計算(Quinn，2004).Nystr?m法求解Lippmann-Schwinger積分方程是在頻率域中實現(xiàn)的，各頻率片之間是獨立計算的，可以通過MPI并行處理技術(shù)實現(xiàn)頻率波場計算并行化.此外，當?shù)叵律⑸洚惓ｓw規(guī)模較大時，剖分網(wǎng)格數(shù)量將會很大，求解單頻率下的線性方程組時通常含有大量稠密矩陣向量乘法運算相關(guān)的for循環(huán)，在單個節(jié)點上使用OpenMP實現(xiàn)for循環(huán)的并行化可進一步提高效率.基于MPI+OpenMP的混合并行模式可以減少MPI的通信開銷，充分利用資源，得到更高的加速比，在有效降低其內(nèi)存消耗的基礎(chǔ)上最大限度地提高積分方程法散射波場數(shù)值模擬的計算效率.

5.1 基于MPI的頻率并行

基于MPI并行實現(xiàn)頻率片并行化的算法設(shè)計過程中，為獲得單炮所有頻率下的波場值，可以采用主從式并行模式，主進程進行消息發(fā)送與接收以及其他輸入輸出任務，從進程負責所有頻率的波場值計算.具體并行實現(xiàn)為主進程進行各從進程頻率數(shù)和偏移量的計算并將消息發(fā)送給從進程，當從進程計算完相應頻率數(shù)的波場值后，主進程接收各從進程發(fā)送的相應頻率數(shù)、偏移量和波場值，最終得到所有接收道和所有頻率波場值生成的二維數(shù)組.基于MPI并行的散射波場數(shù)值模擬實現(xiàn)步驟如下：

(1)調(diào)用MPI_Init函數(shù),并生成通信域，參數(shù)初始化；

(2)根據(jù)集群中的節(jié)點數(shù)量，調(diào)用MPI_Comm_size函數(shù)來確定MPI_COMM_WORLD中的總進程數(shù)量(一個節(jié)點執(zhí)行一個進程)，并調(diào)用MPI_Comm_rank確定進程編號；

(3)主進程根據(jù)設(shè)定的總進程數(shù)計算各從進程需要計算的頻率數(shù)和頻率偏移量并發(fā)送給相應的從進程；

(4)各從進程接收所需計算的頻率數(shù)和頻率偏移量，并計算相應頻率數(shù)下的所有道的波場值，并將計算結(jié)果發(fā)送回主進程；

(5)主進程接收從進程發(fā)送的頻率數(shù)、頻率偏移量和計算結(jié)果，對所有頻率下的波場值進行逆傅氏變換得到時域單炮記錄；

(6)調(diào)用MPI_Finalize，結(jié)束并行環(huán)境.

5.2 基于OpenMP的節(jié)點內(nèi)for循環(huán)并行

L-S方程的Nystr?m近似迭代解法實現(xiàn)過程中涉及以for循環(huán)形式出現(xiàn)的稠密矩陣向量乘法運算，而線程級并行的OpenMP并不需要進行復雜的線程創(chuàng)建、同步、銷毀等工作.基于OpenMP的節(jié)點內(nèi)循環(huán)并行化實現(xiàn)步驟如下：

(1)按照單個計算節(jié)點上的CPU核的數(shù)目，調(diào)用omp_set_num_threads設(shè)定線程數(shù)；

(2)在無數(shù)據(jù)相關(guān)性的for循環(huán)外部調(diào)用#pragma omp parallel for語句，將該循環(huán)并行化；

(3)調(diào)用private(〈variable list〉)和shared(〈variable list〉)設(shè)置并行域內(nèi)的私有變量和共享變量.

6 數(shù)值試驗

下面通過數(shù)值試驗來測試基于Nystr?m法的IE-FFT并行算法在解決地震散射波場正演模擬中的有效性以及計算效率.首先設(shè)計了球形單體散射模型，通過與直接解法得到的數(shù)值結(jié)果進行對比，分析IE-FFT算法在存儲和計算效率上的優(yōu)勢.然后用更為復雜的SEG/EAGE的鹽丘模型對快速算法進行進一步測試.數(shù)值試驗環(huán)境：Linux操作系統(tǒng)，包含11個節(jié)點，每個節(jié)點2個CPU，每個CPU包含8個核(CPU型號：Intel(R) Xeon(R) CPU E5-2620 0 @ 2.00 GHz )，使用基于MPI+OpenMP的并行算法對散射模型進行單炮數(shù)值模擬.

6.1 球體散射模型

球狀散射體速度模型為均勻的背景介質(zhì)中含有一半徑200 m的高速球狀散射體(如圖1).模型大小為3000 m×3000 m，球散射體速度為3000 m×3000 m，背景速度2000 m·s-1，震源為主頻15 Hz的雷克子波，震源位置rs=(1500 m,10 m)，檢波器置于z=10 m的界面上，空間采樣間隔為dx=dz=5 m，時間采樣間隔為4 ms，記錄時間2.4 s.

圖1 球體散射模型Fig.1 Sphere scattering model

基于MPI+OpenMP，開辟11個進程進行基于MPI的頻率并行計算，每個進程開6個線程進行基于OpenMP的for循環(huán)并行計算，分別采用直接法和IE-FFT迭代算法計算球體散射模型正演結(jié)果.對于單球體散射速度模型，僅需要在散射異常體上進行離散.兩種方法得到的單頻波場(f=15 Hz)和單炮記錄結(jié)果如圖2和圖3，IE-FFT迭代快速算法與直接法求解得出的數(shù)值結(jié)果具有較好的一致性.為了更細致的分析，分別從兩種算法得到的單炮記錄中隨機抽取2道進行對比(如圖4)，結(jié)果顯示出，IE-FFT算法得到的數(shù)值結(jié)果在同相軸能量上與直接求解法結(jié)果有極小差別，量級在誤差允許范圍之內(nèi)，兩種方法同相軸相位對應一致.

圖2 球體散射模型直接法與IE-FFT迭代法15Hz單頻波場數(shù)值結(jié)果(a) 直接法單頻波場實部； (b) 直接法單頻波場虛部； (c) IE-FFT迭代法單頻波場實部； (d) IE-FFT迭代法單頻波場虛部.Fig.2 Numerical results of 15Hz single-frequency wave field using the direct method and IE-FFT iterative method on the sphere scattering model(a) The real part of the single-frequency wave field using direct method； (b) The imaginary part of the single-frequency wave field using direct method； (c) The real part of the single-frequency wave field using the IE-FFT iterative method； (d) The imaginary part of the single-frequency wave field using the IE-FFT iterative method.

圖3 球體散射模型直接法(a) 與IE-FFT迭代法; (b) 數(shù)值結(jié)果.Fig.3 Numerical results using the direct method (a) and IE-FFT iterative method (b) on the sphere scattering model

圖4 球體散射模型直接法(實線)與IE-FFT迭代法(點劃線)單道對比(a) 第40道單道對比； (b) 第150道彈道對比.Fig.4 Comparison between the direct method (solid line) and the IE-FFT iterative method (dotted line) for the sphere scattering model(a) Comparison of trace 40； (b) Comparison of trace 150.

球體散射模型直接法與IE-FFT迭代法的CPU時間與內(nèi)存消耗對比見表1.數(shù)值實現(xiàn)過程中，直接法耗費CPU時間為4624.27 s，每個進程內(nèi)存消耗為515 M，IE-FFT法迭代50次耗費CPU計算時間為396.37 s，每個進程內(nèi)存消耗為78 M.IE-FFT迭代算法的數(shù)值結(jié)果可以較好地吻合直接法結(jié)果，而計算速度是直接法的11倍，并且應用格林函數(shù)的Toeplitz性質(zhì)降低了系數(shù)矩陣的內(nèi)存需求，內(nèi)存消耗僅約為直接法的1/7.

6.2 SEG/EAGE鹽丘模型

當?shù)叵陆橘|(zhì)中包含較大尺度的復雜散射體且速度擾動更強時，積分方程法應用于此類介質(zhì)模型時需要在全空間中進行離散計算.而直接法因其離散后的矩陣為NxNz階方陣，普通的微機難以獲取相應的棧內(nèi)存來存儲系數(shù)矩陣且計算效率不滿足實際應用需求.為了驗證IE-FFT迭代算法在復雜強散射介質(zhì)中的適應性，選取SEG/GAGE鹽丘速度模型.模型參數(shù)如圖5所示，模型大小為7800 m×2090 m，背景速度1500 m·s-1，震源子波為主頻15 Hz的雷克子波，源位置rs=(3800 m,10 m)，檢波器置于z=10 m的界面上，迭代法空間采樣間隔為dx=dz=10 m，時間采樣間隔為4 ms，記錄時間4 s.

圖5 SEG/EAGE鹽丘模型Fig.5 SEG/EAGE salt dome model

IE-FFT快速迭代算法得到的單頻波場(f=15 Hz)和單炮地震記錄如圖6和圖7所示.從單炮地震記錄上可以清晰分辨鹽丘頂端以及巖體右側(cè)橫向起伏變化劇烈處的散射同相軸，對強散射體的識別能力較強，并且模型底部測試層的同相軸清晰、能量強，并未受到高速異常體的屏蔽影響.在直接法無法求解大型線性方程組的情況下，IE-FFT快速迭代法迭代100次耗費CPU計算時間和每個進程內(nèi)存消耗為為9466.54 s和321 M(見表1).

圖6 IE-FFT迭代算法15HZ單頻波場數(shù)值結(jié)果(a) 單頻波場實部； (b) 單頻波場虛部.Fig.6 Numerical results of 15Hz single-frequency wave field using the IE-FFT iterative method on the SEG/EAGE salt dome model(a) The real part of the single-frequency wave field； (b) The imaginary part of the single-frequency wave field.

圖7 IE-FFT迭代算法散射波場數(shù)值結(jié)果Fig.7 Numerical result using the IE-FFT iterative method on the SEG/EAGE salt dome model

表1 直接法與IE-FFT迭代法CPU時間與內(nèi)存消耗Table 1 CPU time and memory consumption of the direct method and the IE-FFT iterative method

7 結(jié)論

傳統(tǒng)上，對L-S方程的數(shù)值求解一般采用矩量法.此時，相應線性方程組的系數(shù)矩陣元素都是以內(nèi)積積分的形式出現(xiàn)的.雖然可以利用數(shù)值積分法計算這些內(nèi)積積分且可以獲得較高的精度，但是其計算過程復雜、對計算機資源要求高且效率低.與此相反，基于數(shù)值積分公式的Nystr?m離散可使系數(shù)矩陣只包含對積分核(Green函數(shù))的積分計算.因此，從計算效率方面考慮，對于大尺度地質(zhì)條件下的地震散射波場計算問題，Nystr?m法會更加適合.事實上，通過對L-S方程的變形處理，可使得最終形成的線性方程組的系數(shù)矩陣具有Toeplitz性質(zhì)，進而可將存儲量由對N階矩陣直接進行存儲的O(N×N)降為O(N).此外，利用FFT和IFFT計算離散空間褶積，可以加速線性方程組求解迭代計算過程中的矩矢乘積，使其計算復雜度由O(N×N)降為O(NlogN).再有，根據(jù)各單頻波場在計算上相互獨立的特點，采用了基于MPI+OpenMP的進程級和線程級相結(jié)合的并行化計算方案.數(shù)值結(jié)果表明，相較于傳統(tǒng)的積分方程數(shù)值算法，利用Nystr?m離散改寫后的等價L-S方程且同時采用FFT計算離散褶積和并行計算的數(shù)值求解方案可以得到與直接求解法相當?shù)臄?shù)值結(jié)果，在保證計算精度的同時，內(nèi)存和CPU耗時則顯著降低，有效地克服了傳統(tǒng)數(shù)值算法在求解L-S方程時系數(shù)矩陣存儲和線性方程組求解效率低的問題.