基于低降速帶復雜地形的瑞雷波數(shù)值模擬研究

王道遠，王鑫，李 勤

西安石油大學地球科學與工程學院陜西省油氣成藏地質(zhì)學重點實驗室，陜西西安710065

0 引言

在地震勘探過程中，隨著丘陵、盆地等近地表地區(qū)的深層次研究，發(fā)現(xiàn)近地表地震-地質(zhì)條件下接收到的地震波場成分十分復雜，地震資料的精細處理面臨著巨大挑戰(zhàn)，數(shù)值模擬研究越發(fā)必要。在反射地震勘探中經(jīng)常會接收到一類低頻率、強能量的面波干擾，通常作為地震勘探中的干擾波被壓制掉。但是面波中含有豐富的信息，與近地表低降速帶有直接的關(guān)系。瑞雷面波勘探是一種有效獲取近地表地層結(jié)構(gòu)的途徑，可以有效利用反射地震資料中的瑞雷面波，使其變廢為寶，節(jié)省一筆非常可觀的勘探費用[1]。近年來，瑞雷面波勘探研究正在進行廣泛開展。

地震數(shù)值模擬是研究地震波在復雜介質(zhì)中傳播的有效方法之一[2]。Rayleigh首次提出了瑞雷面波理論，并于1887年證明了瑞雷波的存在[3]。Haskell在最早時提出瑞雷面波頻散曲線正演Haskell算法，此算法易出現(xiàn)高頻數(shù)值不穩(wěn)定現(xiàn)象，而且通常計算速度偏慢[4]。為解決所面臨的問題，眾多學者從數(shù)值角度開展研究，簡化頻散方程的表現(xiàn)形式，解決了頻散方程的不穩(wěn)定性導出了新的頻散方程，有效避免了高頻數(shù)值的不穩(wěn)定，同時計算速度提升[5-8]。

為了通過地震波場更好地分析瑞雷面波在近地表彈性介質(zhì)中的表現(xiàn)特征、響應(yīng)機理及傳播規(guī)律，國內(nèi)學者進行了一系列研究工作。周竹生等進行了基于彈性介質(zhì)中瑞利波有限差分數(shù)值模擬，采用交錯網(wǎng)格有限差分法對多種彈性介質(zhì)模型進行了正演模擬，得到相應(yīng)的波場快照和地震記錄[9]。熊章強等進行了吸收邊界及自由邊界條件研究，正演模擬得到震源在自由界面時的瑞利面波波場，并解釋了波場特征[10]。張大洲等通過對三層介質(zhì)數(shù)值模擬得到波場記錄，揭示了瑞雷波傳播特征[11]。鄒延延等模擬了簡單夾層和巖溶介質(zhì)模型地震波動力特征[12]。眾多學者進行了復雜近地表的地震波波場模擬，得到了復雜近地表波場響應(yīng)特征[13-15]。邵廣周等對比、分析瑞利波的頻散特征和各模式的實際情況，同時也了解了瑞利波各模式在波場中的能量分布情況[16]。蘇傳行等進行了起伏地形瑞雷波模擬研究[17]。

雖然上述方法對瑞雷面波進行了有效的模擬，但對于近地表低降速帶中的波場特征，尤其是瑞雷面波的特征分析不夠全面?；诟唠A交錯網(wǎng)格有限差分原理，對不同復雜近地表地質(zhì)模型進行了數(shù)值模擬，分別從合成地震記錄和波場快照兩個方面分析了不同條件下瑞雷面波的響應(yīng)特征及其傳播規(guī)律。在此基礎(chǔ)上，對面波頻散及數(shù)值模擬邊界條件等問題進行了綜合、全面的分析，從而為實際工程勘探中近地表低降速帶區(qū)域含有軟弱夾層和異常體的情況提供了參考。瑞雷面波正演模擬有助于更準確、更全面地了解瑞雷波的傳播特征，為實際工程勘探提供理論依據(jù)。

1 波動方程數(shù)值模擬基本原理

1.1 各向同性介質(zhì)一階速度-應(yīng)力彈性波方程

在二維情況下，一階速度-應(yīng)力彈性波方程如下所示[18]：

式中，θxx、θzz、θxz為應(yīng)力的三個分量；vx表示速度的水平分量，vz表示速度的垂直分量；λ、μ表示拉梅常數(shù)，λ=ρ（vp2-2vs2），μ=ρvs2，vp表示介質(zhì)中的縱波速度，vs表示橫波速度；ρ表示介質(zhì)密度。

1.2 交錯網(wǎng)格高階有限差分格式

1.2.1 時間上的2M階有限差分近似

通過利用交錯網(wǎng)格求取一階應(yīng)力—速度方程時，速度分量在t+Δt/2時刻進行計算，應(yīng)力分量在t時刻進行計算。將各分量進行Taylor公式展開化簡得：

（2）式和（3）式相減可得2階時間差分近似公式：

式中，Δt為時間步長。

1.2.2 空間上的2N階有限差分近似

空間上采取2N階有限差分近似，使用Taylor公式展開計算得：

式中：為待定的差分系數(shù)，Δx為空間網(wǎng)格間距。對于待定系數(shù)，根據(jù)Taylor公式的2N階展開，有：

可以求得高階空間差分系數(shù)。

1.3 震源設(shè)置

震源函數(shù)可以使用高斯一階導函數(shù)子波、雷克子波等表示，在進行全波場數(shù)值模擬時選取不同的震源函數(shù)對結(jié)果影響很大[19]。

選取雷克子波，頻率fm=30Hz為震源函數(shù)。

1.4 穩(wěn)定性條件

穩(wěn)定性條件在進行數(shù)值計算時較為關(guān)鍵，進行合理的參數(shù)選擇，可以避免出現(xiàn)嚴重頻散、結(jié)果溢出等現(xiàn)象，得出較好的模擬結(jié)果。采用的穩(wěn)定性條件為：

1.5 邊界條件

在面波的數(shù)值模擬過程中，首先要考慮自由邊界情況，面波沿自由表面?zhèn)鞑?，自由邊界條件必然會對面波的傳播有著一定影響。其次是對吸收邊界條件的探討，在地震勘探當中，地震波在一個非均勻各項異性半空間進行傳播，而在計算機模擬時，不可給予如此大的計算區(qū)域，計算是有限的，這時需要引進吸收邊界條件，從而減小吸收邊界條件對正演模擬的影響。

采用理想匹配層（PML）吸收邊界條件處理人工邊界反射，通過在待計算區(qū)域的左右兩側(cè)分別加一個吸收邊界層，將所有的入射波包括縱波和橫波都進行吸收，不產(chǎn)生反射，這樣在邊界處有較好的吸收效果。在近地表研究各種波時，不得不面臨淺層對地震波的影響，地震波激發(fā)后，地震波迅速到達淺層的反射界面，隨后反彈至自由表面，然而在自由表面的處理相對不易，產(chǎn)生向下的反射波的同時，也將會產(chǎn)生不少的轉(zhuǎn)換波。對此本文采用聲學與彈性介質(zhì)邊界方法（AEA）。AEA方法操作簡便，正演模擬效果良好，其二維處理如下[20]：

λ、θzz、ρ表示自由邊界以上的拉梅常數(shù)、正應(yīng)力及密度。圖1為不同吸收邊界厚度的波場快照對比。

圖1 不同吸收邊界厚度的波場快照對比，0層PML（a），10層PML（b），30層PML（c）Fig.1 Comparison of wave field snapshot with different absorption boundary thickness: 0 layer of PML （a）, 10 layers of PML （b）, 30 layers of PML （c）

2 模型試算

為了模擬近地表低降速帶區(qū)域的波場特征，本文設(shè)計了幾種地質(zhì)模型，通過增加模型的復雜程度，模擬近似真實條件下近地表及低降速帶區(qū)域的瑞雷面波波場特征。

2.1 近地表低降速層介質(zhì)模型

圖2為簡化的近地表低降速層介質(zhì)模型，模擬真實情況下地下低速層對地震勘探的影響。設(shè)計模型大小為200 m×200 m，主頻采用30 Hz雷克子波，各層介質(zhì)的參數(shù)如表1所示，模型x方向網(wǎng)格步長為0.5 m，z方向網(wǎng)格步長為0.5 m，炮點坐標選取為（100，0），模擬時間為100 ms，得到的全波場快照如圖2所示。

表1 近地表低降速層介質(zhì)模型參數(shù)Table 1 Model parameters of near surface low velocity layer

圖2（a）為模擬得到的近地表低降速層介質(zhì)模型波場快照，圖2（b）為模擬得到的地震記錄。在圖2中，R表示面波，S表示橫波，P表示縱波，PP表示反射縱波，SS表示反射橫波。從圖2（b）中可以看到各波同相軸，P波斜率較小，傳播速度較快，S波次之，而R波傳播速度較慢。從圖2（a）中可以看出，由于每層介質(zhì)參數(shù)不同，地震波在不同層交界面形成反射，各層的分層情況清晰可見。入射縱波和入射橫波在穿過第一層界面后形成反射波后，能量隨之減弱，形成的透射縱波和透射橫波繼續(xù)穿過第二層界面，形成二次反射波和透射波，由于第二層和第三層分界面波阻抗較強，以致深層地震波能量偏弱。面波在自由表面由縱橫波互相間疊加干涉形成，在地表傳播速度較低，能量極強且遠大于P波。

2.2 傾斜斷層介質(zhì)模型

圖3為簡化的近地表斜面斷層模型，在工程地質(zhì)調(diào)查中常見，其特點是由于受到地質(zhì)作用，兩側(cè)巖塊沿斷裂面發(fā)生位移而導致上下面介質(zhì)參數(shù)不同。設(shè)計模型大小為200 m×200 m，主頻采用30 Hz雷克子波，上層介質(zhì)S1參數(shù)采用：Vp1=830 m/s，Vs1=490 m/s，ρ1=1.4 g/cm3，下層介質(zhì)S2參數(shù)采用：Vp2=1 400 m/s，Vs2=770 m/s，ρ2=2.0 g/cm3，模型x方向網(wǎng)格步長為0.5 m，z方向網(wǎng)格步長為0.5 m，炮點坐標選取為（100，0），模擬時間為110 ms，得到的全波場快照如圖4所示。

圖3 傾斜斷層介質(zhì)模型示意圖Fig.3 Schematic diagram of inclined fault medium model

圖4（a）為模擬得到的傾斜斷層介質(zhì)模型波場快照，圖4（b）為模擬得到的地震記錄。在圖4中，R表示面波，S表示橫波，P表示縱波，PP表示反射縱波，SS表示反射橫波。從圖4（a）中可以明顯看出傾斜斷層面的存在，地震波傳播到斷裂面處時方向發(fā)生偏折，產(chǎn)生反射波和透射波，由于兩側(cè)速度差較大，形成一強波阻抗界面，淺層面波能量最強，深層地 震波能量較少，瑞雷面波在遇到斷裂面后部分能量被反彈，部分能量發(fā)生透射。從圖4（b）中可以看出同相軸不以炮點為對稱軸左右對稱，由于地震波在通過截面時速度減小，同相軸產(chǎn)生一定的角度。

圖4 傾斜斷層介質(zhì)模型波場快照（a）及地震記錄（b）Fig.4 Wave field snapshot （a） and seismic record （b） of inclined fault medium model

2.3 垂直地質(zhì)異常體介質(zhì)模型

在地質(zhì)調(diào)查中，常見較大尺度的地質(zhì)異常構(gòu)造現(xiàn)象。圖5為簡化的垂直地質(zhì)異常體介質(zhì)模型，設(shè)計模型大小為200 m×200 m，主頻采用30Hz雷克子波，上層介質(zhì)S1參數(shù)采用：Vp1=760 m/s，Vs1=450 m/s，ρ1=1.5 g/cm3，下層介質(zhì)S2參數(shù)采用：Vp2=1 600 m/s，Vs2=900 m/s，ρ2=2.1 g/cm3，模型x方向網(wǎng)格步長為0.5 m，z方向網(wǎng)格步長為0.5 m，炮點坐標選?。?00，0），模擬時間為100 ms，得到的全波場快照如圖6。

圖5 垂直地質(zhì)異常體介質(zhì)模型示意圖Fig.5 Schematic diagram of vertical geological anomaly mode

圖6 垂直地質(zhì)異常體介質(zhì)模型波場快照（a）及地震記錄（b）Fig.6 Snapshot of wave field （a） and seismic record （b） of vertical geological anomaly model

圖6（a）為模擬得到的垂直地質(zhì)異常體介質(zhì)模型波場快照，圖6（b）為模擬得到的地震記錄。在圖6中，R表示面波，S表示橫波，P表示縱波，R’表示反射面波。從圖6（a）中可以看出，在地震波傳播過程中，左右波形未與炮點為中心左右對稱，左側(cè)出現(xiàn)較多能量團聚集，并局限在一個區(qū)域內(nèi)，可以確認為一處地質(zhì)異常體。在異常體附近分布反射波以及拐點處產(chǎn)生的繞射波，各種波互相疊加干涉，震相較難準確識別。瑞雷面波在通過地質(zhì)異常體時，部分能量發(fā)生反射，形成反射瑞雷面波，瑞雷面波在異常體內(nèi)震蕩反射，在遇到接觸面后小部分能量發(fā)生透射，大部分能量繼續(xù)進行反射和向外圍散發(fā)，瑞雷面波在通過地質(zhì)異常體后，能量發(fā)生較大衰減，傳播距離落后于均勻介質(zhì)。在圖6（b）中可以看出，同相軸不以炮點為對稱軸左右對稱，發(fā)生折疊現(xiàn)象，能量在垂直異常體內(nèi)震蕩反射聚集，在兩側(cè)形成反射波和透射波。

3 結(jié)論

通過交錯網(wǎng)格有限差分法對近地表低降速層介質(zhì)模型、傾斜斷層介質(zhì)模型及垂直地質(zhì)異常體介質(zhì)模型進行了正演模擬研究，高階差分格式可以明顯削弱數(shù)值頻散的影響。通過對模擬結(jié)果分析發(fā)現(xiàn)，近地表低降速帶區(qū)域相較工程層面尺度較大，地質(zhì)情況復雜，地震波場豐富。在對近地表低降速層介質(zhì)模型模擬時，分層情況清晰可見，波形圖完整，而真實近地表地質(zhì)復雜，傾斜斷層及地質(zhì)異常體廣泛分布，對波場能量分布影響較大。瑞雷面波在通過復雜地形時，會發(fā)生反射、透射現(xiàn)象，能量重新分配，能量團聚集，在地震記錄上出現(xiàn)同相軸的多次偏折。因此，利用對低降速帶復雜地形的數(shù)值模擬，分析了近地表構(gòu)造對波場的影響以及瑞雷面波的響應(yīng)機理和傳播規(guī)律，結(jié)果對實際地質(zhì)勘探具有指導意義。