6 法截弧与大地线

一个很常见的问题：从点 $A$ 到点 $B$，哪条路径是最短的？例如规划飞机或轮船的航线。

我们都知道两点之间，直线段最短。不过那是在平面上，地球不是平面。

如果把地球处理为一个球，那么其表面上连接两点的最短路径是大圆的一条劣弧，也就是过 $A,B$ 和球心 $O$ 的平面与球相交所得到的路径。这也称为大圆航线。可是如果考虑扁率，将地球作为椭球来处理，要如何找出这条线呢？

法截弧

一种很自然的想法是，尝试改用铅垂线构面，截椭球得到弧。换句话说就是取法截线的一部分。

定义 $A$ 到 $B$ 的法截弧为：由椭球面上 $A$ 点的法线与 $B$ 点所确定的法截面，与椭球面相割得到的曲线。

这里有一个问题：用 $A$ 点法线构造且过 $B$ 的法截面，和 $B$ 点法线构造且过 $A$ 的法截面，是同一个吗？

推导过程

利用我们在 第 3 节 和 第 4 节 中推过的式子，有位置向量 $\mathbf{r}$ 和法向量 $\mathbf{n}$：

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

不妨取 $A$ 点位置向量为 ${\mathbf{r}}_1$，$B$ 为 ${\mathbf{r}}_2$。直线 $AB$ 的方向我们可以定义一个 $\mathbf{l}$ 表示：

$$\mathbf{l}=\frac{1}{N}({\mathbf{r}}_1-{\mathbf{r}}_2)=\left[\begin{array}{c}\cos B_1\cos L_1-\cos B_2\cos L_2\\ \cos B_1\sin L_1-\cos B_2\sin L_2\\ (1-e^2)(\sin B_1-\sin B_2)\end{array}\right]$$

如果两个法截面是同一个，意味着 ${\mathbf{n}}_1,{\mathbf{n}}_2,\mathbf{l}$ 共面，三个向量线性相关，有实数 $p,q$ 满足

$$\mathbf{l}=p{\mathbf{n}}_1+q{\mathbf{n}}_2$$

展开，待定系数求 $p,q$：

$$\begin{array}{rl}\cos B_1\cos L_1-\cos B_2\cos L_2& =p\cos B_1\cos L_1+q\cos B_2\cos L_2\\ \cos B_1\sin L_1-\cos B_2\sin L_2& =p\cos B_1\sin L_1+q\cos B_2\sin L_2\\ (1-e^2)\sin B_1-(1-e^2)\sin B_2& =p\sin B_1+q\sin B_2\end{array}$$

两个变量，三个方程。如果这些系数值各不相同，考虑前两个方程，那就应该有 $p=1,q=-1$，可是这样就不符合第三个方程了。这意味着对于一般情况，$p,q$ 是无解的，也就是说三个向量不在一个平面内，两个法截面不重合。

考虑一些特殊情况。如果 $A,B$ 都在赤道面上，$\sin B_1=\sin B_2=0$，第三个方程消失，此时有解。或者 $A,B$ 在同一条子午线上，$L_1=L_2$，前两个方程的经度项全部约去，退化为同一个方程，此时有解。但这明显不是普遍情况。

推导结果告诉我们：$B$ 点的法线，通常不在 $A$ 点这个法截面上。换言之，$A$ 点的法线和 $B$ 点的法线通常是异面直线，找不到一个平面既过 $A$ 点的法线、又过 $B$ 点的法线。

如此一来，$A$ 和 $B$ 之间便可以画出两条法截弧：其中一条过 $A$ 点的法线和 $B$ 点，另一条过 $A$ 点和 $B$ 点的法线。根据上面的定义，它们分别称为「$A$ 到 $B$ 的法截弧」「$B$ 到 $A$ 的法截弧」。这样的两条法截弧称为相对法截弧。

下面是一幅交互式图形，展示了椭球上 $A,B$ 两点各自法向量、法截面，以及截出的两条相对法截弧。可以拖动 $A,B$ 两点，观察法截弧形态的变化。

这样来看，法截弧会是最短路径吗？不能相信直觉。我们接下来从数学手段推导。

变分法基础

WARNING

变分法相关内容教材没有要求，甚至没有提及，直接给了结果公式。也就是说这一部分是拓展内容，可以选择性阅读。

下面的推导过程可能不完全严谨，仅作为变分法入门介绍。不过对我们够用了。

欧拉 - 拉格朗日方程

推导欧拉 - 拉格朗日方程

我们先考虑更简单的情况：平面上两点之间的最短路径。设曲线为 $y=y(x)$，弧长泛函为：

$$J[y]={\int }_a^b\sqrt{1+y^{\prime 2}}\mathrm{d}x$$

这个 $J[y]$ 可以理解为关于函数的函数，参数是函数 $y=y(x)$，结果是一个数。我们的目标是找到最合适的 $y=y(x)$，使得 $S$ 最小。

更抽象地，我们可以表达为：

$$J[y]={\int }_a^{b}F(x,y,y^{\prime })\mathrm{d}x$$

这里 $F$ 是一个已知函数，依赖三个变量：自变量 $x$、要寻找的函数 $y=y(x)$、函数的导数 $y^{\prime }=\mathrm{d}y/\mathrm{d}x$。

现在给曲线加一个很小的扰动：

$$y_{\epsilon }(x)=y(x)+\epsilon \eta (x)$$

其中 $\epsilon$ 是一个很小的实数，$\eta (x)$ 是我们任意选取的光滑函数。因为两端点固定，所以要求 $\eta (a)=\eta (b)=0$。

将新曲线代入积分，得到一个关于普通变量 $\epsilon$ 的函数：

$$\Phi (\epsilon )={\int }_a^{b}F(x,y+\epsilon \eta ,y^{\prime }+\epsilon {\eta }^{\prime })\mathrm{d}x$$

如果原曲线是最短路径，那么至少在它附近，弧长对这种微小扰动的一阶变化应该为零。也就是说有

$${\frac{\partial \Phi }{\partial \epsilon }|}_{\epsilon =0}=0$$

那我们来求这个导数。

$$\begin{array}{rl}{\frac{\partial \Phi }{\partial \epsilon }|}_{\epsilon =0}& ={\int }_a^b\left(\frac{\partial F}{\partial y}\frac{\partial y}{\partial \epsilon }+\frac{\partial F}{\partial y^{\prime }}\frac{\partial y^{\prime }}{\partial \epsilon }\right)\mathrm{d}x\\ & ={\int }_a^b\left(\frac{\partial F}{\partial y}\eta +\frac{\partial F}{\partial y^{\prime }}{\eta }^{\prime }\right)\mathrm{d}x\end{array}$$

第二项的 ${\eta }^{\prime }$ 很麻烦，考虑把它化成只含 $\eta$ 的形式。回忆 分部积分法：

$${\int }_a^{b}u{v}^{\prime }\mathrm{d}x=[uv{]}_a^b-{\int }_a^b{u}^{\prime }v\mathrm{d}x$$

令 $u=\frac{\partial F}{\partial y^{\prime }}$，$v=\eta$，于是

$${\int }_a^b\frac{\partial F}{\partial y^{\prime }}{\eta }^{\prime }\mathrm{d}x={\left[\frac{\partial F}{\partial y^{\prime }}\eta \right]}_a^b-{\int }_a^b\frac{\mathrm{d}}{\mathrm{d}x}\left(\frac{\partial F}{\partial y^{\prime }}\right)\eta \mathrm{d}x$$

因为端点固定，$\eta (a)=\eta (b)=0$，所以第一项就消失了。带回去得到

$${\frac{\partial \Phi }{\partial \epsilon }|}_{\epsilon =0}={\int }_a^b\left[\frac{\partial F}{\partial y}-\frac{\mathrm{d}}{\mathrm{d}x}\left(\frac{\partial F}{\partial y^{\prime }}\right)\right]\eta \mathrm{d}x=0$$

由于 $\eta (x)$ 是任意选取的，整个式子要想恒等于零，则方括号里的部分必须恒为零：

$$\boxed{\frac{\partial F}{\partial y}-\frac{\mathrm{d}}{\mathrm{d}x}\left(\frac{\partial F}{\partial y^{\prime }}\right)=0}$$

这称为 欧拉 - 拉格朗日方程。

我们尝试把 $F=\sqrt{1+y^{\prime 2}}$ 套进来，就有：

$$\frac{\partial F}{\partial y}=0,\frac{\partial F}{\partial y^{\prime }}=\frac{y^{\prime }}{\sqrt{1+y^{\prime 2}}}$$

于是

$$\frac{\mathrm{d}}{\mathrm{d}x}\left(\frac{y^{\prime }}{\sqrt{1+y^{\prime 2}}}\right)=0$$

也就是说 $\frac{y^{\prime }}{\sqrt{1+y^{\prime 2}}}$ 是一个常数，因此 $y^{\prime }$ 是常数，令其为 $k$，则 $y(x)=kx+C$，其中 $C$ 是常数。两点之间确实是直线段最短，太神奇了。

总结一下就是，曲面上的最短路径曲线，满足欧拉 - 拉格朗日方程：

$$\boxed{\frac{\partial F}{\partial y}-\frac{\mathrm{d}}{\mathrm{d}x}\left(\frac{\partial F}{\partial y^{\prime }}\right)=0}$$

我们可以用它推导参考椭球上的最短路径。

推导参考椭球上的最短路径

推导参考椭球上的最短路径

用上一节的 第一基本形式：

$$\mathrm{d}S=\sqrt{M^2\mathrm{d}{B}^2+N^2{\cos }^{2}B\cdot \mathrm{d}{L}^2}$$

用参数 $t$ 表示曲线 $B(t),L(t)$，设其关于 $t$ 的导数分别为 $B^{\prime },L^{\prime }$，就有

$$S=\int \sqrt{M^2{B}^{\prime 2}+N^2{\cos }^2{L}^{\prime 2}}\mathrm{d}t$$

这就是椭球面上的弧长泛函。

直接对这个平方根形式应用变分法当然可以，但计算非常繁琐。这里有一个微分几何中常用的技巧：改为研究能量积分

$$E=\frac{1}{2}\int (M^2{B}^{\prime 2}+N^2{\cos }^2{L}^{\prime 2})\mathrm{d}t$$

对于固定端点、采用适当参数化的曲线，能量的极值曲线与弧长的极值曲线对应；特别是取恒定速率参数时，两者给出相同的测地线轨迹。

这个结论我们这里不证。

不妨设 $G=M^2{B}^{\prime 2}+N^2{\cos }^2{L}^{\prime 2}$。

先考虑纬度方向，使用欧拉 - 拉格朗日方程：

$$\boxed{\frac{\mathrm{d}}{\mathrm{d}t}\frac{\partial G}{\partial B^{\prime }}-\frac{\partial G}{\partial B}=0}$$

其中：

$$\frac{\partial G}{\partial B}=\frac{\partial M^2}{\partial B}{B}^{\prime 2}+\frac{\partial (N^2{\cos }^{2}B)}{\partial B}{L}^{\prime 2}$$

又有：

$$\begin{array}{c}\frac{\partial G}{\partial B^{\prime }}=2M^2{B}^{\prime }\\ \frac{\mathrm{d}}{\mathrm{d}t}(2M^2{B}^{\prime })=2\frac{\mathrm{d}{M}^2}{\mathrm{d}t}{B}^{\prime }+2M^2{B}^{\prime \prime }=2\frac{\partial M^2}{\partial B}{B}^{\prime 2}+2M^2{B}^{\prime \prime }\end{array}$$

代入就是：

$$\begin{array}{rl}2M^2{B}^{\prime \prime }+2\frac{\partial M^2}{\partial B}{B}^{\prime 2}-\frac{\partial M^2}{\partial B}{B}^{\prime 2}-\frac{\partial (N^2{\cos }^{2}B)}{\partial B}{L}^{\prime 2}& =0\\ 2M^2{B}^{\prime \prime }+\frac{\partial M^2}{\partial B}{B}^{\prime 2}-\frac{\partial (N^2{\cos }^{2}B)}{\partial B}{L}^{\prime 2}& =0\\ B^{\prime \prime }+\frac{\partial M^2/\partial B}{2M^2}{B}^{\prime 2}-\frac{\partial (N^2{\cos }^{2}B)/\partial B}{2M^2}{L}^{\prime 2}& =0\end{array}$$

以弧长 $S$ 为参数，即取 $t=S$，整理式子：

$$\boxed{\frac{{\mathrm{d}}^{2}B}{\mathrm{d}{S}^2}+\frac{\partial M^2/\partial B}{2M^2}{\left(\frac{\mathrm{d}B}{\mathrm{d}S}\right)}^2-\frac{\partial (N^2{\cos }^{2}B)/\partial B}{2M^2}{\left(\frac{\mathrm{d}L}{\mathrm{d}S}\right)}^2=0}$$

计算两个系数：

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

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

代入得到

$$\boxed{\frac{{\mathrm{d}}^{2}B}{\mathrm{d}{S}^2}+\frac{3e^2\sin B\cos B}{W^2}{\left(\frac{\mathrm{d}B}{\mathrm{d}S}\right)}^2+\frac{W^2\sin B\cos B}{1-e^2}{\left(\frac{\mathrm{d}L}{\mathrm{d}S}\right)}^2=0}$$

再考虑经度方向上，使用欧拉 - 拉格朗日方程：

$$\boxed{\frac{\mathrm{d}}{\mathrm{d}t}\frac{\partial G}{\partial L^{\prime }}-\frac{\partial G}{\partial L}=0}$$

$G$ 不显含 $L$，所以就是

$$\frac{\mathrm{d}}{\mathrm{d}t}{N}^2{\cos }^{2}B{L}^{\prime }=0$$

同样取弧长为参数，整理一下：

$$\begin{array}{rl}I& =\frac{\mathrm{d}}{\mathrm{d}S}\left(N^2{\cos }^{2}B\frac{\mathrm{d}L}{\mathrm{d}S}\right)\\ & =\frac{\mathrm{d}(N^2{\cos }^{2}B)}{\mathrm{d}S}\frac{\mathrm{d}L}{\mathrm{d}S}+N^2{\cos }^{2}B\frac{{\mathrm{d}}^{2}L}{\mathrm{d}{S}^2}\\ & =\frac{\mathrm{d}(N^2{\cos }^{2}B)}{\mathrm{d}B}\frac{\mathrm{d}B}{\mathrm{d}S}\frac{\mathrm{d}L}{\mathrm{d}S}+N^2{\cos }^{2}B\frac{{\mathrm{d}}^{2}L}{\mathrm{d}{S}^2}\\ & =0\end{array}$$

其中

$$\begin{array}{rl}\frac{\mathrm{d}(N^2{\cos }^{2}B)}{\mathrm{d}B}& =\frac{\mathrm{d}}{\mathrm{d}B}\frac{a^2{\cos }^{2}B}{1-e^2{\sin }^{2}B}\\ & =-\frac{2a^2\cos B\sin B}{1-e^2{\sin }^{2}B}+\frac{2e^2{a}^2{\cos }^{3}B\sin B}{(1-e^2{\sin }^{2}B{)}^2}\\ & =2a^2\cos B\sin B\cdot \frac{e^2{\cos }^{2}B-(1-e^2{\sin }^{2}B)}{(1-e^2{\sin }^{2}B{)}^2}\\ & =-\frac{(1-e^2)2a^2\cos B\sin B}{(1-e^2{\sin }^{2}B{)}^2}\end{array}$$

带回得到

$$-\frac{(1-e^2)2a^2\cos B\sin B}{(1-e^2{\sin }^{2}B{)}^2}\frac{\mathrm{d}B}{\mathrm{d}S}\frac{\mathrm{d}L}{\mathrm{d}S}+N^2{\cos }^{2}B\frac{{\mathrm{d}}^{2}L}{\mathrm{d}{S}^2}=0$$

整理一下就是

$$\boxed{\frac{{\mathrm{d}}^{2}L}{\mathrm{d}{S}^2}-\frac{2(1-e^2)\tan B}{W^2}\frac{\mathrm{d}B}{\mathrm{d}S}\frac{\mathrm{d}L}{\mathrm{d}S}=0}$$

最终我们得到了这样的最短路径表达式：

$$\boxed{\begin{array}{c}\frac{{\mathrm{d}}^{2}B}{\mathrm{d}{S}^2}+\frac{3e^2\sin B\cos B}{W^2}{\left(\frac{\mathrm{d}B}{\mathrm{d}S}\right)}^2+\frac{W^2\sin B\cos B}{1-e^2}{\left(\frac{\mathrm{d}L}{\mathrm{d}S}\right)}^2=0\\ \frac{{\mathrm{d}}^{2}L}{\mathrm{d}{S}^2}-\frac{2(1-e^2)\tan B}{W^2}\frac{\mathrm{d}B}{\mathrm{d}S}\frac{\mathrm{d}L}{\mathrm{d}S}=0\end{array}}$$

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

NOTE

教材上使用了第二偏心率 $e^{\prime }$ 来表达，并使用了一些中间符号，形式上更好看一点：

$$\boxed{\begin{array}{c}\frac{{\mathrm{d}}^{2}B}{\mathrm{d}{S}^2}+\frac{3t{\eta }^2}{V^2}{\left(\frac{\mathrm{d}B}{\mathrm{d}S}\right)}^2+V^{2}t{\cos }^{2}B{\left(\frac{\mathrm{d}L}{\mathrm{d}S}\right)}^2=0\\ \frac{{\mathrm{d}}^{2}L}{\mathrm{d}{S}^2}-\frac{2t}{V^2}\frac{\mathrm{d}B}{\mathrm{d}S}\frac{\mathrm{d}L}{\mathrm{d}S}=0\end{array}}$$

其中 $V=\sqrt{1+e^{\prime 2}{\cos }^{2}B}$，$\eta =e^{\prime }\cos B$，$t=\tan B$。

两种表达方式是等价的。利用 $e$ 和 $e^{\prime }$ 的定义式进行代换可以比较容易地从其中一种推出另一种。这里就不重复推了。

大地线

我们刚才推出的这条最短曲线称为大地线，也称测地线。其英文 geodesic 正来自大地测量学 geodesy。

我们来研究它的性质。

克莱劳定理

将参考椭球的第一基本形式：

$$\mathrm{d}{S}^2=M^2\mathrm{d}{B}^2+N^2{\cos }^{2}B\mathrm{d}{L}^2$$

整理为使用弧长 $S$ 作为参数，有：

$$M^2{\left(\frac{\mathrm{d}B}{\mathrm{d}S}\right)}^2+N^2{\cos }^{2}B{\left(\frac{\mathrm{d}L}{\mathrm{d}S}\right)}^2=1$$

这个 $p^2+q^2=1$ 的形式，不由得让人想起三角函数。而 $\frac{\mathrm{d}B}{\mathrm{d}S}$ 表达北向分量、$\frac{\mathrm{d}L}{\mathrm{d}S}$ 表达东向分量，在实际意义上也十分契合。

我们引入大地方位角 $A$，约定从北方向起算，向东为正，和我们在测量学中使用的方位角定义一致。那么沿大地线前进时，弧长微元在北向和东向的分量分别为：

$$\begin{array}{rl}\cos A& =M\frac{\mathrm{d}B}{\mathrm{d}S}\\ \sin A& =N\cos B\frac{\mathrm{d}L}{\mathrm{d}S}\end{array}$$

在推导大地线方程时，经度方向上应用欧拉 - 拉格朗日公式时有这样一个中间结论：

$$\frac{\mathrm{d}}{\mathrm{d}t}\frac{\partial G}{\partial L^{\prime }}-\frac{\partial G}{\partial L}=0\Rightarrow \frac{\mathrm{d}}{\mathrm{d}S}\left(N^2{\cos }^{2}B\frac{\mathrm{d}L}{\mathrm{d}S}\right)=0$$

也就是说被微分的部分是常数，设为 $C$：

$$N^2{\cos }^{2}B\frac{\mathrm{d}L}{\mathrm{d}S}=N\cos B\left(N\cos B\frac{\mathrm{d}L}{\mathrm{d}S}\right)=C$$

代入方位角关系

$$N\cos B\frac{\mathrm{d}L}{\mathrm{d}S}=\sin A$$

得到了

$$N\cos B\sin A=r_B\sin A=C$$

这称为克莱劳定理或者克莱劳关系，Clairaut 有时也译作克莱罗。其意义是：大地线上各点的平行圈半径，与该点的大地线方位角的正弦之积为常数。同一条大地线上，大地线方位角与平行圈半径的变动相互制约。$C\ne 0$ 即不与子午圈重合时，越向两极移动，方向与赤道夹角越小。

大地线的曲率与测地曲率

前面我们推得的大地线方程，其形式很「隐」，是通过 $B,L$ 对 $S$ 的微分来表达的。我们需要一些更实在的东西，比如位置向量 $\mathbf{r}$ 和单位切向量 $\mathbf{T}=\frac{\mathrm{d}\mathbf{r}}{\mathrm{d}S}$。这里求导能得到直接模为 $1$ 的单位切向量，是由于取弧长 $S$ 作为参数，所以 $\left|\frac{\mathrm{d}\mathbf{r}}{\mathrm{d}S}\right|=1$。

我们此前在求子午圈和卯酉圈曲率的时候，利用了曲率向量 $\mathbf{K}=\frac{\mathrm{d}\mathbf{T}}{\mathrm{d}S}$，曲率 $k=|\mathbf{K}|$。类似地，我们来求大地线的曲率向量。

推导 $\mathbf{K}$

有单位切向量 $\mathbf{T}$ 满足

$$\mathbf{T}=\frac{\mathrm{d}\mathbf{r}}{\mathrm{d}S}=\frac{\mathrm{d}\mathbf{r}}{\mathrm{d}B}\frac{\mathrm{d}B}{\mathrm{d}S}+\frac{\mathrm{d}\mathbf{r}}{\mathrm{d}L}\frac{\mathrm{d}L}{\mathrm{d}S}$$

大地线也在椭球面上，因此满足椭球面的性质，我们 第 4 节 已经推过了

$$\begin{array}{rlrl}\frac{\partial \mathbf{r}}{\partial B}& =M{\mathbf{e}}_N,& {\mathbf{e}}_N& =\left[\begin{array}{c}-\sin B\cos L\\ -\sin B\sin L\\ \cos B\end{array}\right]\\ \frac{\partial \mathbf{r}}{\partial L}& =N\cos B{\mathbf{e}}_E,& {\mathbf{e}}_E& =\left[\begin{array}{c}-\sin L\\ \cos L\\ 0\end{array}\right]\end{array}$$

再利用前面定义的大地方位角 $A$：

$$\cos A=M\frac{\mathrm{d}B}{\mathrm{d}S}=M{B}^{\prime },\sin A=N\cos B\frac{\mathrm{d}L}{\mathrm{d}S}=N{L}^{\prime }\cos B$$

代入就可以把 $\mathbf{T}$ 写成经纬两个方向分量：

$$\mathbf{T}={\mathbf{e}}_N\cos A+{\mathbf{e}}_E\sin A$$

先对两个单位向量求导：

$$\begin{array}{rl}\frac{\partial {\mathbf{e}}_N}{\partial B}& =\left[\begin{array}{c}-\cos B\cos L\\ -\cos B\sin L\\ -\sin B\end{array}\right]=-\mathbf{n}\\ \frac{\partial {\mathbf{e}}_N}{\partial L}& =\left[\begin{array}{c}\sin B\sin L\\ -\sin B\cos L\\ 0\end{array}\right]=-\sin B{\mathbf{e}}_E\\ \frac{\partial {\mathbf{e}}_E}{\partial B}& =0\\ \frac{\partial {\mathbf{e}}_E}{\partial L}& =\left[\begin{array}{c}-\cos L\\ -\sin L\\ 0\end{array}\right]=\sin B{\mathbf{e}}_N-\cos B\mathbf{n}\end{array}$$

其中 $\mathbf{n}$ 是我们之前推过的椭球面单位法向量 $\mathbf{n}=\left[\begin{array}{c}\cos B\cos L\\ \cos B\sin L\\ \sin B\end{array}\right]$。

接下来还是使用 $B^{\prime },L^{\prime }$ 表达 $B,L$ 对 参数 $S$ 的导数。利用链式法则：

$$\begin{array}{rl}\frac{\mathrm{d}{\mathbf{e}}_N}{\mathrm{d}S}& =-B^{\prime }\mathbf{n}-L^{\prime }\sin B{\mathbf{e}}_E\\ \frac{\mathrm{d}{\mathbf{e}}_E}{\mathrm{d}S}& =L^{\prime }\sin B{\mathbf{e}}_N-L^{\prime }\cos B\mathbf{n}\end{array}$$

所以

$$\begin{array}{rl}\mathbf{K}& =\frac{\mathrm{d}}{\mathrm{d}S}({\mathbf{e}}_N\cos A+{\mathbf{e}}_E\sin A)\\ & =-A^{\prime }\sin A{\mathbf{e}}_N+\cos A\frac{\mathrm{d}{\mathbf{e}}_N}{\mathrm{d}S}\\ & +A^{\prime }\cos A{\mathbf{e}}_E+\sin A\frac{\mathrm{d}{\mathbf{e}}_E}{\mathrm{d}S}\end{array}$$

再代入两个基底向量的导数，得到

$$\boxed{\begin{array}{rl}\mathbf{K}& =(-A^{\prime }\sin A+L^{\prime }\sin B\sin A){\mathbf{e}}_N\\ & +(A^{\prime }\cos A-L^{\prime }\sin B\cos A){\mathbf{e}}_E\\ & -(B^{\prime }\cos A+L^{\prime }\cos B\sin A)\mathbf{n}\end{array}}$$

如此我们得到了曲率向量 $\mathbf{K}$ 在北、东、法线三个基底方向上的分量表达式，其中的 $A^{\prime }$ 是 $A$ 对参数 $S$ 的导数。这个 $A^{\prime }$ 比较麻烦，想办法消掉它。不妨设 $s=\sin A$，$c=\cos A$。于是

$$\begin{array}{c}{s}^{\prime }=A^{\prime }\cos A,c^{\prime }=-A^{\prime }\sin A\\ \boxed{c{s}^{\prime }-s{c}^{\prime }=A^{\prime }({\cos }^{2}A+{\sin }^{2}A)=A^{\prime }}\end{array}$$

所以我们只需要分别求 $c^{\prime }$ 和 $s^{\prime }$，不需要对 $A^{\prime }$ 求导。

先求 $s’$。

$$\begin{array}{rl}{s}^{\prime }& =(L^{\prime }\cdot N\cos B{)}^{\prime }=L^{\prime \prime }(N\cos B)+L^{\prime }(N\cos B{)}^{\prime \prime }\\ & =N\cos B\left(L^{\prime \prime }+B^{\prime }{L}^{\prime }\frac{\mathrm{d}(N\cos B)/\mathrm{d}B}{N\cos B}\right)\\ & =N\cos B\left(L^{\prime \prime }+B^{\prime }{L}^{\prime }\frac{\mathrm{d}}{\mathrm{d}B}\ln (N\cos B)\right)\end{array}$$

其中

$$\begin{array}{rl}\frac{\mathrm{d}}{\mathrm{d}B}\ln (N\cos B)& =\frac{\mathrm{d}}{\mathrm{d}B}\ln \frac{a\cos B}{\sqrt{1-e^2{\sin }^{2}B}}\\ & =\frac{\mathrm{d}}{\mathrm{d}B}\left[\ln a+\ln \cos B-\frac{1}{2}\ln (1-e^2{\sin }^{2}B)\right]\\ & =0-\frac{1}{\cos B}\sin B+\frac{1}{2}\frac{1}{1-e^2{\sin }^{2}B}\cdot 2e^2\sin B\cos B\\ & =-\tan B+\frac{e^2\sin B\cos B}{1-e^2{\sin }^{2}B}\\ & =-\tan B+\frac{e^2{\cos }^{2}B\tan B}{1-e^2{\sin }^{2}B}\\ & =\frac{(e^2{\cos }^{2}B-1+e^2{\sin }^{2}B)\tan B}{1-e^2{\sin }^{2}B}\\ & =-\frac{(1-e^2)\tan B}{W^2}\end{array}$$

带回

$$s^{\prime }=N\cos B\left[L^{\prime \prime }-\frac{(1-e^2)\tan B}{W^2}{B}^{\prime }{L}^{\prime }\right]$$

根据第二个大地线方程

$$L^{\prime \prime }-\frac{2(1-e^2)\tan B}{W^2}{B}^{\prime }{L}^{\prime }=0$$

所以

$$s^{\prime }=N\cos B\frac{(1-e^2)\tan B}{W^2}{B}^{\prime }{L}^{\prime }=\boxed{\frac{(1-e^2)\tan B}{W^2}{B}^{\prime }s}$$

接下来再求 $c^{\prime }=M{B}^{\prime }$。

$$c^{\prime }=\frac{\mathrm{d}M}{\mathrm{d}B}{B}^{\prime 2}+M{B}^{\prime \prime }$$

其中

$$\begin{array}{rl}\frac{\mathrm{d}M}{\mathrm{d}B}& =\frac{\mathrm{d}}{\mathrm{d}B}\frac{a(1-e^2)}{(1-e^2{\sin }^{2}B{)}^{3/2}}\\ & =\frac{3}{2}a(1-e^2)(1-e^2{\sin }^{2}B{)}^{-5/2}\cdot 2e^2\sin B\cos B\\ & =\frac{3e^2\sin B\cos B}{W^2}M\end{array}$$

所以

$$\boxed{c^{\prime }=-\frac{M{W}^2\sin B\cos B}{1-e^2}{L}^{\prime 2}}$$

接下来利用 $A^{\prime }=c{s}^{\prime }-s{c}^{\prime }$，有

$$A^{\prime }=sM\left[\frac{(1-e^2)\tan B}{W^2}{B}^{\prime 2}+\frac{W^2\sin B\cos B}{1-e^2}{L}^{\prime 2}\right]$$

现在代入

$$B^{\prime }=\frac{\cos A}{M},L^{\prime }=\frac{\sin A}{N\cos B}$$

于是有

$$\begin{array}{rl}{A}^{\prime }& =sM\left[\frac{(1-e^2)\tan B{\cos }^{2}A}{M^2{W}^2}+\frac{W^2\sin B\cos B{\sin }^{2}A}{(1-e^2)N^2{\cos }^{2}B}\right]\\ & =sM\left[\frac{(1-e^2)\tan B{\cos }^{2}A}{M^2{W}^2}+\frac{W^2\tan B{\sin }^{2}A}{(1-e^2)N^2}\right]\\ & =sM\left[\frac{W^4\tan B{\cos }^{2}A}{a^2(1-e^2)}+\frac{W^4\tan B{\sin }^{2}A}{a^2(1-e^2)}\right]\\ & =sM\frac{W^4\tan B}{a^2(1-e^2)}({\cos }^{2}A+{\sin }^{2}A)\\ & =s\frac{a(1-e^2)}{W^3}\frac{W^4\tan B}{a^2(1-e^2)}\\ & =s\frac{W\tan B}{a}\\ & =\frac{s\tan B}{N}\end{array}$$

再代入 $s=\sin A=N{L}^{\prime }\cos B$，得到

$$\boxed{A^{\prime }=L^{\prime }\sin B}$$

回到前面的曲率向量 $\mathbf{K}$：

$$\begin{array}{rl}\mathbf{K}& =(-A^{\prime }\sin A+L^{\prime }\sin B\sin A){\mathbf{e}}_N\\ & +(A^{\prime }\cos A-L^{\prime }\sin B\cos A){\mathbf{e}}_E\\ & -(B^{\prime }\cos A+L^{\prime }\cos B\sin A)\mathbf{n}\end{array}$$

${\mathbf{e}}_N$ 和 ${\mathbf{e}}_E$ 项的系数都是零！顺势将 $B’,L^{\prime }$ 关于 $A$ 的表达式代入：

$$\begin{array}{rl}\mathbf{K}& =-(B^{\prime }\cos A+L^{\prime }\cos B\sin A)\mathbf{n}\\ & =\boxed{\left(\frac{{\cos }^{2}A}{M}+\frac{{\sin }^{2}A}{N}\right)\mathbf{n}}\end{array}$$

最终我们得到大地线上某个点的曲率向量：

$$\mathbf{K}=\left(\frac{{\cos }^{2}A}{M}+\frac{{\sin }^{2}A}{N}\right)\mathbf{n}$$

曲率向量与该点的椭球法线平行，没有在椭球切向的分量。换言之，大地线始终只沿椭球面法线方向弯曲，切向上没有弯曲。这符合直觉，因为大地线是连接两点的最短曲线，如果在切向上弯曲，相当于「走了弯路」。

我们定义测地曲率 $k_g$ 为曲线曲率向量在切平面上的投影长，法曲率 $k_n$ 为曲率向量在法线上的投影长。那么对于大地线来说，其曲率等于法曲率，而测地曲率为零。也就是

$$\begin{array}{c}k=k_n=\frac{{\cos }^{2}A}{M}+\frac{{\sin }^{2}A}{N}\\ k_g=0\end{array}$$

这是大地线的其中一种定义方式：曲面上测地曲率处处为零的曲线。我们刚才从最短路径出发的另一种定义方式是：曲面上两点的局域最短路径。这两种定义是等价的。

法截线与大地线的关系

回忆一下，我们在 第 4 节求任意方向法截线的曲率半径时，求出的曲率也是 $\frac{{\cos }^{2}A}{M}+\frac{{\sin }^{2}A}{N}$！

但是我们考虑 $A$ 到 $B$ 的法截弧，它是 $A$ 点的法截线，在 $A$ 确实满足这个曲率公式。但一般情况下它不是 $B$ 点的法截线，在 $B$ 点处的方位角与 $B$ 点真正的法截线不同，因而曲率也不同。对于 $B$ 到 $A$ 的法截弧也一样，在 $B$ 点满足、$A$ 点不满足。可是大地线处处满足这个关系。所以我们终于得出结论：法截弧通常不是大地线、不是椭球面上两点的最短路径。

这个「处处满足」关系，可否导出一个更浅显的结论呢？

我们定义：过空间曲线上 $P$ 点的切线和 $P$ 点的邻近一点 $Q$ 可作一平面 ，当 $Q$ 点沿着曲线趋近于 $P$ 时，平面的极限位置为曲线在 $P$ 点的密切平面。

密切平面是该点处切向量 $\mathbf{T}$ 和曲率向量 $\mathbf{K}$ 张成的平面，这在直觉上很好理解。既然大地线的曲率向量平行于法向量，那么大地线密切平面必过法线。换言之，大地线上点的密切平面，始终是该点的法截面。

当两测站的距离在 $600\mathrm{km}$ 以内时，用两点之间的法截弧长度替代大地线长度的误差不超过千分之一毫米（Vanicek, 1986）。因此在小范围时为了计算方便，可以用法截弧近似替代大地线。

大地线的一阶微分方程

我们在推导法截线曲率向量的过程中，求 $A^{\prime }$ 的倒数第二步：

$$A^{\prime }=\frac{s\tan B}{N}$$

由此，我们可以通过 $B,L,A$ 三个量，表达大地线的一阶微分方程：

$$\boxed{\frac{\mathrm{d}B}{\mathrm{d}S}=\frac{\cos A}{M},\frac{\mathrm{d}L}{\mathrm{d}S}=\frac{\sin A}{N\cos B},\frac{\mathrm{d}A}{\mathrm{d}S}=\frac{\tan B}{N}\sin A}$$

小结

大地线的二阶微分方程

$$\begin{array}{c}\frac{{\mathrm{d}}^{2}B}{\mathrm{d}{S}^2}+\frac{3e^2\sin B\cos B}{W^2}{\left(\frac{\mathrm{d}B}{\mathrm{d}S}\right)}^2+\frac{W^2\sin B\cos B}{1-e^2}{\left(\frac{\mathrm{d}L}{\mathrm{d}S}\right)}^2=0\\ \frac{{\mathrm{d}}^{2}L}{\mathrm{d}{S}^2}-\frac{2(1-e^2)\tan B}{W^2}\frac{\mathrm{d}B}{\mathrm{d}S}\frac{\mathrm{d}L}{\mathrm{d}S}=0\end{array}$$

大地线的一阶微分方程

$$\frac{\mathrm{d}B}{\mathrm{d}S}=\frac{\cos A}{M},\frac{\mathrm{d}L}{\mathrm{d}S}=\frac{\sin A}{N\cos B},\frac{\mathrm{d}A}{\mathrm{d}S}=\frac{\tan B}{N}\sin A$$

大地线的曲率

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

大地线的性质

• 大地线上任何点的密切平面就是该点的法截面
• 曲面上连接任何两点的最短直线必为大地线
• 大地线的测地曲率等于零