深層地熱花崗巖體地震波數值模擬及傳播機制

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

摘要 ：深層地熱資源是一種可再生的、蘊藏巨大能量的清潔能源，但其地球物理響應特征不明確，中深層地熱探測成功率較低。為了研究地震波在深層地熱巖體中的傳播規律與波場特征，建立兩種深層地熱花崗巖模型，并使用等效交錯網格有限差分法實現聲波與彈性波的數值模擬。結果表明：地熱花崗巖體速度在溫度的影響下要高于圍巖的，會產生高速屏蔽現象，使得透射波能量變弱，限制了地熱巖體下部地震波傳播能量；相比于聲波，彈性波具有更豐富的波場信息，波型與能量轉換使得彈性波地震記錄也比聲波地震記錄復雜。

關鍵詞 ：深層地熱巖； 地震波傳播機制； 數值模擬； 有限差分模擬

中圖分類號 ：P 631.4 ???文獻標志碼 ：A

引用格式 ：黃建平，劉英輝，李偉，等.深層地熱花崗巖體地震波數值模擬及傳播機制［J］.中國石油大學學報（自然科學版），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

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

1 方法原理

1.1 聲波與彈性波介質數值模擬方法

在聲波介質中，二維一階應力-速度聲波方程 ?［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為質點的振動速度；u為應力；ρ為密度；v ?P 為縱波速度。

在彈性波介質中二維一階應力-速度彈性波方程 ?［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方向上的正應力；σ ?xz 為剪切應力； λ、μ為拉梅系數。

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

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

同理式（2）可以變為常密度彈性波方程為

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）

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

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）

對于二階方程，通常采用規則網格有限差分法，其差分格式 ??［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方向上的網格點位置； Δ x、 Δ z和 Δ t分別為空間和時間差分間隔；n為時間上的網格點位置。

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

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

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 震源加載方式

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

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

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 邊界條件

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

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

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

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

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

2 數值試驗

設計兩個地熱花崗巖的速度模型：一個為高放射性花崗巖侵入體模型，另一個是近代火山模型。兩個模型都具有因地熱環境而產生高溫、高速的火成巖。

2.1 高放射性侵入花崗巖體

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

2.1.1 聲波數值模擬

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

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

2.1.2 彈性波數值模擬

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

與聲波不同的是，彈性波不僅有P波還存在S波，相同傳播時間1.25 s時，圖5中出現了反射S波。另外，P（S）波不僅可以轉換為P（S）波，也可以轉換為S（P）波，因此波場信息比聲波更豐富、復雜，更能表征深層地熱巖體。

2.2 近代火山型

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

2.2.1 聲波數值模擬

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

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

2.2.2 彈性波數值模擬

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

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

3 結 論

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

（2）全彈性波場比聲波波場具有更豐富的波場信息，更能全面地表征深層地熱巖體；深層地熱環境下形成的高溫侵入巖體與圍巖相比具有較高的速度，會形成高速屏蔽現象，使得透射波能量變弱，限制了地熱巖體下部地震波的傳播能量，在波場快照和炮記錄中表現為在高速體附近波形比較模糊。

參考文獻 ：

［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］ ?王丹，魏水建，賈躍瑋，等.地熱資源地震勘探方法綜述［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］ ?藺文靜，劉志明，王婉麗，等.中國地熱資源及其潛力評估［J］.中國地質，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］ ?楊冶，姜志海，岳建華，等.干熱巖勘探過程中地球物理方法技術應用探討［J］. 地球物理學進展， 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］ ?董良國，馬在田，曹景忠，等.一階彈性波方程交錯網格高階差分解法［J］.地球物理學報，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］ ?雍鵬，黃建平，李振春，等.一種時空域優化的高精度交錯網格差分算子與正演模擬［J］.地球物理學報，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］ ?黃建平，李慶洋，雍鵬，等.復雜介質地震波正演模擬方法及優化［M］.北京：科學出版社，2022：1-203.

［11］ ?李慶洋. 復雜介質地震波正演模擬與最小二乘偏移研究［D］.青島：中國石油大學（華東），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網格與標準交錯網格耦合機制下的復雜各向異性正演模擬［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介質中純qP波逆時偏移成像方法［J］.中國石油大學學報（自然科學版），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］ ?黃建平，婁璐烽，彭煒颋，等.一種基于拉格朗日乘子的空間域差分系數優化方法［J］.中國石油大學學報（自然科學版），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］ ?段沛然，李青陽，趙志強，等.等效交錯網格高階有限差分法標量波波場模擬［J］.地球物理學進展，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網格起伏地表正演模擬方法及穩定性分析［J］.中國石油大學學報（自然科學版），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.

（編輯 沈玉英）