经纬高度转xyz坐标

hbghlyj

球上的经纬度
椭球上的经纬度 子午线的横截面是一个椭圆。 纬度是通过观察北极星来测量的，如右图。

坐标转换公式
From geodetic to ECEF coordinates\begin{aligned}X&=\left(N(\phi )+h\right)\cos {\phi }\cos {\lambda }\\Y&=\left(N(\phi )+h\right)\cos {\phi }\sin {\lambda }\\Z&=\left({\frac {b^{2}}{a^{2}}}N(\phi )+h\right)\sin {\phi }\\&=\left((1-e^{2})N(\phi )+h\right)\sin {\phi }\\&=\left((1-f)^{2}N(\phi )+h\right)\sin {\phi }\end{aligned}其中$$N(\phi )={\frac {a^{2}}{\sqrt {a^{2}\cos ^{2}\phi +b^{2}\sin ^{2}\phi }}}={\frac {a}{\sqrt {1-e^{2}\sin ^{2}\phi }}},$$$a$和$b$分别是赤道半径（半长轴）和极半径（半短轴）。
$e$ 是离心率, $ e^{2}=1-{\frac {b^{2}}{a^{2}}} $.
Javascript 经纬高度转xyz坐标

性质
以下以与地心坐标系相同的方式适用于经度：
$$ {\frac {X}{\cos \lambda }}-{\frac {Y}{\sin \lambda }}=0$$
以下适用于纬度：$$ {\frac {p}{\cos \phi }}-{\frac {Z}{\sin \phi }}-e^{2}N(\phi )=0,$$
其中$p={\sqrt{X^{2}+Y^{2}}}$. 为了消去参数$h$，减去
$$ {\frac {p}{\cos \phi }}=N+h $$
和$$ {\frac {Z}{\sin \phi }}={\frac {b^{2}}{a^{2}}}N+h$$
此外：$$ \tan \phi =(Z/p)/(1-e^{2}N/(N+h))$$通过微分确认坐标的正交性：\begin{aligned}{\begin{pmatrix}dX\\dY\\dZ\end{pmatrix}}&={\begin{pmatrix}-\sin \lambda &-\sin \phi \cos \lambda &\cos \phi \cos \lambda \\\cos \lambda &-\sin \phi \sin \lambda &\cos \phi \sin \lambda \\0&\cos \phi &\sin \phi \\\end{pmatrix}}{\begin{pmatrix}dE\\dN\\dU\end{pmatrix}},\\[3pt]{\begin{pmatrix}dE\\dN\\dU\end{pmatrix}}&={\begin{pmatrix}\left(N(\phi )+h\right)\cos \phi &0&0\\0&M(\phi )+h&0\\0&0&1\\\end{pmatrix}}{\begin{pmatrix}d\lambda \\d\phi \\dh\end{pmatrix}},\end{aligned}其中$$ M(\phi )={\frac {a\left(1-e^{2}\right)}{\left(1-e^{2}\sin ^{2}\phi \right )^{\frac {3}{2}}}} $$
（另见“椭球上的子午线弧”）。

hbghlyj

下面的“Bowring's irrational geodetic-latitude equation”，简单地从上述性质推导出来，可以通过 Newton-Raphson 迭代法有效求解：\[\kappa -1-{\frac {e^{2}a\kappa }{\sqrt {p^{2}+\left(1-e^{2}\right)Z^{2}\kappa ^{2}}}}=0,\]where $ \kappa ={\frac {p}{Z}}\tan \phi $ and $ p={\sqrt {X^{2}+Y^{2}}} $ as before. The height is calculated as:
\begin{aligned}h&=e^{-2}\left(\kappa ^{-1}-{\kappa _{0}}^{-1}\right){\sqrt {p^{2}+Z^{2}\kappa ^{2}}},\\\kappa _{0}&=\left(1-e^{2}\right)^{-1}.\end{aligned}The iteration can be transformed into the following calculation:
$$ \kappa _{i+1}={\frac {c_{i}+\left(1-e^{2}\right)Z^{2}\kappa _{i}^{3}}{c_{i}-p^{2}}}=1+{\frac {p^{2}+\left(1-e^{2}\right)Z^{2}\kappa _{i}^{3}}{c_{i}-p^{2}}}, $$where $ c_{i}={\frac {\left(p^{2}+\left(1-e^{2}\right)Z^{2}\kappa _{i}^{2}\right)^{\frac {3}{2}}}{ae^{2}}}. $
The constant $ \,\kappa _{0} $ is a good starter value for the iteration when $ h\approx 0 $. Bowring showed that the single iteration produces a sufficiently accurate solution. He used extra trigonometric functions in his original formulation.

hbghlyj

由上式推导出的$\kappa$的四次方程，可用“Ferrari法”解得：\begin{aligned}\zeta &=\left(1-e^{2}\right){\frac {z^{2}}{a^{2}}},\\[4pt]\rho &={\frac {1}{6}}\left({\frac {p^{2}}{a^{2}}}+\zeta -e^{4}\right),\\[4pt]s&={\frac {e^{4}\zeta p^{2}}{4\rho ^{3}a^{2}}},\\[4pt]t&={\sqrt[{3}]{1+s+{\sqrt {s(s+2)}}}},\\[4pt]u&=\rho \left(t+1+{\frac {1}{t}}\right),\\[4pt]v&={\sqrt {u^{2}+e^{4}\zeta }},\\[4pt]w&=e^{2}{\frac {u+v-\zeta }{2v}},\\[4pt]\kappa &=1+e^{2}{\frac {{\sqrt {u+v+w^{2}}}+w}{u+v}}.\end{aligned}The application of Ferrari's solution
A number of techniques and algorithms are available but the most accurate, according to Zhu, is the following procedure established by Heikkinen, as cited by Zhu. It is assumed that geodetic parameters $ \{a,\,b,\,e\} $ are known
\begin{aligned}a&=6378137.0{\text{ m. Earth Equatorial Radius}}\\[3pt]b&=6356752.3142{\text{ m. Earth Polar Radius}}\\[3pt]e^{2}&={\frac {a^{2}-b^{2}}{a^{2}}}\\[3pt]e'^{2}&={\frac {a^{2}-b^{2}}{b^{2}}}\\[3pt]p&={\sqrt {X^{2}+Y^{2}}}\\[3pt]F&=54b^{2}Z^{2}\\[3pt]G&=p^{2}+\left(1-e^{2}\right)Z^{2}-e^{2}\left(a^{2}-b^{2}\right)\\[3pt]c&={\frac {e^{4}Fp^{2}}{G^{3}}}\\[3pt]s&={\sqrt[{3}]{1+c+{\sqrt {c^{2}+2c}}}}\\[3pt]k&=s+1+{\frac {1}{s}}\\[3pt]P&={\frac {F}{3k^{2}G^{2}}}\\[3pt]Q&={\sqrt {1+2e^{4}P}}\\[3pt]r_{0}&={\frac {-Pe^{2}p}{1+Q}}+{\sqrt {{\frac {1}{2}}a^{2}\left(1+{\frac {1}{Q}}\right)-{\frac {P\left(1-e^{2}\right)Z^{2}}{Q(1+Q)}}-{\frac {1}{2}}Pp^{2}}}\\[3pt]U&={\sqrt {\left(p-e^{2}r_{0}\right)^{2}+Z^{2}}}\\[3pt]V&={\sqrt {\left(p-e^{2}r_{0}\right)^{2}+\left(1-e^{2}\right)Z^{2}}}\\[3pt]z_{0}&={\frac {b^{2}Z}{aV}}\\[3pt]h&=U\left(1-{\frac {b^{2}}{aV}}\right)\\[3pt]\phi &=\arctan \left[{\frac {Z+e'^{2}z_{0}}{p}}\right]\\[3pt]\lambda &=\operatorname {arctan2} [Y,\,X]\end{aligned}Note: arctan2[Y, X] is the four-quadrant inverse tangent function.

hbghlyj

对于小 $e^2$ 幂级数 $ \kappa =\sum _{i\geq 0}\alpha _{i}e^{2i}$ 开始以
\begin{aligned}\alpha _{0}&=1;\\\alpha _{1}&={\frac {a}{\sqrt {Z^{2}+p^{2}}}};\\\alpha _{2}&={\frac {aZ^{2}{\sqrt {Z^{2}+p^{2}}}+2a^{2}p^{2}}{2\left(Z^{2}+p^{2}\right)^{2}}}.\end{aligned}