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

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

汶川震源有限差分?jǐn)?shù)值模擬研究

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

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

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

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

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

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

在二維空間域內(nèi)的X-Z平面內(nèi)，以X軸為水平軸，以Z軸為垂直軸，其二維二階波動(dòng)方程可以表示為：

(1)

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

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

(2)

式中，Vx、Vz分別為質(zhì)點(diǎn)振動(dòng)速度的水平分量和垂直分量，m/s；ρ為介質(zhì)密度,kg/m3；p為壓力,Pa。

1.5 評(píng)價(jià)標(biāo)準(zhǔn) 無(wú)壓紅:受壓部位皮膚與其他皮膚無(wú)變化;輕度壓紅:受壓部位皮膚有紅暈，解除壓力后30 min可以自行消退;瘀紅:受壓部位皮膚瘀紫，解除壓力后也無(wú)法消退，判定為壓瘡1期。

1.1 差分格式

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

(3)

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

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

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

1.3 完全匹配吸收邊界

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

(4)

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

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

(5)

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

圖1 完全匹配層吸收邊界

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

2 初始震源

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

根據(jù)以上擬合的汶川地震震源時(shí)間圖形，通過(guò)計(jì)算，擬合出一條近似的曲線圖來(lái)模擬汶川震源(見(jiàn)圖3)。

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

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

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

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

3.1 利用聲波方程模擬不同時(shí)間的波場(chǎng)快照

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

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

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

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

圖5 以圖3(a)震源為初始震源進(jìn)行的聲波方程模擬時(shí)不同時(shí)間的波場(chǎng)快照

圖6 以圖3(b)震源為初始震源進(jìn)行的聲波方程模擬時(shí)不同時(shí)間的波場(chǎng)快照

圖7 以圖3(c)震源為初始震源進(jìn)行的聲波方程模擬時(shí)不同時(shí)間的波場(chǎng)快照

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

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

3.2 震相模擬

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

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

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

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

4 結(jié) 語(yǔ)

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

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

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

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

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

[2]嚴(yán)珍珍,張懷,楊長(zhǎng)春，等.汶川大地震地震波傳播的譜元法數(shù)值模擬[J].中國(guó)科學(xué)，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