深层地热花岗岩体地震波数值模拟及传播机制

黄建平 刘英辉 李伟 张盟勃 王扬州 杨永红

摘要 ：深層地热资源是一种可再生的、蕴藏巨大能量的清洁能源，但其地球物理响应特征不明确，中深层地热探测成功率较低。为了研究地震波在深层地热岩体中的传播规律与波场特征，建立两种深层地热花岗岩模型，并使用等效交错网格有限差分法实现声波与弹性波的数值模拟。结果表明：地热花岗岩体速度在温度的影响下要高于围岩的，会产生高速屏蔽现象，使得透射波能量变弱，限制了地热岩体下部地震波传播能量；相比于声波，弹性波具有更丰富的波场信息，波型与能量转换使得弹性波地震记录也比声波地震记录复杂。

关键词 ：深层地热岩； 地震波传播机制； 数值模拟； 有限差分模拟

中图分类号 ：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］ 。相比于发达国家，中国对于深层地热资源勘查开发及其技术还处于初级阶段，没有形成系统的勘探开发方法。现阶段主要面临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.

（编辑 沈玉英）