基于自適應平方根濾波的MHD/MEMS信號融合方法*

蔣 愷，李醒飛，劉 帆，拓衛曉

(天津大學精密測試技術及儀器國家重點實驗室，天津 300072)

航天器在軌運行期間，極易受機械振動和空間環境等多種因素的影響，誘發出低幅值、寬頻帶的微角振動[1-2]。這些微角振動源覆蓋頻帶達千Hz，幅值小，難以測量，處于0.01 μrad到幾十μrad量級。隨著高精度航天器和機載系統探測距離逐步提升，微角振動成為限制有效載荷指向精度進一步提高的主要因素之一[3]。因此，實現在軌微角振動特征的測量具有巨大意義。

傳統陀螺儀的測量精度很高，但工作帶寬一般小于20 Hz，微角振動的寬頻帶測量需求(0.1 kHz～1 kHz)會增加陀螺儀的體積或降低陀螺儀的精度[4]。磁流體動力學(Magnetohydrodynamics，MHD)角速度傳感器具有寬頻帶、低噪聲、體積小等優良特性，極其適合測量高頻角振動[2]。但MHD角速度傳感器存在原理極限，低頻時受反電動勢和粘滯力的影響較大，輸出性能不佳，單獨使用時難以滿足微角振動測量的需求。

盡管可以通過進一步處理MHD角速度傳感器的輸出[5]或增加MHD泵等[6]手段改善其非工作頻段的測量能力，拓寬測量頻帶，但要兼顧極寬的測量頻帶和保持原本的靈敏度仍較困難。為了實現這一目標，國內外諸多學者致力于全頻帶慣性傳感器數據融合技術，來彌補MHD角速度傳感器低頻測量的盲區。周新力[7]利用混疊濾波器實現了MHD角速度傳感器和MEMS陀螺儀的數據融合，實現了0～40 Hz頻段內相位波動±1°、幅值波動小于1%的融合效果。Marcelo[8]提出了閉環控制濾波方法，利用高帶寬角位移傳感器(Angular Displacement Sensor，ADS)和低帶寬慣性參考單元(Inertial Reference Unit，IRU)實現了0～100 Hz帶寬內姿態角的測量。毛耀[9]提出了閉環控制濾波的最優校正方法，利用MEMS陀螺儀和MEMS加速度計數據融合實現了1 Hz～100 Hz帶寬內最大幅值誤差2.1 dB，最大相位誤差17°的融合效果，有效提升了MEMS陀螺儀的測量帶寬。拓衛曉[10]利用閉環控制濾波實現了MEMS陀螺儀和MHD角速度傳感器的數據融合，實現了0.1 Hz～60 Hz頻帶內最大幅值畸變小于0.9 dB，最大相位誤差小于34°的融合效果。上述研究均是從能夠直觀評估融合后數據的相位、幅值以及延遲等特性的頻域角度開展的，無法兼顧傳感器的漂移和噪聲等指標。基于時域的融合方法能夠改善傳感器融合中的漂移、噪聲等指標，在傳感器融合領域也得到了學者的廣泛關注。Iwata[11]根據應用維納濾波實現了不同傳感器的信號融合，但維納濾波要求輸入信號統計特性已知，且濾波過程中不能自適應調節。金釗、何世民等人[12-13]建立了MHD角速度傳感器和星敏感器的卡爾曼濾波框架，并對融合算法進行了仿真驗證。紀越[14]對不同頻段的數據采用不同的卡爾曼增益計算式，實現MHD角速度傳感器和MEMS陀螺儀的融合。相比于MHD角速度傳感器，融合輸出在低頻段的信噪比提升了近10 dB。但該算法需要計算多個卡爾曼增益，并需要識別輸入信號的大致范圍，采用對應的卡爾曼增益，對于實現算法的MCU提出了較高的要求。韓佳辰[15]應用Allan方差實時計算測量協方差矩陣，實現了MHD角速度傳感器和MEMS陀螺儀的自適應融合。相較于MHD角速度傳感器，融合后信噪比提高了近10 dB，均方根誤差降低了30%。但上述研究均停留在算法和仿真層面，尚未對具體的硬件實現過程進行研究。具體實現時，主控芯片資源有限，部分算法難以實現，有些則因為計算精度的限制，在實際濾波過程中與計算機仿真結果不一致，導致協方差矩陣易失去正定性，導致濾波發散。

針對上述問題，本文基于Joseph形式的協方差更新計算式，提出了一種平方根卡爾曼濾波算法，實現了MHD角速度傳感器和MEMS陀螺儀的信號融合。算法在MATLAB平臺中通過仿真驗證有效性后，移植到FPGA中，并搭配外圍電路進行了實驗測試，實驗結果表明提出的算法在實現MHD角速度傳感器與MEMS陀螺儀數據融合方面具有良好結果。

1 傳感器復頻域模型的建立

1.1 MHD角速度傳感器頻域模型

MHD角速度傳感器表頭部分由敏感元件與信號檢測電路兩部分構成。圖1為敏感元件原理示意圖，上層圓環為流體通道，內部充滿慣性較大的導電流體；下層圓環為永磁體，與傳感器外殼固連。當載體發生轉動時，導電流體由于慣性作用與慣性空間保持相對靜止，從而與永磁體之間產生相對流速，導電流體切割磁感線，在內外壁間產生與輸入角速度成正比的動生電動勢E。根據非慣性系下的納維斯托克斯方程，建立如下的傳感器敏感元件部分傳遞函數

圖1 MHD角速度傳感器工作原理示意圖

式中:B為磁場強度，ν為導電流體的運動粘度，M為哈特曼數，r為流體環均方根半徑，W為流體環寬度，h為流體環的高度。

敏感元件輸出的弱電壓信號經過調理電路后，作為表頭的輸出，因此MHD角速度傳感器整機的頻響特性應為敏感元件部分傳遞函數與調理電路傳遞函數的級聯。所用的MHD角速度傳感器為實驗室自研的MHD03型號角速度傳感器，帶寬為3 Hz～1 000 Hz，復頻域模型可表示為

式中:所有系數均為常數。

1.2 MEMS陀螺儀的頻域模型

低頻性能良好的陀螺儀種類繁多，如光纖陀螺儀、動力調諧陀螺儀等，但大部分都存在體積大的缺陷，難以在安裝體積要求較高的場合取得應用。MEMS陀螺儀具有體積小、性能可靠的優點，較為適合用于補償MHD角速度傳感器低頻性能。本文選用的MEMS陀螺儀為日本SILICON SENSING公司生產的CRS03型陀螺儀(帶寬0～20 Hz)，其復頻域模型可表示為

式中:所有系數均為常數。

圖2所示為MHD角速度傳感器和MEMS陀螺儀復頻域模型對應的伯德圖，兩種傳感器的幅度曲線在0～1 000 Hz的頻域帶內互補，具備實現組合測量的前提。

圖2 MHD角速度傳感器和MEMS陀螺儀伯德圖

2 基于Allan方差的平方根濾波算法

2.1 基于Allan方差的卡爾曼濾波模型

實現兩個具有不同頻帶特性傳感器融合的卡爾曼濾波模型如圖3所示。圖中，ωtrue為濾波系統真實輸入；ωL，ωH分別為MEMS陀螺儀和MHD角速度傳感器的量測值；AL，BL，AH，BH分別為MEMS陀螺儀和MHD角速度傳感器的系數矩陣；x1為由濾波系統得到的真實角速度；XL＝[x2x3]T，ηL，νL分別為MEMS陀螺儀的狀態變量、系統噪聲和測量噪聲；XH＝[x4x5x6x7]T，ηH，νH分別為MHD角速度傳感器的狀態變量、系統噪聲和測量噪聲。

圖3 角速度測量模型

根據圖3建立以下自適應卡爾曼濾波方程組。

①狀態更新方程

式中:Φk，k-1＝I＋AT，T為所選擇的AD的采集周期。

②時間更新方程

估計測量噪聲:

更新狀態值和系統過程噪聲方差矩陣:

由上式可知，算法在通過測量數據進行濾波遞推時，實時估計系統的測量噪聲，從而達到自適應濾波的目的。由于整個系統階數較高，且系數矩陣Φ陣病態，導致協方差矩陣P陣中數據量級相差較大，由于舍入誤差的影響，協方差矩陣P陣不再保持原本的正定性，致使濾波發散。

2.2 基于Allan方差的平方根濾波模型的建立

根據矩陣理論，當矩陣正定時，可以分解為兩個三角矩陣的乘積，且這兩個矩陣互為轉置。在卡爾曼濾波中，協方差矩陣P滿足該條件，故有

由2.1可知，協方差矩陣P在傳遞過程中失去正定性，從而導致濾波發散，而利用其平方根Δ進行傳遞時，則可以保證協方差矩陣的正定性[16]。

基于序貫濾波的平方根濾波流程串行計算較多，難以矩陣化并行實現，結合矩陣平方根的概念和Joseph形式的協方差更新計算式，提出了平方根濾波模型，并利用Allan方差計算式更新量測噪聲，實現自適應濾波。具體計算過程如下:

①狀態更新方程

式中:第二式表示通過QR分解得到Δk/k-1。

②時間更新方程

測量噪聲計算:

更新狀態值和系統過程噪聲方差矩陣:

式中:第三式與式(8)中第二式意義相同，表示式子左邊的矩陣經過QR分解得到式子右邊的矩陣。

經過上述處理，使用協方差矩陣的平方根替代協方差矩陣完成卡爾曼濾波的迭代計算，保證了濾波過程的穩定進行。

3 硬件系統設計與仿真驗證

3.1 硬件系統設計

算法需搭配相應的硬件方能正常工作，設計的硬件結構如圖4所示，電源系統負責整個電路以及傳感器的供電。在傳感器輸出信號采集之前，由于MHD角速度傳感器的標度因數約為MEMS陀螺儀的一半，需要對兩種傳感器的輸出做一定的調理，具體方案為MHD角速度傳感器的輸出信號進行2倍放大處理，MEMS陀螺儀經過同型號的單倍放大芯片。AD采集后將信號傳輸至FPGA中，在FPGA將高低頻信號融合后傳輸給上位機，在上位機中通過Labview完成進一步分析。

圖4 電路結構設計

圖4中供電模塊負責滿足各個模塊的供電需求，根據各模塊的需求設計得到電源模塊如圖5所示。外部輸入12 V電源，分為兩路，一路經過DC-DC升壓并轉換后得到±15 V電壓，經過LDO降壓得到±12 V電壓，由于經過LDO降壓，此時電源已較為干凈，可用于給傳感器供電；另一路經過DC-DC降壓得到9 V電壓，在經過多個型號的LDO降壓分別得到5VA、5VB、5VC、5VD等多種5 V電壓，5VA和5VB為5 V模擬電壓，其對于電源質量要求較高，其中5VB為傳感器的供電電源，其輸出質量與電源相關，需要采用獨立的LDO轉換；5VC和5VD為數字電壓，5VC為FPGA板供電電源，5VD為采集板供電電源。

圖5 電源結構設計

3.2 FPGA邏輯設計

FPGA中執行的邏輯如圖6所示，當檢測到DRDY端口低電平時，FPGA通過SPI接口接收兩路AD同步采集的MHD角速度傳感器和MEMS陀螺儀的輸出，經過定點轉浮點、平方根濾波模塊后得到融合輸出，將數據轉換為定點數后通過約定好的協議進行編碼發送，在上位機中進行解碼得到融合數據。

圖6 FPGA邏輯設計

考慮協方差矩陣更新需要時間較長，因此在將更新好的數據傳輸出去后再進行計算，以圖減小數據在FPGA中暫存時間，從而縮減算法實現造成的相位延遲。

3.3 融合算法的仿真驗證

本文基于Modelsim平臺對融合算法在FPGA中的執行效果進行仿真驗證，選取0.5 Hz、5 Hz、50 Hz三個頻率點分別進行驗證。仿真得到的融合輸出與兩個傳感器原始數據的對比如圖7～圖9所示，融合輸出在三種頻段內都能很好地跟隨輸入信號。

圖7 0.5 Hz處融合效果

圖8 5 Hz處融合效果

圖9 50 Hz處融合效果

4 實驗驗證

4.1 實驗驗證平臺

實驗驗證平臺如圖10所示，將兩個傳感器通過工裝固連，安裝在轉臺的正中心，通過工裝的平面度和垂直度確保兩個傳感器的敏感軸和轉臺的旋轉軸共線。

圖10 融合算法驗證實驗平臺

測試時，控制轉臺以Ω＝sin(2πft)rad/s的角速度轉動，其中f為指定的測試頻率。安裝在轉臺上的傳感器感應到轉臺轉動，輸出與角速度相關的模擬電壓值，雙通道采集模塊采集兩路傳感器的輸出，通過SPI協議將量化的電壓值發送給FPGA，在片內對兩路傳感器數據進行濾波融合，最后將融合數據通過串口輸出到上位機，在上位機采用FFT計算該頻率下的幅度和相位。

4.2 靜態噪聲功率譜測試

本文采用經典功率譜評價融合輸出的靜態噪聲。測試靜態噪聲時，保持轉臺輸入為零。FPGA串口輸出端僅可輸出一路信號，在測試時控制FPGA依次輸出MHD角速度傳感器輸出、MEMS陀螺儀輸出和融合輸出，通過上位機的程序進行接收。在MATLAB中利用平均加窗周期法(Pwelch函數)對接收到的數據進行處理，得到輸出噪聲的功率譜密度曲線如圖11所示，融合信號的靜態噪聲特性與MHD角速度傳感器持平。

圖11 噪聲功率譜密度

4.3 標度因數測試

幅頻響應特性表征著輸出與輸入的比值隨頻率的變化規律，是傳感器的重要特性，而標度因數正是這一特性的表現。為測試融合輸出標度因數隨頻率變化的情況，控制轉臺以穩定幅度、不同頻率轉動，進行掃頻實驗。具體的掃頻頻率選擇為0.1 Hz～0.9 Hz內，掃頻頻率間隔為0.1 Hz；1 Hz～20 Hz內，掃頻頻率間隔為1 Hz；20 Hz～60 Hz內，掃頻頻率間隔為5 Hz。分別對輸入角速度幅值為20 °/s、30 °/s、40 °/s、50 °/s進行測試，對于不同幅值相同頻率點的標度因數取算術平均值，得到的融合輸出標度因數如表1所示。

表1 標度因數測試

根據表1可知，在測試的頻帶內，經過算法融合的標度因數相對于MHD角速度傳感器和MEMS陀螺儀的標度因數，明顯更為穩定。融合后的標度因數頻帶內最大波動在2.5%以內，低頻段穩定在0.04±0.000 5 V/(°/s)以內，表明算法有效補償了MHD角速度傳感器的低頻性能。

4.4 線性度測試

線性度是傳感器的重要特性之一，為測試系統標度因數隨幅度變化的情況，控制轉臺以不同幅度、相同頻率轉動，進行測試實驗。1 °/s～50 °/s范圍內，幅度間隔為5 °/s。對于每一個頻率點都進行以上幅度點的掃幅測試，頻率點的選取為0.2 Hz、2 Hz、20 Hz、50 Hz，以上頻率點涵蓋了低頻段、中頻段以及高頻段，能夠較好地評價系統的線性情況。取測量數據與擬合直線之間的最大誤差與50 °/s時測量數據的百分比作為評價基準，由于MHD角速度傳感器在0.2 Hz處基本感應不到信號，故在0.2 Hz處不對MHD傳感器的線性度進行評估，計算得到的其余各頻率點線性度如表2所示。從表2可知，融合后的信號保持著MHD角速度傳感器優良的線性度，滿足測量的需求。

表2 線性度測試

5 結論

MHD角速度傳感器具有低噪聲、寬頻帶等優點，但由于反電動勢和粘滯力的影響，具有原理極限，難以滿足微角振動DC-1 kHz的測量需求。MEMS陀螺儀具有體積小、低頻性能良好的優點，在頻帶上與MHD角速度傳感器呈互補特性。本文基于Joseph形式的卡爾曼濾波，提出了一種平方根濾波算法，利用Allan方差實時更新測量協方差矩陣，實現了MHD角速度傳感器和MEMS陀螺儀數據的實時融合。線性度和靜態噪聲測試結果表明融合輸出的噪聲水平和線性度與MHD角速度傳感器相當。掃頻實驗結果表明算法能夠有效補償MHD角速度傳感器的低頻測量盲區。