一种确定地球重力场的六边形网格剖分方法

1.本发明属于物理大地测量技术领域，尤其涉及一种确定地球重力场的六边形网格剖分方法。


背景技术：

2.人类了解地球重力场，从传统的地面重力测量规划、数据统计、数据处理等，都是在地理网格分布方式下进行的。而随着重力场数据量和数据分辨率的提高，传统地理网格分布带来的局限性逐渐凸显，包括网格面积不相等、网格平滑因子复杂、高纬度区域数据冗余、网格角度分辨率低、网格不具有相邻一致性、球谐分析中混频效应显著、积分离散化误差大等。


技术实现要素：

3.本发明针对传统地理网格形状和面积随纬度发生变化的问题，提出一种确定地球重力场的六边形网格剖分方法，不仅有效解决了传统地理网格形状和面积随纬度发生变化的问题，还有效提高了地球重力场模型求解精度，为实现重力测量高效规划、重力数据高效使用、高精度处理和重力场模型构建提供了重要的理论支撑和借鉴意义。
4.为了实现上述目的，本发明采用以下技术方案：
5.本发明提出一种确定地球重力场的六边形网格剖分方法，包括：
6.步骤1，基于二十面体及synder等积投影构建4孔剖分的球面六边形网格系统isea4h，通过isea4h生成不同分辨率下的全球均匀分布的六边形网格，并确定网格分辨率、网格数量、网格中点坐标、网格边界坐标；
7.步骤2，利用现有重力观测资料，基于isea4h生成的六边形网格剖分区域，进行数据统计并生成六边形网格重力异常数据；
8.步骤3，利用不同分辨率下的全球均匀分布的六边形网格重力异常值，根据球谐展开理论公式计算球谐系数；
9.步骤4，利用不同分辨率下的全球均匀分布的六边形网格重力异常值，利用重力场stokes定理，基于离散积分求和得到局部大地水准面。
10.进一步地，所述基于二十面体及synder等积投影构建4孔剖分的球面六边形网格系统isea4h包括：
11.选用二十面体作为基础面，其与地球球体之间的定位关系满足：二十面体的 12个顶点中的两个顶点位于地球南北两极，且还有一个顶点位于起始子午面上；将二十面体展开到平面上，在平面上进行六边形网格的剖分，选用4孔径剖分，生成不同分辨率的多面体表面的六边形网格；通过synder投影建立平面与球面之间的位置对应关系；完成基于二十面体及synder等积投影的4孔六边形网格系统isea4h。
12.进一步地，所述步骤2包括：
13.以isea4h生成的六边形网格边界坐标作为条件，判断离散点观测资料是否位于网
格范围内，统计所有位于网格内部的重力资料，并融合生成代表网格所在区域的重力异常值；
14.利用现有多源重力场资料，包括卫星重力数据、卫星测高数据、地面重力数据、海洋重力数据、航空重力数据，融合多种数据生成全球不同分辨率的六边形网格重力异常数据。
15.进一步地，所述步骤3包括：
16.根据重力场理论，得到球面上重力异常点值与n阶重力场模型扰动位系数之间的数学关系，即球谐分析的数学模型：
[0017][0018]
其中n为重力场模型的最大阶数，n≤n,m≤n，δg是球坐标为(r,θ,λ)的点的空间重力异常，θ表示地心余纬，λ表示地心经度，r表示点至地心的距离， gm为地心引力常数，a是参考半径，表示球谐展开的扰动位系数，为归一化缔合legendre函数。
[0019]
进一步地，所述步骤4包括：
[0020]
在球近似下，基于全球的stokes积分理论实现大地水准面高的求解，其表达式如下：
[0021][0022]
其中，n表示大地水准面高，r是地球平均半径，γ是参考椭球表面的正常重力，(ψ,α)分别指球面角距和方位角，s(ψ)为球面stokes核函数，其表达式为：
[0023][0024]
采用离散化求和代替积分计算，求解大地水准面高。
[0025]
与现有技术相比，本发明具有的有益效果：
[0026]
本发明的一种确定地球重力场的六边形网格剖分方法，相比常用的地理网格剖分，该六边形网格具有全球均匀分布特性以及等积特性，基于该六边形网格剖分下的重力数据在实际工程应用和数据管理、统计中具有优势，能够减小代表误差，其统计重力异常值能够显著代表网格内的重力场情况，利用该六边形网格重力异常值进行球谐分析，不仅大大减少了必须观测数据的数据量，还能够减少混频效应，提升球谐分析的精度，同时相比地理网格，利用该六边形网格重力异常值求解局部大地水准面的离散化误差和中央网格奇异值补偿精度也较高。
附图说明
[0027]
图1为本发明实施例一种确定地球重力场的六边形网格剖分方法的流程图；
[0028]
图2为本发明实施例isea4h网格生成过程示意图；
[0029]
图3为本发明实施例isea4h网格的分布规律和管理规则；
[0030]
图4为本发明实施例四边形格网和六边形网格的积分离散化示意图；
[0031]
图5为本发明实施例全矩阵最小二乘方法恢复的360阶球谐系数的阶误差图；
[0032]
图6为本发明实施例两种网格由离散化采样不足造成的混频效应比较；
[0033]
图7为本发明实施例两种网格由模型截断造成的混频效应比较。
具体实施方式
[0034]
下面结合附图和具体的实施例对本发明做进一步的解释说明：
[0035]
如图1所示，本发明的一种确定地球重力场的六边形网格剖分方法，首先，确定一种球面均匀分布的等面积网格，基于实测重力数据，统计网格平均重力异常，然后基于该网格分布下的重力异常进行球谐分析，确定地球重力场模型位系数，基于该模型系数可实现任意点位处的重力场元快速球谐综合计算。
[0036]
具体步骤如下：
[0037]
步骤1：全球等积六边形网格生成。根据球面离散网格系统构建理论，确定生成该网格系统的五个独立要素(基础多面体、球与多面体的定位、平面剖分方法、投影方法以及点与网格对应关系)，生成一种基于二十面体的、synder等积投影的4孔剖分的球面六边形网格系统isea4h，确定网格分辨率、网格数量、网格中点坐标、网格边界坐标等内容。
[0038]
由于二十面体是五个理想多面体中面数最多的，因此该多面体与球的吻合度也是最好的，从而选用二十面体作为基础面，其与地球球体之间的定位关系满足：二十面体的12个顶点中的两个顶点位于地球南北两极，且还有一个顶点位于起始子午面上，如图2中(a)所示；二十面体可展开到平面上(图2中(b))；在平面上可进行六边形网格的剖分，选用4孔径剖分，剖分的次数越多，生成的六边形网格越多，分辨率越高，从而生成不同分辨率的多面体表面的六边形网格(图 2中(c))；地图投影建立了平面与球面之间的位置对应关系，这一过程是整个球面六边形网格构建的关键环节，采用美国著名的synder投影可实现基于多面体的等积投影，且在保证经纬网连续的同时减小了变形(图2中(d))。经过上述步骤，生成了一种是基于二十面体synder等积投影的4孔六边形网格系统 (icosahedral snyder equal area aperture 4hexagon，isea4h)，如图2中(e)所示。
[0039]
利用上述方法生成的isea4h网格，无论剖分多细，整个系统最终只包含12 个五边形，且其中心位于二十面体的12个顶点处，其余全部为球面六边形。不同分辨率的网格参数如下表所示：
[0040]
表1 isea4h网格的参数统计(球半径r＝6378.136km)
[0041]
level l(分辨率)网格数量网格面积(km2)网格间距(km)0124.260
×
1073.682
×
1031421.217
×
1071.968
×
10321623.156
×
1061.002
×
10336427.963
×
1055.035
×
102425621.995
×
1052.520
×
1025102424.991
×
1041.260
×
1026409621.248
×
1046.303
×
1071638423.120
×
1033.152
×
1086553627.800
×
1021.576
×
10
926214421.950
×
1027.87910104857624.875
×
103.939
[0042]
上述方法构建的六边形网格能够很好的保证所有网格面积近似相等(除12 个五边形外)，更重要的优点是，isea4h网格系统保证了每个网格之间的位置关系简单，可用5个矩阵表示，如图3所示，这有助于使用计算机对网格数据进行管理和检索。
[0043]
步骤2：生成网格重力异常数据。利用现有重力观测资料，基于isea4h生成的六边形网格剖分区域，进行数据统计并生成网格重力异常数据。
[0044]
以isea4h生成的六边形网格边界坐标作为条件，判断离散点观测资料是否位于网格范围内，统计所有位于网格内部的重力资料，并融合生成代表网格所在区域的重力异常值。
[0045]
实践表明，对于同一套实测资料，基于六边形网格的有效网格占比(即包含重力实测点的网格占总网格的比例)比四边形大，六边形网格空间重力异常代表误差相比等面积的四边形网格空间重力异常代表误差小，充分说明了六边形网格在实际工程应用和数据管理、统计中的优势。
[0046]
步骤3：利用全球均匀分布的六边形网格重力异常值，可根据球谐展开理论公式，计算球谐系数，相比传统地理网格下的球谐分析，六边形网格下的球谐分析不仅能够减小必须观测量，还能够有效的减少球谐分析中的混频效应，提升计算精度。
[0047]
球谐分析方法通常使用两种方法：数值积分法与最小二乘方法，数值积分方法因其计算量小及计算近似等特点，早期应用较多，随着计算性能的提高以及多种类型数据联合构建模型的需要，最小二乘方法得到了普遍应用。根据重力场理论，可以得到球面上重力异常点值与n阶重力场模型扰动位系数之间的数学关系，即球谐分析的数学模型(球谐展开理论公式)：
[0048][0049]
其中n为重力场模型的最大阶数，δg是球坐标为(r,θ,λ)的点的空间(或自由)
[0050]
重力异常，θ表示地心余纬，λ表示地心经度，r表示点至地心的距离，gm为地心引力常数，a是参考半径，(gm,a)由位模型所采用的地球参考椭球决定(egm2008模型使用的a为6378137m)，表示球谐展开的扰动位系数，为归一化缔合legendre函数。将上述数学模型表示为以下简写形式
[0051]
l＝f(x)
ꢀꢀ
(8)
[0052]
其中
[0053][0054]
n表示所恢复模型的阶数，num表示球面分布的离散点重力异常总个数，f 函数是关于x的线性方程，根据线性最小二乘平差模型，有
[0055][0056]
v表示观测量残差，a为设计矩阵，设p为重力异常权矩阵，得到法方程
[0057]
[0058]
通常p矩阵取对角阵，最终可得位系数模型的解及相应的协方差阵
[0059][0060][0061]
其中a
t
pa为法矩阵，使用n来表示。法矩阵n的条件数决定了该方程结构的好坏，同时n的稀疏特性(例如块对角结构)决定着超高阶次球谐模型最小二乘求解的快速实现。
[0062]
上述球谐分析解算过程中，当计算n
max
阶次的球谐位系数，根据nyquist采样定律，要求网格最小分辨率为
[0063][0064]
所以，求解n
max
阶位系数模型，至少需要个地理网格重力异常数据，而基于六边形网格计算的球谐位系数，只需要六边形网格重力异常数据大于位系数未知数的个数，即(n
max
1)2即可，大大节约了必须观测值的数量。
[0065]
球谐分析过程中，会出现混频现象。混频效应是对连续信号进行离散化采样过程中，采样数量不足以完全恢复原有信号从而导致高频率的信号被混入低频中而发生混淆，这是信号分析中的概念。一方面是离散化采样不足产生的，当采样密度无穷大的时候，该部分混频效应将会消失；另一方面是因为球谐展开的阶次截断造成的，这是由最小二乘法矩阵的结构决定的。
[0066]
本步骤中，使用六边形网格下的重力异常资料计算球谐位系数与传统地理网格计算球谐系数相比，大大改善了观测点的分布结构，使得法矩阵的结构稳定，不仅节约了必须观测值的数量，还削弱了球谐分析过程中的混频效应，可有效提升重力场解算精度。
[0067]
步骤4：利用球面六边形网格重力异常，利用重力场stokes定理，可基于离散积分求和得到局部大地水准面高数值模型。相比传统地理网格下的球谐分析，六边形网格顾及了stokes核函数的各项同性特性，具有更小的离散化误差和较高的中央奇异性补偿精度，可大大提升大地水准面求解精度。
[0068]
在球近似下，基于全球的stokes积分理论可实现大地水准面高的求解，其表达式如下
[0069][0070]
其中，n表示大地水准面高，r是地球平均半径，γ是参考椭球表面的正常重力，(ψ,α)分别指球面角距和方位角，s(ψ)称为球面stokes核函数，其表达式为
[0071][0072]
实际应用中采用离散化求和代替积分计算，求解大地水准面高。而在分辨率和计算量相同的情况下，六边形网格相比四边形网格，具有更高的计算精度。
[0073]
本步骤中，使用六边形网格下的重力异常资料，采用离散求和计算大地水准面，相比等分辨率、等数据量、等计算量的地理网格的重力异常，离散化误差小，计算精度更高。
[0074]
为验证本发明效果，进行如下实验：
[0075]
(1)六边形网格在重力资料重力异常统计中的优势评估
[0076]
利用表2中的分辨率和数量近似的四边形地理网格和六边形网格对我国重力资料(约100多万点数)的代表误差等统计参数。
[0077]
表2我国重力资料网格划分及对比
[0078][0079]
代表误差指面积内任意一点的重力异常代表该面积平均重力异常所产生的中误差，以表示网格平均重力异常的代表误差，根据统计理论，有
[0080][0081]
其中，n代表网格内的实测重力点数量，δgi为第i个空间重力异常，表示指定面积内所有空间重力异常的平均值，即网格平均重力异常。
[0082]
表3具有实测点的网格统计
[0083][0084]
表4代表误差的统计情况
[0085][0086]
可以看出，由于形状更加接近于圆而结构更优，六边形网格在平均面积近似下，相比四边形网格更加聚中，网格内重力实测点之间的距离更近，因此基于六边形网格的代表误差小于四边形网格，加之六边形网格有效网格的占比高，说明了六边形网格的在实际数据统计方面的优势。
[0087]
(2)六边形网格在球谐分析中的优势评估
[0088]
利用egm2008模型前359阶次系数计算了7层isea共163842个网格重力异常点值和传统地理网格分布下360
×
720即259200个网格重力异常点值，利用全矩阵最小二乘方法恢复同样阶次的系数，可得所恢复系数的阶rms图，如图 5所示。
[0089]
由图可知，传统地理网格分布至少需要360
×
720即259200个重力异常点值方能保
证恢复359阶以下位系数的矩阵不秩亏，而基于isea4h六边形网格分布，只需要163842个点值就能稳定而高精度的恢复359阶次的球谐系数，节省了约 37％的观测量，这是六边形网格相比地理网格的优势。
[0090]
为了说明由采样离散化造成的混频效应，取模型egm2008最高n
max
＝365 阶次计算模型重力异常，利用最小二乘全矩阵恢复相同365阶次的球谐系数，结果与原球谐系数真值相比存在差异，我们认为该误差是由于离散化采样不足造成的，而与模型的截断没有关系，结果见图6。
[0091]
图6说明，相比地理网格，isea4h网格在满足必须观测量的情况下，更能够保证球谐系数的恢复而不失真，而传统奈奎斯特采样准则已经不再适用于该种分布情况，地理网格由于采样频率l＝360，小于模型恢复的最大阶次n
max
＝365，因此其恢复精度太差，可认为是混频效应的结果。
[0092]
为了说明由模型截断造成的混频效应，需要采样率远远高于奈奎斯特准则，即离散化误差足够小，方可认为模型截断是造成混频效应的唯一因素。根据奈奎斯特采样准侧，全球259200个地理网格对应分辨率为30
′
，取egm2008模型前 92阶次系数计算模型重力异常，基于该重力异常，利用全矩阵最小二乘方法恢复阶次n
max
＝90的球谐系数，由于采样点数相对于90阶次来说相对较密，所以认为gsha误差主要为模型截断造成的混频效应，结果见图7。
[0093]
图7中结果显示，两种网格分布下的球谐系数恢复误差均较大，实验过程虽然不符合常规习惯，但是该计算过程能够很好地反映由于模型的截断造成的混频效应大小，结果显示，在isea4h六边形网格中，91～92阶次的信号污染了前90 阶次的系数且量级平均在10-12
，而四边形地理网格中，信号污染主要集中在低次和高阶系数部分且量级平均在10-11
。可见，相比四边形，六边形网格限制高频信号混频入低频信号的能力更强。
[0094]
(3)六边形网格在重力场数值积分中的优势评估
[0095]
根据球面stokes计算公式(15)，如果被积分的重力异常δg为常数时，其可将该常数移至积分号外，得到
[0096][0097]
其中
[0098][0099]
积分区域从球冠半径5
°
到180
°
，利用上式计算值作为真值，统计六边形网格和四边形网格的离散化误差，由于不同分辨率造成中央网格不参与计算引入的系统差异，因此上述差异的平均值不具有统计意义，这里只统计各个积分半径造成误差的标准差std，结果如表5所示。
[0100]
表5 stokes核函数积分离散化误差统计
[0101][0102]
结果可见，在基于四边形网格和六边形网格进行离散积分运算中，相比四边形网格，六边形网格具有更小的积分离散化误差。
[0103]
综上，本发明的一种确定地球重力场的六边形网格剖分方法，相比常用的地理网格剖分，该六边形网格具有全球均匀分布特性以及等积特性，基于该六边形网格剖分下的重力数据在实际工程应用和数据管理、统计中具有优势，能够减小代表误差，其统计重力异常值能够显著代表网格内的重力场情况，利用该六边形网格重力异常值进行球谐分析，不仅大大减少了必须观测数据的数据量，还能够减少混频效应，提升球谐分析的精度，同时相比地理网格，利用该六边形网格重力异常值求解局部大地水准面的离散化误差和中央网格奇异值补偿精度也较高。
[0104]
以上所示仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。