地球形状与曲率半径

\qquad 实际的地球表面是一个凹凸不平、形状十分复杂的物理面，难以准确量化描述。为了方便研究，通常以平静的海平面为基准，并把它向大陆延伸形成一个封闭曲面，称之为大地水准面。它所包含的的几何形状称为大地水准体。地球的实际形状如下图所示：
\qquad 实际应用中总是希望可以使用相对简单的数学形式来近似拟合地球几何形状，常用的近似有三种：
\qquad 1.近似成圆球体， 中心在地球质心，半径6371km，适用于精度要求不高的场景。
\qquad 2.近似成旋转椭球体，地球自转轴与与椭圆短半轴重合，椭圆度$\frac{1}{300}$(长短半轴相差21km)，椭圆绕其短半轴旋转一周构成的旋转椭球体的表面，该描述中地区吃到近似为原型。
\qquad 3.近似成三轴椭球体，其中赤道是椭圆的，赤道椭圆度约为$\frac{1}{100000}$(长短半轴相差约60m)。
\qquad 三轴椭球体虽然描述更精确，但相关计算更加复杂，其赤道的椭圆度并不大，可近似为圆形，因此旋转椭球体的应用最为广泛。
\qquad 包含南北极点的平面称为子午面，子午面与旋转椭球面的交线成为子午圈(经圈)，经过英国格林尼治的经线称为本初子午线(零度经线)。任一经线所在子午面与本初子午面之间的夹角定义为经度$(\varphi )$，夹角方向与地球自转方向相同，取值范围$-180°\sim 180°$。包含旋转椭球中心且垂直自转轴的平面称为赤道面，赤道面与旋转椭球面的交线称为赤道。平行于赤道的平面与旋转椭球面的交线成为纬圈。
\qquad 对于地球旋转椭球体而言，确定其三维形状参数的关键在于确定二维子午圈椭圆。

1.地球几何形状

\qquad 建立地心右手坐标系，即地心地固坐标系$(ECEF)$，坐标原点在地心，与地球固联，跟随地球自转一起相对于惯性空间转动。取子午圈椭圆上任意一点$P$与地心$O$的连线$PO$与$O{X}_e$的夹角成为地心纬度，记为$\varphi$，取值范围$-90°\sim 90°$。南纬为负，北纬为正。过$P$点的法线$PQ$与$O{X}_e$的夹角称为地理纬度。记为$L$，取值范围$-90°\sim 90°$。$PO$称为地心垂线，$PQ$称为地理垂线。
\qquad 椭圆方程为：
$$\frac{x^2}{R_e^2}+\frac{z^2}{R_p^2}=1\text{\hspace{1em}(1.1)}$$ \qquad 其中$R_e,R_p$分别为椭圆的长半轴和短半轴。
\qquad 椭圆扁率(或称椭圆度，flattening)定义为：
$$f=\frac{R_e-R_p}{R_e}\text{\hspace{1em}(1.2)}$$ \qquad 椭圆偏心率定义为：
$$e=\frac{\sqrt{R_e^2-R_p^2}}{R_e}\text{\hspace{1em}(1.3)}$$ \qquad 第二偏心率定义为：
$$e^{\prime }=\frac{\sqrt{R_e^2-R_p^2}}{R_p}\text{\hspace{1em}(1.4)}$$ \qquad 且有$$e^{\prime 2}=\frac{e^2}{1-e^2}\text{\hspace{1em}(1.5)}$$ \qquad 由式(1.2)和式(1.3)分别可得：
$$R_p=(1-f)R_e\text{\hspace{1em}(1.6)}$$ $$R_p=R_e\sqrt{1-e^2}\text{\hspace{1em}(1.7)}$$ \qquad 比较式(1.6)和式(1.7)可得：
$$f=1-\sqrt{1-e^2}\text{\hspace{1em}(1.8)}$$ $$e^2=2f-f^2\text{\hspace{1em}(1.9)}$$ \qquad 将椭圆方程两端对$x$求导，并将式(7)带入得：$$\frac{2x}{R_e^2}+\frac{2z\cdot \frac{dz}{dx}}{R_e^2(1-e^2)}=0\text{\hspace{1em}(1.10)}$$ \qquad 整理得：$$\frac{dz}{dx}=-(1-e^2)\frac{x}{z}\text{\hspace{1em}(1.11)}$$ \qquad 式中$\frac{dz}{dx}$表示椭圆在$P$点的切线$PB$的斜率。切线PB与法线PQ相互垂直，法线PQ的斜率为$tanL$，则有$$\frac{dz}{dx}tanL=-(1-e^2)\frac{x}{z}tanL=-1\text{\hspace{1em}(1.12)}$$ 由式(1.12)可得：$$z=x(1-e^2)tanL\text{\hspace{1em}(1.13)}$$ \qquad 将式(1.7)和式(1.13)带入椭圆方程，可求得以地理纬度为参数的椭圆方程：$$\{\begin{array}{ll}x& =\frac{R_e}{\sqrt{1-e^{2}si{n}^{2}L}}cosL\\ z& =\frac{R_e(1-e^2)}{\sqrt{1-e^{2}si{n}^{2}L}}sinL\end{array}\text{\hspace{1em}(1.14)}$$ \qquad 记线段长度$\overline{PQ}=R_N$，则有：$$x=R_{N}cosL\text{\hspace{1em}(1.15)}$$ \qquad 比较式(1.14)和式(1.15)可得：$$R_N=\frac{R_e}{\sqrt{1-e^{2}si{n}^{2}L}}\text{\hspace{1em}(1.16)}$$ \qquad 式(1.14)可化简为：$$\{\begin{array}{ll}x& =R_{N}cosL\\ z& =R_N(1-e^2)sinL\end{array}\text{\hspace{1em}(1.17)}$$ \qquad 对于同一点的地理纬度和地心纬度是存在偏差的，或者说地理垂线和地心垂线之间是存在偏差的。对于地心纬度有$\tan \phi =\frac{z}{x}$由式(1.13)可得：$$\tan \phi =(1-e^2)tanL\text{\hspace{1em}(1.18)}$$ 地理纬度与地心纬度间的偏差记作$\Delta L=L-\phi$，则有$$\begin{array}{rl}tan\Delta L& =tan(L-\phi )=\frac{tanL-tan\phi }{1+tanLtan\phi }\\ & =\frac{tanL-(1-e^2)tanL}{1+tanL(1-e^2)tanL}=\frac{e^{2}sinLcosL}{1-e^{2}si{n}^{2}L}\end{array}\text{\hspace{1em}(1.19)}$$ \qquad 若将$\Delta L$和$e$视为无穷小量，则式(1.19)可近似为:$$\Delta L\approx e^{2}sinLcosL(1+e^{2}si{n}^{2}L)=\frac{e^2(1+e^{2}si{n}^{2}L)}{2}sin2L\text{\hspace{1em}(1.20)}$$ \qquad 进一步有:$$\Delta L\approx \frac{e^2}{2}sin2L\approx fsin2L\text{\hspace{1em}(1.21)}$$ \qquad 若取椭球扁率$f=\frac{1}{298.257}$，则在纬度$L=45°$处的地理纬度与地心纬度的偏差$\Delta L=11.5^{\prime }$。

2.曲率半径

\qquad 在微分几何中，我们把曲率的倒数称为曲率半径，即$r=\frac{1}{k}$。曲率是指平面上一条曲线在特定点的弯曲程度。直线不弯曲，其各点的曲率均为0，圆上各点的曲率都一样，但小圆比大圆弯曲的厉害。
\qquad 一条曲线$AB$的平均曲率${\overline{K}}_{AB}$定义为动点沿AB从A移动到B时的切线转过的角度$|\Delta \alpha |$与AB的长度$|\Delta s|之比$，即$${\overline{K}}_{AB}=\frac{|\Delta \alpha |}{|\Delta s|}\text{\hspace{1em}(2.1)}$$ \qquad 曲线AB上A点处的曲率$K_A$定义为：$$K_A=\underset{P沿\overline{AB}\to A}{\lim }{\overline{K}}_{AB}=\underset{\Delta s\to 0}{\lim }|\frac{\Delta \alpha }{\Delta s}|\text{\hspace{1em}(2.2)}$$ \qquad 假设函数$y(x)$二阶可导，曲线上一点$(x,y(x))$处的曲率$K=K(x)$。

\qquad 如上图所示，曲线C的方程为$y=y(x)$，点$M(x,y(x))，M^{\prime }(x+\Delta x,y(x+\Delta x))$，点M处的切线倾角$\alpha =\alpha (x)$，则$tan(\alpha )=y(x{)}^{\prime }$。
\qquad 在曲线C上取一点$M_0(x_0,y(x_0))$作为度量弧段的基点，弧$M_{0}M$的长度为$s(x)$，则曲率定义为：$$\begin{array}{rl}K=K(x)& =\underset{M^{\prime }沿C\to M}{\lim }{\overline{K}}_{M{M}^{\prime }}=\underset{M^{\prime }沿C\to M}{\lim }|\frac{\alpha (x+\Delta x)-\alpha (x)}{s(x+\Delta x)-s(x)}|\\ & =\underset{\Delta x\to 0}{\lim }|\frac{[\alpha (x+\Delta x)-\alpha (x)]/\Delta (x)}{[s(x+\Delta x)-s(x)]/\Delta (x)}|\\ & =\frac{|{\alpha }^{\prime }(x)|}{|s^{\prime }(x)|}\end{array}\text{\hspace{1em}(2.3)}$$ \qquad 因为$$tan(x)=y^{\prime }(x)⟹\alpha (x)=arctan(y^{\prime }(x))⟹{\alpha }^{\prime }(x)=\frac{y^{\prime \prime }(x)}{1+y^{\prime }(x)}\text{\hspace{1em}(2.4)}$$ \qquad 将$s^{\prime }(x)=\sqrt{1+y^{\prime }(x{)}^2}$及式(2.4)带入式(2.3)得$y=y(x)$在$(x,y(x))$处的曲率为：$$K(x)=\frac{|{\alpha }^{\prime }(x)|}{|s^{\prime }(x)|}=\frac{|y^{\prime \prime }(x)|}{(1+y^{\prime }(x{)}^2{)}^{\frac{3}{2}}}\text{\hspace{1em}(2.5)}$$ \qquad 或 $$K=\frac{|y^{\prime \prime }|}{(1+y^{\prime 2}{)}^{\frac{3}{2}}}\text{\hspace{1em}(2.6)}$$ \qquad 由于地球近似为一个参考椭球，所以地球表面的曲率半径不是一个常数，而是随位置和方位的不同而发生变化。因此在惯导系统中进行经纬度解算时，必须考虑地球半径变化的情况。在求取地球表面某一点M的纬度时，首先根据北向速度$V_N$求出纬度的变化率$\dot{L}=\frac{V_N}{R_N}$，对其积分得到纬度$L$。
\qquad 根据曲率半径的定义，可得地球子午面内主曲率半径$R_N$为：$$R_N=\frac{1}{K}=\frac{[1+(\frac{d{z}_e}{d{x}_e}{)}^2{]}^{\frac{3}{2}}}{\frac{d^2{z}_e}{d^2{x}_e}}\text{\hspace{1em}(2.7)}$$ \qquad 其中$$\frac{d{z}_e}{d{x}_e}=-cotL\text{\hspace{1em}(2.8)}$$ \qquad 由(2.7)得$$\frac{d^2{z}_e}{d^2{x}_e}=\frac{1}{si{n}^{2}L}\frac{dL}{d{x}_e}\text{\hspace{1em}(2.9)}$$ \qquad 由(1.14)得$$\frac{d{x}_e}{dL}=\frac{(1-e{)}^2{R}_{e}tanLse{c}^{2}L}{[(1-e{)}^2{R}_{e}ta{n}^{2}L+1{]}^{\frac{3}{2}}}\text{\hspace{1em}(2.10)}$$ \qquad 将式(2.9)带入式(2.8)的：$$\frac{d^2{x}_e}{d^{2}L}=\frac{1}{si{n}^{2}L}\frac{-[(1-e{)}^2{R}_{e}ta{n}^{2}L+1{]}^{\frac{3}{2}}}{(1-e{)}^2{R}_{e}ta{n}^{2}L+1}=-\frac{[(1-e{)}^2{R}_{e}ta{n}^{2}L+1{]}^{\frac{3}{2}}}{(1-e{)}^2{R}_{e}ta{n}^{3}L}\text{\hspace{1em}(2.11)}$$ \qquad 展开$(1-e{)}^2$项并略去$e^2$项得:$$R_N\approx \frac{(1-2e)R_e}{(1-2esi{n}^{2}L{)}^{\frac{3}{2}}}\approx R_e(1-2e+3esi{n}^{2}L)\text{\hspace{1em}(2.12)}$$ \qquad 卯酉圈的主曲率半径$R_M$为:$$R_M=\frac{x_e}{cosL}=\frac{R_e}{[(1-e{)}^{2}si{n}^{2}L+co{s}^{2}L{]}^{\frac{1}{2}}}\text{\hspace{1em}(2.13)}$$ \qquad 展开$(1-e{)}^2$项并略去$e^2$项得:$$R_M\approx \frac{R_e}{(1-2esi{n}^{2}L{)}^{\frac{1}{2}}}\approx R_e(1+esi{n}^{2}L)\text{\hspace{1em}(2.14)}$$

参考文献
[1] 同济大学数学系. 高等数学.第7版[M]. 高等教育出版社, 2014.
[2] 严恭敏, 翁浚. 捷联惯导算法与组合导航原理[M]. 西安:西北工业大学出版社, 2020.
[3] 朱家海. 惯性导航[M]. 北京:国防工业出版社,2008.