TI介質中P-SV波傳播的SBP-SAT模擬

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

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

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

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

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

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

1 波動方程SBP-SAT方法構建

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

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

圖1 常見TI介質模型示意圖

(1)

式中:ρ為彈性體密度;σij為應力;vx和vy分別為質點速度在x軸和y軸上的分量。為表述方便僅考慮繞z軸旋轉情形,設極化角為θ,剛度矩陣元素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)

傳統有限差分法在模擬曲線域時大量角點會影響模擬結果的穩定性。將曲線坐標系映射到直角坐標系并劃定合適的網格形式能更好模擬層狀地質結構。

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

圖2 坐標變換

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

(9)

為實現這一映射關系需要使用雅可比矩陣行列式和鏈式法則

(10)

曲線域的參數向量方向導數為

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

(11)

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

Ut=AξUξ+AηUη

(12)

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

(13)

拓展SBP算子得

(14)

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

(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框架來說,設置人工邊界的各種邊界性質只需要選用相應的反射系數即可。

p=Rq

(21)

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

1.3 系統能量分析

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

(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均為非正值,說明系統的能量有界,證明矩陣對稱型彈性波動方程框架滿足穩定性的要求。彈性波動方程的矩陣對稱形式使能量分析過程簡化,并建立起了通用框架。

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

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

(25)

其中,

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

2 算例分析

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

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

圖3 多塊模型示意圖

建立的模型尺寸1 980 m×1 980 m,有3×3個子計算域,每個計算域的網格數為200×200,采用4階精度傳統SBP算子,時間步長Δt取1 ms。定義模型的彈性系數分別為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處。

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

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

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

2.2 分層VTI介質模型SBP-SAT方法模擬

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

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

表1 雙層VTI介質模型參數表格

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

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

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

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

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

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

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

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

圖6 多塊模型示意圖

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

圖7 x方向速度幅值

圖8 y方向速度幅值

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

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

3 準確性驗證和計算效率

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

(27)

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

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

圖9 x方向速度幅值

圖10 x方向的節點幅值差

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

4 結 論

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