曲俊海，夏元清，李 靜，王海穩(wěn)

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

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

針對提升陀螺信號質(zhì)量的問題，國內(nèi)外學(xué)者提出的軟件濾波方法主要分為三種：

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

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

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

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

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

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

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

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

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

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

③ 由新的小波系數(shù)進行重構(gòu)得到去噪后的信號。

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

2 基于Daubechies基的小波濾波

2.1 邊界問題

小波分解和重構(gòu)是小波去噪的關(guān)鍵環(huán)節(jié)，工程中多使用 Daubechies小波基（簡稱為 DbN小波基，N為小波的階數(shù)）來實現(xiàn)小波分解與重構(gòu)。

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

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

Mallat算法實現(xiàn)一次小波重構(gòu)的過程如圖3所示。

圖3 Mallat算法一次小波重構(gòu)實現(xiàn)過程Fig.3 Process of one layer wavelet reconstruction in Mallat

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

設(shè)s(k)為原始信號的L個離散值（1 ≤k≤L），h(k)、g(k)分別為分解低通濾波器和分解高通濾波器，其長度為M。原始信號分別與分解低通濾波器、分解高通濾波器進行卷積得到近似小波系數(shù)c(k)和細節(jié)小波系數(shù)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時刻之后的值參與運算，即右端越界。

重構(gòu)過程卷積運算表達式參見式(3)，可以看到，重構(gòu)過程也面臨數(shù)據(jù)越界現(xiàn)象。

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

從上述小波分解與重構(gòu)的過程中可以看出，基于Daubechies基的小波濾波實現(xiàn)過程中，存在數(shù)據(jù)越界現(xiàn)象，因此實現(xiàn)步驟中包含邊界延拓。

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

① 零延拓

原信號：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),…

為考察邊界延拓方法對數(shù)據(jù)越界問題的效果，進行了仿真實驗。仿真對象采用信號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}。根據(jù)已有試驗經(jīng)驗，在對陀螺信號濾波時，Db3小波基往往優(yōu)于其他Dbn（n≠1）小波基，故本文中所有試驗均選用Db3小波。由于信號分解出的低頻系數(shù)決定了信號的基本輪廓，所以只對采用幾種延拓方法所得到的一級分解低頻系數(shù)進行分析，并與使用真實正弦信號延拓后所得到的低頻系數(shù)進行對比。仿真結(jié)果如表1所示。

表1 低頻系數(shù)對比結(jié)果Tab.1 Comparison on low frequency coefficients

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

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

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

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

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

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

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

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

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

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

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

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

圖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}是一個依緊支撐尺度函數(shù)φ的多分辨率分析，如果s∈L（2R）是如圖5(a)所示的連續(xù)階梯函數(shù)，那么當選擇j，滿足2j=N時（N為信號長度），當φ滿足以下條件：

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

對上述定理展開證明，φ是緊支撐的，其非零集被限定在一閉包[0,1]中，因此等式中x的積分區(qū)間是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)，則連續(xù)階梯函數(shù)在Vj中的投影可最佳逼近原函數(shù)s。

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

3.2 實時濾波算法設(shè)計

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

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

② 若k＜L，由于數(shù)據(jù)量較少，不進行濾波，若kL≥，構(gòu)造濾波器輸入信號序列

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

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

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

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

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

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

4 試驗與結(jié)果

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

試驗一：靜態(tài)陀螺信號濾波

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

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

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

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

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

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

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

試驗二：動態(tài)陀螺信號濾波

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

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

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

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

圖9 動態(tài)陀螺信號濾波結(jié)果Fig.9 The filtered result with dynamic FOG signal

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

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

對分別使用濾波算法前后陀螺信號閉環(huán)時，所對應(yīng)的速度響應(yīng)均方差進行統(tǒng)計，結(jié)果如表3所示。試驗結(jié)果說明，在對陀螺原始信號進行濾波，并以其作為速度反饋后，由于反饋信號質(zhì)量的提升，光纖陀螺噪聲對系統(tǒng)帶寬的制約降低，從而可在最優(yōu)PID控制器的基礎(chǔ)上進一步上調(diào)比例系數(shù)，最終測試平臺的穩(wěn)定精度提高約14%。

表3 濾波前后的平臺速度響應(yīng)均方差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 基于濾波后信號閉環(huán)的速度響應(yīng)Fig.11 Speed response of platform corresponding using filtered signal in closed-loop calculation

5 結(jié) 論

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