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

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

0 引 言

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

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

本文針對氣象無人機的測風特點和Kalman 濾波在數(shù)據處理中的局限性，將強跟蹤Kalman 濾波和抗野值算法應用于無人機探測數(shù)據處理中，達到抑制濾波發(fā)散，提高濾波精度的目的。

1 濾波算法

1.1 Kalman 濾波算法

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

狀態(tài)方程

量測方程

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

時間更新方程:

狀態(tài)一步預測

一步驗前誤差方差陣

量測更新方程:

增益陣

濾波方程

驗后誤差方差陣

新息序列

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

1.2 改進的強跟蹤Kalman 濾波算法

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

其中

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

式中 各項加權系數(shù)均為1/k，但風場估計中應該強調新進信息的作用，因此，可以改變因子，使和式中的各項乘以不同的加權系數(shù)。應用指數(shù)加權方法，新建加權系數(shù)μi，并使其滿足

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

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

1.3 抗野值修正算法

帶有野值的數(shù)據樣本往往會使Kalman 濾波器對系統(tǒng)的狀態(tài)預報修正錯誤，使濾波結果發(fā)生偏移，以往常采用3δ 準則進行野值判別，但其判別尺度把握不靈活且修正值的準確度不高。為了提高抗野值能力，保證濾波精度，文獻［7］提出了一種利用新息序列加權的方法來消除野值影響的方法。當量測數(shù)據不包含野值時，濾波器能夠有效利用新息，同時對量測噪聲進行在線估計;當量測數(shù)據中存在野值時，能克服其不利影響，盡可能還原系統(tǒng)真實狀態(tài)?？挂爸邓惴ㄈ缦?/p>

其中

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

2 抗野值自適應濾波算法

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

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

3 測風實驗與結果分析

3.1 風場模型建立

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

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

3.2 模擬數(shù)據處理

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

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

圖2 仿真模型風速估計對比圖Fig 2 Comparison of wind speed estimation in simulation model

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

3.3 實測數(shù)據處理

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

圖3 風速估計對比圖Fig 3 Comparison of wind speed estimation

圖4 風速估計局部對比圖Fig 4 Local comparison of wind speed estimation

4 結 論

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

［1］ 張曉芳，嚴 衛(wèi).中高層大氣探測技術的研究進展［J］.氣象科學，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］ 胡 鋒，孫國基.Kalman 濾波的抗野值修正［J］.自動化學報，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］ 柯 晶，錢積新.基于邏輯轉換的改進強跟蹤卡爾曼濾波器［J］.電子學報，2002，30(6):925－927.

［7］ 高 寧，周躍慶，楊 曄，等.抗野值自適應卡爾曼濾波方法的研究［J］.中國慣性技術學報，2003，11(3):25－28.

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

［9］ 蔡 崧，產竹旺.基于自適應濾波的風場測量仿真試驗平臺［J］.計算機工程，2003，29(18):192－194.

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