一种高温介质地震波传播模拟方法

摘要

本发明属于地震波模拟领域，具体提供一种高温介质地震波传播模拟方法。本发明旨在解决传统的地震波正演模型为常温介质导致的模拟结果不准确的问题。本发明的模拟方法包括以下步骤：利用边界元法将热弹性波动方程转化为边界元形式方程；利用边界元形式方程求得地震波在均匀岩石物理模型区域的边界的应力值和位移值；将边界的应力值和位移值代入边界元形式方程求得地震波在均匀岩石物理模型内区域的应力值和位移值；根据均匀岩石物理模型区域的边界和内区域的全部应力值和位移值，确定地震波的波形记录。本发明可以准确考察高温对地震波传播的影响，对深层及超深层油气勘探领域拓展提供理论基础，对地球深部观测数据异常现象提供合理的解释。

权利要求

1.一种高温介质地震波传播模拟方法，其特征在于，包括以下步骤：步骤S10、利用边界元法将热弹性波动方程转化为边界元形式方程；

步骤S20、利用边界元形式方程求得地震波在均匀岩石物理模型区域的边界的应力值和位移值；

步骤S30、将所述边界的应力值和位移值代入边界元形式方程求得地震波在均匀岩石物理模型内区域的应力值和位移值；

步骤S40、根据均匀岩石物理模型区域的边界和内区域的全部应力值和位移值确定地震波的波形记录；

其中，所述步骤S10具体包括以下步骤：

步骤S101、将热弹性波动方程写成频率域形式并转化为二维的情形如下：公式(1)中，

其中，Δ代表拉普拉斯算子，i为虚数单位，λ和μ代表拉梅常数，u代表质点的位移，ux代表质点在x方向的位移，uz代表质点在z方向的位移，γ代表应力的热力系数，θ代表介质的相对温度，ρ代表介质的密度，ω代表角频率，Xx和Xz分别代表x方向和z方向的体力，κ代表热扩散系数，Q代表内热源，k代表导热系数，T0代表介质的初始温度，T代表介质的实时温度，αt代表热膨胀系数，η为中间变量，c3为单位初始体积常应变比热；

步骤S102、将上述热弹性波动方程转化为如下边界元形式方程：公式(2)中，Σ*和V*是热弹性波动方程的应力基本解和位移的基本解，u和p是待求的位移值和应力值，i和j为第i个或者第j个离散单元，N为离散的总单元数，Γ为离散单元的边界；

其中，所述步骤S20进一步包括：

步骤S201、给出均匀岩石物理模型参数、观测系统中各检波点及各炮点的位置坐标、边界元法网格单元的各个控制点坐标及区域的划分、震源子波的信息以及均匀岩石物理模型区域边界；

步骤S202、将步骤S201中给出的所有信息代入边界元形式方程中，求得地震波在均匀岩石物理模型区域的边界的应力值和位移值。

2.根据权利要求1所述的高温介质地震波传播模拟方法，其特征在于，在步骤S202之前，还包括以下步骤：利用傅里叶变换将时间域的热弹性波动方程转化为频率域的热弹性波动方程，以便在频率域进行边界元形式方程的计算。

3.根据权利要求2所述的高温介质地震波传播模拟方法，其特征在于，所述步骤S202进一步包括：步骤S2021、利用形函数对所述均匀岩石物理模型的边界进行离散，计算离散点坐标，同时离散外边界吸收层数；

步骤S2022、利用贝塞尔函数计算所述边界元形式方程的积分核函数，得出第一边界矩阵系数；

步骤S2023、利用ASSEM函数对内边界和外边界进行处理，即将所述第一边界矩阵系数进行组装；

步骤S2024、采用无限元方法进行人工截断边界处理得出第二边界矩阵系数，并将所述第二边界矩阵系数组合到所述第一边界矩阵系数中；

步骤S2025、利用块状高斯消元法对代入步骤S2024中组合后的第一边界矩阵系数的边界元形式方程进行求解，得出地震波在均匀岩石物理模型区域的边界的应力值和位移值。

4.根据权利要求3所述的高温介质地震波传播模拟方法，其特征在于，所述步骤S30进一步包括：步骤S301、将地震波在均匀岩石物理模型区域的边界的应力值和位移值代入边界元形式方程；

步骤S302、利用贝塞尔函数计算上述边界元形式方程的积分核函数，得出第三边界矩阵系数；

步骤S303、利用块状高斯消元法对代入第三边界矩阵系数的边界元形式方程进行求解，得出地震波在均匀岩石物理模型的内区域的应力值和位移值。

5.根据权利要求4所述的高温介质地震波传播模拟方法，其特征在于，所述步骤S40进一步包括：利用傅里叶变换将地震波在均匀岩石物理模型的边界的应力值和位移值以及内区域的应力值和位移值从频率域转化为时间域；

根据时间域下的地震波在均匀岩石物理模型的边界的应力值和位移值以及内区域的应力值和位移值获得地震波的波形记录。

6.根据权利要求5所述的高温介质地震波传播模拟方法，其特征在于，所述均匀岩石物理模型参数包括介质的密度ρ、拉梅常数λ和μ、热扩散系数κ、热膨胀系数αt、单位初始体积常应变比热cε。

7.根据权利要求6所述的高温介质地震波传播模拟方法，其特征在于，所述震源子波为高斯子波或者雷克子波，所述震源子波的信息包括震源子波的起始频率、中心频率和终止频率。

说明书

一种高温介质地震波传播模拟方法

技术领域

[0001] 本发明属于地球物理勘探研究领域，具体提供一种高温介质地震波传播模拟方法。

背景技术

[0002] 地震波波场正演模拟越来越受到地球物理学者们的重视,它的地位也显得越来越重要，正演模拟技术的日趋成熟也使得反演技术在地震勘探中应用的准确性提高。而传统的正演模型往往是常温介质，并没有考虑影响弹性介质物理属性的温度因素，特别是高温状态下介质的弹性属性变化。

[0003] 众所周知，在不同的温度下，其岩石密度、弹性模量、波速会有较大的差异，尤其对于有裂缝存在的情况下。随着温度的增加，原始的缝洞开裂密度逐渐变大，其必然会引起储层岩性物性参数的变化。所以在实际的深层—超深层地震勘探中，必须考虑到温度的影响，才能更加准确的地反映真实地下构造，故高温介质中地震波传播模拟方法至关重要。

[0004] 相应地，本领域需要一种高温介质地震波传播模拟方法来解决上述问题。

发明内容

[0005] 为了解决现有技术中的上述问题，即为了解决传统的地震波正演模型为常温介质导致的模拟结果不准确的问题，本发明提供了一种高温介质地震波传播模拟方法，所述模拟方法包括以下步骤：

[0006] 步骤S10、利用边界元法将热弹性波动方程转化为边界元形式方程；

[0007] 步骤S20、利用边界元形式方程求得地震波在均匀岩石物理模型区域的边界的应力值和位移值；

[0008] 步骤S30、将所述边界的应力值和位移值代入边界元形式方程求得地震波在均匀岩石物理模型内区域的应力值和位移值；

[0009] 步骤S40、根据均匀岩石物理模型区域的边界和内区域的全部应力值和位移值确定地震波的波形记录。

[0010] 在上述高温介质地震波传播模拟方法的优选技术方案中，所述步骤S10具体包括以下步骤：

[0011] 步骤S101、将热弹性波动方程写成频率域形式并转化为二维的情形如下：

[0012]

[0013] 公式(1)中，

[0014] θ＝T-T0,γ＝(3λ+2μ)αt，

[0015] 其中，Δ代表拉普拉斯算子，i为虚数单位，λ和μ代表拉梅常数，u代表质点的位移，ux代表质点在x方向的位移，uz代表质点在z方向的位移，γ代表应力的热力系数，θ代表介质的相对温度，ρ代表介质的密度，ω代表角频率，Xx和Xz分别代表x方向和z方向的体力，i代表第i个离散单元，κ代表热扩散系数，Q代表内热源，k代表导热系数，T0代表介质的初始温度，T代表介质的实时温度，αt代表线膨胀系数，η为中间变量。

[0016] 步骤S102、将上述热弹性波动方程转化为如下边界元形式方程：

[0017]

[0018] 公式(2)中，u和p是公式(2)待求的位移值和应力值，i和j为第i个或者第j个离散单元，N为离散的总单元数，Σ*和V*是热弹性波动方程的应力基本解和位移的基本解，Γ为离散单元的边界。

[0019] 在上述高温介质地震波传播模拟方法的优选技术方案中，所述步骤S20进一步包括：

[0020] 步骤S201、给出均匀岩石物理模型参数、观测系统中各检波点及各炮点的位置坐标、边界元法网格单元的各个控制点坐标及区域的划分、震源子波的信息以及均匀岩石物理模型区域边界；

[0021] 步骤S202、将上述所有信息代入边界元形式方程中，求得地震波在均匀岩石物理模型区域的边界的应力值和位移值。

[0022] 在上述高温介质地震波传播模拟方法的优选技术方案中，在步骤S202之前，所述模拟方法还包括以下步骤：

[0023] 利用傅里叶变换将所述的时间域的热弹性波动方程转化为频率域的热弹性波动方程，以便在频率域进行边界元形式方程的计算。

[0024] 在上述高温介质地震波传播模拟方法的优选技术方案中，所述步骤S202进一步包括：

[0025] 步骤S2021、利用形函数对所述均匀岩石物理模型的边界进行离散，计算离散点坐标，同时离散外边界吸收层数；

[0026] 步骤S2022、利用贝塞尔函数计算所述边界元形式方程的积分核函数，得出第一边界矩阵系数；

[0027] 步骤S2023、利用ASSEM函数对内边界和外边界处理，即将所述第一边界矩阵系数进行组装；

[0028] 步骤S2024、采用无限元方法进行人工截断边界处理得出第二边界矩阵系数，并将所述第二边界矩阵系数组合到所述第一边界矩阵系数中；

[0029] 步骤S2025、利用块状高斯消元法对代入步骤S2024中组合后的第一边界矩阵系数的边界元形式方式进行求解，得出地震波在均匀岩石物理模型区域的边界的应力值和位移值。

[0030] 在上述高温介质地震波传播模拟方法的优选技术方案中，所述步骤S30进一步包括：

[0031] 步骤S301、将地震波在均匀岩石物理模型区域的边界的应力值和位移值代入边界元形式方程；

[0032] 步骤S302、利用贝塞尔函数计算上述边界元形式方程的积分核函数，得出第三边界矩阵系数；

[0033] 步骤S303、利用块状高斯消元法对代入第三边界矩阵系数的边界元形式方式进行求解，得出地震波在均匀岩石物理模型的内区域的应力值和位移值。

[0034] 在上述高温介质地震波传播模拟方法的优选技术方案中，所述步骤S40进一步包括：利用傅里叶变换将地震波在均匀岩石物理模型的边界的应力值和位移值以及内区域的应力值和位移值从频率域转化为时间域；根据时间域下的地震波在均匀岩石物理模型的边界的应力值和位移值以及内区域的应力值和位移值获得地震波的波形记录。

[0035] 在上述高温介质地震波传播模拟方法的优选技术方案中，所述均匀岩石物理模型参数包括密度ρ、拉梅常数λ和μ、热扩散系数κ、热膨胀系数αt、单位初始体积常应变比热cε。

[0036] 在上述高温介质地震波传播模拟方法的优选技术方案中，所述震源子波为高斯子波或者雷克子波，所述震源子波信息包括震源子波的起始频率、中心频率和终止频率。

[0037] 本领域技术人员能够理解的是，在本发明的技术方案中，高温介质地震波传播模拟方法包括以下步骤：步骤S10、利用边界元法将热弹性波动方程转化为边界元形式方程；步骤S20、利用边界元形式方程求得地震波在均匀岩石物理模型区域的边界的应力值和位移值；步骤S30、将所述边界的应力值和位移值代入边界元形式方程求得地震波在均匀岩石物理模型内区域的应力值和位移值；步骤S40、根据均匀岩石物理模型区域的边界和内区域的全部应力值和位移值确定地震波的波形记录。

[0038] 通过这样的设置，有如下有益效果：1、由于边界元法的理论基础是边界积分方程理论，其特点是降维、解析与离散相结合、能自动满足无穷远处边界条件，因此该方法只需在区域的边界上划分单元，不必在无穷远边界上划分单元，且计算精度高，适于无限域问题。2、引入温度参量的热弹性波动方程，计及到温度的影响，使结果更加贴切真实情况。本发明实现了基于热弹性波动方程的高温介质地震波传播模拟方法，开展热弹理论研究，可以准确考察高温对地震波传播的影响，对深层及超深层油气勘探领域拓展提供理论基础，对解释地球深部观测数据异常现象提供合理的解释。

附图说明

[0039] 下面参照附图来描述本发明的优选实施方式，附图中：

[0040] 图1是本发明一种实施例的高温介质地震波传播模拟方法的主要步骤流程图；

[0041] 图2是本发明一种实施例的高温介质地震波传播模拟方法的具体步骤流程图；

[0042] 图3是本发明一种实施例的地震波水平方向的波形图；

[0043] 图4是本发明一种实施例的地震波垂直方向的波形图。

具体实施方式

[0044] 为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

[0045] 下面参照附图来描述本发明的优选实施方式。本领域技术人员应当理解的是，这些实施方式仅仅用于解释本发明的技术原理，并非旨在限制本发明的保护范围。

[0046] 地震波波场正演模拟方法的种类很多，一般将其分为波动方程法、积分方程法、几何射线法三大类。我们基于耦合热弹性理论建立热弹性波动方程，利用边界元数值方法进行高温下地震波传播正演模拟。热弹性耦合理论由热力学基本定律、介质的本构理论和Helmholtz自由能表达式导出的热弹性材料的热传导方程中，除了待定的温度场函数外，还含有应变率。它表明物体上的温度场不仅取决于热源，以及各有关的热力学物性系数和换热边界条件，而且还受到弹性变形应变率的影响，或者说，弹性变形的体积应变率将在一定的程度上改变物体上的热量的传递，因此热传导方程和热弹性运动方程必须耦合求解。

[0047] 边界元法是对边界积分方程离散求解的现代数值分析方法，属于近似解法之一。边界积分方程从定解问题的控制方程转换而得，其最大的特点就是降低了求解问题的维数，将二维问题化为其边界线上的一维问题，它只以边界变量为基本变量，域内位置可以在需要时根据边界变量求出。基于热弹性理论及边界元法，我们发展一种热介质中弹性波边界元法正演模拟。

[0048] 如图1所示，图1是本发明一种实施例的高温介质地震波传播模拟方法的主要步骤流程图。参照图1，本发明的高温介质地震波传播模拟方法主要包括以下步骤：

[0049] 步骤S10、利用边界元法将热弹性波动方程转化为边界元形式方程；

[0050] 优选地，步骤S10进一步包括：步骤S101、将热弹性波动方程写成频率域形式并转化为二维的情形如下：

[0051]

[0052] 公式(1)中，

[0053] θ＝T-T0,γ＝(3λ+2μ)αt,

[0054] 其中，Δ代表拉普拉斯算子，i为虚数单位，λ和μ代表拉梅常数，u代表质点的位移，ux代表质点在x方向的位移，uz代表质点在z方向的位移，γ代表应力的热力系数，θ代表介质的相对温度，ρ代表介质的密度，ω代表角频率，Xx和Xz分别代表x方向和z方向的体力，i代表第i个离散单元，κ代表热扩散系数，Q代表内热源，k代表导热系数，T0代表介质的初始温度，T代表介质的实时温度，αt代表线膨胀系数，η为中间变量。

[0055] 如下图所示，在本发明的实施例中常数和一些参数的取值如下：

[0056] ρ(kg·m-3) λ(N·m-2) μ(N·m-2) κ(m2/s) αt(1/K) cε(J/m3/K)2.65 4×109 6×109 0.089 0.33×10-5 1.17×102

[0057] 步骤S102、将上述热弹性波动方程转化为如下边界元形式方程：

[0058]

[0059] 公式(2)中，Σ*和V*是热弹性波动方程的应力基本解和位移的基本解，u和p是公式(2)待求的位移值和应力值，i和j为第i个或者第j个离散单元，N为离散的总单元数，Γ为离散单元的边界。

[0060] 一方面，通过引入温度参量的热弹性波动方程，计及到温度的影响，使模拟结果更加贴切实际情况。另一方面，由于边界元法的理论基础是边界积分方程理论，其特点是降维、解析与离散相结合、能自动满足无穷远处边界条件，因此该方法只需在区域的边界上划分单元，不必在无穷远边界上划分单元，且计算精度高，适于无限域问题。

[0061] 步骤S20、利用边界元形式方程求得地震波在均匀岩石物理模型区域的边界的应力值和位移值。

[0062] 在本发明的一种实施例中，步骤S20进一步包括：步骤S201、给出均匀岩石物理模型参数、观测系统中各检波点及各炮点的位置坐标、边界元法网格单元的各个控制点坐标及区域的划分、震源子波的信息以及均匀岩石物理模型区域边界。

[0063] 其中，均匀岩石物理模型参数具体包括：密度ρ、拉梅常数λ和μ、热扩散系数κ、热膨胀系数αt、单位初始体积常应变比热cε共六个参数。这六个参数的信息需要从靶区的原位岩石获取，对靶区的原位岩石进行高温原位热岩石物理实验测量，使求得的参数更加准确、真实可靠。其中，震源子波的信息包括起止频率、中心频率、终止频率等信息。

[0064] 步骤S202、将上述所有信息代入边界元形式方程中，求得地震波在均匀岩石物理模型区域的边界的应力值和位移值。

[0065] 需要说明的是，由于在时间域进行边界元形式方程的计算比较复杂，因此为了方便计算，在步骤S202之前，所述方法还包括：利用傅里叶变换将时间域的热弹性波动方程转化为频率域的热弹性波动方程，即公式(1)所示，以便在频率域进行边界元形式方程的计算。

[0066] 在本发明的一种实施例中，步骤S202进一步包括：

[0067] 步骤S2021、利用形函数对均匀岩石物理模型的边界进行离散，计算离散点坐标，同时离散外边界吸收层数；

[0068] 如步骤S201已经给出边界元法网格单元的各个控制点坐标及区域的划分，利用形函数及区域边界的控制点，对区域边界中进行插值，从而计算离散点坐标及离散外边界吸收层数；

[0069] 步骤S2022、利用贝塞尔函数计算所述边界元形式方程的积分核函数，得出第一边界矩阵系数；

[0070] 步骤S2023、利用ASSEM函数对内边界和外边界处理，即将所述第一边界矩阵系数进行组装；

[0071] 步骤S2024、采用无限元方法进行人工截断边界处理得出第二边界矩阵系数，并将所述第二边界矩阵系数组合到所述第一边界矩阵系数中；即在得到第一边界矩阵系数的基础上，再增加扩充矩阵边界系数。

[0072] 步骤S2025、利用块状高斯消元法对代入步骤S2024中组合后的第一边界矩阵系数的边界元形式方式进行求解，得出地震波在均匀岩石物理模型区域的边界的应力值和位移值。

[0073] 通过求出地震波在均匀岩石物理模型区域的边界的应力值和位移值，来进行地震波在均匀岩石物理模型内区域的应力值和位移值的求解。

[0074] 步骤S30、将边界的应力值和位移值代入边界元形式方程求得地震波在均匀岩石物理模型内区域的应力值和位移值；

[0075] 在本发明的一种实施例中，步骤S30进一步包括：

[0076] 步骤S301、将地震波在均匀岩石物理模型区域的边界的应力值和位移值代入边界元形式方程；

[0077] 步骤S302、利用贝塞尔函数计算上述边界元形式方程的积分核函数，得出第三边界矩阵系数；

[0078] 步骤S303、利用块状高斯消元法对代入第三边界矩阵系数的边界元形式方式进行求解，得出地震波在均匀岩石物理模型的内区域的应力值和位移值。

[0079] 可以理解的是，要得出多组地震波在均匀岩石物理模型的边界以及内区域的应力值和位移值，需要输入震源子波的某个频率类型的不同频率，得出地震波在震源子波的不同频率下的在均匀岩石物理模型的边界以及内区域的应力值和位移值。如改变给定的震源子波的起始频率：30赫兹、40赫兹、50赫兹等等，计算出地震波在相对应的赫兹下的在均匀岩石物理模型的边界以及内区域的应力值和位移值。

[0080] 步骤S40、根据均匀岩石物理模型区域的边界和内区域的全部应力值和位移值确定地震波的波形记录。

[0081] 可以理解的是，由于计算出来的应力值和位移值是在频率下，因此步骤S40进一步包括：利用傅里叶变换将地震波在均匀岩石物理模型的边界的应力值和位移值以及内区域的应力值和位移值从频率域转化为时间域；根据时间域下的地震波在均匀岩石物理模型的边界的应力值和位移值以及内区域的应力值和位移值获得地震波的波形记录。

[0082] 需要说明的是，上述实施例是在模拟的岩石物理模型为均匀介质下进行的，在岩石物理模型为多层不均匀时，需要引入ANGLE函数对每层地程的倾斜角度进行计算，然后再进行步骤S20、S30和S40。

[0083] 在一种可能的实施例中，如图2所示，图2是本发明一种实施例的高温介质地震波传播模拟方法的具体步骤流程图。参照图2，高温介质地震波传播模拟方法包括：

[0084] 1)给定输入的均匀岩石物理模型，其模型包括：密度ρ、拉梅常数λ和μ、热扩散系数κ、热膨胀系数αt、单位初始体积常应变比热cε共六个参数的信息；2)给出观测系统中各检波点及各炮点的位置坐标；3)给出边界元法网格单元的各个控制点坐标及区域的划分；4)给出震源子波的信息，包括起止频率、中心频率、终止频率等信息；5)根据前面4个步骤得到的所有输入信息，输入子模块EDITP中，同时子模块EDITP包含模拟的均匀岩石物理模型区域边界以及模拟震源子波的类型(程序中可选定高斯子波或者雷克子波)；我们选用形函数对均匀岩石物理模型的边界进行离散处理，从而计算边界上的离散点坐标，同时离散外界吸收边界层数；利用傅里叶变换(FFT)对时间域的热弹性波动方程转化为频率域的热弹性波动方程，从而在频率域进行边界元形式方程的离散；调用ASSEM函数对第一边界矩阵系数进行组装，其中第一边界矩阵系数为通过调用贝塞尔函数计算边界元形式方程的积分核函数得出；本程序中采用无限元方法进行人工截断边界处理对第一边界矩阵系数进行重新组装；然后通过块状高斯消元求解，得出地震波在均匀岩石物理模型区域的边界的应力值和位移值；最后将地震波在均匀岩石物理模型区域的边界的应力值和位移值代入边界元形式方程中，通过贝塞尔函数和块状高斯消元求解检波点波场(即地震波在均匀岩石物理模型的区域内的应力值和位移值)；6)将频率域的结果转换为时间域，得出全部节点的位移和应力(即地震波在均匀岩石物理模型的区域内和边界的应力值和位移值)，为后续的反演做准备，从而进行后续的分析。

[0085] 如图3和图4所示，图3是本发明一种实施例的地震波水平方向的波形图；图4是本发明一种实施例的地震波垂直方向的波形图。参照图3和图4，从特定的输入参数中，我们得到地震波水平和垂直方向的记录，图中我们可以看出，高温介质下的地震波传播记录中一共存在三种波形，即第一纵波、第二纵波(热波)及横波。在不考虑温度效应中，即传统的均匀各项同性介质中的仅仅存在两种记录即横波和纵波记录，我们利用弹性理论模拟地震波在热介质中的传播有新的波形出现，为地震勘探提供新思路。

[0086] 通过上述描述可以看出，本发明的高温介质地震波传播模拟方法具有以下优点：1、输入的模型岩石物理参数准确，其中密度ρ、拉梅常数λ和μ、热扩散系数κ、热膨胀系数αt、单位初始体积常应变比热cε共六个参数的信息共六个参数的信息需要靶区的原位岩石，从而进行热岩石物理实验测量，求得的参数更加真实可靠。2、利用边界元法，边界元法的理论基础是边界积分方程理论，其特点是降维、解析与离散相结合、能自动满足无穷远处边界条件。该方法只需在区域的边界上划分单元，不必在无穷远边界上划分单元，计算精度高，适于无限域问题。3、引入温度参量的热弹性波动方程，计及到温度的影响，使结果更加贴切真实情况。本发明实现了基于热弹性波动方程的高温介质地震波传播模拟方法，开展热弹理论研究，可以准确考察高温对地震波传播的影响，对深层及超深层油气勘探领域拓展提供理论基础，对解释地球深部观测数据异常现象提供合理的解释。

[0087] 至此，已经结合附图所示的优选实施方式描述了本发明的技术方案，但是，本领域技术人员容易理解的是，本发明的保护范围显然不局限于这些具体实施方式。在不偏离本发明的原理的前提下，本领域技术人员可以对相关技术特征作出等同的更改或替换，这些更改或替换之后的技术方案都将落入本发明的保护范围之内。