CN114529688B - 一种基于六边形网格索引下的多分辨率多面体计算全球地形改正的方法 - Google Patents

Abstract

本发明公开一种基于六边形网格索引下的多分辨率多面体计算全球地形改正的方法，包括：基于二十面体构建7孔剖分的球面六边形网格系统H3；基于H3及各网格参数构建H3下的地形分区策略；利用六边形网格与三角面之间的对应关系，确定多分辨率地形的多面体模型，顾及海水质量密度亏损，利用多面体积分方法，确定待计算点处所受的地形改正影响；利用构建的多面体模型及其对应海水质量密度，使用多面体积分方法，确定待计算点处的各多面体模型对应的引力值，多组引力值之和即为最终该待计算点处的地形改正值。本发明具有更小的地形改正的长波长误差，同时还解决了利用多分辨地形模型实现局部区域地形改正的高精度计算。

Description

一种基于六边形网格索引下的多分辨率多面体计算全球地形 改正的方法

技术领域

本发明属于大地测量与测绘工程、固体地球物理技术领域，尤其涉及一种基于六边形网格索引下的多分辨率多面体计算全球地形改正的方法。

背景技术

地形对于相关地球重力场元(例如重力异常、扰动重力、大地水准面高、垂线偏差、扰动位的二阶梯度等)有重要影响，该影响的精确计算在大地测量的应用和大地水准面的构建中发挥着重要作用，对于地球物理学中反演地球内部密度和地质构造方面也具有重要作用。地形改正多采用棱柱和Tesseroid方法，随着全球地形模型分辨率不断提高至优于30m，计算耗时和效率直接影响了大地水准面精化和重力场资料处理的效率。

发明内容

本发明针对大地水准面精化和地球重力场确定过程中，计算耗时和效率直接影响了大地水准面精化和重力场资料处理的效率的问题，提出一种基于六边形网格索引下的多分辨率多面体计算全球地形改正的方法，不仅有效顾及了地球曲率、全球地形质量影响及海水质量密度亏损，还解决了利用多分辨地形模型实现局部区域地形改正的高精度计算，相比以往棱柱积分方法、点质量、线质量、Tesseroid方法，计算效率明显提高。

为了实现上述目的，本发明采用以下技术方案：

本发明提出一种基于六边形网格索引下的多分辨率多面体计算全球地形改正的方法，包括：

步骤1：基于二十面体构建7孔剖分的球面六边形网格系统H3，利用该网格系统确定网格分辨率、网格数量、网格中点坐标、网格边界坐标、网格相互关系；

步骤2：基于六边形网格系统H3及该系统中各网格参数构建六边形网格系统H3下的地形分区策略；

步骤3：利用六边形网格与三角面之间的对应关系，确定多分辨率地形的多面体模型，顾及海水质量密度亏损，利用多面体积分方法，确定待计算点处所受的地形改正影响；

步骤4：利用构建的多分辨率地形的多面体模型及其对应海水质量密度，使用多面体积分方法，确定待计算点处的各多面体模型对应的引力值，多组引力值之和即为最终该待计算点处的地形改正值。

进一步地，所述六边形网格系统H3相邻分辨率中大的六边形网格包含7个小的六边形网格，大六边形和小六边形的方向存在19.1°的转角，利用该网格系统可以进行不同分辨率六边形网格之间的快速索引。

进一步地，所述步骤2包括：

基于六边形网格系统H3及该系统中各网格参数、按照如下方式构建六边形网格系统H3下的地形分区策略：

区域1：区域半径为60km，网格分辨率为0.174km；

区域2：区域半径为158km，网格分辨率为0.46km；

区域3：区域半径为418km，网格分辨率为1.2km；

区域4：区域半径为1107km，网格分辨率为3.2km；

区域5：区域半径为全球，网格分辨率为8.5km。

进一步地，所述步骤3包括：

根据六边形网格剖分与三角形剖分之间的对偶关系，利用六边形网格对三角形进行检索和管理，重建六边形网格分布下的规则三角网地形，基于多面体积分方法求解该地形质量产生的引力作用；利用H3六边形网格索引系统，生成全球分布的六边形网格，每个六边形网格的中心点与6个边界点构成6个近似等边三角形，从而形成六边形网格与三角网格之间的对应关系；

根据上述关系，利用多分辨率六边形网格剖分，得到以三角面组成的多分辨率地形的多面体模型，以利用多面体积分方法求解地形改正：当求解顾及全球地形质量以及海水质量密度亏损下的地形改正时，将陆地平均地壳密度设为2.67g/cm³，海水密度为1.03g/cm³，利用上述三角形面组成的点位数据，结合全球1′的地形数据ETOPO1以及30m的局部区域地形数据ASTER，构建出3套不同海水密度的多面体模型，分别为，以大地水准面高加上全球地形模型确定的地形大地高为边界、密度为2.67g/cm³的多面体，陆地区域以大地水准面为外边界、海洋区域以实际地形为边界、密度为-1.03g/cm³的多面体，以整个大地水准面为外边界的密度为-1.64g/cm³的多面体。

进一步地，所述步骤4包括：

设多面体共由n个三角面ABC组成，则待计算点Q所受地形改正的引力影响表达式为：

其中

其中d_As为奇异项改正值，表示第p个三角面的外法线方向单位矢量，/>是以Q点在第p个三角面上的垂足点Q′为原点建立的用于描述待求解引力矢量三分量坐标系的三个轴单位矢量，h_p是计算点Q到第p个三角面的距离，h_pq是计算点Q在第p个三角面上的垂足Q′到第p个三角面上第q条边的距离；σ_pq表示Q点的垂足点Q′在直线AB外法线一侧为1，另一侧为-1；G表示万有引力常数，ρ表示多面体质量的密度；S₁、S₂、l₁、l₂分别表示AQ″、BQ″、AQ、BQ的距离，Q″表示Q′在第p个三角面第q条边上的垂足点。

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

本发明的一种基于六边形网格索引下的多分辨率多面体计算全球地形改正的方法，利用多面体方法实现顾及全球高分辨率地形的地形改正计算，并利用六边形网格与三角形之间的对偶关系，实现多种分辨率的多面体构建，从而高效的实现局部区域地形改正的高精度计算。相比常用的棱柱积分方法与Tesseroid方法，本发明不仅能够顾及地球曲率、全球地形质量影响及海水质量密度亏损，因此具有更小的地形改正的长波长误差，同时还解决了利用多分辨地形模型实现局部区域地形改正的高精度计算。

附图说明

图1为本发明实施例一种基于六边形网格索引下的多分辨率多面体计算全球地形改正的方法的流程图；

图2为本发明实施例7孔径球面六边形网格系统构建示意图；

图3为本发明实施例六边形网格多分辨率组合示意图(全球分辨率Level 2)；

图4为本发明实施例H3六边形网格与三角形剖分之间的对应关系(L＝0,1,2)；

图5为本发明实施例基于全球地形生成3套不同密度的多面体模型示意图；

图6为本发明实施例多面体积分公式变量示意图；

图7为本发明实施例多面体从表面积分到线积分转换中的奇异情况；

图8为本发明实施例嵩山区域实测重力点分布示例图；

图9为本发明实施例嵩山试验区地形多面体模型示意图(全球网格分辨率L＝3)；

图10为本发明实施例地形质量的引力效应与空间重力异常之间的关系。

具体实施方式

下面结合附图和具体的实施例对本发明做进一步的解释说明：

如图1所示，本发明的一种基于六边形网格索引下的多分辨率多面体计算全球地形改正的方法，首先，利用球面六边形网格系统对地球进行多分辨率剖分，得到不同分辨率的六边形网格，实现六边形网格中点坐标、边界点坐标以及个网格之间的位置关系索引；其次，依据地形改正计算中多分辨率地形的分区策略，根据区域的位置和区域大小，结合对地形改正计算精度的要求，确定六边形网格下的分区参数；然后，根据六边形网格与三角面之间的对应关系，确定多分辨率多面体地形模型，顾及海水质量密度亏损，利用多面体积分方法，确定待计算点处所受的地形改正影响。

具体步骤如下：

步骤1：全球六边形离散网格系统剖分。根据球面离散网格系统构建理论，生成一种基于二十面体的7孔剖分的球面六边形网格系统H3，利用该网格系统，方便快捷的确定网格分辨率、网格数量、网格中点坐标、网格边界坐标、网格相互关系等内容。

利用球面六边形网格构建理论生成全球7孔径六边形网格H3系统，如图2所示，该网格系统能够实现相邻分辨率中大的六边形网格近似包含7个小的六边形网格，大六边形和小六边形的方向存在19.1°的转角，该网格系统实现了不同分辨率六边形网格之间的快速索引。不同分辨率的网格参数如下表所示：

表1 H3网格的参数统计(球半径R＝6378.136km) 引自(https://h3geo.org/)

步骤2：确定六边形网格多分辨率地形分区策略。结合地形改正的理论和计算方法，参考目前使用最为广泛的地理网格下地形多分辨率分区策略，构建H3六边形网格系统下的地形分区策略。

计算顾及全球地形影响的地形改正过程中，由于随着地形分辨率的升高，计算量迅速增加，导致存在计算难题。为解决解算效率的难题，通常采用计算点周边区域高分辨率的地形数据，远区采用较低分辨率的地形数据，从而在保证精度较少损失的情况下，实现地形改正的高效计算。目前普遍使用的地理网格的地形多分辨率分区策略是将全球地形分为以下5个区域，区域大小及地形分辨率如下表所示：

表2地形改正的地理网格多分辨率分区策略

基于本发明使用的H3六边形网格系统，参考表1中对应的网格参数，确定如下六边形网格下全球多分辨率地形分区策略：

表3地形改正的六边形网格多分辨率分区策略

表3中的分区策略确定的地形模型能够接近于表2中的地形分辨率和各分辨率地形的范围。根据上述分区策略，构建了以珠穆朗玛峰区域计算点为中心的多分辨率多面体地形模型，为了增强展示性，此处所有六边形网格的分辨率均降低3层，即全球采用L＝2分辨率的六边形网格，其详细划分如图3所示。

步骤3：利用六边形网格与三角面之间的对应关系，确定多分辨率地形的多面体模型，顾及海水质量密度亏损，利用多面体积分方法，确定待计算点处所受的地形改正影响。

根据六边形网格剖分与三角形剖分之间的对偶关系，可以利用六边形网格对三角形进行检索和管理，重建六边形网格分布下的规则三角网地形，基于多面体Polyhedron积分方法求解该地形质量产生的引力作用。利用H3六边形网格索引系统，可生成全球分布的六边形网格，包括12个五边形，每个六边形网格的中心点与6个边界点可构成6个近似等边三角形，从而形成了六边形网格与三角网格之间的对应关系，如图4所示。

根据上述关系，利用多分辨率六边形网格剖分，可得到以三角面组成的多分辨率地形的多面体模型，从而方便利用多面体积分方法求解地形改正。当顾及全球地形质量以及海水质量密度亏损下的地形改正时，将陆地平均地壳密度设为2.67g/cm³，相比于该地壳密度值，海水密度为1.03g/cm³，因此存在-1.64g/cm³的补偿密度。因此，利用上述三角形面组成的点位数据，结合全球1′的地形数据ETOPO1以及30m的局部区域地形数据ASTER，可构建出3套不同密度的多面体模型，这3套不同密度的多面体所综合的影响等价于陆地区域2.67g/cm³以及海洋区域-1.64g/cm³产生的地形改正效应，如图5所示，其中“多面体模型1”是以大地水准面高加上全球地形模型确定的地形大地高为边界的、密度为2.67g/cm³的多面体，“多面体模型2”是陆地区域以大地水准面为外边界，海洋区域以实际地形为边界的、密度-1.03g/cm³的多面体，“多面体模型3”是以整个大地水准面为外边界的密度为-1.64g/cm³的多面体。其中大地水准面可利用EGM2008地球重力场模型计算得到，其精度完全满足地形改正的需要。

步骤4：利用前述构建的三组多面体模型及其对应海水密度，使用多面体(Polyhedron)积分方法，确定待计算点处的3组引力值，3组引力值之和，即为最终该待计算点处的地形改正值。

无论是对于数字高程模型(DEM)的规则网格高度值，还是对于描述某些埋藏质量异常的点云，通用多面体都是一种高效而灵活的建模方法，该方法是通过先将体积积分转换为表面积分，然后再转换为线积分来实现求解的，所以它能够应用于更为复杂形状的分布，其主要的特性是，不仅可用于形状复杂的多面体计算，还可以仅使用两个不同的线积分来表示所有引力场量，包括引力位、3个引力矢量分量、6个引力位二阶梯度张量。

如图6所示，整个多面体模型能够分解为n个三角面ABC分别对计算点Q的引力，然后累加求和即为整个均质多面体对Q点的引力值，利用三角面坐标和计算点Q的坐标，可确定图中l₁、l₂、h_p、h_pq以及与e₁′、e₂′、e₃′之间的夹角。

设多面体共由n个三角面ABC组成，那么Q点所受地形改正的引力影响表达式为

其中

其中d_As为奇异项改正值，表示第p个三角面的外法线方向单位矢量，/>是以Q点在第p个三角面上的垂足点Q′为原点建立的用于描述待求解引力矢量三分量坐标系的三个轴单位矢量，h_p是计算点Q到第p个三角面的距离，h_pq是计算点Q在第p个三角面上的垂足Q′到第p个三角面上第q条边的距离；σ_pq表示Q点的垂足点Q′在直线AB外法线一侧为1，另一侧为-1；G表示万有引力常数，ρ表示多面体质量的密度；式中S₁、S₂、l₁、l₂表示图6中AQ″、BQ″、AQ、BQ的距离，Q″表示Q′在第p个三角面第q条边上的垂足点。

多面体积分应用中，当计算点在平面上的投影点在区域中心、边界边和边界点上的时候，其不连续性使得散度定理无法应用，因此，当投影点分别在内部、边线上和边界点上的时候(图7)，分别需要加入奇异项改正。

对应图7的三种分布情况，(1)式中所加的奇异项改正d_As分别为

d_As1＝-2πh_p d_As2＝-πh_p d_As3＝-θh_p (3)

其中π为圆周率，θ为Q点在第p个三角面上的垂足Q′所在三角面的两条边的夹角大小。

本步骤中，使用多面体积分方法计算全球地形的地形改正，相比棱柱积分方法与Tesseroid方法，不仅三角面的多面体模型与真实地形的吻合度相比离散的地理网格要高，而且多面体积分计算的效率要普遍高于前面两种方法。

为验证本发明效果，本发明以计算我国嵩山地区152个高密度的加密重力测量点处的地形改正为例(图8)，采用上述多分辨率分区策略，利用该区域周边1″分辨率的ASTER地形数据以及全球ETOPO1地形数据，构建了适用于该区域的多面体模型，示意图如图9所示。

图10显示了基于L＝4和L＝5两个多面体模型计算的地形引力大小与实测重力资料计算得到的空间重力异常数据。两种方法计算结果之间的系统误差是由多方面原因造成的，包括计算地形影响所使用的参考球或者椭球半径、远区海洋质量亏损的影响、棱柱积分没考虑远区地形效应的影响等，而从其局部的变化来说，两者具有非常好的一致性。

综上，本发明的一种基于六边形网格索引下的多分辨率多面体计算全球地形改正的方法，利用多面体方法实现顾及全球高分辨率地形的地形改正计算，并利用六边形网格与三角形之间的对偶关系，实现多种分辨率的多面体构建，从而高效的实现局部区域地形改正的高精度计算。相比常用的棱柱积分方法与Tesseroid方法，本发明不仅能够顾及地球曲率、全球地形质量影响及海水质量密度亏损，因此具有更小的地形改正的长波长误差，同时还解决了利用多分辨地形模型实现局部区域地形改正的高精度计算。

以上所示仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。

Claims (4)

1.一种基于六边形网格索引下的多分辨率多面体计算全球地形改正的方法，其特征在于，包括：

步骤1：基于二十面体构建7孔剖分的球面六边形网格系统H3，利用该网格系统确定网格分辨率、网格数量、网格中点坐标、网格边界坐标、网格相互关系；

步骤2：基于六边形网格系统H3及该系统中各网格参数构建六边形网格系统H3下的地形分区策略；

步骤3：利用六边形网格与三角面之间的对应关系，确定多分辨率地形的多面体模型，顾及海水质量密度亏损，利用多面体积分方法，确定待计算点处所受的地形改正影响；

所述步骤3包括：

根据六边形网格剖分与三角形剖分之间的对偶关系，利用六边形网格对三角形进行检索和管理，重建六边形网格分布下的规则三角网地形，基于多面体积分方法求解该地形质量产生的引力作用；利用H3六边形网格索引系统，生成全球分布的六边形网格，每个六边形网格的中心点与6个边界点构成6个近似等边三角形，从而形成六边形网格与三角网格之间的对应关系；

根据上述关系，利用多分辨率六边形网格剖分，得到以三角面组成的多分辨率地形的多面体模型，以利用多面体积分方法求解地形改正：当求解顾及全球地形质量以及海水质量密度亏损下的地形改正时，将陆地平均地壳密度设为2.67g/cm³，海水密度为1.03g/cm³，利用上述三角形面组成的点位数据，结合全球1′的地形数据ETOPO1以及30m的局部区域地形数据ASTER，构建出3套不同海水密度的多面体模型，分别为，以大地水准面高加上全球地形模型确定的地形大地高为边界、密度为2.67g/cm³的多面体，陆地区域以大地水准面为外边界、海洋区域以实际地形为边界、密度为-1.03g/cm³的多面体，以整个大地水准面为外边界的密度为-1.64g/cm³的多面体；

步骤4：利用构建的多分辨率地形的多面体模型及其对应海水质量密度，使用多面体积分方法，确定待计算点处的各多面体模型对应的引力值，多组引力值之和即为最终该待计算点处的地形改正值。

2.根据权利要求1所述的一种基于六边形网格索引下的多分辨率多面体计算全球地形改正的方法，其特征在于，所述六边形网格系统H3相邻分辨率中大的六边形网格包含7个小的六边形网格，大六边形和小六边形的方向存在19.1°的转角，利用该网格系统可以进行不同分辨率六边形网格之间的快速索引。

3.根据权利要求1所述的一种基于六边形网格索引下的多分辨率多面体计算全球地形改正的方法，其特征在于，所述步骤2包括：

基于六边形网格系统H3及该系统中各网格参数、按照如下方式构建六边形网格系统H3下的地形分区策略：

区域1：区域半径为60km，网格分辨率为0.174km；

区域2：区域半径为158km，网格分辨率为0.46km；

区域3：区域半径为418km，网格分辨率为1.2km；

区域4：区域半径为1107km，网格分辨率为3.2km；

区域5：区域半径为全球，网格分辨率为8.5km。

4.根据权利要求1所述的一种基于六边形网格索引下的多分辨率多面体计算全球地形改正的方法，其特征在于，所述步骤4包括：

设多面体共由n个三角面ABC组成，则待计算点Q所受地形改正的引力影响表达式为：

其中

其中d_As为奇异项改正值，表示第p个三角面的外法线方向单位矢量，/>是以Q点在第p个三角面上的垂足点Q′为原点建立的用于描述待求解引力矢量三分量坐标系的三个轴单位矢量，h_p是计算点Q到第p个三角面的距离，h_pq是计算点Q在第p个三角面上的垂足Q′到第p个三角面上第q条边的距离；σ_pq表示Q点的垂足点Q′在直线AB外法线一侧为1，另一侧为-1；G表示万有引力常数，ρ表示多面体质量的密度；S₁、S₂、l₁、l₂分别表示AQ″、BQ″、AQ、BQ的距离，Q″表示Q′在第p个三角面第q条边上的垂足点。