汶川震源有限差分數(shù)值模擬研究

干大勇，趙晨露，張勁超，孟祥琦

(成都理工大學地球物理學院，四川 成都 610059)

汶川震源有限差分數(shù)值模擬研究

干大勇，趙晨露，張勁超，孟祥琦

(成都理工大學地球物理學院，四川 成都 610059)

2008年5月12日在四川汶川發(fā)生8.0級地震，震中位于龍門山斷裂帶上。地震激發(fā)的地震波可以用來研究震源破裂過程及地球內部的層次結構。對地震臺站接收到的地震記錄的擬合圖重新進行了擬合，得到一個極易用數(shù)學式表達的震源時間函數(shù)，并結合龍門山的實際情況，劃分為碎屑巖層、碳酸鹽層及基底層模擬區(qū)域，并運用交錯網(wǎng)格的有限差分法對聲波方程進行數(shù)值模擬。研究表明，利用3個不同初始震源都能得到較好的波場快照圖，而初始震源持續(xù)時間短的波場快照精度更高；不同初始震源不同地層的單點P波形圖都能與根據(jù)地震臺站接收記錄擬合的P波形圖相符合。因此，通過數(shù)值模擬能夠反映汶川地震波傳播的相關特征。

有限差分；聲波；數(shù)值模擬；汶川震源

有限差分法長期以來是求取微分方程數(shù)值解的一個主要方法，也是目前地震勘探資料數(shù)字處理中采用的求解波動方程的主要方法。波動方程有限差分法正演模擬，對認識地震波傳播規(guī)律、進行地震屬性研究、地震資料地質解釋、儲層評價等，均具有重要的理論和實際意義。以波動方程為基礎的有限差分法能夠準確的描述地震波在地下復雜介質中運動學和動力學特征，對地震波的波場特征和傳播規(guī)律能夠進行完整的描述。筆者主要以二維聲波方程為基礎，分析有限差分數(shù)值模擬中的頻散和邊界條件，對汶川地震激發(fā)的地震波傳播進行數(shù)值模擬。

1 有限差分聲波數(shù)值模擬

在二維空間域內的X-Z平面內，以X軸為水平軸，以Z軸為垂直軸，其二維二階波動方程可以表示為：

(1)

式中，c0為傳播速度，m/s；u為位移函數(shù)，m。

根據(jù)均勻彈性各項同性介質速度、應力、位移3個場變量以及場變量關系方程，不考慮剪應力的影響并令其為零，將二階方程轉變?yōu)橐浑A方程：

(2)

式中，Vx、Vz分別為質點振動速度的水平分量和垂直分量，m/s；ρ為介質密度,kg/m3；p為壓力,Pa。

1.5 評價標準 無壓紅:受壓部位皮膚與其他皮膚無變化;輕度壓紅:受壓部位皮膚有紅暈，解除壓力后30 min可以自行消退;瘀紅:受壓部位皮膚瘀紫，解除壓力后也無法消退，判定為壓瘡1期。

1.1 差分格式

采用經(jīng)典的二階差分格式對式(2)進行離散：

(3)

式中，k表示時間采樣間隔,s；i表示X軸采樣間隔,m；j表示Z軸采樣間隔,m；Δx和Δz分別表示空間步長,m；Δt表示時間步長,s。

1.2 差分網(wǎng)格頻散

不同頻率的波分量以不同相速度傳播時，會在聲波方程有限差分數(shù)值解中存在差分網(wǎng)格頻散。因此，在差分模擬中若對空間和時間步長取值不當會造成頻散。對一種數(shù)值解法,需要知道計算穩(wěn)定的離散參數(shù)區(qū)域,即分析解法的穩(wěn)定性，據(jù)此筆者采用Courant的穩(wěn)定性條件[3]:

1.3 完全匹配吸收邊界

對式(2)引入人工邊界條件，得完全匹配層方程為：

(4)

式中，dx、dz分別為x、z方向的衰減系數(shù)。

波長按傳播距離進行衰減時衰減速度很快。當衰減吸收系數(shù)隨空間位置變化時，不會在介質中產(chǎn)生任何反射。衰減吸收函數(shù)采用Collino[4]與Berenger J[5]等提出的通式：

(5)

式中，i代表從模型鑲邊后的外邊界到有效模型區(qū)域內邊界的網(wǎng)格點序號；PML代表完全匹配吸收層沿坐標軸方向所占的網(wǎng)格點個數(shù)；m代表吸收衰減系數(shù)的階數(shù)；B為衰減幅度因子。

圖1 完全匹配層吸收邊界

完全匹配層作為吸收邊界的基本原理是在研究區(qū)域內四周加入人工邊界。圖1所示為在加入的完全匹配層中分為3個區(qū)域(在區(qū)域1令dx=0,dz≠0；在區(qū)域2令dz=0,dx≠0； 在區(qū)域3令dx≠0,dz≠0)。當波傳播到完全匹配層的邊界時按距離的指數(shù)衰減，數(shù)值也越來越近似為零，這就不會產(chǎn)生邊界反射。

2 初始震源

聲波方程數(shù)值模擬初始震源設置為與汶川地震震源相類似的震源。圖2所示為汶川地震發(fā)生時某地震臺站接收到的地震記錄擬合的時間持續(xù)圖。

根據(jù)以上擬合的汶川地震震源時間圖形，通過計算，擬合出一條近似的曲線圖來模擬汶川震源(見圖3)。

3 數(shù)值計算與分析

圖2 汶川地震震源時間函數(shù)

數(shù)值模擬中盡量與汶川地震實際情況相符合，震源設置在13.5km處，模擬區(qū)域為24km×24km，對整個區(qū)域分為3個不同介質模型，分別為碎屑巖層0～5km、碳酸鹽巖層5～9.6km和基底9.6～24km(見圖4)。

模擬模型大小為600×600個網(wǎng)絡，網(wǎng)格間距空間Δx=Δy=40m，時間步長Δt=0.01s,分別用3種不同初始震源對聲波方程進行數(shù)值模擬。

3.1 利用聲波方程模擬不同時間的波場快照

1)以圖3(a)震源為初始震源 以圖3(a)震源為初始震源進行的聲波方程模擬時不同時間的波場快照如圖5所示。當t=2s時，地震波已從基底傳播到碳酸鹽巖層(見圖5(a))；t=4s時地震波已到達碎屑巖層，并且在碎屑巖與碳酸鹽巖的分界面存在界面反射(見圖5(b))；t=5.6s時地震波傳播到地表(見圖5(c))。由于加入了邊界細搜條件，人工邊界反射不是很明顯。

圖3 擬合近似的地震震源時間函數(shù)

圖4 數(shù)值模擬區(qū)域

2)以圖3(b)震源為初始震源 以圖3(b)震源為初始震源進行的聲波方程模擬時不同時間的波場快照如圖6所示。地震波由一介質傳入另一介質時都存在界面反射，由于地震波在6s內是連續(xù)傳播的，所以在短時間內界面反射現(xiàn)象不是很明顯。由于存在邊界反射，已經(jīng)影響了地震波傳播，因而在6s后會看到明顯的邊界反射效果(見圖6(a)和圖6(b))。在t=5.6s時已經(jīng)開始有明顯的邊界反射出現(xiàn)，并且地震波的前波已到達地表(見圖6(c))。當?shù)卣鸩ㄖ鞑ǖ竭_地表時，在碳酸鹽巖層和基底出現(xiàn)明顯的邊界反射(見圖6(d))。

圖5 以圖3(a)震源為初始震源進行的聲波方程模擬時不同時間的波場快照

圖6 以圖3(b)震源為初始震源進行的聲波方程模擬時不同時間的波場快照

圖7 以圖3(c)震源為初始震源進行的聲波方程模擬時不同時間的波場快照

圖8 根據(jù)地震臺站接收到的記錄擬合的P波形圖

3)以圖3(c)震源為初始震源 以圖3(c)震源為初始震源進行的聲波方程模擬時不同時間的波場快照如圖7所示。地震波前波到達地表(見圖7(a))，地震波主波已經(jīng)開始出現(xiàn)(見圖7(b))，第3次波出現(xiàn)(見圖7(c))，第4波出現(xiàn)(見圖7(d))。在對以上數(shù)值模擬的基礎上，選擇3個不同介質的點，給出其P波形圖，并與地震臺站接收到的記錄擬合的P波形圖(見圖8)進行比較。

3.2 震相模擬

1)以圖3(a)為初始震源 以圖3(a)為初始震源模擬得到的各點P波形圖如圖9所示。由圖9可知，當持續(xù)時間為1s時，由于時間短精度高，在短時內單點P波形圖效果較好，因而能更好地描述地震波的傳播特征。當波傳播到各界面出現(xiàn)邊界反射導致負值，隨著時間推移振動逐漸停止。

2)以圖3(b)為初始震源 以圖3(b)為初始震源模擬得到的各點P波形圖如圖10所示。由圖10可知，當持續(xù)時間為10s時，單點P波形圖與地震臺站接收到的地震記錄擬合的P波形圖最為相似,說明時間和空間的精度影響的頻散最小,從而更能反映地震波在地層的傳播形態(tài)，且隨著時間的推移該圖形會趨于平穩(wěn)。

3)以圖3(c)為初始震源 以圖3(c)為初始震源模擬得到的各點P波形圖如圖11所示。由圖11可知，當持續(xù)時間為100s時，雖然計算效率降低了，但更能結合實際情況，很好地反映了地震波的特征。

總之，通過計算模擬得到的P波形圖與地震臺站接收到的地震記錄擬合的P波形圖大致相似，說明上述方法能夠較準確地反映地震波在地下的反射、透射、散射等傳播形態(tài)和各種復雜的運動學和動力學特征。

4 結 語

根據(jù)3層地質模型，通過穩(wěn)定性條件與邊界條件的約束，利用3種不同持續(xù)時間的汶川初始震源進行聲波方程數(shù)值模擬，并將模擬結果與根據(jù)地震臺站接收的地震記錄擬合的P波形圖對比。研究表明，利用3個不同初始震源都能得到較好的波場快照圖，而初始震源持續(xù)時間短的波場快照精度更高；不同初始震源不同地層的單點P波形圖都能與根據(jù)地震臺站接收記錄擬合的P波形圖相符合。因此，通過聲波方程的有限差分數(shù)值模擬能夠反映地震波在地層的傳播特征。

圖9 以圖3(a)為初始震源模擬得到的各點P波形圖

圖10 以圖3(b)為初始震源模擬得到的各點P波形圖

圖11 以圖3(c)為初始震源得到的各點P波形圖

[1]王衛(wèi)民,趙連峰,李娟，等.四川汶川8.0級地震震源過程[J].地球物理學報，2008，51(5):1403-1410.

[2]嚴珍珍,張懷,楊長春，等.汶川大地震地震波傳播的譜元法數(shù)值模擬[J].中國科學，2009，39(4):393-402.

[3]Collina F,Tsogka. Application of the perfectly matched absorbing layer model to the linear elastodynamic problem in anisotropic heterogeneous media[J].Geophysics,2001,66(1):294-307.

[4]Zakaria A, Penrose J. The Two Dimensional Numerical Modeling Of Acoustic Wave Propagation in Shallow Wate[J].Austrialian Acoustical Society Conference,2000,15:1-5.

[5]Berenger J.P A perfectly matched layer for the absorption of electromagnetic waves[J].Journal of Computational Physics，1994,114(2):185-200.

[編輯] 李啟棟

10.3969/j.issn.1673-1409(N).2012.10.019

P315.3

A

1673-1409(2012)10-N060-05