一種利用Nystr?m離散與FFT快速褶積的散射地震波并行計(jì)算方法

徐楊楊， 孫建國(guó) ， 商耀達(dá)

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

0 引言

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

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

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

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

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

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

1 等價(jià)與歸一化等價(jià)L-S方程

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

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

(1)

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

(2)

(3)

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

2 等價(jià)L-S方程的Nystr?m離散

令Pn是插值算子，{rj,j=1,2,…,n}?Ω是插值基點(diǎ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的中心.此時(shí)，需要計(jì)算弱奇異積分：

(11)

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

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

Lippmann-Schwinger積分方程的常元Nystr?m近似生成的系數(shù)矩陣為稠密矩陣，并且每一個(gè)矩陣元都必須計(jì)算格林函數(shù)在剖分元上的弱奇異積分.當(dāng)散射體規(guī)模較大時(shí)，網(wǎng)格剖分?jǐn)?shù)量也很大，在二維情況下，系數(shù)矩陣由N2個(gè)元素組成(N=NxNz，Nx、Nz分別為x、z方向上的網(wǎng)格剖分?jǐn)?shù)量).采用高斯消元法或LU分解等直接法求解時(shí)需要存儲(chǔ)整個(gè)矩陣及其逆矩陣，內(nèi)存需求為O(N2)并且矩陣求逆計(jì)算復(fù)雜度為O(N3).而迭代法無(wú)需存儲(chǔ)逆矩陣，其復(fù)雜度下降為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ái)求解L-S方程離散生成的線性系統(tǒng).

(12)

(13)

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

(14)

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

(15a)

(15b)

其中：

(15c)

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

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

(16)

(17)

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

(18)

從而有:

(19)

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

(20)

其中:

0≤i≤Nx-1，

(21)

5 頻率域散射波場(chǎng)正演的并行化

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

5.1 基于MPI的頻率并行

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

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

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

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

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

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

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

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

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

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

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

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

6 數(shù)值試驗(yàn)

下面通過(guò)數(shù)值試驗(yàn)來(lái)測(cè)試基于Nystr?m法的IE-FFT并行算法在解決地震散射波場(chǎng)正演模擬中的有效性以及計(jì)算效率.首先設(shè)計(jì)了球形單體散射模型，通過(guò)與直接解法得到的數(shù)值結(jié)果進(jìn)行對(duì)比，分析IE-FFT算法在存儲(chǔ)和計(jì)算效率上的優(yōu)勢(shì).然后用更為復(fù)雜的SEG/EAGE的鹽丘模型對(duì)快速算法進(jìn)行進(jìn)一步測(cè)試.數(shù)值試驗(yàn)環(huán)境：Linux操作系統(tǒng)，包含11個(gè)節(jié)點(diǎn)，每個(gè)節(jié)點(diǎn)2個(gè)CPU，每個(gè)CPU包含8個(gè)核(CPU型號(hào)：Intel(R) Xeon(R) CPU E5-2620 0 @ 2.00 GHz )，使用基于MPI+OpenMP的并行算法對(duì)散射模型進(jìn)行單炮數(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，時(shí)間采樣間隔為4 ms，記錄時(shí)間2.4 s.

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

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

圖2 球體散射模型直接法與IE-FFT迭代法15Hz單頻波場(chǎng)數(shù)值結(jié)果(a) 直接法單頻波場(chǎng)實(shí)部； (b) 直接法單頻波場(chǎng)虛部； (c) IE-FFT迭代法單頻波場(chǎng)實(shí)部； (d) IE-FFT迭代法單頻波場(chǎng)虛部.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 球體散射模型直接法(實(shí)線)與IE-FFT迭代法(點(diǎn)劃線)單道對(duì)比(a) 第40道單道對(duì)比； (b) 第150道彈道對(duì)比.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時(shí)間與內(nèi)存消耗對(duì)比見(jiàn)表1.數(shù)值實(shí)現(xiàn)過(guò)程中，直接法耗費(fèi)CPU時(shí)間為4624.27 s，每個(gè)進(jìn)程內(nèi)存消耗為515 M，IE-FFT法迭代50次耗費(fèi)CPU計(jì)算時(shí)間為396.37 s，每個(gè)進(jìn)程內(nèi)存消耗為78 M.IE-FFT迭代算法的數(shù)值結(jié)果可以較好地吻合直接法結(jié)果，而計(jì)算速度是直接法的11倍，并且應(yīng)用格林函數(shù)的Toeplitz性質(zhì)降低了系數(shù)矩陣的內(nèi)存需求，內(nèi)存消耗僅約為直接法的1/7.

6.2 SEG/EAGE鹽丘模型

當(dāng)?shù)叵陆橘|(zhì)中包含較大尺度的復(fù)雜散射體且速度擾動(dòng)更強(qiáng)時(shí)，積分方程法應(yīng)用于此類介質(zhì)模型時(shí)需要在全空間中進(jìn)行離散計(jì)算.而直接法因其離散后的矩陣為NxNz階方陣，普通的微機(jī)難以獲取相應(yīng)的棧內(nèi)存來(lái)存儲(chǔ)系數(shù)矩陣且計(jì)算效率不滿足實(shí)際應(yīng)用需求.為了驗(yàn)證IE-FFT迭代算法在復(fù)雜強(qiáng)散射介質(zhì)中的適應(yīng)性，選取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，時(shí)間采樣間隔為4 ms，記錄時(shí)間4 s.

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

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

圖6 IE-FFT迭代算法15HZ單頻波場(chǎng)數(shù)值結(jié)果(a) 單頻波場(chǎng)實(shí)部； (b) 單頻波場(chǎng)虛部.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迭代算法散射波場(chǎng)數(shù)值結(jié)果Fig.7 Numerical result using the IE-FFT iterative method on the SEG/EAGE salt dome model

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

7 結(jié)論

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