黃建平 劉英輝 李偉 張盟勃 王揚州 楊永紅

摘要 ：深層地?zé)豳Y源是一種可再生的、蘊藏巨大能量的清潔能源，但其地球物理響應(yīng)特征不明確，中深層地?zé)崽綔y成功率較低。為了研究地震波在深層地?zé)釒r體中的傳播規(guī)律與波場特征，建立兩種深層地?zé)峄◢弾r模型，并使用等效交錯網(wǎng)格有限差分法實現(xiàn)聲波與彈性波的數(shù)值模擬。結(jié)果表明：地?zé)峄◢弾r體速度在溫度的影響下要高于圍巖的，會產(chǎn)生高速屏蔽現(xiàn)象，使得透射波能量變?nèi)?，限制了地?zé)釒r體下部地震波傳播能量；相比于聲波，彈性波具有更豐富的波場信息，波型與能量轉(zhuǎn)換使得彈性波地震記錄也比聲波地震記錄復(fù)雜。

關(guān)鍵詞 ：深層地?zé)釒r； 地震波傳播機制； 數(shù)值模擬； 有限差分模擬

中圖分類號 ：P 631.4 ???文獻標(biāo)志碼 ：A

引用格式 ：黃建平，劉英輝，李偉，等.深層地?zé)峄◢弾r體地震波數(shù)值模擬及傳播機制［J］.中國石油大學(xué)學(xué)報（自然科學(xué)版），2024，48（1）：63-69.

HUANG Jianping， LIU Yinghui， LI Wei， et al. Numerical simulation of seismic wave in deep geothermal rock mass［J］. Journal of China University of Petroleum （Edition of Natural Science），2024，48（1）：63-69.

Numerical simulation of seismic wave in deep geothermal rock mass

HUANG Jianping 1， LIU Yinghui 1， LI Wei 2， ZHANG Mengbo 3， WANG Yangzhou 4， YANG Yonghong 5

（1.School of Geosciences in China University of Petroleum（East China）， Qingdao 266580， China;

2.Shandong Energy Group， Jinan 250101， China;

3.PetroChina Changqing Oilfield， Xi an 710018， China;

4.Shandong Energy Group South America Company Limited， Qingdao 266000， China;

5.SINOPEC Shengli Oilfield Exploration and Development Research Institute， Dongying 257029， China）

Abstract ： Deep geothermal resources represent a substantial form of renewable and clean energy.However， their geophysical response characteristics remain poorly understood， leading to a low success rate in exploring these resources. This study aims to ?investigate the propagation theory and wavefield characteristics of seismic waves within deep geothermal rock masses. To achieve this， two deep geothermal granite models were established， and ?numerical simulations of acoustic and elastic waves were implemented using the equivalent staggered grid finite difference method. The numerical simulation results reveal that the velocity of geothermal granite is much higher than that of surrounding rock due to temperature influences， which results in a phenomenon known as high-speed shielding and weakens the transmitted wave energy， consequently limiting the propagation of seismic waves within the lower part of the geothermal rock mass. In addition， compared with acoustic waves， elastic waves offer richer information about the wavefield. The conversion of wave mode and energy makes seismic records of elastic waves more complex than acoustic waves.

Keywords ： deep geothermal rock; seismic wave propagation mechanism; numerical simulation; finite difference modeling

中深層地?zé)豳Y源包括常規(guī)水熱型地?zé)嵯到y(tǒng)（深度在0.2～3 km）和增強型地?zé)嵯到y(tǒng)，其中增強型地?zé)嵯到y(tǒng)又叫干熱巖。干熱巖溫度大于200 ℃，埋深超過3 km，內(nèi)部不存在流體或僅存在少量流體，其蘊藏的能量是煤炭、石油和天然氣總數(shù)的30倍 ?［1- 2］ 。干熱巖是一種可再生的清潔能源，在供暖、工業(yè)發(fā)電和農(nóng)業(yè)種植等領(lǐng)域都有廣泛的應(yīng)用 ?［3］ 。因此開發(fā)、利用深層地?zé)豳Y源對于國家減少二氧化碳排放量，堅定實施可持續(xù)發(fā)展戰(zhàn)略具有重要意義 ?［4-5］ 。相比于發(fā)達國家，中國對于深層地?zé)豳Y源勘查開發(fā)及其技術(shù)還處于初級階段，沒有形成系統(tǒng)的勘探開發(fā)方法。現(xiàn)階段主要面臨3大挑戰(zhàn) ?［6］ ：第一，中深層（3～6 km）資源儲量豐富，成儲與成藏機制認(rèn)識淺；第二，復(fù)雜地質(zhì)體反射信號雜亂，高速屏蔽作用強，信噪比低，地震資料品質(zhì)差，高精度成像困難；第三，深層地?zé)醿拥刭|(zhì)-地球物理特征復(fù)雜，地球物理響應(yīng)特征不明確，中深層地?zé)崽綔y成功率低。這些問題都是對于深層地?zé)釒r體下地震波的傳播規(guī)律和波場特征認(rèn)識較少。因此，研究地震波在深層地?zé)釒r體中的傳播規(guī)律及波場特征對于深層地?zé)豳Y源的勘探開發(fā)至關(guān)重要。有限差分波動方程數(shù)值模擬是研究地震波傳播規(guī)律的主要方法 ?［7］ 。董良國等 ?［8］ 應(yīng)用交錯網(wǎng)格有限差分法實現(xiàn)了一階應(yīng)力-速度彈性波方程數(shù)值模擬。雍鵬等 ?［9］ 使用共軛梯度法優(yōu)化時空差分算子后實現(xiàn)了一階應(yīng)力-速度聲波方程數(shù)值模擬。黃建平等 ?［10］ 、李慶洋 ?［11］ 、李娜等 ?［12］ 、慕鑫茹等 ?［13］ 對復(fù)雜介質(zhì)中地震波的正演模擬理論、方法及其策略進行了系統(tǒng)總結(jié)。這些研究大都是數(shù)值模擬方法及算法優(yōu)化方面，對于深層地?zé)釒r體環(huán)境下地震波傳播規(guī)律的研究較少。筆者建立兩種深層地?zé)峄◢弾r的速度模型，并對該模型的聲波和彈性波方程進行數(shù)值模擬；在數(shù)值模擬過程中使用等效交錯網(wǎng)格進行剖分，在保證與交錯網(wǎng)格同等精度的同時，減少內(nèi)存的占用量。

1 方法原理

1.1 聲波與彈性波介質(zhì)數(shù)值模擬方法

在聲波介質(zhì)中，二維一階應(yīng)力-速度聲波方程 ?［9］ 為

u ?t =ρv 2 ?P ???v x ?x + ?v z ?z ?，

v x ?t = 1 ρ ??u ?x ?，

v z ?t = 1 ρ ??u ?z ?. ?（1）

式中，v x和v z為質(zhì)點的振動速度；u為應(yīng)力；ρ為密度；v ?P 為縱波速度。

在彈性波介質(zhì)中二維一階應(yīng)力-速度彈性波方程 ?［8］ 為

v x ?t = 1 ρ ???σ ?xx ??x + ?σ ?xz ??z ?，

v z ?t = 1 ρ ???σ ?xz ??x + ?σ ?zz ??z ?，

σ ?xx ??t =（λ+2μ） ?v x ?x +λ ?v z ?z ?，

σ ?zz ??t =（λ+2μ） ?v z ?z +λ ?v x ?x ?，

σ ?xz ??t =μ ??v z ?x + ?v x ?z ?. ??（2）

式中， σ ?xx 和σ ?zz 分別為x、z方向上的正應(yīng)力；σ ?xz 為剪切應(yīng)力； λ、μ為拉梅系數(shù)。

對式（1）中第一式進行時間求導(dǎo)，并把第二、三式代入其中，可以得到二階的常密度聲波方程為

2u ?t 2 =v 2 ?P ????2u ?x 2 + ??2u ?z 2 ?. （3）

同理式（2）可以變?yōu)槌Ｃ芏葟椥圆ǚ匠虨?/p>

2σ ?xx ??t 2 =（λ+2μ） ??2σ ?xx ??x 2 +（λ+μ） ??2σ ?zz ??x z +μ ??2σ ?xx ??z 2 ?，

2σ ?zz ??t 2 =（λ+2μ） ??2σ ?zz ??z 2 +（λ+μ） ??2σ ?xx ??x z +μ ??2σ ?zz ??x 2 ?. ?（4）

由于波動方程的形式不同，數(shù)值模擬過程中采用的網(wǎng)格剖分格式也不同，對于式（1）和（2）中的一階方程通常采用交錯網(wǎng)格有限差分法，其差分格式（以聲波方程為例）為

u ?n+1 ??i，j =u n ?i，j +ρv 2 ?P ??Δ t ?Δ x ∑ L ?m=1 ?a m v ?n+ 1 2 ???x ?i+ 2m-1 2 ，j ?-v ?n+ 1 2 ???x ?i- 2m-1 2 ，j ??+

ρv 2 ?P ??Δ t ?Δ z ∑ L ?m=1 ?a m v ?n+ 1 2 ???z ?i，j+ 2m-1 2 ??-v ?n+ 1 2 ???z ?i，j- 2m-1 2 ???，

v ?n+ 1 2 ???x ?i+ 1 2 ，j ?=

v ?n- 1 2 ???x ?i+ 1 2 ，j ?+ ?Δ t ρ Δ x ∑ L ?m=1 ?a m（u n ?i+m，j -u n ?i-m+1，j ），

v ?n+ 1 2 ???z ?i，j+ 1 2 ??=v ?n- 1 2 ???z ?i，j+ 1 2 ??+ ?Δ t ρ Δ z ∑ L ?m=1 ?a m（u n ?i，j+m -u n ?i，j-m+1 ）. ?（5）

對于二階方程，通常采用規(guī)則網(wǎng)格有限差分法，其差分格式 ??［14］ ?為

u ?n+1 ??i，j =2u n ?i，j -u ?n-1 ??i，j +v 2 ?P ??Δ t 2 ?Δ x 2 ∑ L ?m=1 ?a m（u n ?i+m，j -u n ?i-m，j ）+

v 2 ?P ??Δ t 2 ?Δ z 2 ∑ L ?m=1 ?a m（u n ?i，j+m -u n ?i，j-m ）. （6）

式中，i，j分別為空間x和z方向上的網(wǎng)格點位置； Δ x、 Δ z和 Δ t分別為空間和時間差分間隔；n為時間上的網(wǎng)格點位置。

從差分格式中可以看出，在數(shù)值模擬過程中，交錯網(wǎng)格需要存儲2個應(yīng)力場和4個速度場，而規(guī)則網(wǎng)格只需存儲3個應(yīng)力場。相較于交錯網(wǎng)格，規(guī)則網(wǎng)格可以減少內(nèi)存的使用；但由于二階方程是常密度假設(shè)，忽略了密度的影響，當(dāng)密度變化較大時，規(guī)則網(wǎng)格的精度會降低，即交錯網(wǎng)格精度高于規(guī)則網(wǎng)格。

為了在保證精度的前提下減小存儲量，使用等效交錯網(wǎng)格 ?［7，15-16］ 進行剖分，即保證了與交錯網(wǎng)格同等的精度與穩(wěn)定性，也減少了內(nèi)存使用量，如表1所示。在大型的實際資料處理以及三維介質(zhì)的成像或者反演中具有較大的應(yīng)用前景。以聲波方程為例，給出了等效交錯網(wǎng)格的差分格式為

u ?n+1 ??i，j =2u n ?i，j -u ?n-1 ??i，j +v 2 ?P ??Δ t 2 ?Δ x 2 ∑ L ?m=1 ?∑ N ?l=1 ?a ma l［（u n ?i+m+l-1，j -u n ?i-m+l，j ）-（u n ?i+m-l，j -u n ?i-m-l+1，j ）］+ v 2 ?P ??Δ t 2 ?Δ z 2 ∑ L ?m=1 ?∑ N ?l=1 ?b mb l［（u n ?i，j+m+l-1 -u n ?i，j-m+l ）-（u n ?i，j+m-l - u n ?i，j-m-l+1 ）］. （7）

1.2 震源加載方式

加載震源是在有限差分網(wǎng)格點上施加特定的震源，地震勘探中通常使用炸藥作為震源產(chǎn)生地震波。爆炸震源是在有限差分網(wǎng)格點上施加作用力，作用是在球腔內(nèi)產(chǎn)生對稱的徑向作用力。均勻各向同性中，爆炸震源相當(dāng)于純 P波震源，只產(chǎn)生純P波波場，其彈性波波場如圖1所示。 只存在P波，如果模型發(fā)生變化，將在速度分界面上產(chǎn)生反射S波及透射S波。

地震波數(shù)值模擬過程中采用合適的震源函數(shù)至關(guān)重要，雷克子波是常用的震源函數(shù)之一，其函數(shù)解析式為

f（t）=［1-2 π ?2f 2 ??avg ??（t-t 0） ?2］ exp［-π ?2f 2 ??avg ??（t-t 0） ?2］. （8）

式中，f ??avg ?為子波主頻；t 0為時間延遲項。

本文中模型測試中雷克子波選取的主頻均為15 ?Hz。

1.3 邊界條件

描述地震波在地層中的傳播規(guī)律和特征的波動方程是基于無限空間介質(zhì)建立的，但由于使用計算機進行數(shù)值模擬運算時，計算范圍有限，必然會人為地限定計算區(qū)域。在數(shù)值模擬過程中，當(dāng)?shù)卣鸩ㄈ肷涞竭@種人為的截斷邊界時會產(chǎn)生反射波。這種邊界反射在實際中是不存在的，會嚴(yán)重干擾有效信號，因此為了消除這種人工邊界反射必須施加合適的邊界條件，吸收或者衰減這種邊界反射波的能量。數(shù)值模擬中常用的邊界條件有：①衰減邊界條件；②完全匹配層邊界條件 （perfectly matched layer，PML） ?［17］ 。

模型測試中選用的邊界條件均為衰減邊界條件 ?［18］ ，通過在人工邊界處設(shè)置衰減帶消除邊界反射。在每一時間步長計算中衰減帶的每個網(wǎng)格點上的波場振幅都需要乘以一個衰減因子表示為

G（i）= exp ［-c 2（N-i） 2］. （9）

式中，i為計算點到邊界的距離；c為衰減系數(shù)，模型試算中取0.12；N為衰減帶的網(wǎng)格數(shù)，模型試算中取25。

這種邊界條件的設(shè)置方法比較簡單，適用性強，是地震勘探數(shù)值模擬中常用的消除人工邊界反射的方法之一。

2 數(shù)值試驗

設(shè)計兩個地?zé)峄◢弾r的速度模型：一個為高放射性花崗巖侵入體模型，另一個是近代火山模型。兩個模型都具有因地?zé)岘h(huán)境而產(chǎn)生高溫、高速的火成巖。

2.1 高放射性侵入花崗巖體

圖2為高放射性侵入花崗巖體模型。 模型大小為5 km×5 km，其特點是在模型左下部分有一塊高放射性花崗巖侵入體，其速度在6.1～6.6 km/s內(nèi)隨機變化。分別采用聲波方程和彈性波方程等效交錯網(wǎng)格有限差分法對該模型進行數(shù)值模擬。

2.1.1 聲波數(shù)值模擬

聲波方程數(shù)值模擬中，時間采樣間隔為0.000 5 s， ?空間差分間隔在 x和z 方向上均取10 m；主頻為15 Hz的ricker子波作為震源設(shè)置在模型中（2 500 m，10 m）位置處。圖3為不同時刻的聲波波場快照，圖4為聲波模擬地震記錄。

由圖3可以看出，當(dāng)聲波傳播時間為1.25 s時，聲波傳播到兩個速度交界的凸界面，開始出現(xiàn)反射波；當(dāng)傳播時間為1.5 s時，遇到第二個速度交界面，界面下出現(xiàn)高速的花崗巖侵入體，會發(fā)生高速屏蔽現(xiàn)象，透射波能量變?nèi)?，限制了地震波的穿透能力，在波場快照中表現(xiàn)為波場比較模糊（圖3（c））。

2.1.2 彈性波數(shù)值模擬

彈性波方程數(shù)值模擬中，時間采樣間隔為0000 5 s， 空間差分間隔在 x和z 方向上均取10 m；縱波速度與橫波速度之比為1.732；主頻為15 Hz的ricker子波作為震源設(shè)置在模型中（2 500 m，10 m）位置處。圖5為不同時刻彈性波波場快照（左、右圖分別為彈性波的水平和垂直分量），圖6為彈性波模擬地震記錄。

與聲波不同的是，彈性波不僅有P波還存在S波，相同傳播時間1.25 s時，圖5中出現(xiàn)了反射S波。另外，P（S）波不僅可以轉(zhuǎn)換為P（S）波，也可以轉(zhuǎn)換為S（P）波，因此波場信息比聲波更豐富、復(fù)雜，更能表征深層地?zé)釒r體。

2.2 近代火山型

圖7為近代火山型地?zé)釒r體模型，模型大小為7 km×5 km，由于近代火山噴發(fā)結(jié)束后保持休眠的時間過短，其形成的火山口離地表較近，巖漿凝固后 形成的火成巖穿透不同的地層到達近地表，并且具有較高的速度。分別采用聲波方程和彈性波方程等效交錯網(wǎng)格有限差分法對該模型進行了數(shù)值模擬。

2.2.1 聲波數(shù)值模擬

聲波方程數(shù)值模擬中，時間采樣間隔為0.000 5 s，空間差分間隔在 x和z 方向上均取10 m；主頻為15 Hz的ricker子波作為震源設(shè)置在模型中（3 500 m，10 m）位置處。圖8為不同時刻聲波波場快照，圖9為聲波模擬地震記錄。

從圖8、 9中可以看出，近代火山模型的高速屏蔽現(xiàn)象比高放射性侵入花崗巖模型較嚴(yán)重，在水平方向4 km附近的火山通道往下的波場能量都比較弱，使得地震記錄的信噪比較低。

2.2.2 彈性波數(shù)值模擬

彈性波方程數(shù)值模擬中，時間采樣間隔為0000 5 s，空間差分間隔在 x和z 方向上均取10 m；縱波速度與橫波速度比為1.732；主頻為15 Hz的ricker子波作為震源設(shè)置在模型中（3 500 m，10 m）位置處。圖10為不同時刻的彈性波波場快照（左側(cè)為水平分量，右側(cè)為垂直分量）。圖11為彈性波模擬地震記錄。

除了高速屏蔽現(xiàn)象外，模型的構(gòu)造起伏也使得反射同相軸發(fā)生畸變，隆起和凹陷構(gòu)造會使地震波能量發(fā)散和聚集，對反射波、直達波的同相軸形狀影響很大，在成像過程中容易形成構(gòu)造假象。

3 結(jié) 論

（1）等效交錯網(wǎng)格有限差分法實現(xiàn)聲波與彈性波方程正演模擬，可以在保證精度的同時減少內(nèi)存的使用量，在處理大型的實際資料及三維成像和反演中具有較好的應(yīng)用前景。

（2）全彈性波場比聲波波場具有更豐富的波場信息，更能全面地表征深層地?zé)釒r體；深層地?zé)岘h(huán)境下形成的高溫侵入巖體與圍巖相比具有較高的速度，會形成高速屏蔽現(xiàn)象，使得透射波能量變?nèi)?，限制了地?zé)釒r體下部地震波的傳播能量，在波場快照和炮記錄中表現(xiàn)為在高速體附近波形比較模糊。

參考文獻 ：

［1］ ?JAYA M S， SHAPIRO S A， KRISTINSDOTTIR L H， et al. Temperature dependence of seismic properties in geothermal rocks at reservoir conditions［J］. Geothermics， 2010，39（1）：115-123.

［2］ ?王丹，魏水建，賈躍瑋，等.地?zé)豳Y源地震勘探方法綜述［J］.物探與化探，2015，39（2）：253-261.

WANG Dan， WEI Shuijian， JIA Yuewei， et al. An overview of methods for geothermal seismic exploration［J］. Geophysical and Geochemical Exploration， 2015，39（2）：253-261.

［3］ ?SCHMELZBACH C， GREENHALGH S， REISER F， et al. Advanced seismic processing/imaging techniques and their potential for geothermal exploration［J］. Interpretation， 2016，4（4）：SR1-SR18.

［4］ ?藺文靜，劉志明，王婉麗，等.中國地?zé)豳Y源及其潛力評估［J］.中國地質(zhì)，2013，40（1）：312-321.

LING Wenjing， LIU Zhiming， WANG Wanli， et al. The assessment of geothermal resources potential of China［J］. Geology in China， 2013，40（1）：312-321.

［5］ ?BIROL F， OLERJARNIK P. Will. China lead the world into a clean-energy future［J］. Economics of Energy and Environmental Policy， 2012，1（1）：5-9.

［6］ ?楊冶，姜志海，岳建華，等.干熱巖勘探過程中地球物理方法技術(shù)應(yīng)用探討［J］. 地球物理學(xué)進展， 2019，34（4）：1556-1567.

YANG Ye， JIANG Zhihai， YUE Jianhua， et al. Discussion on application of geophysical methods in hot dry rock（HDR）exploration［J］. Progress in Geophysics， 2019，34（4）：1556-1567.

［7］ ?YONG Peng， HUANG Jianping， LI Zhenchun， et al. Optimized equivalent staggered-grid FD method for elastic wave modeling based on plane wave solutions［J］. Geophysical Supplements to the Monthly Notices of the Royal Astronomical Society， 2016，208（2）：1157-1172.

［8］ ?董良國，馬在田，曹景忠，等.一階彈性波方程交錯網(wǎng)格高階差分解法［J］.地球物理學(xué)報，2000，43（3）：411-419.

DONG Liangguo， MA Zaitian， CAO Jingzhong， et al. The staggered-grid high-order difference method of first-order elastic wave equation［J］. Chinese Journal of Geophysics， 2000，43（3）：411-419.

［9］ ?雍鵬，黃建平，李振春，等.一種時空域優(yōu)化的高精度交錯網(wǎng)格差分算子與正演模擬［J］.地球物理學(xué)報，2016，59（11）：4223-4233.

YONG Peng， HUANG Jianping， LI Zhenchun， et al. Optimized staggered-grid finite-difference method in time-space domain based on exact time evolution schemes［J］. Chinese Journal of Geophysics， 2016，59（11）：4223-4233.

［10］ ?黃建平，李慶洋，雍鵬，等.復(fù)雜介質(zhì)地震波正演模擬方法及優(yōu)化［M］.北京：科學(xué)出版社，2022：1-203.

［11］ ?李慶洋. 復(fù)雜介質(zhì)地震波正演模擬與最小二乘偏移研究［D］.青島：中國石油大學(xué)（華東），2017.

LI Qingyang. Study on seismic forward modeling and least-squares migration in complex media［D］.Qingdao：China University of Petroleum （East China）， 2017.

［12］ ?李娜，李振春，黃建平，等.Lebedev網(wǎng)格與標(biāo)準(zhǔn)交錯網(wǎng)格耦合機制下的復(fù)雜各向異性正演模擬［J］.石油地球物理勘探，2014，49（1）：121-131.

LI Na， LI Zhenchun， HUANG Jianping， et al. Complex anisotropic forward modeling under the coupling mechanism of Lebedev grid and standard staggered grid［J］.Oil Geophysical Prospecting， 2014，49（1）：121-131.

［13］ ?慕鑫茹，黃建平，黎國龍，等.黏聲TTI介質(zhì)中純qP波逆時偏移成像方法［J］.中國石油大學(xué)學(xué)報（自然科學(xué)版），2023，47（2）：44-52.

MU Xinru， HUANG Jianping， LI Guolong， et al. Reverse time migration in viscoacoustic TTI media based on pure qP wave equation［J］. Journal of China University of Petroleum （Edition of Natural Science）， 2023，47（2）：44-52.

［14］ ?黃建平，婁璐烽，彭煒颋，等.一種基于拉格朗日乘子的空間域差分系數(shù)優(yōu)化方法［J］.中國石油大學(xué)學(xué)報（自然科學(xué)版），2022，46（4）：30-40.

HUANG Jianping， LOU Lufeng， PENG Weiting， et al. An optimization method of finite difference coefficient in spatial domain based on Lagrange multipliers［J］. Journal of China University of Petroleum （Edition of Natural Science）， 2022，46（4）：30-40.

［15］ ?段沛然，李青陽，趙志強，等.等效交錯網(wǎng)格高階有限差分法標(biāo)量波波場模擬［J］.地球物理學(xué)進展，2019，34（3）：1032-1040.

DUAN Peiran， LI Qingyang， ZHAO Zhiqiang， et al. High-order finite-difference method based on equivalent staggered grid scheme for scalar wavefield simulation［J］. Progress in Geophysics， 2019，34（3）：1032-1040.

［16］ ?ZOU Qiang， HUANG Jianping， YONG Peng， et al. 3D elastic waveform modeling with an optimized equivalent staggered-grid finite-difference method［J］. Petroleum Science， 2020，17：967-989.

［17］ ?黃建平，楊宇，李振春，等.基于M－PML邊界的Lebedev網(wǎng)格起伏地表正演模擬方法及穩(wěn)定性分析［J］.中國石油大學(xué)學(xué)報（自然科學(xué)版），2016，40（4）：47-56.

HUANG Jianping， YANG Yu， LI Zhenchun， et al. Lebedev grid finite-difference modeling for irregular free-surface and stability analysis based on M-PML boundary condition［J］. Journal of China University of Petroleum （Edition of Natural Science）， 2016，40（4）：47-56.

［18］ ?CERJAN C， KOSLOFF D， KOSLOFF R， et al. A non-reflecting boundary condition for discrete acoustic and elastic wave equations［J］. Geophysics， 1985，50（4）：705-708.

（編輯 沈玉英）