一種基于粗差檢驗整體平差的周跳探測與修復方法

趙祥偉，陳正宇

(中國能源建設集團江蘇省電力設計院有限公司，江蘇 南京 211102)

0 引言

GNSS動態測量中，載波相位觀測值對于定位結果的影響較大，同時是系統性與長時間的影響[1-3]。為了得到高精度的定位結果，必須要準確地修復觀測值中可能存在的周跳。但是，因為GNSS工作的機動性較強，觀測環境變化較大，受到對流程、電離層、多路徑等因素的影響也各不相同，這對載波觀測值的質量造成了多方面的影響[4]。因而，周跳的現象經常發生。如何快速高效地探測修復周跳成為了測繪研究者們急需解決的問題。

隨著研究的深入，對于周跳探測修復工作的研究成果也越來越多。方法更是層出不窮，目前已有的檢測周跳的方法主要有多項式擬合法[5-6]、差分法[7]、多普勒頻移法[3,8]、相位/偽距組合法[9]、電離層殘差法[2]、濾波法[10]和站間或站星間差分法[11]。

本文選取現有的相位/偽距組合法，結合粗差探測[12-14]方法，對周跳探測方法進行適當的改進，通過實例分析得出該方法能夠探測出 2周以上的周跳值，且探測方法簡單，探測速度較快，探測結果可靠。

1 周跳檢驗量構成

1.1 偽距/載波組合法

在周跳探測的諸多方法中，一般常利用模糊度參數進行周跳的探測修復。在信號觀測方程中，偽距觀測方程為[15]：

式(1)中，偽距觀測方程中不包括模糊度項，因而認為偽距觀測值不會受到周跳的影響。相應的載波相位觀測方程可表示為：

綜合式(1)和(2)，為分離出模糊度參數值，通常需要對其進行組合，將相位觀測值直接減去偽距觀測值，消除接收機鐘差和衛星鐘差等公共誤差，忽略次要因素的影響，如多路徑效應、噪聲項等。相位偽距觀測方程作差，表示為：

式中：λ表示載波的波長；?表示載波相位觀測值；P表示偽距觀測值；I表示電離層延遲量；a表示整周模糊度；m?和mP分別表示載波相位觀測值和偽距觀測值的觀測值噪聲；δ?和δP分別表示載波相位觀測值和偽距觀測值殘差。

1.2 構成周跳檢驗量

如式(3)所示，經過相位偽距觀測值作差后，方程里剩余有電離層殘差、多路徑殘差以及觀測噪聲等因素影響。當歷元采樣率較小時，電離層變化值不會太大，而多路徑殘差以及觀測噪聲等影響同樣變化不大，理論上通過歷元間作差可以消除絕大部分的值，假設觀測值中含有周跳，便可以通過殘差序列反映出來。此時，周跳檢驗量D表示為：

式中：t表示對應的檢驗時間；Δt為數據的采樣間隔。

如果信號采樣間隔較短，因而電離層殘差、多路徑效應在時間上相關性較大，歷元間求差后，其對周跳檢驗量的影響較小，可以忽略。由于一般民用用戶不能獲得雙頻P碼觀測值，通常是可獲得L1頻率的C/A碼偽距和L2頻率的P2碼數據，這時可將P2碼偽距代替C/A碼偽距計算檢驗量D，以便獲得相當精度的D檢驗量。

此時，周跳檢驗量D的精度主要受偽距測量誤差的影響。偽距測量的測距精度對于P碼約為29 cm；而載波相位相應的波長為λ1=19.03 cm，λ2=24.42 cm。對 式(4)運 用 誤 差傳播定律得：①對于L1：mD≈1.8(周)；②對于L2：mD≈2.2(周)。

因此在能獲得L2的P碼觀測值的情形下，以三倍檢驗量方差為限，D檢驗量能探測出大于7～8周的周跳。

2 粗差檢驗與整體平差

對于上文構成的周跳探測檢驗量D，其基本服從正態分布。對其中出現的異常值，通常作以下處理：

令Li=λ?(ti)-P(ti),i=1, 2, …,n。n為歷元個數，則相應的檢驗量殘差陣可表示為：

根據式(7)可知，系數矩陣B為(n-1)×n的矩陣，無法根據最小二乘方法獲得方程的唯一解或最優解。

2.1 粗差檢驗

對式(7)進行粗差定位，結合檢驗量構成的誤差方程，利用粗差探測方法，對發生周跳的歷元進行定位，具體方案流程如下。

通 過式(7)中的常 數項l可知，E(l)=E(L)-L0。觀測值的數學期望與近似觀測值的差是一個很小的數，一般規定小于三倍的中誤差，即為觀測值的閉合差。通過對常數項l進行檢驗可以確定發生周跳的具體位置，具體做法為：①首先結合誤差方程計算出常數項l，給出每個常數項l的具體初值。②根據誤差傳播定律，求出l的中誤差，其中L0的中誤差初值根據先驗經驗值給出，以三倍中誤差m為限，給出具體常數項的分類范圍，分為大于m和小于m兩個數據集合。③用小于m的li并集的補和大于m的li的集合分別求交集，假設結果是單個元素，那么可以確定該元素值為粗差，此刻標記上粗差的具體位置i。④假設集合中出現兩個或兩個以上的元素，那么需要重新建立誤差方程式進行甄別，循環求解。

2.2 整體平差修復周跳

此時，假定在觀測的歷元i處發生了周跳，那么從1～(i-1)歷元處的觀測方程系數值均為0，此時，方程呈現出嚴重病態，可以不予考慮，即無周跳的歷元不予處理。在i～n之間的所有歷元均受到周跳的影響，此時需要進行處理得到發生周跳值的大小。B矩陣的劃分如下：

結合式(7)和式(8)可得，最小二乘表達式轉化為：

此時，根據最小二乘原理即可以求得周跳改正的最優解，循環計算直至i=n，此時便可修復觀測值中所有的周跳。

根據上述方法便可對任意觀測文件進行周跳的探測與修復工作，該方法簡單易行，易于編程實現。

3 實例分析

為驗證上述方法的可靠性與實用性，選取2019年7月1日上海佘山站當天的觀測數據。截取2019年7月1日1時06分0秒～2019年7月1日2時13分0秒的22號衛星觀測數據，選取觀測數據，根據式(4)組成檢驗量，其殘差陣序列如圖 1所示。

圖 1 無周跳下檢驗量D殘差序列

圖 1中，由偽距/載波構成的殘差期望值分布在±0.5周左右，故而取三倍中誤差作為檢驗范圍，閾值線如圖1所示。故可以判定該方法在不加入任何其他的改正手段后，其周跳探測范圍為±2周左右。

為驗證本文方法對周跳探測的效果，人為在觀測值中加入周跳。加入方式：在10歷元處加入1周周跳；40歷元處加入5周周跳；80歷元處加入8周周跳，加入周跳后的檢驗量殘差序列圖如圖2所示。

圖 2 加入周跳后檢驗量殘差分布圖

圖 2中，加在10歷元處的1周的周跳無法通過閾值線反映出來，周跳的值與噪聲值相當，淹沒在噪聲中，故而會出現探測的誤差；而加在40歷元處的5周的周跳與80歷元處的8周的周跳在圖中很明顯可以辨別出來，周跳的值在相應的閾值線之上，此時通過閾值可以判別出具體哪個歷元發生了周跳。

在本文方法中，根據殘差陣進行粗差探測，探測出殘差陣中異常值發生的具體位置，探測流程按照2.1節中的步驟進行。通過粗差探測，可以得到只有在i=40和i=80時，用于求交集的兩個誤差方程常數項集合的交集才為單個元素，表示在該地方存在粗差；而i=1處的粗差因為偽距/載波組合法的周跳探測能力限制而無法進行粗差辨別，出現探測 盲區。

當i=40時，進行第一次循環平差計算時，結果見表1所列。

表1 第一次循環平差周跳修復結果

由表1可以看出，在40歷元處，經過第一次整體平差，修正周跳后的載波觀測值與真值相比相差了1周的誤差，這是因為加在10歷元處的1周周跳沒有正確修復的緣故，故而存在 1周的修復誤差，后面的歷元類似；在80歷元處修正后的值與真值差值為9周，其原因是人為加在上面的8周周跳以及1周周跳粗差。第一次平差計算結束后，i指向80歷元處，進行第二次平差計算。

當i=80時，進行第二次循環平差計算時，結果見表2所列。

表2 第二次循環平差周跳修復結果

由表2的平差結果可以看出，經過第二次平差后，加入周跳后觀測值的周跳值得到了修復，40和80歷元處的較大周跳得到了正確修復，而加在10歷元處的周跳因為探測方法的局限性，無法得到正確辨別，因而存在一定的探測盲區。

4 結論

本文對經典方法偽距/載波組合法的原理進行了詳細闡述，并結合粗差探測的方法對周跳進行探測修復。從原理上來說，該方法可行而且易于實現編程。結合實例，驗證了本文方法的優缺點。

通過仿真算例實例驗證得出，本文方法對于2周以上的周跳具有較強的探測修復能力，且探測速度較快，不需要復雜的判斷迭代過程，探測準確，修復精度高；但是對于一些小周跳，由于初始方法的周跳探測能力限制，使其對小周跳的探測無能為力。因此下一步的工作主要是如何利用更為精確的周跳探測方法結合更加優化的粗差探測方法，同時可以考慮多頻觀測條件下的周跳探測方法，進一步優化本文提出的方法，提升該方法的周跳探測適用 范圍。