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

于佳慧 張世濤 張巖

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

關(guān)鍵詞： MATLAB;GNSS;周跳的探測與修復(fù);抗差多項式

中圖分類號： TP3 ? ?文獻(xiàn)標(biāo)識碼： A ? ?DOI：10.3969/j.issn.1003-6970.2019.08.041

本文著錄格式：于佳慧，張世濤，張巖. 基于MATLAB的抗差多項式擬合方法在GNSS周跳探測中的應(yīng)用研究[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具有使用觀測精度高、測量時間短、全天候、技術(shù)自動化程度高、觀測范圍大、能提供三維坐標(biāo)測量成果等特點[1]。在GPS數(shù)據(jù)處理過程中，周跳的存在會使觀測值中出現(xiàn)一個偏差，這會使觀測值失真，若不修復(fù)周跳，剔除這些錯誤成果，會對測繪工作造成嚴(yán)重后果。因此，要想提高定位精度就必須能準(zhǔn)確地探測周跳并且修復(fù)載波相位。目前可用于接收機(jī)的周跳探測方法主要有高次差法、電離層求差法、偽距和載波相位觀測值組合、多項式擬合法來探測和修復(fù)周跳等。

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

1 ?MATLAB簡介

1.1 ?什么是MATLAB

MATLAB軟件是由美國Math Works公司推出的一種交互式、面向?qū)ο蟮某绦蛟O(shè)計語言，廣泛用于工業(yè)界和學(xué)術(shù)界。不僅能很好的處理數(shù)值問題，同時還能解決符號問題，繪制圖形。

1.2 ?MATLAB所具有的優(yōu)勢

MATLAB系統(tǒng)包括以下五個部分：

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

（2） MATLAB具有完整工作界面。包括MATLAB編程和調(diào)試的環(huán)境，對工作空間中的變量進(jìn)行管理及進(jìn)行輸入輸出的數(shù)據(jù)[4]。

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

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

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

2 ?抗差多項式擬合法

2.1 ?周跳概念和產(chǎn)生的原因

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

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

2.2 ?抗差估計求解擬合系數(shù)

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

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

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

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

本部分首先對高次差法的原理和適用范圍優(yōu)缺點進(jìn)行了介紹，其次介紹了利用傳統(tǒng)多項式探測周跳原理，在此基礎(chǔ)上進(jìn)行時間標(biāo)準(zhǔn)化、歷元間求差、重新定權(quán)等改進(jìn)，從而建立抗差多項式模型，利用抗差多項式探測和修復(fù)周跳。

3 ?編制周跳探測與修復(fù)系統(tǒng)程序

3.1 ?程序的設(shè)計框架

研究利用多項式[8]擬合法周跳探測和修復(fù)方法，對傳統(tǒng)的多項式擬合法周跳探測和修復(fù)方法進(jìn)行改進(jìn)，并確定多項式擬合的階數(shù)。

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

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

1. 解決多項式的基礎(chǔ)問題

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

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

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

2. 高次差法檢驗觀測數(shù)據(jù)質(zhì)量

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

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

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

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

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

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

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

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

4 ?結(jié)論

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

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

參考文獻(xiàn)

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

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

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

[4] 蔡旭輝， 劉衛(wèi)國， 蔡立艷. MATLAB基礎(chǔ)與應(yīng)用教程[M].人民郵電出版社， 2009： 1-2.

[5] 關(guān)昊. GPS載波相位周跳探測與修復(fù)方法的研究[D]. 遼寧工程技術(shù)大學(xué)， 2011.

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

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

[8] 胡夢英， 賀祖國. 基于重開始共軛思想的改進(jìn)多項式插值法[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] 聶敬云， 李春青， 李威威， 等. 關(guān)于遺傳算法優(yōu)化的最小二乘支持向量機(jī)在MBR仿真預(yù)測中的研究[J]. 軟件， 2015， 36（5）： 40-44.

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

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

[13] 史小玲. 集團(tuán)網(wǎng)網(wǎng)絡(luò)可靠性的探討及實例研究[J]. 軟件， 2016， 37（01）： 117-119.