李晉斐,趙冬青,王棟民,蔡聰聰,賈曉雪,張樂添
(信息工程大學 地理空間信息學院,鄭州 450001)
在多傳感器融合定位中,由于慣性導航系統能夠持續輸出姿態和位置信息,一般多以慣性導航誤差模型建立狀態方程。慣性元件隨機誤差對慣性導航系統產生的影響較大,因此需有效降低慣性元件的隨機誤差,達到提高慣性導航系統測量精度的目的。小波變換具有在時域和頻域同時表征局部特征的能力[1],采用多分辨率特性,使得去除噪聲的過程中保持原始信息,能夠很好地刻畫信號的非平穩特征,且不需要建立精確的隨機誤差模型,即可分析處理非平穩信號,故非常適用于處理陀螺的狀態信息[2]。
所有小波去噪算法都有共同的結構,但其有效性非常依賴相關參數的選擇,主要涉及小波基函數、閾值選取規則、閾值函數、分解尺度等[3]。該過程中,很多學者關注了如何改進閾值函數來提高去噪能力,并取得了一定效果,而對于其他參數的選擇卻十分困難[4],更多的是以經驗或直接給定參數來得到最終結果,很少分析這些參數對最終效果的影響,也沒有明確的小波基函數與分解層數等最優參數的確定方法。例如,對于慣性元件隨機誤差的小波去噪參數,文獻[5]直接選定db4 小波基,3 尺度分解;文獻[6]直接選定db6 小波基,5 尺度分解。然而,根據多次實驗對比發現,采用小波去噪時,去噪效果隨著每個尺度下所選擇的小波變化而變化,很難直接確定分解層數;而分解層數與小波基函數的選取對去噪效果的影響不可忽略。要使去噪效果最好,必須選擇最佳的小波參數,通常選擇信噪比(signal-noise ratio,SNR)、相關系數 R、均方根誤差(root mean square error,RMSE)、平滑度r等參數作為指標來衡量小波去噪效果[7]。在實際應用中,由于信號真值是未知的,上述單一的評價標準通常無法準確判斷去噪效果,具有很大的局限性。因此,陶珂和朱建軍[8]通過融合去噪過程中信號RMSE、SNR 和平滑度的變化量構造具有收斂特性的新統計量,通過識別拐點判斷最佳分解尺度。朱建軍等[9]通過變異系數法進行定權,選擇信號的平滑度與RMSE 2 種指標進行組合,進而指導小波去噪相關參數的選擇。王旭和王昶[10]利用熵權法將RMSE 與平滑度變化量線性組合得到復合評價指標,利用其收斂特性確定最優分解層數。實際計算表明,在真值未知時,這些方法都有較好效果,但在準確度上存在不同程度的缺陷,且在判斷收斂性方面缺乏相應的理論依據。因此,本文構建了一種新的評價指標,并通過對慣性元件隨機誤差的分析,驗證其在小波降噪過程中對相關參數選擇的有效性。
小波閾值去噪算法的實質是抑制信號的無用部分,最大限度保留信號的有用部分。利用小波分析的多尺度特性,把含有噪聲的信號分解到多個尺度上,再在每一尺度上分別設定對應的閾值,其中,將噪聲信號歸納為小于該閾值的小波系數,將有用信號歸納為大于該閾值的小波系數,分別去除各尺度上的噪聲信號,同時保留有用信號,即可重構出去噪后的信號[11]。可以看出,同閾值函數一樣,小波基函數和分解層數的選擇對小波閾值的去噪結果也有很大影響。
一個含噪聲的一維信號可表示為[12]
式中:f(t)為 含有噪聲的信號;s(t) 為 純凈信號;σ為噪聲標準差;w(t)為噪聲信號。
有用信號一般為低頻信號,而噪聲通常為高頻信號。例如,對原始信號進行3 層分解:
式中:cAi為 分解的近似部分,cDi為分解的細節部分,i =1,2,3,噪聲部分通常包含在c D3、cD2、cD1中。
總結去噪過程,主要包括以下3 個步驟[13]:
1)對原始信號進行小波分解。選擇最佳的小波基和分解層數,將原始信號進行多尺度分解得到相對應的小波系數。
2)小波分解閾值量化。選擇相對應的閾值及閾值函數,對小波系數進行閾值量化處理。
3)小波系數重構。將閾值化處理后的小波系數進行重構,得到去噪后的信號。
從上可知,最優的小波閾值去噪方案,不僅與閾值函數的選取有關,小波基函數、分解層數的選擇都直接影響信號去噪的質量,而對于小波基函數與分解層數的選擇大多基于經驗或直接給定得到,未充分分析其對去噪效果的影響,如果能選擇最優的小波基函數和分解層數,對于提高去噪效果將會有明顯的幫助[14]。
小波基函數及分解層數都直接影響小波去噪效果,而使用不同的小波基與分解層數,去噪后的效果也不相同。如果分解尺度過小,則信號中仍存在較多的噪聲數據,無法獲得理想的去噪結果;如果分解尺度過大,則會將信號中的部分細節信息當做噪聲刪除,容易造成信號失真,同時增加計算復雜度,而去噪效果卻沒有顯著改善[15]。因此,在去噪效果和分解層數之間找到平衡非常重要,需要構建一種全面的質量評價指標來衡量去噪效果。
當前,RMSE、SNR、平滑度r 等指標常被用于分析小波去噪的效果[16]。RMSE 指去噪信號與原始信號之間方差的平方根,其值越小說明去噪效果越優;SNR 是原始信號功率與噪聲信號功率之間的比值,其值越大說明去噪效果越優[17];平滑度是去噪信號一階差分與原始信號一階差分之間方差根比值,代表信號整體的變異信息,其值越小說明去噪效果越優。
為了評估這些指標的有效性,本文采用一組真值已知、小波去噪專用的測試數據即多普勒信號和一組真值未知的陀螺輸出數據來進行小波分析,同陀螺信號類似,多普勒信號是一種非平穩信號,在原始信號上加入8 dB 的高斯白噪聲,以及系數為0.1 的隨機游走,信號長度為1 024 個點,采樣頻率為1 Hz。由于已知未加噪聲的純凈信號,可求得重構信號與理論純凈信號的RMSE 及SNR,當RMSE最小且SNR 最大時,所配參數即為最佳的去噪參數。真實信號為采集的加拿大諾瓦泰公司SPANISA-100C 型高精度光纖慣導的靜態陀螺數據,時長為5 h,采樣頻率為200 Hz。利用上述評價指標來判斷最優分解層數,對于2 類信號,都以sym 4 作為小波基函數,分解層數為2~10。圖1 為真值已知情況下多普勒信號小波去噪的質量評價結果,從RMSE 和SNR 結果中可以看出,最優分解層數為6 層,平滑度無法判定。圖2 為真值未知情況下陀螺輸出數據小波去噪的質量評價結果,此時無法找出最優分解層數,這是因為所測信號的真值是未知的,參與計算的是含有噪聲的原始信號與去噪信號,如果未去掉任何噪聲或去噪效果差,但得到的RMSE 值卻很小,SNR 值很大,這與評價依據不一致,縱然平滑度是減小的,但卻無法判斷最優層數。通過分析可得在真值未知情況下3 種評價指標特點如表1 所示。

表1 真值未知時評價指標特點Table 1 Characteristics of evaluation indexes when truth value is unknown

圖1 真值已知的多普勒信號的單一評價指標趨勢Fig.1 Trends of single evaluation index of Doppler signal with known truth value

圖2 真值未知的陀螺輸出數據的單一評價指標趨勢Fig.2 Trends of single evaluation index of gyro output data with unknown truth value
因此,單一的評價指標具有一定的欺騙性,難以滿足慣性元件誤差分析的需求。若將多類評價指標進行賦權融合,以更加全面的視角來構建復合評價指標,從而滿足在真值未知情況下小波去噪參數的選擇。
考慮到模型需簡單且具有可行性,根據RMSE、SNR 及平滑度的幾何與物理意義,選取RMSE 與平滑度2 個指標來構建復合評價指標。主要原因在于:RMSE 和SNR 都是關注信號的細節信息,而平滑度關注的是信號近似信息,為反映去噪信號的整體特征,應使用平滑度指標,同時,考慮到第1 節實驗中SNR 值與平滑度值隨分解層數的增加都減小,而RMSE 值隨分解層數的增加而增大,故用平滑度與RMSE 構建復合評價指標。隨著分解層數不斷增加,必定出現一個極值,即去噪信號的近似信息與細節信息達到了最優的組合,對應的是最優分解層數。
由于RMSE 與平滑度值變化的范圍不同,先將其歸一化[18],使2 個指標都處于[0,1]區間。歸一化算法為
式中:PRMSE和 Pr分別為歸一化后RMSE 和平滑度的值;m in和 max分別為最小值和最大值函數。
由于這2 個指標對信號特征的描述程度不同,在復合過程中權重也不一致,需要進行賦權處理。目前,有關評價指標權系數的確定方法很多,主要分為主觀賦權法、客觀賦權法、主客觀集成賦權法3 類[19],其中客觀賦權法主要以其強大的數學理論做支撐,相比其他方法更加可靠。本文基于信息熵權法和變異系數法2 種客觀賦權思想,提出了一種基于熵權與變異系數組合賦權的復合評價指標。
首先,以RMSE 為例,給出熵權法定權的計算公式為

式中:σPRMSE和 μPRMSE分別為RMSE 的標準差和均值;σPr和 μPr分 別為平滑度的標準差和均值;CVPRMSE和CVPr分 別為RMSE 變異系數和平滑度變異系數;WPRMSE和 WPr分別為RMSE 和平滑度的權值。
得到上述2 種賦權方法所計算的權值后,為了彌補各自不足,最大限度地減少有用信息的損失,將熵權法與變異系數法相結合,構建組合賦權法,從而使指標的權值更加客觀合理,得到最終權值,則有

在本文中,由于熵權法和變異系數法沒有明確的可靠性比較,設為0.5 將組合賦權得到的RMSE 與平滑度的權值分別與歸一化得到的RMSE 與平滑度值相乘并作和,得到復合評價指標 H,其表達式為
根據RMSE 與平滑度的定義及定權過程可知,隨著分解層數的增加,復合評價指標 H將會出現一個極值,且這個極值為最小值,所對應的分解層數即為最優分解層數。
為了能夠準確比較驗證本文所提出的復合評價指標在判斷最優參數方面的有效性,采用第1 節真值已知的仿真信號進行分析。
分析時采用控制變量法,即在小波基不變的情況下,采用不同分解層數對含有噪聲的信號進行處理,并選用不同的小波基分別對含噪信號進行處理。其他參數保持不變,閾值估計方法采用基于史坦的無偏似然估計原理(rigrsure)及軟閾值處理函數。為驗證方法的適用性,本文選取了文獻[9]和文獻[10] 2 種較為成熟的方法進行對比。
以sym4 小波基及db5 小波基為例,先計算多普勒信號在真值已知情況下(原始信號為無噪聲的純凈信號)的RMSE 及SNR,再計算其在真值未知情況下(原始信號為加噪聲的信號)的RMSE 與平滑度,并計算本文的復合評價指標 H、采用文獻[9]的指標 T 和文獻[10]的指標S,結果如表2 和表3 所示。
從表2 可以看出,選用sym4 小波基進行去噪處理時,在真值已知情況下,當分解層數為6 時,RMSE 最小,同時SNR 最大,說明最優分解層數為6;在真值未知情況下,隨著分解層數的增加,RMSE不斷增加,平滑度r 不斷降低,這與前文判斷一致,無法判斷最優分解層數。根據復合評價指標 H,當分解層數為6 時出現最小值,說明最優分解層數為6,這與實際情況相符;文獻[9]的指標 T也有類似的結論;而文獻[10]的評價指標 S從第5 層開始變化率趨于平緩,判定最優分解層數為5,與實際情況不符,同時對于判斷變化率是否趨于平緩具有一定程度的主觀性。

表2 不同分解層數利用sym4 小波基處理的評價指標Table 2 Evaluation indexes of sym4 wavelet basis processing for different decomposition layers
從表3 可以看出,選用db5 小波基進行去噪處理時,在真值已知情況下最優分解層數為7,而在真值未知情況下,利用指標 H得到的最優分解層數為7,與實際情況相符;指標 T得到的最優分解層數為6,與實際情況不符;指標 S得到的最優分解層數為7,與實際情況相符。由此可以看到,復合評價指標H相比現有方法精度更高,具有更高的可信度。

表3 不同分解層數利用db5 小波基處理的評價指標Table 3 Evaluation indexes of db5 wavelet basis processing for different decomposition layers
為了更加全面地分析本文方法的可靠性,避免單一小波基函數對實驗結果的誤導,在真值未知情況下分別選擇不同類型的小波基進行去噪分析,得到復合評價指標,從而判斷最優分解層數,并以真值已知情況下得到的最優分解層數為參考,對比本文方法與現有方法,結果如表4 所示。

表4 不同小波基對應的最優分解層數Table 4 Number of optimal decomposition layers corresponding to different wavelet bases
從表4 可以得到,本文方法的正確率為80%,文獻[9]和文獻[10]所用方法的正確率分別為53%和47%,由此可看出,相比于現有方法,本文所提方法具有更高的準確率及可信度,而且可廣泛應用于各種類型的小波基,是一種可靠的小波去噪質量評價方法。
為驗證本文方法對慣性元件隨機誤差降噪處理的適用性,分別利用加拿大諾瓦泰公司SPAN-ISA-100C 型高精度光纖慣導和邁普時空研發的MPM39 型MEMS 慣性測量單元進行數據的采集。將2 款不同精度慣導設備安裝在靜態環境中,并預熱30 min 后開始工作。為了使2 款設備的微小誤差能夠表現出來,需盡可能采集多的數據,為了達到光纖慣導采樣時間,將2 款設備同樣采集5 h 實驗數據[20],記錄得到本次實驗測試數據樣本,采樣頻率設為200 Hz。
首先,利用Allan 方差分析誤差信號的特點,Allan方差的最大特點就是能夠對各種誤差源進行單獨辨識分析。2 款陀螺儀的誤差分析結果如表5 和表6 所示。

表5 SPAN-ISA-100C 陀螺儀Allan 方差分析結果Table 5 Allan variance analysis results of SPAN-ISA-100C gyroscope

表6 MP-M39 陀螺儀Allan 方差分析結果Table 6 Allan variance analysis results of MP-M39 gyroscope
由表5 和表6 可知,陀螺儀的主要誤差來自于角度隨機游走、零偏不穩定性及角速率隨機游走,隨著陀螺精度的降低,3 軸各項誤差源也隨之變得不穩定,同時誤差系數變大。
在小波去噪之前,確保其他參數一致,采用rigrsure 閾值選取準則、軟閾值處理函數,選取不同的小波基函數,分解層數為2~10 層,采用控制變量法,選定sym6 小波基函數,分別利用本文復合評價指標 H 及文獻[9]的指標 T進行最優分解層數的判定。2 款慣導設備的3 軸陀螺數據所對應的最優分解層數列于表7。其中,x1、y1、z1 分別表示SPAN-ISA-100C 陀螺儀的3 軸數據,x2、y2、z2 分別表示MP-M39 陀螺儀的3 軸數據。
在實驗中同時發現,對于這2 款慣導設備共6 組數據,每組數據選取不同小波基函數得到的最優分解層數是相同的,SPAN-ISA-100C 的x 軸輸出數據最優分解層數為4,其余數據得到的最優分解層數為5。由此,可以將對應最優分解層數設置為表7 中的分解層數,其他參數保持不變,根據復合評價指標 H判定最優小波基函數,得到對應最優小波基函數如表8 所示。

表7 陀螺數據對應的最優分解層數Table 7 Number of optimal decomposition layers corresponding to gyro data

表8 陀螺數據對應的最佳小波基函數Table 8 Optimal wavelet basis functions corresponding to gyro data
可知,復合評價指標 H可以將小波閾值去噪參數中的分解層數與小波基函數的選擇統一起來,對SPAN-ISA-100C 的x 軸輸出數據選用sym6 小波基函數,4 層分解層數;y 軸輸出數據選用sym4 小波基函數,5 層分解層數;z 軸輸出數據選用db8 小波基函數,5 層分解層數,取得最優的去噪效果。對MP-M39 的3 軸輸出數據皆選用db7 小波基函數,5 層分解層數,取得最好的去噪效果。另外,對于不同的傳感器所選擇的最優參數往往不同,且對于同一傳感器,由于其所處環境不同,所選用小波去噪參數亦不一致?;诮涷炈x擇的參數缺乏理論指導,并不一定準確可靠。
同時,從表7 可以看出,利用2 種方法進行小波去噪時,SPAN-ISA-100C 的x 軸數據對應分解層數是相同的,但其他最優分解層數有所不同。為了直觀判定2 種方法所得結果的準確性,以SPAN-ISA-100C 的y 軸輸出為例,分別利用2 種方法得到降噪信號的時域曲線與頻譜曲線,結果如圖3 和圖4所示。

圖3 原始信號與降噪信號時域曲線Fig.3 Time-domain curves of original signal and denoising signal

圖4 原始信號與降噪信號頻譜曲線Fig.4 Spectrum curves of original signal and denoising signal
圖3 中,小波去噪后的信號都保留了原始信號的變化趨勢,相比于4 層分解,5 層分解降噪信號峰值曲線更加光滑,波形圖更加平穩,消除噪聲更加徹底。圖4 為信號的單邊振幅譜圖,去噪主要剔除高頻無用細節信息,同時保留了低頻有用近似信息,這也與小波去噪原理相符。5 層分解信號比4 層分解信號更有效地剔除了高頻噪聲,同時也說明5 層分解降噪信號的SNR 更高。由此可得,選取5 層分解降噪在此狀態下效果最優,這也從另一個側面說明了復合評價指標 H的有效性。
為了更加清晰地判斷小波去噪效果,圖5 給出了原始信號與分解層數為4 和5 的去噪信號的功率譜密度估計??梢钥闯?,保留了低頻有用信息,剔除了高頻無用信息,同時,對中間的過渡信息適當進行了保留,相比于4 層分解降噪信號,5 層分解降噪信號剔除了較為微弱的噪聲信號,也證明了5 層分解層數是最優分解層數。由于篇幅所限,本文不再分析其他輸出數據去噪效果對比結果。

圖5 原始信號與降噪信號功率譜密度圖Fig.5 Power spectral density diagram of original signal and denoising signal
針對傳統方法中經驗模型缺乏理論依據的不足,本文提出了一種基于組合賦權法的小波去噪質量評價方法,得出以下結論:
1)通過實驗發現,在真值未知情況下,單一評價指標無法指導小波閾值去噪過程中相關參數的選擇,需構建更加全面的質量評價指標。
2)選取RMSE、平滑度2 個指標作為評價依據,對其進行初始化,通過引入信息熵權與變異系數2 種客觀賦權法進行組合賦權,將歸一化值與權值線性組合形成了一種新的小波去噪質量評價指標,并得出了相關計算公式。
3)通過對真值已知的模擬信號和真值未知的實測信號進行處理,證明該方法能夠指導小波閾值去噪過程中相關參數的選擇,具有更高的準確性與廣泛適用性。