于佳慧 張世濤 張巖
摘 ?要: 為了提高GNSS載波相位的定位精度,本文采用歷元間求差與抗差估計(jì)相結(jié)合的方法對(duì)過(guò)去常常使用的多項(xiàng)式進(jìn)行改良,最終提出用抗差多項(xiàng)式擬合來(lái)探測(cè)并修復(fù)周跳。本文基于Matlab軟件編制程序,并對(duì)GNSS載波相位周跳進(jìn)行了試算。實(shí)驗(yàn)結(jié)果表明:此方法可以做到對(duì)GNSS周跳探測(cè)與修復(fù)功能,具有文件讀寫、衛(wèi)星觀測(cè)值提取、周跳探測(cè)與修復(fù)、繪制擬合的誤差結(jié)果圖形等功能。通過(guò)這種方法可以探測(cè)并修復(fù)1周以上的周跳。得出結(jié)論:通過(guò)抗差多項(xiàng)式擬合方法可以探測(cè)并修復(fù)GNSS周跳。
關(guān)鍵詞: MATLAB;GNSS;周跳的探測(cè)與修復(fù);抗差多項(xiàng)式
中圖分類號(hào): TP3 ? ?文獻(xiàn)標(biāo)識(shí)碼: A ? ?DOI:10.3969/j.issn.1003-6970.2019.08.041
本文著錄格式:于佳慧,張世濤,張巖. 基于MATLAB的抗差多項(xiàng)式擬合方法在GNSS周跳探測(cè)中的應(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具有使用觀測(cè)精度高、測(cè)量時(shí)間短、全天候、技術(shù)自動(dòng)化程度高、觀測(cè)范圍大、能提供三維坐標(biāo)測(cè)量成果等特點(diǎn)[1]。在GPS數(shù)據(jù)處理過(guò)程中,周跳的存在會(huì)使觀測(cè)值中出現(xiàn)一個(gè)偏差,這會(huì)使觀測(cè)值失真,若不修復(fù)周跳,剔除這些錯(cuò)誤成果,會(huì)對(duì)測(cè)繪工作造成嚴(yán)重后果。因此,要想提高定位精度就必須能準(zhǔn)確地探測(cè)周跳并且修復(fù)載波相位。目前可用于接收機(jī)的周跳探測(cè)方法主要有高次差法、電離層求差法、偽距和載波相位觀測(cè)值組合、多項(xiàng)式擬合法來(lái)探測(cè)和修復(fù)周跳等。
總之,有很多方法可以探測(cè)和修復(fù)周跳,但不同的周跳探測(cè)與修復(fù)方法對(duì)提高計(jì)算效率及可靠性都有所不同。常規(guī)方法對(duì)載波相位觀測(cè)值周跳探測(cè)和修復(fù)雖然都有其優(yōu)點(diǎn),但也存在其局限性,鑒于此,本文將以多項(xiàng)式擬合法為基礎(chǔ),先進(jìn)行時(shí)間標(biāo)準(zhǔn)化,然后采用歷元間求差與抗差估計(jì)相結(jié)合的方法對(duì)傳統(tǒng)算法中存在的問題進(jìn)行改進(jìn),并利用GNSS/GLONASS實(shí)測(cè)數(shù)據(jù)驗(yàn)證改進(jìn)后的算法的有效性[2]。
1 ?MATLAB簡(jiǎn)介
1.1 ?什么是MATLAB
MATLAB軟件是由美國(guó)Math Works公司推出的一種交互式、面向?qū)ο蟮某绦蛟O(shè)計(jì)語(yǔ)言,廣泛用于工業(yè)界和學(xué)術(shù)界。不僅能很好的處理數(shù)值問題,同時(shí)還能解決符號(hào)問題,繪制圖形。
1.2 ?MATLAB所具有的優(yōu)勢(shì)
MATLAB系統(tǒng)包括以下五個(gè)部分:
(1)MATLAB具有多種語(yǔ)言體系。在MATLAB可用的編程語(yǔ)言,對(duì)應(yīng)于工具箱中六個(gè)目錄,分別為ops(操作符與特殊字符)、lang(編程語(yǔ)言結(jié)構(gòu))、strfun(字符串操作)、iofun(文件的輸入輸出)、time fun(時(shí)間和日期)和data type(數(shù)據(jù)類型和結(jié)構(gòu))[3]。
(2) MATLAB具有完整工作界面。包括MATLAB編程和調(diào)試的環(huán)境,對(duì)工作空間中的變量進(jìn)行管理及進(jìn)行輸入輸出的數(shù)據(jù)[4]。
(3)MATLAB具有強(qiáng)大的圖形處理能力。其中包括二、三維圖形可視化,也包括定義圖形外觀的操作,最終可以為MATLAB應(yīng)用程序創(chuàng)建整套的用戶界面。
(4)數(shù)學(xué)函數(shù)庫(kù)。包括初等函數(shù);也包括更復(fù)雜的功能,如矩陣功能(矩陣求逆、求矩陣的特征值)、快速傅里葉變換等[3]。
(5)MATLAB應(yīng)用程序接口(API)。提供接口程序使程序員可以編寫C/C6、Fortran程序調(diào)用MATLAB作為計(jì)算引擎,并用MATLAB調(diào)用外部應(yīng)用程序或者讀寫MATLAB文件[4]。MATLAB中的M文件的語(yǔ)法與其他的高級(jí)語(yǔ)言類似,但語(yǔ)法較其他一般的高級(jí)語(yǔ)言來(lái)說(shuō)要簡(jiǎn)單,它只是一個(gè)ASCII碼簡(jiǎn)單的文本文件,其語(yǔ)法編制與數(shù)學(xué)語(yǔ)言非常接近。它不僅是一種程序化的編程語(yǔ)言,也是一種解釋性的編程語(yǔ)言,可逐行解釋和運(yùn)行程序,程序更易調(diào)試[4]。
2 ?抗差多項(xiàng)式擬合法
2.1 ?周跳概念和產(chǎn)生的原因
載波相位可以用在高精度衛(wèi)星導(dǎo)航定位。完整且連續(xù)的周期通過(guò)使用多普勒計(jì)數(shù)完成,而載波相位在導(dǎo)航定位中主要是測(cè)量不完整的周期。信噪比低、信號(hào)質(zhì)量差以及接收機(jī)接收到的質(zhì)量差等方面都可能引起完整的周期部分突變,叫作整周跳變,簡(jiǎn)稱為周跳(cycle slip)。
載波相位是由 、 、 這些部分組成的。雖然可以以極高的精度測(cè)定 ,但也要正確地測(cè)定 和 的情況下,才能確定載波相位值。在載波相位測(cè)量的實(shí)際觀測(cè)值中,小于一周的部分稱為觀測(cè)時(shí)刻的觀測(cè)值。衛(wèi)星載波信號(hào)與接收機(jī)參考振蕩信號(hào)形成的差頻信號(hào)不足一周。接收機(jī)載波跟蹤回路中的鑒相器測(cè)量可以實(shí)現(xiàn)對(duì)載波信號(hào)的測(cè)量。只有接收機(jī)的載波信號(hào)和參考振蕩信號(hào)能在觀測(cè)時(shí)段產(chǎn)生差頻信號(hào),才可以得到正確的觀測(cè)值。整個(gè)星期的計(jì)數(shù)不一樣。它是計(jì)數(shù)器從第一個(gè)觀測(cè)時(shí)間到當(dāng)前觀測(cè)時(shí)間積累的差分頻率信號(hào)中的整數(shù)頻帶數(shù)。如果由于外界因素,通過(guò)使用多普勒計(jì)數(shù)法的計(jì)數(shù)器在兩個(gè)觀測(cè)周期之間的某個(gè)時(shí)段停止計(jì)數(shù),從而使整個(gè)周期的計(jì)數(shù)比實(shí)際值少n周,那么當(dāng)計(jì)數(shù)器恢復(fù)正常運(yùn)行時(shí),所有載波相位波段將包含相同的誤差值,該誤差值小于整個(gè)周期的正常值(如圖1)。這樣一個(gè)系統(tǒng)性的誤差,在一個(gè)周期內(nèi)計(jì)數(shù)和不足一周的部分保持正確計(jì)數(shù)的情況,稱為整周跳變,簡(jiǎn)稱周跳[5]。
2.2 ?抗差估計(jì)求解擬合系數(shù)
傳統(tǒng)的多項(xiàng)式擬合方法一般計(jì)算擬合系數(shù)是采用最小二乘法,但當(dāng)觀測(cè)值不服從正態(tài)分布的假設(shè)時(shí),即m個(gè)擬合觀測(cè)值出現(xiàn)粗差或周跳時(shí),最小二乘法估計(jì)就不具有抗干擾性,錯(cuò)誤的估計(jì)結(jié)果可能由單個(gè)觀測(cè)值的偏差導(dǎo)致。而利用抗差估計(jì)法與利用最小二乘估計(jì)法不同,它不像最小二乘那樣過(guò)分的追求有效性與無(wú)偏性等內(nèi)部性質(zhì),抗差估計(jì)則更加重視估值的實(shí)際可靠性和抗差性[6]。從第一個(gè)歷元開始,先一次取m個(gè)歷元 和對(duì)應(yīng)的載波相位觀測(cè)值 代入式(1),通過(guò)抗差估計(jì)法解算其擬合系數(shù)。假設(shè)初始權(quán)矩陣為單位陣,即 ,之后進(jìn)行平差計(jì)算。根據(jù)第一次平差后得到的改正數(shù)陣 和單位權(quán)中誤差 ,利用下式(2)
對(duì)各相位觀測(cè)值重新定權(quán),如果 則認(rèn)為m個(gè)歷元觀測(cè)值中不存在周跳,最后的平差結(jié)果與第一次平差結(jié)果相同。若有偏差,可利用權(quán)陣 在再次進(jìn)行平差計(jì)算,直至 時(shí)停止迭代,取第j次迭代結(jié)果作為最終結(jié)果。這樣可以確保參加擬合的m個(gè)相位觀測(cè)值為不含周跳的值,從而確保得到的擬合系數(shù)更加準(zhǔn)確和可靠。
2.3 ?利用抗差多項(xiàng)式探測(cè)和修復(fù)周跳
根據(jù)求得的擬合系數(shù)擬合,并帶入m+1個(gè)歷元值,來(lái)估計(jì)第m+1個(gè)歷元的相位值,將利用抗差多項(xiàng)式擬合出來(lái)的擬合值與對(duì)應(yīng)歷元的實(shí)際觀測(cè)作差,若差值 時(shí),則認(rèn)為 不含周跳,去掉本次擬合的第一個(gè)觀測(cè)值,加入第m+1個(gè)歷元觀測(cè)值繼續(xù)進(jìn)行擬合;當(dāng)差值 時(shí),則認(rèn)為 含有周跳,采用外推的整周計(jì)數(shù)去取代有周跳的實(shí)際觀測(cè)值中的整周計(jì)數(shù),但不足一周的部分仍然保持不動(dòng),然后繼續(xù)上述過(guò)程,直到最后一個(gè)相位觀測(cè)值為止[7]。
本部分首先對(duì)高次差法的原理和適用范圍優(yōu)缺點(diǎn)進(jìn)行了介紹,其次介紹了利用傳統(tǒng)多項(xiàng)式探測(cè)周跳原理,在此基礎(chǔ)上進(jìn)行時(shí)間標(biāo)準(zhǔn)化、歷元間求差、重新定權(quán)等改進(jìn),從而建立抗差多項(xiàng)式模型,利用抗差多項(xiàng)式探測(cè)和修復(fù)周跳。
3 ?編制周跳探測(cè)與修復(fù)系統(tǒng)程序
3.1 ?程序的設(shè)計(jì)框架
研究利用多項(xiàng)式[8]擬合法周跳探測(cè)和修復(fù)方法,對(duì)傳統(tǒng)的多項(xiàng)式擬合法周跳探測(cè)和修復(fù)方法進(jìn)行改進(jìn),并確定多項(xiàng)式擬合的階數(shù)。
在無(wú)周跳“干凈”觀測(cè)數(shù)據(jù)中加入周跳,通過(guò)差分技術(shù)對(duì)相鄰歷元之間的觀測(cè)數(shù)據(jù)進(jìn)行差分,對(duì)差分?jǐn)?shù)據(jù)進(jìn)行多項(xiàng)式建模擬合與預(yù)報(bào),探測(cè)周跳;利用抗差多項(xiàng)式建模探測(cè)周跳,即編制抗差多項(xiàng)式程序設(shè)計(jì)框架(如圖2)。
3.2 ?周跳探測(cè)程序算例分析
1. 解決多項(xiàng)式的基礎(chǔ)問題
衛(wèi)星到地面的距離對(duì)時(shí)間的四階以上導(dǎo)數(shù)已趨于零,其變化規(guī)律是隨機(jī)的,無(wú)法再用多項(xiàng)式進(jìn)行擬合。故多項(xiàng)式擬合的階數(shù)n取3~4階,在本實(shí)驗(yàn)中n=3的條件下既能保證精度又保證了運(yùn)算速度。
擬合窗口m選取時(shí),當(dāng)擬合窗口寬度越大時(shí),外推值越準(zhǔn)確,擬合中誤差越小,但過(guò)大又會(huì)增加計(jì)算量。當(dāng)擬合窗口寬度越小時(shí),外推值精度越差。所以通過(guò)取不同的值進(jìn)行實(shí)驗(yàn),根據(jù)實(shí)驗(yàn)取m=10比較合適。
門檻系數(shù)k的取值是隨不同的情況改變的。當(dāng)k較大時(shí),表示只有當(dāng)觀測(cè)值偏離外推值很大的時(shí)候才認(rèn)為它是異常的;當(dāng)k較小時(shí),表示當(dāng)觀測(cè)值偏離外推值較小的時(shí)候就認(rèn)為它是異常值了。一般取k在2到9之間的數(shù)。經(jīng)過(guò)多次取值實(shí)驗(yàn),結(jié)果表明門檻系數(shù)取k=4時(shí)效果最好。
2. 高次差法檢驗(yàn)觀測(cè)數(shù)據(jù)質(zhì)量
算例基于2014年1月14日國(guó)內(nèi)的GPS/GLONASS的實(shí)測(cè)數(shù)據(jù),觀測(cè)時(shí)間為10分鐘,采樣間隔為1秒。使用高階差分的方法先對(duì)16時(shí)45分0秒到16時(shí)55分0秒GPS的10號(hào)衛(wèi)星共600個(gè)L1載波的據(jù)進(jìn)行周跳探測(cè),對(duì)觀測(cè)值求四階差分的時(shí)間序列圖(如圖3),從(如圖3)可以看出,實(shí)驗(yàn)結(jié)果中沒有周跳[9]。
3. 多項(xiàng)式擬合法算例分析
本文選取上述無(wú)周跳的數(shù)據(jù)中從16時(shí)45分0秒到16時(shí)45分60秒的實(shí)測(cè)GPS數(shù)據(jù)進(jìn)行多項(xiàng)式擬合實(shí)驗(yàn),以第3秒(即第4個(gè)歷元)到第12秒(即第13個(gè)歷元)的數(shù)據(jù)作為擬合多項(xiàng)式的數(shù)據(jù),外推出從第13秒到第60秒的數(shù)據(jù),并設(shè)計(jì)實(shí)驗(yàn)方案:方案1,傳統(tǒng)的多項(xiàng)式擬合法;方案2,歷元間求差并采用最小二乘[10]平差的多項(xiàng)式擬合法;方案3,歷元間求差并采用抗差估計(jì)的多項(xiàng)式擬合法。
在無(wú)周跳情況下(如圖4),三種方案擬合能力的比較,前一段時(shí)間曲線變化比較平緩,較符合實(shí)際情況;隨著觀測(cè)歷元數(shù)的增大,方案一與方案二曲線的波動(dòng)越來(lái)越大,出現(xiàn)發(fā)散的情況,有可能造成周跳判斷失誤。方案三可明顯看出采用抗差估計(jì)后每個(gè)歷元的擬合中誤差都有較大改善,GPS觀測(cè)量擬合值與實(shí)際值之差在0.01周以內(nèi)。
人為在第6秒(即第7個(gè)觀測(cè)歷元)上加上1周周跳(如圖5)三種方案擬合效果的比較,前一段時(shí)間3條曲線變化比較平緩;隨著觀測(cè)歷元數(shù)的增大,方案一與方案二擬合值與實(shí)際值之差越來(lái)越大,且方案一中擬合值與實(shí)際值之差最大,甚至達(dá)到了800周,這都是由于在參加擬合的m個(gè)觀測(cè)值中存在周跳影響的,而方案三可明顯看出采用抗差多項(xiàng)式擬合法后每個(gè)歷元的擬合值與實(shí)際值之差都較小且在0.01周以內(nèi),這是因?yàn)榉桨溉蓪?duì)參加擬合的數(shù)據(jù)進(jìn)行篩選,對(duì)存在周跳的數(shù)據(jù)不用事先利用高次差法判定后,予以剔除,這樣可以確保參加擬合的m個(gè)觀測(cè)值為“干凈”的,得到的可視化[11]擬合系數(shù)也更加準(zhǔn)確和可靠。?
人為在非參加擬合的多項(xiàng)式數(shù)據(jù)中第14秒(即第15個(gè)觀測(cè)歷元)上加上1周周跳(如圖6)。
三種方案擬合能力的比較,前一段時(shí)間3條曲線變化比較平緩,較符合實(shí)際情況,并在第15個(gè)觀測(cè)歷元上發(fā)現(xiàn)周跳,探測(cè)結(jié)果(如表1);隨著觀測(cè)歷元數(shù)的增大,方案一與方案二擬合中誤差越來(lái)越大,然而方案三的擬合誤差在0.01周,這表明:利用抗差多項(xiàng)式法即可探測(cè)出周跳的位置和大小,又可達(dá)到修復(fù)周跳的目的。
人為在第14秒(即第15個(gè)觀測(cè)歷元)上加上1周周跳、在第20秒上加上5周周跳、在第30秒上加上10周周跳(如圖7),三種方案擬合能力的比較,前一段時(shí)間3條曲線變化較符合實(shí)際情況,并在第14、20、30秒觀測(cè)時(shí)刻發(fā)現(xiàn)周跳。隨著觀測(cè)歷元數(shù)的增大,方案一與方案二擬合中誤差越來(lái)越大,雖然有連續(xù)周跳卻不影響其擬合精度,然而方案三的擬合誤差在有連續(xù)周跳情況下仍然控制在0.01周內(nèi),這表明:利用抗差多項(xiàng)式法也可探測(cè)出連續(xù)周跳的位置和大小,見下(如表2)并對(duì)其周跳進(jìn)行修復(fù)[12]。
4 ?結(jié)論
本文介紹了周跳產(chǎn)生的原因和周跳對(duì)定位精度的影響,在此基礎(chǔ)上,以抗差多項(xiàng)式擬合法周跳探測(cè)與修復(fù)為研究對(duì)象,得出以下結(jié)論:采用抗差估計(jì)多項(xiàng)式擬合法,由于對(duì)傳統(tǒng)多項(xiàng)式擬合進(jìn)行了重新定權(quán)的改進(jìn),將含有周跳的擬合數(shù)據(jù)剔除,所以外推出來(lái)的擬合值與實(shí)際值之差會(huì)很小。在參加擬合的數(shù)據(jù)不存在周跳的情況下,顯然采用抗差多項(xiàng)式擬合法探測(cè)出周跳的能力更好;由于參加擬合的數(shù)據(jù)不存在周跳,所以擬合值不受加入周跳的影響,而采用抗差多項(xiàng)式擬合法能將周跳修復(fù)。
采用抗差估計(jì)多項(xiàng)式擬合法對(duì)傳統(tǒng)多項(xiàng)式擬合進(jìn)行了改進(jìn),解決了傳統(tǒng)多項(xiàng)式擬合中擬合值與實(shí)際觀測(cè)值之差發(fā)散的問題,結(jié)合抗差估計(jì)保證了估值的實(shí)際抗差性和可靠性[13]。通過(guò)實(shí)測(cè)數(shù)據(jù)進(jìn)行實(shí)驗(yàn),結(jié)果驗(yàn)證了該方法探測(cè)與修復(fù)周跳的效果比傳統(tǒng)的多項(xiàng)式擬合法和在傳統(tǒng)多項(xiàng)式擬合基礎(chǔ)上只進(jìn)行歷元間求差的方法更準(zhǔn)確可靠。
參考文獻(xiàn)
[1] 周忠謨, 易杰軍, 周琪編著. GPS衛(wèi)星測(cè)量原理與應(yīng)用[M].湖北武漢. 測(cè)繪出版社, 1995: 90-91
[2] 李明, 高星偉, 徐愛功. 一種改進(jìn)的周跳多項(xiàng)式擬合方法[J]. 測(cè)繪科學(xué), 2008, 33(4): 82-83.
[3] 孫曉, 王勇, 王寶山. MATLAB軟件在測(cè)量平差中的應(yīng)用[J]. 焦作工學(xué)院學(xué)(自然科學(xué)版), 2002(5).
[4] 蔡旭輝, 劉衛(wèi)國(guó), 蔡立艷. MATLAB基礎(chǔ)與應(yīng)用教程[M].人民郵電出版社, 2009: 1-2.
[5] 關(guān)昊. GPS載波相位周跳探測(cè)與修復(fù)方法的研究[D]. 遼寧工程技術(shù)大學(xué), 2011.
[6] 王福麗, 成英燕, 韋鋮, 王曉明. 利用抗差多項(xiàng)式擬合法探測(cè)修GNSS周跳[J]. 大地測(cè)量與地球動(dòng)力學(xué), 2013(3): 129-132.
[7] 邢會(huì)敏, 楊福芹. 基于改進(jìn)多項(xiàng)式擬合法的周跳探測(cè)與修復(fù)[J]. 科技信息, 2011, (22): 108-109.
[8] 胡夢(mèng)英, 賀祖國(guó). 基于重開始共軛思想的改進(jìn)多項(xiàng)式插值法[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ù)測(cè)中的研究[J]. 軟件, 2015, 36(5): 40-44.
[11] 孫彥超, 章堅(jiān)民. 基于MATLAB的DG選址定容可視化系統(tǒng)的設(shè)計(jì)與實(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ò)可靠性的探討及實(shí)例研究[J]. 軟件, 2016, 37(01): 117-119.