CN113109883A - 基于等参变换全球离散网格球坐标下卫星重力场正演方法 - Google Patents

Abstract

本发明公开了一种基于等参变换全球离散网格(dggrid)的球坐标下卫星重力场正演方法，首先选择一定经纬度范围作为观测网范围，并确定观测高度和地球物理模型，利用等面积投影方法，离散地球物理模型，从而获得一系列dggrid网格；此外，基于六边形棱柱体积分公式，计算获得dggrid网格插值形函数、积分点和积分权系数，继而利用等参变换计算所有dggrid网格于各观测点的异常响应。相比于传统经纬度网格的非均匀性和极点奇异性缺点(在极点附近四边形会退化成三角形)，dggrid网格具备一致相邻性特点。与任意四面体解析解对比，基于等参变换全球离散网格(dggrid)的球坐标下卫星重力场正演精度较高。

Description

基于等参变换全球离散网格球坐标下卫星重力场正演方法

技术领域

本发明属于地球物理勘探技术领域，具体涉及基于等参变换全球离散网格球坐标下卫星重力场正演方法。

背景技术

快速高效正演计算是大规模位场反演的基础。当前，重力及其梯度正演方法的有效性、时效性是制约着大规模重力资料反演解释的重要原因。对于局部地区的重力资料而言，可将重力观测面近似为平面,其数据处理问题可在直角坐标系中进行理论推导和计算。然而，对于区域性的,乃至全球尺度问题，重力观测面已不能近似为平面，而是一个球面问题，需要利用球坐标进行表达与处理。于球坐标系下，计算重力位方法可大致分为两类：基于频率域的球谐函数方法和基于空间域的直接积分方法。对于前者，球谐函数方法在中低阶(N<360)相对于数值积分方法计算效率高，但在(超)高阶(N>～2700)会受到数值不稳定性的影响。且球谐函数方法不能适当地考虑给定质量分布的垂直延伸，实际计算结果幅值稍大，不适用于精确地全局或局部重力场建模。

大尺度条件下，地球表面是一个曲面，其反演模型空间的离散化形式也与笛卡尔坐标系下有差别。在球坐标系下，物性反演沿纬度、经度及径向等三个方向分别按照等间隔方式划分为一个个单元块体，这类单元称之为经纬度网格的Tesseroid单元体(Anderson，1976)。经纬度网格的最大优点是网格编码过程简单，且与三维空间坐标的转换易于实现。但随着纬度的增加，从赤道到两极会产生退化现象，在极点附近四边形会退化成三角形，产生大量数据冗余，因此所占内存比较多。

基于以上描述，亟需一种卫星重力场正演方法,以解决现有技术计算过程中产生大量冗余数据，且计算结果不精确的问题。

发明内容

本发明目的在于提供一种基于等参变换全球离散网格球坐标下卫星重力场正演方法,通过该方法可快速获得更高精度的重力、重力矢量及重力梯度张量正演计算结果。

为解决上述技术问题，本发明的一种基于等参变换全球离散网格球坐标下卫星重力场正演方法的具体技术方案如下：

一种基于等参变换全球离散网格球坐标下卫星重力场正演方法，包括以下步骤：

S1、构造具有n个节点的正六边形棱柱；

S2、确定正六边形棱柱在局部坐标系中的空间位置和积分点个数；

S3、构建待分析形函数N：

S4、将n节点正六边形棱柱的顶点坐标代入形函数N，构建如下矩阵形式方程组：

当矩阵C的秩小于n_t时，分析各v_i的相关性，重新遴选

直至C的秩等于n_t，则此时

S5、获得各积分点的积分权重；

S6、获得dggrid网格插值形函数、积分点和积分权系数；

S7、利用等参变换计算所有dggrid网格与各观测点的异常响应。

作为优选，在步骤S1中，正六边形棱柱是由任意n个节点六边形棱柱等参变换而来。

作为优选，在步骤S2中，正六边形棱柱在局部坐标系中的空间位置和积分点个数是根据六边形棱柱的对称性确定的。

作为优选，在步骤S3中，构建待分析形函数N的方法如下：

利用帕斯卡三角形来确定一待用多项式列表V：

V＝{1,x,y,z,xy,…,x⁵…,y⁵…,z⁵…}，

从三维待用多项式列表V选择n_t个多项式

然后构建待分析形函数N：

为便于表述，利用

将N改写为：

作为优选，在步骤S5中，各积分点的积分权重是根据等参变换中雅可比矩阵的定义计算得来的，该矩阵的行列式等于六边形棱柱的体积，即利用各形函数N_i基于六边形棱柱的体积分等于六边形棱柱的体积，构建线性方程组获得。

作为优选，在步骤S7中，异常响应计算方法如下：

球坐标系中任意地质体对其外部空间产生的引力位均可由牛顿万有引力积分公式表达：

其中，G为万有引力常数；Ω地质体所处的积分域；dm与dv为积分质量元与体积元；r为观测点与积分单元的距离；u(x,y,z)为观测点处的引力位；ρ(x',y',z')为地质体的密度或异常密度分布；

根据上述的积分公式对径向r求导，即得到重力或重力异常：

通过u对x和y求偏导得到重力、重力向量或重力梯度张量的各个分量。

有益效果：

本发明充分利用等参变换可将系统坐标系中复杂形体转换为局部坐标系中规则形体的优点，将球壳六边形棱柱的积分区间通过等参变换转换至正六边形棱柱的积分区间，具有以下优点：

1.球壳六边形棱柱具有极强的对称性，这为核函数的重复使用提供了可能，可极大的降低正演核矩阵内存占用。

2.就物性网格离散而言，各球壳六边形棱柱间具有一致相邻性，这避免了传统tesseroid网格因网格位置的不同，而tesseroid网格大小存在显著差异的问题。

3.避免实现传统tesseroid网格积分时需要获得积分根的问题，从而加快正演计算速度。

4.通过本发明提供的计算流程，易于获得更高次积分型函数和积分点，从而获得更高精度的重力、重力矢量及重力梯度张量正演计算结果。

附图说明

图1为本发明的方法流程图；

图2为本发明用于对比验证的四面体剖分立方体示意图；

图3为本发明用于对比验证的立方体正演获得重力矢量和重力梯度张量结果；

图4为本发明用于对比验证的四面体剖分立方体后正演获得重力矢量和重力梯度张量结果；

图5为本发明用于对比验证的立方体及四面体剖分立方体后正演获得重力矢量和重力梯度张量结果的差；

图6为本发明的全球离散网格示意图(线框图)；

图7为本发明的全球离散网格示意图(渲染图)；

图8为本发明的观测示意图及四面体剖分单一全球离散网格示意图；

图9为本发明的用于对比验证四面体剖分单一全球离散网格重力梯度张量分量结果对比图；

图10为本发明的用于对比验证四面体剖分单一全球离散网格重力梯度矢量分量结果对比图。

具体实施方式

为了更好地了解本发明的方法及步骤，下面结合附图，对本发明一种基于等参变换全球离散网格球坐标下卫星重力场正演方法做进一步详细的描述。

首先，引入任意四面体网格正演计算技术，以验证本发明所提正演方法的有效性和可靠性，在此之前，需要验证任意四面体网格正演计算方法的准确性。为此，将与具有精确解析解的长方体单元作为验证任意四面体网格单元正演计算精度的参考，可以分如下几步实现该目的。第一步，以任意四面体网格单元为例，通过三维Hammer积分，获得其重力矢量和重力梯度张量的解析解；第二步，如图2所示，利用Delaunay三角化技术，将长方体单元离散为多个任意四面体单元；第三步，分别计算上述两种正演方法的重力矢量和重力梯度张量值，通过对比图3-图5,可以发现基于三维Hammer积分重力矢量/重力梯度张量解析解具有很高的计算精度。

其次，在此基础上，利用任意四面体网格正演计算方法验证本发明所提正演方法的有效性和可靠性，涉及如下几个步骤：

步骤一，如图6和图7所示，确定研究区域球壳厚度H，沿着球壳半径方向将其等分为n_z层，根据全球离散网格定义，在球面上构建其中的第一级网格(亦可以进一步在球面上构建其中的第n级网格)，其中五边形全球离散网格12个，六边形全球离散网格30个，并利用这些球面上的网格，沿着半径方向逐层离散研究区域球壳，即每一层球壳被离散为五边形棱柱的全球离散网格12个，六边形棱柱的全球离散网格30个，这里将其称之为全球离散网格球壳单元。

步骤二，获得本发明基于等参变换全球离散网格球坐标下卫星重力场正演计算公式，具体方法如下：

球坐标系中任意地质体对其外部空间产生的引力位均可由牛顿万有引力积分公式为：

其中，G为万有引力常数；Ω地质体所处的积分域；dm与dv分别为积分质量元与体积元；r为观测点与积分单元的距离；u(x,y,z)为观测点处的引力位；ρ(x',y',z')为地质体的密度或异常密度分布。

根据上述的积分公式对径向r求导，即得到重力或重力异常：

通过u对x和y求偏导得到重力、重力向量或重力梯度张量的各个分量。

为实现上述公式的积分，首先选择一定经纬度范围作为观测网范围，并确定观测高度和地球物理模型，在等参变换下，构造自n个节点六边形棱柱至n个节点正六边形棱柱的投影。

根据六边形棱柱的对称性确定等参变换的与局部坐标系中的空间位置和积分点个数；

利用帕斯卡三角形确定一待用多项式列表V：

V＝{1,x,y,z,xy,…,x⁵…,y⁵…,z⁵…}；

自三维待用多项式列表V选择n_t个多项式

然后构建待分析形函数N：

为便于表述，利用

将N改写为：

将n节点正六边形棱柱的顶点坐标代入形函数，构建如下矩阵形式方程组：

当矩阵C的秩小于n_t时，分析各v_i的相关性，重新遴选

直至C的秩等于n_t，则此时

根据等参变换中雅可比矩阵的定义，该矩阵的行列式等于六边形棱柱的体积，即利用各形函数N_i基于六边形棱柱的体积分等于六边形棱柱的体积，构建线性方程组，以获得各积分点的积分权重。

基于六边形棱柱体积分公式获得dggrid网格插值形函数、积分点和积分权系数。

步骤三，选择任意一个或多个全球离散网格球壳单元，利用随机数生成函数，对所选全球离散网格球壳单元，赋予随机剩余密度值，以模拟地壳中的复杂地质构造，利用等参变换计算所有dggrid网格于各观测点的重力矢量和重力梯度张量异常响应。

再者，如图8所示，利用Delaunay三角化技术，将上述对所选全球离散网格球壳单元，离散一系列的任意四面体网格单元，进而利用基于三维Hammer积分重力矢量和重力梯度张量解析解，计算上述对所选全球离散网格球壳单元的重力矢量和重力梯度张量异常响应。

可选地，在Delaunay三角化过程中，可限制四面体网格的奇强度不小于1.4，以保证正演计算精度。

最后，通过对比图9和图10,可以发现本发明的重力矢量/重力梯度张量正演结果具有很高的计算精度。

可以理解，本领域技术人员知悉的，在不脱离本发明的精神和范围的情况下，可以对该方法进行各种改变或等效替换。另外，在本发明的教导下，可以对该方法进行修改以适应具体的情况而不会脱离本发明的精神和范围。因此，本发明不受此处所公开的具体实施例的限制，所有落入本申请的权利要求范围内的实施例都属于本发明所保护的范围内。

Claims (6)

1.一种基于等参变换全球离散网格球坐标下卫星重力场正演方法，其特征在于，包括以下步骤：

S1、构造具有n个节点的正六边形棱柱；

S2、确定正六边形棱柱在局部坐标系中的空间位置和积分点个数；

S3、构建待分析形函数N：

S4、将n节点正六边形棱柱的顶点坐标代入形函数N，构建如下矩阵形式方程组：

当矩阵C的秩小于n_t时，分析各v_i的相关性，然后重新遴选

直至C的秩等于n，则此时

S5、获得各积分点的积分权重；

S6、获得dggrid网格插值形函数、积分点和积分权系数；

S7、利用等参变换计算所有dggrid网格与各观测点的异常响应。

2.根据权利要求1所述的基于等参变换全球离散网格球坐标下卫星重力场正演方法，其特征在于，在步骤S1中，正六边形棱柱是由任意n个节点六边形棱柱等参变换而来。

3.根据权利要求1所述的基于等参变换全球离散网格球坐标下卫星重力场正演方法，其特征在于，在步骤S2中，正六边形棱柱在局部坐标系中的空间位置和积分点个数是根据六边形棱柱的对称性确定的。

4.根据权利要求1所述的基于等参变换全球离散网格球坐标下卫星重力场正演方法，其特征在于，在步骤S3中，构建待分析形函数N的方法如下：

利用帕斯卡三角形确定一待用多项式列表V：

V＝{1,x,y,z,xy,…,x⁵…,y⁵…,z⁵…}，

从三维待用多项式列表V中选择n_t个多项式

然后构建如下待分析形函数N：

为便于表述，利用

将N改写为：

5.根据权利要求1所述的基于等参变换全球离散网格球坐标下卫星重力场正演方法，其特征在于，在步骤S5中，各积分点的积分权重是根据等参变换中雅可比矩阵的定义计算得来的，该矩阵的行列式等于六边形棱柱的体积，即利用各形函数N_i基于六边形棱柱的体积分等于六边形棱柱的体积，构建线性方程组获得。

6.根据权利要求1所述的基于等参变换全球离散网格球坐标下卫星重力场正演方法，其特征在于，在步骤S7中，异常响应计算方法如下：

球坐标系中任意地质体对其外部空间产生的引力位均可由牛顿万有引力积分公式表达：

其中，G为万有引力常数；Ω地质体所处的积分域；dm与dv分别为积分质量元和体积元；r为观测点与积分单元的距离；u(x,y,z)为观测点处的引力位；ρ(x',y',z')为地质体的密度或异常密度分布；

根据上述的积分公式对径向r求导，即得到重力或重力异常：

通过u对x和y求偏导得到重力、重力向量或重力梯度张量的各个分量。

Images