面向深層干熱巖體的全波形速度反演建模方法

楊繼東 于由財 劉朋 高建明 黃建平 楊永紅

摘要 ：隨著勘探深度增加，地震信號衰減嚴重，常規速度建模方法精度較低，無法滿足深層干熱巖體高精度勘探的需求。為提高速度建模的精度，采用全波形反演方法，但受局部優化算法的限制，其應用存在收斂速度慢、反演深度淺以及容易陷入局部極值等問題。為此，對梯度應用預條件處理和平滑處理，并使用共軛梯度優化算法，以解決地下照明不均勻問題，提高反演深度、精度和收斂速度。同時，為緩解局部極小值問題，引入局部相似性全波形反演方法來更新初始模型，并對構建的干熱巖模型進行數值測試。結果表明：即使在初始模型極不準確的情況下，該方法仍能避免周波跳躍的不利影響，實現穩健的迭代更新；該方法能夠顯著增加反演深度和精度，并且對高速巖體有較好的刻畫，最終能獲得淺、中、深全層系高精度速度模型；提出的全波形反演方法為深層高溫花崗巖體的勘探提供了一個切實可行的建模流程。

關鍵詞 ：干熱巖； 高溫花崗巖體； 全波形反演； 深層速度建模； 局部相似性

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

引用格式 ：楊繼東，于由財，劉朋，等.面向深層干熱巖體的全波形速度反演建模方法［J］.中國石油大學學報（自然科學版），2024，48（1）：70-76.

YANG Jidong， YU Youcai， LIU Peng， et al. Velocity building using full waveform inversion for deep hot dry rock［J］. Journal of China University of Petroleum（Edition of Natural Science）， 2024，48（1）：70-76.

Velocity building using full waveform inversion for deep hot dry rock

YANG Jidong 1， YU Youcai 1， LIU Peng 2， GAO Jianming 2， ?HUANG Jianping 1， YANG Yonghong 3

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

2.Shandong Energy Group South America Limited， Qingdao 266580， China;

3.Shengli Oilfield Exploration and Development Research Institute， ?SINOPEC， Dongying 257000， China）

Abstract ：Severe seismic signal attenuation with increasing exploration depth renders conventional velocity building methods less accurate for high-precision exploration of deep hot dry rock （HDR） bodies. To address this challenge， the full waveform inversion （FWI） method is used to enhance the accuracy of velocity building. However， FWI encounters issues like slow convergence， limited inversion depth， and susceptibility to local optimal solutions due to constraints imposed by local optimization algorithms. To mitigate these challenges， we introduce a novel approach that utilizes preconditioning and smoothing gradients in conjunction with the conjugate gradient optimization algorithm. This strategy aims to rectify unbalanced illumination issues and speed up convergence. Additionally， to alleviate the issue of local minima， a local coherence misfit is integrated into FWI to update the velocity model. Numerical simulations conducted on typical HDR models show that the proposed method is capable of circumventing the negative effects of cycle skipping issues， ensuring stable iterative updating even with highly inaccurate initial model. Furthermore， the method significantly improves the depth and accuracy of inversion， providing a more precise depiction of high-velocity rock bodies， eventually obtaining high-accuracy velocity models in both shallow and deep layers. The FWI method proposed in this study provides a practical and efficient velocity building tool for detecting deep high-temperature granite bodies.

Keywords ： hot dry rock; high-temperature granite body; full waveform inversion; velocity building for deep structure; local coherence

“碳達峰”“碳中和”概念的提出，標志著中國將繼續堅定不移地推進能源革命，加快構建清潔低碳、安全高效的能源體系。干熱巖是深層地熱資源重要的賦存形式之一，具有資源量大、分布廣、無污染、安全、利用率高等優點，是今后能源開發的重要方向 ?［1］ 。王鈞 ?［2］ 通過對地溫場分布的研究，闡明了地溫場的形成原因和控制因素。汪集旸等 ?［3］ 對中國大陸干熱巖地熱資源的潛力進行了評估，并圈定出了有利靶區。甘浩男等 ?［4］ 將中國干熱巖資源的賦存類型分為高放射性產熱型、沉積盆地型、近代火山型和強烈構造活動帶型，并分別對其成因機制進行了研究。由于干熱巖儲層埋藏較深（通常大于3 km），溫度高（普遍高于200 ℃），硬度大，因此開發風險大，成本高，這決定了地球物理勘探的重要性 ?［5-6］ 。常用于地熱勘探的地球物理方法有重力勘探 ?［7-8］ 、磁法勘探 ?［9-10］ 、電法勘探 ?［11-12］ 、地震勘探 ?［13-15］ 和遙感 ?［16-18］ 等。其中地震勘探具有高精度、高分辨率、大探測深度的優勢，可以精確推測斷層位置、埋深和產狀，進一步可以通過波速分布圈定地熱資源的位置。然而，隨著探測深度的增加以及對成像精度的要求越來越高，傳統速度建模方法如疊加速度分析、偏移速度分析和層析成像等已經難以滿足地震資料處理和解釋的需求。全波形反演可以綜合利用運動學和動力學信息，通過不斷使合成數據擬合匹配觀測數據而更新優化速度模型，來獲得可以準確描述地下速度分布的高精度模型。Lailly ?［19］ 和Tarantola ?［20］ 在聲學介質中直接將觀測與合成數據進行擬合，革新了地球內部速度結構的反演。Tromp ?［21］ 總結形成了伴隨層析成像的框架，可以自由選擇包括相位、振幅、波形在內的目標函數。之后許多專家學者在目標函數、梯度、尋優方法上進行了一些改進 ?［22-26］ ，在實際資料應用中也取得了很好的效果。祝賀君等 ?［27］ 將地震學全波形反演進展進行了分類與總結歸納。筆者針對深層高溫干熱巖體建模問題，提出一種面向深層干熱巖體的全波形速度反演建模方法。首先，根據中國地殼結構背景和前人研究，構建兩種干熱巖速度模型。然后，在低頻段采用局部相似性全波形反演 ?［28］ 對低精度的初始模型進行更新，以緩解常規全波形反演中可能出現的周波跳躍問題。接下來，采用分頻多尺度常規全波形反演策略細化速度結構。最終，得到高精度的深層速度模型。

1 方法原理

1.1 常規全波形反演

在常規聲波全波形反演中，觀測數據與合成數據差異的L2范數通常用作目標泛函：

f（v）= 1 2 ∑ r ∫［d ??syn ?（x r，t）-d ??obs ?（x r，t）］ 2 d t.

（1）

式中，f（v）為目標泛函；v為模型速度參數；r為檢波器；d ??syn ?（x r，t）與d ??obs ?（x r，t）分別為接收位置x r處記錄時間為t的合成數據與觀測數據。

梯度可以通過求取目標泛函相對模型速度的導數來推導：

f（v） ?v =∑ r 0 ∫ ?d ??syn ?（x r，t） ?v d ??adj ?（x r，t） d t.

（2）

式中，d ??adj ?（x r，t）為常規全波形反演的伴隨源，為合成數據與觀測數據的殘差，表示為

d ??adj ?（x r，t）=d ??syn ?（x r，t）-d ??obs ?（x r，t）. （3）

1.2 局部相似性波形反演

為了緩解因初始模型精度低而導致的“周波跳躍”問題，在低頻段反演時采用局部相似性目標泛函的全波形反演。其目標泛函表示為

f（v）=- 1 N ∑ r ∫C（x r，t） d t. （4）

式中，N為所有數據的總采樣點數；C（x r，t）為由局部合成和觀測記錄計算的相關系數，其取值范圍是［-1，1］。由于全波形反演的優化算法是最小化函數值，而此時的目標是使相關系數最大化，因此將相關系數求和并取均值后再與負一相乘，這樣就可以直接替換常規目標泛函應用在全波形反演框架中。相關系數C（x r，t）的計算方法如下：

C（x r，t）=

syn ?（k，h，x r，t） ???obs ?（k，h，x r，t） d k d h

syn ?（k，h，x r，t） 2 d k d h ? ???obs ?（k，h，x r，t） 2 d k d h ??. ?（5）

式中， ???syn ?（k，h，x r，t）和 ???obs ?（k，h，x r，t）分別為去除均值后的局部尺度數據。局部尺度數據D ??syn ?（k，h，x r，t）與D ??obs ?（k，h，x r，t）是用二維局部窗口獲得，k和h分別是局部窗口的長和寬，獲取局部數據表示如下：

D ??syn ?（k，h，x r，t）=d ??syn ?（x r，t）g（k，h），

D ??obs ?（k，h，x r，t）=d ??obs ?（x r，t）g（k，h）. ?（6）

式中，g（k，h）為二維 Gabor 正變換。為了將局部尺度數據恢復成真實數據尺度，相應的逆變換表示為

d ??syn ?（x r，t）=D ??syn ?（k，h，x r，t）g（k，h） d k d h，

d ??obs ?（x r，t）=D ??obs ?（k，h，x r，t）g（k，h） d k d h. ?（7）

對應的梯度可以表示為

f（v） ?v =- 1 N ∑ r ∫ ?d ??syn ?（x r，t） ?v d ??adj ?（x r，t） d t. （8）

式（8）中伴隨源d ??adj ?（x r，t）表示為

d ??adj ?（x r，t）=? ? adj ?（k，h，x r，t）g（k，h） d k d h. （9）

式中， ???adj ?（k，h，x r，t）為去除均值后的局部伴隨源。局部伴隨源D ??adj ?（k，h，x r，t）表示為

D ??adj ?（k，h，x r，t）=

2 ??syn ??（k，h，x r，t） d k d h ??? ?2 ??obs ??（k，h，x r，t） d k d h ?????obs ?（k，h，x r，t）- C（x r，t） ???syn ?（k，h，x r，t）. （10）

1.3 預條件共軛梯度法

全波形反演的優化算法是通過從初始模型逐步迭代優化來實現的，迭代的序列可以表示為

v ?k+1 =v k+α kδv k. （11）

式中，v ?k+1 和v k分別為第k+1和k次迭代的模型速度參數；α k和δv k分別為第k次用拋物線擬合法估計的最佳步長 ?［29］ 和下降方向。在全波形反演中常用的基于梯度的算法是共軛梯度法，其下降方向為

δv 0=-g 0=- ?f（v 0） ?v 0 ?，

δv ?k+1 =-g ?k+1 +β ?k+1 δv k.

（12）

式中，g ?k+1 為第k+1次迭代的梯度；β ?k+1 為第k+1次的加權因子，可以采用不同的方式來計算，筆者使用的計算方式可表示為

β ?k+1 = ?g ??T ??k+1 （ g ??k+1 - g ?k） ??g ??T ?k ?g ?k ?. （13）

對于預條件共軛梯度法，使用對角漢森矩陣對下降方向進行修改：

δv 0=- H ??-1 ?0 g ?0，

δv ?k+1 =- H ??-1 ??k+1

g ??k+1 +β ?k+1 δv k. ?（14）

式中， H ??-1 ??k+1 為第k+1次迭代的逆漢森矩陣。相應地將β ?k+1 修改為

β ?k+1 = （ H ??-1 ??k+1 ?g ??k+1 ） ?T （ H ??-1 ??k+1 ?g ??k+1 - H ??-1 ?k g ?k） （ H ??-1 ?k g ?k） ?T （ H

-1 ?k g ?k） ?. （15）

2 數值模擬

建立并采用兩個干熱巖體模型（近代火山型與高放射性產熱型），從不精確的一維初始速度模型開始，使用所提出的方法得到了高分辨率的反演結果，驗證了此方法的優點和適應性。

2.1 近代火山型干熱巖模型

近代火山型干熱巖是一種與火山活動密切相關的特殊巖石類型，巖性多為花崗巖，具有高溫度、高硬度、高密度和高熔點等特點。近代火山型干熱巖模型如圖1（a）所示。在該模型中，伸張應力導致形成多個正斷層，為地下熱流體提供了從深部到地表的上升通道。

巖漿從地幔上涌，形成了巖漿房，這些巖漿房內未冷卻的巖漿為附近巖石提供持續熱源，使巖基被持續加熱，形成一個完整、高效的地熱系統。對應的一維初始模型如圖1（b）所示。模型離散為281×651網格點，網格間距為25 m。地震記錄采用交錯網格有限差分方法計算，55炮均勻分布在地表上，使用全接收觀測系統，記錄時間為9.2 s，時間采樣間隔為2 ms。

觀測記錄和合成記錄分別如圖2（a）、（b）所示。本文中在低頻段（0～3 Hz）使用局部相似性全波形反演更新初始模型，來緩解因初始模型不準確而可能出現的周波跳躍問題。第一次迭代的伴隨源如圖2（d）所示，可以看出其不再是觀測數據和合成數據的簡單相減（圖2（c）），而是根據局部相似性自適應的平衡波形，首先擬合早至波和大偏移距回轉波部分，并隨著模型的恢復，不斷自動增加擬合范圍，實現自適應從低波數擬合到高波數擬合，從而有效緩解周波跳躍問題。此外，該方法可以有效補償大偏移距范圍下的能量損失，有助于恢復深部速度結構。

圖3顯示了常規的梯度以及使用對角漢森預條件和平滑濾波處理后的梯度。經過處理的梯度與常規梯度相比，壓制了高頻分量使得整體更加平滑，這有利于緩解周波跳躍問題，而且具有更深的照明范圍，可以有效增加反演的深度，加快模型的更新和收斂。

常規全波形反演在0～3 Hz頻段上進行30次迭代的結果如圖4（a）所示，該結果顯示模型進行了不連續的點狀、條狀更新，這會使之后的高頻段反演變得極不穩定。為了克服這一問題，本文中首先采用局部相似性全波形反演對初始模型進行更新，在0～3 Hz頻段上進行了60次迭代，成功地恢復了模型的宏觀結構，結果如圖4（b）所示。為進一步提高模型精度，本文中使用分頻多尺度全波形反演繼續更新圖4（b）所示的速度模型，最終獲得了圖4（c）所示的反演結果。該結果已經有效地恢復了淺、中、深層以及巖漿高速體的速度結構，表明本文方法對巖漿高速侵入體的速度建模具有很好的適應性。

2.2 高放射性產熱型干熱巖模型

高放射性產熱型干熱巖是一種富含鈾、釷等放射性元素的礦物資源，具有高輻射、高溫度和高產熱等特點。該類巖石多為酸性花崗巖，主要由鉀長石等礦物組成。高放射性產熱型干熱巖模型如圖5（a）所示。在此模型中地表呈現地塹構造特征，而其深部是明顯的背斜結構。背斜促進了深源巖漿的上升，當這些巖漿接觸到冷的地殼巖體時，它們不僅冷卻和結晶，還導致了熱和化學成分的交換。這種交互過程加強了系統的地熱特性。最關鍵的熱源來自于巖漿中豐富的放射性元素，它們通過衰變持續地釋放熱量。對應的一維初始模型如圖5（b）所示。模型離散為281×601網格點，橫、縱向網格間距都為25 m。地震記錄采用交錯網格有限差分方法計算，51炮均勻分布在地表上，使用全接收觀測系統，記錄時間為7.5 s，時間采樣間隔為1.5 ms。

常規全波形反演在0～3 Hz頻段進行30次迭代更新的結果如圖6（a）所示。由于初始模型不準確，導致模型的更新陷入了周波跳躍問題，這會對進一步的高頻段反演造成嚴重的干擾，最終導致錯誤的結果和解釋。因此本文中仍采用局部相似性全波形反演方法對初始模型進行更新，在低頻段（0～3 Hz）進行50次迭代更新，結果如圖6（b）所示。該方法成功恢復了模型深部和中部目標層的速度輪廓。接下來使用分頻多尺度全波形反演策略對圖6（b）所示速度模型做進一步優化，得到如圖6（c）所示的最終反演結果。結果表現出較高的精度，能夠非常有效地揭示地下高放射性花崗巖體的輪廓，同時也證明了該方法對深層花崗巖體的速度建模有很好的適應性。

3 結束語

首先采用對周波跳躍問題不敏感的局部相似性全波形反演方法對初始模型進行更新，以使其宏觀結構與真實模型相匹配，并做進一步細化以使接下來的分頻多尺度常規全波形反演可以穩健進行。在優化算法方面，本研究采用預條件共軛梯度算法。該算法不但可以增加反演的深度和精度，還能夠提高反演的收斂速度并節省計算成本。本文中提出的方法在兩個干熱巖模型的算例中進行了驗證，其結果表明，即使在初始模型極不準確的情況下，該方法仍能避免周波跳躍的不利影響，最終獲得高精度的深層速度模型。該方法在干熱巖勘探方面具有廣闊的應用前景。

參考文獻 ：

［1］ ?藺文靜，王貴玲，邵景力，等.我國干熱巖資源分布及勘探：進展與啟示［J］.地質學報，2021，95（5）：1366-1381.

LIN Wenjing， WANG Guiling， SHAO Jingli， et al. Distribution and exploration of hot dry rock resources in China： progress and inspiration［J］. Acta Geologica Sinica， 2021，95（5）：1366-1381.

［2］ ?王鈞.東南沿海地區地溫場的形成及其分布規律［J］.地震地質，1985，7（1）：49-58.

WANG Jun. Distribution and formation of the geothermal field along the coast， southeastern China［J］. Seismology and Geology， 1985，7（1）：49-58.

［3］ ?汪集旸，胡圣標，龐忠和，等.中國大陸干熱巖地熱資源潛力評估［J］.科技導報，2012，30（32）：25-31.

WANG Jiyang， HU Shengbiao， PANG Zhonghe， et al.Estimate of geothermal resources potential for hot dry rock in the continental area of China［J］. Science and Technology Review， 2012，30（32）：25-31.

［4］ ?甘浩男，王貴玲，藺文靜，等.中國干熱巖資源主要賦存類型與成因模式［J］.科技導報，2015，33（19）：22-27.

GAN Haonan， WANG Guiling， LIN Wenjing， et al. Research on the occurrence types and genetic models of hot dry rock resources in China ［J］. Science and Technology Review， 2015，33（19）：22-27.

［5］ ?FU G， PENG S， WANG R， et al. Seismic prediction and evaluation techniques for hot dry rock exploration and development［J］. Journal of Geophysics and Engineering， 2022，19（4）：694-705.

［6］ ?曾昭發，陳雄，李靜，等.地熱地球物理勘探新進展［J］.地熱能，2012（5）：3-12.

ZENG Zhaofa， CHEN Xiong， LI Jing， et al. Advancement of geothermal geophysics exploration［J］. Geothermal Energy， 2012（5）：3-12.

［7］ ?羅國平.重力勘探在城市區地熱調查中的應用［J］.中國煤田地質，2000，12（4）：68-72.

LUO Guoping. The application of gravity prospecting in geothermal survey at city and its suburbs［J］. Coal Geology of China， 2000，12（4）：68-72.

［8］ ?劉振華，李世峰，楊特波，等.綜合物探技術在邯鄲地熱田勘查中的應用［J］.工程地球物理學報，2013，10（1）：111-116.

LIU Zhenhua， LI Shifeng， YANG Tebo， et al. The application of integral geophysical survey technology in geothermal exploration of Handan area［J］. Journal of Engineering Geophysics， 2013，10（1）：111-116.

［9］ ?SOENGKONO S， HOCHSTEIN M P. Magnetic anomalies over the Wairakei geothermal field， central north island， New Zealand［J］. Geothermal Resources Council Transactions， 1992，16：273-278.

［10］ ?孫二虎，陳愛珍，劉鴻福，等.測氡法-磁法在祁縣王村地熱勘測中的應用［J］. 太原理工大學學報， 2004，35（3）：307-310.

SUN Erhu， CHEN Aizhen， LIU Hongfu， et al. The application of radon and magnetic survey in geothermal exploration in Wangcun， Qixian County［J］. Journal of Taiyuan University of Technology， 2004，35（3）：307-310.

［11］ ?SINGH S B， DROLIA R K， SHARMA S R， et al. Application of resistivity surveying to geothermal exploration in the Puga Valley， India［J］. Geoexploration， 1983，21（1）：1-11.

［12］ ?郭守鋆，李百祥，周小波，等.電法在甘子河斷裂對流型地熱資源勘查中的應用［J］.物探與化探，2013，37（2）：229-232.

GUO Shoujun， LI Baixiang， ZHOU Xiaobo， et al. The application of the remote-sensing and electrical methods to geothermal water exploration［J］. Geophysical and Geochemical Exploration， 2013，37（2）：229-232.

［13］ ?BARIA R， MICHELET S， BAUMGRTNER J， et al. Microseismic monitoring of the world s largest potential HDR reservoir［C/OL］//Proceedings of Twenty-Ninth Workshop on Geothermal Reservoir Engineering. Stanford University， Stanford， California， January 26-28， 2004［2023-03-12］. https：//pangea.stanford.edu/ERE/pdf/IGAstandard/SGW/2004/Baria.pdf.

［14］ ?KARASTATHIS V K， PAPOULIA J， Di FIORE B， et al. Deep structure investigations of the geothermal field of the North Euboean Gulf， Greece， using 3-D local earthquake tomography and curie point depth analysis［J］. Journal of Volcanology and Geothermal Research， 2011，206（3/4）：106-120.

［15］ ?孫黨生，雷煒，李洪濤，等.高分辨率地震勘探在地熱資源勘查中的應用［J］.勘察科學技術，2002（6）：55-59.

SUN Dangsheng， LEI Wei， LI Hongtao， et al. Application of high-resolution on seismic exploration method to the prospecting of geothermal resources［J］. Exploration Science and Technology， 2002（6）：55-59.

［16］ ?許軍強，白朝軍，劉嘉宜.遙感技術支持下的長白山火山區地熱預測研究［J］.自然資源遙感，2009，20（1）：68-71.

XU Junqiang， BAI Chaojun， LIU Jiayi. Geothermalresource prognosis based on remote sensing technology in Changbaishan volcanic area［J］. Remote Sensing for Land and Resources， 2009，20（1）：68-71.

［17］ ?van der MEER F， HECKER C， van RUITENBEEK F， et al. Geologic remote sensing for geothermal exploration：a review［J］. International Journal of Applied Earth Observation and Geoinformation， 2014，33：255-269.

［18］ ?NORINI G， GROPPELLI G， SULPIZIO R， et al. Structural analysis and thermal remote sensing of the Los Humeros Volcanic Complex： implications for volcano structure and geothermal exploration［J］. Journal of Volcanology and Geothermal Research， 2015，301：221-237.

［19］ ?LAILLY P， SANTOSA F. Migration methods： partial but efficient solutions to the seismic inverse problem［C/OL］// Inverse problems of acoustic and elastic waves， Philadelphia， Pennsylvania， US， Jan， 1984［2023-01-22］.https：//searchworks.stanford.edu/view/1622537.

［20］ ?TARANTOLA A. Inversion of seismic reflection data in the acoustic approximation［J］. Geophysics， 1984，49（8）：1259-1266.

［21］ ?TROMP J， TAPE C， LIU Q. Seismic tomography， adjoint methods， time reversal and banana-doughnut kernels［J］. Geophysical Journal International， 2005，160（1）：195-216.

［22］ ?CHOI Y， ALKHALIFAH T. Application of multi-source waveform inversion to marine streamer data using the global correlation norm［J］. Geophysical Prospecting， 2012，60：748-758.

［23］ ?WARNER M， GUASCH L. Adaptive waveform inversion：theory［J］. Geophysics， 2016，81（6）：R429-R445.

［24］ ?YANG J， ZHU H， LI X， et al. Estimating P wave velocity and attenuation structures using full waveform inversion based on a time domain complex-valued viscoacoustic wave equation： the method［J］. Journal of Geophysical Research： Solid Earth， 2020，125（6）：e2019JB019129.

［25］ ?HU Y， HAN L G， XU Z， et al. Demodulation envelope multi-scale full waveform inversion based on precise seismic source function［J］. Chinese Journal of Geophysics， 2017，60（3）：1088-1105.

［26］ ?YANG J， XU J， HUANG J， et al. The connection of velocity and impedance sensitivity kernels with scattering-angle filtering and its application in full waveform inversion［J］. Frontiers in Earth Science， 2022，10：961750.

［27］ ?祝賀君，劉沁雅，楊繼東.地震學全波形反演進展［J］.地球與行星物理論評，2023，54（3）：287-317.

ZHU Hejun， LIU Qinya， YANG Jidong. Recent progress on full waveform inversion［J］. Progress in Earth and Planetary Science， 2023，54（3）：287-317.

［28］ ?YU Y， YANG J， HUANG J， et al. Full waveform inversion using a high-dimensional local-coherence misfit function［J］. IEEE Transactions on Geoscience and Remote Sensing， 2023，61：1-8.

［29］ ?NOCEDAL J， WRIGHT S J. Numerical optimization［M］.2nd edi. New York：Springer， 2006：101-134.

（編輯 修榮榮）