重力局部地形改正量的计算方法比较

欧阳明达，周 巍，张敏利

(1.地理信息工程国家重点实验室，陕西 西安 710054;2.西安测绘总站，陕西 西安 710054)

地形改正在大地水准面计算、重力勘探、地壳密度结构研究等方面有重要应用[1-2]。现有大量文献对地形改正计算提出了一系列解决方案，且通常将地形分为内区和外区。传统的平面计算方案，外区被视为无限布格平板，内区被近似为一有限范围且以计算点为中心的平面[3-4]。FFT方法的出现提高了地形改正值的计算效率，但需要提供规则化格网形式的地形高输入数据[5-7]。随着卫星重力观测技术的发展，以及全球性大尺度海量观测数据的获取，学者们开始了对球坐标系下重力场元的研究，该方法将地表近似为贴合的球面，外区为球面布格层，其量级约为平面无限布格层的两倍，内区则采用球面积分方法进行计算[8-10]。

本文重点关注局部重力地形改正，即内区地形改正值的获取，给出了平面积分方法、平面FFT方法、球面积分方法的计算公式，以西部某山区为例开展试算，对其结果进行了比较，并得出了一些有益的结论。

1 平面积分

根据定义，不规则地形起伏相对于计算点P(x,y,z)产生的引力位为

(1)

式中，G为牛顿引力常数，ρ为地壳密度，σ代表积分区域，x、y、z、h为流动点的坐标及高程，xp、yp、zp、hp为计算点的坐标及高程，式(1)的负垂直方向分量即为重力的地形改正：

(2)

假设DTM以N×M的规则格网数据表示，每个网格内的地形用一个质量均匀分布而高程不变的棱柱体代替，则得到用质量棱柱地形模型表示的地形改正式：

(3)

令

(4)

(5)

上式右侧可改写为[1+x]-1/2的形式，其中，x=[(hp-h)/l]2，级数展开为

(6)

从式(6)可知，地形改正计算时，流动点与计算点连线的坡度不能超过45°，很明显，这种情况非常少见，即便在山区也与大多数地形地貌严重不符。将式(6)引入式(5)，则积分式可改写为

(7)

即

δg≈C1+C2+C3.

(8)

本文取其1阶近似，得到

(9)

对其进行FFT变换，改写成谱形式：

(10)

采用FFT方法计算时，在计算点处，核函数出现奇异，通常采用引入核函数项增加常数因子(α)的方法解决这一问题，即使用新的水平距离函数代替原来的函数，具体实现过程参考文献[6]。

2 球面积分

平面积分忽视了地球曲率影响，对较远区来说存在误差。将积分区域视为球近似，则重力的局部地形改正可表示为

(11)

式中，λ′、φ′、r′为流动点的球坐标，L-1为牛顿积分核函数(即流动点到计算点空间距离L的倒数)，且

(12)

cosψ=cosφcosφ′+sinφsinφ′cos(λ-λ′).

(13)

r(3t2-1)ln|r′-rt+L|+C.

(14)

其中，t=cosψ，C为积分常数。

计算地形改正的奇异积分问题，采用“梯度法”处理[11]。以计算点为中心，取一足够小的球冠，将其视为平面，半径为Δφ0=Δλ0，表达式为

(15)

其中，

(16)

[h(φp,λp+Δλ0)-h(φp,λp-Δλ0)].

(17)

3 算例与结果分析

3.1 数据准备与方案设计

选择我国西部山区 35°～38°N，95°～98°E范围进行计算，高程数据来自于SRTM模型，经平滑得到分辨率为30″的高精度地形数据，如图1所示，横坐标、纵坐标为经纬度范围，数值单位为m。由于计算需要，地形模型需向外分别扩充2°范围。试验区地形高程的最大值为5 592 m，最小值为2 675 m，平均值为3 695 m。

图1 地形模型

3.2 结果分析与比较

根据相关文献结论，选择积分半径为100 km[12]。图2表示出了采用3种方案得到的局部重力地形改正值。表1给出了其相关统计信息，可以看出，地形改正值结果多为正值，不同方法具有显著差异；计算点周边地形起伏大，其改正值较大；中部盆地受较远地区复杂地形的影响较小，改正值较小，多位于0 mGal附近；局部地形改正量受山区、平原、盆地等的地貌影响较小(见图2)，数值单位为mGal，尽管南部为山区地貌，但其上地形起伏较为平缓，因而地形改正量并不大。

图3给出了3种方案的计算结果对比，表2给出了其相关统计信息。从表2和图3可以看出，不同方案计算结果存在显著差异，平面方法与球面方法差异较大，差值结果的标准差分别为1.59 mGal和2.29 mGal，主要原因是两种方法对地形体的近似模型不同。平面积分方法和平面FFT方法的差值位于-1～6 mGal，标准差为0.78 mGal，不存在较大差别。中部地区差异较小，在0 mGal左右，南部和北部山区计算结果的差异化较大，说明计算点周边的复杂地形对改正量造成显著影响。后续将进一步研究改进、提高精度的空间：一是进一步提高输入的地形模型分辨率和精度；二是尽可能扩大积分半径，对较远区地形影响可以采用粗/细格网相结合的方式进行计算提高效率，改善精度；三是采用严密计算方式，提高平面FFT方法计算阶数。

图2 局部重力地形改正值

图3 局部地形改正值结果比较

表1 局部地形改正值结果统计信息 mGal

表2 局部地形改正值比较结果统计信息 mGal

4 结束语

本文给出了局部地形改正的平面积分方法、平面FFT方法、球面积分方法的具体计算表达式，以西部山区为例，计算了重力局部地形改正值，比较了不同结果的差异。结果表明：平面积分方法和平面FFT方法计算结果接近；球面积分计算结果与平面积分方法和平面FFT方法具有较大差异。三种方法存在差异的本质原因是近似面分别作为球面和平面产生的系统性影响不同，在对大面积地形开展积分时，球面更加能贴合地球形状，精度较高；局部地形改正值大小主要与计算点近区地形起伏有关，若近区地形起伏较大会对改正值产生很大影响，若近区地形较为平缓则其改正值多在0mGal附近，重力值几乎不受影响。计算重力局部地形改正时应根据地形起伏、分辨率，计算效率等要求选择合适的计算方法。在平原地区，地形起伏变化不大，无论采用球面积分或是平面积分，近区地形在垂直方向上没有较大差异，积分结果主要受算式中常数因子影响，这种情况考虑平面方法即可。鉴于FFT方法的速度优势，针对大面积的地形改正计算，此方法既能保证精度和球面积分近似，又能提高计算效率；在山区，球面积分对平面积分有很强的改善作用，为精度考虑，采用球面积分方法为宜。