A Slight Diffs on Coordinate

最近在迁移一个处理 GPS 数据的 pipeline 到 BigQuery 上。期间在 Google Cloud 的犄角旮旯里遇到了各种各样稀奇古怪的错误，让人防不胜防。我都怀疑我是 Google Cloud 找来帮他们修 bug 的免费劳力。其中最让人心烦的是 BigQuery 的 GEOGRAPHY 类型仅支持 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$$

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

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

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

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

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

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

$$\nabla^2 V = \frac{\partial^2 V}{\partial x^2} + \frac{\partial^2 V}{\partial y^2} + \frac{\partial^2 V}{\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}$ 的散度

   $$\nabla \cdot (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}$$

   两项都是 0：前一项因为 $\frac{\partial \hat{\phi}}{\partial \theta} = 0$；后一项因为 $\hat{\phi} \cdot \frac{\partial \hat{\phi}}{\partial \phi} = \hat{\phi} \cdot \left( -(\sin\theta \hat{r} + \cos\theta \hat{\theta}) \right) = 0$，$\partial_\phi \hat{\phi}$ 落在 $\hat{r}$-$\hat{\theta}$ 平面里，和 $\hat{\phi}$ 正交。所以 $\nabla \cdot \hat{\phi} = 0$。 最后可以得到方位向分量 $A_\phi \hat{\phi}$ 的散度

   $$\nabla \cdot (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}$$

勒让德多项式的引入

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

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

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

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

然后我们利用我们设的引力位的分离解 $V = 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$ 就会回到原点）。这也和地球上的实际情形对得上：$\phi$ 就是经度，沿着 $\phi$ 增大的方向（也就是沿着一条纬线）一直走，最后一定绕回原点。

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

$$\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, \dots$) 并且 $|m| \le l$ 时，这个方程在闭区间 $x \in [-1, 1]$ 上（也就是包含了南北两极）才有有限的解（因为物理上的势场不能在极点发散）。这个解叫连带勒让德函数 (Associated Legendre Function)，记作 $P_l^m(x)$。

也就是说，对给定的一对 $(l, m)$，$\theta$ 方向的解就是

$$\Theta(\theta) = P_l^m(\cos\theta)$$

它本身可以从普通的勒让德多项式 $P_l$ 求 $m$ 阶导数得到，后一个式子叫 Rodrigues 公式：

+ 通常物理教材里的定义通常还会在前面乘一个 $(-1)^m$（Condon–Shortley 相位），大地测量学的惯例是不带，所以配套的球谐系数符号也一样有区别。 $$P_l(x) = \frac{1}{2^l\,l!}\frac{d^l}{dx^l}\left(x^2 - 1\right)^l, \qquad P_l^m(x) = \left(1 - x^2\right)^{m/2}\frac{d^m}{dx^m}P_l(x)$$

$m = 0$ 的时候连带勒让德函数退化成普通的勒让德多项式，这时可以直接写成一个有限次多项式：

$$P_l(x) = \frac{1}{2^l} \sum_{k=0}^{\lfloor l/2 \rfloor} (-1)^k \frac{(2l-2k)!}{k!\,(l-k)!\,(l-2k)!}\, x^{l-2k}$$

这里求和的 $k$ 只是展开的项号，和描述经度频率的那个 $m$ 毫无关系。这两个下标长得实在太像，是整段推导里最容易看串行的地方。

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

$$R(r)=A_lr^l + B_l\frac{1}{r^{l+1}}$$ $$\Theta(\theta) = P_l^m(\cos\theta)$$ $$\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 V = 0$ 在球坐标系下的通解就是径向解和角向解以及方位角向解的乘积，并对所有可能的 $l$ 和 $m$ 进行叠加：

$$V(r, \theta, \phi) = \sum_{l=0}^{\infty} \sum_{m=0}^{l} \frac{P_l^m(\cos\theta)}{r^{l+1}} \left( C_{lm}\cos(m\phi)+ S_{lm}\sin(m\phi) \right)$$

这里 $m$ 只从 $0$ 数到 $l$ 就够了：$\cos(m\phi)$ 和 $\sin(m\phi)$ 两项加在一起已经把 $\pm m$ 都装进去了，再从 $-l$ 数一遍就是重复计数。两个系数 $C_{lm}$、$S_{lm}$ 则把前面 $R(r)$ 和 $\Phi(\phi)$ 里那些待定常数一并吸收了进去。

不过现在还有一个量纲混乱的问题， 当 $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$ 来消除量纲。公式变成了这样：

$$V(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)，意味着引力位 $V$ 对 $\phi$ 的导数必须为 0：

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

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

$$\cos(0 \cdot \phi) = 1, \qquad \sin(0 \cdot \phi) = 0$$

经度方向被彻底抹掉了。剩下的 $m = 0$ 项叫带谐项（Zonal Harmonics），因为它们只随纬度变化，等值线是一条条平行于赤道的「带」。于是双重求和塌缩成单重的：

+ $P_l(\cos\theta)$ 就是 $P_l^0(\cos\theta)$，即 $m=0$ 时退化成的普通勒让德多项式。 $$V(r, \theta) = \frac{GM}{r} \sum_{l=0}^{\infty} \left( \frac{a}{r} \right)^{l} C_{l0}\, P_l(\cos\theta)$$

接着还能靠三个条件继续删项。

$l = 0$ 项就是牛顿的点质量解。 $P_0 = 1$，取 $C_{00} = 1$，这一项给出 $GM/r$，相当于把地球当成一个质点时的引力位。

$l = 1$ 的三项恒为零。 球谐展开里的一次项正比于地球质心相对于坐标原点的偏移量。既然我们把坐标原点取在地球质心上（现代基准的定义如此），这个偏移量就是零，一次项自动消失。

所有奇数次项也为零。 旋转椭球关于赤道面是南北对称的 真实的地球并不完全南北对称（南半球稍微「胖」一点，俗称梨形），所以实测的 $J_3$ 不为零，只是量级比 $J_2$ 小了三个数量级。不过毕竟 GRS 80 是个理想模型，它直接定义了严格南北对称的。 ，也就是 $V(r, \theta) = V(r, \pi - \theta)$。而勒让德多项式满足 $P_l(-x) = (-1)^l P_l(x)$，奇数 $l$ 的 $P_l$ 是奇函数，代入 $x \to -x$ 会变号。要同时满足南北对称，奇数次项的系数只能全取零。

删到最后只剩下偶数次带谐项。大地测量学习惯把系数换成 $J_l = -C_{l0}$ 来写，负号是这个定义带来的：

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

其中 $J_2$ 被称为动力形状因子，量级大概是 $10^{-3}$。这个参数对建模出来的地球的“扁形”贡献最大。

到这里还没完。上面整套球谐展开解的是引力位 $V$ 的拉普拉斯方程。前面说了因为地球一直在自转，我们脚下所谓的“重力”里还有离心力的那一份贡献。

地球上一点到自转轴的距离是 $d = r\sin\theta$，根据高中物理可知离心力为：

$$Z = \frac{1}{2}\omega^2 d^2 = \frac{1}{2}\omega^2 r^2 \sin^2\theta$$

两者相加，才得到 GRS 80 的正常重力位：

+ 顺带一提这个重力位表达式里 $\nabla^2 V = 0$，而 $\nabla^2 U = 2\omega^2$。 $$U(r, \theta) = V + Z = \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$$

到这一步终于有了等势面的完整的表达式。然后我们就可以用 $a$、$GM$、$J_2$、$\omega$ 这四个参数定义一个椭球。

+ 查公开资料可知到 GRS 80 的四个定义常数：$a = 6378137\ \text{m}$（长半轴）、$GM = 3.986005 \times 10^{14}\ \text{m}^3/\text{s}^2$（地心引力常数）、$J_2 = 1.08263 \times 10^{-3}$（动力形状因子）、$\omega = 7.292115 \times 10^{-5}\ \text{rad/s}$（自转角速度）。至于扁率 $1/f = 298.257222101$ 是从前面四个数推算出来的。

最后在这个正常重力场里挑出一个特定的等势面 $U = U_0$，让它尽可能贴近真实的平均海面，这个面就是 GRS 80 的参考椭球面。

从重力势回到高度

+ 这一套记号是大地测量学的通行惯例。$V$ 是引力位、$Z$ 是离心力位、$U = V + Z$ 是正常重力位（理想椭球模型的）、$W$ 是真实地球的重力位、$T = W - U$ 是扰动位。

到这里依然还没结束。真实地球的重力位 $W$ 和上面建模出来的正常重力位 $U$ 并不一致。它们之间的差叫扰动位 $T = W - U$。真实的大地水准面和参考椭球面之间的高度差，可以由 Bruns 公式给出：

$$N = \frac{T}{\gamma}$$

其中 $\gamma$ 是椭球面上的正常重力。这个 $N$ 叫大地水准面差距（Geoid Undulation，日语里叫「ジオイド高」）。它在全球范围内大约在 $-106\ \text{m}$ 到 $+85\ \text{m}$ 之间起伏。前面图 2 里面那几个凹凸不平的大土豆，其实就是对这种情况的一种夸张。

所以准确来说「高度」其实由两个独立的量组成：

• 椭球高 $h$：从参考椭球面沿法线量到你的距离。GNSS 直接解算出来的就是这个。
• 正高 / 标高 $H$：从大地水准面沿铅垂线量到你的距离。地图上的海拔、水准点的成果、你家门口的标高牌都是这个。

它们之间差的就是大地水准面差距： 严格来说这个式子是近似的：铅垂线和椭球面法线并不平行，两者的夹角叫垂线偏差（图 5 里 $O'$ 和 $O''$ 的分岔），在日本一般是几秒到十几秒。不过它带进 $h = H + N$ 里的误差在毫米级，可以忘掉。

$$h = H + N$$

日本的大地水准面差距大概在 30 到 40 米之间，东京湾一带约 36 米，山区能超过 40 米。国土地理院为此专门提供了格网模型 这里参考的是「日本のジオイド2011」。2025 年 4 月起换成了「ジオイド2024」。 。也就是说理论上 GPS 传感器直接测出来的高度和地图数据里的道路标高，天生就差了三四十米。

经纬度其实是一个方程的解

□ 图7. 大地纬度、经度、椭球高都是相对于「准拠椭球」这个人为规定的曲面定义的, 出典：国土地理院ウェブサイト

大地纬度 $\varphi$ 从定义上讲是过该点的椭球面法线与赤道面的夹角。这里面你用哪一种「椭球面」是人为规定的，相应的「法线」也是相对这个椭球面才有意义。所以我们平时描述坐标“东经 xxx 度，北纬 xxx 度”，其实并不严谨，因为我们没有说明自己用的是哪一个参考系。如果不在同一个坐标系上的话，那完全就是跨服聊天了。

把经纬度和椭球高还原成一个真正客观的、地心直角坐标 $(X, Y, Z)$，需要用到椭球的参数：

$$e^2 = \frac{a^2 - b^2}{a^2} = 2f - f^2$$ $$\qquad N(\varphi) = \frac{a}{\sqrt{1 - e^2\sin^2\varphi}}$$

值得一提的是 $N(\varphi)$ 是 卯酉圈曲率半径 （Radius of Curvature in the Prime Vertical），跟上面提到的 大地水准面差距 $N$ 不是一东西。 符号撞车是大地测量学的传统艺能。 然后可以求得地心直角坐标 $(X, Y, Z)$

$$\begin{cases} X = \left(N(\varphi) + h\right)\cos\varphi\cos\lambda \\ Y = \left(N(\varphi) + h\right)\cos\varphi\sin\lambda \\ Z = \left(N(\varphi)(1 - e^2) + h\right)\sin\varphi \end{cases}$$

这也可以从数学上说明，只有在给定了 $a$、$f$、以及这个椭球在空间中的摆放方式之后，才对应 $(\varphi, \lambda, h)$ 到一个确定的物理位置。换一个椭球，或者把同一个椭球平移一米，同一个物理点的经纬度数值就变了。

顺便看一个反面例子。如果你不用大地纬度而用地心纬度 $\psi$（就是真正从地心连线量出来的角度），两者的关系是：

$$\tan\psi = (1 - e^2)\tan\varphi$$

在 45 度附近两者相差最大，约 $f$ 弧度，也就是 $11.5'$、$0.19^\circ$，换算到地面上是差不多得有 21 公里。所以经纬度这个东西，脱离了它的定义体系就是三个没有意义的浮点数。

基准：确定椭球如何摆放

GRS 80 只是对地球形状做了数学建模，但没说椭球到底该摆在地球的什么位置上、椭球的三个轴该指向哪里。为此我们需要定义 大地基准（Geodetic Datum）。

老式的基准是「局部」的：挑一个天文点，实测它的经纬度和方位，然后让椭球贴合本地的大地水准面。日本的旧「日本測地系」（Tokyo Datum，epsg:4301）就是这么来的。这个基准采用的是 Bessel 1841 椭球（长半轴 6377397.155 米，比 GRS 80 小了 740 米），原点在东京天文台。它在日本本土贴得很好，但椭球中心离地球质心差了几百米，所以整体上旧日本測地系和世界測地系之间的区别很大，离谱的地方大概能有纬度 $+12''$、经度 $-12''$ 左右的系统性偏移，直接差上 400 多米。

现代的基准是「地心」的，靠 VLBI、SLR、GNSS 全球联测把椭球中心钉在地球质心上。这里的关键角色是 ITRF（International Terrestrial Reference Frame，国际地球参考框架）。

要注意 ITRS 和 ITRF 的区别：ITRS 是一套定义（原点在地球质心、尺度用 SI 米、朝向如何演化）；而 ITRF 是这套定义的一个实现，也就是具体的数值。他们会列出全球几百个观测站在某个参考历元的坐标和它们的速度。从 ITRF94、ITRF2008 到 ITRF2014、ITRF2020，每一版都是重新平差出来的一张新表。

WGS 84 则是美国国防部自己维护的另一套实现，它随着 GPS 监测站的重新解算也在不断更新，历史上有 G730（1994）、G873（1997）、G1150（2002）、G1674（2012）、G1762（2013）、G2139（2021）这些版本，版本号是解算所用的 GPS 周。每一版都在向同期的 ITRF 靠近，现在两者的一致性在厘米级。

于是就有了一个非常容易踩进去的坑：epsg:4326 并不是上面任何一个具体实现。EPSG 把它定义成一个「基准集合」（Datum Ensemble），意思是「WGS 84 的某个实现，我不告诉你是哪个」，它的标称定位精度是 2 米。所以当我说「我的 GPS 数据是 epsg:4326」的时候，我实际上说的是「这份数据大概落在真值 2 米以内」。

地球在动：元期与今期

□ 图9. 元期是什么：上图是某一瞬间测定并冻结下来的「测量成果」，规整的格网；下图是若干年后的现实，地壳变动把这张网拽变形了。坐标成果没变，地面变了, 出典：国土地理院ウェブサイト

要确定一个大地基准，需要定义四样东西：

1. 一个椭球体（$a$、$f$）
2. 原点（椭球中心和地球质心的关系）
3. 三个坐标轴的指向（$Z$ 轴对准哪个极、$X$ 轴对准哪条本初子午线）
4. 一个历元（Epoch）

前三个前面都解释过了，也比较好理解。但历元这玩意儿就很抽象了。

现在把日本的三个坐标系摆在一起看：

坐标系	EPSG	椭球	基于的框架	元期
旧 日本測地系	4301	Bessel 1841	东京天文台（局部）	—
JGD 2000	4612	GRS 80	ITRF94	1997.0
JGD 2011	6668	GRS 80	ITRF2008（东日本）/ 继承 JGD 2000	2011.4 / 1997.0
WGS 84 ensemble	4326	WGS 84	WGS 84（G730…G2139）	观测时刻

现代坐标系采取的椭球几乎完全相同。GRS 80 和 WGS 84 的椭球长半轴都是 6378137 米，扁率倒数分别是 298.257222101 和 298.257223563。小数点后第 6 位才有差别，换算到短半轴大约是 0.1 毫米。所以椭球体几乎不是误差来源。真正的差别在于元期。日本列岛在太平洋板块和菲律宾海板块的挤压下，相对于 ITRF 每年移动好几厘米。JGD 2000 的元期是 1997.0，也就是说它记录的是 1997 年年初那一瞬间日本各个基准点的位置。到今天已经过去快三十年，光是稳态的板块运动就累积了几十厘米。

所谓 元期 （reference epoch）是坐标成果所对应的那个时刻； 今期 （observation epoch）是你实际观测的时刻。GNSS 实测得到的永远是今期坐标，而地籍、地图、公共测量成果要求的是元期坐标。国土地理院为此提供「セミ・ダイナミック補正」的年度参数，把今期实测值折算回元期。简而言之元期是坐标成果是一张在某个瞬间按下快门的照片，而地面此后还在继续动。

然后众所周知的 2011 年 3 月 11 日，M9.0。牡鹿的观测点向东移动了约 5.3 米，同时下沉了约 1.2 米。这个量级已经不是「修正参数」能糊过去的了，国土地理院只能重新测量、重新发布成果，这就是 JGD 2011。

有意思的是 JGD 2011 是个「拼接」出来的基准：东日本和北陆地区用 ITRF2008、元期 2011.4 重新算过，其余地区直接继承了 JGD 2000 的成果。所以 JGD 2000 和 JGD 2011 在西日本和北海道几乎相同，在东北却能差好几米。

+ 追记：JGD 2024 的 EPSG code 在 2026 年 4 月的 EPSG Dataset v12.055 里补上了。因为水平位置没动，经纬度还是 6668、平面直角座標系还是 6669–6687，一个都没变；新增的只有垂直基准 JGD2024 (vertical) EPSG:11317，以及把它和各系水平坐标拼起来的复合 CRS。原来的 EPSG:6697（JGD2011 + JGD2011 (vertical) height）这类复合 CRS 相应被取代。所以只有带高度的数据才需要改 code，只存经纬度的表完全不受影响。

再往后，2025 年 4 月 1 日起日本又切到了 JGD 2024（測地成果2024）。这一版的水平位置直接继承 JGD 2011 没有变，改的是高度体系：标高不再靠水准网逐点传递，而是由 GNSS 的椭球高减去新的重力法大地水准面模型「ジオイド2024」得到，局部标高最大变动约 60 厘米。

两个坐标系之间的换算

知道了差别在哪，换算的思路也就清楚了。基准之间的转换本质上分两类。

第一类是刚体变换。把两边都化成地心直角坐标，然后套七参数的 Helmert 变换（Bursa–Wolf 模型）：三个平移、三个微小旋转、一个尺度因子。

$$\begin{pmatrix} X \\ Y \\ Z \end{pmatrix}_{B} = \begin{pmatrix} t_x \\ t_y \\ t_z \end{pmatrix} + (1 + s) \begin{pmatrix} 1 & -r_z & r_y \\ r_z & 1 & -r_x \\ -r_y & r_x & 1 \end{pmatrix} \begin{pmatrix} X \\ Y \\ Z \end{pmatrix}_{A}$$

因为旋转量是微角度，$\sin r \approx r$、$\cos r \approx 1$，旋转矩阵可以线性化成上面这个形式。ITRF 各版本之间的官方转换参数则是十四个：上面七个，外加它们各自的时间变化率，用的时候先按历元差把参数演化一遍：

$$p(t) = p(t_0) + \dot{p}\,(t - t_0)$$

这也是为什么谈 ITRF 之间的转换必须同时给出历元，不给历元的转换参数是没有意义的。

第二类是格网插值。地壳变动、地震位移这种东西不是刚体运动，一个七参数模型拟合不出来，只能实测一堆点然后做格网。

日本的 TKY2JGD.par（旧測地系到 JGD 2000）、touhokutaiheiyouoki2011.par（东北地震的 PatchJGD 修正）都是这种东西，用的时候在格网上做双线性插值。

实际写代码的时候这些通常都由 PROJ 处理，不过这里有一个坑：

from pyproj import Transformer

# always_xy=True：强制输入输出都按 (lon, lat) 排列
t = Transformer.from_crs("EPSG:6668", "EPSG:4326", always_xy=True)
print(t.description)

epsg:6668 到 epsg:4326 这一对，PROJ 会告诉你它选的是 Ballpark geographic offset from JGD2011 to WGS 84。「Ballpark」是 PROJ 的行话，意思是“它找不到任何有依据的转换，于是干脆什么都没做，只是把数值原样搬过去了，精度未知”。这是因为 epsg:4326 是个 2 米精度的 ensemble，而 JGD 2011 和它的差异本来就在这个 2 米以内，没有哪个转换是「正确」的。

所以这里必须根据使用场景做判断：如果你的容差是米级，Ballpark 的答案就是对的，把 6668 当 4326 用完全没问题。但如果你要的是厘米级，那就得放弃 epsg:4326 这个含糊的写法，指明具体的实现（比如 EPSG:9755，WGS 84 (G2139)）和历元，让 PROJ 走真正的时变转换。

另外那个 always_xy=True 很重要。epsg:4326 在 EPSG 的官方定义里轴序是纬度在前、经度在后，而 GeoJSON、WKT 实践、BigQuery、绝大多数 web 地图用的都是经度在前。pyproj 默认尊重 EPSG 的定义，于是不加这个参数就会悄悄地把你的坐标转置。

+ 讲道理轴序搞错在日本其实还行：日本的纬度在 24–46 度、经度在 122–154 度之间，一交换纬度就超出 $[-90, 90]$，ST_GEOGPOINT 会直接报错。在纬度和经度都小于 90 的地方（比如地中海沿岸），坐标就会在你不知不觉地时候跑到另一个半球。

写到这里我心里就直冒火，这种傻逼问题到处都是：我以前升级 Spark 版本，本来想顺带升级一下依赖的 org.datasyslab 的 GeoSpark 库，然后发现他们直接另起炉灶搞了 Sedona。我折腾半天把 Sedona 调通之后发现计算结果完全不对！我还寻思是怎么回事是不是我有问题，最后发现这个经纬度前后顺序它给我默认颠倒了。

回到最初的问题

绕回文章开头那三份数据：

• WGS 84（4326）：一个 ensemble，跟着 GPS 星历走，本质上是「今期」坐标，标称精度 2 米。
• JGD 2000（4612）：GRS 80 + ITRF94，冻结在 1997.0。
• JGD 2011（6668）：GRS 80 + ITRF2008 冻结在 2011.4（东日本），其余继承 1997.0。

三者两两之间的水平差异，在西日本是几十厘米量级，在东日本因为地震重测过反而更接近现在的真值。

至于让 BigQuery 上的结果和我本地跑出来的不一样的是定义的区别。BigQuery 的 GEOGRAPHY 类型定义里说：点是 WGS 84 椭球面上的位置，用经度和大地纬度表示，经度在前；而边是球面上的最短路径，即球面测地线。这里有个微妙的不一致：点定义在椭球面上，边却定义在球面上。关键是 GEOGRAPHY 里没有「直线」这个概念。

GeoJSON 规定的边是平面直线，BigQuery 在导入 GeoJSON 时必须做曲面化（tessellation），沿着原来的边加插值点。最后折线与原始边的偏差最大能有 **10 米。

而且你在平面上画的一条直线道路，转成经纬度再进 BigQuery，它就不是原来那条线了。短边无所谓，长边会偏得很离谱：从 (139°E, 35°N) 到 (140°E, 36°N) 这条 143 公里的斜线，「经纬度平面上的直线」和「球面测地线」在中间最多相差约 460 米。

所以 PostGIS 的 geometry 是平面的、geography 是球面的，ST_Distance 在两者上的语义不同。我本地用 geometry 跑通的逻辑，搬到 BigQuery 上语义就变了。而我当时以为是坐标系转错了，于是往完全错误的方向查了很久。关键是这种差异几乎无法体现在具体坐标数值上，但等到算距离、算面积、判断相交的时候悄悄给你整个大的。真他妈的傻逼啊。

总结

把这篇文章里出现过的所有误差来源按数量级排一下：

误差来源	量级
经纬度轴序搞反	$10^3$ km
把旧日本測地系当 WGS 84 用	~400 m
把平面上的直线当成球面测地线（或反之）	数百 m（143 km 的边约 460 m）
椭球高 vs 标高（大地水准面差距）	30–40 m（垂直方向）
GeoJSON 平面边 → 球面边的曲面化容差	≤ 10 m
消费级 GPS 在城市峡谷里的定位误差	5–50 m
epsg:4326 ensemble 的标称精度	2 m
JGD 2000/2011 元期与今期之间的地壳变动	几十 cm（震后局部数 m）
GRS 80 与 WGS 84 椭球体之差	0.1 mm

所以落到实践上需要注意这些地方：

1. 要根据使用场景决定是否进行参考系转换。比如做地图匹配、GPS 本身就有 5–10 米误差的时候，几十厘米的基准差异可以直接忽略；但如果要把轨迹 snap 到 3.5 米宽的车道，或者判断车辆是否进入了某个几十米的围栏，几十厘米就足以改变结论。
2. 如果要讨论高精度定位，每份数据都必须带 CRS 和历元。万一数据里没写可以自己验：拿几个已知地物对一下，看偏移是系统性的 400 米还是随机的 10 米，基本就能猜出来源。
3. 写入数据库之前要统一参考系。可以全部转成 epsg:4326。
4. 经纬度要分成两列分开存 (lon, lat)。而且使用地理信息库的时候 一定要注意经纬度顺序。
5. 椭球高和标高也要分成两列存。
6. 需要面积、缓冲、平行线偏移这些需要严格平面几何的时候，先投影到平面直角座標系再算；BigQuery 里只能处理球面几何的问题。