基于Shi-Tomasi角點檢測算法的彈道目標平動補償

李玉璽 馮存前 許旭光 韓立珣 李 江

(1. 空軍工程大學研究生院, 陜西西安 710051; 2. 空軍工程大學防空反導學院, 陜西西安 710051)

1 引言

當今世界,彈道導彈技術實現了快速地發展[1],對各國產生了重大影響。為了應對彈道導彈的威脅,各國都在抓緊研究如何對彈道導彈進行攔截。對于彈道導彈防御這項技術來說,首要的是對彈道導彈進行準確識別[2],識別出真假彈頭。然而,隨著分導式多彈頭以及導彈誘餌技術的日益成熟[3],在飛行過程中,彈道導彈會釋放假彈頭以及誘餌,同時自身也會產生許多碎片,對雷達形成干擾。因此,對于依靠傳統特征進行識別的方法來說,很難從群目標中準確識別出真彈頭,從而使得很難依靠傳統的識別技術獲得彈頭的相關參數。然而,物體的微小運動是獨一無二的,因此可以利用該特性對彈道導彈進行識別[4]。Victor C. Chen首先提出了微多普勒的概念[5],并將其引入到目標識別領域中。之后,利用微多普勒信息識別目標成為新的研究熱點。

一般情況下,彈道導彈在飛行中,其運動由平動和微動組成。其中,平動會導致散射點回波多普勒曲線出現傾斜,影響彈道導彈微動信息的提取。因此,為了提取微動信息,應該進行平動補償。目前,對平動分量的補償已經進行了大量的研究并取得了一定的成果。文獻[6]采用高階模糊函數,搜索圖像峰值估計出平動參數,實現了目標的平動補償,但求解需要相對較長的觀測時間。文獻[7]通過Radon變換,對目標整體微多普勒曲線進行投影變換,進而實現平動補償。文獻[8]提出了一種延遲共軛相乘的方法,消除微動保留平動,對目標進行平動補償。文獻[9]利用Viterbi算法,根據多普勒率與微多普勒關系得到微動參數,進而實現平動補償。文獻[10]通過提取交叉點信息,得到平動參數,實現了對目標的平動補償。

基于上述的研究現狀,本文提出了Shi-Tomasi角點檢測算法[11],進行平動補償。同時,針對檢測時間較長的問題,本文對算法進行了優化,更好的估計出平動參數,實現了較高精度的平動補償。

2 錐體彈道目標的進動模型

在飛行的過程中,彈道目標的飛行姿態受到多種因素的影響。為了簡化分析,本文只考慮重力對飛行姿態的影響。若重力作用在導彈的質心上,則飛行姿態僅僅為自旋;若重力沒有作用在導彈的質心上,則彈道目標的飛行姿態為一方面繞著自身的對稱軸進行自旋,另外一方面繞著空間某一定軸做錐旋運動。自旋和錐旋的疊加,最終形成了進動的飛行姿態[12]。

圖1 進動錐體彈道目標的散射模型Fig.1 Scattering model of precession cone ballistic target

如圖1所示,坐標系O-XYZ為錐體目標的參考坐標系,其原點O為錐體彈道目標的質心,坐標系中的Z軸為錐體運動的進動軸。在目標運動的起始時刻,雷達視線LOS在參考坐標系O-XYZ中的方位角和仰角分別為α和β,與自旋軸的夾角為φ,錐旋角與自旋角的夾角,即進動角為θ,目標質心與雷達的初始距離為R0(這里沒有在圖中標明)。目標的錐旋角速度為ωc,自旋角速度為ωs,平動速度為ν。質心到散射點A的距離為l1,到底面圓圓心距離為l2。

根據文獻[13],可以得到,散射點的微動距離變化公式如下所示,即

rA=l1cosβ(t)
rB=-l2cosβ(t)-rsinβ(t)
rC=-l2cosβ(t)+rsinβ(t)

(1)

其中:

cosβ(t)=cosθcosα+sinθsinαsin(ωct-φ0)

φ0為初始時刻t=0時,目標對稱軸在面XOY內的投影與OX軸的夾角。由于夾角φ0僅僅影響微多普勒的初相,對周期以及幅值沒有影響,因此我們可以取任意一數值,不失一般性,可以令φ0=0[14]。

當雷達波束照射到目標的時候,目標會形成三個典型散射點,錐頂散射點A點,底面邊緣滑動散射點B點、C點,如圖1所示。在實際應用中,由于存在遮擋效應,三個散射點不一定同時被雷達觀測到。根據文獻[15],我們可以得到相關散射點的可見角度,如下表1所示,其中ε為錐體目標的半錐角。

表1 散射點的可見角度

本文對角度在0<β<(ε∩(π-ε))的范圍內進行分析,此時三個散射點均可產生回波,對每個散射點產生的回波進行疊加,形成一個總的回波[16]。下面對散射點的頻率與運動之間的關系進行分析。其他情況均可照此分析,這里不在贅述。

圖2 雷達與錐體目標散射點之間距離示意圖Fig.2 Schematic diagram of the distance between the radar and the cone target scattering point

本文以散射點A為例。如圖2所示,雷達位于O′點位置。在實際的飛行過程中,彈道目標距離雷達較遠,因此O′O與O′A的夾角比較小,此時我們可以近似得到|O′A|≈|O′H|、|O′A|≈|O′O|+|OA|,即|O′H|=RA=R+rA。

多普勒頻率可以由下式得到

(2)

式中:fd為多普勒頻移,νr為散射點沿雷達視線方向的徑向速度,λ為載波波長。

當目標向著雷達運動時,νr>0,回波載頻比發射的載頻大,反之νr<0,回波載頻比發射的載頻小。因此,雷達只要能夠測出散射點徑向距離的變化率,即可得到散射點的頻率[17]。

(3)

式中:νrA散射點A沿雷達視線方向的徑向速度,ft和fm分別為目標平動與微動導致的多普勒頻率。

在目標的運動過程中,彈道目標的平動分量數學表達式可以近似用多項式來描述[9],即

(4)

若時間比較短,采用低階的多項式數學模型就可以較好的描述平動分量。本文采用3階多項式數學模型來對平動分量進行描述,即:

R=R0+a1t+a2t2/2+a3t3/6

(5)

式中:R0表示錐體目標距雷達的初始距離,a1表示錐體目標沿雷達視線的徑向速度,a2表示錐體目標沿雷達視線的徑向加速度,a3表示錐體目標沿雷達視線的2階徑向加速度。

因此,錐體目標運動產生的頻率為

(6)

假設雷達發射載頻為f0的單頻信號,即

st(t)=exp(j 2πf0t)

(7)

則雷達接收到的各個散射點回波為

(8)

式中:Ai為散射點回波的系數,Ri(t)為t時刻散射點與雷達間的距離。

接收的回波信號包含載頻信號exp(j 2πf0t),在采樣之前,可以通過正交解調予以去除[18]。解調后的散射點的基帶信號可以表示為

(9)

根據上式,可以得到散射點的時頻曲線圖。通過觀察目標散射點的時頻曲線圖,可以看出在時頻曲線圖上存在多個交點。

在實際中,彈道目標的平動速度比較大,產生的多普勒頻率要遠遠大于雷達脈沖重復頻率,得到的回波頻率遠大于基帶[-PRF/2,PRF/2],從而會發生模糊或折疊[10]。若僅僅考慮多普勒模糊情況,可以將目標多普勒頻率fi拆分如下式所示

fi=n·PRF+fn

(10)

式中:n·PRF為模糊項,n=round(fi/PRF)為模糊數,fn為剩余項。根據式(6)和(10)得到

(11)

(12)

此時,錐頂散射點A點,底面邊緣滑動散射點B點、C點具有相同的頻率值,也即時頻曲線圖的交點。在這些交點處,由式(12)可知,頻率值僅僅與平動項有關。因此,可以通過對交點信息的提取得到平動分量,進而實現對散射點的平動補償[10]。

3 時頻方法

由于時頻對仿真實驗存在相關影響,因此對時頻方法進行分析。

3.1 短時傅里葉變換

短時傅里葉變換(short-time Fourier transform, STFT)的基本原理是用時寬很窄的窗函數提取信號,這時可以將取出的信號看作是平穩的,用傅里葉變換進行分析,得到頻率信息[19]。

信號s(t)的短時傅里葉變換定義為

(13)

式中g(t)為一個時間寬度很短的窗函數。

3.2 Gabor變換

Gabor展開可以表示為

(14)

式中:amn稱為Gabor展開系數,g(t)為Gabor基函數,ΔT和Δf分別表示時間和頻率采樣間隔。

amn可以用下面式子求出

(15)

式中:γ(t)是與g(t)雙正交對偶的基函數[20]。

3.3 Cohen類時頻分布

Cohen類時頻分布可以統一寫成

(16)

Cohen類時頻分布還可以用模糊函數表示,即

(17)

式中:F是∏(t,f)的二維傅里葉變換[21]。

Rihaczek分布定義為:

Rx(t,f)=x(t)x*(f)e-j2πft

(18)

它對應于在Cohen類時頻分布中取F(τ,ν)=ejπ ντ[22]。

Page分布定義為:

(19)

它表示信號在時間t之前的能量譜密度的微分,也是Cohen類時頻分布的一種,對應取F(τ,ν)=e-jπν|τ|[22]。

4 提取角點的圖像特征

4.1 角點檢測算法

根據第2節所述,為了實現平動補償,需要確定時頻曲線中交點的位置,為此引入角點檢測算法。角點一般可以定義為兩條邊的交點,或者可以理解為在角點的周圍內具有兩個不同區域的不同方向的邊界。角點檢測算法(Corner Detection Algorithm)是用來獲得圖像中特征點的一種方法,也稱為特征點檢測算法[23]。算法的基本原理是選取某一區域,在圖像上任意移動選取的區域。通過對比移動前后,區域中像素點灰度值變化來確定是否存在角點。任意移動所選取的區域,如果圖像的灰度值產生較大變化,那么可以認為在該區域中存在角點[24]。

Moravec角點檢測算法是最早應用的算法。首先,該算法以某一像素點為中心確定一個區域,每個像素點都會形成一個區域。算法通過計算相鄰區域的平方差之和,檢測相鄰區域的相關性,進而確定是否存在角點。然而,該算法有明顯的缺陷。算法只是單純計算了幾個確定方向的像素值,未實現各向同性。同時由于所取的范圍簡單的正方形區域,很容易產生噪聲。為了解決上述相關問題,提出了Harris算法。首先,Harris算法計算矩陣行列式值與跡的差值,求出角點響應函數,之后進行閾值處理,得到局部極大值。但是該算法不具有尺度不變性,同時檢測所需的時間較長,實際應用較差。

Shi-Tomasi算法是Harris算法的改進。相對于Harris算法來說,Shi-Tomasi算法通過計算自相關矩陣的特征值,用所得特征值與設定的閾值進行比較。若特征值均大于閾值,則可以確定角點。Shi-Tomasi算法具有良好的可重復性,比較簡單,檢測效率較高,同時解決了Moravec算法無法對圖像進行各向性計算等問題。在很多情況下,可以得到更好的結果。

如圖3所示,窗口在平面上移動。若在光滑平面上,如圖3(a),窗口在任意方向上沒有發生任何變化;若在邊緣區域上,如圖3(b),窗口在邊緣所對應的方向沒有發生變化;若在角點區域上,如圖3(c),窗口在各個方向均發生變化。根據圖3,Shi-Tomasi算法利用窗口移動判斷在任意方向上是否發生變化,對圖像中的角點進行檢測。

(2)該橋上部結構車行道實心板梁的部分鉸縫失效，個別1.5 m車行道板、1.0 m車行道板存在單梁受力的跡象，因此車行道實心板梁的檢算均按單梁受力狀態進行，同時給出鉸縫完好狀態下的檢算結果。

圖3 窗口運動示意圖Fig.3 Schematic diagram of window movement

在圖像(x,y)處,對窗口進行平移,平移量為[u,ν],則移動前后的圖像灰度值變化為

(20)

式中:w(x,y)為加權函數,位于(x,y)處的窗口,I(x,y)為(x,y)處像素點灰度值大小,I(x+u,y+ν)為(x+u,y+ν)處像素點灰度值大小。

根據全微分公式可以得到:

I(x+u,y+ν)=I(x,y)+Ixu+Iyν+o(u2,ν2)

(21)

帶入式(20)中得到

(22)

將上式改為矩陣形式,如下式所示

(23)

記:

(24)

M為2×2維矩陣,則式(23)可以化簡為

(25)

設矩陣M的特征值為λ1、λ2,定義

R=min(λ1,λ2)

(26)

R>threshold

(27)

即可提取R的局部最大值。根據局部最大值R進行非極大值抑制,實現對角點的檢測。

根據上述所述的原理,Shi-Tomasi算法對圖像中所有的像素點進行檢測,得到大量角點,但是多數角點并不是我們需要的,同時導致檢測時間較長。為此對該算法進行改進,剔除不需要的角點,縮短檢測時間。

具體主要是在計算響應函數之前,對圖像中的像素點進行初步篩選,確定候選的角點。算法將檢測區域內中心像素點與邊緣像素點的灰度值進行比較,確定兩者的相似度。若兩者的灰度值之差在一個給定的范圍之內,那么這兩個像素點是相似的,反之,是不相似的。計算每個邊緣像素點與中心像素點灰度值之差,可以得出檢測區域內與中心像素點相似的個數n,根據n的值確定是否為候選角點。

通過該方法,可以將大量無關的像素點剔除,大大縮短后續程序所需要的時間,提高了算法的效率。具體算法步驟如下:

步驟3確定候選角點。對窗口內邊緣像素點與中心像素點灰度值進行比較,判斷是否在所給出的閾值之內,得到相似像素點的個數,進而確定是否為待選角點。

步驟4利用公式(26)計算像素點的響應函數R。

步驟5進行非極大值限制,判斷像素點是否存在角點。

4.2 回歸分析

在前文的敘述中,平動分量利用3階多項式模型描述即可,即

R=R0+a1t+a2t2/2+a3t3/6

(28)

由公式(2)可知多普勒頻率fd與目標的徑向速度νr正相關,即

fd∝νr

(29)

而νr=-dR/dt,因此,在時頻曲線圖曲線中的交點應該成拋物線分布。根據角點的參數,對上述的交點進行回歸分析,得出平動分量的相關參數。

假設在一直角坐標系中,有n個點,其坐標分別為(x1,y1),(x2,y2),…,(xn,yn)?，F在對這n個點進行k階多項式擬合,即

f(x)=a0+a1x+a2x2+...+akxk

(30)

將上述點的橫坐標代入k階多項式中得

(31)

為了讓擬合效果更好,我們需要使多項式計算出來的f(xi)與真實yi之間最小。這里可以利用下式進行評價:

(32)

即

(33)

此時,利用偏導求出最小值。下面將z分別對a0、a1、a2、…、ak求偏導,得到下面的式子。令式子均為0,求解可得到a0、a1、a2、…、ak的值。

(34)

化簡整理得到如下多個方程式

(35)

將其轉換成矩陣形式,令

(36)

(37)

則矩陣形式如下所示

XTXa=XTY

(38)

解得

a=(XTX)-1XTY

(39)

5 仿真實驗

5.1 時頻方法對檢測效果的影響

利用短時傅里葉變換對信號進行分析,得到圖4(a)。利用Gabor變換,對信號進行時頻分析得到散射點時頻曲線圖,如圖4(b)所示。利用Rihaczek分布,對信號進行時頻分析得到散射點時頻曲線圖,如圖4(c)所示。利用Page分布,對信號進行時頻分析得到散射點時頻曲線圖,如圖4(d)所示。

從圖4(b)可以看出時頻效果較差。這是由于在實際的臨界采樣情況下,窗函數不能構成標架,不能保證其時域和頻域的緊支撐性能,無法反映信號在時頻平面能量的分布[20]。從圖4(c)、(d)可以看到,交叉項比較嚴重,時頻曲線圖中交點無法清晰反映出來,后續檢測會產生較大誤差。而根據圖4(a)可以看出,對于短時傅里葉變換,交叉項的影響較小,可以較好提取出角點進行檢測。因此本文選擇短時傅里葉變換進行仿真實驗。

圖4(a)是利用短時傅里葉變換仿真的錐體目標在飛行的過程中三個散射點產生回波的時頻曲線圖,從圖中我們可以看到是時頻曲線圖不是水平的,而是帶有一定傾斜角度的,說明這里面存在著平動分量。同時還可以觀察到,三條曲線有多個相交點,根據前文所述,這些交點只包含平動分量信息,因此可以利用這些交點進行平動分量的補償。

圖4 散射點時頻曲線圖Fig.4 Time-frequency curve of scattering points

采用高斯平滑濾波,對圖4(a)進行平滑處理得到二值化圖像,如圖5所示。之后再對圖5進行骨架抽取,如圖6所示。

圖5 二值化圖像Fig.5 Binarized image

圖6 骨架提取圖像Fig.6 Skeleton extraction image

根據上文所述的角點檢測算法,對圖6的時頻曲線圖進行角點檢測,提取相關角點,如圖7所示。觀察圖7,算法檢測出大量的角點。利用2.1小節的改進方法對算法進行改進,剔除不需要的角點,得到新的角點分布圖,如圖8所示。通過觀察,曲線有9個交點,其中第4個、第6個、第8個交點檢測存在較大誤差,其余6個交點均可通過算法準確檢測出來,確定出相關交點的位置信息。由于交點只包含平動分量信息,根據所得的交點位置信息可以得出平動分量信息。根據圖像行列數與時頻圖坐標對應關系,將確定的位置轉化成時間-頻率信息。圖5～圖8是進行一次仿真得到的仿真圖。實驗中進行多次蒙特卡羅仿真。通過實驗條件仿真觀察時頻圖,每個離散時間對應一個角點,這為利用多項式回歸分析估計平動分量提供了可能。根據式(6)與式(39)可以估計得到徑向速度a1=-2.0009 m/s,徑向加速度a2=1.97915 m/s2,2階徑向加速度a3=0.9976 m/s3。對比假設的數據,計算相對誤差分別為εa1=0.045%、εa2=1.0425%、εa3=0.24%。

圖7 角點圖Fig.7 Corner diagram

圖8 圖像中顯示特征點Fig.8 Feature points displayed in the image

同時,在進行提取的過程中,比較改進前后的Shi-Tomasi算法計算響應函數R以及確定角點所耗費時間。在改進之前,其所耗費的時間為0.728964 s,在改進之后,其所耗費的時間為0.022917 s。通過對比,可以看到改進之后的算法所需時間大約為原算法的3.14%,大大縮短了相關的計算時間,更有利于在實際中進行應用。

上述實驗均是在同一臺計算機上進行的,計算機的硬件配置參數:處理器Intel i5-10210U、內存16G、顯卡型號MX250。

根據公式(6),減去估計出的平動分量,得到散射點的微動分量。圖9為減去平動分量之后的錐體目標散射點的時頻曲線圖。根據圖像所示,平動分量已經去掉,得到了較好的微動分量時頻曲線圖。驗證了該算法的有效性。

圖9 平動補償后的散射點時頻曲線圖Fig.9 Time-frequency curve of scattering points after translation compensation

5.2 不同參數下算法的檢測效果

下面在不同信噪比下,對本文所提出的算法的檢測效果進行進一步分析。接下來在不同的信噪比條件下進行多次蒙特卡羅仿真,其他仿真參數與上面設置一致。在不同的信噪比條件下仿真,每個離散時間對應一個角點,利用多項式擬合得到的仿真結果如表2所示。同時,對平動參數a1、a2、a3進行整體的誤差分析,可以得到參數的聯合相對誤差,如表2所示。

表2 不同信噪比下的相對誤差

根據表2所得到的數據可以看出,在不同的信噪比下,平動分量所產生的相對誤差,基本保持穩定,因此該算法具有良好的魯棒性。表中有極個別數據波動較大,可能是由于提取的角點偏差較大。

同時,相較于文獻[6]的方法,可以看出本文所提的算法,參數估計的精度較好,補償的效果更好。

針對不同的彈頭特性及系統參數對圖像灰度的影響,通過改變彈頭的自旋頻率以及雷達發射信號的脈沖重復頻率進行仿真實驗,其結果如表3、表4所示。

表3 不同自旋頻率下的相對誤差

表4 不同脈沖重復頻率下的相對誤差

根據仿真實驗結果可以得到,增大自旋頻率以及脈沖重復頻率可以提高實驗結果的檢測精度。

6 結論

彈道目標的雷達時頻曲線圖是由多個散射點疊加而成,但由于目標存在平動運動,因此為了實現對彈道目標的識別,分析相關特性,需要進行平動補償。為此,本文提出了Shi-Tomasi角點檢測算法,通過分析角點的特性對平動分量進行了補償。仿真實驗驗證了該算法的有效性、魯棒性。該算法實現了對圖像的各向性進行檢測,比較簡單,同時又對算法進行了優化,縮短了運行時間,提高了在實際應用中的效率。

根據表2可以看出,由于提取的角點不精確,進而導致在估計平動補償的相關參數時產生較大的誤差。因此,接下來可以繼續對算法進行改進,使提取的角點更加精確。表3針對不同的彈頭特性對圖像灰度的影響進行了仿真,隨著自旋頻率的增大,檢測的點數隨之增加,角點檢測精度提高,誤差變小。而表4針對系統參數對圖像灰度的影響進行了仿真,隨著脈沖重復頻率的增大,時頻效果更好,角點檢測精度提高,誤差變小。

在實際中,目標的散射點特性都比較復雜,同時在目標的飛行過程中,伴隨著目標飛行的行程中,還存在大量的彈頭碎片以及誘餌。因此,下一步,可以對在復雜散射點情況下對群目標的微動特性進行研究。

同時,對中段彈道目標進行平動補償之后,下一步工作就是對得到的信號進行分離。通過對散射點時頻曲線的研究,得到彈道目標的特性。