周艷青, 薛河儒, 姜新華, 王思宇, 王 靜

(內(nèi)蒙古農(nóng)業(yè)大學(xué) 計算機(jī)與信息工程學(xué)院,呼和浩特 010018)

錫林河流域?qū)儆诘湫托圆菰土饔?地處內(nèi)蒙古高原中部,流域總面積約為10 800 km2. 是錫林郭勒草原上的主要河流,因此對它展開研究具有重要的意義[1].IPCC第四次評估報告指出,近50年來全球地表增溫速度明顯加快,若氣候變化超出生態(tài)系統(tǒng)的彈性閾值,將嚴(yán)重破壞生態(tài)系統(tǒng)的結(jié)構(gòu)和功能,后果不堪設(shè)想[2].而氣象站的設(shè)定可以獲得有關(guān)氣候變化的氣象數(shù)據(jù).但是氣象數(shù)據(jù)是以時間序列而采集,通常以分鐘為單位,采集的密度濃且指標(biāo)多. 隨著時間的累積將獲得大量監(jiān)測的氣象數(shù)據(jù),增加了數(shù)據(jù)的復(fù)雜性. 而數(shù)據(jù)融合技術(shù)則是大數(shù)據(jù)處理中處理數(shù)據(jù)的一種方式. 通常的數(shù)據(jù)融合算法針對多個傳感器在同一時間不同空間的特征數(shù)據(jù)的融合. 而在實際應(yīng)用中,一個氣象站節(jié)點(diǎn)上不止一類傳感器,況且每種傳感器都用有自己的特性,例如溫度傳感器所采集的溫度變化緩慢,如果按照固定的采樣頻率,將會產(chǎn)生大量的重復(fù)冗余的數(shù)據(jù),并增加網(wǎng)絡(luò)的傳輸量,所以數(shù)據(jù)融合的研究成為人們關(guān)注的焦點(diǎn)[3]. 例如,劉衛(wèi)萍等人將數(shù)據(jù)融合技術(shù)應(yīng)用于環(huán)境測量中,對多傳感器監(jiān)測的環(huán)境數(shù)據(jù)進(jìn)行轉(zhuǎn)換、相關(guān)性分析、融合,降低數(shù)據(jù)的規(guī)模[4]. 張鵬鵬等對礦井安全監(jiān)測數(shù)據(jù)實施融合,提高數(shù)據(jù)采集的準(zhǔn)確性[5].

數(shù)據(jù)融合是充分利用不同時間與空間的多傳感器數(shù)據(jù)資源,采用計算機(jī)技術(shù)按時間序列獲得傳感器的觀測數(shù)據(jù),在一定準(zhǔn)則下進(jìn)行分析. 當(dāng)前,主要的數(shù)據(jù)融合方法有加權(quán)平均、卡爾曼濾波法、貝葉斯估計法、證據(jù)推理等. 其中貝葉斯估計法和證據(jù)推理是在靜態(tài)環(huán)境適用于高層的數(shù)據(jù)融合; 加權(quán)平均和卡爾曼濾波法適用于動態(tài)環(huán)境的低層數(shù)據(jù)的融合. 卡爾曼濾波是目前應(yīng)用最廣泛的濾波方法. 劉超云等人提出卡爾曼濾波算法對多個滑坡體位移監(jiān)測數(shù)據(jù)進(jìn)行濾波融合,從而對其穩(wěn)定狀態(tài)和變化趨勢做出預(yù)測[6]. 鄒波等人提出基于EKF的改進(jìn)非線性定姿估計方法,通過多傳感器互補(bǔ),利用儀器的進(jìn)行誤差修正,得到準(zhǔn)確的姿態(tài)角[7]. Ou等人通過多傳感器互補(bǔ),采用連續(xù)EKF進(jìn)行姿態(tài)估計,提高姿態(tài)角估計的動態(tài)和穩(wěn)態(tài)性能[8]. 龍慧提出基于Consensus濾波的分布式卡爾曼融合算法,引入一致濾波算法用于計算節(jié)點(diǎn)平均觀測數(shù)據(jù)和平均你協(xié)方差,從而獲得各節(jié)點(diǎn)的分布式狀態(tài)估計[9].Subhro等人提出一種單一時間尺度的分布式卡爾曼濾波算法,獲得不穩(wěn)定系統(tǒng)的有節(jié)MSE的無偏估計[10].Lvanjko等通過采用擴(kuò)展卡爾曼濾波和無跡卡爾曼濾波的融合實現(xiàn)移動機(jī)器人姿態(tài)的跟蹤[11]. 吳勇等提出一種收縮無跡卡爾曼濾波器,并應(yīng)用于SLAM問題中,通過設(shè)置收縮參數(shù)降低計算復(fù)雜度[12].

卡爾曼濾波是線性無偏最小方差估計,EKF是將非線性系統(tǒng)線性化,然后進(jìn)行卡爾曼濾波,但是線性化處理時需要用雅克比矩陣,其繁瑣的計算過程導(dǎo)致該方法實現(xiàn)相對困難. UKF是針對非線性系統(tǒng),但是計算量大. 與EKF相比較,它有更高的估計精度和更強(qiáng)的魯棒性及穩(wěn)定性,但當(dāng)周圍環(huán)境發(fā)生較大變化,其精度和穩(wěn)定性都會大大降低[13].

已提出的各種不同卡爾曼濾波算法,在融合的精度上和計算的復(fù)雜性有所提高,但是都是針對不同空間的多個傳感器的數(shù)據(jù)融合算法,屬于橫向融合. 而實際采集的氣象數(shù)據(jù)僅來源于相距較遠(yuǎn)的兩個氣象站,采集數(shù)據(jù)密度密,所以將上述算法應(yīng)用于縱向基于時間序列數(shù)據(jù)的融合,計算復(fù)雜性相對高. 因此,本文將卡爾曼濾波算法應(yīng)用于縱向基于時間序列數(shù)據(jù)的融合.因此本文以典型草原流域錫林河流域為研究對象,引入分布圖法,利用傳統(tǒng)的卡爾曼濾波算法對氣象數(shù)據(jù)中變化緩慢的空氣溫度指標(biāo)進(jìn)行同一空間不同時間序列的融合,以減少數(shù)據(jù)的傳輸量,提高數(shù)據(jù)的準(zhǔn)確性與科學(xué)性,同時方便以后的計算. 通常探測數(shù)據(jù)采用平均值作為瞬時值,如果測量的數(shù)據(jù)出現(xiàn)不完整性或者存在異常,將會導(dǎo)致最終結(jié)果的不準(zhǔn)確. 所以處理異常或者缺失數(shù)據(jù)進(jìn)行融合將是本研究的重點(diǎn).

1 融合算法

1.1 卡爾曼濾波算法

卡爾曼濾波算法采用遞歸方法來解決線性濾波問題,主要用于動態(tài)環(huán)境中冗余信息的融合,根據(jù)上一狀態(tài)的估計值和當(dāng)前狀態(tài)的觀測值推出當(dāng)前狀態(tài)的估計的濾波方法[14,15]. 首先引入一個離散控制過程的系統(tǒng)和系統(tǒng)的測量值,用公式(1)和(2)描述：

其中,Xk是k時刻對系統(tǒng)狀態(tài)變量,Uk是k時刻對系統(tǒng)的控制量,A和B是系統(tǒng)的參數(shù),Zk是k時刻的測量值,H為測量系統(tǒng)的參數(shù),Wk和Vk是過程和測量噪聲,假設(shè)其均為高斯白噪聲,噪聲協(xié)方差分別用q,r表示.

由公式(1)和(2),可得出卡爾曼濾波器的時間更新方程,包含時間遞推狀態(tài)變量計算和向前推算誤差協(xié)方差的計算：

基于現(xiàn)在時刻的觀測值和前一時刻的估計值,可得卡爾曼濾波器的狀態(tài)更新方程,包含觀測量的更新估計、卡爾曼增益計算和更新誤差協(xié)方差：

卡爾曼濾波算法通過不斷的預(yù)測和校正獲得最終的準(zhǔn)確的測量值.

1.2 卡爾曼濾波的參數(shù)設(shè)置

在實際濾波時,需要根據(jù)測量數(shù)據(jù)的實際情況設(shè)置初值X0,初始誤差協(xié)方差P0(P0≠0),過程噪聲協(xié)方差Q,和測量噪聲協(xié)方差R,使得濾波收斂速度和效果最佳. 由于卡爾曼濾波算法是一個最優(yōu)化自回歸的算法,所以X0可以取0或者測量的初始值. 通常P0越小表示初始的估計較好. 如果測量的環(huán)境相對穩(wěn)定,Q可設(shè)置為一個確定的值,Q的取值越接近0,融合的曲線越光滑,但不宜特別小.R與測量儀器的精度有關(guān),R取值與Q相類似,R越小,濾波效果好,收斂快[16].

1.3 改進(jìn)的卡爾曼濾波算法

1.3.1 分布圖法

對于一個相對穩(wěn)定的環(huán)境,傳統(tǒng)的卡爾曼濾波算法可以獲得較好的融合結(jié)果. 當(dāng)是當(dāng)數(shù)據(jù)出現(xiàn)異常或者缺失,傳統(tǒng)的卡爾曼濾波算法的融合結(jié)果將會出現(xiàn)波動. 針對該問題,對傳統(tǒng)的卡爾曼濾波算法實施改進(jìn),引入了分布圖法[17],對測量數(shù)據(jù)進(jìn)行處理,在利用卡爾曼濾波算法融合,得到可靠的融合結(jié)果. 分布圖法通過計算判斷區(qū)間[ρ1,ρ2]來排除50%的離異值干擾. 而且中位值和四分位離散度的選擇與極值點(diǎn)的大小無關(guān),僅取決于數(shù)據(jù)的分布位置,有效區(qū)間的獲得與疏失數(shù)據(jù)的關(guān)系不大,而且利用分布圖法獲得的數(shù)據(jù)不受數(shù)據(jù)分布的限制[18]. 所以分布圖法消除疏失數(shù)據(jù)具有運(yùn)算量小、魯棒性強(qiáng)、實時性好的優(yōu)點(diǎn)[19]. 為了保證數(shù)據(jù)的維數(shù)不變,將采用估計值來代替疏失數(shù)據(jù),這樣可以減少維數(shù)的判斷,提高計算的效率.

首先引入分布圖法的參數(shù),將測量的氣象溫度數(shù)據(jù)按從小到大的順序排列,設(shè)為T1,T2,…,TN,則中位數(shù)Tm按公式(8)的定義,中位數(shù)將數(shù)據(jù)序列分為上四分位FU和下四分位FL,四分位離散度dF為FU和FL的差值. 如果氣象溫度數(shù)據(jù)與中位數(shù)的距離大于βdF,則稱該數(shù)據(jù)為無效數(shù)據(jù),β為常數(shù).

假設(shè)有效數(shù)據(jù)的判斷區(qū)間為[ρ1,ρ2],則不包含在這個區(qū)間內(nèi)的數(shù)據(jù)認(rèn)為異常數(shù)據(jù). 其中,

當(dāng)判斷出異常數(shù)據(jù),為了使得數(shù)據(jù)維數(shù)保持一致,則通過對中間溫度數(shù)據(jù)求平均值,來代替當(dāng)前時刻的估計值,如表1.

表1 中間數(shù)據(jù)的平均值計算

1.3.2 算法步驟

改進(jìn)的卡爾曼濾波算法的具體步驟如下：

1)讀入融合的測量數(shù)據(jù),并對數(shù)據(jù)進(jìn)行排序,利用分布圖法確定異?；蚴枋?shù)據(jù),并用平均值來代替測量數(shù)據(jù)中的異常數(shù)據(jù);

2)初始化卡爾曼濾波算法的參數(shù),包括參數(shù)A,B,H,以及P0和X0的初始值,并計算過程和測量的噪聲的協(xié)方差q,r;

3)利用卡爾曼濾波算法中的公式,通過循環(huán)迭代計算,對測量數(shù)據(jù)實施融合.

2 實驗及結(jié)果分析

選取錫林河流域內(nèi)石門景區(qū)氣象站2016年1月1日每隔5分鐘所采集的空氣溫度數(shù)據(jù),所以1小時之內(nèi)可以采集到12個空氣溫室數(shù)據(jù),24小時,共288個采用數(shù)據(jù),部分?jǐn)?shù)據(jù)如表2所示,利用改進(jìn)的卡爾曼濾波算法按照小時對數(shù)據(jù)實施融合,將融合為12個數(shù)據(jù).

由于空氣溫度傳感器采集過程相對穩(wěn)定,且采集的數(shù)據(jù)與溫度是直接對應(yīng)的,所以采用一維線性的離散系統(tǒng),將變量A和H均設(shè)置為1. 噪聲的來源主要源于環(huán)境和采集儀器,所以Wk和Vk且均為零均值的獨(dú)立的高斯白噪聲,方差分別為q=0.04和r=0.2.

1)原始數(shù)據(jù)的融合

表2中的空氣溫度數(shù)據(jù)如圖1所示. 可以明顯的看出一天溫度變化趨勢,中午兩點(diǎn)空氣溫度已經(jīng)達(dá)到了極值,將其使用平均值、卡爾曼濾波和改進(jìn)的卡爾曼濾波融合后的結(jié)果如圖1(b)所示.

從圖1結(jié)果對比中可以看出,三種方法的融合趨勢大致相同. 在上午6點(diǎn),卡爾曼濾波算法的融合值存在微小的突變; 凌晨2點(diǎn)到3點(diǎn),平均值算法的融合結(jié)果也存在突變,所以與平均值算法和卡爾曼濾波算法相比,改進(jìn)的卡爾曼濾波算法的融合曲線更光滑,符合空氣溫度的緩慢變化規(guī)律

表 2 空氣溫度采樣數(shù)據(jù) (單位：°C)

圖1 原始數(shù)據(jù)及融合結(jié)果

2)添加擾動數(shù)據(jù)和突變數(shù)據(jù)

為了驗證改進(jìn)的算法的性能,在原始的空氣數(shù)據(jù)中設(shè)置了四個擾動數(shù)據(jù)和兩個突變數(shù)據(jù),其中突變數(shù)據(jù)是將原始數(shù)據(jù)設(shè)置為最大或者是0,0點(diǎn)和9點(diǎn)設(shè)置減小的擾動數(shù)據(jù),13點(diǎn)和22點(diǎn)設(shè)置增大的擾動數(shù)據(jù).設(shè)置的位置如表3所示.

表3 添加的擾動和突變的數(shù)據(jù)(℃)

由于設(shè)置了6個疏失數(shù)據(jù),利用分布圖法判斷該數(shù)據(jù),并賦予估計值,如表4所示.

從表中可知,估計值與原始數(shù)據(jù)值相差較小,誤差介于-0.4～0.28.

表4 疏失數(shù)據(jù)的估計值(℃)

對擾動數(shù)據(jù)和突變數(shù)據(jù)分別用平均值、卡爾曼濾波和改進(jìn)的卡爾波濾波算法融合,其融合結(jié)果顯示在圖2中,相應(yīng)的改變的值顯示在表5和表6中.

依據(jù)圖2(a)和表5可知,在0點(diǎn)設(shè)置的擾動數(shù)據(jù),平均值法融合結(jié)果顯然從-12.74℃下降到-13.76℃,不能避免擾動數(shù)據(jù)的干擾. 卡爾曼濾波算法的融合結(jié)果僅在-0.33到0.02的小范圍內(nèi)波動. 而改進(jìn)的算法的融合結(jié)果的波動范圍更小,為-0.1到0.09. 同時,與卡爾曼濾波算法,改進(jìn)方法的融合結(jié)果與設(shè)置擾動數(shù)據(jù)的增減性一致. 因此,改進(jìn)的算法不受擾動數(shù)據(jù)的影響.

圖2 設(shè)置異常數(shù)據(jù)的融合結(jié)果

表5 三種方法對擾動數(shù)據(jù)的融合結(jié)果對比(℃)

表6 三種方法對突變數(shù)據(jù)的融合結(jié)果對比(℃)

從圖2(b)和表6可知,平均值法對突變數(shù)據(jù)的融合結(jié)果影響大于其他兩種方法的融合結(jié)果. 對于改進(jìn)的融合算法在16點(diǎn)時,添加突變數(shù)據(jù)和原始數(shù)據(jù)的融合結(jié)果是一樣的,是由于改進(jìn)的融合算法將突變數(shù)據(jù)剔除,用表1中計算的平均值代替突變數(shù)據(jù). 因此改進(jìn)的方法能避免突變數(shù)據(jù). 綜上所述,利用分布圖法判斷疏失數(shù)據(jù),并用估計值來代替疏失數(shù)據(jù),在利用卡爾曼濾波算法進(jìn)行融合,能夠獲得光滑的融合結(jié)果,具有較強(qiáng)的穩(wěn)定性、健壯性.

對于時間跨度比較的大的溫度數(shù)據(jù),可以利用本文的算法對數(shù)據(jù)進(jìn)行融合,從而了解溫度的變化趨勢.選擇2016年1月1日到10日的數(shù)據(jù),每隔五分鐘采集,則共有2880條記錄. 按小時融合后的數(shù)據(jù)共240條,顯示在圖3中; 按天進(jìn)行融合的數(shù)據(jù)共10條,如圖4所示. 從中可以清楚的看到,溫度在1月6日和1月7日出現(xiàn)極值點(diǎn),且1月7日的溫度最低.

圖3 按小時的融合結(jié)果

圖4 按天的融合結(jié)果

3 結(jié)束語

因為氣象數(shù)據(jù)中的空氣溫度變化緩慢,并且受到噪聲干擾小,提出了一種線性模型的改進(jìn)的卡爾曼濾波融合算法. 它根據(jù)數(shù)據(jù)緩慢變化趨勢,利用上一時刻的估計值和當(dāng)前時刻的觀測值推出當(dāng)前狀態(tài)進(jìn)行融合.通過仿真實驗,改進(jìn)的算法能夠避免擾動數(shù)據(jù)和突變數(shù)據(jù)的干擾,為同一空間的時間序列數(shù)據(jù)融合提供新的方法,為下一步的氣象數(shù)據(jù)的研究奠定基礎(chǔ).

1 張雪峰,牛建明,張慶,等. 內(nèi)蒙古錫林河流域草地生態(tài)系統(tǒng)土壤保持功能及其空間分布. 草業(yè)學(xué)報,2015,24(1)：12-20. [doi：10.11686/cyxb20150103]

2 宋小園,朱仲元,張圣微,等. 錫林河流域氣候變化特征診斷分析. 干旱區(qū)資源與環(huán)境,2016,30(4)：151-158.

3 李翀,王沁,李磊,等. 一種基于時間序列的節(jié)點(diǎn)級數(shù)據(jù)融合方法. 計算機(jī)科學(xué),2008,35(11A)：307-311.

4 劉衛(wèi)萍,王寧,周曉磊,等. 數(shù)據(jù)融合技術(shù)在環(huán)境監(jiān)測領(lǐng)域的應(yīng)用. 計算機(jī)系統(tǒng)應(yīng)用,2016,25(6)：88-93.

5 張鵬鵬,俞阿龍,孫詩裕,等. 多傳感器數(shù)據(jù)融合在礦井安全監(jiān)測中的應(yīng)用. 工礦自動化,2015,41(12)：5-8.

6 劉超云,尹小波,張彬. 基于Kalman濾波數(shù)據(jù)融合技術(shù)的滑坡變形分析與預(yù)測. 中國地質(zhì)災(zāi)害與防治學(xué)報,2015,26(4)：30-35,42.

7 鄒波,張華,姜軍. 多傳感信息融合的改進(jìn)擴(kuò)展卡爾曼濾波定姿. 計算機(jī)應(yīng)用研究,2014,31(4)：1035-1038,1042.

8 Ou Y,Xia YQ,Fu MY. A modified method of nonlinear attitude estimation based on EKF. Proceedings of the 12th International Conference on Control Automation Robotics &Vision (ICARCV). Guangzhou,China. 2012. 901-906.

9 龍慧. 基于Consensus濾波的分布式卡爾曼信息融合方法.物聯(lián)網(wǎng)技術(shù),2011,(3)：61-64.

10 Das S,Moura JMF. Distributed Kalman filtering with dynamic observations consensus. IEEE Transactions on Signal Processing,2015,63(17)：4458-4473. [doi：10.1109/TSP.2015.2424205]

11 Ivanjko E,Vasak M,Petrovic I. Kalman filter theory based mobile robot pose tracking using occupancy grid maps.Proceedings of International Conference on Control and Automation. Budapest,Hungary. 2005,2. 869-874.

12 吳勇,關(guān)勝曉. 基于無跡卡爾曼濾波器的改進(jìn)SLAM問題求 解 方 法. 計 算 機(jī) 系 統(tǒng) 應(yīng) 用,2017,26(3)：30-36. [doi：10.15888/j.cnki.csa.005674]

13 周衛(wèi)東,喬相偉,吉宇人,等. 基于新息和殘差的自適應(yīng)UKF 算法. 宇航學(xué)報,2010,31(7)：1798-1804.

14 蔡小慶,魯小利,張偉娟,等. 基于卡爾曼濾波數(shù)據(jù)融合的溫室監(jiān)控系統(tǒng). 電子測試,2016,(6)：59-60.

15 Chen YK,Si XC,Li ZG. Research on Kalman-filter based multisensor data fusion. Journal of Systems Engineering and Electronics,2007,18(3)：497-502. [doi：10.1016/S1004-4132(07)60119-4]

16 Wang J,Xue HR,Jiang XH. Application of Kalman filtering algorithm in greenhouse environment monitoring.Proceedings of the 2013 2nd International Symposium on Instrumentation and Measurement,Sensor Network and Automation (IMSNA). Toronto,ON,Canada. 2013.539-544.

17 范滿紅,馬勝前,陳彥,等. 基于多傳感器數(shù)據(jù)融合的溫濕度監(jiān)測系統(tǒng). 壓電與聲光,2012,34(3)：459-462,465. [doi：10.11977/j.issn.1004-2474.2012.03.037]

18 夏卓君. 分布圖法在疏失誤差處理中的應(yīng)用. 實用測試技術(shù),2002,28(2)：33-34,8.

19 滕召勝. 基于多傳感器數(shù)據(jù)融合的熱處理爐溫度測量方法. 計量學(xué)報,2000,21(2)：148-152.