一种密度为多项式的任意多面体的重力场正演方法

摘要

本发明公开了一种密度为多项式的任意多面体的重力场正演方法，包括以下步骤：步骤1、确定任意多面体的点和面的坐标、密度多项式表达式和观测点坐标，得到包含体积分的重力位和重力场计算表达式；步骤2、建立局部坐标系，通过矢量恒等式、散度定理和梯度定理，建立体积分的迭代计算方式，将体积分转换为面积分；步骤3、建立局部极坐标系，建立面积分的迭代计算方式，并通过引入固体角，消除面积分中的奇异性，将面积分转换为线积分；步骤4、建立几何关系，通过迭代计算求解线积分，从而得到重力场。本发明能够高精度、快速计算密度随水平方向和垂直方向呈任意阶多项式变化的任意多面体产生的重力位和重力场。

说明书

技术领域

本发明涉及一种地球物理领域的重力正演方法，特别是航空、地面和井中重力勘探等的重力高精度解析计算方法。针对密度为任意高阶多项式的任意多面体在任意观测点处产生的重力位和重力场高精度正演方法。

背景技术

重力勘探是地球勘探方法中的一种重要的、基本的方法，被广泛用于资源勘探、深部地质构造勘探、莫霍界面探测和水文环境物探中。重力勘探包括航空重力勘探、地面重力勘探和井中重力勘探。测量的重力数据和地下地质体的真实密度分布相关，通过进行解析信号分析或者直接借助反演手段，可以得到地下密度分布情况，从而达到解决特定的地质问题的目的。

目前，重力正演常常采用单元积分后进行叠加求和的方式得到重力场，重力解析计算主要针对规则的六面体、多面体和其他具有特殊几何特征的模型，如球体，圆柱体等。常用的网格剖分技术是规则的六面体，并假设每个单元的密度为常数。由于六面体具有简单的几何特性和较容易得到其高精度的解析计算表达式，因此长期以来一直受到人们的青睐。当地下地质体密度分布比较复杂时，六面体的不足日益显著，尤其是无法很好逼近具有复杂几何模型的地质体，。因此，多面体模型逐渐被应用到正演和反演计算中，采用多面体进行网格剖分，最为突出的是四面体单元。同时，由于地下地质体密度分布的不均匀性，如果简单的采用密度为常数的多面体去拟合，那么需要将地质体剖分成众多的小的多面体才能够达到有效逼近真实地质体，无疑加大了计算量，降低了计算效率。因此，众多学者采用密度随水平和垂直方向呈一次、二次和三次多项式变化的多面体进行有效逼近地下复杂地质体，从而达到减少网格剖分个数和提高计算效率的目的。同时，当观测点靠近地质体时，特别是在地质体内部，重力位和重力场计算存在奇异性。因此，有必要设计一种无奇异性的、高精度的密度为多项式的任意多面体的重力场正演方法。

发明内容

本发明所解决的技术问题是，针对现有技术的不足，提供了一种密度为多项式的任意多面体的重力场解析计算方法，能够对多面体在任意观测点处产生的重力位和重力场进行无奇异的高精度、快速的计算求解，提高了计算的精度，为反演提供核心动力。

本发明的技术方案为：

一种密度为多项式的任意多面体的重力场正演方法，包括以下步骤：首先确定任意多面体的点和面的坐标、密度多项式表达式和观测点坐标，得到包含体积分的重力位和重力场计算表达式；然后建立局部坐标系，通过矢量恒等式、散度定理和梯度定理，建立体积分的迭代计算方式，将体积分转换为面积分；再建立极坐标系，建立面积分的迭代计算方式，并通过引入固体角，消除面积分中的奇异性，将面积分转换为线积分；最后建立几何关系，通过迭代的方式计算求解线积分，从而得到重力位和重力场。

具体步骤为：

步骤1、确定任意多面体H的点和面的坐标、密度多项式表达式λ(r)和观测点坐标，将多面体引起的重力位和重力场计算表达式转换为简单的体积分：

如图1所示，对于给定的任意多面体H，其密度λ(r)沿水平方向和垂直方向呈任意阶多项式变化，具体表达式如下所示：

其中，a_pqt为常数，n为密度多项式表达式的阶数，为任意自然数；整数p,q,t分别为在λ(r)在x,y,z方向的阶数，满足0≤p≤n，0≤q≤n和0≤t≤n；那么观测点r'处的重力位φ可以表示为：

其中，G＝6.673×10^-11m³kg^-1s^-2为引力常数，R＝|r-r'|为观测点r'到积分点r的距离，积分点r在多面体H内部，有r∈H；

将方程(2.1)带入方程(2.2)得：

其中，φ_pqt表示如下积分：

观测点r'处的重力场g可以通过对式(2.2)进行求梯度得到：

这里，使用了等式表示求重力位在观测点r'处的梯度。g_pqt表示如下积分：

步骤2、将上述计算表达式中的体积分转换为面积分，从而得到只包含面积分的重力位和重力场计算公式：

首先，建立局部坐标系，设观测点r'为坐标原点，即r'＝(0,0,0)，那么，有：

其中，和分别为x,y,z方向的单位向量；

设多面体H由N个面(多边形)组成，即：

其中，为组成多面体H的面，为第i个面；同时，设第i个面包含M条边，即：

其中，表示构成面的边，C_ij为第i个面的第j条边；

然后，借助矢量恒等式、散度定理和梯度定理，分别对φ_pqt和g_pqt进行处理，得到其迭代计算表达式。

对于φ_pqt：

如果p＞0，借助矢量恒等式，有：

其中表示散度，表示梯度。

进一步地，借助散度定理，当p≥2时，方程(2.4)中的积分φ_pqt可以表示为：

其中，为面的外法线方向，为第i个面的外法线方向。

当p＝1时，φ_1qt可以表示为：

类似地，如果q＞0，那么基于公式(2.8)，借助矢量恒等式和散度定理，当q≥2时，有：

当q＝1时，φ_p1t可以表示为：

类似地，如果t＞0，那么基于公式(2.9)，借助矢量恒等式和散度定理，当t≥2时，有：

当t＝1时，φ_pq1可以表示为：

如果p＝q＝t＝0，由公式(4)可知，φ₀₀₀可以表示为：

对于g_pqt，有：

借助矢量恒等式：

和梯度定理，式(2.6)中的积分可以表示为：

如果p＝0，如果q＝0，如果t＝0，

现在，积分φ_pqt和积分g_pqt分别通过式(2.13)-(2.19)和式(2.21)转换成多个面积分和体积分。将方程(2.13)-(2.19)和方程(2.21)中的体积分写成统一的形式：

I_v＝∫∫∫_Hx^ay^bz^cR^wdv(2.22)

其中，阶数a,b,c为整数，且满足0≤a≤p，0≤b≤q和0≤c≤t，并且w值为±1。

借助假设r'＝(0,0,0)，式(2.22)中的积分核函数可以表示为：

类似地，如果a＞0，则借助矢量恒等式，当a≥2时，有：

当a＝1时，方程(2.22)变为：

类似地，如果b＞0，则借助矢量恒等式，当b≥2时，有：

当b＝1时，有：

如果c＞0，则借助矢量恒等式，当c≥2时，

当c＝1时，有：

当a＝b＝c＝0，方程(2.22)可以表示为：

I_v＝∫∫∫_HR^wdv(2.31)

从式(2.26)可知，当a≥2时，积分I_v转换为对各个面的面积分和一个新的体积分。同时，x的阶数降低2，R的阶数升高2，即w+2。由于w的值为±1，因此，在每一个迭代步骤中，R的阶数一直为奇数。对于任意自然数a≥2，通过次迭代，x的阶数降低为1或者0，其中表示对向下取整；当a＝1时，通过式(2.27)，体积分I_v转换为面积分。同理，利用式(2.28)和式(2.29)分别降低将y和z的阶数降低为1或者0。当a＝b＝c＝0，体积分I_v为一个简单的积分核函数为R^w的体积分。因此，通过循环调用式上述迭代公式计算体积分，降低xyz的阶数【为了减少计算量，优先对x、y和z中阶数较小且阶数不等于零的变量进行降维】，最终将体积分转换为面积分和一个相对简单的体积分，该体积分的积分核函数为R^δ，其中δ≥-1，且为奇数，表示为：

K_v＝∫∫∫_HR^δdv,δ≥-1,(2.32)

计算式(2.32)中的体积分。考虑到矢量恒等式：

结合散度定理，有：

其中，为观测点r'到面的距离。需要注意的是，h_i可以为正数，也可以为负数，取决于观测点r'和面的相对位置。通过式(2.34)，核函数为R^δ的体积分被转换为核函数为R^δ的面积分。

从而，将所有的体积分转换为面积分。

将式(2.13)-(2.19)和式(2.21)中的剩余的面积分、式(2.26)-(2.31)中新的面积分以及式(2.34)中的面积分写成统一的形式：

其中，δ为奇数。式(2.34)中的面积分可以被看为式(2.35)的面积分的一种特殊情况，即α＝β＝γ＝0。

步骤3、面积分转换为线积分

接下来，计算式(2.35)中的面积分。为了便于表述，使用符号Γ表示面如果α＞0，那么，借助式(2.23)中的矢量恒等式，当α≥2时，有：

其中，为边C_ij的外法线方向；

当α＝1时，面积分I_s变成：

类似的，如果β＞0，则运用式(2.24)矢量恒等式，当β≥2时，有：

当β＝1时，有：

如果γ＞0，则运用式(2.24)矢量恒等式，当γ≥2时，有：

当γ＝1时，有：

当α＝β＝γ＝0时，面积分I_s变成：

I_s＝∫∫_ΓR^δds.(2.42)

从式(2.36)可知，对于α≥2，面积分I_s被转化为对各条边的线积分和其他面积分。在第一、第三和第四个面积分中，x的阶数降低1，在第二个面积分中，x的阶数降低2；同时，第三和第四个面积分中，若β或γ等于0，则对应的面积分为0，若β和γ不等于0，则直接计算梯度后得到一个降维的面积分。同理，利用式(2.38)和式(2.40)分别降低将y和z的阶数降低为1或者0。通过循环调用式上述迭代公式，降低xyz的阶数【为了减少计算量，优先对x、y和z中阶数较小且阶数不等于零的变量进行降维】，面积分最终被转换为一系列的线积分和面积分，线积分的积分核函数为x^μy^νz^λR^ξ(0≤μ≤α，0≤ν≤β，0≤λ≤γ，1≤ξ≤γ)，面积分的积分核函数为R^ζ(ζ≥-1，且为奇数)。该面积分形式与式(2.34)和式(2.42)中的面积分形式一致，因此可以同时一次处理。将上述面积分和线积分写成统一的形式，有：

K_s＝∫∫_ΓR^ζds,ζ≥-1(2.43)

其中，ζ和ξ为奇数。

首先，推导式(2.43)中的面积分。如果ζ＝-1，当R→0，即观测点r'靠近多面体H，积分∫∫_ΓR^-1ds存在奇异性，因此，需要特殊处理。

为了消除∫∫_ΓR^-1ds的奇异性，高精度地计算面积分，在各个面上建立极坐标系如图1所示，为面的外法线方向，点o为观测点r'在面的投影点，为点o在面上的固体角，其中β(o)_ij为边C_ij对应的固体角，即点o和边C_ij的两个顶点构成的角。和分别为边C_ij的外法线方向和切线方向，满足边C_ij使用变量s_ij参数化，即v_ij0和v_ij1分别为边C_ij的两个顶点。R_ij0和R_ij1分别为观测点r'到两个顶点的距离。点r到点o的距离用ρ＝|r-o|表示。

因此，观测点r'到点r的距离可以表示为那么，存在如下矢量恒等式：

其中，为面散度符号。

那么，式(2.43)中的面积分可以表示为：

从而将所有的面积分转换为线积分。

其中，在边C_ij上为一常数。ρ＝|r-o|表示点r到点o的距离。β(o)为固体角，当点o在面内部时，β(o)＝2π；当点o在边上时，β(o)＝π；当点o在顶点上时，β(o)＝θ，其中θ为该点所在的两条边组成的角的度数；当点o面的外部时，β(o)＝0。此时，∫∫_ΓR^ζds不存在奇异性，这是因为当β(o)≠0时，积分存在奇异性，然而奇异性被考虑到固体角中。固体角可以通过如下公式进行计算：

其中，为点o到点r_⊥的单位向量。点r_⊥为观测点r'在边C_ij的投影点。

步骤4、计算观测点重力位和重力场：

采用迭代的方式计算步骤3中式(2.46)中的线积分：

为了迭代计算式(2.48)，需要初始积分可以通过如下表达式进行计算：

并且，式(2.48)中的线积分可以采用迭代公式进行计算求解：

其中，初始线积分为：

最后，计算式(2.44)中的线积分。如图1所示，存在如下几何关系：

其中，w₁＝x,w₂＝y,w₃＝z和那么，线积分可以表示为：

通过二项展开，式(2.53)中的线积分包含一系列的简单的线积分。忽略各线积分的系数，可以将之写成统一的形式，并通过迭代方式进行计算求解：

其中，1≤u≤μ+ν+λ，u为s_ij的阶数。方程(2.54)表明每一次迭代，s_ij的阶数降低2。因此，如果s_ij为奇数，经过式(2.54)迭代计算，线积分最后变成核函数为s_ijR^ξ的线积分，可以通过如下计算表达式进行计算：

如果s_ij为偶数，那么核函数为线积分最后变成核函数为R^ξ的线积分，可以通过式(2.50)进行迭代求解。

因此，可以通过上述计算公式，本发明给出了密度为任意阶多项式的任意多面体的重力场高精度正演方法。

有益效果：

本发明通过建立局部坐标系，借助散度定理和梯度定理，将重力位和重力场中包含的体积分转换为面积分，再通过建立极坐标系，消除了积分中存在的奇异性，并将面积分转化为线积分，最后推导了各个积分的迭代计算表达式。本发明能够高精度计算密度随水平方向和垂直方向呈任意阶多项式变化的任意多面体产生的重力位和重力场，为快速重力正演和反演提供核心动力。

本发明既可以计算在地质体外部的重力位和重力场，也可以用于计算在异常体表面和内部任意点出的重力位和重力场。为航空、地面和井中重力计算提供了很好的计算方法，为正演和反演提供了核心动力。

附图说明

图1为任意多面体示意图。r'为观测点，H为多面体，为第i个面，C_ij为第i个面的第j条边，h_i为观测点r'到面的距离，点o为观测点r'在面上的投影点，为面的外法线方向，和分别为边C_ij的外法线方向和切线方向，点r_⊥为观测点r'在边C_ij的投影点，v_ij0和v_ij1分别为边C_ij的两个顶点，R_ij0和R_ij1分别为观测点r'到两个顶点的距离。

图2为正十二体示意图。正十二面体的边长为虚线表示测线。

图3为正十二面体在测线上产生的重力场。实线为采用高斯数值积分计算的结果，圆圈表示采用本发明的方法计算的结果。图3(a)为g_x分量对比结果；图3(b)为g_y分量对比结果；图3(c)为g_z分量对比结果；图3(d)为二者的相对误差。

具体实施方式

以下结合附图和具体实施方式对本发明作进一步的说明。

本发明涉及的重力计算方法包括以下步骤：

步骤1、多面体引起的重力位和重力场计算表达式转换为简单的体积分。

步骤2、体积分转换为面积分：建立局部坐标系，通过矢量恒等式、散度定理和梯度定理，推导任意高阶多项式的体积分的迭代计算方式，并将之转换为面积分。

步骤3、面积分转换为线积分：、建立局部极坐标系，推导面积分的迭代计算求解方式，并通过引入固体角，消除面积分中的奇异性，将面积分转换为线积分。

步骤4、计算观测点重力位和重力场：建立几何关系，通过迭代的方式计算求解线积分，从而正演得到密度为多项式的任意多面体的重力场。

以下为本发明计算一个正十二面体重力场的实例。

如图2所示，设单一正十二面体的二十个顶点坐标为和其中，所有数值的单位为km。正十二面体的密度随水平方向和垂直方程呈多项式变化：

λ(r)＝50+2.6x+1.8y+3.2z+26.27z²+1.6z³(2.56)

其中，密度的单位为kg/m³。50个观测点均匀分布在x＝25km到x＝1250km测线上(y＝z＝0km)。分别采用本发明方法和高阶高斯数值积分方法(64×64×64点)计算重力场。计算结果如图3所示。实线为采用高斯数值积分计算的结果，圆圈表示采用本发明的方法计算的结果。图(a)为g_x分量对比结果，图(b)为g_y分量对比结果，图(c)为g_z分量对比结果，二者的结果十分一致，表明本发明方法的正确性。图(d)为二者的相对误差，表明在观测点距离目标体的距离是目标体维度的200倍范围内，本发明方法可以高精度计算目标体在观测点出产生的重力场。