基于隨機化張量算法的紅外弱小目標檢測

蹇 淵,黃自力,王 詢

(西南技術物理研究所,成都 610041,中國)

0 引 言

紅外弱小目標檢測技術一直是制導預警系統、機載、陸載、星載監視系統的關鍵技術。俄烏沖突中,小型無人機、高超音速彈藥的大量使用,顛覆了傳統戰爭模式,弱小目標檢測技術的需求變得越來越迫切。許多國家在該技術的研發上投入了大量科研力量。弱小目標檢測技術的難點在于:目標在圖像中通常僅占幾十甚至幾個像素,缺乏明顯的形狀、紋理、顏色等信息,對比度低,受周圍環境的干擾,目標容易被淹沒在背景中?；谀繕藦娂y理特征的算法適用性較差。工程中使用廣泛弱小目標檢測方法是基于單幀圖像的局部信息算法,具有代表性的算法包括:最大均值和最大中值濾波器[1]、頂帽(top-hat)變換[2]、局部對比度方法[3-4]等。這些算法通常認為目標灰度與周圍鄰域存在著較明顯的差異,并通過背景抑制、形態學估計等方法提取目標。此類算法的優點主要體現在:算法原理簡單、計算量通常較小、易于硬件平臺實現、實時性高;缺點是:算法基于的假設前提是目標在局部區域必須是最突出的,而在實際應用場景中,目標可能會出現在亮度較高的背景周圍(如云層、道路、建筑邊界等),使用局部信息檢測出目標將變得十分困難[5]。

隨著壓縮感知技術的興起,基于低秩稀疏分解模型[6]的紅外弱小目標檢測方法憑借其數學意義明確、解釋性強,受到了研究者們的青睞。該方法認為圖像的背景具有一定的非局部自相關的性質,相似性較高,表現出低秩的性質,而前景(弱小目標)在整幅圖像中僅占較少的像素,表現出稀疏的性質。因此,將原始圖像代入某種優化模型進行計算,便能實現圖像的前景與背景分離。2013年,GAO等人提出了紅外圖像塊(infrared patch-image,IPI)方法[7],該方法使用滑動窗口將單幀圖像拆分成大小相同的若干小塊,把獲得的圖像塊拉直成向量,按順序拼接成矩陣,使用優化模型提取矩陣的稀疏與低秩部分,再按原順序還原為圖像,從而獲得目標。由于該算法思路新穎、檢測性能較好,基于圖像塊的方法受到了研究者們的認可[8-9]?；趶埩康娜跣∧繕藱z測方法是基于圖像塊方法的擴展。提出這一方法的研究者通常認為:圖像塊經向量化后拼接為矩陣的方式破壞了圖像中的空間相關聯信息,這種信息的破壞導致了弱小目標檢測性能的下降。2017年,DAI等人首先將張量模型引入到弱小目標檢測算法中,即重加權紅外塊張量(reweighted infrared patch-tensor,RIPT)算法[10],該方法將圖像分塊并堆疊,構造為3階張量,并使用張量低秩稀疏分解模型進行處理,最終獲得分離后的弱小目標;模型中將3階張量直接按各指標展開成矩陣的操作導致了算法需要計算較大矩陣的奇異值分解(singular value decomposition,SVD)。為了提升算法性能,基于張量積運算的張量奇異值分解[11](tensor-singular value decomposition,t-SVD)逐漸在張量優化模型的求解中流行了起來,如:基于張量核范數[12]的張量魯棒性主成分分析方法[13]、基于張量核范數的部分和的方法[14]、基于張量加權核范數的方法[15]等,都應用到了t-SVD。

近幾年,隨著圖像探測器幀率的不斷提升,在基于張量的弱小目標檢測方法中,使用多幀圖像構造張量的方法成為了一個新的研究趨勢[16-18]。多幀圖像構造的張量又被稱為時空張量。提出這一方法的研究者通常認為,單幀圖像構造張量的方法僅使用了圖像的空間信息,而多幀圖像構造張量的方法能夠將時間信息與空間信息結合,帶來更好的目標檢測效果。

無論是基于單幀或多幀的張量弱小目標檢測方法,問題的求解通常被歸結于求解一個低秩稀疏分解優化模型。對于大多數模型,其求解步驟中通常需要計算t-SVD(即一系列矩陣的SVD),這一步驟通常是模型求解過程中計算復雜度最高的部分。為了降低該步驟的計算復雜度,一些使用隨機化方法的t-SVD替代算法被提出,例如:基于隨機化奇異值分解[19]的隨機化t-SVD[20]、基于矩陣骨骼分解的張量骨骼分解[21]等,這些算法的步驟是:首先使用隨機采樣、投影等方法對原始張量進行預處理,獲得體積較小的張量;再使用處理后的張量進行相應的分解,從而實現降低算法復雜度的目的。隨機化算法由于其令人滿意的可靠性與計算的有效性吸引了大量的算法研究者[22]。

基于圖像時空張量以及隨機化算法的思想,本文作者結合張量核范數與圖像結構張量[23]提出了一種基于時空張量的弱小目標檢測方法,并在算法求解過程中引入隨機化張量算法,有效地提升了算法的計算效率。

1 圖像時空張量

本文中所指的圖像時空張量是由圖像序列分塊、堆疊構造而來。分塊是使用固定尺寸的窗口將圖像序列中的每一幅圖像分成若干大小相同的小塊,堆疊即將分割后的圖像塊按照相應的位置進行堆疊,從而獲得一系列的圖像塊序列。將上述構造時空張量的方法展示為圖1。圖像堆疊的順序為圖中箭頭指向,將圖1右側任一個獨立的圖像塊序列稱為圖像時空張量。

圖1 圖像時空張量構造示意圖Fig.1 Illustration of spatial-temporal tensor construction

圖像塊序列記為張量A∈Cn1×n2×n3(C為復數域,n1×n2為圖像塊的大小,n3表示圖像塊序列中圖像塊的個數)。為了便于表示,將A(:,:,k)表示A的第k個正向切面(符號“:”為該維度的所有元素),即對應于圖像塊序列中的第k個圖像塊,顯然A(:,:,k)為一個n1×n2矩陣;A(i,j,k)表示A的第(i,j,k)元素,即圖像序列中第k個圖像塊中第i行、第j列處的圖像灰度值。

2 算法原理

一幅紅外圖像通?？梢员硎緸楸尘?、前景、噪聲三部分之和。在簡化的情況下,圖像被認為由前景和背景兩部分之和組成?；趶埩康募t外弱小目標檢測方法的原理是:通過某種操作將圖像構造為張量,并代入預先設計的張量優化模型,通過算法計算分離出目標。以下分別從張量低秩稀疏分解優化模型、模型的求解介紹本文作者提出的紅外弱小目標檢測算法。

2.1 基于張量低秩稀疏分解的紅外弱小目標檢測方法

本文中所指的基于張量的紅外弱小目標檢測方法是通過張量低秩稀疏分解模型的求解,將紅外弱小目標(稀疏部分)從原始的張量中分離出來。圖2展示了原始張量與低秩、稀疏張量間的關系。即原始張量A經模型計算后得到低秩張量B、稀疏張量T,滿足B+T=A。稀疏張量T即為包含弱小目標的前景張量。

圖2 張量稀疏、低秩分解示意圖Fig.2 Illustration of low-rank and sparse tensor decomposition

給定一個張量A∈Rn1×n2×n3(R表示實數域),一般形式的張量低秩稀疏分解模型如下:

(B+T=A)

(1)

即計算滿足條件B+T=A的張量B和T使得目標函數trank(B)+λ‖T‖0達到最小值。模型中trank(B)為B的秩,‖T‖0表示T非零元素的個數,λ>0為一設定的常數。麻煩的是:式(1)是一個非確定性多項式(nondeterministic polynomially,NP)問題,無法直接求解。通常的做法是將式(1)進行轉化,例如使用張量核范數[12]、張量核范數的部分和[24]等代替trank(B),使用T的l1范數[25]‖T‖1代替‖T‖0,以及在目標函數中添加其它的正則項抑制圖像背景中邊緣或噪聲等?？紤]到模型的簡潔,本文作者在張量穩健主成分分析(tensor robust principal component analysis,TRPCA)模型[12]的基礎上引入權重張量獲得如下的優化模型:

(B+T=A)

(2)

式中:K為A的權重張量[24];張量K⊙T表示K與T的哈達馬達積[25];‖K⊙T‖1為張量K⊙T的l1范數;‖B‖*為B的核范數(tensor nuclear norm,TNN)[12]。與參考文獻[14]中的優化模型不同之處在于式(2)使用了TNN對低秩張量B進行約束,這使得式(2)為一個凸優化問題。

2.2 模型的求解

式(2)的拉格朗日函數為:

L(B,T,Y,μ)=‖B‖*+λ‖K⊙T‖1+

(3)

式中:Y∈Rn1×n2×n3為拉格朗日乘子;μ>0為懲罰因子;〈Y,B+T-A〉為張量Y與張量B+T-A的內積;‖B+T-A‖F為張量B+T-A的Frobenius范數[25]。式(3)可使用交錯方向法(alternating direction multiplier method,ADMM)[26]通過迭代求解。根據參考文獻[12]和[14]中模型的求解方法,設第l步已獲得B,T,Y,μ,K的迭代值,使用下標l,分別記為Bl,Tl,Yl,μl,Kl,則第l+1步迭代為:

(4)

(5)

Yl+1=Yl+μl(Tl+1+Bl+1-A)

(6)

μl+1=min(ρμl,μmax)

(7)

(8)

式中:ρ>1為增長系數;μmax為設定的懲罰因子的最大值;c/(|Tl|+ε)為3階張量;|·|為取模運算;ε和c為大于0的常數;除法與加法運算均為元素運算;argmin(·)函數表示其括號中算式取得最小值時變量的取值。式(8)中的重加權方式見相關文獻[27],其作用是提高問題的收斂速度。

求解上述迭代問題的關鍵在于計算式(4)、式(5)的解。幸運的是式(4)、式(5)均有閉式解[12]。首先給出式(5)的閉式解表達式。

給定權重張量K∈Rn1×n2×n3張量以及常數τ>0,收縮算子SτK(·)(下標τK表示τ與張量K的各個元素相乘)的定義如下,對于任意的3階張量E∈Rn1×n2×n3,張量H=SτK(E)的元素滿足:

H(i,j,k)=sign(E(i,j,k))×

max(0,|E(i,j,k)|-τK(i,j,k))

(9)

式中:i=1,2,…,n1;j=1,2,…,n2;k=1,2,…,n3;sign(·)為符號函數;max(·,·)表示取兩元素的最大值運算。因此式(5)的閉式解為:

(10)

式中:Sμl-1λKl(·)是關于μl-1λKl的收縮算子,下標表示μl-1λ和張量Kl的各個元素相乘。

在TNN的意義下,式(4)的閉式解表達式可由張量奇異值閾值算子(tensor-singular value thresholding,t-SVT)[12]給出,t-SVT的實質是計算輸入張量的快速傅里葉(逆)變換與一系列矩陣的SVD。由于計算SVD復雜度通常較高(對于一個m×n矩陣,其SVD的時間復雜度為O(mn2)),這導致式(4)的求解是每一次迭代中計算復雜度最高的部分。為了降低算法的計算復雜度,本文中使用隨機化奇異值閾值算子(randomized tensor singular value thresholding,rt-SVT)Dτ(·)近似替代t-SVT算子,這里τ>0為常數。Dτ(·)的計算步驟在下節中給出,值得一提的是,參考文獻[28]和[29]中提到了類似的算法,但在構造正交投影張量時與本文中的算法略有差異。

根據上述算法,式(4)的近似解可表示為:

(11)

式中:D1/μl(·)為隨機化奇異值閾值算子。

2.2.2 使用ADMM方法求解 給定A∈Rn1×n2×n3以及大于0的常數λ,μ0,μmax,ρ,c和ε,迭代終止閾值η>0,最大迭代步數Nmax。使用ADMM求解式(2)的步驟如下:(a)將張量T,B,Y的各元素初始化為0,根據參考文獻[24]中的方法初始化A的權重張量K,即K0;(b)開始迭代,令迭代次數l=0,進行以下計算;(c)根據2.2.1中算法計算Bl+1=D1/μl(A-Tl-Yl/μl);(d)計算Tl+1=Sμl-1λKl(A-Yl/μl-Bl+1);(e)計算Yl+1=Yl+μl(Tl+1+Bl+1-A);(f)計算μl+1=min(ρμl,μmax),Kl+1=K0⊙(c/(|Tl|+ε));(g)取l=l+1,重復步驟(c)～(f),直到滿足如下條件之一終止迭代‖Bl-Bl+1‖∞≤η,‖Tl-Tl+1‖∞≤η,‖A-Bl+1-Tl+1‖∞≤η,l>Nmax(這里‖Bl-Bl+1‖∞為張量Bl-Bl+1的無窮范數[12]);(h)記B=Bl,T=Tl,并輸出B,T。當該算法中的輸入張量為由紅外圖像構造的時空張量時,則計算輸出的張量T即為包含弱小目標的前景張量。將其按原時空張量的構造順序還原,便得到了前景(目標)圖像序列,這一目標檢測過程如圖3所示。值得注意的是,紅外圖像序列構成的時空張量由多個大小為n1×n2×n3的張量構成,因此計算時需要將這些張量依次輸入本節中算法進行計算。

圖3 本文中的算法步驟Fig.3 Procedure of the proposed method in this paper

3 實驗與分析

3.1 實驗數據與評價指標

為了客觀地對算法的性能進行評價,這里使用信雜比增益(signal-to-clutter ratio gain,SCRG)、背景抑制因子(background suppression factor,BSF)對算法背景抑制的能力和目標增強的性能進行分析。在介紹SCRG前,首先引入目標局部信雜比(signal-to-clutter ratio,SCR)的定義。圖4展示了目標與背景區域的關系。圖中目標大小為a×b,局部背景大小為(a+2d)×(b+2d),本文中取常數d=20。設μtarget與μlocal分別表示目標區域與目標的局部鄰域灰度平均值;σlocal為目標局部區域灰度值的標準差,則目標的局部信雜比RSCR定義如下[10]:

圖4 目標與局部背景的關系示意圖Fig.4 Relationship between the target and local background

(12)

目標的局部信雜比反映了圖像中目標的檢測難易程度,通常局部信雜比越低,目標的檢測難度越大。根據目標局部信雜比的定義,設Rin、Rout分別為輸入圖像與輸出圖像(經算法處理后圖像)的目標局部信雜比;σin、σout分別表示輸入圖像與輸出圖像目標局部背景區域灰度的標準差,則目標的信雜比增益RSCRG與背景抑制因子RBSF分別定義如下[10]:

(13)

這兩個值反映了算法對目標局部背景的抑制程度,值越大,則說明算法對背景的抑制能力越強。

本實驗中使用了4組不同場景的紅外圖像序列,這些序列均來源于公開數據集[31]。圖5a～圖5d分別展示了這4組不同場景的紅外圖像。為了便于觀察,圖中使用箭頭指向的方框標出了弱小目標的位置。根據RSCR的定義,4組圖像中目標的平均RSCR、圖像的尺寸以及目標的平均尺寸分別見表1。觀察表1與圖4能夠看出,實驗中選取的圖像背景復雜,且目標具有尺寸小、信雜比低的特點。

表1 實驗數據的目標特性/pixelTable 1 Target properties of different experimental data /pixel

圖5 4種不同場景的原始數據紅外圖像Fig.5 Original infrared images with four different scenes

3.2 實驗結果

表2 不同算法的背景抑制性能比較Table 2 Comparison of background suppression performance of different algorithms

圖6 4種算法對應圖5中4種場景的處理結果Fig.6 Processing results of the four algorithms correspond to the four scenes in Fig.5

表3是上述不同算法分別處理4組數據中單幅圖像的平均用時(本實驗中使用的計算機配置為:Intel i5-8600K 6核3.6 GHz,主存24 GB)。根據表中數據,本文作者提出算法的計算速度明顯快于基于低秩稀疏分解方法的IPT與RIPT。雖然top-hat速度最快,但從表3中反映出的背景抑制能力來看,top-hat的背景抑制能力最弱,而本文中提出的算法較好地平衡了計算的效率與目標的檢測性能之間的關系,因此是有效的。

表3 不同算法處理單幅圖像的平均用時Table 3 Average processing time per frame of different algorithms

為了進一步探討本文中提出的隨機化張量算法與確定型算法的差異,將2.2.2節中的rt-SVT替換為t-SVT[12],即直接使用確定型算法求解式(2)。在保持與本文中隨機化算法參數一致的情況下,使用該確定型算法分別處理圖5a～圖5d中4組數據得到的處理結果如圖7所示。圖中方框標出的區域為目標真實區域。對比圖7與圖6發現,隨機化算法與確定型算法處理上述圖像序列,獲得的目標檢測結果幾乎一致。為了量化這種差異,將該確定型算法處理結果的SCRG值與BSF值進行比較,如表4所示。比較表4與表2中的數據發現,隨機化算法與確定型算法的背景抑制性能具有較小的差異,在圖5a、圖5d序列中,隨機化算法表現出了更好的背景抑制性能。

表4 確定型算法的背景抑制性能Table 4 Background suppression performance of the deterministic algorithm

圖7 確定型算法的處理結果展示Fig.7 Results of the deterministic algorithm

使用與上述實驗相同的計算機與軟件,獲得了確定型算法處理圖5a～圖5d 4組數據中單幀圖像的平均用時,如表5所示。通過對比表3與表5中的數據能夠發現,本文中提出的隨機化算法計算速度更快。上述對比實驗進一步表明,使用隨機化方法對確定型算法加速,不僅能夠較好地繼承原確定型算法的目標檢測性能,而且能夠有效地提升算法的計算速度。

表5 確定型算法處理單幅圖像的平均用時Table 5 Average processing time per frame of the deterministic algorithm

4 結 論

紅外弱小目標檢測[22]是具有挑戰性的問題,在工程應用方面,一些基于形態學濾波的算法使用廣泛,這些算法可以通過流水線處理進行加速,通常表現出速度快、實時性高等特點。與上述方法不同,基于張量低秩稀疏分解的方法是使用優化模型求解圖像問題,相比于傳統的算法,優化模型的方法通常表現出更好的圖像處理效果。然而基于優化的算法求解的過程需要進行多次迭代,每一步迭代通常都需要計算一系列張量(矩陣)分解,如SVD等,通常難以通過流水線操作加速,需要使用多(單)核處理平臺如:數字信號處理芯片等來計算,這些平臺提供了與線性代數相關的函數庫,為算法的實現提供了有效的工具。雖然這些庫已經獲得了較好的優化,但從原理上來講,計算SVD仍然是一個需要迭代的過程[30],其計算的復雜度與矩陣的尺寸密切相關。因此,減小參與SVD計算的張量(矩陣)尺寸是降低算法復雜度和減少硬件資源開銷的關鍵。

本文中將隨機化張量算法引入到張量優化問題的求解中,通過隨機投影方法降低了計算步驟中參與分解的張量尺寸,有效地提升了基于張量低秩稀疏分解的弱小目標檢測算法的計算效率。在圖像張量的構造方面,本文中使用了時空張量的構造方法,使得構造的張量較好地利用了圖像序列的空間與時間信息。雖然通過數值實驗驗證,本文作者提出的算法表現出了較好的性能,但仍然存在一些有待解決的問題,例如:不同場景的圖像數據如何自適應地設置隨機化算法中采樣參數、構造圖像時空張量時一些參數的設置等。這些問題的解決將有助于算法綜合性能的進一步提升。