-
Stokes理论和Molodensky理论是以重力异常为基本物理输入量、使用正高/正常高高程系统的第三大地边值问题的两大理论体系。在第三边值问题中，重力异常为大地水准面上的重力与参考椭球面上的正常重力之差(Stokes理论)或地表重力与近似地形面上的正常重力之差(Molodensky理论)。1969年，英国学者Hotine提出了Hotine积分公式，提出了以扰动重力为边值条件的第二大地边值问题。由于扰动重力是同一点的重力与正常重力之差，因此扰动重力较重力异常具有更明确的物理意义。
过去由于无法获得大地高数据，因此第二大地边值理论的发展和应用受到限制。现阶段，SRTM(shuttle radar topography mission)数据在确定高精度的(似)大地水准面中获得越来越广泛的应用[1-4]。由于SRTM数据以EGM96(earth gravitational model 96)大地水准面起伏为高程基准，因此通过EGM96重力场模型可方便地将SRTM数据转化为大地高。与此同时，卫星重力技术的发展使重力场模型的精度不断提高，因而可借助高精度的重力场模型将离散重力点的正常高数据也转化为大地高。SRTM数据和高精度重力场模型的出现，使得大地高数据能够间接地获得，因此也初步具备了研究与应用第二大地边值问题的条件。
在第二边值问题的理论层面，Moritz最先提出GPS边值问题及其一般解式。文献[5]分析了重力异常和扰动重力的特点，指出以扰动重力作为边值条件更加简单、严密。文献[6]详细论述了以地心参考椭球面为边界面的第二大地边值问题的理论意义和实际意义，进而引出了该问题的原理和内涵。文献[7-12]给出了GPS边值问题的定义、数学表达，并通过类似于Molodensky理论中的分析方法，得到了第二边值问题以地表为参考面的实用公式。在上述研究成果的基础上，本文研究将第三边值问题中发展和应用均比较成熟的Molodensky理论方法应用于第二边值问题。为此，本文推导了扰动重力梯度与重力异常梯度的关系、Hotine算子[13]与梯度算子的关系，给出了针对第二边值问题的Molodensky算法，并通过实验验证在第二边值问题中应用Molodensky算法的可行性。
HTML
-
传统第三边值问题中，用于解算似大地水准面的Molodensky算法公式为：
式中，dΩ′为球面微分单元；ζA、γA、hA、ΔgA分别为计算点A处的高程异常、正常重力、高程、重力异常；R为地球平均半径；ψ为积分点和计算点间的角距；Δg为Stokes积分中流动积分点的重力异常；C为局部地形改正；S(ψ)为Stokes积分的核函数，定义为：
其中，t=sin(ψ/2)；Pn(cosψ)为勒让德函数；δhA2是hA2(计算点高程的平方)的零阶项与一阶项之和，通常采用下式计算：
其中，φA、λA分别为计算点的球面纬度和经度。
-
类似于Stokes积分的Moritz解析延拓解的形式，李斐等[7-9]给出了以地面为参考面的Hotine积分的实用公式：
式中，δg为Hotine积分中流动积分点的扰动重力；h为Hotine积分中流动积分点的高程；H(ψ)为Hotine积分的核函数；L表示梯度算子；L(δg)为地面扰动重力沿地心向径的梯度，其计算公式为：
式中，δg′为梯度算子中流动积分点的扰动重力；r为计算点的地心向径；R为地球平均半径；ℓ0为计算点与积分点的平面距离，计算公式为ℓ0=2Rsin(ψ/2)。
Hotine积分核函数H(ψ)在实际应用中有多种形式，其中包含零阶、一阶项的标准形式为：
式(6)所示的标准Hotine核函数的Legendre级数包含了零阶、一阶项。由于通常认为扰动位不含零阶、一阶项，因此标准Hotine核函数也可以是如下不含零阶、一阶项的形式：
此外还可将Hotine核函数除去与参考的重力场模型阶数一致的Legendre级数项，可称为截断Hotine核函数：
式中，L为移去Legendre级数的最高阶数。
-
下面给出第二边值问题中的Molodensky算法的推导过程。
球近似下，重力异常和扰动重力与扰动位T的关系分别为：
由式(9)可知Δg=δg-2T/R，进而可推得扰动重力梯度与重力异常梯度存在下列关系式：
式中，2δg/R量级很小，可忽略(梯度算子式(5)也包含2δg/R项，该项在梯度算子中通常也被忽略)，即有L(δg)≈L(Δg)。在传统的Molodensky理论中，通过假设重力异常和高程呈线性关系，得到如下公式[14]：
因此有hL(δg)≈-C+πGρL(h2)，将其代入式(4)得：
下面推导Hotine算子与梯度算子的关系式。对于球面上的连续函数f，设g为其梯度，则有：
gn为g的n阶球谐项，根据球谐展开原理，有[14]:
式中，fn为f的n阶球谐项。根据球函数的展开理论，L(f)的Hotine积分可进行如下推导：
式中，H表示包含零阶、一阶项的标准Hotine算子(见式(6))。式(15)的结果说明，当Hotine积分核为包含零阶、一阶项的标准形式时，Hotine算子与梯度算子严格满足关系式HL(f)=-γ-1I(f)(I表示单位算子，I(f)=f)。当Hotine核函数为式(7)中不含零阶、一阶项的形式时，根据式(15)可很容易得出：
式中，f0、f1分别为f的零阶和一阶球谐项。
在传统的Molodensky理论中推导了梯度算子与Stokes算子间的关系式为：
式中，S表示Stokes算子，该式等号右端最后一项量级很小，式(1)实际上省略了该项的计算。比较式(15)、(16)与式(17)可知，Hotine算子与梯度算子的关系式无需忽略微小项，因而比传统Molodensky理论中采用的Stokes算子与梯度算子的近似关系式SL(f)=-γ-1I(f-f0-f1)更为严格准确。
将式(15)代入式(12)，得到Hotine积分核包括零阶、一阶项时的Molodensky算法为：
将式(16)代入式(12)，即为Hotine积分核不包括零阶、一阶项时的Molodensky算法：
式中，δhA2的计算参见式(3)。类似于传统的Molodensky理论，式(19)的结果也是基于重力异常不包含零阶和一阶项这一前提。
-
SRTM数据以EGM96地球重力场模型计算的大地水准面为高程基准面。为了分别进行第二、第三边值实验，本实验将SRTM3数据通过以下方式分别转化为大地高和正常高：①将SRTM3正高数据与EGM96大地水准面高相加以得到大地高数据；②通过EGM96模型的高程异常——大地水准面高改正模型(改正模型的球谐系数[15]与EGM96模型同时发布)，将SRTM3正高数据转化为正常高数据。
实验区实测重力离散点共有70 379个，范围为105.5°~116.5°E，25.5°~34.5°N。将离散重力点格网化为106°~116°E、26°~34°N区域范围的1.5′×1.5′分辨率的扰动重力和重力异常数据。该区域地形高范围为0~3 370 m，平均515 m，图 1为区域地形分布图。
Figure 1. Topographic Map of the Experiment Zone
表 1对格网扰动重力和重力异常的最小值、最大值、平均值及均方根(root mean square, RMS)进行了统计，其分布如图 2、图 3所示。
统计量 min max mean RMS 扰动重力 -141.70 194.74 -28.21 41.97 重力异常 -134.69 194.58 -21.44 36.89 Table 1. Statistics of the Grid Gravity Disturbances and Gravity Anomalies/mGal
Figure 2. Distribution of the Gravity Disturbances
Figure 3. Distribution of the Gravity Anomalies
表 1中扰动重力和重力异常的统计结果相近，反映了该两类物理量具有相近的频谱特性。图 2扰动重力的分布与图 3重力异常的分布比较相似，说明本文推导的扰动重力梯度和重力异常梯度近似相等的关系式与实际情况是比较符合的。
下面以Eigen-6C4重力场模型的前360阶为参考重力场进行边值解算实验。生成的重力似大地水准面范围为108°~114°E、28°~32°N。该区域GPS水准检核点共计68个，其分布如图 4所示。
Figure 4. Distribution of the GPS Levelling Points
首先利用Moritz提出的梯度算子计算Hotine积分，Hotine积分核分别采用以下3种形式：①包含零阶和一阶项的标准形式；②不含零阶和一阶项的标准形式；③截断形式，Hotine积分半径为2°。计算的似大地水准面精度统计见表 2，其中包含RMS和标准差(standard deriation, SD)两项精度指标。
积分核 min max mean RMS SD ① -1.537 -0.169 -0.770 0.870 0.407 ② -1.505 -0.164 -0.749 0.848 0.400 ③ -0.120 0.296 0.100 0.130 0.083 Table 2. Accuracy of the Gravimetric Quasi-geoid by the Hotine Integral Using Gradient Operator/m
从表 2中可看出，两种形式的标准Hotine核函数计算的似大地水准面的精度较差。因此当使用梯度算子直接计算似大地水准面时，应使用截断形式的核函数，此时解算的似大地水准面精度为13.0 cm(均方根)、8.3 cm(标准差)。
下面分别利用传统第三边值问题中的Molodensky理论方法和本文推导的针对第二边值问题的Molodensky算法求解似大地水准面。其中局部地形改正的内区积分半径为0.5°，采用3″地形计算，外区积分半径为1°，采用9″地形计算。表 3统计了传统第三边值的Molodensky算法与第二边值问题的Molodensky算法计算的重力似大地水准面的精度，其中Hotine积分核采用两种形式：①包含零阶、一阶项；②不含零阶、一阶项，Hotine积分半径为2°。
积分核 min max mean RMS SD Stokes -0.217 0.108 -0.062 0.087 0.062 Hotine① -0.251 0.131 -0.042 0.086 0.076 Hotine② -0.210 0.080 -0.079 0.100 0.062 Table 3. Accuracy of the Gravimetric Quasi-geoid by Molodensky Algorithms/m
从表 3可以看出，第二边值问题的Molodensky算法与第三边值问题的Molodensky算法计算的重力似大地水准面达到同等精度水平，且高于使用梯度算子计算的重力似大地水准面精度，说明了Molodensky理论方法可以适用于第二边值问题。图 5所示为第二边值问题的Molodensky算法计算的重力似大地水准面。
Figure 5. The Gravimetric Quasi-geoid
-
本文在李斐等人给出的Hotine积分实用公式的基础上，借鉴第三边值问题中计算似大地水准面的Molodensky理论方法，给出了第二边值问题的Molodensky解算方法。结论如下：
1) 重力异常的梯度与扰动重力的梯度相差2δg/R，量级很小，因此可近似认为重力异常的梯度与扰动重力的梯度相等。
2) Stokes算子与梯度算子间的关系式SL(f)=-γ-1I(f-f0-f1)是种近似关系；包含零阶、一阶项的Hotine算子与梯度算子严格满足HL(f)=-γ-1I(f)；不包含零阶、一阶项的Hotine算子与梯度算子的关系式为HL(f)=-γ-1I(f-f0-f1)。
3) 通过实验得出第二边值问题的Molodensky算法计算的重力似大地水准面精度与传统第三边值问题中Molodensky算法计算的重力似大地水准面精度相当，且高于直接用梯度算子计算的重力似大地水准面的精度，这说明基于Molodensky理论对第二边值问题也是适用的。
本文借助于SRTM数据和高精度的重力场模型，证明了Molodensky算法在第二大地边值问题解算中的适用性。现阶段第二边值和第三边值已经能够以同样的精度水平计算似大地水准面，随着未来GPS对大地高测量精度的提高，第二边值必将具有更广阔的发展前景。
DownLoad:
