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

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

1 引言

近年來(lái),隨著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)測(cè)(Receiver Autonomous Integrity Monitoring,RAIM)作為提高定位可靠性的一種重要手段,其在PPP中的應(yīng)用也受到越來(lái)越多的關(guān)注[2],[3]。

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

鑒于此,本文利用三頻消電離層組合觀測(cè)值提取三頻非差組合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)組合觀測(cè)值消除電離層延遲誤差,并利用擴(kuò)展卡爾曼濾波進(jìn)行定位解算[16],狀態(tài)方程和觀測(cè)方程分別為

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ī)游走過(guò)程描述。假定各頻段上碼觀測(cè)值和相位觀測(cè)值噪聲水平分別一致,則IF組合觀測(cè)噪聲方差與單個(gè)頻段的噪聲方差關(guān)系為

(3)

(4)

(5)

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

(6)

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

PPP定位解算的EKF方程為

(7)

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

(8)

式中,PL0為無(wú)故障情況下的保護(hù)級(jí);PL1為有故障情況下的保護(hù)級(jí);k0和k1為放大系數(shù),由完好性風(fēng)險(xiǎn)和連續(xù)性風(fēng)險(xiǎn)確定;σp為定位誤差標(biāo)準(zhǔn)差,由Pk主對(duì)角線元素得到;λ為非中心化參數(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ì)算過(guò)程中,特征斜率和定位誤差標(biāo)準(zhǔn)差均直接或間接與Rk相關(guān),因此載波噪聲模型的精確與否直接影響完好性監(jiān)測(cè)算法的可用性。

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

3.1 載波觀測(cè)噪聲序列提取

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

(11)

(12)

(13)

(14)

式中,f5為L(zhǎng)5載波頻率。

在無(wú)周跳情況下,BDIF為常值。通過(guò)分段區(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),通過(guò)求解三頻非差組合PMN序列的噪聲分布模型可以間接獲得IF組合的觀測(cè)噪聲模型,從而實(shí)現(xiàn)保護(hù)級(jí)計(jì)算。

3.2 三頻組合PMN分布模型

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

(17)

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

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

(18)

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

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

本節(jié)主要在該模型基礎(chǔ)上,利用三頻組合PMN序列,求解針對(duì)載波噪聲的模型參數(shù)。與偽距噪聲類似,三頻組合PMN由接收機(jī)噪聲和多路徑誤差組成,其中接收機(jī)噪聲與高度角無(wú)關(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(·)表示對(duì)應(yīng)高度角區(qū)間中PMN的統(tǒng)計(jì)數(shù)。

所有高度角區(qū)間對(duì)應(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)化求解問(wèn)題為

(25)

滿足式(25)的σ1,σ2,a,η即為理論上的模型參數(shù)組合。該目標(biāo)函數(shù)中含有正態(tài)分布函數(shù),該問(wèn)題是典型的多維非線性優(yōu)化問(wèn)題,采用一般線性優(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)化算法,通過(guò)模擬灰狼群體在捕食過(guò)程中的覓食行為實(shí)現(xiàn)目標(biāo)優(yōu)化。在GWO算法中,狼群由高到低排序?yàn)棣?β,δ,ω4個(gè)等級(jí),其中α,β,δ灰狼等級(jí)最高,其位置向量表示3組帕累托最優(yōu)解,ω灰狼位置根據(jù)α,β,δ更新,狼群通過(guò)等級(jí)更替實(shí)現(xiàn)最優(yōu)解的搜索,更新過(guò)程如下

(26)

(27)

式中,t為當(dāng)前迭代次數(shù);Xα,Xβ,Xδ分別為α,β,δ狼個(gè)體位置向量;X表示ω狼個(gè)體的位置向量;D表示狼個(gè)體與最優(yōu)解的距離;C和H為權(quán)重系數(shù)。在迭代過(guò)程中,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)中,通過(guò)對(duì)4個(gè)模型參數(shù)的快速尋優(yōu),以達(dá)到提高模型求解準(zhǔn)確性和高效性的目的?；贕WO算法的模型參數(shù)求解流程如圖1所示。具體實(shí)現(xiàn)過(guò)程包括以下步驟:

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

1)對(duì)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日)的觀測(cè)數(shù)據(jù),接收機(jī)型號(hào)為Trimble NetR9,天線為“TRM55971.00”,采樣率為30s。由于GPS衛(wèi)星目前只有Block ⅡF類型衛(wèi)星具有三頻觀測(cè)數(shù)據(jù),因此只提取該類型衛(wèi)星的PMN序列,觀測(cè)數(shù)據(jù)總量為1.04×106,可以用于高度角閾值和模型參數(shù)的確定。

4.1 高度角閾值確定

PMN在不同高度角范圍具有不同的噪聲水平,利用Q-Q Plot方法、峰度和偏度系數(shù)對(duì)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序列的對(duì)稱性和“厚尾”情況,分別計(jì)算衛(wèi)星在不同高度角區(qū)間的偏度值和峰度值,結(jié)果如圖3所示?？梢钥吹?各高度角區(qū)間中載波噪聲對(duì)應(yīng)的偏度始終在零值附近波動(dòng),說(shuō)明PMN序列的分布具有對(duì)稱性,高度角的變化對(duì)PMN偏度影響較小。在0～20°高度角區(qū)間內(nèi),PMN峰度值顯著大于3,PMN序列的分布具有尖頂和“厚尾”特征;在大于20°高度角區(qū)間中,PMN峰度值趨于一致,近似為3,說(shuō)明基本服從高斯分布。

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

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

4.2 模型參數(shù)求解

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

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

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

表2 PMN模型參數(shù)

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

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

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

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

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

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

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

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

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

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

圖8 衛(wèi)星天空?qǐng)D

圖9 載體飛行軌跡

動(dòng)態(tài)定位解算模式采用IF組合PPP浮點(diǎn)解,衛(wèi)星軌道和鐘差采用IGS改正產(chǎn)品改正,對(duì)流層誤差等采用模型改正,衛(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分布模型對(duì)應(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%,說(shuō)明后者模型更加精確,克服了經(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í)具有保守性的問(wèn)題,本文基于GWO算法求解了一種載波噪聲分布模型。仿真結(jié)果表明,相比于經(jīng)驗(yàn)觀測(cè)噪聲模型,文中所提PMN分布模型對(duì)于載波噪聲的分布描述更為精確,使精密單點(diǎn)定位的垂向保護(hù)級(jí)降低了38.1%,水平保護(hù)級(jí)降低了27.1%,提高了完好性監(jiān)測(cè)的可用性。