劉九龍，曹月玲，胡小工

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

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

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

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

1 偽距測量的多路徑誤差

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

GNSS多頻偽距/相位觀測方程如下:

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

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

(3)

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

(4)

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

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

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

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

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

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

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

2 MHM算法和CNMC算法

2.1 MHM算法

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

下面驗證多路徑誤差的空間重復(fù)性。多路徑誤差的時間重復(fù)性是從空間重復(fù)性推導(dǎo)出來的。MHM網(wǎng)格是基于衛(wèi)星在空間位置的殘差構(gòu)建的。所有的網(wǎng)格在地心坐標(biāo)系中被定義為高度角和方位角。本文分別使用2019年5月21日至23日和2019年5月24日至26日北京監(jiān)測站B1C頻點(diǎn)的偽距殘差構(gòu)建了兩個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中測量方位角從以北向為起點(diǎn)順時針旋轉(zhuǎn)360°，高度角以地平為起點(diǎn)指向天頂，中心表示高度角為90°，最大圓表示高度角為0°。從圖2可以看出，MHM圖有非常相似的空間分布。較大的多路徑誤差分布在5°～25°高度角的環(huán)區(qū)間。

2.2 CNMC算法

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

由式(3)可得：

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

(5)

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

(6)

(7)

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

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

(8)

實時偽距和多路徑改正算法如下。

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

(9)

后續(xù)遞歸為：

(10)

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

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

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

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

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

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

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

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

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

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

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

(d) 39號衛(wèi)星B1I頻點(diǎn)偽距殘差隨高度角變化圖(d) Figure of pseudorange residuals varies with altitude at B1I frequency of SATSCID 39圖4 30和39號衛(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)計

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

圖7為2019年5月25日至28日偽距定位誤差結(jié)果，由圖7可以得出，經(jīng)過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)星偽距測量多路徑誤差，本文利用實測數(shù)據(jù)分析了偽距測量殘差，MEO衛(wèi)星在低高度角，即出入境時，偽距多路徑誤差十分明顯，甚至可以達(dá)到米級。然后利用衛(wèi)星在空域的重復(fù)性建立MHM模型，但是MHM模型對于偽距多路徑誤差的改正并不明顯。同時利用CNMC算法對偽距觀測量進(jìn)行處理。實際數(shù)據(jù)分析表明，該算法在高度角大于30°時可以將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%。