基于卡爾曼濾波的含噪聲小樣本數據處理方法

劉芬范洪強呂濤李謙錢權

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

在數據采集過程中噪聲不可避免.這些噪聲數據使得進行機器學習、數據挖掘等任務時產生不準確的輸出,因此對噪聲數據進行處理十分必要.傳統的含噪聲數據處理方法主要是通過一系列算法檢測數據當中的噪聲數據并刪除[1].這種直接刪除噪聲數據的方法會損失該條記錄的所有信息,因此更適用于數據量較大的數據集.當遇到小樣本數據集時,由于數據本身很寶貴,因此需要設計新的處理噪聲數據的方法,在充分利用數據信息的同時對噪聲數據進行處理.

本研究提出了一種基于卡爾曼濾波的含噪聲小樣本數據處理方法.首先,結合領域知識,利用數據集背后的物理模型或經驗公式,設計系統轉移矩陣和系統控制矩陣,并建立系統模型;然后,利用得到的系統模型預測數據,并利用卡爾曼增益結合系統模型的預測結果和觀測數據生成修正數據;最后,使用修正后的數據進行數據挖掘或其他數據分析任務.本方法基于卡爾曼濾波對數據進行修正,達到含噪聲數據處理的效果,由于不需要刪除數據,因此較適合小樣本數據集.

1 相關研究

含噪聲數據的處理方法主要有基于分箱的方法、基于模型的方法以及基于機器學習的方法.

1.1 基于分箱的方法

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

除了上述基本分箱法以外,還可以根據任務自定義分箱方法.Kerber[2]提出了一種ChiMerge分箱方法,用χ2統計量將連續屬性值進行離散化.傅濤等[4]結合了基于分箱的Fuzzy C Means算法用于檢測網絡入侵數據.該方法用分箱方式劃分數據集,解決了原Fuzzy C Means聚類算法頻繁更換聚類中心的問題.

1.2 基于模型的方法

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

1.3 基于機器學習的方法

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

現有的噪聲數據處理方法還存在一些問題.基于分箱的方法要求處理的數據具有一定的規模,若數據樣本過少,則數據分箱之后取平均值、眾數或中位數等方式的平滑效果就會大打折扣.這是因為在數據量較少的情況下,很難通過統計的方式得到數據的分布規律.基于模型的方法和基于機器學習的方法本質上是通過一定的手段檢測出噪聲數據,然后將噪聲數據去除以達到降噪的目的.但直接去除整條噪聲記錄會損失大量信息.尤其是面對小樣本噪聲數據時,若直接去除某條噪聲數據會損失很大比例的信息資源.因此,針對小樣本噪聲數據需要設計新的去噪方法.

本工作提出了一種基于卡爾曼濾波的噪聲數據處理方法,結合了模型方法中的回歸以及分箱方法中的平滑技術.卡爾曼濾波處理方法將數據視為系統的狀態,通過數據的經驗公式建立系統的狀態轉移模型和狀態控制模型,以此來預測系統的狀態,得到系統預測數據.最后,通過觀測數據修正系統預測數據,得到最后的修正數據.該修正數據即為降噪之后的數據.卡爾曼濾波處理方法通過修正數據的方式來達到降噪的效果,不用刪除數據,更加適用于小樣本數據.

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

2.1 卡爾曼濾波

卡爾曼濾波是一種線性二次估計方法[16].該方法利用系統的動態模型,通過系統的控制輸入數據以及觀測數據來對系統的狀態量進行估計.該方法相比只使用動態模型獲得的估計值更準確,是一種常用的模型數據和觀測數據融合的算法.卡爾曼濾波的基本思想如圖1所示.

首先,卡爾曼濾波器根據k-1時刻得到的系統狀態的最終估計^xk-1(見圖1中的黑色曲線),產生k時刻系統狀態的預測估計^x-k(見圖1中的綠色曲線).然后，根據k時刻系統的測量估計yk(見圖1中的藍色直線),將預測估計以及測量估計根據相應的“權重”進行融合，產生k時刻系統狀態的最終估計(見圖1中的紅色曲線).這里的“權重”是由預測系統的不確定性和觀測系統的不確定性綜合得出.以上過程在每一個時間步重復，得到每一個時間步系統狀態的最終估計.最終狀態位于預測狀態和測量狀態之間，比任何一個單獨狀態都具有更好的不確定性估計.

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

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

式中:X(k)代表k時刻的系統狀態向量;Z(k)代表觀測向量;u(k)代表輸入向量;v(k)代表m維測量噪聲向量;A(k+1,k)代表從k時刻到k+1時刻的系統狀態轉移矩陣;B(k+1,k)代表從k時刻到k+1時刻的系統控制矩陣;H(k+1)代表k+1時刻的預測輸出轉移矩陣.式(1)被稱為狀態方程,式(2)被稱為測量方程.

卡爾曼濾波算法包括預測和更新兩個階段,其流程如圖2所示.在預測階段,卡爾曼濾波器根據式(3)產生預測模型對當前狀態變量的估計,根據式(4)產生當前系統噪聲的協方差,即

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

當系統觀察到下一個觀測數據(觀測數據會有一定的誤差,如隨機噪聲等)時,進入到更新階段.卡爾曼增益的計算公式為

式中:K表示卡爾曼增益;H表示預測輸出轉移矩陣,即從預測估計到測量估計的計算方式;R代表測量誤差協方差.

根據式(6)更新當前系統狀態量的估計,根據式(7)計算下一時刻需要用到的估計協方差矩陣,即

式中,z代表觀測數據.

卡爾曼算法是迭代的,不斷重復預測和更新步驟.它可以實時運行,僅使用當前輸入測量值和先前計算的狀態及其不確定性矩陣,不需要額外的歷史信息.

2.2 擴展卡爾曼濾波

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

在擴展卡爾曼濾波中,系統的狀態轉移模型和觀測模型如式(8)和(9)所示,二者只要是可微函數即可,并不要求是線性函數.

式中:函數f用于從過去的估計值中計算預測的狀態;函數h用于以預測的狀態計算預測的測量值.f和h不能直接計算協方差,需要計算它們的偏導矩陣(Jacobian矩陣).狀態轉移模型和觀測模型的Jacobian矩陣計算如下:

在解決了f和h不能直接計算協方差的問題后,擴展卡爾曼濾波的迭代方程如式(12)～(16)所示.在預測階段,根據式(12)計算預測模型對當前狀態變量的估計,根據式(13)計算當前系統噪聲的協方差,即當系統觀察到下一個觀測數據時,進入到更新階段.根據式(14)計算卡爾曼增益,根據式(15)修正模型得到的數據,根據式(16)更新系統噪聲的協方差P,為下一步迭代做準備.與卡爾曼濾波一樣,擴展卡爾曼算法也是迭代的,不斷重復預測和更新步驟.

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

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

本工作提出的基于卡爾曼濾波的含噪聲小樣本數據處理方法的流程如圖3所示.首先,對原數據進行處理,用線性模型對數據進行擬合,得到初步的線性模型參數.其次,根據線性模型參數,建立卡爾曼濾波模型中的系統狀態轉移矩陣和系統狀態控制矩陣.再次,根據得到的系統狀態轉移矩陣和系統狀態控制矩陣建立系統模型,結合系統模型得出的結果與觀測數據,對原數據進行修正.最后,根據修正得到的數據進行后續的數據挖掘應用.

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

3 含噪聲材料數據處理實驗

3.1 數據集

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

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

3.2 實驗設置

耐候鋼的大氣腐蝕動力學[25-28]可以表示為

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

對式(17)等號兩邊各取自然對數,推導出式(18).從式(18)可以看出,對耐候鋼的大氣腐蝕增重ΔW和腐蝕時間N分別取對數后,二者呈線性關系,即

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

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

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

由于本實驗數據集中觀測變量即為系統狀態,因此將預測輸出轉移矩陣H設置成值為1的單元素矩陣.考慮到測量會造成噪聲干擾,因此引入觀測噪聲σz,系統的測量方程為

式中,zk為觀測變量,噪聲服從均值為0的正態分布.

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

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

3.3 實驗結果與分析

3.3.1 腐蝕動力學方程建模及卡爾曼濾波效果

本實驗展示了通過卡爾曼濾波處理含噪聲材料數據,來建立腐蝕動力學方程模型的結果,以及數據經過卡爾曼濾波處理后的修正效果.

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

圖4 數據集Data1、Data2和Data3的實驗結果Fig.4 Experimental results of Data1,Data2,and Data3

表1 耐候鋼數據腐蝕動力學方程參數擬合結果Table 1 Parameters estimation results of corrosion kinetics for weathering steel

3.3.2 不同系統噪聲協方差對卡爾曼濾波效果的影響

卡爾曼濾波模型中有一個非常重要的參數Q(系統噪聲協方差),用來表示狀態轉換矩陣與實際過程之間的誤差.Q越小代表整個系統更加相信模型的結果,Q越大則代表整個系統更加相信觀測數據的結果.圖5顯示了不同Q矩陣下數據集Data1的卡爾曼濾波結果.可以看出:當Q越小時整個系統更加偏向于模型的結果,因此橙色線更加接近藍色線;當Q越大時整個系統更加偏向于觀測數據,此時橙色線更加接近綠色線.不同的Q矩陣對數據修正的結果影響較大,因此需要選擇合適的Q矩陣.

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

3.3.3 卡爾曼濾波處理噪聲材料數據后的建模實驗

為了驗證本工作提出的方法對噪聲數據的處理效果,分別測試了未經處理的數據和經過卡爾曼濾波處理的數據的實驗效果.將時間當作輸入,分別用ARIMA模型和RF模型來預測該時間下材料的腐蝕增重.

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

表2 經卡爾曼濾波處理含噪聲數據后的ARIMA回歸和RF回歸結果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的觀測數據具有更高的波動性(見圖4),數據噪聲較大,因此模型預測的R2較小,效果更差;“Kalman-Raw”列中Data1和Data3的R2要大于Data2,Data1的R2由0.863 0增加到0.943 5,Data3的R2由0.851 1增加到0.945 7,這說明經過卡爾曼濾波處理后,Data1和Data3數據中的噪聲被平滑,模型的結果變好.Data2數據經過卡爾曼濾波處理后,其模型的R2反而變小.這是因為這組數據原本的噪聲就比較小,而經卡爾曼濾波處理后反而引入了新噪聲,對原數據產生了影響,模型結果變差.因此,卡爾曼濾波處理對本身就存在較大噪聲的數據才會產生較為積極的影響.在RF模型中,“Kalman-Raw”列中的R2均大于“Raw-Raw”列.這進一步驗證了卡爾曼濾波處理對含噪聲數據會產生較為積極的影響(兩個模型的R2平均提升了6.4%).對比表2中的數據可以看出,卡爾曼濾波處理對時間序列回歸模型的影響較大,這是因為時序模型對數據噪聲更為敏感.

3.3.4 擴展卡爾曼濾波處理噪聲材料數據后的建模實驗

本實驗測試了擴展卡爾曼濾波對于噪聲數據的處理效果.首先按照3.2節中所述方法建模得到了式(17)中F、n的值.假設式(17)的函數關系為f,則可得到

將腐蝕增重ΔW視為系統狀態,在k時刻,擴展卡爾曼濾波把系統按照k-1時刻的方程進行線性化以適應非線性系統.將f(N)按照泰勒公式在Nk-1處展開,則系統模型描述為

式中:xk為系統狀態;N代表時間.系統的測量方程為

式中:zk為觀測變量,噪聲服從均值為0的正態分布.

由式(29)和(30)可知,系統狀態轉移模型和觀測模型分別為非線性模型和線性模型,因此相應的Jacobian矩陣分別為f′(Nk-1)和值為1的單元素矩陣.已知系統的狀態轉移模型、觀測模型以及二者相應的Jacobian矩陣,然后按照式(12)～(16)進行擴展卡爾曼濾波迭代,對原數據進行擴展卡爾曼濾波處理,用于同3.3.3節的下游任務.

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

表3 經擴展卡爾曼濾波處理含噪聲數據后的ARIMA回歸和RF回歸結果Table 3 ARIMA regression and RF regression results of noisy data processed by extended Kalman filter

4 結束語

本工作提出了一種基于卡爾曼濾波的含噪聲小樣本數據處理方法.首先,結合領域知識,用數據集背后的物理模型或經驗公式擬合初步模型,得到系統狀態轉移矩陣和系統控制矩陣.然后,根據得到的系統狀態轉移矩陣和系統狀態控制矩陣建立系統模型.最后,綜合系統模型得出的結果與觀測數據,得到最后的修正數據.

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

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