将GPS高程测量测得的大地高转换为工程上应用的正常高，关键在于高程异常的求解^{[1, 2]}。由美国地理空间情报局(National Geospatial Intelligence Agency,NGA)发布的空间分辨率为5′(约为9 km)的EGM2008重力场模型能够精确表达高程异常的长波项^{[3, 4, 5, 6]}。高程异常的短波项由地形变化引起^[7]，在地形起伏大的地区，地形变化是影响高程异常的重要因素，而EGM2008重力场模型不能很好地解决地形变化产生的短波项的影响。在实际工程应用中求解高程异常时，常利用数值逼近方法考虑地形因素的影响^{[8, 9]}，而忽略高程异常的长波部分；或只考虑长波项的影响，不考虑地形因素的影响^{[10, 11, 12, 13]}；亦或结合EGM2008重力场模型并考虑地形影响，求解地形影响时采用重力引力位法、棱柱积分法等方法^{[14, 15]}，计算过程复杂。

利用“移去-拟合-恢复”法的思想，本文提出了一种基于EGM2008模型求解高程异常长波项，利用数字高程模型(digital elevation model,DEM)格网数据求解地形起伏变化引起的高程异常短波项，并结合二次曲面拟合法拟合残余高程异常的GPS高程拟合方法。通过港珠澳大桥的实例数据，将仅考虑地形因素的二次曲面拟合，仅结合EGM2008模型的二次曲面拟合方法与本文提出方法的GPS高程拟合精度进行比较。

本文选择“移去-拟合-恢复”法求解高程异常^[16]，该方法综合利用EGM2008模型与数学拟合方法，无须实测重力数据，只须联测少量水准点，尤其是在地形起伏较大、水准联测点较少的地区，能方便地进行GPS高程拟合。根据物理大地测量学理论，高程异常ζ可以表示为^[12]：

式中，ζ_EGM是由EGM2008模型求得的高程异常长波部分；ζDEM是地形起伏变化引起的高程异常短波部分；ζres是残余高程异常。“移去-拟合-恢复”法思路为，从已知点的ζ中移去ζ_EGM与ζ_DEM，获得ζ_res进行模型拟合；基于拟合模型计算待求点的ζ_res，再加上ζDEM与ζEGM，就可得到待求点的高程异常ζ，具体处理流程见图 1。

图 1 “移去-拟合-恢复”法求解高程异常的流程图 Fig. 1 Flowchart of the Height Anomaly Determination Based on the Remove-fit-restore Approach

图选项

EGM2008模型是由NGA发布的全球超高阶重力场模型，它以PGM2007B为参考，综合利用GRACE卫星重力数据、全球5′×5′重力异常数据、TOPEX卫星测高数据、地形数据以及地面重力数据^{[17, 18]}。该模型的阶次完全至2 159^[19](另外球谐系数的阶扩展至2 190，次为2 159)，先利用EGM2008模型求得扰动位T^[20]，再根据Burns公式计算地面上任一点的高程异常长波项：

式中，ρ、φ和λ分别为计算点地心向径、余纬和经度；GM为地心引力常数；a为椭球长半径；C_nm、S_nm为完全规格化位系数;P_nm(sinφ)为完全规格化缔合函数；γ为计算点的正常重力值。

由地形起伏变化引起的高程异常变化量δ可表示为^[16]：

式中，G为万有引力常数;ρ为地球质量密度;h_i为流动单元平均高程;h为待求点高程;γ为待求点正常重力;l为待求点到流动单元的距离;dxdy为流动单元面积。

高程异常短波项ζ_DEM的求解步骤如下。

1)获取所在区域的DEM数据。目前美国国家航空和宇宙航行局(The National Aeronautics and Space Administration，NASA)公开的DEM在中国境内的分辨率为3″(90 m)，根据拟合区域的坐标最大值与最小值选择对应的DEM格网数据。

2)划分流动单元。

3)求解流动单元的平均高程，即单元4个格网点的平均高程，利用式(3)可以求得每一个格网点的高程异常。

4)按双线性内插法计算待求点的高程异常短波项：

式中，a、b、c、d为流动单元的4个格网点(依次为左下角、左上角、右上角、右下角);l为流动单元边长;Δx、Δy为待求点与流动单元左下角格网点的坐标差。

残余高程异常主要由模型误差、数据误差等引起，可通过数学函数对其规律性进行拟合。常用的数学函数有二次曲面函数、平面函数、反距离加权内插函数等^{[11, 12]}。在实际工程中，大多采用二次曲面函数拟合残余高程异常^{[14, 15, 21]}，其函数表达式为：

式中，ζ_res表示残余高程异常;x、y表示点的重心化坐标;a₀、…、a₅表示拟合系数。联测水准点的个数应不少于6个，使用最小二乘法进行求解。

本文利用某跨海大桥工程项目的实测数据进行高程拟合方法分析。该工程数据点分布在长约53 km、宽约28 km的范围内。使用GPS接收机观测了39个点，并用精密水准测量方法联测了其中的13个点。将7个联测点作为已知点(高程均值33.036 7 m，高程最大值106.174 6 m，最小值1.036 4 m)，6个联测点作为检核点(高程均值37.890 9 m，高程最大值89.489 3 m，最小值-0.228 6 m)，点位平面分布情况见图 2。图 2显示，用于拟合模型的7个已知点基本覆盖工程范围，6个检核点均在拟合范围之内。

图 2 GPS点位分布图 Fig. 2 Site Distribution of GPS Stations

图选项

分别使用表 1中的三种方法(分别为仅考虑地形因素的二次曲面拟合,仅结合EGM2008模型的二次曲面拟合方法与本文方法)进行GPS高程拟合。表 1中“是”/“否”表示考虑/不考虑相关影响，在考虑地形因素时，本文将流动单元划分为270 m×270 m的格网来求解由地形起伏变化引起的高程异常变化量。

对以上三种方法按式(6)、(7)进行内外符合精度的评定：

式中，STD表示内符合精度;ε为已知点的高程异常残差;n为已知点的个数；RMS表示外符合精度;Δ为检核点的高程异常残差;m为检核点的个数。内符合精度反映了利用已知点进行建模的模型质量，外符合精度反映了利用所建模型求解检核点高程异常的整体质量。

图 3和图 4分别是三种方法拟合已知点和检核点高程异常的残差比较。从图 3可以看出，7个已知点在本文提出的方法下拟合残差均最小，相对方法1和方法2均有mm级的改进精度；本文方法的高程异常残差最大值为0.26 cm，小于方法1的0.70 cm和方法2的0.72 cm，其残差均值为-0.01 cm，也小于方法1的0.13 cm和方法2的0.02 cm；本文提出方法的内符合精度为0.16 cm，优于方法1的0.43 cm和方法2的0.41 cm。实验结果表明，利用本文提出的方法能够提高建模精度。

图 3 不同方法已知点高程异常残差比较 Fig. 3 Comparison of the Height Anomaly Residuals of Known Points by Different Methods

图选项

图 4 不同方法检核点高程异常残差比较 Fig. 4 Comparison of the Height Anomaly Residuals of Checked Points by Different Methods

图选项

从图 4可以看出，6个检核点在本文方法中高程异常残差最小，精度较方法1和方法2均有mm级改进；本文方法的高程异常残差最大值为1.72 cm，小于方法1的1.96 cm和方法2的1.82 cm；本文方法的外符合精度达到了0.96 cm，优于方法1的1.61 cm和方法2的1.20 cm。

改变已知点与检核点的个数，用上述三种方法进行验证，其拟合效果见表 2。需要注意的是，当用1个检核点进行检核时，为防止因已知点的选择而造成拟合偏差过大，本文用交叉验证法^[22]进行处理。

表 2显示，在改变已知点与检核点的个数后，本文提出方法的内、外符合精度较方法1、方法2均有mm级的改进；已知点与检核点的个数虽然发生变化，但本文提出方法的内符合精度均达到mm级，外符合精度均达到cm级。实验结果表明,本文提出的方法能够提高GPS高程拟合的精度。

利用“移去-拟合-恢复”的思想，本文提出了一种基于EGM2008模型，并考虑地形因素影响的二次曲面拟合GPS高程的方法。相比以往的高程异常求解方法，本文提出的方法精度较高。通过某跨海大桥工程的实测数据，比较分析了仅考虑地形因素的二次曲面拟合法、仅结合EGM2008模型的二次曲面拟合法和本文提出的基于EGM2008模型和地形改正的二次曲面拟合法的拟合效果。结果表明，本文提出的方法精度较高，其内符合精度达到0.16 cm，外符合精度达到0.96 cm，能够提高GPS高程拟合精度。此外，由于本文的研究工作是基于工程项目建设的实际地形特征完成的，提出的方法对于其他地形区域的适用性还需要进一步分析和验证。结合多种地形特征，进一步完善计算模型，提高模型的通用性是下一步研究的重点。