卡爾曼濾波算法在MEMS 熱式風速風向傳感器中的應用?

趙 天田 潤易真翔秦 明

(東南大學電子科學與工程學院，江蘇 南京210097)

風速風向是氣象測量中重要參數之一，能夠直觀地表征空氣流動及氣壓分布的情況，因此，風速風向的測量在氣象預測、航天、軍事等領域均有重要意義[1]。 傳統的風速風向傳感器主要包括風杯式和超聲式，但是此類傳感器具有維護成本高、小型化困難的缺點[2]。 MEMS 技術的發展為解決上述問題提供了新的技術方案。 MEMS 熱式風速風向傳感器主要由加熱元件和溫度敏感元件組成，基于強迫對流原理，通過監測芯片表面氣體流動引起的熱場梯度變化實現風速和風向測量[3-5]。 1974 年，荷蘭代爾夫特理工大學的Van Putten 和Middelhoek 首次采用標準的硅加工技術成功研制了熱流量傳感器[6]。自此之后，世界各大高校紛紛開始了MEMS 風速風向傳感器的研究[7-11]，通過近半個世紀的理論和技術積累，該傳感器已逐漸可應用于實際生產生活中[12-16]。 相較于傳統的風速風向儀，MEMS 熱式風速風向傳感器具有體積小、功耗低、精度高、便于集成和可批量化制備的優點[17]。

但是，在實際的長期應用中，本征噪聲的存在容易造成MEMS 風速風向傳感器輸出不穩定等性能退化的現象。 由于本征噪聲源自環境或傳感器本身存在的分子的隨機運動，屬于白噪聲。 這部分噪聲是傳感器本征固有的，無法通過電路設計解決這一問題[18]。 在此之前，實驗室研究的MEMS 風速風向傳感器系統均采用加權滑動均值濾波算法以減小本征白噪聲[18，19]，但該算法存在濾波效果差、內存消耗多等問題。 因此，為了解決上述問題，本文首次將Kalman 濾波算法引入到MEMS 風速風向傳感器中，并通過仿真及實驗證明通過該數字濾波算法可以有效提高響應速度，降低測量誤差。

1 Kalman 濾波原理及算法

1.1 算法簡介

Kalman R E 在20 世紀60 年代初提出了Kalman濾波算法[20]。 Kalman 濾波算法是一種時域算法，其將信號處理的過程視為一個線性系統，將本征噪聲視為系統的一個激勵，用狀態方程來描述輸入和輸出的關系，以此形成濾波算法。 Kalman 濾波只涉及時域變量，可以對平穩的一維隨機過程如風速傳感器的輸出進行估計。 此外，由于算法的遞推性，可以在計算機上實現即時處理。 自卡爾曼濾波算法被發明以來，其已在各類工程實踐如慣性導航，制導系統，全球定位系統，通信與信號處理，航空航天等領域被廣泛應用[21]。

Kalman 濾波算法的核心原理是綜合考慮上一時刻的狀態最優估計Xk-1以及當前時刻的測量結果Zk，并根據系統噪聲和測量噪聲的統計特點對當前時刻的狀態的最優估計Xk進行預測，具體的推導過程利用到了線性系統原理以及射影定理等，較為復雜，在此不再贅述[21]。 基于以上理論，可以給出以下狀態空間模型，如式(1)所示:

在該模型中，Xk和Xk-1均是n×1 維的狀態向量；φk/k-1是n×n維的狀態轉移矩陣；Γk/k-1是n×l維的系統噪聲分配矩陣；Zk是m×1 維的量測向量；Hk是m×n維的量測矩陣；此外Wk-1和Vk分別是l×1維的系統噪聲向量和m×1 維的測量噪聲向量，且這兩個向量均為期望為零的白噪聲序列，彼此互不相關，即滿足式(2):

Kalman 濾波算法主要分為以下五步，如式(3)～式(7)所示:

①計算狀態一步預測:

②計算狀態一步預測方差:

③計算濾波增益:

④計算狀態估計:

⑤計算新狀態估計均方誤差:

通過以上五步循環迭代即可實現濾波。 具體的計算過程可以用圖1 來表示。

圖1 Kalman 濾波迭代框圖[22]

可以看到對于這個系統來說，量測信息Zk是主要的輸入，并且對于非時變系統，系統參數包括φk/k-1，Γk/k-1，Hk以及噪聲參數Qk-1，Rk均為定值，這樣在實際計算中只需要存儲Xk，Pk兩個變量，輸入也僅有Zk，對于單片機來說內存消耗小，計算較簡單。 以下將介紹本系統中的具體參數設置。

1.2 參數設置

對于本文實驗中所涉及的風速風向傳感器而言，其系統狀態可以設置為

其中Xk表示系統狀態，Tk表示當前時刻測溫元件處的溫度，Ak表示測溫元件處溫度的變化速率，因此對于下一時刻的狀態預測，可以認為狀態保持不變，即:

因此，可以得到狀態轉移矩陣:

其中Δt為采樣間隔或采樣周期。

本文中涉及的熱式風速傳感器是通過熱溫差的方式實現風速測量的，具體原理為，當有風流動時會使傳感器表面的熱場發生改變，導致傳感器上下游出現溫差ΔT，通過一對熱電堆可以對該溫差進行測量，測量是通過Seebeck 效應實現的，Seebeck 效應又稱第一熱電效應，是指半導體兩端的溫度差異而引起電壓差的熱電現象，對于同一種材料有簡單的表達式如式(11)所示[23]:

根據以上原理，給出式(12)所示測量方法:

即有量測矩陣:

以上兩式中α均為Seebeck 系數，此外在實際應用中根據測量原理不同，這里的溫度T可以指當前溫度，也可以指當前溫度差或是其他適用于風速測量的物理參數。

考慮到利用Seebeck 效應的計算其實是一個線性過程，因此完全可以直接將熱電堆的輸出電壓Vk和熱電堆輸出電壓的變化速度Dk視為系統狀態的組成，這樣Seebeck 系數α就不會出現在表達式中，可以顯著簡化運算，如式(14)所示:

至于噪聲參數，測量噪聲V可以通過對傳感器噪聲測量得到，對于確定的傳感器，其具有與傳感器物理結構相關的本征噪聲。 為了測量該部分噪聲的大小，可將傳感器置于風洞中，不經過數字濾波直接測量給定風速下傳感器輸出電壓的變化范圍，將其峰峰值作為測量噪聲V的取值。

對于系統噪聲W，在卡爾曼濾波算法中，該參數與狀態一步預測結果的可信度有關，如果該值較大，則濾波結果將更接近于測量值，如果該值較小，則濾波結果將更趨近于預測值。 在實際使用中可以通過調整該參數大小以取得更好的濾波效果。

對于采樣間隔Δt這個參數值的選擇有幾個需要考慮的地方。 首先采樣過程相當于一個低通濾波器，采樣間隔Δt越大，濾波器的截止頻率越低。 如果增大采樣間隔，很多噪聲會在采樣過程中被消除，但是如果采樣間隔過大，就會導致信息的丟失。 其次，選取采樣間隔時還應考慮單片機本身的性能，采樣間隔至少要大于一次Kalman 循環所需要的計算時間，所以Δt必須根據實際需求進行權衡。

2 Kalman 濾波算法軟件仿真

本節通過軟件仿真驗證Kalman 濾波算法的有效性，所選仿真軟件為MATLAB 2016a。

2.1 濾波效果仿真

首先，利用示波器對8.0 m/s 風速下傳感器的輸出信號進行采集，去除有效信號提取出噪聲信號，分別采用加權滑動均值濾波算法和Kalman 濾波算法，對比兩種算法的濾波效果。 本次實驗中共采集了50 s 的數據，采樣間隔為0.000 8 s，因此共有62 500 個數據。 所得結果如圖2 所示。

圖2 Kalman 濾波及加權滑動均值濾波兩算法濾波效果仿真結果對比

如圖2 所示，原始噪聲數據的峰峰值高達9.97 mV，使用加權滑動均值濾波后的峰峰值減小至1.39 mV，而使用Kalman 濾波后的峰峰值為4.72 mV。可以看到兩種濾波算法都有不錯的濾波效果，其中加權滑動均值濾波效果要更好一些，但是實際應用中，不能只看濾波效果，還應考慮算法的響應速度。

2.2 響應速度仿真

為了比較兩種算法的響應速度，在25 s 時刻添加了一個大小為1 V 的階躍信號，之后分別進行加權滑動均值濾波和Kalman 濾波，觀察其響應速度，其結果如圖3 所示。

圖3 Kalman 濾波算法及加權滑動均值濾波算法響應速度仿真結果對比

此處選取首次到達目標電壓的時間作為響應時間。 可以看到，加權滑動均值濾波的響應時間為7.5 s，而Kalman 濾波的響應時間僅為1.64 s，Kalman 濾波的響應速度明顯占優。 實際上，根據國家標準，有效的風向風速測量系統需要3 s 內回傳一組數據，即響應時間應不小于3 s。 因此，按照上述仿真的參數設置，加權滑動均值濾波的響應速度并未達到國家標準的要求，故無法實際使用。 與之形成對比，Kalman 濾波雖然濾波效果稍差，但響應速度上仍有很多余量，考慮到在數字濾波算法中，響應速度和濾波效果是一對相互矛盾的變量，因此在實際使用中可以通過調整參數，犧牲部分響應速度以加強濾波效果。

通過仿真結果可以發現Kalman 濾波雖然濾波效果遜于加權滑動均值濾波，但相比之下響應速度快很多，在實際使用中，需要針對實際指標要求判斷兩種濾波算法的優劣。

3 Kalman 濾波算法實驗結果

3.1 實驗條件

為了驗證實際使用中Kalman 濾波的效果，將Kalman 濾波算法及加權滑動均值濾波算法移植到8位單片機C8051F410 上，在風速8 m/s，濕度25%，環境溫度18 ℃的條件下進行測試。 實驗裝置及測控電路原理如圖4 所示。

實驗過程中，將風速傳感器與支架固定后放入風洞內，當風洞工作時，風速傳感器表面的熱場會根據風速大小發生不同程度改變，風速信號會首先被轉化為熱信號，再通過Seebeck 效應將熱信號轉化為電壓信號，之后通過PCB 上的數據采集電路(包括濾波、放大電路) 將數據傳入單片機并進行Kalman 濾波、數據處理得到風速風向信息，最后通過UART 串口發送至PC 上顯示，此外PCB 板上還有一塊加熱電路使芯片表面電路略高于環境溫度，便于測量溫場變化。

圖5 為實際測試中所使用的MEMS 熱式風速風向傳感器，PCB 電路及測試用風洞的實物照片。

圖4 實驗裝置及測控電路原理示意圖

圖5 MEMS 熱式風速風向傳感器，PCB 電路及測試用風洞的實物照片

在本次實驗中，根據1.2 節中參數設置的規則，參數值設置如下:

系統噪聲方差:

測量噪聲方差:

該值為風速8 m/s 時傳感器低頻噪聲(0.1 Hz～10.0 Hz)的測量值。

測量間隔:

3.2 實驗結果

根據上述參數，進行實驗，得到輸出電壓，結果如圖6 所示。

圖6 Kalman 濾波算法及加權滑動均值濾波算法實際濾波效果對比

首先考慮濾波效果，可以看到原始信號噪聲的峰峰值高達56.2 mV，加權滑動均值濾波后的峰峰值為19.7 mV，而Kalman 濾波后的峰峰值僅為9.7 mV，相較于原始信號降低了82.7%，相較于加權滑動均值濾波的結果降低了50.8%。

為了測試兩種算法的響應速度，先將傳感器啟動后放置于風洞中，設置風速為0，放置15 min 待其輸出穩定后將風速設置為8 m/s，等待5 min 后，收集數據，結果如圖7 所示。

圖7 Kalman 濾波算法及加權滑動均值濾波算法實際響應速度對比

不難看出，加權滑動均值濾波的響應時間為5.76 s，而Kalman 濾波的響應時間僅為2.56 s，相較于加權滑動濾波快了55.6%，滿足3 s 回傳時間的國家標準。 值得注意的是，圖7 中Kalman 濾波結果存在一段過沖，但考慮到過過沖電壓較小，為0.04 V，滿足本項目中指標要求(在該風速段中換算到電壓值為0.048 V)，因此，仍取首次到達目標電壓值的時間為響應時間。

之后為計算濾波算法對測量誤差的影響，考慮到無法直接將電壓值與風洞設置的準確值進行對比，需要先將輸出電壓轉化為風速數據。 由于表征輸出電壓與風速關系的數學表達式較復雜[24]，因此在嵌入式代碼中直接使用線性插值的方法將電壓與風速一一對應起來，這樣就可以對測量誤差進行表征。

圖8(a)和8(b)分別為不同風速下經過濾波、插值后所得風速結果較風洞輸入值的絕對誤差和相對誤差，很明顯，風速測量的絕對誤差隨風速增大而增大，相對誤差則隨風速增大而減小。 此外，應用Kalman 濾波算法后， 最大風速絕對誤差從±0.34 m/s 縮小至±0.25 m/s，減小了27%。 絕對誤差在風速為11.2 m/s 時減小最多，由±0.20 m/s 縮小為±0.08 m/s，減幅高達60%。 相對誤差在風速為0.8 m/s 時減小最多，由2.7%減小為1.42%。

可以看出，在上述參數設置條件下，Kalman 濾波無論在濾波效果、響應速度和誤差上都明顯優于加權滑動均值濾波，而且響應時間尚有一定余量，如有必要仍可進一步強化濾波效果。 因此，Kalman 濾波算法可以有效應用于本文所涉及的MEMS 熱式風速風向傳感器。

圖8 Kalman 濾波算法及加權滑動均值濾波算法誤差對比

4 總結與展望

本論文就Kalman 濾波應用于MEMS 熱式風速風向傳感器的可行性進行了討論。 首先，介紹了Kalman 濾波算法的原理，逐步分析了Kalman 濾波算法的流程，詳細給出了實際應用時的參數設置方法和注意事項。 然后，利用MATLAB 對Kalman 濾波算法進行仿真，并與傳統的加權滑動均值濾波算法進行對比，從仿真結果可以看出，Kalman 濾波算法在濾波效果上稍遜于加權滑動均值濾波，但響應速度有顯著優勢。 最后，將Kalman 濾波算法寫入單片機，直接采集傳感器輸出信號進行濾波，觀察濾波結果并與加權滑動均值濾波的結果進行對比，通過調整參數可以使Kalman 濾波的濾波效果、響應時間以及誤差都優于加權滑動均值濾波，由此確定Kalman 濾波算法可以應用于MEMS 熱式風速風向傳感器，該算法能有效抑制測量中出現的低頻白噪聲，提高傳感器的精確度和穩定性。