基于改進(jìn)卡爾曼濾波的CO2泄漏監(jiān)測數(shù)據(jù)融合算法

楊兆發(fā)， 胡業(yè)林

(安徽理工大學(xué)電氣與信息工程學(xué)院 安徽 淮南 232001)

0 引言

在應(yīng)對全球氣候變化的CCUS(Carbon Capture, Utilization and Sequestration)技術(shù)[1]得到開創(chuàng)性應(yīng)用嘗試的背景下，對CO2封存泄漏[2]的研究和監(jiān)測成為值得關(guān)注的熱點。關(guān)于CCUS封存環(huán)境地面CO2泄漏監(jiān)測，一方面，監(jiān)測環(huán)境處于有人類活動干擾的自然環(huán)境下，導(dǎo)致監(jiān)測到的數(shù)據(jù)存在一定的誤差[3]。另一方面，地質(zhì)封存的CO2產(chǎn)生泄漏的過程是緩慢的，一般在幾十年甚至上百年的時間尺度[4]。針對該情況下的監(jiān)測數(shù)據(jù)處理問題，一個新的解決思路是采用數(shù)據(jù)濾波融合算法來提高監(jiān)測的準(zhǔn)確度。

卡爾曼濾波算法是在動態(tài)數(shù)據(jù)濾波領(lǐng)域應(yīng)用最廣泛的算法[5-8]。段杰等[5]利用卡爾曼濾波算法對農(nóng)業(yè)大棚信息監(jiān)測系統(tǒng)采集終端傳感器采集到的大量數(shù)據(jù)進(jìn)行融合處理，提高了監(jiān)測精度，得到了更加平穩(wěn)的數(shù)據(jù)值。吳勇等[7]研究SLAM問題過程中提出一種收縮無跡卡爾曼濾波器，通過設(shè)置收縮參數(shù)降低了計算復(fù)雜度。周艷青等[8]研究錫林河流域空氣溫度數(shù)據(jù)的規(guī)律時，提出了一種引入分布圖法改進(jìn)的卡爾曼濾波的空氣溫度數(shù)據(jù)融合算法，獲得了較強的抗干擾性和穩(wěn)定性，提高了氣象數(shù)據(jù)處理的準(zhǔn)確性。

考慮到CCUS封存環(huán)境下的地面CO2泄漏監(jiān)測數(shù)據(jù)處理分析在超長時間尺度下減少干擾和降低冗余的需要，采用在監(jiān)測時間序列上分段預(yù)處理，將預(yù)處理算法融入傳統(tǒng)卡爾曼濾波融合算法，形成改進(jìn)的卡爾曼濾波算法。對長時間尺度的監(jiān)測數(shù)據(jù)基于時間序列預(yù)處理段長度值的合理取值分析是本研究的重點。研究結(jié)果效用的評判采用相關(guān)度、信噪比和均方根誤差作為評價指標(biāo)。

1 卡爾曼濾波融合算法基本原理

卡爾曼濾波數(shù)據(jù)融合算法是一種采用自回歸的遞推運算來實現(xiàn)由兩個系統(tǒng)組成的觀測目標(biāo)進(jìn)行數(shù)據(jù)濾波融合的算法，適用于動態(tài)監(jiān)測數(shù)據(jù)的濾波與融合。其濾波的基本原理表述如下。

觀測對象本身過程控制系統(tǒng)的狀態(tài)方程為：

Xk=AXk-1+BUk+Wk

(1)

式中，Xk為k時刻的系統(tǒng)狀態(tài)變量；Xk-1為k-1時刻的系統(tǒng)狀態(tài)變量；Uk為k時刻的系統(tǒng)控制量；A、B為系統(tǒng)參數(shù)；Wk表示過程控制系統(tǒng)固有擾動，一般為均值高斯白噪聲，其噪聲協(xié)方差矩陣記為Q。

對觀測對象進(jìn)行觀測的量測方程為：

Zk=HXk+Vk

(2)

式中，Zk為k時刻的測量值；H為測量系統(tǒng)的參數(shù)；Vk表示傳感器量測系統(tǒng)固有偏差，其噪聲協(xié)方差矩陣記為R。

(3)

(4)

(5)

(6)

(7)

通過以上5個核心公式循環(huán)遞推可以實現(xiàn)動態(tài)監(jiān)測目標(biāo)信息的數(shù)據(jù)濾波，能夠?qū)ΜF(xiàn)場采集的數(shù)據(jù)進(jìn)行實時的更新和處理。

2 改進(jìn)卡爾曼濾波融合算法

2.1 基于時間序列的分段預(yù)處理

為了對傳統(tǒng)卡爾曼濾波算法進(jìn)行改進(jìn)，首先對采集到的數(shù)據(jù)序列進(jìn)行處理。假設(shè)采樣時間間隔為ΔT，那么采樣時間序列記為

[ΔT,2ΔT,…,NΔT,(N+1)ΔT,…,2NΔT,(2N+1)ΔT,…]選取其中一個時間序列預(yù)處理段長度值為N的采樣區(qū)間，其時間序列記為

[ΔT,2ΔT,…,NΔT]

其對應(yīng)的數(shù)據(jù)序列為

[x1,x2,…,xn]

將這個長度為N的采樣區(qū)間的數(shù)據(jù)序列進(jìn)行從小到大的排序，得到新的數(shù)據(jù)序列為

[y1,y2,…,yn]

其中y1≤y2≤…≤yn。

在這里的處理中，選取其中的半數(shù)中間數(shù)據(jù)取算術(shù)平均值作為該長度為N的采樣區(qū)間的數(shù)據(jù)序列的一個代表值。由于這個長度N是奇偶皆可的自然數(shù)，這個半數(shù)中位數(shù)的計算公式設(shè)計為較強適用性的公式，其具體的計算公式為

(8)

其中[]為向下取整。

經(jīng)過以上預(yù)處理，原來采樣時間間隔為ΔT的采樣時間序列為

[ΔT,2ΔT,…,NΔT,(N+1)ΔT,…,2NΔT,

(2N+1)ΔT,…]

擴大后的采樣時間間隔為

[Δt,2Δt,3Δt,…]

那么其相對應(yīng)的數(shù)據(jù)序列則為

[Y1,Y2,Y3,…]

從使得改進(jìn)算法設(shè)計的適用性更廣泛的角度而言，就低改進(jìn)區(qū)間N值的情況展開進(jìn)一步的分析討論，最終確定算法的改進(jìn)預(yù)處理廣泛適用公式為

(9)

其中[]為向下取整。

到此，對傳統(tǒng)卡爾曼濾波監(jiān)測數(shù)據(jù)序列的改進(jìn)達(dá)到了兩個方面的目的：(1)在這個處理上，選取的是半數(shù)中位數(shù)的算術(shù)平均值作為替代，這個操作舍棄了處理段數(shù)據(jù)序列中的極大值和極小值，即極大地減小了粗大誤差的影響；(2)將長度為N采樣區(qū)間上的N個數(shù)據(jù)用一個半數(shù)中位數(shù)的算術(shù)平均值取代，做到了把N個數(shù)據(jù)壓縮為一個數(shù)據(jù)的效果，即降低了數(shù)據(jù)冗余度。

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

通過預(yù)處理后，進(jìn)行卡爾曼濾波，需要對相關(guān)變量進(jìn)行同步更新。其需要更新的變量如下：(1)原來的k時刻經(jīng)N長度變換后變成了k＇時刻；(2)對應(yīng)的系統(tǒng)狀態(tài)變量和量測變量分別由Xk變?yōu)閄k′和Zk變?yōu)閆k′，其系統(tǒng)控制量由Uk變成Uk′；(3)預(yù)處理操作不影響過程控制系統(tǒng)固有擾動和量測系統(tǒng)固有偏差，即過程控制系統(tǒng)噪聲協(xié)方差矩陣Q和量測系統(tǒng)噪聲協(xié)方差矩陣R不變，即Q′=Q，R′=R。

以上對應(yīng)變量更新之后，便可以通過卡爾曼濾波算法的5個標(biāo)準(zhǔn)循環(huán)遞推公式(3)至公式(7)進(jìn)行數(shù)據(jù)濾波融合運算。

2.3 預(yù)處理段長度的確定與算法性能評價

為驗證改進(jìn)算法的合理性、有效性和進(jìn)步性，結(jié)合兩個方面的指標(biāo)：濾波融合數(shù)據(jù)與原始數(shù)據(jù)的相關(guān)度和信噪比，通過這兩個因素的值歸一化后的融合值fN來對N的合適取值進(jìn)行評估。其相關(guān)處理原理如下。

數(shù)據(jù)濾波結(jié)果與理想真值之間的相關(guān)度

(10)

式中，X和Y分別表示濾波結(jié)果數(shù)據(jù)和理想真值數(shù)據(jù)，E和D分別表示數(shù)學(xué)中的期望和方差。通常R值越接近1，表明變量之間的相關(guān)性越強，其處理效果越好；R值越接近0，表明變量之間的相關(guān)性越弱，其處理效果越差。

數(shù)據(jù)濾波結(jié)果與理想真值的信噪比(SNR)

(11)

式中，k表示采樣點的數(shù)目，N表示信號的長度，xk表示理想數(shù)據(jù)和，yk表示濾波后的數(shù)據(jù)。

將信噪比歸一化處理后和相關(guān)度進(jìn)行融合成為一個評價指標(biāo)，用來評估N值選擇的改進(jìn)效果，其公式定義為

(12)

(13)

為進(jìn)一步驗證該算法的有效性和優(yōu)異性，引入了數(shù)據(jù)濾波性能評價指標(biāo)均方根誤差(RMSE)，其定義為

(14)

3 驗證與分析

3.1 實驗數(shù)據(jù)獲取與初始化設(shè)置

這里僅針對數(shù)據(jù)濾波融合算法在CCUS封存環(huán)境下CO2泄漏監(jiān)測數(shù)據(jù)處理領(lǐng)域的可行性進(jìn)行研究，不對CO2泄漏的具體形成原因和擴散機理進(jìn)行探討。為獲取可供算法仿真驗證的實驗數(shù)據(jù)，搭建了如圖1所示地面CO2泄漏擴散試驗平臺來模擬CCUS封存環(huán)境下地面CO2泄漏的情況：設(shè)置一個尺寸為3.0 m×1.8 m×2.6 m(長×寬×高)的實驗密封艙用于進(jìn)行CO2擴散實驗，實驗艙預(yù)留一個0.7 m×1.0 m通風(fēng)口來模擬自然通風(fēng)的室內(nèi)環(huán)境，通風(fēng)口尺寸可調(diào)節(jié)。設(shè)計示意圖如圖1(a)所示。設(shè)計了如圖1(b)所示的1、2、3、4共4個監(jiān)測點，二氧化碳濃度傳感器采用COZIR-WX-20%，量程為0%～20%，上電穩(wěn)定時間小于10 s，在數(shù)據(jù)流模式下每秒采集并上傳2次數(shù)據(jù)，精度為±70 ppm，輸出模式為TTL電平輸出，各項參數(shù)均能滿足實驗要求。

圖1 數(shù)據(jù)獲取實驗圖

由于泄漏擴散試驗的試驗環(huán)境是設(shè)置在相對穩(wěn)定的“理想”環(huán)境下進(jìn)行的，因此，采用泄漏擴散試驗泄漏口為“DN100-0.6-大矩形”下監(jiān)測點1的監(jiān)測數(shù)據(jù)用于算法仿真驗證。在相同實驗條件下，通過改變尺寸可調(diào)節(jié)的通風(fēng)口通風(fēng)情況控制干擾條件(模擬環(huán)境監(jiān)測噪聲)形成含有干擾噪聲的數(shù)據(jù)，對其采用所提改進(jìn)算法進(jìn)行數(shù)據(jù)濾波，將處理結(jié)果與“理想”實驗數(shù)據(jù)進(jìn)行對比分析?？柭鼮V波算法的初始化設(shè)置，將變量A和H均設(shè)置為1；數(shù)據(jù)的噪聲干擾的來源主要源于環(huán)境和采集儀器，所以，根據(jù)傳感器的特性，將其測量數(shù)據(jù)的誤差協(xié)方差矩陣設(shè)置為r=0.07，過程控制的誤差協(xié)方差矩陣設(shè)置為q=0.04。

3.2 對不同長度N值的合理取值實驗分析

為了獲悉該算法的時間序列預(yù)處理段長度值N的合適值，采用的2.3節(jié)所述的算法對其進(jìn)行了仿真研究，其實驗結(jié)果如圖2所示。

從圖2(a)不同N值變化下的濾波融合結(jié)果信噪比中可以看出：N=1為經(jīng)典卡爾曼濾波算法處理下的數(shù)據(jù)濾波融合結(jié)果信噪比；當(dāng)N<18時，隨著N的增加，其信噪比也隨之增加，即處理效果在變好，改進(jìn)效果在顯現(xiàn)，這個增加的勢頭不會一直持續(xù)下去；當(dāng)N>18時，隨著N的繼續(xù)增加，就出現(xiàn)了一個下行的趨勢，這意味著這個改進(jìn)呈現(xiàn)了適得其反的效果；從工程學(xué)最優(yōu)化的角度考慮，選擇這個最優(yōu)的峰值N=18作為預(yù)處理長度段的設(shè)定值可以達(dá)到改進(jìn)的最佳效果。從圖2(b)不同N變化下濾波結(jié)果與理想真值的相關(guān)度情況可以看出，其呈現(xiàn)的趨勢的是：隨著N的增加，相關(guān)度有明顯提升；當(dāng)提升到一定的程度(N>7)之后，其相關(guān)度的變化便不再明顯。將圖2(a)和(b)的情況融合之后的結(jié)果為圖2(c)，通過圖2(c)的情況可以確定N的合適取值。

圖2 對N合理取值的探討

通過對以上實驗仿真結(jié)果的觀察，可以直觀地發(fā)現(xiàn)其改進(jìn)效果為：在經(jīng)典卡爾曼濾波的基礎(chǔ)上，隨著N的增加，其改進(jìn)效果有所提升；當(dāng)獲得一個峰值后，N繼續(xù)增加，其效果出現(xiàn)下降趨勢。因此，在預(yù)設(shè)允許壓縮的尺度內(nèi)有一個最優(yōu)可改進(jìn)預(yù)處理段長度值，在該實驗研究中，這個最優(yōu)的N值取18。

3.3 不同預(yù)處理長度N改進(jìn)算法與傳統(tǒng)算法的對比實驗

通過上面的仿真實驗可知，對于N取值為18時，其效果最佳。接下來，將選取幾個不同預(yù)處理長度值N的情況運用該算法進(jìn)行處理，以獲得更健全的效果對比信息。代表性N值的選取說明：由于N=1為未經(jīng)預(yù)處理改進(jìn)的經(jīng)典卡爾曼濾波算法，N=18為該實驗研究中的最佳改進(jìn)值，為使得加以對比的N值具有代表性，選取1到18的中間值N=9，以及從9到18之后向外延伸一倍單位距離的值N=27作為對比仿真實驗的N值。其仿真實驗結(jié)果如圖3所示。

為使實驗結(jié)果效果看起來更加清晰可辨，在實驗結(jié)果中對每一種N值的實驗仿真結(jié)果作了全局的展現(xiàn)，同時，將局部就行放大處理以更好地觀察對比效果；另外，還將其處理效果的量化評判數(shù)據(jù)匯總與表1中。從圖3和表1中可以看出：N=1為經(jīng)典卡爾曼濾波算法，在預(yù)處理階段的結(jié)果等于原含噪量測值，在卡爾曼濾波后的信噪比僅從21.97 dB提升至26.89 dB；N=9較N=1處理結(jié)果的信噪比有了較大的提升，約提升了29%，均方根誤差大大降低；N=18為該仿真對比實驗下最佳的預(yù)處理預(yù)設(shè)值，其相對于N=1和N=9的情況，在信噪比方面都有巨大提升，最終處理結(jié)果的信噪比已達(dá)40.82 dB；N=27是在最佳預(yù)處理值N=18之后的選值，其處理效果與3.2節(jié)證實的結(jié)論一致，其相對于N=18時的情況，處理結(jié)果的信噪比已經(jīng)下降了2.79 dB，在均方根誤差方面沒有出現(xiàn)相對較大的變化。

圖3 不同N取值下的數(shù)據(jù)濾波結(jié)果對比

表1 結(jié)果數(shù)據(jù)

4 結(jié)論

通過對數(shù)據(jù)濾波融合算法和CCUS封存環(huán)境CO2泄漏監(jiān)測數(shù)據(jù)處理的研究，得出了以下結(jié)論：

(1)研究了應(yīng)用于CCUS封存環(huán)境下地面CO2泄漏監(jiān)測的數(shù)據(jù)融合算法，提出了一種通過時間序列上的分段預(yù)處理去除原始數(shù)據(jù)的局部粗大誤差、并降低其冗余，然后將處理后的數(shù)據(jù)進(jìn)行卡爾曼濾波得到最佳監(jiān)測數(shù)據(jù)的改進(jìn)卡爾曼濾波算法，可為相關(guān)領(lǐng)域的數(shù)據(jù)融合算法提供理論支持；

(2)針對提出的算法，深入研究了不同N取值的確定，結(jié)果表明，隨著N的增加，其改進(jìn)效果有所提升，當(dāng)獲得一個峰值后，N繼續(xù)增加，其效果出現(xiàn)下降趨勢，因此，在預(yù)設(shè)的時間序列預(yù)處理段長度內(nèi)有一個最佳長度值；

(3)對不同N值的預(yù)處理改進(jìn)算法進(jìn)行了仿真實驗，對比分析表明，所選取的最佳預(yù)處理值相較于其他值有明顯的數(shù)據(jù)濾波處理效果，證明了改進(jìn)算法的有效性和優(yōu)異性。