一種基于Haar小波的光纖陀螺信號實時濾波方法

曲俊海，夏元清，李 靜，王海穩

（1.北京理工大學 自動化學院，北京 100081；2.北方自動控制技術研究所，太原 030006；3.太原科技大學 電子信息工程學院，太原 030024）

慣性穩定平臺控制系統多以光纖陀螺作為敏感相對慣性空間旋轉角速度的測量元件，利用反方向補償該角速率的變化，實時閉環控制實現穩定及跟蹤功能[1]。但是，光纖陀螺輸出的角速度信號中含有噪聲，該輸出噪聲會直接影響到慣性穩定平臺的控制精度。尤其當處于低速工況時，未經濾波處理的陀螺輸出噪聲分量的標準偏差值甚至高于陀螺的期望輸出值，系統輸出除了跟隨輸入信號外，還要跟隨陀螺的輸出噪聲，控制效果受噪聲影響很大。因此，光纖陀螺的性能是影響系統控制精度的關鍵因素。

針對提升陀螺信號質量的問題，國內外學者提出的軟件濾波方法主要分為三種：

一是利用卡爾曼濾波對陀螺信號進行濾波[3-4]。這類濾波方法在隨機噪聲模型建立準確的情況下，具有較好的濾波效果，然而在實際應用中很難得到準確的模型，因此濾波效果往往受到模型精確性的較大限制。二是平滑濾波，例如滑動平均，最小二乘回歸或者Savitzky-Golay回歸等[5-6]。這類方法不依賴于先驗知識，而且運算量小，適合在線進行，然而此類方法的魯棒性較差，當陀螺信號出現有用的突變信號或者異常的野值時，濾波效果會受到很大影響。

三是變換域濾波，該類方法通常根據有用信號和噪聲的頻率分布特性[7]，設計相應的低通或帶阻濾波器對噪聲進行濾除，此類方法由于實現簡單且效果明顯，曾在工程上被廣泛應用。

然而事實上，光纖陀螺的噪聲分布于各個頻帶，頻域濾波法無法徹底濾除低頻噪聲或表征信號的突變點。小波變換是一種在時-頻方面都具有表征信號局部特征能力的方法，利用正交小波的多分辨率特性對信號進行處理，相當于多個帶寬不同的濾波器對信號進行處理，根據信號和噪聲的不同特征，將其區分開來。

目前已有文獻[8-10]利用小波方法對陀螺信號進行濾波，多為基于Daubechies小波基的閾值小波去噪方法。這類小波去噪方法又大可分為兩類：一類是事后濾波方法，即將所有陀螺數據一次性作為輸入，整體得到濾波后的值，此類方法不能滿足實時性應用需求；另一類是滑動窗濾波方法，設置一定長度的滑動窗，以最新一段陀螺數據為輸入，或者對滑動窗截取的最新一段陀螺數據進行一定方法的擴充，讓當前時刻信號采樣值位于待濾波數據的非邊緣段，以此作為濾波輸入數據，得到濾波數據后取對應時刻的輸出，此類方法確可滿足實時濾波需求，但精度較差。因此，提出了對光纖陀螺輸出信號進行小波實時濾波的技術需求。

本文在進行大量仿真實驗的基礎上，說明了影響現有小波方法實時濾波效果的根源，即邊界問題；進而從理論出發，對邊界問題產生的原因進行深入剖析，得出其源于對支撐集的要求，在指出Harr小波基的支撐集符合要求后，進一步證明了Harr小波基可實現連續階梯信號的逼近；最后提出一種基于Haar小波的實時光纖陀螺信號濾波方法。該方法既利用了小波優越的信號-噪聲分離能力，又可滿足實時應用背景的要求，有效地實現了光纖陀螺信號噪聲的實時濾除。

1 基于小波理論的去噪方法

小波閾值方法由 Donoho提出，該方法認為信號對應的小波系數包含信號的重要信息，其幅值較大，但數目較少，而噪聲對應的小波系數是一致分布的，個數較多，但幅值小。閾值法由于實現簡單、計算量小，在工程中得到了廣泛的應用。

利用小波閾值法去噪一般分為3個步驟：

① 對信號進行小波分解，得到各尺度下的小波系數，以三層分解為例，小波分解結構如圖1所示，cn表示第n層的近似小波系數（也稱低頻系數），dn表示第n層的細節小波系數（也稱為高頻系數）；

② 根據噪聲能量及分布為每個尺度的小波系數選擇合適的閾值，對細節小波系數進行閾值操作得到新的小波系數組合；

③ 由新的小波系數進行重構得到去噪后的信號。

圖1 小波分解示意圖Fig.1 Diagram of wavelet decomposition

2 基于Daubechies基的小波濾波

2.1 邊界問題

小波分解和重構是小波去噪的關鍵環節，工程中多使用 Daubechies小波基（簡稱為 DbN小波基，N為小波的階數）來實現小波分解與重構。

實際工程中的小波濾波過程的實現都依賴于Mallat塔式多分辨分解與重構算法，該算法實現一次小波分解的過程如圖2所示。圖中：()s k為原信號；在基于 Daubechies基的小波濾波中，()h k、 ()g k是Daubechies小波基的分解低通濾波器和分解高通濾波器，原信號通過與分解低通和分解高通卷積得到低頻系數和高頻系數； ()c k和 ()d k分別為分解后的低頻系數和高頻系數；()c k′和()d k′為經過下采樣后的低頻和高頻系數。

圖2 Mallat算法一次小波分解實現過程Fig.2 Process of one layer wavelet decomposition in Mallat

Mallat算法實現一次小波重構的過程如圖3所示。

圖3 Mallat算法一次小波重構實現過程Fig.3 Process of one layer wavelet reconstruction in Mallat

Mallat算法是在假定原始信號無限長的基礎上進行的，但實際工程中，我們通常都是對有限長信號進行濾波。

設s(k)為原始信號的L個離散值（1 ≤k≤L），h(k)、g(k)分別為分解低通濾波器和分解高通濾波器，其長度為M。原始信號分別與分解低通濾波器、分解高通濾波器進行卷積得到近似小波系數c(k)和細節小波系數d(k)，見式(1)和式(2)。

可以看到，當計算c(1),…,c(M-1)和d(1),…,d(M- 1)時，需要原信號s(k)初始時刻之前的值參與運算，即左端越界；當計算c(L+ 1),…,c(L+M-1)和d(L+ 1),…,d(L+M-1)時，需要原信號s(k)在L時刻之后的值參與運算，即右端越界。

重構過程卷積運算表達式參見式(3)，可以看到，重構過程也面臨數據越界現象。

2.2 邊界問題對實時濾波性能影響分析

從上述小波分解與重構的過程中可以看出，基于Daubechies基的小波濾波實現過程中，存在數據越界現象，因此實現步驟中包含邊界延拓。

延拓方法主要有零延拓、周期延拓、對稱延拓三種，其具體實現方式為：

① 零延拓

原信號：s(1),s(2),…,s(k)

延拓后信號：0,…,0,s(1),s(2),…,s(k),0,…,0

② 對稱延拓

原信號：s(1),s(2),…,s(k)

延拓后信號：

…,s(3),s(2),s(1),s(1),s(2),

…,s(k),s(k),s(k- 1),s(k-2),…

③ 周期延拓

原信號：s(1),s(2),…,s(k)

延拓后信號：

…,s(k- 2),s(k-1),s(k),s(1),s(2),

…,s(k),s(1),s(2),s(3),…

為考察邊界延拓方法對數據越界問題的效果，進行了仿真實驗。仿真對象采用信號s=sin(0.2πt)，原信號取 10個采樣點，采樣序列為{0.31038,0.00159,-0.30735,-0.58624,-0.80779,-0.95037,-0.99999,-0.951 84,-0.810 61,-0.590 10}。根據已有試驗經驗，在對陀螺信號濾波時，Db3小波基往往優于其他Dbn（n≠1）小波基，故本文中所有試驗均選用Db3小波。由于信號分解出的低頻系數決定了信號的基本輪廓，所以只對采用幾種延拓方法所得到的一級分解低頻系數進行分析，并與使用真實正弦信號延拓后所得到的低頻系數進行對比。仿真結果如表1所示。

表1 低頻系數對比結果Tab.1 Comparison on low frequency coefficients

從表1中可以看出，與使用真實信號延拓相比，其他無論何種延拓方式，中間的低頻系數（系數 3、系數4、系數5）相同，兩邊的低頻系數（系數1、系數2、系數6、系數7）都出現了偏差。

為更直觀地說明邊界問題的影響，模擬生成10 (°)/s的速度信號，并加入高斯白噪聲作為陀螺仿真原信號，在所有原信號生成后，將其作為輸入進行一次性整體濾波（本文稱為離線濾波），使用Db3對輸入信號進行閾值去噪，延拓方法為零延拓，濾波效果如圖4所示。

圖4 基于Db3小波的離線濾波結果Fig.4 Offline filtered result based on Db3

從圖4中可以看到，邊界問題對原信號左邊界和右邊界的濾波結果都產生了較大的影響，使得濾波值的準確度明顯變低。而在實時小波濾波時，我們獲取的當前時刻信號正位于待濾波信號段的最右端，其濾波結果的準確度取決于右邊界濾波結果的準確度，因此，在使用Daubechies小波基時，克服邊界問題對實時光纖陀螺濾波至關重要。

3 基于Haar小波的實時陀螺信號濾波方法

3.1 Harr小波適用性分析與證明

由第2節可知，當采用Db3小波進行光纖陀螺濾波時，在數據邊界濾波結果的準確度將大為降低，而在實時控制系統中，用于閉環運算的數據正是當前的實時數據（即在邊界的數據），因此邊界問題處理不好會對系統控制產生不利的影響。分析引起此問題的原因是由于塔式分解計算的過程中需要用到采樣點以外的數據，而在實際系統中無法準確預知未來的數據。為克服邊界問題，從小波理論出發，進一步分析可得到引發此問題的原因。由圖5(a) Db3小波尺度函數圖可以看出，Db3小波函數的支撐區間為0 ≤t≤ 5，因此當進行下一層分解時其支撐區間為0 ≤t≤ 2.5，需要6個系數才能運算得到，如圖5(b)所示。

而按照塔式分解的方法，在邊界分解時超過2個小波系數就需要用到原信號的未來數據了。因此帶來了邊界問題。

圖5 (a) Db3小波尺度函數圖Fig.5 (a) Wavelet scale function of Db3

圖5 (b) 第一層分解的小波尺度函數Fig.5 (b) Wavelet scale function when doing first layer decomposition

解決此問題的根本方法是使小波函數的支撐集不大于1，即0 ≤t≤ 1，而按照Db小波的構造辦法可知，只有Db1即Harr小波的支撐集可限制在0 ≤t≤1的范圍內。但 Harr小波缺點是逼近連續信號的能力差，因此還需要證明Harr小波用于光纖陀螺濾波的適用性問題。

本文研究的對象是光纖陀螺信號，實際工程中，都是以一定時間間隔對光纖陀螺信號進行數據錄取。某靜置光纖陀螺輸出信號如圖6(a)所示，對采樣區間[9,10]的陀螺信號進行局部放大，得到圖6(b)。可以看到，待處理的光纖陀螺信號是連續階梯信號。

圖6 (a) 某光纖陀螺輸出信號Fig.6 (a) Original signal of FOG

圖6 (b) 陀螺信號局部放大圖Fig.6 (b) Partial amplification of the original signal

小波分解時，首先要確定近似空間Vj，使其能最佳地反映原信號s的各種信息，然后選擇sj∈Vj，最佳地逼近s。由于 2j/2φ(2jx-k)是Vj空間的標準正交基，因此s在Vj上的投影Pj s可表示為：

其中，

定理：若 {Vj,j∈Z}是一個依緊支撐尺度函數φ的多分辨率分析，如果s∈L（2R）是如圖5(a)所示的連續階梯函數，那么當選擇j，滿足2j=N時（N為信號長度），當φ滿足以下條件：

ⅱ）φ支撐集為[0,1]，可得式(1)中，

對上述定理展開證明，φ是緊支撐的，其非零集被限定在一閉包[0,1]中，因此等式中x的積分區間是0 ≤ 2jx-k≤ 1。進行變量替換t= 2jx-k，可得：

如圖6(b)所示有s(2-jt+ 2-jk) =s(2-jk)，t∈ [ 0,1]，因此可用s(2-jk)代替s(2-jt+ 2-jk)，因此：

由以上分析可知，當采用 Harr小波基時，如果2j=N，取=s(k/2j)，則連續階梯函數在Vj中的投影可最佳逼近原函數s。

通過上述證明可發現，對于連續階梯信號，采用Harr小波基能夠很好地逼近。因此，對于光纖陀螺信號實時濾波的問題，采用Harr小波基既可以完美地逼近原函數，又不存在邊界問題，是最佳的選擇。

3.2 實時濾波算法設計

為實現實時濾波功能，采取滑動窗濾波的方法，在采集到最新時刻陀螺數據后，利用最新時刻及之前時刻的陀螺數據作為待濾波數據。設滑動窗寬度為L，k時刻的光纖陀螺采樣值為x(k)，k時刻的實時濾波結果用表示，實時濾波算法的詳細實現步驟如下：

① 對窗寬度L（2的冪級數）、分解總層級N進行設置。

② 若k＜L，由于數據量較少，不進行濾波，若kL≥，構造濾波器輸入信號序列

③ 使用Harr小波高、低通分解系數對輸入信號進行第n級分解，抽二下采樣后分別得到細節系數和近似系數

其中：若n=1，則待分解信號為輸入信號；若n＞ 1，則待分解信號為n-1層所得近似系數；i的總數取決于每層級待分解信號的長度。

④ 按照一定的閾值處理規則，對多層分解結束后的細節系數進行閾值處理，而近似系數保留。

⑤ 利用Harr小波高、低通重構系數完成對原信號的重構，得到濾波輸出序列其中，即為k時刻的實時濾波結果，即完成一次單點的實時濾波。

⑥ 當采樣時刻k遞增為 +1k時，讀入下一個時刻的陀螺信號值，返回步驟②進行下一次的濾波處理。

從實現步驟也可以發現，Harr小波分解的高、低通分解系數的個數都為 2，抽二下采樣后可完全消除輸入信號外延數據的影響。

4 試驗與結果

為對基于 Haar小波基的陀螺信號濾波方法的性能進行評估，分別利用靜態陀螺數據和動態陀螺數據進行兩組濾波試驗。綜合考慮濾波效率和性能，濾波窗寬度設置為L= 24=16。窗和窗之間有重疊地向后平移，迭代輸出實時濾波結果。

試驗一：靜態陀螺信號濾波

將光纖陀螺放在靜止基座上，令其測量軸位于水平面內并且向東，測得一組靜態陀螺數據作為待濾波原信號。分別利用已有Db3小波方法和本文所提基于Harr小波的濾波方法對陀螺實測數據進行濾波。分解總層級設為3級，使用固定閾值法確定閾值，閾值處理方法為軟閾值法，濾波結果如圖7所示。

從圖7(a)可以看到，兩種方法都能很好地對陀螺的隨機噪聲進行抑制，說明小波方法應用于陀螺去噪是可行的。進一步觀察圖7(b)所示的濾波結果局部放大圖，可以看到，基于Haar小波的陀螺信號濾波結果的隨機噪聲要更小一些。為進一步說明基于Haar小波的濾波方法的優越性，采取3組實測靜態數據作為原信號，以標準差作為評判濾波方法性能的標準，對各信號的隨機誤差進行定量分析，結果如表2所示。

圖7 (a) 靜態陀螺信號濾波結果Fig.7 (a) The filtered results with static FOG signal

圖7 (b) 靜態陀螺信號濾波結果局部放大圖Fig.7 (b) Partial amplification of filtered results with static FOG signal

表2 基于靜態陀螺信號試驗的標準差統計結果Tab.2 Statistical results of STD experimented with static signals (°/s)

由于實測數據是靜態陀螺數據，假設真值為0的情況下，濾波結果標準差越接近零值越精確，由表2數據可得基于Haar小波的信號濾波方法可更好地抑制隨機噪聲。

由3.2節算法步驟已知，不同于離線濾波，本文算法是實時進行的。該算法基于TI公司的TMS320F2812芯片實現，代碼占用空間為4970byte，滿足硬件資源存儲需求；運行時長約為106 μs，遠小于采樣時長1 ms，運算實時性可以保證。為更好地觀察陀螺實時濾波效果，對動態陀螺信號展開濾波試驗。

試驗二：動態陀螺信號濾波

以某慣性穩定平臺控制系統（以下簡稱平臺）為對象展開試驗，該平臺配有三軸光纖陀螺，并安裝于電動六自由度搖擺臺上。搖擺試驗時，平臺工作于速度環，其速度響應（即陀螺輸出）反映了平臺的抗干擾能力，故搖擺試驗中常使用速度均方差來衡量平臺的穩定精度。

本試驗中，在完成對平臺的最優PID控制器設計后，以5°/0.5Hz的幅值/頻率對平臺俯仰向進行搖擺試驗，方位向和橫滾向無擾動，采用原始陀螺數據（不作任何濾波處理）作反饋，對俯仰向陀螺進行采樣，得到平臺速度響應曲線如圖8所示。受系統響應特性的影響，在每個搖擺周期換向點（即擾動加速度最大點）處，平臺會出現較大的穩定誤差，從而使陀螺輸出如圖8所示的頻率約為1 Hz的尖峰信號。

圖8 基于濾波前陀螺信號閉環的速度響應Fig.8 Speed response of platform without using filtered signal in closed-loop calculation

對如圖8所示的原始陀螺數據進行基于本文算法的半實物仿真濾波試驗，為保證通用性，參數設置同靜態陀螺濾波試驗，濾波結果如圖9所示，濾波前后信號頻譜圖如圖10所示。

圖9 動態陀螺信號濾波結果Fig.9 The filtered result with dynamic FOG signal

從圖9所示濾波結果可以得出，對于動態陀螺信號，本文方法仍可以有效地降低隨機噪聲，進一步說明本文方法用于實時陀螺濾波的有效性。圖10所示濾波前后信號頻率成分類似，說明濾波的實施未損失原信號的帶寬。濾波后信號與原信號包絡對比，存在約5 ms時間的滯后，對于速度閉環系統（帶寬約幾十赫茲），此滯后時間完全可以接受。

為定量地評估實時濾波方法對實際平臺綜合性能的影響，直接使用在線實時濾波后信號閉環，重新進行最優PID控制器設計并開展同條件搖擺試驗，新的速度響應信號如圖11所示。

對分別使用濾波算法前后陀螺信號閉環時，所對應的速度響應均方差進行統計，結果如表3所示。試驗結果說明，在對陀螺原始信號進行濾波，并以其作為速度反饋后，由于反饋信號質量的提升，光纖陀螺噪聲對系統帶寬的制約降低，從而可在最優PID控制器的基礎上進一步上調比例系數，最終測試平臺的穩定精度提高約14%。

表3 濾波前后的平臺速度響應均方差Tab.3 Mean square errors of speed response before and after using filtered signal

圖10 濾波前后信號頻譜圖Fig.10 Spectrum diagram of signals before and after being filtered

圖11 基于濾波后信號閉環的速度響應Fig.11 Speed response of platform corresponding using filtered signal in closed-loop calculation

5 結 論

本文首先從Mallat算法步驟的角度出發，對數據越界現象進行描述。然后通過仿真說明，基于Daubechies類小波方法存在嚴重的邊界問題，對于實時陀螺濾波并不適用。進而通過理論分析，詳細證明了Harr小波在本文背景下的適用性相關問題。最后提出一種基于Haar小波的實時陀螺信號濾波方法，并進行了兩組基于靜態和動態數據的濾波試驗，最終驗證本文所提方法對于提升實際慣性穩定平臺控制系統的積極作用。