利用矩谐分析的局部大地水准面确定

路媛琦，章传银，蒋 涛

(1. 青岛市勘察测绘研究院，山东 青岛 266032； 2. 中国测绘科学研究院，北京 100039； 3. 山东科技大学，山东 青岛 266590)

大地水准面是一个与静止的平均海水面重合并延伸到大陆内部的、包围整个地球的、封闭的重力位水准面，可以作为定义全球或区域高程系统的基准面[1]。确定大地水准面是建立高程基准统一的有效途径，并可结合GNSS实现大地高到正高的转换，实现卫星三维空间定位。

全球大地水准面的确定通常采用球谐展开模型，但确定局部大地水准面所利用的局部重力数据无法满足球谐函数的正交性，球谐展开模型并不适用于表达区域大地水准面。区域大地水准面确定的常用方法一种是基于Stokes或Molodensky理论求解大地测量边值问题，我国第二代似大地水准面模型(CQG2000)即是采用此类传统方法构建。另一种被许多学者应用的解析法为球冠谐分析；李建成等首次提出将球冠谐分析应用于区域重力场建模，证明了球冠谐展开表达区域重力场的可行性。球冠谐分析在大区域范围内建立大地水准面是一种理想方法，但应用于较小范围的大地水准面确定时，存在数值计算上的困难，精度较差。

确定局部大地水准面的谱方法除球冠谐分析外，还有一种方法是矩谐分析。矩谐分析在地磁学研究中应用广泛；Alldredge最早使用矩谐分析用于区域地磁场建模。1992年边少锋[1]将矩谐分析引入地球重力场逼近，由于矩谐分析的基函数是三角函数，球冠谐分析的基函数是非整阶缔合Legendre函数，展开到较高阶次时前者的计算远比后者稳定、快捷。

本文利用EGM-2008重力位模型，加入标准差为2 mGal的高斯白噪声，模拟陆地重力及飞行高度为4 km的航空重力值；并在研究区域内建立直角坐标系，求解地球引力位的Laplace方程，得到由矩谐系数表达的重力场参数表达式；将模拟重力数值代入矩谐展开式中求得矩谐系数、残余重力扰动及残余大地水准面；最后恢复大地水准面并与真实值进行对比，以此说明矩谐分析确定大地水准面的可靠性、准确性。

1 矩谐分析

已知地球外部空间某一区域或地球表面的重力观测值(重力异常或重力扰动)，求定该区域内地球重力场的矩谐系数，这一过程叫作矩谐分析[2]。

1.1 矩谐展开模型

矩谐分析(rectangular harmonic analysis，RHA)的基本原理是对地球某一地区构建区域大小的矩形面，以此近似球面；以矩形面的中心点为原点建立局部直角坐标系(如图1所示)；并在其中求解Laplace方程，在选定一定的截断水平后，引力场和引力位是包含有限个待定系数的已知函数，可以用所研究区域内一组观测点上的引力场值确定这些系数，并用这些系数表示不同波长的重力场信息，便于求解大地水准面[3-6]。

地球引力位满足Laplace方程[6-7]

(1)

Laplace方程在直角坐标系中表示为

(2)

采用分离变量法，引力位可写为

V(x，y，z)=f(x)g(y)h(z)

(3)

将式(3)代入式(2)得局部直角坐标系中Laplace方程的解[3]为

(4)

式中，Cnm为矩谐系数，即引力位系数；ψ(x)，ψ(y)为基函数，可分别表示为

(5)

其中

(6)

τnm定义为

(7)

式中，Dx、Dy分别表示矩形计算区域在东西(x)和南北(y)方向上的距离。

图1 矩谐分析使用的局部直角坐标系

实用上矩谐级数Cnm不可能展开到无穷阶数，应截断至某个最高阶数N和次数M，则扣除参考椭球的正常引位之后，扰动位可展开为如下有限阶次的矩谐级数

(8)

由Molodensky边值条件得重力异常与扰动位满足如下近似关系

(9)

于是联合式(8)和式(9)可推出重力异常的矩谐展开式为

(10)

式中，r为地心距离。

同理，重力扰动、大地水准面差距的矩谐展开式分别为

(11)

(12)

式中，γ0为正常重力。

1.2 截断阶数

在实际应用中，无论球谐级数还是矩谐级数都不可能展开到无穷多次，因此对级数表达式选择适当的截断阶数非常重要，过低的截断阶数会丢失有价值的信息，而过高的截断阶数会导致计算量增加，从而使结果不稳定，出现所谓的“龙格”现象[8-9]。因此需要根据数据的数量和精度要求[10]，确定模型的截断阶数Nmax。最优截断阶数的确定仍旧困难，通常得到最高截断阶数后不断进行尝试性的计算比较，以获取最优截断阶数。

最高截断阶数公式[8]为

最佳截断阶数公式[11]为

1.3 边界效应

对于区域大地水准面而言，其研究范围的局部性限制了重力场信号的完全周期性；但矩谐展开模型采用的周期性函数Fourier级数成立的前提是假定待求信号在计算区域内为周期函数[12]。式(8)在计算区域内收敛于真实地球扰动位，在计算区域边界处近似等于两边界扰动位的平均值，因此在边界处会产生Gibbs振荡现象。因此可扩展所求区域，以此降低Gibbs震荡现象的影响，使得计算区域的重力场信号具有周期性[13-15]。

举例说明，若研究范围大小的矩形长、宽分别为为Dx、Dy，数据点范围长、宽同为dx、dy；可将计算区域的平面范围扩展为Dx+dx、Dy+dy，此时待估重力场信号在区间D+d上满足周期性条件，此时式(7)为

当扩展参数增大到一定大小时可以有效降低边缘效应，但过大的扩展参数会使求解过程存在病态，同样影响带球扰动位系数的估算精度[3]。

2 数值计算与分析

2.1 基于地面重力数据的算例与分析

为验证矩谐分析基于陆地重力试验数据的可靠性，设计如下试验：获取纬度为[29°，34°]、经度为[111°，116°]围成的5°×5°区域，Dx=504.37 km，Dy=601.34 km，格网间隔为2.5′×2.5′，共计14 400个计算点。利用EGM2008重力位模型的2～2190阶系数计算得到重力异常的模拟观测值，并加入数学期望为0、标准差为2 mGal的高斯白噪声；移去基于GOCO05S模型的2～200阶系数计算的参考重力异常，得到残余重力异常；根据残余重力异常的矩谐展开式得到矩谐系数；根据大地水准面差距的矩谐表达式求得残余大地水准面，恢复参考大地水准面得到最终的大地水准面；最后利用EGM2008模型2～2190阶位系数计算的大地水准面高作为真实值，比较基于矩谐分析的结果与真实值间差距，说明矩谐分析的可靠性与精度。

根据上文所述可知，矩谐分析受周期延拓边界效应的制约，并成为影响其精度的关键因素。本试验利用扩展参数法以减小边界效应的影响，图2中虚线范围为3°×3°中心区域。由于在不同截断阶数下大地水准面误差变化规律相同，因此以截断阶数30为例进行边界效应的研究，并绘制误差分布图。

选取扩展参数为截断阶数30、扩展参数为50 km 绘制大地水准面对比图，如图3、图4所示。

2.2 基于航空重力数据的算例与分析

由于航空重力测量获取的是飞行高度处的重力扰动值，但所需值为大地水准面上的扰动值，而将飞行高度处的数值向下延拓至大地水准面的过程中，噪声会被放大，重力场信息因此失真，这是一个病态求解的过程。但航空重力测量对于获取高频重力场信息有不可替代的作用，因此如何能在延拓过程中将重力场信息高度还原，获得稳定的重力场解是航空重力数据处理的难点问题。

为验证矩谐分析的全面性、可靠性，本文基于航空重力数据设计如下试验。获取纬度为[29°，34°]、经度为[111°，116°]围成的5°×5°区域，Dx=504.37 km，Dy=601.34 km，格网间隔为2.5′×2.5′，共计14 400个计算点。利用EGM2008重力位模型的2～2190阶系数分别计算飞行高度为4 km的重力扰动的模拟观测值，并加入数学期望为0、标准差为2 mGal的高斯白噪声；移去基于GOCO05S模型的2～200阶系数计算的参考重力扰动，得到残余重力扰动；根据残余重力扰动的矩谐展开式得到矩谐系数；根据大地水准面差距的矩谐表达式和重力扰动的矩谐展开式求得残余大地水准面；最后恢复参考大地水准面及参考重力扰动得到最终的大地水准面；最后利用EGM2008模型2～2190阶位系数计算的大地水准面作为真实值，比较基于矩谐分析的计算结果与真实值间的差距。

图2 截断阶数30时取不同扩展参数矩谐分析所得大地水准面的误差

图3 真实大地水准面

图4 由重力异常值所得大地水准面

与重力异常处理方法相同，航空重力值也选用扩展参数法以减小边界效应的影响。由于不同截断阶数时的大地水准面误差变化规律相同，因此选取截断阶数30为例绘制误差分布图，如图5所示。

图5 截断阶数30时取不同扩展参数矩谐分析所得大地水准面的误差

选取扩展参数为截断阶数30、扩展参数为50 km 绘制大地水准面对比图，如图6、图7所示。

图6 真实大地水准面值

图7 由重力扰动值所得大地水准面

3 分析与结论

如图2、图5所示，当截断阶数固定时，针对每一截断阶数，分别选取0、30、50、70 km的扩展参数计算其大地水准面的误差分布。由于各截断阶数的变化规律一致，因此选取截断阶数30为案例进行分析。表1、表2分别给出了不同研究范围内取不同扩展参数时矩谐分析所得的大地水准面的误差统计信息。

表15°×5°范围时取不同扩展参数时矩谐分析所得大地水准面的误差统计m

范围扩展参数Δ最小值最大值平均值均方差5°×5°0-0.3380.311-0.0150.06730-0.3590.299-0.0150.06750-0.5460.273-0.0150.06370-0.5540.269-0.0140.063

对于全图5°×5°区域范围，在4个角点处误差出现极大值，边界效应明显，对于本研究区域尤以左上角部分最为严重，误差最大值达到3 dm；在扩展参数为0、30 km时均方误差基本不变；当扩展参数不断增大时，整体误差逐步降低，当扩展参数为50 km 时，均方误差产生小型跳变，最大值降低到2.7 cm，均方误差为6.3 cm；继续增大扩展参数，达到70 km时，并不能减小均方误差，且会增长计算时间，影响计算效率。

表23°×3°范围时取不同扩展参数时矩谐分析所得大地水准面的误差统计m

范围扩展参数Δ最小值最大值平均值均方差3°×3°0-0.1570.1190.0150.04830-0.1510.1160.0150.04750-0.1510.1150.0150.04270-0.1510.1150.0140.042

可以明显看出图示中心3°×3°区域在没有加入扩展参数时相较于周边区域误差较小，但依然受到长波系统误差的影响；当加入扩展参数时，随着扩展参数的增大，中心3°区域大地水准面的误差逐步明显降低，与5°区域变化趋势一致；在扩展参数增大到50 km时误差值发生一个小型跳变，降低到4.2 cm；继续增大扩展参数(Δ=70 km)均方误差没有明显变化，并不能提高大地水准面的精度。

按同样的方式对航空重力数据进行统计，误差统计见表3、表4。在不同范围的不同扩展参数下其变化趋势与重力异常数据一致。

综上所述，边界效用确实影响矩谐分析的精度，但采用扩展参数方法后，精度上升明显。扩展参数为50 km时矩谐分析得到的大地水准面的精度最高，扩展参数过小(如0、30 km)对精度的提高无法达到最佳，当扩展参数过大时(如70 km)，并不能继续提高精度，相反会大大增加解算时间以及复杂程度；因此实际计算中宜选取扩展参数为50 km进行解算。

表35°×5°范围时取不同扩展参数时矩谐分析所得大地水准面的误差统计m

范围扩展参数Δ最小值最大值平均值均方差5°×5°0-0.3340.273-0.0140.06830-0.3540.300-0.0150.06650-0.5460.275-0.0150.06370-0.5490.271-0.0140.062

表43°×3°范围时取不同扩展参数时矩谐分析所得大地水准面的误差统计m

范围扩展参数Δ最小值最大值平均值均方差3°×3°0-0.1580.104-0.0220.04630-0.1560.107-0.0200.04450-0.1530.109-0.0180.04370-0.1540.110-0.0170.043

限制矩谐分析精度的另一大难题是截断阶数的选取。为获得针对矩谐分析的适当的截断阶数，避免选取截断阶数较大时产生“龙格”现象，截断阶数较小时无法获取有效的数据值，本文针对扩展参数为50的不同截断阶数的矩谐展开结果进行分析，误差分布如图8所示。数值统计分析见表5。

表5取不同展开阶数时矩谐分析所得大地水准面的误差统计m

范围截断阶数(N=M)最小值最大值平均值均方差3°×3°20-0.1560.119-0.0190.05925-0.1520.117-0.0160.04830-0.1530.119-0.0180.04335-0.1510.115-0.0160.04340-0.1520.116-0.0160.043

可以看出，当截断阶数由20增加至40时误差逐步减小，均方误差值由5.9 cm减小至4.3 cm，误差减小了27%；其中尤以截断阶数从20增加至25时变化最为明显，均方误差减少1.1 cm；当截断阶数达到30时，误差变化趋于稳定，继续增大截断阶数并不能减小误差。相较于截断到35阶时，截断阶数为30阶时误差分布更为均匀，最大值与最小值差距较小。而截断阶数达到40时虽然均方误差值不变，但计算时间大大增加，不利于进行大量数据点的解算。故在实际计算中可选取截断阶数30进行矩谐分析。

图8 选取扩展参数为50 km时不同截断阶数下矩谐分析所得大地水准面的误差

4 结 语

针对球谐分析及球冠谐分析在求解小区域大地水准面时的解算难点，本文利用更稳当的矩谐分析求解大地水准面。为验证理论的可靠有效性，设计完成了基于EGM2008、基于陆地重力数据和航空重力数据的试验。试验结果表明，利用矩谐分析求解大地水准面，基于陆地观测数据精度可达到4 cm，基于航空重力观测数据精度可达到4 cm。相比于应用广泛的stokes理论求解大地测量边值问题所得结果，精度提高明显。因此，矩谐分析是确定区域大地水准面的理想方法，值得应用推广。