4. 正常重力场计算与纬度校正

布格异常的本质,是把实测重力值减去“同一点在理想地球模型下应有的重力值”,再扣除高程与中间层影响。其中“理想地球模型下应有的重力值”就是正常重力。它随纬度变化,因此纬度校正并不是一个可选的修饰项,而是布格异常计算链条中的基础环节。本节围绕正常重力公式的选择、不同纬度下的取值、纬度校正的物理含义与实现,以及编程中的误差控制展开。

一、正常重力公式:从椭球模型到1980公式

正常重力场由旋转椭球体(参考椭球)产生,其表面正常重力是纬度的函数。国际上常用的是1980年大地测量参考系统(GRS80)给出的正常重力公式:

γ(φ) = 9.780327 × (1 + 0.0053024·sin²φ − 0.0000058·sin²2φ)   [m/s²]
      = 978032.7 × (1 + 0.0053024·sin²φ − 0.0000058·sin²2φ)   [mGal]

式中 φ 为地理纬度,γ 为椭球面正常重力。第一项反映赤道与两极的重力差(约5.3 Gal),第二项是椭球扁率带来的高阶修正。工程中还有1930年卡西尼斯公式(Cassinis)等旧式公式,二者在数值上可差数毫伽,跨系统对比数据时必须统一公式版本,否则会引入系统性偏差。

需要注意,该公式给出的是椭球面(大地水准面附近)的正常重力。实际测点在海拔 H 处,还需把正常重力归算到测点高度,即“正常重力高度改正”,通常近似为每升高1 m重力减少0.3086 mGal。因此完整的正常场值为:

γ(φ, H) = γ(φ) − 0.3086·H      (H 以 m 计,结果以 mGal 计)

二、不同纬度下的正常重力值

用1980公式计算几个典型纬度(椭球面,单位 mGal):

纬度 φsin²φγ(φ) / mGal
0°(赤道)0978032.7
30°0.2500979321.6
45°0.5000980624.2
60°0.7500981921.9
90°(极点)1983218.6

从赤道到极点,正常重力增大约5.19 Gal,平均每度约57.7 mGal,中纬度处每度约81 mGal。这意味着:若测点纬度存在0.01°(约1.1 km)的定位误差,正常重力就会偏差约0.8 mGal,足以影响中小比例尺重力资料的解释。因此高精度重力测量对点位坐标的要求远高于一般地形测量。

三、纬度校正的物理意义

纬度校正(又称正常场校正)的目的是消除地球椭球形状与自转离心力随纬度变化所引起的重力背景差异,使不同纬度的测点归算到同一比较基准,从而凸显地下密度不均匀产生的局部异常。其物理来源有二:一是地球扁率使两极离地心更近,引力增大;二是离心力随纬度降低而增大(赤道最大、两极为零),使赤道重力减小。两者叠加,形成上述约5.2 Gal的纬向梯度。

纬度校正不是“减去一个常数”,而是逐点减去该点纬度对应的正常重力。校正后的布格异常中,长波纬向趋势被压制,短波密度异常得以保留。

四、实现步骤与编程示例

标准实现流程为:读取测点经纬度与高程 → 按统一椭球公式计算 γ(φ) → 施加高度改正得到 γ(φ,H) → 用实测重力减去该值,得到纬度校正后的重力(再进入布格改正)。下面给出可直接运行的 Python 片段:

import numpy as np

def normal_gravity(lat_deg, h_m=0.0):
    """GRS80 正常重力,返回 mGal;lat_deg 为地理纬度,h_m 为海拔"""
    phi = np.radians(lat_deg)
    s2 = np.sin(phi)**2
    g = 978032.7 * (1 + 0.0053024*s2 - 0.0000058*np.sin(2*phi)**2)
    return g - 0.3086 * h_m

# 示例:北纬30°、海拔500 m 的测点
print(normal_gravity(30.0, 500.0))   # 约 979167.3 mGal

实际工程中应把公式封装为批量函数,对整条测线一次性计算,并保留公式版本号与椭球参数,便于资料归档与复算。

五、误差分析与常见陷阱

纬度校正的误差主要来自三方面:定位误差(纬度精度直接放大为重力误差,0.001°约0.08 mGal)、公式版本不一致(新旧公式差可达数毫伽)、高度改正系数近似(0.3086 mGal/m 是平均密度假设下的近似,高海拔或高精度测量应改用更严格的正常重力高度公式)。此外,若误将大地纬度与地心纬度混用,在高纬度地区会带来明显偏差。建议在程序中显式注明纬度类型,并用已知标准点(如某基准站)做一次闭环校验,确认计算结果与权威值一致后再批量处理。

正常重力随纬度变化与纬度校正流程 γ(φ) 纬向变化(椭球面) 纬度 φ / ° γ / mGal 0 30 45 60 90 978000 980600 983200 978033 980624 983219 赤道→极点增大约 5.19 Gal 纬度校正处理链 测点纬度 φ、高程 H γ(φ)=978032.7(1+0.0053024sin²φ−…) GRS80 椭球面正常重力 γ(φ,H)=γ(φ)−0.3086·H 高度改正(mGal) g_纬校 = g_实测 − γ(φ,H) 进入布格改正 误差源:定位 0.001°≈0.08 mGal;公式版本差数 mGal

把上述环节串起来看,纬度校正的核心是“逐点、按统一椭球、含高度改正”地计算正常重力并相减。只要公式版本、纬度类型与高度改正系数三者保持一致,纬度校正就能稳定地把纬向背景从布格异常中剥离,为后续密度反演提供干净的输入。