基于小波分解與恒星日濾波的BDS多路徑削減方法

楊威，張秋昭，馬曉東

（1.中國礦業大學 環境與測繪學院，江蘇 徐州 221116；2.中鐵大橋勘測設計院集團有限公司，湖北 武漢 430051）

0 引 言

隨著科學技術與工程技術發展，國內大型復雜結構建筑數量不斷增加，傳統變形監測方法局限嚴重，應用全球衛星導航系統（global navigation satellite system，GNSS）觀測技術進行變形監測已成為保證建筑結構健康安全的重要渠道，然而，受軌道誤差、電離層延遲、對流層延遲、相位中心偏移與多路徑誤差等影響，GNSS載波相位測量技術精度一般為厘米級，要達到變形監測的毫米、亞毫米級，必須嚴格消除誤差項。短基線（＜3 km）情況下，通過載波相位雙差技術可基本消除多路徑誤差以外的誤差項。

國內外學者提出了多種方法降低測站周圍環境、衛星測站幾何結構對多路徑誤差的影響。除改進天線和硬件外，還改進了數據后處理方法，常用的數據后處理方法包括信噪比法［1］和恒星日濾波法［2］。然而，信噪比與多路徑效應影響之間存在時間差且無法量化［3］，具有一定局限性。文獻［4］發現在測站恒定不動且周圍環境不變時，衛星與地面測站及周圍反射面的幾何構型具有周期重復性。J.F.Genrich等［5］依據這一特征首先提出恒星日濾波法（sidereal filtering，SF），由此逐漸發展出基于觀測值域和坐標域的SF［6-8］。利用坐標域SF對首次觀測解算坐標提取多路徑誤差序列，并以此對重復觀測定位結果進行修正，容易實現且計算簡單。但該方法采用平均星下點周期，忽略了衛星間的差異性，無法有效應用于具有3種不同類型軌道的北斗衛星導航系統（Beidou satellites navigation system，BDS）多路徑緩解中［9］。觀測值域SF通過對首次觀測值提取單差殘差或雙差殘差，再對每顆衛星殘差進行濾波處理得到各衛星的多路徑誤差序列，然后對每顆衛星觀測值域進行改正［10］。該方法考慮了每顆衛星重復周期不完全相同的特點，適應性強，計算相對復雜。此外，在應用SF提取多路徑誤差序列時，還要采用合理的濾波方法，常用的濾波法有小波濾波［11］、EMD濾波［12］、Vondrak濾波［13］等。

由于BDS獨特的星座設計，眾多學者發現BDS衛星的多路徑效應特征顯著不同于其他系統的，因而進行了大量研究。WU X L等［14］提出一種緩解IGSO衛星和GEO衛星偽距多路徑的方法；MA X等［15］發現BDS系統GEO衛星觀測值的多路徑效應呈低頻變化特征，而IGSO衛星和MEO衛星主要呈高頻變化特征；YE S等［9］研究了BDS衛星多路徑重復周期，并提出了BDS系統的載波相位多路徑消除方法；DONG D N等［16］提出了多路徑半天球圖（multipath hemispherical map，MHM），實現了多路徑效應的實時在線解算和改正；DAI W等［17］對MHM算法和恒星日濾波算法進行了對比分析，提出了多路徑誤差參數化改進模型。

本文針對BDS系統星座特性，提出一種恒星日濾波與小波結合的多路徑削減方法，該方法基于重構后的載波相位和偽距單差殘差提取多路徑誤差序列，對每顆BDS衛星分別進行多路徑誤差改正，采用單頻單歷元LAMBDA算法解算基線處理實測數據。結果表明該方法能有效削弱多路徑誤差，提高定位精度。

1 多路徑誤差理論

多路徑效應引起的路程差［10］可表示為

式中：D為接收機天線與反射體間的水平距離；θ為反射信號的入射角；λ為信號波長；α為測站周圍反射體的反射系數；Δφm為多路徑效應引起的相位延遲量。

由式（1）可知，多路徑誤差的大小受反射體反射系數、反射信號入射角和反射體到接收機天線的距離等因素影響。靜態基線下，反射系數與距離可認為不變，入射角則與衛星相對地面的運行軌跡密切相關。而衛星相對地面的運行軌跡是周期重復的，所以多路徑誤差也應隨衛星運行軌跡周期變化。

2 BDS衛星星下點重復周期

由于衛星多路徑誤差變化具有周期性，研究中可利用首次觀測提取出的多路徑誤差序列改正重復觀測數據，衛星的多路徑周期應與其星下點周期一致，因此必須計算出每顆衛星的星下點周期。對于BDS系統的3種衛星，GEO衛星和IGSO衛星的星下點軌跡周期約為1 d，而MEO衛星的星下點軌跡周期約為7 d。由于周期近似于整數天，需要求出每顆衛星重復周期的時間提前量。根據廣播星歷給出的軌道半徑與角速度改正，求出該時刻的衛星運行角速度

對于GEO衛星和IGSO衛星，時間提前量

對于MEO衛星，時間提前量

以2017年1月和2月的BDS星歷數據為基礎，從中提取所有歷元時刻的相關軌道參數，參照式（3）～（4），求出所有觀測到的衛星在各個星歷時刻所對應星下點周期的時間提前量，見圖1～3。

圖1 GEO衛星時間提前量Fig.1 Time advance of GEO satellite

圖2 IGSO衛星時間提前量Fig.2 Time advance of IGSO satellite

圖3 MEO衛星時間提前量Fig.3 Time advance of MEO satellite

由圖1～3可知，3種不同類型的衛星對應的平均時間提前量各不相同，同一類型的衛星對應的時間提前量也不一樣，即便是同一顆衛星在不同時間段的時間提前量也不同，所以修正多路徑誤差不僅要精確到單顆衛星，還要精確到每顆衛星的對應時刻。

3 單差殘差重構

短基線雙差偽距和載波相位觀測方程為

設靜態模型下可觀測衛星顆數為m，則單頻單歷元可得到2（m-1）個雙差方程，將近似值歸入觀測值并線性化可得

其中，x和y分別為雙差模糊度和雙差基線向量，A和B分別為二者的系數陣，l為偽距和載波組合雙差觀測值，dd為雙差殘差項。

短基線雙差基本可以消除鐘差、電離層延遲、對流層延遲和星歷誤差，所以雙差殘差項殘余主要為多路徑誤差。由于雙差過程中基準星為最大高度角衛星，并非恒定衛星，造成部分衛星的部分歷元雙差殘差項缺失，可以通過反求單差殘差規避。由單差到雙差的轉換過程為

其中，sd為單差殘差，D為轉換矩陣，且為奇異矩陣，無法求逆。

為解決D無法求逆這一問題，采用C.Albe等［18］在2000年提出零均值假設，即，使矩陣D變為非奇異矩陣，則式（7）可以寫為

其中，wi為控制條件方程質量的關鍵因子，即觀測值的權，該值受觀測值自身質量影響。本文采用高度角定權的方式確定權陣大小，即wi＝sin2（θ）。則單差殘差可重構為

4 單頻單歷元BDS多路徑的SF算法

考慮到BDS衛星重復周期的差異性，應基于觀測值域進行恒星日濾波。對于單頻單歷元載波相位解算，觀測方程必須包括偽距方程，否則會造成秩虧，因此應同時對載波和偽距方程進行去多路徑處理。其主要步驟如下。

（1）以解算結果的均值x作為基線向量的已知解，將其帶入第一次觀測數據解算，反算整周模糊度并取整得y，將x和y回帶到式（2）得到雙差載波和偽距觀測值殘差v。

（2）利用式（6）將雙差觀測值殘差轉換為單差觀測值殘差。

（3）使用小波濾波方法分解單差殘差項，提取每顆衛星的載波和偽距多路徑誤差序列。

（4）基于第一次觀測數據的廣播星歷計算每顆衛星的時間提前量，根據時間提前量調整單差殘差多路徑誤差時間序列并轉化為雙差序列。

（5）從第二次觀測數據雙差方程的殘差項中減去雙差殘差多路徑序列，用LAMBDA算法單頻單歷元解算基線。

圖4 單差殘差恒星日濾波算法流程Fig.4 Flow of single difference residual sidereal filter algorithm

5 多路徑誤差改進試驗

采用一組實測單頻數據驗證本方法，數據采集時間為2016年9月9日至26日，每日采集4 h BDS數據，采樣頻率為1 s，地點為江蘇省徐州市中國礦業大學環測樓樓頂，接收機型號為天寶R10，主站安置于無遮擋無多路徑環境并開啟接收機多路徑抑制功能，移動站周圍模擬多路徑環境，并關閉接收機多路徑抑制功能。以其中17，23和24日的數據為基礎進行改正和展示。

5.1 計算重復觀測提前量

觀測時期C01～C05為GEO衛星，C06～C10，C15為IGSO衛星，C11，C12和C14為MEO衛星。根據DBS衛星星下點重復周期可知，9月24日1～10號和15號衛星的時間提前量應為9月23日對應時間段的廣播星歷計算的平均值，9月24日11，12，14號衛星的時間提前量應為9月17日對應時間段的廣播星歷計算的平均值，其計算結果如表1所示。

表1 9月24日BDS衛星重復觀測提前量Tab.1 Advances of repeated observation of BDSsatellite on September 24

5.2 單差殘差

采用單差殘差重構步驟，提取17，23和24日每顆BDS衛星的單差殘差序列，受篇幅所限圖5～7僅給出C01，C06和C14衛星的載波相位單差殘差，圖8為C14衛星的偽距單差殘差，采用相關系數法分別求取每顆衛星17與24日、23與24日的單差殘差序列的相關性與錯位時間的關系曲線（即遍歷計算兩個序列錯位［-10 000～1 000］s的相關系數值），其最大相關值點為重復觀測提前量。表2給出示例衛星分別在17—24日和23—24日多路徑誤差曲線的最大相關值及其對應的重復觀測時間提前量。

表2 C01，C06，C14衛星單差殘差序列相關性與時間提前量Tab.2 Correlations and advances of C01，C06，C14 single difference residual sequence

圖5 1號GEO衛星載波單差殘差序列及相關系數Fig.5 Carrier phase single difference residual sequence and correlation coefficient of GEO 01

圖8 14號MEO衛星偽距單差殘差序列及相關性Fig.8 Code single difference residual sequence and correlation coefficient of MEO 14

根據BDS衛星星下點重復周期可知，GEO衛星和IGSO衛星的星下點軌跡周期均為近1 d，故C01和C06衛星在23—24日的單差殘差序列中表現出良好的相關性。C01衛星載波和偽距最大相關系數分別為0.712 7和0.755 7，最大相關處時間提前量分別為224和218 s，與通過星歷計算的226 s十分接近。同時C06衛星相位和偽距最大相關系數分別為0.832 2和0.947 4，最大相關處時間提前量分別為242和255 s，與244 s理論值相符。MEO衛星的星下點軌跡周期均約7 d，故C14衛星在17日和24日的單差殘差序列中表現出良好的相關性，最大相關系數分別為0.839和0.822 9，最大相關處時間提前量為分別為1 701和1 702 s，與理論值1 704 s相符；23—24日的單差殘差序列則并沒有明顯的同趨性，最大相關系數僅為0.401。

圖6 6號IGSO衛星載波單差殘差序列及相關性Fig.6 Carrier phase single difference residual sequence and correlation coefficient of IGSO 06

圖7 14號MEO衛星載波單差殘差序列及相關性Fig.7 Carrier phase single difference residual sequence and correlation coefficient of MEO 14

5.3 sym6小波提取多路徑誤差

小波濾波器通過伸縮平移變換將信號分解成不同的頻率層［19］。文獻［20］指出多路徑誤差主要存在于低頻層，而高頻層主要包括白噪聲或建筑物的高頻振動信號，因此需去掉高頻部分而保留低頻部分，重構原始單差殘差序列。對比不同小波基函數對模擬噪聲的去噪效果，選取sym6小波提取多路徑誤差序列，分解級數為4層，選取兩層低頻的進行重構。

采用sym6小波方法對17日MEO衛星和23日GEO衛星、IGSO衛星的載波和偽距單差殘差序列進行分解，將提取的低頻部分視為多路徑誤差。圖9為C01，C06和C14衛星的單差殘差多路徑時間序列。

圖9 單差殘差提取的多路徑誤差序列Fig.9 Multipath error sequences of single difference residual extraction

5.4 改正多路徑誤差

根據時間提前量調整17日的單差載波和偽距多路徑誤差序列，使之與24日觀測數據的多路徑誤差在時間上保持一致。將提取的單差多路徑誤差轉化為雙差形式后帶入雙差方程求解，其轉換方程為式（8），相當于雙差轉單差的逆變換。圖10為使用本文方法改進多路徑誤差前后的E，N，U偏差圖，該圖以基線解算結果的平均值為參照解，以誤差值的標準差為評價標準，表3為E，N，U改進效果。

圖10 多路徑誤差改進效果Fig.10 Improvement effects of multipath error mitigation

表3 基線解算E，N，U偏差改進效果Tab.3 Improvement effect of baseline solution deviation of E，N，U directions

如圖10所示，左圖為改進前的多路徑誤差基線E，N，U偏差圖，右圖為改進后的偏差圖，改進后的偏差平滑效果明顯。由表3可知，改進前的E，N，U偏差分別為2.87，3.59和8.05 mm，改進后相應削弱為1.72，2.03和5.00 mm，精度相應提升了40.1%，43.5%和37.9%。

6 結 語

本文系統闡釋了BDS系統3種衛星的星下點重復周期特性，針對由此造成的多路徑誤差特性，提出了一種應用于BDS衛星多路徑削弱的數據后處理方法。該方法基于觀測值域重構每顆衛星的相位和偽距單差殘差序列，應用sym6小波分解，從中提取多路徑誤差序列，并將其按提前量調整序列后改正到重復觀測數據的相位和偽距方程中，然后使用LAMBDA方法單頻單歷元解算基線向量。最后，對一組實測數據實施上述方法，結果表明，該方法能有效削弱多路徑誤差，提高了定位精度，使E，N和U精度分別提升40.1%，43.5%和37.9%。