CN106855904A - 一种二度体重力异常计算方法 - Google Patents

Abstract

本发明提出了一种二度体重力异常计算方法，本发明通过复杂二度体模型表示和矩形模型重力异常计算步骤，其中矩形模型重力异常计算包括加权系数计算、一维离散卷积计算、重力异常值合成步骤。本发明实现了重力异常计算在效率和精度上的统一。本发明解决了现有二度体重力异常计算方法不能同时保证计算效率和计算精度，无法满足大规模重力异常数据密度精细反演、人机交互建模和解释的需求的问题。

Description

一种二度体重力异常计算方法

技术领域

本发明涉及一种二度体重力异常计算方法，特别是任意截面形状、任意密度分布的二度体重力异常快速、高精度计算方法。

背景技术

重力勘探以地球重力场为根据、密度差异为物理基础，是一种研究地球构造及寻找矿产资源的地球物理勘探方法，该方法效率高、操作简单、成本低、勘探深度大、实施过程没有过多条件限制。重力异常正演计算是指根据密度分布计算重力异常，反演是指根据观测重力值计算密度分布。正演是反演的基础，正演计算的效率直接影响反演计算的效率，而正演计算精度直接影响反演结果的质量。一直以来，学者们都非常重视重力异常的正演计算。实际情况中，常见这样一类线性地质体：其走向方向的尺度要远比垂直其走向方向的尺度大，诸如此类地质体，实际场源分布便可以用走向方向无限延伸的二度体代替。

伴随测绘技术和仪器的发展，重力测量在测量精度、空间分辨率和测量范围上都有了显著提高，为重力反演提供了大规模高精度、高分辨率数据。同时，随着计算机软硬件水平的不断提高，人机交互建模、解释也日益得到人们的重视，在地球物理勘探中发挥着越来越重要作用。人机交互建模能够使人们通过直观的方式对地质体进行建模，更容易融合地下结构的先验信息。反演方法与人机交互建模、解释方法相辅相成，将极大提高人们对地球内部结构的认识。

对于任意截面形状、任意密度分布的二度体重力异常计算，一般采用数值方法。文献(张岭,郝天珧.基于Delaunay剖分的二维非规则重力建模及重力计算.地球物理学报,2006,49(3):877-884.)公开了一种Delaunay剖分方式，将截面分为若干三角形，将二度体分解为若干三角棱柱的组合，利用解析公式计算变密度三角棱柱重力异常，并将其累加，这种方法保证了计算精度，但计算效率低；文献(朱自强,曾思红,鲁光银等.二度体的重力张量有限元正演模拟.物探与化探,2010,34(5):668-671.)公开了一种计算二度体重力梯度张量的有限元方法；文献(Reeder K,Louie J,Kent G.Efficient 2D finite elementgravity modeling using convolution.SEG Technical Program Expanded Abstracts,2014:1254-1258.)公开了一种有限元方法和卷积算法相结合的二度体重力异常计算方法。采用有限元方法计算二度体重力异常，能够较精确刻画复杂二度体，计算精度高但计算效率很低。

剖分方式和计算方法共同决定了重力异常计算的效率和精度。计算效率和计算精度是一对矛盾体，目前已有的方法存在的最大问题是不能同时保证计算效率和计算精度，无法满足大规模重力异常数据密度精细反演、人机交互建模和解释的需求。因此，寻找一种计算效率高、同时能保证计算精度的计算方法具有重要的现实意义。

发明内容

本发明针对目前现有二度体重力异常计算方法不能同时保证计算效率和计算精度，无法满足大规模重力异常密度精细反演成像、人机交互建模和解释的需求问题，本发明其目的在于提出一种二度体重力异常计算方法，其是一种任意截面形状、任意密度分布的二度体重力异常快速、高精度计算方法。

本发明的技术方案是：

一种二度体重力异常计算方法，包括以下步骤：

第一步复杂二度体模型表示：

确定观测点，观测点坐标(x_m,z₀)，基于观测点建立包含所有目标区域的矩形模型，确定矩形模型在x，z方向的起始位置，使得包含起伏地形的目标区域完全嵌入在该矩形模型中；

然后将该矩形模型均匀划分成若干个规则小矩形，确定小矩形的几何尺寸Δx，Δz；

最后根据目标区域的密度分布，对每个小矩形密度进行赋值，每个小矩形密度为常值，不同矩形密度取值不同，以此刻画任意截面形状、任意密度分布的二度体；将位于空气部分的小矩形的密度值设为零，以此刻画起伏地形；

第二步矩形模型重力异常计算；

第一步中给出的矩形模型，其重力异常计算公式为

式中，(x_m,z₀)表示观测点坐标，z₀为常值；L表示矩形模型在z方向剖分小矩形的个数；M表示矩形模型在x方向剖分小矩形的个数；(ξ_p,ζ_r)表示编号为(p,r)的小矩形几何中心坐标即在矩形模型中x方向为第p个且在z方向为第r个的小矩形的几何中心坐标；ρ(ξ_p,ζ_r)表示编号为(p,r)的小矩形的密度值；h(x_m-ξ_p,z₀-ζ_r)表示加权系数。

进一步地，第二步中，矩形模型其重力异常计算公式的解算方法如下：

a，根据观测点坐标(x_m,z₀)和小矩形几何中心坐标(ξ_p,ζ_r)，计算加权系数h(x_m-ξ_p,z₀-ζ_r)，其计算公式为

式中，γ表示万有引力常数，arctan()表示反余切函数运算符，ln()表示自然对数运算符；其它符号含义如下

X₁＝ξ_p-0.5Δx-x_m，X₂＝ξ_p+0.5Δx-x_m，Z₁＝ζ_r-0.5Δz-z₀，Z₂＝ζ_r+0.5Δz-z₀

Δx，Δz表示小矩形几何尺寸。

b，在矩形模型中在x方向上位于同一行的所有小矩形为一层。采用一维离散卷积快速计算方法来计算一层矩形模型重力异常，其计算公式为

式中，表示第r层(r＝1,2,…,L)矩形模型在测线z₀产生的重力异常；(x_m,z₀)表示观测点坐标；

c，将各层矩形模型重力异常进行累加，得到整个矩形模型的重力异常，即

在第二步中的步骤b中，采用一维离散卷积快速计算方法来计算一层矩形模型重力异常，其步骤为：

(1)将加权系数h(x₁-ξ_p,z₀-ζ_r)排列成向量t，记为

t＝[t₀,t₁,t₂,…,t_M-1,0,t_M-1,t_M-2,…,t₂,t₁]^T (5)

式中，矩阵元素t_i与加权系数h(x₁-ξ_p,z₀-ζ_r)存在关系

t_i＝h(x₁-ξ_i+1,z₀-ζ_r) (6)

(2)将第r层密度值ρ(ξ_p,ζ_r)(p＝1,2,…,M)排列成向量ρ，向量元素ρ_i与密度值存在关系

ρ_i＝ρ(ξ_i,ζ_r) (7)

将向量ρ补零扩展成向量ρ^ext

式中，0_M×1表示M×1零向量；

(3)计算

式中，fft()表示一维快速傅里叶变换；

(4)计算

式中，“.*”表示对应元素相乘运算；

(5)计算

式中，ifft()表示一维快速傅里叶反变换；

(6)提取矩阵g^ext的前M行元素，构成向量g，即为一维离散卷积计算结果，向量g中元素就是所求

与现有技术相比，本发明具有以下优点：

(1)模型表示方法简单、灵活，很容易刻画任意截面形状、任意密度分布复杂二度体以及起伏地形；

(2)能够实现任意截面形状、任意密度情况下复杂二度体重力异常的快速、高精度计算，可以满足大规模重力密度精细反演、人机交互建模和解释的需求；

(3)大规模正演计算时，算法不但计算效率和计算精度高，并且所需计算机内存小。

附图说明

图1为本发明的计算流程图；

图2为复杂二度体模型表示；

图3为截面为圆形二度体模型图；

图4为重力异常计算值和理论值对比图；

图5为计算值与理论值的相对误差；

图中符号说明如下：

ρ：表示密度。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明实施方式作进一步地详细描述。

参照图1，为本实施例一种二度体重力异常计算方法的流程图，包括以下步骤：

1、复杂二度体模型表示：

首先，确定观测点，观测点坐标(x_m,z₀)，基于观测点建立包含所有目标区域的规则矩形模型，确定矩形模型在x，z方向的起始位置，使得目标区域(包含起伏地形)完全嵌入在该矩形模型中；

其次，根据实际问题需求即根据预定的针对目标区域确定的小矩形剖分个数，将矩形模型均匀划分成多个规则小矩形(如图2所示)，确定小矩形的几何尺寸Δx，Δz；

最后，根据已知的目标区域的密度分布，对每个小矩形密度进行赋值，根据小矩形和二度体截面的几何关系，当小矩形中心位于二度体截面内时，赋值给所在位置的二度体密度，每个小矩形密度为常值，不同小矩形密度取值不同，以此刻画任意截面形状、任意密度分布的二度体。位于空气部分的小矩形，其密度值设为零，以此刻画起伏地形。

2、矩形模型重力异常计算：

步骤一中给出的矩形模型，其重力异常计算公式为

式中，(x_m,z₀)表示观测点坐标，z₀为常值；L表示矩形模型在z方向剖分小矩形的个数；M表示矩形模型在x方向剖分小矩形的个数；(ξ_p,ζ_r)表示编号为(p,r)的小矩形几何中心坐标即在矩形模型中x方向为第p个且在z方向为第r个的小矩形的几何中心坐标；ρ(ξ_p,ζ_r)表示编号为(p,r)的小矩形的密度值；h(x_m-ξ_p,z₀-ζ_r)表示加权系数。

实现式(1)的计算，分为三个环节：

首先，根据观测点坐标(x_m,z₀)和小矩形几何中心坐标(ξ_p,ζ_r)，计算加权系数h(x_m-ξ_p,z₀-ζ_r)，其计算公式为

式中，γ表示万有引力常数，arctan()表示反余切函数运算符，ln()表示自然对数运算符；其它符号含义如下

X₁＝ξ_p-0.5Δx-x_m，X₂＝ξ_p+0.5Δx-x_m，Z₁＝ζ_r-0.5Δz-z₀，Z₂＝ζ_r+0.5Δz-z₀

Δx，Δz表示小矩形几何尺寸。

其次，在矩形模型中在x方向上位于同一行的所有小矩形为一层，采用一维离散卷积快速计算方法来计算一层矩形模型重力异常，其计算公式为

式中，表示第r层(r＝1,2,…,L)矩形模型在测线z₀产生的重力异常；(x_m,z₀)表示观测点坐标；

最后，将各层矩形模型重力异常进行累加，得到整个矩形模型的重力异常，即

步骤二中所述的一维离散卷积快速计算方法，其步骤为：

(1)将加权系数h(x₁-ξ_p,z₀-ζ_r)排列成向量t，记为

t＝[t₀,t₁,t₂,…,t_M-1,0,t_M-1,t_M-2,…,t₂,t₁]^T (5)

式中，矩阵元素t_i与加权系数h(x₁-ξ_p,z₀-ζ_r)存在关系

t_i＝h(x₁-ξ_i+1,z₀-ζ_r) (6)

(2)将第r层密度值ρ(ξ_p,ζ_r)(p＝1,2,…,M)排列成向量ρ，向量元素ρ_i与密度值存在关系

ρ_i＝ρ(ξ_i,ζ_r) (7)

将向量ρ补零扩展成向量ρ^ext

式中，0_M×1表示M×1零向量；

(3)计算

式中，fft()表示一维快速傅里叶变换；

(4)计算

式中，“.*”表示对应元素相乘运算；

(5)计算

式中，ifft()表示一维快速傅里叶反变换；

(6)提取矩阵g^ext的前M行元素，构成向量g，即为一维离散卷积计算结果，向量g中元素就是所求

下面对本发明方法的效果进行检验。

为了说明本发明所提出的方法用于计算任意截面形状、任意密度分布情况下复杂二度体重力异常时的效率和精度，设计了如下二度体组合模型(图3所示)：

密度为常值的矩形区域内嵌一个密度均匀的圆形，圆心与矩形中心重合。矩形范围为：x方向从-1000m到1000m，z方向从0m到1000m(z轴向下为正)；圆形半径为400m。矩形密度为0g/cm³，圆形的密度为1g/cm³。将矩形剖分成10000×10000个大小相同的小矩形，计算高度为-50m测线(图3中虚线所示)上的重力异常，计算点个数为10000。

重力异常算法利用Fortran语言编程实现，运行程序所用的个人台式机配置为：CPU为i7-2620，主频为2.7GHz，内存为32GB，四核八线程。运行所需时间约为6秒，由此可见算法效率很高。重力异常计算值和理论值如图4所示，从形态上看，两者是一致的。相对误差由理论值减去计算值得到差值的绝对值除以理论值得到(图5所示)，对相对误差进行统计，统计结果由表1给出，可知算法精度很高。

表1重力异常理论值和计算值相对误差统计

本发明是一个有机整体，即在特定的模型表示方式条件下，建立矩形二度体重力异常叠加模型，根据一种特殊的加权系数计算公式，采用一维离散卷积快速计算方法，实现了重力异常计算在效率和精度上的统一。

以上包含了本发明优选实施例的说明，这是为了详细说明本发明的技术特征，并不是想要将发明内容限制在实施例所描述的具体形式中，依据本发明内容主旨进行的其他修改和变型也受本专利保护。本发明内容的主旨是由权利要求书所界定，而非由实施例的具体描述所界定。

Claims (3)

1.一种二度体重力异常计算方法，其特征在于，包括以下步骤：

第一步 复杂二度体模型表示：

确定观测点，观测点坐标(x_m,z₀)，基于观测点建立包含所有目标区域的矩形模型，确定矩形模型在x，z方向的起始位置，使得包含起伏地形的目标区域完全嵌入在该矩形模型中；

将该矩形模型均匀划分成若干个小矩形，确定小矩形的几何尺寸Δx，Δz；

最后，根据目标区域的密度分布，对每个小矩形密度进行赋值，每个小矩形密度为常值，不同小矩形密度取值不同，以此刻画任意截面形状、任意密度分布的二度体；将位于空气部分的小矩形的密度值设为零，以此刻画起伏地形；

第二步 矩形模型重力异常计算；

第一步中给出的矩形模型，其重力异常计算公式为

$g_z(x_m,z_0)=\underset{r=1}{\overset{L}{\Σ}}\underset{p=1}{\overset{M}{\Σ}}\mathrm{\ρ}({\mathrm{\ξ}}_p,{\mathrm{\ζ}}_r)h(x_m-{\mathrm{\ξ}}_p,z_0-{\mathrm{\ζ}}_r)---\left(1\right)$

式中，(x_m,z₀)表示观测点坐标，z₀为常值；L表示矩形模型在z方向剖分小矩形的个数；M表示矩形模型在x方向剖分小矩形的个数；(ξ_p,ζ_r)表示编号为(p,r)的小矩形几何中心坐标；ρ(ξ_p,ζ_r)表示编号为(p,r)的小矩形的密度值；h(x_m-ξ_p,z₀-ζ_r)表示加权系数。

2.根据权利要求1所述的二度体重力异常计算方法，其特征在于：第二步中，矩形模型其重力异常计算公式的解算方法如下：

a，根据观测点坐标(x_m,z₀)和小矩形几何中心坐标(ξ_p,ζ_r)，计算加权系数h(x_m-ξ_p,z₀-ζ_r)，其计算公式为

$h(x_m-{\mathrm{\ξ}}_p,z_0-{\mathrm{\ζ}}_r)=-\mathrm{\γ}\left[\begin{array}{c}{X}_{2}ln\left(\frac{X_2^2+Z_2^2}{X_2^2+Z_1^2}\right)+X_{1}ln\left(\frac{X_1^2+Z_1^2}{X_1^2+Z_2^2}\right)\\ +2\left|Z_2\right|(\arctan \frac{X_2}{\left|Z_2\right|}-\arctan \frac{X_1}{\left|Z_2\right|})\\ +2\left|Z_1\right|(arctan\frac{X_1}{\left|Z_1\right|}-\arctan \frac{X_2}{\left|Z_1\right|})\end{array}\right]---\left(2\right)$

式中，γ表示万有引力常数，arctan()表示反余切函数运算符，ln()表示自然对数运算符；其它符号含义如下

X₁＝ξ_p-0.5Δx-x_m，X₂＝ξ_p+0.5Δx-x_m，Z₁＝ζ_r-0.5Δz-z₀，Z₂＝ζ_r+0.5Δz-z₀

Δx，Δz表示小矩形几何尺寸；

b，在矩形模型中在x方向上位于同一行的所有小矩形为一层，采用一维离散卷积快速计算方法来计算一层矩形模型重力异常，其计算公式为

$g_z^r(x_m,z_0)=\underset{p=1}{\overset{M}{\Σ}}\mathrm{\ρ}({\mathrm{\ξ}}_p,{\mathrm{\ζ}}_r)h(x_m-{\mathrm{\ξ}}_p,z_0-{\mathrm{\ζ}}_r)---\left(3\right)$

式中，表示第r层矩形模型在测线z₀产生的重力异常，其中r＝1,2,…,L；(x_m,z₀)表示观测点坐标；

c，将各层矩形模型重力异常进行累加，得到整个矩形模型的重力异常，即

$g_z(x_m,z_0)=\underset{r=1}{\overset{L}{\Σ}}{g}_z^r(x_m,z_0)---\left(4\right).$

3.根据权利要求2所述的二度体重力异常计算方法，其特征在于：在第二步中的步骤b中，采用一维离散卷积快速计算方法来计算一层矩形模型重力异常，其步骤为：

(1)将加权系数h(x₁-ξ_p,z₀-ζ_r)排列成向量t，记为

t＝[t₀,t₁,t₂,…,t_M-1,0,t_M-1,t_M-2,…,t₂,t₁]^T (5)

式中，矩阵元素t_i与加权系数h(x₁-ξ_p,z₀-ζ_r)存在关系

t_i＝h(x₁-ξ_i+1,z₀-ζ_r) (6)

(2)将第r层密度值ρ(ξ_p,ζ_r)(p＝1,2,…,M)排列成向量ρ，向量元素ρ_i与密度值存在关系

ρ_i＝ρ(ξ_i,ζ_r) (7)

将向量ρ补零扩展成向量ρ^ext

${\mathrm{\ρ}}^{ext}=\left[\begin{array}{c}\mathrm{\ρ}\\ 0_{M\×1}\end{array}\right]---\left(8\right)$

式中，0_M×1表示M×1零向量；

(3)计算

式中，fft()表示一维快速傅里叶变换；

(4)计算

式中，“.*”表示对应元素相乘运算；

(5)计算

式中，ifft()表示一维快速傅里叶反变换；

(6)提取矩阵g^ext的前M行元素，构成向量g，即为一维离散卷积计算结果，向量g中元素就是所求