面向深部地熱巖體的彈性波逆時偏移成像方法

黃建平 楊秀金 張鑫 王揚州 陳亮 高建明

摘要 ：地熱能作為一種清潔能源正受到全世界的日益關注，因此對于深部地熱巖體的高精度成像格外重要。傳統的彈性波逆時偏移方法是利用耦合縱橫波直接成像，從而得到耦合波場的成像結果;然而這種方法可能會出現串擾假象，導致成像結果不夠清晰；為了解決這個問題，采用基于解耦的彈性波方程，實現縱橫波場的分離；通過利用內積成像條件，對兩個典型的干熱巖模型進行數值測試。結果表明，相較于耦合波場成像結果，基于解耦方程分離得到的彈性波成像剖面具有更清晰的同相軸，深部能量更加均衡，該方法能夠實現對深部地熱巖體的高質量成像。

關鍵詞 ：深部地熱巖體； 彈性波成像； 縱橫波分離； 逆時偏移

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

引用格式 ：黃建平，楊秀金，張鑫，等.面向深部地熱巖體的彈性波逆時偏移成像方法［J］.中國石油大學學報（自然科學版），2024，48（1）：55-62.

HUANG Jianping， YANG Xiujin， ZHANG Xin， et al. Elastic wave reverse time migration imaging method for deep geothermal rock mass［J］. Journal of China University of Petroleum （Edition of Natural Science），2024，48（1）：55-62.

Elastic wave reverse time migration imaging method

for deep geothermal rock mass

HUANG Jianping 1， YANG Xiujin 1， ZHANG Xin 2， WANG Yangzhou 2， CHEN Liang 1， GAO Jianming 2

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

2.Shandong Energy Group South America Company Limited， Qingdao 266000， China）

Abstract ： Geothermal energy is gaining global attention as a clean energy source， underscoring the significance of precise imaging of deep geothermal reservoirs. The conventional elastic reverse time migration （E-RTM） technique relies on direct imaging of coupled P- and S- waves， which often leads to ?imaging artifacts and unclear images. In this study， we propose an approach based on the decoupled elastic wave equation to independently handle ?P- and S- wavefields.We apply this approach to two typical hot dry rock models using inner product imaging conditions. The numerical simulations demonstrate that the decoupled elastic wave imaging profiles ?provide clearer event representations and a more evenly distributed energy profile in the deeper regions compared to imaging results obtained from coupled wavefields. This method showcases the potential for achieving high-quality imaging of deep geothermal reservoirs.

Keywords ： deep geothermal rock mass; elastic wave imaging; separation of P- and S-waves; reverse time migration

多波多分量勘探技術是地層巖性識別、陡傾界面成像的有效手段，能夠提供更詳細的地下介質信息，在油氣和地熱等能源探測領域具有重要意義。相比于傳統的聲波逆時偏移，彈性波逆時偏移考慮了地震波在地下傳播時的多波耦合特征，可以較為充分地利用彈性波多分量地震數據中豐富的縱橫波場信息，實現對復雜地質結構區塊的高精度成像 ?［1］ 。傳統彈性波逆時偏移由于成像時未進行縱橫波場分離，導致耦合波場間有串擾的產生，嚴重影響成像質量 ?［2］ 。為了去除串擾偽影和提高成像質量，國內外學者做了許多研究 ?［3-5］ ，提出采用波場分離的方法可以進行有效壓制。常用的波場分離方法有基于 Helmholtz 分解的波場分離以及基于解耦方程的波場分離。Helmholtz分解自提出以來，就被廣泛地應用于彈性波場分解中 ?［6-8］ 。Sun等 ??［9］ 在進行彈性波場外推時，使用散度及旋度算子分離出地表附近的縱波場和橫波場，再分別進行PP和轉換PS的兩次聲波逆時偏移。Dellinger等 ??［10］ 基于Helmholtz理論提出了波數域的類散度算子和類旋度算子，并利用Christoffel方程實現了各向異性均勻彈性介質中的波場分解。Yan等 ??［11］ 將分離算子從波數域進一步轉換到空間域，并制定了一個濾波器分離具有垂直對稱軸的橫向各向同性介質中的波場。但是以上基于Helmholtz 分解得到的縱波和橫波波場通常被理解為標量勢與矢量勢，得到的振幅以速度的倒數衰減，相位推遲π/2 ?［1］ ，不能準確地完成波場分離，需要進行額外的保幅處理和相位校正。為了探尋一種波場分離方法以更貼近于彈性波的傳播規律，基于解耦方程進行波場分離的研究也取得了許多進展。馬德堂等 ??［12］ 根據 Helmholtz 分解理論推導出縱波和橫波分離的二階位移方程，通過偽譜法實現了在波場外推過程中完成縱波和橫波的分離，這為后續研究提供了一種新思路。李振春等 ?［13］ 和Zhang等 ?［14］ 進一步推導出一階速度-應力彈性波場的分離方程，并利用交錯網格高階有限差分算法驗證了二維方程的可行性；Xiao 等 ?［15］ 引入輔助 P 波和 S 波應力場進行波場分解，實現了VSP局部彈性逆時偏移，使得基于解耦方程的波場分離更加具有矢量傳播的物理意義； Wang等 ??［16］ 證明 Zhang 和 Xiao的公式是形式不一樣的等價方程，并給出了新縱波和橫波分離的彈性波方程?；诮怦罘匠痰玫降氖鞘噶靠v波和矢量橫波，若仍然采取互相關成像條件進行成像，則得到的反射率模型與傳統意義上的不一樣 ?［17］ ，會給解釋工作人員帶來困難。為此適用于矢量波場成像的內積成像條件被進一步發展 ?［18-19］ ，這種成像條件能較為準確地提取矢量波場各位置的成像值。然而逆時偏移是一種構造成像方法，對地熱體的準確模擬和成像還可以考慮溫度和物質形態等因素變化 ?［20-22］ ，進而發展出更適用的方法。筆者首先介紹波場分離的基本公式和所應用的成像條件；其次將上述提到的公式和成像條件對層狀模型進行波場模擬，通過波場快照和單道分析說明其可行性；最后將矢量波場分離公式應用到以侵入花崗型和近代火山型為代表的深部干熱巖體模型進行逆時偏移，以檢測本方法的優越性、適用性和準確性。

1 方法原理

本文中彈性波逆時偏移可分為3步：①構建震源波場并保存各時刻波場值；②構建檢波波場并保存各時刻波場值；③使用震源歸一化內積成像條件進行成像。

1.1 二維波動方程

波動方程反映的是受到激勵時地下介質每時刻的運動情況，它能表示地震波的傳播過程。二維情況下，通常假定地震波沒有 y分量，即 ???y =0，所以可給出各向同性介質縱橫波耦合傳播的一階速度-應力方程 ??［23］ 為

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

v z ?t = 1 ρ ???τ ?zx ??x + ?τ ?zz ??z ?，

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

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

τ ?xz ??t =μ ??v x ?z + ?v z ?x ?. ??（1）

式中，v x和v z分別為地震波的橫向速度和縱向速度；τ ?xx 、τ ?zz 和τ ?xz 分別為x方向正應力、z方向正應力和xz方向的剪切應力；λ和μ為拉梅系數；ρ為介質密度。

1.2 彈性波場分離

傳統的波場分離方法是基于 ?Helmholtz 分解對總波場取散度得到縱波場，對波場取旋度得到橫波場。此類波場分離方法對原始矢量波場施加了散度和旋度算子，相當于對其求空間導數，導致分離得到的縱波和橫波分量產生90°相移，振幅分別產生 ??1 v ?P ?、 1 v ?S ?的改變。Wang等 ??［16］ 通過引入一個輔助縱波應力，進一步發展了一種波場矢量分解公式，能較好地解決這種現象，有效地保存了彈性波場的振幅和相位信息。

τ ?P ??t =（λ+2μ） ??v x ?x + ?v z ?z ?，

v ?x P ???t = 1 ρ ??τ ?P ??x ?，

v ?z P ???t = 1 ρ ??τ ?P ??z ?，

v ?x S ?=v x-v ?x P ?，

v ?z S ?=v z-v ?z P ?. ?（2）

式中，τ ?P 為引入的縱波應力；v ?x P ?和v ?z P ?分別為縱波在x和z方向上的速度分量；v ?x S ?和v ?z S ?分別為橫波在x和z方向上的速度分量。

為了實現震源波場延拓與檢波波場延拓，采用時間二階、空間八階的規則網格有限差分格式對式（1）和（2）進行離散。在進行波場模擬時設定完全匹配層（PML）吸收邊界，減弱由于人工截斷邊界而導致的邊界反射現象，更有效地壓制邊界反射和提高目標模擬區域的成像精度。

1.3 內積成像條件

考慮到基于解耦波動方程得出的縱波和橫波是矢量形式，若仍采用傳統的互相關成像條件來成像，不能反映出反射系數的物理意義，內積成像條件是解決此問題的有效手段 ?［19］ 。復雜模型逆時偏移剖面的深部成像通常存在照明不足等問題，使用震源歸一化技術可以有效地補償地震波的能量損失，有利于增強模型深部的能量。為此給出一種震源歸一化內積成像條件為

I ??PP ?（ x ）= ∫ S ??P （ x ，t）· R ??P （ x ，t） d t ∫ S ??P （ x ，t）· S ??P （ x ，t） d t ?，

I ??PS ?（ x ）= ∫ S ??P （ x ，t）· R ??S （ x ，t） d t ∫ S ??P （ x ，t）· S ??P （ x ，t） d t ?，

I ??SP ?（ x ）= ∫ S ??S （ x ，t）· R ??P （ x ，t） d t ∫ S ??S （ x ，t）· S ??S （ x ，t） d t ?，

I ??SS ?（ x ）= ∫ S ??S （ x ，t）· R ??S （ x ，t） d t ∫ S ??S （ x ，t）· S ??S （ x ，t） d t ?. ?（3）

式中，I（ x ）為各空間位置成像值；S為震源波場；R為檢波波場；下標 P和S 分別表示縱波和橫波；符號“·”為內積算子。

例如，I ??PS ?（ x ）表示震源縱波與檢波橫波在成像點 x 處的內積成像值，內積計算格式由 S ??P （ x ，t）· R ??S （ x ，t）= S ?P （ x ，t） ?R ?S （ x ，t） ?cos 〈 S ??P （ x ，t）， ?R ??S （ x ，t）〉 得出，其中S ?P （ x ，t）表示震源縱波場在t時刻、 x 位置處的波場值，R ?S （ x ，t）表示檢波橫波場在t時刻、 x 位置處的波場值。

通常情況下地震波的主要能量入射角小于40°，然而內積成像條件發生極性反轉的條件是入射角大于45°，因此采用內積成像條件不會出現極性反轉現象 ?［17］ 。針對使用雙程波動方程進行波場模擬而導致成像結果存在低頻噪音的問題，采用Laplacian濾波對所有逆時偏移剖面進行處理，有效地壓制了此現象。

2 數值模擬

為了驗證本文提到的波場分離方法的有效性以及對深部高溫地熱巖體構造成像的適用性，通過對層狀模型、侵入花崗型和近代火山型干熱巖模型進行測試。層狀模型用于驗證基于解耦方程進行波場分離的準確性，侵入花崗型和近代火山型用于檢驗基于解耦方程對深部地熱巖體等復雜地質結構成像的優越性和適用性。本文中均使用時間二階、空間八階的規則網格有限差分格式完成波場延拓過程，采用集中力震源完成激發（震源增益為1 000）和完全匹配層邊界條件進行吸收。橫波速度根據泊松體的縱橫波比求得，即 V ?S = V ?P 1.732 ，介質密度根據Gardner 公式求取，即ρ=0.31V ?0.25 ??P。

2.1 層狀模型

為了更為直觀地認識基于解耦方程的波場分離方法，采用一個三層模型（圖1）測試其有效性和準確性，該模型假定縱波速度分別為3 000、4 000、3 000 m/s，各層厚度分別為600、1 800和600 m。模型的網格點數為401×301，網格間距為10 m，時間采樣率為1 ms。震源位于模型正中央（201，151），震源函數為25 Hz的雷克子波。

從0.45 s時刻記錄到的波場快照中可以看到，基于解耦方程的波場分離效果較好，可以將純縱波、純橫波和轉換波從耦合波場中實現分離，各波型關于震源分別上下、左右對稱，符合地震波傳播的基本規律。

為了更清楚展示波場分離的振幅和相位特征，在所有波場快照的1.5 km處（圖2（a）中白線示意位置）抽取單道。圖3（a）～（c）分別為耦合波場、縱波場和橫波場水平分量單道信息，圖3（d）～（f）分別為耦合波場、縱波場和橫波場垂向分量單道信息。通過對比可以看出，分離得到的縱波和橫波振幅沒有衰減為耦合波場的 1 ?V ?P ?， 1 V ?S ，相位也沒有延遲，均基本上保留了與耦合波場一樣的振幅和相位特性。

由此可見，基于解耦方程的彈性波場可以較為準確地分離出縱波和橫波，在傳播過程中可以保持其振幅和相位基本不變，無需再額外進行振幅補償和相位校正。相比于傳統基于Helmholtz分解對波場取散度和旋度的方法，該方法更簡易、涉及的步驟更少。

2.2 侵入花崗型和近代火山型干熱巖

侵入花崗型和近代火山型作為干熱巖的兩種常見類型，具有較大的研究價值，發展一種高精度成像方法，對地熱能的探測和開發具有重要意義。已經對基于解耦方程的彈性波場進行了分析，通過數值模擬論證了其優點，然后使用該方程分別對面積為10 km×7 km的侵入花崗模型和近代火山模型目標區域進行逆時偏移成像研究，所采用的偏移速度模型均由真實速度模型進行5點平滑后所得。

2.2.1 侵入花崗型

侵入花崗型干熱巖是由以巖漿巖深層侵入為主的強烈構造活動導致。圖4為侵入花崗模型的縱波速度分布。該模型含2條正斷層，地層間多為不整合接觸關系，較多位置出現尖滅現象，構造比較復雜。模型中最小速度為2 000 m/s，最大速度為5 500 m/s?；◢弾r基巖位于5 km深度以下，有3條主支呈無規則形狀向上侵入5套地層可至3 km。

對模型采樣1 001×701個，橫縱空間采樣間隔均為10 m，時間采樣率為0.3 ms，時間采樣21 000個，雷克子波主頻為20 Hz，以200 m的炮間距在地 表共布設51個炮點。采用全接收的觀測系統，檢波點共1 001個，檢波點間隔為10 m。

圖5（a）、（b）為耦合波場成像結果，圖5（c）～（f）為PP、PS、SP、SS成像結果。整體來看，圖5（c）～（f） 的同相軸能量分布更均衡，偏移假象更少，成像效果較好，可以較為清晰地實現侵入花崗巖體的基本構造成像。

為了方便分析解耦方程對成像結果精度的影響，截取圖5中黑色標識框圈定的區域進行分析，將其左側小黑框圈定區域記為局部1（圖6），中間黑框圈定的區域記為局部2（圖7）。

圖6（a）、（b）和圖7（a）、（b）為耦合波場在水平方向和垂直方向的局部成像結果，在白色箭頭所指向的位置處，可以看到圖6（a）、（b）;圖7（a）中出現目標體同相軸明顯錯段、能量微弱等現象，圖7（b）中出現明顯的同相軸合并，偏移假象等現象。都沒有很好地對目標體成像，這可能會給后續的解釋工作帶來困難。而在圖6（c）～（f）和圖7（c）～（f）的同位置處，尤其是圖6（f）和圖7（f），可明顯地看到同相軸不連續位置更少、同相軸能量更均衡、偏移假象更少、同相軸分辨率更高，可更準確地對目標體構造成像結果進行分析。

2.2.2 近代火山型

近代火山型干熱巖主要是由深層高溫巖漿向上涌出形成。圖8為近代火山型的縱波速度分布情況。該模型左側出現地塹構造，產生一條正斷層和一條逆斷層，含多套不規則形狀地層，構造復雜。最小速度為2 000 m/s，最大速度為5 500 m/s。高溫火山基巖位于5 km至深處，有2條主支左右入侵其他地層，巖漿房緊挨于火山基巖右旁，其中巖漿向上涌出至1 km處。

對該模型采樣也是1 001×701個，橫縱空間采樣間隔均為10 m，時間采樣率為0.3 ms，時間采樣 20 000 個，子波主頻為20 Hz，以200 m的炮間距在地表共布設51個炮點，同樣采用的是檢波點間隔為10 m的全接收觀測系統。

由于巖株的速度相比其附近地層的要大很多，所以在1 km深度附近的巖株頂部發生了強烈反射，使得圖9在該位置表現出同相軸的能量較強。相應地，透射波的能量相對減少，該區域至深處地震波能量較弱，目標體同相軸較為不明顯。為了更好地分析基于解耦方程進行逆時偏移的成像效果，選擇的區域最好是目標體同相軸能量較強的區域，故選擇圖9中黑色框圈定區域進行成像結果的局部分析，以進一步說明成像精度。

圖10（a）、（b）為P-S波耦合的成像結果的局部，在箭頭所指向的區域，目標體的同相軸能量較弱，有些地方幾乎沒有出現同相軸，成像結果不好，這是耦合波成像的常見缺陷。圖10（c）～（f）為PP、PS、SP和SS波成像結果局部，相對于圖10（a）、（b），同位置成像效果明顯變好，同相軸有了較好的呈現，分辨率提高很多。

3 結 論

（1）基于解耦方程的波場分解避免了傳統Helmholtz分解對波場取散度和旋度導致的振幅減小、相位推遲的問題，無需后續額外校正，可在實現高精度逆時偏移成像的同時，將地震資料處理步驟減少。

（2）侵入花崗型和近代火山型兩個典型干熱巖模型的PP、PS、SP和SS成像結果說明，較于耦合波場成像，基于解耦方程分離縱橫波得到的彈性波成像結果普遍更清晰，同相軸錯段位置減少和信息更豐富，深部能量更均衡，分辨率更高。

參考文獻 ：

［1］ ?劉偉，李慧，張偉.基于兩種波場分離方法的彈性波逆時偏移對比研究［J］.地球物理學進展，2020，35（3）：1000-1009.

LIU Wei， LI Hui， ZHANG Wei. Comparative study on elastic reverse time migration based on two methods of waves field separation［J］. Progress in Geophysics， 2020，35（3）：1000-1009.

［2］ ?GU B， LI Z， MA X， et al. Multi-component elastic reverse time migration based on the P-and S-wave separated velocity–stress equations ［J］. Journal of Applied Geophysics， 2015，112：62-78.

［3］ ?DU Q， ZHAO Q， LI Q， et al. A new decoupling and elastic propagator for efficient elastic reverse time migration［J］. Geophysics， 2020，85（5）：A31-A36.

［4］ ?GU B， HUANG J， HAN J， et al. Least-squares inversion-based elastic reverse time migration with PP-and PS-angle-domain common-imaging gathers［J］. Geophysics， 2021，86（1）：S29-S44.

［5］ ?CHEN Y， FOMEL S. 3D true-amplitude elastic wave-vector decomposition in heterogeneous anisotropic media［J］. Geophysics， 2023，88（3）：1-45.

［6］ ?ZHU H. Elastic wavefield separation based on the Helmholtz decomposition elastic wave field separation［J］. Geophysics， 2017，82（2）：S173-S183.

［7］ ?YANG J， ZHU H， HUANG J， et al. 2D isotropic elastic Gaussian-beam migration for common-shot multicomponent records［J］. Geophysics， 2018，83（2）：S127-S140.

［8］ ?GU B， REN Z， LI Q， et al. An application of vector wavefield decomposition to 3D elastic reverse time migration and field data test［J］. Computers & Geosciences， 2019，131：112-131.

［9］ ?SUN R， MCMECHAN G A. Scalar reverse-time depth migration of prestack elastic seismic data［J］. Geophysics， 2001，66（5）：1519-1527.

［10］ ?DELLINGER J， ETGEN J. Wave-field separation in two-dimensional anisotropic media［J］. Geophysics， 1990，55（7）：914-919.

［11］ ?YAN J， SAVA P. Elastic wave-mode separation for VTI media［J］. Geophysics， 2009，74（5）：WB19-WB32.

［12］ ?馬德堂，朱光明.彈性波波場P波和S波分解的數值模擬［J］.石油地球物理勘探，2003，38（5）：482-486.

MA Detang， ZHU Guangming. P- and S-wave separated elastic wave equation numerical modeling［J］. Oil Geophysical Prospecting， 2003，38（5）：482-486.

［13］ ?李振春，張華，劉慶敏，等.彈性波交錯網格高階有限差分法波場分離數值模擬［J］.石油地球物理勘探，2007，42（5）：510-515.

LI Zhenchun， ZHANG Hua， LIU Qingmin， et al. Numeric simulation of elastic wavefield separation by staggering grid high-order finite-difference algorithm［J］. Oil Geophysical Prospecting， 2007，42（5）：510-515.

［14］ ?ZHANG J， TIAN Z， WANG C. P-and S-wave separated elastic wave equation numerical modeling using 2D staggered-grid［R］. SEG 2104，2007.

［15］ ?XIAO X， LEANEY W S. Local vertical seismic profiling （VSP） elastic reverse-time migration and migration resolution： salt-flank imaging with transmitted P-to-S waves［J］. Geophysics， 2010，75（2）：S35-S49.

［16］ ??WANG W， MCMECHAN G A， ZHANG Q. Comparison of two algorithms for isotropic elastic P and S vector decomposition［J］. Geophysics， 2015，80（4）：T147-T160.

［17］ ?ZHANG W， SHI Y. Imaging conditions for elastic reverse time migration［J］. Geophysics， 2019，84（2）：S95-S111.

［18］ ?李振春，雍鵬，黃建平，等.基于矢量波場分離彈性波逆時偏移成像［J］.中國石油大學學報（自然科學版），2016，40（1）：42-48.

LI Zhenchun， YONG Peng， HUANG Jianping， et al. Elastic wave reverse time migration based on vector wavefield separation ［J］. Journal of China University of Petroleum （Edition of Natural Science）， 2016，40（1）：42-48.

［19］ ?YONG P， HUANG J， LI Z， et al. Elastic-wave reverse-time migration based on decoupled elastic-wave equations and inner-product imaging condition ［J］. Journal of Geophysics and Engineering， 2016，13（6）：953-963.

［20］ ?孫金聲，廖波，王金堂，等.溫度場擾動對天然氣水合物相態的影響機制［J］.中國石油大學學報（自然科學版），2023，47（5）：115-121.

SUN Jinsheng， LIAO Bo， WANG Jintang， et al. Mechanism on influence of temperature field perturbation on phase state of natural gas hydrate［J］. Journal of China University of Petroleum （Edition of Natural Science），2023，47（5）：115-121.

［21］ ?李玉坤，于文廣，李玉星，等.超聲臨界折射縱波測量應力的溫度影響［J］.中國石油大學學報（自然科學版），2021，45（2）：134-140.

LI Yukun， YU Wenguang， LI Yuxing， et al. Effect of temperature on measurement of stress by ultrasonic critical refractive longitudinal wave［J］. Journal of China University of Petroleum （Edition of Natural Science）， 2021，45（2）：134-140.

［22］ ?王瑩瑩，劉書杰，張崇，等.水下采油樹過流通道溫變規律和計算方法［J］.中國石油大學學報（自然科學版），2022，46（3）：148-157.

WANG Yingying， LIU Shujie， ZHANG Chong， et al. Law of temperature change and calculation method in subsea christmastree internal overflow channel［J］. Journal of China University of Petroleum（Edition of Natural Science）， 2022，46（3）：148-157.

［23］ ?VIRIEUX J. P-SV wave propagation in heterogeneous media： velocity-stress finite-difference method［J］. Geophysics， 1986，51（4）：889-901.

（編輯 沈玉英）