周建科，印興耀，曹丹平

(中國(guó)石油大學(xué)(華東)地球科學(xué)與技術(shù)學(xué)院，山東青島266580)

波動(dòng)方程正演模擬對(duì)認(rèn)識(shí)地震波傳播規(guī)律和指導(dǎo)地震資料解釋等具有重要意義[1-2]，要想精確地模擬地震波在地下介質(zhì)中的傳播，一方面要采用較為接近實(shí)際地層的地球物理模型，另一方面要采用高精度的數(shù)值模擬算法。有限差分法(FDM)和有限元法(FEM)是地震波場(chǎng)正演模擬中非常重要的兩種方法，有限差分法具有編程簡(jiǎn)單，計(jì)算效率高等特點(diǎn)，在地震模擬中得到廣泛的研究和應(yīng)用，但難以處理起伏地表和自由邊界條件[3-4];有限元法能夠適應(yīng)劇烈的物性變化，自然地滿足自由邊界條件，因此具有精度高、可模擬任意復(fù)雜結(jié)構(gòu)、易于進(jìn)行邊界處理等優(yōu)點(diǎn)[5-8]，其缺點(diǎn)是計(jì)算量大和內(nèi)存占用高。

有限元法具有多種網(wǎng)格剖分方式，能夠?qū)θ我鈴?fù)雜的邊界進(jìn)行有效剖分，因而能夠準(zhǔn)確模擬起伏地形、復(fù)雜構(gòu)造和復(fù)雜介質(zhì)條件下的地震波場(chǎng)，種類繁多的插值函數(shù)可以提供不同精度的數(shù)值計(jì)算結(jié)果。當(dāng)單元采用線性函數(shù)進(jìn)行插值時(shí)，得到的數(shù)值結(jié)果往往精度較低，達(dá)不到預(yù)期目標(biāo)，當(dāng)采用高次函數(shù)進(jìn)行插值時(shí)，具有較高的數(shù)值模擬精度，但計(jì)算量以及內(nèi)存的占用也會(huì)增加。史明娟等[9]以及劉云等[10]采用雙二次插值有限元法(以下簡(jiǎn)稱雙二次插值法)實(shí)現(xiàn)了大地電磁的數(shù)值模擬，得到高精度的數(shù)值結(jié)果;戴前偉等[11]實(shí)現(xiàn)了雙二次插值的探地雷達(dá)有限元數(shù)值模擬，與雙線性插值有限元法(以下簡(jiǎn)稱雙線性插值法)相比，雙二次插值法的數(shù)值模擬結(jié)果更能突出異常體的響應(yīng)。

應(yīng)用有限元法解波傳播問(wèn)題時(shí)需占用巨量的機(jī)時(shí)和內(nèi)存，使得它在地球物理中的應(yīng)用受到極大限制。郭建[12]和劉有山等[8]根據(jù)結(jié)構(gòu)剛度矩陣稀疏性的特點(diǎn)，采用稀疏矩陣的壓縮存儲(chǔ)行(Compressed Sparse Row，CSR)格式。在此不僅考慮了結(jié)構(gòu)剛度矩陣的稀疏性，還考慮其對(duì)稱性，采用緊湊存儲(chǔ)法只存儲(chǔ)結(jié)構(gòu)剛度矩陣下三角部分的非零元素，同時(shí)，在波場(chǎng)的遞推過(guò)程中，只有結(jié)構(gòu)剛度矩陣每一行的非零元素與對(duì)應(yīng)的節(jié)點(diǎn)位移進(jìn)行相乘;根據(jù)質(zhì)量守恒原則，采用對(duì)角的集中質(zhì)量矩陣代替一致質(zhì)量矩陣，避免顯式有限元法在每個(gè)時(shí)間步長(zhǎng)上進(jìn)行大型稀疏矩陣的求逆運(yùn)算。

我們?cè)诰匦尉W(wǎng)格剖分情況下，采用雙二次插值法實(shí)現(xiàn)了二維聲波方程的數(shù)值模擬，采用對(duì)角的集中質(zhì)量矩陣和緊湊存儲(chǔ)法提高計(jì)算效率、減少內(nèi)存的占用，為實(shí)現(xiàn)有限元法的高精度數(shù)值模擬提供一定的指導(dǎo)。

1 方法原理

1.1 有限元法求解二維聲波方程

使用伽勒金法求解二維標(biāo)量波動(dòng)方程得到有限元方程組[13]

(1)

式中：U(t)代表t時(shí)刻結(jié)構(gòu)節(jié)點(diǎn)位移列向量;M和K分別代表結(jié)構(gòu)的質(zhì)量矩陣和剛度矩陣：

(2)

式中：Me和Ke分別為單元的質(zhì)量矩陣和剛度矩陣;N為結(jié)構(gòu)中單元的個(gè)數(shù);c代表單元內(nèi)介質(zhì)的聲波速度;N為形函數(shù)矩陣;∑代表單元質(zhì)量(剛度)矩陣組裝為結(jié)構(gòu)質(zhì)量(剛度)矩陣的法則。

對(duì)于(1)式，時(shí)間因子采用二階中心差分格式進(jìn)行求解，并考慮震源的加載，得到

(3)

式中：F(t)為t時(shí)刻作用在節(jié)點(diǎn)上的等效載荷;Δt為時(shí)間采樣間隔。

進(jìn)行波場(chǎng)的遞推時(shí)，需要給定初始條件：

(4)

1.2 單元內(nèi)雙二次插值

首先將求解區(qū)域Ω剖分成矩形單元，在單元內(nèi)進(jìn)行雙二次插值。在每個(gè)單元上取8個(gè)節(jié)點(diǎn)：4個(gè)角點(diǎn)和4條邊的中點(diǎn)，8個(gè)節(jié)點(diǎn)的編號(hào)次序及坐標(biāo)如圖1所示。母單元與子單元間的坐標(biāo)變換關(guān)系為

(5)

其中，x0,y0是子單元中點(diǎn)的坐標(biāo)，a與b分別為子單元的長(zhǎng)和寬。

采用雙二次插值，母單元內(nèi)任意一點(diǎn)的位移u可以表示為

(6)

式中：a1,…,a8為待定常數(shù);ξ,η為母單元上的坐標(biāo)。

將母單元8個(gè)節(jié)點(diǎn)的坐標(biāo)以及相應(yīng)的位移代入(6)式，整理可得

(7)

圖1 單元剖分示意a 子單元; b 母單元

(8)

將(5)式、(8)式代入(2)式，采用等參單元和高斯數(shù)值積分[14]，得到單元?jiǎng)偠染仃嘖e和一致質(zhì)量矩陣Me。

兩矩陣計(jì)算式分別為

1.3 集中質(zhì)量矩陣

在地震波場(chǎng)數(shù)值模擬中，一般采用對(duì)角的集中質(zhì)量矩陣代替一致質(zhì)量矩陣，以避免時(shí)間循環(huán)過(guò)程中的求逆運(yùn)算。對(duì)線性插值單元集中質(zhì)量矩陣的求取，一般將一致質(zhì)量矩陣每一行所有元素疊加到對(duì)角線上，然后令非對(duì)角線元素為零。但對(duì)高次插值單元而言，這樣的做法不可取，例如將(10)式中的矩陣的每一行元素疊加到對(duì)角線上后，4個(gè)角節(jié)點(diǎn)分配到的質(zhì)量為負(fù)值，這在力學(xué)上是不合理的，同時(shí)也將對(duì)計(jì)算精度產(chǎn)生不良影響[15]。

對(duì)于矩形雙二次插值單元，采用質(zhì)量守恒原理來(lái)確定單元集中質(zhì)量矩陣。設(shè)每個(gè)角節(jié)點(diǎn)分配到的單元質(zhì)量為m1，每條邊的中點(diǎn)分配到的單元質(zhì)量為m2，單元的總質(zhì)量為m。定義質(zhì)量分布參數(shù)β滿足

(11)

根據(jù)質(zhì)量守恒定理，有

(12)

將(12)式代入(11)式整理得到：

(13)

(14)

矩形單元雙二次插值的集中質(zhì)量矩陣的形成不是唯一的，例如以單元一致質(zhì)量矩陣的對(duì)角線元素為權(quán)重時(shí)有β=3/76=0.0395，此外取β=1/36=0.0278在實(shí)際中也運(yùn)用得較多[15]。本文進(jìn)行數(shù)值模擬時(shí)，β的取值為3/76，即β=0.0395。

1.4 結(jié)構(gòu)剛度矩陣的緊湊存儲(chǔ)法

根據(jù)結(jié)構(gòu)剛度矩陣的集成法則(附錄A)可知，結(jié)構(gòu)剛度矩陣每一行只有少量的非零元素，且具有對(duì)稱性，每個(gè)節(jié)點(diǎn)的位移只受到與其共點(diǎn)單元節(jié)點(diǎn)的影響。圖2是計(jì)算區(qū)域中的任一共點(diǎn)單元，可見(jiàn)，對(duì)于非邊界角節(jié)點(diǎn)(如節(jié)點(diǎn)L3)而言，對(duì)其位移產(chǎn)生影響的節(jié)點(diǎn)只有21個(gè)(包括該節(jié)點(diǎn)本身)，而對(duì)非邊界中點(diǎn)的位移產(chǎn)生影響的節(jié)點(diǎn)只有13個(gè)，與邊界節(jié)點(diǎn)位移相關(guān)的節(jié)點(diǎn)則更少。因此由雙二次插值法得到的結(jié)構(gòu)剛度矩陣每一行非零元素最多為21個(gè)，再根據(jù)其對(duì)稱性，只需存儲(chǔ)結(jié)構(gòu)剛度矩陣下三角部分的非零元素即可。根據(jù)以上分析，本文采用緊湊存儲(chǔ)法來(lái)存儲(chǔ)結(jié)構(gòu)剛度矩陣，需要3個(gè)一維數(shù)組：浮點(diǎn)型數(shù)組CS按行存儲(chǔ)結(jié)構(gòu)剛度矩陣下三角部分的非零元素，整型數(shù)組GA存放CS中對(duì)應(yīng)元素的列號(hào)，整型數(shù)組ID存放下三角部分每一行非零元素的個(gè)數(shù)。

以水平方向有nx個(gè)單元，垂直方向有nz個(gè)單元的均勻剖分結(jié)構(gòu)為例，結(jié)構(gòu)節(jié)點(diǎn)與單元編號(hào)規(guī)則為：先從左到右，然后由上至下，單元節(jié)點(diǎn)的局部編號(hào)如圖1所示。下面以節(jié)點(diǎn)L3為例來(lái)說(shuō)明緊湊存儲(chǔ)法的原理。節(jié)點(diǎn)L3決定結(jié)構(gòu)剛度矩陣的第L3行，該行非零元素的列號(hào)為與其共點(diǎn)單元節(jié)點(diǎn)的結(jié)構(gòu)編號(hào)，再根據(jù)結(jié)構(gòu)剛度矩陣的對(duì)稱性，第L3行所需存儲(chǔ)的元素僅為11個(gè)，其列號(hào)分別為：L1-2，L1-1，L1，L1+1，L1+2，L2-1，L2，L2+1，L3-2，L3-1，L3，這11個(gè)數(shù)用數(shù)組GA來(lái)儲(chǔ)存，以便后續(xù)使用。一個(gè)具有nx×nz個(gè)單元的結(jié)構(gòu)，結(jié)構(gòu)剛度矩陣中元素個(gè)數(shù)為(3nxnz+2nx+2nz+1)2，而采用緊湊存儲(chǔ)法后，結(jié)構(gòu)剛度矩陣中需要存儲(chǔ)元素的個(gè)數(shù)僅為：25nxnz+5nx+5nz+1，輔助數(shù)組GA與ID中元素的個(gè)數(shù)分別為25nxnz+5nx+5nz+1與3nxnz+2nx+2nz+1。

圖2 共點(diǎn)單元示意

由(3)式可知，在時(shí)間循環(huán)過(guò)程中，涉及到結(jié)構(gòu)剛度矩陣與結(jié)構(gòu)節(jié)點(diǎn)位移列向量的矩陣乘法，如果結(jié)構(gòu)剛度矩陣中的全部元素參與運(yùn)算，這將非常耗時(shí)。眾所周知，矩陣中的零元素對(duì)矩陣乘法是沒(méi)有貢獻(xiàn)的，利用數(shù)組ID和GA可以找到非零元素在結(jié)構(gòu)剛度矩陣中的具體位置，因此在時(shí)間循環(huán)過(guò)程中，只需結(jié)構(gòu)剛度矩陣每一行的非零元素與對(duì)應(yīng)的節(jié)點(diǎn)位移進(jìn)行相乘。變帶寬存儲(chǔ)法[14]雖然也只存儲(chǔ)了結(jié)構(gòu)剛度矩陣下三角部分的元素，但其把半帶寬中的零元素也儲(chǔ)存了，同時(shí)在計(jì)算中，半帶寬中的零元素也參與了運(yùn)算，對(duì)大型結(jié)構(gòu)不再適用(節(jié)點(diǎn)數(shù)越多，同一單元節(jié)點(diǎn)編號(hào)的最大值與最小值之差也就越大，半帶越寬，半帶寬中的零元素也就越多)。

在內(nèi)存為4GB，處理器為Intel(R)Core(TM)2 DUO E8400 @3.00GHz的計(jì)算機(jī)上對(duì)比分析了緊湊存儲(chǔ)法與變帶寬儲(chǔ)存法的計(jì)算效率以及所占內(nèi)存，具體結(jié)果見(jiàn)表1。當(dāng)模型網(wǎng)格數(shù)為100×100時(shí)，變帶寬存儲(chǔ)法進(jìn)行300次時(shí)間循環(huán)的耗時(shí)和結(jié)構(gòu)剛度矩陣內(nèi)存占用分別是緊湊存儲(chǔ)法的19.20倍和14.47倍，當(dāng)網(wǎng)格數(shù)為400×400時(shí)，相應(yīng)的倍數(shù)變?yōu)?1.90和56.90。

表1 不同存儲(chǔ)方法的CPU耗時(shí)以及內(nèi)存占用統(tǒng)計(jì)

2 模型試算

采用C語(yǔ)言編制雙二次插值法二維聲波方程地震波場(chǎng)數(shù)值模擬程序，運(yùn)行在內(nèi)存4GB，處理器Intel(R)Core(TM)2 DUO E8400 @3.00GHz的計(jì)算機(jī)上。震源選用定點(diǎn)上延遲的雷克子波，峰值頻率fp=40Hz。

2.1 計(jì)算效率比較

為了說(shuō)明雙二次插值法的優(yōu)缺點(diǎn)(計(jì)算效率及內(nèi)存占用情況)，分別采用雙線性插值法和雙二次插值法進(jìn)行地震波場(chǎng)的數(shù)值模擬，在無(wú)可見(jiàn)數(shù)值頻散的情況下，比較二者的計(jì)算耗時(shí)以及存儲(chǔ)需求量。在本算例中選取縱波波速c=1900m/s的均勻介質(zhì)模型，計(jì)算區(qū)域Ω={(x,z)|0≤x≤2000m,0≤z≤2000m}，時(shí)間采樣間隔Δt=0.4ms，震源位于模型中央。采用2種算法得到的T=0.4s時(shí)刻波場(chǎng)快照如圖3所示。其中，圖3a至圖3c是在不同尺寸單元網(wǎng)格下采用雙線性插值法得到的波場(chǎng)快照，當(dāng)單元網(wǎng)格大小為2.5m×2.5m(圖3a)和2.0m×2.0m(圖3b)時(shí)，波場(chǎng)快照可見(jiàn)微弱的數(shù)值頻散，當(dāng)單元網(wǎng)格減小到1.7m×1.7m時(shí)(圖3c)無(wú)可見(jiàn)數(shù)值頻散，此時(shí)結(jié)構(gòu)中的網(wǎng)格數(shù)為1176×1176;圖3d至圖3f則為在不同尺寸單元網(wǎng)格下采用雙二次插值法得到的波場(chǎng)快照，在6.0m×6.0m(圖3d)以及5.5m×5.5m(圖3e)的單元網(wǎng)格下，波場(chǎng)快照具有微弱的數(shù)值頻散，當(dāng)單元網(wǎng)格減小到5.0m×5.0m時(shí)，數(shù)值頻散得以消除，此時(shí)網(wǎng)格數(shù)為400×400。也就是說(shuō)，當(dāng)震源主頻為40Hz、介質(zhì)速度為1900m/s時(shí)，在保證無(wú)可見(jiàn)數(shù)值頻散情況下，雙線性插值法的單元網(wǎng)格最大尺寸為1.7m×1.7m，雙二次插值法的單元網(wǎng)格最大尺寸為5.0m×5.0m。

表2給出了本算例在無(wú)可見(jiàn)數(shù)值頻散下，采用緊湊存儲(chǔ)法時(shí)2種數(shù)值模擬方法的CPU耗時(shí)(1001次時(shí)間循環(huán))以及內(nèi)存占用情況。可見(jiàn)，為使兩者的數(shù)值模擬結(jié)果均無(wú)可見(jiàn)的數(shù)值頻散，雙線性插值法占用的內(nèi)存是雙二次插值法的一倍，且雙二次插值法的單步耗時(shí)也比雙線性插值法少。在這里值得指出的是，雙二次插值法的穩(wěn)定性條件比雙線性插值法更為苛刻，因此雙線性插值法的時(shí)間采樣間隔可以取得比雙二次插值法大，當(dāng)兩者計(jì)算時(shí)間相同時(shí)，雙線性插值法可以選取更大的時(shí)間步長(zhǎng)，時(shí)間循環(huán)次數(shù)也就變少，有可能使得計(jì)算耗時(shí)比雙二次插值法少，因此本文得出的結(jié)論是雙二次插值法的單步耗時(shí)比雙線性插值法少。

表2 無(wú)可見(jiàn)數(shù)值頻散情況下兩種數(shù)值模擬方法的CPU耗時(shí)以及內(nèi)存占用量

圖3 均勻介質(zhì)模型T=0.4s時(shí)刻波場(chǎng)快照a 雙線性插值法,單元網(wǎng)格大小2.5m×2.5m; b 雙線性插值法,單元網(wǎng)格大小2m×2m; c 雙線性插值法,單元網(wǎng)格大小1.7m×1.7m; d 雙二次插值法,單元網(wǎng)格大小6.0m×6.0m; e 雙二次插值法,單元網(wǎng)格大小5.5m×5.5m; f 雙二次插值法,單元網(wǎng)格大小5.0m×5.0m

2.2 復(fù)雜構(gòu)造模型

為了進(jìn)一步驗(yàn)證雙二次插值法對(duì)復(fù)雜構(gòu)造的模擬精度，設(shè)計(jì)了如圖4所示的模型，模型包括一個(gè)垂直斷層和一個(gè)背斜構(gòu)造，長(zhǎng)為3600m，深2000m。模型中的階梯狀界面上有4個(gè)拐點(diǎn)，其中點(diǎn)1與頂界面以及點(diǎn)2的垂直距離均為400m，點(diǎn)2與點(diǎn)3和4的水平距離分別為400,2000m，點(diǎn)2到其下水平界面的垂直距離為400m，點(diǎn)4與右邊界的水平距離為800m，背斜弧高500m，由上至下速度依次為2000，3000，4000m/s。震源位于(1800m，200m)處，四周邊界均采用改進(jìn)的Sarma吸收邊界條件[16]，厚度為400m。單元網(wǎng)格大小為4m×4m，考慮到穩(wěn)定性條件，時(shí)間采樣間隔取0.2ms，計(jì)算時(shí)長(zhǎng)T=1.6s。

圖5給出了0.5，0.6，0.7，0.8，0.9，1.0s時(shí)刻的波場(chǎng)快照(單元中心點(diǎn)的波場(chǎng)值由單元8個(gè)節(jié)點(diǎn)的波場(chǎng)值插值得到)。從圖5中可以清楚地看到地震波在4個(gè)拐點(diǎn)處產(chǎn)生的繞射波、界面的反射波以及繞射波形成的多次反射波，說(shuō)明雙二次插值法在模擬復(fù)雜構(gòu)造時(shí)具有較好的計(jì)算精度。單炮記錄如圖6所示，由于離散背斜界面采用矩形網(wǎng)格，因此可以觀察到由階梯狀界面產(chǎn)生的非物理繞射波。

在該數(shù)值算例中，實(shí)際參與計(jì)算(包含邊界條件)的網(wǎng)格數(shù)為1100×700，進(jìn)行8001次時(shí)間循環(huán)耗時(shí)4045s，平均單步耗時(shí)0.5056s，采用緊湊存儲(chǔ)法存儲(chǔ)結(jié)構(gòu)剛度矩陣占用的內(nèi)存為155.76MB，結(jié)構(gòu)集中質(zhì)量矩陣占用的內(nèi)存為8.83MB。

圖4 復(fù)雜構(gòu)造速度模型

圖5 復(fù)雜構(gòu)造模型不同時(shí)刻的波場(chǎng)快照a 0.5s; b 0.6s; c 0.7s; d 0.8s; e 0.9s; f 1.0s

圖6 復(fù)雜構(gòu)造模型單炮記錄

3 結(jié)束語(yǔ)

我們采用雙二次插值法實(shí)現(xiàn)了二維聲波方程的有限元法數(shù)值模擬，采用質(zhì)量守恒定律得到對(duì)角的集中質(zhì)量矩陣，結(jié)構(gòu)剛度矩陣采用緊湊存儲(chǔ)。模型試算結(jié)果表明，雙二次插值法能夠提高地震波場(chǎng)數(shù)值模擬的精度，采用集中質(zhì)量矩陣和剛度矩陣的緊湊存儲(chǔ)后，內(nèi)存消耗以及計(jì)算速度能夠滿足實(shí)際生產(chǎn)要求。雙二次插值法和雙線性插值法的對(duì)比實(shí)驗(yàn)表明，在無(wú)可見(jiàn)數(shù)值頻散情況下，雙二次插值法消耗的內(nèi)存更少，單步耗時(shí)更短?；陔p二次插值法的以上優(yōu)越性，該算法可望進(jìn)一步推廣到彈性波數(shù)值模擬以及偏移成像等領(lǐng)域。

附錄A 雙二次插值法結(jié)構(gòu)剛度矩陣的組裝

根據(jù)單元?jiǎng)偠染仃囅禂?shù)與結(jié)構(gòu)剛度矩陣系數(shù)的對(duì)應(yīng)關(guān)系(具體對(duì)應(yīng)關(guān)系請(qǐng)參考文獻(xiàn)[12])，結(jié)構(gòu)剛度矩陣第L3(圖2)行需要進(jìn)行疊加的元素為

參 考 文 獻(xiàn)

[1] 吳國(guó)忱,王華忠.波場(chǎng)模擬中的數(shù)值頻散分析與校正策略[J].地球物理學(xué)進(jìn)展,2005,20(1):58-65

Wu G C,Wang H Z.Analysis of numerical dispersion in wave-field simulation[J].Progress in Geophysics,2005,20(1):58-65

[2] 孫林潔,印興耀.基于PML邊界條件的高倍可變網(wǎng)格有限差分?jǐn)?shù)值模擬方法[J].地球物理學(xué)報(bào),2011,54(6):1614-1623

Sun L J,Yin X Y.A finite-difference scheme based on PML boundary condition with high power grid step variation[J].Chinese Journal of Geophysics,2011,54(6):1614-1623

[3] Alford R M,Kelly K R,Boore D M.Accuracy of finite-difference modeling of the acoustic wave equation[J].Geophysics,1974,39(6):834-842

[4] 殷文,印興耀,吳國(guó)忱,等.高精度頻率域彈性波方程有限差分方法及波場(chǎng)模擬[J].地球物理學(xué)報(bào),2006,49(2):561-568

Yin W,Yin X Y,Wu G C,et al.The method of finite difference of high precision elastic wave equations in the frequency domain and wave-field simulation[J].Chinese Journal of Geophysics,2006,49(2):561-568

[5] Marfurt K J.Accuracy of finite-difference and finite-element modeling of the scalar and elastic wave equations[J]. Geophysics,1984,49(5):533-549

[6] 薛東川,王尚旭,焦淑靜.起伏地表復(fù)雜介質(zhì)波動(dòng)方程有限元數(shù)值模擬方法[J].地球物理學(xué)進(jìn)展,2007,22(2):522-529

Xue D C,Wang S X,Jiao S J.Wave equation finite-element modeling including rugged topography and complicated medium[J]. Progress in Geophysics,2007,22(2):522-529

[7] 薛東川，王尚旭.波動(dòng)方程有限元疊前逆時(shí)偏移[J].石油地球物理勘探,2008,43(1):17-21

Xue D C,Wang S X.Wave-equation finite element prestack reverse-time migration[J].Oil Geophysical Prospecting,2008,43(1):17-21

[8] 劉有山,滕吉文,劉少林,等.稀疏存儲(chǔ)的顯式有限元三角網(wǎng)格地震波數(shù)值模擬及其PML吸收邊界條件[J].地球物理學(xué)報(bào),2013,56(9):3085-3099

Liu Y S，Teng J W，Liu S L，et al.Explicit finite element method with triangle meshes stored by sparse format and its perfectly matched layers absorbing boundary condition[J].Chinese Journal of Geophysics,2013,56(9):3085-3099

[9] 史明娟,徐世浙,劉斌.大地電磁二次函數(shù)插值的有限元法正演模擬[J].地球物理學(xué)報(bào),1997,40(3):421-430

Shi M J,Xu S Z,Liu B.Finite element method using quadratic element in MT forward modeling[J]. Chinese Journal of Geophysics,1997,40(3):421-430

[10] 劉云,王緒本.大地電磁二維自適應(yīng)地形有限元正演模擬[J].地震地質(zhì),2010,32(3):382-391

Liu Y,Wang X B.FEM using adaptive topography in 2D MT forward modeling[J].Seismology and Geology,2010,32(3):382-391

[11] 戴前偉,王洪華,馮德山,等.基于雙二次插值的探地雷達(dá)有限元數(shù)值模擬[J].地球物理學(xué)進(jìn)展,2012,27(2):736-743

Dai Q W,Wang H H,Fei D S,et al.Finite element numerical simulation for GPR based on quadratic interpolation[J].Progress in Geophysics,2012,27(2):

736-743

[12] 郭建.一種有限元快速算法[J].石油物探,1991,30(2):36-43

Guo J.A kind of fast finite element algorithm[J].Geophysical Prospecting for Petroleum,1991,30(2):36-43

[13] 杜世通.變速不均勻介質(zhì)中波動(dòng)方程的有限元法數(shù)值解[J].華東石油學(xué)院學(xué)報(bào),1982,6(2):1-20

Du S T.Finite element numerical solution of wave propagation in non-homogeneous medium with variable velocities[J].Journal of East China Petroleum Institute,1982,6(2):1-20

[14] 徐世浙.地球物理中的有限元法[M].北京:科技出版社,1994:1-308

Xu S Z.Finite Element method for geophysics[M].Beijing:Science Press,1994:1-308

[15] 王勖成.有限單元法[M].北京:清華大學(xué)出版社,2003:472-475

Wang M C.Finite element method[M].Beijing:Tsinghua University Press,2003:472-475

[16] 王月英,宋建國(guó).波場(chǎng)正演模擬中Sarma邊界條件的改進(jìn)[J].石油物探,2007,46(4):359-362

Wang Y Y，Song J G.Improvement of Sarma boundary condition in wavefield forward modeling[J].Geophysical Prospecting for Petroleum,2007,46(4):359-362