近区地形直接与间接影响的棱柱模型算法

马 健，魏子卿

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

应用Stokes/Hotine边值理论解算(似)大地水准面时要求边界面外不存在地形质量[1-2]，对此可通过Helmert第二压缩法将边界面外的地形压缩成覆盖在边界面上的质量薄层，由此产生了Stokes-Helmert大地边值理论[3]并逐步应用到实践中[4-7]。按照Helmert第二压缩法进行地形压缩时，地形质量的移动对重力产生直接影响，对大地水准面起伏/高程异常产生间接影响。在基于Helmert第二压缩法的边值问题中，直接影响将重力数据从真实空间转化到Helmert空间，而间接影响将Stokes/Hotine积分得到的Helmert高程异常/大地水准面高从Helmert空间恢复到真实空间。由于直接、间接影响直接参与(似)大地水准面解算，因此提高直接、间接影响的计算精度对于生成高精度的(似)大地水准面有着重要意义。

目前国外对地形直接、间接影响算法的研究相对较多，其中，文献[8—9]给出了地形直接、间接影响的球谐表达式，文献[10—11]提出综合严格积分方法和球谐方法计算地形影响的算法。国内在解算大地水准面中地形的研究主要为局部地形改正[12-16]、层间改正算法[17]等方面，对Helmert压缩理论中地形直接、间接影响研究相对较少。随着国内专家学者对基于Helmert第二压缩法的边值解算的应用与研究[4，18]的深化，对地形直接、间接影响展开研究是非常有必要的。

为了兼顾计算精度和计算效率，近区地形影响通常采用积分算法计算，远区则使用Molodensky截断理论或FFT技术结合高程球谐模型计算[19-21]。本文主要研究近区地形影响的计算。目前使用的地形直接、间接影响算法仍是二重积分形式，该算法在中央区(与计算点位于同一地心向径的积分点所在的地形格网，也可称为中央网格)存在奇异性。实践中通常采用圆柱近似或借助球面布格层等方法单独计算中央区地形影响[22]，增大了计算的复杂性。

目前地形质量产生的引力位和引力的棱柱模型算法已经比较成熟[23-24]，而局部地形改正的棱柱算法在实践中已经获得广泛应用[14]，但尚未出现关于地形直接、间接影响的棱柱模型算法的研究。根据定义，地形改正是计算点周围地形起伏对计算点产生的引力作用，而地形直接、间接影响是边界面外地形压缩到边界面上产生的引力和引力位，其本质均是地形质量(及其压缩形式的变体)产生的引力位的函数[25]，因此理论上应可得到直接、间接影响的棱柱模型算法。鉴于此，本文通过一系列推导给出了近区地形直接、间接影响的棱柱模型公式。为了避免棱柱模型算法本身存在平面近似误差，可进一步根据文献[24]提出的方法方便地转化为顾及地球曲率的棱柱模型算法。地形影响的棱柱模型算法不仅提高了地形影响的计算精度，而且在中央区不存在奇异性，无需单独计算中央区地形影响，从而使地形影响传统的“中央区影响+近区影响+远区影响”计算模式转变为“近区影响+远区影响”计算模式。

1 地形影响的传统积分公式

使用Helmert第二压缩法进行地形压缩时，需依据一定的压缩准则。本文采用质量守恒的压缩准则，即在球近似下压缩面密度σ与积分网格的地形高度h的关系为[26]

(1)

式中，ρ为地形密度；R为地球平均半径。将大地水准面外真实地形产生的引力位表示为Vt，压缩地形产生的引力位表示为Vc，其算法分别为

(2)

直接影响的积分算法公式为[8,19,27]

(3)

式中，μ=Gρ，计算时ρ通常取地壳平均密度。

地形压缩产生的残余位(真实地形产生的引力位与压缩地形产生的引力位之差)的积分算法公式为[9,19,22,27]

(4)

根据Bruns公式，当计算点位于地面上时，得到地形压缩对高程异常的间接影响δζ为[9,19,27]

δζ=(Vt-Vc)/γ

(5)

当计算点位于边界面上时(即hP=0)，得到地形压缩对大地水准面高的间接影响δN为[9,19,27]

δN=(Vt-Vc)/γ0

(6)

式中，γ为地面点正常高处的正常重力；γ0为参考椭球面上的正常重力。

2 近区地形影响的棱柱模型算法推导

鉴于目前地形影响积分公式还是二重积分的形式，以下对地形影响的棱柱模型算法展开推导。通常地形是以网格形式存储的，在空间直角坐标系下，任一地形网格Π与计算点P的相对位置如图1所示。

图中P为计算点，坐标系原点O为计算点在边界面(大地水准面或参考椭球面)上的投影，因此有xP=0、yP=0。x轴指向北、y轴指向东，x1、x2分别为棱柱南、北边界的x坐标值，y1、y2分别为棱柱西、东边界的y坐标值。hP为计算点高程，h为地形网格(实际计算中代表流动积分点)的高程。

图1 地形网格与计算点相对位置示意图Fig.1 The schematic diagram of a topographic grid away from the computing point

2.1 直接影响

根据直接影响的定义，直角坐标系下直接影响的公式为

(7)

式(7)可化为下列二重积分

(8)

已知积分方程式

(9)

式中，C为不定积分式中的任意常数。根据式(9)可得式(8)的一重积分形式为

(10)

进一步运用积分式

(11)

可将式(10)转化为

(12)

若f(x)为变量x的函数，g(y)为变量y的函数，根据式(13)

(13)

可对式(12)进一步简化得到直接影响的棱柱模型算法公式

(14)

2.2 间接影响

由于间接影响可由残余位通过Bruns公式得到，因此这里首先对残余位的棱柱模型算法进行推导。直角坐标系下残余地形位的基本公式为

(15)

运用式(9)可将式(15)转化为下列二重积分

(16)

根据式(9)、式(11)，将式(16)转化为一重积分

(17)

由于当ab<1时，反正切函数具有下列性质

(18)

从而可得出

(19)

另外，根据分部积分法可以得到下列不定积分成立

(20)

根据式(19)、式(20)及式(13)，由式(17)得到残余地形位的棱柱模型公式为

(21)

将式(21)分别代入式(5)、式(6)即为地形压缩对高程异常、大地水准面高的间接影响。

地球是一个曲面，因此地形直接、间接影响的棱柱模型算法将近区范围视为平面，存在平面近似误差，当近区范围较大时，需顾及地球曲率的影响。根据文献[24]提出的算法，当顾及地球曲率的影响时需对计算点高程进行如下改化

(22)

下面从理论上对传统积分算法存在的误差进行说明。由于角距ψ是纬度θ、经度λ的函数，因此可将式(3)、式(4)的直接、间接影响积分算法简化为下列形式

(23)

当求取二重积分式(23)的单个网格的离散积分结果时，计算中采用的近似公式为

(24)

式中，Δθ=θ2-θ1；Δλ=λ2-λ1。这一算法实际上是以网格中心点处积分核的大小作为网格积分核的平均值。对于一元积分

(25)

令f(θ)=sinθ，上下界θ2=π/2、θ1=0，则式(25)左端等于1/2，右端等于π/4，显然是不成立的。当θ1、θ2差值减小时，式(25)两端的差异(误差)也将随之缩小。需要特别说明的是，当f为kθ/cosθ(k为常数，此时积分核与自变量为线性关系)时，式(25)是恒成立的。上述对一元积分式(25)的分析同样适用于二重积分式(24)。

传统积分算法是二重积分形式，由于地形影响积分核函数与几何位置并非简单的线性关系。通过以上分析可知，传统积分算法计算地形影响时以网格中心点处积分核的大小作为积分单元积分核的平均值，在一定程度上引入了近似误差。由于棱柱模型公式不存在此项误差，因而可提高地形影响的计算精度。此外，棱柱模型公式在中央区不存在奇异性，实现了中央网格与其他网格算法的统一。

需要说明的是，同地形改正的棱柱模型算法一样，当x1、x2、y1、y2中任意一个值接近0时，棱柱算法存在奇异性，但当满足这一条件时传统积分算法是非奇异的，因此对于满足这一条件的个别网格，计算地形影响时使用传统积分算法即可。

3 试 验

为了定量说明本文推导的棱柱模型算法的精度，本文在我国中部地区进行了地形影响试验，区域范围为104°E～108°E，32°N～36°N，海拔最高3941 m，最低415 m，平均1562 m。由于顾及地球曲率的棱柱模型算法既避免了未顾及地球曲率的棱柱算法的平面近似误差，又克服了传统积分算法以格网中心点处积分核大小作为格网平均值存在的近似误差，因此该算法计算的地形影响更加准确、可靠。本文试验中将顾及地球曲率的棱柱模型计算的地形影响视为地形影响的准确值。

利用顾及地球曲率的棱柱算法分别计算了2.5′×2.5′分辨率下0.5°×0.5°、1°×1°、3°×3°共3种积分半径的直接、间接影响，其中间接影响为地形压缩对高程异常的间接影响(下同)。

表1中3种积分半径下直接、间接影响各自的差异较小，说明近区地形影响主要受计算点附近地形的影响，而积分半径对其影响不明显(1 mGal=1×10-5m/s2)。为了说明不顾及地球曲率的棱柱算法的平面近似误差和传统积分算法以网格中点积分核大小作为网格积分核平均值的近似误差的量级，将两种算法计算的地形影响分别与表1中顾及地球曲率的棱柱模型计算的地形影响作差。表2对直接影响的差异进行了统计。

表1顾及曲率的棱柱模型计算的不同积分半径的地形影响

Tab.1Topographiceffectscalculatedbyprismmodelsconsideringcurvatureunderdifferentintegralradius

地形影响算法最小值最大值平均值RMS直接影响/(mGal)0.5°×0.5°-64.68061.066-2.06910.7671°×1°-63.96264.690-0.67210.5483°×3°-63.45266.8960.43910.545间接影响/m0.5°×0.5°0.0080.7510.1520.1861°×1°0.0070.7560.1510.1863°×3°0.0010.7410.1410.177

表2两种算法与顾及曲率的棱柱模型计算的直接影响差异

Tab.2 The direct effect differences between the two algorithms and the prism model considering curvature mGal

由于顾及曲率的棱柱算法计算的地形影响近似作为地形影响的真实值，因此表2中未顾及曲率与顾及曲率的棱柱算法差异反映了未顾及曲率的棱柱算法的平面近似误差，而积分算法与顾及曲率的棱柱算法差异反映了积分算法以网格中点的积分核大小作为网格平均值的计算模式引入的近似误差。从表2可以看出，两种误差随积分半径变化不明显。3种积分半径下，未顾及地球曲率的直接影响棱柱算法优于直接影响的传统积分算法：未顾及曲率的棱柱模型计算的直接影响的平面近似误差很小，仅为0.06 mGal，而传统积分算法计算的直接影响误差达3.36 mGal。表3统计了两种算法计算的间接影响与顾及曲率的棱柱模型计算的间接影响之间的差异。

表3两种算法与顾及曲率的棱柱模型计算的间接影响差异

Tab.3 The indirect effect differences between the two algorithms and the prism model considering curvature cm

由于直接、间接影响具有不同的频谱特性，两种算法的直接、间接影响误差也呈现不同的特点。表3说明了未顾及曲率的棱柱算法计算的间接影响误差随积分半径的增大而增大：0.5°积分半径下误差仅为0.08 cm，比积分算法的间接影响误差约优一个数量级；1°积分半径下误差达到0.17 cm；当积分半径为3°时，未顾及曲率的棱柱模型算法计算的间接影响误差明显增大，已超过积分算法计算的间接影响误差，说明当间接影响的积分半径较大时，应使用顾及地球曲率的棱柱算法。

为了反映地形影响在不同分辨率下的情况，表4统计了顾及地球曲率的棱柱算法计算的30″×30″与3″×3″分辨率的直接、间接影响，积分半径为1°。

表4顾及曲率的棱柱算法计算的不同分辨率的地形影响

Tab.4Thetopographiceffectsbytheprismalgorithmsconsideringcurvaturewithdifferentresolutions

地形影响分辨率最小值最大值平均值RMS直接影响/(mGal)30″×30″-54.34290.3658.24014.5183″×3″-53.12892.03210.24415.694间接影响/m30″×30″0.0070.7640.1510.1873″×3″0.0070.7660.1510.187

根据Nyquist采样定律，2.5′×2.5′网格含有0～4320阶地形信息，30″×30″网格包含0～21 600频段的地形信息，而3″×3″网格包含0～216 000频段的地形信息。相同积分半径(1°)下，表1、表4中30″×30″分辨率的直接影响与2.5′×2.5′分辨率的直接影响差异较为明显，说明直接影响在4320～21 600频率波段的信息含量较多。表1、表4中1°半径的3种分辨率的间接影响差异很小，说明间接影响在超过4320阶的频率区间的信息较少。为了说明30″×30″、3″×3″分辨率下传统积分算法计算的地形影响的误差情况，将积分算法与顾及曲率的棱柱算法计算的地形影响差异统计于表5。

表5传统积分算法与顾及曲率的棱柱算法计算的地形影响差异

Tab.5Thedifferencesofthetopographiceffectsbetweenthetraditionalintegrationalgorithmsandtheprismalgorithms

地形影响分辨率最小差值最大差值平均差值RMS直接影响/(mGal)30″×30″-0.5685.3043.8374.0403″×3″-0.158-0.015-0.1060.110间接影响/cm30″×30″-0.1870.6860.0480.1613″×3″-0.964-0.060-0.2910.341

从表5可以看出，30″×30″分辨率的直接影响误差为4.04 mGal，3″×3″误差仅为0.11 mGal，说明分辨率为3″×3″时传统积分算法计算的直接影响误差基本可以忽略。但对于间接影响，试验区30″×30″分辨率的误差为0.16 cm，而3″×3″分辨率的误差为0.34 cm，此时相应的最大误差分别为0.68 cm、0.96 cm(绝对值)，因此当构建高精度的(似)大地水准面时应使用顾及曲率的棱柱模型算法计算间接影响。

4 结 论

使用Helmert第二压缩法对边界面外地形进行压缩时，对重力数据和(似)大地水准面分别产生直接、间接影响。近区地形影响的传统算法是积分算法，由于传统积分算法为二重积分形式，计算时以网格中点积分核的大小作为网格积分核平均值存在一定的近似误差。此外，传统积分算法在中央区存在奇异性，需单独计算中央网格地形影响，增加了地形影响计算的复杂度。为此，本文推导了近区地形直接、间接影响的棱柱模型算法，为进一步避免棱柱模型存在的平面近似误差，可使用顾及地球曲率的棱柱模型算法计算地形影响。由于棱柱模型算法在中央区不存在奇异性，因此地形影响的计算模式由传统的“中央区影响+近区影响+远区影响”转变为“近区影响+远区影响”，实现了算法的简化。

当试验区地形影响取1°积分半径时，对于2.5′×2.5′分辨率的地形，未顾及曲率的棱柱模型计算的直接影响的平面近似误差仅为0.06 mGal，而传统积分算法的误差达3.36 mGal，此时两种算法计算的间接影响误差分别为0.17 cm与0.74 cm；当使用30″×30″与3″×3″分辨率地形时，传统积分算法计算的直接影响误差分别为4.04 mGal与0.11 mGal，间接影响误差分别为0.16 cm与0.34 cm。综上所述，对于直接影响，未顾及曲率的棱柱模型的平面近似误差很小，在地形分辨率较高时传统积分算法的误差也是可忽略的，但对于间接影响，在高精度要求的大地水准面计算中应尽量使用顾及地球曲率的棱柱模型算法。