進動彈道目標平動補償與分離

韓立珣，田波，馮存前，2，*，賀思三

（1.空軍工程大學 防空反導學院，西安710051； 2.信息感知技術協同創新中心，西安710077）

彈道導彈中段防御一直是國際政治和軍事領域關注的重點。彈道目標識別是一種典型的非合作式目標識別，目標的形狀、結構、表面材料電磁參數和常規運動特性等特征對先驗信息要求較高，而攻擊方彈道導彈參數通常很難獲取，這些特征在彈道導彈目標識別中的實用性受到限制，而且彈道目標通常會在中段釋放子彈頭誘餌，此時傳統的目標識別方式難以從中分辨出真彈頭，而基于難以模仿的目標固有屬性微動特性的應用為其提供了新的突破口[1-3]。

目前，國內外在彈道目標微動特性方面已經做了大量的研究，其中，文獻[4]得出多點和噪聲對微多普勒提取影響不大這一結論。文獻[5]提出可以將中段彈道目標的平動近似為多項式表達，利用最小二乘法估計出了平動參數，但該方法只適用于單散射源或目標含有一個強散射源和若干個弱散射源的情況。文獻[6]利用形態學中的骨架提取方法得到了清晰的微多普勒曲線。文獻[7]利用Viterbi算法構造時頻濾波器從而提取出時頻圖中各個信號分量，但該方法只能用于時頻曲線交疊不明顯的情況，不具備普遍適用性。

本文針對以上問題，提出了一種利用圖像處理領域中的角點檢測進行平動補償的方法，可以對含有多散射點的目標實現較好效果的平動補償。然后利用基于交點信息的分段Viterbi算法對回波信號實現了分離處理，較好地解決了傳統Viterbi算法的頻率跳變及錯誤關聯問題。

1 進動錐體彈道目標微動特性

錐體彈道目標的進動模型如圖1所示。以彈頭的進動軸為z軸，彈頭質心為坐標原點O，進動軸與彈頭自旋軸為yOz平面，x軸滿足右手螺旋定理。設彈頭質心O距彈頭頂端A的距離為h1，距彈頭底部圓心的距離為h2，彈頭總高度設為h，彈頭底面圓的半徑為r，半錐角為ε。彈頭進動軸與自旋軸夾角為θ，初始時刻的方位角為φ0，進動角速度為ωc，自旋角速度為ωs。雷達視線（LOS）入射方向的俯仰角α，方位角為φ，與自旋軸的夾角為 β。

圖1 錐體彈道目標進動模型Fig.1 Cone ballistic target's precession model

理論計算和測量試驗均表明，錐體彈道目標的散射中心通常表現為入射線和目標對稱軸所構成的平面與目標不連續處邊緣的交點[8]。由此可知，當雷達波束照射到處于進動狀態的錐體彈道目標時，會在錐頂A，錐底邊緣與電磁波入射平面的交點B、C處形成散射中心。

在不考慮遮擋效應的情況下，文獻[8]分別對3個散射點進行分析，可以得到在t時刻各個散射點的微多普勒公式分別為

式中：F（t）=sinθsinαsin（ωct+φ0）+cosθcosα，λ為信號波長；ω為角速速。

可以看出，對于理想散射點A，該點的運動僅僅表現為自旋，運動方式較為簡單，符合正弦調制規律；對于滑動散射點B、C，這兩點的運動表現為自旋和錐旋的合成，運動形式比較復雜，不再滿足簡單的正弦調制規律。

此時假設雷達發射波長為λ的信號，則接收到的雷達基頻回波信號可以表示為

式中：Rk（t）為t時刻k散射點與雷達的徑向距離。

物體處于高速平動時會出現多普勒折疊現象，彈道目標在中段飛行時，彈道是相對平穩的，此時雷達會處于“粗跟”狀態。可以利用某個觀測時間內的脈沖測得的速度vi對回波信號進行重構補償，從而完成對速度的粗補償[9]。粗補償后的徑向距離可以近似表示為

式中：r0、v0、a和rk（t）分別為t時刻散射點k的徑向初始距離、速度、加速度和微動距離。

在實際雷達回波中要考慮遮擋效應帶來的影響，通過分析可以得出各個散射點的可見條件如表1所示。

表1 散射點可見角度Table 1 Visible angle of scattering points

可以設一個條件函數 γ（t），當滿足表1時，γ（t）的函數值為1，其余時刻為0。

綜上所述，基頻回波信號可以表示為

基頻信號相位的一階導數即為頻率信息，因此可以推導出進動錐體彈道目標的頻率為

式中：ft和fi分別為目標平動與微動導致的頻率。

2 基于角點檢測的平動補償

通過式（1）和式（5）可以看出，當各個散射點的頻率相同時（表現為時頻圖上的曲線交點），應當滿足fi=0這一條件，即時頻圖上曲線的交點所對應的頻率與目標散射點的微動特性無關，完全受平動控制，所以可以通過交點信息反推平動信息，從而進行精度較高的平動補償。

在這里引進圖像處理領域中的角點檢測[10]算法提取時頻圖的交點信息。

現有的角點檢測算法一般分為3類[11]：基于結構邊緣輪廓的角點檢測[12]；基于圖像灰度強度的角點檢測[13]；基于模板匹配的角點檢測。常見的角點檢測算法有Harris算法、Harris-Laplace算法、He&Yung[14]算法等。其中基于Harris算法的角點檢測方法被公認為是效果較好的角點檢測方法之一。但Harris算法存在對尺度變換敏感且提取出的角點是像素級別的缺陷，Harris-Laplace算法存在極值、定位精度以及冗余檢測等問題。如將上述方法直接運用到時頻圖的交點檢測中，檢測的效果不是很理想。

可以分別在以下2個方面加以改進。

1）對時頻圖進行預處理。采用數學形態學圖像處理技術對時頻圖進行加工。通過這一處理可以較好地消除噪聲的影響、提高分析精度。

2）采用改進的Harris角點檢測方法。針對Harris角點檢測方法的不足，采用一種基于雙邊濾波器的Harris角點檢測方法。具體流程如下：

從而建立雙邊濾波函數w（x，y）為

式中：u為歸一化常數；σk和 σh分別為空間距離標準差和灰度距離標準差。

步驟2 利用雙邊濾波函數得到自相關函數矩陣M。

步驟3 計算角點檢測算子R。

式中：det M 表示M 矩陣的行列式；tr M 表示矩陣M 的跡；q為一常數，這里取0.06。當R值為像素點在3×3范圍內最大值且大于設定的閥值時即可認為該點為角點。

在檢測出角點后，提取出角點的坐標（x，y），數值x即為時頻圖交點的時間點取值，數值y即為時頻圖交點的頻率取值。原則上只需2個交點信息即可求解出平動項參數v0、a的信息，在這里為保證精度減少誤差，可以將交點兩兩組合求解出一組平動項參數信息再取平均即可估計出平動參數。

在利用時頻圖交點求解出平動參數即v′0、a′后，即可利用平動補償函數s′（t）對原來的回波進行平動補償。

3 基于交點信息的曲線分離

Viterbi算法[15]是一種以信號能量大小為依據，對多分量信號進行抽取的算法。這種算法假設瞬時頻率曲線是一條相對平滑的曲線，對瞬時頻率估計就是為了使下式的估計路徑最小化：

式中：K為時頻分布中所有可能路徑的集合；N為采樣點個數；u（x）為定義在GD（n，k（n））函數上的代 價 函 數，表 征 點 的 權 值；g（x，y）表 示的代價函數，在這里依據下式取值：

式中：ξ為相鄰2個瞬時頻率點頻率變化的期望值，一般取決于時頻變換的頻率分辨率[16]。

在經過平動補償后，時頻曲線將會聚焦在零頻附近，此時如果采用傳統的Viterbi算法進行曲線分離，將會在時頻曲線交點處產生錯誤關聯從而無法精確分離各個曲線[17]，而且往往會因為交點處的復雜情況導致計算時間大大增加。在這里提出一種可以充分利用之前所提取出交點信息的分段Viterbi算法對交叉程度較高的彈道目標回波時頻圖進行分離，具體步驟如下：

步驟1 根據交點信息將曲線分段。

步驟2 利用Viterbi算法對每段圖像分別進行曲線分離。

步驟3 分別計算每段圖像交點附近不同曲線的斜率，利用斜率大小對曲線進行編號。

步驟4 將相鄰兩段圖像編號相同的曲線合并，從而實現整體曲線的分離。

4 仿真實驗

進行以下仿真實驗驗證本文算法的有效性。

仿真參數設置：設雷達發射載頻為6 GHz，脈沖重復頻率為500 Hz，帶寬為5 MHz的單頻信號，觀測時間為2 s，信噪比為2 d B。錐體彈道目標的錐體高度h1=2 m、h2=0.5 m，底面半徑r=0.5 m，進動角 θ=10°，進動角速度 ωc=4πrad/s，雷達視線入射方向的俯仰角α=，與自旋軸的夾角β=。經過平動預補償后的v=-2 m/s、a=2.5 m/s2。

圖2仿真了在經過平動預補償后，錐體彈道目標的3個散射中心造成的回波多普勒曲線，可以看出在觀測時間內3條曲線共有8個共同的交點，且3條曲線隨時間作同方向的近似線性傾斜運動。

隨后對圖2中的時頻曲線進行加工。建立一個25×25像素的高斯空間掩模對圖像進行平滑處理，再進一步將其轉化為二值化圖像，然后提取圖像骨架得到圖3。

圖2 錐體彈道目標回波時頻圖Fig.2 Time-frequency image of cone ballistic target echo

圖3 時頻曲線骨架Fig.3 Time-frequency curve skeleton

對于圖3可以采用本文提及的角點檢測算法得到角點坐標。因角點檢測算法是對整個圖像進行檢測，圖像的坐標軸和單位等部分被當成檢驗對象后，會出現大量位于坐標軸上的多余角點并增加檢驗時間，從而影響整體檢測的結果，故這里在角點檢測前先將圖像的坐標隱去，只對所需要的時頻曲線進行檢測。

綜上，分析兩罪犯罪構成要件的不同之處，我們很容易將二者區分，并清晰探知嫖宿幼女罪在刑法分則中所處位置的意義，它與強奸罪有重合部分，但又各司其職，屬于特別法與一般法的關系。

將檢測到的角點繪制在圖4中。可以看出此時檢測出12個角點，其中第2、5、6、7個曲線交點處有多個角點聚集，第3個曲線交點未檢測出角點。由這12個曲線交點坐標估計出平動參數v=-1.958 8 m/s，a=2.460 5 m/s2。

利用估計出的平動參數結合式（11）對回波信號時頻圖進行平動補償可以得到圖5。從圖5中可以看出，平動補償很好地消除了平動項對時頻圖的影響，此時時頻曲線近似分布在零頻附近，從而驗證了本文算法的有效性。

再在圖5的基礎上對時頻曲線進行加工并提取出角點，可以得到圖6。從圖6中可以看出，總共檢測出9個角點，這9個角點其中的8個是曲線的交點，有一個是曲線的最高點。利用這9個交點坐標對圖5進行分段處理，將時頻曲線分為9段可以得到圖7。再利用本文提出的分段Viterbi算法對每段圖像進行處理，對每一段的時頻曲線進行分離，結果如圖8所示。

圖4 角點檢測圖Fig.4 Corner detection image

圖5 平動補償后的時頻圖Fig.5 Time-frequency image after translation compensation

圖6 平動補償后的角點檢測圖Fig.6 Corner detection image after translation compensation

將相鄰兩段圖像斜率相同的曲線用相同顏色標識，然后將相同顏色的曲線合并即可得到完整的分離曲線，如圖9所示。

圖10為采用文獻[7]中的Viterbi算法直接對時頻圖中曲線進行分離的結果。

從圖9、圖10的對比中可以看出，本文提出的分段Viterbi算法較傳統的Viterbi算法在交點處分離曲線的準確度有很大提升。原因是本文提出的分段Viterbi算法利用時頻曲線中的交點對圖像進行處理，充分利用了圖像本身的特性，將復雜的處理過程簡潔化，從而避免了Viterbi算法帶來的頻率跳變及錯誤關聯問題。

為進一步分析本文提出算法在不同信噪比情況下的分離效果，接下來在不同信噪比條件下各進行100次的蒙特卡羅仿真，定義歸一化均方根誤差為

圖7 分段處理效果圖Fig.7 Periods processing effect image

圖8 分段Viterbi算法分離處理Fig.8 Separation processing by segmentation Viterbi algorithm

圖9 分段Viterbi算法抽取的曲線Fig.9 Curves extracted by segmentation Viterbi algorithm

圖10 Viterbi算法抽取的曲線Fig.10 Curves extracted by Viterbi algorithm

圖11 本文算法與Viterbi算法對比Fig.11 Comparison between proposed algorithm and Viterbi algorithm

隨著信噪比的增加，歸一化均方根誤差整體呈下降趨勢，且在信噪比大于6 d B后，曲線趨于平穩，這是由于骨架提取的固有缺陷造成的。從圖中可以明顯看出，本文所提算法的歸一化均方根誤差要小于Viterbi算法。

5 結 論

本文通過分析進動錐體彈道目標微動特性，提出了一種利用時頻曲線交點信息進行平動補償與曲線分離的方法。

1）彈道目標的時頻圖是多個散射點微動分量疊加的結果，想要直接進行平動補償往往很困難。利用交點處微多普勒與散射點微動特性完全受平動效應支配這一結論，本文先采用改進的Harris角點檢測方法提取出時頻曲線交點坐標，再根據交點進行平動補償，仿真實驗表明該方法具有良好的平動補償效果。

2）針對傳統Viterbi算法在交點處出現的頻率跳變及錯誤關聯現象，本文充分利用交點信息，先對時頻圖進行分段處理，再利用Viterbi算法對每一段時頻圖進行分離處理，最后再進行曲線合并得到完整的分離曲線。仿真實驗表明該理念的正確性。

3）對算法進行了蒙特卡羅仿真實驗，仿真結果表明在不同信噪比情況下本文算法的歸一化均方根誤差均小于Viterbi算法，具有良好的抗噪性與較高的準確度。

4）在仿真中發現，由于骨架提取的固有缺陷，通過本文算法分離出的曲線與原始數據始終有一定的誤差，如何減少這一誤差，提高曲線分離準確度將是下一步研究的重點。