ECEF和大地坐标系的相互转化

在阅读 RTKLIB的源码时，发现了ECEF和大地坐标系的相互转换的函数，大地坐标系（φ,λ,h）转成ECEF（X,Y,Z）与所看书籍（GPS原理与接收机，谢刚，电子工业出版社）的公式是一样的，而ECEF转成大地坐标系的公式则与上述书籍和 RTKLIB的使用手册（P135）中的公式都不一样。简而言之，遇到的问题如下图所示。

一开始想着能不能直接在某个文献中直接找到与程序算法一致的处理流程，但找了很久都没能如愿。最后，只能想着彻底从平面几何上理解这个数学推导过程，从而再尝试理解为什么程序会那么写。

在查找介绍大地坐标系的文献中，总能看到酉卯圆曲率半径这个名词。百度百科上给的定义是：过椭球面上一点的法线，可作无限个法截面，其中一个与该点子午面相垂直的法截面同椭球面相截形成的闭合的圈称为卯酉圈，如下图中的PEE′所示。

图 1

投影到二维平面中，就是下面这幅图。其中，Pn就称为酉卯圆曲率半径，PT则为椭圆的切线，角B则是酉卯圆曲率半径与长轴的夹角，称为纬度（也就是说通常意义上某个点的纬度并不是该点到地心的连线与长轴的夹角！）。

图 2

由椭圆切线的斜率公式，可知

                                                                                                                                                                                       (1)

联立上面两式，可得

                                                                                                                                                                                                    (2)

其中，

。再将上式带入到椭圆的标准方程

中，即可解得

                                                                                                                                                                    (3)

另外，在△PnA中通过简单的几何关系可知，

。所以有

                                                                                                                                                                                             (4)

又因为

，所以有

                                                                                                                                                                                    (5)

注意，上面这个关于PQ、Qn的公式在后面推导ECEF和大地坐标系的相互转换公式时，是很重要的。

为了说明坐标转换的一般性，这里的待转换点P是在地球外部的，示意图如下。图 3在理解（φ,λ,h）→XYZ中的X、Y坐标的变换公式时较为重要。

图 3

而要想真正简单、清晰地理解大地坐标系（φ,λ,h）与ECEF（X,Y,Z）如何转换，还是要从二维平面图中着手。下图可以看作是椭球体的主视图，XY平面压缩成了水平X轴。

图 4

1. 大地坐标系（φ,λ,h）转成ECEF（X,Y,Z）。

OD是点P在XY平面上的投影向量，其长度为

                                                                                                                                                                        (6)

再结合图3中XY平面中经度角λ的示意图，可得

                                                                                                                                                                                       (7)

而关于z坐标，则有

                                                                                                                                                           (8)

1. ECEF（X,Y,Z）转成大地坐标系（φ,λ,h）。

首先是经度角λ，在图3的XY平面中可以很清晰地看到

                                                                                                                                                                                              (9)

对于纬度角φ和高度h，可以在在图4中分析其几何关系。

                                                                                                                                                                             (10)

其中，

                                                                                                                                                                      (11)

从上述公式中可以看到，在计算φ时还是会用到纬度角φ的，所以不能直接用上述公式来计算纬度角。这里给出RTKLIB中ecf2pos函数中的算法步骤：

①　假设PD=PE，计算出夹角φ’。

②　根据式(4)和式(11)再计算N’和PE’。

③　对比②中的PE’与①中的PE之间的差别是否小于截断因子。不符合条件，将②中的PE’代入到①中，继续循环计算；符合条件，则说明此时φ’=φ，可以跳出循环。

④　将最终的PE代入到式(10)中，从而计算出纬度角φ和高度h。

至此，关于RTKLIB的源码中的pos2ecf和ecf2pos函数的计算过程的理论公式推导过程，就证明完了。

这个事情证明，不会推公式的算法程序员是不合格的。找到的现成理论依据在复杂编程问题时可能会不够高效，或者当算法流程稍作改变时就会看不懂别人写的程序了。