時間相關(guān)觀測的幾種動態(tài)數(shù)據(jù)處理方法

李博峰，章浙濤

同濟(jì)大學(xué)測繪與地理信息學(xué)院，上海 200092

在現(xiàn)代大地測量動態(tài)數(shù)據(jù)處理中，只要觀測值數(shù)量大于未知參數(shù)個數(shù)，即可建立數(shù)學(xué)模型并采用逐點單歷元的平差方法求得所需的動態(tài)解[1-2]。數(shù)學(xué)模型包括函數(shù)模型和隨機(jī)模型兩部分，前者表征的是觀測值和未知參數(shù)之間的關(guān)系，而后者描述的是觀測值的精度以及相互之間的關(guān)系，兩者同等重要且缺一不可[3-4]。

通常，可以認(rèn)為觀測誤差為高斯白噪聲，即假設(shè)觀測值之間不存在時間相關(guān)性。然而在實際觀測數(shù)據(jù)中，由于非模型化誤差的存在[5-6]，導(dǎo)致觀測值中往往包含具有時間相關(guān)性的有色噪聲，從而影響了已建立數(shù)學(xué)模型的正確性[7-10]。事實上，在實際應(yīng)用中如忽略觀測值時間相關(guān)性往往會得到不可靠的參數(shù)解及其精度[11-13]。文獻(xiàn)[14]詳細(xì)分析了忽略觀測值時間相關(guān)性得到的參數(shù)解精度往往過于樂觀且不符合實際。文獻(xiàn)[15]指出觀測值時間相關(guān)性將引起GPS長基線定位結(jié)果差異達(dá)2 cm，在高精度定位領(lǐng)域如變形監(jiān)測中不可忽略。文獻(xiàn)[16]研究了顧及時間相關(guān)性可以獲得符合實際的最小可探測粗差，從而幫助用戶在質(zhì)量控制應(yīng)用中得到更正確的結(jié)果。此外，在GNSS應(yīng)用中，文獻(xiàn)[17]證明了當(dāng)考慮觀測值時間相關(guān)性時，可以提高模糊度固定的成功率，在單頻觀測值中尤其明顯。

然而，在動態(tài)數(shù)據(jù)處理中，顧及觀測值時間相關(guān)性的平差方法雖然得到了一定程度的研究[1,18-19]，尤其是在GNSS導(dǎo)航領(lǐng)域[20-22]，但在統(tǒng)一性和實用性兩個問題上尚未完全解決。首先，當(dāng)前應(yīng)用最為普遍的是通過時間差分變換方式，將原始相關(guān)觀測值變換為獨(dú)立觀測值，從而達(dá)到消除觀測值時間相關(guān)性的目的[18,20,23]。事實上，類似的還有去相關(guān)變換法，如利用LDLT分解對觀測值的方差-協(xié)方差陣進(jìn)行變換[15]。遺憾的是，這些方法的特點以及相互之間的關(guān)系目前并沒有得到詳細(xì)論證，因此容易讓人產(chǎn)生困惑，造成混亂。其次，由于時間相關(guān)觀測的動態(tài)數(shù)據(jù)處理中的方差-協(xié)方差矩陣不再是一個分塊對角矩陣，單歷元或者法方程疊加的方法不再適用，直接求逆又會引起巨大的計算量[24]，有色噪聲模型往往又很難實時準(zhǔn)確地獲取[25-27]，因此在實際應(yīng)用時，需要有簡化的方法在能夠保證解算效率的同時又能有效顧及時間相關(guān)性。

因此，本文將深入研究幾種時間相關(guān)觀測的動態(tài)數(shù)據(jù)處理方法。首先，對時間相關(guān)觀測的動態(tài)數(shù)據(jù)處理方法進(jìn)行了擴(kuò)展和統(tǒng)一。利用極大驗后(MAP)估計原理，提出了顧及觀測值時間相關(guān)性動態(tài)解的第3種思路，并針對無歷元公用參數(shù)和有歷元公用參數(shù)兩類應(yīng)用情景，詳細(xì)推導(dǎo)了去相關(guān)變換法、差分變換法和MAP估計法這3種顧及時間相關(guān)性的動態(tài)解的解析公式，并重點探討了它們的適用范圍以及等價性。其次，導(dǎo)出了一種實用的時間相關(guān)觀測的動態(tài)數(shù)據(jù)處理方法，即針對嚴(yán)格考慮時間相關(guān)性會導(dǎo)致計算復(fù)雜且效率低的客觀問題，采用自相關(guān)函數(shù)(ACF)實時估計時間相關(guān)性，從而發(fā)展了一種簡化且能有效表征觀測值時間相關(guān)性影響的退化形式。最后，采用雙頻GPS試驗數(shù)據(jù)，對本文理論進(jìn)行了驗證。

1 無歷元公用參數(shù)的時間相關(guān)觀測模型

首先探討無歷元公用參數(shù)的時間相關(guān)觀測模型，例如采用GNSS偽距或者模糊度固定后的載波相位進(jìn)行動態(tài)定位。設(shè)連續(xù)K個歷元觀測值時間相關(guān)，對應(yīng)的觀測方程及其方差-協(xié)方差陣為

y=Ax+ε

(1a)

(1b)

1.1 去相關(guān)變換法

利用LDLT分解法對方差-協(xié)方差陣Qyy進(jìn)行變換，令分解形式為Qyy=UDUT，其中，U為單位下三角矩陣，對應(yīng)的遞歸公式為[28]

(2)

式中，Uij表示單位下三角矩陣的第i行第j列子矩陣；D=blkdiag(D1,D2,…,DK)為分塊對角矩陣。令LQyyLT=D，則

(3)

去相關(guān)變換法的基本思想是，通過對相關(guān)觀測值改造，得到一組新的獨(dú)立觀測量。具體的，對觀測方程式(1a)左乘矩陣L，得到新的觀測方程

(4a)

(4b)

(5)

(6)

(7)

1.2 差分變換法

將觀測方程式(1a)按前i-1個歷元和第i個歷元分塊表達(dá)，則前i-1個歷元的觀測方程為

(8)

(9)

(10)

(11)

(12)

(13)

(14)

1.3 極大驗后估計法

(15)

(16)

(17a)

(17b)

(18)

(19)

(20)

(21)

將式(20)和式(21)代入式(18)得式(14)，即證明了MAP法與差分變換法的動態(tài)解是等價的。類似的，MAP估計法不需要前后歷元的觀測值類型和維數(shù)相同，而差分變換法則需要。

2 有歷元公用參數(shù)的時間相關(guān)觀測模型

在某些應(yīng)用場景下，歷元間往往存在公用參數(shù)，例如，在GNSS長基線精密定位中，除了隨歷元變化的位置參數(shù)外，還含有未固定的模糊度參數(shù)和天頂對流層參數(shù)。此時，即使歷元間觀測值相互獨(dú)立時，也不能再用單歷元觀測值單獨(dú)平差的方法，而需要多個歷元聯(lián)合進(jìn)行平差。但傳統(tǒng)的多歷元平差方法因其需要不同歷元的觀測值，容易引起解算效率低的問題。因此，為了能夠保證解算結(jié)果不變的情況下，盡最大可能地提高解算效率，本節(jié)將探討有歷元公用參數(shù)的時間相關(guān)觀測模型的動態(tài)數(shù)據(jù)處理方法。設(shè)連續(xù)K個歷元的觀測值存在時間相關(guān)，除了隨歷元變化的參數(shù)x外，歷元間還存在公用參數(shù)ξ，則K個歷元的觀測方程為

y=Ax+Cξ+ε

(22)

對觀測方程式(22)左乘矩陣L，得到類似于式(4a)的去相關(guān)后的觀測方程

(23)

(24)

(25)

式中，將式(25)與式(24)的第2個方程進(jìn)行最小二乘準(zhǔn)則下融合，得法方程

(26)

采用最小二乘法方程約化且根據(jù)矩陣反演公式，推導(dǎo)得到第i個歷元的參數(shù)最小二乘估值為

(27)

3 基于ACF的退化形式

雖然上述解算公式嚴(yán)密準(zhǔn)確，但涉及大量的矩陣運(yùn)算，計算復(fù)雜且效率低。折中的方法是采用盡量簡潔且能有效刻畫時間相關(guān)性的協(xié)方差函數(shù)，從而在滿足計算效率的同時能符合實際時間相關(guān)誤差特征。假設(shè)各歷元觀測值的方差-協(xié)方差陣相同，利用ACF估計該時間相關(guān)系數(shù)，且相鄰歷元的時間相關(guān)系數(shù)為ρ，則原始方差-協(xié)方差陣為

(28)

(29)

4 算例與分析

為驗證上述理論公式推導(dǎo)的正確性和實用性，選取了同一測區(qū)同一時間段內(nèi)的兩條基線長度分別為11.4和42.5 km的雙頻相位GPS數(shù)據(jù)，分別命名為No.1和No.2。采樣間隔為1 s，時間長度為1 h，模糊度已事先固定。首先利用傳統(tǒng)的單歷元單獨(dú)平差方法，獲得定位結(jié)果，并采用ACF計算觀測值殘差的時間相關(guān)性，通過得到的時間相關(guān)系數(shù)可以判斷相應(yīng)觀測值是否存在時間相關(guān)性。圖1是兩條基線的雙頻相位觀測值時間相關(guān)系數(shù)，其中每條曲線代表一個雙差衛(wèi)星對。顯然，對于這兩條基線，時間相關(guān)性在一段時間內(nèi)都是顯著存在的。因此，采用顧及觀測值時間相關(guān)性的動態(tài)數(shù)據(jù)處理方法是必要的。

圖1 GPS觀測值時間相關(guān)系數(shù)Fig.1 Time correlation coefficients of GPS observations

為分析觀測值時間相關(guān)性的影響因素，統(tǒng)計了這兩條基線所有雙頻相位觀測值時間相關(guān)性小于0.3以及接近0時的平均時間間隔，這里0.3和0分別認(rèn)為時間相關(guān)性不顯著以及不存在。計算可得，當(dāng)基線長度為11.4 km時，觀測值時間相關(guān)性不顯著以及不存在的時間間隔平均為86和181 s；而當(dāng)基線長度變?yōu)?2.5 km時，相應(yīng)時間間隔平均為220和396 s。即在200～400 s之后，觀測值時間相關(guān)性才可以完全忽略。事實上，之所以不同基線長度的觀測值可以忽略時間相關(guān)性的時間間隔不同，是由于上述時間相關(guān)性主要是因為觀測值中存在一些難以通過經(jīng)驗?zāi)Ｐ透恼皡?shù)化吸收等的非模型化誤差所引起的[5-6,30]。顯然，當(dāng)基線長度越長，所受到的大氣誤差越嚴(yán)重，因而非模型化誤差也越嚴(yán)重，即造成了觀測值時間相關(guān)性也越強(qiáng)的現(xiàn)象。

接著，采用上述時間相關(guān)觀測模型中的3種方法分別進(jìn)行計算，即A、B和C分別表示去相關(guān)變換法、差分變換法和MAP估計法。統(tǒng)計了利用這些方法對No.1和No.2這兩條基線進(jìn)行解算后的定位誤差的均值和標(biāo)準(zhǔn)差，可以發(fā)現(xiàn)這3種方法的定位誤差都在厘米級，且他們的定位結(jié)果等價，具體統(tǒng)計結(jié)果如表1所示。因此，驗證了本文公式推導(dǎo)的正確性及3種方法等價性的論斷。

表1 兩組數(shù)據(jù)的定位統(tǒng)計結(jié)果

下面將以GNSS模糊度固定為例來研究考慮時間相關(guān)性的作用?；贚AMBDA方法，分別采用L1、L2和L1+L2進(jìn)行單歷元模糊度固定，并與真值比較分析模糊度固定成功率。以No.1為例，表2顯示的是不顧及和顧及時間相關(guān)性的模糊度固定成功率。當(dāng)考慮時間相關(guān)性時，模糊度固定的成功率顯著提高，平均提高約7%。其原因是考慮時間相關(guān)性的方差-協(xié)方差陣更加符合實際觀測數(shù)據(jù)情況，導(dǎo)致模糊度搜索空間與浮點解相容[17]。

圖2 No.1數(shù)據(jù)U方向上不顧及和顧及時間相關(guān)性的定位精度比較Fig.2 Positioning precisions of ignoring and considering time correlations in U direction of dataset No.1

表2No.1數(shù)據(jù)不顧及和顧及時間相關(guān)性的模糊度固定成功率

Tab.2AmbiguityresolutionsuccessratesofignoringandconsideringtimecorrelationsofdatasetNo.1(%)

5 結(jié) 論

本文圍繞解決現(xiàn)代大地測量數(shù)據(jù)處理中觀測值時間相關(guān)的問題，深入研究并提出了幾種時間相關(guān)觀測的動態(tài)解法，包括無歷元公用參數(shù)和有歷元公用參數(shù)兩類模型下3種方法的推導(dǎo)，即去相關(guān)變換法、差分變換法和MAP估計法，同時結(jié)合ACF，導(dǎo)出了它們的等價退化形式，得出以下結(jié)論：

(1) 推導(dǎo)出的3種顧及觀測值時間相關(guān)性的動態(tài)解的結(jié)果是等價的且各有特點，因此在實際使用中可選擇任意一種易于實現(xiàn)的方法。

(2) 導(dǎo)出的等價退化形式理論上簡化了歷元之間的關(guān)系，且有效提高了計算效率，可以應(yīng)用于實時觀測數(shù)據(jù)，便于編程計算。

(3) 算例表明，實際觀測值中往往存在時間相關(guān)性，利用顧及觀測值時間相關(guān)性的動態(tài)解能得到更為可靠和準(zhǔn)確的計算結(jié)果，如可以提高模糊度固定成功率。