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

（1.空軍工程大學(xué) 防空反導(dǎo)學(xué)院，西安710051； 2.信息感知技術(shù)協(xié)同創(chuàng)新中心，西安710077）

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

目前，國(guó)內(nèi)外在彈道目標(biāo)微動(dòng)特性方面已經(jīng)做了大量的研究，其中，文獻(xiàn)[4]得出多點(diǎn)和噪聲對(duì)微多普勒提取影響不大這一結(jié)論。文獻(xiàn)[5]提出可以將中段彈道目標(biāo)的平動(dòng)近似為多項(xiàng)式表達(dá)，利用最小二乘法估計(jì)出了平動(dòng)參數(shù)，但該方法只適用于單散射源或目標(biāo)含有一個(gè)強(qiáng)散射源和若干個(gè)弱散射源的情況。文獻(xiàn)[6]利用形態(tài)學(xué)中的骨架提取方法得到了清晰的微多普勒曲線。文獻(xiàn)[7]利用Viterbi算法構(gòu)造時(shí)頻濾波器從而提取出時(shí)頻圖中各個(gè)信號(hào)分量，但該方法只能用于時(shí)頻曲線交疊不明顯的情況，不具備普遍適用性。

本文針對(duì)以上問(wèn)題，提出了一種利用圖像處理領(lǐng)域中的角點(diǎn)檢測(cè)進(jìn)行平動(dòng)補(bǔ)償?shù)姆椒?，可以?duì)含有多散射點(diǎn)的目標(biāo)實(shí)現(xiàn)較好效果的平動(dòng)補(bǔ)償。然后利用基于交點(diǎn)信息的分段Viterbi算法對(duì)回波信號(hào)實(shí)現(xiàn)了分離處理，較好地解決了傳統(tǒng)Viterbi算法的頻率跳變及錯(cuò)誤關(guān)聯(lián)問(wèn)題。

1 進(jìn)動(dòng)錐體彈道目標(biāo)微動(dòng)特性

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

圖1 錐體彈道目標(biāo)進(jìn)動(dòng)模型Fig.1 Cone ballistic target's precession model

理論計(jì)算和測(cè)量試驗(yàn)均表明，錐體彈道目標(biāo)的散射中心通常表現(xiàn)為入射線和目標(biāo)對(duì)稱軸所構(gòu)成的平面與目標(biāo)不連續(xù)處邊緣的交點(diǎn)[8]。由此可知，當(dāng)雷達(dá)波束照射到處于進(jìn)動(dòng)狀態(tài)的錐體彈道目標(biāo)時(shí)，會(huì)在錐頂A，錐底邊緣與電磁波入射平面的交點(diǎn)B、C處形成散射中心。

在不考慮遮擋效應(yīng)的情況下，文獻(xiàn)[8]分別對(duì)3個(gè)散射點(diǎn)進(jìn)行分析，可以得到在t時(shí)刻各個(gè)散射點(diǎn)的微多普勒公式分別為

式中：F（t）=sinθsinαsin（ωct+φ0）+cosθcosα，λ為信號(hào)波長(zhǎng)；ω為角速速。

可以看出，對(duì)于理想散射點(diǎn)A，該點(diǎn)的運(yùn)動(dòng)僅僅表現(xiàn)為自旋，運(yùn)動(dòng)方式較為簡(jiǎn)單，符合正弦調(diào)制規(guī)律；對(duì)于滑動(dòng)散射點(diǎn)B、C，這兩點(diǎn)的運(yùn)動(dòng)表現(xiàn)為自旋和錐旋的合成，運(yùn)動(dòng)形式比較復(fù)雜，不再滿足簡(jiǎn)單的正弦調(diào)制規(guī)律。

此時(shí)假設(shè)雷達(dá)發(fā)射波長(zhǎng)為λ的信號(hào)，則接收到的雷達(dá)基頻回波信號(hào)可以表示為

式中：Rk（t）為t時(shí)刻k散射點(diǎn)與雷達(dá)的徑向距離。

物體處于高速平動(dòng)時(shí)會(huì)出現(xiàn)多普勒折疊現(xiàn)象，彈道目標(biāo)在中段飛行時(shí)，彈道是相對(duì)平穩(wěn)的，此時(shí)雷達(dá)會(huì)處于“粗跟”狀態(tài)?？梢岳媚硞€(gè)觀測(cè)時(shí)間內(nèi)的脈沖測(cè)得的速度vi對(duì)回波信號(hào)進(jìn)行重構(gòu)補(bǔ)償，從而完成對(duì)速度的粗補(bǔ)償[9]。粗補(bǔ)償后的徑向距離可以近似表示為

式中：r0、v0、a和rk（t）分別為t時(shí)刻散射點(diǎn)k的徑向初始距離、速度、加速度和微動(dòng)距離。

在實(shí)際雷達(dá)回波中要考慮遮擋效應(yīng)帶來(lái)的影響，通過(guò)分析可以得出各個(gè)散射點(diǎn)的可見條件如表1所示。

表1 散射點(diǎn)可見角度Table 1 Visible angle of scattering points

可以設(shè)一個(gè)條件函數(shù) γ（t），當(dāng)滿足表1時(shí)，γ（t）的函數(shù)值為1，其余時(shí)刻為0。

綜上所述，基頻回波信號(hào)可以表示為

基頻信號(hào)相位的一階導(dǎo)數(shù)即為頻率信息，因此可以推導(dǎo)出進(jìn)動(dòng)錐體彈道目標(biāo)的頻率為

式中：ft和fi分別為目標(biāo)平動(dòng)與微動(dòng)導(dǎo)致的頻率。

2 基于角點(diǎn)檢測(cè)的平動(dòng)補(bǔ)償

通過(guò)式（1）和式（5）可以看出，當(dāng)各個(gè)散射點(diǎn)的頻率相同時(shí)（表現(xiàn)為時(shí)頻圖上的曲線交點(diǎn)），應(yīng)當(dāng)滿足fi=0這一條件，即時(shí)頻圖上曲線的交點(diǎn)所對(duì)應(yīng)的頻率與目標(biāo)散射點(diǎn)的微動(dòng)特性無(wú)關(guān)，完全受平動(dòng)控制，所以可以通過(guò)交點(diǎn)信息反推平動(dòng)信息，從而進(jìn)行精度較高的平動(dòng)補(bǔ)償。

在這里引進(jìn)圖像處理領(lǐng)域中的角點(diǎn)檢測(cè)[10]算法提取時(shí)頻圖的交點(diǎn)信息。

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

可以分別在以下2個(gè)方面加以改進(jìn)。

1）對(duì)時(shí)頻圖進(jìn)行預(yù)處理。采用數(shù)學(xué)形態(tài)學(xué)圖像處理技術(shù)對(duì)時(shí)頻圖進(jìn)行加工。通過(guò)這一處理可以較好地消除噪聲的影響、提高分析精度。

2）采用改進(jìn)的Harris角點(diǎn)檢測(cè)方法。針對(duì)Harris角點(diǎn)檢測(cè)方法的不足，采用一種基于雙邊濾波器的Harris角點(diǎn)檢測(cè)方法。具體流程如下：

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

式中：u為歸一化常數(shù)；σk和 σh分別為空間距離標(biāo)準(zhǔn)差和灰度距離標(biāo)準(zhǔn)差。

步驟2 利用雙邊濾波函數(shù)得到自相關(guān)函數(shù)矩陣M。

步驟3 計(jì)算角點(diǎn)檢測(cè)算子R。

式中：det M 表示M 矩陣的行列式；tr M 表示矩陣M 的跡；q為一常數(shù)，這里取0.06。當(dāng)R值為像素點(diǎn)在3×3范圍內(nèi)最大值且大于設(shè)定的閥值時(shí)即可認(rèn)為該點(diǎn)為角點(diǎn)。

在檢測(cè)出角點(diǎn)后，提取出角點(diǎn)的坐標(biāo)（x，y），數(shù)值x即為時(shí)頻圖交點(diǎn)的時(shí)間點(diǎn)取值，數(shù)值y即為時(shí)頻圖交點(diǎn)的頻率取值。原則上只需2個(gè)交點(diǎn)信息即可求解出平動(dòng)項(xiàng)參數(shù)v0、a的信息，在這里為保證精度減少誤差，可以將交點(diǎn)兩兩組合求解出一組平動(dòng)項(xiàng)參數(shù)信息再取平均即可估計(jì)出平動(dòng)參數(shù)。

在利用時(shí)頻圖交點(diǎn)求解出平動(dòng)參數(shù)即v′0、a′后，即可利用平動(dòng)補(bǔ)償函數(shù)s′（t）對(duì)原來(lái)的回波進(jìn)行平動(dòng)補(bǔ)償。

3 基于交點(diǎn)信息的曲線分離

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

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

式中：ξ為相鄰2個(gè)瞬時(shí)頻率點(diǎn)頻率變化的期望值，一般取決于時(shí)頻變換的頻率分辨率[16]。

在經(jīng)過(guò)平動(dòng)補(bǔ)償后，時(shí)頻曲線將會(huì)聚焦在零頻附近，此時(shí)如果采用傳統(tǒng)的Viterbi算法進(jìn)行曲線分離，將會(huì)在時(shí)頻曲線交點(diǎn)處產(chǎn)生錯(cuò)誤關(guān)聯(lián)從而無(wú)法精確分離各個(gè)曲線[17]，而且往往會(huì)因?yàn)榻稽c(diǎn)處的復(fù)雜情況導(dǎo)致計(jì)算時(shí)間大大增加。在這里提出一種可以充分利用之前所提取出交點(diǎn)信息的分段Viterbi算法對(duì)交叉程度較高的彈道目標(biāo)回波時(shí)頻圖進(jìn)行分離，具體步驟如下：

步驟1 根據(jù)交點(diǎn)信息將曲線分段。

步驟2 利用Viterbi算法對(duì)每段圖像分別進(jìn)行曲線分離。

步驟3 分別計(jì)算每段圖像交點(diǎn)附近不同曲線的斜率，利用斜率大小對(duì)曲線進(jìn)行編號(hào)。

步驟4 將相鄰兩段圖像編號(hào)相同的曲線合并，從而實(shí)現(xiàn)整體曲線的分離。

4 仿真實(shí)驗(yàn)

進(jìn)行以下仿真實(shí)驗(yàn)驗(yàn)證本文算法的有效性。

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

圖2仿真了在經(jīng)過(guò)平動(dòng)預(yù)補(bǔ)償后，錐體彈道目標(biāo)的3個(gè)散射中心造成的回波多普勒曲線，可以看出在觀測(cè)時(shí)間內(nèi)3條曲線共有8個(gè)共同的交點(diǎn)，且3條曲線隨時(shí)間作同方向的近似線性傾斜運(yùn)動(dòng)。

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

圖2 錐體彈道目標(biāo)回波時(shí)頻圖Fig.2 Time-frequency image of cone ballistic target echo

圖3 時(shí)頻曲線骨架Fig.3 Time-frequency curve skeleton

對(duì)于圖3可以采用本文提及的角點(diǎn)檢測(cè)算法得到角點(diǎn)坐標(biāo)。因角點(diǎn)檢測(cè)算法是對(duì)整個(gè)圖像進(jìn)行檢測(cè)，圖像的坐標(biāo)軸和單位等部分被當(dāng)成檢驗(yàn)對(duì)象后，會(huì)出現(xiàn)大量位于坐標(biāo)軸上的多余角點(diǎn)并增加檢驗(yàn)時(shí)間，從而影響整體檢測(cè)的結(jié)果，故這里在角點(diǎn)檢測(cè)前先將圖像的坐標(biāo)隱去，只對(duì)所需要的時(shí)頻曲線進(jìn)行檢測(cè)。

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

將檢測(cè)到的角點(diǎn)繪制在圖4中?？梢钥闯龃藭r(shí)檢測(cè)出12個(gè)角點(diǎn)，其中第2、5、6、7個(gè)曲線交點(diǎn)處有多個(gè)角點(diǎn)聚集，第3個(gè)曲線交點(diǎn)未檢測(cè)出角點(diǎn)。由這12個(gè)曲線交點(diǎn)坐標(biāo)估計(jì)出平動(dòng)參數(shù)v=-1.958 8 m/s，a=2.460 5 m/s2。

利用估計(jì)出的平動(dòng)參數(shù)結(jié)合式（11）對(duì)回波信號(hào)時(shí)頻圖進(jìn)行平動(dòng)補(bǔ)償可以得到圖5。從圖5中可以看出，平動(dòng)補(bǔ)償很好地消除了平動(dòng)項(xiàng)對(duì)時(shí)頻圖的影響，此時(shí)時(shí)頻曲線近似分布在零頻附近，從而驗(yàn)證了本文算法的有效性。

再在圖5的基礎(chǔ)上對(duì)時(shí)頻曲線進(jìn)行加工并提取出角點(diǎn)，可以得到圖6。從圖6中可以看出，總共檢測(cè)出9個(gè)角點(diǎn)，這9個(gè)角點(diǎn)其中的8個(gè)是曲線的交點(diǎn)，有一個(gè)是曲線的最高點(diǎn)。利用這9個(gè)交點(diǎn)坐標(biāo)對(duì)圖5進(jìn)行分段處理，將時(shí)頻曲線分為9段可以得到圖7。再利用本文提出的分段Viterbi算法對(duì)每段圖像進(jìn)行處理，對(duì)每一段的時(shí)頻曲線進(jìn)行分離，結(jié)果如圖8所示。

圖4 角點(diǎn)檢測(cè)圖Fig.4 Corner detection image

圖5 平動(dòng)補(bǔ)償后的時(shí)頻圖Fig.5 Time-frequency image after translation compensation

圖6 平動(dòng)補(bǔ)償后的角點(diǎn)檢測(cè)圖Fig.6 Corner detection image after translation compensation

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

圖10為采用文獻(xiàn)[7]中的Viterbi算法直接對(duì)時(shí)頻圖中曲線進(jìn)行分離的結(jié)果。

從圖9、圖10的對(duì)比中可以看出，本文提出的分段Viterbi算法較傳統(tǒng)的Viterbi算法在交點(diǎn)處分離曲線的準(zhǔn)確度有很大提升。原因是本文提出的分段Viterbi算法利用時(shí)頻曲線中的交點(diǎn)對(duì)圖像進(jìn)行處理，充分利用了圖像本身的特性，將復(fù)雜的處理過(guò)程簡(jiǎn)潔化，從而避免了Viterbi算法帶來(lái)的頻率跳變及錯(cuò)誤關(guān)聯(lián)問(wèn)題。

為進(jìn)一步分析本文提出算法在不同信噪比情況下的分離效果，接下來(lái)在不同信噪比條件下各進(jìn)行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算法對(duì)比Fig.11 Comparison between proposed algorithm and Viterbi algorithm

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

5 結(jié) 論

本文通過(guò)分析進(jìn)動(dòng)錐體彈道目標(biāo)微動(dòng)特性，提出了一種利用時(shí)頻曲線交點(diǎn)信息進(jìn)行平動(dòng)補(bǔ)償與曲線分離的方法。

1）彈道目標(biāo)的時(shí)頻圖是多個(gè)散射點(diǎn)微動(dòng)分量疊加的結(jié)果，想要直接進(jìn)行平動(dòng)補(bǔ)償往往很困難。利用交點(diǎn)處微多普勒與散射點(diǎn)微動(dòng)特性完全受平動(dòng)效應(yīng)支配這一結(jié)論，本文先采用改進(jìn)的Harris角點(diǎn)檢測(cè)方法提取出時(shí)頻曲線交點(diǎn)坐標(biāo)，再根據(jù)交點(diǎn)進(jìn)行平動(dòng)補(bǔ)償，仿真實(shí)驗(yàn)表明該方法具有良好的平動(dòng)補(bǔ)償效果。

2）針對(duì)傳統(tǒng)Viterbi算法在交點(diǎn)處出現(xiàn)的頻率跳變及錯(cuò)誤關(guān)聯(lián)現(xiàn)象，本文充分利用交點(diǎn)信息，先對(duì)時(shí)頻圖進(jìn)行分段處理，再利用Viterbi算法對(duì)每一段時(shí)頻圖進(jìn)行分離處理，最后再進(jìn)行曲線合并得到完整的分離曲線。仿真實(shí)驗(yàn)表明該理念的正確性。

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

4）在仿真中發(fā)現(xiàn)，由于骨架提取的固有缺陷，通過(guò)本文算法分離出的曲線與原始數(shù)據(jù)始終有一定的誤差，如何減少這一誤差，提高曲線分離準(zhǔn)確度將是下一步研究的重點(diǎn)。