劉芬范洪強(qiáng)呂濤李謙錢(qián)權(quán)

(1.上海大學(xué)計(jì)算機(jī)工程與科學(xué)學(xué)院,上海200444;2.上海大學(xué)材料科學(xué)與工程學(xué)院,上海200444;3.上海大學(xué)材料基因組工程研究院材料信息與數(shù)據(jù)科學(xué)中心,上海200444;4.重慶大學(xué)材料科學(xué)與工程學(xué)院國(guó)家鎂合金材料工程技術(shù)研究中心,重慶400044;5.之江實(shí)驗(yàn)室,浙江杭州311100)

在數(shù)據(jù)采集過(guò)程中噪聲不可避免.這些噪聲數(shù)據(jù)使得進(jìn)行機(jī)器學(xué)習(xí)、數(shù)據(jù)挖掘等任務(wù)時(shí)產(chǎn)生不準(zhǔn)確的輸出,因此對(duì)噪聲數(shù)據(jù)進(jìn)行處理十分必要.傳統(tǒng)的含噪聲數(shù)據(jù)處理方法主要是通過(guò)一系列算法檢測(cè)數(shù)據(jù)當(dāng)中的噪聲數(shù)據(jù)并刪除[1].這種直接刪除噪聲數(shù)據(jù)的方法會(huì)損失該條記錄的所有信息,因此更適用于數(shù)據(jù)量較大的數(shù)據(jù)集.當(dāng)遇到小樣本數(shù)據(jù)集時(shí),由于數(shù)據(jù)本身很寶貴,因此需要設(shè)計(jì)新的處理噪聲數(shù)據(jù)的方法,在充分利用數(shù)據(jù)信息的同時(shí)對(duì)噪聲數(shù)據(jù)進(jìn)行處理.

本研究提出了一種基于卡爾曼濾波的含噪聲小樣本數(shù)據(jù)處理方法.首先,結(jié)合領(lǐng)域知識(shí),利用數(shù)據(jù)集背后的物理模型或經(jīng)驗(yàn)公式,設(shè)計(jì)系統(tǒng)轉(zhuǎn)移矩陣和系統(tǒng)控制矩陣,并建立系統(tǒng)模型;然后,利用得到的系統(tǒng)模型預(yù)測(cè)數(shù)據(jù),并利用卡爾曼增益結(jié)合系統(tǒng)模型的預(yù)測(cè)結(jié)果和觀測(cè)數(shù)據(jù)生成修正數(shù)據(jù);最后,使用修正后的數(shù)據(jù)進(jìn)行數(shù)據(jù)挖掘或其他數(shù)據(jù)分析任務(wù).本方法基于卡爾曼濾波對(duì)數(shù)據(jù)進(jìn)行修正,達(dá)到含噪聲數(shù)據(jù)處理的效果,由于不需要?jiǎng)h除數(shù)據(jù),因此較適合小樣本數(shù)據(jù)集.

1 相關(guān)研究

含噪聲數(shù)據(jù)的處理方法主要有基于分箱的方法、基于模型的方法以及基于機(jī)器學(xué)習(xí)的方法.

1.1 基于分箱的方法

分箱方法是一種局部平滑法.將數(shù)據(jù)分成一個(gè)個(gè)“箱子”,然后將同一個(gè)箱子里面的數(shù)據(jù)進(jìn)行取平均值、取眾數(shù)、取中位數(shù)等方式進(jìn)行平滑,以此達(dá)到數(shù)據(jù)去噪的效果,獲得穩(wěn)定的數(shù)據(jù)特征.常見(jiàn)的分箱方法有等深分箱法和等寬分箱法[2-3].等深分箱法是將數(shù)據(jù)平均分成n個(gè)箱子,每個(gè)箱子中的數(shù)據(jù)數(shù)目是相等的.等寬分箱法是將數(shù)據(jù)的值區(qū)間等分,每個(gè)箱子中的數(shù)據(jù)區(qū)間大小相同.

除了上述基本分箱法以外,還可以根據(jù)任務(wù)自定義分箱方法.Kerber[2]提出了一種ChiMerge分箱方法,用χ2統(tǒng)計(jì)量將連續(xù)屬性值進(jìn)行離散化.傅濤等[4]結(jié)合了基于分箱的Fuzzy C Means算法用于檢測(cè)網(wǎng)絡(luò)入侵?jǐn)?shù)據(jù).該方法用分箱方式劃分?jǐn)?shù)據(jù)集,解決了原Fuzzy C Means聚類(lèi)算法頻繁更換聚類(lèi)中心的問(wèn)題.

1.2 基于模型的方法

基于模型的方法是使用原數(shù)據(jù)訓(xùn)練一個(gè)回歸模型,然后用訓(xùn)練好的模型對(duì)目標(biāo)變量進(jìn)行預(yù)測(cè).若觀測(cè)值與預(yù)測(cè)值相差過(guò)大,則被視為噪聲數(shù)據(jù).對(duì)于時(shí)序數(shù)據(jù)主要有自回歸(autoregressive,AR)模型、自回歸移動(dòng)平均(autoregressive integrated moving average,ARIMA)模型[5]等方法.Li[6]提出了一種基于支持向量回歸模型的圖像去噪方法,在帶有真實(shí)標(biāo)簽的含噪聲圖像上進(jìn)行訓(xùn)練,得到了多個(gè)支持向量機(jī)模型及其權(quán)重.然后,利用訓(xùn)練好的支持向量機(jī)及其權(quán)值,對(duì)受到不同程度隨機(jī)噪聲影響的圖像進(jìn)行逐像素點(diǎn)去噪.L′opez-Rubio等[7]提出了一種基于核回歸的磁共振三維圖像去噪算法,改進(jìn)了傳統(tǒng)核回歸方法,使其適應(yīng)符合Rician分布的噪聲.Huang等[8]提出了基于低秩判別的最小二乘回歸模型,用來(lái)消除標(biāo)簽空間的噪聲.Bai等[9]利用圖像序列數(shù)據(jù)中不同像素的振幅變化特征,對(duì)每個(gè)像素的相應(yīng)時(shí)間序列施加高斯過(guò)程(Gaussian process,GP)模型,并通過(guò)GP回歸進(jìn)行圖像序列數(shù)據(jù)去噪.

1.3 基于機(jī)器學(xué)習(xí)的方法

基于機(jī)器學(xué)習(xí)方法是用機(jī)器學(xué)習(xí)模型對(duì)數(shù)據(jù)進(jìn)行處理,然后通過(guò)某種方式檢測(cè)噪聲數(shù)據(jù).何添等[10]提出了基于維諾圖(Voronoi diagram,VOD)和K-Means算法的噪聲數(shù)據(jù)處理方法,通過(guò)K-Means算法對(duì)原數(shù)據(jù)進(jìn)行聚類(lèi),再通過(guò)VOD方法檢測(cè)噪聲數(shù)據(jù).Dong等[11]提出了一種基于聚類(lèi)的圖像數(shù)據(jù)的稀疏表示去噪算法.通過(guò)變分方法結(jié)合局部圖像模型和非局部圖像模型各自的優(yōu)點(diǎn):一方面通過(guò)構(gòu)建一個(gè)基本函數(shù)字典提高稀疏性;另一方面通過(guò)聚類(lèi)將稀疏性與圖像源的自相似性聯(lián)系起來(lái).胡嬌嬌[12]提出了基于卷積神經(jīng)網(wǎng)絡(luò)的噪聲處理模型,將噪聲數(shù)據(jù)處理轉(zhuǎn)化為分類(lèi)或預(yù)測(cè)問(wèn)題.Zhang等[13]提出了一種基于神經(jīng)網(wǎng)絡(luò)的小波域圖像去噪方法.首先對(duì)帶噪圖像進(jìn)行小波變換得到子帶,然后用訓(xùn)練好的分層神經(jīng)網(wǎng)絡(luò)對(duì)每個(gè)子帶進(jìn)行小波變換,得到去噪后的小波系數(shù),最后對(duì)去噪后的小波系數(shù)進(jìn)行反變換得到去噪后的圖像.Mozhde等[14]提出了一種基于離散小波變換和人工神經(jīng)網(wǎng)絡(luò)(discrete wavelet transform-artificial neural network,DWT-ANN)的自適應(yīng)算法,過(guò)濾嘈雜環(huán)境中的肺聲信號(hào)數(shù)據(jù).該算法結(jié)合離散小波變換的多分辨率特性與人工神經(jīng)網(wǎng)絡(luò)的非線性特性,設(shè)計(jì)了不同信噪比下的獨(dú)立模型.通過(guò)集合這些獨(dú)立模型來(lái)消除輸入信號(hào)的信噪比信息,便于后續(xù)去噪操作.Yao等[15]提出了一種多尺度殘差融合網(wǎng)絡(luò)(multi-scale residual net,MRF-Net)用于圖像數(shù)據(jù)的去噪.該網(wǎng)絡(luò)綜合了圖像的空間和上下文信息,首先通過(guò)擴(kuò)張卷積層擴(kuò)大網(wǎng)絡(luò)的接受域,然后進(jìn)行多尺度特征提取形成多級(jí)特征圖,最后將多級(jí)特征圖進(jìn)行殘差融合生成殘差圖像,用于去除圖像噪聲.

現(xiàn)有的噪聲數(shù)據(jù)處理方法還存在一些問(wèn)題.基于分箱的方法要求處理的數(shù)據(jù)具有一定的規(guī)模,若數(shù)據(jù)樣本過(guò)少,則數(shù)據(jù)分箱之后取平均值、眾數(shù)或中位數(shù)等方式的平滑效果就會(huì)大打折扣.這是因?yàn)樵跀?shù)據(jù)量較少的情況下,很難通過(guò)統(tǒng)計(jì)的方式得到數(shù)據(jù)的分布規(guī)律.基于模型的方法和基于機(jī)器學(xué)習(xí)的方法本質(zhì)上是通過(guò)一定的手段檢測(cè)出噪聲數(shù)據(jù),然后將噪聲數(shù)據(jù)去除以達(dá)到降噪的目的.但直接去除整條噪聲記錄會(huì)損失大量信息.尤其是面對(duì)小樣本噪聲數(shù)據(jù)時(shí),若直接去除某條噪聲數(shù)據(jù)會(huì)損失很大比例的信息資源.因此,針對(duì)小樣本噪聲數(shù)據(jù)需要設(shè)計(jì)新的去噪方法.

本工作提出了一種基于卡爾曼濾波的噪聲數(shù)據(jù)處理方法,結(jié)合了模型方法中的回歸以及分箱方法中的平滑技術(shù).卡爾曼濾波處理方法將數(shù)據(jù)視為系統(tǒng)的狀態(tài),通過(guò)數(shù)據(jù)的經(jīng)驗(yàn)公式建立系統(tǒng)的狀態(tài)轉(zhuǎn)移模型和狀態(tài)控制模型,以此來(lái)預(yù)測(cè)系統(tǒng)的狀態(tài),得到系統(tǒng)預(yù)測(cè)數(shù)據(jù).最后,通過(guò)觀測(cè)數(shù)據(jù)修正系統(tǒng)預(yù)測(cè)數(shù)據(jù),得到最后的修正數(shù)據(jù).該修正數(shù)據(jù)即為降噪之后的數(shù)據(jù).卡爾曼濾波處理方法通過(guò)修正數(shù)據(jù)的方式來(lái)達(dá)到降噪的效果,不用刪除數(shù)據(jù),更加適用于小樣本數(shù)據(jù).

2 基于卡爾曼濾波的含噪聲小樣本數(shù)據(jù)處理模型

2.1 卡爾曼濾波

卡爾曼濾波是一種線性二次估計(jì)方法[16].該方法利用系統(tǒng)的動(dòng)態(tài)模型,通過(guò)系統(tǒng)的控制輸入數(shù)據(jù)以及觀測(cè)數(shù)據(jù)來(lái)對(duì)系統(tǒng)的狀態(tài)量進(jìn)行估計(jì).該方法相比只使用動(dòng)態(tài)模型獲得的估計(jì)值更準(zhǔn)確,是一種常用的模型數(shù)據(jù)和觀測(cè)數(shù)據(jù)融合的算法.卡爾曼濾波的基本思想如圖1所示.

首先,卡爾曼濾波器根據(jù)k-1時(shí)刻得到的系統(tǒng)狀態(tài)的最終估計(jì)^xk-1(見(jiàn)圖1中的黑色曲線),產(chǎn)生k時(shí)刻系統(tǒng)狀態(tài)的預(yù)測(cè)估計(jì)^x-k(見(jiàn)圖1中的綠色曲線).然后，根據(jù)k時(shí)刻系統(tǒng)的測(cè)量估計(jì)yk(見(jiàn)圖1中的藍(lán)色直線),將預(yù)測(cè)估計(jì)以及測(cè)量估計(jì)根據(jù)相應(yīng)的“權(quán)重”進(jìn)行融合，產(chǎn)生k時(shí)刻系統(tǒng)狀態(tài)的最終估計(jì)(見(jiàn)圖1中的紅色曲線).這里的“權(quán)重”是由預(yù)測(cè)系統(tǒng)的不確定性和觀測(cè)系統(tǒng)的不確定性綜合得出.以上過(guò)程在每一個(gè)時(shí)間步重復(fù)，得到每一個(gè)時(shí)間步系統(tǒng)狀態(tài)的最終估計(jì).最終狀態(tài)位于預(yù)測(cè)狀態(tài)和測(cè)量狀態(tài)之間，比任何一個(gè)單獨(dú)狀態(tài)都具有更好的不確定性估計(jì).

圖1 卡爾曼濾波示意圖Fig.1 Diagram of Kalman filter

在卡爾曼濾波中，一般線性離散系統(tǒng)可以表示為

式中:X(k)代表k時(shí)刻的系統(tǒng)狀態(tài)向量;Z(k)代表觀測(cè)向量;u(k)代表輸入向量;v(k)代表m維測(cè)量噪聲向量;A(k+1,k)代表從k時(shí)刻到k+1時(shí)刻的系統(tǒng)狀態(tài)轉(zhuǎn)移矩陣;B(k+1,k)代表從k時(shí)刻到k+1時(shí)刻的系統(tǒng)控制矩陣;H(k+1)代表k+1時(shí)刻的預(yù)測(cè)輸出轉(zhuǎn)移矩陣.式(1)被稱為狀態(tài)方程,式(2)被稱為測(cè)量方程.

卡爾曼濾波算法包括預(yù)測(cè)和更新兩個(gè)階段,其流程如圖2所示.在預(yù)測(cè)階段,卡爾曼濾波器根據(jù)式(3)產(chǎn)生預(yù)測(cè)模型對(duì)當(dāng)前狀態(tài)變量的估計(jì),根據(jù)式(4)產(chǎn)生當(dāng)前系統(tǒng)噪聲的協(xié)方差,即

圖2 卡爾曼濾波算法迭代圖Fig.2 Iteration diagram of Kalman filter algorithm

當(dāng)系統(tǒng)觀察到下一個(gè)觀測(cè)數(shù)據(jù)(觀測(cè)數(shù)據(jù)會(huì)有一定的誤差,如隨機(jī)噪聲等)時(shí),進(jìn)入到更新階段.卡爾曼增益的計(jì)算公式為

式中:K表示卡爾曼增益;H表示預(yù)測(cè)輸出轉(zhuǎn)移矩陣,即從預(yù)測(cè)估計(jì)到測(cè)量估計(jì)的計(jì)算方式;R代表測(cè)量誤差協(xié)方差.

根據(jù)式(6)更新當(dāng)前系統(tǒng)狀態(tài)量的估計(jì),根據(jù)式(7)計(jì)算下一時(shí)刻需要用到的估計(jì)協(xié)方差矩陣,即

式中,z代表觀測(cè)數(shù)據(jù).

卡爾曼算法是迭代的,不斷重復(fù)預(yù)測(cè)和更新步驟.它可以實(shí)時(shí)運(yùn)行,僅使用當(dāng)前輸入測(cè)量值和先前計(jì)算的狀態(tài)及其不確定性矩陣,不需要額外的歷史信息.

2.2 擴(kuò)展卡爾曼濾波

基本卡爾曼濾波要求必須是線性系統(tǒng),但大部分都是非線性系統(tǒng),其中的“非線性性質(zhì)”既可能存在于過(guò)程模型(process model)中,又可能存在于觀測(cè)模型(observation model)中,或者二者都有.當(dāng)處理非線性系統(tǒng)時(shí),常用擴(kuò)展卡爾曼濾波[17].

在擴(kuò)展卡爾曼濾波中,系統(tǒng)的狀態(tài)轉(zhuǎn)移模型和觀測(cè)模型如式(8)和(9)所示,二者只要是可微函數(shù)即可,并不要求是線性函數(shù).

式中:函數(shù)f用于從過(guò)去的估計(jì)值中計(jì)算預(yù)測(cè)的狀態(tài);函數(shù)h用于以預(yù)測(cè)的狀態(tài)計(jì)算預(yù)測(cè)的測(cè)量值.f和h不能直接計(jì)算協(xié)方差,需要計(jì)算它們的偏導(dǎo)矩陣(Jacobian矩陣).狀態(tài)轉(zhuǎn)移模型和觀測(cè)模型的Jacobian矩陣計(jì)算如下:

在解決了f和h不能直接計(jì)算協(xié)方差的問(wèn)題后,擴(kuò)展卡爾曼濾波的迭代方程如式(12)～(16)所示.在預(yù)測(cè)階段,根據(jù)式(12)計(jì)算預(yù)測(cè)模型對(duì)當(dāng)前狀態(tài)變量的估計(jì),根據(jù)式(13)計(jì)算當(dāng)前系統(tǒng)噪聲的協(xié)方差,即當(dāng)系統(tǒng)觀察到下一個(gè)觀測(cè)數(shù)據(jù)時(shí),進(jìn)入到更新階段.根據(jù)式(14)計(jì)算卡爾曼增益,根據(jù)式(15)修正模型得到的數(shù)據(jù),根據(jù)式(16)更新系統(tǒng)噪聲的協(xié)方差P,為下一步迭代做準(zhǔn)備.與卡爾曼濾波一樣,擴(kuò)展卡爾曼算法也是迭代的,不斷重復(fù)預(yù)測(cè)和更新步驟.

觀察卡爾曼濾波的計(jì)算公式(式(3)～(7))以及擴(kuò)展卡爾曼濾波的計(jì)算公式(式(12)～(16))可以看出,二者幾乎是一樣的,只是由于非線性系統(tǒng)無(wú)法直接計(jì)算狀態(tài)轉(zhuǎn)移模型和觀測(cè)模型的協(xié)方差矩陣,因此將卡爾曼濾波計(jì)算過(guò)程中的狀態(tài)轉(zhuǎn)移矩陣A和狀態(tài)控制矩陣B替換為相應(yīng)的偏導(dǎo)矩陣Fk和Hk.觀察Fk和Hk的計(jì)算公式可以看出,擴(kuò)展卡爾曼濾波相當(dāng)于把系統(tǒng)在k時(shí)刻按照k-1時(shí)刻的方程進(jìn)行線性化,即系統(tǒng)當(dāng)前時(shí)刻的值等于前一時(shí)刻的值加上前后兩時(shí)刻值之間的變化量.實(shí)現(xiàn)方法為不斷迭代計(jì)算系統(tǒng)在前一時(shí)刻的偏導(dǎo)矩陣值.

2.3 基于卡爾曼濾波的含噪聲小樣本數(shù)據(jù)處理模型

本工作提出的基于卡爾曼濾波的含噪聲小樣本數(shù)據(jù)處理方法的流程如圖3所示.首先,對(duì)原數(shù)據(jù)進(jìn)行處理,用線性模型對(duì)數(shù)據(jù)進(jìn)行擬合,得到初步的線性模型參數(shù).其次,根據(jù)線性模型參數(shù),建立卡爾曼濾波模型中的系統(tǒng)狀態(tài)轉(zhuǎn)移矩陣和系統(tǒng)狀態(tài)控制矩陣.再次,根據(jù)得到的系統(tǒng)狀態(tài)轉(zhuǎn)移矩陣和系統(tǒng)狀態(tài)控制矩陣建立系統(tǒng)模型,結(jié)合系統(tǒng)模型得出的結(jié)果與觀測(cè)數(shù)據(jù),對(duì)原數(shù)據(jù)進(jìn)行修正.最后,根據(jù)修正得到的數(shù)據(jù)進(jìn)行后續(xù)的數(shù)據(jù)挖掘應(yīng)用.

圖3 基于卡爾曼濾波的含噪聲小樣本數(shù)據(jù)處理流程Fig.3 Flow chart of noisy sample data processing based on Kalman filter

3 含噪聲材料數(shù)據(jù)處理實(shí)驗(yàn)

3.1 數(shù)據(jù)集

金屬腐蝕在自然界中廣泛存在,研究耐腐蝕金屬已成為材料研究領(lǐng)域的熱點(diǎn)[18-20].耐候鋼是一種通過(guò)在普通鋼中添加一定量合金元素所制成的低合金鋼,具有較強(qiáng)的耐腐蝕性[21].耐候鋼在初期和普通碳鋼一樣也會(huì)銹蝕.但隨著合金表面的腐蝕發(fā)展,可以在表面形成一層致密的非晶態(tài)銹層組織,這層致密銹層是其抗大氣腐蝕的關(guān)鍵[22-23].低成本、高強(qiáng)度、高耐腐蝕性的耐候鋼開(kāi)發(fā)已成為熱點(diǎn)[24],因此對(duì)耐候鋼腐蝕失效數(shù)據(jù)的研究十分必要.

本工作以實(shí)驗(yàn)室腐蝕加速模擬所獲得的BC500耐候鋼腐蝕失效數(shù)據(jù)集為例,提出了一種基于卡爾曼濾波的含噪聲小樣本數(shù)據(jù)處理方法,用于處理耐候鋼腐蝕失效數(shù)據(jù)中的數(shù)據(jù)噪聲,便于后續(xù)的腐蝕失效數(shù)據(jù)擬合等任務(wù).耐候鋼腐蝕失效數(shù)據(jù)集包含90條記錄,每30條記錄為一組,共分為3組,編號(hào)為Data1、Data2和Data3,分別代表一次耐候鋼腐蝕實(shí)驗(yàn).數(shù)據(jù)集共含有7維特征:溫度、濕度、工業(yè)污染物(Na2SO4、NaHCO3和NaNO3)濃度、光照、腐蝕時(shí)間.目標(biāo)變量為耐候鋼的腐蝕速率.任務(wù)目標(biāo)為處理該數(shù)據(jù)集中的噪聲,為后續(xù)擬合耐候鋼的腐蝕速率做準(zhǔn)備.

3.2 實(shí)驗(yàn)設(shè)置

耐候鋼的大氣腐蝕動(dòng)力學(xué)[25-28]可以表示為

式中:ΔW為腐蝕增重(mg);N為腐蝕時(shí)間(h),當(dāng)腐蝕時(shí)間N=1時(shí),F=ΔW.F和n為常數(shù),其中F是作為測(cè)試樣品初始腐蝕速率的度量;指數(shù)n反映了腐蝕動(dòng)力學(xué)的特性.n＜1表示腐蝕進(jìn)程受到抑制;n＞1表示腐蝕進(jìn)程受到加速;n=1表示腐蝕速率恒定.n值越小,表明該階段的腐蝕強(qiáng)度越小.n的取值與耐候鋼所處的環(huán)境密切相關(guān).

對(duì)式(17)等號(hào)兩邊各取自然對(duì)數(shù),推導(dǎo)出式(18).從式(18)可以看出,對(duì)耐候鋼的大氣腐蝕增重ΔW和腐蝕時(shí)間N分別取對(duì)數(shù)后,二者呈線性關(guān)系,即

在本實(shí)驗(yàn)數(shù)據(jù)集所涉及到的腐蝕因素(特征)中,溫度、濕度、工業(yè)污染物濃度、光照這些特征會(huì)影響材料的腐蝕速率,即影響式(17)中n的大小.將數(shù)據(jù)中的腐蝕速率與對(duì)應(yīng)的時(shí)間相乘,可以得到式(17)中的腐蝕增重ΔW.根據(jù)式(18)所示的lnΔW和lnN的關(guān)系,可以用線性模型擬合數(shù)據(jù)初步建立腐蝕動(dòng)力學(xué)方程.假設(shè)線性模型的擬合結(jié)果為

式中:x代表輸入(lnN);y代表輸出(lnΔW);k和b代表線性方程的斜率和截距.根據(jù)線性擬合結(jié)果,可計(jì)算得出F和n的值,即

綜上,根據(jù)線性擬合結(jié)果,已建立了一個(gè)耐候鋼腐蝕增重的系統(tǒng)模型,其中大氣腐蝕增重ΔW被視為系統(tǒng)狀態(tài);系統(tǒng)的狀態(tài)轉(zhuǎn)移矩陣A設(shè)置成值為1的單元素矩陣;系統(tǒng)的控制轉(zhuǎn)移矩陣B設(shè)置成值為F的單元素矩陣,F的值由式(21)給出.由于實(shí)驗(yàn)所用的數(shù)據(jù)集為每12 h采集一次,因此系統(tǒng)的輸入設(shè)置為12n,n的值由式(20)給出.系統(tǒng)模型可以描述為

由于本實(shí)驗(yàn)數(shù)據(jù)集中觀測(cè)變量即為系統(tǒng)狀態(tài),因此將預(yù)測(cè)輸出轉(zhuǎn)移矩陣H設(shè)置成值為1的單元素矩陣.考慮到測(cè)量會(huì)造成噪聲干擾,因此引入觀測(cè)噪聲σz,系統(tǒng)的測(cè)量方程為

式中,zk為觀測(cè)變量,噪聲服從均值為0的正態(tài)分布.

初始化估計(jì)協(xié)方差矩陣P和預(yù)測(cè)模型的誤差協(xié)方差矩陣Q.根據(jù)系統(tǒng)模型公式,本工作所用的P和Q矩陣均為單元素矩陣.按照上述方法迭代進(jìn)行卡爾曼濾波,其預(yù)測(cè)階段的計(jì)算方法如式(22)和(24)所示,更新階段的計(jì)算方法如式(25)～(27)所示.

經(jīng)過(guò)卡爾曼濾波處理之后,原數(shù)據(jù)已經(jīng)過(guò)修正,可用于下游任務(wù)模型.本工作所用的下游任務(wù)模型為差分整合移動(dòng)平均自回歸(autoregressive integrated moving average,ARIMA)模型和隨機(jī)森林(random forest,RF)模型,用來(lái)預(yù)測(cè)材料數(shù)據(jù)的腐蝕增重.將時(shí)間當(dāng)作輸入,預(yù)測(cè)該時(shí)間下材料的腐蝕增重.

3.3 實(shí)驗(yàn)結(jié)果與分析

3.3.1 腐蝕動(dòng)力學(xué)方程建模及卡爾曼濾波效果

本實(shí)驗(yàn)展示了通過(guò)卡爾曼濾波處理含噪聲材料數(shù)據(jù),來(lái)建立腐蝕動(dòng)力學(xué)方程模型的結(jié)果,以及數(shù)據(jù)經(jīng)過(guò)卡爾曼濾波處理后的修正效果.

表1展示了Data1、Data2和Data3共3組數(shù)據(jù)擬合的腐蝕動(dòng)力學(xué)方程中的各項(xiàng)參數(shù),其中b和k是線性擬合結(jié)果,n和A是系統(tǒng)狀態(tài)轉(zhuǎn)移矩陣及系統(tǒng)狀態(tài)控制矩陣的相關(guān)參數(shù).可以看出,3組數(shù)據(jù)的n指數(shù)均小于1,可見(jiàn)耐候鋼在腐蝕進(jìn)程中其腐蝕速率逐步受到抑制,腐蝕強(qiáng)度逐漸變小.圖4為Data1、Data2和Data3擬合的腐蝕動(dòng)力學(xué)方程的可視化結(jié)果以及經(jīng)卡爾曼濾波處理后的效果.可以看出,卡爾曼濾波處理對(duì)含噪聲數(shù)據(jù)具有一定的平滑效果,修正后的數(shù)據(jù)(橙色線)處于模型數(shù)據(jù)(藍(lán)色線)和觀測(cè)數(shù)據(jù)(綠色線)之間,且在一定程度上平滑了原數(shù)據(jù)中起伏過(guò)大(實(shí)驗(yàn)中引入的誤差)的地方.

圖4 數(shù)據(jù)集Data1、Data2和Data3的實(shí)驗(yàn)結(jié)果Fig.4 Experimental results of Data1,Data2,and Data3

表1 耐候鋼數(shù)據(jù)腐蝕動(dòng)力學(xué)方程參數(shù)擬合結(jié)果Table 1 Parameters estimation results of corrosion kinetics for weathering steel

3.3.2 不同系統(tǒng)噪聲協(xié)方差對(duì)卡爾曼濾波效果的影響

卡爾曼濾波模型中有一個(gè)非常重要的參數(shù)Q(系統(tǒng)噪聲協(xié)方差),用來(lái)表示狀態(tài)轉(zhuǎn)換矩陣與實(shí)際過(guò)程之間的誤差.Q越小代表整個(gè)系統(tǒng)更加相信模型的結(jié)果,Q越大則代表整個(gè)系統(tǒng)更加相信觀測(cè)數(shù)據(jù)的結(jié)果.圖5顯示了不同Q矩陣下數(shù)據(jù)集Data1的卡爾曼濾波結(jié)果.可以看出:當(dāng)Q越小時(shí)整個(gè)系統(tǒng)更加偏向于模型的結(jié)果,因此橙色線更加接近藍(lán)色線;當(dāng)Q越大時(shí)整個(gè)系統(tǒng)更加偏向于觀測(cè)數(shù)據(jù),此時(shí)橙色線更加接近綠色線.不同的Q矩陣對(duì)數(shù)據(jù)修正的結(jié)果影響較大,因此需要選擇合適的Q矩陣.

圖5 數(shù)據(jù)集Data1在不同Q矩陣值下的卡爾曼濾波結(jié)果比較Fig.5 Data processing results of Data1 data under different Q matrix of Kalman filter

3.3.3 卡爾曼濾波處理噪聲材料數(shù)據(jù)后的建模實(shí)驗(yàn)

為了驗(yàn)證本工作提出的方法對(duì)噪聲數(shù)據(jù)的處理效果,分別測(cè)試了未經(jīng)處理的數(shù)據(jù)和經(jīng)過(guò)卡爾曼濾波處理的數(shù)據(jù)的實(shí)驗(yàn)效果.將時(shí)間當(dāng)作輸入,分別用ARIMA模型和RF模型來(lái)預(yù)測(cè)該時(shí)間下材料的腐蝕增重.

在ARIMA任務(wù)模型中,將每組數(shù)據(jù)(Data1、Data2和Data3)的前20條作為訓(xùn)練集,后10條作為測(cè)試集,用訓(xùn)練好的ARIMA模型預(yù)測(cè)后10個(gè)腐蝕增重,并計(jì)算結(jié)果的決定系數(shù)R2.在RF模型中,將每組數(shù)據(jù)的前25條作為訓(xùn)練集,后5條作為測(cè)試集,用訓(xùn)練好的RF模型預(yù)測(cè)后5個(gè)腐蝕增重,并計(jì)算結(jié)果的R2.ARIMA模型和RF模型的實(shí)驗(yàn)效果如表2所示,其中“Raw-Raw”列表示用未經(jīng)過(guò)卡爾曼濾波處理的數(shù)據(jù)訓(xùn)練的模型,預(yù)測(cè)未經(jīng)過(guò)卡爾曼濾波處理的數(shù)據(jù)的R2,“Kalman-Raw”列表示用經(jīng)過(guò)卡爾曼濾波處理的數(shù)據(jù)訓(xùn)練的模型,預(yù)測(cè)未經(jīng)過(guò)卡爾曼濾波處理的數(shù)據(jù)的R2.

表2 經(jīng)卡爾曼濾波處理含噪聲數(shù)據(jù)后的ARIMA回歸和RF回歸結(jié)果Table 2 ARIMA regression and RF regression results of noisy data processed by Kalman filter

由表2可以看出:在ARIMA模型中,“Raw-Raw”列中Data1和Data3的R2要小于Data2,這是由于Data1和Data3的觀測(cè)數(shù)據(jù)具有更高的波動(dòng)性(見(jiàn)圖4),數(shù)據(jù)噪聲較大,因此模型預(yù)測(cè)的R2較小,效果更差;“Kalman-Raw”列中Data1和Data3的R2要大于Data2,Data1的R2由0.863 0增加到0.943 5,Data3的R2由0.851 1增加到0.945 7,這說(shuō)明經(jīng)過(guò)卡爾曼濾波處理后,Data1和Data3數(shù)據(jù)中的噪聲被平滑,模型的結(jié)果變好.Data2數(shù)據(jù)經(jīng)過(guò)卡爾曼濾波處理后,其模型的R2反而變小.這是因?yàn)檫@組數(shù)據(jù)原本的噪聲就比較小,而經(jīng)卡爾曼濾波處理后反而引入了新噪聲,對(duì)原數(shù)據(jù)產(chǎn)生了影響,模型結(jié)果變差.因此,卡爾曼濾波處理對(duì)本身就存在較大噪聲的數(shù)據(jù)才會(huì)產(chǎn)生較為積極的影響.在RF模型中,“Kalman-Raw”列中的R2均大于“Raw-Raw”列.這進(jìn)一步驗(yàn)證了卡爾曼濾波處理對(duì)含噪聲數(shù)據(jù)會(huì)產(chǎn)生較為積極的影響(兩個(gè)模型的R2平均提升了6.4%).對(duì)比表2中的數(shù)據(jù)可以看出,卡爾曼濾波處理對(duì)時(shí)間序列回歸模型的影響較大,這是因?yàn)闀r(shí)序模型對(duì)數(shù)據(jù)噪聲更為敏感.

3.3.4 擴(kuò)展卡爾曼濾波處理噪聲材料數(shù)據(jù)后的建模實(shí)驗(yàn)

本實(shí)驗(yàn)測(cè)試了擴(kuò)展卡爾曼濾波對(duì)于噪聲數(shù)據(jù)的處理效果.首先按照3.2節(jié)中所述方法建模得到了式(17)中F、n的值.假設(shè)式(17)的函數(shù)關(guān)系為f,則可得到

將腐蝕增重ΔW視為系統(tǒng)狀態(tài),在k時(shí)刻,擴(kuò)展卡爾曼濾波把系統(tǒng)按照k-1時(shí)刻的方程進(jìn)行線性化以適應(yīng)非線性系統(tǒng).將f(N)按照泰勒公式在Nk-1處展開(kāi),則系統(tǒng)模型描述為

式中:xk為系統(tǒng)狀態(tài);N代表時(shí)間.系統(tǒng)的測(cè)量方程為

式中:zk為觀測(cè)變量,噪聲服從均值為0的正態(tài)分布.

由式(29)和(30)可知,系統(tǒng)狀態(tài)轉(zhuǎn)移模型和觀測(cè)模型分別為非線性模型和線性模型,因此相應(yīng)的Jacobian矩陣分別為f′(Nk-1)和值為1的單元素矩陣.已知系統(tǒng)的狀態(tài)轉(zhuǎn)移模型、觀測(cè)模型以及二者相應(yīng)的Jacobian矩陣,然后按照式(12)～(16)進(jìn)行擴(kuò)展卡爾曼濾波迭代,對(duì)原數(shù)據(jù)進(jìn)行擴(kuò)展卡爾曼濾波處理,用于同3.3.3節(jié)的下游任務(wù).

分別測(cè)試了經(jīng)過(guò)擴(kuò)展卡爾曼濾波處理過(guò)后的數(shù)據(jù)在ARIMA模型和RF模型上的實(shí)驗(yàn)結(jié)果,如表3所示.可以看出,雖然經(jīng)過(guò)擴(kuò)展卡爾曼濾波處理之后的含噪聲數(shù)據(jù),其R2也有所增加,但對(duì)比3.3.3節(jié)的實(shí)驗(yàn)結(jié)果,擴(kuò)展卡爾曼濾波處理的效果略差于卡爾曼濾波處理,但差別不大(兩個(gè)模型的R2平均提升了4.9%).這是由于本實(shí)驗(yàn)中的材料數(shù)據(jù)取自然對(duì)數(shù)后符合線性關(guān)系,因此經(jīng)擴(kuò)展卡爾曼濾波處理后效果不明顯.

表3 經(jīng)擴(kuò)展卡爾曼濾波處理含噪聲數(shù)據(jù)后的ARIMA回歸和RF回歸結(jié)果Table 3 ARIMA regression and RF regression results of noisy data processed by extended Kalman filter

4 結(jié)束語(yǔ)

本工作提出了一種基于卡爾曼濾波的含噪聲小樣本數(shù)據(jù)處理方法.首先,結(jié)合領(lǐng)域知識(shí),用數(shù)據(jù)集背后的物理模型或經(jīng)驗(yàn)公式擬合初步模型,得到系統(tǒng)狀態(tài)轉(zhuǎn)移矩陣和系統(tǒng)控制矩陣.然后,根據(jù)得到的系統(tǒng)狀態(tài)轉(zhuǎn)移矩陣和系統(tǒng)狀態(tài)控制矩陣建立系統(tǒng)模型.最后,綜合系統(tǒng)模型得出的結(jié)果與觀測(cè)數(shù)據(jù),得到最后的修正數(shù)據(jù).

與傳統(tǒng)的噪聲數(shù)據(jù)處理方法不同,卡爾曼濾波處理方法不會(huì)直接刪掉噪聲數(shù)據(jù),而是通過(guò)系統(tǒng)模型的計(jì)算結(jié)果與觀測(cè)數(shù)據(jù)融合,對(duì)原數(shù)據(jù)進(jìn)行平滑,從而達(dá)到去噪的效果.該方法對(duì)小樣本數(shù)據(jù)更加友好,通過(guò)修正數(shù)據(jù)既達(dá)到了降噪的效果,又不會(huì)因?yàn)閯h去噪聲而損失寶貴的材料數(shù)據(jù)樣本.BC500腐蝕數(shù)據(jù)的實(shí)驗(yàn)結(jié)果表明,卡爾曼濾波處理含噪聲數(shù)據(jù)取得了較好的效果,Data1預(yù)測(cè)的R2由0.863 0增加到0.943 5,Data3的R2由0.851 1增加到0.945 7.實(shí)驗(yàn)結(jié)果還表明,卡爾曼濾波處理對(duì)時(shí)間序列模型影響更大,在時(shí)間序列模型上取得了更好的效果.

Q矩陣代表的是系統(tǒng)噪聲協(xié)方差.該超參數(shù)被用來(lái)表示狀態(tài)轉(zhuǎn)換矩陣與實(shí)際過(guò)程之間的誤差.Q越小代表整個(gè)系統(tǒng)更加相信模型的結(jié)果,Q越大則表示整個(gè)系統(tǒng)更加相信觀測(cè)數(shù)據(jù).因此,選擇合適的Q矩陣對(duì)實(shí)驗(yàn)結(jié)果的影響較大.由于難以獲得噪聲協(xié)方差矩陣Q的良好估計(jì),本工作僅通過(guò)實(shí)驗(yàn)不同Q值對(duì)結(jié)果的影響找到了一個(gè)較為合適的Q值,但如何從數(shù)據(jù)中估計(jì)協(xié)方差,可作進(jìn)一步研究.例如,采用自協(xié)方差最小二乘[29-30]、場(chǎng)卡爾曼濾波器[31]等來(lái)自動(dòng)估計(jì)系統(tǒng)噪聲協(xié)方差.