基于改進橢球模型的航線算法＊

李方能 李厚樸 吳延坤

（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，比考慮高程時的航程還要大，顯然是不合理的，究其原因是該式將歸化緯度和地心緯度的概念混淆了。

6 結語

目前，機載區域導航計算機在計算航跡時采用的是大圓航線模型，為滿足遠距離高精度導航的需要，本文建立了顧及飛機飛行高度的改進橢球模型，系統地討論了大橢圓航線要素的計算問題。研究表明：

1）利用改進橢球模型的對稱性，推導出了計算大橢圓航線初始航線角的直接公式，給出了象限判斷原則，與以往文獻中的公式相比，該式形式簡潔明了，便于應用。

2）前人給出的航程計算公式在理論上是不嚴密的，本文采用空間向量分析方法，借助 Mathematica代數系統，推導出了嚴密的計算公式，該式系數為關于大橢圓偏心率的符號形式，可解決不同參考橢球下航程的計算問題。算例分析表明在進行遠距離航行時，飛行高度對航程的影響比較明顯，需要加以考慮。

3）采用本文給出的計算方法改進飛機的區域導航計算軟件后，可以得到更為準確可靠的初始航線角和航程，從而可以減小飛機的顯示誤差，降低飛行員的心理負荷，進一步提高飛機的安全性和經濟性。

［1］胡慧玲.RNP和RNAV概念淺析［J］.空中交通管理，2007（3）：24－25.

［2］馬志剛，王宏建.必備導航性能／區域導航［J］.中國民航飛行學院學報，2001，12（2）：33－36.

［3］劉渡輝，帥斌，王大海，等.區域導航航路和進離場程序設計分析與仿真［J］.交通運輸工程與信息學報，2006，4（3）：123－127.

［4］方學東.基于大圓航線的RNAV航路規劃［J］.中國民航飛行學院學報，2006，17（2）：3－5.

［5］周其煥，陸學軍.FANS航線和最短航線算法問題［J］.民航經濟與技術，1997（188）：51－55.

［6］方學東，黃潤秋，魏光興.基于橢球模型的FANS航線算法［J］.西南交通大學學報，2006，41（4）：512－516.

［7］方學東.考慮地球扁率的RNAV航路規劃［J］.測繪科學，2008，33（6）：149－150.

［8］王新龍，馬閃.高空長航時無人機高精度自主定位方法［J］.航空學報，2008，29（增刊）：39－45.

［9］邊少鋒，柴洪州，金際航.大地坐標系與大地基準［M］.北京：國防工業出版社，2005：7－34.

［10］孔祥元，郭際明，劉宗泉.大地測量學基礎［M］.武漢：武漢大學出版社，2005：97－105.

［11］邊少鋒，許江寧.計算機代數系統與大地測量數學分析［M］.北京：國防工業出版社，2004：10－53.