椭球面系列---大地坐标和笛卡尔坐标的相互转换

椭球体下，尤其是地球的旋转椭球体下，大地坐标和笛卡尔坐标的相互转换是最基础的算法了。本章给出两种坐标系下相互转换的原理及相应的转换公式，供参考。

大地坐标

大地坐标（Geodetic coordinate）是大地测量中以参考椭球面为基准面的坐标，点P的位置用大地经度λ、大地纬度φ和大地高H表示。

大地坐标多应用于大地测量学，测绘学等。具体为：

• 大地经度
  大地经度是通过该点的大地子午面与起始大地子午面（通过格林尼治天文台的子午面）之间的夹角。规定以起始子午面起算，向东由0°至180°称为东经；向西由0°至180°称为西经。
• 大地纬度
  大地纬度是P点在椭球面的投影点的法线与赤道面的夹角，规定由赤道面起算，由赤道面向北从0°至90°称为北纬；向南从0°到90°称为南纬。P点位于椭球面的投影点的法线方向上。
• 大地高度
  大地高是地面点沿法线到参考椭球面的距离。

注意大地纬度与地心纬度的区别！大地坐标的示意图如下。

笛卡尔坐标

点P的笛卡尔坐标即为参考椭球体中心直角坐标系下的坐标，使用$(x,y,z)$表示。

大地坐标$(\lambda ,\phi ,h)$转换为笛卡尔坐标$(x,y,z)$

大地坐标向笛卡尔坐标的转换为直接转换，不需要迭代。

见上图，已知$P$点的大地坐标$(\lambda ,\phi ,h)$，$P^{\prime }$为$P$点在椭球面上的投影点，即$P^{\prime }P$为$P^{\prime }$点的法线方向(定义为$\text{n}$)。

法线向量$\text{n}$很容易由大地经度和大地纬度计算得到：
$$\begin{array}{c}\text{n}=\left[\begin{array}{c}cos(\phi )cos(\lambda )\\ cos(\phi )sin(\lambda )\\ sin(\phi )\end{array}\right]\end{array}\text{\hspace{1em}?}$$
在椭球面系列—基本性质一文中我们知道，若已知椭球面上点$P^{\prime }$的法线向量$\text{n}$，则可求解$k$:
$$\begin{array}{c}{k}^2={\text{n}}^T{\text{C}}^{?1}\text{n}\end{array}\text{\hspace{1em}?}$$
从而可以直接计算得到椭球面上点$P^{\prime }$的笛卡尔坐标$(x^{\prime },y^{\prime },z^{\prime }{)}^T$：
$$\begin{array}{c}{P}^{\prime }=(x^{\prime },y^{\prime },z^{\prime }{)}^T=\frac{1}{k}{\text{C}}^{?1}\text{n}\end{array}\text{\hspace{1em}?}$$
得到点$P^{\prime }$后，则由法向量$\text{n}$和高度$h$可直接得到点$P$的笛卡尔坐标$(x,y,z{)}^T$。
$$\begin{array}{c}P=\left[\begin{array}{c}x\\ y\\ z\end{array}\right]=\left[\begin{array}{c}{x}^{\prime }\\ y^{\prime }\\ z^{\prime }\end{array}\right]+\frac{\text{n}}{|\text{n}|}h\end{array}\text{\hspace{1em}?}$$

笛卡尔坐标$(x,y,z)$转换为大地坐标$(\lambda ,\phi ,h)$

笛卡尔坐标向大地坐标的转换为隐式转换，需要迭代求解。

仍然见上图，$P^{\prime }$点的法线方向(定义为$\text{n}$)由其笛卡尔坐标$(x^{\prime },y^{\prime },z^{\prime }{)}^T$表示（这里我们不知道大地坐标，所以不能像式(1)方式来表示）则有：
$$\begin{array}{c}\text{n}=\left[\begin{array}{c}{x}^{\prime }/a^2\\ y^{\prime }/b^2\\ z^{\prime }/c^2\end{array}\right]\end{array}\text{\hspace{1em}?}$$
很明显，此法线$\text{n}$为梯度向量的$1/2$。

$P^{\prime }$点和$P$点的关系可表示为（与式(4)形式相同）：
$$\begin{array}{c}\left[\begin{array}{c}x\\ y\\ z\end{array}\right]=\left[\begin{array}{c}{x}^{\prime }\\ y^{\prime }\\ z^{\prime }\end{array}\right]+t\left[\begin{array}{c}{x}^{\prime }/a^2\\ y^{\prime }/b^2\\ z^{\prime }/c^2\end{array}\right]\end{array}\text{\hspace{1em}?}$$
上式中，$t$为法向量的系数，选取合适的$t$，则可使得上式成立。

$P^{\prime }$点和$P$点的关系为（与式(4)形式相同）：
$$\begin{array}{c}\left[\begin{array}{c}x\\ y\\ z\end{array}\right]=\left[\begin{array}{c}{x}^{\prime }\\ y^{\prime }\\ z^{\prime }\end{array}\right]+t\left[\begin{array}{c}{x}^{\prime }/a^2\\ y^{\prime }/b^2\\ z^{\prime }/c^2\end{array}\right]\end{array}\text{\hspace{1em}?}$$
将上式形式改变为：
$$\begin{array}{c}\left[\begin{array}{c}{x}^{\prime }\\ y^{\prime }\\ z^{\prime }\end{array}\right]=\left[\begin{array}{c}x/(1+t/a^2)\\ y/(1+t/b^2)\\ z/(1+t/c^2)\end{array}\right]\end{array}\text{\hspace{1em}?}$$
由于$P^{\prime }$为椭球面上的点，因此满足椭球面方程：
$$\begin{array}{c}\frac{x^{\prime 2}}{a^2}+\frac{y^{\prime 2}}{b^2}+\frac{z^{\prime 2}}{c^2}=1\end{array}\text{\hspace{1em}?}$$
将式(8)带入上式，并定义函数$f(t)$:
$$\begin{array}{c}f(t)=\frac{x^2}{a^2(1+t/a^2{)}^2}+\frac{y^2}{b^2(1+t/b^2{)}^2}+\frac{z^2}{c^2(1+t/c^2{)}^2}?1\end{array}\text{\hspace{1em}?}$$
由于$P$点坐标$(x,y,z)$已知，则变为求解方程$f(t)=0$的根。

函数$f(t)$为一元函数，且为隐式函数，因此求根一般采用牛顿迭代法(使用一阶导数$f^{\prime }(t)$)，即：
$$\begin{array}{c}f(t)=f(t_0)+f^{\prime }(t_0).\Delta t=0\end{array}\text{\hspace{1em}?}$$
迭代求解时，每一步$t$可由前一次数值给出：
$$\begin{array}{c}{t}_{n+1}=t_n?f(t_n)/f^{\prime }(t_n)\end{array}\text{\hspace{1em}?}$$
由式(10)可得一阶导数$f^{\prime }(t)$:
$$\begin{array}{c}{f}^{\prime }(t)=\frac{?2x^2}{a^4(1+t/a^2{)}^3}?\frac{2y^2}{b^4(1+t/b^2{)}^3}?\frac{2z^2}{c^4(1+t/c^2{)}^3}\end{array}\text{\hspace{1em}?}$$

下面给出迭代参数$t$的初值$t_0$的求解。

由上图可知，$OP$与椭球面的交点$P_0$与$P^{\prime }$点很近，因此首先求解$P_0$点的坐标$(x_0,y_0,z_0)$：
$$\begin{array}{c}\left[\begin{array}{c}{x}_0\\ y_0\\ z_0\end{array}\right]=\frac{1}{d}.\left[\begin{array}{c}x\\ y\\ z\end{array}\right]\end{array}\text{\hspace{1em}?}$$
上式中$d$可由$P$点的坐标$(x,y,z)$求得：
$$\begin{array}{c}\frac{x^2}{a^2}+\frac{y^2}{b^2}+\frac{z^2}{c^2}=d^2\end{array}\text{\hspace{1em}?}$$

迭代开始时，将$P_0$当作$P^{\prime }$即可。则根据上图几何关系和式(7)，有：
$$h=(1?1/d)|\text{P}|=t|\text{n}|$$
因此$t_0$为：
$$\begin{array}{c}t=(1?1/d)|\text{P}|/|\text{n}|\end{array}\text{\hspace{1em}?}$$
上式中：
$$|\text{P}|=\sqrt{x^2+y^2+z^2}$$
$$|\text{n}|=\sqrt{x_0^2/a^4+y_0^2/b^4+z_0^2/c^4}$$
使用牛顿迭代法求得$t$后，则带入式(8)可求得椭球面投影点于$P^{\prime }$的坐标$(x^{\prime },y^{\prime },z^{\prime }{)}^T$。

有了$P^{\prime }$点坐标，则椭球面法向量$\text{n}$可由式(5)给出。

最终，大地坐标为：
$$\begin{array}{c}?\begin{array}{l}\lambda =ta{n}^{?1}(n_y,n_x)\\ \phi =si{n}^{?1}(n_z/|\text{n}|)\\ h=t/|\text{n}|\end{array}\end{array}\text{\hspace{1em}?}$$
注意，上式中的$|\text{n}|$为$\sqrt{x^{\prime 2}/a^4+y^{\prime 2}/b^4+z^{\prime 2}/c^4}$。

参考:
[1]: 椭球面系列—基本性质
[2]: [学习内容：求一个点到椭球面的距离（下）