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

王道遠，王鑫，李 勤

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

0 引言

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

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

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

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

1 波動方程數值模擬基本原理

1.1 各向同性介質一階速度-應力彈性波方程

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

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

1.2 交錯網格高階有限差分格式

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

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

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

式中，Δt為時間步長。

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

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

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

可以求得高階空間差分系數。

1.3 震源設置

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

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

1.4 穩定性條件

穩定性條件在進行數值計算時較為關鍵，進行合理的參數選擇，可以避免出現嚴重頻散、結果溢出等現象，得出較好的模擬結果。采用的穩定性條件為：

1.5 邊界條件

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

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

λ、θzz、ρ表示自由邊界以上的拉梅常數、正應力及密度。圖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 模型試算

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

2.1 近地表低降速層介質模型

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

表1 近地表低降速層介質模型參數Table 1 Model parameters of near surface low velocity layer

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

2.2 傾斜斷層介質模型

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

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

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

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

2.3 垂直地質異常體介質模型

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

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

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

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

3 結論

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