面向深层干热岩体的全波形速度反演建模方法

杨继东 于由财 刘朋 高建明 黄建平 杨永红

摘要 ：隨着勘探深度增加，地震信号衰减严重，常规速度建模方法精度较低，无法满足深层干热岩体高精度勘探的需求。为提高速度建模的精度，采用全波形反演方法，但受局部优化算法的限制，其应用存在收敛速度慢、反演深度浅以及容易陷入局部极值等问题。为此，对梯度应用预条件处理和平滑处理，并使用共轭梯度优化算法，以解决地下照明不均匀问题，提高反演深度、精度和收敛速度。同时，为缓解局部极小值问题，引入局部相似性全波形反演方法来更新初始模型，并对构建的干热岩模型进行数值测试。结果表明：即使在初始模型极不准确的情况下，该方法仍能避免周波跳跃的不利影响，实现稳健的迭代更新；该方法能够显著增加反演深度和精度，并且对高速岩体有较好的刻画，最终能获得浅、中、深全层系高精度速度模型；提出的全波形反演方法为深层高温花岗岩体的勘探提供了一个切实可行的建模流程。

关键词 ：干热岩； 高温花岗岩体； 全波形反演； 深层速度建模； 局部相似性

中图分类号 ：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.

（编辑 修荣荣）