Laurenfrost's Blog.

今朝有酒今朝醉，明日愁来明日愁。

A Slight Diffs on Coordinate

聊聊地理坐标系

Oct 7, 2025 2 min read
地理信息系统

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

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

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

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

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

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 $$

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

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

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

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

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

那怎么办呢，只能人为规定一个比较合理而且好算的一个等势面当作真实地球的重力基准面了。我们称之为大地水准面 (Geoid)。而 GRS 80 就是在这一整套思想指导下的地球等势面椭球模型。

想象一下站在这个旋转的椭球表面上。你感受到两个力：万有引力 + 离心力。

为了让椭球表面成为一个等位势面，重力方向必须处处垂直于椭球表面。如果不是垂直的，物体就会有一个沿着表面的分力，从而“滑向”位势更低的地方，这就不是一个平衡的等位势面了。

这里就引入了克莱罗定理 (Clairaut’s Theorem)。定理指出，对于一个处于流体静力平衡状态（可以想象成一个液态的、旋转的星球，最终形成的稳定形状）的旋转天体，它的几何形状的扁平程度与其表面重力变化的扁平程度之间的关系为：

$$ f + \beta \approx \frac{5}{2}m $$

$f$ (几何扁率 Geometric Flattening) 描述椭球在几何上有多扁。 $$ f = \frac{a-b}{a} $$ 其中 $a$ 是赤道半径，$b$ 是极半径。
$\beta$ (重力扁率 Gravity Flattening) 这个参数描述了表面重力值从赤道到两极的变化程度。 $$ \beta = \frac{g_p - g_e}{g_e} $$ 其中 $g_p$ 是在极点的重力加速度，$g_e$ 是在赤道的重力加速度。由于离心力在两极为零、在赤道最大，同时两极离地心更近，所以极点的重力 $g_p$ 会大于赤道的重力 $g_e$。
$m$ (自转影响因子) 代表在赤道上，离心力与引力（或重力）的比值。它量化了“自转”这个物理过程的影响力。 $$ m \approx \frac{\omega^2 a}{g_e} \approx \frac{\omega^2 a^3}{GM} $$

这揭示了一个深刻的内在联系：一个星球的形状、自转和表面重力场这三者是紧密耦合、相互制约的。你不能随意改变其中一个而不影响其他两个。如果一个星球转得更快（$m$ 变大），那么为了保持平衡，它的几何形状和重力场就必须变得更扁（$f$ 和 $β$ 相应增大）。

观察地球上某一点 P，设 P 点当前的纬度为 $\phi$。那么根据克莱罗定理可以得到推导公式索米里安纳公式 (Somigliana’s formula)：

$$ g_{\phi} = g_e \frac{1+k\sin^2{\phi}}{\sqrt{1 - e^2 \sin^2{\phi}}} $$

$g_\phi$ 是在纬度 $\phi$ 的正常重力值。
$g_e$ 是在赤道的正常重力值。
$e^2 = 2f - f^2$ 是椭球第一偏心率的平方。
$k$ 是一个由椭球基本参数（$a$, $b$, $g_e$, $g_p$）决定的常数，称为正常重力常数。

ITRF 94 坐标系是国际大地测量学和地球物理学联合会（IUGG）、国际大地测量学协会（IAG）、国际天文学会（IAU）决定建立的。并且由国际地球自转局 IERS（International Earth Rotation Service）准确及时提供自转参数，保持这个地球参考框架。

が構築した３次元直交座標系です。この座標系では、地球の重心に原点を置き、X軸を本初子午線と赤道との交点方向に、Y軸を東経90度の方向に、Z軸を北極（地軸の北端）の方向にとって、空間上の位置をX、Y、Zの数字の組で表現します。 ※国際地球基準座標系（ITRF）ではIERS基準子午線を本初子午線としています。IERS基準子午線は、グリニッジ子午線（英国のグリニッジ天文台のエアリー子午環を通る子午線）の102mほど東を通過します。

　Geodetic Reference System（測地基準系）1980の略語です。GRS80は、IAG（International Association of Geodesy：国際測地学協会）及びIUGG（International Union of Geodesy and Geophysics：国際測地学及び地球物理学連合）が1979年に採択した、地球の形状、重力定数、角速度等地球の物理学的な定数及び計算式です。GRS80では、楕円体の形状や軸の方向及び地球重心を楕円体の原点とすることも定められています。この楕円体をGRS80楕円体といいます。

国际上广泛使用的 WGS 84。

SELECT
  ST_INTERSECTS(
    ST_GeomFromWKT('POLYGON((139.85625 37.6, 139.85625 37.6041666666667, 139.85 37.6041666666667, 139.85 37.6, 139.85625 37.6))'),
    ST_GeomFromWKT('MULTILINESTRING ((139.85187100237266 37.59986200415802, 139.85148000514474 37.600000007507305))')
  ) AS do_they_intersect;