結(jié)合總體最小二乘的自適應(yīng)抗差濾波算法

劉 青 崔先強(qiáng) 宋迎春

(中南大學(xué)地球科學(xué)與信息物理學(xué)院，湖南 長沙410083)

0 引言

在動(dòng)態(tài)導(dǎo)航定位中，由于受到電離層、對(duì)流層、衛(wèi)星軌道以及載體受大氣摩擦和本身機(jī)動(dòng)等因素的影響，觀測(cè)噪聲和系統(tǒng)噪聲一般都為時(shí)間相關(guān)的有色噪聲。而標(biāo)準(zhǔn)Kalman濾波則要求觀測(cè)噪聲和系統(tǒng)噪聲為互不相關(guān)的高斯白噪聲，如在有色噪聲條件下繼續(xù)使用標(biāo)準(zhǔn)Kalman濾波必將會(huì)大大降低狀態(tài)參數(shù)估值的精度和可靠性，嚴(yán)重時(shí)可能會(huì)導(dǎo)致濾波發(fā)散。為了解決這一問題，國內(nèi)外學(xué)者提出了多種控制有色噪聲影響的算法，主要有矩陣增廣[1-3]和函數(shù)模型擬合[4-5]，后來，楊元喜等人又提出了一系列的自適應(yīng)抗差濾波算法[6-10]。以上這些算法都能有效地減弱或控制有色噪聲的影響，提高動(dòng)態(tài)導(dǎo)航定位的精度和可靠性，尤其是自適應(yīng)抗差濾波，不僅能夠控制有色噪聲的影響，而且能夠同事抵制觀測(cè)和動(dòng)力學(xué)模型異常的影響。

在實(shí)際的數(shù)據(jù)處理過程中發(fā)現(xiàn)，使用狀態(tài)預(yù)報(bào)值計(jì)算出的觀測(cè)方程系數(shù)矩陣往往會(huì)存在一定的偏差，這是因?yàn)檩d體的運(yùn)動(dòng)狀態(tài)總是在不斷變化中，使得預(yù)先給出的動(dòng)力學(xué)模型無法精確地描述載體的實(shí)際運(yùn)動(dòng)情況。此時(shí)，如果直接忽略觀測(cè)系數(shù)矩陣偏差的影響，必將會(huì)在一定程度上降低濾波結(jié)果的精度和可靠性。針對(duì)這種觀測(cè)方程系數(shù)矩陣存在偏差的情況，20世紀(jì)80年代總體最小二乘(total least squares,TLS)方法得以提出，該方法是一種能同時(shí)顧及觀測(cè)值誤差和模型系數(shù)矩陣偏差的數(shù)學(xué)方法，是最小二乘的擴(kuò)展，在系數(shù)矩陣存在偏差時(shí)，總體最小二乘解比最小二乘解更為真實(shí)可靠。關(guān)于TLS的研究最早可以追溯到19世紀(jì)，1877年Adcock[11]研究單變量直線擬合問題。隨后，Pearson[1 2]、 Koopmans[13]、Madansky[14]和 York[15]等人均對(duì)其進(jìn)行了研究和拓展。到1980年Golub和Van Loan才從數(shù)值分析的觀點(diǎn)首次對(duì)這種方法進(jìn)行了整體分析，并正式稱之為整體最小二乘[16]。此后，由于奇異值分解（singular value decomposition，SVD）的引入以及越來越多的應(yīng)用領(lǐng)域需要顧及系數(shù)矩陣偏差，許多學(xué)者對(duì)總體最小二乘進(jìn)行了更為深入的研究[17-20]。近幾年來，測(cè)繪學(xué)科內(nèi)的許多學(xué)者也開始了對(duì)總體最小二乘的研究，并在測(cè)量數(shù)據(jù)處理中得到了成功應(yīng)用[21-24]。

但是，現(xiàn)有的研究成果基本上都是針對(duì)單一存在有色噪聲或觀測(cè)系數(shù)矩陣偏差而進(jìn)行的，而在動(dòng)態(tài)導(dǎo)航定位中有色噪聲和系數(shù)矩陣偏差往往同時(shí)存在?；诖?，本文嘗試將現(xiàn)有的總體最小二乘方法和自適應(yīng)抗差濾波算法相結(jié)合，以便能夠同時(shí)減弱或控制有色噪聲和觀測(cè)系數(shù)矩陣偏差的影響，提高動(dòng)態(tài)導(dǎo)航定位的精度和可靠性。

1 自適應(yīng)抗差濾波

在動(dòng)態(tài)導(dǎo)航定位中，標(biāo)準(zhǔn)Kalman濾波的函數(shù)模型表示如下：

其中，Lk為觀測(cè)向量；Ak為觀測(cè)設(shè)計(jì)矩陣；Xk為歷元k時(shí)刻的狀態(tài)向量；Φk，k-1為狀態(tài)轉(zhuǎn)移矩陣；vk，ωk分別為觀測(cè)噪聲和狀態(tài)噪聲。

式中X?k-1指的是第tk-1歷元的狀態(tài)參數(shù)向量估值。因此，在tk歷元狀態(tài)預(yù)測(cè)信息向量和觀測(cè)向量的誤差方程可寫為：

式中Vk為觀測(cè)殘差向量，X?k為tk第歷元的狀態(tài)參數(shù)估值，Σk和Pk分別為觀測(cè)向量的協(xié)方差矩陣和權(quán)矩陣，為狀態(tài)預(yù)測(cè)向量的權(quán)矩陣。

自適應(yīng)抗差濾波原則為[25]：

自適應(yīng)抗差濾波的解向量為：

其中，c0和c1的取值范圍分別為1.0～1.5和3.0～4.5。統(tǒng)計(jì)量為：

式中的c0和c1取值可與式（9）中相同。

2 總體最小二乘算法求解系數(shù)矩陣偏差

為了能在自適應(yīng)抗差濾波的基礎(chǔ)上進(jìn)一步提高狀態(tài)參數(shù)估值的精度，必須先將觀測(cè)系數(shù)矩陣偏差求解出來。因此，我們先用總體最小二乘算法求解出觀測(cè)系數(shù)矩陣偏差EA，并對(duì)觀測(cè)方程系數(shù)矩陣進(jìn)行修正。

總體最小二乘問題有多種解法，本文中采用一種常見的迭代解法[28]。

當(dāng)系數(shù)矩陣Ak存在偏差時(shí)，觀測(cè)方程應(yīng)改寫為[29]：

若觀測(cè)向量殘差為V?k，同時(shí)系數(shù)矩陣 Ak的偏差為E?A，則上式的誤差方程可寫為：

式中，In和 It分別為 n 和 t階單位矩陣，vec（EA）是將矩陣 EA按列拉直得到的列向量化函數(shù)，vec（EA）∈R（n×t）×1；?表示克羅內(nèi)克(Kronecker)積。

總體最小二乘估計(jì)準(zhǔn)則為：

其中 K 為拉格朗日乘數(shù)。為求 Ф 的極小值，將其分別對(duì)V?k、vec（E?A）和X?k求一階導(dǎo)數(shù)，并令其為零：

由上式可見，γ是最小特征值，參數(shù)的總體最小二乘解是在γ=V?Tk V?k+vec（E?A）Tvec（E?A）=min≥0 時(shí)的特征向量[30]。 因此，由式（25）可得：

上式即為參數(shù)的總體最小二乘解，采用迭代求解。

計(jì)算步驟為：

由以上迭代計(jì)算得出參數(shù)估值X?k后，回代（22）式即可計(jì)算出系數(shù)矩陣偏差EA的值。

3 結(jié)合總體最小二乘的自適應(yīng)抗差濾波算法

在GNSS導(dǎo)航定位中，我們常常會(huì)遇到觀測(cè)噪聲和動(dòng)力學(xué)模型噪聲為有色噪聲且觀測(cè)系數(shù)矩陣存在偏差的情況，而現(xiàn)有的濾波算法并不能同時(shí)顧及有色噪聲和觀測(cè)系數(shù)矩陣偏差影響，這在一定程度上降低了導(dǎo)航定位結(jié)果的精度和可靠性，在實(shí)際的導(dǎo)航定位數(shù)據(jù)處理中不容忽視。這里，我們將把總體最小二乘和自適應(yīng)抗差濾波相結(jié)合，以便能同時(shí)減弱或消除有色噪聲和觀測(cè)系數(shù)矩陣偏差的影響。

式中各符號(hào)意義與（1）和（2）式相同。

在組成了新的系數(shù)矩陣和觀測(cè)方程后，我們可以用新的函數(shù)模型按自適應(yīng)抗差濾波算法計(jì)算狀態(tài)參數(shù)估值，這樣就能夠同時(shí)控制有色噪聲和觀測(cè)系數(shù)矩陣偏差的影響了。

4 計(jì)算與比較

本文對(duì)一組機(jī)載動(dòng)態(tài)GPS觀測(cè)數(shù)據(jù)進(jìn)行了處理，該組數(shù)據(jù)于1996年9月20日由兩臺(tái)Trimble4000SSE型接收機(jī)采集，其中一臺(tái)固定在參考站，一臺(tái)安置于飛機(jī)上。飛機(jī)的初始位置距離固定參考站約1公里，經(jīng)過10分鐘的初始化后起飛，飛行時(shí)間持續(xù)約90分鐘。

我們使用高精度的載波相位雙差解作為參考值，以便對(duì)濾波結(jié)果進(jìn)行比較和分析，使用C/A碼雙差觀測(cè)值進(jìn)行濾波結(jié)算，其結(jié)果與參考值進(jìn)行比較。濾波模型采用常速度模型，位置、速度和C/A碼的初始方差分別為 0.2 m2,、10-4m2·s-2和 1m2，速度的譜密度取為 5×10-3m2·s-3。狀態(tài)模型方差—協(xié)方差與文獻(xiàn)[31]中相同。

采用以下四種方案進(jìn)行比較分析：

方案1：標(biāo)準(zhǔn)Kalman濾波；

方案2：結(jié)合總體最小二乘的標(biāo)準(zhǔn)Kalman濾波；

方案3：自適應(yīng)抗差濾波；

方案4：結(jié)合總體最小二乘的自適應(yīng)抗差濾波。

計(jì)算結(jié)果如圖1～4和表1，圖中只給出了X軸方向的結(jié)果，Y軸和Z軸方向的結(jié)果與X軸方向類似。

圖1 標(biāo)準(zhǔn)Kalman濾波Fig.1 Standard Kalman filtering

圖2 結(jié)合總體最小二乘的標(biāo)準(zhǔn)Kalman濾波Fig.2 Standard Kalman filtering combined with Total Least Squares

圖3 自適應(yīng)抗差濾波Fig.3 Adaptively robust filtering

圖4 結(jié)合總體最小二乘的自適應(yīng)抗差濾波Fig.4 Adaptively robust filtering combined with Total Least Squares

表1 四種方案均方根誤差（RMS）/mTab.1 Root Mean Squares Error of Four Plans

由以上圖表可得出以下結(jié)論：

（1）由圖1、2可知，在飛機(jī)的飛行過程中有兩處明顯的非平穩(wěn)擾動(dòng)狀態(tài)，一處在800歷元附近，另一處在3000～3500歷元之間，即飛機(jī)起飛和轉(zhuǎn)彎時(shí)。這兩處比較大的擾動(dòng)對(duì)標(biāo)準(zhǔn)Kalman濾波和結(jié)合總體最小二乘的Kalman濾波都有很顯著的影響，這也說明結(jié)合總體最小二乘的Kalman濾波并不具有抵制動(dòng)態(tài)擾動(dòng)異常的能力。

（2）比較表1中方案1和方案2的結(jié)果可知，經(jīng)過總體最小二乘方法改正觀測(cè)系數(shù)矩陣偏差后，標(biāo)準(zhǔn)Kalman濾波的精度有一定的提高。這就說明觀測(cè)系數(shù)矩陣確實(shí)是存在偏差的，并且使用總體最小二乘方法能夠較好地改正這種偏差。

（3）將方案3和方案2的結(jié)果相比較可以發(fā)現(xiàn)，自適應(yīng)抗差濾波的結(jié)果比結(jié)合總體最小二乘的標(biāo)準(zhǔn)Kalman濾波的結(jié)果要好，而且整個(gè)濾波圖形比較平穩(wěn)，這是因?yàn)樽赃m應(yīng)抗差濾波不僅通過觀測(cè)等價(jià)權(quán)控制了觀測(cè)異常的影響，而且通過自適應(yīng)因子調(diào)整了動(dòng)力學(xué)模型信息對(duì)狀態(tài)參數(shù)估值的貢獻(xiàn)。同時(shí)也說明了，本算例中有色噪聲和異常對(duì)濾波精度的影響要遠(yuǎn)遠(yuǎn)大于觀測(cè)系數(shù)矩陣偏差。

（4）表1中方案4的結(jié)果要稍好于方案3，這說明在經(jīng)過總體最小二乘改正觀測(cè)系數(shù)矩陣偏差后，自適應(yīng)抗差濾波結(jié)果的精度仍然有一定的提高。這說明結(jié)合總體最小二乘的自適應(yīng)抗差濾波不僅能改正觀測(cè)系數(shù)矩陣偏差，而且能夠控制有色噪聲及異常的影響，新方法的可行性和有效性得到了驗(yàn)證。

5 結(jié)束語

在動(dòng)態(tài)導(dǎo)航定位中，常常會(huì)遇到觀測(cè)噪聲和動(dòng)力學(xué)模型噪聲都為有色噪聲，且觀測(cè)方程系數(shù)矩陣存在偏差的情況，數(shù)據(jù)處理時(shí)無論是忽略有色噪聲還是觀測(cè)系數(shù)矩陣偏差的影響，都將會(huì)降低導(dǎo)航定位的精度和可靠性。因此，實(shí)際的數(shù)據(jù)處理中有色噪聲和觀測(cè)系數(shù)矩陣偏差都不容忽視。本文提出的結(jié)合總體最小二乘的自適應(yīng)抗差濾波算法，不僅使用總體最小二乘方法改正了觀測(cè)系數(shù)矩陣偏差，而且通過自適應(yīng)抗差濾波控制有色噪聲和異常的影響，能夠進(jìn)一步提高動(dòng)態(tài)導(dǎo)航定位的精度和可靠性。

［1］張金槐，蔡洪.飛行器試驗(yàn)統(tǒng)計(jì)學(xué)[M].國防科技大學(xué)出版社,1995

［2］趙長勝.有色噪聲濾波理論與算法[M].測(cè)繪出版社,2011

［3］趙長勝，陶本藻.有色噪聲作用下的卡爾曼濾波[J].武漢大學(xué)學(xué)報(bào):信息科學(xué)版,2008,33(2):180-182

［4］楊元喜,崔先強(qiáng).動(dòng)態(tài)定位有色噪聲影響函數(shù)——以一階AR模型為例[J].測(cè)繪學(xué)報(bào),2003,32(1):6-10.

［5］崔先強(qiáng),楊元喜,高為廣.多種有色噪聲自適應(yīng)濾波算法的比較[J].武漢大學(xué)學(xué)報(bào):信息科學(xué)版,2006,31(8):731-735.

［6］崔先強(qiáng)，楊元喜.分類因子自適應(yīng)抗差濾波[J].自然科學(xué)進(jìn)展,2006,16(4):490-494.

［7］崔先強(qiáng),楊元喜,張曉東.導(dǎo)航解算中的有色噪聲及其協(xié)方差矩陣自適應(yīng)擬合[J].空間科學(xué)學(xué)報(bào),2012,32(4):592-597.

［8］楊元喜,何海波，徐天河.論動(dòng)態(tài)自適應(yīng)濾波[J].測(cè)繪學(xué)報(bào),2001,30(4):293-298.

［9］楊元喜,徐天河.基于移動(dòng)開窗法協(xié)方差估計(jì)和方差分量估計(jì)的自適應(yīng)濾波[J].武漢大學(xué)學(xué)報(bào):信息科學(xué)版，2003,28(6):714-718.

［10］楊元喜.動(dòng)態(tài)定位自適應(yīng)濾波解的性質(zhì)[J].測(cè)繪學(xué)報(bào),2003,32(3):189-192.

［11］Adcock R.Note on the Method of Least Squares[J].Analyst,1877(4):183-184.

［12］Pearson K.On Lines and Planes of Closest Fit to Points in Space[J].Philos Mag,1901(2):559-572.

［13］Koopmans T C.Linear Regression Analysis of Economic Time Series[M].De Erven F,Bohn N V,1937.

［14］Madansky A.The Fitting of Straight Lines when both Variables are Subject to Error[J].J Amer Statist Assoc,1959,54:173-205.

［15］York D.Least Squares Fitting of a Straight Line[J].Can J Phys,1966,44:1079-1086.

［16］Golub G H,Van Loan C F.An Analysis of the Total Least Squares problem[J].SIAM J Number Anal,1980,17:883-893.

［17］Gleser L J.Estimation in a Multivariate Errors in Variables Regression Model:Large Sample Results[J].Ann Statist,1989,9(1):24-44.

［18］Van Huffel S,Vandewalle J.Algebra Connection Between the Least Squares and Total Least Squares Problems[J].Numerical Math,1989,55:431-449.

［19］Stoicaa P,Vibergb M.Weighted LS and TLS Approaches Yiedl Asymptotically Equivalent Results[J].Signal Processing,1995,45:255-259.

［20］Schuermans M,Markovsky I,Wentzell P D,et al.On the Equivalance Between Total Least Squares and Maximum Likelihood PCA[J].Anal Chim Acta,2005,544:254-267.

［21］俞錦成.關(guān)于整體最小二乘問題的可解性[J].南京師大學(xué)報(bào):自然科學(xué)版,1996,19(1):13-16

［22］魯鐵定,陶本藻,周世健.基于整體最小二乘法的線性回歸建模和解法[J].武漢大學(xué)學(xué)報(bào):信息科學(xué)版,2008,33(5):504-507.

［23］孔建,姚宜斌,吳寒.整體最小二乘的迭代解法[J].武漢大學(xué)學(xué)報(bào):信息科學(xué)版,2010,35(6):711-714.

［24］丁克良.整體最小二乘理論及其在測(cè)量數(shù)據(jù)處理中的若干應(yīng)用研究[D].武漢:中國科學(xué)院測(cè)量與地球物理研究所,2006.

［25］楊元喜.自適應(yīng)動(dòng)態(tài)導(dǎo)航定位[M].測(cè)繪出版社,2006.

［26］Y Yang,H He,G Xu.Adaptively Robust Filtering for Kinematic Geodetic Positioning[J].Journal of Geodesy,2001,75:109-116

［27］楊元喜,任夏,許艷.自適應(yīng)抗差濾波理論及其應(yīng)用的主要進(jìn)展[J].導(dǎo)航定位學(xué)報(bào),2013,1(1):9-15.

［28］沈云中,陶本藻.實(shí)用測(cè)量數(shù)據(jù)處理方法[M].測(cè)繪出版社,2012.

［29］魯鐵定,周世健.總體最小二乘的迭代解法[J].武漢大學(xué)學(xué)報(bào):信息科學(xué)版,2010,35(11):1051-1054.

［30］魯鐵定.總體最小二乘平差理論及其在測(cè)繪數(shù)據(jù)處理中的應(yīng)用[D].武漢大學(xué),2010.

［31］Schwarz K P,Cannon M E,Wong RVC.A Comparison of GPS Kinematic Models for Determination of Position and Velocity along a Trajectory[J].Manuscr Geod,1989,14:345-353.