對流層時間序列的Kalman濾波組合方法

李 敏，徐天河

對流層時間序列的Kalman濾波組合方法

李 敏1，徐天河2，3，4

(1.長安大學 地質(zhì)工程與測繪工程學院，西安 710054；2.地理信息工程國家重點實驗室，西安 710054； 3.宇航動力學國家重點實驗室，西安 710043；4.西安測繪研究所，西安 710054)

由于目前IGS提供的對流層最終產(chǎn)品滯后約一星期到四星期，即使是超快速產(chǎn)品也滯后約3小時，這嚴重影響了目前實時氣象的應用與需求，而且對流層產(chǎn)品的連續(xù)和實時組合的相關研究較少，大多集中在歐洲的試點項目?；诖?，本文主要研究了一種對流層時間序列的短期高精度連續(xù)組合方法，利用各個分析中心給出的對流層時間序列，利用抗差Kalman濾波和最小二乘方差分量估計的原理，進行GPS或VLBI的短期時間序列的實時的或事后的連續(xù)組合。先估計各個分析中心產(chǎn)品偏差，同時計算出其權因子，然后利用Kalman濾波技術進行產(chǎn)品組合可以獲得滯后的或?qū)崟r的組合解及其標準差，組合解的平均精度達到0.85 mm。

對流層；時間序列；組合；抗差Kalman濾波；最小二乘方差分量估計

0 引言

對流層延遲是全球衛(wèi)星導航系統(tǒng)(global navigation satellite system，GNSS)定位的主要誤差源之一，其隨著衛(wèi)星高度角的降低而增大，在低高度角的情況(10°以下)，中低緯度地區(qū)對流層延遲誤差可達20～40 m[1]。目前國際GNSS服務組織(International GNSS Service，IGS)和國際VLBI服務組織(International VLBI Service，IVS)等機構均提供了對流層天頂相位延遲產(chǎn)品ZTD，另外歐洲大地測量參考框架網(wǎng)[2](The International Association of Geodesy Reference Frame Sub-Commission for Europe，EUREF)和歐洲GPS氣象觀測方案[2](targeting optimal use of GPS humidity measurements in meteorology，TOUGH)也提供了滯后的和實時的對流層產(chǎn)品，這些產(chǎn)品均來自他們各自分析中心的法方程層面或參數(shù)層面的組合解。

IGS提供的對流層產(chǎn)品有兩種：最終和超快速對流層產(chǎn)品，分別由各個分析中心采用最終的和超快速的軌道和鐘差產(chǎn)品計算而來，最初由Gendt等人通過組合各家分析中心的SINEX文件得到最終解[3]。EUREF永久跟蹤網(wǎng)有超過100個站，且分布密集、均勻，其觀測數(shù)據(jù)由12家分析中心解算生成SINEX形式的站坐標和對流層產(chǎn)品，再綜合成最終的周解EUREF產(chǎn)品。這兩種情況下的事后解已無法滿足目前實時氣象應用的需求。2005年，國際上COST-716項目通過實時處理EUREF觀測網(wǎng)絡的模式，獲得了間隔為1 h 45 min的對流層實時產(chǎn)品。針對連續(xù)、實時的對流層產(chǎn)品綜合，目前國內(nèi)外研究涉及較少，本章詳細研究了一種對流層時間序列的高精度連續(xù)、實時組合方法，利用抗差Kalman濾波和最小二乘方差分量估計的原理，進行對流層短期時間序列的組合。本文以EUREF產(chǎn)品實例驗證該方法的有效性。

1 算法描述

1.1 對流層時間序列偏差估計

假設各個分析中心的解是連續(xù)的，將解化成以天或周為單位的模塊，而每天又包含小時為單位的對流層解。每個分析中心選取的站、采用的軟件、解的策略、高度截止角、對流層總延遲間隔、采樣間隔等不盡相同，給出的對流層產(chǎn)品之間不可避免地存在系統(tǒng)誤差，這可通過估計不同分析中心的時間序列偏差加以解決。而采樣間隔的差異可通過高精度插值方法來進行統(tǒng)一。

(1)

未知量為組合的時間序列解和組合序列與單個時間序列的偏差，顯然單個時間序列估值與這些未知量是線性相關的，因此，對于第k個時間序列有：

(2)

從(2)式可以看出時間序列偏差可以理解為不同分析中心處理策略不一致引起的系統(tǒng)誤差，這附加系統(tǒng)參數(shù)的平差模型的合理性將在后面得到檢驗。

假設y所包含的觀測噪聲向量是隨機的(白噪聲)，則所有分析中心的形如(1)式的方程表達成線性化的Gauss-Markov模型為：

(3)

(4)

式(3)和(4)中，E和D分別表示期望和方差向量；Y是N×1維向量的組合解；b是K×1維向量的時間序列偏差值；AY和Ab分別是Y和b的(K·N)×N維和(K·N)×K維的設計矩陣，寫成表達式的形式有：

(5)

式(5)中，INN是N×N維的單位矩陣，eN是N維全為1的列向量，Ab中省略元素均為0。

利用Kalman濾波來進行時間序列的連續(xù)組合，用來建模時間序列偏差的隨機游走模型[4]表示為：

(6)

在上一步組合中的時間序列偏差估值可當作下一次組合的先驗信息，寫成表達式：

bi,i-1=bi-1

(7)

式(7)中，bi,i-1為第i步預測的時間序列偏差向量。bi的協(xié)方差矩陣Qbi的增值過程描述如下：

Qbi,i-1=Qbi-1+Q

(8)

式(8)中：

(9)

顧及式(7)與式(8)，可以用Kalman濾波來計算組合ZTD解Yi和在第i步的時間序列偏差bi的濾波估計值：

(10)

(11)

式(10)和式(11)的時間序列偏差的觀測更新方程為

(12)

(13)

式(13)中，增益矩陣為

(14)

殘差向量為

vi=yi-AYYi-Abbi,i-1

(15)

這樣，迭代開始，最初的K個先驗時間序列偏差全設為0，迭代時，檢查每一步的vi，由于引入了系統(tǒng)誤差參數(shù)(時間序列偏差)，該殘差可看作額外殘差，而額外殘差與粗差有很大的關系[5]，因此以2.5σ為標準來檢驗額外殘差從而探測超限值。

注意從式(13)、式(14)可以看出要成功得出某一歷元的時間序列偏差，至少需要兩個分析中心和至少兩個共同歷元的解，一旦一個塊中某個歷元出現(xiàn)不到兩個分析中心的殘差滿足要求(不超限)，將無法組合下去，這時可利用的就只能是前面歷元獲得的先驗信息對當前結果進行預報。

額外殘差的產(chǎn)生是由于各個分析中心所用軟件以及處理策略的不同引起的。將其標準化有：

(16)

(17)

式(17)中，k0、k1在本文分別取1.5和2.5。在這里當殘差在2.5σ以內(nèi)時，才利用構造的等價權重新確定驗后標準差。超限時解和標準差全部舍去。

對于式(10)、式(11)，其并不是對所有時間序列組合塊都是有效的，也就是說可能有模型誤差的產(chǎn)生，這時有必要進行平差模型的顯著性檢驗。

對于殘差向量vi，其協(xié)方差陣：

(18)

(19)

以此為統(tǒng)計量可用來檢驗原假設和備擇假設：

(20)

1.2 對流層時間序列權估計

在2.1節(jié)的組合中，各個分析中心的權看作是一樣的。然而，前面已經(jīng)提到，每個分析中心處理策略的不一致將導致解之間產(chǎn)生相對偏差。因此，有必要對各個分析中心重新定權。

這里用到最小二乘方差分量估計的思想[8-10]。假設第k個分析中心權為wk，這些權將作為觀測值的協(xié)方差矩陣中的額外未知因子。

將向量Y和b組合在一起變成(K+N)×1維向量x，那么Ab和AY變成(K·N)×(K+N)維矩陣A，則：

(21)

如果把權因子看作協(xié)方差矩陣Qy的方差分量，那么觀測值的協(xié)方差矩陣的隨機模型可寫作：

(22)

式(22)中，Qk是含有N個對角元素的(K·N)×(K·N)維對角矩陣，對角元素為第k個時間序列ZTD估值的方差，其余元素為0，表達式如下：

Qk=

(23)

因此，如果將Qk看作協(xié)因數(shù)矩陣，式(21)可重寫為：

(24)

這樣可用最小二乘方差分量估計原理解出wk，根據(jù)文獻[10]，解是無偏的且具有最小方差分量的，限于篇幅，解算過程不再介紹。式(24)的最小二乘解為[10]：

(25)

式(25)中，

(26)

(27)

同時間序列偏差一樣，各個塊之間權的獲取也是連續(xù)的。對于s個時間序列塊，則有：

(28)

(29)

再次應用Kalman濾波，第i-1和i次(塊)組合的增值方程為

(30)

權的預測值為：wi,i-1=wi-1，Qwi-1的增值過程為：

(31)

那么類似于式(25)新的權因子的法方程為：

(32)

式(32)中，

(33)

(34)

式(32)不同式(25)的是，它考慮了上一步組合步驟，即考慮了先驗信息，因此它采用的是序貫或濾波組合方式計算出時間序列權因子。

2 步驟總結

對流層時間序列的連續(xù)組合步驟可總結如下：

(1)檢查文件和設置先驗信息。檢查各個時間序列文件，剔除解的標準差大于10mm的觀測值，各個對流層時間序列先驗偏差向量全設為0，先驗權設為1。

(2)時間序列組合和偏差估計，利用權因子修正觀測值的協(xié)方差矩陣，并作為當前濾波的先驗信息。在每次迭代時檢查殘差向量，根據(jù)式(17)來重新計算驗后標準差，并剔除超限值，計算時間序列偏差值和組合解。

(3)時間序列權估計，利用先驗信息和殘差向量來估計時間序列的權因子。

(4)時間序列偏差向量和權及其相應的協(xié)方差陣用來作為下一次連續(xù)組合步驟的先驗信息，重復(2)、(3)步驟，直到滿足收斂條件。

(5)結果輸出。

3 算例分析

由于只能獲取個別IGS分析中心(2013年少于4個)的對流層解，因此本文選取EUREF的分析中心解來驗證本文方法的正確性，筆者自編程序組合計算。

EUREF永久跟蹤網(wǎng)提供了每小時的對流層解，由各個分析中心解的時間序列組合得到。為了驗證本文方法的正確性，采取它的各個分析中心的時間序列進行組合，限于篇幅，本文只選取一個站MATE(Italy)的解來計算。時間為GPS周1 722～1 724，分析中心有ASI、MUT、BKG、BEK、RGA、UPA[11](別的分析中心沒有該站ZTD解)。它們的統(tǒng)計信息如下：

表1 各個AC統(tǒng)計信息表

注：AC表示分析中心；GIP和BSW表示GIPSY和Bernese軟件；NET表示網(wǎng)解；wN表示濕的Niell模型

從表1可以看出，本文選取的分析中心處理策略基本一致，測站個數(shù)略有差異，因此各個時間序列存在相關性。需要指出的是：除了ASI與別的分析中心在多數(shù)個方面有所差別外，其它分析中心都大致相同。

從2007年IGS開始使用新的ZTD產(chǎn)品，表現(xiàn)在GPS天線相位中心標準的改變[2]。由于采用更高精度的最終軌道、鐘差等產(chǎn)品，ZTD解的精度有所提高，六個分析中心所有解標準差均小于10 mm。組合時每天解為一個塊，所有塊連續(xù)組合。由于各個分析中心各天解精度基本一致，可以認為每天的時間序列偏差在其平均值周圍變化很小，在這里設σb和σw分別為0.3m和10-5，這樣設置能迭代4～5步后收斂。時間序列偏差和兩類權因子的估計結果如圖1-3所示：

圖1 時間序列偏差

圖2 時間序列權

圖3 不考慮先驗信息的時間序列權

圖1反映了時間序列偏差的變化，各個時間序列偏差值隨時間變化不大，基本都在某一值附近波動，但它們之間差別較大，ASI的偏差值最大，在1.5～2mm的比較多，MUT的最小，基本在±1mm，同時可以看出各個分析中心之間是存在一定的系統(tǒng)誤差的，需進行改正，不然對組合解會產(chǎn)生破壞性的影響。圖2反映了時間序列權的變化，由于是連續(xù)組合，經(jīng)過一段時間后，時間序列權排列錯落有致，在一定程度上反映了它們解的精度的總體變化；圖3是不考慮先驗信息的權的變化，和圖2比較可以看出他們是基本一致的，權因子的平均變化趨勢是，ASI的權偏大，RGA的偏小，這和它們的驗后標準差以及殘差有關，說明該時間序列權能總體反映解的長期精度。

從圖1和圖2看出時間序列偏差與權并不是完全一致，因此評定時間序列的精度要綜合考慮殘差和標準差等多方面信息。

組合解和EUREF解的差值以及組合解標準差如圖4、5所示：

圖4 組合解和EUREF解差值

圖5 組合解的標準差

較差和組合解的標準差的統(tǒng)計信息如表2所示：從圖4和表2可以看出，本文組合解和EUREF解差別在±1.2mm以內(nèi)，且已沒有系統(tǒng)差，這主要是因為都采用了時間序列組合。組合解的標準差基本在1mm。

表2 較差和標準差統(tǒng)計表

圖6 模型有效性檢驗

本文的組合模型并不是在每塊都是有效的，以5%的顯著水平計算超限值kα=270。從圖5可以看出在第三個塊有模型誤差產(chǎn)生，從第四塊之后模型趨于穩(wěn)定，且沒有模型誤差產(chǎn)生，說明該模型生效需要一定的連續(xù)組合步驟。

4 結束語

本文詳細探討了各個分析中心ZTD時間序列的Kalman濾波組合問題，該方法可進行事后的或?qū)崟rZTD組合，核心在于利用抗差Kalman濾波實時估計各個時間序列的偏差。由于各個分析中心解的精度參差不齊，也存在明顯的系統(tǒng)誤差，因此有必要對這些解進行偏差補償和加權組合。本文通過對EUREF的六家分析中心驗證計算，得到的組合解的平均標準差在0.85mm，和EUREF最終解的的一致性水平在0.32mm。本文提出的方法也可拓展到GPS和VLBI技術間的對流層產(chǎn)品組合。實時ZTD產(chǎn)品估計及組合將對氣象方面的研究有重要的參考和借鑒價值[12]。

[1] 張雙成，張鵬飛，范鵬飛.GPS對流層改正模型的最新進展及對比分析[J].大地測量與地球動力學，2012，32(2)：91-95.

[2]SUNGH，BYUN，YOAZE，etal.ANewTypeofTroposphereZenithPathDelayProductoftheInternationalGNSSService[J].JournalofGeodesy，2009，83(7):367-373.

[3]GENDTG.IGSCombinationofTroposphericEstimates[EB/OL].[2015-04-28].http://www.oso.chalmers.se/users/kge/cost/dwl/igs_wg_99.pdf.

[4]GELBA.AppliedOptimalEstimation[M].Cambridge，Massachusetts，USA:TheMITPress，1974.

[5] 張勤，張菊清.近代測量數(shù)據(jù)處理與應用[M].北京：測繪出版社，2010.

[6] 楊元喜.抗差估計理論及其應用[M].北京：八一出版社，1993.

[7]TIBERIUSC.RecursiveDataProcessingforKinematicGPSSurveying[D].Delft，Netherland:DelftUniversityofTechnology，1998.

[8]KOCHK.ParameterEstimationandHypothesisTestinginLinearModels[M].BerlinHeidelberg:Springer，1999.

[9]TEUNISSENPJG.Least-squaresVarianceComponentEstimation[J].JournalofGeodesy，2008，82(11):65-82.

[10]TEUNISSENPJG.TowardsaLeast-squaresFrameworkforAdjustingandTestingofbothFunctionalandStochasticModels[EB/OL].[2015-04-28].http://saegnss1.curtin.edu.au/Publications/2004/Teunissen2004Towards.pdf.

[11]HABRICHH，GOLTZM，WIESENSARTERE.GPS-GLONASS-GalileoTrackingDataandProducts[DB/OL].(2013-02-11)[2015-04-28].http://igs.bkg.bund.de/.

[12]KRüGELM，THALLERD，TESMERV，etal.TroposphericParameters:CombinationStudiesBasedonHomogeneousVLBIandGPSData[J].JournalofGeodesy，2007，81(6):515-527.

Research on the Combination of Troposphere Time Series with Kalman Filtering

LIMin1，XUTian-he2，3，4

(1.School of Geology Engineering and Surveying，Chang’an University，Xi’an 710054，China； 2.State Key Laboratory of Geo-information Engineering，Xi’an 710054，China； 3.State Key Laboratory of Astronautic Dynamics，Xi’an 710043,China； 4.Xi’an Research Institute of Surveying and Mapping，Xi’an 710054，China)

as the final troposphere products provided by IGS were delayed for at least one week，even though the ultro-rapid products also delayed for three hours，they all affected the real-time meteorological application，and also it’s rare to see the research papers on sequential or real-time combination for troposphere time series，most were focused on the experimental projects in Europe.Thus，on this paper the method of high precision and sequential combination for short-term troposphere time series provided by several analysis centers was discussed in detail，with the technique of robust Kalman filter and least-squares variances component estimation，a real-time or lagging sequential combination of GPS or VLBI short-term time series can be manipulated.First the bias and weight factor for each analysis center solution were estimated，and then the combined solution and its standard deviation can be obtained by Kalman filtering，the average accuracy can approach 0.85mm.

troposphere；time series；combination；robust Kalman filtering；least-squares variances component estimation

李敏，徐天河.對流層時間序列的Kalman濾波組合方法[J].導航定位學報,2015,3(3):79-84.(LI Min, XU Tian-he.Research on the Combination of Troposphere Time Series with Kalman Filtering[J].Journal of Navigation and Positioning,2015,3(3):79-84.)

10.16547/j.cnki.10-1096.20150316.

2015-05-18

國家自然科學基金(41174008)、宇航動力學國家重點實驗室開放基金(2014ADL-DW0101)。

李敏(1989—)，男，安徽桐城人，碩士生，現(xiàn)主要從事GNSS精密定位研究。

P228

A

2095-4999(2015)-03-0079-06