马 健，魏子卿，任红飞

1. 地理信息工程国家重点实验室，陕西 西安 710054； 2. 西安测绘研究所，陕西 西安 710054； 3. 信息工程大学地理空间信息学院，河南 郑州 450001

Stokes/Hotine积分是计算(似)大地水准面的理论工具[1-6]。Stokes/Hotine积分理论上需要全球积分，但由于全球重力数据难以获取，且全球积分计算量较大，因此实际计算通常使用移去恢复技术，此时边值解算仅需进行近区Stokes/Hotine积分。为了减弱远区影响、提高(似)大地水准面的解算精度，出现了众多核函数的修正方法。核函数修正方法可分为确定性修正与随机性修正两类。确定性修正理论包括有Meissl方法[7]、Molodensky方法、Sjöberg方法[8]、Wong和Gore方法等[9-10]，随机性修正需将模型重力数据和实测重力数据的误差阶方差作为先验信息，而实测数据的误差阶方差信息很难获取，限制了随机性修正方法的应用[11]。

使用移去恢复技术时需选定一个模型重力场作为参考场，因此Stokes/Hotine积分中包含了实测重力和模型重力两类重力数据。重力测量、水准测量、地形归算[12-13]、格网化[14]等使Stokes/Hotine积分中的实测重力数据存在一定的长波误差[15]，而重力场模型的长波精度较高(satellite-only卫星重力场模型的前20阶已达到非常高的精度[16-17])，因此Stokes/Hotine积分时，应选择对实测重力数据具有高通滤波属性的积分核函数。截断核函数(the spheroidal kernel，由标准核函数除去一定阶数的Legendre多项式得到)是一种非常有效的高通滤波器[18]，可有效地避免实测重力数据的长波误差传播到(似)大地水准面中，因此在基于移去恢复的边值解算中应用较为广泛。但模型重力数据的累积误差随着模型阶数的增加而增大，而实测重力数据的中高频误差相对较小，因此为了获得高精度的(似)大地水准面不宜过滤过多频段的实测重力信息。

传统积分核函数的研究主要为了减小远区截断误差，但在移去恢复模式下，积分核对边值解算精度的影响方式发生了转变：在高频波段(大于模型重力场截止阶数的频率波段)，核函数决定了高频远区截断误差的大小以及谱泄露的程度；在低频波段(小于模型重力场截止阶数的频率波段)，核函数决定了实测与模型两类重力数据对低阶高程异常/大地水准面高的贡献率。因此，在移去恢复模式中核函数的研究应分高频和低频两个频段展开。

Stokes积分核与Hotine积分核[19-21]具有相似的频谱特性，Stokes核函数理论可方便地应用于Hotine积分核，反之亦然。根据Stokes-Helmert理论方法构建的加拿大重力大地水准面CGG2010[22]使用了高低阶均修正的截断Stokes核函数。本文以Hotine积分核为例在其基础上提出了更加高效的低阶修正的截断核函数，不仅包括余弦低阶修正核函数，还构造了一种线型低阶修正核函数。本文首先对扩展Hotine核函数展开频谱分析，然后根据Hotine-Helmert理论方法使用不同的Hotine核函数解算了试验区的似大地水准面，以验证不同Hotine核函数的应用效果。

1 Hotine截断核函数

在移去恢复模式下，使用传统截断Hotine核函数计算高程异常(高程异常为似大地水准面到参考椭球面的距离)的算法公式为

(1)

式中，ζ为高程异常；ζM为模型高程异常；C1为Hotine的近区积分范围；R为地球平均半径；r为计算点的地心向径；γ为地面点的正常重力；ψ为流动点与计算点间的球心角距；δgT为地面重力扰动；δgM为模型重力扰动；*表示向下延拓过程[23-26]，通过向下延拓可将地面重力值转化为边界面重力值；Hbl(r,ψ)为扩展形式的Hotine截断核函数，为方便表述下文省略“扩展”。

目前广泛使用的Hotine截断核函数有两种。一种是高低阶均截断的核函数

(2)

式中，L为低阶截断频率(阶数)；M为由地形分辨率决定的高阶截断频率；Pn(cosψ)为Legendre函数。另一种为仅低阶截断的核函数

(3)

式中，Hstd为标准形式的Hotine积分核

(4)

2 修正的Hotine截断核函数

修正的截断核函数是通过对传统截断核函数增加修正因子得到的，即

(5)

式中,αn为修正因子，也可称为平滑因子。

2.1 高低阶均修正

文献[22]提出了余弦函数构建的高低阶均修正的截断核函数，其修正因子为

(6)

式中，μ为低阶修正带宽；υ为高阶修正带宽，修正带宽表示修正的阶数区间。从式(6)可以看出该修正因子在L-μ阶到L阶的频率波段从0逐渐增大至1，在M阶到M+υ阶的频率波段由1逐渐减少至0。当修正区间(L-μ,L)和(M,M+υ)的各阶修正因子均取0时，高低阶均修正的截断核函数转化为式(2)所示的传统高低阶均截断的核函数。

由于边值解算采用的重力数据、地形数据均为离散的网格数据，因而需对Hotine积分进行离散化。当对中央网格(计算点地心向径穿过的网格)进行离散化时，通常将中央网格近似为圆形区域。采用高低阶均修正的Hotine截断核函数时，通过推导可得出中央网格重力扰动对高程异常贡献ζ0的近似公式为

(7)

2.2 低阶修正

本文在高低阶均修正的截断核函数的基础上提出了仅低阶修正的截断核函数，两种修正核函数的差别在于修正因子的不同。低阶修正的截断核函数的修正因子αn为

(8)

从式(8)可看出，低阶修正核函数的修正因子在L-μ阶到L阶的频率区间由0逐步增大至1。当低阶修正频率区间(L-μ,L)的各阶修正因子均取0时，低阶修正核函数转化为传统仅低阶截断的核函数。式(8)所示的低阶修正核函数的修正因子是由余弦函数构造的。本文进一步提出一种线型低阶修正的核函数，其修正因子为

(9)

式(9)所示的线型低阶修正核函数的修正因子同样具有在L-μ阶至L阶的频段由0逐步增大至1的特点。图1为两种低阶修正核函数的修正因子的取值。

由图1可知，对于余弦修正核函数，修正因子在修正初始频段和修正末尾频段变化相对较慢，而在修正的中间频段变化相对较快；对于线型低阶修正核函数，修正因子在修正频段保持固定的变化速度。修正因子的类型及参数设置与实测重力数据和模型重力数据的相对精度有关。当使用低阶修正的Hotine截断核函数时，通过推导可得出中央网格重力扰动对高程异常贡献ζ0为

(10)

式中，lP0(ψ0)的计算公式为

3 移去恢复模式下Hotine积分的谱分解式

在移去恢复模式下，Hotine积分仅需近区积分。根据勒让德函数的球面积分算法和Molodensky截断理论，近区Hotine积分可表示为如下谱形式

(11)

(12)

式(11)、式(12)不仅适用于修正形式的Hotine截断核函数，对于传统Hotine截断核函数也是适用的，此时αn只有0和1两种取值(传统Hotine截断核函数为修正Hotine截断核函数的特例)。根据式(11)所示的谱表示方法可对式(1)进行如下改化

(13)

4 试验与分析

4.1 频谱分析

根据上文分析，在移去恢复模式下核函数在高、低频波段对高程异常的影响方式不同，因此对核函数的频谱分析应分高、低频两个波段进行。本节频谱分析过程中将计算点高程设为1000 m，Hotine积分半径设为1°，低、高阶截断频率(阶数)分别取为360、5400。图2为使用不同的Hotine核函数时实测重力数据对高程异常的高阶贡献率，其中修正因子采用余弦函数形式，低、高阶修正带宽均为180阶。

根据图2，在高频波段的初始频段，高低阶均截断与低阶截断的Hotine核函数对高程异常的高阶贡献率差异很小，而高低阶均修正与仅低阶修正的Hotine核函数的高阶贡献率的差异也非常小。贡献率越接近1，重力数据越能够有效地传播到解算的高程异常中。图2中传统截断核函数在高频波段的初始频段贡献率较低，说明传统截断核函数存在较严重的谱泄露现象。修正的截断核函数能够有效地控制谱泄露现象，提高了实测数据在高频波段初始频段的贡献率。在高阶截断频率附近，几种核函数中高低阶均截断的核函数存在比较明显的谱泄露现象。图3为使用不同的Hotine核函数时模型和实测数据对高程异常的低阶贡献率。

图1 两种修正因子的取值Fig.1 Values of the two modification factors

图2 实测数据的高阶贡献率Fig.2 High-degree contribution rates of the measured data

图3 模型和实测数据的低阶贡献率Fig.3 Low-degree contribution rates of the measured data and model data

从图3可以看出，在低频波段，两种截断Hotine核函数的低阶贡献率差异很小，两种修正Hotine核函数的低阶贡献率差异也很小。在低频波段的初始频段，使用截断和修正形式的核函数时模型数据的贡献率均接近1，说明在低频波段的初始频段高程异常主要来自模型数据的贡献。在修正频段(L-μ～L，此处为180～360)，使用修正核函数时实测数据对高程异常的贡献多于使用传统截断核函数时实测数据对高程异常的贡献。

下面分析不同修正带宽下实测和模型数据对高程异常的贡献率的变化，使用的核函数分别为余弦低阶修正和线型低阶修正的Hotine截断核函数。图4所示为不同修正带宽下实测数据的高阶贡献率。

图4 不同修正带宽下实测数据的高阶贡献率Fig.4 High-degree contribution rates of the measured data using different modification bandwidths

由图4可知，在高频波段的初始频段，对于两种(余弦与线型)低阶修正的Hotine核函数，低阶修正带宽为90阶时实测重力数据的贡献率均明显小于修正带宽为180、270、360阶时的贡献率，说明90阶修正带宽下两种低阶修正核函数对谱泄露现象的改善程度不及180、270、360阶修正带宽下修正核函数对谱泄露的改善程度。图4中，使用余弦与线型低阶修正核函数时实测重力数据的高阶贡献率相差不大，但当修正带宽为360阶时，使用线型低阶修正核函数时实测重力数据的高阶贡献率更接近1，说明此时线型低阶修正核函数稍优于余弦低阶修正核函数。图5所示为不同修正带宽下模型和实测重力数据的低阶贡献率。

图5 不同修正带宽下实测和模型数据的低阶贡献率Fig.5 Low-degree contribution rates of the measured data and model data using different modification bandwidths

从图5可以看出，使用两种低阶修正核函数时实测与模型重力数据对低阶高程异常的贡献率相差不大。在未修正的频率区间(0～L-μ)，低阶高程异常的贡献主要来自于模型重力扰动；在修正的频率区间(L-μ～L)，实测重力扰动的贡献逐渐增大，而模型重力扰动的贡献逐渐减少。从图5中还可看出，修正带宽越大，模型重力数据的贡献越少而实测重力数据的贡献越多，因此可通过调整修正带宽控制实测和模型重力数据对低阶高程异常贡献的权重。

4.2 应用试验

为了说明不同Hotine核函数的应用效果，针对不同形式的Hotine核函数解算似大地水准面进行了试验。似大地水准面的区域范围为108°E—114°E，28°N—32°N，该区域海拔最高2700 m，平均607 m。本试验收集了105.5°E—116.5°E，25.5°N—34.5°N范围内的70 822个地面离散重力点，剔除粗差后，将剩余的70 379个实测重力点作为边值解算的基础重力数据。本试验边值解算采用基于移去恢复的Hotine-Helmert理论方法[27]，其中地形直接、间接影响采用文献[28]给出的算法公式。试验区共有68个GPS水准点，利用试验区GPS水准点对截断至360阶和2190阶的EIGEN-6C4、EGM2008重力场模型分别进行精度检核，其结果统计于表1。

表1 截断到不同阶数的重力场模型精度

Tab.1 Accuracies of the gravity field models with different cutoff degreesm

表1中，截断到360和2190阶的EIGEN-6C4模型精度(标准差)分别为±28.2 cm与±7.8 cm，而截断到360和2190阶的EGM2008模型精度分别为±29.5 cm与±10.6 cm，由此可看出试验区EIGEN-6C4模型精度优于EGM2008模型。

下面首先利用两种截断(高低阶均截断、低阶截断)Hotine核函数计算似大地水准面，为了提高计算速度，本试验采用2′×2′的格网分辨率。表2统计了解算的重力似大地水准面的精度与Hotine积分所需时间，其中Hotine积分半径为1°，高阶截断阶数取5400阶(根据Nyquist采样定律可知2′×2′分辨率格网对应5400阶的高阶截断频率)，低阶截断阶数取360，参考模型取EIGEN-6C4模型的前360阶。

表2 两种传统截断核函数计算的重力似大地水准面精度

Tab.2 Accuracies of the gravimetric quasi-geoids using the two traditional spheroidal kernels

核函数最小值/m最大值/m平均值/m均方差/m标准差/m耗时/min高低阶均截断-0.0580.4270.1260.1520.086150.0低阶截断-0.0580.4260.1260.1520.08710.5

从表2可以看出，两种传统截断核函数计算的似大地水准面精度基本没有差别，但图2中两种传统截断核函数在高阶截断频率附近差异较大，这反映了高程异常对高频信息不敏感的频谱特性。根据表2，在解算精度一致的前提下，使用低阶截断核函数时Hotine积分所需的计算时间较高低阶均截断的核函数明显缩短，因此对于两种截断核函数，低阶截断的核函数比高低阶均截断的核函数更适于边值解算。

表3比较了两种修正(高低阶均修正、低阶修正)Hotine核函数在边值解算精度与计算时间上的差别，其中低阶修正带宽取180阶，高阶修正带宽分别取180、540阶，修正核函数采用余弦修正形式，其他参数设置与表2相同。

表3 两种修正核函数计算的重力似大地水准面精度

Tab.3 Accuracies of the gravimetric quasi-geoids using the two modified spheroidal kernels

核函数高阶修正带宽/阶最小值/m最大值/m平均值/m均方差/m标准差/m耗时/min高低阶均修正1800.0080.4080.1500.1660.072156.75400.0080.4080.1500.1660.072166.8低阶修正—0.0080.4070.1500.1660.07212.3

从表3可得出，当低阶修正带宽一定时，两种修正核函数解算的似大地水准面的精度基本相同。另外还可看出，高阶修正带宽对似大地水准面的精度的影响非常小。由于低阶修正核函数较高低阶均修正核函数大大缩短了计算时间，而二者的计算精度相同，因此低阶修正的核函数比高低阶均修正的核函数更适于边值解算。此外，对比表2、表3可看出，低阶修正核函数计算的似大地水准面精度优于低阶截断核函数的计算精度，二者在计算效率方面的差别也不大，因此低阶修正的截断核函数的应用效果优于其他核函数。

为了说明低阶修正带宽对边值解算精度的影响，表4统计了使用不同的低阶修正带宽时求解的重力似大地水准面的精度，其中低阶修正核函数采用余弦修正形式，地形分辨率为1.5′×1.5′，Hotine积分半径为1°，同样以EIGEN-6C4重力场模型的前360阶作为参考模型。

表4 余弦低阶修正核函数计算的重力似大地水准面的精度

Tab.4 Accuracies of the gravimetric quasi-geoids using the cosine low-degree modified kernels

低阶修正带宽/阶最小值/m最大值/m平均值/m均方差/m标准差/m0(未修正)-0.0690.3490.1230.1470.08290-0.0360.3330.1310.1480.071180-0.0070.3090.1430.1570.066270 0.0050.2620.1410.1530.059360-0.0220.1940.1040.1150.048

低阶修正带宽取0时低阶修正核函数实际即为传统的低阶截断核函数。从表4可知，低阶修正带宽对似大地水准面的精度具有较大的影响。在本文试验区360阶低阶修正带宽解算的似大地水准面精度最高(±4.8 cm)，而传统截断核函数计算的似大地水准面精度仅为±8.2 cm。似大地水准面精度的提高一方面是由于低阶修正核函数有效地控制了传统截断核函数存在的谱泄露问题，另一方面是由于低阶修正核函数增大了实测重力数据在低阶修正频段对高程异常的贡献。

表5统计了不同低阶修正带宽下线型低阶修正Hotine核函数计算的重力似大地水准面精度，其参数设置同表4。

表5 线型低阶修正核函数计算的重力似大地水准面的精度

Tab.5 Accuracies of the gravimetric quasi-geoids using the linear low-degree modified kernels

低阶修正带宽/阶最小值/m最大值/m平均值/m均方差/m标准差/m90-0.0360.3330.1310.1490.071180-0.0100.3060.1420.1560.065270-0.0060.2460.1300.1410.056360-0.0700.1580.0730.0870.048

比较表4、表5可知，在相同的低阶修正带宽下，线型低阶修正与余弦低阶修正核函数计算的似大地水准面精度基本一致。将表4、表5的结果与表2进行比较可以看出低阶修正核函数的应用效果较好。

由于不同重力场模型的精度不同，而修正带宽起到调整模型和实测重力数据在高程异常中贡献的比重的作用，因此使用不同的参考模型时适宜的低阶修正带宽也存在差异。为此，表6统计了以EIGEN-6C4模型与EGM2008模型的前360阶作为参考模型时不同低阶截断阶数下适宜的低阶修正带宽(即使用该低阶修正带宽时边值解算精度相对较高)及相应的似大地水准面精度，其中核函数采用余弦低阶修正核函数。

表6 以重力场模型的前360阶作为参考模型时解算的重力似大地水准面精度

Tab.6 Accuracies of the gravimetric quasi-geoids with the first 360 degrees of the gravity field model as the reference field

重力场模型低阶截断阶数/阶适宜的低阶修正带宽/阶最小值/m最大值/m平均值/m均方差/m标准差/mEIGEN-6C4模型360360-0.022 0.194 0.104 0.115 0.048240240-0.1070.1810.0600.0790.0511200-0.0770.2550.0910.1050.053EGM2008模型360360-0.604 -0.092 -0.304 0.314 0.081240240-0.646-0.144-0.3490.3560.07112060-0.687-0.189-0.3830.3900.074

表6中，当以EIGEN-6C4模型的前360阶作为参考模型时，低阶截断阶数与低阶修正带宽均取360解算的似大地水准面精度最高(±4.8 cm)，而当低阶截断阶数取120时，核函数不进行修正时边值解算精度相对较高。当以EGM2008模型的前360阶作为参考模型时，低阶截断阶数与低阶修正带宽均取240解算的似大地水准面精度最高(±7.1 cm)，但不及以EIGEN-6C4模型的前360阶作为参考模型的边值解算精度。

表7统计了采用2190阶的EIGEN-6C4模型与EGM2008模型作为参考模型时，不同低阶截断阶数下适宜的低阶修正带宽及相应的似大地水准面精度。

表7中，当以2190阶的EIGEN-6C4模型与EGM2008模型作为参考模型时，低阶截断阶数取2190阶计算的似大地水准面精度并不理想。根据上文频谱分析可知，当低阶截断阶数取2190阶时解算的高程异常主要来自于模型数据的贡献，因此不能有效地利用实测数据的信息。当2190阶EIGEN-6C4模型作为参考场时边值解算精度(最高达±4.5 cm)优于2190阶EGM2008模型作为参考场时边值解算精度(最高达±6.8 cm)。对比表6可以看出，以2190阶EIGEN-6C4模型作为参考模型较以360阶EIGEN-6C4模型作为参考模型求解的似大地水准面的精度提高不明显(分别为±4.5 cm与±4.8 cm)，这反映了远区高程异常的高频(361～2190阶频段)信息非常少。

表7 以2190阶重力场模型作为参考模型解算的似大地水准面精度

Tab.7 Accuracies of the gravimetric quasi-geoids with the 2190 degree gravity field model as the referencefield

重力场模型低阶截断阶数/阶适宜的低阶修正带宽/阶最小值/m最大值/m平均值/m均方差/m标准差/mEIGEN-6C4模型2190360-0.1100.3240.1270.1520.084360360-0.0160.1910.1040.1140.047240240-0.0900.1550.0610.0770.0471200-0.0420.2040.0910.1010.045EGM2008模型2190360-0.646-0.049-0.2790.3000.113360360-0.604-0.092-0.3040.3140.081240240-0.629-0.170-0.3480.3540.06812060-0.687-0.189-0.3830.3900.074

5 结 论

利用基于移去恢复技术的Stokes-Helmert边值理论或Hotine-Helmert边值理论求解(似)大地水准面时，截断形式的Stokes/Hotine核函数是一种有效的高通滤波器。但传统的截断核函数存在谱泄露现象，影响了(似)大地水准面的解算精度，因此本文引入并进一步发展了一种低阶修正的截断核函数。本文的创新点与结论主要有：

(1) 在余弦函数构造的高低阶均修正的核函数的基础上提出了余弦低阶修正核函数，进一步提出了线型低阶修正核函数，并给出了使用高低阶均修正核函数与低阶修正核函数时Hotine积分的中央区算法。

(2) 对移去恢复模式下核函数在解算似大地水准面中的作用进行了频谱分析，说明了在移去恢复模式下核函数研究的内容及内涵已发生转变：在高频波段，核函数影响着高频远区截断误差的大小以及谱泄露的程度，但在低频波段，核函数决定了实测重力数据与模型重力数据对高程异常的贡献率。通过对不同核函数的频谱分析得出，修正核函数能够有效地控制谱泄露现象，并且增大了实测数据在修正频段对高程异常的贡献率。余弦低阶修正与线型低阶修正核函数的频谱特性比较接近。

(3) 试验结果表明，采用相同的低阶修正带宽时，线型和余弦低阶修正的核函数计算的似大地水准面精度与高低阶均修正的核函数的计算精度一致，均优于传统截断核函数计算的似大地水准面精度。在计算效率上，低阶修正的核函数明显优于高低阶均修正的核函数。

(4) 通过本文的试验可看出，在基于Helmert第二压缩法的边值问题中低阶修正的截断核函数能够有效地改善边值解算的精度，具体应用时应结合参考模型精度、试验区实测重力数据精度、模型截止阶数等因素确定适宜的低阶截断阶数以及低阶修正带宽。

本文的研究是对文献[29]提出的厘米级精度要求下重新研究核函数的改进和构造截断函数问题进行的理论探索与实践，在Stokes-Helmert/Hotine-Helmert边值理论的研究和应用方面具有一定的参考价值。