楊在林, 魏夕杰, 孫 鋮, 楊 勇

(1.哈爾濱工程大學 航天與建筑工程學院,哈爾濱 150001; 2.哈爾濱工程大學 先進船舶材料與力學工業(yè)與信息化部重點實驗室,哈爾濱 150001; 3.廣西民族大學 建筑工程學院,南寧 530006)

地震模擬方法作為正演過程中的有力工具和策略,也是地震勘探過程中十分重要的一環(huán),它立足于不同介質(zhì)的參數(shù)特性,通過數(shù)值計算的方式來模擬波在介質(zhì)中的傳播過程,最終可得到波的傳播規(guī)律。使用各向異性模型表示具有空間變化的斷裂方向和密度、應力分布、復雜層理等的巖石具有特別重要的意義。波動方程數(shù)值直接解法以網(wǎng)格化的模型作為基礎,將復雜的地形地質(zhì)條件轉(zhuǎn)化成可以計算求解的形式。它包含了地震波傳播過程中的很多信息,因此得到的地震記錄往往比較全面真實,模擬精度相對較高。地震波動方程的數(shù)值解法主要包括:有限差分法、有限元法、邊界元法、譜元法、偽譜法等。

自Alterman等[1]開始使用其進行地震波場模擬至今,有限差分法已發(fā)展成為地震數(shù)值模擬領(lǐng)域最廣泛采用的數(shù)值解法之一,具有計算速度快、層狀介質(zhì)模擬能力好、編程簡便等優(yōu)點。有限差分方法也存在需要繼續(xù)完善改進的地方,比如在模擬波在復雜形狀模型邊界和內(nèi)部傳播時易產(chǎn)生穩(wěn)定性問題,時空域長期穩(wěn)定保證難度較大。這些都是我們期望優(yōu)化的方向?；谑褂媚芰糠ㄗC明穩(wěn)定性的目的,Kreiss等[2]提出了一階導數(shù)的SBP算子基本理論;Carpenter等[3]推導出具有最小模板的二階SBP算子形式,其內(nèi)部算子具有與一階導數(shù)內(nèi)部算子相同的節(jié)點數(shù)量。Mattsson[4]對SBP算子形式進行修改,使其更加簡潔。在處理帶邊界或界面問題時若使用強施加邊界條件,將會破壞SBP算子的性質(zhì)。得益于Funaro等[5]使用懲戒項的啟示,Carpenter等[6]提出了通過SAT方法弱施加邊界條件。Mattsson[7]系統(tǒng)地對比了不同方法施加邊界條件后對離散系統(tǒng)穩(wěn)定性的影響,確定采用SAT方法的離散系統(tǒng)更具優(yōu)勢,SBP-SAT離散框架由此確立起來。國內(nèi)對該方法逐漸重視,Sun等[8-10]在SBP-SAT方法離散框架的基礎上搭配完全匹配層作為人工邊界,有效壓制了邊界反射波的產(chǎn)生,獲得了穩(wěn)定的模擬結(jié)果。楊在林等[11]推導出了彈性波動方程SBP-SAT離散方程,進一步探究了該方法在波動領(lǐng)域的應用。

經(jīng)典反射地震勘探通常要求地表水平且地下介質(zhì)水平層狀,但隨著勘探領(lǐng)域向更復雜地質(zhì)區(qū)域深入,這些假設條件變得不再適用[12-16]。Igel等[17]將交錯網(wǎng)格有限差分法應用于各向異性介質(zhì),并且實現(xiàn)了柱坐標系和球坐標系下的數(shù)值模擬。不同數(shù)值解方法在橫向各向同性(transversely isotropy,TI)介質(zhì)運動學屬性研究上各有特色[18-20],對于SBP-SAT方法,在曲線坐標系下對角范數(shù)建立的SBP算子能夠滿足分部求和的性質(zhì),進一步的研究由此展開[21]。Carpenter等[22]和Nordstr?m等[23]將SAT方法拓展到曲線坐標下線性問題的多種邊界條件及塊狀界面條件。Virta等[24]建立了適用于復雜形狀與非均勻介質(zhì)的聲波SBP-SAT框架,可用于模擬介質(zhì)不連續(xù)模型。孫鋮[25]建立了曲線網(wǎng)格矩陣對稱型波動方程的多塊SBP-SAT方法,并建立了聲波SMF體系。SBP-SAT方法與適當?shù)奈諏踊蜻吔鐥l件結(jié)合并在曲線域中模擬波的傳播,公式推導和程序編譯得到簡化,仿真過程穩(wěn)定性強,應用前景非常廣泛。

本文基于速度-應力形式的彈性波動方程,采用SBP-SAT方法進行了矩陣對稱形式的波動方程離散,使用能量法證明其穩(wěn)定性。對幾種TI介質(zhì)模型進行數(shù)值模擬,得到波場信息來分析其中P-SV波的傳播規(guī)律,并對比COMSOL結(jié)果驗證SBP-SAT方法的正確性。

1 波動方程SBP-SAT方法構(gòu)建

1.1 彈性波動方程的矩陣對稱化

常見TI介質(zhì)模型示意圖,如圖1所示。傾斜橫向各向同性介質(zhì)速度-應力形式的P-SV波波動方程(體力為零)為

圖1 常見TI介質(zhì)模型示意圖

(1)

式中:ρ為彈性體密度;σij為應力;vx和vy分別為質(zhì)點速度在x軸和y軸上的分量。為表述方便僅考慮繞z軸旋轉(zhuǎn)情形,設極化角為θ,剛度矩陣元素dij為

d11=(C11cos2θ+C13sin2θ)cos2θ+

(C13cos2θ+C33sin2θ)sin2θ+C44sin2(2θ)

d13=(C11cos2θ+C13sin2θ)sin2θ+

(C13cos2θ+C33sin2θ)cos2θ-C44sin2(2θ)

d15=0.5(C11cos2θ+C13sin2θ)sin(2θ)-

0.5(C13cos2θ+C33sin2θ)sin(2θ)-

C44sin(2θ)(cos2θ-sin2θ)

d33=(C11sin2θ+C13cos2θ)sin2θ+

(C13sin2θ+C33cos2θ)cos2θ+C44sin2(2θ)

d35=0.5(C11sin2θ+C13cos2θ)sin(2θ)-

0.5(C13sin2θ+C33cos2θ)sin(2θ)+

C44sin(2θ)(cos2θ-sin2θ)

d55=0.25[C11sin(2θ)-C13sin(2θ)]sin(2θ)-

0.25[C13sin(2θ)-C33sin2θ)]sin(2θ)+

C44(cos2θ-sin2θ)2

(2)

將波動方程表示為矩陣形式

(3)

(4)

M=PTΛ-1/2P,M-1=PTΛ1/2P

(5)

代入(3)中得到

(6)

進一步,得

(7)

至此我們得到了SMF的彈性波動方程

(8)

傳統(tǒng)有限差分法在模擬曲線域時大量角點會影響模擬結(jié)果的穩(wěn)定性。將曲線坐標系映射到直角坐標系并劃定合適的網(wǎng)格形式能更好模擬層狀地質(zhì)結(jié)構(gòu)。

將曲線域Ω∈(x,y)轉(zhuǎn)換成直角域Ω′∈(ξ,η),如圖2所示。坐標映射關(guān)系為

圖2 坐標變換

[x(ξ,η),y(ξ,η)]?[ξ(x,y),η(x,y)]

(9)

為實現(xiàn)這一映射關(guān)系需要使用雅可比矩陣行列式和鏈式法則

(10)

曲線域的參數(shù)向量方向?qū)?shù)為

JUx=Uξyη+UηyξJUy=Uξxη+Uηxξ

(11)

利用式(10)可以得到彈性波動方程式(8)在曲線域的矩陣對稱形式為

Ut=AξUξ+AηUη

(12)

由于Ax和Ay為對稱矩陣,所以Aξ和Aη也為對稱矩陣。參數(shù)矩陣為

(13)

拓展SBP算子得

(14)

Iξ和Iη分別為(Nξ+1)階矩陣和(Nη+1)階矩陣,拓展系數(shù)矩陣得

(15)

邊界向量為

EW=e1ξ?Iη?I5EE=erξ?Iη?I5ES=Iξ?elη?I5EN=I1ξ?erη?I5

(16)

求解域邊界定義為W,E,S,N,得到SBP-SAT方法建立的矩陣對稱型彈性波動方程的離散形式為

Ut=AξDξU+AηDηU+SATW+SATE+SATS-SATN

(17)

式中,SATW,SATE,SATS,SATN的表達式分別為

(18)

1.2 適定的邊界條件與半離散近似

通過特征邊界條件的推導能夠得到傳入和傳出求解域的物理量。波動方程變形得到

(19)

(20)

對于SMF框架來說,設置人工邊界的各種邊界性質(zhì)只需要選用相應的反射系數(shù)即可。

p=Rq

(21)

反射系數(shù)的定義域為[-1,1],取值為0時即無反射邊界條件,取值為1時即自由邊界條件。傳入和傳出的特征變量關(guān)系可以定義不同的邊界條件,也意味著更大的適用范圍。邊界條件的通用框架建立起來。

1.3 系統(tǒng)能量分析

結(jié)合1.2節(jié)中得到的邊界條件式(20),對彈性波動方程在連續(xù)域進行能量分析。左乘UT后與其轉(zhuǎn)置相加,并在求解域上進行積分,得到

(22)

式(22)等價于

(23)

邊界項BTW,BTE,BTS,BTN分別為

BTW=-?W[UT(Ax)U]dSBTE=?E[UT(Ax)U]dSBTS=-?S[UT(Ay)U]dSBTN=?N[UT(Ay)U]dS

(24)

邊界條件為無反射邊界條件時,BTW,BTE,BTS,BTN均為非正值,說明系統(tǒng)的能量有界,證明矩陣對稱型彈性波動方程框架滿足穩(wěn)定性的要求。彈性波動方程的矩陣對稱形式使能量分析過程簡化,并建立起了通用框架。

離散化的彈性波動方程的證明過程與上述過程類似,當滿足邊界條件式(20)時,式(18)恒等于0。

將式(17)左乘UTHξη再與該式子轉(zhuǎn)置相加后,得

(25)

其中,

邊界項值均為非正,說明矩陣對稱型彈性波動方程的離散框架滿足穩(wěn)定性要求。

2 算例分析

2.1 傾斜橫向各向同性介質(zhì)模型SBP-SAT方法模擬

傾斜橫向各向同性(tilted transverse isotropic ,TTI)介質(zhì)是具有傾斜角度對稱軸的TI介質(zhì)形式,其剛度矩陣可對垂直橫向各向同性(transverse isotropy with a vertical axis of symmetry,VTI)介質(zhì)彈性系數(shù)矩陣進行Bond變換而得到,如圖3所示。

圖3 多塊模型示意圖

建立的模型尺寸1 980 m×1 980 m,有3×3個子計算域,每個計算域的網(wǎng)格數(shù)為200×200,采用4階精度傳統(tǒng)SBP算子,時間步長Δt取1 ms。定義模型的彈性系數(shù)分別為C11=28.00 GPa,C13=8.00 GPa,C33=15.00 GPa,C44=3.00 GPa,極化角取45°。初始條件設置為

g(x,y)=exp{-α[(x-x0)2+(y-y0)2]}

(26)

α影響波前面,取0.002目的是使波場效果更清晰明顯。單炮記錄位置選擇在深度825 m處。

獲得的結(jié)果包括波場快照和單炮記錄,在圖4的波場中存在波和波。同時波場沒有明顯的數(shù)值頻散現(xiàn)象,展現(xiàn)出SBP-SAT方法數(shù)值模擬的保真性和精準性。

圖4 極化角45°時的單炮記錄和x,y方向速度幅值

在TTI介質(zhì)中由于波的偏振方向與傳播方向存在夾角,除特定傳播方向外不存在純橫波或純縱波,需要使用克里斯托弗方程求解偏振向量,相關(guān)波場分離可以作為下一步的研究方向。分析圖4可知,介質(zhì)中波的傳播速度受到介質(zhì)參數(shù)影響不能保持相同,波場快照中外側(cè)的為波,內(nèi)側(cè)的為波,并且波的波前面相較于算例2.2對應位置更圓,波波前面角度變化并出現(xiàn)波面尖角現(xiàn)象。橫縱波場是耦合在一起的。

2.2 分層VTI介質(zhì)模型SBP-SAT方法模擬

模型的尺寸、子計算域數(shù)量和網(wǎng)格數(shù)與2.1節(jié)中相同,采用6階精度傳統(tǒng)SBP算子,時間步長取1 ms。多塊模型分上下兩層,上層0～1 320 m,對應子計算域2、3、5、6、8、9;下層1 320～1 980 m,對應子計算域1、4、7。極化角均取0°,構(gòu)成分層的VTI介質(zhì)模型。

初始條件和單炮記錄位置均不變。

表1 雙層VTI介質(zhì)模型參數(shù)表格

本算例主要為分析在分層界面處P-SV波的轉(zhuǎn)換和傳播規(guī)律。觀察圖5可以看到在同一介質(zhì)內(nèi)部,波前面清晰連貫,但在介質(zhì)分界處有明顯的波前面斷裂和速度差異,且上下兩部分形成明顯完整的反射波、透射波和轉(zhuǎn)換波。

圖5 單炮記錄和x,y方向速度幅值

為了震相的簡單標注,在此規(guī)定:或為波或波;角標1、2分別為上、下兩層;上標為上行波;下標為下行波。分析圖5可以得到數(shù)字:1為直達qP波;2為直達qSV波;3qP1qp1為反射波;4qP1qS1為反射轉(zhuǎn)換波;5qSV1qSV1為反射波;6qSV1qSV2為透射波;7qP1qSV2為透射轉(zhuǎn)換波;8qP1qP2為透射波。

從圖5可知,隨時間增大,波前面也在增大,波和波耦合在一起傳播,只在某些特殊方向上獨立傳播。當波遇到波阻抗分界面時會發(fā)生轉(zhuǎn)換,波不僅會產(chǎn)生反射波和透射波,還會產(chǎn)生反射轉(zhuǎn)換波和透射轉(zhuǎn)換波。波也有同樣的現(xiàn)象,其波場變化比較明顯,由于介質(zhì)的速度各向異性導致波的速度隨傳播方向變化而變化,產(chǎn)生了波面尖角。

另外波前面到達邊界后沒有明顯的反射,證明無反射邊界的效果良好。

2.3 含有裂縫的VTI介質(zhì)模型SBP-SAT方法模擬

不同地質(zhì)年代沉積形成的地層可近似為TI介質(zhì),而其中裂縫、斷層是普遍存在的,使用數(shù)值模擬方法對含裂縫的VTI介質(zhì)開展地震動模擬研究具有豐富的工程背景。

建立的模型尺寸、子計算域數(shù)量和每個計算域的網(wǎng)格數(shù)均不變,采用6階精度傳統(tǒng)SBP算子,時間步長Δt取1 ms。裂縫位于子計算域4～5,裂縫曲線設置為正弦函數(shù),子計算域5為曲線域,如圖6所示。定義彈性系數(shù)取C11=28 GPa,C13=8 GPa,C33=15 GPa,C44=3 GPa,極化角取0°。

圖6 多塊模型示意圖

初始條件施加在方向上,其他條件不變,如圖7、圖8所示。

圖7 x方向速度幅值

圖8 y方向速度幅值

SAT弱施加的特點使得邊界設置富有多樣性,在裂縫、空腔等類型邊界條件建立上,不必采用介質(zhì)填充等相對繁瑣的方式,而是可以通過設置反射參數(shù)這種簡單的方式來完成。本算例中子計算域4和5邊界的反射系數(shù)均設置為1,即設為自由邊界。觀察圖7,曲線域中波的傳播連續(xù)清晰,波場效果穩(wěn)定,同時裂縫邊界的反射作用也很明顯。

進一步,對相鄰兩子計算域網(wǎng)格、節(jié)點位置不對應的情形,可以采取子計算域各自建立更優(yōu)的節(jié)點分布的策略,這大大拓寬了方法的應用范圍。

3 準確性驗證和計算效率

為求簡明,建立均勻各向同性介質(zhì)模型,尺寸設置為3 000 m×3 000 m,橫波波速為2 000 m/s,縱波波速為1 100 m/s,密度為1 100 kg/m3。初始條件為

(27)

采用4階精度傳統(tǒng)SBP算子,時間步長為0.001 s,子計算域網(wǎng)格劃分為200×200。通過節(jié)點(1 500 m,1 000 m)的x方向求解速度幅值,與COMSOL軟件使用間斷伽遼金法計算的結(jié)果進行比較,觀察幅值差隨時間變化。

如圖9和圖10所示,結(jié)果吻合度較好,能夠驗證結(jié)果的準確性。

圖9 x方向速度幅值

圖10 x方向的節(jié)點幅值差

采用三種網(wǎng)格大小來計算CPU時間,模型采用2.1節(jié)的模型,網(wǎng)格數(shù)為150×150,300×300,600×600,迭代500次,分別得到的CPU計算時間為6.564 s,21.368 s,110.048 s。

4 結(jié) 論

本文使用SBP-SAT方法,對不同類型的橫向各向同性介質(zhì)模型進行了波場數(shù)值模擬研究。結(jié)果表明,SBP-SAT方法表現(xiàn)出了很好的適用性。在模擬分層介質(zhì)和含裂縫的情形時,既能夠保證計算精度又能夠保證穩(wěn)定性,進一步優(yōu)化子計算域劃分方式下有繼續(xù)提高計算效率的潛力。SBP-SAT方法同樣適用于其他介質(zhì)中波動方程的求解,推導過程與本文討論的方法相近。