路 陽,何 偉,王 樂
(中國電子科技集團公司第二十研究所,西安 710068)
GBAS 在飛機起降的應用,屬于高端應用,能有效地提高飛機起降引導精度,為飛機航行保駕護航。但在實際應用中,由于新增衛(wèi)星的偽距精度差,采用最小二乘解算,會引起定位解算結果的變化,從而造成高度跳變。在本文中,采用有效的抗差估計方法,能有效地提高偽距收斂精度,縮短收斂時間,提升GBAS 差分定位在飛機航行引導中的可靠性。
在GBAS 差分定位中,地面站偽距差分修正量的正確與否直接影響到機載端接收機的定位精度。實際應用中差分量生成流程如圖1所示。
圖1 差分量生成流程
生成步驟如下:
(1)原始殘差生成
一顆衛(wèi)星(編號為i)通過星歷外推可以得到t時刻衛(wèi)星位置為(xi,yi,zi),差分衛(wèi)星地面站基準接收機天線的位置為(xr,yr,zr),從差分衛(wèi)星地面站接收機天線r到衛(wèi)星i的幾何距離(真距)rr(i)為:
對4 個接收機同時收星,衛(wèi)星(編號為i)在t時刻的偽距測量值為ρr(i)(m),則含有殘差的偽距差分修正量為ρc(i)(m):
式中,i代表衛(wèi)星號,i=1,2,3…32;m=1,2,3,4,代表4個接收機號。
(2)原始校正差分量的生成
由于式(2)中計算公式所得的差分量含有本地鐘差,雖然對機載差分定位沒有影響,但為了減少空地傳輸字節(jié)量,一般采用如下公式,估算本地鐘差并消除:
(3)融合差分量的生成
對4 臺校正本地鐘差后的差分修正量進行融合處理,生成最終的差分修正量:
機載接收機(編號為u)對共視衛(wèi)星i的偽距測量值為ρu(i),則修正后的機載關于這顆星的偽距測量值為:
機載和地面有N顆共視衛(wèi)星,機載接收機修正后的偽距定位方程為:
通過解算式(7)可以計算得到機載接收機的差分定位結果(x,y,z)。再通過坐標轉換得到大地坐標系(Lon′ ,Lat′ ,Height′)。
當i號低仰角衛(wèi)星剛剛跟蹤上,偽距的跟蹤精度還沒穩(wěn)定時,偽距測量值為ρr(i)中包含的有誤差分量 Δρr(i),Δρr(i)精度收斂圖如圖2所示。在接收機剛跟蹤上衛(wèi)星時,產生的偽距誤差達到3 m 以上,隨著時間的推移,偽距誤差逐漸減小,由實踐經驗可知通常在40 s 以后,能夠穩(wěn)定在厘米級。圖3所示為地面?zhèn)尉嗾`差引起對應的機載端高度變化圖。
圖2 偽距精度收斂圖
圖3 機載端高度變化圖
i號衛(wèi)星偽距未達穩(wěn)定時,對4 個接收機會產生錯誤的差分修正量ρc′(i)為:
采用上述方法估算本地鐘差并消除:
對4 臺校正本地鐘差后的差分修正量進行融合處理,生成最終的差分修正量為:
機載接收機收到錯誤的差分修正量后,修正共視衛(wèi)星i的偽距測量值,則得到帶有偏差的修正后的偽距為:
機載和地面有N顆共視衛(wèi)星,其中包含一顆含有錯誤偽距修正值的機載偽距值,機載接收機的偽距定位方程:
式(13)中,由于第一個方程式中的錯誤偽距,導致在解算方程組時得到誤差較大的機載差分定位(x′,y′,z′),坐標轉換得到錯誤的大地坐標系(Lon′ ,Lat′ ,Height′)。因此,當一顆低仰角衛(wèi)星被接收機接收到,偽距測量值的跟蹤精度還沒到達穩(wěn)定時,會產生錯誤的機載修正偽距,導致式(13)解算出錯,產生高度跳變。從圖2~圖3也可以看出,只有當地面接收機的偽距跟蹤精度穩(wěn)定后,解算上述方程組才會得到正確的定位結果,從而機載的差分定位結果恢復到正確的定位值[1]。
由此看出,只有偽距精度達到穩(wěn)定狀態(tài)后,才能生成正確的偽距差分修正量。傳統(tǒng)方案根據實踐經驗,當新增衛(wèi)星出現60 s 后,再讓其參與差分運算,但此方法效率不高,較為簡單原始,缺乏理論依據。而本文采用的自適應抗差估計算法能夠對新增衛(wèi)星偽距精度作出靈活自主判斷,將其應用到新增衛(wèi)星偽距測量值的預處理中,使問題得到更好的解決。
2.1.1 傳統(tǒng)最小二乘估計
設觀測向量 L= [L1L2…Ln]T,Li(i=1,2,…,n)獨立,未知參數向 X= [X1X2…Xm]T,其誤差向量Δ,線性化后的觀測方程為:
誤差方程為:
經典最小二乘(LS)估計要求:
由式(17)可得:
式中,∑X為X 的協方差矩陣;為方差因子估值。
我們以偽距雙差作為對偽距精度的估計,并假設雙差值各歷元間相互獨立,由此可得式(18)~式(20)中權陣P 為n階單位陣,系數矩陣A 是各元素值均為1 的n×1 維列向量,故:
2.1.2 極大似然性估計(M 估計)
設觀測樣本L1,L2,…,Ln,測值Li分布密度為f(X-Li),極大似然估計要求參數估值滿足:
用ρ(·)代替 -Inf(·),則極大似然原理要求:
即M 估計轉化為求上述極小化問題的解,對式(22)求導后可得:
式中,ψ是ρ的導函數。
對于M 估計,將式(25)改寫,令:
將式(26)代入式(25)可得:
式中,iW取Tukey 的ψ函數所相應的值,即:
式中,C=6.0,
2.2.1 相關最小二乘估計
設各觀測值L1,L2, …,Ln均服從正態(tài)分布,故其聯合概率密度函數為:
式中,L 為n維測值向量;∑L為L 的協方差矩陣;為測值協方差矩陣的行列式。式(30)可簡寫為:
實際求解時真實誤差 Δ = L -E(L)無法直接求得,故只能求測值殘差 V = L -,其中為L的估計值,因此式(32)可轉換為:
將式(15)代入式(34)后求導可得:
表示成矩陣形式,可得估計值[2]為:
2.2.2 相關IGG3 方案
相關IGG3 方案基于相關最小二乘估計,且根據測值的不同情況分別作出處理,即當測值偏差較小時(測值主部)采用效率高的LS 估計;當測值超出一定范圍時采用降權估計;當個別測值明顯過大時,采用淘汰法(零權估計)。由此可得如下相關等價權[3]:
式中,k0可取值為1.0~1.5,k1可取值為2.5~3.0。
偽距精度采用偽距雙差進行估計,偽距雙差變量X的近似真值為0 m,對其進行了10 次觀測得到:L=[0.099268,0.071994,0.027277,0.099270,0.108213,0.072887,0.077807,0.061709,0.094797,0.028172]。將L的首個歷元測值L(1)加入大約0.3 m 的粗差后變?yōu)長′,然后對L′分別用M 估計、相關最小二乘估計、相關IGG3 方案這三種抗差估計算法進行處理,并使各算法最終收斂時,每種抗差方法的最終收斂結果、迭代次數及最終收斂精度如表1所示。
表1 三種抗差算法收斂結果對照表
各抗差算法收斂結果對比圖如圖4所示。
圖4 三種抗差算法收斂效果對比圖
由以上試驗結果可得如下結論:從最終收斂結果方面,相關IGG3 方案最優(yōu),其次為M 估計,最后為相關最小二乘估計;從迭代次數方面,相關IGG3 方案和相關最小二乘估計均為2 次,算法效率最高,而M 估計迭代5 次,效率最低;從最終收斂精度方面,IGG3 方案最大,說明其收斂速度最快,而M 估計的最終收斂精度最小,說明其收斂速度最慢,相關最小二乘估計收斂速度介于兩者之間。
表2為某地面接收機靜態(tài)下的連續(xù)28 個觀測歷元的偽距雙差觀測值[4]。從數據變化上看,各歷元觀測值圍繞偽距雙差的近似真值0 上下小幅波動,我們依舊將觀測值首個歷元加入大約0.3 m 的粗差后,分別采用三種抗差算法進行估計,試驗結果如圖5所示。
表2 接收機連續(xù)28個觀測歷元的偽距雙差觀測值
由圖5可以明顯得出如下結論:三種抗差算法均能較好剔除較大的偽距粗差(首個歷元跳點被“過濾”掉),從估計值準確度方面,起始階段具有獨立性的M 估計相比于其他兩個具有相關性的估計算法效果較好,但隨著歷元數不斷增多,后期(從14 歷元起)兩個相關估計算法占優(yōu),具有與真值更為接近的估計,相比較而言,IGG3 方案比相關LS 更能有效剔除觀測值明顯異常的點,并對部分測值較高點采用降權估計,因此采用IGG3 方案更為合適;從估計算法的連續(xù)性方面,具有獨立性的M 估計連續(xù)性最好,而兩個相關估計算法由于需要運用到權陣P,但在某些連續(xù)歷元相關性較差的時刻權陣P 無法生成,因此會導致這些時刻無法生成估計值,故連續(xù)性較差。
圖5 三種方案抗差結果圖
通過上述分析可知,任何算法都不會十全十美,正所謂“魚與熊掌不可兼得”,必須在算法估值準確度和連續(xù)性之間做出取舍,因此實際應用中,將M 估計和IGG3 方案綜合使用效果更好,即在連續(xù)歷元相關性較好的時刻,權陣P 存在,故此時采用相關IGG3 方案,而在連續(xù)歷元相關性較差的時刻,權陣P 不存在,故此時采用獨立M 估計,具體試驗結果如圖6所示。
圖6 M 與IGG3 組合抗差結果圖
由圖6可知,原始偽距雙差在34 s 以后精度達到0.1 m 以下,而采用M 與IGG3 組合抗差估計算法后,偽距雙差在13 s 以后就能達標,由此可見,該算法適用于接收機初始跟蹤鎖定階段的偽距抗差校正,能極大縮短接收機偽距的收斂時間,并能適當提高偽距精度[5]。
通過以上分析可知,采用M 與IGG3 組合抗差估計算法并將其應用到新增衛(wèi)星偽距測量值的預處理中,可極大提高接收機的收斂速度和精度。