周艷青, 薛河儒, 姜新華, 王思宇, 王 靜
(內(nèi)蒙古農(nóng)業(yè)大學(xué) 計(jì)算機(jī)與信息工程學(xué)院,呼和浩特 010018)
錫林河流域?qū)儆诘湫托圆菰土饔?地處內(nèi)蒙古高原中部,流域總面積約為10 800 km2. 是錫林郭勒草原上的主要河流,因此對(duì)它展開(kāi)研究具有重要的意義[1].IPCC第四次評(píng)估報(bào)告指出,近50年來(lái)全球地表增溫速度明顯加快,若氣候變化超出生態(tài)系統(tǒng)的彈性閾值,將嚴(yán)重破壞生態(tài)系統(tǒng)的結(jié)構(gòu)和功能,后果不堪設(shè)想[2].而氣象站的設(shè)定可以獲得有關(guān)氣候變化的氣象數(shù)據(jù).但是氣象數(shù)據(jù)是以時(shí)間序列而采集,通常以分鐘為單位,采集的密度濃且指標(biāo)多. 隨著時(shí)間的累積將獲得大量監(jiān)測(cè)的氣象數(shù)據(jù),增加了數(shù)據(jù)的復(fù)雜性. 而數(shù)據(jù)融合技術(shù)則是大數(shù)據(jù)處理中處理數(shù)據(jù)的一種方式. 通常的數(shù)據(jù)融合算法針對(duì)多個(gè)傳感器在同一時(shí)間不同空間的特征數(shù)據(jù)的融合. 而在實(shí)際應(yīng)用中,一個(gè)氣象站節(jié)點(diǎn)上不止一類(lèi)傳感器,況且每種傳感器都用有自己的特性,例如溫度傳感器所采集的溫度變化緩慢,如果按照固定的采樣頻率,將會(huì)產(chǎn)生大量的重復(fù)冗余的數(shù)據(jù),并增加網(wǎng)絡(luò)的傳輸量,所以數(shù)據(jù)融合的研究成為人們關(guān)注的焦點(diǎn)[3]. 例如,劉衛(wèi)萍等人將數(shù)據(jù)融合技術(shù)應(yīng)用于環(huán)境測(cè)量中,對(duì)多傳感器監(jiān)測(cè)的環(huán)境數(shù)據(jù)進(jìn)行轉(zhuǎn)換、相關(guān)性分析、融合,降低數(shù)據(jù)的規(guī)模[4]. 張鵬鵬等對(duì)礦井安全監(jiān)測(cè)數(shù)據(jù)實(shí)施融合,提高數(shù)據(jù)采集的準(zhǔn)確性[5].
數(shù)據(jù)融合是充分利用不同時(shí)間與空間的多傳感器數(shù)據(jù)資源,采用計(jì)算機(jī)技術(shù)按時(shí)間序列獲得傳感器的觀測(cè)數(shù)據(jù),在一定準(zhǔn)則下進(jìn)行分析. 當(dāng)前,主要的數(shù)據(jù)融合方法有加權(quán)平均、卡爾曼濾波法、貝葉斯估計(jì)法、證據(jù)推理等. 其中貝葉斯估計(jì)法和證據(jù)推理是在靜態(tài)環(huán)境適用于高層的數(shù)據(jù)融合; 加權(quán)平均和卡爾曼濾波法適用于動(dòng)態(tài)環(huán)境的低層數(shù)據(jù)的融合. 卡爾曼濾波是目前應(yīng)用最廣泛的濾波方法. 劉超云等人提出卡爾曼濾波算法對(duì)多個(gè)滑坡體位移監(jiān)測(cè)數(shù)據(jù)進(jìn)行濾波融合,從而對(duì)其穩(wěn)定狀態(tài)和變化趨勢(shì)做出預(yù)測(cè)[6]. 鄒波等人提出基于EKF的改進(jìn)非線性定姿估計(jì)方法,通過(guò)多傳感器互補(bǔ),利用儀器的進(jìn)行誤差修正,得到準(zhǔn)確的姿態(tài)角[7]. Ou等人通過(guò)多傳感器互補(bǔ),采用連續(xù)EKF進(jìn)行姿態(tài)估計(jì),提高姿態(tài)角估計(jì)的動(dòng)態(tài)和穩(wěn)態(tài)性能[8]. 龍慧提出基于Consensus濾波的分布式卡爾曼融合算法,引入一致濾波算法用于計(jì)算節(jié)點(diǎn)平均觀測(cè)數(shù)據(jù)和平均你協(xié)方差,從而獲得各節(jié)點(diǎn)的分布式狀態(tài)估計(jì)[9].Subhro等人提出一種單一時(shí)間尺度的分布式卡爾曼濾波算法,獲得不穩(wěn)定系統(tǒng)的有節(jié)MSE的無(wú)偏估計(jì)[10].Lvanjko等通過(guò)采用擴(kuò)展卡爾曼濾波和無(wú)跡卡爾曼濾波的融合實(shí)現(xiàn)移動(dòng)機(jī)器人姿態(tài)的跟蹤[11]. 吳勇等提出一種收縮無(wú)跡卡爾曼濾波器,并應(yīng)用于SLAM問(wèn)題中,通過(guò)設(shè)置收縮參數(shù)降低計(jì)算復(fù)雜度[12].
卡爾曼濾波是線性無(wú)偏最小方差估計(jì),EKF是將非線性系統(tǒng)線性化,然后進(jìn)行卡爾曼濾波,但是線性化處理時(shí)需要用雅克比矩陣,其繁瑣的計(jì)算過(guò)程導(dǎo)致該方法實(shí)現(xiàn)相對(duì)困難. UKF是針對(duì)非線性系統(tǒng),但是計(jì)算量大. 與EKF相比較,它有更高的估計(jì)精度和更強(qiáng)的魯棒性及穩(wěn)定性,但當(dāng)周?chē)h(huán)境發(fā)生較大變化,其精度和穩(wěn)定性都會(huì)大大降低[13].
已提出的各種不同卡爾曼濾波算法,在融合的精度上和計(jì)算的復(fù)雜性有所提高,但是都是針對(duì)不同空間的多個(gè)傳感器的數(shù)據(jù)融合算法,屬于橫向融合. 而實(shí)際采集的氣象數(shù)據(jù)僅來(lái)源于相距較遠(yuǎn)的兩個(gè)氣象站,采集數(shù)據(jù)密度密,所以將上述算法應(yīng)用于縱向基于時(shí)間序列數(shù)據(jù)的融合,計(jì)算復(fù)雜性相對(duì)高. 因此,本文將卡爾曼濾波算法應(yīng)用于縱向基于時(shí)間序列數(shù)據(jù)的融合.因此本文以典型草原流域錫林河流域?yàn)檠芯繉?duì)象,引入分布圖法,利用傳統(tǒng)的卡爾曼濾波算法對(duì)氣象數(shù)據(jù)中變化緩慢的空氣溫度指標(biāo)進(jìn)行同一空間不同時(shí)間序列的融合,以減少數(shù)據(jù)的傳輸量,提高數(shù)據(jù)的準(zhǔn)確性與科學(xué)性,同時(shí)方便以后的計(jì)算. 通常探測(cè)數(shù)據(jù)采用平均值作為瞬時(shí)值,如果測(cè)量的數(shù)據(jù)出現(xiàn)不完整性或者存在異常,將會(huì)導(dǎo)致最終結(jié)果的不準(zhǔn)確. 所以處理異常或者缺失數(shù)據(jù)進(jìn)行融合將是本研究的重點(diǎn).
卡爾曼濾波算法采用遞歸方法來(lái)解決線性濾波問(wèn)題,主要用于動(dòng)態(tài)環(huán)境中冗余信息的融合,根據(jù)上一狀態(tài)的估計(jì)值和當(dāng)前狀態(tài)的觀測(cè)值推出當(dāng)前狀態(tài)的估計(jì)的濾波方法[14,15]. 首先引入一個(gè)離散控制過(guò)程的系統(tǒng)和系統(tǒng)的測(cè)量值,用公式(1)和(2)描述:

其中,Xk是k時(shí)刻對(duì)系統(tǒng)狀態(tài)變量,Uk是k時(shí)刻對(duì)系統(tǒng)的控制量,A和B是系統(tǒng)的參數(shù),Zk是k時(shí)刻的測(cè)量值,H為測(cè)量系統(tǒng)的參數(shù),Wk和Vk是過(guò)程和測(cè)量噪聲,假設(shè)其均為高斯白噪聲,噪聲協(xié)方差分別用q,r表示.
由公式(1)和(2),可得出卡爾曼濾波器的時(shí)間更新方程,包含時(shí)間遞推狀態(tài)變量計(jì)算和向前推算誤差協(xié)方差的計(jì)算:

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

卡爾曼濾波算法通過(guò)不斷的預(yù)測(cè)和校正獲得最終的準(zhǔn)確的測(cè)量值.
在實(shí)際濾波時(shí),需要根據(jù)測(cè)量數(shù)據(jù)的實(shí)際情況設(shè)置初值X0,初始誤差協(xié)方差P0(P0≠0),過(guò)程噪聲協(xié)方差Q,和測(cè)量噪聲協(xié)方差R,使得濾波收斂速度和效果最佳. 由于卡爾曼濾波算法是一個(gè)最優(yōu)化自回歸的算法,所以X0可以取0或者測(cè)量的初始值. 通常P0越小表示初始的估計(jì)較好. 如果測(cè)量的環(huán)境相對(duì)穩(wěn)定,Q可設(shè)置為一個(gè)確定的值,Q的取值越接近0,融合的曲線越光滑,但不宜特別小.R與測(cè)量?jī)x器的精度有關(guān),R取值與Q相類(lèi)似,R越小,濾波效果好,收斂快[16].
對(duì)于一個(gè)相對(duì)穩(wěn)定的環(huán)境,傳統(tǒng)的卡爾曼濾波算法可以獲得較好的融合結(jié)果. 當(dāng)是當(dāng)數(shù)據(jù)出現(xiàn)異常或者缺失,傳統(tǒng)的卡爾曼濾波算法的融合結(jié)果將會(huì)出現(xiàn)波動(dòng). 針對(duì)該問(wèn)題,對(duì)傳統(tǒng)的卡爾曼濾波算法實(shí)施改進(jìn),引入了分布圖法[17],對(duì)測(cè)量數(shù)據(jù)進(jìn)行處理,在利用卡爾曼濾波算法融合,得到可靠的融合結(jié)果. 分布圖法通過(guò)計(jì)算判斷區(qū)間[ρ1,ρ2]來(lái)排除50%的離異值干擾. 而且中位值和四分位離散度的選擇與極值點(diǎn)的大小無(wú)關(guān),僅取決于數(shù)據(jù)的分布位置,有效區(qū)間的獲得與疏失數(shù)據(jù)的關(guān)系不大,而且利用分布圖法獲得的數(shù)據(jù)不受數(shù)據(jù)分布的限制[18]. 所以分布圖法消除疏失數(shù)據(jù)具有運(yùn)算量小、魯棒性強(qiáng)、實(shí)時(shí)性好的優(yōu)點(diǎn)[19]. 為了保證數(shù)據(jù)的維數(shù)不變,將采用估計(jì)值來(lái)代替疏失數(shù)據(jù),這樣可以減少維數(shù)的判斷,提高計(jì)算的效率.
首先引入分布圖法的參數(shù),將測(cè)量的氣象溫度數(shù)據(jù)按從小到大的順序排列,設(shè)為T(mén)1,T2,…,TN,則中位數(shù)Tm按公式(8)的定義,中位數(shù)將數(shù)據(jù)序列分為上四分位FU和下四分位FL,四分位離散度dF為FU和FL的差值. 如果氣象溫度數(shù)據(jù)與中位數(shù)的距離大于βdF,則稱(chēng)該數(shù)據(jù)為無(wú)效數(shù)據(jù),β為常數(shù).

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

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

表1 中間數(shù)據(jù)的平均值計(jì)算
改進(jìn)的卡爾曼濾波算法的具體步驟如下:
1)讀入融合的測(cè)量數(shù)據(jù),并對(duì)數(shù)據(jù)進(jìn)行排序,利用分布圖法確定異常或疏失數(shù)據(jù),并用平均值來(lái)代替測(cè)量數(shù)據(jù)中的異常數(shù)據(jù);
2)初始化卡爾曼濾波算法的參數(shù),包括參數(shù)A,B,H,以及P0和X0的初始值,并計(jì)算過(guò)程和測(cè)量的噪聲的協(xié)方差q,r;
3)利用卡爾曼濾波算法中的公式,通過(guò)循環(huán)迭代計(jì)算,對(duì)測(cè)量數(shù)據(jù)實(shí)施融合.
選取錫林河流域內(nèi)石門(mén)景區(qū)氣象站2016年1月1日每隔5分鐘所采集的空氣溫度數(shù)據(jù),所以1小時(shí)之內(nèi)可以采集到12個(gè)空氣溫室數(shù)據(jù),24小時(shí),共288個(gè)采用數(shù)據(jù),部分?jǐn)?shù)據(jù)如表2所示,利用改進(jìn)的卡爾曼濾波算法按照小時(shí)對(duì)數(shù)據(jù)實(shí)施融合,將融合為12個(gè)數(shù)據(jù).
由于空氣溫度傳感器采集過(guò)程相對(duì)穩(wěn)定,且采集的數(shù)據(jù)與溫度是直接對(duì)應(yīng)的,所以采用一維線性的離散系統(tǒng),將變量A和H均設(shè)置為1. 噪聲的來(lái)源主要源于環(huán)境和采集儀器,所以Wk和Vk且均為零均值的獨(dú)立的高斯白噪聲,方差分別為q=0.04和r=0.2.
1)原始數(shù)據(jù)的融合
表2中的空氣溫度數(shù)據(jù)如圖1所示. 可以明顯的看出一天溫度變化趨勢(shì),中午兩點(diǎn)空氣溫度已經(jīng)達(dá)到了極值,將其使用平均值、卡爾曼濾波和改進(jìn)的卡爾曼濾波融合后的結(jié)果如圖1(b)所示.
從圖1結(jié)果對(duì)比中可以看出,三種方法的融合趨勢(shì)大致相同. 在上午6點(diǎn),卡爾曼濾波算法的融合值存在微小的突變; 凌晨2點(diǎn)到3點(diǎn),平均值算法的融合結(jié)果也存在突變,所以與平均值算法和卡爾曼濾波算法相比,改進(jìn)的卡爾曼濾波算法的融合曲線更光滑,符合空氣溫度的緩慢變化規(guī)律

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

圖1 原始數(shù)據(jù)及融合結(jié)果
2)添加擾動(dòng)數(shù)據(jù)和突變數(shù)據(jù)
為了驗(yàn)證改進(jìn)的算法的性能,在原始的空氣數(shù)據(jù)中設(shè)置了四個(gè)擾動(dòng)數(shù)據(jù)和兩個(gè)突變數(shù)據(jù),其中突變數(shù)據(jù)是將原始數(shù)據(jù)設(shè)置為最大或者是0,0點(diǎn)和9點(diǎn)設(shè)置減小的擾動(dòng)數(shù)據(jù),13點(diǎn)和22點(diǎn)設(shè)置增大的擾動(dòng)數(shù)據(jù).設(shè)置的位置如表3所示.

表3 添加的擾動(dòng)和突變的數(shù)據(jù)(℃)
由于設(shè)置了6個(gè)疏失數(shù)據(jù),利用分布圖法判斷該數(shù)據(jù),并賦予估計(jì)值,如表4所示.
從表中可知,估計(jì)值與原始數(shù)據(jù)值相差較小,誤差介于-0.4~0.28.

表4 疏失數(shù)據(jù)的估計(jì)值(℃)
對(duì)擾動(dòng)數(shù)據(jù)和突變數(shù)據(jù)分別用平均值、卡爾曼濾波和改進(jìn)的卡爾波濾波算法融合,其融合結(jié)果顯示在圖2中,相應(yīng)的改變的值顯示在表5和表6中.
依據(jù)圖2(a)和表5可知,在0點(diǎn)設(shè)置的擾動(dòng)數(shù)據(jù),平均值法融合結(jié)果顯然從-12.74℃下降到-13.76℃,不能避免擾動(dòng)數(shù)據(jù)的干擾. 卡爾曼濾波算法的融合結(jié)果僅在-0.33到0.02的小范圍內(nèi)波動(dòng). 而改進(jìn)的算法的融合結(jié)果的波動(dòng)范圍更小,為-0.1到0.09. 同時(shí),與卡爾曼濾波算法,改進(jìn)方法的融合結(jié)果與設(shè)置擾動(dòng)數(shù)據(jù)的增減性一致. 因此,改進(jìn)的算法不受擾動(dòng)數(shù)據(jù)的影響.

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

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

表6 三種方法對(duì)突變數(shù)據(jù)的融合結(jié)果對(duì)比(℃)
從圖2(b)和表6可知,平均值法對(duì)突變數(shù)據(jù)的融合結(jié)果影響大于其他兩種方法的融合結(jié)果. 對(duì)于改進(jìn)的融合算法在16點(diǎn)時(shí),添加突變數(shù)據(jù)和原始數(shù)據(jù)的融合結(jié)果是一樣的,是由于改進(jìn)的融合算法將突變數(shù)據(jù)剔除,用表1中計(jì)算的平均值代替突變數(shù)據(jù). 因此改進(jìn)的方法能避免突變數(shù)據(jù). 綜上所述,利用分布圖法判斷疏失數(shù)據(jù),并用估計(jì)值來(lái)代替疏失數(shù)據(jù),在利用卡爾曼濾波算法進(jìn)行融合,能夠獲得光滑的融合結(jié)果,具有較強(qiáng)的穩(wěn)定性、健壯性.
對(duì)于時(shí)間跨度比較的大的溫度數(shù)據(jù),可以利用本文的算法對(duì)數(shù)據(jù)進(jìn)行融合,從而了解溫度的變化趨勢(shì).選擇2016年1月1日到10日的數(shù)據(jù),每隔五分鐘采集,則共有2880條記錄. 按小時(shí)融合后的數(shù)據(jù)共240條,顯示在圖3中; 按天進(jìn)行融合的數(shù)據(jù)共10條,如圖4所示. 從中可以清楚的看到,溫度在1月6日和1月7日出現(xiàn)極值點(diǎn),且1月7日的溫度最低.

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

圖4 按天的融合結(jié)果
因?yàn)闅庀髷?shù)據(jù)中的空氣溫度變化緩慢,并且受到噪聲干擾小,提出了一種線性模型的改進(jìn)的卡爾曼濾波融合算法. 它根據(jù)數(shù)據(jù)緩慢變化趨勢(shì),利用上一時(shí)刻的估計(jì)值和當(dāng)前時(shí)刻的觀測(cè)值推出當(dāng)前狀態(tài)進(jìn)行融合.通過(guò)仿真實(shí)驗(yàn),改進(jìn)的算法能夠避免擾動(dòng)數(shù)據(jù)和突變數(shù)據(jù)的干擾,為同一空間的時(shí)間序列數(shù)據(jù)融合提供新的方法,為下一步的氣象數(shù)據(jù)的研究奠定基礎(chǔ).
1 張雪峰,牛建明,張慶,等. 內(nèi)蒙古錫林河流域草地生態(tài)系統(tǒng)土壤保持功能及其空間分布. 草業(yè)學(xué)報(bào),2015,24(1):12-20. [doi:10.11686/cyxb20150103]
2 宋小園,朱仲元,張圣微,等. 錫林河流域氣候變化特征診斷分析. 干旱區(qū)資源與環(huán)境,2016,30(4):151-158.
3 李翀,王沁,李磊,等. 一種基于時(shí)間序列的節(jié)點(diǎn)級(jí)數(shù)據(jù)融合方法. 計(jì)算機(jī)科學(xué),2008,35(11A):307-311.
4 劉衛(wèi)萍,王寧,周曉磊,等. 數(shù)據(jù)融合技術(shù)在環(huán)境監(jiān)測(cè)領(lǐng)域的應(yīng)用. 計(jì)算機(jī)系統(tǒng)應(yīng)用,2016,25(6):88-93.
5 張鵬鵬,俞阿龍,孫詩(shī)裕,等. 多傳感器數(shù)據(jù)融合在礦井安全監(jiān)測(cè)中的應(yīng)用. 工礦自動(dòng)化,2015,41(12):5-8.
6 劉超云,尹小波,張彬. 基于Kalman濾波數(shù)據(jù)融合技術(shù)的滑坡變形分析與預(yù)測(cè). 中國(guó)地質(zhì)災(zāi)害與防治學(xué)報(bào),2015,26(4):30-35,42.
7 鄒波,張華,姜軍. 多傳感信息融合的改進(jìn)擴(kuò)展卡爾曼濾波定姿. 計(jì)算機(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)勝曉. 基于無(wú)跡卡爾曼濾波器的改進(jìn)SLAM問(wèn)題求 解 方 法. 計(jì) 算 機(jī) 系 統(tǒng) 應(yīng) 用,2017,26(3):30-36. [doi:10.15888/j.cnki.csa.005674]
13 周衛(wèi)東,喬相偉,吉宇人,等. 基于新息和殘差的自適應(yīng)UKF 算法. 宇航學(xué)報(bào),2010,31(7):1798-1804.
14 蔡小慶,魯小利,張偉娟,等. 基于卡爾曼濾波數(shù)據(jù)融合的溫室監(jiān)控系統(tǒng). 電子測(cè)試,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 范滿(mǎn)紅,馬勝前,陳彥,等. 基于多傳感器數(shù)據(jù)融合的溫濕度監(jiān)測(cè)系統(tǒng). 壓電與聲光,2012,34(3):459-462,465. [doi:10.11977/j.issn.1004-2474.2012.03.037]
18 夏卓君. 分布圖法在疏失誤差處理中的應(yīng)用. 實(shí)用測(cè)試技術(shù),2002,28(2):33-34,8.
19 滕召勝. 基于多傳感器數(shù)據(jù)融合的熱處理爐溫度測(cè)量方法. 計(jì)量學(xué)報(bào),2000,21(2):148-152.