最近在迁移一个处理 GPS 数据的 pipeline 到 BigQuery 上。期间在 Google Cloud 的犄角旮旯里遇到了各种各样稀奇古怪的错误，让人防不胜防。我都怀疑我是 Google Cloud 找来帮他们修 bug 的免费劳力。其中最让人心烦的是 BigQuery 的 geometry 类型仅支持 WGS 84（World Geodetic System 1984）坐标系（当然这也是国际上比较通用的坐标系），而且它的计算默认是考虑球面坐标的，这就导致了一系列坐标误差带来的结果不同。

但一开始我并没有意识到这一点。我手头上有几个不同来源的数据：

1. 需要处理的汽车行驶 GPS 记录 epsg:4326 (WGS 84)
2. 从地图提供商获取的地图道路数据 epsg:4612 (JGD 2000)
3. 从日本国土地理院获取的公开数据 epsg:6668 (JGD 2011)

理论上这几个坐标系是非常近似的。而这也正是 JGD 2000 标准设定的初衷：基于全球统一的世界坐标系 国际地球参考框架 ITRF 94 和 GRS 80 设定的椭球体。

那么 ITRF 94 和 GRS 80 究竟是啥玩意儿呢？

^* 图1. 一个理想的标准椭球体

GRS 叫 大地参考系统（Geodetic Reference System）是一套理论上的地球数学模型，它试图用一套参数来描述一个理想的、均质的一个参考椭球体。国际大地测量学和地球物理学联合会（IUGG）在 1979 年采纳的这个计算方法叫 GRS 80。

我们都知道对于一个椭球体来说，设长轴半径为 $a$，短轴半径为 $b$，那么这个椭球在直角坐标系 $O_{xyz}$ 中可表示为：

$$ \frac{x^2}{a^2} + \frac{y^2}{a^2} + \frac{z^2}{b^2} = 1 $$

但地球在物理性质上肯定不是一个绝对刚体，随着它自己的自转公转地球都会发生一定程度的形变。加之占据了地球大部分表面的海洋，会因为地月地日的潮汐效应不断变动，整个地球就是个扭来扭去的凹凸不平的大土豆。

图2. 一个夸张 10000 倍的不平整的地球

虽说地球扭来扭去，潮汐升升落落，但它总不是在乱扭。海平面升降转换之间总该有一个中间位置，海浪在这个中间位置上下浮动。假设一下，如果地球在一个完全静止的空间中，它自己也不动的话，海洋自然也是平静的。那这个时候的平静海面不就是地球的表面了吗？

图3. 等势面与大地水准面

从这个思路出发，我们认识到，其实在一个绝对平静、不受任何干扰的全球海洋所形成的表面，其实不就是一个重力势能完全相等的一个曲面吗？那么我们只要计算出一个重力的等位势面出来，不就可以得到地球表面了吗？这个理想中的表面也叫“大地水准面”。

^* 图4. 大地水准面也是不规则的

想象你在爬山。通常说的海拔高度，字面意思指的是距离海平面的高度。在这个高度你积累的重力势能，按理说，和你从海平面的位置爬上相同高度的位置的重力势能应该是一样的。把一块石头从等高线上的一点移动到另一点，只要高度不变，重力对它不做功。所以我们说的地图上的等高线，在物理学上就应该是“等势线”。

不过问题还没结束，地球的重力场也不是完美的对称球体。因为地球在自转，它会产生离心力影响势能。而且地球是扁的，质量其实分布不均，重力所指向的“下”并不一定完全指向地球球心。比如说在山上，山的质量产生的引力会使得在该处的重力偏向山。相同的“铅垂线”一旦拉长实际上成了一个弧形。这导致地球周围的等位势面变得非常复杂。最麻烦的是它也不是一个椭球形，我们计算的时候没有办法简便地表示它。

那怎么办呢，只能人为规定一个比较合理而且好算的一个等势面当作真实地球的重力基准面了。观察地球的重力场，我们可以理想化地认为在地球外部的真空空间中没有质量分布，这样引力势 $U$ 就成了一个没有源的势场，满足拉普拉斯方程：$\nabla^2 U = 0$，其中 $\nabla^2$ 是拉普拉斯算子。它在笛卡尔坐标系（$x,y,z$）下的标准定义为：

$$ \nabla^2 U = \frac{\partial^2 U}{\partial x^2} + \frac{\partial^2 U}{\partial y^2} + \frac{\partial^2 U}{\partial z^2} = 0 $$

^* 图6. 在笛卡尔坐标系上构建球坐标系

具体到一个球坐标系 $(r, \theta, \phi)$ 上到时候，我们有

• $r$：地心距离。
• $\theta$：余纬 (Colatitude)，即从北极点 ($z$轴) 向下的角度 ($0 \le \theta \le \pi$)。
• $\phi$：经度 (Longitude) ($0 \le \phi \le 2\pi$)。

笛卡尔坐标 ($x, y, z$) 与球坐标 ($r, \theta, \phi$) 的变换关系为： $$ \begin{cases} x = r \sin\theta \cos\phi \ y = r \sin\theta \sin\phi \ z = r \cos\theta \end{cases} $$

而且

$$ r^2 = x^2 + y^2 + z^2 $$

接下来我们需要把拉普拉斯算子从笛卡尔坐标系转换到球坐标系上。

拉普拉斯算子在球坐标系下的推导

为了体现这个过程的一般性，也是顺带复习一下大一学的高数，这里重新推导一下。

首先我们定义一个一般的标量场 $\psi$。在空间中移动微小距离 $d\mathbf{l}$，对应的标量场 $\psi$ 的变化量 $d\psi$ 为：

$$d\psi = \nabla \psi \cdot d\mathbf{l}$$

在球坐标系中，微小位移矢量 $d\mathbf{l}$ 由三段弧长组成：

$$d\mathbf{l} = (dr)\hat{r} + (r d\theta)\hat{\theta} + (r\sin\theta d\phi)\hat{phi}$$

• 径向移动距离是 $dr$
• 角向方向移动距离是弧长 $r d\theta$
• 方位向方向移动距离是弧长 $r\sin\theta d\phi$

为了使 $\nabla \psi \cdot d\mathbf{l}$ 等于全微分 $d\psi = \frac{\partial \psi}{\partial r}dr + \frac{\partial \psi}{\partial \theta}d\theta + \frac{\partial \psi}{\partial \phi}d\phi$，$\nabla \psi$ 在各个基矢量方向的分量必须抵消掉弧长系数。

因此，球坐标梯度算子为：

$$ \nabla = \hat{r}\frac{\partial}{\partial r} + \hat{\theta}\frac{1}{r}\frac{\partial}{\partial \theta} + \hat{\phi}\frac{1}{r\sin\theta}\frac{\partial}{\partial \phi} $$

为了进一步计算拉普拉斯算子 $\nabla^2 \psi = \nabla \cdot (\nabla \psi)$，这里引入基矢量的微分关系。

在笛卡尔坐标系中，$\hat{i}, \hat{j}, \hat{k}$ 是常矢量，导数为0。但在球坐标系中，$\hat{r}, \hat{\theta}, \hat{\phi}$ 的方向随位置改变，它们对坐标的偏导数不为零。它们与笛卡尔基矢量的关系为：

$$ \begin{aligned} \hat{r} &= \sin\theta\cos\phi \hat{i} + \sin\theta\sin\phi \hat{j} + \cos\theta \hat{k} \ \hat{\theta} &= \cos\theta\cos\phi \hat{i} + \cos\theta\sin\phi \hat{j} - \sin\theta \hat{k} \ \hat{\phi} &= -\sin\phi \hat{i} + \cos\phi \hat{j} \end{aligned} $$

接下来对球坐标系的基矢量求偏导：

对 $\partial r$ 求导：

$$\frac{\partial \hat{r}}{\partial r} = 0, \quad \frac{\partial \hat{\theta}}{\partial r} = 0, \quad \frac{\partial \hat{\phi}}{\partial r} = 0$$

对 $\partial \theta$ 求导：

$$ \frac{\partial \hat{r}}{\partial \theta} = \hat{\theta}, \quad \frac{\partial \hat{\theta}}{\partial \theta} = -\hat{r}, \quad \frac{\partial \hat{\phi}}{\partial \theta} = 0 $$

对 $\partial \phi$ 求导：

$$ \frac{\partial \hat{r}}{\partial \phi} = \sin\theta \hat{\phi}, \quad \frac{\partial \hat{\theta}}{\partial \phi} = \cos\theta \hat{\phi}, \quad \frac{\partial \hat{\phi}}{\partial \phi} = -(\sin\theta \hat{r} + \cos\theta \hat{\theta}) $$

现在我们把梯度算子作用在梯度向量上。设 $\nabla \psi = \mathbf{A} = A_r \hat{r} + A_\theta \hat{\theta} + A_\phi \hat{\phi}$，其中：

$$ A_r = \frac{\partial \psi}{\partial r}, \quad A_\theta = \frac{1}{r}\frac{\partial \psi}{\partial \theta}, \quad A_\phi = \frac{1}{r\sin\theta}\frac{\partial \psi}{\partial \phi} $$

拉普拉斯算子即为在前面求得的球坐标梯度 $\nabla = \hat{r}\frac{\partial}{\partial r} + \hat{\theta}\frac{1}{r}\frac{\partial}{\partial \theta} + \hat{\phi}\frac{1}{r\sin\theta}\frac{\partial}{\partial \phi}$ 的基础上再求散度：

$$ \nabla^2 \psi = \nabla \cdot \mathbf{A} = \left( \hat{r}\frac{\partial}{\partial r} + \frac{\hat{\theta}}{r}\frac{\partial}{\partial \theta} + \frac{\hat{\phi}}{r\sin\theta}\frac{\partial}{\partial \phi} \right) \cdot (A_r \hat{r} + A_\theta \hat{\theta} + A_\phi \hat{\phi}) $$

这里分成三部分来求和：

1. 径向分量 $A_r \hat{r}$ 的散度 $$ \nabla \cdot (A_r \hat{r}) = (\nabla A_r) \cdot \hat{r} + A_r (\nabla \cdot \hat{r}) $$ 第一项，因为只有 $\hat{r}$ 方向的导数非零，可以直接得到 $$ (\nabla A_r) \cdot \hat{r} = \frac{\partial A_r}{\partial r} $$ 第二项中 $\nabla \cdot \hat{r}$ 可以展开为 $$ \nabla \cdot \hat{r} = \frac{1}{r}\hat{\theta} \cdot \frac{\partial \hat{r}}{\partial \theta} + \frac{1}{r\sin\theta}\hat{\phi} \cdot \frac{\partial \hat{r}}{\partial \phi} $$ 代入微分关系可以进一步得到 $$ \nabla \cdot \hat{r} = \frac{1}{r}\hat{\theta}\cdot\hat{\theta} + \frac{1}{r\sin\theta}\hat{\phi}\cdot(\sin\theta\hat{\phi}) = \frac{1}{r} + \frac{1}{r} = \frac{2}{r} $$ 加起来得到径向分量 $A_r \hat{r}$ 的散度 $$ \nabla \cdot (A_r \hat{r}) = \frac{\partial A_r}{\partial r} + \frac{2}{r}A_r $$

2. 角向分量 $A_\theta \hat{\theta}$ 的散度

   $$ \nabla \cdot (A_\theta \hat{\theta}) = (\nabla A_\theta) \cdot \hat{\theta}

   • A_\theta (\nabla \cdot \hat{\theta}) $$

   第一项只有 $\hat{\theta}$ 方向导数贡献，可以直接得到 $$(\nabla A_\theta) \cdot \hat{\theta} = \frac{1}{r}\frac{\partial A_\theta}{\partial \theta}$$ 第二项中 $\nabla \cdot \hat{\theta}$ 可以展开为 $$ \nabla \cdot \hat{\theta} = \frac{1}{r}\hat{\theta} \cdot \frac{\partial \hat{\theta}}{\partial \theta}

   • \frac{1}{r\sin\theta}\hat{\phi} \cdot \frac{\partial \hat{\theta}}{\partial \phi} $$ 代入微分关系可以进一步得到 $$ \nabla \cdot \hat{\theta} = \frac{1}{r}\hat{\theta}\cdot(-\hat{r}) + \frac{1}{r\sin\theta}\hat{\phi}\cdot(\cos\theta\hat{\phi}) = 0 + \frac{\cos\theta}{r\sin\theta} = \frac{\cot\theta}{r} $$ 加起来得到角向分量 $A_\theta \hat{\theta}$ 的散度 $$ A_\theta \hat{\theta} = \frac{1}{r}\frac{\partial A_\theta}{\partial \theta} + \frac{\cot\theta}{r}A_\theta $$

3. 方位向分量 $A_\phi \hat{\phi}$ 的散度 $$ \nabla \cdot (A_\phi \hat{\phi}) = (\nabla A_\phi) \cdot \hat{\phi}

   • A_\phi (\nabla \cdot \hat{\phi}) $$ 第一项只有 $\hat{\phi}$ 方向导数贡献，可以直接得到 $$ \frac{1}{r\sin\theta}\frac{\partial A_\phi}{\partial \phi} $$ 第二项：$\nabla \cdot \hat{\phi}$ 可以展开为 $$ \nabla \cdot \hat{\phi} = \frac{1}{r}\hat{\theta} \cdot \frac{\partial \hat{\phi}}{\partial \theta}
   • \frac{1}{r\sin\theta}\hat{\phi} \cdot \frac{\partial \hat{\phi}}{\partial \phi} $$ 利用 $\frac{\partial \hat{\phi}}{\partial \phi}$ 的结果，点积均为0。$\hat{\phi}$ 垂直于 $\hat{r}$ 和 $\hat{\theta}$ 的变化方向，且 $\partial_\phi \hat{\phi}$ 垂直于 $\hat{\phi}$。所以散度为 0。 最后可以得到方位向分量 $A_\phi \hat{\phi}$ 的散度 $$A_\phi \hat{\phi} = \frac{1}{r\sin\theta}\frac{\partial A_\phi}{\partial \phi}$$

把这三个方向的分量加起来就是：

$$ \nabla^2 \psi = \nabla \cdot \mathbf{A} = (\frac{\partial A_r}{\partial r} + \frac{2}{r}A_r)

• (\frac{1}{r}\frac{\partial A_\theta}{\partial \theta} + \frac{\cot\theta}{r}A_\theta)
• (\frac{1}{r\sin\theta}\frac{\partial A_\phi}{\partial \phi}) $$

将 $A_r, A_\theta, A_\phi$ 的定义代回：

$$ \left[\frac{\partial}{\partial r}\left(\frac{\partial \psi}{\partial r}\right) + \frac{2}{r}\left(\frac{\partial \psi}{\partial r}\right)\right] + \left[\frac{1}{r}\frac{\partial}{\partial \theta}\left( \frac{1}{r}\frac{\partial \psi}{\partial \theta} \right) + \frac{\cot\theta}{r}\left( \frac{1}{r}\frac{\partial \psi}{\partial \theta} \right)\right] + \left[\frac{1}{r\sin\theta}\frac{\partial}{\partial \phi}\left( \frac{1}{r\sin\theta}\frac{\partial \psi}{\partial \phi} \right)\right] $$

进一步化简，我们就得到了球坐标系下的拉普拉斯算子：

$$\nabla^2 \psi = \frac{1}{r^2}\frac{\partial}{\partial r}\left( r^2 \frac{\partial \psi}{\partial r} \right) + \frac{1}{r^2\sin\theta}\frac{\partial}{\partial \theta}\left( \sin\theta \frac{\partial \psi}{\partial \theta} \right) + \frac{1}{r^2\sin^2\theta}\frac{\partial^2 \psi}{\partial \phi^2}$$

勒让德多项式的引入

回到地球的重力场，前面说了我们可以理想化地认为地球是一个没有源的势场 $U$，满足拉普拉斯方程：$\nabla^2 U = 0$，所以我们有：

$$ \nabla^2 U = \frac{1}{r^2}\frac{\partial}{\partial r}\left( r^2 \frac{\partial U}{\partial r} \right) + \frac{1}{r^2\sin\theta}\frac{\partial}{\partial \theta}\left( \sin\theta \frac{\partial U}{\partial \theta} \right) + \frac{1}{r^2\sin^2\theta}\frac{\partial^2 U}{\partial \phi^2} = 0 $$

前面的拉普拉斯算子经过在球坐标系下的推导，在一般的标量场 $\psi$ 下可以独立拆分成三部分：径向 ($r$)，角向 ($\theta$)，还有方位向 ($\phi$)。那么我们也可以认为在地球势场 $U$ 下也可以有类似的情况，所以这里设 $U(r, \theta, \phi)$ 也有变量分离的形式解：

$$ U(r, \theta, \phi) = R(r) \cdot \Theta(\theta) \cdot \Phi(\phi) $$

然后我们利用我们设的重力场的分离解 $U = R \cdot \Theta \cdot \Phi$ 代入拉普拉斯方程：

$$ \frac{\Theta \Phi}{r^2} \frac{d}{dr} \left(r^2 \frac{dR}{dr} \right) + \frac{\Phi R}{r^2 \sin\theta}\frac{d}{d\theta} \left(\sin\theta\frac{d\Theta}{d\theta} \right) + \frac{\Theta R}{r^2 \sin^2\theta}\frac{d^2\Phi}{d\phi^2} = 0 $$

两边同时乘以 $\frac{r^2}{\Theta\Phi R}$，并移项，得到：

$$ \frac{1}{R}\frac{d}{dr}\left(r^2\frac{dR}{dr}\right) = -\frac{1}{\Theta\sin\theta}\frac{d}{d\theta}\left(\sin\theta\frac{d\Theta}{d\theta}\right) -\frac{1}{\Phi \sin^2\theta}\frac{d^2\Phi}{d\phi^2} $$

现在左边只与 $r$ 有关，右边只与角度有关。为了后续求解方便，我们将这个常数设为 $l(l+1)$，其中 $l$ 是分离常数，一般为实数。于是我们就有：

$$ \frac{1}{R}\frac{d}{dr}\left(r^2\frac{dR}{dr}\right) = l(l+1) $$

$$ \frac{1}{\Theta\sin\theta}\frac{d}{d\theta}\left(\sin\theta\frac{d\Theta}{d\theta}\right) +\frac{1}{\Phi \sin^2\theta}\frac{d^2\Phi}{d\phi^2} =-l(l+1) $$

这里的第一个式子，我们就分离出了 $r$ 的方程：

$$ \frac{1}{R}\frac{d}{dr}\left(r^2\frac{dR}{dr}\right) = l(l+1) $$

它的通解为 $R(r)=A_lr^l + B_l\frac{1}{r^{l+1}}$。

给第二个式子两边同时乘以 $\sin^2\theta$，并移项，得到：

$$ \frac{1}{\Theta}\sin\theta\frac{d}{d\theta}\left(\sin\theta\frac{d\Theta}{d\theta}\right) +l(l+1)\sin^2\theta =-\frac{1}{\Phi}\frac{d^2\Phi}{d\phi^2} $$

现在左边只与 $\theta$ 有关，右边只与 $\phi$ 有关。为了后续求解方便，我们再次引入常数 $m^2$。于是我们就有：

$$ -\frac{1}{\Phi}\frac{d^2\Phi}{d\phi^2} = m^2 $$

$$ \frac{1}{\Theta}\sin\theta\frac{d}{d\theta}\left(\sin\theta\frac{d\Theta}{d\theta}\right) +l(l+1)\sin^2\theta =m^2 $$

到这一步我们已经分离出了 $\phi$ 的方程：

$$ \frac{d^2 \Phi}{d\phi^2} = -m^2 \Phi $$

这是一个简谐振动方程，解是 $\Phi(\phi) = B_1\cos(m\phi)+ B_2\sin(m\phi)$。这里的 $m$ 必须是整数（因为转一圈 $2\pi$ 就会回到原点）。这就和绕地球纬度方向的运动一致了：沿着纬度方向最后就会绕地一周。

针对第二个式子，我们继续整理后如下：

$$ \frac{1}{\sin\theta} \frac{d}{d\theta} \left( \sin\theta \frac{d\Theta}{d\theta} \right) + \left[ l(l+1) - \frac{m^2}{\sin^2\theta} \right] \Theta = 0 $$

现在里面充满了三角函数 $\sin\theta$。为了把它变成多项式方程，我们必须消灭三角函数。我们引入一个新的自变量 $x$，定义为：

$$x = \cos\theta$$

它的取值范围是：当 $\theta$ 从 $0$ 到 $\pi$ 时，$x$ 从 $1$ 到 $-1$。

现在我们需要把对 $\theta$ 的导数变成对 $x$ 的导数：

$$ \frac{d}{d\theta} = \frac{dx}{d\theta} \frac{d}{dx} = (-\sin\theta) \frac{d}{dx} = -\sqrt{1-x^2} \frac{d}{dx} $$

然后可以求得

$$ \sin\theta \frac{d}{d\theta} = \sin\theta (-\sin\theta \frac{d}{dx}) = -\sin^2\theta \frac{d}{dx} = -(1-x^2) \frac{d}{dx} $$

继续代换到上面的 $\theta$ 方程中，就可以得到一个不含任何三角函数的形式：

$$ \frac{d}{dx} \left[ (1-x^2) \frac{d\Theta}{dx} \right] + \left[ l(l+1) - \frac{m^2}{1-x^2} \right] \Theta = 0 $$

或者展开写成标准形式：

$$ (1-x^2) \frac{d^2\Theta}{dx^2} - 2x \frac{d\Theta}{dx} + \left[ l(l+1) - \frac{m^2}{1-x^2} \right] \Theta = 0 $$

这得到了 “连带勒让德微分方程 (Associated Legendre Differential Equation)”。

连带勒让德微分方程的特点就是，只有当 $l$ 是整数 ($0, 1, 2…$) 时，这个方程在区间 $x \in [-1, 1]$ 上才有有限的解（这很重要，因为物理上的势场不能无穷大）。这种形式被称为连带勒让德多项式 (Associated Legendre Polynomials)，这个对 $\theta$ 的解就是 $\Theta(\theta) = \sum_{l=0}^\infty \sum_{m=0}^M P_l^m(x) = \sum_{l=0}^\infty P_l^m(\cos\theta)$，其中 $P_l^m(x)$ 展开为：

$$ P_l^m(x) = \frac{1}{2^l} (-1)^m \frac{(2l-2m)!}{m!(l-m)!(l-2m)!}x^{l-2m} $$

以上我们就得到了三个方向上的解。

$$R(r)=A_lr^l + B_l\frac{1}{r^{l+1}}$$ $$\Theta(\theta) = \sum_{l=0}^\infty \sum_{m=0}^M P_l^m(x)$$ $$\Phi(\phi) = B_1\cos(m\phi)+ B_2\sin(m\phi)$$

但是对于径向上的解，我们的求解区域是地球表面以外的空间：$r > R_{\text{earth}}$。

首先看 $A_l r^l$，当 $r \to \infty$（飞向宇宙深处）时，对于 $l = 0$：$A_0$ 是常数。对于 $l \ge 1$：$A_l r^l$ 会趋向于无穷大。也就是说如果 $l \ge 1$，意味着你离地球越远，地球的引力势就越强。这显然违反了物理常识。因此，为了满足物理现实，所以这一部分的 $A_l = 0$。

对于另外一部分 $B_l \frac{1}{r^{l+1}}$，当 $r \to \infty$ 时，$\frac{1}{r^{l+1}}$ 趋向于 0。这完全符合我们的预期：离地球越远，引力影响越微弱，最终消失。所以，我们保留这一项。

这样一来，径向解 $R(r)$ 就变成：

$$R(r) = B_l\frac{1}{r^{l+1}} = \frac{B_{lm}}{r^{l+1}}$$

因此，拉普拉斯方程 $\nabla^2 U = 0$ 在球坐标系下的通解就是径向解和角向解以及方位角向解的乘积，并对所有可能的 $l$ 和 $m$ 进行叠加：

$$ U(r, \theta, \phi) = R(r) \cdot \Theta(\theta) \cdot \Phi(\phi) = \sum_{l=0}^{\infty} \sum_{m=-l}^{l} \frac{B_{lm}}{r^{l+1}} P_l^m(x) \left( B_1\cos(m\phi)+ B_2\sin(m\phi) \right) $$

不过现在还有一个量纲混乱的问题， 当 $l=0$ 时，项是 $1/r$，所以 $B_0$ 的单位是 $m^3/s^2$。 当 $l=2$ 时，项是 $1/r^3$，所以 $B_2$ 的单位是 $m^5/s^2$。

我们强制提取出公因子 $\frac{GM}{r}$，并引入地球半径 $a$ 来消除量纲。公式变成了这样：

$$ U(r, \theta, \phi) = \frac{GM}{r} \sum_{l=0}^{\infty} \sum_{m=0}^{l} \left( \frac{a}{r} \right)^l \bar{P}{lm}(\cos\theta) [ \bar{C}{lm} \cos(m\phi) + \bar{S}_{lm} \sin(m\phi) ] $$

GRS 80 的定义

GRS80 定义了一个“理想的旋转椭球体”。我们取一个二维的椭圆（子午圈）。以短轴（Z轴）为中心，将其旋转 360 度。这就形成了一个旋转椭球体。它绕 Z 轴（自转轴）完全对称。

这样一来，只要你的纬度 ($\theta$) 和高度 ($r$) 是一样的，你脚下的地球形状就是完全一样的，地下的质量分布也是完全一样的。既然质量分布与经度无关，那么产生的引力势自然也与经度无关。

我们在前面推导过，通用的球谐函数项是：

$$ P_l^m(\cos\theta) \cdot [C_{lm} \cos(m\phi) + S_{lm} \sin(m\phi)] $$

这里的 $m$ (Order, 次数) 专门负责描述随经度 $\phi$ 的变化频率。在 GRS80 中，由于定义了旋转对称性 (Rotational Symmetry)，意味着势场 $U$ 对 $\phi$ 的导数必须为 0：

$$ \frac{\partial U}{\partial \phi} = 0 $$

为了满足这个条件，数学上唯一的选择就是强制所有 $m \neq 0$ 的系数全部为零。只保留 $m=0$ 的项。

$$\cos(0 \cdot \phi) + \sin(0 \cdot \phi) = 1 + 0 = 1$$

到最后 GRS80 的重力势公式就简化成了：

$$ U(r, \theta) = \frac{GM}{r} \left[ 1 - \sum_{n=1}^{\infty} J_{2n} \left( \frac{a}{r} \right)^{2n} P_{2n}(\cos\theta) \right] + \frac{1}{2}\omega^2 r^2 \sin^2\theta $$

这里用到的 $P_{2n}(\cos\theta)$ 其实就是 $P_{2n}^0(\cos\theta)$。