基于超高阶地球重力场模型的GNSS 高程转化方法

赵保成，徐健，徐坚

(1.长江科学院 空间信息技术应用研究所,武汉 430010；2.武汉市智慧流域工程技术研究中心,武汉 430010)

0 引言

近年来，随着卫星导航定位技术和移动网络通讯技术的快速发展，基于连续运行参考基站(CORS)的网络实时动态(RTK)技术发展势头迅猛，千寻知寸FindCM、中国移动CORS、北斗CORS、厘清/Locate-CM 以及各省的CORS 服务等极大地丰富了测绘技术人员的技术手段.借助于CORS 服务，在移动网络覆盖区域，单人使用单台全球卫星导航系统(GNSS)接收机便能够瞬时获取在指定参考坐标系下的厘米级三维(3D)定位结果，相比于传统GNSS 应用的电台“1+N”模式，该服务使测绘作业效率得到了极大地提升.但不可忽略的是，无论采用何种CORS 服务，卫星接收机直接获取的高程信息始终是基于参考椭球面的大地高，而在我国的高程体系中使用的却是基于似大地水准面的正常高，二者高程系统不一致导致GNSS 在高程测量中的应用受到了一定的限制.为了使GNSS 测得的大地高直接应用于工程建设中，必须要进行高程数据转换[1-3].

GNSS 高程转换方法主要有两种.一种是GNSS水准高程法，其原理为利用均匀分布于测区的同名点的大地高和水准高(正常高)成果通过数学函数模型确定局部似大地水准面，由于一般测区内很少布设有高精度的GNSS 水准点，如果从外部引测，时间及人工成本高，因此该方法的适用范围有限；另一种方法是GNSS 重力高程法，其原理是利用地面或者卫星重力资料求解待测点的高程异常值，结合GNSS 接收机获取高精度大地高，进而求得待测点正常高[4].

随着全球卫星重力测量技术的不断发展，利用地面重力数据结合卫星重力数据获取的超高阶地球重力场模型的精度和分辨率在不断提高.基于此，众多行业学者尝试利用超高阶地球重力场模型进行了诸多系统GNSS 高程转换方法研究.文献[5]利用EGM2008 地球重力场模型去除高程异常中的长波项，进而使用加权组合模型对高程异常的剩余项拟合，拟合效率得到了一定的提升；文献[6]基于二次曲面和EGM2008 地球重力场模型，构建了多种移去-恢复型函数拟合模型，探究了各种模型的适用性；文献[7]利用矿区水准观测数据和EGM2008 地球重力场模型，综合物理和几何方法建立了矿区局部似大地水准面模型；文献[8]利用EIGEN-6C4 地球重力场模型进行GNSS 高程拟合，拟合精度得到了进一步提升；文献[9]详细阐述了利用地球重力场模型位系数计算高程异常的数学模型，同时实现了全球高程异常的可视化；文献[10]利用EIGEN-6C4、EIGEN-6C2、EGM2008 地球重力场模型结合移去-恢复法进行了高程转换方法研究，得出了在我国顾及EIGEN-6C4地球重力场模型的GNSS 高程转换精度更高的结论；文献[11]通过GNSS 高程转换精度分析，得出了XGM2019 重力场模型与我国似大地水准面符合程度更高的结论；文献[12]讨论了在不同地形条件下利用EGM2008 地球重力场模型转换大地高为正常高的方法；文献[13]利用常规水准限差评定了基于EGM2008 地球重力场模型的GNSS 拟合高程等级及适用性.

综上，目前GNSS 高程转换方法的研究存在两方面问题，一是研究虽引入了高阶地球重力场模型，但多是针对单种地球重力场模型并结合不同插值拟合方法，鲜有对利用目前已存的各种超高阶地球重力场模型进行高程转换方法效果的综合评估；二是插值拟合方法多是从拟合函数的数学模型角度出发，而对于参与插值拟合的特征点的选取方法及原则并未做深入研究.本文在总结了多位行业学者研究成果的基础上，针对截至2022年已存的5 种超高阶地球重力场模型，结合改进的“移去-拟合-恢复”法的GNSS 高程转换方法开展了研究，本研究针对在缺乏高程基础控制资料及地面重力资料测区的高程系统的建立工作上提供解决方案及思路.

1 高程系统简介

目前，常用的高程系统有三种，分别是以参考椭球面为基准面的大地高系统、以大地水准面为基准面的正高系统及以似大地水准面为基准面的正常高系统.大地高H为地面点沿法线至参考椭球面的距离；正高h为地面点沿铅垂线至大地水准面的距离；正常高hr为地面点沿铅垂线至似大地水准面的距离[14].在不考虑垂线偏差影响情况下，三种高程系统之间的几何关系如图1 所示，数学关系由式(1)～(2)所示，式中ζ 表示高程异常，N表示大地水准面差距.我国法定使用的1985 国家高程基准的高程系统即是正常高系统.

图1 高程系统关系图

2 原理与方法

2.1 超高阶地球重力场模型

1)地球重力场模型(EGM2008)

EGM2008 地球重力场模型是由美国国家地理空间情报局通过先进建模技术与算法，以PGM2007B为参考，结合Topex 卫星测高数据、GRACE 卫星重力数据、地面重力数据、地形数据，得到拓展至2 190阶球谐系数的全球超高阶地球重力场模型.EGM2008地球重力场模型的地球覆盖率高达83.80%，模型的基本分辨率为5'×5'，最高分辨率为1'×1'.

2)利用新技术改进的欧洲地球重力模型(EIGEN-6C4)

EIGEN-6C4 地球重力场模型于2014年发布，其球谐系数的阶数和次数均为2 190.数据源主要为1985—2010年的LAGEO 数据、2003—2012年的GRACERL03GRGS数据、完整的GOCE-SGG 数据、DTU10 地面重力数据.

3)SGG-UGM-2 地球重力场模型

SGG-UGM-2 是由我国梁伟等学者于2020年采用球谐分析方法构建并发布的一个新的2 190 阶地球重力场模型，使用的数据包括卫星重力观测数据、卫星测高数据和EGM2008 模型重力异常数据.

4)XGM2019e_2159 地球重力场模型

该模型于2019年由慕尼黑技术大学天文和物理大地测量学研究所发布，其阶次完全至2 159(球谐系数的阶扩展至2 190)，空间分辨率5'×5'，数据来源包括GOCO 06 卫星重力模型、陆地和海洋重力异常和地面上的地形导出重力.

5)GECO 地球重力场模型

该模型于2015年发布，其球谐系数的阶为2 190，主要采用简单谱组合原理将GO_CONS_GCF_2_TIM_R5 模型和EGM2008 模型结合组成.

2.2 改进的移去-拟合-恢复法

根据物理大地测量学理论[15]，地球任意点高程异常通常由2 部分构成，分别为中长波项与短波项，如式(3)所示，高程异常的中长波分量变化较小，趋势相对平缓，短波分量主要由地形起伏引起，变化相对较大[16].

“移去-拟合-恢复法”是GNSS 高程转换的常规方法，其基本思想为：首先将参与拟合点的高程异常中变化较平稳的中长波分量移除，只对变化较大的短波项进行拟合，然后利用空间插值计算出待求点上的短波项，最后在待定点上恢复中长波分量.

然而，在拟合的步骤中，参与拟合点的选取往往具有一定的盲目性，极易造成拟合曲面关键信息的损失，最终导致高程转换精度降低.为了解决上述问题，本文提出了一种改进的“移去-拟合-恢复”法，其具体操作流程如图2 所示.

图2 改进的移去-拟合-恢复法流程图

1)确定测区范围，将测区网格化处理，格网尺寸根据项目成果要求的比例尺精度需求确定，利用全球重力场模型计算测区格网点的高程异常值，并将其作为真实高程异常值的中长波分量，进而构建基于高程异常中长波分量的似大地水准面模型并将其栅格化显示.

2)利用ArcMap 软件中的空间分析模块对上一步得到的似大地水准面曲面模型(基于高程异常中长波分量)进行地形表面及3D 分析，计算其高程异常曲面的极大值、极小值、坡向、坡度、曲率、粗糙度等参数，进而选择特征点作为GNSS 水准联测点参与后续的曲面拟合计算.对于特征点选取有以下原则：所选点位在似大地水准面上趋势变化较大；点位的个数需保证曲面拟合参数方程有解；点位的位置需覆盖整个测区.

3)利用GNSS 水准联测法获取特征点的高精度大地高、正常高，求得真实高程异常，进而移去拟合点上的高程异常中长波分量，将参与拟合点的平面坐标作为输入值，高程异常短波项作为输出值，构建二次曲面拟合模型.

4)利用地球重力场模型计算出待求点的高程异常中长波分量，再加上拟合步骤中该点高程异常短波项拟合值，最终得到待求点的高程异常值.

2.3 利用地球重力场模型求解高程异常中长波项

可以利用2.1 节介绍的5 种超高阶地球重力场模型获取地表任意点重力势的泛函值，例如重力异常、高程异常、垂线偏差等[15].利用Bruns 公式求解地球上任意点的高程异常中长波项ζGM，公式为

式中：ρ为计算点地心向径(计算点到地心的距离)；Ψ为余纬(纬度90°)；λ为大地经度；GM为包含大气层在内的地心引力常数；a为地球椭球长半径；、nm为完全规格化位系数；(sin Ψ)为完全规格化缔合函数；γ为计算点的正常重力位.

2.4 二次多项式曲面拟合

高程异常短波项的曲面拟合方法较多，常用的方法有二次曲面拟合、多面函数法、曲面样条函数法、移动曲面法等.在实际工程应用中，为方便计算，大多采用二次曲面函数拟合高程异常短波项，其数学表达式为

式中：ζ短波为特征点的高程异常短波项；(xi,yi)为特征点在某参考系统内的平面直角坐标；(a0,a1,a2,a3,a4,a5)为模型待定系数.

根据曲面拟合函数模型可知，共有6 个未知参数，因此为求解模型系数，至少需要列立6 个方程.当已知点个数大于6 时，表达式为

将上式变换为误差方程形式如式(7)所示，其中ε为模型的拟合残差：

根据最小二乘法原理，令模型拟合残差平方和最小(i≥6)，在此条件下，可以求得二次曲面模型待定系数X，如式(13)所示，其中，P为权矩阵，在实际应用中，常使用单位权矩阵代入：

3 算例验证

3.1 实验区概况

实验区位于湖北省某县城城区，此区域内相对高差较小，区内均匀地设置有20 个等精度GNSS 水准联测点，点位平面按照《全球定位系统(GPS)测量规范》中D 级要求施测，高程值采用《国家三、四等水准测量规范》进行三等水准联测.具体的点位分布如图3 所示.

图3 实验点位分布图

3.2 实验方案设计及拟合精度评定

按照2.2 节中介绍的拟合点选点方法选取图3中用红色三角形表示的6 个点作为拟合点，选取黑色圆形表示的14 个点作为检核点，开展高程异常拟合实验.实验方法采用改进的“移去-拟合-恢复”法，分别采用2.1 节介绍的5 种超高阶地球重力场模型进行高程异常中长波项的求解，拟合方法均采用二次曲面拟合法.

拟合效果一般通过对比拟合值与真值的符合程度来评估，主要指标有中误差、极差、偏度值、峰度值等.中误差对一组测量值中的异常值非常敏感，能较好地反映测量结果波动大小，因此本文主要使用中误差评定高程异常拟合精度，极差、偏度值、峰度值作为辅助参考值.中误差计算式为

式中：mζ为高程异常拟合中误差；Δ为在检核点上曲面拟合模型插值计算的高程异常值与GNSS 水准计算得出的高程异常值(真值)的差值；n为检核点的个数.

3.3 实验对比分析

如图4 所示，本实验区内使用5 种地球重力场模型求解的高程异常均为负值，基本与真实高程异常值的变化趋势一致，其中SGG-UGM-2 模型计算得到的高程异常与真实高程异常变化趋势符合程度更高.

图4 模型计算高程异常图

由表1 所示，5 种地球重力场模型计算的高程异常值与真值残差均较小，中误差均不超过±0.012 0 m，残差最大值为GECO 模型计算得出的结果，残差值为0.094 3 m.其中，SGG-UGM-2 地球重力场模型计算得出的高程异常值与真值符合程度更高，中误差为±0.009 3 m.

表1 模型高程异常精度评定表 m

由图5 及表2 可知，5 种方案的拟合中误差分别为±5.785 4 mm、±5.810 0 mm、±5.815 3 mm、±4.786 6 mm、±7.622 5 mm，拟合效果均比较优秀，中误差均不超过±8.000 0 mm，在此表明了基于超高阶地球重力场模型结合改进的“移去-拟合-恢复”法进行GNSS 高程转换的正确性与优越性.其中，利用XGM2019e_2159重力场模型的方案4 在各方案中的中误差、极差、偏度值、峰度值均为最小，分别为±4.786 6 mm、18.875 7 mm、-0.648 8、0.887 8，从而证明了利用该模型计算结果相对准确和可靠.而模型EGM2008、EIGEN-6C4和SGG-UGM-2 拟合精度大体一致，模型GECO 拟合精度相对最差.

图5 模型计算高程异常拟合残差图

表2 5 种方案的拟合精度评定表

4 结束语

基于5 种不同的超高阶地球重力场模型，本文在湖北省某县城城区开展了GNSS 高程转换方法实验研究，通过对实验结果对比分析，表明了基于超高阶地球重力场模型结合改进的“移去-拟合-恢复”法的GNSS 高程转换方法的可行性与可靠性.在不考虑我国高程基准与全球高程基准之间高差的影响下，在实验测区内，SGG-UGM-2 地球重力场模型计算得出的高程异常值与真值符合程度在5 种模型中最高.在改进的“移去-拟合-恢复”法实验中，基于XGM2019e_2159地球重力场模型的高程转换效果最好，高程异常拟合中误差、极差、偏度值和峰度值最小.

本文提出了基于地球重力场模型改进的 “移去-拟合-恢复”方法，相比于仅依赖数学模型的拟合方法更具有实际的物理意义.除此之外，该方法还在曲面拟合特征点的选择上更具有科学性和合理性，拟合残差更小.

在绝大多数的工程建设中，高程控制测量往往采用水准法或三角高程法，虽然保证了精度，但是上述两种方法的时间及人力成本较高.对于面积广阔、处于偏远地区、严重缺乏已有高程资料的测区的高程测量工作，本文的研究成果或将具有一定的参考指导意义.