本文内容来源于《测绘学报》2024年第8期（审图号GS京(2024)1748号）
基于ico_HEALPix网格的超高阶地球重力场建模方法
张展鹏
,1, 李新星,1,2,3, 刘长建1, 范昊鹏1, 裴宪勇1
1.信息工程大学地理空间信息学院，河南 郑州 450001
2.地理信息工程国家重点实验室，陕西 西安 710054
3.智慧地球重点实验室，北京 100020
摘要：
国家自然科学基金(U23A2028)(41404020)(42174007)(42074013); 地理信息工程国家重点实验室基金(SKLGIE2023-Z-1-1)
作者简介：
张展鹏（2000—），男，硕士生，研究方向为测量数据处理理论与方法。E-mail：[email protected]通讯作者: 李新星 E-mail：[email protected]
摘要：
本文针对传统地理网格数据剖分的重力场模型构建在高纬度区域出现数据冗余的问题，将分层等积等纬度像素化（HEALPix）网格结构引入地球重力场解算中，提出了利用二十面体HEALPix（ico_HEALPix）网格的超高阶地球重力场建模方法，实现了全球3600阶次球谐位系数的高效构建，同时针对ico_HEALPix网格在球谐分析过程中法矩阵不是严格块对角化结构的问题，设计了迭代算法，有效提高模型构建的精度。试验表明，ico_HEALPix网格数据在数据量小于地理网格500万的前提下，通过迭代方法构建的全球地球重力场模型精度可达到优于地理网格的效果，球谐位系数误差阶RMS提升1~2个数量级，还解决了地理网格南北极点畸变和数据冗余的问题，提高了网格的数据利用率。
关键词
HEALPix网格; ico_HEALPix网格; 超高阶地球重力场模型; 迭代算法; 最小二乘法; XGM2019e
本文引用格式
张展鹏, 李新星, 刘长建, 范昊鹏, 裴宪勇. 基于ico_HEALPix网格的超高阶地球重力场建模方法[J]. 测绘学报, 2024, 53(8): 1531-1539 doi:10.11947/j.AGCS.2024.20230151
ZHANG Zhanpeng, LI Xinxing, LIU Changjian, FAN Haopeng, PEI Xianyong. Ultra-high-order Earth gravity field modeling method based on ico_HEALPix grid[J]. Acta Geodaetica et Cartographica Sinica, 2024, 53(8): 1531-1539 doi:10.11947/j.AGCS.2024.20230151
阅读全文
http://xb.chinasmp.com/article/2024/1001-1595/1001-1595-2024-08-1531.shtml
基于传统等角地理网格构建的地球重力场模型面临着许多无法克服的局限性，包括网格面积不等、观测数据冗余、重力场平滑因子复杂、积分离散化误差大等，以此建立和管理全球多分辨率数据时，在两极区域会产生较大的形变，并且高纬度区域具有大的数据冗余，相应导致球谐分析过程的解算效能较差[1-5]。针对球谐分析解算中传统地理网格所面临的问题，已有多名学者采用新型网格剖分方式进行相关解算并取得较好结果。文献[6]求解地球重力场模型的法矩阵结构的好坏受到了球面网格分布的影响，球面均匀分布的网格结构对于球谐分析和综合问题非常重要；文献[7]采用快速傅里叶变换技术(fast fourier transform，FFT)技术验证了基于八面体和二十面体的近似全球均匀分布点位相比地理网格，在球谐分析和综合方面解算效率和精度明显提高。天文学中，分析天球的宇宙微波背景辐射(cosmic microwave background，CMB)数据可以更好地理解宇宙的几何形状和组成[8]，同时为解决天球面上离散数据的可视化和数据快速分析的需求，文献[9-10]设计了分层等积等纬度像素化(hierarchical equal area isolatitude pixelation，HEALPix)网格结构，因该网格优良的特性及其同名软件的发布[11]，使得HEALPix剖分方案在天体物理学[12]、天文学[13-14]、地球物理学[15]、行星科学[16]、核工程[17]和计算机视觉[18]等应用中得到推广。文献[19]利用非均匀FFT、双傅里叶球面法和Slevinsky快速球面谐波变换的新组合对基于HEALPix这一等纬度等面积的网格类型进行球谐分析，结果表明提升球谐分析精度的同时，收敛速度至少高出两倍以上。针对地理网格重力数据构建地球重力场模型时存在的数据分布不均匀、高纬区域数据冗余等问题，本文将HEALPix网格剖分方案引入球谐分析，证实常规HEALPix网格在球谐分析中因高纬度区域数据经度方向采样不足造成较大误差，相应提出改进的ico_HEALPix(icosahedral HEALPix)网格，利用块对角最小二乘(least-squares，LS)快速计算方法[20-22]，针对法矩阵非完全块对角化问题，设计迭代算法，实现3600阶次的地球重力场模型构建。
1 建模方法
1.1 HEALPix网格及其推广
HEALPix网格系统是一种多层次复合网格模型，在两极地区采用球面三角投影，非两极地区采用圆柱等面积投影。具有3个核心特征[23]:①数据按照分辨率分层构建；②球面的网格面积相等；③网格在球面上等纬度分布。该模型的空间网格和平面网格如图1所示。
图1
图1 空间网格[19]和平面网格
Fig.1 Spatial grid[19]and planar grid
HEALPix网格系统包括12块基本网格，从北到南分为3行，每行均匀分布4个基本网格，每个基本网格可以将边长细分为Nside=2k份，其中k=1,2,3,…,n表示网格剖分的层次，所以每一个基本网格划分为
块等面积的子网格。ico_HEALPix(icosahedral HEALPix)网格是本文对常规HEALPix网格的改化，将HEALPix网格3×4的12个基础网格扩充到4×5共20个基本网格，其平面展开如图2(b)所示。从北到南共4行，每行均匀分布5个基本网格。为便于重力场解算，将4×5基本网格的南半球整体平移π/5，可得到关于赤道南北对称的网格，平移后的网格平面分布如图2(c)所示。ico_HEALPix保留了HEALPix网格系统依分辨率分层构建、网格等面积等纬度分布等优良特性。
图2
图2 HEALPix网格及其改进网格平面展开
Fig.2 HEALPix grid and its improved grid plane expansion
如图3所示，表示ico_HEALPix网格在k=1时所对应的平面展开图和全球分布图。其中，边线代表了所剖分网格的边界，实体点表示每一小块地理网格所对应的中心点的位置。
图3
图3 ico_HEALPix网格平面展开和全球分布
Fig.3 Plane expansion and global distribution of ico_HEALPix grid
式(1)给出了从北到南不同纬度圈i所对应的经度分辨率大小λi
(1)
ico_HEALPix网格共有等面积但形状不同的子网格
块，这些网格的中心点分别位于5×Nside-1个纬度圈上，单位球面上每一网格像素大小均为
，每个像素可利用网格纬度圈索引值(Nside,i,j)计算其在球面上的坐标(θ,λ)，其中i表示从北到南数第i层纬度圈，j表示对应纬度圈所在的列号，θ为球面余纬，在北极为0，在南极为π,λ为经度，范围[0,2π]。当1≤i<Nside时定义为北极区域，当4Nside<i<5Nside时定义为南极区域，北极区域与南极区域之间定义为赤道区域，此时Nside≤i≤4Nside。设z=cos θ，网格中心点坐标计算公式如下。北极区域，索引值(i,j)满足1≤i<Nside,1≤j≤5i
(2)
(3)
(4)
式中，s表示在纬度圈索引体系中，索引值j相对纬度为0处的偏移量。
赤道区域，索引值(i,j)满足Nside≤i≤4Nside,1≤j≤5Nside
(5)
(6)
(7)
南极区域，索引值(i,j)满足4Nside<i<5Nside,1≤j≤25Nside-5i
(8)
(9)
(10)
通过上述公式，可以由当前索引值(i,j)计算网格中心点在球面上的球坐标(θ,λ),s取值为0或1。
为充分说明ico_HEALPix的优势，在保证地理网格和ico_HEALPix网格数大致相等的条件下，对比两网格的结构形状。选取分辨率18°×18°的地理网格和剖分层次k=1的ico_HEALPix网格，两网格的点数对比见表1。
表1 地理网格和ico_HEALPix网格数量对比
Tab.1 Comparison of the number of geographic grids and ico_HEALPix grids
| 网格 | 使用模型阶次 | 网格数量 |
|---|---|---|
| 地理网格 | 7 | 98 |
| ico_HEALPix网格 | 7 | 80 |
新窗口打开| 下载CSV
图4为地理网格和改进后的HEALPix网格在描述北极地区地球重力场的分布。
图4
图4 地理网格和ico_HEALPix网格北极地区结构对比
Fig.4 Comparison of Arctic structure between geographic gird and ico_HEALPix grid
如图4(a)所示，基于经纬度剖分的地理网格，从两极到赤道，网格的面积不相等，不同纬度圈的网格面积差异巨大，甚至不在同一量级，存在较大的畸变问题；如图4(b)所示，ico_HEALPix网格等面积等纬度分布，南北极地区的面积畸变问题得到了有效的解决，网格均匀无缝分布。为说明ico_HEALPix网格在赤道地区相比地理网格的优势，如图5所示，分别截取17°E到110°E,30°N到30°S的网格范围进行对比分析。
图5
图5 地理网格和ico_HEALPix网格赤道地区对比
Fig.5 Comparison of geographic grid and ico_HEALPix grid equatorial region
由图5可知，地理网格赤道地区网格形状过大，等面积剖分的ico_HEALPix网格使得赤道地区网格面积小于地理网格，赤道地区分辨率得到了提高。整体结构形状对比如图6所示。
图6
图6 地理网格和ico_HEALPix网格整体结构对比
Fig.6 Comparison of the overall structure between geographic grid and ico_HEALPix gird
在网格数量近似相等的条件下，相比于地理网格，ico_HEALPix网格有以下3方面优势：①ico_HEALPix网格的等面积剖分解决了地理网格极区网格过密，数据利用率低的问题；②ico_HEALPix网格的等面积剖分使得网格在赤道地区数据分辨率优于地理网格，网格结构构造更为合理；③ico_HEALPix网格近似球面均匀分布，各向同性相比地理网格更强。
1.2 重力异常和位系数的关系
重力异常的球谐展开公式为[24-25]
(11)
式中，Δg为空间重力异常观测量；(r,θ,λ)表示重力异常点的球坐标；
为所求的一组重力场模型扰动位系数，即待求量；a为正常椭球的长半径；
(cos θ)为完全归一化缔合Legendre函数；Deg是球谐模型的最高阶次。式(11)是重力异常值Δg与位系数
之间联系的桥梁，也是通过最小二乘法求解位系数的基础。将观测式(11)表示为最小二乘通用表达式
(12)
式中，
(13)
F函数是关于球谐系数Xa的线性方程，根据最小二乘平差模型，有法方程[26]
(14)
式中，N=ATPA称为法方程的系数矩阵，简称为法矩阵；U=ATPLb为法方程自由项向量。矩阵A的各个元素计算式为
(15)
(16)
式中，(θi,λj)为重力异常网格数据索引值(i,j)对应的球心余纬和球心经度。其解为
(17)
上述求解过程需要满足重力异常Δg的观测值个数大于待求量球谐系数
的个数，即对于ico_HEALPix和HEALPix网格，只需要满足总网格点数不少于待求位系数个数，而对于地理网格分布来说dθ=dλ，则需要满足Nyquist采样定律，即Deg≤π/dθ，采用最小二乘方法时不能取等号，即Deg<π/dθ[21]。表2给出了不同分辨率ico_HEALPix网格可恢复的模型阶次及对应阶次模型所需要的最少地理网格点数。
表2 网格之间所需点数对应关系对比
Tab.2 Comparison of the required number of points between grids
| Deg | k | Nside | 位系数个数 | ico_HEALPix点数 | HEALPix点数 | 地理网格点数 |
|---|---|---|---|---|---|---|
| 180 | 6 | 64 | 32 757 | 81 920 | 49 152 | 64 800 |
| 360 | 7 | 128 | 130 317 | 327 680 | 196 608 | 259 200 |
| 720 | 8 | 256 | 519 837 | 1 310 720 | 786 432 | 1 036 800 |
| 1440 | 9 | 512 | 2 076 477 | 5 242 880 | 3 145 728 | 4 147 200 |
| 1800 | 9 | 512 | 3 243 597 | 5 242 880 | — | 6 480 000 |
| 1800 | 10 | 1024 | 3 243 597 | — | 12 582 912 | 6 480 000 |
| 2160 | 10 | 1024 | 4 669 917 | 20 971 520 | 12 582 912 | 9 331 200 |
| 3600 | 10 | 1024 | 12 967 197 | 20 971 520 | — | 25 920 000 |
| 3600 | 11 | 2048 | 12 967 197 | — | 50 331 648 | 25 920 000 |
| 5400 | 11 | 2048 | 29 170 797 | 83 886 080 | 50 331 648 | 58 320 000 |
注：表2中“—”表示网格点数少于位系数个数，不能够用于相应阶次模型球谐系数的求解。
新窗口打开| 下载CSV
为了说明本文所提ico_HEALPix网格在计算球谐系数方面的优势，选择k=11的HEALPix网格数据和k=10的ico_HEALPix网格数据分别恢复5400阶和3600阶模型位系数。
1.3 块对角最小二乘迭代方法
XGM2019e是慕尼黑技术大学于2019年发布的一组5540阶次组合全球重力场模型，对应空间分辨率为2'[27]。本文采用XGM2019e前6阶的模型系数计算k=1的ico_HEALPix网格中点的重力异常值，然后利用最小二乘方法求解球谐系数，其法矩阵N系数如图7所示，其中白色部分代表系数值接近零，黑色部分代表非零值，可知，基于ico_HEALPix网格分布构建的6阶地球重力场模型的法矩阵为非完全块对角矩阵，块对角元素之外存在非零值，因此对于大型的法逆矩阵N-1的求解是困难的。
图7
图7 k=1的ico_HEALPix网格法矩阵
Fig.7 The ico_HEALPix grid method matrix of k=1
因为N为主对角占优的对称正定矩阵，对法矩阵N中的“块”逐个求逆，忽略非块对角矩阵参数，则可以得到N-1的一种逼近N*，利用该矩阵，可以得到循环迭代求解公式
(18)
式中，r(k+1)是第k+1次迭代的残余矢量。将式(18)结合得到新的表达式
(19)
该形式满足松弛因子为1的Richardson迭代形式，当矩阵(I-N*N)的谱半径小于1，则上述迭代方法收敛，最终迭代的结果一致逼近于真值。
针对该问题，本文设计了迭代算法，具体迭代流程如图8所示。
图8
图8 块对角最小二乘迭代流程
Fig.8 Block diagonal least squares iterative flow chart
2 数值模拟试验及分析
2.1 HEALPix网格与传统地理网格的比较
本文提出ico_HEALPix网格而摒弃使用HEALPix网格，是因为利用HEALPix分布方式进行重力异常采样并恢复地球重力场模型时存在一定的局限性。为了说明该一点，以k=6的HEALPix网格分布作为采样点，利用采样点处XGM2019模型180阶次系数仿真的模型重力异常值，采用块对角最小二乘方法恢复180阶次模型位系数，再正算采样点处的模型点重力异常，与原始仿真值相比，差异分布如图9所示。
图9
图9 HEALPix网格恢复球谐系数造成的差异
Fig.9 The difference caused by HEALPix grid restoring spherical harmonic coefficient
由图9可知，在中低纬度区域，所求模型的重力异常与原始仿真值高度吻合；而在高纬度区域，所求模型的重力异常与原始仿真值之间存在数个mGal的差异。进一步分析认为，这一差异源于HEALPix网格在高纬度区域经度方向的采样过于稀疏。为充分说明HEALPix网格存在的问题，选取XGM2019e前5400阶的模型系数作为输入数据，仿真分辨率2'×2'的地理网格和剖分层次k=11的HEALPix网格中心点的重力异常值，所需的网格点数见表3。
表3 地理网格和HEALPix网格数量对比
Tab.3 Comparison of the number of geographic grids and HEALPix grids
| 网格 | 阶次 | 网格数量 |
|---|---|---|
| 地理网格 | 5400 | 58 320 000 |
| HEALPix网格 | 5400 | 50 331 648 |
新窗口打开| 下载CSV
仿真得到地理网格和HEALPix网格中心点重力异常数据，然后采用块对角最小二乘迭代方法反算模型系数，与XGM2019e模型真值进行比较，计算所求模型位系数的阶误差RMS，其计算公式为[28]
(20)
利用图8所示的流程对HEALPix网格构建的超高阶次地球重力场模型进行迭代运算，统计未迭代、迭代4次、迭代9次的阶误差RMS情况，并与地理网格求解的模型精度结果进行比较，对比如图10所示。因为地理网格法矩阵严格满足块对角形式，所以不存在迭代的情况。
图10
图10 误差阶RMS比较
Fig.10 Error order RMS comparison
HEALPix网格未迭代、HEALPix网格迭代9次和地理网格的误差谱图如图11所示。
图11
图11 模型球谐系数误差谱
Fig.11 Model spherical harmonic coefficient error spectrum
由图10、图11知，经过迭代求解，基于HEALPix网格构建的地球重力场模型中低阶次精度得到了一定的提高，但模型高阶次部分精度恢复较差，分析原因在于HEALPix网格在高纬度区域经度方向采样率低。
2.2 ico_HEALPix与传统地理网格的比较
ico_HEALPix网格在经度方向的采样率是HEALPix的1.25倍，能够较好地改善HEALPix网格在高纬度区域经度方向采样过低的问题，并采用仿真试验进行验证。选取ico_HEALPix网格点数少于地理网格的XGM2019e前3600阶的模型系数作为输入数据，仿真分辨率3'×3'的地理网格和剖分层次k=10的ico_HEALPix网格中点处的重力异常值，以此作为观测数据，采用块对角最小二乘迭代方法恢复前3600阶次的模型位系数。使用XGM2019e作为标准模型，地理网格和ico_HEALPix网格分别恢复3600阶次的模型位系数，所需的网格数见表4。
表4 ico_HEALPix网格和地理网格数量对比
Tab.4 Comparison of the number of geographic grids and ico_HEALPix grids
| 网格 | 阶次 | 网格数量 |
|---|---|---|
| 地理网格 | 3600 | 25 920 000 |
| ico_HEALPix网格 | 3600 | 20 971 520 |
新窗口打开| 下载CSV
利用图8所示的流程对ico_HEALPix网格构建的超高阶次地球重力场模型进行迭代运算，统计未迭代近似解、迭代4次、迭代9次的阶误差RMS情况，并与地理网格求解的模型精度结果进行比较，对比如图12所示。
图12
图12 误差阶RMS比较
Fig.12 Error order RMS comparison
ico_HEALPix网格未迭代、ico_HEALPix网格迭代9次和地理网格的误差谱图如图13所示。
图13
图13 模型球谐系数误差谱
Fig.13 Model spherical harmonic coefficient error spectrum
由图12、图13可知，ico_HEALPix网格数据在数据量小于地理网格约500万的前提下，通过迭代方法构建的全球地球重力场模型精度可达到优于地理网格的效果，系数误差阶RMS提升1~2个数量级。
3 结论
针对地理网格剖分在重力场球谐分析过程中存在的局限性，包括模型求解过程中的混频效应明显、求解精度受限、Nyquist采样准则与数据分辨率存在固有矛盾等，本文在模型求解中引入球面等面积ico_HEALPix网格剖分，设计了块最小二乘迭代算法，以损失一定计算效率为代价，实现了基于ico_HEALPix网格的超高阶地球重力场模型高精度构建。以XGM2019e模型作为输入真值，经过球谐综合与球谐分析试算结果表明：经过迭代，基于ico_HEALPix网格构建的超高阶次地球重力场模型精度可达到优于地理网格的效果，同时，该方案提高了原模型可存储的数据分辨率精度，减小了同等情况下的数据存储容量，解决了南北极点附近的畸变问题，提高了网格数据的利用率和可靠性。后续可探索其他新的均匀网格剖分方式来避免法矩阵的迭代求解以提升解算效率及ico_HEALPix网格在削弱混频效应等方面的应用优势研究。本文所构建的模型可为提升球谐系数模型构建精度，丰富构建地球重力场模型理论，完善球谐系数模型理论体系提供借鉴与参考。
初审：张艳玲复审：宋启凡
终审：金 君
资讯
○ 中南大学石岩教授：时空异常探测：从数据驱动到知识驱动的内涵转变与实现路径|《测绘学报》2024年53卷第8期