基于載波高斯包絡模型的精密單點定位保護級

王吉濤,許承東,趙 靖,蘭曉偉

(1.北京理工大學宇航學院,北京 100081;2. 中國交通通信信息中心,北京 100011)

1 引言

近年來,隨著GNSS導航服務性能的不斷提高,精密單點定位(Precise Point Positioning,PPP)這一高精度定位技術逐漸成為關注熱點[1],正逐步應用于智能交通系統、鐵路和航海等領域。接收機自主完好性監測(Receiver Autonomous Integrity Monitoring,RAIM)作為提高定位可靠性的一種重要手段,其在PPP中的應用也受到越來越多的關注[2],[3]。

RAIM算法一般通過經驗噪聲模型估計觀測噪聲的分布情況,進而計算定位誤差限值即保護級(Protection Level,PL)來評估完好性風險[4],[5]。常用的經驗模型主要包括衛星高度角法[6]、信噪比方法[7]和高度角信噪比聯合法[8][9]等,上述方法大都根據經驗選取方差因子,利用觀測噪聲水平和高度角、信噪比等因素的關系建立隨機模型,再在此基礎上進行定位解算和完好性監測。但是在PPP實際應用中,觀測噪聲易受多路徑效應影響,出現“厚尾”現象。針對“厚尾”現象,大多采用方差膨脹方法對觀測噪聲尾部包絡[10][11],但是該方法對于觀測噪聲的分布特性估計過于保守,易造成PPP完好性監測算法的低可用性。此外,目前研究主要利用大量觀測數據針對偽距噪聲[12][13]或雙差載波噪聲[14][15]隨機模型開展研究,由于載波多路徑和接收機噪聲(Phase Multipath and Noise,PMN)提取相對困難,直接針對非差PMN的隨機模型研究相對較少。

鑒于此,本文利用三頻消電離層組合觀測值提取三頻非差組合PMN序列,建立真實環境下的PMN分布模型,并基于灰狼優化算法求解模型參數,使得保護級計算結果更為嚴密,以期為PPP的應用推廣提供參考。

2 PPP保護級計算框架

PPP模型通常采用消電離層(Iono-Free,IF)組合觀測值消除電離層延遲誤差,并利用擴展卡爾曼濾波進行定位解算[16],狀態方程和觀測方程分別為

Xk=Φk/k-1Xk-1+wk

(1)

Zk=HkXk+vk

(2)

通常wk和vk為零均值高斯白噪聲且互不相關,即wk～N(0,Qk),vk～N(0,Rk),Qk和Rk分別為wk和vk的協方差矩陣。Qk可由隨機游走過程描述。假定各頻段上碼觀測值和相位觀測值噪聲水平分別一致,則IF組合觀測噪聲方差與單個頻段的噪聲方差關系為

(3)

(4)

(5)

觀測噪聲協方差矩陣可以表示為

(6)

式中,diag[·]表示對角矩陣,n為可見星數目。

PPP定位解算的EKF方程為

(7)

完好性監測中常用保護級包絡定位誤差,基于特征斜率的保護級計算通常表述為[16]

(8)

式中,PL0為無故障情況下的保護級;PL1為有故障情況下的保護級;k0和k1為放大系數,由完好性風險和連續性風險確定;σp為定位誤差標準差,由Pk主對角線元素得到;λ為非中心化參數,由連續性風險、可見星數目和衛星故障概率確定;Slopei為特征斜率,用于描述檢驗統計量和定位誤差之間的線性關系。特征斜率的計算方法為[17]

(9)

式中,VSlopei為衛星i的垂向特征斜率;HSlopei為衛星i的水平特征斜率;矩陣下標ij表示矩陣的第i行第j列元素索引。

由此可以看到,在保護級的計算過程中,特征斜率和定位誤差標準差均直接或間接與Rk相關,因此載波噪聲模型的精確與否直接影響完好性監測算法的可用性。

3 PMN高斯包絡模型

3.1 載波觀測噪聲序列提取

三頻載波觀測值兩兩組合得到兩組IF組合觀測值,通過作差消除載波噪聲以外的誤差項,進而得到三頻組合PMN,具體表示為[18]

LDIF=(c1L1+c2L2)-(c3L1+c4L5)

=(c1-c3)L1+c2L2-c4L5

=(c1-c3)εL1+c2εL2-c4εL5+BDIF

(10)

式中,LDIF為三頻載波觀測值組合;Lj(j=1,2,5)為載波觀測值;BDIF為模糊度項和硬件延遲偏差;εLj((j=1,2,5)為頻段j的PMN;ci(1,2,3,4)為組合系數,表達式為

(11)

(12)

(13)

(14)

式中,f5為L5載波頻率。

在無周跳情況下,BDIF為常值。通過分段區間去均值處理得到三頻組合PMN序列

εPMN=(c1-c3)εL1+c2εL2-c4εL5

(15)

根據誤差傳播原理并基于各頻段載波噪聲水平一致假設,三頻非差組合PMN序列噪聲方差與單頻載波噪聲方差關系為

εPMN=(c1-c3)εL1+c2εL2-c4εL5

(16)

結合式(3)和式(16),通過求解三頻非差組合PMN序列的噪聲分布模型可以間接獲得IF組合的觀測噪聲模型,從而實現保護級計算。

3.2 三頻組合PMN分布模型

基于高度角的經驗噪聲模型[19]為

(17)

式中,σ0為經驗噪聲方差,θ為衛星高度角。對于偽距噪聲,σ0=0.3m;對于載波噪聲,σ0=0.003m。

該經驗模型對良好環境中的觀測噪聲描述程度較好,但是估計較為保守,導致計算的保護級結果偏大,不符合實際情況。文獻[12]提出一種可用于非高斯觀測噪聲分布估計的觀測噪聲分布模型,即

(18)

對于偽距噪聲,模型中各參數取值見表1。

表1 偽距噪聲模型參數[12]

本節主要在該模型基礎上,利用三頻組合PMN序列,求解針對載波噪聲的模型參數。與偽距噪聲類似,三頻組合PMN由接收機噪聲和多路徑誤差組成,其中接收機噪聲與高度角無關,多路徑誤差與高度角指數相關,當衛星高度角高于某閾值時,PMN由高斯模型描述;當高度角低于某閾值時,PMN由高斯膨脹模型描述。PMN仍服從零均值高斯分布,其概率分布函數(CDF)為

(19)

式中,f(x;0,σ2)表示均值為0,方差為σ2的正態分布概率密度函數。

3.3 模型參數求解

(20)

這樣,不同高度角區間內PMN的真實CDF可以用頻數近似描述為

(21)

式中num(·)表示對應高度角區間中PMN的統計數。

所有高度角區間對應PMN的實際CDF用矩陣形式描述為

(22)

類似地,由式(19)確定的PMN的理論CDF用矩陣形式描述為

(23)

由4個模型參數確定的理論CDF和實際CDF應當保持一致,即

ΔF=-F=0

(24)

由此建立優化求解問題為

(25)

滿足式(25)的σ1,σ2,a,η即為理論上的模型參數組合。該目標函數中含有正態分布函數,該問題是典型的多維非線性優化問題,采用一般線性優化方法很難保證求解結果的正確性和計算效率。因此,為快速獲得準確的PMN模型參數,本文基于灰狼優化算法[20](Grey Wolf Optimizer,GWO)并結合蒙特卡洛方法實現參數優化求解。GWO算法是一種新型群智能優化算法,通過模擬灰狼群體在捕食過程中的覓食行為實現目標優化。在GWO算法中,狼群由高到低排序為α,β,δ,ω4個等級,其中α,β,δ灰狼等級最高,其位置向量表示3組帕累托最優解,ω灰狼位置根據α,β,δ更新,狼群通過等級更替實現最優解的搜索,更新過程如下

(26)

(27)

式中,t為當前迭代次數;Xα,Xβ,Xδ分別為α,β,δ狼個體位置向量;X表示ω狼個體的位置向量;D表示狼個體與最優解的距離;C和H為權重系數。在迭代過程中,C和H根據如下規則更新

(28)

式中,T為最大迭代次數;r1,r2為[0,1]隨機數。

本文中,X為PMN分布模型參數,即

X=[σ1,σ2,a,η]

(29)

本文將GWO算法應用到PMN分布模型參數的求解當中,通過對4個模型參數的快速尋優,以達到提高模型求解準確性和高效性的目的?；贕WO算法的模型參數求解流程如圖1所示。具體實現過程包括以下步驟:

圖1 基于GWO的模型參數求解流程

1)對PMN序列進行高度角和數值區間劃分,統計各區間PMN樣本頻數,建立實際CDF矩陣。

2)初始化灰狼種群數量M和最大迭代次數T等參數,設置參數σ1,σ2,a和η的取值范圍,選擇式(25)作為適應度函數。

3)計算狼群個體位置及其適應度。

4)遍歷灰狼種群,計算所有個體適應度函數,將適應度函數值排序前3位的灰狼個體標記為α,β,δ。

5)根據式(26)(27)計算ω狼與α,β,δ狼的距離,并更新α,β,δ狼的位置。

6)根據式(28)更新算法中參數C和H。

7)判斷是否達到最大迭代次數T,若達到則保存最優組合優化解Xα,即模型參數的最優值,否則返回步驟4)。

4 仿真分析

選取JFNG站100天(2020年1月1日至4月9日)的觀測數據,接收機型號為Trimble NetR9,天線為“TRM55971.00”,采樣率為30s。由于GPS衛星目前只有Block ⅡF類型衛星具有三頻觀測數據,因此只提取該類型衛星的PMN序列,觀測數據總量為1.04×106,可以用于高度角閾值和模型參數的確定。

4.1 高度角閾值確定

PMN在不同高度角范圍具有不同的噪聲水平,利用Q-Q Plot方法、峰度和偏度系數對PMN的高斯特性進行分析。以G32衛星為例,按照高度角間隔Δθ=5°將PMN序列從0到90°劃分為18組,繪制PMN噪聲序列在10°～15°,20°～25°,45°～50°和60°～65°四個高度角區間內的Q-Q Plot圖,如圖2所示。

圖2 載波噪聲Q-Q Plot

從圖2中可以看到,GPS衛星在高度角較大情況下的PMN序列完全符合高斯分布,但在10°～15°小高度角區間下,PMN序列基本符合高斯分布,但是在“尾部”與高斯分布存在偏離。

為檢驗PMN序列的對稱性和“厚尾”情況,分別計算衛星在不同高度角區間的偏度值和峰度值,結果如圖3所示?？梢钥吹?各高度角區間中載波噪聲對應的偏度始終在零值附近波動,說明PMN序列的分布具有對稱性,高度角的變化對PMN偏度影響較小。在0～20°高度角區間內,PMN峰度值顯著大于3,PMN序列的分布具有尖頂和“厚尾”特征;在大于20°高度角區間中,PMN峰度值趨于一致,近似為3,說明基本服從高斯分布。

圖3 載波噪聲的偏度和峰度值

不同高度角下的載波噪聲QQ-Plot、偏度和峰度值表明,應分段描述不同高度角區間的噪聲分布特性,這進一步驗證了PMN分布模型的合理性,并可確定高度角閾值為20°。當高度角大于20°時,高斯分布對PMN的分布描述效果較好。但當高度角小于20°時,PMN分布呈現厚尾特性,若仍用高斯分布描述其分布特性,與實際情況不符,因此采用膨脹系數η進行高斯包絡。

4.2 模型參數求解

將每組高度角區間中的PMN序列從小到大排列后劃分為200個子區間,高度角閾值取為20°。為削弱GWO算法不同初值選取對參數優化結果的影響,采用蒙特卡洛方法進行多次計算,蒙特卡洛次數設為50,取最優估計結果,得到各衛星模型參數估計值,如圖4所示。

圖4 模型參數估計結果

通過比較可以發現,各衛星的PMN分布模型參數基本一致,取各衛星模型參數的平均值作為最終結果,得到三頻組合PMN的分布模型參數,并由式(16)得到單頻PMN模型參數,具體結果見表2。

表2 PMN模型參數

4.3 觀測噪聲包絡分析

為驗證所求模型參數的有效性,采集一組實測數據并提取三頻非差組合PMN序列,觀測數據信息見表3。分別利用經驗模型和PMN分布模型對載波噪聲進行3σ包絡分析,結果如圖5所示。

圖5 兩種模型噪聲包絡情況

表3 觀測數據信息

提取所有衛星在5°～20°和45°～60°高度角區間中的三頻組合PMN序列,對比經驗模型和PMN分布模型在低高度角和大高度角兩種情況下對載波噪聲概率密度的包絡情況,結果如圖6和圖7所示。

圖6 低高度角噪聲模型概率密度包絡

圖7 大高度角噪聲模型概率密度包絡

結合圖5、圖6和圖7結果可以發現,兩種模型均能對載波噪聲進行包絡,但是經驗模型對PMN的噪聲分布情況估計明顯過于保守,與經驗模型相比,PMN分布模型的包絡曲線與真實PMN貼合更加緊密,對載波噪聲包絡效果更好。

4.4 保護級結果評估

為進一步驗證所提載波噪聲模型在PPP完好性監測中的性能,選取機載動態觀測數據進行保護級計算。觀測時長約為30分鐘,采樣率為1 s,接收機類型為“NovAtel OEM7”,衛星天空圖和飛行軌跡分別如圖8和圖9所示。

圖8 衛星天空圖

圖9 載體飛行軌跡

動態定位解算模式采用IF組合PPP浮點解,衛星軌道和鐘差采用IGS改正產品改正,對流層誤差等采用模型改正,衛星截止高度角設置為7°。兩種隨機模型方案包括:①式(17)所示的經驗噪聲模型②本文所求解的載波噪聲模型。用于計算保護級的完好性參數設置見表4。

表4 完好性參數設置

圖10和圖11展示了兩種隨機模型方案的定位誤差和保護級計算結果?？梢钥吹絇MN分布模型對應的VPL和HPL始終小于經驗噪聲模型,VPL值最大相差1.13 m(第1794歷元),HPL值最大相差0.45 m(第1789歷元)。

圖10 垂向保護級和定位誤差

圖11 水平保護級和定位誤差

兩種保護級收斂后的統計信息見表5。與經驗噪聲模型相比,PMN分布模型的VPL平均降低了38.1%,HPL平均降低了27.1%,說明后者模型更加精確,克服了經驗模型保護級計算結果偏保守的缺點。

表5 保護級的統計信息

5 結論

為解決經驗噪聲模型計算精密單點定位保護級具有保守性的問題,本文基于GWO算法求解了一種載波噪聲分布模型。仿真結果表明,相比于經驗觀測噪聲模型,文中所提PMN分布模型對于載波噪聲的分布描述更為精確,使精密單點定位的垂向保護級降低了38.1%,水平保護級降低了27.1%,提高了完好性監測的可用性。