基于MATLAB的抗差多項式擬合方法在GNSS周跳探測中的應用研究

于佳慧 張世濤 張巖

摘 ?要： 為了提高GNSS載波相位的定位精度，本文采用歷元間求差與抗差估計相結合的方法對過去常常使用的多項式進行改良，最終提出用抗差多項式擬合來探測并修復周跳。本文基于Matlab軟件編制程序，并對GNSS載波相位周跳進行了試算。實驗結果表明：此方法可以做到對GNSS周跳探測與修復功能，具有文件讀寫、衛星觀測值提取、周跳探測與修復、繪制擬合的誤差結果圖形等功能。通過這種方法可以探測并修復1周以上的周跳。得出結論：通過抗差多項式擬合方法可以探測并修復GNSS周跳。

關鍵詞： MATLAB;GNSS;周跳的探測與修復;抗差多項式

中圖分類號： TP3 ? ?文獻標識碼： A ? ?DOI：10.3969/j.issn.1003-6970.2019.08.041

本文著錄格式：于佳慧，張世濤，張巖. 基于MATLAB的抗差多項式擬合方法在GNSS周跳探測中的應用研究[J]. 軟件，2019，40（8）：175180

【Abstract】： In order to improve the positioning accuracy of GNSS carrier phase， this paper uses the method of inter-calendar difference and robust estimation to improve the traditional multinomial， and finally proposes a method of using robust multinomial fitting to detect and repair cycle slip. In this paper， based on Matlab software， the phase cycle jump of GNSS carrier is calculated. The experimental results show that this method can detect and repair GNSS cycle slip， and has the functions of file reading and writing， satellite observation data extraction， cycle slip detection and repair， drawing fitting error result figure and so on. Through this program， the cycle jump can be detected and repaired for more than 1 week. Draw Conclusion： the robust polynomial fitting method can be used in the detection and repair of GNSS cycle slip.

【Key words】： MATLAB; GNSS; Detection and repair; Robust Polynomials

0 ?引言

GNSS具有使用觀測精度高、測量時間短、全天候、技術自動化程度高、觀測范圍大、能提供三維坐標測量成果等特點[1]。在GPS數據處理過程中，周跳的存在會使觀測值中出現一個偏差，這會使觀測值失真，若不修復周跳，剔除這些錯誤成果，會對測繪工作造成嚴重后果。因此，要想提高定位精度就必須能準確地探測周跳并且修復載波相位。目前可用于接收機的周跳探測方法主要有高次差法、電離層求差法、偽距和載波相位觀測值組合、多項式擬合法來探測和修復周跳等。

總之，有很多方法可以探測和修復周跳，但不同的周跳探測與修復方法對提高計算效率及可靠性都有所不同。常規方法對載波相位觀測值周跳探測和修復雖然都有其優點，但也存在其局限性，鑒于此，本文將以多項式擬合法為基礎，先進行時間標準化，然后采用歷元間求差與抗差估計相結合的方法對傳統算法中存在的問題進行改進，并利用GNSS/GLONASS實測數據驗證改進后的算法的有效性[2]。

1 ?MATLAB簡介

1.1 ?什么是MATLAB

MATLAB軟件是由美國Math Works公司推出的一種交互式、面向對象的程序設計語言，廣泛用于工業界和學術界。不僅能很好的處理數值問題，同時還能解決符號問題，繪制圖形。

1.2 ?MATLAB所具有的優勢

MATLAB系統包括以下五個部分：

（1）MATLAB具有多種語言體系。在MATLAB可用的編程語言，對應于工具箱中六個目錄，分別為ops（操作符與特殊字符）、lang（編程語言結構）、strfun（字符串操作）、iofun（文件的輸入輸出）、time fun（時間和日期）和data type（數據類型和結構）[3]。

（2） MATLAB具有完整工作界面。包括MATLAB編程和調試的環境，對工作空間中的變量進行管理及進行輸入輸出的數據[4]。

（3）MATLAB具有強大的圖形處理能力。其中包括二、三維圖形可視化，也包括定義圖形外觀的操作，最終可以為MATLAB應用程序創建整套的用戶界面。

（4）數學函數庫。包括初等函數;也包括更復雜的功能，如矩陣功能（矩陣求逆、求矩陣的特征值）、快速傅里葉變換等[3]。

（5）MATLAB應用程序接口（API）。提供接口程序使程序員可以編寫C/C6、Fortran程序調用MATLAB作為計算引擎，并用MATLAB調用外部應用程序或者讀寫MATLAB文件[4]。MATLAB中的M文件的語法與其他的高級語言類似，但語法較其他一般的高級語言來說要簡單，它只是一個ASCII碼簡單的文本文件，其語法編制與數學語言非常接近。它不僅是一種程序化的編程語言，也是一種解釋性的編程語言，可逐行解釋和運行程序，程序更易調試[4]。

2 ?抗差多項式擬合法

2.1 ?周跳概念和產生的原因

載波相位可以用在高精度衛星導航定位。完整且連續的周期通過使用多普勒計數完成，而載波相位在導航定位中主要是測量不完整的周期。信噪比低、信號質量差以及接收機接收到的質量差等方面都可能引起完整的周期部分突變，叫作整周跳變，簡稱為周跳（cycle slip）。

載波相位是由 、 、 這些部分組成的。雖然可以以極高的精度測定 ，但也要正確地測定 和 的情況下，才能確定載波相位值。在載波相位測量的實際觀測值中，小于一周的部分稱為觀測時刻的觀測值。衛星載波信號與接收機參考振蕩信號形成的差頻信號不足一周。接收機載波跟蹤回路中的鑒相器測量可以實現對載波信號的測量。只有接收機的載波信號和參考振蕩信號能在觀測時段產生差頻信號，才可以得到正確的觀測值。整個星期的計數不一樣。它是計數器從第一個觀測時間到當前觀測時間積累的差分頻率信號中的整數頻帶數。如果由于外界因素，通過使用多普勒計數法的計數器在兩個觀測周期之間的某個時段停止計數，從而使整個周期的計數比實際值少n周，那么當計數器恢復正常運行時，所有載波相位波段將包含相同的誤差值，該誤差值小于整個周期的正常值（如圖1）。這樣一個系統性的誤差，在一個周期內計數和不足一周的部分保持正確計數的情況，稱為整周跳變，簡稱周跳[5]。

2.2 ?抗差估計求解擬合系數

傳統的多項式擬合方法一般計算擬合系數是采用最小二乘法，但當觀測值不服從正態分布的假設時，即m個擬合觀測值出現粗差或周跳時，最小二乘法估計就不具有抗干擾性，錯誤的估計結果可能由單個觀測值的偏差導致。而利用抗差估計法與利用最小二乘估計法不同，它不像最小二乘那樣過分的追求有效性與無偏性等內部性質，抗差估計則更加重視估值的實際可靠性和抗差性[6]。從第一個歷元開始，先一次取m個歷元 和對應的載波相位觀測值 代入式（1），通過抗差估計法解算其擬合系數。假設初始權矩陣為單位陣，即 ，之后進行平差計算。根據第一次平差后得到的改正數陣 和單位權中誤差 ，利用下式（2）

對各相位觀測值重新定權，如果 則認為m個歷元觀測值中不存在周跳，最后的平差結果與第一次平差結果相同。若有偏差，可利用權陣 在再次進行平差計算，直至 時停止迭代，取第j次迭代結果作為最終結果。這樣可以確保參加擬合的m個相位觀測值為不含周跳的值，從而確保得到的擬合系數更加準確和可靠。

2.3 ?利用抗差多項式探測和修復周跳

根據求得的擬合系數擬合，并帶入m+1個歷元值，來估計第m+1個歷元的相位值，將利用抗差多項式擬合出來的擬合值與對應歷元的實際觀測作差，若差值 時，則認為 不含周跳，去掉本次擬合的第一個觀測值，加入第m+1個歷元觀測值繼續進行擬合;當差值 時，則認為 含有周跳，采用外推的整周計數去取代有周跳的實際觀測值中的整周計數，但不足一周的部分仍然保持不動，然后繼續上述過程，直到最后一個相位觀測值為止[7]。

本部分首先對高次差法的原理和適用范圍優缺點進行了介紹，其次介紹了利用傳統多項式探測周跳原理，在此基礎上進行時間標準化、歷元間求差、重新定權等改進，從而建立抗差多項式模型，利用抗差多項式探測和修復周跳。

3 ?編制周跳探測與修復系統程序

3.1 ?程序的設計框架

研究利用多項式[8]擬合法周跳探測和修復方法，對傳統的多項式擬合法周跳探測和修復方法進行改進，并確定多項式擬合的階數。

在無周跳“干凈”觀測數據中加入周跳，通過差分技術對相鄰歷元之間的觀測數據進行差分，對差分數據進行多項式建模擬合與預報，探測周跳;利用抗差多項式建模探測周跳，即編制抗差多項式程序設計框架（如圖2）。

3.2 ?周跳探測程序算例分析

1. 解決多項式的基礎問題

衛星到地面的距離對時間的四階以上導數已趨于零，其變化規律是隨機的，無法再用多項式進行擬合。故多項式擬合的階數n取3～4階，在本實驗中n=3的條件下既能保證精度又保證了運算速度。

擬合窗口m選取時，當擬合窗口寬度越大時，外推值越準確，擬合中誤差越小，但過大又會增加計算量。當擬合窗口寬度越小時，外推值精度越差。所以通過取不同的值進行實驗，根據實驗取m=10比較合適。

門檻系數k的取值是隨不同的情況改變的。當k較大時，表示只有當觀測值偏離外推值很大的時候才認為它是異常的;當k較小時，表示當觀測值偏離外推值較小的時候就認為它是異常值了。一般取k在2到9之間的數。經過多次取值實驗，結果表明門檻系數取k=4時效果最好。

2. 高次差法檢驗觀測數據質量

算例基于2014年1月14日國內的GPS/GLONASS的實測數據，觀測時間為10分鐘，采樣間隔為1秒。使用高階差分的方法先對16時45分0秒到16時55分0秒GPS的10號衛星共600個L1載波的據進行周跳探測，對觀測值求四階差分的時間序列圖（如圖3），從（如圖3）可以看出，實驗結果中沒有周跳[9]。

3. 多項式擬合法算例分析

本文選取上述無周跳的數據中從16時45分0秒到16時45分60秒的實測GPS數據進行多項式擬合實驗，以第3秒（即第4個歷元）到第12秒（即第13個歷元）的數據作為擬合多項式的數據，外推出從第13秒到第60秒的數據，并設計實驗方案：方案1，傳統的多項式擬合法;方案2，歷元間求差并采用最小二乘[10]平差的多項式擬合法;方案3，歷元間求差并采用抗差估計的多項式擬合法。

在無周跳情況下（如圖4），三種方案擬合能力的比較，前一段時間曲線變化比較平緩，較符合實際情況;隨著觀測歷元數的增大，方案一與方案二曲線的波動越來越大，出現發散的情況，有可能造成周跳判斷失誤。方案三可明顯看出采用抗差估計后每個歷元的擬合中誤差都有較大改善，GPS觀測量擬合值與實際值之差在0.01周以內。

人為在第6秒（即第7個觀測歷元）上加上1周周跳（如圖5）三種方案擬合效果的比較，前一段時間3條曲線變化比較平緩;隨著觀測歷元數的增大，方案一與方案二擬合值與實際值之差越來越大，且方案一中擬合值與實際值之差最大，甚至達到了800周，這都是由于在參加擬合的m個觀測值中存在周跳影響的，而方案三可明顯看出采用抗差多項式擬合法后每個歷元的擬合值與實際值之差都較小且在0.01周以內，這是因為方案三可對參加擬合的數據進行篩選，對存在周跳的數據不用事先利用高次差法判定后，予以剔除，這樣可以確保參加擬合的m個觀測值為“干凈”的，得到的可視化[11]擬合系數也更加準確和可靠。?

人為在非參加擬合的多項式數據中第14秒（即第15個觀測歷元）上加上1周周跳（如圖6）。

三種方案擬合能力的比較，前一段時間3條曲線變化比較平緩，較符合實際情況，并在第15個觀測歷元上發現周跳，探測結果（如表1）;隨著觀測歷元數的增大，方案一與方案二擬合中誤差越來越大，然而方案三的擬合誤差在0.01周，這表明：利用抗差多項式法即可探測出周跳的位置和大小，又可達到修復周跳的目的。

人為在第14秒（即第15個觀測歷元）上加上1周周跳、在第20秒上加上5周周跳、在第30秒上加上10周周跳（如圖7），三種方案擬合能力的比較，前一段時間3條曲線變化較符合實際情況，并在第14、20、30秒觀測時刻發現周跳。隨著觀測歷元數的增大，方案一與方案二擬合中誤差越來越大，雖然有連續周跳卻不影響其擬合精度，然而方案三的擬合誤差在有連續周跳情況下仍然控制在0.01周內，這表明：利用抗差多項式法也可探測出連續周跳的位置和大小，見下（如表2）并對其周跳進行修復[12]。

4 ?結論

本文介紹了周跳產生的原因和周跳對定位精度的影響，在此基礎上，以抗差多項式擬合法周跳探測與修復為研究對象，得出以下結論：采用抗差估計多項式擬合法，由于對傳統多項式擬合進行了重新定權的改進，將含有周跳的擬合數據剔除，所以外推出來的擬合值與實際值之差會很小。在參加擬合的數據不存在周跳的情況下，顯然采用抗差多項式擬合法探測出周跳的能力更好;由于參加擬合的數據不存在周跳，所以擬合值不受加入周跳的影響，而采用抗差多項式擬合法能將周跳修復。

采用抗差估計多項式擬合法對傳統多項式擬合進行了改進，解決了傳統多項式擬合中擬合值與實際觀測值之差發散的問題，結合抗差估計保證了估值的實際抗差性和可靠性[13]。通過實測數據進行實驗，結果驗證了該方法探測與修復周跳的效果比傳統的多項式擬合法和在傳統多項式擬合基礎上只進行歷元間求差的方法更準確可靠。

參考文獻

[1] 周忠謨， 易杰軍， 周琪編著. GPS衛星測量原理與應用[M].湖北武漢. 測繪出版社， 1995： 90-91

[2] 李明， 高星偉， 徐愛功. 一種改進的周跳多項式擬合方法[J]. 測繪科學， 2008， 33（4）： 82-83.

[3] 孫曉， 王勇， 王寶山. MATLAB軟件在測量平差中的應用[J]. 焦作工學院學（自然科學版）， 2002（5）.

[4] 蔡旭輝， 劉衛國， 蔡立艷. MATLAB基礎與應用教程[M].人民郵電出版社， 2009： 1-2.

[5] 關昊. GPS載波相位周跳探測與修復方法的研究[D]. 遼寧工程技術大學， 2011.

[6] 王福麗， 成英燕， 韋鋮， 王曉明. 利用抗差多項式擬合法探測修GNSS周跳[J]. 大地測量與地球動力學， 2013（3）： 129-132.

[7] 邢會敏， 楊福芹. 基于改進多項式擬合法的周跳探測與修復[J]. 科技信息， 2011， （22）： 108-109.

[8] 胡夢英， 賀祖國. 基于重開始共軛思想的改進多項式插值法[J]. 軟件， 2015， 36（11）： 48-51

[9] BISNATH S B. Efficient automated cycle-slip correction of dual-frequency kinem atic GPS data[A]. Proceedings of ION GPS 2000. the l3th International Technical Meeting of The Institute of Navigation[C]， Salt Lake City， Utah， 2000， 145- 154.

[10] 聶敬云， 李春青， 李威威， 等. 關于遺傳算法優化的最小二乘支持向量機在MBR仿真預測中的研究[J]. 軟件， 2015， 36（5）： 40-44.

[11] 孫彥超， 章堅民. 基于MATLAB的DG選址定容可視化系統的設計與實現[J]. 軟件， 2018， 39（4）： 142-147.

[12] 胡云康， 姜蘇， 吳志榮， 等. 基于改進的紋理合成圖像修復算法[J]. 軟件， 2016， 37（4）： 60-63.

[13] 史小玲. 集團網網絡可靠性的探討及實例研究[J]. 軟件， 2016， 37（01）： 117-119.