0 引言

地形改正是消除观测点附近高出或低于观测点水平面的地形质量对观测重力的影响而加的改正.对于局部地区的重力资料而言，可以将重力观测面近似为平面，所有数据处理问题都可以在直角坐标系中进行理论推导和计算，然而，对于区域性的，乃至全球尺度问题，重力观测面已不能近似为平面，而是一个球面问题(杜劲松等，2012).在空间域，一般将全球地形剖分成单元体，计算每个单元体产生的引力效应，然后累加求和即得到全球地形的重力效应(杜劲松等，2012).目前已知的解决办法主要有Tesseroid单元体泰勒展开式方法(Wild-Pfeiffer，2008)和球冠体积分方法(杜劲松等，2012)等，Tesseroid单元体泰勒展开式方法不存在严格的解析解，球冠体积分方法则需要对地形模型进行重构，降低了计算效率，而相应的球坐标系下的三维球谐谱公式也是展开到三次项的近似表达(Tsoulis，2001).

本文研究基于Tenzer等(2007)给出的极坐标下的牛顿引力公式，但该公式中存在除积分元素之外的变量，本文给出推导过程解决了这一问题，对公式进一步完善.选取中国及其邻域进行了实验，将其划分为1°×1°经纬网格，计算地形改正并绘制成图，并且和使用球谐谱方法得到的结果对比.

1 计算方法介绍 1.1 计算公式

全球地形的重力效应由每个单元体产生的重力效应累加而成，每个单元体在观测点产生的引力位公式为(MacMillan，1958)：

(1)

(1) 式存在严格的解析解，这里，δV(r, Ω)是Ω在观测点产生的引力位，Ω代表单元体，ψ是球心角，α是球面方位角，r′₁、r′₂是单元体半径上下限，α₁、α₂是球面方位角上下限，ψ₁、ψ₂是球心角上下限，r是观测点至球心的距离，取值6388 km，r′是单元体上的点至球心的距离，是欧式距离，G是万有引力常数，取值0.0000000000667 m³·kg^-1·s^-2，ρ是密度，取值2670 kg·m^-3(部分元素的几何意义如图 1所示).

单元体在观测点产生的重力异常在各个方向上的分量表达式为：

(1) 在纬度方向上：

(2)

(2) 在经度方向上：

(3)

(3) 在半径方向上：

(4)

式中相关表达式详见文献(Tenzer et al., 2007)，其中有一表达式为

(5)

因公式(5)中含有除积分元素之外的未知量λ′，需将其消去，公式推导如下：由球面三角公式(6)、(7)(海斯卡涅和莫里兹，1979)计算得出α、ψ，如图 1所示，可得：

(6)

(7)

将(6)式变形得：

(8)

将(8)式代入(7)式化简得：

(9)

(9) 式除以(8)式化简得：

(10)

再利用球面三角公式得：

(11)

将(11)式代入(10)式得：

(12)

将(12)式代入(5)式中消除掉λ′.

1.2 单元体角元素积分上下限的确定办法

本文中的地形改正算法是基于极坐标的，单元体半径的积分上下限是在球上直接确定的，地球半径(取值6378 km)加上高程即可，很容易求得.而角元素的积分上下限是在极坐标下由观测点和积分单元体的相对位置确定.如图 2所示，O为球心，A点为观测点在单位球面上的投影，球面四边形BCDE为单元体在单位球面上的投影，为了使积分结果更加准确，球面四边形BCDE以OB为轴旋转至BC′D′E′位置，则可以计算得：

(13)

(14)

然后，依据α₁、ψ₁在其所在弧线方向延长估算α₂、ψ₂则有:

(15)

(16)

这种角元素积分上下限的确定方法省去了地形模型重构这一步骤，提高了效率，适用于全球地形产生的重力效应的计算，这里，λ_B、ϕ_B分别是B点的经度和纬度，ϕ_C是C点的纬度.

1.3 计算流程

(1) 如图 3a所示，过观测点作平行于赤道面的平面，与经线所在的大圆面将全球划为①、②、③、④四个部分.

这里，规定①④在观测点处经度方向的重力分量为正，②③为负.

(2) 如图 3b所示，过观测点作大圆面，且此大圆面与赤道面的交线垂直于观测点经线所在的大圆面，将全球分为①②两个部分.

这里，规定①在观测点处纬度方向的重力分量为正，②为负.

(3) 依据每个单元体所在的部分，计算每个单元体在观测点处纬度方向的分量δg_ϕ(r, Ω)，经度方向的分量δg_λ(r, Ω)，半径方向上的分量δg_r(r, Ω).

(4) 观测点处纬度方向的分量总和记为∑_ϕ，经线方向的分量总和记为∑_λ，半径方向的总和记为∑_r，总值用∑表示，具体求法为：

(1) ∑_ϕ=∑δg_ϕ(r, Ω_i)，

(2) ∑_λ=∑δg_λ(r, Ω_i)，

(3) ∑_r=∑δg_r(r, Ω_i)，

则.

2 实验及结果分析

实验区域为纬度0°~60°N，经度60°~150°E的中国及其邻近区域，综合考虑，最终选择以1°×1°分辨率验证本算法的可行性和准确性.实验中，观测点高度为10 km，高程数据来自ETOPO1，ETOPO1是分辨率为1弧分的全球地形起伏模型，其包含了陆地地形和海洋水深的数据，是目前已知的分辨率较高的地形起伏数据，对其进行稀疏提取，得到1°×1°分辨率的地形起伏模型，本文中高程为负的取0计算.

图 4a是使用极坐标下球面方法计算的中国及其邻域的地形改正结果分布图，从图可见，相较其他区域，青藏高原地区的地形改正较大，与周边区域很明显的差异，向南越过喜马拉雅山脉过渡到印度大平原，阿尔金—祁连山脉向南与昆仑山脉将柴达木盆地包围、向北与天山山脉将塔里木盆地包围，形成很明显的分界线，与中国的重力梯级带相符合. 图 4b是使用球谐谱方法计算的中国及其邻域的地形改正结果分布图，图 4c是中国及其邻域的DEM图(这里为了直观将高程取反).图 4中a与c的相似度为99.7%，有很强的相关性，验证了本方法的正确性. 图 4a和图 4b中纬度5°~20°N，经度85°~105°E范围内的地形改正有较大的区别，图 4a中该区域的地形改正更加平滑，而图 4b中则略显粗糙，究其原因可能是周围陆地地形在两种方法下的影响不同.

图 5是中国及其邻域使用本文算法与球谐谱方法得到的地形改正之差的分布图.从图中可以看出，在局部地区差值的大小和地形的起伏呈现较强的相关性，这为提高本文算法地形改正计算精度提供了思路.

3 结语

本文将基于极坐标的牛顿引力公式用于地形改正的计算，不同于以往局部和区域性的地形改正方法，本文算法有效地避免了大区域乃至全球范围计算时将球面看成平面带来的误差影响.从中国及其邻域地形改正结果图分析可得，准确地反映了地形的起伏变化，是一种可行的算法.由于本算法不需要进行坐标旋转，适用于重力三维反演.以后，将进行进一步的研究，将其应用到更多方面.

致谢 感谢审稿专家和编辑部的热情支持和帮助！