基于碼載偏離度的改進(jìn)自適應(yīng)Hatch濾波算法

胡杰， 章林， 朱倚嫻， 周玲

(1.中國電子科技集團(tuán)公司第二十八研究所， 江蘇 南京 210007； 2.空中交通管理系統(tǒng)與技術(shù)國家重點(diǎn)實(shí)驗(yàn)室， 江蘇 南京 210007；3.南通大學(xué) 機(jī)械工程學(xué)院， 江蘇 南通 226019； 4.運(yùn)城學(xué)院 物理與電子工程系， 山西 運(yùn)城 044000)

0 引言

基于標(biāo)準(zhǔn)衛(wèi)星信號(hào)的導(dǎo)航系統(tǒng)無法滿足終端區(qū)精密進(jìn)近著陸導(dǎo)航性能需求[1-2]，因此，國際民航組織(ICAO)提出在終端區(qū)對(duì)衛(wèi)星信號(hào)進(jìn)行增強(qiáng)。地基增強(qiáng)系統(tǒng)(GBAS)通過地面站子系統(tǒng)播發(fā)偽距校正值和完好性參數(shù)提高了機(jī)載導(dǎo)航性能[3]?；谌蚨ㄎ幌到y(tǒng)(GPS)L1 C/A碼GBAS已經(jīng)獲得ICAO I類精密進(jìn)近著陸許可認(rèn)證[4]。為提高終端區(qū)進(jìn)場(chǎng)線路設(shè)計(jì)靈活性，ICAO提出未來將利用GBAS代替?zhèn)鹘y(tǒng)陸基導(dǎo)航系統(tǒng)，用于支持飛機(jī)實(shí)現(xiàn)I類甚至更高精度等級(jí)精密進(jìn)近著陸。

載波相位測(cè)距精確高，且受多路徑影響小，因此在進(jìn)行導(dǎo)航數(shù)據(jù)處理前需要進(jìn)行平滑預(yù)處理，以提高系統(tǒng)精度[5]。文獻(xiàn)[6]給出了Hatch濾波器算法推導(dǎo)過程，其本質(zhì)是利用算術(shù)平均獲得偽距測(cè)量值低頻部分，利用時(shí)間差分算子獲得載波相位測(cè)量值高頻部分，以提高偽距測(cè)量精度，目前在衛(wèi)星導(dǎo)航數(shù)據(jù)處理領(lǐng)域得到廣泛應(yīng)用。由算法原理知，濾波器平滑時(shí)間越大，其帶寬越小，則噪聲抑制效果越好。但是衛(wèi)星信號(hào)觀測(cè)值中的偽碼- 載波偏離度會(huì)導(dǎo)致Hatch濾波器輸出一項(xiàng)額外的常值偏置誤差。根據(jù)分析知，該項(xiàng)誤差與濾波器平滑時(shí)間呈正比，因此必須選擇一個(gè)適合的平滑時(shí)間[7]。美國聯(lián)邦航空局(FAA)研究報(bào)告指出，GBAS中Hatch濾波器平滑時(shí)間宜設(shè)置為100 s固定值[8]，但是相關(guān)學(xué)者研究表明，對(duì)于所有可見衛(wèi)星在全時(shí)間段內(nèi)使用固定平滑時(shí)間缺乏合理性，例如，當(dāng)衛(wèi)星仰角較低時(shí)偽距所受多徑效應(yīng)較為明顯，此時(shí)宜提高濾波器平滑時(shí)間以抑制高頻噪聲。由于GPS L1是目前唯一獲得ICAO認(rèn)證的用于衛(wèi)星信號(hào)增強(qiáng)頻點(diǎn)[9]，為提高系統(tǒng)精度與完好性，許多學(xué)者進(jìn)行了相關(guān)技術(shù)攻關(guān)。文獻(xiàn)[10]為抑制算法發(fā)散對(duì)Hatch濾波器系數(shù)進(jìn)行歸一化，但是其最優(yōu)系數(shù)設(shè)置受限較大。文獻(xiàn)[11]建立了電離層延時(shí)分段線性模型，其相比常規(guī)線性模型更加符合實(shí)際特性，但該方法計(jì)算量相對(duì)較大，影響了算法的實(shí)時(shí)性。文獻(xiàn)[12]建立了平滑時(shí)間與電離層風(fēng)暴梯度數(shù)學(xué)模型，給出了濾波器最優(yōu)平滑時(shí)間估計(jì)方法，當(dāng)電離層異常時(shí)該方法假設(shè)平滑時(shí)間為無窮大，該假設(shè)在實(shí)際中并不成立，且在計(jì)算最優(yōu)平滑時(shí)間時(shí)并未考慮偽距測(cè)量噪聲影響。文獻(xiàn)[13]利用Klobuchar模型對(duì)電離層延時(shí)進(jìn)行建模并根據(jù)可見衛(wèi)星仰角大小自適應(yīng)確定Hatch濾波器平滑時(shí)間，分析該文獻(xiàn)可知，電離層延時(shí)建模結(jié)果與實(shí)際值之間存在較大偏差，因此導(dǎo)致濾波器平滑精度不能保證。文獻(xiàn)[14]提出一種電離層異常檢測(cè)算法，并利用低通濾波器抑制偽距測(cè)量噪聲以及實(shí)現(xiàn)電離層延時(shí)變化率在線估計(jì)，分析該文獻(xiàn)可知，電離層延時(shí)檢測(cè)靈敏度與低通濾波器時(shí)間常數(shù)呈反比，當(dāng)模擬數(shù)值大小為0.012 m/s的電離層延時(shí)變化率時(shí)，低通濾波器響應(yīng)時(shí)間為94 s，其估計(jì)電離層延時(shí)變化率的靈敏度低，且其在計(jì)算最優(yōu)平滑時(shí)間時(shí)未考慮偽距噪聲標(biāo)準(zhǔn)差對(duì)平滑時(shí)間估計(jì)影響。

本文提出一種適用于單頻GBAS的自適應(yīng)Hatch濾波算法，該算法根據(jù)衛(wèi)星信號(hào)偽碼- 載波偏離度在線估計(jì)電離層延時(shí)變化率和偽距測(cè)量噪聲標(biāo)準(zhǔn)差，并推導(dǎo)得到平滑后偽距誤差均方根與電離層延時(shí)變化率、偽距測(cè)量噪聲標(biāo)準(zhǔn)差以及濾波器平滑時(shí)間三者之間的函數(shù)模型，進(jìn)一步根據(jù)該函數(shù)表達(dá)式計(jì)算得到濾波器最優(yōu)平滑時(shí)間。

1 碼減載波建模及電離層延時(shí)變化率估計(jì)

衛(wèi)星信號(hào)偽距和載波相位測(cè)距觀測(cè)方程[15]可表示為

(1)

式中：ρ、φ分別表示接收機(jī)輸出偽距和載波相位測(cè)量值；r為一項(xiàng)等效距離，由衛(wèi)星與接收機(jī)之間幾何距離、接收機(jī)和衛(wèi)星鐘差以及對(duì)流層延時(shí)等組成；ι表示電離層延時(shí)；N表示載波相位整周模糊度；Mρ、Mφ分別表示偽距和載波相位測(cè)量噪聲。

任意t時(shí)刻利用偽距和載波相位測(cè)量值做差可得該顆衛(wèi)星的碼減載波，即

yt=ρt-φt=2ιt-N+Mt,ρ-Mt,φ，

(2)

式中：yt表示t時(shí)刻偽距和載波相位測(cè)量差值；ρt、φt分別表示t時(shí)刻偽距和載波相位測(cè)量值；ιt表示t時(shí)刻電離層延時(shí)；Mt,ρ、Mt,φ分別表示t時(shí)刻偽距和載波相位測(cè)量噪聲。

由(2)式可知，碼減載波值中包含2倍電離層延時(shí)、載波相位整周模糊度、偽距和載波相位測(cè)量噪聲等。電離層延時(shí)具有隨時(shí)間變化特性，但是在短時(shí)間內(nèi)一般可認(rèn)為其具有線性變化特性，因此，(2)式又可表示為

yt=b+2t·a+Mt,ρ-Mt,φ≈b+2t·a+Mt,ρ，

(3)

式中：b表示常值項(xiàng)；a表示電離層延時(shí)變化率，由于載波相位測(cè)量噪聲相比偽距測(cè)量噪聲要小很多，其可以忽略不計(jì)，將t1～tn時(shí)間段內(nèi)的碼減載波寫成矩陣形式，n表示大于1的自然數(shù)，其中tn-t1取值為60 s，可表示為

(4)

即

Y=Xβ+M.

(5)

(6)

圖1所示為一組GPS衛(wèi)星觀測(cè)數(shù)據(jù)中10號(hào)衛(wèi)星7 200 s內(nèi)碼減載波隨時(shí)間變化曲線，由于碼減載波中包含載波相位整周模糊度，設(shè)置碼減載波初始值為零。計(jì)算得到的碼減載波噪聲較大，本文利用滑動(dòng)平均濾波器進(jìn)行平滑預(yù)處理，平均滑口大小設(shè)置為800 s，經(jīng)過滑動(dòng)平均濾波后的碼減載波反映了當(dāng)前衛(wèi)星所受電離層延時(shí)變化大小和趨勢(shì)。

圖1 10號(hào)衛(wèi)星碼減載波曲線Fig.1 Code-minus-carrier curve of Satellite 10

進(jìn)一步利用最小二乘法可計(jì)算得到電離層延時(shí)變化率a，其中tn取值為60 s，因此系統(tǒng)每間隔60 s更新一次電離層延時(shí)變化率計(jì)算值。當(dāng)衛(wèi)星截止角設(shè)置為10°時(shí)，10號(hào)衛(wèi)星在天頂方向可見總時(shí)長(zhǎng)為20 254 s，因此共進(jìn)行了337次最小二乘求解電離層延時(shí)變化率，如圖2所示。為了驗(yàn)證最小二乘法估計(jì)得到的電離層延時(shí)變化率的準(zhǔn)確性，本文以雙頻載波相位測(cè)量值計(jì)算得到的電離層延時(shí)變化率作為參考比較對(duì)象，其計(jì)算過程為

(7)

ΔId,t=SL1,t-SL1,t-1，

(8)

式中：SL1,t、ΔId,t分別表示t時(shí)刻根據(jù)接收機(jī)輸出L1和L2頻點(diǎn)載波相位測(cè)量值計(jì)算得到的電離層延時(shí)及其變化率；fL1、fL2分別表示接收機(jī)輸出L1和L2測(cè)量頻率；φL1、φL2分別表示L1和L2頻點(diǎn)載波相位測(cè)量值；NL1、NL2分別表示L1和L2頻點(diǎn)載波相位整周模糊度。

圖2 電離層延時(shí)變化率估計(jì)曲線Fig.2 Estimated curve of ionospheric delay rate

對(duì)比圖2中藍(lán)色和紅色曲線可以看出，本文建立的碼減載波線性模型能夠較好地估計(jì)得到電離層延時(shí)變化率，在此基礎(chǔ)上可以分離得到偽距測(cè)量噪聲，并進(jìn)一步計(jì)算其標(biāo)準(zhǔn)差。

2 基于碼載偏離度的平滑時(shí)間自適應(yīng)確定方法

2.1 Hatch濾波器誤差分析

任意k時(shí)刻Hatch濾波器離散遞推形式[16]可表示為

(9)

式中：ρs,k、ρs,k-1分別表示k和k-1時(shí)刻平滑后偽距測(cè)量值；K表示濾波器平滑時(shí)間；ρk表示k時(shí)刻偽距測(cè)量值；φk、φk-1分別表示k和k-1時(shí)刻載波相位測(cè)量值。經(jīng)Hatch濾波器平滑后的偽距誤差εk可表示為

εk=ρs,k-ρk.

(10)

(9)式代入(10)式，可得平滑后偽距誤差表達(dá)式為

(11)

式中：Mρ，k表示k時(shí)刻偽距測(cè)量噪聲。

假設(shè)ε初始值為0，電離層延時(shí)變化率為a，則由(11)式可得平滑后偽距誤差遞推解為

(12)

(13)

由(13)式可知，平滑后偽距誤差包括系統(tǒng)誤差和隨機(jī)誤差兩部分，其中，λ表示系統(tǒng)誤差，Δ表示隨機(jī)誤差。系統(tǒng)誤差由電離層延時(shí)變化率和濾波器平滑時(shí)間構(gòu)成，隨機(jī)誤差則由偽距測(cè)量噪聲和濾波器平滑時(shí)間構(gòu)成。穩(wěn)態(tài)誤差εss可以表示為電離層延時(shí)變化率、偽距測(cè)量噪聲和濾波器平滑時(shí)間的函數(shù)。當(dāng)電離層延時(shí)變化率不等于零時(shí)，經(jīng)Hatch濾波器平滑后偽距會(huì)產(chǎn)生一項(xiàng)偏置誤差，其大小與濾波器平滑時(shí)間呈正比，例如，當(dāng)a=5 mm/s、濾波器平滑時(shí)間設(shè)置為100 s時(shí)，該誤差項(xiàng)最大為1 m，因此濾波器平滑時(shí)間不宜設(shè)置過大。同時(shí)，由偽距平滑后隨機(jī)誤差可知，濾波器平滑時(shí)間K值越大，則偽距噪聲抑制效果越好，為了提高偽距平滑精度需要增大K值，實(shí)際應(yīng)用中需要根據(jù)衛(wèi)星信號(hào)傳播環(huán)境自適應(yīng)確定濾波器平滑時(shí)間。

2.2 計(jì)算最優(yōu)平滑時(shí)間

由(13)式可知，平滑后偽距誤差中系統(tǒng)誤差和隨機(jī)誤差可表示為

λ=-2(K-1)a，

(14)

(15)

進(jìn)一步根據(jù)(14)式和(15)式，可計(jì)算得到平滑后偽距誤差均方根RMS為

(16)

(17)

(18)

λ2=4(K-1)2a2.

(19)

(18)式和(19)式代入(16)式中，有

(20)

由(20)式可知，平滑后偽距誤差均方根大小由Hatch濾波器平滑時(shí)間、電離層延時(shí)變化率以及偽距測(cè)量噪聲標(biāo)準(zhǔn)差決定，如(21)式所示：

RMS(εss)=f(K,a,σMρ).

(21)

分析(21)式可知，當(dāng)a和σMρ已知時(shí)，函數(shù)f(K,a,σMρ)最小值時(shí)所對(duì)應(yīng)的K值即為最優(yōu)平滑時(shí)間，對(duì) (21) 式求偏導(dǎo)數(shù)，可得

(22)

(22)式進(jìn)一步可整理簡(jiǎn)化為

(23)

本文利用Cardano法實(shí)現(xiàn)(23)式解析求解，其中Cardano法求解三次方程過程如(24)式～(27)式所示，假設(shè)三次方程表示為

x3+αx2+βx+γ=0，

(24)

式中：x表示三次方程待求解未知根；α，β，γ∈R.

利用Cardano方法可得(24)式的解析解為

(25)

由于濾波器平滑時(shí)間為正整數(shù)，(23)式的正實(shí)根取整即為所求最優(yōu)平滑時(shí)間Ko：

Ko=312(-q+q2+4p3)+312(-q-q2+4p3)+13,

(26)

圖3為最優(yōu)平滑時(shí)間Ko隨電離層延時(shí)變化率a和偽距測(cè)量噪聲標(biāo)準(zhǔn)差σMρ變化三維曲面圖，a取值范圍為0.1～8 mm/s，σMρ取值范圍為0.01～6 m. 由圖3可以看出，Ko取值與σMρ呈正比，與a呈反比，即σMρ值越大，a值越小，則對(duì)應(yīng)的Ko值越大。

圖3 最優(yōu)平滑時(shí)間Ko隨a和σMρ變化三維曲面圖Fig.3 Three-dimensional surface graph of the optimal smoothing time Ko changing with a and σMρ

綜上所述，本文所提出的自適應(yīng)Hatch濾波算法流程如圖4所示。

圖4 自適應(yīng)Hatch濾波算法流程Fig.4 Flowchart of adaptive Hatch filter algorithm

首先，依次計(jì)算所有可見衛(wèi)星的碼減載波，建立其1階線性模型，并利用最小二乘法估計(jì)電離層延時(shí)變化率a和偽距測(cè)量噪聲標(biāo)準(zhǔn)差σMρ；然后，利用1階線性時(shí)不變低通濾波器進(jìn)行濾波處理，以平滑高頻噪聲；最后，根據(jù)計(jì)算得到的a和σMρ以及最優(yōu)平滑時(shí)間計(jì)算公式確定得到Ko值，該值即為Hatch濾波器最優(yōu)平滑時(shí)間。

3 驗(yàn)證實(shí)驗(yàn)

3.1 實(shí)驗(yàn)環(huán)境

利用實(shí)驗(yàn)室研制GBAS原理樣機(jī)進(jìn)行驗(yàn)證實(shí)驗(yàn)，該系統(tǒng)實(shí)時(shí)接收GPS衛(wèi)星信號(hào)后經(jīng)過數(shù)據(jù)處理中心后輸出偽距差分校正值和完好性參數(shù)等。GBAS地面系統(tǒng)由基準(zhǔn)站、數(shù)據(jù)處理中心、VDB發(fā)射電臺(tái)等組成，模擬機(jī)載用戶(位置域)由機(jī)載多模接收機(jī)、機(jī)載數(shù)據(jù)處理中心以及VDB接收電臺(tái)等組成，GPS衛(wèi)星接收天線架設(shè)在實(shí)驗(yàn)室樓頂，如圖5所示。

圖5 GBAS地面系統(tǒng)GPS天線Fig.5 GPS antennas of GBAS ground subsystem

3.2 結(jié)果分析

以一組用戶站接收機(jī)輸出GPS衛(wèi)星信號(hào)偽距和載波相位測(cè)量數(shù)據(jù)為例，對(duì)本文所提算法進(jìn)行驗(yàn)證，其中，碼減載波滑動(dòng)平滑窗口長(zhǎng)度設(shè)置為800 s，a和σMρ每間隔60 s更新一次，因此最優(yōu)平滑時(shí)間Ko也是每間隔60 s更新計(jì)算一次。圖6為其中1號(hào)衛(wèi)星估計(jì)得到的a和σMρ隨時(shí)間變化曲線，圖7為利用估計(jì)得到的a和σMρ計(jì)算得到的Ko隨時(shí)間變化曲線，圖8為1號(hào)衛(wèi)星仰角隨時(shí)間變化曲線。

圖6 1號(hào)衛(wèi)星a和σMρ隨時(shí)間變化曲線Fig.6 a and σMρ of Satellite 1 changing with time

圖7 1號(hào)衛(wèi)星Ko隨時(shí)間變化曲線Fig.7 Ko of Satellite 1 changing with time

圖8 1號(hào)衛(wèi)星仰角隨時(shí)間變化曲線Fig.8 Elevation angle of Satellite 1 changing with time

由圖6和圖8可以看出：σMρ隨著衛(wèi)星仰角的增大而減小，表明衛(wèi)星仰角越低，偽距所受多路徑越大；a值變化較為平穩(wěn)，也進(jìn)一步表明電離層延時(shí)在較短時(shí)間內(nèi)可認(rèn)為具有線性特性。由圖7可以看出，當(dāng)電離層變化平穩(wěn)時(shí)，Ko隨σMρ值的增大而變大，即衛(wèi)星仰角較低時(shí)，宜增大Ko值，以抑制偽距測(cè)量隨機(jī)誤差。

圖9給出了模擬機(jī)載用戶3 600 s內(nèi)差分定位誤差曲線。為了凸顯本文所提算法的有效性，分別利用固定平滑時(shí)間Hatch濾波器和本文所提自適應(yīng)Hatch濾波算法對(duì)偽距進(jìn)行平滑，δx、δy和δz分別表示x軸、y軸和z軸3個(gè)方向位置誤差。圖10給出了該時(shí)間段內(nèi)GPS可見衛(wèi)星仰角隨時(shí)間變化曲線。

由圖9可以看出，自適應(yīng)濾波算法位置解算精度較高，特別是在1 200～2 000 s之間，自適應(yīng)濾波算法精度優(yōu)勢(shì)更加明顯，根據(jù)圖10可知，在1 200～2 000 s之間最多有4顆衛(wèi)星仰角低于15°，即20號(hào)、22號(hào)、29號(hào)以及12號(hào)衛(wèi)星，因此這一時(shí)間段內(nèi)相應(yīng)偽距所受的多路徑也較大，利用自適應(yīng)濾波算法能夠更好地起到平滑濾波效果，位置解算精度也相對(duì)較高。

為進(jìn)一步驗(yàn)證本文所提出算法的有效性，分不同時(shí)間段(早上8:00—9:00，中午11:00—12:00，晚上19:00—20:00)進(jìn)行3組驗(yàn)證實(shí)驗(yàn)，每組持續(xù)時(shí)間為3 600 s，表1為兩種濾波方法對(duì)應(yīng)的位置誤差均方根值對(duì)比統(tǒng)計(jì)結(jié)果。

由表1可以看出，中午太陽光輻射強(qiáng)時(shí)，衛(wèi)星信號(hào)所受電離層延時(shí)最大，因此其對(duì)應(yīng)差分定位精度也相對(duì)最差。同時(shí)可以看出，當(dāng)采用本文所提出的自適應(yīng)Hatch濾波算法時(shí)3個(gè)方向位置解算精度相比傳統(tǒng)100 s固定平滑時(shí)間算法都得到了提升，其中，實(shí)驗(yàn)2中x軸方向位置精度最大提高了53.19%，由此進(jìn)一步驗(yàn)證了本文所提方法的正確性。

圖9 不同濾波方法位置誤差對(duì)比曲線Fig.9 Comparison of position errors of different filtering methods

圖10 可見衛(wèi)星仰角變化曲線Fig.10 Changing curves of elevation angles of all visible satellites

表1 位置誤差均方根值對(duì)比統(tǒng)計(jì)Tab.1 Statistical comparison of RMSs of position errors

4 結(jié)論

為提高單頻GBAS中載波相位平滑偽距精度，本文針對(duì)Hatch濾波器平滑精度易受電離層延時(shí)和偽距測(cè)量噪聲影響的問題進(jìn)行了分析和研究。得到以下主要結(jié)論：

1) 電離層延時(shí)短時(shí)間內(nèi)具有線性變化特性，碼減載波值反映了2倍電離層延時(shí)，通過建立碼減載波1階線性模型，在此基礎(chǔ)上利用最小二乘法能夠估計(jì)得到電離層延時(shí)變化率。

2) 根據(jù)Hatch濾波算法推導(dǎo)得到平滑后偽距誤差均方根函數(shù)模型，該函數(shù)模型大小由電離層延時(shí)變化率、偽距測(cè)量噪聲標(biāo)準(zhǔn)差以及平滑時(shí)間決定，由此進(jìn)一步可得濾波器最優(yōu)平滑時(shí)間。

3) 進(jìn)行了系統(tǒng)驗(yàn)證實(shí)驗(yàn)，結(jié)果表明本文所提算法能夠根據(jù)衛(wèi)星信號(hào)實(shí)際環(huán)境自適應(yīng)確定最優(yōu)平滑時(shí)間；模擬機(jī)載用戶(位置域)差分定位結(jié)果表明本文所提出的自適應(yīng)濾波算法3個(gè)方向位置解算精度相比固定100 s平滑時(shí)間都有提高，表明了本文所提算法的正確性。