基于聚類求解TVAR模型的目標(biāo)微多普勒分析

祿曉飛, 靳碩靜, 洪 靈, 戴奉周,*

(1. 中國人民解放軍63620部隊, 甘肅 蘭州 732750; 2. 西安電子科技大學(xué)電子工程學(xué)院,陜西 西安 710071; 3. 陜西師范大學(xué)計算機科學(xué)學(xué)院, 陜西 西安 710061)

0 引 言

目標(biāo)的微動是目標(biāo)或其部分結(jié)構(gòu)相對于雷達在徑向上的非勻速運動,如振動、轉(zhuǎn)動、加速運動等。微動在自然界中十分常見,如行人手腳的擺動,電機、直升機旋翼的旋轉(zhuǎn),橋梁的振動,彈道導(dǎo)彈目標(biāo)的自旋、進動、章動等[1-4]。由微動引起的相位變化通常是非線性的,如目標(biāo)自旋的微多普勒為正弦類曲線形式,目標(biāo)進動的微多普勒為諧波函數(shù)形式,目標(biāo)章動的微多普勒為更加復(fù)雜的時變函數(shù)形式。由雷達目標(biāo)微動產(chǎn)生的微多普勒是時變的,這是由于在運動過程中，目標(biāo)結(jié)構(gòu)部件和其主體會產(chǎn)生相互作用[5-7]。由空間錐體目標(biāo)微動產(chǎn)生的微多普勒特征具有獨特性,能夠反映出目標(biāo)的精細特征[8-9],為雷達目標(biāo)的檢測與識別提供重要依據(jù),具有重要的研究價值和意義。

時頻分析方法作為微多普勒分析的一種有效手段,已被廣泛應(yīng)用。時頻分析方法包括參數(shù)化時頻分析和非參數(shù)化時頻分析[10],在非參數(shù)化時頻分析中,線性時頻分析短時傅里葉變換(short-time fourier transform, STFT)方法的時頻分辨率與窗長有關(guān),而窗函數(shù)受Heisenberg不確定準(zhǔn)則的限制,時間-頻率分辨率相互制約,導(dǎo)致分辨率不高;雙線性時頻分布Wigner-Ville分布方法具有對稱性、時移性、組合性等特征,可克服時頻分辨率的Hisenbergur界,分辨率較高,但其對多分量信號或具有調(diào)制規(guī)律的信號會產(chǎn)生嚴(yán)重的交叉干擾,這是二次型時頻分布的固有結(jié)果[11-12]。短時分?jǐn)?shù)階傅里葉變換[13]緩解了線性時頻分析方法不確定性原理的局限性,計算簡單,但在克服交叉項的干擾時,無法處理復(fù)雜信號。采用Hilbert-Huang變換獲得了比傳統(tǒng)Cohen類時頻分析方法更好的性能,但其對噪聲非常敏感。為了解決非參數(shù)化時頻分析方法時頻分辨率低、對噪聲缺乏魯棒性或交叉項嚴(yán)重等不足,有學(xué)者提出了基于參數(shù)化的時頻分析方法。其中,塊稀疏前后向時變自回歸(block-sparse forward-backward time-varying autoregressive, BS-FBTVAR)模型利用模型參數(shù)中的塊稀疏性對剛體目標(biāo)進行微多普勒分析,利用了相對短時間窗內(nèi)良好的頻率分辨率,并將塊稀疏貝葉斯學(xué)習(xí)(block-sparse Bayesian learning, BSBL)算法應(yīng)用于模型的參數(shù)求解[14-16],不存在交叉項的影響。

然而，使用BSBL算法對BS-FBTVAR模型的求解并沒有考慮到BS-FBTVAR模型中的塊稀疏時不變系數(shù)與其鄰域系數(shù)之間的相關(guān)性,與理想的時頻分辨率有一定差距。在本文中,使用改進的BSBL(extended BSBL, EBSBL)算法[17]來求解BS-FBTVAR模型,通過對EBSBL模型中由人工構(gòu)造的增廣向量的積分形式進行分析,得到一個新的結(jié)構(gòu)化稀疏先驗，通過適當(dāng)處理鄰域的超參數(shù)來促進相鄰稀疏系數(shù)之間的相關(guān)性。這種聚類先驗采用了與(sparse Bayesian learning, SBL)理論中的模式耦合先驗相似的公式[18-19]。剛體目標(biāo)時變自回歸(time-varying autoregre-ssive, TVAR)模型的時不變塊稀疏系數(shù)的塊邊界是已知的,通過結(jié)合該先驗信息,能夠加快算法迭代的收斂速度。實驗表明,所提算法的性能優(yōu)于傳統(tǒng)算法,不僅時頻分辨率更高,抗噪聲性能也較強。

1 BS-FBTVAR模型

本節(jié)介紹參數(shù)化時頻分析方法中的BS-FBTVAR模型[16]。

前后向預(yù)測的TVAR模型用于分析非平穩(wěn)隨機信號的瞬時頻率特性,彈道導(dǎo)彈目標(biāo)的雷達回波數(shù)據(jù)為非平穩(wěn)隨機信號。設(shè)目標(biāo)回波數(shù)據(jù)為x(n),可由p階前后向TVAR模型表示。前后向預(yù)測的TVAR模型可分別表示為如下形式:

n=p,p+1,…,Nx-1

(1)

n=0,1,…,Nx-p-1

(2)

(3)

(4)

式中：“^ ”表示取估計值。

為了盡可能地利用觀測數(shù)據(jù),使前后向模型預(yù)測誤差平均值最小,從而得到模型參數(shù)的解。根據(jù)式(1)和式(2)，可得前、后向預(yù)測誤差分別為

n=p,p+1,…,Nx-1

(5)

n=0,1,…,Nx-p-1

(6)

進而,前、后向預(yù)測誤差總和為

(7)

求解TVAR模型參數(shù)ai(n)(i=1,2,…,p)，需滿足ξ最小。時變系數(shù)ai(n)不易直接求解,因此將前后向TVAR模型的時變系數(shù)用基函數(shù)與時不變常系數(shù)的線性組合進行表示,將TVAR模型中時變系數(shù)的估計轉(zhuǎn)化為時不變塊稀疏常系數(shù)向量的估計[20],進而使用壓縮感知理論[21]對時不變常系數(shù)向量進行重構(gòu)。將ai(n)由一組離散余弦基fj(n)(j=1,2,…,q)[22-25]展開為

(8)

式中:q為基函數(shù)的維數(shù);aij(i=1,2,…,p；j=1,2,…,q)均為常數(shù)?；瘮?shù)的選擇有多種方式,例如離散余弦(discrete cosine transform, DCT)基、離散傅里葉(discrete Fourier transform, DFT)基、Chebyshev基、多項式基等。由于在所使用的雷達目標(biāo)回波數(shù)據(jù)中，TVAR系數(shù)在離散余弦基函數(shù)下展開的稀疏性最好,本文選用離散余弦基函數(shù)。將式(8)代入式(1)和式(2),可得

n=p,p+1,…,Nx-1

(9)

n=0,1,…,Nx-p-1

(10)

將式(9)和式(10)改寫為矩陣形式:

Yf=-Zfβ+wf

(11)

Yb=-Zbβ*+wb

(12)

式中:Yf=[x(p),x(p+1),…,x(Nx-1)]T,Yb=[x(0),x(1),…,x(Nx-p-1)]T分別為TVAR模型的前向和后向觀測數(shù)據(jù)。wf和wb分別為前、后向TVAR模型的高斯白噪聲。Zf和Zb分別為TVAR模型的前向和后向觀測矩陣

(13)

Zf=[A1,A2,…,Aq]

Zb=[B1,B2,…,Bq]

(14)

(15)

2 基于聚類先驗求解BS-FBTVAR模型

本節(jié)利用文獻[16]中改進的EBSBL算法以及剛體目標(biāo)塊邊界已知的先驗條件來求解BS-FBTVAR模型。

2.1 β的聚類結(jié)構(gòu)先驗

EBSBL模型[17]用來解決未知塊邊界的塊稀疏問題,該模型引入了一個擴展的隱藏塊集zi,zi具有固定的塊大小h,i=1∶g,塊的數(shù)量g=Nβ-h+1,β的長度Nβ=pq，以及一個線性變換,將β表示為如下形式:

(16)

式中:Ei∈RNβ×h是一個除了第i行到第i+h-1行的部分被單位矩陣替代的零矩陣。分塊向量z=[z1，z2，…，zg]是具有塊大小h的等分割塊稀疏信號。將式(16)代入式(11)和式(12),得到BS-FBTVAR模型的矩陣形式，如下所示:

(17)

式中:j∈{f,b1},f表示前向TVAR模型中用到的數(shù)據(jù),b1表示后向TVAR模型中用到的數(shù)據(jù)。式(17)是一個具有已知塊劃分的塊稀疏恢復(fù)問題,可由BSBL算法求解。由BSBL算法中假設(shè)的分塊向量z服從p(z;{γi}i,B)=N(0,Σ0),Σ0=diag (γ1B,γ2B,…,γgB)來施加塊稀疏性。每個塊滿足多元高斯分布:

p(zi;γi,B)=N(0,γi,B)

(18)

由于在本文的問題中,塊稀疏系數(shù)向量β是塊邊界已知的,以p為塊大小進行劃分,因此h=p,則zi=[zi(0),zi(1),zi(p-1)]。根據(jù)式(16)對時不變系數(shù)向量β進行計算,則β的第i個元素為當(dāng)前塊zi和相鄰的(p-1)個塊中元素的和,考慮TVAR模型參數(shù)p的取值為奇、偶數(shù)的不同情況,總結(jié)得到

(19)

式(19)與模式耦合先驗公式[18]相似。

2.2 求解BS-FBTVAR模型

利用β在式(18)中的聚類結(jié)構(gòu)先驗,本文的稀疏恢復(fù)算法[24]可以在貝葉斯框架下表述。考慮噪聲為均值為0、方差為γ0的高斯白噪聲,則式(11)和式(12)中信號模型的似然函數(shù)可以寫為

(20)

(21)

(22)

超參數(shù)γ和γ0通過在系數(shù)β上邊緣化，然后執(zhí)行最大似然優(yōu)化，從數(shù)據(jù)Yj中進行估計。邊際分布可以由下式給出:

(23)

(24)

通過分析式(24)中的代價函數(shù)，并遵循文獻[26]中的步驟,很容易證明式(24)的所有局部最小值都是作為稀疏解實現(xiàn)的。

為了獲得超參數(shù)γ和γ0的估計,本文執(zhí)行最大期望(expectation-maximization, EM)算法來等效地最大化聯(lián)合分布p(Yj,β;γj,0,γj),可以將似然性和式(16)中的聚類結(jié)構(gòu)先驗組合為

p(Yj,β;γj,0,γj)=p(Yj|β;γj,0)p(β;γj)

(25)

在EM算法中,β可以被視為隱藏變量,可以最大化Eβ|Yj,γj,0,γj{p(Yj,β;γj,0,γj)}來找到γj,0和γj的估計值。其中，Eβ|Yj,γj,0,γj{·}表示關(guān)于β的后驗分布的期望。

將式(24)代入Eβ|Yj,γj,0,γj{p(Yj,β;γj,0,γj)},本文可以通過最大化Eβ|Yj,γj{lnp(β;γj)}等效地估計γj:

Eβ|Yj,γj{ lnp(β;γj)}=

(26)

(27)

γj,0的更新形式如下:

(28)

由于前后向預(yù)測都是在恢復(fù)β,需要做如下約束:

(29)

(30)

(31)

將目標(biāo)的時不變系數(shù)向量β估計出來后,通過離散余弦基的線性組合公式將β轉(zhuǎn)換為時變系數(shù)ai(n),利用瞬時功率公式計算得到目標(biāo)的瞬時功率譜,由式(3)給出。

2.3 模型參數(shù)的選取準(zhǔn)則

BS-FBTVAR模型中的模型階次p、基函數(shù)個數(shù)q與雷達回波數(shù)據(jù)的散射點個數(shù)和噪聲環(huán)境有關(guān),p和q的選取需要采用一定的準(zhǔn)則,共包含p×q個未知量。本文使用最小描述長度(minimum description length, MDL)準(zhǔn)則來實現(xiàn)p與q的選擇[29]:

MDL (p,q)=

(32)

3 實驗結(jié)果

本節(jié)分別在電磁仿真工具軟件產(chǎn)生雷達目標(biāo)回波的情況下和使用目標(biāo)實測回波數(shù)據(jù)的情況下對目標(biāo)進行微多普勒分析。

3.1 仿真產(chǎn)生雷達回波數(shù)據(jù)的微多普勒分析結(jié)果

為了驗證所提出的基于聚類先驗的BS-FBTVAR模型求解算法的有效性,本文采用電磁仿真軟件CST Studio Suite 2020產(chǎn)生空間錐體目標(biāo)的雷達回波。目標(biāo)使用錐體模型,錐頂為半球體,錐體材質(zhì)為理想良導(dǎo)體。錐體參數(shù)設(shè)置如下:高度為1.5 m,底面半徑為0.9 m,在距離錐體底面0.5 m處的錐體表面開出等間隔的半球體凹槽,半球體半徑為0.05 m。仿真雷達目標(biāo)回波所用的雷達參數(shù)設(shè)置如下:載頻為10 GHz,帶寬為2 GHz,脈沖重復(fù)頻率(pulse repetition frequency, PRF)為800 Hz,駐留時間為1 s。目標(biāo)的自旋頻率為1 rad/s,進動頻率為1 rad/s,雷達入射方向與錐體中心軸線垂直,模型圖如圖1所示。

圖1 目標(biāo)模型示意圖Fig.1 Schematic diagram of the target model

圖2所示為仿真雷達回波數(shù)據(jù)所對應(yīng)的真實的時變功率譜。其中，兩條正弦曲線對應(yīng)目標(biāo)的頂點A和頂點B,A為點散射中心,B為滑動散射中心。C點只有在雷達視線范圍內(nèi)才能被觀測到,在圖中是不可見的。中間的3條正弦曲線對應(yīng)3個凹點D、E、F,即點散射中心。由于目標(biāo)在旋轉(zhuǎn),因此空間錐體目標(biāo)等間隔的3個凹槽會被雷達視線輪流觀測,在時頻圖中呈現(xiàn)出不完整的曲線。

圖2 真實的微多普勒時頻圖示意圖Fig.2 Schematic diagram of real micro-Doppler time-frequency diagram

在實驗中,將超參數(shù)設(shè)置為a0=1e-2,b0=1e-3,c0=1e-2,d0=1e-3,通過MDL準(zhǔn)則得到p和q。在噪聲較大時,p的取值會大于目標(biāo)含有的散射點的個數(shù)。將收斂閾值設(shè)置為1e -3,對目標(biāo)在無噪聲、20 dB噪聲、10 dB噪聲環(huán)境中分別使用STFT時頻分析方法、基于BSBL求解BS-FBTVAR模型的時頻分析方法、本文算法進行微多普勒分析,得到目標(biāo)回波的時變功率譜估計結(jié)果,設(shè)置功率譜幅度的動態(tài)范圍一致,均為120 dB,如圖3所示。

圖3 仿真目標(biāo)回波在無噪聲、20 dB噪聲以及10 dB噪聲環(huán)境中使用3種時頻分析算法的時變功率譜估計示意圖Fig.3 Schematic diagram of time-varying power spectrum estimation of the use of three time-frequency analysis algorithms in noise-free and noise environment (20 dB and 10 dB)

從圖2與圖3可以看出,本文算法較為準(zhǔn)確地估計出了目標(biāo)回波的時變功率譜。STFT時頻分析方法整體分辨率較低,對于目標(biāo)模型中的D、E、F點散射中心而言,當(dāng)處于即將被遮擋的狀態(tài)時,回波信號能量較低,產(chǎn)生的微多普勒分辨率更低。對于由散射中心A、B微動導(dǎo)致的微多普勒特征模糊的問題,無法同時分辨較近的兩條交叉曲線。基于BSBL求解BS-FBTVAR模型的時頻分析方法的時頻分辨率精度有所提高,能夠?qū)⒂赡繕?biāo)模型滑動散射中心和點散射中心自旋和進動產(chǎn)生的微多普勒頻率完整地分析出來,但是在對點散射中心D、E、F的微動進行分析時,在正弦曲線的各個周期的上頂點和下頂點處均出現(xiàn)了錯誤分量,時頻圖較為渾濁。本文算法能夠?qū)⒂煽臻g錐體目標(biāo)的滑動散射中心和點散射中心微動產(chǎn)生的微多普勒頻率完整地提取出來,明顯分辨出距離較近的兩條交叉曲線,時頻圖更為清晰,時頻分辨率更高。

由圖3可以看出,在目標(biāo)散射點處于即將被遮擋的狀態(tài)或者低信噪比情況下,STFT時頻分析方法時頻能量較弱。而基于TVAR模型的時頻分析方法時頻能量較強,清晰可見,這是因為STFT方法對回波數(shù)據(jù)進行滑窗處理,只處理短時窗內(nèi)的部分回波數(shù)據(jù),而基于TVAR模型的時頻分析方法對全部回波數(shù)據(jù)進行處理。

采用Renyi熵[30]對3種時頻分析算法的時頻分辨性能進行進一步量化分析。Renyi熵是描述系統(tǒng)混亂度的一種衡量指標(biāo),熵值越低，代表時頻譜圖的能量聚集性越好,其計算公式為

(33)

式中：P(t,f)為瞬時功率譜值;α為Renyi熵的次序。本文使用3次Renyi熵對時頻聚集性進行評價,對噪聲環(huán)境為0～20 dB的目標(biāo)回波分別使用3種時頻分析方法計算時頻譜圖的Renyi熵,取100次實驗的平均結(jié)果如圖4所示。由圖4可以看出,基于BSBL求解BS-FBTVAR模型算法與本文算法的時頻能量聚集性都明顯高于STFT時頻分析方法,本文算法的Renyi熵值在0～5 dB噪聲環(huán)境下與基于BSBL的BS-FBTVAR的時頻分析算法基本一致,在5～20 dB噪聲情況下最低,時頻聚集性更好,抗噪性能良好。

圖4 目標(biāo)回波的3種時頻譜圖的Renyi熵值Fig.4 Renyi entropy values of three time-frequency spectrograms of target echoes

3.2 利用實測數(shù)據(jù)對目標(biāo)進行微多普勒分析

雷達參數(shù)設(shè)置如下:載頻為15 GHz,帶寬為2 GHz,PRF為200 Hz。目標(biāo)使用光滑錐體,高為0.85 m,底面半徑為0.225 m,中部和底部分別貼有2個半徑為0.02 m的半球體,兩個半球體對稱分布,如圖5所示。

圖5 目標(biāo)示意圖Fig.5 Target schematic diagram

錐體目標(biāo)自旋頻率為1.5 Hz,雷達發(fā)射電磁波照射錐體目標(biāo)并接收其回波,對其進行解線調(diào)、快時間傅里葉變換以及去除剩余視頻項和包絡(luò)斜置項處理,得到目標(biāo)的高分辨距離像(high resolution range profile, HRRP)數(shù)據(jù)。取HRRP數(shù)據(jù)中錐頂處、中部貼有半球體處以及錐體底部距離單元內(nèi)的數(shù)據(jù),分別使用STFT、BSBL求解BS-FBTVAR模型算法和本文算法，進行微多普勒分析,得到目標(biāo)回波的時變功率譜估計,設(shè)置功率譜幅度的動態(tài)范圍一致,均為120 dB,如圖6所示。圖6(a)～圖6(c)為對目標(biāo)的錐頂距離單元進行微多普勒分析的結(jié)果。該距離單元含有1個點散射中心,在目標(biāo)進行自旋時,此散射中心不會產(chǎn)生微多普勒頻率,時頻圖為零頻。對比3種算法可以看出，STFT算法和基于BSBL求解BS-FBTVAR模型算法都在零頻附近的較大范圍內(nèi)出現(xiàn)了頻譜能量,而本文算法在零頻的能量聚集較為集中,時頻分辨率最高。圖6(b)～圖6(d)為對目標(biāo)的中部貼有半球體處的距離單元進行微多普勒分析的結(jié)果。該距離單元包含2個對稱分布的點散射中心,自旋產(chǎn)生的微多普勒為兩條正弦曲線。圖6(g)～圖6(i)為對目標(biāo)的底部距離單元進行微多普勒分析的結(jié)果,包含1個滑動散射中心和2個點散射中心,2個點散射中心產(chǎn)生的微多普勒為兩條正弦曲線,滑動散射中心不會產(chǎn)生微多普勒頻率,在時頻圖中表現(xiàn)為零頻[31]。

圖6 目標(biāo)實測回波數(shù)據(jù)的時變功率譜估計示意圖Fig.6 Schematic diagram of time-varying power spectrum estimation of target measured echo data

由圖6可以看出,在對目標(biāo)的實測回波數(shù)據(jù)進行微多普勒分析時,本文算法能夠準(zhǔn)確估計出目標(biāo)微動產(chǎn)生的微多普勒頻率,并且較STFT算法和基于BSBL求解BS-FBTVAR模型的時頻分析算法，時頻圖清晰度更高。

由Renyi熵值也可以表明本文算法的優(yōu)越性,表1是對實測回波數(shù)據(jù)3種算法的熵值對比。由表1可以看出,本文算法Renyi熵值最低,時頻聚集性最好。

表1 實測數(shù)據(jù)3種算法Renyi熵定量分析對比

4 結(jié) 論

本文提出了一種新的求解BS-FBTVAR模型的方法,進行目標(biāo)的微多普勒分析,利用改進EBSBL算法,通過人工構(gòu)造增廣向量模型的潛在關(guān)系,得到一個結(jié)構(gòu)化稀疏先驗,然后考慮其超參數(shù)和其相鄰超參數(shù)之間的相互作用,得到聚類結(jié)構(gòu)先驗,并施加剛體目標(biāo)塊邊界已知的先驗信息來求解BS-FBTVAR模型,從而對目標(biāo)進行時變功率譜估計。實驗表明所提算法性能優(yōu)于傳統(tǒng)的時頻分析算法,給出的時頻圖時頻分辨率更高,多普勒曲線的聚集性也更好,并且抗噪聲性能良好。