中国地质大学（武汉）周东权：子午线弧长正解展开通式及数学分析 |《测绘学报》2024年53卷第9期

本文内容来源于《测绘学报》2024年第9期（审图号GS京(2024)1896号）

周东权,1,2, 边少锋,2, 黄晓颖3

1.中国地质大学（武汉）地理与信息工程学院，湖北　武汉　430078

2.中国地质大学（武汉）地质探测与评估教育部重点实验室，湖北　武汉　430074

3.海军工程大学作战运筹与规划系，湖北　武汉　430033

基金项目

 国家自然科学基金(42342024); 中央高校基金(GLAB2024ZR05)

作者简介

作者简介：周东权（1999—），男，硕士，研究方向为椭球大地测量。E-mail：1431645283@qq.com

通讯作者: 边少锋　E-mail：sfbian@sina.com

摘要

子午线弧长公式通常表示为大地纬度B的级数展开式，其系数以第一偏心率e为参数，本文利用第三扁率n对公式进行重新推导，并将其表示为三角函数的3种形式：倍角形式、指数形式和二倍角形式。重新推导的公式中各系数中的分母值都明显变小，个别高阶项系数消失，结构简单、形式简洁。基于此，对3种表示形式的8阶展开式和10阶展开式进行截断误差分析，结果表明3种表示形式的8阶展开式精度最低也有毫米级，能满足日常的使用场景，而10阶展开式精度提高了至少2个数量级，能满足高精度的使用场景。在高精度的使用场景下，不同表示形式的子午线弧长公式具有不同的实用性，分析表明在0°N—30°N低纬度地区推荐使用三角函数的指数形式，在30°N—55°N中纬度地区推荐使用三角函数的二倍角形式，而在55°N—90°N高纬度地区推荐使用三角函数的倍角形式。

关键词

子午线弧长; 第三扁率; 数学分析; 区域分析

本文引用格式

周东权, 边少锋, 黄晓颖. 子午线弧长正解展开通式及数学分析[J]. 测绘学报, 2024, 53(9): 1790-1798 doi:10.11947/j.AGCS.2024.20230546

ZHOU Dongquan, BIAN Shaofeng, HUANG Xiaoying. General series expansion and mathematical analysis of the direct solution of meridian arc length[J]. Acta Geodaetica et Cartographica Sinica, 2024, 53(9): 1790-1798 doi:10.11947/j.AGCS.2024.20230546

子午线弧长的计算是椭球大地测量学中最基本的数学问题[1-4]，是高斯投影经差展开公式的必算项。同时，子午线弧长计算本质上是第二类椭圆积分[5-9]，设椭球长半轴为a，第一偏心率为e，则椭球面上由赤道到大地纬度B处的子午线弧长X为

(1)

过去的研究中，学者对子午线弧长进行了大量的推导，其正解展开式通常表示为三角函数的倍角形式[10-16]。同时子午线弧长公式系数大多使用第一偏心率e作为变量，在此形式下，收敛速度慢且系数形式复杂。针对这一问题，国内外已有多位学者使用第三扁率n代替第一偏心率e重新推导公式系数，其定义如下[17]

(2)

式中，a和b分别为椭球的长、短半轴。

根据第三扁率的定义式可得第三扁率与第一偏心率的关系式

(3)

以第三扁率为参数的级数展开式往往具有更快的收敛速度[18]。过去学者们对于第三扁率的研究较少，在子午线弧长研究领域，文献[19-20]曾使用第三扁率重新表示了子午线弧长正解公式。近年来，第三扁率开始被广泛应用于大地测量学的各项研究中，国外大地测量学家在椭球大地测量学的相关研究中均引入了第三扁率[21-23]，结果表明，所得公式系数结构更加简洁。国内学者也使用第三扁率对收敛速度和计算精度进行了系统研究[24-27]。由此可见，第三扁率的引入能够有效简化公式系数且提高收敛速度，但是对于子午线弧长的展开式，目前还停留在对系数的改化上，并未引入第三扁率对公式进行重新推导并形成常用的通式形式。

虽然子午线弧长是一个传统的研究问题，但是借助现代强大的计算机代数工具[28]，本文从数学分析的角度系统剖析子午线弧长常用表达式的简化表示，在第三扁率展开的基础上进一步简化公式，并数学分析3种形式表达式的适用性，最后以高精度使用场景为标准，通过对不同形式子午线弧长的正解公式进行严谨的数学分析，为实际的大地测量学应用提供帮助与建议。

1　子午线弧长关于大地纬度的3种表示形式

1.1　三角函数的倍角形式

子午线弧长往往可以表示成大地纬度B的三角函数的倍角形式的展开式，考虑大地测量学高精度的计算要求，本文通过计算机代数对其展开式从8阶进一步拓展至10阶，可得

(4)

式中，系数为

(5)

使用第三扁率n代替第一偏心率e，进行重新推导，得到以第三扁率n表示的子午线弧长正解展开公式

(6)

式中，系数为

(7)

可以明显看出，第三扁率n为参数的子午线正解公式中的各系数的分母值都明显变小，个别高阶项系数消失，形式非常简洁。因此后续将以第三扁率n为参数进行推导。

1.2　三角函数的指数形式

为了避免复杂的多倍角计算，可将三角函数倍角计算转换成三角函数指数计算。由此可得三角函数的指数形式的子午线弧长正解展开式

(8)

式中，系数为

(9)

1.3　三角函数的二倍角形式

为了避免复杂的多倍角计算，除了1.2节中使用的将三角函数的倍角形式转换成三角函数的指数形式计算的方法外，还可改为三角函数的二倍角形式

(10)

式中，系数为

(11)

可以发现，式(11)第三扁率n重新推导后的子午线弧长正解公式系数与传统公式相比更加简单，个别高阶项系数消失，形式简洁，这是由第三扁率的定义决定的，第三扁率值更小，收敛更快。在此基础上，将对以第三扁率n重新推导后的子午线弧长正解公式进行系统的数学分析。

2　3种表示形式的数学分析

2.1　子午线弧长正解展开式通式

为了系统地数学分析，需要建立不同形式的计算通式，结果见表1。表1为各种形式展开式的计算通式，以大地纬度为变量，在以第三扁率进行重新推导时，保持了前半部分公式表达一致，后半部分的级数展开存在差别，因此误差由后半部分产生，为后面的误差分析提供了很好的思路。

表1   子午线正解展开通式

Tab.1  The general series expansion of direct solution of meridian arc length

表示形式	正解展开通式
三角函数的倍角形式
三角函数的指数形式
三角函数的二倍角形式

注：式中t和M为自然数。

新窗口打开| 下载CSV

2.2　子午线弧长正解展开式截断误差通式

子午线弧长展开式理论上可以无限展开，但在实际大地测量应用场景中，一般会在保证精度的情况下，尽可能地让级数展开式简洁，由此会产生一定的计算误差，该误差称为(高阶)截断误差。因为第三扁率的级数收敛很快，为简单起见，只将截断的第一项作为截断误差。

通常来说，子午线弧长展开式级数展开至8阶或许就可以满足普通使用场景的精度要求，但对于高精度的使用场景，往往可能需要展开至10阶以上。基于此，在适用性分析前，需要先探究以第三扁率n重新推导后的子午线弧长正解公式在8阶展开式和10阶展开式的精度情况。设子午线弧长正解公式的截断误差为ε，结合表1可得子午线弧长正解展开式截断误差的计算通式见表2。

表2   子午线弧长正解展开式截断误差计算通式

Tab.2  The general formula for calculating the truncation error of the series expansion of the direct solution of meridian arc length

表示形式	8阶截断误差计算通式	10阶截断误差计算通式
三角函数的倍角形式	ε11≈a(1－n)b10sin(10B)	ε12≈a(1－n)b12sin(12B)
三角函数的指数形式	ε21≈a(1－n)cos B×c9sin9B	ε22≈a(1－n)cos B×c11sin11B
三角函数的二倍角形式	ε31≈a(1－n)sin(2B)×d5cos4(2B)	ε32≈a(1－n)sin(2B)×d6cos5(2B)

新窗口打开| 下载CSV

结合式(7)、式(9)和式(11)可计算8阶、10阶截断误差计算通式的系数

(12)

2.3　截断误差极值分析

对子午线弧长公式的截断误差进行极值分析需要对通式进行求导，可得导数ε'见表3。

表3   子午线弧长正解展开式截断误差计算通式一阶导数

Tab.3  The first-order derivative of the general formula for calculating the truncation error in the series expansion of the direct solution of meridian arc length

表示形式	截断误差计算通式
三角函数的倍角形式
三角函数的指数形式
三角函数的二倍角形式

新窗口打开| 下载CSV

令ε'=0，可得极值点通式(表4)。

表4   子午线弧长正解展开式截断误差极值点通式

Tab.4  The general formula for the extreme points of the truncation error in the series expansion of the direct solution of meridian arc length

表示形式	8阶截断误差极值点通式	10阶截断误差极值点通式
三角函数的倍角形式
三角函数的指数形式
三角函数的二倍角形式

新窗口打开| 下载CSV

当一阶导数为0、二阶导数不为0时，该点为极值点。表4为3种表示形式子午线弧长正解展开式8阶截断误差和10阶截断误差的极值点通式，顾及大地纬度0≤B≤π/2，并通过计算二阶导数值，可以判断是否为极值点，最终得到子午线弧长正解展开式截断误差极值统计见表5。

表5   子午线弧长正解展开式截断误差极值统计

Tab.5  The extreme value statistics of truncation error in the series expansion of the direct solution of meridian arc length

表示形式	阶数	极值点/rad	极值/m
三角函数的倍角形式	8阶		±4.244 18×10－7
	10阶		±6.977 73×10－11
三角函数的指数形式	8阶		－0.000 046 386
	10阶		2.556 33×10－8
三角函数的二倍角形式	8阶		±2.107 85×10－8
	10阶		±5.778 77×10－10

新窗口打开| 下载CSV

由表5可以看出8阶展开式下，三角函数的倍角形式展开式存在5个极值点，极值为±4.244 18×10-7 m；三角函数的指数形式展开式则只存在一个极小值点，极值为-0.000 046 386 m；三角函数的二倍角形式存在两个极值点，极值为2.107 85×10-8 m。而在10阶展开式下，三角函数的倍角形式展开式存在6个极值点，极值为±6.977 73×10-11 m；三角函数的指数形式展开式依旧存在一个极小值点，极值为2.556 33×10-8 m；三角函数的二倍角形式依旧存在两个极值点，极值为±5.778 77×10-10 m。

截断误差分析表明，在8阶展开式下，3种形式最大误差都能在毫米级以下，满足大部分日常使用的场景，但是对于高精度的使用场景来说，三角函数的指数形式的级数展开式往往需要在更高阶的级数展开；在10阶展开式时，三角函数的倍角形式和指数形式误差极值精度提高了3个数量级，三角函数的二倍角形式误差极值精度提高了2个数量级，可以看出，已经可以满足高精度的使用要求。在此基础上，以高精度的使用场景为标准，探究子午线弧长公式在10阶展开式下不同区域各种形式表达式的实用性。

以北半球为研究区域，对子午线弧长公式在10阶展开式下截断误差空间分布进行了分析，如图1所示，可以看出，三角函数的倍角形式下子午线弧长正解展开式波动较为频繁，呈周期性，但波峰不高，波动不大；三角函数的指数形式下，存在一个比较高的波峰，误差存在波动较大的情况；三角函数的二倍角形式下，存在两个波峰，但是波峰较低，总体变化较为平缓。

图1

图1   子午线弧长正解展开式截断误差空间分布

Fig.1   The spatial distribution of truncation error in the series expansion of the direct solution of meridian arc length

在此基础上，对3种形式子午线弧长正解展开式截断误差进行组合分析，对不同区域子午线正解展开式进行分析，如图2—图4所示，为方便比较，对靠近x轴的曲线进行了放大处理。由图2可以看出，三角函数的指数形式在0°N—30°N低纬度地区精度非常好，基本趋于0，这是由于B≤30°时候sinB是一个较小的值，三角函数的指数形式下的值更加小，因此误差非常小，同时发现三角函数的指数形式子午线弧长公式比倍角形式误差几乎小了3倍以上，适用性更强。由图3可以看出，在30°N—55°N中纬度地区，三角函数的二倍角形式和倍角形式子午线正解展开式精度更高，但三角函数的二倍角形式更加平稳，这是由于在B→45°时cos(2B)→0，同时发现三角函数的二倍角形式比倍角形式误差也几乎小了2倍以上，适用性更强。由图4可以看出，在55°N—90°N高纬度地区，三角函数的倍角形式精度则更加平稳，这是由于在B→90°时sin(2nB)→0，同时发现三角函数的倍角形式比二倍角形式误差几乎小了6倍以上，适用性更强。

图2

图2   0°N—30°N低纬度地区子午线弧长正解展开式截断误差分析

Fig.2   The analysis of truncation error in the series expansion of the direct solution of meridian arc lengthat low latitudes region of 0°N—30°N

图3

图3   30°N—55°N中纬度地区子午线弧长正解展开式截断误差分析

Fig.3   The analysis of truncation error in the series expansion of the direct solution of meridian arc length at mid-latitudes region of 30°N—55°N

图4

图4   55°N—90°N高纬度地区子午线弧长正解展开式截断误差分析

Fig.4   The analysis of truncation error in the series expansion of the direct solution of meridian arc length at high latitudes region of 55°N—90°N

我国是一个南北跨度较大的国家，最北境位于漠河以北的黑龙江主航道的中心线上(53°34'N)，最南境在南沙群岛的曾母暗沙附近(3°51'N)，根据研究可以发现在我国建立区域性的地图投影，利用三角函数的指数形式和二倍角形式进行子午线弧长的正解展开式能实现更加精确的投影表示。

3　结论

子午线弧长的正解公式是第二类椭圆积分，无法用常规的积分方法直接求解，在求解过程中往往使用级数展开的方法进行积分，在此基础上，针对以第一偏心率展开的公式结构复杂、系数冗长等问题，以第三扁率对子午线弧长公式进行重新推导，并将其表示成三角函数的倍角形式、指数形式和二倍角形式3种形式，最终以系统的数学分析对其进行了实用性分析，结果表明：

(1)以第三扁率n为参数的级数展开收敛速度快，更加适合级数展开，利用第三扁率n重新推导的子午线弧长公式各系数中的分母值都明显变小，个别高阶项系数消失，结构和形式简单。

(2)对子午线弧长正解公式的3种表示形式进行截断误差分析，发现子午线弧长正解公式8阶展开式能满足普通场景的使用情况，而10阶展开式截断误差精度提升明显，三角函数的倍角形式和指数形式精度提高了3个数量级，三角函数的二倍角形式精度提高了2个数量级，更适用于高精度场景。

(3)对子午线弧长正解公式8阶、10阶展开式进行了实用性分析，发现在0—30°低纬度地区，三角函数的指数形式子午线弧长公式比倍角形式误差几乎小了3倍以上，适用性更强；在30°—55°中纬度地区，三角函数的二倍角形式比倍角形式误差也几乎小了3倍以上，适用性更强；而在55°—90°高纬度地区，三角函数的倍角形式比二倍角形式误差几乎小了6倍以上，适用性更强。

初审：张艳玲

复审：宋启凡

终审：金   君