基于改进椭球模型的航线算法＊

2014-07-11李方能李厚朴吴延坤

舰船电子工程 2014年7期

李方能李厚朴吴延坤

（1.海军工程大学导航工程系武汉 430033）（2.中石化地球物理公司中原分公司濮阳 257001）

1 引言

传统的航路结构是基于连接地基导航设备限定的各个固定航路点的航线纵横交错而成的，通常按照导航台—导航台的原则进行设计，每个航段的航线角和航程根据等角航线的模型进行标注，飞机只能在以导航台为中心的辐射线或弧线上向台或背台飞行。20世纪70年代后期，随着全球航空业的飞速发展，空中流量急剧增加，传统的航路结构难以满足航班量增加的要求，延误和空域拥挤问题已成为航空界关注的焦点；同时，随着全球导航卫星系统等远程导航设备的发展，导航精度越来越高，传统的航路结构就空域使用和机载设备能力的应用而言显得越来越缺乏经济性、灵活性和高效性［1］。为此，国际民航组织于1991年提出了区域导航（Area Navigation，RNAV）的概念。区域导航是一种导航方法，允许飞机在地基导航设备的基准台覆盖范围内或自主导航设备能力限度内，或两者配合下按任何希望的路径飞行［1～2］。区域导航的优势是可以建立路程更短的径直航线以节约飞行成本，实现任意两点间的直线飞行，建立平行或双线航路以增加交通流量，提高空域资源利用率，缩短飞行间隔，进一步提高空中交通管制的灵活性［2～3］。区域导航是目前国际航空界的研究热点，代表未来导航技术的发展方向。

现代客机普遍装载有惯性导航系统（INS）和全球导航卫星系统（GNSS）等远程导航设备，可以实现国际民航组织倡导的区域导航逐点飞行［4］。目前的机载区域导航计算机在导航计算时采用的是大圆航线模型，该模型是建立在理想地球正球体基础上的。然而地球并不是一个标准的球体，而是近似于旋转椭球体，以大圆航线为模型的计算结果和真实数据存在一定误差，这对飞机的正常飞行和安全具有极大的影响。已有学者注意到了这一问题并进行了研究。周其焕等指出应当选择考虑地球扁率影响的大椭圆航线模型进行航迹计算［5］；方学东等研究了基于大椭圆航线模型的区域导航航路规划问题，推导出了初始航线角及航程的计算公式［6～7］。以往文献在推导航线角计算式时没有利用旋转椭球体的对称性，给出的表达式复杂繁琐；导出的航程计算公式在理论上不严密，若直接用于机载区域导航计算机软件，必将产生偏差；并且忽略了飞机飞行高度对航线角和航程计算的影响。由于当今的民航客机飞行高度最大可达15km，而在军事领域中日益受到关注的高空长航时无人机的飞行高度一般在20km以上［8］，因此为使机载区域导航计算机计算出更为准确的航迹，引导飞机沿航线精确飞行，建立考虑飞行高度的航线模型并给出正确的航线要素计算公式是十分必要的。鉴于此，本文对考虑飞行高度的航线算法问题进行了深入研究，推导出了初始航线角的直接计算公式以及航程的严密计算式，从而为机载区域导航计算机确定航线要素的提供了理论基础。

2 地球模型及导航坐标系

地球的自然表面高低起伏不平，是一个非常复杂而又不规则的曲面，无数学规律而言。为便于计算，需要选定一个与地球形状极为相近的、可用简单数学公式表示的、且能确定其与地球相关位置的表面作为基准面。研究表明，地球的赤道半径比短半径约大21km，地球的真实形状接近于南北稍扁的旋转椭球。

表1 若干参考椭球的相关参数

参考椭球是与大地体吻合最好的旋转椭球，它代表地球的数学模型，可由长半径a和扁率α来描述［9～10］。17世纪以来，许多国家的学者和机构根据不同地区、不同年代的观测资料推算出了若干地球参考椭球，具体如表1所示。

在航路规划、制作导航数据库及引导飞机飞行的过程中，存在两种坐标系，一种是以参考椭球的中心为坐标原点的空间直角坐标系，另一种是大地坐标系。空间直角坐标系中某点位置可用该点在坐标系各个轴上的投影（X，Y，Z）来表示；大地坐标系则采用大地纬度B、大地经度L和大地高程H来描述空间位置。大地坐标转换至空间直角坐标的数学关系为［9～10］

3 考虑飞行高度的改进椭球模型

图1 子午面直角坐标系

如图1所示，P点表示飞机当前所在的航路点，其大地坐标（BP，LP，HP）可由机载远程导航设备实测得到，由式（1）可得其空间直角坐标（XP，YP，ZP），图1中VP＝。在保持地球参考椭球中心O和短半轴方向不变的情况下，作一新椭球，使新椭球通过P点，并保持P在新椭球和原椭球中的大地经纬度不变，该新椭球即本文建立的考虑飞行高度的改进椭球。设新椭球的长半轴为an，第一偏心率为en，则由大地测量学可知［10］：

式中OK表示过P点的新、旧椭球的法线PK与短轴的交点K至椭球中心O的距离，Nn、N分别表示新、旧椭球中P点处的卯酉圈曲率半径，则有：

由式（2）并顾及式（3）可得：

至此，新椭球的第一偏心率en和长半轴an完全确定。

4 大椭圆航线要素计算

建立新椭球后，此时飞机航迹为新椭球上的大椭圆航线，即通过航路起点、终点和椭球中心所作的平面切割椭球得到的一个截面椭圆。飞机沿该航线航行时，必须随时确定飞机的大地坐标，并结合航路点大地坐标计算出初始航线角和航程。

4.1 大椭圆航线的初始航线角

图2 大椭圆航线示意图

设航行起点的大地坐标为P1（B1，L1，H1），可按上节方法建立考虑高度H1的新椭球，如图2所示，新椭球的长半轴an和第一偏心率en可按上节给出的公式计算，只是需将（BP，LP，HP）替换为（B1，L1，H1）。航行起点P1在新椭球上的大地坐标变为（B1，L1，0），设航行终点P2在新椭球上的大地坐标为（B2，L2，0），经纬度坐标可由导航数据库提取。顾及椭球的对称性，可设P1（B1，0，0）和P2（B2，ΔL，0），ΔL＝L2－L1，并对可能出现ΔL＞180°或ΔL＜－180°的情况进行规范化处理，即如果ΔL＞180°，则ΔL＝ΔL－360°，直至ΔL≤180°；如果ΔL＜－180°，则ΔL＝ΔL＋360°，直至ΔL≥－180°。P0QO是通过P1、P2两点和椭球中心O的截面椭圆，易知该截面椭圆的横轴必与新椭球的赤道圆的直径相重合，因而截面椭圆的长半轴等于该新椭球的长半轴an。Q点为该截面椭圆与椭球面某一子午线正交处，为截面椭圆的极，该点为截面椭圆上纬度最高处，OQ的长度即截面椭圆的短半径¯b。初始航线角A是指P1所在子午线北向顺时针旋转至椭圆P0QO在该点处的切线（指向P2）所经过的角度，变化范围为（0°，360°）。欲求初始航线角，首先需要确定子午圈所在平面NP1O与截面椭圆P0QO的夹角A1。令式（1）中H＝0，可得新椭球中大地坐标转换为空间直角坐标的公式为

上式展开后得到：dX＋fY＋gZ＝0，式中d＝－Y2Z1，f＝X2Z1－X1Z2，g＝X1Y2。通过P1、N、O的空间平面方程为：Y＝0。平面NP1O与截面椭圆P0QO的夹角A1亦即两平面法线矢量间的夹角，由向量点积可得：

表2 初始航线角A与A1的关系

4.2 大椭圆航线的航程

如图2所示，大椭圆航线的航程S即截面椭圆中劣弧P1P2的长度。为求航程，首先需要确定截面椭圆P0QO的短半径¯b和偏心率¯e。

该截面椭圆的极点Q处的地心纬度Φ0即赤道平面与截面椭圆P1P2O所在平面的夹角。考虑到赤道面方程为Z＝0，由两平面法线向量点积可得：

子午椭圆NOE的直角坐标方程式为

式中bn为新椭球的短半轴。直角坐标V、Z与极坐标ρ、Φ之间的关系为

将式（11）代入式（10），可得：

将由式（9）确定的Φ0代入上式，即得截面椭圆P0QO的短半径：

至此，可确定出该截面椭圆的偏心率：

该截面椭圆与赤道面的交线P0P′0的方程为

即：dX＋fY＝0，该直线与赤道圆的交点坐标可由式（16）确定：

易知这两个点的坐标分别为

由向量点积可得：

整理后得到：

根据文献［6］的说明，该式中的变量φ即截面椭圆P0QO中的夹角变量。事实上，利用该公式来计算航程从理论上讲是不严密的。因为在截面P0QO中，夹角变量实际上是该椭圆的地心纬度变量，而式（22）中的自变量φ则是该椭圆的归化纬度变量（此时截面椭圆的参数方程为：x＝ancosφ，y＝¯bsinφ），地心纬度与归化纬度之间的关系详见文献［9］和文献［10］，因此利用式（22）计算航程显然是不严密的。本节将给出严密的计算公式。

根据式（12）的有关推导，截面椭圆P0QO的极径¯ρ可以表示为

式中φ为该椭圆的地心纬度变量。该椭圆的参数方程为

在截面椭圆P0QO中，从赤道至φ处的弧长S（φ）可由对弧长的曲线积分表达式得到：

式（25）为椭圆积分，无法直接积出。注意到¯e很小，因此可将被积函数展开为¯e的幂级数形式，之后再进行积分，以上过程可利用计算机代数系统Mathematica［11］完成。在Mathematica中输入以下命令：

可以得到：

以上结果可进一步整理为

式中系数为

将θ1、θ2代入式（26），可依次得到自P0至P1、P2处的弧长S（θ1）、S（θ2）。至此，大椭圆航线的航程S为

考虑到在某些情况下θ1有可能大于θ2，S（θ2）与S（θ1）的差为负，因此这里取S（θ2）与S（θ1）差的绝对值作为航程。

特别，如果飞机的初始航线角A＝90°或A＝270°，即飞机沿等纬圈航行时，式（6）及式（26）将不再适用，此时航程S的计算公式应为

5 算例分析

设飞机从北京P1（40°N，116°E）起飞，飞行高度约为10km，终点为底特律P2（43°N，83°W），选用表1中的WGS84椭球参数，利用本文给出的有关公式确定其初始航线角和航程。

1）初始航线角A的计算

WGS84椭球的第一偏心率为

由式（3）知参考椭球中起点P1处的卯酉圈曲率半径为：N1＝6386976.172，新椭球中该点处的卯酉圈曲率半径为：Nn＝6386976.172＋10000＝6396976.172，由式（4）和式（5）可得新椭球的第一偏心率和长半径：en＝0.081755244960998，an＝6388137.010。根据新椭球的对称性，可设P1（40°N，0°E）和P2（43°N，ΔL），ΔL＝－83°－116°＋360°＝161°E。根据式（6）可得新椭球中航行起点和终点的空间直角坐标分别为

由式（7）可得：d＝－6.22228×1012，f＝－3.93106×1013，g＝7.46532×1012。根据式（8）可得：A1＝13.8863°，由表2知，初始航线角A＝A1。

2）航程S的计算

由式（13）和式（14）可得截面椭圆P0QO的短半径和第一偏心率为：¯b＝6367475.599，¯e＝0.080363052034348；根据式（20）和式（21），可得：θ1＝0.709457，θ2＝2.377986分别代入式（25）得到：S（θ1）＝4.52989×106S（θ2）＝1.51612×107，因此航程为

作者对不考虑飞机飞行高度的情况进行了计算，得到此时初始航线角为 13.8864°，航程为1.06147×107m，与考虑高度时的航程相差16600m，可见飞行高度对航程的影响比较明显。将此时的θ1、θ2代入文献［6］和文献［7］给出的式（22），得到的航程为1.06352×107m，比考虑高程时的航程还要大，显然是不合理的，究其原因是该式将归化纬度和地心纬度的概念混淆了。