基于雙二次插值的地震波場有限元法數值模擬

周建科，印興耀，曹丹平

(中國石油大學(華東)地球科學與技術學院，山東青島266580)

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

有限元法具有多種網格剖分方式，能夠對任意復雜的邊界進行有效剖分，因而能夠準確模擬起伏地形、復雜構造和復雜介質條件下的地震波場，種類繁多的插值函數可以提供不同精度的數值計算結果。當單元采用線性函數進行插值時，得到的數值結果往往精度較低，達不到預期目標，當采用高次函數進行插值時，具有較高的數值模擬精度，但計算量以及內存的占用也會增加。史明娟等[9]以及劉云等[10]采用雙二次插值有限元法(以下簡稱雙二次插值法)實現了大地電磁的數值模擬，得到高精度的數值結果;戴前偉等[11]實現了雙二次插值的探地雷達有限元數值模擬，與雙線性插值有限元法(以下簡稱雙線性插值法)相比，雙二次插值法的數值模擬結果更能突出異常體的響應。

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

我們在矩形網格剖分情況下，采用雙二次插值法實現了二維聲波方程的數值模擬，采用對角的集中質量矩陣和緊湊存儲法提高計算效率、減少內存的占用，為實現有限元法的高精度數值模擬提供一定的指導。

1 方法原理

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

使用伽勒金法求解二維標量波動方程得到有限元方程組[13]

(1)

式中：U(t)代表t時刻結構節(jié)點位移列向量;M和K分別代表結構的質量矩陣和剛度矩陣：

(2)

式中：Me和Ke分別為單元的質量矩陣和剛度矩陣;N為結構中單元的個數;c代表單元內介質的聲波速度;N為形函數矩陣;∑代表單元質量(剛度)矩陣組裝為結構質量(剛度)矩陣的法則。

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

(3)

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

進行波場的遞推時，需要給定初始條件：

(4)

1.2 單元內雙二次插值

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

(5)

其中，x0,y0是子單元中點的坐標，a與b分別為子單元的長和寬。

采用雙二次插值，母單元內任意一點的位移u可以表示為

(6)

式中：a1,…,a8為待定常數;ξ,η為母單元上的坐標。

將母單元8個節(jié)點的坐標以及相應的位移代入(6)式，整理可得

(7)

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

(8)

將(5)式、(8)式代入(2)式，采用等參單元和高斯數值積分[14]，得到單元剛度矩陣Ke和一致質量矩陣Me。

兩矩陣計算式分別為

1.3 集中質量矩陣

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

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

(11)

根據質量守恒定理，有

(12)

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

(13)

(14)

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

1.4 結構剛度矩陣的緊湊存儲法

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

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

圖2 共點單元示意

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

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

表1 不同存儲方法的CPU耗時以及內存占用統(tǒng)計

2 模型試算

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

2.1 計算效率比較

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

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

表2 無可見數值頻散情況下兩種數值模擬方法的CPU耗時以及內存占用量

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

2.2 復雜構造模型

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

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

在該數值算例中，實際參與計算(包含邊界條件)的網格數為1100×700，進行8001次時間循環(huán)耗時4045s，平均單步耗時0.5056s，采用緊湊存儲法存儲結構剛度矩陣占用的內存為155.76MB，結構集中質量矩陣占用的內存為8.83MB。

圖4 復雜構造速度模型

圖5 復雜構造模型不同時刻的波場快照a 0.5s; b 0.6s; c 0.7s; d 0.8s; e 0.9s; f 1.0s

圖6 復雜構造模型單炮記錄

3 結束語

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

附錄A 雙二次插值法結構剛度矩陣的組裝

根據單元剛度矩陣系數與結構剛度矩陣系數的對應關系(具體對應關系請參考文獻[12])，結構剛度矩陣第L3(圖2)行需要進行疊加的元素為

參 考 文 獻

[1] 吳國忱,王華忠.波場模擬中的數值頻散分析與校正策略[J].地球物理學進展,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邊界條件的高倍可變網格有限差分數值模擬方法[J].地球物理學報,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] 殷文,印興耀,吳國忱,等.高精度頻率域彈性波方程有限差分方法及波場模擬[J].地球物理學報,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] 薛東川,王尚旭,焦淑靜.起伏地表復雜介質波動方程有限元數值模擬方法[J].地球物理學進展,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] 薛東川，王尚旭.波動方程有限元疊前逆時偏移[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] 劉有山,滕吉文,劉少林,等.稀疏存儲的顯式有限元三角網格地震波數值模擬及其PML吸收邊界條件[J].地球物理學報,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] 史明娟,徐世浙,劉斌.大地電磁二次函數插值的有限元法正演模擬[J].地球物理學報,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] 劉云,王緒本.大地電磁二維自適應地形有限元正演模擬[J].地震地質,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] 戴前偉,王洪華,馮德山,等.基于雙二次插值的探地雷達有限元數值模擬[J].地球物理學進展,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] 杜世通.變速不均勻介質中波動方程的有限元法數值解[J].華東石油學院學報,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].北京:清華大學出版社,2003:472-475

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

[16] 王月英,宋建國.波場正演模擬中Sarma邊界條件的改進[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