張晨陽, 張 亞, 李世中
(中北大學 機電工程學院, 太原 030051)
引信仿真測試是評價引信綜合性能的有效途徑,對于提升硬目標侵徹武器的攻擊能力具有重要意義。開展引信仿真測試系統的研究,能夠縮短引信的研制時間、節省研制費用、提高研制質量,為引信的總體設計和建立完善的引信評價體系提供技術保障[1]。在侵徹過程中,彈體與彈靶的作用過程十分復雜,有多種振動信號作用在彈體上,獲得的侵徹過載信號不僅包含彈體自身的剛體加速度信號,還疊加了各種噪聲信號。因此,將彈體的剛體加速度信號從侵徹過載信號中分離出來,是引信仿真測試的關鍵一步。
對于侵徹過載信號的處理,已經展開了多方面的研究。張兵等[2]在壓電傳感器的兩端加入三種減震片,通過彈體沖擊試驗發現,加入減震片后傳感器輸出信號的頻率范圍有了明顯縮小,并且三種不同材質的減震片分別可以濾除大于4 kHz、7 kHz、10 kHz的噪聲信號。黃娟等[3]提出一種小波去噪與HHT變換相結合的故障信號提取方法,先對時域信號進行小波去噪,再采用HHT變換進行時頻分析,根據分析得到一系列本征模態函數判斷軸承的故障情況。趙海峰等[4]提出一種基于奇異值分解的侵徹過載信號降噪方法,通過對侵徹過載信號進行奇異值分解,由主體奇異值分量穩定原則確定信號重構矩陣,從而去除侵徹過載信號中的噪聲分量,提取出彈體的剛體加速度信號。毋文峰等[5]采用經驗模態分解對單通道信號進行盲分離,提取含有不同故障信息的本征模態函數與原單通道信號組成多維信號,利用奇異值分解重構多通道混合信號,對多通道混合信號進行盲源分離,但由于經驗模態分解在分解過程中存在欠包絡、過包絡、模態混疊、端點效應等問題,使得信號處理結果差強人意。集合經驗模態分解[6]雖然在一定程度上解決了模態混疊問題,但同時加大了算法的計算量,且端點效應問題仍未解決。
本文提出一種基于變分模態分解的侵徹過載信號盲源分離算法。為了將單通道信號轉化為多維信號,首先對信號進行變分模態分解,由本征模態函數和單通道信號組成多維信號;然后對多維信號的自相關矩陣進行奇異值分解估計源信號數目,并計算各本征模態函數與單通道信號的相關系數,由相關系數的大小選擇相應的本征模態函數與單通道信號重構多通道信號;最后利用盲源分離中特征矩陣聯合近似對角化算法(joint approximate diagonalization of eigenmatrices, JADE)[7]對多通道信號進行盲源分離。
盲源分離是一種在信號傳輸路徑的各項參數與信號源均未知的情況下,根據輸入源信號的統計特性,由混合觀測信號還原出源信號獨立成分的信號分離算法。根據源信號數目小于、等于和大于觀測信號數目,可將盲源分離分為超定、正定和欠定盲源分離三種[8]。傳統的盲源分離算法主要針對源信號信息充足的超定和正定盲源分離,而對于欠定盲源分離,由于其觀測信息缺乏,現有的盲源分離算法不能有效地求解出源信號的估計信號[9]。
假設源信號S(t)由N個未知的獨立信號組成,經未知的信道傳送后,由M個傳感器采集獲得觀測信號X(t),盲源分離的數學模型見式(1)
X(t)=AS(t)
(1)
式中:X(t)=[X1(t),X2(t),…,XM(t)]T為M維觀測信號;S(t)=[S1(t),S2(t),…,SN(t)]T為N維源信號;A為M×N階混合矩陣。
盲源分離就是在S(t)和A都未知的情況下尋求一個N×M階的矩陣W,通過該矩陣和觀測信號X(t)分離出源信號S(t)的近似估計信號Y(t)[10],盲源分離過程如圖1所示。
Y(t)=WX(t)
(2)

圖1 盲源分離過程圖
變分模態分解(variational mode decomposition, VMD)是由Dragomiretskiy等[11]根據維納濾波和變分思想提出的一種尺度可變的信號分解算法,通過將信號自適應地分解為若干個本征模態函數,采用交替方向乘子法迭代求解各模態的中心頻率,從而將各模態解調至對應的基頻帶。與經驗模態分解的遞歸式分解不同,VMD是在變分模型的約束下進行信號的分解過程,具有完備的理論支撐,在復雜信號分析領域里具有廣泛應用[12]。
變分模態分解將源信號分解為有限個具有限帶寬的本征模態函數,每個本征模態函數圍繞在其中心頻率附近,以各本征模態函數的帶寬之和最小為約束條件,建立約束變分模型(見式(3)),尋求各個信號分量對應的本征模態函數,從而將源信號自適應地分解為各本征模態函數的線性組合。
(3)
式中:{uk}={u1,u2,…,uK}為K個本征模態函數;{ωk}={ω1,ω2,…,ωK}為K個本征模態函數對應的中心頻率。
為了求解約束變分模型,引入二次懲罰因子α和拉格朗日懲罰算子λ(t),將變分模型最優解的約束變分問題轉化為非約束變分問題,通過交替方向乘子法在傅里葉域內求解該非約束變分問題,表達式見式(4)

(4)
通過不斷更新各本征模態分量的中心頻率和帶寬,并根據源信號的頻域特性劃分頻帶,從而將源信號自適應地分解為一系列本征模態函數,其中各變量的更新見式(5)~式(8),算法求解流程如圖2所示。
(5)

圖2 VMD算法流程圖
(6)
(7)
(8)
式中:n=n+1;k=1,2,…,K;τ為預設誤差。
在盲源分離算法中,要求觀測信號數目大于或等于源信號數目,而在實際情況中觀測信號數目往往小于源信號數目,且有時只能采集到單通道觀測信號,為了解決單通道觀測信號盲源分離問題,提出一種基于VMD的盲源分離算法。首先由VMD將單通道觀測信號自適應地分解為若干個本征模態函數,并將單通道觀測信號與各本征模態函數組成多維信號;然后利用奇異值分解估計單通道觀測信號的源信號數目[13],并計算各本征模態函數與單通道觀測信號的相關系數,根據源信號數目,選擇相關系數較大的本征模態函數與單通道觀測信號重構多通道觀測信號;最后采用盲源分離中的JADE算法對多通道觀測信號進行盲源分離,算法流程如圖3所示。

圖3 基于VMD的盲源分離算法流程圖
首先選擇適宜的采樣頻率和本征模態函數個數K,通過對單通道觀測信號X(t)進行變分模態分解,從而將X(t)自適應地分解為K個本征模態函數X(t)IMF=[IMF1,IMF2,…,IMFK];然后將X(t)與K個本征模態函數組成多維觀測信號y(t)=[X(t),IMF1,…,IMFK],即解決了盲源分離中觀測信號數目小于源信號數目的問題。對y(t)的自相關矩陣進行奇異值分解,y(t)的自相關矩陣Ryy見式(9),對Ryy進行奇異值分解見式(10)
Ryy=E[y(t)×yH(t)]
(9)
Ryy=UΛVT
(10)
式中:U為左奇異向量;V為右奇異向量;Λ=diag{λ1,λ2,…,λn}為奇異值向量對應的特征值;yH(t)為多維觀測信號y(t)的共軛轉置。
Λ中前幾個特征值的數值較大,已經包含信號中的大部分信息,通過計算自相關矩陣Ryy的特征值的數值占比(見式(11)),即可估算出單通道觀測信號的數目。
(11)
當m值取i(1≤i 確定源信號的數目后,通過計算各本征模態函數與X(t)IMF的相關系數(見式(12)),選擇相關系數較大的前m-1個本征模態函數,與X(t)重構多通道觀測信號H(t)=[X(t),IMF1,…,IMFm-1]。 ρ(X(t),IMFi)= (12) 采用盲源分離中的JADE算法對多通道觀測信號H(t)進行分析處理。首先對H(t)進行去中心化和白化處理,求得白化矩陣W,由盲源分離模型可得白化信號Z(t)見式(13) Z(t)=WX(t)=WAS(t)=US(t) (13) 式中,U為一個N×N階的酉矩陣。 從而將求解混合矩陣A轉化為計算酉矩陣U。取一任意的N×N階矩陣M,則白化信號Z(t)的四階累積量矩陣的第a行第b列元素為見式(14) (14) 式中:QZ(M)為Z(t)的四階累積量矩陣;Pabcd(Z)為Z(t)的第a,b,c,d分量的四階累計量;mcd為矩陣M的第c行第d列元素。 將UTQZ(Mi)U中對角元素的平方和(見式(15))作為矩陣對角化的目標函數,通過最小化目標函數F(U)求解出酉矩陣U,即可求出源信號的估計信號Y(t)見式(16)。 (15) Y(t)=UTWX(t) (16) 采用彈載測試系統進行侵徹試驗以獲得侵徹過載信號時間歷程曲線。試驗過程如下:侵徹試驗前在試驗彈體底部安裝彈載測試系統,在彈體侵徹過程中由彈載測試系統中的傳感器模塊采集彈體與彈靶產生的侵徹過載信號;侵徹過載信號經信號調理模塊升壓、整流、放大、濾波后,傳送到數據采集存儲模塊中進行模數轉換并存儲到數據存儲器中;侵徹試驗結束后,通過彈載測試系統的數據接口讀取試驗數據。測試系統架構如圖4所示,試驗結束后,通過測試系統的數據接口將數據讀取出來。試驗測得的彈體侵徹過載信號曲線如圖5所示,試驗測得彈體侵徹深度為103.52 cm,彈體觸靶速度為478.3 m/s。 圖4 測試系統架構圖 圖5 侵徹過載信號時間歷程曲線 以1 kHz的采樣頻率對試驗獲得的侵徹過載信號進行變分模態分解。通過對K取2~10,將各本征模態函數的中心頻率按從低到高排列,觀察其中心頻率的變化情況,如表1所示。由表1可知,當K取8時,中心頻率在454 Hz左右的本征模態函數信息丟失;K取10時,第8個分量的中心頻率738 Hz和第9個分量的中心頻率787 Hz較為接近,容易出現模態混疊現象。因此將本征模態函數的分解個數K設置為9,各本征模態函數圖像如圖6所示。 表1 不同模態個數時各模態的中心頻率 (a) IMF1~IMF5函數圖像 將9個本征模態函數與侵徹過載信號組成多維觀測信號y(t),對y(t)的自相關矩陣Ryy進行奇異值分解,得到各奇異值向量對應的特征值,并計算各特征值及前n個特征值在全體特征值中所占比例,如表2所示。 表2 特征值占比及前n個特征值占比 由表2可知,當n=3時,前3個特征值在總體特征值中占比為98.59%;當n>3時各特征值占比均小于1%,表明前3個特征值已經包含多維觀測信號y(t)中的大部分信息,因此源信號的估計數目為3。然后計算各本征模態函數與侵徹過載信號的相關系數,選擇相關系數最大的2個本征模態函數與侵徹過載信號重構多通道觀測信號,計算結果如表3所示。 表3 各本征模態函數與侵徹過載信號相關系數 由表3可知,η5與η7在9個相關系數中數值較大,因此選擇IMF5、IMF7與侵徹過載加信號X(t)重構多通道觀測信號H(t)=[X(t),IMF5,IMF7],采用JADE對多通道觀測信號H(t)進行盲源分離,分離出的信號分量如圖7所示。 (a) 分量1 由圖7可知,H(t)經本文方法處理后分離出3個信號分量。其中,分量1的加速度幅值曲線較為光滑,由式(14)計算分量1與侵徹過載信號的相關系數ρ=0.952 7,可認為分量1包含了侵徹過載信號曲線的大部分信號特征;分量2和分量3為高頻噪聲振動信號。對分量2和分量3進行頻譜分析如圖8所示,侵徹過載信號頻譜如圖9所示。 圖8 分量2和分量3頻譜圖 圖9 侵徹過載信號頻譜圖 由圖8可知,分量2和分量3的頻率集中點主要在1 072.7 Hz和3 875.4 Hz附近,與圖9中的頻率集中點1 061.1 Hz和3 863.9 Hz較為接近。因此,可認為侵徹過載信號經本文方法處理后,分離出的信號分量1與原信號的相似度約為95%,反映了彈體在侵徹過程中的剛體加速度特征;同時還分離出兩個頻率主要集中在1 072.7 Hz和3 875.4 Hz的高頻振動信號分量2和分量3,由于彈體在侵徹過程中存在復雜的應力作用,因此這兩個頻率集中點可能為彈體的兩次或多次諧振頻率。 為了檢驗本文方法處理侵徹過載信號的可靠性和有效性,采用小波分解法和文獻[5]中經驗模態分解與獨立分量分析法處理侵徹試驗獲得的侵徹過載信號,提取出的侵徹過載信號曲線如圖10所示。 圖10 三種方法提取出的侵徹過載信號曲線 由圖10可知,本文方法分離的侵徹過載信號的曲線光滑度更高;文獻[5]方法提取的侵徹過載信號在曲線快速上升和快速下降過程中存在曲線尖點問題;小波分解法提取的侵徹過載信號缺失部分細節特征,且該方法需人工選擇小波基函數、分解層數和閾值,算法本身缺乏自適應性。為了進一步比較三種方法的侵徹過載信號特征分離效果,先對每條侵徹過載信號曲線分別在時域上進行一次積分,得到三條侵徹速度變化曲線,如圖11所示;然后對每條侵徹速度變化曲線分別在時域上再進行一次積分,得到三條侵徹深度變化曲線,如圖12所示。 圖11 侵徹速度變化曲線 圖12 侵徹深度變化曲線 由圖11、12可知,本文方法得到的彈體觸靶速度為485.9 m/s,彈體侵徹位移為100.61 cm;文獻[5]中的方法得到的彈體觸靶速度為458.7 m/s,彈體侵徹位移為97.18 cm;小波分解法得到的彈體觸靶速度為465.5 m/s,彈體侵徹位移為98.44 cm。為了更直觀地對比三種處理方法的處理效果,將三種處理方法得到的試驗結果與實際侵徹數據進行對比,如表4所示。 表4 三種處理方法的試驗結果對比 由表4可知,與文獻[5]方法和小波分解法相比,本文方法得到的觸靶速度和侵徹深度與侵徹試驗數據的誤差最小,所需處理時間最短,信號處理效率最高,因此本文方法具有更優的處理效果以及更高的可靠性和可行性。 變分模態分解能夠自適應地將輸入信號分解為若干個本征模態函數的線性組合,在保留信號局部特征規律的同時,有效去除了噪聲信號,為單通道信號的盲源分離奠定基礎。通過對本征模態函數的自相關矩陣進行奇異值分解并計算各本征模態函數與單通道觀測信號的相關系數,將單通道觀測信號轉化為多通道觀測信號,解決了欠定盲源分離中觀測信號數目小于源信號數目的問題,然后由盲源分離中的JADE算法對重構的多通道觀測信號進行分析處理。通過對分離的侵徹過載信號進行進一步分析處理,檢驗了本文方法的可靠性和可行性。與傳統信號處理方法相比,本文方法具有以下優點: (1) 本文方法結合時頻分析在非平穩信號分析中的優點,能夠有效分離侵徹過程中作用在彈體上的多種頻率響應,分離效果良好; (2) 盲源分離在混疊信號分離方面具有獨特的優勢,能夠深入研究侵徹過程中作用在彈體上的復雜振動響應,挖掘出蘊含在復雜波形中的深層信息; (3) 本文方法很好地解決了欠定盲源分離中的源數估計問題,實現了侵徹過載信號的盲源分離。通過分析侵徹過載信號,得到侵徹過載信號的數據統計特征,挖掘其中蘊含的深層信息,從而分析彈體的侵徹過程,可為引信結構設計、彈體強度設計、防御工事的材料選擇和結構設計提供重要參考依據。2.2 多通道觀測信號盲源分離
3 侵徹試驗分析
3.1 侵徹試驗


3.2 侵徹過載信號盲源分離







3.3 試驗結果分析




4 結 論