湯 云,金元剛,付 平

(1.解放軍92941部隊,遼寧葫蘆島125001;2.中海油能源發(fā)展有限公司,天津300450)

0 引 言

在動態(tài)測量數(shù)據(jù)處理中,系統(tǒng)誤差是很棘手的問題,研究系統(tǒng)誤差的有效處理方法正受到越來越多的關(guān)注。系統(tǒng)誤差的特性非常復(fù)雜,既具有延續(xù)性,又有偶然性,因而很難用統(tǒng)一的、確定性的模型來準(zhǔn)確描述它。就GPS動態(tài)測量來說,信號受到多種誤差源的誤差干擾,有的誤差(例如電離層折射影響)可以利用雙頻接收機的方式得到較好補償,而多路徑影響卻很難找到有效的改正模型。系統(tǒng)誤差雖經(jīng)部分模型改正,或差分改正,仍有部分殘余留在觀測量中,或多或少影響到定位結(jié)果的精度[1]。

在很多GPS定位技術(shù)的文獻(xiàn)資料上,定位誤差都按如下方法進(jìn)行估計[2]

式中:ei是第i次采樣的誤差;n為樣本容量;ˉx為系統(tǒng)誤差;σ為標(biāo)準(zhǔn)差。

上述誤差估計方法的誤差模型是用隨機變量來描述的。對于GPS靜態(tài)定位或靜態(tài)測量時是正確的,但將這種方法用于GPS動態(tài)定位或動態(tài)測量中一般是不正確的,因為只有當(dāng)系統(tǒng)誤差是常量時,該方法才可適用。而實際上GPS動態(tài)測量的系統(tǒng)誤差一般并不是常量,其誤差曲線并不水平。對于GPS動態(tài)測量系統(tǒng)的誤差估計,應(yīng)力求提高誤差估計結(jié)果的精度或可信度。為了科學(xué)客觀地評定一個GPS動態(tài)測量系統(tǒng)的性能,有必要研究適合GPS動態(tài)測量的數(shù)據(jù)處理方法。

1 GPS動態(tài)測量的誤差特性

GPS動態(tài)測量的誤差是隨機過程或經(jīng)采樣以后形成的隨機序列,當(dāng)用戶臺處于動態(tài)時,這是很顯然的。在GPS動態(tài)測量系統(tǒng)的鑒定試驗中,可以采用靜態(tài)目標(biāo)試驗動態(tài)性能,這種方法盡管是采用靜止目標(biāo),但本質(zhì)上不同于靜態(tài)測量。因為噪聲引起的誤差、傳播介質(zhì)產(chǎn)生的誤差(包括電離層、對流層折射和多徑效應(yīng)引起的誤差)、衛(wèi)星鐘差、衛(wèi)星星歷誤差等經(jīng)差分后的殘余誤差都是隨機過程,采用單歷元解算,這種特性被保留下來并離散化為隨機時間序列。

GPS測量的誤差中,傳播介質(zhì)產(chǎn)生的誤差、衛(wèi)星鐘差、衛(wèi)星星歷誤差等存在隨時間慢變的成分;不同歷元可能有不同導(dǎo)航衛(wèi)星的信息參與運算,其DOP值一般情況下也是不同的,由此產(chǎn)生的誤差也是隨時間變化的;在載波相位差分定位中,周跳引起的誤差在不同歷元也是變化的。這些隨時間慢變的誤差成分形成了系統(tǒng)誤差,它們并不表現(xiàn)為常值,而是時間的隨機函數(shù)。除上述慢變成分外,它們的快變部分以及噪聲引起的誤差則成為噪聲狀態(tài)的隨機誤差部分。

GPS動態(tài)測量系統(tǒng)的誤差是由熱噪聲、量化噪聲、電離層和對流層折射產(chǎn)生的誤差、衛(wèi)星鐘差、衛(wèi)星星歷誤差、基準(zhǔn)站誤差、模糊度解算誤差等許多誤差分量綜合的結(jié)果。在接收機和衛(wèi)星系統(tǒng)正常工作的情況下,各種誤差對總誤差的影響都不顯著。根據(jù)中心極限定理,系統(tǒng)測量誤差的概率分布為正態(tài)分布。

綜合以上分析,GPS動態(tài)測量系統(tǒng)的動態(tài)測量誤差是非平穩(wěn)的正態(tài)隨機序列。誤差的這一特性是研究數(shù)據(jù)處理方法的出發(fā)點。

2 處理方法的理論基礎(chǔ)

2.1 非平穩(wěn)隨機過程誤差處理方法

動態(tài)測量誤差是非平穩(wěn)的正態(tài)隨機過程或經(jīng)采樣得到的隨機序列,用x(t)表示。其中有變化緩慢的誤差成分y(t)和快變的誤差成分η(t),即

式中:慢變誤差y(t)稱為系統(tǒng)誤差,而把零均值的η(t)稱作隨機誤差。

根據(jù)隨機過程的理論,估計動態(tài)測量誤差 x(t)的均值和方差需得到 x(t)在ti時刻的抽樣值x1(ti),x2(ti),…,xn(ti),它們是來自x(ti)的樣本點。其中xj(t)是指來自x(t)的第j個樣本記錄,對應(yīng)于系統(tǒng)在相同狀態(tài)下(以保證它們來自同一過程x(t))所做的第j次試驗(第j航次、架次)。用隨機變量樣本值的統(tǒng)計方法(即式(1)和式(2))得到ti時刻的估計。對記錄的每一時刻都做如此處理,就可得到均值函數(shù)和方差函數(shù)[3]。

但是上述處理方法對GPS動態(tài)測量誤差估計存在很大困難。需要做許多次試驗才能保證誤差估計的精度。這多次試驗要保持系統(tǒng)狀態(tài)相同可以說是無法實現(xiàn)的。比如,不同次數(shù)的試驗要求跟蹤相同的衛(wèi)星,要求DOP值相等,這些都是不可能的。因此動態(tài)誤差的估計需要另想辦法。

2.2 各態(tài)歷經(jīng)隨機過程誤差處理方法

按照隨機過程理論,如果過程具有各態(tài)歷經(jīng)性,則上述估計等價于對一個記錄的相應(yīng)時間均值的估計,即可通過對一個記錄時間均值的計算得到x(t)的均值和方差的估計[4]。

經(jīng)過論證有如下結(jié)論:對于只求正態(tài)隨機過程的一、二階矩,如果過程是平穩(wěn)的,則具有各態(tài)歷經(jīng)。但是如前所述x(t)一般并不平穩(wěn)。

對于式(3)中的系統(tǒng)誤差y(t),并不同于雜亂無章的噪聲,它的出現(xiàn)具有一定規(guī)律性。只要試驗時間足夠長,即可用一個現(xiàn)實(一次試驗)數(shù)據(jù)得到。能否用時間均值作為動態(tài)誤差的估計,關(guān)鍵是看噪聲η(t)是否具有各態(tài)歷經(jīng)。將x(t)中的系統(tǒng)誤差y(t)分離以后,所剩η(t)如果平穩(wěn),理論上就可用一個現(xiàn)實的數(shù)據(jù)完成誤差統(tǒng)計。

3 系統(tǒng)誤差分離

系統(tǒng)誤差是具有規(guī)律性的慢變誤差,可以用一個解析函數(shù)逼近。由函數(shù)逼近理論,如果函數(shù) f(t)在區(qū)間[a,b]上連續(xù),則在該區(qū)間上總可用一個代數(shù)多項式P(t)去一致逼近。用n階多項式a0P0(t)+a1P1(t)+…+anPn(t)作為系統(tǒng)誤差y(t)的模型,由此得

其離散序列為

式中:k=1,2,…,N。N為序列長度。令

參數(shù)A的最優(yōu)估計為

式中:R是Q的方差陣,或者是η(t)的相關(guān)函數(shù)陣。

假如η(t)二階平穩(wěn)且不相關(guān),則可令R-1=I,而實際上η(t)是相關(guān)的,并且有時是強相關(guān),尤其在定位處理中采用了卡爾曼濾波后更是如此。為了避免計算R的困難,常用的方法是獨立抽樣。獨立抽樣的問題存在相關(guān)時間計算是否較為準(zhǔn)確,而主要問題是獨立抽樣引起y(t)的失真。通過頻域研究發(fā)現(xiàn)獨立抽樣引起y(t)的失真幾乎是不可避免的,這種方法不能得到系統(tǒng)誤差的最優(yōu)估計。

求A的最優(yōu)估計可考慮如下算法

①令R=I,得到初始估計A0;

②用A0分離系統(tǒng)誤差,其殘差作為η(k),求方差陣R0。由于A0不是最優(yōu)估計,R0只是R的初始估計;

③令R=R0,用式(5)估計A得到A1,A1是比A0更好的估計;

④用A1分離系統(tǒng)誤差后再次估計R得到R1,R1是比R0更好的估計;

⑤再令R=R1重新估計A,依此循環(huán)直到滿足計算精度。

這種方法實際上會產(chǎn)生很大的計算誤差。GPS動態(tài)測量系統(tǒng)采樣率一般不低于20 Hz,為了充分考核動態(tài)測量精度,一次試驗要連續(xù)記錄30 min以上,誤差序列將有幾萬組以上的采樣值。x(t)的方差陣將達(dá)幾萬階以上,很難保證其正定性。即使能保證正定,這樣高階矩陣的求逆要經(jīng)過許多步的計算,使計算誤差大量累積。

通過研究和實際計算形成了求解系統(tǒng)誤差參數(shù)的如下方法。

零均值平穩(wěn)隨機序列可用ARMA模型描述,而ARMA可用足夠高階的AR模型模擬,因此,零均值平穩(wěn)隨機序列可用AR模型建模

式中:ε(k)為白噪聲;φ(B)為算子多項式

將式(6)代入誤差模型式(4),有

即

令

式(7)為

如果得到N+q個樣本值,有

式中:

由于ε(k)是零均值白噪聲,因此ai的最優(yōu)估計為

可以進(jìn)一步推導(dǎo)出AR模型參數(shù)的估計

式中:Φ=[φ1 φ2 … φL]T,R(L)是L維隨機變量η(k)的方差陣,而 Σ(L)是L維隨機變量η(k)的協(xié)方差組成的列向量。

Φ還可用最小二乘估計得到

式中:

式(11)、式(12)中的R-1(L)和(BTB)-1都只是L階矩陣求逆,而L不會超過10,顯然比幾萬階矩陣求逆誤差小得多。

計算方法與前述逐次求R矩陣的方法相同。首先令φ(B)=1,得到初始估計A0,進(jìn)而得到初始估計 Φ0。利用 Φ0重新計算得到A1和 Φ1,直至達(dá)到計算精度。

4 實測GPS數(shù)據(jù)誤差分析

為了驗證算法的有效性,采用2009年3月20日某無人機飛行記錄的GPS測量數(shù)據(jù)進(jìn)行多項式擬合及誤差分離。無人機飛行了大約1 h,機載GPS的數(shù)據(jù)更新率是20 Hz,為了作圖清晰直觀,采用提取1 Hz的實測數(shù)據(jù)進(jìn)行處理,數(shù)據(jù)長度為1000 s。圖1～圖6是應(yīng)用前面所述的方法對實測數(shù)據(jù)的高度和速度進(jìn)行多項式擬合并進(jìn)行誤差分離的結(jié)果。

圖1 高度測量值及其擬合曲線

圖3 高度誤差自相關(guān)函數(shù)

圖1、2、3表示的是GPS高度數(shù)據(jù)分析結(jié)果。從圖1可以看出,多項式能較好擬合無人機飛行高度曲線;圖2、3顯示高度誤差具有明顯的周期性和隨機性,其自相關(guān)函數(shù)緩慢衰減且周期振蕩。圖4、5、6則表示對與高度記錄相同時段所對應(yīng)的飛行速度曲線分析結(jié)果。從分離出的速度誤差來看,它與高度誤差一樣,具有周期性和隨機性,其自相關(guān)函數(shù)同樣表現(xiàn)出振蕩衰減性質(zhì)。圖2、圖5中的出現(xiàn)數(shù)據(jù)跳變,是因為在此期間無人機處于俯沖階段,姿態(tài)變化較大,引起短暫衛(wèi)星失鎖,使整周模糊度發(fā)生跳變所致。對該時段的其他的GPS數(shù)據(jù)采用上述處理方法,可以得到類似的結(jié)論。

圖5 速度誤差分離曲線

圖6 速度誤差自相關(guān)函數(shù)

5 結(jié) 論

經(jīng)過理論分析及實例驗證,為GPS動態(tài)測量誤差的統(tǒng)計分析確立了完整的分析處理方法。除上述誤差特性、理論基礎(chǔ)、系統(tǒng)誤差分離以外,還包括數(shù)據(jù)的預(yù)處理(數(shù)據(jù)分析、野點判別及修復(fù)等)、多項式和AR模型階數(shù)的判別、非線性模型、平穩(wěn)性檢驗、隨機誤差估計的性質(zhì)及無偏修正、AR模型建模的遞推算法、非線性建模方法、利用BTB矩陣性質(zhì)構(gòu)造多項式(非正交)參數(shù)估計的遞推算法等。限于篇幅,這些問題將在后續(xù)的研究中再補充討論。

[1] 歐吉坤,王振杰.GPS精密測量中系統(tǒng)誤差的分離方法[J].數(shù)據(jù)采集與處理,2003,12(4):365-368

[2] 李洪濤,許國昌,薛鴻印,等.GPS應(yīng)用程序設(shè)計[M].北京:科學(xué)出版社,1999:256-262

[3] 費業(yè)泰.誤差理論與數(shù)據(jù)處理.第5版[M].北京:機械工業(yè)出版社,2005:162-174.

[4] 梁晉文,陳林才,何 貢,等.誤差理論與數(shù)據(jù)處理.第2版[M].北京:中國計量出版社,2001:136-162.