北斗三號(hào)系統(tǒng)中軌道地球衛(wèi)星偽距多路徑誤差分析和改正

劉九龍，曹月玲，胡小工

(1. 中國(guó)科學(xué)院上海天文臺(tái) 天文地球動(dòng)力學(xué)中心， 上海 200030；2. 中國(guó)科學(xué)院大學(xué)， 北京 100049； 3. 上海市空間導(dǎo)航與定位技術(shù)重點(diǎn)實(shí)驗(yàn)室， 上海 200030)

隨著北斗三號(hào)系統(tǒng)(BeiDou Navigation Satellite System Phase Ⅲ，BDS-3)逐步完成全球組網(wǎng)并提供導(dǎo)航定位服務(wù)，偽距多路徑已成為影響定位的主要誤差來(lái)源。當(dāng)直射導(dǎo)航信號(hào)到達(dá)接收器天線時(shí)，反射和折射的信號(hào)也同時(shí)到達(dá)接收器天線并污染了直射信號(hào)。這種現(xiàn)象被稱為多徑效應(yīng)，而這些間接信號(hào)引起的測(cè)量誤差就是多徑誤差。針對(duì)固定站或靜止站接收機(jī)常用改正多路徑誤差的方法除了在測(cè)站選址上減少接收機(jī)天線周圍的反射物體外，還可以從基于硬件和基于軟件算法兩方面考慮?；谟布夹g(shù)改善接收機(jī)天線增益模式可以消除部分多路徑誤差的影響。基于軟件算法的方法包括對(duì)多路徑誤差進(jìn)行建模[1]；使用信噪比(Signal-to-Noise Ratio，SNR)對(duì)全球?qū)Ш叫l(wèi)星系統(tǒng)(Global Navigation Satellite System, GNSS)測(cè)量數(shù)據(jù)進(jìn)行加權(quán)處理[2]以及利用各種濾波技術(shù)在頻域、空域或時(shí)域改正多路徑誤差[3]，如卡爾曼濾波[4]、小波分析[5]、恒星日濾波(Sidereal Filter, SF)[6]；奇異譜分析的北斗恒星日濾波算法[7]；利用北斗衛(wèi)星三頻信號(hào)組合[8]；利用多線性擬合和多項(xiàng)式擬合方法[9]等。

多路徑半天球圖(Multipath Hemispherical Map，MHM)模型方法就是利用導(dǎo)航衛(wèi)星在空域的重復(fù)性建立多路徑模型，從而對(duì)多路徑誤差進(jìn)行改正。如利用MHM模型對(duì)GPS短基線偽距和相位測(cè)量中的多路徑誤差進(jìn)行改正[10]，利用MHM模型對(duì)BDS偽距和相位雙差測(cè)量中的多路徑誤差進(jìn)行改正[11]。本文利用雙頻單點(diǎn)偽距和載波相位觀測(cè)數(shù)據(jù)，建立偽距MHM模型，從而對(duì)偽距多路徑誤差進(jìn)行改正。

偽距噪聲和多路徑校正(Code Noise and Multipath Correction，CNMC)法是利用載波相位平滑偽距的方法，對(duì)偽距多路徑誤差進(jìn)行改正。針對(duì)BDS-2的地球靜止軌道(Geostationary Earth Orbit, GEO)和傾斜地球同步軌道(Inclined Geosynchronous Orbit, IGSO)衛(wèi)星，取得了很好的改正效果[12-13]。BDS-3衛(wèi)星在播發(fā)新體制信號(hào)B1C和B2a的基礎(chǔ)上，同時(shí)繼續(xù)播發(fā)平穩(wěn)過(guò)渡信號(hào)B1I和B3I信號(hào)。本文利用CNMC算法針對(duì)BDS-3地球中軌道(Medium Earth Orbit, MEO)衛(wèi)星B1C和B1I信號(hào)的偽距多路徑誤差進(jìn)行改正，利用新體制信號(hào)比對(duì)多路徑誤差改正前后偽距的定位精度。

1 偽距測(cè)量的多路徑誤差

根據(jù)基本導(dǎo)航原理，如式(1)所示，偽距測(cè)量值包含衛(wèi)星軌道信息、衛(wèi)星鐘差、接收機(jī)坐標(biāo)、接收機(jī)鐘差、電離層延遲、中性大氣延遲、多路徑誤差和隨機(jī)噪聲。除了整周模糊度，載波相位測(cè)量和偽距測(cè)量包含相同的衛(wèi)星軌道信息、衛(wèi)星鐘差、接收機(jī)坐標(biāo)、接收機(jī)鐘差和中性大氣延遲。如果忽略高階電離層延遲，對(duì)于相同的觀測(cè)頻率，偽距的電離層延遲與載波相位的電離層延遲大小相等、方向相反。相比于偽距的多路徑誤差，載波相位測(cè)量的多路徑誤差可以忽略不計(jì)[12]。

GNSS多頻偽距/相位觀測(cè)方程如下:

Pi=|Rsat-Rrcv|+Δtrcvlk-Δtsatclk+Δttrop+Δtiono,i+

Δtcor+τsat+τrcv+MPi+εPi

(1)

Li=|Rsat-Rrcv|+Δtrcvlk-Δtsatclk+Δttrop-

Δtiono,i+Δtcor+τsat+τrcv+MLi+λi·Ni+εLi

(2)

其中：Pi為i頻點(diǎn)偽距觀測(cè)量；Li為i頻點(diǎn)相位觀測(cè)量；Rsat是衛(wèi)星位置坐標(biāo)；Rrcv是接收機(jī)位置坐標(biāo)；Δtcor是觀測(cè)誤差改正，包括衛(wèi)星和接收機(jī)的天線相位中心改正、測(cè)站潮汐改正；Δtrcvlk是接收機(jī)鐘差；Δtsatclk是衛(wèi)星鐘差；Δttrop是中性大氣時(shí)延；Δtiono,i是i頻點(diǎn)電離層時(shí)延；τsat和τrcv是衛(wèi)星和接收機(jī)的硬件時(shí)延；Ni是整周模糊度；MPi和εPi分別是偽距測(cè)量多路徑誤差和隨機(jī)誤差；MLi和εLi分別是載波相位測(cè)量多路徑誤差和隨機(jī)誤差。兩者做差可得：

Pi-Li=2Δtiono,i+MPi-λi·Ni

(3)

其中，隨機(jī)誤差包含在偽距多路徑誤差MPi中。對(duì)于雙頻BDS-3信號(hào)，i對(duì)應(yīng)B1C和B2a兩個(gè)頻點(diǎn)信號(hào)。利用電離層延遲與信號(hào)頻率平方成反比的特點(diǎn)可以得到包含偽距殘差的偽距多路徑表達(dá)式：

(4)

其中，f1和f2分別對(duì)應(yīng)B1C和B2a的頻率。式(4)忽略了載波相位觀測(cè)的多路徑誤差，因?yàn)檩d波相位多路徑誤差遠(yuǎn)小于偽距多路徑誤差?？梢酝ㄟ^(guò)對(duì)式(4)取一段時(shí)間窗的平均值作為整周模糊度扣除，同時(shí)由于隨機(jī)誤差符合正態(tài)分布，可以假定其均值為0，即可得到對(duì)應(yīng)頻點(diǎn)的偽距多路徑值。

利用北京監(jiān)測(cè)站2019年5月21日至24日的雙頻偽距/相位觀測(cè)數(shù)據(jù)，計(jì)算偽距殘差，如圖1所示。

(a) 30號(hào)衛(wèi)星B1C頻點(diǎn)偽距多路徑誤差隨時(shí)間變化圖(a) Figure of pseudorange multipath error varies with time at B1C frequency of SATSCID 30

(b) 30號(hào)衛(wèi)星B1C頻點(diǎn)偽距多路徑誤差隨高度角變化圖(b) Figure of pseudorange multipath error varies with altitude at B1C frequency of SATSCID 30

(c) 39號(hào)衛(wèi)星B1C頻點(diǎn)偽距多路徑誤差隨時(shí)間變化圖(c) Figure of pseudorange multipath error varies with time at B1C frequency of SATSCID 39

(d) 39號(hào)衛(wèi)星B1C頻點(diǎn)偽距多路徑誤差隨高度角變化圖(d) Figure of pseudorange multipath error varies with altitude at B1C frequency of SATSCID 39圖1 30號(hào)和39號(hào)衛(wèi)星B1C頻點(diǎn)偽距多路徑結(jié)果Fig.1 Results of pseudorange multipath error at B1C frequency of SATSCID 30 and 39

圖1為30號(hào)和39號(hào)衛(wèi)星B1C頻點(diǎn)偽距多路徑結(jié)果，從圖1可以看出，BDS-3 MEO衛(wèi)星伴隨出入境(即低高度角時(shí))，偽距多路徑誤差會(huì)有大幅增加。當(dāng)高度角小于15°時(shí)，偽距多路徑影響可以達(dá)到米級(jí)，隨著高度角的增加而多路徑影響逐步縮??；當(dāng)高度角大于30°時(shí)，偽距多路徑誤差的影響會(huì)小于0.5 m。由此可見(jiàn)，在偽距測(cè)量中，多路徑誤差的影響十分顯著，這體現(xiàn)了偽距多路徑誤差改正算法的重要性。

2 MHM算法和CNMC算法

2.1 MHM算法

如果衛(wèi)星是來(lái)自相同的GNSS類別(GPS、Galileo或BDS，但是GLONASS除外，因?yàn)槊恳活wGLONASS都有不同的頻率)，并且衛(wèi)星的星上偽距多路徑可以忽略，那么多路徑誤差的影響只和測(cè)站周圍環(huán)境有關(guān)。如果測(cè)站環(huán)境保持不變，那么衛(wèi)星信號(hào)從相同的天空半球面?zhèn)鱽?lái)，無(wú)論信號(hào)來(lái)自哪一顆衛(wèi)星，對(duì)接收機(jī)的多路徑誤差都應(yīng)該是相同的。這種多路徑誤差的空間重復(fù)性是建立MHM模型的基礎(chǔ)。MHM方法通過(guò)累積幾天的偽距殘差來(lái)估計(jì)以高度角和方位角為坐標(biāo)的格網(wǎng)點(diǎn)平均值。然后對(duì)實(shí)時(shí)數(shù)據(jù)進(jìn)行處理，扣除所在格網(wǎng)的偽距多路徑估值。MHM模型是和觀測(cè)頻率相關(guān)的，對(duì)于不同的觀測(cè)頻率，應(yīng)建立不同的MHM模型。

下面驗(yàn)證多路徑誤差的空間重復(fù)性。多路徑誤差的時(shí)間重復(fù)性是從空間重復(fù)性推導(dǎo)出來(lái)的。MHM網(wǎng)格是基于衛(wèi)星在空間位置的殘差構(gòu)建的。所有的網(wǎng)格在地心坐標(biāo)系中被定義為高度角和方位角。本文分別使用2019年5月21日至23日和2019年5月24日至26日北京監(jiān)測(cè)站B1C頻點(diǎn)的偽距殘差構(gòu)建了兩個(gè)1°×1°的MHM網(wǎng)格模型(如圖2所示)。

(a) 2019年5月21日至23日B1C偽距多路徑誤差MHM圖(a) MHM graph of B1C pseudorange multipath error on May 21 to 23, 2019

(b) 2019年5月24日至26日B1C偽距多路徑誤差MHM圖(b) MHM graph of B1C pseudorange multipath error on May 24 to 26, 2019圖2 MHM圖Fig.2 MHM graph

圖2中測(cè)量方位角從以北向?yàn)槠瘘c(diǎn)順時(shí)針旋轉(zhuǎn)360°，高度角以地平為起點(diǎn)指向天頂，中心表示高度角為90°，最大圓表示高度角為0°。從圖2可以看出，MHM圖有非常相似的空間分布。較大的多路徑誤差分布在5°～25°高度角的環(huán)區(qū)間。

2.2 CNMC算法

根據(jù)偽距和載波相位觀測(cè)方程式(1)～(2)，可以容易地驗(yàn)證，同歷元偽距和載波相位觀測(cè)量之差(Code Minus Carrier, CMC)包括雙倍電離層延遲改正、模糊度和通道時(shí)延組合值、偽距多路徑誤差和噪聲信息。根據(jù)雙頻載波相位觀測(cè)量可計(jì)算當(dāng)前歷元(帶模糊度的)電離層延遲改正，如果能估計(jì)出模糊度和通道時(shí)延組合值，就可以分離出偽距多路徑和噪聲。CNMC方法通過(guò)遞推估計(jì)模糊度和通道時(shí)延組合值。具體方法如下：

由式(3)可得：

CMC=2Δtiono,i+MPi-λi·Ni

(5)

利用B1C和B2a雙頻相位觀測(cè)數(shù)據(jù)可以計(jì)算B1C和B2a頻點(diǎn)電離層延遲改正：

(6)

(7)

式中，b包含了載波相位的頻率間偏差和雙頻整周模糊度之差，在后續(xù)的估計(jì)中假定為常數(shù)。偽距和載波相位觀測(cè)量之差扣除雙倍電離層延遲改正后，根據(jù)遞推公式估計(jì)模糊度和通道時(shí)延的組合值，即偽距和載波相位之差(Code Minus Carrier Bias， CMCB)。偽距多路徑誤差(包含隨機(jī)誤差)如下。

MPi=Pi(t1)-Li(t1)-2Δtiono,i-ΔCMCBi

(8)

實(shí)時(shí)偽距和多路徑改正算法如下。

這里假定偽距多路徑誤差和隨機(jī)誤差初始值為0，則初始化

(9)

后續(xù)遞歸為：

(10)

式中，T為遞推時(shí)間間隔，n為遞推個(gè)數(shù)，δtiono,i(t)表示t時(shí)刻雙倍電離層延遲。利用該遞歸公式就可以得到偽距多路徑改正。

3 結(jié)果統(tǒng)計(jì)和分析

利用2019年5月21日至24日數(shù)據(jù)建立MHM模型，改正5月25日至26日的偽距多路徑誤差。同時(shí)與利用CNMC算法改正的偽距殘差值進(jìn)行比對(duì)，結(jié)果如圖3所示。

(a) 30號(hào)衛(wèi)星B1C頻點(diǎn)偽距殘差隨時(shí)間變化圖(a) Figure of pseudorange residuals varies with time at B1C frequency of SATSCID 30

(b) 30號(hào)衛(wèi)星B1C頻點(diǎn)偽距殘差隨高度角變化圖(b) Figure of pseudorange residuals varies with altitude at B1C frequency of SATSCID 30

(c) 39號(hào)衛(wèi)星B1C頻點(diǎn)偽距殘差隨時(shí)間變化圖(c) Figure of pseudorange residuals varies with time at B1C frequency of SATSCID 39

(d) 39號(hào)衛(wèi)星B1C頻點(diǎn)偽距殘差隨高度角變化圖(d) Figure of pseudorange residuals varies with altitude at B1C frequency of SATSCID 39圖3 30和39號(hào)衛(wèi)星B1C頻點(diǎn)偽距多路徑結(jié)果Fig.3 Results of pseudorange multipath error at B1C frequency of SATSCID 30 and 39

圖3為30和39號(hào)衛(wèi)星B1C頻點(diǎn)偽距多路徑結(jié)果。利用同樣的方法，對(duì)B1I頻點(diǎn)的偽距多路徑誤差進(jìn)行計(jì)算和改正。圖4為30號(hào)和39號(hào)衛(wèi)星B1I頻點(diǎn)偽距多路徑結(jié)果。圖5和圖6為偽距殘差的均方根誤差(Root Mean Square, RMS)。

由圖3、 圖5和圖4、圖6可以得出，BDS-3 MEO衛(wèi)星B1C和B1I頻點(diǎn)偽距多路徑誤差在分米量級(jí)，是不可忽視的誤差來(lái)源，其中E為衛(wèi)星高度角。在不進(jìn)行多路徑誤差改正時(shí)，高度角小于30°時(shí)，偽距殘差在幾十厘米。MHM方法對(duì)于偽距多路徑誤差的降低并不顯著，更適合用于短基線高精度測(cè)量。CNMC方法對(duì)于偽距多路徑誤差改正效果十分明顯，在高度角大于15°時(shí)，可以將偽距殘差降低至厘米量級(jí)。B1C和B1I頻點(diǎn)不同高度角偽距殘差RMS統(tǒng)計(jì)結(jié)果分別如表1和表2所示。結(jié)合表1和表2分析，CNMC算法在高度角大于15°時(shí)可以將B1C和B1I頻點(diǎn)偽距多路徑誤差減小60%以上，使偽距精度得到了大幅提升。同時(shí)經(jīng)過(guò)改正的偽距數(shù)據(jù)，其偽距殘差時(shí)間序列變得更加平穩(wěn)。與BDS-2衛(wèi)星B1I頻點(diǎn)30 cm大小的偽距多路徑值相比[14]，BDS-3衛(wèi)星B1I頻點(diǎn)的偽距多路徑誤差減小60%以上。

(a) 30號(hào)衛(wèi)星B1I頻點(diǎn)偽距殘差隨時(shí)間變化圖(a) Figure of pseudorange residuals varies with time at B1I frequency of SATSCID 30

(b) 30號(hào)衛(wèi)星B1I頻點(diǎn)偽距殘差隨高度角變化圖(b) Figure of pseudorange residuals varies with altitude at B1I frequency of SATSCID 30

(c) 39號(hào)衛(wèi)星B1I頻點(diǎn)偽距殘差隨時(shí)間變化圖(c) Figure of pseudorange residuals varies with time at B1I frequency of SATSCID 39

(d) 39號(hào)衛(wèi)星B1I頻點(diǎn)偽距殘差隨高度角變化圖(d) Figure of pseudorange residuals varies with altitude at B1I frequency of SATSCID 39圖4 30和39號(hào)衛(wèi)星B1I頻點(diǎn)偽距多路徑結(jié)果Fig.4 Results of pseudorange multipath error at B1I frequency of SATSCID 30 and 39

(a) 5°

(b) 15°

(c) 30°

(d) 5°

(a) 5°

(b) 15°

(c) 30°

(d) 5°

表1 B1C頻點(diǎn)不同高度角偽距殘差RMS統(tǒng)計(jì)

表2 B1I頻點(diǎn)不同高度角偽距殘差RMS統(tǒng)計(jì)

圖7為2019年5月25日至28日偽距定位誤差結(jié)果，由圖7可以得出，經(jīng)過(guò)CNMC改正后可以使偽距定位精度提高約10%。

(a) 偽距北向定位誤差(a) Pseudorange positioning error in the north direction

(b) 偽距東向定位誤差(b) Pseudorange positioning error in the east direction

(c) 偽距天頂向定位誤差(c) Pseudorange positioning error in the upward direction

(d) 偽距三維定位誤差(d) Pseudorange positioning error in three dimensions圖7 2019年5月25日至28日偽距定位誤差Fig.7 Pseudorange positioning error on May 25 to 28, 2019

4 結(jié)論

為了改正BDS-3 MEO衛(wèi)星偽距測(cè)量多路徑誤差，本文利用實(shí)測(cè)數(shù)據(jù)分析了偽距測(cè)量殘差，MEO衛(wèi)星在低高度角，即出入境時(shí)，偽距多路徑誤差十分明顯，甚至可以達(dá)到米級(jí)。然后利用衛(wèi)星在空域的重復(fù)性建立MHM模型，但是MHM模型對(duì)于偽距多路徑誤差的改正并不明顯。同時(shí)利用CNMC算法對(duì)偽距觀測(cè)量進(jìn)行處理。實(shí)際數(shù)據(jù)分析表明，該算法在高度角大于30°時(shí)可以將B1C頻點(diǎn)多路徑誤差水平從0.11 m降低至0.03 m，使B1I頻點(diǎn)偽距多路徑誤差從0.074 m降低至0.024 m，使MEO衛(wèi)星的B1C和B1I頻點(diǎn)多路徑誤差降低了60%以上，使B1C頻點(diǎn)偽距定位精度提高約10%。