c++高斯投影正反算_论文推荐 | 李松林：常用海图投影平面上大椭圆航线的表象与曲率分析...

《测绘学报》

构建与学术的桥梁 拉近与权威的距离

常用海图投影平面上大椭圆航线的表象与曲率分析

李松林, 陈成, 边少锋, 李厚朴, 刘强

海军工程大学导航工程系, 湖北 武汉 430033

收稿日期：2018-07-19；修回日期：2019-06-03

基金项目：国家自然科学基金(Nos.41631072；41571441；41604010)

第一作者简介：李松林(1991-), 男, 博士生, 研究方向为地图投影。E-mail:songlin829@126.com

通信作者：陈成, E-mail: chenchengzyq@163.com

摘要：利用空间矢量方法推导出了椭球面上只与起止点地理坐标有关的大椭圆航线方程，代入4种常用海图投影的正解公式，得到不同投影平面上的大椭圆参数方程；利用上述参数方程进而推导出了不同投影面上大椭圆航线的曲率与曲率半径公式。选取伦敦到纽约的大椭圆航线为例，通过绘制不同投影面上的大椭圆航线并分析其曲率、曲率半径变化曲线可知，大椭圆航线在日晷投影上的表象为曲率处处为0的直线，而在其他投影面上的表象为曲率较小但不断变化的曲线。利用推导的曲率半径公式可以计算各类大椭圆航线上任意位置的“代曲直距”，方便在不同比例尺的海图上对大椭圆航线进行量测和绘制，提高作图效率。

关键词：大椭圆航线 海图投影 曲率 曲率半径 代曲直距

Representation and curvature analysis of great ellipse on common chart projection plane

LI Songlin, CHEN Cheng, BIAN Shaofeng, LI Houpu, LIU Qiang

Department of Navigation, Naval University of Engineering, Wuhan 430033, China

Foundation support: The National Natural Science Foundation of China(Nos. 41631072; 41571441; 41604010)

First author: LI Songlin(1991—), male, PhD candidate, majors in map projection. E-mail: songlin829@126.com.

Corresponding author: CHEN Cheng, E-mail: chenchengzyq@163.com.

Abstract: By using the space vector method, the equation of great ellipse on ellipsoidal surface only related to the geographical coordinates of the starting and ending points was derived, the parameter equations of great ellipse routes on various projection planes were obtained from the great ellipse equation and the positive solution formulas of the four kinds of common projection. And then the curvature and radius of curvature of great ellipse routes on the four kinds of projection planes were derived. The great ellipse route from London to New York was taken as an example, by drawing the great ellipse route on different projection planes and analyzing the curve of curvature and curvature radius, the conclusion is drawn that the representation of great ellipse route on gnomonic projection plane is straight line, while which of great ellipse route on the other projection planes are curves with curvature changing slightly. The formulas of curvature radius derived in this paper can be used to calculate the "substitution distance" of great ellipse routes, which is convenient to measure and draw the great ellipse routes on the nautical chart with any scale, and improve the efficiency of nautical drawing.

Key words: great ellipse route nautical chart projection curvature radius of curvature substitution distance

测地线是指曲面上测地曲率处处为零的曲线，也是曲面上任意两点之间距离最短的曲线^[1]，因其具有短程性而广泛用于大地测量、航海导航等领域。将地球看成球体的情况下，参考球面上的大圆线即为测地线，相关参数可以方便求出；而地球实际上是一个近似的旋转椭球体，针对地球椭球面上的测地线方程或弧长却很难计算^[2]，又由于地球椭球面上的大椭圆与大地线在距离上非常接近，15 000 km的远程距离只有几米之差，因此在航海导航应用过程中，一般用大椭圆线代替测地线作为最短航线，以实现经济高效航行^[3]。

针对大椭圆航线及其相关应用问题，国内外学者进行了深入的研究并取得丰硕的成果^[4-21]。文献[4—8]分别推导了大椭圆地理坐标解算公式，探讨了大椭圆航线与等角航线的关系并比较了二者在远洋航行中优劣，分析了大椭圆航法与大圆航法的航程计算误差，并结合地图投影理论，研究了大椭圆表象为直线的日晷投影以及墨卡托海图上大地线曲率问题；文献[9—12]采用位置矢量的方法，研究了大椭圆航线航程、方位角和正反解问题；文献[13—15]则采用单位速度矢量的方法，研究了大椭圆航线航程、方位角和正反解问题；文献[16—17]借助计算机代数系统Mathematica，研究了大椭圆航行方法及其相关的导航参数计算，并导出了大椭圆航法航程符号形式的计算公式；文献[18—19]研究了基于大椭圆航线的RNAV航路规划问题；文献[20]在总结前人成果的基础上，采用空间矢量代数方法直接求解大椭圆顶点，提出了基于Newton-Raphson(N-R)的等距离航线设计算法。

总结上述研究成果，大多围绕大椭圆航线的地理坐标、航程和方位角等参数的解算，基于大椭圆航线的航路规划与算法设计等问题展开，针对不同类型投影平面下大椭圆航线的描写，及其曲率和代曲直距问题缺乏系统的研究。而在航线规划与航海导航应用中，大椭圆航线的曲率和代曲直距问题关乎作图的精确性、便捷性和航行的安全性、经济性，因此本文分别针对4种常用投影(墨卡托投影、高斯投影、极球面投影和日晷投影)平面下大椭圆航线的表象、曲率和代曲直距进行分析研究，以解决航线绘制过程中的“以直代曲”问题，为航海人员进行大椭圆航线的量测、设计和绘制提供理论与计算分析依据。

1 常用投影平面上的大椭圆航线参数方程

以旋转椭球体球心为原点，建立空间直角坐标系如图 1所示。i、j、k为X、Y、Z坐标轴单位向量。过椭球面上两点P₁(B₁,l₁)、P₂(B₂,l₂)的球心截面OP₁P₂与椭球面的交线上两点之间的劣弧即为大椭圆航线，记为P₁P₂。根据文献[20]中的推导，截面椭圆OP₁P₂所在平面的法向量n的矢量表达式为

(1)

图 1 大椭圆航线示意图Fig. 1 Sketch of great ellipse route

图选项

式中

(2)

设点P₁(B₁,l₁)到P₂(B₂,l₂)间的大椭圆航线上的任意一点为P(B,l)，根据大地坐标与空间直角坐标的关系可得P点的坐标矢量P=NcosBcosl·i+NcosBsinl·j+N(1-e²)sinB·k。由于P是截面椭圆OP₁P₂与椭球面交线上的点，因此向量n与P垂直，即n·P=0，据此可得P点坐标方程，亦即大椭圆P₁P₂的方程为

(3)

即

(4)

等价于

(5)

式(4)说明，大椭圆航线上任意一点大地坐标与椭球参数(长短半轴、偏心率)无关，只与起止点大地坐标有关。将式(5)大椭圆航线方程写成B(l)=

，并分别代入椭球墨卡托投影、高斯投影、极球面投影和日晷投影4种常用投影方程中，可以得到4种常用投影平面上的大椭圆线的参数方程。

其中，墨卡托投影平面上的大椭圆线参数方程为^[22-23]

(6)

式中，a为椭球长半轴；q为等量纬度，其表达式为

(7)

高斯投影平面上大椭圆线参数方程为^[22-25]

(8)

式中，系数αi见文献[23]；u、v为中间变量，表达式为

(9)

极球面投影平面上大椭圆线参数方程为^[22-23]

(10)

式中，Rφ为等角球半径，表达式为Rφ=

。

日晷投影平面上大椭圆线参数方程为^[22-23]

(11)

式中，B₀为椭球面与地心纬度球面切点处纬度，相应的地心纬度为φ₀，表达式为φ₀=arctan[(1-e²)tanB₀]；Rg为切点处的椭球动径，表达式为Rg=

。为便于比较分析，上述引用的投影公式比例系数均取为1。

2 常用投影平面上大椭圆航线的曲率与曲率半径公式推导

由大椭圆航线方程式(5)，求得大地纬度B对经差l的一阶、二阶导数为

(12)

(13)

由式(7)求得等量纬度q对大地纬度B的一阶导数为

(14)

进而得到等量纬度q对经差l的一阶、二阶导数为

(15)

(16)

2.1 墨卡托投影平面上大椭圆航线曲率公式

联系椭球面墨卡托投影公式(6)，可得大椭圆航线平面坐标的一阶、二阶导数

(17)

(18)

最终得到大椭圆航线在墨卡托投影平面上的曲率公式为

(19)

曲率半径为M₁=

。与经线重合的大椭圆航线在墨卡托投影平面上表现为直线，曲率为0，曲率半径为无穷大。

2.2 高斯投影平面上大椭圆航线曲率公式

由式(9)首先得到中间变量u、v对经差l的一阶、二阶导数为

(20)

(21)

然后由式(8)得到投影平面坐标对经差l的一阶、二阶导数为

(22)

(23)

最终求得曲率κ₂=M₂=

大椭圆航线与经线重合时，大椭圆航线方程与l无关，此时航线的表象就是经线的表象，可令经度为定值，以纬度为曲线参数进行类似的计算。设此时的大椭圆航线l≡lGE，以等量纬度为曲线参数，由式(9)得中间变量u、v对等量纬度q的一阶导数和二阶导数为

(24)

(25)

然后只需把式(20)—式(23)中的经度参数改为等量纬度参数，即以

替换

、

，即可计算大椭圆航线的曲率等参数。

2.3 极球面投影平面上大椭圆航线曲率公式

联系椭球面极球面投影方程式(10)，可得一阶、二阶导数为

(26)

(27)

最终得到大椭圆航线在极球面投影平面上的曲率公式为(极球面投影平面y为纵坐标，x为横坐标)

(28)

曲率半径为M₃=

2.4 日晷投影平面上大椭圆航线曲率公式

联系椭球面日晷投影公式(11)，得到大椭圆航线上任意一点P(B,l)与起始点P₁(B₁,l₁)的日晷投影坐标差为

(29)

由式(5)有

(30)

因此，有

(31)

由式(31)可知，大椭圆在日晷投影平面上表象为一直线，曲率为0，曲率半径为无穷大。与经线重合的大椭圆航线在日晷投影平面上显然也表现为直线，曲率也为0。

3 算例分析

选取纽约(40°43′N, 74°00′W)至伦敦(51°30′N, 00°05′E)的大椭圆航线为例，取WGS 84参考椭球参数，a=6 378 137 m，e=0.081 819 190 843，根据曲率与曲率半径公式，分别绘出墨卡托投影、高斯投影、极球面投影和日晷投影平面上的曲率随大椭圆航线经度变化的曲线如图 2所示，曲率半径变化曲线如图 3所示。

图 2 常用投影平面上大椭圆航线曲率变化曲线Fig. 2 Curves of curvature of great ellipse routes on common projection planes

图选项

图 3 常用投影平面上大椭圆航线曲率半径变化曲线Fig. 3 Curves of curvature radius of great ellipse routes on common projection planes

图选项

结合图 2、图 3分析：随着经度变化，日晷投影平面上的大椭圆航线曲率恒为零，说明该投影平面上的大椭圆航线确实为一条直线，而其他3种投影平面上大椭圆航线曲率不为零，说明其投影平面上大椭圆航线的表象为曲线。其中，极球面和高斯投影平面上大椭圆航线的曲率数量级为10^-8，而墨卡托投影平面上数量级为10^-7，说明墨卡托投影平面上大椭圆航线弯曲更明显，这是由于墨卡托投影在中高纬度地区的变形相对较大导致的。同时还可以看出，极球面投影平面上大椭圆航线的曲率变化非常小(看起来为一条直线，但将坐标尺度变小后发现仍为一条曲线)，说明大椭圆航线在该投影平面上的表象近似为一圆弧(半径足够大)。综上，在日晷投影平面上进行大椭圆航线的量测和绘算最为方便，而在墨卡托投影平面上针对中高纬度地区大椭圆航线的量测误差较大。

借助计算机代数系统Mathematica，采用GeoGraphics绘图命令，绘出4种常用投影平面上纽约至伦敦的大椭圆航线如图 4所示(其中图 4(b)用横墨卡托投影来代替高斯投影)。由图 4可看出大椭圆航线在日晷投影平面上为一条直线，在墨卡托投影平面上为一条弯曲较为明显的曲线，这充分验证了上述分析结论的准确性，也验证了前文推导公式的正确性。

图 4 常用投影平面上的大椭圆航线Fig. 4 Great ellipse route on common projection plane

图选项

通过上述分析可知，大椭圆航线在海图上绘制时，一般都为曲线，难以直接精确绘制，因此可将海图上的大椭圆航线进行分段处理，相邻两点用直线连接，“以直代曲”，绘制精度满足《中国航海图编绘规范》中规定的位置偏差(“以直代曲”时的拱高)不大于制图允许误差Δ=0.01 cm即可。根据文献[8, 26]，海图投影平面上一般曲线的代曲直距满足如下不等式

(32)

式中，k=89.44；S以厘米为单位(图上距离)，若以实地距离计算，为SC₀/100=0.894 4C₀为比例尺分母，M为曲率半径(以千米为单位)。

以墨卡托投影平面上纽约至伦敦的大椭圆航线为例，通过计算得到在经度l=-23°1.5′处，曲率绝对值达到极大值点，相应的曲率为κ=-1.262 4×10^-7m^-1，在该处曲率半径达到极小值点，其值为M=7.921 4×10⁶m。因此，在1：500 000墨卡托航海图上，大椭圆航线代曲直距为

。同理，在其他不同投影平面上，可以计算得到大椭圆航线的曲率、曲率半径随经度变化的极值点，从而在选定的比例尺下计算最大代曲直距，结果见表 1。同样的，可以利用推导出的曲率半径公式计算大椭圆航线上任意位置点的代曲直距。

表 1 常用投影平面上大椭圆航线的最大代曲直距Tab. 1 Maximum of substitution distance of great ellipse route on common projection plane

表选项

由于日晷投影平面上大椭圆航线为曲率处处为0的直线，因此代曲直距可取任意值，这是日晷投影最重要的优点，为航线量测提供了方便。

4 结论

本文围绕墨卡托投影、高斯投影、极球面投影和日晷投影4种常用海图投影平面上的大椭圆航线的表象和曲率进行研究，推导出了不同投影平面上的大椭圆航线参数方程、曲率和曲率半径公式，并根据上述方程，绘制了不同投影面上纽约至伦敦大椭圆航线的曲率和曲率半径变化曲线。结果表明，大椭圆航线在日晷投影平面上的表象为一条直线，在其他3种投影面上为曲率较小但曲率不断变化的曲线，其中相比高斯投影和极球面投影，墨卡托投影平面上的大椭圆航线曲率和曲率变化率较大。因大椭圆航线的经济性及其在日晷投影海图上量测和绘算的便利性，使得日晷投影在舰船远洋航行中具有重要应用价值。但由于大椭圆航线在其他常用海图上绘制时一般都为曲率不断变化的曲线，难以直接精确绘制，因此一般将海图上的大椭圆航线进行分段处理，“以直代曲”。利用本文推导的公式可以方便地计算相应航线的曲率半径，进而计算大椭圆航线上任意位置的代曲直距，以解决航线绘制过程中的“以直代曲”问题。本文的研究工作可为航海人员进行大椭圆航线的量测、设计和绘制提供便利参考。

【引文格式】李松林, 陈成, 边少锋, 等. 常用海图投影平面上大椭圆航线的表象与曲率分析. 测绘学报，2019，48(10)：1331-1338. DOI: 10.11947/j.AGCS.2019.20180348

权威 | 专业 | 学术 | 前沿微信投稿邮箱 | song_qi_fan@163.com

欢迎加入《测绘学报》作者QQ群： 297834524

进群请备注：姓名+单位+稿件编号