韋楠楠 張興敢

（南京大學(xué)電子科學(xué)與工程學(xué)院，江蘇南京 210023）

1 引言

彈道導(dǎo)彈的飛行過程可分為助推段、中段和再入段，中段是彈道導(dǎo)彈飛行時間最長的一段，被認為是防御的關(guān)鍵階段。構(gòu)成彈道中段目標的物體有彈頭、誘餌、彈體碎片、金屬箔條、干擾機等，這些目標的形狀較為簡單，彈頭和誘餌一般是錐、柱、錐臺、球冠的組合體，碎片一般是不規(guī)則薄片；這些目標的運動形式也較為特定，是平動與微動（進動、章動、翻滾等）的復(fù)合運動。彈頭及誘餌等目標由于受到質(zhì)量及質(zhì)量分布等因素影響，它們的進動頻率、章動角等運動特性存在明顯差異，因此常被用于彈道目標識別的研究中。

隨著雷達技術(shù)的發(fā)展，彈道目標的運動細節(jié)能夠更精細地刻畫出來，由于微動能夠?qū)﹄姶挪ㄟM行額外的調(diào)制，導(dǎo)致多普勒頻移現(xiàn)象的產(chǎn)生，其表現(xiàn)在時頻圖上，則是頻率周期性的變化，許多學(xué)者［1-3］對此問題進行了廣泛研究。國內(nèi)外對彈道目標的微動特性研究主要有兩種思路：一是研究學(xué)者們［3-6］利用窄帶回波信號獲得的RCS 序列或時頻圖提取出目標微動參數(shù)，但是目標多散射點的回波信號在一個波束中會相互交疊，從而產(chǎn)生微多普勒折疊的現(xiàn)象。另一種研究思路則是通過分析寬帶雷達的一維距離像（high range resolution profile，HRRP）序列的調(diào)制特征來提取微動參數(shù)，文獻［7］利用相鄰HRRP之間的差分值形成一維差分序列，再估計出目標進動頻率，文獻［8］利用逆Radon變換和偽Zernike矩方法提取出HRRP序列的微動特征，文獻［9］通過正弦曲線擬合的方法提取出目標散射中心的微動參數(shù)。寬帶雷達的分辨率較高，可以在距離維上分辨出目標對應(yīng)的多個散射點，同時還能解決多分量微多普勒信號混疊問題。無論是利用窄帶信號或是寬帶信號提取微動特征，由于誘餌與彈頭的相似度越來越高，并且反導(dǎo)作戰(zhàn)對系統(tǒng)實時性要求越來越高，而微動特征提取與目標成像往往依賴于長時間觀測，因此如何在數(shù)據(jù)率、積累時長、帶寬等條件的限制下提取出穩(wěn)定的、可分性高的微動特征是需要考慮的問題。

本文的研究則是通過寬帶雷達獲得的HRRP序列提取目標進動頻率。本文以錐柱狀目標作為分析對象，首先通過電磁仿真獲得錐柱狀目標的轉(zhuǎn)臺特性，分析其多個散射點在半空間范圍內(nèi)的可見范圍；其次，通過構(gòu)建錐柱狀目標的進動模型分析HRRP 序列周期性變化的過程；建立錐柱體彈頭滑動型散射中心的位置模型，在此基礎(chǔ)上推導(dǎo)滑動型散射中心的微多普勒的數(shù)學(xué)表達式?；谝陨戏治?，本文提出了自適應(yīng)閾值法獲取目標進動頻率的算法，該算法提取HRRP的強散射中心單元，對該單元進行時頻處理獲得時頻圖，通過自適應(yīng)閾值法提取出時頻圖多個散射點的時頻脊線，據(jù)此估計目標微動頻率。本文實驗采用了衛(wèi)星工具箱（STK）和電磁仿真工具（FEKO）聯(lián)合仿真，生成彈道導(dǎo)彈寬帶雷達回波數(shù)據(jù)，通過蒙特卡洛實驗和對比實驗證明了算法具有較好的穩(wěn)定性和抗噪性。

2 彈道目標進動模型及散射特性

2.1 錐柱體目標的散射特性

旋轉(zhuǎn)體目標的散射場，主要包括尖頂散射、錐面行波、曲面不連續(xù)處散射、柱面行波和角頂散射，這五類散射中，行波是較弱的一類散射場。雷達獲得的回波信號可以描述為各散射中心的回波合成，當雷達發(fā)射電磁波照射在目標上，目標旋轉(zhuǎn)一定角度后，如果電磁波照亮部分的形狀、介質(zhì)等影響回波因素未發(fā)生變化，那么目標的散射場也基本不變。

建立錐柱狀彈頭模型如圖1 所示，雷達發(fā)射帶寬為8.5 GHz～9.5 GHz 的步進頻信號，通過FEKO計算得到的目標一維距離像在0°～180°范圍內(nèi)的變化情況如圖2 所示。五個散射中心的可見范圍如下：散射中心A為0～(π -γ)，散射中心B、C在(0～π)范圍內(nèi)均可見，散射中心D 在范圍內(nèi)可見，在散射中心E的可見范圍為(0～γ)。通常防御雷達位于彈頭的前下方，雷達視線和目標主軸的夾角為50°～80°。

2.2 錐柱體目標的進動模型分析

建立錐柱體彈頭目標的進動模型如圖3 所示，假設(shè)以目標質(zhì)心O 點為坐標原點建立參考坐標系OXYZ。目標圍繞自身對稱軸OA 旋轉(zhuǎn)的運動稱為自旋，自旋角速度設(shè)為ωZ；目標繞著某一旋轉(zhuǎn)軸（非幾何對稱軸）旋轉(zhuǎn)的運動稱為錐旋，錐旋軸設(shè)為OZ，目標錐旋角速度設(shè)為ωp；OZ 與OA 的夾角θ，稱為進動角。假設(shè)雷達視線（radar light of sight，RLOS）方向位于YOZ 平面內(nèi)，與OZ 的夾角α稱為雷達視線角；RLOS 與OA 軸的夾角ρ稱為姿態(tài)角。假設(shè)在t=0時刻，自旋軸在YOZ平面的初始方位角為φ0，則姿態(tài)角ρ隨時間的變化公式為

式（1）中雷達視線角α?xí)S著彈道目標飛行時間的推移而改變，彈道導(dǎo)彈中段的飛行時間一般為十幾分鐘，而進動周期為秒級，α在一個進動周期內(nèi)的變化幾乎可以忽略，所以在求進動參數(shù)時α可看作固定值。進動角θ對ρ(t)的影響與α等效。式（1）表明了目標進動的周期性會導(dǎo)致姿態(tài)角ρ周期性變化。對于旋轉(zhuǎn)對稱體而言，自旋不會導(dǎo)致散射場變化，但錐旋和雷達視線的變化會影響散射場，所以目標周期性進動會導(dǎo)致一維距離像上的散射中心位置和幅度發(fā)生周期性變化。

2.3 錐柱體目標散射中心的微多普勒模型

假設(shè)錐柱體目標平動經(jīng)過補償，只考慮進動帶來的頻率調(diào)制。進動目標的散射場變化的影響因素主要是錐旋，因此僅考慮目標錐旋運動產(chǎn)生的微多普勒特性。目標模型如圖3 所示，假設(shè)目標質(zhì)心點O 到錐柱結(jié)合的平面的距離為d1，與柱體底面的距離為d2，柱體底面半徑為r。目標對稱軸OA 在XOY 平面內(nèi)的投影為OV，OV 與OX 的夾角為φ，在彈體上取一個固連其上的坐標系OVWZ，OV、OW、OZ 構(gòu)成右手系。記參考坐標系OXYZ 的單位向量基為[eX，eY，eZ]，連體坐標系OVWZ 的單位向量基為[eV，eW，eZ]，則目標體錐旋角速度向量可定義為：ω=φ?eZ=ωp?eZ，其中φ=φ0+ωpt，φ0為OA 的初始位置在XOY 平面內(nèi)的投影與OX 軸的夾角，由于φ0僅改變微多普勒的初相，不影響微多普勒的周期及幅度，所以令φ0=0，則φ=ωpt。

錐柱體目標的散射中心可分為錐頂局部散射中心（圖1中的點A）、邊緣（棱線）型散射中心（圖1中的B、C、D、E）。根據(jù)電磁散射理論，彈頭頂部A通常會形成一個較強的散射中心，由于該散射中心位于彈體對稱軸上，其運動規(guī)律和彈頭運動規(guī)律相同，即以角速度ωp繞OZ做圓周運動。根據(jù)文獻［10］的分析，錐頂散射點A 在OVWZ 坐標系的矢量rA=(sinθ?eV+cosθ?eW)?LA，LA為OA長度，則點A的速度可表示為

雷達視線方向為RLOS=[0，sinα，cosα]T，雷達發(fā)射電磁波的波長為λ，那么錐頂散射中心A 在雷達視線方向上引起的微多普勒可表示為：

由式（3）可以看出，錐頂散射中心A 的微多普勒表達式是正弦形式。

邊緣（棱線）型散射中心的微多普勒，除了受彈體微運動的影響外，還和雷達視線與錐軸所構(gòu)成的平面有關(guān)，當旋轉(zhuǎn)對稱彈頭在自由空間進動時，散射中心B、C、D、E 的微運動規(guī)律和彈頭本身的進動規(guī)律并不完全一致，因此，該處散射中心的微多普勒將不再是簡單的正弦形式。以散射中心B、E 為例，根據(jù)文獻［11-12］的分析，散射中心B、E 在雷達視線方向上引起的微多普勒可表示為：

其中：F(t) ?sinθsinαsinωpt+cosθcosα，i=B，E。從式（4）可以看出，滑動型散射中心的微動引起的微多普勒在余弦形式之外，增加了的調(diào)制項。對式（4）進行泰勒展開，并忽略較小的高階項，可以得到散射中心B、E的微多普勒近似表達式為

將式（5）化簡為式（6）形式

底部邊緣形成的散射中心C、D 的微多普勒形式也如上式（4）所給出，只是d1替換為坐標原點O 到底部所在平面的距離d2，同時只取加號即可。從式（5）、（6）來看，滑動型散射模型的微多普勒呈現(xiàn)顯著的非單頻性，但其微多普勒曲線表現(xiàn)出周期性，這種周期性與目標的錐旋頻率相對應(yīng)。滑動型散射中心的微多普勒在錐旋頻率的基礎(chǔ)上增加了倍頻分量，但并未改變微多普勒的周期性。

以圖1目標模型進行仿真分析，利用式（1）生成不同雷達視線角下的目標姿態(tài)角數(shù)據(jù)，并結(jié)合目標電磁仿真數(shù)據(jù)，產(chǎn)生不同視線角下的目標HRRP 序列。設(shè)目標進動頻率為1.5 Hz，進動角度為10°，雷達工作參數(shù)同2.1，脈沖重復(fù)頻率為500 Hz，觀察時間設(shè)為2 s。時頻分析方法采用100 點的漢明窗和128點的短時傅里葉變化處理。

圖4、圖5給出了目標分別在雷達視線角為50°、80°下的微多普勒時頻分析結(jié)果。圖4（a）的HRRP序列顯示50°視線角下目標散射點出現(xiàn)越距離單元走動，而圖5（a）顯示的80°視線角下的目標散射點出現(xiàn)重合，這一現(xiàn)象與圖2所顯示的半空間HRRP散射點變化一致，由于目標進動的周期性，一維距離像的幅度與能量強弱產(chǎn)生周期性變化。各自取不同距離單元的回波進行時頻分析，當雷達視線角為50°時，圖4 顯示目標回波在第32、33 個距離單元的時頻圖可以觀察到區(qū)分明顯的時頻脊線，而在第34、35、36個距離單元的時頻圖出現(xiàn)了不同程度的交疊與重合；圖5 顯示目標在80°的雷達視線角下，其回波第32個距離單元的時頻圖出現(xiàn)嚴重交疊現(xiàn)象。

旋轉(zhuǎn)體目標的進動導(dǎo)致各散射點的微動變化曲線為正弦函數(shù)或正弦函數(shù)的疊加形式，而在某些視角下，可能會出現(xiàn)散射點無法分離、時頻脊線交疊重合的現(xiàn)象，但我們發(fā)現(xiàn)時頻脊線的振幅與能量強弱均具有周期性，該周期性與進動頻率相對應(yīng)。在本文的研究中，試圖從圖像處理技術(shù)的角度從時頻圖中提取出不同散射點的時頻脊線，從而估計出目標的進動頻率。

3 基于HRRP序列的進動頻率估計方法

3.1 提取HRRP的強散射中心單元

假設(shè)雷達回波信號經(jīng)過平動補償，收到Np幀HRRP 序列，其中每幅距離像的距離單元數(shù)Ns。為進行進動頻率的提取，將其構(gòu)建為一維距離像數(shù)據(jù)矩陣：

3.2 微多普勒瞬時頻率估計

進動頻率提取和分離的前提，首先通過對H0r的Nst個強散射中心單元依次進行STFT處理，得到散射點的時頻分布圖，由于HRRP在經(jīng)過平動補償后，時頻曲線會聚焦在零頻附近，在時頻曲線交點處，容易產(chǎn)生錯誤關(guān)聯(lián)而無法分離各個曲線。在這里提出利用自適應(yīng)閾值的方法從時頻圖中提取時頻脊線。

STFT的定義為：

x(t)為待處理信號，g(t)為窗函數(shù)。假設(shè)時頻圖中存在p條時頻脊線，通過遍歷每一時刻的頻譜提取時頻脊線，設(shè)t時刻的頻譜為：

Ei為第i個頻點的幅值，ωi為第i個頻點的頻率值。首先對ft(i)進行歸一化：操作步驟如下：

步驟1首先，設(shè)立能量幅度閾值e為0.15，判斷滿足ft(ωi)＞e條件的頻點值個數(shù)，假設(shè)滿足該條件的頻點值有k個，即ωi1，ωi2，...，ωik；若k＞p，則令e=e+0.01，通過提高閾值，去掉旁瓣；若k=0，則令e=e-0.005；若k＜p，則令ωi[(k+1)：p]=ωik，直至k=p，則退出循環(huán)。此時，在t時刻可提取出p個頻點值。

步驟2提取完所有時刻點的頻譜值后，此時獲得一個p×m維矩陣表示第k個散射點，在第tm時刻提取的頻率值大小。由于目標短時間運動過程中，可認為回波信號是平穩(wěn)的，其時頻變化曲線的幅度大小相對平穩(wěn)。利用此特點，對Wp×m矩陣進行重排，每一列按照頻率幅值由小到大排序，，重排后的矩陣為。

步驟3利用卡爾曼濾波方法對的每一條時頻脊線進行平滑處理。

需要說明的是，步驟1 中的閾值e的設(shè)定沒有嚴格的要求，由于能量幅度已經(jīng)過歸一化處理，因此需要滿足e∈(0，1)，該步驟的最終目的是提取出p個頻點值，e的值越大，能提取到的頻點值越少，e的值越小，則提取到的頻點值越精確。

3.3 進動頻率估計

在2.3 節(jié)的分析中，進動時目標的微多普勒頻率曲線為周期函數(shù)，其頻率分量一定是進動頻率的倍數(shù)，因此這為進動頻率的提取提供了途徑。該曲線的振動頻率等于進動頻率。通過3.2 的三個步驟，提取了散射中心的時頻線Xk(tm)

其中Xk(tm)為第k個散射點的微多普勒曲線。時頻脊線去除直流分量的影響，利用快速傅里葉變換分別對p條時頻脊線曲線進行頻譜分析，

獲得相應(yīng)的頻譜曲線Fk(fn)，fn=f1，...，fN，fN為采樣頻率，取能量較大的頻率分量作為進動頻率估計值

為第k個散射點進動頻率估計值。最后取p條時頻曲線的進動頻率的平均值作為最后估計的進動頻率為：

在3.1 的預(yù)處理中，提取出了Nst個強散射中心單元，對每個強散射中心單元重復(fù)以上步驟，獲得Nst個進動頻率估計值，對其求均值，得到最終的頻率估計值。

4 仿真數(shù)據(jù)的分析與處理

實驗數(shù)據(jù)由仿真獲得，仿真流程如圖6所示，首先采用衛(wèi)星工具箱（STK）生成彈道數(shù)據(jù)，再利用電磁散射仿真工具（FEKO）生成電磁仿真數(shù)據(jù)。

STK軟件負責對戰(zhàn)場環(huán)境、彈頭平動微動、地球自轉(zhuǎn)、引力等因素進行綜合仿真，并生成任意時刻下雷達視線角、飛行距離、速度等數(shù)據(jù)。實驗所用彈道發(fā)射地經(jīng)緯度為(43.32，63.11)，落點經(jīng)緯度為(116.40，39.89)，關(guān)機點位置(82.76，60.52)，取彈道中段的飛行時間進行實驗，觀測時間總長為297 s。在觀測段中雷達視線與導(dǎo)彈運動軌跡切線方向夾角（即雷達視線角）變化如圖7所示。

目標散射模型由FEKO 軟件進行仿真。電磁仿真參數(shù)設(shè)置如下：雷達信號為X波段的步進頻信號，信號載頻起始頻率為8.5 GHz，終止頻率為9.5 GHz，工作頻率步長為15.675 MHz，電磁計算的求解方式為PO算法。目標尺寸如圖1所示，目標運動參數(shù)見表1，三個目標的目標尺寸相同，微動參數(shù)不同。

表1 目標微動參數(shù)Tab.1 Simulation parameters

通過FEKO 軟件計算生成電磁數(shù)據(jù)包含在.ffe文件中，其包括了各入射點和反射點的散射電場E(n)。目標頻域響應(yīng)序列Sr(t，n)與E(n)的對應(yīng)關(guān)系如式（17）所示為：

對序列Sr(t，n)進行逆傅里葉變換即獲得HRRP。

4.1 進動頻率估計

圖8 是雷達前1000 次回波得到的彈頭及誘餌的一維距離像序列。提取彈頭和誘餌的一維距離像序列的強散射中心單元，并應(yīng)用STFT 獲得時頻分布圖如圖9 所示，從圖中可觀察到時頻脊線的幅度強弱存在周期性變化，并且也不是嚴格意義上的正弦曲線。

接下來對時頻圖提取時頻脊線，效果如圖10所示，從提取效果來看，本文所提的數(shù)據(jù)關(guān)聯(lián)的方法，能夠準確分離出時頻脊線。

對獲得的時頻脊線去直流分量后分別進行FFT，如圖11 所示，T1 目標的P1、P2散射點時頻脊線的頻譜最大值分別為1.4648 Hz、1.4648 Hz，T2 目標的P1、P2散射點時頻脊線的FFT 頻率最大值分別為1.9531 Hz、1.9531 Hz。

最后獲得T1 的估計值為1.4726 Hz，T2 的估計值為2.0215 Hz，與真實值誤差較小。

4.2 性能分析

為了驗證本文算法在不同噪聲影響下的穩(wěn)定性，本文對仿真的信號加入高斯白噪聲，信噪比從SNR=-5 以1 dB 的步進遞增到20 dB，不同信噪比下進行100次蒙特卡洛實驗獲得進動參數(shù)的估計誤差。取T2目標進行進動頻率估計實驗，雷達脈沖重復(fù)頻率設(shè)為500 Hz，回波積累時長為2 s。結(jié)果如圖12所示（圖12（b）的均方誤差結(jié)果利用20 log（）處理成dB 單位），從結(jié)果來看，當信噪比大于1 dB 時，進動頻率的估計均方根誤差趨于平穩(wěn)，與文獻［7］提出的序列差分方法進行比較，本文算法在信噪比較低的情況下，依然具有良好的提取效果。

4.3 彈道中段目標進動特征提取

為了驗證本文算法在彈道目標整個觀測段的提取效果，對彈道回波疊加信噪比為12 dB 的高斯白噪聲，回波積累窗長為2 s，窗口步進1 s，共有297 個統(tǒng)計窗口。雷達工作方式不變，PRF 設(shè)為500 Hz，目標參數(shù)不變，對每個統(tǒng)計窗中的HRRP 序列提取目標進動頻率，最后對297個估計值求平均，結(jié)果如表2 所示，三個目標的估計結(jié)果均靠近真實值。

表2 彈道中段目標進動頻率估計結(jié)果Tab.2 Precession frequency estimation of ballistic missiles in midcourse phase

5 結(jié)論

本文通過分析彈道導(dǎo)彈的錐柱體目標的進動模型與散射點模型，利用散射點在HRRP 序列上的周期特性，運用時頻分析方法獲得時頻圖，提取目標的微多普勒頻率，并估計目標進動參數(shù)。仿真結(jié)果表明，該方法利用圖像處理方式能夠在一定程度上消除噪聲的影響，并且能解決不同視角下時頻脊線重合的問題，有效地估計出進動參數(shù)，具有較好的穩(wěn)定性和抗噪性。

中段彈頭目標結(jié)構(gòu)一般較為簡單，散射中心位置要滿足在觀測期間有一個強散射中心的條件。文中實驗只分析了兩個散射中心的情況，實際中，目標可能存在兩個或兩個以上的散射中心，本文所提的方法依然適用，并且散射點對應(yīng)位置關(guān)聯(lián)的方法也相對簡單。