時間相關觀測的幾種動態數據處理方法
李博峰 , 章浙濤
同濟大學測繪與地理信息學院, 上海 200092
收稿日期：2018-04-27；修回日期：2018-09-07
基金項目：國家自然科學基金（41574031；41622401）；上海市科技委員會科技創新行動計劃（17511109501；17DZ1100802；17DZ1100902）
第一作者簡介：李博峰(1983-), 男, 博士, 教授, 研究方向為多頻多模GNSS數據處理理論及應用新技術。E-mail:bofeng_li@tongji.edu.cn
摘要：現代大地測量數據中往往存在時間相關性，忽略觀測值時間相關性會影響參數估值的精度和可靠性。因此，本文研究了幾種時間相關觀測的動態數據處理方法。首先，對觀測值時間相關性的處理方法進行了擴展和統一。利用極大驗後估計原理，提出了第3種思路，並在無曆元公用參數和有曆元公用參數兩類模型下，詳細推導了去相關變換法、差分變換法和極大驗後估計法這3種思路的動態解，並對它們的特點和等價性進行了論證。其次，為了平衡實際應用中的解算效率且同時有效顧及時間相關性，利用自相關函數表達時間相關性，導出了相應退化形式。最後，通過GPS實測數據驗證了本文理論公式的正確性和實用性。
關鍵詞：時間相關 動態解 去相關變換 差分變換 極大驗後估計
Several Kinematic Data Processing Methods for Time-correlated Observations
LI Bofen , ZHANG Zhetao
Abstract: Time correlations always exist in modern geodetic data, and ignoring these time correlations will affect the precision and reliability of solutions. In this paper, several kinematic data processing methods for time-correlated observations are studied. Firstly, the method for processing the time-correlated observations is expanded and unified. Based on the theory of maximum a posteriori estimation, the third idea is proposed. Two types of situations with and without common parameters are both investigated by using the decorrelation transformation, differential transformation and maximum a posteriori estimation solutions. Besides, the characteristics and equivalence of above three methods are studied. Secondly, in order to balance the computational efficiency in real applications and meantime effectively capture the time correlations, the corresponding reduced forms based on the autocorrelation function are deduced. Finally, with GPS real data, the correctness and practicability of derived formulae are evaluated.
Key words: time correlation kinematic solution decorrelation transformation differential transformation maximum a posteriori estimation
在現代大地測量動態數據處理中，只要觀測值數量大於未知參數個數，即可建立數學模型並採用逐點單曆元的平差方法求得所需的動態解[1-2]。數學模型包括函數模型和隨機模型兩部分，前者表徵的是觀測值和未知參數之間的關係，而後者描述的是觀測值的精度以及相互之間的關係，兩者同等重要且缺一不可[3-4]。
通常，可以認為觀測誤差為高斯白噪聲，即假設觀測值之間不存在時間相關性。然而在實際觀測數據中，由於非模型化誤差的存在[5-6]，導致觀測值中往往包含具有時間相關性的有色噪聲，從而影響了已建立數學模型的正確性[7-10]。事實上，在實際應用中如忽略觀測值時間相關性往往會得到不可靠的參數解及其精度[11-13]。文獻[14]詳細分析了忽略觀測值時間相關性得到的參數解精度往往過於樂觀且不符合實際。文獻[15]指出觀測值時間相關性將引起GPS長基線定位結果差異達2cm，在高精度定位領域如變形監測中不可忽略。文獻[16]研究了顧及時間相關性可以獲得符合實際的最小可探測粗差，從而幫助用戶在質量控制應用中得到更正確的結果。此外，在GNSS應用中，文獻[17]證明了當考慮觀測值時間相關性時，可以提高模糊度固定的成功率，在單頻觀測值中尤其明顯。
然而，在動態數據處理中，顧及觀測值時間相關性的平差方法雖然得到了一定程度的研究[1, 18-19]，尤其是在GNSS導航領域[20-22]，但在統一性和實用性兩個問題上尚未完全解決。首先，當前應用最為普遍的是通過時間差分變換方式，將原始相關觀測值變換為獨立觀測值，從而達到消除觀測值時間相關性的目的[18, 20, 23]。事實上，類似的還有去相關變換法，如利用LDLT分解對觀測值的方差-協方差陣進行變換[15]。遺憾的是，這些方法的特點以及相互之間的關係目前並沒有得到詳細論證，因此容易讓人產生困惑，造成混亂。其次，由於時間相關觀測的動態數據處理中的方差-協方差矩陣不再是一個分塊對角矩陣，單曆元或者法方程疊加的方法不再適用，直接求逆又會引起巨大的計算量[24]，有色噪聲模型往往又很難實時準確地獲取[25-27]，因此在實際應用時，需要有簡化的方法在能夠保證解算效率的同時又能有效顧及時間相關性。
因此，本文將深入研究幾種時間相關觀測的動態數據處理方法。首先，對時間相關觀測的動態數據處理方法進行了擴展和統一。利用極大驗後(MAP)估計原理，提出了顧及觀測值時間相關性動態解的第3種思路，並針對無曆元公用參數和有曆元公用參數兩類應用情景，詳細推導了去相關變換法、差分變換法和MAP估計法這3種顧及時間相關性的動態解的解析公式，並重點探討了它們的適用範圍以及等價性。其次，導出了一種實用的時間相關觀測的動態數據處理方法，即針對嚴格考慮時間相關性會導致計算複雜且效率低的客觀問題，採用自相關函數(ACF)實時估計時間相關性，從而發展了一種簡化且能有效表徵觀測值時間相關性影響的退化形式。最後，採用雙頻GPS試驗數據，對本文理論進行了驗證。
1 無曆元公用參數的時間相關觀測模型
首先探討無曆元公用參數的時間相關觀測模型，例如採用GNSS偽距或者模糊度固定後的載波相位進行動態定位。設連續K個曆元觀測值時間相關，對應的觀測方程及其方差-協方差陣為
(1a)
(1b)
式中，y=[y1T, y2T, …, yKT]T；A=blkdiag(A1, A2, …, AK)；x=[x1T, x2T, …, xKT]T；ε=[ε1T, ε2T, …, εKT]T。當曆元間觀測值相互獨立時，即Qyiyj=0(i≠j)，可採用單曆元觀測值單獨平差的方法進行求解；而當曆元間觀測值時間相關時，即Qyiyj≠0，需要聯立多曆元觀測值聯合平差的方法進行求解。本文將分別採用去相關變換法、差分變換法和MAP估計法推導時間相關觀測的動態數據處理解法，並著重論證它們的等價性。
1.1 去相關變換法
利用LDLT分解法對方差-協方差陣Qyy進行變換，令分解形式為Qyy=UDUT，其中，U為單位下三角矩陣，對應的遞歸公式為[28]
(2)
式中，Uij表示單位下三角矩陣的第i行第j列子矩陣；D=blkdiag(D1, D2, …, DK)為分塊對角矩陣。令LQyyLT=D，則
(3)
去相關變換法的基本思想是，通過對相關觀測值改造，得到一組新的獨立觀測量。具體的，對觀測方程式(1a)左乘矩陣L，得到新的觀測方程
(4a)
式中， 。變換後觀測值的方差-協方差陣為
(4b)
式中，。由於D為分塊對角矩陣，顯然變換後各曆元觀測值相互獨立。令 表示所有i-1個曆元參數構成的向量，則所有i個曆元的觀測方程及其方差-協方差陣可表示為
(5)
式中， 。注意，此時變換後的第i個曆元觀測方程中，包含了所有i個曆元的參數。對應的最小二乘法方程為
(6)
式中， 。採用最小二乘法方程約化消去參數，並根據矩陣反演公式，推導出第i個曆元參數的最小二乘估值為
(7)
式中，為以參數估值為條件的觀測值的方差-協方差陣。
1.2 差分變換法
將觀測方程式(1a)按前i-1個曆元和第i個曆元分塊表達，則前i-1個曆元的觀測方程為
(8)
式中， 。方差-協方差陣式(1b)相應地表示為
(9)
式中，。
將第i個曆元觀測值改造為，則新觀測方程為
(10)
式中，。根據協方差傳播律，前i-1個曆元觀測值與改造後的第i個曆元觀測值的方差-協方差陣為
(11)
顯然，當時，觀測值與li相互獨立。此時，前i個曆元觀測方程及其方差-協方差陣為
(12)
式中，。對應的最小二乘法方程為
(13)
式中， 。採用最小二乘法方程約化消去參數，並根據矩陣反演公式，推導出第i個曆元參數的最小二乘估值為
(14)
式中，為以參數估值為條件的觀測值的方差-協方差陣。
與去相關變換法比較，可以證明 ，即差分變換法得到的參數估值與去相關變換法的參數估值等價。但差分變換法要求差分曆元的觀測值類型與維數相同，而去相關變換法不需要該條件。
1.3 極大驗後估計法
MAP估計理論的本質是用觀測值y對參數先驗信息進行改造，其改造程度由觀測值方差及其與參數的相關程度共同決定[4]。首先給出MAP估計的標準形式，令隨機變量和的先驗統計信息為
(15)
式中，E為期望算子。根據MAP理論，當隨機變量的採樣(實際觀測)y存在時，隨機變量對應的估值為[4]
(16)
式中，帶下劃線的變量和表示隨機變量；不帶下劃線的變量y表示隨機變量的採樣值。
本節採用MAP估計推導時間相關觀測的動態解，利用差分變換法中的分塊觀測方程，按照MAP標準形式進行變量對應，即 ，則改造後的第i個曆元觀測值為
(17a)
(17b)
式(17a)的含義不難理解，由於觀測值li與yi相關，當已知li的估值時，將必然對觀測值yi產生更新(改造)。此時，改造後的觀測值事實上已經顧及了觀測值之間的相關性，因此直接採用改造後的觀測值能得到與整體平差解等價的第i個曆元參數估值，即
(18)
由於上文已證明了去相關變換法與差分變換法等價，這裡只證明MAP法的動態解方程式(18)與差分變換法方程式(14)等價。可導出 ，代入式(17b)得
(19)
顧及差分變換法中且，則
(20)
此外，由式(17a)得，同時顧及差分變換法中，則
(21)
將式(20)和式(21)代入式(18)得式(14)，即證明了MAP法與差分變換法的動態解是等價的。類似的，MAP估計法不需要前後曆元的觀測值類型和維數相同，而差分變換法則需要。
2 有曆元公用參數的時間相關觀測模型
在某些應用場景下，曆元間往往存在公用參數，例如，在GNSS長基線精密定位中，除了隨曆元變化的位置參數外，還含有未固定的模糊度參數和天頂對流層參數。此時，即使曆元間觀測值相互獨立時，也不能再用單曆元觀測值單獨平差的方法，而需要多個曆元聯合進行平差。但傳統的多曆元平差方法因其需要不同曆元的觀測值，容易引起解算效率低的問題。因此，為了能夠保證解算結果不變的情況下，盡最大可能地提高解算效率，本節將探討有曆元公用參數的時間相關觀測模型的動態數據處理方法。設連續K個曆元的觀測值存在時間相關，除了隨曆元變化的參數x外，曆元間還存在公用參數ξ，則K個曆元的觀測方程為
(22)
式中，C=[C1T, C2T, …, CKT]T，對應的方差-協方差陣形式與式(1b)相同。由於上節已證明了去相關變換法、差分變換法和MAP估計法的等價性，本文利用去相關變換法進行推導。
對觀測方程式(22)左乘矩陣L，得到類似於式(4a)的去相關後的觀測方程
(23)
式中，顯然，變換後各曆元觀測值相互獨立，變換後的第i個曆元觀測方程包含了所有i個曆元的參數。為進行法方程約化，對i個曆元的觀測方程分塊及其方差-協方差陣重整為
(24)
式中，，其餘變量與式(5)中定義相同。顯然，採用前i-1個曆元觀測值可計算參數的估值，其中ξi-1表示採用前i-1個曆元觀測值計算的參數ξ的估值。令採用前i-1個曆元觀測值求得的參數估值為
(25)
式中，將式(25)與式(24)的第2個方程進行最小二乘準則下融合，得法方程
(26)
式中，。注意不同，表示採用所有i個曆元觀測值求得的參數的估值。
採用最小二乘法方程約化且根據矩陣反演公式，推導得到第i個曆元的參數最小二乘估值為
(27)
式中，。
3 基於ACF的退化形式
雖然上述解算公式嚴密準確，但涉及大量的矩陣運算，計算複雜且效率低。折中的方法是採用儘量簡潔且能有效刻畫時間相關性的協方差函數，從而在滿足計算效率的同時能符合實際時間相關誤差特徵。假設各曆元觀測值的方差-協方差陣相同，利用ACF估計該時間相關係數，且相鄰曆元的時間相關係數為ρ，則原始方差-協方差陣為
(28)
式中，符號「⊗」代表克羅內克積，即滿足A⊗B= 為矩陣A中的元素。上述方差-協方差陣形式也可以稱之為部分延續模式[29]。對應的分解LRLT=D，其中L= ρ2])。
以觀測方程式(1a)為例，利用去相關變換法進行推導。具體的，對觀測方程式(1a)左乘矩陣L⊗I。對應去相關後的觀測方程式(4a)中的變量為 。變換後觀測值的方差-協方差陣為，其中。對應觀測方程式(5)中，Ei=[0, …，0, -ρAi-1]，將相應的變量代入式(7)，得
(29)
式中，為以參數估值為條件的觀測值的方差-協方差陣。顯然，基於ACF的退化形式極其簡潔，即Ei、GTCi和可以被利用ACF計算的ρ簡化，而且退化為，簡化了曆元之間的關係。因此，該退化形式既能符合實際時間相關誤差特徵又顯著提高了計算效率。
4 算例與分析
為驗證上述理論公式推導的正確性和實用性，選取了同一測區同一時間段內的兩條基線長度分別為11.4和42.5km的雙頻相位GPS數據，分別命名為No.1和No.2。採樣間隔為1s，時間長度為1h，模糊度已事先固定。首先利用傳統的單曆元單獨平差方法，獲得定位結果，並採用ACF計算觀測值殘差的時間相關性，通過得到的時間相關係數可以判斷相應觀測值是否存在時間相關性。圖 1是兩條基線的雙頻相位觀測值時間相關係數，其中每條曲線代表一個雙差衛星對。顯然，對於這兩條基線，時間相關性在一段時間內都是顯著存在的。因此，採用顧及觀測值時間相關性的動態數據處理方法是必要的。
圖 1 GPS觀測值時間相關係數Fig. 1 Time correlation coefficients of GPS observations
為分析觀測值時間相關性的影響因素，統計了這兩條基線所有雙頻相位觀測值時間相關性小於0.3以及接近0時的平均時間間隔，這裡0.3和0分別認為時間相關性不顯著以及不存在。計算可得，當基線長度為11.4km時，觀測值時間相關性不顯著以及不存在的時間間隔平均為86和181s；而當基線長度變為42.5km時，相應時間間隔平均為220和396s。即在200~400s之後，觀測值時間相關性才可以完全忽略。事實上，之所以不同基線長度的觀測值可以忽略時間相關性的時間間隔不同，是由於上述時間相關性主要是因為觀測值中存在一些難以通過經驗模型改正及參數化吸收等的非模型化誤差所引起的[5-6, 30]。顯然，當基線長度越長，所受到的大氣誤差越嚴重，因而非模型化誤差也越嚴重，即造成了觀測值時間相關性也越強的現象。
接著，採用上述時間相關觀測模型中的3種方法分別進行計算，即A、B和C分別表示去相關變換法、差分變換法和MAP估計法。統計了利用這些方法對No.1和No.2這兩條基線進行解算後的定位誤差的均值和標準差，可以發現這3種方法的定位誤差都在釐米級，且他們的定位結果等價，具體統計結果如表 1所示。因此，驗證了本文公式推導的正確性及3種方法等價性的論斷。
表 1 兩組數據的定位統計結果Tab. 1 Positioning statistical results of the two datasets
mm數據均值為進一步證明採用時間相關觀測動態解的必要性，下面將比較不顧及時間相關性和顧及時間相關性的平差結果。因3種顧及時間相關性的方法A、B和C等價，因此以方法A為例與不顧及時間相關性的方法作比較。由於兩類方法的定位誤差結果不是特別大，這與其他相關研究的結論是相同的[14, 30-31]，因此沒有展示它們的定位結果。但對一個完整的平差結果來說，參數估計固然重要，但精度評定同樣不可或缺，因此筆者著重分析這兩類方法對定位精度的影響。推導可得K個曆元的基線解方差-協方差陣的精度評定公式 (K-1)，其中為K個曆元內的基線解，R定義同式(28)，eK為K行全1列向量。顯然，當忽略時間相關性時，R變為K階單位矩陣，此時滿足 (K-1)[14]。以No.1為例，設曆元窗口為60s，通過計算可得U方向上不顧及和顧及時間相關性的定位精度，如圖 2所示。可以發現，當忽略觀測值時間相關性時，其基線解的精度值小於顧及時間相關性的精度值。這表明了當忽略觀測值時間相關性時，本文得到的平差結果是過於樂觀且不符合實際的[6, 14]。因此，為了獲得更可靠的基線解，需要考慮觀測值時間相關性。
圖 2 No.1數據U方向上不顧及和顧及時間相關性的定位精度比較Fig. 2 Positioning precisions of ignoring and considering time correlations in U direction of dataset No.1下面將以GNSS模糊度固定為例來研究考慮時間相關性的作用。基於LAMBDA方法，分別採用L1、L2和L1+L2進行單曆元模糊度固定，並與真值比較分析模糊度固定成功率。以No.1為例，表 2顯示的是不顧及和顧及時間相關性的模糊度固定成功率。當考慮時間相關性時，模糊度固定的成功率顯著提高，平均提高約7%。其原因是考慮時間相關性的方差-協方差陣更加符合實際觀測數據情況，導致模糊度搜索空間與浮點解相容[17]。
表 2 No.1數據不顧及和顧及時間相關性的模糊度固定成功率Tab. 2 Ambiguity resolution success rates of ignoring and considering time correlations of dataset No.1
(%)模糊度固定策略L1L2L1+L2不顧及時間相關性47.258.385.7顧及時間相關性55.065.492.2
5 結論
本文圍繞解決現代大地測量數據處理中觀測值時間相關的問題，深入研究並提出了幾種時間相關觀測的動態解法，包括無曆元公用參數和有曆元公用參數兩類模型下3種方法的推導，即去相關變換法、差分變換法和MAP估計法，同時結合ACF，導出了它們的等價退化形式，得出以下結論：
(1) 推導出的3種顧及觀測值時間相關性的動態解的結果是等價的且各有特點，因此在實際使用中可選擇任意一種易於實現的方法。
(2) 導出的等價退化形式理論上簡化了曆元之間的關係，且有效提高了計算效率，可以應用於實時觀測數據，便於編程計算。
(3) 算例表明，實際觀測值中往往存在時間相關性，利用顧及觀測值時間相關性的動態解能得到更為可靠和準確的計算結果，如可以提高模糊度固定成功率。
【引文格式】李博峰, 章浙濤. 時間相關觀測的幾種動態數據處理方法. 測繪學報，2018，47(12)：1563-1570. DOI: 10.11947/j.AGCS.2018.20180192
權威 | 專業 | 學術 | 前沿
微信投稿郵箱 | song_qi_fan@163.com
微信公眾號中搜索「測繪學報」，關注我們，長按上圖二維碼，關注學術前沿動態。
歡迎加入《測繪學報》作者QQ群： 297834524
進群請備註：姓名+單位+稿件編號