民用航空慣性基準系統(tǒng)重力補償方法

王鑫耀,李榮冰,劉 剛,劉建業(yè)

（1.南京航空航天大學導航研究中心,南京 211106;2.航空工業(yè)西安飛行自動控制研究所,西安 710065)

0 引言

民用航空在國家發(fā)展中具有重要的戰(zhàn)略地位,為促進民用航空的發(fā)展,國際民航組織（International Civil Aviation Organization,ICAO)提出了所需導航性能（Required Navigation Performance,RNP)的概念。在飛行的各個階段,RNP運行標準用RNP x表示,x為所需滿足的精度性能指標,表示總體系統(tǒng)誤差,參考航空無線電技術(shù)委員會（Radio Technical Commission for Aeronautics,RTCA)DO-236標準中的精度性能指標分配方法,導航系統(tǒng)的定位精度應優(yōu)于0.3xnmile。RNP進近（Approach,APCH)是飛行難度最高的階段,所需精度一般為0.1nmile～0.3nmile。

慣性基準系統(tǒng)在民用航空中是位置、速度、姿態(tài)等信息的核心信息源,即使是在衛(wèi)星失效的情況下,也能在短期內(nèi)獨立完成導航工作并保持可靠的性能[1]。由于民用航空應用的安全關鍵性,導航系統(tǒng)必須經(jīng)過嚴格的測試,以確保其滿足RNP各個階段的指標要求。在導航系統(tǒng)的安全認證過程中,必須妥善處理所有的系統(tǒng)誤差,其中的重力模型誤差是需要關注的一項,尤其是在純慣性運行期間。

RTCA制定的全球定位系統(tǒng)/廣域增強系統(tǒng)機載設備的最低運行性能標準（DO-229)的附錄R（緊組合系統(tǒng)的要求和測試程序)是緊組合系統(tǒng)研制的行業(yè)標準。附錄R中指出,重力模型誤差會導致導航系統(tǒng)誤差,特別是在衛(wèi)星導航失效后的純慣性工作狀態(tài)下會產(chǎn)生問題。因此,附錄R要求導航系統(tǒng)必須應用某種形式的重力補償方法,使系統(tǒng)即使在衛(wèi)星失效的純慣性工作狀態(tài)下運行于重力模型誤差較大的地區(qū)時也能繼續(xù)約束導航誤差。重力補償能約束由重力模型不準確造成的導航誤差,但無法彌補慣導系統(tǒng)本身固有的誤差累積。

文獻[2]分析了重力誤差模型對慣性導航系統(tǒng)精度的影響,文獻[3]推導了水平姿態(tài)誤差與水平重力擾動間的傳遞函數(shù),但均未對重力擾動進行補償。文獻[4]假設重力擾動符合均勻分布,其重力補償模型與實際不符。文獻[5]使用重力網(wǎng)格數(shù)據(jù)進行插值補償,會受到區(qū)域限制。文獻[6]使用重力梯度儀實時補償重力擾動,但重力梯度儀的主要用途是在水下,空中的重力梯度變化量不足以支持慣導的補償。近年來,學者們研究了高階重力場球諧模型,文獻[7]、文獻[8]利用EGM2008直接計算重力擾動,但沒有涉及地球自轉(zhuǎn),并未計算真實重力,對重力概念的闡述也不夠完整。

EGM2008是美國國家地理空間情報局（National Geospatial-intelligence Agency,NGA)發(fā)布的一個地球引力場的球諧模型,由EGM2008所計算得到的在美國和澳大利亞上空的垂線偏差與獨立的天體測量值的差異在1.1″～1.3″[9]。基于 EGM2008模型,可以獲得高精度的地球引力數(shù)據(jù)。但由于高階重力場球諧模型非常復雜,實際應用中存在困難,現(xiàn)有的研究主要集中在重力擾動補償?shù)姆抡骝炞C分析階段。

針對現(xiàn)有研究中存在的重力補償模型不準確、重力概念闡述不完整的問題,本文闡述了基于高階重力場球諧模型EGM2008的完整重力矢量計算方法,分析了慣性基準系統(tǒng)在沒有衛(wèi)星輔助的純慣性運行期間由重力模型誤差引起的定位誤差及重力補償后的改進效果。針對高階重力場模型在應用場景約束下使用困難的問題,分析了重力模型階次與性能優(yōu)化的關系。

1 重力模型誤差對慣性基準系統(tǒng)的影響分析

慣性基準系統(tǒng)在導航解算過程中需要載體相對于導航系的運動加速度,而加速度計無法區(qū)分運動加速度和重力加速度,加速度計輸出的是比力值。由比力提取載體運動加速度的過程稱為比力變換,比力方程為

在傳統(tǒng)的重力模型下,重力矢量總是垂直于當?shù)厮矫?在水平面內(nèi)沒有分量,稱之為正常重力矢量[10]。在當?shù)氐乩碜鴺讼迪?正常重力矢量表示為

垂直分量γ的大小由Somigliana函數(shù)計算得到,其是緯度L的函數(shù)

式（3)中,e為橢球偏心率。

由于地球內(nèi)部質(zhì)量分布不均勻,真實重力矢量與正常重力矢量間存在重力模型誤差,稱為重力擾動。重力擾動分為垂線偏斜和重力值異常:垂線偏斜定義為重力矢量與垂直方向的偏差,重力值異常定義為重力矢量大小的幅度誤差。重力擾動不是隨機誤差,而是在局部地區(qū)的每次飛行中都會出現(xiàn)的相同的系統(tǒng)誤差。

慣性基準系統(tǒng)使用正常重力矢量代替真實重力矢量進行慣導解算時,會因為重力擾動而使運動加速度計算結(jié)果存在系統(tǒng)誤差并隨積分運算影響導航解算結(jié)果。速度誤差方程為

式（4)中,Φn為姿態(tài)誤差向量的反對稱矩陣,fn和Δn分別為加速度計輸出的比力值和零偏在導航系的投影,δg為重力擾動矢量。

相比于重力值異常,垂線偏斜是影響慣性基準系統(tǒng)性能的主要因素,全球垂線偏斜最大在100″數(shù)量級,而1″的垂線偏斜相當于5μg的水平重力分量。在慣性器件精度不高的系統(tǒng)中,重力擾動相對于器件誤差可以忽略。而隨著慣性器件精度的提升,重力擾動與器件誤差處在同一數(shù)量級甚至高于器件誤差,這種情況下的重力擾動不可忽略,有必要采用更高精度的重力模型。重力補償下的慣性基準系統(tǒng)框圖如圖1所示。

圖1 重力補償下的慣性基準系統(tǒng)原理圖Fig.1 Schematic diagram of inertial reference system under gravity compensation

重力補償?shù)男畔⒘鞒虨?慣性基準系統(tǒng)向重力場模型提供位置信息,根據(jù)位置信息計算真實重力矢量,利用該重力矢量進行比力變換,解算后得到導航信息。需要注意的是,由慣性器件等誤差源造成的位置偏差也會引起重力誤差,但重力場變化比較緩慢,一般采用網(wǎng)格補償?shù)男问?所以由位置偏差引起的重力誤差相比重力擾動可忽略,本文研究的重力補償以重力擾動為主。

2 基于EGM2008模型的重力計算方法

地球引力與離心力的合力稱為重力,重力矢量可以通過計算重力勢能的梯度獲得。重力勢能由引力勢能和離心勢能組成

球坐標系下的引力勢能可表示為

式（6)中,G為萬有引力常數(shù),M為地球質(zhì)量,r為計算點到球心的距離,θ為極距或稱為球心余緯度,λ為經(jīng)度,a為參考橢球長半軸,Cnm和Snm為球諧模型的系數(shù),n為球諧模型的階,m為球諧模型的次,Pnm（cosθ)為 n階 m 次 Legendre函數(shù)項。

球坐標系下的離心勢能可表示為

式（7)中,ω為地球自轉(zhuǎn)角速度。

計算引力勢能的梯度

計算離心勢能的梯度

將計算得到的引力勢能梯度與離心勢能梯度相加得到重力矢量,根據(jù)球坐標（λ,θ,r)與球面直角坐標（x,y,z)的轉(zhuǎn)換關系

將重力矢量從球坐標系轉(zhuǎn)換到球面直角坐標系下

導航解算通常在當?shù)氐乩碜鴺讼迪峦瓿?使用的是大地緯度,而球面直角坐標系使用的是地心余緯度,因此需要將球面直角坐標系下的重力矢量轉(zhuǎn)換到當?shù)氐乩碜鴺讼迪?坐標系的轉(zhuǎn)換關系如圖2所示。將球面直角坐標系的y軸反向為y′, （x,y′,z)繞 x（xn)軸轉(zhuǎn)過 - α 角即可轉(zhuǎn)換為當?shù)氐乩碜鴺讼担▁n,yn,zn)。

圖2 坐標系轉(zhuǎn)換示意圖Fig.2 Schematic diagram of coordinate system conversion

旋轉(zhuǎn)矩陣為

綜上,當?shù)氐乩碜鴺讼迪碌闹亓κ噶繛?/p>

式（14)中,R（- β)=diag（1,-1,1)。

需要注意的是,基于球諧模型的重力計算中用到的是完全正常化的締合Legendre函數(shù),通過遞推計算得到。完全正常化的締合Legendre函數(shù)及其導數(shù)的遞推關系見文獻[11]。

為驗證本文利用球諧函數(shù)計算重力矢量方法的有效性,對全球的真實重力矢量進行了計算,采用WGS-84參考橢球模型,網(wǎng)格分辨率為1°×1°。

計算得到的東向、北向、天向的真實重力矢量分量如圖3～圖5所示。其中,東向和北向的重力分量較小,采用毫伽（mGal)為單位（1mGal=1×10-5m/s2)。

圖3 東向重力分量Fig.3 Diagram of eastward gravity component

圖4 北向重力分量Fig.4 Diagram of northward gravity component

將計算得到的重力矢量大小與國際全球地球模型中心（International Centre for Global Earth Models,ICGEM)提供的重力值進行對比,得到的殘差均值為-1.3×10-5mGal,標準差為6.1×10-2mGal,驗證了本文重力矢量計算方法的有效性。

本章節(jié)完整闡述了基于EGM2008球諧模型的重力矢量計算方法,與文獻[7]、文獻[8]中的計算方法相比,本文方法補充了對離心加速度和坐標系轉(zhuǎn)換的說明,使重力概念更完整,計算得到的重力矢量可直接用于地理系內(nèi)的導航解算,不需要再增加額外的計算。

3 純慣性運行狀態(tài)重力補償效果的仿真驗證

民用航空對安全性有極高的要求,由于衛(wèi)星無線電信號容易受到干擾,在衛(wèi)星信號不可用的情況下,慣性基準系統(tǒng)應當能夠獨立工作,在短時間內(nèi)繼續(xù)保持可靠的性能。在純慣性運行期間,重力模型的誤差不可忽略。

為驗證球諧重力模型對純慣性工作狀態(tài)下的慣性基準系統(tǒng)補償?shù)挠行?對某航跡下的導航系統(tǒng)性能進行數(shù)值仿真,仿真條件如下:三軸陀螺儀的常值漂移為0.005（°)/h,三軸加速度計的零偏 為 20μg, 航 跡 的 起 點 設 置 為 （30°N,0°E,5000m),飛行器沿30°N緯線以250m/s的速度向東飛行。通常情況下,民用航空要求導航系統(tǒng)有20min的純慣性運行性能,因此設置仿真時間為20min,飛行距離為162nmile。整個飛行過程中,航向不變無機動,飛機保持平飛,陀螺儀和加速度計的信息由數(shù)值解法給出。

根據(jù)本文的重力矢量計算方法計算飛行航跡上的重力擾動,重力擾動情況及統(tǒng)計信息如圖6和表1所示。

圖6 航跡重力擾動情況Fig.6 Diagram of trajectory gravity disturbance

表1 航跡重力擾動數(shù)據(jù)統(tǒng)計Table 1 Data statistics of trajectory gravity disturbance

在無重力補償?shù)那闆r下,設定航跡下的仿真位置誤差如圖7所示。

圖7 無重力補償?shù)奈恢谜`差Fig.7 Position errors with no gravity compensation

相同仿真條件下,采用5′×5′分辨率的重力網(wǎng)格數(shù)據(jù)對系統(tǒng)進行補償,每飛過一個網(wǎng)格重新計算一次重力矢量,補償后的位置誤差如圖8所示。

圖8 5′×5′重力補償?shù)奈恢谜`差Fig.8 Position errors with 5′×5′gravity compensation

從仿真結(jié)果可以直觀看出,位置誤差經(jīng)過重力補償后顯著減小,導航系統(tǒng)的精度得到有效提高。為進一步分析重力補償?shù)男Ч?對仿真結(jié)果進行統(tǒng)計,統(tǒng)計結(jié)果如表2所示。

表2 重力擾動補償前后的位置誤差Table 2 Position errors before and after gravity disturbance compensation

由表2可知,在民用航空要求的20min純慣性運行期間,針對特定航線和慣性器件精度水平,重力補償前的水平位置誤差為210.6m,該誤差由重力模型誤差和器件誤差共同引起;重力補償后的水平位置誤差為83.1m,該誤差主要由器件誤差引起。因此,重力模型誤差導致的定位誤差為127.5m,占到了水平定位誤差的61%,由此可見重力補償?shù)谋匾浴?/p>

重力模型誤差引起的定位誤差為確定性誤差而非隨機誤差,該誤差在多數(shù)航路階段是微不足道的,但在進近階段卻不能忽略。以RNP0.3的進近為例,分配給導航系統(tǒng)的定位精度應優(yōu)于0.09nmile（166.7m),該誤差占到了導航定位精度要求的76%,因此面對一定的RNP要求,重力模型誤差引起的定位誤差不可忽略,這限制了慣性基準系統(tǒng)在民用航空中支撐高RNP要求的能力。在采用5′×5′分辨率的重力網(wǎng)格數(shù)據(jù)補償后,導航定位精度顯著提升,有助于支撐慣性基準系統(tǒng)達到更高的RNP要求。

4 面向應用場景約束的重力模型階次與性能優(yōu)化分析

EGM2008模型的階數(shù)為2190,實時的重力計算將對實際應用場景下導航計算機的CPU性能和內(nèi)存占用帶來巨大的挑戰(zhàn)。在導航計算機性能不足的情況下,可以犧牲部分重力信息,對模型進行適當裁剪,以達到實時計算的要求。

為驗證重力模型裁剪的合理性,采用與本文上一章節(jié)相同的航跡,降低重力模型的階數(shù),分析使用裁剪后重力模型的補償效果。選擇的模型階數(shù)為1600、1000、600、360、180、90、20, 計算航跡上的水平重力擾動,以2190階的重力擾動為基準,計算降階后的重力擾動殘差,結(jié)果如圖9所示。

圖9 水平重力擾動的殘差對比Fig.9 Residuals comparison of horizontal gravity disturbance

由圖9可知,360階以下隨著階數(shù)的增加,重力擾動殘差有明顯收斂的趨勢,360階的殘差在5mGal以內(nèi),隨著階數(shù)繼續(xù)提高殘差收斂的趨勢并不明顯。

（1)重力場模型階次與補償效果關系分析

使用裁剪的重力模型補償導航系統(tǒng),忽略了器件誤差以避免器件誤差的隨機性影響效果對比,定位誤差的補償效果如表3所示。

表3 裁剪重力模型的補償效果對比Table 3 Comparison of compensation effect for cropped gravity model

由表3可知,360階及以上的重力模型對導航定位誤差的補償效果在90%以上。隨著模型階數(shù)的增加,補償效果并不能獲得線性提升,過高的階數(shù)對性能的提升有限,反而增加了導航計算機的負擔。

（2)重力場模型階次與計算機資源占用關系分析

為適應應用場景約束下的導航計算機的性能條件,給實際中的模型選擇提供參考,對重力模型的內(nèi)存占用和計算耗時進行分析,結(jié)果如表4所示。模型數(shù)據(jù)存儲在一維數(shù)組中,采用雙精度浮點數(shù)據(jù)類型,單次耗時取10次計算耗時的均值。本次計算在PC平臺下完成,使用的CPU型號為Intel i5-9300H,由于實際應用場景下的計算性能存在差異,本文計算的耗時僅為相對結(jié)果,可為實際應用場景下的重力模型集成提供參考。

表4 裁剪重力模型的內(nèi)存占用及耗時對比Table 4 Comparison of memory usage and time consumption for cropped gravity model

重力計算的時間和空間復雜度為O（n2),隨著階數(shù)的降低,內(nèi)存占用和耗時明顯減少。實際應用中,應根據(jù)導航計算機的內(nèi)存容量和CPU性能,結(jié)合精度要求綜合選用模型階數(shù)。

此外,由于重力場變化較為緩慢,采用5′×5′分辨率的重力網(wǎng)格數(shù)據(jù)足以取得較好的補償效果,因此不必在慣導解算的每個周期都進行重力矢量的更新,只需在飛行器飛出網(wǎng)格后完成更新即可。對于支持多線程的導航計算機,重力計算可以作為一個單獨的線程與慣導解算并發(fā)進行,將重力計算的時間分攤到多個慣導解算周期內(nèi),從而進一步降低計算機的計算壓力。

5 結(jié)論

本文研究了慣性基準系統(tǒng)重力補償技術(shù),闡明了基于EGM2008球諧模型的重力計算方法,并對真實重力矢量進行了計算,驗證了計算方法的有效性。通過仿真對EGM2008重力模型應用下的慣性基準系統(tǒng)的性能進行了分析,結(jié)果表明:針對特定航線和慣性器件精度水平,在純慣性運行期間,重力模型誤差是一定RNP要求飛行階段的重要誤差源,采用基于EGM2008模型的重力進行慣導解算可使定位誤差顯著減小,水平定位精度提升達到61%,對滿足一定RNP要求的慣性基準系統(tǒng)研制具有支撐作用。最后,面向受約束的實際應用場景,分析了模型階次與性能優(yōu)化的關系,仿真表明:在特定航線上使用裁剪模型,在保證精度性能的同時能帶來97%的計算效率優(yōu)化,為實際應用場景下的模型選擇提供了參考。