简介
基于地球重力模型(Earth Gravitational Model) 进行一系列解算, 包括大地水准面差距(Geoid Height) 解算以及重力加速度异常解算.最后, 针对 EGM2008, 我会给出这些算法的鲁棒解及其 C 语言实现.
为什么要使用 EGM
在我们需要用到高程信息和重力加速度信息时, 如果精度要求并不是很高, 我们一般通过 WGS84 模型来计算高程和重力加速度. 但在某些高精度场合(例如惯性导航系统)中, 我们需要对 WGS84 模型进行补偿. 这些补偿由 EGM 提供, 包括高程异常和重力异常.EGM 通过卫星对地球各个地点的重力进行测定, 然后给出一系列球谐展开系数. 基于这些系数, 我们可以计算大地水准面差距, 重力加速度异常等.
现在常用的两个 EGM 为 EGM96(1996) 和 EGM2008(2008). 前者的阶数和次数都为 360; 后者的阶数和次数分别为 2190 和 2160.
基本原理
理想的重力势(Gravitational potential, 单位: $m^2s^{-2}$ )如果将地球看作密度均匀的球体, 其重力势 V 可以通过三重积分计算:
$$V=k\int\int\int\frac{\rho}{l}\,dv$$
式中,
k 表示引力常数;
ρ 表示地球密度;
l 表示被吸引点(attracted point)距离;
Geoid 的重力势
在 Geoid 中, 地球的实际重力势通过 n 次和 m 阶的位势(geopotential)系数 $C_{nm}$ 和 $S_{nm}$ 给出.
这些系数以正交级数以及伴随 Legendre 函数(associated Legendre functions, ALF) $\bar{P}_{\text{nm}}$ 表示完全归一化伴随 Legendre 函数(fully normalised associated Legendre functions):
$$V(\Psi,\lambda,r)=\frac{\text{kM}}{r}\left(1+\sum_{n=2}^{n_{\max}}\left(\frac{a}{r}\right)^n\sum_{m=0}^n(\bar{C}_{\text{nm}}\cos(m\lambda)+\bar{S}_{\text{nm}}\sin(m\lambda))\bar{P}_{\text{nm}}(\sin\Psi)\right)$$
式中,
k 表示引力常数;
M 表示地球质量;
r 表示地理半径;
ψ 和 λ 分别表示经度和纬度;
坐标
需要注意的是, 在上面的等式里, 被吸引点以地理坐标(φ, λ, H)表示; 而 Geoid 中的级数展开则是基于球坐标(spherical coordinates), (ψ, λ, r), 所以需要将地理坐标转换为球坐标.
坐标准换是这样实现的: 首先将地理坐标转换为 3-D 笛卡尔坐标, 然后再将其转换为球坐标. 在地理坐标->3D笛卡尔坐标的准换中, 我们需要用到 WGS84 中的一些常数. 这些细节将会在后续章节中讨论.
大地水准面差距和重力加速度异常的计算
基于上面重力势的等式, 我们可以计算大地水准面差距和重力加速度异常.标准高(Normal height, 单位: m)
标准高是在做标准换中实现的(地理坐标->3D 笛卡尔坐标->球坐标), 特性由 WGS84 的常数决定.
标准重力(Normal gravity, 单位: m/s^2)
在纬度 φ 的标准重力 γ(φ) 由下式给出:
$$\gamma(\phi)=\frac{\gamma_e\cos^2\phi+(1-f)\gamma_p\sin^2\phi}{\sqrt{\cos^2\phi+(1-f)^2\sin^2\phi}} $$
其中 γe 和 γp 分别是位于赤道和地球极部的标准重力; f 是 WGS84 中的扁率.
大地水准面差距(Geoid height, 单位: m)
大地水准面差距 ζ 由下面的公式给出:
$$\zeta(\Psi,\lambda,r)=\frac{\text{kM}}{\gamma(\phi)r}\left(\sum_{n=2}^{n_{\max}}\left(\frac{a}{r}\right)^n\sum_{m=0}^n(\bar{C}_{\text{nm}}\cos(m\lambda))+\bar{S}_{\text{nm}}\sin(m\lambda)\bar{P}_{\text{nm}}(\sin\Psi)\right)$$
其中 γ(φ) 是在 φ 的标准重力(normal gravity).
重力加速度异常(Gravitational acceleration anomalies, 单位: m/s^2)
重力加速度异常 Δg 由下面的公式给出:
$$\text{$\Delta$g}(\Psi,\lambda,r)=\frac{\text{kM}}{r^2}\left(\sum_{n=2}^{n_{\max}}(n-1)\left(\frac{a}{r}\right)^n\sum_{m=0}^n\left(\bar{C}_{\text{nm}}\cos(m\lambda)+\bar{S}_{\text{nm}}\sin(m\lambda)\bar{P}_{\text{nm}}(\sin\Psi)\right)\right)$$
看起来很简单
好像只需要下面的步骤:1. 计算球坐标(地理坐标->3D笛卡尔坐标->球坐标);
2. 计算重力势 V(ψ, λ, r);
3. 通过 V, 根据实际需求(大地水准面差距 / 重力加速度)计算对应的值.
理论上是这样的. 可是在 ALF 的经典算法中, 双阶乘的存在使得这些看似简单的步骤在高阶数时变得难以计算. 在高阶数时, 运算将非常难以进行, 并且会产生数值不稳定等问题.
现在我们已经可以对 EGM 进行解算. 但是在实际应用中, 我们需要更加鲁棒的算法. 接下来将介绍这些算法并且用 C 语言实现他们.
关于 ALF, 可以参考 这里.
具体实现
坐标转换
就像上面提到的, 在进行实际的 EGM 解算前, 我们必须先进行坐标转换.WGS84 到 3D 笛卡尔的转换
我们选用的地理坐标系是 WGS84 坐标系, 所以我们需要先把 WGS84 坐标系下的地理坐标转换到 3D 笛卡尔坐标. 这个过程的 C 语言实现如下:
void geoCoor2Cart3d(
const float a, const float finv,
const float phi, const float lambda, const float h,
float *x, float *y, float *z
) {
/* 计算偏心率的平方 */
const float esq = finv < 1.0e-8 ? 0.0 : (2.0 - 1.0 / finv) / finv;
const float sinPhi = sinf(DEG2RAD(phi));
/* N: 卯酉圈曲率半径 */
/* N = a / sqrt(1 - e^2 sin^2 phi) */
const float N = a / sqrtf(1.0 - esq * sinPhi * sinPhi);
const float p = (N + h) * cosf(DEG2RAD(phi));
*x = p * cosf(DEG2RAD(lambda));
*y = p * sinf(DEG2RAD(lambda));
*z = (N * (1.0 - esq) + h) * sinPhi;
}
在 WGS84 中, a = 6378137.0; finv = 298.257223563.
经纬度以角度格式输入, DEG2RAD 宏将角度转为弧度.
3D 笛卡尔坐标到球坐标的转换
这个转换简单得多, 我们只需要用 C 标准库中提供的 hypot 函数和 atan2 函数就可以实现了.
void cart3d2Sph(
const float x, const float y, const float z,
float *az, float *elev, float *r
) {
const float hypotxy = hypotf(x, y);
*r = hypotf(hypotxy, z);
*elev = atan2f(z, hypotxy);
*az = atan2f(y, x);
}
完全归一化伴随 Legendre 函数
上面提到, ALF 的经典算法在高阶数时难以计算, 我们使用 ALF 的递归解, 完全归一化伴随 Legendre 函数的递归解:$$\bar{P}_{n,m}(t)=a_{n,m}\bar{P}_{n-1,m}(t)-b_{n,m}\bar{P}_{n-2,m}(t)$$
其中,
$$a_{n,m}=\sqrt{\frac{(2n-1)(2n+1)}{(n-m)(n+m)}}, b_{n,m}=\sqrt{\frac{(2n+1)(n+m-1)(n-m-1)}{(n-m)(n+m)(2n-3)}}.$$
当 n = m 时, $\bar{P}_{n,m}(t)$ 作为递归式的种子值, 并且 $\bar{P}_{0,0}(t)(t)=1.$ 它的值由下式给出:
$$\bar{P}_{m,m}(t)=u\sqrt{\frac{2m+1}{2m}}\bar{P}_{m-1,m-1}(t).$$
其中,
$$u=\sqrt{1-\sin^2 t}.$$