李 陽(yáng)，孔 毅，趙現(xiàn)斌

(1.解放軍理工大學(xué) 氣象海洋學(xué)院，江蘇 南京211101;2.中國(guó)人民解放軍94994 部隊(duì)，江蘇 南京210017)

0 引 言

隨著氣象無(wú)人機(jī)技術(shù)的發(fā)展，利用無(wú)人機(jī)搭載皮托靜壓管進(jìn)行空中風(fēng)［1］探測(cè)日益引起國(guó)內(nèi)外關(guān)注。皮托靜壓管測(cè)風(fēng)是把無(wú)人機(jī)作為探測(cè)平臺(tái)［2］，搭載機(jī)載氣象傳感器和數(shù)據(jù)采集裝置進(jìn)行風(fēng)速測(cè)量。這種大區(qū)域、長(zhǎng)時(shí)間、連續(xù)風(fēng)場(chǎng)探測(cè)方法在航天器發(fā)射與返回、重要武器試驗(yàn)、戰(zhàn)場(chǎng)氣象測(cè)量、惡劣天氣監(jiān)測(cè)、龍卷風(fēng)近距環(huán)境探測(cè)監(jiān)視等應(yīng)用中具有獨(dú)特的作用和優(yōu)勢(shì)。

在高空風(fēng)場(chǎng)探測(cè)過(guò)程中，由于受傳感器工作狀態(tài)、環(huán)境雜波等因素的影響，量測(cè)數(shù)據(jù)中經(jīng)常存在一些異常值。因此，風(fēng)場(chǎng)數(shù)據(jù)處理方法是無(wú)人機(jī)測(cè)風(fēng)的一項(xiàng)重要內(nèi)容，對(duì)提高測(cè)風(fēng)精度具有重要意義。當(dāng)前常使用Kalman濾波算法對(duì)風(fēng)場(chǎng)數(shù)據(jù)進(jìn)行處理，但其過(guò)于依賴系統(tǒng)噪聲和量測(cè)噪聲的統(tǒng)計(jì)特性，且缺乏對(duì)測(cè)量數(shù)據(jù)過(guò)失誤差的抗擾性和對(duì)傳感器突發(fā)性故障的容錯(cuò)能力［3］，降低了濾波精度和抗野值能力，容易造成系統(tǒng)的不穩(wěn)定。

本文針對(duì)氣象無(wú)人機(jī)的測(cè)風(fēng)特點(diǎn)和Kalman 濾波在數(shù)據(jù)處理中的局限性，將強(qiáng)跟蹤Kalman 濾波和抗野值算法應(yīng)用于無(wú)人機(jī)探測(cè)數(shù)據(jù)處理中，達(dá)到抑制濾波發(fā)散，提高濾波精度的目的。

1 濾波算法

1.1 Kalman 濾波算法

Kalman 濾波是典型的最小方差(MMSE)估計(jì)方法，采用 遞歸技術(shù)，利用k－1時(shí)刻狀態(tài)值給出k的預(yù)測(cè)值，并保證該均方誤差最小。建立離散系統(tǒng)模型

狀態(tài)方程

量測(cè)方程

上述模型中，Xk為狀態(tài)向量，Zk為觀測(cè)向量，φk|k－1為k－1 時(shí)刻到k 時(shí)刻的狀態(tài)轉(zhuǎn)移矩陣，Hk為量測(cè)矩陣，Wk－1和Vk分別是均值為零，方差為Qk－1和Rk且相互獨(dú)立的過(guò)程噪聲和量測(cè)噪聲。由于Kalman 濾波是一個(gè)帶回饋的估計(jì)方法，可將其分為時(shí)間更新和量測(cè)更新兩個(gè)階段。

時(shí)間更新方程:

狀態(tài)一步預(yù)測(cè)

一步驗(yàn)前誤差方差陣

量測(cè)更新方程:

增益陣

濾波方程

驗(yàn)后誤差方差陣

新息序列

上述方程未提供相關(guān)噪聲的估計(jì)模型，且無(wú)人機(jī)測(cè)風(fēng)受噪聲影響較大，因此，需對(duì)系統(tǒng)噪聲和量測(cè)噪聲進(jìn)行分析。由于系統(tǒng)干擾存在穩(wěn)定性，可設(shè)系統(tǒng)噪聲方差Q 為定值;量測(cè)噪聲受周?chē)h(huán)境的影響較為明顯，需對(duì)量測(cè)噪聲方差R 在線估計(jì)。利用Sage-Husa 自適應(yīng)濾波［4］，得到量測(cè)噪聲方差

1.2 改進(jìn)的強(qiáng)跟蹤Kalman 濾波算法

強(qiáng)跟蹤Kalman 濾波［5］是在Kalman 濾波理論的基礎(chǔ)上提出的，為保證濾波器的可靠收斂，在一步驗(yàn)前誤差方差陣中引入可在線計(jì)算的時(shí)變漸消矩陣［6］。該方法在風(fēng)場(chǎng)突變時(shí)具有很強(qiáng)的跟蹤能力，同時(shí)對(duì)噪聲統(tǒng)計(jì)特性的敏感性也較低。因此，將式(4)修改為

其中

應(yīng)用統(tǒng)計(jì)方法，可用算術(shù)平均值將v0(k)近似表示

式中 各項(xiàng)加權(quán)系數(shù)均為1/k，但風(fēng)場(chǎng)估計(jì)中應(yīng)該強(qiáng)調(diào)新進(jìn)信息的作用，因此，可以改變因子，使和式中的各項(xiàng)乘以不同的加權(quán)系數(shù)。應(yīng)用指數(shù)加權(quán)方法，新建加權(quán)系數(shù)μi，并使其滿足

由等比數(shù)列求和公式可推導(dǎo)出

其中，a 為遺忘因子，式(15)中的和式各項(xiàng)乘以加權(quán)系數(shù)μk－i，代替原來(lái)的1/k，得到v0(k)近似估計(jì)算法

1.3 抗野值修正算法

帶有野值的數(shù)據(jù)樣本往往會(huì)使Kalman 濾波器對(duì)系統(tǒng)的狀態(tài)預(yù)報(bào)修正錯(cuò)誤，使濾波結(jié)果發(fā)生偏移，以往常采用3δ 準(zhǔn)則進(jìn)行野值判別，但其判別尺度把握不靈活且修正值的準(zhǔn)確度不高。為了提高抗野值能力，保證濾波精度，文獻(xiàn)［7］提出了一種利用新息序列加權(quán)的方法來(lái)消除野值影響的方法。當(dāng)量測(cè)數(shù)據(jù)不包含野值時(shí)，濾波器能夠有效利用新息，同時(shí)對(duì)量測(cè)噪聲進(jìn)行在線估計(jì);當(dāng)量測(cè)數(shù)據(jù)中存在野值時(shí)，能克服其不利影響，盡可能還原系統(tǒng)真實(shí)狀態(tài)。抗野值算法如下

其中

式(18)為壓縮影響函數(shù);式(19)為權(quán)矩陣;式(20)為選取的門(mén)限常數(shù)序列。此外，λk為矩陣的最大特征值，為以χ2(m)分布置信度為(1－α) ×100%的上分位點(diǎn)，通常情況下α 取0.05 或0.025，m 為量測(cè)數(shù)據(jù)Zk的維數(shù)。上述模型抗野值原理是通過(guò)引入壓縮影響函數(shù)，改變野值的影響權(quán)重，抑制野值帶來(lái)的異常新息，達(dá)到抗野值的目的。

2 抗野值自適應(yīng)濾波算法

自然風(fēng)大致可分解為定常風(fēng)、紊流、風(fēng)切變和突風(fēng)四個(gè)分量［8］。通常情況下，在風(fēng)場(chǎng)探測(cè)系統(tǒng)中，將風(fēng)速表示成定常風(fēng)和紊流之和，并在此基礎(chǔ)上疊加量測(cè)噪聲，一旦量測(cè)值中出現(xiàn)野值或?yàn)V波發(fā)散時(shí)，誤差協(xié)方差陣都將無(wú)界，此時(shí)實(shí)際的新息往往比理論預(yù)計(jì)的大很多倍。實(shí)際的新息用vk表示，理論預(yù)計(jì)誤差的新息用新息序列協(xié)方差陣E［vk來(lái)描述，。因此，野值是否存在或?yàn)V波是否發(fā)散的依據(jù)為

式中 κ 為判別因子，其大小往往取3 或4。當(dāng)其成立時(shí)，量測(cè)值Zk為正常量測(cè)值，直接采用引入Sage-Husa 的Kalman 濾波進(jìn)行最優(yōu)估計(jì);反之，說(shuō)明量測(cè)值Zk為野值或?yàn)V波發(fā)散。設(shè)置連續(xù)野值門(mén)限值，當(dāng)連續(xù)的野值數(shù)多于該門(mén)限值或?yàn)V波發(fā)散時(shí)采用改進(jìn)的強(qiáng)跟蹤Kalman 濾波算法;當(dāng)存在少數(shù)野值時(shí)采用抗野值修正算法。抗野值強(qiáng)跟蹤Kalman 濾波算法流程如圖1。

3 測(cè)風(fēng)實(shí)驗(yàn)與結(jié)果分析

3.1 風(fēng)場(chǎng)模型建立

飛行仿真中，定常風(fēng)用低頻風(fēng)模型［9］描述，紊流場(chǎng)由基于空間相關(guān)函數(shù)的方法［10］建立，風(fēng)切邊、突風(fēng)等在內(nèi)的大氣擾動(dòng)分量直接疊加到量測(cè)噪聲上，本文不予單獨(dú)考慮，建立風(fēng)場(chǎng)模型

式中 VD，VW，y 分別為定常風(fēng)速、紊流風(fēng)速和量測(cè)風(fēng)速;a為紊流場(chǎng)衰減系數(shù);vD，vW，vO為互不相關(guān)的零均值高斯白噪聲，方差分別為QD，QW，RO。a 和QW由采樣步長(zhǎng)h、紊流尺度L 和紊流強(qiáng)度σ 決定

3.2 模擬數(shù)據(jù)處理

在模擬風(fēng)場(chǎng)中，加入服從高斯分布的量測(cè)噪聲，設(shè)置濾波各參數(shù)初值如下

由于風(fēng)是一種大范圍分布的空氣運(yùn)動(dòng)，受周?chē)h(huán)境和測(cè)量?jī)x器的影響導(dǎo)致數(shù)據(jù)中摻雜野值，因此，仿真中在多處加入孤立野值和連續(xù)野值。分別使用Kalman 濾波算法和本文提出的抗野值強(qiáng)跟蹤Kalman 濾波算法進(jìn)行濾波處理，圖2 是不同濾波算法的風(fēng)速估計(jì)比較圖。

圖2 仿真模型風(fēng)速估計(jì)對(duì)比圖Fig 2 Comparison of wind speed estimation in simulation model

由圖2 可以看出:Kalman 濾波具有一定的風(fēng)場(chǎng)跟蹤能力，但并未取得理想的濾波效果，受成片野值影響，結(jié)果出現(xiàn)大幅度跳變，濾波出現(xiàn)發(fā)散;而抗野值強(qiáng)跟蹤Kalman 濾波在保證濾波精度的同時(shí)，能夠有效抑制孤立野值或連續(xù)的影響，且精度較高，具有良好的魯棒性。因此，抗野值強(qiáng)跟蹤Kalman 濾波算法無(wú)論在濾波的精度還是收斂性方面均遠(yuǎn)優(yōu)于Kalman 濾波，適合工程應(yīng)用。

3.3 實(shí)測(cè)數(shù)據(jù)處理

以某型氣象無(wú)人機(jī)2008 年4 月17 日飛行探測(cè)數(shù)據(jù)為例，將本文方法用于實(shí)測(cè)數(shù)據(jù)處理中，如圖3 所示。圖中可以看到無(wú)人機(jī)初期和末期處于爬升和下降階段，受氣流擾動(dòng)較為嚴(yán)重，風(fēng)速值出現(xiàn)一定幅度的波動(dòng)，甚至出現(xiàn)野值，兩種數(shù)據(jù)處理方法都具有一定的跟蹤能力，但Kalman 規(guī)避連續(xù)野值的效果不佳，且濾波精度不高，而抗野值強(qiáng)跟蹤Kalman 濾波成功處理了連續(xù)野值，表現(xiàn)出較好的抗野值能力。圖4 是圖3 的局部放大圖。由圖4(a)可知，抗野值強(qiáng)跟蹤Kalman 濾波算法較Kalman 濾波算法具有更好的風(fēng)場(chǎng)跟蹤能力;由圖4(b)可知，抗野值強(qiáng)跟蹤Kalman 濾波算法具有更優(yōu)秀的抗野值能力。

圖3 風(fēng)速估計(jì)對(duì)比圖Fig 3 Comparison of wind speed estimation

圖4 風(fēng)速估計(jì)局部對(duì)比圖Fig 4 Local comparison of wind speed estimation

4 結(jié) 論

濾波的精度和穩(wěn)定性是無(wú)人機(jī)測(cè)風(fēng)數(shù)據(jù)處理的關(guān)鍵，常用的Kalman 濾波算法的跟蹤能力和抗野值能力欠佳。因此，本文結(jié)合Sage-Husa 濾波算法、強(qiáng)跟蹤濾波算法、抗野值修正算法的特點(diǎn)，根據(jù)實(shí)際新息和理論預(yù)計(jì)的判別關(guān)系，自適應(yīng)選擇濾波算法。最后利用仿真和實(shí)測(cè)風(fēng)場(chǎng)數(shù)據(jù)證明了該方法不僅克服了孤立野值和成片野值對(duì)濾波帶來(lái)的不利影響，同時(shí)也保證了濾波精度，具有較強(qiáng)的魯棒性，可作為無(wú)人機(jī)測(cè)風(fēng)系統(tǒng)有效的數(shù)據(jù)處理方法。

［1］ 張曉芳，嚴(yán) 衛(wèi).中高層大氣探測(cè)技術(shù)的研究進(jìn)展［J］.氣象科學(xué)，2007，27(4):457－463.

［2］ Kang H，Yang Q，Butler C，et al.Optimization of sensor locations for measurement of flue gas flow in industrial ducts and stacks using neural networks［J］.IEEE Transactions on Instrumentation and Measurement，2000，49(2):228－233.

［3］ 胡 鋒，孫國(guó)基.Kalman 濾波的抗野值修正［J］.自動(dòng)化學(xué)報(bào)，1999，25(5):692－696.

［4］ Efe M，Bather J A，Atherton D P.An adaptive Kalman filter with sequential rescaling of process noise［C］∥American Control Conference，1999:3913－3917.

［5］ Ristic B，Arulampalam M S.Tracking a manoeuvring target using angle-only measurements algorithms and performance［J］.Signal Processing，2003，83:1223－1238.

［6］ 柯 晶，錢(qián)積新.基于邏輯轉(zhuǎn)換的改進(jìn)強(qiáng)跟蹤卡爾曼濾波器［J］.電子學(xué)報(bào)，2002，30(6):925－927.

［7］ 高 寧，周躍慶，楊 曄，等.抗野值自適應(yīng)卡爾曼濾波方法的研究［J］.中國(guó)慣性技術(shù)學(xué)報(bào)，2003，11(3):25－28.

［8］ 肖亞倫，金長(zhǎng)江.大氣擾動(dòng)中的飛行原理［M］.北京:國(guó)防工業(yè)出版社，1993.

［9］ 蔡 崧，產(chǎn)竹旺.基于自適應(yīng)濾波的風(fēng)場(chǎng)測(cè)量仿真試驗(yàn)平臺(tái)［J］.計(jì)算機(jī)工程，2003，29(18):192－194.

［10］陸宇平，胡亞海.基于空間相關(guān)函數(shù)的二維紊流場(chǎng)數(shù)值生成法［J］.南京航空航天大學(xué)學(xué)報(bào)，1999，31(2):139－145.