CN113267830A - 基于非结构网格下二维重力梯度与地震数据联合反演方法 - Google Patents

Abstract

本发明公开了一种基于非结构网格下二维重力梯度与地震数据联合反演方法，包括如下过程：获得重力梯度异常数据；将地下半空间剖分成个若干三角形单元；实现非规则网格下二维重力梯度及地震声波数据的正演计算；然后分别获得不规则网格下二维重力梯度、地震声波数据单独反演结果；将交叉梯度约束函数引入非规则网格中，实现非规则网格下的二维重力梯度与地震声波数据的联合反演；进行多次迭代计算，检验反演拟合程度，直至得到高精度的反演结果。本发明对二维离散空间实现不规则网格剖分，更加拟合非规则目标地质体，将交叉梯度约束函数引入非结构网格中，实现非规则网格下的重力梯度与地震数据联合反演。

Description

基于非结构网格下二维重力梯度与地震数据联合反演方法

技术领域

本发明属于地球物理学技术领域，尤其涉及一种基于非结构网格下二维重力梯度与地震数据联合反演方法。

背景技术

重力梯度勘探以勘探目标与围岩间密度差异为物性基础，通过测量重力梯度异常来研究地下空间地质构造特征。因其具有经济、勘探深度大及快速获得面积上信息等优点，在地球深部构造探测、区域地质构造单元划分、沉积盆地圈定、固体矿产及油气资源勘查等领域获得广泛应用。地震信号记录地震波穿过地下不同介质的响应，对介质的波速参数敏感，且在射线覆盖密集的区域具有较高的垂向分辨率。近年来，随着硬件设备的升级和计算机技术的快速发展，勘探装备精度与效率都有了很大提高，同时数据处理与解释方法也由传统的定性解释逐步向定量解释发展。目前，我国矿产资源接替基地面临的主要找矿难题是：老矿山深部和各类隐伏区的探矿难度大，急需先进、高效的理论和技术方法指导深部找矿。

随着卫星、航空等勘探手段的综合应用，以及数据精度和勘探精度要求的不断提高，开展高精度联合反演方法用于深部找矿是十分有必要的。

地球物理反演是探索地下结构的最佳途径之一。单一地球物理反演方法具有多解性和不确定性，因此联合多种地球物理数据进行联合反演可以实现相互约束、相互补充，使模型和数据拟合度达到最佳。密度反演是根据观测的异常数据对地下地质体密度分布的定量计算，进而可以用来推测地质体的空间分布体积。但是随着深度的增加，重力梯度异常对地下目标体的敏感度衰减，缺少垂向分辨率；地震信号记录地震波穿过地下不同介质的响应，对介质的波速参数敏感，且在射线覆盖密集的区域具有较高的垂向分辨率。因此将重力梯度和地震数据联合，可以发挥各数据在分辨率和敏感度上的优势。相比于单独反演，重力梯度与地震联合反演在减少反演多解性、提高模型质量和精细综合地球物理解释等方面都发挥出重要作用。

现今重力梯度与地震数据联合反演的离散方式大多是采用结构化网格，非结构网格下重力梯度与地震数据联合反演采用的方式依据速度和密度转换关系式，这种方式更多依赖先验转换关系式。由于非规则网格可以较高精度地拟合任何结构，对于非规则目标具有很好的适用性。因此开展一种非结构网格下的重力梯度与地震数据联合反演方法是有必要的。

发明内容

本发明将交叉梯度约束函数引入非规则网格下的二维重力梯度与地震联合反演，对于规则模型、非规则模型均具有很好的适用性，可以实现高精度的联合反演。同时，引入不同物性参数间的耦合方式—基于交叉梯度函数约束，不再依赖于先验速度和密度转换关系式。

具体的，本发明是通过以下技术方案实现的：

提供一种基于非结构网格下二维重力梯度与地震数据联合反演方法，包括以下步骤：

步骤1：获得重力梯度异常数据；

步骤2：将地下半空间剖分成个若干三角形单元；

步骤3：实现非规则网格下二维重力梯度及地震声波数据的正演计算，然后分别获得不规则网格下二维重力梯度、地震声波数据单独反演结果；

步骤4：将交叉梯度约束函数引入非规则网格中，实现非规则网格下的二维重力梯度与地震声波数据的联合反演；

步骤5：进行多次迭代计算，检验反演拟合程度，直至得到高精度的反演结果。

作为本发明的进一步说明，所述实现非规则网格下二维重力梯度及地震声波数据的正演计算，然后分别获得不规则网格下二维重力梯度、地震声波数据单独反演结果具体为：

所述三角形的重力异常g可以转换为各边贡献的和，表示为

G是牛顿的引力常数,ρ代表每一个单元的密度值。x_i、z_i为三角形单元格中每个点的x坐标和z坐标。根据格林定理，可以表示为：

逆时针旋转坐标轴(XOZ)直到旋转坐标系(XOZ)Z轴平行于该边L_ij的外法向量为止；故关于一条边L_ij的重力异常g_ij、重力梯度异常g_zij表示为

式中G为牛顿引力常数，ρ为每一个单元的密度，θ代表逆时针旋转的角度；

整个三角形的重力梯度异常计算如下:

g_z＝g_zij+g_zjk+g_zki

地震声波方程如下：

p(x,t)代表波场值；ρ(x)和v(x)是地下密度和速度分布；s(x_s,t)代表震源函数；

格子法的核心是积分平衡的微分方程弱形式，在k点邻域对所述地震声波方程作面积分并应用所述格林定理可得到：

式中m是围绕节点k的三角单元数，S_kl是节点k周围第l个三角形单元中的虚线段，a和β是虚线段包络线的外法线方向余弦，

是震源函数的面积分；

利用动力学计算中的集中质量模型及三角单元线性插值方法，可得到上式的空间离散形式：

其中

式中下标r指代三角单元中的节点，A是三角单元面积，b_k和a_k是三角单元的几何参数，用二阶中心差分来离散公式左侧项中的时间导数，即可实现波场值的迭代更新。

作为本发明的进一步说明，所述将交叉梯度约束函数引入非规则网格中，实现非规则网格下的二维重力梯度与地震声波数据的联合反演具体为：

其中：

式中下标r指代三角单元中的节点，A是三角单元面积，b_k和a_k是三角单元的几何参数；

m⁽¹⁾和m⁽²⁾代表密度和速度参数d_obs ⁽¹⁾和d_obs ⁽²⁾分别代表观测重力梯度数据和地震数据；α₁和α₂代表数据权，用来平衡重力梯度数据和地震波形数据的拟合差；τ₁和τ₂代表正则化系数；

和

分别代表密度和速度部分的模型正则化算子；Φ为密度模型与速度模型的交叉梯度约束函数。

与现有技术相比，本发明具有以下有益的技术效果：

本发明对二维离散空间实现不规则网格剖分，非结构网格剖分更加精细化，更加拟合非规则目标地质体，并将交叉梯度约束函数引入非结构网格中，实现非规则网格下的重力梯度与地震数据联合反演。

附图说明

图1是本发明提供的基于非结构网格下二维重力梯度与地震数据联合反演方法的流程图；

图2是地下半空间网格剖分示意图；

图3是坐标系旋转示意图；

图4是局部非规则网格剖分示意图；

图5是联合反演密度结果图；

图6是联合反演速度结果图。

具体实施方式

为了能够更清楚地理解本发明的上述目的、特征和优点，下面结合附图和具体实施例对本发明进行详细描述。需要说明的是，在不冲突的情况下，本发明的实施例及实施例中的特征可以相互组合。

除非另有定义，本文所使用的所有的技术和科学术语与属于本发明的技术领域的技术人员通常理解的含义相同。本文中在本发明的说明书中所使用的术语只是为了描述具体的实施例的目的，不是旨在于限制本发明。

图1为基于非结构网格下二维重力梯度与地震数据联合反演方法的流程图，该方法包括如下步骤：

步骤1：获得重力梯度异常数据(Δg_z(m×1))；

步骤2：将地下半空间剖分成个若干三角形单元(如图2所示)；

步骤3：实现非规则网格下二维重力梯度及地震声波数据的正演计算，然后分别获得不规则网格下二维重力梯度、地震声波数据单独反演结果；

对于二维情况，通常对异常源进行假设，认为其在沿某个坐标轴方向是没有变化的，因此通常对其界面进行分析便可以知道二度体的重力异常分布情况。

三角形的重力异常g可以转换为各边贡献的和，可以表示为

G是牛顿的引力常数，ρ代表每一个单元的密度值。x_i、z_i为三角形单元格中每个点的x坐标和z坐标。根据格林定理，可以表示为：

根据积分的定义，坐标系的转动不会影响积分的计算，只会引起坐标的变化。逆时针旋转坐标轴(XOZ)直到旋转坐标系(XOZ)Z轴平行于该边L_ij的外法向量为止(如图3所示)。所以关于一条边L_ij的重力异常g_ij、重力梯度异常g_zij可以表示为

式中G为牛顿引力常数，ρ为每一个单元的密度。θ代表逆时针旋转的角度。

整个三角形的重力梯度异常可以计算如下:

g_z＝g_zij+g_zjk+g_zki

地震声波方程如下：

p(x,t)代表波场值。ρ(x)和v(x)是地下密度和速度分布。s(x_s,t)代表震源函数。

格子法的核心是积分平衡的微分方程弱形式，如图2所示，在k点邻域对声波方程作面积分并应用格林定理可得到：

式中m是围绕节点k的三角单元数,S_kl是节点k周围第l个三角形单元中的虚线段，a和β是虚线段包络线的外法线方向余弦，

是震源函数的面积分。利用动力学计算中的集中质量模型及三角单元线性插值方法，可得到上式的空间离散形式：

其中

式中下标r指代三角单元中的节点，A是三角单元面积，b_k和a_k是三角单元的几何参数，以图4中三角形单元ijk为例，b_k＝(z_i-z_j)/2，a_k＝(x_i-x_j)/2。用二阶中心差分来离散公式左侧项中的时间导数，即可实现波场值的迭代更新。

步骤4：将交叉梯度约束函数引入非规则网格中，实现非规则网格下的二维重力梯度与地震数据的联合反演；

其中：

式中下标r指代三角单元中的节点，A是三角单元面积，b_k和a_k是三角单元的几何参数，以图4中三角形单元ijk为例，b_k＝(z_i-z_j)/2，a_k＝(x_i-x_j)/2。

m⁽¹⁾和m⁽²⁾代表密度和速度参数d_obs ⁽¹⁾和d_obs ⁽²⁾分别代表观测重力梯度数据和地震数据。α₁和α₂代表数据权，用来平衡重力梯度数据和地震波形数据的拟合差。τ₁和τ₂代表正则化系数。

和

分别代表密度和速度部分的模型正则化算子。Φ为密度模型与速度模型的交叉梯度约束函数。

步骤5：进行多次迭代计算，检验反演拟合程度，最终得到的联合反演密度结果如图5所示，联合反演速度结果如图6所示。

最后应说明的是，以上实施例仅用以说明本发明的技术方案而非限制，尽管参照较佳实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，可以对本发明的技术方案进行修改或等同替换，而不脱离本发明技术方案的精神和范围。

Claims (3)

1.一种基于非结构网格下二维重力梯度与地震数据联合反演方法，其特征在于，包括以下步骤：

步骤1：获得重力梯度异常数据；

步骤2：将地下半空间剖分成个若干三角形单元；

步骤3：实现非规则网格下二维重力梯度及地震声波数据的正演计算，然后分别获得不规则网格下二维重力梯度、地震声波数据单独反演结果；

步骤4：将交叉梯度约束函数引入非规则网格中，实现非规则网格下的二维重力梯度与地震声波数据的联合反演；

步骤5：进行多次迭代计算，检验反演拟合程度，直至得到高精度的反演结果。

2.根据权利要求1所述的基于非结构网格下二维重力梯度与地震数据联合反演方法，其特征在于，所述实现非规则网格下二维重力梯度及地震声波数据的正演计算，然后分别获得不规则网格下二维重力梯度、地震声波数据单独反演结果具体为：

所述三角形的重力异常g可以转换为各边贡献的和，表示为

G是牛顿的引力常数,ρ代表每一个单元的密度值。x_i、z_i为三角形单元格中每个点的x坐标和z坐标。根据格林定理，可以表示为：

逆时针旋转坐标轴(XOZ)直到旋转坐标系(XOZ)Z轴平行于该边L_ij的外法向量为止；故关于一条边L_ij的重力异常g_ij、重力梯度异常g_zij表示为

式中G为牛顿引力常数，ρ为每一个单元的密度，θ代表逆时针旋转的角度；

整个三角形的重力梯度异常计算如下:

g_z＝g_zij+g_zjk+g_zki

地震声波方程如下：

p(x,t)代表波场值；ρ(x)和v(x)是地下密度和速度分布；s(x_s,t)代表震源函数；

格子法的核心是积分平衡的微分方程弱形式，在k点邻域对所述地震声波方程作面积分并应用所述格林定理可得到：

式中m是围绕节点k的三角单元数，S_kl是节点k周围第l个三角形单元中的虚线段，a和β是虚线段包络线的外法线方向余弦，

是震源函数的面积分；

利用动力学计算中的集中质量模型及三角单元线性插值方法，可得到上式的空间离散形式：

其中

式中下标r指代三角单元中的节点，A是三角单元面积，b_k和a_k是三角单元的几何参数，用二阶中心差分来离散公式左侧项中的时间导数，即可实现波场值的迭代更新。

3.根据权利要求1所述的基于非结构网格下二维重力梯度与地震数据联合反演方法，其特征在于，所述将交叉梯度约束函数引入非规则网格中，实现非规则网格下的二维重力梯度与地震声波数据的联合反演具体为：

其中：

式中下标r指代三角单元中的节点，A是三角单元面积，b_k和a_k是三角单元的几何参数；

m⁽¹⁾和m⁽²⁾代表密度和速度参数d_obs ⁽¹⁾和d_obs ⁽²⁾分别代表观测重力梯度数据和地震数据；α₁和α₂代表数据权，用来平衡重力梯度数据和地震波形数据的拟合差；τ₁和τ₂代表正则化系数；

和

分别代表密度和速度部分的模型正则化算子；Φ为密度模型与速度模型的交叉梯度约束函数。

Images