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}$$

这个 $N$ 是一个有意义的量。我们后续会提到它。

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

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

$$\{\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}}$

这个 $M$ 是一个有意义的量。我们后续会提到它。

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

$${\mathbf{T}}_{\text{meridian}}=\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$ 处的平行圈单位切向量：

$${\mathbf{T}}_{\text{parallel}}\left[\begin{array}{c}-\sin L\\ \cos L\\ 0\end{array}\right]$$

其模长为 $1$。

椭球面法向量

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

$$\mathbf{n}=\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$。

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

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

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

椭球法截线的曲率及曲率半径

空间曲线的曲率及曲率半径

高数课程中介绍过曲率以及曲率半径。我们简单回顾一下。

曲线的曲率 $k$ 用于反映曲线的弯曲程度。一段弧 $\stackrel{⌢}{M{M}_0}$ 的平均曲率定义为 $M$ 点切向量与 $M^{\prime }$ 点切向量的夹角，与这段弧长之比。$M,M^{\prime }$ 两点无限接近时，定义这段曲线的曲率为 $M$ 处的曲率。定义曲率的倒数为曲率半径。曲率越大或曲率半径越小，则曲线的完全程度越高。

则曲线上一点 $C$ 的曲率记为该点切向量的一阶导数向量的模，或为该点位置向量 $\mathbf{r}(s)$ 二阶导数向量的模，即有

$$k(s){|}_{s=s_c}=\left|{\frac{\mathrm{d}\mathbf{T}(s)}{\mathrm{d}s}|}_{s=s_c}\right|=\left|{\frac{{\mathrm{d}}^2\mathbf{r}(s)}{\mathrm{d}{s}^2}|}_{s=s_c}\right|$$

一种几何的理解方式

首先，若曲线以弧长 $s$ 为参数，则 $|\mathrm{d}\mathbf{r}|=\mathrm{d}s$，$|\mathbf{T}|=\left|\frac{\mathrm{d}\mathbf{r}}{\mathrm{d}s}\right|=1$。

假设有一段弧 $\mathrm{d}s$，前后位置向量关于 $s$ 的导数分别为 $\mathbf{T},{\mathbf{T}}^{\prime }$。$\mathbf{T},{\mathbf{T}}^{\prime }$ 在几何上是切向量，与 $\mathrm{d}\mathbf{T}$ 组成三角形。由于 $\mathrm{d}\mathbf{T}/\mathbf{T}\to 0$，二者的夹角 $\mathrm{d}\alpha =\sin \mathrm{d}\alpha =\tan \mathrm{d}\alpha =|\mathrm{d}\mathbf{T}|/|\mathbf{T}|=|\mathrm{d}\mathbf{T}|$。

所以 $k=\frac{\mathrm{d}\alpha }{\mathrm{d}s}=\left|\frac{\mathrm{d}\mathbf{T}}{\mathrm{d}s}\right|$。

若曲线不以弧长 $s$ 为参数，则由 $\frac{\mathrm{d}\mathbf{r}(s)}{\mathrm{d}s}=\frac{\mathrm{d}\mathbf{r}(s)}{\mathrm{d}t}\frac{\mathrm{d}t}{\mathrm{d}s}$ 再推及二阶导数式，即可得到在任意参数 $t$ 下，曲率的计算公式为

$$k(t)=\frac{\left|\frac{\mathrm{d}\mathbf{r}}{\mathrm{d}t}\times \frac{{\mathrm{d}}^2\mathbf{r}}{\mathrm{d}{t}^2}\right|}{{\left|\frac{\mathrm{d}\mathbf{r}}{\mathrm{d}t}\right|}^3}$$

一种物理的理解方式

我们知道曲率半径的引入让我们可以将更一般的曲线运动当作圆周运动处理。回忆一下高中物理中介绍的圆周运动向心加速度公式：

$$a=\frac{v^2}{r}$$

把半径 $r$ 改成曲率半径 $\rho$，整理一下：

$$\begin{array}{c}|{\mathbf{a}}_{向心}|=\frac{|\mathbf{v}{|}^2}{\rho }\\ k=\frac{1}{\rho }=\frac{|{\mathbf{a}}_{向心}|}{|\mathbf{v}{|}^2}\end{array}$$

这里的 ${\mathbf{a}}_{向心}$ 是向心加速度，也就是正交于速度 $\mathbf{v}$ 方向上的加速度，也就是有：

$$|{\mathbf{a}}_{向心}|=|{\mathbf{a}}_{总}|\sin ⟨\mathbf{v},{\mathbf{a}}_{总}⟩=\frac{|\mathbf{v}\times {\mathbf{a}}_{总}|}{|\mathbf{v}|}$$

带回原式就有

$$k=\frac{|\mathbf{v}\times {\mathbf{a}}_{总}|}{|\mathbf{v}{|}^3}$$

取参数为时间 $t$，那么就有 $\mathbf{v}=\frac{\mathrm{d}\mathbf{r}}{\mathrm{d}t}$，${\mathbf{a}}_{总}=\frac{{\mathrm{d}}^2\mathbf{r}}{\mathrm{d}{t}^2}$，这就得到了：

$$k=\frac{\left|\frac{\mathrm{d}\mathbf{r}}{\mathrm{d}t}\times \frac{{\mathrm{d}}^2\mathbf{r}}{\mathrm{d}{t}^2}\right|}{{\left|\frac{\mathrm{d}\mathbf{r}}{\mathrm{d}t}\right|}^3}$$

法截线的定义

过椭球面上一点的椭球面法线所作的平面称为法截面，法截面同椭球的交线则称为法截线（也称法截弧）。上一节的天球部分中提到的垂直圈（地平经圈）的概念与此类似。

包含椭球面上一点的法线可以构成无数个法截面，相应的也有无数条法截线，它们都是平面曲线，而且曲率各不相同。

由于旋转曲面的特性，曲率的大小与经度无关。

子午圈曲率半径

子午线是特殊的法截线。

不失一般性，取 $XOZ$ 坐标面上的起始子午线，推求其上任一点 $M$ 的子午线曲率半径。为简化计算，使用归化纬度 $u$ 作为参数，则有 $M$ 的位置向量 $\mathbf{r}(u)$

$$\mathbf{r}(u)=\left[\begin{array}{c}X(u)\\ Y(u)\\ Z(u)\end{array}\right]=\left[\begin{array}{c}a\cos u\\ 0\\ b\sin u\end{array}\right]$$

求一阶和二阶导数

$$\begin{array}{rl}\frac{\mathrm{d}\mathbf{r}(u)}{\mathrm{d}u}& =\left[\begin{array}{c}-a\sin u\\ 0\\ b\cos u\end{array}\right]\\ \frac{{\mathrm{d}}^2\mathbf{r}(u)}{\mathrm{d}{u}^2}& =\left[\begin{array}{c}-a\cos u\\ 0\\ -b\sin u\end{array}\right]\end{array}$$

于是得到

$$\begin{array}{c}\frac{\mathrm{d}\mathbf{r}(u)}{\mathrm{d}u}\times \frac{{\mathrm{d}}^2\mathbf{r}(u)}{\mathrm{d}{u}^2}={\left[\begin{array}{ccc}0& -ab& 0\end{array}\right]}^{\mathrm{T}}\\ \left|\frac{\mathrm{d}\mathbf{r}(u)}{\mathrm{d}u}\times \frac{{\mathrm{d}}^2\mathbf{r}(u)}{\mathrm{d}{u}^2}\right|=ab\\ \left|\frac{\mathrm{d}\mathbf{r}(u)}{\mathrm{d}u}\right|=\sqrt{b^2{\cos }^{2}u+a^2{\sin }^{2}u}\end{array}$$

代入前面的曲率半径即可得到

$$k(u)=\frac{ab}{(b^2{\cos }^{2}u+a^2{\sin }^{2}u{)}^{3/2}}$$

利用我们前面提到过的 $u$ 与 $B$ 的关系，可以推出以大地纬度 $B$ 表达的曲率半径 $k(B)$：

$$k(B)=\frac{W^3}{a(1-e^2)}$$

推导过程

我们之前在推导 以大地纬度与大地经度为参数的椭球参数方程 时得到过如下式子：

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

设上面 $k(u)$ 分母中的表达式 $I=b^2{\cos }^{2}u+a^2{\sin }^{2}u$，有

$$\begin{array}{rl}I& =\frac{b^2{\cos }^{2}B}{W^2}+\frac{b^2{\sin }^{2}B}{W^2}\\ & =\frac{b^2}{W^2}\end{array}$$

带回 $k(u)$ 得到

$$\begin{array}{rl}k(B)& =\frac{ab}{I^{3/2}}=ab{\left(\frac{W^2}{b^2}\right)}^{3/2}=\frac{a{W}^3}{b^2}\\ & =\frac{W^3}{a(1-e^2)}\end{array}$$

通常使用字母 $M$ 代表子午线曲率半径：

$$M=\frac{1}{k(B)}=\frac{a(1-e^2)}{W^3}$$

这个 $M$ 就是我们前面在 子午圈切向量 那里遇到的那个 $M$。

卯酉圈曲率半径

卯酉圈在某处的曲率 $k_n$ 满足：

$$k_n=k_r\cos B$$

其中 $k_r$ 是平行圈在该处的曲率。此结论涉及微分几何中的麦尼厄定理，这里不证。

平行圈是圆。所以其曲率半径即为圆的半径 $r_B$。我们已经有了大地纬度与大地经度表达的椭球参数方程：

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

定 $B$，则 $L$ 为参数画出一个圆，半径为 $r_B=N\cos B$。所以有

$$\begin{array}{c}{k}_r=\frac{1}{r_B}=\frac{1}{N\cos B}\\ k_n=k_r\cos B=\frac{1}{N\cos B}\cos B=\frac{1}{N}\end{array}$$

所以我们此前提到的 $N$ 就是某点卯酉圈在该点的曲率半径。

任意方向法截线的曲率半径

根据微分几何中的欧拉公式，大地方位角 $A$ 方向上法截线的曲率与子午、卯酉曲率半径的关系为：

$$k_A=\frac{{\cos }^{2}A}{M}+\frac{{\sin }^{2}A}{N}$$

因此，该方向上的曲率半径为

$$R_A=\frac{1}{k_A}=\frac{MN}{N{\cos }^{2}A+M{\sin }^{2}A}$$

当 $A=0,\frac{\pi }{2},\pi ,\frac{3\pi }{2}$ 时，取得极值。

平均曲率半径

平均曲率定义为某点所有方向法截线曲率半径的平均值。

$$R=\frac{2}{\pi }{\int }_0^{\pi /2}\frac{MN}{N{\cos }^{2}A+M{\sin }^{2}A}\mathrm{d}A$$

三角换元，设 $t=\sqrt{\frac{M}{N}}\tan A$，则有 $\mathrm{d}t=\sqrt{\frac{M}{N}}\frac{1}{{\cos }^{2}A}\mathrm{d}A$

代入得到

$$\begin{array}{rl}R& =\frac{2}{\pi }\sqrt{MN}{\int }_0^{\infty }\frac{\mathrm{d}t}{1+t^2}\\ & =\frac{2}{\pi }\sqrt{MN}\arctan t{|}_0^{\infty }\\ & =\sqrt{MN}=\frac{a\sqrt{1-e^2}}{W}\end{array}$$

公式汇总

一般方程：

• 椭球方程：$\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]$

曲率半径：

• 子午圈曲率半径：$M=\frac{a(1-e^2)}{W^3}$
• 卯酉圈曲率半径：$N=\frac{a}{W}$
• 平均曲率半径：$R=\sqrt{MN}=\frac{a\sqrt{1-e^2}}{W^2}$

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