基于離散小波分解的火工沖擊數(shù)據(jù)有效性分析與校正方法

王錫雄， 秦朝燁， 丁繼鋒， 褚福磊

(1.清華大學(xué) 機械工程系，北京　100084；2. 北京空間飛行器設(shè)計總體部， 北京　100094)

?

基于離散小波分解的火工沖擊數(shù)據(jù)有效性分析與校正方法

王錫雄1， 秦朝燁1， 丁繼鋒2， 褚福磊1

(1.清華大學(xué) 機械工程系，北京100084；2. 北京空間飛行器設(shè)計總體部， 北京100094)

摘要：由于火工沖擊試驗條件的限制，測試所得沖擊加速度信號中通常含有低頻趨勢項，導(dǎo)致試驗數(shù)據(jù)不準確。因此需要對信號進行有效性分析，校正信號，消除積分零漂?，F(xiàn)有的校正方法大都是基于經(jīng)驗的定性方法，具有一定的局限性?；陔x散小波分解，提出了一種沖擊數(shù)據(jù)有效性分析方法，該方法通過分析信號的頻譜與相關(guān)性，能夠定量確定小波分解層數(shù)的范圍，還討論了Db小波基函數(shù)的消失矩對數(shù)據(jù)校正結(jié)果的影響。通過對某分離沖擊試驗數(shù)據(jù)進行有效性分析與校正，證明該方法能夠有效抑制積分零漂，消除趨勢項對沖擊譜的影響，提高試驗數(shù)據(jù)的準確性和有效性。

關(guān)鍵詞：火工沖擊；離散小波分解；沖擊響應(yīng)譜；分解層數(shù)；消失矩

由于火工沖擊高頻、瞬態(tài)、高量級的力學(xué)環(huán)境特點，沖擊試驗對采集系統(tǒng)要求很高，一般要求數(shù)據(jù)帶寬[0,10 000]Hz以上、采樣頻率為數(shù)據(jù)最高頻率的10倍以上、模擬反混疊濾波器衰減率達到60 dB/倍頻程，這給信號測試帶來了較大困難和一定的隨機性[1]。傳感器本身的諧振、數(shù)據(jù)采集系統(tǒng)性能不足、噪聲干擾等現(xiàn)象都可能導(dǎo)致測得的信號與真實信號之間存在趨勢項[2-3]。這種趨勢項在采用卷積方法求解沖擊響應(yīng)譜的過程中會被放大，產(chǎn)生零位漂移(Zero shift)，影響試驗數(shù)據(jù)的有效性。Chu[4]發(fā)現(xiàn)沖擊加速度零漂以來，目前的沖擊試驗水平還無法杜絕零漂現(xiàn)象，因此有效性的分析與校正，是沖擊數(shù)據(jù)處理的重要步驟。

沖擊數(shù)據(jù)有效性的判定方法已經(jīng)比較成熟，美軍標與NASA等相關(guān)標準[5-7]都建議采用零漂法和正負沖擊譜法進行判定。而沖擊數(shù)據(jù)的校正主要是從所測信號中分離低頻趨勢項的過程，主要方法包括最小二乘法、低通濾波法、離散小波變換(Discrete Wavelet Transform，DWT)和經(jīng)驗?zāi)Ｊ椒纸?Empirical Mode Decomposition，EMD)等[8-14]。前兩種方法不適合處理非平穩(wěn)的沖擊信號，且需要對信號具有一定的先驗知識[8]。EMD方法具有自適應(yīng)性，但分解所得本征模函數(shù)(IMF)均值不為0，因此要保證信號消除積分零漂， EMD通常結(jié)合最小二乘法使用[10-11]。小波基函數(shù)具有多項式正交性，因此DWT能夠有效提取低階多項式趨勢項，是一種簡單通用，無需先驗知識的方法。DWT方法最早由Smallwood[12]引入，用于校正積分零漂的沖擊試驗數(shù)據(jù)，隨后由Edward[13]給出了試驗驗，袁宏杰等[14]也采用DWT進行了沖擊數(shù)據(jù)分析與處理。以上針對DWT的研究主要是根據(jù)經(jīng)驗進行小波系數(shù)的選擇，沒有統(tǒng)一定量的選擇依據(jù)，需要花費大量的時間與精力逐一驗證校正結(jié)果。

為了改進這一缺點，本文基于離散小波分解，提出了一種沖擊數(shù)據(jù)有效性分析方法，該方法通過分析信號的頻譜與相關(guān)性，能夠定量確定小波分解層數(shù)范圍。最后本文采用該方法對某分離沖擊試驗數(shù)據(jù)進行了有效性分析與校正，驗證了該方法是正確有效的。

1火工沖擊數(shù)據(jù)有效性

常用的沖擊數(shù)據(jù)有效性判別方法有零漂法和正負沖擊譜法，兩種方法具有一致性。零漂法如圖1所示，將實測沖擊加速度信號S1進行時域積分，得到的速度與位移若不收斂于零，則表明存在零位漂移。

圖1　零漂現(xiàn)象Fig.1 Zero shift phenomenon

計算沖擊響應(yīng)譜(Shock Response Spectrum， SRS)時[15]，以含趨勢項的S1作為激勵，使用Duhamel積分[16]計算單自由度系統(tǒng)的響應(yīng)。由于趨勢項的影響，計算所得的響應(yīng)也會產(chǎn)生零漂，進而產(chǎn)生如圖2所示正負沖擊響應(yīng)譜不一致的現(xiàn)象。

圖2　正負譜不一致Fig.2 Differences between positive and negative SRS

在數(shù)據(jù)有效性檢驗過程中，若發(fā)生以上兩種現(xiàn)象，則說明數(shù)據(jù)有效性存在問題，數(shù)據(jù)需要經(jīng)過校正才可使用。

2數(shù)據(jù)校正方法-DWT

離散小波變換是將小波因子離散化標準化的小波變換，在信號處理、圖像處理和故障診斷等領(lǐng)域有著廣泛應(yīng)用[17-19]。下面基于DWT方法，從小波理論出發(fā)，分析DWT各參數(shù)的意義以及小波參數(shù)與數(shù)據(jù)有效性的關(guān)系，提出定量選取小波參數(shù)的方法。

2.1DWT理論

對連續(xù)小波變換的尺度因子和時間因子進行離散化，得到離散小波變換[17]：

(1)

求得小波系數(shù)后，可以將原信號f(x)分解為不同尺度因子與時間因子的小波基函數(shù)的線性組合：

(2)

(3)

根據(jù)Mallat快速算法(Fast Wavelet Algorithm, FWA)，小波尺度因子反映了信號的頻率信息，離散小波分解可以看作一組高低通濾波器[18]，其分解后頻帶如圖3所示。對于沖擊信號，高頻部分Dj為信號的主要特征，趨勢項包含在最低頻率AN中。根據(jù)奈奎斯特頻率freN與分解層數(shù)N，可估算趨勢項AN的頻率范圍fre(AN)[17]：

fre(AN)∈[0,freN/2N]

(4)

圖3　DWT頻帶Fig.3 DWT frequency band

2.2小波關(guān)鍵參數(shù)的選取

如式(1)所示，離散小波的尺度平移系數(shù)都作了離散化規(guī)范，分解效果主要取決于分解層數(shù)N和小波基函數(shù)。為了避免截斷誤差，DWT要求采用正交緊支的小波基函數(shù)[17]，一般采用Db小波(Daubechies Wavelet)，因此基函數(shù)的選取可轉(zhuǎn)化為Db小波消失矩(Vanishing Moments)的確定。

(1) 分解層數(shù)的確定

確定小波分解層數(shù)的本質(zhì)在于，找到一個合適的正整數(shù)N，使原始信號進行N次DWT分解后滿足：① 校正后的信號fd中包含信號的中高頻主要特征；② 去除的低頻信號fa中包含導(dǎo)致積分零漂的誤差趨勢項。

首先根據(jù)信號頻譜確定主要頻率范圍，如圖4所示，沖擊信號S1的頻率成分很廣，且頻率成分主要集中在中高頻。信號在1 000 Hz之后頻率成分能量較大，可初步判斷信號的趨勢項頻帶為[0,1 000]Hz，通過式(4)可估算對應(yīng)的分解層數(shù)N>5。

圖4　信號頻譜與分解頻帶Fig.4 S1 frequency spectrum and the trend frequency range

為了進一步定量確定分解層次，引入fd與原始信號的Pearson相關(guān)系數(shù)為校正相關(guān)系數(shù)Xorfa，fa與原始信號的相關(guān)系數(shù)為殘余相關(guān)系數(shù)Xorfd，以便評價校后信號與原信號的差別：

(5)

式中：cov(x,y)為變量x與y的協(xié)方差，σx為變量x的標準差。

如圖5所示，隨著分解層數(shù)的增加，fa的頻帶變窄，與原信號的相關(guān)性性逐漸降低，而fd與原信號的相關(guān)性性逐漸增加。為了保證fa中包含信號的主要特征，設(shè)定閾值確保Xorfa>0.9，為強相關(guān)；同樣設(shè)定Xorfd<0.1，使分離剔除的近似分量與原始信號弱相關(guān)，通過這種方法可以確定小波分解層數(shù)的下界。

圖5　殘余相關(guān)系數(shù)和校正相關(guān)系數(shù)Fig.5 Residual correlation coefficient and corrected correlation coefficient

令校正前信號f的兩次積分與fa的兩次積分的互相關(guān)系數(shù)為位移零漂系數(shù)XorZ，這一參數(shù)用于判斷積分后的零漂趨勢是否一致，反映了fd中對趨勢項的包含程度：

(6)

由于母小波均值為0，其重構(gòu)的細節(jié)分量fd均值為0，主要的趨勢項都包含在父小波重構(gòu)的fa中，因此如圖6所示，位移零漂系數(shù)都很大。經(jīng)過兩次積分，信號的高頻振蕩被抑制，體現(xiàn)的主要是趨勢項的零漂特性。取XorZ開始減小時的分解層數(shù)即可以確定分解的上界。

圖6　位移零漂系數(shù)Fig.6 Zero shift correlation coefficient of displacement

確定了小波分解層數(shù)的上下界范圍之后，結(jié)合信號的頻率范圍，對比范圍內(nèi)的分解效果，選用校正結(jié)果最好的層級對信號進行校正。

(2) 消失矩對校正的影響

如式(7)所示，具有M階消失矩的小波基函數(shù)，與k次線性多項式函數(shù)正交，?k∈Z,0

(7)

因此為了在細節(jié)分量fd中更好地保留高階變化，消失矩不宜過小，一般取k≥3[12-13]。如圖7所示，使用不同消失矩的Db小波對信號進行7次分解提取的趨勢項fa。Db2分解的fa中依然存在許多未被識別的高階變化分量，這些變化分量都在校正后被舍棄，可能會影響校正效果。而Db5與Db7的結(jié)果就比較平滑。當消失矩取值在合理的范圍時，fd中已經(jīng)能夠消去大多的多項式趨勢項，再增加消失矩對校正結(jié)果影響不大。

圖7　不同消失矩Db小波分解提取的faFig.7 Trend decomposed by DWT of Db2, Db5 and Db7

3沖擊試驗數(shù)據(jù)校正實例

本節(jié)以某星箭分離沖擊試驗實測數(shù)據(jù)為例，進行有效性分析與校正，驗證本文方法的有效性。該試驗為單機級沖擊響應(yīng)驗證試驗，采用真實火工品作為沖擊源；采集系統(tǒng)使用了壓電式加速度傳感器，采樣頻率為100 kHz，奈奎斯特頻率為50 kHz。圖8所示為試驗實測沖擊加速度信號S2，對S2進行兩次數(shù)值積分，發(fā)生了明顯的零漂現(xiàn)象。觀察圖9所示信號S2的頻譜，可以判斷到低頻趨勢項的大致范圍為[0,1 000]Hz。

圖8　分離沖擊實測信號S2的積分零漂Fig.8 Zero shift of S2

圖9　信號S2的頻譜與分解頻帶Fig.9 S2 frequency spectrum and the trend frequency range

然后根據(jù)2.2節(jié)的參數(shù)選擇方法，進一步通過3個相關(guān)系數(shù)確定分解層數(shù)范圍為[5,10]，經(jīng)過對比分析，最終確定分解層數(shù)為6。使用小波基函數(shù)Db7，將S2進行6次分解重構(gòu)得到的各尺度分量如圖10所示，最終分解得到的近似分量A6為一個較平滑的瞬態(tài)沖擊趨勢項。

圖10　S2的6層分解結(jié)果Fig.10 Level 6 decomposition of S2

圖11　S2的DWT校正結(jié)果Fig.11 DWT correction of S2

對趨勢項fa進行兩次數(shù)值積分得到的結(jié)果如圖12所示，可以看出：

(1) 趨勢項fa幅值范圍為[-150,150]m/s2，對比原始信號的幅值[-10 000,10 000]m/s2很小，殘余相關(guān)系數(shù)Xorfa=0.033 04。這說明分離出的趨勢項與原始信號為弱相關(guān)，fa中未包含原始信號的主要能量成分。

(2) 將fa與S2時域積分后，無論是速度還是位移，零漂趨勢幾乎重合，位移零漂相關(guān)系數(shù)XorZ=1，二者完全正相關(guān)。說明最終分解得到的fa中包含導(dǎo)致積分零漂產(chǎn)生的低頻趨勢項部分。

圖12　趨勢項fa與原信號S2對比Fig.12 Comparison of fa and S2 in time domain

如圖13所示,時域上的校正前后信號對比，fd與S2的加速度曲線幾乎重合，校正相關(guān)系數(shù)Xorfd=0.998 8。說明校正前后信號為強相關(guān)，fd中包含了原信號S2中的主要能量成分。對校正后的信號fd進行兩次時域積分，積分后無論是速度還是加速度都收斂于0，零漂現(xiàn)象已經(jīng)消除。

圖13　校正前后信號對比Fig.13 Comparison of fa and S2 in time domain

對比圖14所示校正前后的頻譜，可知信號高于分解頻率的成分幾乎沒有能量損失，而低于分解頻率的趨勢項成分得到了很好的抑制。

圖14　校正前后信號頻譜對比Fig.14 Comparison of fa and S2 in frequency domain

計算修正后的正負沖擊譜,觀察圖15可知響應(yīng)的正負最值已基本一致，響應(yīng)的零漂現(xiàn)象已經(jīng)消除。

圖15　校正后信號的正負沖擊響應(yīng)譜Fig.15 Corrected positive and negative SRS

校正前后的沖擊響應(yīng)譜如圖16所示，對比可知，在[1,400]Hz的頻率范圍內(nèi)，校正前后的SRS相差了1～2個數(shù)量級。沖擊信號的數(shù)據(jù)趨勢項嚴重干擾了沖擊響應(yīng)譜對沖擊力學(xué)環(huán)境的預(yù)計。若不進行數(shù)據(jù)有效性分析，必定造成預(yù)示結(jié)果偏大，影響沖擊試驗的真實性。本文提出的方法能夠確定小波分解的關(guān)鍵參數(shù)，有效提取信號中的低頻趨勢項，實現(xiàn)信號校正。

圖16　校正前后信號的沖擊響應(yīng)譜對比Fig.16 Comparison of original and corrected SRS

4結(jié)論

本文基于傳統(tǒng)的沖擊分離數(shù)據(jù)處理技術(shù)，提出了一種基于DWT的有效性分析與校正方法，引入了3個相關(guān)系數(shù)，能夠定量確定DWT分解層數(shù)。最后本文采用該方法對實測星箭分離沖擊數(shù)據(jù)進行了分析，得到了以下結(jié)論：

(1) 火工沖擊屬于高頻、瞬態(tài)、高量級的嚴酷力學(xué)環(huán)境，對沖擊信號的采集造成了一定的困難，容易造成所測信號的不準確，因此有必要對沖擊信號進行有效性分析與校正。

(2) DWT分解層數(shù)是控制信號校正的主要參數(shù)，可以結(jié)合頻譜和3個相關(guān)系數(shù)確定其范圍。

(3) 小波基函數(shù)的消失矩不宜取得過小，通常取k≥3，但只要在合理的范圍內(nèi)，消失矩對校正結(jié)果影響較小。

(4) 本文提出的基于DWT的校正方法，具有通用性強，無需先驗的特點，適合沖擊信號校正。經(jīng)試驗驗證，該方法能夠有效消除積分零漂，修正沖擊響應(yīng)譜，保證沖擊試驗數(shù)據(jù)正確有效。

參 考 文 獻

[1] Lee J R, Chia C C, Kong C W. Review of pyroshock wave measurement and simulation for space systems[J]. Measurement, 2012, 45(4): 631-642.

[2] Bateman V I, Merritt R. Validation of pyroshock data[J]. Journal of the IEST, 2012, 55(1): 40-56.

[3] Irvine T. An introduction to the shock response spectrum[DB/OL]. http://www.vibrationdata.com/tutorials2/srs_intr.pdf, 2012.

[4] Chu A. Zeroshift of piezoelectric accelerometers in pyroshock measurements[C]//Bulletin of the 57th Shock and Vibration. Washington DC,USA:SAVIAC,1987:71-80.

[5] IEST-RP-DTE032.2. Recommended practice for pyroshock testing[S]. Institute of Environmental Sciences and Technology, 2009.

[6] Mulville D R. Pyroshock test criteria, NASA technical standard[S]. Report NASA-STD-7003A, 2010.

[7] MIL-STD-810G Committee. Environmental test methods and engineering guides[S]. 2009.

[8] 龍源, 謝全民, 鐘明壽, 等. 爆破震動測試信號預(yù)處理分析中趨勢項去除方法研究[J]. 工程力學(xué), 2012, 29(10): 63-68.

LONG Yuan, XIE Quan-min, ZHONG Ming-shou, et al. Research on trend removing methods in preprocessing analysis of blasting vibration monitoring signals[J]. Eegineering Mechanics, 2012, 29(10): 63-68.

[9] Wu Z, Huang N E, Long S R, et al. On the trend, detrending, and variability of nonlinear and nonstationary time series[J]. Proceedings of the National Academy of Sciences, 2007, 104(38): 14889-14894.

[10] 陸凡東, 方向, 郭濤, 等. EMD與SLS法在爆破振動加速度信號時域積分中的應(yīng)用[J]. 振動與沖擊,2012,31(9): 90-93.

LU Fan-dong, FANG Xiang, GUO Tao, et al. Application of EMD and SLS in time integration of blasting vibration acceleration[J].Journal of Vibration and Shock,2012,31(9): 90-93.

[11] 胡燦陽, 陳清軍. 基于 EMD 和最小二乘法的基線飄移研究[J]. 振動與沖擊, 2010, 29(3): 162-167.

HU Can-yang, CHEN Qing-jun. Research on base line drift using least-square and EMD[J]. Journal of Vibration and Shock, 2010,29(3):162-167.

[12] Smallwood D O, Cap J S. Salvaging pyrotechnic data with minor overloads and offsets[J]. Journal of the IEST, 1999, 42(3): 27-35.

[13] Edwards T S. An improved wavelet correction for zero shifted accelerometer data[J]. Shock and Vibration, 2003, 10(3): 159-167.

[14] 袁宏杰, 姜同敏. 實測爆炸分離沖擊數(shù)據(jù)的分析和處理[J]. 固體火箭技術(shù), 2006, 29(1): 72-74.

YUAN Hong-jie, JIANG Tong-min. Analysis and treatment of measured pyrotechnic shock data[J]. Journal of Solid Rocket Technology, 2006, 29(1): 72-74.

[15] Smallwood D O. An improved recursive formula for calculating shock response spectra[J]. Shock and Vibration Bulletin, 1981, 51(2): 211-217.

[16] Thomson W. Theory of vibration with applications[M]. Boca Raton,FL:CRC Press, 1996.

[17] Daubechies I. Ten lectures on wavelets[M]. Philadelphia: Society for Industrial and Applied Mathematics, 1992.

[18] Fugal D L. Conceptual wavelets in digital signal processing[M]. San Diego,CA:Space & Signals Technical Publishing, 2009.

[19] 褚福磊, 彭志科, 馮志鵬, 等. 機械故障診斷中的現(xiàn)代信號處理方法[M]. 北京: 科學(xué)出版社, 2009.

基金項目：國家自然科學(xué)基金(11102096)資助

收稿日期：2015-06-08修改稿收到日期：2015-08-11

通信作者褚福磊 男，博士，教授，博士生導(dǎo)師，1959年生

中圖分類號:V416

文獻標志碼：A

DOI:10.13465/j.cnki.jvs.2016.14.001

Validation and correction of pyroshock data based on discrete wavelet decomposition

WANG Xi-xiong1, QIN Zhao-ye1, DING Ji-feng2, CHU Fu-lei1

(1. Department of Mechanical Engineering, Tsinghua University, Beijing 100084, China；2. Beijing Institute of Spacecraft System Engineering, Beijing 100094, China)

Abstract:Acceleration data measured in pyroshock experiment usually contain a low-frequency trend term, which makes the signal inaccuracy. Therefore the pyroshock data require validation and correction to avoid the integral zero shift. Traditional correction methods are mostly qualitative methods based on experience, making some limitation in application. A method based on the discrete wavelet decomposition was proposed. In the method wavelet correction parameters were chosen quantitatively via analying the frequency spectrum and correlation coefficients. The influence of vanishing moments on the result of correction was discussed. The proposed method was applied to the validation and correction of experimental data of pyroshock. It is revealed that the method can remove integral zero shift, correct the shock spectrum and improve the accuracy and effectiveness of pyroshock test systems.

Key words:pyroshock; discrete wavelet transform; shock response spectrum; decomposition level; vanishing moments

第一作者 王錫雄 男，博士生，1988年生