3 地球椭球

椭球的几何元素

无论是参考椭球抑或是平均地球椭球，从几何上来说都是旋转椭球，亦即由长半轴为 $a$、短半轴为 $b$ 的椭圆绕其短轴而成的几何形体。

以椭球的中心 $O$ 为原点，旋转轴为 $Z$ 轴，赤道面为 $XOY$ 平面，构建空间直角坐标系，则旋转椭球面可以描述为

$$\frac{X^2}{a^2}+\frac{Y^2}{a^2}+\frac{Z^2}{b^2}=1$$

由椭球的长短半轴 $a,b$ 即可确定椭球的形状和大小。通常仅采用一个长度元素 $a$ 以确定椭球的大小，而用扁率 $\alpha$（有时也用 $f$）或第一偏心率 $e$ 来反映椭球的形状（扁平程度）。

定义扁率 $\alpha$ 与第一偏心率 $e$：

$$\begin{array}{c}\alpha =\frac{a-b}{a}\\ e=\frac{\sqrt{a^2-b^2}}{a},1-e^2=\frac{b^2}{a^2}\end{array}$$

二者的关系为：

$$e^2=2\alpha -{\alpha }^2$$

通常采用数值接近于零的 $\alpha$ 或 $e$ 取代 $b$，以便于作级数展开等运算。扁率 $\alpha$ 的倒数 $1/\alpha$ 常称作反扁率。

近似看作圆球面时，$\alpha =e=0$。

有时也使用第二偏心率 $e^{\prime }=\frac{\sqrt{a^2-b^2}}{b}$，满足 $1+e^{\prime 2}=\frac{a^2}{b^2}$。

椭球上的点和线

经线和纬线的曲线方程

我们在高等数学中学过，空间曲线通常使用两个面的交线表达。

以起始子午面为 $XOZ$ 坐标面建系，则起始子午线（$0^{\circ }$ 经线）是 $XOZ$ 面与椭球的交线。其方程为：

$$\frac{X^2}{a^2}+\frac{Z^2}{b^2}=1,Y=0$$

设 $M_0(X_0,Y_0,Z_0)$ 为该起始经线上一点，其三维坐标满足上式。

设 $M_0$ 旋转一定角度后生成一个点 $M(X,Y,Z)$，则称这个角度为 $M$ 处的大地经度，通常用字母 $L$ 表示。大地经度定义为所在子午面与起始子午面间二面角的平面角，规定向东（$Z$ 轴的右手螺旋方向）为正、向西为负。

经线圈方程可以表达为所在子午面截椭圆：

$$\frac{X^2}{a^2}+\frac{Z^2}{b^2}=1,\frac{Y}{X}=\tan L$$

点 $M_0$ 在旋转过程中，产生了旋转椭球面上的一个纬圈（平行圈）。纬圈是一个圆，所在平面与 $XOY$ 坐标面（赤道面）平行。设该平面与 $Z$ 轴交于 $A$，纬圈半径 $r$ 满足：

$$r=AM=A{M}_0=X_0=\sqrt{X^2+Y^2}$$

纬线圈方程可以表达为所在平面截椭圆：

$$\frac{X^2}{a^2}+\frac{Z^2}{b^2}=1,Z=Z_0$$

椭球面法线与大地纬度

过点 $M$，作平行圈的切线 $MR$，作子午圈的切线 $MT$，则 $MR$ 与 $MT$ 构成在 $M$ 椭球面的切平面。过点 $M$ 作切平面的垂线 $MP$，称为 $M$ 处的椭球面法线。

椭球面法线交旋转轴于 $P$，交赤道面于 $Q$。定义椭球面法线与赤道面的夹角为 $M$ 点处的大地纬度，通常用 $B$ 表示。规定北半球为正，南半球为负。

在微分几何的曲线理论中，利用曲线上任一点处的三个互相正交的单位向量来构成一个三维直角坐标系，称之为 Frenet 活动标架。在椭球面上，这里的 $MP,MR,MT$ 三个方向可以定义此处的 Frenet 标架。若将原点由 $M$ 移至测站中心，即可构建站心地平坐标系。我们将在后续章节详细讨论。

归化纬度与地心纬度

值得一提的是，椭球面法线不是点与椭球中心的连线，只有南北两极的椭球面法线会过椭球中心。

在大地测量中除了用大地纬度 $B$ 外，还采用归化纬度 $u$ 和地心纬度 $\phi$ 来表示椭圆上点的位置。

作中心与椭球中心 $O$ 重合、半径等于椭球长半轴的球。过点 $M$ 作赤道面的垂线，交该球于点 $M^{\prime }$。连接 $M^{\prime }$ 与椭球中心 $O$。定义 $O{M}^{\prime }$ 与赤道面的夹角为 $M$ 点的归化纬度，通常用字母 $\phi$ 表示。

连接点 $M$ 与椭球中心 $O$。定义 $OM$ 与赤道面的夹角为 $M$ 点的地心纬度，通常用字母 $\phi$ 表示。

下面是一幅交互式示意图。$M$ 点可以拖动，也可以点击左下角播放暂停按钮控制自动移动。

不难证明，地心纬度 $\phi$、归化纬度 $u$ 和大地纬度 $B$ 存在如下关系：

$$\tan \phi =\frac{b}{a}\tan u=\frac{b^2}{a^2}\tan B$$

证明过程

为熟悉起见，取过点 $M$ 的子午面重新建立平面直角坐标系 $xOy$。设截出的椭圆方程为 $\frac{x^2}{a^2}+\frac{y^2}{b^2}=1$，截出的圆方程为 $x^2+y^2=a^2$，点 $M$ 坐标为 $(x_M,y_M)$。

最简单的是地心纬度：$\tan \phi =k_{OM}=\frac{y_M}{x_M}$。

接下来考虑归化纬度。过点 $M$ 垂直于赤道面的直线为 $x=x_M$，即 $x_{M^{\prime }}=x_M$。代入圆方程求 $y_{M^{\prime }}$：

$$\{\begin{array}{l}{x}_{M^{\prime }}^2+y_{M^{\prime }}^2=a^2\\ x_{M^{\prime }}=x_M\end{array}\Rightarrow y_{M^{\prime }}^2=a^2-x_M^2$$

由于 $M$ 坐标符合椭圆方程，因此有

$$\begin{array}{rl}\frac{x_M^2}{a^2}+\frac{y_M^2}{b^2}& =1\\ x_M^2+\frac{a^2{y}_M^2}{b^2}& =a^2\\ a^2-x_M^2& =\frac{a^2{y}_M^2}{b^2}\end{array}$$

代入 $y_{M^{\prime }}$ 表达式得到

$$\begin{array}{c}{y}_{M^{\prime }}^2=a^2-x_M^2=\frac{a^2{y}_M^2}{b^2}\\ \Rightarrow y_{M^{\prime }}=\frac{a{y}_M}{b}\end{array}$$

故有

$$\tan u=k_{O{M}^{\prime }}=\frac{y_{M^{\prime }}}{x_{M^{\prime }}}=\frac{a{y}_M}{b}\cdot \frac{1}{x_M}=\frac{a}{b}\cdot \frac{y_M}{x_M}$$

值得一提的是，这里 $y$ 方向上的比例关系实际上对应了椭圆与圆在 $y$ 方向 $b/a$ 的缩放比例关系。我们后面还会碰到这个关系。

考虑大地纬度。设椭圆过点 $M$ 的切线为 $l$。要先求 $l$ 的斜率 $k_l$。

在初等数学中，我们可以用 $k_l$ 表达出 $l$ 的方程 $y-y_M=k(x-x_M)$，与椭圆方程联立并代入消 $y$ 后得到一个关于 $x$ 的二次方程，令判别式 $\Delta =0$ 解方程得到 $k_l$。不过我们现在已经是大学生了，不妨使用高等数学中介绍过的 隐函数求导 直接求解。

已知 $x,y$ 的关系由 $F(x,y)=\frac{x^2}{a^2}+\frac{y^2}{b^2}-1=0$ 确定。由偏导数

$$\begin{array}{c}{F}_x=\frac{\partial F}{\partial x}=\frac{2x}{a^2}\\ F_y=\frac{\partial F}{\partial y}=\frac{2y}{b^2}\end{array}$$

则有

$$\frac{\mathrm{d}y}{\mathrm{d}x}=-\frac{F_x}{F_y}=-\frac{2x}{a^2}\cdot \frac{b^2}{2y}=-\frac{b^2}{a^2}\cdot \frac{x}{y}$$

故

$$k_l={\frac{\mathrm{d}y}{\mathrm{d}x}|}_{\begin{array}{c}x=x_M\\ y=y_M\end{array}}=-\frac{b^2}{a^2}\cdot \frac{x_M}{y_M}$$

则有

$$\tan B=k_{QM}=-\frac{1}{k_l}=\frac{a^2}{b^2}\cdot \frac{y_M}{x_M}$$

总结一下，有：

$$\begin{array}{rl}\tan \phi & =\frac{y_M}{x_M}\\ \tan u& =\frac{a}{b}\cdot \frac{y_M}{x_M}\\ \tan B& =\frac{a^2}{b^2}\cdot \frac{y_M}{x_M}\end{array}$$

因此有

$$\tan \phi =\frac{b}{a}\tan u=\frac{b^2}{a^2}\tan B$$

证毕。

椭球的参数方程

正则曲面上的点有两个自由度，因此曲面上点位的坐标 $(X,Y,Z)$ 可表示成两个参数的单值函数，即曲面的参数方程。而曲线上的点只有一个自由度，可表示为单参数的函数。

以归化纬度与大地经度为参数

我们前面提到了归化纬度 $u$。为何要定义一个以长半轴为半径的球，导出这样一个量？我们知道，在平面上，圆与椭圆具有仿射变换关系。椭圆可以看作是圆在短轴上缩放 $b/a$ 倍的结果。可以利用这一点书写参数方程。

取首子午面（$XOZ$ 坐标面）研究。以 $a$ 为半径的圆，以 $u$ 为参数，参数方程可以写作：

$$\{\begin{array}{l}X=a\cos u\\ Z=a\sin u\end{array}$$

在短轴（$Z$ 轴）上缩放 $b/a$ 倍，即把 $Z$ 坐标乘以 $b/a$，即可得到首子午圈的曲线参数方程：

$$\{\begin{array}{l}X=a\cos u\\ Z=b\sin u\end{array}$$

首子午圈绕 $Z$ 轴旋转形成椭球面。取旋转角（也就是大地经度）$L$ 为第二个参数，即可写出参考椭球的参数方程：

$$\{\begin{array}{l}X=a\cos u\cos L\\ Y=a\cos u\sin L\\ Z=b\sin u\end{array}$$

三角函数将旋转语义表达得很直观。

参数方程里，固定大地经度 $L$ 即可得到对应经圈（子午圈）关于 $u$ 的参数方程；固定归化纬度 $u$ 即可得到对应纬圈（平行圈）关于 $L$ 的参数方程。

以大地纬度与大地经度为参数

我们前面已经提到大地纬度与归化纬度之间的转化关系：

$$\begin{array}{c}\tan \phi =\frac{b}{a}\tan u=\frac{b^2}{a^2}\tan B\\ \Rightarrow \tan u=\frac{b}{a}\tan B\end{array}$$

因此，只需要按此式推出 $\phi$ 与 $B$ 正弦、余弦值的关系，代入即可得到以大地纬度与大地经度为参数的椭球面参数方程了。

推导过程

$$\begin{array}{rl}{\cos }^{2}u& =\frac{1}{1+{\tan }^{2}u}\\ & =\frac{1}{1+\frac{b^2}{a^2}{\tan }^{2}B}\\ & =\frac{{\cos }^{2}B}{{\cos }^{2}B+\frac{b^2}{a^2}{\sin }^{2}B}\end{array}$$

其中的分母：

$$\begin{array}{rl}\text{Denominator}& ={\cos }^{2}B+\frac{b^2}{a^2}{\sin }^{2}B\\ & =1-{\sin }^{2}B+\frac{b^2}{a^2}{\sin }^{2}B\\ & =1-\frac{a^2-b^2}{a^2}{\sin }^{2}B\\ & =1-e^2{\sin }^{2}B\end{array}$$

带回 ${\cos }^{2}u$ 得

$$\begin{array}{rl}{\cos }^{2}u& =\frac{{\cos }^{2}B}{1-e^2{\sin }^{2}B}\\ \Rightarrow \cos u& =\frac{\cos B}{\sqrt{1-e^2{\sin }^{2}B}}\end{array}$$

由几何关系可得，$u$ 与 $B$ 必为同象限角，余弦值正负性相同。故舍去负根。

又有

$$\begin{array}{rl}{\sin }^{2}u& =1-{\cos }^{2}u\\ & =1-\frac{{\cos }^{2}B}{1-e^2{\sin }^{2}B}\\ & =\frac{1-e^2{\sin }^{2}B-{\cos }^{2}B}{1-e^2{\sin }^{2}B}\end{array}$$

其中的分子：

$$\begin{array}{rl}\text{Numerator}& =1-e^2{\sin }^{2}B-{\cos }^{2}B\\ & ={\sin }^{2}B-e^2{\sin }^{2}B\\ & =(1-e^2){\sin }^{2}B\\ & =\frac{b^2}{a^2}{\sin }^{2}B\end{array}$$

带回 ${\sin }^{2}u$ 得

$$\begin{array}{rl}{\sin }^{2}u& =\frac{b^2}{a^2}\frac{{\sin }^{2}B}{1-e^2{\sin }^{2}B}\\ \Rightarrow \sin u& =\frac{b\sin B}{a\sqrt{1-e^2{\sin }^{2}B}}\end{array}$$

由几何关系可得，$u$ 与 $B$ 必为同象限角，正弦值正负性相同。故舍去负根。

综上，有

$$\{\begin{array}{l}\cos u=\frac{\cos B}{\sqrt{1-e^2{\sin }^{2}B}}\\ \sin u=\frac{b\sin B}{a\sqrt{1-e^2{\sin }^{2}B}}\end{array}$$

又已知参考椭球关于 $u,L$ 的参数方程

$$\{\begin{array}{l}X=a\cos u\cos L\\ Y=a\cos u\sin L\\ Z=b\sin u\end{array}$$

代入即可得到关于 $B,L$ 的参数方程。

原始方程：

$$\{\begin{array}{l}X=\frac{a\cos B\cos L}{\sqrt{1-e^2{\sin }^{2}B}}\\ Y=\frac{a\cos B\sin L}{\sqrt{1-e^2{\sin }^{2}B}}\\ Z=\frac{b^2\sin B}{a\sqrt{1-e^2{\sin }^{2}B}}\end{array}$$

采用辅助符号：

$$W=\sqrt{1-e^2{\sin }^{2}B},N=\frac{a}{W}$$

代入则有：

$$\{\begin{array}{l}X=N\cos B\cos L\\ Y=N\cos B\sin L\\ Z=N(1-e^2)\sin B\end{array}$$

以球心纬度与大地经度为参数

推导过程类似大地纬度，从略。

$$\{\begin{array}{l}X=a\cos \phi \cos L\sqrt{\frac{1-e^2}{1-e^2{\cos }^2\phi }}\\ Y=a\cos \phi \sin L\sqrt{\frac{1-e^2}{1-e^2{\cos }^2\phi }}\\ Z=a\sin \phi \sqrt{\frac{1-e^2}{1-e^2{\cos }^2\phi }}\end{array}$$

椭球的几何性质

旋转椭球面具有下述性质：

• 对称性：椭球关于中心点对称、关于赤道面对称、关于任意子午面对称；
• 有界性：$|X|\le a,|Y|\le a,|Z|\le b$。也就是说，椭球面被限定在六个平面 $X=\pm a,Y=\pm a,Z=\pm b$ 所围成的长方体内；
• 正则性：椭球面是一个连续、封闭的正则曲面，每个曲面点上都有唯一确定的非零法向量；
• 不可展性：旋转椭球面不可展开为平面，若强行拉平则会出现重叠和裂口。（柱面与锥面是可展曲面）

椭球的切向量与法向量

子午圈切向量

我们可以把椭球上的点用位置矢量 $\mathbf{r}$ 表示：

$$\mathbf{r}=\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]=\left[\begin{array}{c}N\cos B\cos L\\ N\cos B\sin L\\ N(1-e^2)\sin B\end{array}\right]$$

其中 $W=\sqrt{1-e^2{\sin }^{2}B}$，$N=\frac{a}{W}=\frac{a}{\sqrt{1-e^2{\sin }^{2}B}}$

这样，$\frac{\partial \mathbf{r}}{\partial B}$ 与 $\frac{\partial \mathbf{r}}{\partial L}$ 就分别是椭球面上 $\mathbf{r}$ 处的经线和纬线方向的切向量，也就是子午圈与平行圈的切向量。

$\frac{\partial \mathbf{r}}{\partial B}$ 推导过程

考虑 $X$：

$$\frac{\partial X}{\partial B}=\frac{\partial N}{\partial B}\cos B\cos L-N\sin B\cos L$$

其中：

$$\begin{array}{rl}\frac{\partial N}{\partial B}& =\frac{\partial }{\partial B}\frac{a}{\sqrt{1-e^2{\sin }^{2}B}}\\ & =-\frac{a}{2(1-e^2{\sin }^{2}B{)}^{3/2}}\cdot \frac{\partial }{\partial B}(1-e^2{\sin }^{2}B)\\ & =-\frac{a}{2(1-e^2{\sin }^{2}B{)}^{3/2}}\cdot 2e^2\sin B\cos B\\ & =\frac{a{e}^2\sin B\cos B}{(1-e^2{\sin }^{2}B{)}^{3/2}}\\ & =\frac{a{e}^2\sin B\cos B}{W^3}\end{array}$$

带回 $\frac{\partial X}{\partial B}$ 得到：

$$\begin{array}{rl}\frac{\partial X}{\partial B}& =\frac{a{e}^2\sin B\cos B}{W^3}\cos B\cos L-\frac{a}{W}\sin B\cos L\\ & =\frac{a\sin B\cos L}{W^3}(e^2{\cos }^{2}B-W^2)\\ & =\frac{a\sin B\cos L}{W^3}(e^2{\cos }^{2}B-1+e^2{\sin }^{2}B)\\ & =\frac{a\sin B\cos L}{W^3}(e^2-1)\\ & =-\frac{a(1-e^2)}{W^3}\sin B\cos L\end{array}$$

令 $M=\frac{a(1-e^2)}{W^3}=\frac{a(1-e^2)}{(1-e^2{\sin }^{2}B{)}^{3/2}}$，则有：

$$\frac{\partial X}{\partial B}=-M\sin B\cos L$$

同理有：

$$\frac{\partial Y}{\partial B}=-M\sin B\sin L$$

下面考虑 $Z$：

$$\begin{array}{rl}\frac{\partial Z}{\partial B}& =\frac{\partial }{\partial B}N(1-e^2)\sin B\\ & =(1-e^2)\left(\frac{\partial N}{\partial B}\sin B+N\cos B\right)\\ & =(1-e^2)\left(\frac{a{e}^2\sin B\cos B}{W^3}\sin B+\frac{a}{W}\cos B\right)\\ & =\frac{a(1-e^2)\cos B}{W^3}(e^2\sin B+W^2)\\ & =\frac{a(1-e^2)\cos B}{W^3}(\cancel{e^2\sin B}+1-\cancel{e^2\sin B})\\ & =\frac{a(1-e^2)}{W^3}\cos B\\ & =M\cos B\end{array}$$

综上，有

$$\begin{array}{rl}\frac{\partial X}{\partial B}& =-M\sin B\cos L\\ \frac{\partial Y}{\partial B}& =-M\sin B\sin L\\ \frac{\partial Z}{\partial B}& =M\cos B\end{array}$$

可以得到：

$$\frac{\partial \mathbf{r}}{\partial B}=\left[\begin{array}{c}\partial X/\partial B\\ \partial Y/\partial B\\ \partial Z/\partial B\end{array}\right]=M\left[\begin{array}{c}-\sin B\cos L\\ -\sin B\sin L\\ \cos B\end{array}\right]$$

其中 $M=\frac{a(1-e^2)}{W^3}=\frac{a(1-e^2)}{(1-e^2{\sin }^{2}B{)}^{3/2}}$

所以有 $B,L$ 处的子午圈单位切向量：

$$\left[\begin{array}{c}-\sin B\cos L\\ -\sin B\sin L\\ \cos B\end{array}\right]$$

其模长为 $1$。

平行圈切向量

$\frac{\partial \mathbf{r}}{\partial L}$ 推导过程

$$\begin{array}{rl}\frac{\partial X}{\partial L}& =\frac{\partial }{\partial L}N\cos B\cos L=-N\cos B\sin L\\ \frac{\partial Y}{\partial L}& =\frac{\partial }{\partial L}N\cos B\sin L=N\cos B\cos L\\ \frac{\partial Z}{\partial L}& =\frac{\partial }{\partial L}N(1-e^2)\sin B=0\end{array}$$

可以得到：

$$\frac{\partial \mathbf{r}}{\partial N}=\left[\begin{array}{c}\partial X/\partial N\\ \partial Y/\partial N\\ \partial Z/\partial N\end{array}\right]=N\cos B\left[\begin{array}{c}-\sin L\\ \cos L\\ 0\end{array}\right]$$

所以有 $B,L$ 处的平行圈单位切向量：

$$\left[\begin{array}{c}-\sin L\\ \cos L\\ 0\end{array}\right]$$

其模长为 $1$。

椭球面法向量

将子午圈与平行圈的单位切向量作叉积，即可得到单位法向量：

$$\left[\begin{array}{c}-\sin B\cos L\\ -\sin B\sin L\\ \cos B\end{array}\right]\times \left[\begin{array}{c}-\sin L\\ \cos L\\ 0\end{array}\right]=\left[\begin{array}{c}\cos B\cos L\\ \cos B\sin L\\ \sin B\end{array}\right]$$

其模长为 $1$。

忘记叉积怎么算了？高数笔记传送门

所以，法向量处处存在且不等于零。椭球面是正则的。

椭球面法向量的方向也就是铅垂线的方向。

公式汇总

一般方程：

• 椭球方程：$\frac{X^2}{a^2}+\frac{Y^2}{a^2}+\frac{Z^2}{b^2}=1$
• 子午圈方程：与 $Y=X\tan L$ 联立
• 平行圈方程：与 $Z=Z_0$ 联立

三个纬度关系：

$$\tan \phi =\frac{b}{a}\tan u=\frac{b^2}{a^2}\tan B$$

参数方程：

$$\begin{array}{rl}\boxed{u+L}& \{\begin{array}{l}X=a\cos u\cos L\\ Y=a\cos u\sin L\\ Z=b\sin u\end{array}\\ \boxed{B+L}& \{\begin{array}{l}X=N\cos B\cos L\\ Y=N\cos B\sin L\\ Z=N(1-e^2)\sin B\end{array}N=\frac{a}{\sqrt{1-e^2{\sin }^{2}B}}\\ \boxed{\phi +L}& \{\begin{array}{l}X=a\cos \phi \cos L\sqrt{\frac{1-e^2}{1-e^2{\cos }^2\phi }}\\ Y=a\cos \phi \sin L\sqrt{\frac{1-e^2}{1-e^2{\cos }^2\phi }}\\ Z=a\sin \phi \sqrt{\frac{1-e^2}{1-e^2{\cos }^2\phi }}\end{array}\end{array}$$

切向量与法向量

• 子午圈单位切向量 $\left[\begin{array}{c}-\sin B\cos L\\ -\sin B\sin L\\ \cos B\end{array}\right]$
• 平行圈单位切向量 $\left[\begin{array}{c}-\sin L\\ \cos L\\ 0\end{array}\right]$
• 单位法向量 $\left[\begin{array}{c}\cos B\cos L\\ \cos B\sin L\\ \sin B\end{array}\right]$