坐标系与矩阵(4)：球心坐标与NEU坐标系

前三篇介绍了坐标系和矩阵的数学知识，从本篇开始，我们试图运用这些知识来解决实际问题。

如上图，模拟了一个以球心为原点的固定坐标系，该坐标系有一个名称地心地固坐标系(ECEF)，对应我们之前介绍的坐标系

，而平面场景在我们生活中更为直观，上北下南，左东右西，对应上图中绿色的切平面，简称NEU坐标系，对应之前介绍的坐标系

。于是，给定一点

，我们需要计算一个矩阵

，实现两个坐标系的转换。

这里对应两个环节，(1)球心坐标系的单位换算， 从经纬度

到米单位的笛卡尔坐标

；(2)从ECEF到NEU，从全球坐标系

到本地坐标系

。

整体来看，默认初始时

，从F到M大概需要三个过程：(1)沿

逆时针旋转

，如上图，橙色对应的是

，红色对应的是

，方向均向内；(2)沿着新坐标系中的红轴逆时针旋转

；（3）沿新坐标系的

方向平移到绿色坐标系的原点。

前两个旋转矩阵对应的是：

这样，只要知道平移

对应的三个分量，可以轻松的得到最终的矩阵。然后这里就有一个问题，地球是椭球而不是圆球，如下图所示，维度是

，但准确的切平面对应的法线是

，前者是geocentric normal，后者是geodetic normal，通常情况下，我们都是指的

。

椭球下经纬度和笛卡尔坐标转换的问题暂时先不做处理，假设已知某点对应球心的位置

，椭球在X，Y，Z方向的半径分别为a, b，c。则，up方向也就是该点对应椭球切面的法线(geodetic normal)：

又因为地球对应的椭球体中，a==b。North轴指向北极，而East轴和Up轴，North轴正交，可得East轴必须垂直于

和

：

而

，因此，我们可以获取ENU坐标系三个轴的向量

，这样，对应的转换公式为：

这样，我们在ENU本地坐标系上的一点

，对应球心坐标系上的点

，满足：

如上，我们实现了ECEF和ENU之间的转化，下面，我们讲一下经纬度到ECEF之间的转换，该问题可以抽象为已知经纬度+高度

，这里的

对应ECEF坐标系下的

，求其在笛卡尔坐标系下对应的位置

，这块并不涉及到坐标系的概念，不感兴趣的可以略过。

不失一般性，我们认为地球的椭球为：

因为

是geodetic normal，所以，在点

处的法线为

，同样，由椭圆的属性可知，该点的法线（未归一化）同样可以是

，则可得如下关系：

将

三个未知变量放到等式左侧：

同时，因为

在椭球体上，带入椭球体的公式，可得：

此时，我们可以得出

三个值，也就是

点的未知，最终就可以得到绿色点的位置：

接下来，自然就是逆运算了，已知

，如何求

。

如上图，一旦我们知道

：

已知

，可得：

同样的思路，将

三个未知变量放到等式左侧：

这里，我们定义一个函数S：

这是求解一元四次方程，我觉得应该可以直接公式计算解析解吧？不过标准方式是采用的牛顿迭代法。初始值就是以geocentric和

连线与椭球体的焦点，此时

：

这样，我们可以计算函数S相对于

的导数，用牛顿迭代法不断逼近，找到满足自己要求的

值：

我个人对这个迭代对应的几何意义理解大概如下，迭代调整α，不断的缩小alpha*n_s的距离：

Iteration 1：初始值，得到

Iteration 2：以

的geodetic normal方向进行scale，得到某个未知，连接r，得到

Iteration 3：同理，不断逼近，得到满足容限的近似解

本文到此结束，主要介绍了球心坐标系转换的相关内容，因为涉及到椭球而变得有些复杂，但在坐标系转换这方面并不复杂。下一篇介绍运动学中常用的Denavit-Hartenberg Algorithm。

参考资料：Cesium

3D Engine Design for Virtual Globes第二章