基于載波高斯包絡(luò)模型的精密單點(diǎn)定位保護(hù)級(jí)

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

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

1 引言

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

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

鑒于此,本文利用三頻消電離層組合觀測值提取三頻非差組合PMN序列,建立真實(shí)環(huán)境下的PMN分布模型,并基于灰狼優(yōu)化算法求解模型參數(shù),使得保護(hù)級(jí)計(jì)算結(jié)果更為嚴(yán)密,以期為PPP的應(yīng)用推廣提供參考。

2 PPP保護(hù)級(jí)計(jì)算框架

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

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

(1)

Zk=HkXk+vk

(2)

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

(3)

(4)

(5)

觀測噪聲協(xié)方差矩陣可以表示為

(6)

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

PPP定位解算的EKF方程為

(7)

完好性監(jiān)測中常用保護(hù)級(jí)包絡(luò)定位誤差,基于特征斜率的保護(hù)級(jí)計(jì)算通常表述為[16]

(8)

式中,PL0為無故障情況下的保護(hù)級(jí);PL1為有故障情況下的保護(hù)級(jí);k0和k1為放大系數(shù),由完好性風(fēng)險(xiǎn)和連續(xù)性風(fēng)險(xiǎn)確定;σp為定位誤差標(biāo)準(zhǔn)差,由Pk主對角線元素得到;λ為非中心化參數(shù),由連續(xù)性風(fēng)險(xiǎn)、可見星數(shù)目和衛(wèi)星故障概率確定;Slopei為特征斜率,用于描述檢驗(yàn)統(tǒng)計(jì)量和定位誤差之間的線性關(guān)系。特征斜率的計(jì)算方法為[17]

(9)

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

由此可以看到,在保護(hù)級(jí)的計(jì)算過程中,特征斜率和定位誤差標(biāo)準(zhǔn)差均直接或間接與Rk相關(guān),因此載波噪聲模型的精確與否直接影響完好性監(jiān)測算法的可用性。

3 PMN高斯包絡(luò)模型

3.1 載波觀測噪聲序列提取

三頻載波觀測值兩兩組合得到兩組IF組合觀測值,通過作差消除載波噪聲以外的誤差項(xiàng),進(jìn)而得到三頻組合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為模糊度項(xiàng)和硬件延遲偏差;εLj((j=1,2,5)為頻段j的PMN;ci(1,2,3,4)為組合系數(shù),表達(dá)式為

(11)

(12)

(13)

(14)

式中,f5為L5載波頻率。

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

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

(15)

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

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

(16)

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

3.2 三頻組合PMN分布模型

基于高度角的經(jīng)驗(yàn)噪聲模型[19]為

(17)

式中,σ0為經(jīng)驗(yàn)噪聲方差,θ為衛(wèi)星高度角。對于偽距噪聲,σ0=0.3m;對于載波噪聲,σ0=0.003m。

該經(jīng)驗(yàn)?zāi)Ｐ蛯α己铆h(huán)境中的觀測噪聲描述程度較好,但是估計(jì)較為保守,導(dǎo)致計(jì)算的保護(hù)級(jí)結(jié)果偏大,不符合實(shí)際情況。文獻(xiàn)[12]提出一種可用于非高斯觀測噪聲分布估計(jì)的觀測噪聲分布模型,即

(18)

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

表1 偽距噪聲模型參數(shù)[12]

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

(19)

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

3.3 模型參數(shù)求解

(20)

這樣,不同高度角區(qū)間內(nèi)PMN的真實(shí)CDF可以用頻數(shù)近似描述為

(21)

式中num(·)表示對應(yīng)高度角區(qū)間中PMN的統(tǒng)計(jì)數(shù)。

所有高度角區(qū)間對應(yīng)PMN的實(shí)際CDF用矩陣形式描述為

(22)

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

(23)

由4個(gè)模型參數(shù)確定的理論CDF和實(shí)際CDF應(yīng)當(dāng)保持一致,即

ΔF=-F=0

(24)

由此建立優(yōu)化求解問題為

(25)

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

(26)

(27)

式中,t為當(dāng)前迭代次數(shù);Xα,Xβ,Xδ分別為α,β,δ狼個(gè)體位置向量;X表示ω狼個(gè)體的位置向量;D表示狼個(gè)體與最優(yōu)解的距離;C和H為權(quán)重系數(shù)。在迭代過程中,C和H根據(jù)如下規(guī)則更新

(28)

式中,T為最大迭代次數(shù);r1,r2為[0,1]隨機(jī)數(shù)。

本文中,X為PMN分布模型參數(shù),即

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

(29)

本文將GWO算法應(yīng)用到PMN分布模型參數(shù)的求解當(dāng)中,通過對4個(gè)模型參數(shù)的快速尋優(yōu),以達(dá)到提高模型求解準(zhǔn)確性和高效性的目的。基于GWO算法的模型參數(shù)求解流程如圖1所示。具體實(shí)現(xiàn)過程包括以下步驟:

圖1 基于GWO的模型參數(shù)求解流程

1)對PMN序列進(jìn)行高度角和數(shù)值區(qū)間劃分,統(tǒng)計(jì)各區(qū)間PMN樣本頻數(shù),建立實(shí)際CDF矩陣。

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

3)計(jì)算狼群個(gè)體位置及其適應(yīng)度。

4)遍歷灰狼種群,計(jì)算所有個(gè)體適應(yīng)度函數(shù),將適應(yīng)度函數(shù)值排序前3位的灰狼個(gè)體標(biāo)記為α,β,δ。

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

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

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

4 仿真分析

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

4.1 高度角閾值確定

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

圖2 載波噪聲Q-Q Plot

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

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

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

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

4.2 模型參數(shù)求解

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

圖4 模型參數(shù)估計(jì)結(jié)果

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

表2 PMN模型參數(shù)

4.3 觀測噪聲包絡(luò)分析

為驗(yàn)證所求模型參數(shù)的有效性,采集一組實(shí)測數(shù)據(jù)并提取三頻非差組合PMN序列,觀測數(shù)據(jù)信息見表3。分別利用經(jīng)驗(yàn)?zāi)Ｐ秃蚉MN分布模型對載波噪聲進(jìn)行3σ包絡(luò)分析,結(jié)果如圖5所示。

圖5 兩種模型噪聲包絡(luò)情況

表3 觀測數(shù)據(jù)信息

提取所有衛(wèi)星在5°～20°和45°～60°高度角區(qū)間中的三頻組合PMN序列,對比經(jīng)驗(yàn)?zāi)Ｐ秃蚉MN分布模型在低高度角和大高度角兩種情況下對載波噪聲概率密度的包絡(luò)情況,結(jié)果如圖6和圖7所示。

圖6 低高度角噪聲模型概率密度包絡(luò)

圖7 大高度角噪聲模型概率密度包絡(luò)

結(jié)合圖5、圖6和圖7結(jié)果可以發(fā)現(xiàn),兩種模型均能對載波噪聲進(jìn)行包絡(luò),但是經(jīng)驗(yàn)?zāi)Ｐ蛯MN的噪聲分布情況估計(jì)明顯過于保守,與經(jīng)驗(yàn)?zāi)Ｐ拖啾?PMN分布模型的包絡(luò)曲線與真實(shí)PMN貼合更加緊密,對載波噪聲包絡(luò)效果更好。

4.4 保護(hù)級(jí)結(jié)果評(píng)估

為進(jìn)一步驗(yàn)證所提載波噪聲模型在PPP完好性監(jiān)測中的性能,選取機(jī)載動(dòng)態(tài)觀測數(shù)據(jù)進(jìn)行保護(hù)級(jí)計(jì)算。觀測時(shí)長約為30分鐘,采樣率為1 s,接收機(jī)類型為“NovAtel OEM7”,衛(wèi)星天空圖和飛行軌跡分別如圖8和圖9所示。

圖8 衛(wèi)星天空圖

圖9 載體飛行軌跡

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

表4 完好性參數(shù)設(shè)置

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

圖10 垂向保護(hù)級(jí)和定位誤差

圖11 水平保護(hù)級(jí)和定位誤差

兩種保護(hù)級(jí)收斂后的統(tǒng)計(jì)信息見表5。與經(jīng)驗(yàn)噪聲模型相比,PMN分布模型的VPL平均降低了38.1%,HPL平均降低了27.1%,說明后者模型更加精確,克服了經(jīng)驗(yàn)?zāi)Ｐ捅Ｗo(hù)級(jí)計(jì)算結(jié)果偏保守的缺點(diǎn)。

表5 保護(hù)級(jí)的統(tǒng)計(jì)信息

5 結(jié)論

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