中国地质大学(武汉)地质探测与评估教育部重点实验室，湖北 武汉

收稿日期：2023年2月22日；录用日期：2023年4月5日；发布日期：2023年4月12日

摘要

针对传统方法大地主题问题解算用第一偏心率e和K为参数的幂级数展开收敛速度慢、形式复杂效率低下的问题，以第三扁率n代替以往参数将其级数展开式进行重新推导改化。结果表明，基于n的大地线问题幂级数展开式更为简洁，形式上更加简单，更加便于分析与应用；在中小距离的大地主题反解上，新的表达式精度最高达到0.0001 m，满足高精度的要求；在长距离大地主题反解上，对传统人工推导的结果进行了新的计算，依然保持精度相对一致性。同时结果发现：经过转换后，在保持了上述精度的同时，展开式从原来的10阶降到了3阶，形式上更为简单。与传统的贝塞尔大地反解问题相比，在长距离大地线的解算上也存在着一定优势。

关键词

大地主题，第三扁率，幂级数展开式

The Third Flattening Expansion of Power Series Commonly Used in Geodetic Problem

Dongquan Zhou, Shaofeng Bian, Siyuan Zhang, Zhenghan Li

Key Laboratory of Geological Survey and Evaluation of Ministry of Education, China University of Geosciences (Wuhan), Wuhan Hubei

Received: Feb. 22^nd, 2023; accepted: Apr. 5^th, 2023; published: Apr. 12^th, 2023

ABSTRACT

In view of the problems of slow convergence and low efficiency of the power series expansion with the first flattening e and K as the parameters of the traditional method for solving the geodetic problem, the third flattening n was used to replace the previous parameters and the series expansion was reduced and modified. The results show that the power series expansion of the geodetic problem based on n is more concise, simpler in form, and more convenient for analysis and application. For the inverse solution of the earth theme in small and medium distance, the accuracy of the new expression is up to 0.0001 m, which meets the requirement of high precision. In the inverse solution of long distance geodetic theme, a new calculation is made to the result of traditional artificial derivation, and the accuracy is still consistent. At the same time, it is found that after the conversion, while maintaining the above accuracy, the expansion formula is reduced from the original order 10 to the order 3, which is simpler in form. Compared with the traditional Bessel geodetic inverse problem, it also has some advantages in the calculation of long distance geodetic line.

Keywords:Geodetic Problem, The Third Flattening, Power Series Expansion

Copyright © 2023 by author(s) and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY 4.0).

http://creativecommons.org/licenses/by/4.0/

1. 引言

大地线是椭球面上两点的最短曲线，其空间复杂性导致了大地主题问题的复杂性，一般来说大地主题问题包括正解和反解，正解问题是已知大地线长、起始点坐标和正方位角求终点坐标和反方位角，而反解问题是已知两点的大地坐标求大地线长和方位角 [1] 。图1展示了传统贝塞尔方法 [2] 的投影法则，其核心是将两点间的大地线投影到球面为大圆弧，使得投影前后的方位角不变，并利用大圆弧展开的以偏心率为自变量的幂级数进行大地线的求解，因此该方法与距离无关，在长距离的解算中较为适宜。

图1. 贝塞尔投影的几何要素

传统贝塞尔方法在具体的工程问题上有着较为广泛的应用场景，如海上航行、精确打击等方面。吴祖新等 [3] 参考贝塞尔方法的基本原理，利用智能船舶强大的计算能力和数据获取能力，构建了新型的大地线航法，在GPS的辅助下，实现了船舶航线的智能规划。同时，贝塞尔方程也存在需要改进的地方，传统贝塞尔方法忽略了高程对大地主题问题的影响，针对这种问题，国内外学者从顾及拱高误差 [4] 和考虑高程带来的误差 [5] 等方面出发，建立了新的参考椭球的算法。同时针对传统贝塞尔方法极点问题的求解，也有不少学者提出了新的思路，如大地线极点的归化纬度求解问题 [6] 。

近年来，随着计算机代数分析系统的不断发展，不仅提高了传统大地测量学的精度，而且革新了以往大地测量学的许多数学公式，这就为传统贝塞尔方法的改进准备了条件。边少锋等在《大地测量计算机代数分析》 [7] [8] 一书，对计算机代数系统在大地测量学上的新理论和新思路做出了系统的阐述，在大地线板块不仅对勒让德级数展开和高斯平均级数展开进行了创新讨论，还对传统的贝塞尔大地问题解算做出了改进，用 $\bar{e}$ 代替以往的K。针对传统的贝塞尔大地反解存在的奇异问题，由于繁琐的象限判定影响了解算的效率，就此问题周江华等 [9] [10] 提出了高效率的解决办法，并解决了传统贝塞尔的奇异问题。

对于贝塞尔大地方程，前人的研究解决了奇异问题，并提高了解算的效率，并且通常情况下，级数展开后的方程，展开式的项数越多，精度越高，同时也会带来了效率降低，计算复杂等问题，在实际用途中受到了限制 [11] 。在实际应用比如导弹弹道计算过程需要快速精准的解算结果，需要在满足精度的情况下，降低展开的阶数，提高解算效率。大地线主题还常常利用克莱罗方程进行计算，在椭球大地测量学中有着重要意义，但是同时存在着大地线常数C的求解问题 [12] [13] [14] 。如果得到一个已知的克莱罗常数将可以计算任意点上纬度、方位角和经度之间的精确关系，但是目前尚未有较好的解决办法，大部分大地主题问题依然采用传统的贝塞尔方程来进行求解。

综上所述，贝塞尔方程依然是大地主题研究的重要对象，而第三扁率的出现 [15] [16] [17] 可以进一步提高大地主题问题的计算效率。因此，本文在贝塞尔方程改进的基础上进行公式改化，改写为以第三扁率n为自变量的表达形式，简化其公式，同时保持其计算精度。为了探究公式改化后在精度上是否满足要求，本文与前人研究做出综合比对，对简化后公式的误差进行定量评估，最终给出贝塞尔方程新的实用公式。

2. 第三扁率n的定义

大地纬度B和归化纬度u存在以下关系：

$B=u+n\sin 2u+\frac{1}{2}{n}^2\sin 4u+\frac{1}{3}{n}^3\sin 6u-L$ (1)

其中

$n=\frac{1-\sqrt{1-e^2}}{1+\sqrt{1-e^2}}$ (2)

以上即为第三扁率n的定义，n的定义和拉格朗日共轭级数的引入可以提高式(1)的收敛速度，简化了大地纬度和规划纬度之间的相互转化。

同时，由此定义可知，第一偏心率e和第二偏心率 $e^{\prime }$ 可以分别表示为：

$e=\sqrt{\frac{4n}{{\left(1+n\right)}^2}}$ (3)

$e^{\prime }=\sqrt{\frac{4n}{{\left(1-n\right)}^2}}$ (4)

3. 大地线公式改化的数字描述

3.1. 大地线长与球面大圆弧长的微分方程迭代解

由前人的研究可知，大地线长S与归化纬度u的微分方程为 $\text{d}S=a\sqrt{1-e^2{\cos }^{2}u}\text{d}\sigma$ ， $\sigma$ 为大圆弧变量同时根据球面归化纬度与球面大圆弧长的数学关系式 $\sin u=\cos A_0\sin \sigma$ ，代入dS中，整理可得：

$\text{d}S=b\sqrt{1+{e^{\prime }}^2{\cos }^2{A}_0{\sin }^2\sigma }\text{d}\sigma$ (5)

设 ${\bar{e}}^2={e^{\prime }}^2{\cos }^2{A}_0$ ，其中 $e^{\prime }$ 为参考椭球的第二偏心率，则有

$\text{d}S=b\sqrt{1+{\bar{e}}^2{\sin }^2\sigma }\text{d}\sigma$ (6)

牛顿二项式展开定理的推广中，对于 ${\left(1+x\right)}^a$ 可以推广到非整数指数的情况，即：

${\left(1+x\right)}^a={\displaystyle \underset{k=0}{\overset{\infty }{\sum }}{C}_{\alpha }^k{x}^k}={\displaystyle \underset{k=0}{\overset{\infty }{\sum }}\frac{{\left(\alpha \right)}_k}{k!}{x}^k}$ (7)

其中

${\left(\alpha \right)}_k={\displaystyle \underset{i=\frac{1}{2}-k+1}{\overset{\alpha }{\prod }}i}$ (8)

因此将(7)代入(6)可展得(此处只展至 ${\bar{e}}^{10}$ )

$\text{d}S=b\left(1+\frac{1}{2}{\bar{e}}^2\sin {\left[x\right]}^2-\frac{1}{8}{\bar{e}}^4\sin {\left[x\right]}^4+\frac{1}{16}{\bar{e}}^6\sin {\left[x\right]}^6-\frac{5}{128}{\bar{e}}^8\sin {\left[x\right]}^8+\frac{7}{256}{\bar{e}}^{10}\sin \left[x\right]\right)$ (9)

上式积分并化为三角函数的倍角形式得

$\begin{array}{c}S=b{\displaystyle {\int }_{{\sigma }_1}^{{\sigma }_2}\sqrt{1+{\bar{e}}^2{\sin }^2\sigma }\text{d}\sigma }\\ =b{\left[\alpha \sigma +\beta \sin 2\sigma +\gamma \sin 4\sigma +\delta \sin 6\sigma +\epsilon \sin 8\sigma +ϵ\sin 10\sigma \right]|}_{{\sigma }_1}^{{\sigma }_2}\end{array}$ (10)

式中

$\{\begin{array}{l}\alpha =1+\frac{{\bar{e}}^2}{4}-\frac{3{\bar{e}}^4}{64}+\frac{5{\bar{e}}^6}{256}-\frac{175{\bar{e}}^8}{16384}+\frac{441{\bar{e}}^{10}}{65536}\\ \beta =-\frac{{\bar{e}}^2}{8}+\frac{{\bar{e}}^4}{32}-\frac{15{\bar{e}}^6}{1024}+\frac{35{\bar{e}}^8}{4096}-\frac{735{\bar{e}}^{10}}{131072}\\ \gamma =-\frac{{\bar{e}}^4}{256}+\frac{3{\bar{e}}^6}{1024}-\frac{35{\bar{e}}^8}{16384}+\frac{105{\bar{e}}^{10}}{65536}\\ \delta =-\frac{{\bar{e}}^6}{3072}+\frac{5{\bar{e}}^8}{12288}-\frac{105{\bar{e}}^{10}}{262144}\\ \epsilon =-\frac{5{\bar{e}}^8}{131072}+\frac{35{\bar{e}}^{10}}{524288}\\ ϵ=-\frac{7{\bar{e}}^{10}}{1310720}\end{array}$ (11)

根据第三扁率的定义，可知

$\{\begin{array}{l}{e}^2=\frac{4n}{{\left(1+n\right)}^2}\\ {e^{\prime }}^2=\frac{4n}{{\left(1-n\right)}^2}\end{array}$ (12)

因此可知

${\bar{e}}^2=\frac{4n{\cos }^2{A}_0}{{\left(1-n\right)}^2}$ (13)

代入微分方程 $dS$ 对其展开并积分可得：

$S=b\ast A\ast B$ (14)

其中

$A=\left[\begin{array}{cccc}1& n+2n^2+3n^3& -\frac{3n^2}{4}-3n^3& \frac{5n^3}{4}\\ 0& -\frac{1}{2}n-n^2-\frac{3}{2}{n}^3& \frac{n^2}{2}+2n^3& -\frac{15n^3}{16}\\ 0& 0& -\frac{n^2}{16}-\frac{n^3}{4}& \frac{3n^3}{16}\\ 0& 0& 0& -\frac{n^3}{48}\end{array}\right]\left[\begin{array}{c}1\\ \cos {\left[A_0\right]}^2\\ \cos {\left[A_0\right]}^4\\ \cos {\left[A_0\right]}^6\end{array}\right]$ (15)

$B=\left[\begin{array}{cccc}{\sigma }_2-{\sigma }_1& \sin \left[2{\sigma }_2\right]-\sin \left[2{\sigma }_1\right]& \cdots & \sin \left[6{\sigma }_2\right]-\sin \left[6{\sigma }_1\right]\end{array}\right]$ (16)

对比式(15)和(11)可以发现，以第三扁率n代替原来的公式，让展开式更加紧凑简洁，同时计算效率得到了提高，经检验，只需要展开至n³即可满足中短距离的大地线计算，很好地适用于短距离，

在2000 km及以上的大地线计算中，也比以往的参考结果更加精准，因此该公式的适用范围得到了提高。

3.2. 大地线长与球面大圆弧长微分方程的非迭代解

以上导出了 $\sigma$ 与S的关系，但是由于现实问题中往往是已知S求 $\sigma$ ，上式仍然需要进行迭代运算，为克服迭代运算带来的运算复杂等问题，需要导出S与 $\sigma$ 的无迭代的直接解法，对此，首先定义

$\phi =\frac{S}{\alpha b}=\frac{S}{b\left(1+\frac{{\bar{e}}^2}{4}-\frac{3{\bar{e}}^4}{64}+\frac{5{\bar{e}}^6}{256}-\frac{175{\bar{e}}^8}{16384}+\frac{441{\bar{e}}^{10}}{65536}\right)}$ (17)

式中 $\alpha$ 由上文所求参考系数

同时，可设 $\sigma$ 与S的直接解公式的形式为

$\sigma =\phi +c_2\sin 2\phi +c_4\sin 4\phi +\cdots +c_{10}\sin 10\phi$ (18)

求式(18)的待定系数需要确定5个插值条件，进行多次求导得：

$\{\begin{array}{l}\frac{\text{d}\sigma }{\text{d}\phi }\left(0\right)=1+2c_2+\cdots +10c_{10}\\ \frac{{\text{d}}^3\sigma }{\text{d}{\phi }^3}\left(0\right)=-8c_2-64c_4-\cdots -1000c_{10}\\ \frac{{\text{d}}^5\sigma }{\text{d}{\phi }^5}\left(0\right)=32c_2+1024c_4+\cdots +100000c_{10}\\ \frac{{\text{d}}^7\sigma }{\text{d}{\phi }^7}\left(0\right)=-128c_2-16384c_4-\cdots -10000000c_{10}\\ \frac{{\text{d}}^9\sigma }{\text{d}{\phi }^9}\left(0\right)=512c_2+262144c_4+\cdots +1000000000c_{10}\end{array}$ (19)

同时根据式(17)和式(6)，将 $\sigma$ 对 $\phi$ 求导，即

$\frac{\text{d}\sigma }{\text{d}\phi }=\frac{ab}{dS}\ast \text{d}\sigma =\frac{ab}{b\sqrt{1+{\bar{e}}^2{\sin }^2\sigma }\text{d}\sigma }\text{d}\sigma =\frac{a}{\sqrt{1+{\bar{e}}^2{\sin }^2\sigma }}$ (20)

由式(18)当 $\phi =0$ 时，易知 $\sigma =0$ ，式(20)根据隐函数的求导法则可对 $\sigma$ 进行多次求导得

$\{\begin{array}{l}\frac{\text{d}\sigma }{\text{d}\phi }\left(0\right)=\alpha \\ \frac{{\text{d}}^3\sigma }{\text{d}{\phi }^3}\left(0\right)=-e^2{\alpha }^3\\ \frac{{\text{d}}^5\sigma }{\text{d}{\phi }^5}\left(0\right)=e^2\left(4+13e^2\right){\alpha }^5\\ \frac{{\text{d}}^7\sigma }{\text{d}{\phi }^7}\left(0\right)=-e^2\left(16+\cdots +493e^4\right){\alpha }^7\\ \frac{{\text{d}}^9\sigma }{\text{d}{\phi }^9}\left(0\right)=e^2\left(64+4944e^2+\cdots +37369e^6\right){\alpha }^9\end{array}$ (21)

对式(19)和式(21)其联立矩阵方程，解矩阵方程可得

$\left(\begin{array}{c}{c}_2\\ c_4\\ c_6\\ c_8\\ c_{10}\end{array}\right)={\left(\begin{array}{ccccc}2& \cdots & \cdots & \cdots & 10\\ -8& \cdots & \cdots & \cdots & -1000\\ 32& \cdots & \cdots & \cdots & 100000\\ -128& \cdots & \cdots & \cdots & -10000000\\ 512& \cdots & \cdots & \cdots & 1000000000\end{array}\right)}^{-1}\left(\begin{array}{c}\alpha -1\\ -e^2{\alpha }^3\\ e^2\left(4+13e^2\right){\alpha }^5\\ -e^2\left(16+\cdots +493e^4\right){\alpha }^7\\ e^2\left(64+\cdots +37369e^6\right){\alpha }^9\end{array}\right)$ (22)

式(17)中对 $\alpha$ 给出了定义，求出矩阵的解并代入 $\alpha$ ，通过级数展开的方式化简，可得最终结果如下：

$\left(\begin{array}{c}{c}_2\\ c_4\\ c_6\\ c_8\\ c_{10}\end{array}\right)=\left(\begin{array}{c}\frac{{\bar{e}}^2}{8}-\frac{{\bar{e}}^4}{16}+\frac{71{\bar{e}}^6}{2048}-\frac{85{\bar{e}}^8}{4096}+\frac{20797{\bar{e}}^{10}}{1572864}\\ \frac{5{\bar{e}}^4}{256}-\frac{5{\bar{e}}^6}{256}+\frac{383{\bar{e}}^8}{24576}-\frac{143{\bar{e}}^{10}}{12288}\\ \frac{29{\bar{e}}^6}{6144}-\frac{29{\bar{e}}^8}{4096}+\frac{969{\bar{e}}^{10}}{131072}\\ \frac{539{\bar{e}}^8}{393216}-\frac{539{\bar{e}}^{10}}{196608}\\ \frac{3467e^{10}}{7864320}\end{array}\right)$ (23)

对于式(23)将其代入式(18)可求球面大圆弧的直接解公式，现在式(23)进行改化，表示为第三扁率n的形式的幂级数形式为：

$\left(\begin{array}{c}{c}_2\\ c_4\\ c_6\\ c_8\\ c_{10}\end{array}\right)=M\ast N$

$M=\left[\begin{array}{ccccc}\frac{n}{2}+n^2+\frac{3n^3}{2}+2n^4+\frac{5n^5}{2}& -n^2-4n^3-10n^4-20n^5& \frac{71n^3}{32}+\frac{213n^4}{16}+\frac{1491n^5}{32}& -\frac{20797n^4}{6144}-\frac{20797n^5}{768}& \frac{697613n^5}{61440}\\ 0& \frac{5n^2}{16}+\frac{5n^3}{4}+\frac{25n^4}{8}+\frac{25n^5}{4}& -\frac{5n^3}{4}-\frac{15n^4}{2}-\frac{105n^5}{4}& \frac{383n^4}{96}+\frac{383n^5}{12}& -\frac{143n^5}{12}\\ 0& 0& \frac{29n^3}{96}+\frac{29n^4}{16}+\frac{203n^5}{32}& -\frac{29n^4}{16}-\frac{29n^5}{2}& \frac{969n^5}{128}\\ 0& 0& 0& \frac{539n^4}{1536}+\frac{539n^5}{192}& -\frac{539n^5}{192}\\ 0& 0& 0& 0& \frac{3467n^5}{7680}\end{array}\right]$ (24)

$N=\left[\begin{array}{c}\cos {\left[A_0\right]}^2\\ \cos {\left[A_0\right]}^4\\ \cos {\left[A_0\right]}^6\\ \cos {\left[A_0\right]}^8\\ \cos {\left[A_0\right]}^{10}\end{array}\right]$ (25)

根据式(24)与式(23)相比可以看出，以第三扁率n来进行幂级展开，以矩阵的形式在保持精度的情况，使其系数更小，指数更低，使得形式更加紧凑，同时以克拉索夫斯基椭球为例，n约等于0.0015，而 $\cos \left[A_0\right]$ 则是−1到1的分数，为在实际应用上随着指数的提高是一个很小的数，在实际的应用中，在精度要求不是特高的大地线计算如无线电定位中，还可以化为对角矩阵，进一步对系数进行简化，即式(26)：

$M^{\prime }=\left[\begin{array}{ccccc}\frac{n}{2}+n^2+\frac{3n^3}{2}+2n^4+\frac{5n^5}{2}& 0& 0& 0& 0\\ 0& \frac{5n^2}{16}+\frac{5n^3}{4}+\frac{25n^4}{8}+\frac{25n^5}{4}& 0& 0& 0\\ 0& 0& \frac{29n^3}{96}+\frac{29n^4}{16}+\frac{203n^5}{32}& 0& 0\\ 0& 0& 0& \frac{539n^4}{1536}+\frac{539n^5}{192}& 0\\ 0& 0& 0& 0& \frac{3467n^5}{7680}\end{array}\right]$ (26)

式(26)与式(25)以克拉索夫斯基椭球为参数计算可估算最大误差为式(27)，由此可见，在大部分的应用场景中，式(26)以简略的形式即可满足要求。

$\left(\begin{array}{c}\Delta c_2\\ \Delta c_4\\ \Delta c_6\\ \Delta c_8\\ \Delta c_{10}\end{array}\right)=\left(\begin{array}{c}-2.26964\times {10}^{-6}\\ -4.27501\times {10}^{-9}\\ -9.34022\times {10}^{-14}\\ -2.16408\times {10}^{-14}\\ 0\end{array}\right)$ (27)

3.3. 大椭圆法解大地问题

由于在一个象限内，大椭圆法与贝塞尔法的精度差别不大，因此在一些特殊情况往往也会采用大椭圆法来进行大地问题求解。

图2. 大椭圆法的几何要素

图2为大椭圆法求大地主题问题时的几何要素，其设立经过椭球面和椭球中心一个象限的截面椭圆，长半轴与地球椭球相等，设截面椭圆的短半轴为 $\bar{b}$ ，第一偏心率为 $\bar{e}$ ，以大地线最高点为极轴，则 $\theta$ 为该点到赤道面的极角， $\rho$ 为该点的极径，首先设截面椭圆以 $\left(\rho ,\theta \right)$ 为变量的极坐标方程：

$\rho =\frac{\bar{b}}{\sqrt{1-{\bar{e}}^2{\cos }^2\theta }}$ (28)

对其求导可得：

${\rho }^{\prime }=-\frac{\bar{b}{\bar{e}}^2\sin \theta \cos \theta }{{\left(1-{\bar{e}}^2{\cos }^2\theta \right)}^{\frac{3}{2}}}$ (29)

其中

$\bar{b}=\frac{b}{\sqrt{1-e^2{\cos }^2{\varphi }_n}}=\frac{b}{\sqrt{1-e^2{\cos }^2{\varphi }_1{\sin }^2{A}_1}}$ (30)

$\bar{e}=\sqrt{\frac{a^2-{\bar{b}}^2}{a^2}}=e\sqrt{\frac{1-{\cos }^2{\varphi }_1{\sin }^2{A}_1}{1-e^2{\cos }^2{\varphi }_1{\sin }^2{A}_1}}$ (31)

式中 $\varphi$ 为地心纬度

对其进行椭圆积分，可得大地线

$\begin{array}{c}S={\displaystyle {\int }_{{\theta }_1}^{{\theta }_2}\text{d}S}={\displaystyle {\int }_{{\theta }_1}^{{\theta }_2}\sqrt{{\rho }^{\prime }+{\rho }^2}\text{d}\theta }\\ =\bar{b}{\left[\alpha \theta +\beta \sin 2\theta +\gamma \sin 4\theta +\delta \sin 6\theta +\epsilon \sin 8\theta +ϵ\sin 10\theta \right]|}_{{\theta }_1}^{{\theta }_2}\end{array}$ (32)

式中

$\{\begin{array}{l}\alpha =1+\frac{{\bar{e}}^2}{4}+\frac{13{\bar{e}}^4}{64}+\frac{45{\bar{e}}^6}{256}+\frac{2577{\bar{e}}^8}{16384}+\frac{9417{\bar{e}}^{10}}{65536}\\ \beta =\frac{1{\bar{e}}^2}{8}+\frac{3{\bar{e}}^4}{32}+\frac{95{\bar{e}}^6}{1024}+\frac{385{\bar{e}}^8}{4096}+\frac{12351{\bar{e}}^{10}}{131072}\\ \gamma =-\frac{{\bar{e}}^4}{256}-\frac{5{\bar{e}}^6}{1024}-\frac{19{\bar{e}}^8}{16384}+\frac{249{\bar{e}}^{10}}{65536}\\ \delta =-\frac{5{\bar{e}}^6}{1024}-\frac{35{\bar{e}}^8}{4096}-\frac{2551{\bar{e}}^{10}}{262144}\\ \epsilon =-\frac{261{\bar{e}}^8}{131072}-\frac{2349{\bar{e}}^{10}}{524288}\\ ϵ=-\frac{921{\bar{e}}^{10}}{1310720}\end{array}$ (33)

同时，大地线两点在截面椭圆上的极角 ${\theta }_2$ 和 ${\theta }_1$ 分别为：

$\{\begin{array}{l}\cot {\theta }_1=\frac{\cos A_1}{\tan {\varphi }_1}\\ \cot {\theta }_2=\frac{\cos A_2}{\tan {\varphi }_2}\end{array}$ (34)

按上述理论，同样可将式(33)进行改化，将其表示为第三扁率n的幂级数形式：

$S=\bar{b}\ast C\ast D$

其中

$C=\left[\begin{array}{cccc}1& n+2n^2+3n^3& \frac{13n^2}{4}+13n^3& \frac{45n^3}{4}\\ 0& \frac{n}{2}+n^2+\frac{3n^3}{2}& \frac{3n^2}{2}+6n^3& \frac{95n^3}{16}\\ 0& 0& -\frac{n^2}{16}-\frac{n^3}{4}& -\frac{5n^3}{16}\\ 0& 0& 0& -\frac{n^3}{16}\end{array}\right]\left[\begin{array}{c}1\\ \cos {\left[A_0\right]}^2\\ \cos {\left[A_0\right]}^4\\ \cos {\left[A_0\right]}^6\end{array}\right]$ (35)

$D=\left[\begin{array}{ccc}{\theta }_2-{\theta }_1& \cdots & \sin \left[{\theta }_2\right]-\sin \left[6{\theta }_1\right]\end{array}\right]$ (36)

在式(31)的基础上，对其以第三扁率n代替原来的公式可以让原来的展开式更加紧凑简洁，同时提高计算效率，形成4 × 4的上三角系数矩阵。以上推导可以看出，由于 $\bar{b}$ 受到三个不确定参数的影响，与贝塞尔方法相比，适用场景较少，但是同样能够满足无线电双曲线导航定位的要求，在长距离计算中，在米级、十米级和百米级都能保持较好的精度。

3.4. 大椭圆法与贝塞尔方程精度比较

大椭圆法常常用于一个象限的大地线问题求解，与贝塞尔方程相比，精度有所差异，以式(35)和式(15)做比较，以克拉索夫斯基椭球为参数对双方误差作最大估计得式(37)，由式(37)可知，由于在一个象限内， $\theta$ 与 $\sigma$ 的差别很小，同时 $\bar{b}$ 受地心纬度和方位角影响，误差在一个象限内也受两点位置的影响，在忽略 $\theta$ 与 $\sigma$ 、 $\bar{b}$ 和b的误差的情况下，大椭圆法与贝塞尔方程相比，误差在厘米级，当由于不定因素过多，当两点距离变大时，误差可达米级，当两点超过一个象限的距离时，与贝塞尔方程比，将不再适用。因此下文将以贝塞尔方程为主，对其精度进行数学分析。

$\Delta =\left(\bar{b}-b\right)\ast \left(\begin{array}{c}9.14283\times {10}^{-6}\\ 0.151135\times {10}^{-2}\\ 1.70279\times {10}^{-9}\\ 1.41899\times {10}^{-10}\end{array}\right)\ast \left(D-B\right)$ (37)

4. 大地主题解算计算实例

在上述公式的数学描述中，可以看出贝塞尔方法在大椭圆法中的优势，同时精度也优于大椭圆法，因此本文对改化后的贝塞尔大地问题解算进行互相一致性的精度检验。为了进一步验算本文算法的精度和可靠性，本文结合3组具有代表性的算例，在克拉索夫斯基椭球的参数下采用双精度运算，对大地主题反算的实例进行互相一致性检验。

表1. 大地主题计算模型验证

本文在参考前人科学研究的基础上，结合已有的参考算例，在Mathematica软件上进行算法实现和结果计算：

表1中参考结果1来源于文献 [9] 和《大地测量学基础》与《大地坐标系与大地基准》两本书籍，参考结果2为将《大地测量计算机代数分析》中贝塞尔方程以e为变量的幂级展开式进一步推导到十阶并计算的结果。

算例1可以看出在短距离的计算中，可以保持非常好的精度水平，同时算例1的参考结果1计算的是法截长，在距离很小的时候法截线与大地线的差别也非常小。在50 km的计算精度可以达到0.001 m，在2000 km的距离上，精度最高可以实现达到0.0001 m。算例1还可以看出，在计算机代数系统下，距离很接近的两个点，法截线与大地线差别不大，同时对比分析发现总体精度能保持在0.001 m以内，达到毫米级，对于长距离的大地线计算，引入大地方位角的修正项，在计算机代数系统辅助下，本文研究结果可信度高，同时精度最高能达到0.0001 m，能满足各种场景的实用，特别是高精度需求的领域，或许也是一个重要参考。

5. 结论

1) 本文通过引用第三扁率n对现有的大地线算法进行了改化，在保持精度的同时，简化了公式结构，使得计算效率得到了较大的提高，在大地主题反解大地线的计算中，本文参照经典算例，进行了互相一致性检验，结果表明，新的改化公式可以更好地满足各种场景的应用，同时精度总体达到0.001 m，最高可达0.0001 m。

2) 本文不仅基于贝塞尔公式进行了改化，同时对大椭圆法求解大地线主题问题进行了初步探究，公式表明大椭圆法在精度上表现不优，但是在十米级或百米级的定位误差应用中，如无线电双曲线导航定位中，依旧有着独特优势。

基金项目

国家自然科学基金项目(42074010)。

文章引用

周东权,边少锋,张思远,李正汉. 大地主题问题常用幂级数的第三扁率展开
The Third Flattening Expansion of Power Se-ries Commonly Used in Geodetic Problem[J]. 测绘科学技术, 2023, 11(02): 98-108. https://doi.org/10.12677/GST.2023.112011

参考文献