張 健，婁樹理，任建存

(1．海軍航空工程學(xué)院控制工程系，山東 煙臺264001;2．中國人民解放軍91055部隊，浙江 臺州318050)

1 引言

本文研究的星圖是利用觀測相機(jī)拍攝的以深空為背景的圖像，圖像尺寸為512 pixels×512 pixels，圖像位數(shù)16 bits。星圖中主要包含恒星、空間目標(biāo)(主要指衛(wèi)星)和背景噪聲。恒星和空間目標(biāo)距離相機(jī)較遠(yuǎn)，其所成像僅占據(jù)有限個像素，為點(diǎn)狀分布;星空背景噪聲表現(xiàn)為變化緩慢的低灰度區(qū)域，理論上認(rèn)為星空背景噪聲符合高斯［1］或者泊松分布［2］。星圖中的運(yùn)動目標(biāo)檢測主要的難點(diǎn)包括:宇宙射線、高能粒子會在圖像中產(chǎn)生大量虛假目標(biāo);雜散光將使得圖像背景強(qiáng)弱分布不均勻，難以直接分割出目標(biāo);數(shù)量繁多的駐留物體導(dǎo)致視場內(nèi)出現(xiàn)多個目標(biāo);最后，星圖中恒星數(shù)量巨大，目標(biāo)在運(yùn)動過程中連續(xù)穿越恒星導(dǎo)致恒星對目標(biāo)產(chǎn)生遮擋。

為了對星圖中的運(yùn)動目標(biāo)進(jìn)行檢測，必須利用多幀圖像的運(yùn)動信息。對星圖中的運(yùn)動目標(biāo)進(jìn)行檢測，方法主要有跟蹤前檢測(Detection before Tracking，DBT)算法和檢測前跟蹤(Tracking before Detection，TBD)算法兩類。

基于DBT的空間目標(biāo)檢測算法必須要解決目標(biāo)分割和航跡關(guān)聯(lián)這兩個難題。文獻(xiàn)［3］～［6］中給出了比較有代表性的基于DBT的檢測算法。

基于TBD的空間目標(biāo)檢測算法主要有動態(tài)規(guī)劃算法［7］、極大似然檢測算法［1］、粒子濾波算法［8］以及時序多幀投影算法［2，9，10］等。

理論上，文獻(xiàn)［11］中提出的全維匹配濾波的檢測性能最高，但其計算量太大而得不到實(shí)際應(yīng)用，為解決這一問題，文獻(xiàn)［2］中引入投影算法的思想，提出MTI算法，以損失一定信噪比為代價大幅度降低了算法的復(fù)雜度。美國SBV計劃中的“空間中段試驗(yàn)”衛(wèi)星上就使用了MTI算法，使其成為成功應(yīng)用于天基空間目標(biāo)監(jiān)視系統(tǒng)的在軌檢測算法。

MTI算法是一類典型的TBD檢測算法，該算法首先利用時序多幀最大值投影進(jìn)行背景抑制，然后利用速度濾波檢測目標(biāo)運(yùn)動軌跡。本文將MTI算法應(yīng)用于目標(biāo)檢測的圖像預(yù)處理階段，由于MTI算法要求序列圖像精確配準(zhǔn)。因此，本文首先研究星圖的精確配準(zhǔn)方法。

2 基于SURF算法的星圖精確配準(zhǔn)

近來，基于特征點(diǎn)的圖像配準(zhǔn)引起了眾多研究者的重視。Lowe［12］在前人研究的基礎(chǔ)上，提出了SIFT(Scale Invariant Feature Transform)算法。Bay等人在SIFT算法的基礎(chǔ)上提出了SURF［13］算法。SURF算法將DoH(Detemination of Hessian)中的高斯二階微分模板進(jìn)行了近似簡化，使得模板對圖像的濾波只需要進(jìn)行幾個簡單的加減法運(yùn)算，并且這種運(yùn)算與濾波模板的尺寸無關(guān)。由于SURF算法在運(yùn)算速度上優(yōu)于SIFT算法，因此本文選擇SURF算法進(jìn)行圖像配準(zhǔn)。

由于恒星距離相機(jī)遙遠(yuǎn)，因此相鄰兩幀星圖可以看作服從剛體變換模型。為得到剛體變換模型參數(shù)，至少需要兩幀圖像中3對以上的匹配特征點(diǎn)，本文以20幀圖像為1組進(jìn)行配準(zhǔn)，以每組第1幀圖像作為基準(zhǔn)圖像，后續(xù)圖像都與第1幀配準(zhǔn)。將SURF算法檢測閾值設(shè)為10－12，分解層數(shù)設(shè)為1，第1和第2幀中檢測到的匹配特征點(diǎn)有18對，如圖1所示。

圖1 SURF算法特征點(diǎn)提取匹配結(jié)果Fig．1 Feature points matching results with SURF algorithm

根據(jù)式(1)利用最小二乘法求解待配準(zhǔn)圖像的旋轉(zhuǎn)角度β和平移參數(shù)(dx，dy)，并利用參數(shù)對待配準(zhǔn)圖像進(jìn)行雙三次灰度插值。

式中，F(xiàn)為匹配特征點(diǎn)對數(shù)。

本文利用RMSE衡量星點(diǎn)質(zhì)心的配準(zhǔn)精度，星點(diǎn)質(zhì)心計算采用一階矩法。質(zhì)心的RMSE定義如下［14］:

其中，N為觀測圖像提取的星點(diǎn)個數(shù);(xi，yi)為參考圖像中星點(diǎn)質(zhì)心位置;(x'i，y'i)為配準(zhǔn)后圖像中星點(diǎn)質(zhì)心位置，為兩質(zhì)心間的歐氏距離。各幀中檢測到的匹配恒星數(shù)量如圖2所示。

圖2 各幀中檢測到的恒星數(shù)量Fig．2 The number of s detected tars in each frame

各幀中恒星質(zhì)心坐標(biāo)的RMSE如圖3所示。采用本文算法配準(zhǔn)后，星圖質(zhì)心的RMSE最小達(dá)到0．3269 pixel，平均值為0．5441 pixel，完全滿足后續(xù)時序多幀投影對配準(zhǔn)精度的要求，同時配準(zhǔn)精度要高于文獻(xiàn)［14］中的結(jié)果。

圖3 星象質(zhì)心的RMSE曲線Fig．3 RMSE curve of star centroid

3 基于時序多幀投影的星圖預(yù)處理

本文研究的運(yùn)動目標(biāo)幀間的運(yùn)動速度通常大于10個像素，導(dǎo)致在投影圖像中運(yùn)動軌跡不連續(xù)。本文將MTI算法作適當(dāng)改進(jìn)，用來對星圖背景進(jìn)行抑制。假設(shè)序列圖像中像素灰度表示為r(x，y，t)，其中，(x，y)表示空間坐標(biāo)，x=1，…，m，y=1，…，n。t表示時間序列，t=1，…，N，則本文改進(jìn)的MTI算法的數(shù)學(xué)計算過程如下:

(1)得到序列圖像最大值投影:

(2)得到中值投影:

(3)得到標(biāo)準(zhǔn)差投影:

其中，max2［r(x，y)］為序列圖像第二大值。

(4)從最大值中減去中值:

(5)利用標(biāo)準(zhǔn)差歸一化:

對序列星圖最大值投影歸一化后的結(jié)果如圖4所示。

圖4 本文算法序列星圖最大值投影結(jié)果Fig．4 Maximum value projection of the star images with improved MTI algorithm

(6)成對二值化:

相機(jī)光學(xué)系統(tǒng)一般都進(jìn)行了散焦處理，使得目標(biāo)成像的尺寸最小為2 pixels×2 pixels［15］。根據(jù)這一結(jié)論，本文對MTI算法中的成對二值化方法進(jìn)行改進(jìn):

其中，閾值T由恒虛警原理確定。圖5給出了本文算法對最大值投影的二值化結(jié)果。

圖5 最大值投影的二值化結(jié)果Fig．5 Binarization of maximum value projection map

(7)解幀:

最大值投影時會產(chǎn)生一個標(biāo)記幀F(xiàn)rame，F(xiàn)rame(x，y)記錄坐標(biāo)(x，y)處最大值對應(yīng)的幀序號。利用Frame將b(x，y)中屬于同一幀的所有“1”像素提取出來，組成一幅新的二值圖像T(x，y，t)，t=1，…，N。

本文對MTI算法的改進(jìn)在于:一是將步驟(2)中的均值投影改為中值投影;二是將步驟(6)的成對二值化方法改進(jìn)以減少虛假目標(biāo)數(shù)量;三是對最大值投影進(jìn)行解幀，方便后續(xù)的速度濾波。

4 目標(biāo)運(yùn)動軌跡檢測

4．1 目標(biāo)初始運(yùn)動狀態(tài)的建立

從二值圖像T(x，y，k)，t=1，…，N中可以提取所有疑似目標(biāo)的信息，包括星點(diǎn)的質(zhì)心坐標(biāo)，星象區(qū)域長寬和像素數(shù)，星象的信噪比等。本文后續(xù)所有字符的上標(biāo)代表圖像幀序號，將第k幀中所有疑似目標(biāo)的集合標(biāo)記為Uk。假定Uk中包含少量目標(biāo)點(diǎn)和大量噪聲點(diǎn)，在圖像連續(xù)幀中，可以認(rèn)為空間目標(biāo)作勻速直線運(yùn)動，而噪聲出現(xiàn)的位置是隨機(jī)的，因此利用目標(biāo)運(yùn)動的連續(xù)性和一致性檢測目標(biāo)。假設(shè)Uk+1、Uk+2、Uk+3中的3個星點(diǎn)的坐標(biāo)依次為Pk+1i=，這3個星點(diǎn)的運(yùn)動可描述為:

建立目標(biāo)初始運(yùn)動狀態(tài)的步驟如下:

(2)在U3中找一個星點(diǎn)P，根據(jù)式(9)計算V23和A23。

(3)若|V12－V23|和|A12－A23|均較小，表明上述3個星點(diǎn)在幀間運(yùn)動距離基本相等，運(yùn)動方向基本一致，可以認(rèn)為是目標(biāo)點(diǎn)，把V12和A12作為的運(yùn)動狀態(tài)記錄下來。將除去，加入目標(biāo)集合After1。

(4)重復(fù)步驟(1)～步驟(3)，完成第一幀的后向關(guān)聯(lián)，得到目標(biāo)集合After1。

(5)利用步驟(1)～步驟(4)的相同方法進(jìn)行前向關(guān)聯(lián)，得到目標(biāo)集合Before1，取After1和Before1中目標(biāo)個數(shù)較多的集合作為真實(shí)目標(biāo)集合T1。

4．2 軌跡檢測及目標(biāo)插值

目標(biāo)初始運(yùn)動狀態(tài)建立之后，圖像中仍然存在一定虛假目標(biāo)。本文利用一種二元累積(Binary Integration，BI)的方法進(jìn)行速度濾波。本文BI檢測算法中首先建立目標(biāo)鏈，再將滿足條件的目標(biāo)加入目標(biāo)鏈，目標(biāo)鏈的長度為M，目標(biāo)鏈的起始幀為k，結(jié)束幀為s，s－k+1=N，最后根據(jù)N與M的關(guān)系識別該目標(biāo)鏈。限于篇幅，速度濾波具體流程見圖6中斜線標(biāo)注的區(qū)域。

由于目標(biāo)自身灰度值相比高亮恒星要低，當(dāng)其穿越恒星時，恒星會對目標(biāo)產(chǎn)生遮擋效應(yīng)，在MTI背景抑制過程中會將目標(biāo)屏蔽掉，造成目標(biāo)丟幀。以往檢測算法能應(yīng)對短暫丟幀現(xiàn)象，但當(dāng)恒星數(shù)量巨大，弱小目標(biāo)連續(xù)穿越密集恒星群，丟幀嚴(yán)重時，算法性能下降。BI檢測確認(rèn)目標(biāo)鏈之后，如果在第s幀中該目標(biāo)預(yù)測位置附近沒有其他目標(biāo)出現(xiàn)，說明該目標(biāo)在第s幀丟幀，此時需要根據(jù)第(s－1)幀中目標(biāo)的運(yùn)動信息，對第s幀目標(biāo)進(jìn)行插值，插值流程見圖6中細(xì)點(diǎn)標(biāo)注的區(qū)域。

圖6 速度濾波以及目標(biāo)插值流程圖Fig．6 Flow chart of velocity filtering and interpolation

4．3 檢測結(jié)果

最終經(jīng)本文算法檢測到的20幀序列圖像中目標(biāo)的運(yùn)動軌跡由圖7給出。

圖7 目標(biāo)的運(yùn)動軌跡Fig．7 Trajectory of targets

5 實(shí)驗(yàn)結(jié)果

5．1 算法的檢測概率

時序多幀投影檢測算法的檢測概率PD與投影幀數(shù)N，信噪比S，式(8)中的二值化閾值T，以及目標(biāo)鏈長度M等有關(guān)。檢測算法采用恒虛警原理，當(dāng)虛警概率PFA恒定時，選取T和M使PD最大化。對每一個可能的M值(M=0，…，N)，首先根據(jù)PFA得到T，再結(jié)合S，得到使PD最大化的(T，M)組合。

令PFA=10－2，可得到前述的(T，M)組合，結(jié)果由表1給出。

表1 N=20時，BI檢測算法的檢測閾值TTab．1 The detection threshold T of BI algorithm with N=20

根據(jù)表1的結(jié)果，對不同的(T，M)組合計算BI算法的檢測概率，由計算結(jié)果可知，當(dāng)N=20時，M=8可以取得最高的檢測概率。相應(yīng)地，根據(jù)表1，閾值T=1．9。圖8給出了N=20，T=1．9時，不同信噪比目標(biāo)檢測概率與M的關(guān)系曲線。

圖8 T=1．9時，不同信噪比目標(biāo)的檢測概率曲線Fig．8 T=1．9，different SNR，detection probability curves

利用蒙特卡洛方法分別進(jìn)行105次實(shí)驗(yàn)，得到本文算法與文獻(xiàn)［2］中SBV算法的檢測概率曲線，如圖9中所示。由實(shí)驗(yàn)結(jié)果可見在低信噪比條件下，本文算法的檢測概率要高于SBV算法。

圖9 本文算法與SBV算法檢測概率對比Fig．9 Comparison of detection probability of our method with SBV method

5．2 算法的檢測速度

檢測算法的軟件開發(fā)平臺為MATLAB 2008a，操作系統(tǒng)為Windows 7，硬件環(huán)境為Core i3四核CPU 2．5 GHz，內(nèi)存為4G的PC電腦。表2給出了分別利用文獻(xiàn)［8］中算法和本文算法檢測4組實(shí)拍觀測圖像的平均每幀耗時情況。由結(jié)果可見本文算法的檢測速度明顯快于文獻(xiàn)［8］中算法，主要原因在于本文利用SURF算法檢測圖像特征點(diǎn)，避免了遍歷恒星檢測控制點(diǎn)這個耗時的過程。SURF算法是速度得到極大提升的局部特征檢測方法，用其對觀測圖像進(jìn)行配準(zhǔn)精度高而且速度快。

表2 算法耗時情況Tab．2 Time consuming result of algorithm

6 結(jié)論

由于DBT檢測算法需要事先將圖像配準(zhǔn)，而利用恒星控制點(diǎn)進(jìn)行配準(zhǔn)需要遍歷所有恒星，算法耗時且精度低。為解決這一問題，論文提出了一種基于SURF算法精確配準(zhǔn)和時序多幀投影的TBD檢測算法。算法首先利用SURF描述子將圖像精確配準(zhǔn)并將配準(zhǔn)后的序列圖像在空間上對齊;然后利用改進(jìn)的MTI算法對圖像進(jìn)行最大值投影，并對投影結(jié)果解幀，提取單幀中的疑似目標(biāo);最后利用基于恒虛警的BI檢測器對目標(biāo)軌跡進(jìn)行檢測識別，并對目標(biāo)丟幀位置進(jìn)行插值處理。該算法解決了大尺寸觀測圖像高精度、快速配準(zhǔn)問題，同時解決了SBV中MTI算法無法檢測非連續(xù)運(yùn)動軌跡的問題。

［1］ Li Z Z，Li R X，Wei H G．Research of Image denoising based on high precision of star sensors［C］．6th International Symposium on Advanced Optical Manufacturing and Testing Technologies:Design，Manufacturing，and Testing of Smart Structures，Xiamen，2012:1－6．

［2］ Chu P L．Efficient detection of small moving objects［R］．Lexington，Massachusetts Institute of Technology，Lincoln Laboratory，1992:1－70．

［3］ WABG Weihua，HE Yan，CHEN Zengping．Real—time algorithm for small moving Target detection in photoelectric image sequences［J］．Opto-Electronic Engineering，2006，33(4):14－18．(in Chinese)王衛(wèi)華，何艷，陳曾平．光電圖像序列運(yùn)動弱目標(biāo)實(shí)時檢測算法［J］．光電工程，2006，33(4):14－18．

［4］ZHANG Chunhua，ZHOU Xiaodong，CHEN Weizhen．Target trace acquisition method of star images based on background elimination［J］．Infrared and Laser Engineering，2008，(37):143－146．(in Chinese)張春華，周曉東，陳維真．基于背景抑制的星空圖像目標(biāo)運(yùn)動軌跡提?。跩］．紅外與激光工程，2008，(37):143－146．

［5］ SUN Jinqiu，ZHANGYanning，JIANGLei，et al．Detection of a dim and small moving target using a temporal-spatial fusion filtering algorithm［J］．Mechanical Science and Technology for Aerospace Engineering，2009，28(1):20－24．(in Chinese)孫瑾秋，張艷寧，姜磊，等．基于時空域融合濾波的弱小運(yùn)動目標(biāo)檢測算法［J］．機(jī)械科學(xué)與技術(shù)，2009，28(1):20－24．

［6］CHENG Jun，ZHANG Wei，CONG Mingyu，et al．Research of detecting algorithm for space object based on star map recognition［J］．Optical Technique，2010，36(3):439－444．(in Chinese)程軍，張偉，叢明煜，等．基于星圖識別的空間目標(biāo)檢測算法研究［J］．光學(xué)技術(shù)，2010，36(3):439－444．

［7］ ZHANG Yuye，WANG Chunxin．Space small targets detection based on improved DPA［J］．Acta Electronica Sinica，2010，38(3):556－560．(in Chinese)張玉葉，王春歆．基于改進(jìn)DPA的空間小目標(biāo)檢測算法［J］．電子學(xué)報，2010，38(3):556－560．

［8］ LIU Zhiyong，LILei，YAO Rui，et al．Small and dim space object detection method based on particle filter［J］．Chinese Journal of Stereology and Image Analysis，2011，16(2):113－117．(in Chinese)劉志勇，李磊，姚睿，等．基于粒子濾波的空間弱小目標(biāo)檢測方法［J］．中國體視學(xué)與圖像分析，2011，16(2):113－117．

［9］ CONG Mingyu，HE Wenjia，LU Lihong，et al．Trace extraction of moving point targets in complex background images［J］．Optics and Precision Engineering，2012，20(7):1619－1625．(in Chinese)叢明煜，何文家，逯力紅，等．復(fù)雜背景成像條件下運(yùn)動點(diǎn)目標(biāo)的軌跡提?。跩］．光學(xué) 精密工程，2012，20(7):1619－1625．

［10］CONG Mingyu，HE Wenjia，BAO Wenzhuo，et al．Suppression of cluttered cloud image background by time series multi-frame projection［J］．Optics and Precision Engineering，2012，20(4):826－834．(in Chinese)叢明煜，何文家，鮑文卓，等．云雜波成像背景的時序多幀投影抑制［J］．光學(xué) 精密工程，2012，20(4):826－834．

［11］Carbone CP，Kay S M．A novel normalization algorithm based on the three-dimensional minimum variance spectral estimation［J］．IEEE Transactions on Aerospace and Electronic System，2012，48(1):430－448．

［12］Lowe D G．Distinctive image features from scale-invariant keypoints［J］．International Journal of Computer Vision，2004，60(2):91－110．

［13］Bay H，Ess A，Tuytelaars T，et al．Speeded-up robust features(SURF)［J］．Computer Vision and Image Understanding，2008，110(3):346－359．

［14］JIAO Jichao，ZHAO Baojun，TANG Linbo．Space image registration algorithm based on nonsubsampled contourlet transform and MLESAC［J］．Systems Engineering and Electronics，2010，32(12):2686－2690．(in Chinese)焦繼超，趙保軍，唐林波．一種基于非降采樣Contourlet變換和MLESAC的星空圖像配準(zhǔn)算法［J］．系統(tǒng)工程與電子技術(shù)，2010，32(12):2686－2690．

［15］Sun T，Xing F，You Z．Optical system error analysis and calibration method of high-accuracy star trackers［J］．Sensors，2013，13(4):4598－4623．