大地水准面的确定目前仍基于求解经典大地测量边值问题，最具有代表性的是Stokes理论和Molodesnky理论。为保证大地水准面外部扰动位为谐函数，Stokes理论要求大地水准面外部无质量。将大地水准面外部的地形质量移去，再以一定的方式对移去的地形质量进行补偿，是解算Stokes边值问题处理地形影响的基本方法。然而，无论用何种方式移去地形质量都将使大地水准面发生变化，产生间接影响。

Molodensky级数不必考虑地球质量的影响，无需作与地球密度相关的各种改正，它摆脱了Stokes方法需要已知地形密度分布的困难。但事实上地面及其外部扰动位又必然因地形起伏而变，因此其导出的级数解包含了地形影响项G_n(n=1，2，…)。计算G_n过于复杂耗时，其级数各项正负相间，收敛性目前也尚无定论，这就大大降低了Molodensky级数的实用效果。

近年来，发展了一种新型大地边值问题，即Stokes-Helmert边值问题，即在Stokes理论和计算模型的基础上，结合Helmert地形处理模式，同时兼容Molodensky理论，解算大地水准面外部扰动位的边值问题，又称Stokes双边值问题。Vaníček等于1995年介绍了Stokes-Helmert边值问题的基本原理。Sjöberg于2000年讨论了利用Stokes-Helmert法计算大地水准面过程中的地形影响，指出经典Stokes-Helmert算法已不能满足厘米级大地水准面精度的要求。Nahavandchi于2000年推导了基于Stokes-Helmert法计算大地水准面过程中直接地形影响的计算公式。罗志才于2003年比较了利用经典Stokes-Helmert方法与Sjöberg方法计算地形对大地水准面的影响，指出经典Stokes-Helmert方法存在明显的偏差，应采用更加严密的计算方法和模型。

利用Stokes-Helmert边值问题精化似大地水准面对地形影响处理技术提出了更高要求。目前大地水准面计算中地形影响的计算通常将参考面近似为平面，所有数据处理问题都可以在直角坐标系中进行理论推导和计算。然而，对于大区域乃至全球性的尺度问题，需要利用球坐标进行表达与处理，此时建立在直角坐标系中的数据处理方法已不再适用。因此，本文在已有工作的基础上，推导球坐标的Stokes-Helmert公式和严密地形影响计算公式，实现临沂厘米级大地水准面精化。首先基于确定地球外部扰动重力场的经典理论，借鉴国内外精密确定地表及其外部扰动引力位的严密解算理论，给出Stokes-Helmert边值问题的定义及Helmert重力场概念；然后根据实际算例，精确计算第2类Helmert凝聚法中由地形和凝聚层质量所产生的直接影响，以及地形及相应引力变化的间接影响，构建临沂高精度的似大地水准面模型。

1 Stokes-Helmert边值问题基本原理

Helmert重力场中的Stokes边值问题可表示为

式中，S_t和S_g分别表示地球自然表面和大地水准面；▽为梯度算子，▽²=Δ为Laplace二阶梯度算子；U为正常位，为残差地形位。用V^t和V^c分别表示地形质量和大地水准面S_g上的凝聚质量所产生的引力位，则

在式(1) 中, 式(2) 是一个非线性边界条件，为向量，| |表示▽W的模, 需要在球近似下线性化，即

式中，g为向量的纯量积，将式(3) 中的第3项略去，作泰勒级数展开，保留线性项，则有

略去二阶和高阶项，取球近似，略去地球扁率，并近似取γ与的方向重合，则有

式(6) 在S_g外部成立，用于S_t上的P点有

式中，r为地心距离，δA为对重力的直接地形影响。设点Q为点P沿正常重力线在正常椭球面S_e上的对应点，则γ_P可表示为

式中，N表示大地水准面高；H表示地形高，略去H³等高次项。引进空间改正，令

则式(7) 可写为

顾及，将式(9) 代入式(6)，并用P_g表示地面点P在大地水准面S_g上的对应点，则式(1) 中关于g的边界条件球近似线性化形式为

式中

Δg^h称为Helmert重力异常。点P和点Q分别为地面和正常椭球面上的点，为物质移动对正常重力的间接影响，是将大地水准面上的重力异常归化到Helmert重力场中Δg^h的改正项。δA为地形物质移动对g_P的直接影响。为确定T_Pg^h还需满足约束性边界U_Pg+T_Pg^h+δV_Pg=W₀。其中，U_Pe可由正常位模型表达为

式中，∂/∂h是沿正常重力线方向的导数；N为大地水准面高；∂γ/∂h近似为。略去第3项，且按规定U_Q=W₀，则N可表示为

式(13) 为Helmert重力场的布隆斯公式，其中是对Helmert重力场大地水准面高(经地形调整的大地水准面)δN=T^h/γ_Q的基本间接地形影响。从上述原理可以看出，基于Stokes-Helmert法计算大地水准面的关键问题是精确确定地形的直接影响和间接影响。

2 地形影响 2.1 直接地形影响

移去大地水准面外部地形质量的影响称为地形的直接影响。计算直接影响的经典公式为

式中，G为重力常数；ρ为地壳密度；l₀=2rsin ψ/2；H_p和H分别为地面点P和流动点的正常高，R=r。从严格意义上来讲，式(15) 只适用于远区积分，即l₀≥H。Sjöberg于1995年将直接影响球谐展开至H2，Nahavandchi等1998年将其扩展到H3，Sjöberg于2000年将直接影响写成曲面积分模式。

对于大区域的乃至全球尺度问题，发展椭球坐标系内高精度的地形改正方法尤其重要。本文基于国际通用的GRS80椭球采用Tesseroid单元体引力效应计算地形改正的方法。对于以地心为圆心，高度H=0、H=h₁的两个球面，其半径可近似表示为：r₁=Rradii+0，r₂=Rradii+h₁，Rradii表示所选等质圆的半径(这里选地球的平均半径)。则球面扇形柱体的重力位可表示为

式中，，为计算点P和积分点Q之间的距离；ψ是点P和点Q之间的夹角，即

由于椭球积分困难，扇形柱体的重力位不能解析计算，但是可通过数值计算得到近似解，式(15) 可以进行级数展开。

2.2 间接地形影响

恢复地形质量带来的影响称为间接影响，计算间接影响的经典计算公式为

式(17) 右端第一项称为主项，第二项称为次项。Sjöberg于1995年将间接影响球谐展开至H2，Nahavandchi等于1998年将其扩展到H3。然而由于公式需要对整个地球进行积分，这在计算中是不实际的。这里将给出间接改正的严密球近似公式。假设任意一点P的地形引力可写为

式中，为计算点和流动点之间的夹角。除以γ，可以得到地形对大地水准面的间接影响，即

式(19) 核函数f (H, t)可表示为

式中，

3 数据准备及离散重力数据格网化 3.1 数据准备

本文研究区域范围为34.3°N—36.3°N, 117.4°E—119.2°E，地形数据采用SRTM3模型。图 1为所用的SRTM3地形数据等值线图，高程最大值为1125 m，平均高程为151.7 m。

临沂市境内相对重力数据为995个，其点位分布如图 2圆点所示。为了检验地面重力数据的精度，采用该区18个长期观测的绝对重力观测数据(其点位分布如图 2中星号点所示)对相对重力测量数据进行检核。经分析，测区相对重力测量数据标准差为0.02 mGal。

目前应用最为广泛的全球地球重力场模型是EGM2008地球重力位系数模型。利用位系数计算参考重力场的公式为

式中，C_nm、S_nm为完全规格化重力位系数；GM为地心引力常数；a为参考椭球长半径；P_nmcos θ为完全规格化缔合Legendre函数；n_max为重力场模型的最大阶数2190阶。为了验证EGM2008模型在临沂地区的精度，本文采用实测地面重力数据对EGM2008模型进行外部测试。考虑到地面重力测量成果使用无潮汐基准, 模型重力场计算时也采用无潮汐基准模型。经过分析，模型数据与实测空间重力网格数据的标准差为5.66 mGal。

本区共有165个GPS/水准点，分布如图 3所示。其坐标系统采用2000国家大地坐标系，高程基准采用1985国家高程基准。任取其中142个点(圆点标记)作为高一级的控制点计算大地水准面，其余23个点(星号点标记)作为外部检验点，以检验求得的大地水准面的精度。

3.2 离散重力数据格网化

重力测量的结果是一些分布不规则的离散点重力值。为了便于使用快速傅里叶变换(FFT)等科学计算方法，离散重力数据需要进行格网化。然而由于空间重力异常的变化较大，直接对空间重力异常进行格网化将产生较大的误差，因此，通常在进行格网化之前对重力异常进行归算。本文重力归算采用布格重力异常进行格网化，具体计算步骤如图 4所示。

采用上述格网化步骤对地面离散数据进行格网化，格网化范围为34.3°N—36.3°N, 117.4°E—119.2°E。临沂市境内采用地面观测得到的空间异常，临沂市外空白区域采用EGM2008模型重力异常进行填充。格网化过程中的地形数据采用SRTM3模型。按图 4所述的方法，将离散重力数据格网化成2.5′×2.5′分辨率的网格点。图 5给出了格网化后的空间异常等值线图，其最大值为50.74 mGal，最小值为-6.61 mGal，平均值为9.36 mGal。

4 临沂大地水准面模型精化 4.1 临沂重力似大地水准面计算

采用2.5′×2.5′格网空间重力异常作为输入数据，以EGM2008作为参考重力场，再利用第2类Helmert凝聚法计算大地水准面中的各类地形位及地形引力的影响，即牛顿地形质量引力位和凝聚层位间的残差地形位的间接影响，以及Helmert重力异常由地形质量引力位和凝聚层位所产生的引力影响。图 6给出了由Stokes-Helmert法计算的重力似大地水准面，单位为m。

为了评定重力似大地水准面的精度，采用171个GPS/水准资料与其进行比较。结果表明，重力似大地水准面与GPS/水准高程异常比较，差值最大值为11.43 cm，标准差为4.52 cm。同样，采用171个GPS/水准资料与EGM2008模型计算的大地水准面高进行比较，差值最大为21.04 cm，标准差为6.67 cm。由结果可见，由Stokes-Helmert方法计算的重力似大地水准面比模型结果相比，精度有显著提高。

4.2 GPS水准数据与重力似大地水准面高拟合

利用高精度GPS/水准数据对上述建立的重力似大地水准面进行控制拟合，求得最终的大地水准面模型。由GPS水准测量得到的高程异常及重力似大地水准面起伏之间的拟合模型通常可表达为

式中，H为GPS测量的大地测高；h为水准测量得到的正常高；N为重力似大地水准面起伏；x为基准不一致的拟合参数；A为系数矩阵，取决于拟合模型。本文计算中拟合模型Ax选用多面函数拟合模型，其高程异常函数可表示为

式中，α_i为待定系数；Q (x, y；x_i, y_i)为核函数；(x, y)为未测点坐标；(x_i, y_i)为已测点坐标。核函数一般可取

式中，δ称为平滑因子，为了得到较好的逼近效果，应作一定的试算后确定，一般可取作0；b为一个可供选择的非零实数，一般取1/2，即核函数为正双曲形。

本文采用的GPS/水准点共165个。任意选取142个均匀分布的GPS水准点作为控制点，其他23个点作为检验点。控制点和检核点的点位分布如图 3所示。采用曲面函数拟合法拟合大地水准面的外符合精度，最大值为3.18 cm，最小值为-2.79 cm，标准差为1.61 cm。如图 7所示。

5 结语

本文从讨论地形数据在大地水准面确定中的作用出发，分析比较了Stokes和Molodensky两种经典方法处理地形影响存在的问题，给出了Stokes-Helmert边值问题的定义及数学描述，以及地形直接影响、间接影响的理论表达式，介绍了解算Stokes-Helmert边值问题的详细步骤和流程。

以格网化后的重力场模型为基础，采用Stokes-Helmert边值问题，建立临沂地区2.5′×2.5′分辨率的重力似大地水准面模型。利用GPS水准高程异常对重力似大地水准面模型进行控制拟合，求得最终的大地水准面模型。经检核点检验可知，拟合外符合精度达到1.6 cm。