一种重力梯度的自适应剖分反演方法与流程

1.本发明属于地球物理学技术领域，尤其涉及一种重力梯度的自适应剖分反演方法。


背景技术：

2.近年来随着现代科技水平、地球物理仪器和测试技术的发展，重力梯度测量方法获得了广泛应用。相比于传统重力异常数据，重力梯度数据对浅层异常和微小密度差引起的异常更敏感，精度较高，包含更多的地质信息，在对地质体的边界、形状识别上有较好的应用效果。
3.我国矿产开采深度在500m左右，当前矿产资源形势短缺，为实现矿产资源的可持续性发展，满足社会经济发展的需求，我们应提升反演的效率及精度，加快找矿速度。
4.地下包含着大量的地质信息，为了让反演信息更加精确，可以将剖分网格缩小，但同时反演效率降低，增大cpu内存。


技术实现要素：

5.本发明提出一种重力梯度的自适应剖分反演方法，在边界识别范围内进行网格剖分，使其能够减少剖分网格数量，将剖分区域划分在矿体边界内，该种方法既没有损失精度，又有效的提升了反演效率，减少cpu内存；且考虑到实际地质情况，本发明采用四面体剖分来拟合起伏地形和不规则地质体。
6.具体的，本发明是通过以下技术方案实现的：
7.提供一种重力梯度的自适应剖分反演方法，包括以下步骤：
8.步骤1：获得重力梯度异常数据或正演出模型的重力梯度异常数据；
9.步骤2：采用倾斜角边界识别法，在结果大于等于0的范围内进行四面体剖分；
10.步骤3：构建非全局核函数矩阵，并实现非全局四面体剖分反演计算，得到空间密度分布图。
11.作为本发明的进一步说明，所述倾斜角边界识别法，公式如下：
[0012][0013]
式中：
[0014][0015]
其中g
z
为重力梯度异常。
[0016]
作为本发明的进一步说明，所述构建非全局核函数矩阵的过程具体为：
[0017]
将任意多面体重力异常表达式应用于所述四面体，任意观测点(x
n
,y
n
,z
n
)处的重力异常计算公式如下：
[0018][0019]
重力梯度为重力异常的一阶垂直导：
[0020]
其中：g为核函数矩阵；γ为万有引力常数6.672
×
10
‑
11
m3/(kg
·
s2)；ρ为四面体密度；φ
i
为四面体第i个三角形与z轴的夹角；j
j,i
由如下公式得到：
[0021][0022]
其中：i为三角形序号，j为第i个三角形中的边序号；
[0023][0024][0025][0026]
θ
i
为四面体第i个面外法向在水平方向的投影与x方向的夹角，式中0≤θ
i
≤2π，0≤φ
i
≤π。
[0027]
作为本发明的进一步说明，在实现非全局四面体剖分反演计算时，所求解的方程解为欠定方程的解，在求解欠定方程时，构造下式目标函数：
[0028][0029]
m是密度参数的向量；m
ref
是参考模型，一般取值为0；λ为正则化参数，起到平衡数据拟合函数与权函数的作用；w
d
是数据加权矩阵；为深度权，用于约束反演结果使结果更加符合实际克服趋肤效应，且我们在四面体剖分反演中加上与体积有关的权：v代表四面体的体积，z0是由深度权函数与观测点下的核函数进行匹配得到，α是一个经验值，通常取0.5，β是一个经验值。
[0030]
与现有技术相比，本发明具有以下有益的技术效果：
[0031]
本发明通过边界识别圈定出目标矿体的范围，在该范围内进行网格剖分；并采用非全局四面体剖分来拟合起伏地形和不规则地质体，更加贴合地形形状；相比传统的物性反演方法，本发明能够实现在同样的区域内通过缩小剖分范围来减少剖分网格，以达到加快反演速度，提高反演精度的目的，同时降低对计算机性能的要求。
附图说明
[0032]
图1是本发明提供的重力梯度的自适应剖分反演方法的流程图；
[0033]
图2是三维四面体剖分示意图；
[0034]
图3是传统全局四面体剖分反演结果图；
[0035]
图4是本发明提供的非全局四面体剖分反演结果图。
具体实施方式
[0036]
为了能够更清楚地理解本发明的上述目的、特征和优点，下面结合附图和具体实施例对本发明进行详细描述。需要说明的是，在不冲突的情况下，本发明的实施例及实施例中的特征可以相互组合。
[0037]
除非另有定义，本文所使用的所有的技术和科学术语与属于本发明的技术领域的技术人员通常理解的含义相同。本文中在本发明的说明书中所使用的术语只是为了描述具体的实施例的目的，不是旨在于限制本发明。
[0038]
图1为重力梯度的自适应剖分反演方法的流程图，该方法包括如下步骤：
[0039]
步骤1：获得重力梯度异常数据或正演出模型的重力梯度异常数据；
[0040]
实际测量过程中，我们测量所得到的数据比地下地质体范围更大，因此在反演时，我们将地质体剖分范围用边界识别来界定，在边界识别结果大于等于0的范围内进行网格剖分，可有效提升反演效率和反演精度。
[0041]
步骤2：采用倾斜角边界识别法，在结果大于等于0的范围内进行四面体剖分：
[0042]
边界识别方法采用倾斜角法，公式如下：
[0043][0044]
式中：
[0045][0046]
其中g
z
为重力梯度异常。
[0047]
步骤3：构建非全局核函数矩阵，并实现非全局四面体剖分反演计算，得到空间密度分布图：
[0048]
针对实际地质体的非规则性，我们采用四面体剖分来更好的拟合；okabe(1979)推导出任意多面体重力异常表达式，将其应用于四面体，任意观测点(x
n
,y
n
,z
n
)处的重力异常计算公式如下：
[0049][0050]
重力梯度为重力异常的一阶垂直导：
[0051]
其中：g为核函数矩阵；γ为万有引力常数6.672
×
10
‑
11
m3/(kg
·
s2)；ρ为四面体密度；φ
i
为四面体第i个三角形与z轴的夹角；j
j,i
由如下公式得到：
[0052][0053]
其中：i为三角形序号，j为第i个三角形中的边序号。
[0054][0055][0056][0057]
θ
i
为四面体第i个面外法向在水平方向的投影与x方向的夹角，注意：式中0≤θ
i
≤2π，0≤φ
i
≤π；
[0058]
反演求解的是地下空间密度分布，由于观测点数远小于地下剖分块体数，因此方程解为欠定方程的解。为求解欠定方程，构造下式目标函数：
[0059][0060]
m是密度参数的向量；m
ref
是参考模型，一般取值为0；λ为正则化参数，起到平衡数据拟合函数与权函数的作用；w
d
是数据加权矩阵；为深度权，用于约束反演结果使结果更加符合实际克服趋肤效应，因每个四面体体积大小不一，而体积的大小会影响深度权的作用，因此在四面体剖分反演中加上与体积有关的权：v代表四面体的体积，z0是由深度权函数与观测点下的核函数进行匹配得到，α是一个经验值，通常取0.5，β是一个经验值。
[0061]
下表为本发明与传统物性反演方法，在采用同等反演精度下(相同采样点距、相同成图精度、相同的加权函数)，反演相同数据时采用同一的计算机的同一软件反演时的运行时间以及占用内存对比表。
[0062] 运行时间占用内存本发明270.56s62978kb传统方法747.45s165539kb
[0063]
本发明通过边界识别圈定出目标矿体的范围，在该范围内进行网格剖分；并采用非全局四面体剖分来拟合起伏地形和不规则地质体，更加贴合地形形状；相比传统的物性反演方法，本发明能够实现在同样的区域内通过缩小剖分范围来减少剖分网格，以达到加快反演速度，提高反演精度的目的，同时降低对计算机性能的要求。
[0064]
最后应说明的是，以上实施例仅用以说明本发明的技术方案而非限制，尽管参照
较佳实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，可以对本发明的技术方案进行修改或等同替换，而不脱离本发明技术方案的精神和范围。