CN110967778B - 一种动态坐标系多面体剖分重力布格校正方法 - Google Patents

Abstract

本发明提供一种动力坐标系多面体剖分重力布格校正方法，包括：坐标系的转换：将大地水准上以上的质量体转化到以校正点为原点，以东方向为x方向，以北方向为y方向，以背离地心方向为z方向的动态直角坐标系中；质量多面体剖分：将步骤1中转换后的质量体利用矩形网格切割成多个细小的质量体，根据每一个网格切割成的4点节点的高程特点，可以分为六种情况，对6种情况进行二次剖分剖分成多个四面体；多面体重力值校正量计算：首先确定积分平面的外法线方向，再利用高斯原理坐标转换公式或矢量计算公式计算多面体重力值；验证和实例：采用模拟数据和实际数据进行验证，确定进行布格重力异常校正的正确性。本发明提高了布格校正的精度。

Description

一种动态坐标系多面体剖分重力布格校正方法

技术领域

本发明属于重力布格异常校正技术领域，具体涉及一种动态坐标系多面体剖分重力布格校正方法。

背景技术

布格重力异常是圈定能源矿产、判定断裂系统、确定密度界面起伏、反演三维密度结构等信息，分析构造几何学特征，进而论述运动学、动力学和构造演化的重要数据资料，同时布格重力异常相对正常重力场值而言极小，为确保重力资料获取的地质信息结果准确可靠，保证布格重力数据的精度就显得尤为重要，因此从野外原始数据到布格重力异常的处理中要尽量提高每一步处理的精度。

野外实测重力值受仪器零漂、时间、高度、地形、纬度因素影响，需要移去这些因素从而获取剩余质量所引起的布格重力异常。仪器、时间、高度和纬度的影响处理相对简单，而地形影响较为复杂。Hoyford and Bowie(1912)首先提出要进行起伏地形的校正。Swick(1942)提出经典布格校正(地形影响)分为三部分：中间层校正(布格校正A)、曲率校正(布格校正B)、地形校正(布格校正C)。Cogbill(1979)和LaFehr(1991b)分别提出曲率校正公式。Hammer(1939)把测点周围区域划分为一系列环形域，然后再计算每一小区域对测点影响。Bott(1959)指出计算机可以提高矩形网格的计算效率；Nagy(1966)建立了重力值的直角四棱柱计算公式，现在计算机计算已较易实现地形校正(布格校正C)。地形是重力重要影响因素，地形影响包括中间层校正、曲率校正、地形校正，这些经典处理方法中采用了较多的近似计算，因此影响布格重力异常的精度。

发明内容

本发明提供一种动态坐标系多面体剖分重力布格校正方法，采用一步法，直接把大地水准上以上的质量体转化到动态直角坐标系中，剖分成多个四面体，以解析解计算每一部分的重力异常值进行地形影响校正量计算，以此来减少处理中的近似计算提高精度，提高了布格校正的精度。

本发明的技术方案是：一种动态坐标系多面体剖分重力布格校正方法，包括：

1.坐标系的转换：将大地水准上以上的质量体转化到以校正点为原点，以东方向为x方向，以北方向为y方向，以背离地心方向为z方向的动态直角坐标系中；

2.质量多面体剖分：将步骤1中转换后的质量体利用矩形网格切割成多个细小的质量体，根据每一个网格切割成的4点节点的高程特点，可以分为六种情况，4点全部高程为正，3点为正1点为负，2点为正2点为负(两种)，3点为负1点为正，4点均为负，对6种情况进行二次剖分剖分成多个四面体；

3.多面体重力值校正量计算：首先确定积分平面的外法线方向，再利用高斯原理坐标转换公式或矢量计算公式计算多面体重力值；

4.验证和实例：采用模拟数据和实际数据进行验证，确定进行布格重力异常校正的正确性，所述模拟数据采用简单的凹凸面进行计算，易于判定计算结果的正确性，所述实际数据选取进行过经典校正计算的模型进行验证。

方案进一步地，上述步骤1中大地水准上以上的质量体转化坐标的方程是：

其中，(X_B,Y_B,Z_B)是转换后B点坐标，(x’_A,y’_A,z’_A)，(x’_B,y’_B,z’_B)分别是转换前A点B点的坐标，A点为校正点，θ为校正点的经度，

为校正点的纬度，R为地球半径，h为校正点高程，B为地表上任意一点，θ_B为其经度，

为其纬度，h_B为其高程。

方案进一步地，上述步骤3中积分平面的外法线方向的确定方法是，确定四面体的四个顶点坐标值，利用矢量运算和矩阵运算进行计算最终确定法线方向。

方案进一步地，上述步骤3中高斯原理坐标转换公式的确定过程是，首先把要计算的面进行旋转，x轴和z轴分别旋转α(外法线分量在xoy上的投影与x轴夹角)和β(外法线分量与z轴夹角)，转换后平面三点顶点坐标Z值相同，Z为xyz坐标系转化为XYZ坐标系后，一个面上三角形顶点的Z值，然后面积分转化为线积分，再进行一次坐标旋转，Y轴旋转Ψ角，新的坐标系两轴用η和ξ表示，转化为一边的两顶点的η值相同，Ψ为XOY平面Y轴旋转到η轴的旋转角度，ξ_j+1和ξ_j为一条边两顶点的坐标值，J_j(i)为一边的积分结果，I_i为一个三角形三边积分和，对线进行积分得到高斯原理坐标转换公式是：

方案进一步地，上述步骤3中矢量计算公式为：

I_i＝∑c_ij

其中，投影面为S_i，计算点P投影在投影面上P’，投影面上计算边12，P’在计算边12的投影点为0,0点到1点的距离是l_1ij，0点到2点的距离是l_2ij，v_i为计算点的高程，h_ij为点P’到计算边12的距离，r_0ij、r_1ij、r_2ij分别是计算点到点0、1、2的矢量。

方案进一步地，上述步骤3还包括将两种计算方法进行舍入误差分析、奇异值处理、计算范围分析。

方案进一步地，所述计算范围分析包括根据网格间距确定校正点的计算范围，根据海拔高度，确定计算范围内质量体是否参与计算，若高度在计算点的计算量小于0.01mGal，将不参与计算。

本发明的优点是：本发明采用动态坐标系转换方式，运用角点分割方法，采用多面体矢量公式进行布格校正值的叠加计算，讨论了计算中的舍入误差、奇异值处理、外法线方向确定、计算范围等问题，并通过模拟数据和实际数据的计算，验证了该方法的可行性，减少处理中的近似计算提高精度，提高了布格校正的精度。

附图说明

图1是本发明测量坐标系与动态坐标系转化示意图；

图2是本发明大地水准面以上质量体动态坐标系转换结果图，a所示是一组数据的转换结果图，图b是网格测点把大地水准面以上的质量体切分成了许多小的质量体；

图3是是本发明质量体剖分方法示意图；

图4是本发明外法线向量求取示意图；

图5是本发明重力值求取坐标转换示意图；

图6是本发明矢量计算各量示意图；

图7是本发明舍入误差影响分析的结果图，a是高斯原理坐标转换公式计算结果，b是矢量公式计算结果；

图8是本发明三种公式直角棱柱计算结果对比图，a为直角棱柱计算公式的平面特征和剖面特征，b为基于高斯原理坐标转换计算公式的平面特征和剖面特征，c为矢量公式的平面特征和剖面特征，d为计算模型，e为三剖面叠合结果；

图9是本发明模拟数据验证结果图，a是凸面高程图，b是凸面校正量图，c是凹面高程图，d是凹面校正量图；

图10是本发明选取鄂尔多斯及其周缘重力数据进行了处理的结果图，a是鄂尔多斯及其周缘高程图，b是鄂尔多斯及其周缘自由空间重力异常图，c是鄂尔多斯及其周缘布格重力异常图。

具体实施方式

下面结合附图对本发明做清楚完整的描述，以使本领域的技术人员在不需要作出创造性劳动的条件下，能够充分实施本发明。

本发明的具体实施方式是：一种动态坐标系多面体剖分重力布格校正方法，包括：

1.坐标系的转换：如图1所示，将大地水准上以上的质量体转化到以校正点为原点，以东方向为x方向，以北方向为y方向，以背离地心方向为z方向的动态直角坐标系中，因校正点的变化，因此此直角坐标是动态变化的，随着校正点的不同，坐标原点也不同；

质量体转化坐标的方程是：

其中，θ为校正点的经度，

为校正点的纬度，R为地球半径，h为校正点高程，B为大地水准上以上的质量体，θ_B为其经度，

为其纬度，h_B为其高程；

如图2中a所示是一组数据的转换结果图，校正点中心为(E73,N33)，转化了校正中心周围2度的数据，最底部的平面为大地水准面，上层部分为地表起伏面，中间部分为两面之间的质量体，从图中可见大地水准面以上的质量在转换到为A点为原点的动态坐标系中后，仍然呈现原始球面上的状态，因此没有产生密度分布差异，只是坐标原点成为了校正点A，这对于计算而言，更加简单易操作。从图b可以看出，网格测点把大地水准面以上的质量体切分成了许多小的质量体。

2.质量多面体剖分：将步骤1中转换后的质量体利用矩形网格切割成多个细小的质量体，如图2b所示，根据每一个网格切割成的4点节点的高程特点，可以分为六种情况如图3所示，4点全部高程为正(图3A)，3点为正1点为负(图3B)，2点为正2点为负(图3C和图3D)，3点为负1点为正(图3E)，4点均为负(图3F)，这6种情况的可以适应复杂的地形变化，如陆地、滨海、海域地形特点，对6种情况进行二次剖分剖分成多个四面体，图3A四棱柱可以剖分为两个三棱柱，两个三棱柱可以继续剖分为6个四面体，图3B剖分可剖分为两个三棱柱和一个四面体，图3C可以剖分为两个三棱柱，图3D可以剖分为两个三棱柱和两个四面体，图3E可以两个三棱柱和一个四面体，图3F可以剖分为两个三棱柱，如图3G所示四棱柱可以剖分为两三棱柱，三棱柱可以剖分三个四面体，因此所有质量引起的地形影响，就可以看作所有剖分的四面体引起的引力分量总和；

3.多面体重力值校正量计算：首先确定积分平面的外法线方向，再利用高斯原理坐标转换公式或矢量计算公式计算多面体重力值；

4.验证和实例：采用模拟数据和实际数据进行验证，确定进行布格重力异常校正的正确性，所述模拟数据采用简单的凹凸面进行计算，易于判定计算结果的正确性，所述实际数据选取进行过经典校正计算的模型进行验证。

进一步地，上述步骤3中积分平面的外法线方向的确定方法是，确定四面体的四个顶点坐标值，利用矢量运算和矩阵运算进行计算最终确定法线方向，如图4所示，四面体的四个顶点坐标值为A(x₁,y₁,z₁)，B(x₂,y₂,z₂)，C(x₃,y₃,z₃)，D(x₄,y₄,z₄)。

求取面ABC的外法线方向时，面ABC的上两矢量BA，CA为：

则BA,CA两矢量的叉乘表示面ABC的法线方向，则有：

则面ABC平面方程和过D点垂直面ABC的直线DF方程为:

a(x-x₁)+b(y-y₁)+c(z-z₁)＝0

可以确定线性方程组如下：

为了保证求取面ABC和线DF交点的计算的正确性，比较a、b、c绝对值大小，选取包含绝对值大的组成三元一次线性方程组，求取F点坐标F(x_j,y_j,z_j),则往法线方向n为：

进一步地，上述步骤3中高斯原理坐标转换公式的确定过程是，首先把要计算的面进行旋转，x轴和z轴分别旋转α(外法线分量在xoy上的投影与x轴夹角)和β(外法线分量与z轴夹角)，转换后平面三点顶点坐标Z值相同，Z为xyz坐标系转化为XYZ坐标系后，一个面上三角形顶点的Z值，然后面积分转化为线积分，再进行一次坐标旋转，Y轴旋转Ψ角，新的坐标系两轴用η和ξ表示，转化为一边的两顶点的η值相同，Ψ为XOY平面Y轴旋转到η轴的旋转角度，ξ_j+1和ξ_j为一条边两顶点的坐标值，J_j(i)为一边的积分结果，I_i为一个三角形三边积分和，如图5所示是重力值求取坐标转换示意图，对线进行积分得到高斯原理坐标转换公式是：

进一步地，上述步骤3中矢量计算公式为：

I_i＝∑c_ij

其中，投影面为S_i，计算点P投影在投影面上P’，投影面上计算边12，P’在计算边12的投影点为0,0点到1点的距离是l_1ij，0点到2点的距离是l_2ij，v_i为计算点的高程，h_ij为点P’到计算边12的距离，r_0ij、r_1ij、r_2ij分别是计算点到点0、1、2的矢量，如图6所示是矢量计算各量示意图。

进一步地，上述步骤3还包括将两种计算方法进行舍入误差分析、奇异值处理、计算范围分析。

舍入误差分析的结果如图7所示，是网格间距为0.001度，长四棱柱异常体8顶点经纬度高程值为(72.550,32.551,1000m,72.550,32.550,1000,72.550,32.551,1000,72.550,32.551,100,72.550,32.551,1000m,72.550,32.550,1000,72.550,32.551,1000,72.550,32.551,100)，密度为2730kg/m3,沿32.550纬度的计算结果，图7a所示，随着距离的增大，计算结果不再稳定，甚至超过真值，分析高斯原理坐标转换公式，总共有四项，第一、第二项在对全部积分面积分后和为零，第四项为反正切函数，应用高斯和斯托克斯公式，进行两次坐标旋转后，对于每一条积分线两端点的Z值和η值是相同，第三、第四项问题是相同，都是一个常数乘以一个单调函数的形势，以第三项为例进行说明，第三项是η和ln函数相乘结果，不难看出当两端点相减时变成：

上式中r为计算点到源点积分线端点的距离，当计算点远离源点时，r逐渐增大，ξ_j+1-ξ_j为源的一条积分线两端点的距离，r_j+1-r_j为计算点到源的一条积分线两端点距离差，距离足够大时有(ξ_j+1-ξ_j)＜(r_j+1-r_j)＜＜r_j，因此随着离源的距离越大，η值越大，ln函数值越趋0，设此时值为ω，因计算机有舍入误差(设为σ)，因此计算结果变为η(ω+σ)，因计算距离增大，η增大，ω减小，某一距离时有|ησ|＝ηω，距离继续增大时甚至有|ησ|>ηω。

图7中b是相同模型下运用矢量变形公式计算结果，虽然远离源时，仍然存在舍入误差影响，但误差大小不会超过真值，数值结果是稳定的，公式为：

sign(l₁)＝sign(l₂)：

sign(l₁)≠sign(l₂)：

利用三种公式直角棱柱计算对比，结果如图8所示，d为计算模型，立方体一个角点为坐标原点，各棱边分别平行于坐标轴，顶点坐标如图中所示，平面计算间隔为5m，d异常体平面上的4点恰好为平面计算点，当用矢量公式计算时，r0＝0，v＝0，h＝0，a为直角棱柱计算公式的平面特征和剖面特征，b为基于高斯原理坐标转换计算公式的平面特征和剖面特征，c为矢量公式的平面特征和剖面特征，从平面特征看，除了远距离的舍入误差影响，其他形态完全相同，从平面特征和e三剖面叠合结果看，除了远距的舍入误差影响，其他计算值是基本一致的，在xoy平面内的四个角点恰好是计算点，其r₀，v，h值均为0，此时设c_ij＝0,从计算结果来看，与其它两个公式的结果是相同，说明当r₀＝0，v＝0，h＝0时，c_ij＝0是合理的。

进一步地，所述计算范围分析包括根据网格间距确定校正点的计算范围，根据海拔高度，确定计算范围内质量体是否参与计算，若高度在计算点的计算量小于0.01mGal，将不参与计算，以中纬度地区0.01度网格可以确定计算范围如表1所示，当网格间隔为0.01度，10m高质量体只有在质量体正上方时才参与计算，100m高质量体在校正点0.02度范围内进行计算，500m高质量体在校正点0.07度范围内进行计算，其它以此类推。

表1

如图9所示是通过程序运算获取的数据验证结果，a所示为一个凸面，通过程序运算获取的各点校正量如图b所示，图b同样是中间校正量大、周围小的一系列同心圆，说明校正量与地形起伏成正比关系，与理论相符。图c为一个凹面，通过程序运算获取的各校正量如图d所示，是中间校正小、周围大的一系列同心圆，说明校正量与地形起伏成正关系，与理论相符，通过图9验证了方法的正确性。

为了检验方法对实际数据的处理能力，选取了鄂尔多斯及其周缘重力数据进行了处理，如图10所示，a为该区域的高程数据，通过程序运算获取了该区域的校正量，通过图b的自由空间重力异常与校正量进行运算，获取了该区域的布格重力异常图c所示，把该区域重力布格异常与已有结果进行对比，整体趋势相符，盆山关系吻合好。

以上对本发明的较佳实施例进行了描述，需要指出的是，本发明并不局限于上述特定实施方式，其中未尽详细描述的设备和结构应该理解为用本领域中的普通方式予以实施；任何熟悉本领域的技术人员，在不脱离本发明技术方案范围情况下，依据本发明的技术实质对以上实施例所做的任何简单修改、等同变化及修饰，均仍属于本发明技术方案保护的范围内。

Claims (4)

1.一种动态坐标系多面体剖分重力布格校正方法，其特征在于，包括：

步骤1.坐标系的转换：将大地水准上以上的质量体转化到以校正点为原点，以东方向为x方向，以北方向为y方向，以背离地心方向为z方向的动态直角坐标系中，转化坐标的方程是：

其中，(X_B,Y_B,Z_B)是转换后B点坐标，(x’_A,y’_A,z’_A)，(x’_B,y’_B,z’_B)分别是转换前A点B点的坐标，A点为校正点,θ为校正点的经度，

为校正点的纬度，R为地球半径，h为校正点高程，B为任意地表上一点，θ_B为其经度，

为其纬度，h_B为其高程；

步骤2.质量多面体剖分：将步骤1中转换后的质量体利用矩形网格切割成多个细小的质量体，根据每一个网格切割成的4点节点的高程特点，可以分为六种情况，4点全部高程为正，3点为正1点为负，2点为正2点为负存在两种情况，3点为负1点为正，4点均为负，对6种情况进行二次剖分剖分成多个四面体；

步骤3.多面体重力值校正量计算：首先确定积分平面的外法线方向，再利用高斯原理坐标转换公式或矢量计算公式计算多面体重力值，积分平面的外法线方向的确定方法是，确定四面体的四个顶点坐标值，利用矢量运算和矩阵运算进行计算最终确定法线方向，高斯原理坐标转换公式的确定过程是，首先把要计算的面进行旋转，x轴和z轴分别旋转α和β，α为外法线分量在xoy上的投影与x轴夹角，β为外法线分量与z轴夹角，转换后平面三点顶点坐标Z值相同，Z为xyz坐标系转化为XYZ坐标系后，一个面上三角形顶点的Z值，然后面积分转化为线积分，再进行一次坐标旋转，Y轴旋转Ψ角，新的坐标系两轴用η和ξ表示，转化为一边的两顶点的η值相同，Ψ为XOY平面Y轴旋转到η轴的旋转角度，ξ_j+1和ξ_j为一条边两顶点的坐标值，J_j(i)为一边的积分结果，I_i为一个三角形三边积分和，对线进行积分得到高斯原理坐标转换公式是：

步骤4.验证和实例：采用模拟数据和实际数据进行验证，确定进行布格重力异常校正的正确性，所述模拟数据采用简单的凹凸面进行计算，易于判定计算结果的正确性，所述实际数据选取进行过经典校正计算的模型进行验证。

2.根据权利要求1所述的一种动态坐标系多面体剖分重力布格校正方法，其特征在于，上述步骤3中矢量计算公式为：

其中，投影面为S_i，计算点P投影在投影面上P’，投影面上计算边12，P’在计算边12的投影点为0，0点到1点的距离是l_1ij，0点到2点的距离是l_2ij，v_i为计算点的高程，h_ij为点P’到计算边12的距离，r_0ij、r_1ij、r_2ij分别是计算点到点0、1、2的矢量。

3.根据权利要求2所述的一种动态坐标系多面体剖分重力布格校正方法，其特征在于，上述步骤3还包括将两种计算方法进行舍入误差分析、奇异值处理、计算范围分析。

4.根据权利要求3所述的一种动态坐标系多面体剖分重力布格校正方法，其特征在于，所述计算范围分析包括根据网格间距确定校正点的计算范围，根据海拔高度，确定计算范围内质量体是否参与计算，若高度在计算点的计算量小于0.01mGal，将不参与计算。

Images