Infrared small target detection based on randomized tensor algorithm

红外弱小目标检测技术一直是制导预警系统、机载、陆载、星载监视系统的关键技术。俄乌冲突中，小型无人机、高超音速弹药的大量使用，颠覆了传统战争模式，弱小目标检测技术的需求变得越来越迫切。许多国家在该技术的研发上投入了大量科研力量。弱小目标检测技术的难点在于: 目标在图像中通常仅占几十甚至几个像素，缺乏明显的形状、纹理、颜色等信息，对比度低，受周围环境的干扰，目标容易被淹没在背景中。基于目标强纹理特征的算法适用性较差。工程中使用广泛弱小目标检测方法是基于单帧图像的局部信息算法，具有代表性的算法包括：最大均值和最大中值滤波器^[1]、顶帽(top-hat)变换^[2]、局部对比度方法^[3-4]等。这些算法通常认为目标灰度与周围邻域存在着较明显的差异，并通过背景抑制、形态学估计等方法提取目标。此类算法的优点主要体现在：算法原理简单、计算量通常较小、易于硬件平台实现、实时性高；缺点是：算法基于的假设前提是目标在局部区域必须是最突出的，而在实际应用场景中，目标可能会出现在亮度较高的背景周围(如云层、道路、建筑边界等)，使用局部信息检测出目标将变得十分困难^[5]。

随着压缩感知技术的兴起，基于低秩稀疏分解模型^[6]的红外弱小目标检测方法凭借其数学意义明确、解释性强，受到了研究者们的青睐。该方法认为图像的背景具有一定的非局部自相关的性质，相似性较高，表现出低秩的性质，而前景(弱小目标)在整幅图像中仅占较少的像素，表现出稀疏的性质。因此，将原始图像代入某种优化模型进行计算，便能实现图像的前景与背景分离。2013年，GAO等人提出了红外图像块(infrared patch-image，IPI)方法^[7], 该方法使用滑动窗口将单帧图像拆分成大小相同的若干小块，把获得的图像块拉直成向量，按顺序拼接成矩阵，使用优化模型提取矩阵的稀疏与低秩部分，再按原顺序还原为图像，从而获得目标。由于该算法思路新颖、检测性能较好，基于图像块的方法受到了研究者们的认可^[8-9]。基于张量的弱小目标检测方法是基于图像块方法的扩展。提出这一方法的研究者通常认为: 图像块经向量化后拼接为矩阵的方式破坏了图像中的空间相关联信息，这种信息的破坏导致了弱小目标检测性能的下降。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右侧任一个独立的图像块序列称为图像时空张量。

Figure 1.  Illustration of spatial-temporal tensor construction

图像块序列记为张量$\boldsymbol{A} \in \mathbb{C}^{n_1 \times n_2 \times n_3}$($\mathbb{C}$为复数域，n₁×n₂为图像块的大小，n₃表示图像块序列中图像块的个数)。为了便于表示，将 A (: , : , k)表示 A的第k个正向切面(符号“：”为该维度的所有元素)，即对应于图像块序列中的第k个图像块，显然 A (: , : , k)为一个n₁×n₂矩阵; A (i, j, k)表示 A的第(i, j, k)元素，即图像序列中第k个图像块中第i行、第j列处的图像灰度值。

一幅红外图像通常可以表示为背景、前景、噪声三部分之和。在简化的情况下，图像被认为由前景和背景两部分之和组成。基于张量的红外弱小目标检测方法的原理是：通过某种操作将图像构造为张量，并代入预先设计的张量优化模型，通过算法计算分离出目标。以下分别从张量低秩稀疏分解优化模型、模型的求解介绍本文作者提出的红外弱小目标检测算法。

本文中所指的基于张量的红外弱小目标检测方法是通过张量低秩稀疏分解模型的求解, 将红外弱小目标(稀疏部分)从原始的张量中分离出来。图 2展示了原始张量与低秩、稀疏张量间的关系。即原始张量A经模型计算后得到低秩张量B、稀疏张量T，满足B+T=A。稀疏张量T即为包含弱小目标的前景张量。

Figure 2.  Illustration of low-rank and sparse tensor decomposition

给定一个张量$\boldsymbol{A} \in \mathbb{R}^{n_1 \times n_2 \times n_3}$($\mathbb{R}$表示实数域)，一般形式的张量低秩稀疏分解模型如下：

即计算满足条件B+T=A的张量B和T使得目标函数trank(B)+λ‖T‖₀达到最小值。模型中trank(B)为B的秩，‖T‖₀表示T非零元素的个数，λ>0为一设定的常数。麻烦的是：式(1)是一个非确定性多项式(nondeterministic polynomially，NP)问题，无法直接求解。通常的做法是将式(1)进行转化，例如使用张量核范数^[12]、张量核范数的部分和^[24]等代替trank(B)，使用T的l₁范数^[25]‖T‖₁代替‖T‖₀，以及在目标函数中添加其它的正则项抑制图像背景中边缘或噪声等。考虑到模型的简洁，本文作者在张量稳健主成分分析(tensor robust principal component analysis, TRPCA)模型^[12]的基础上引入权重张量获得如下的优化模型：

式中: K为A的权重张量^[24]；张量K⊙T表示K与T的哈达马达积^[25]; ‖K⊙T‖₁为张量K⊙T的l₁范数; ‖B‖_*为B的核范数(tensor nuclear norm, TNN)^[12]。与参考文献[14]中的优化模型不同之处在于式(2)使用了TNN对低秩张量B进行约束, 这使得式(2)为一个凸优化问题。

式(2)的拉格朗日函数为:

式中: $\boldsymbol{Y} \in \mathbb{R}^{n_1 \times n_2 \times n_3}$为拉格朗日乘子; μ>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，分别记为B_l, T_l, Y_l, μ_l, K_l，则第l+1步迭代为：

式中：ρ>1为增长系数; μ_max为设定的惩罚因子的最大值; c/(|T_l|+ε)为3阶张量; |·|为取模运算；ε和c为大于0的常数; 除法与加法运算均为元素运算; argmin(·)函数表示其括号中算式取得最小值时变量的取值。式(8)中的重加权方式见相关文献^[27]，其作用是提高问题的收敛速度。

求解上述迭代问题的关键在于计算式(4)、式(5)的解。幸运的是式(4)、式(5)均有闭式解^[12]。首先给出式(5)的闭式解表达式。

给定权重张量$\boldsymbol{K} \in \mathbb{R}^{n_1 \times n_2 \times n_3}$张量以及常数τ>0，收缩算子S_τK(·)(下标τK表示τ与张量K的各个元素相乘)的定义如下，对于任意的3阶张量$\boldsymbol{E} \in \mathbb{R}^{n_1 \times n_2 \times n_3}$，张量H=S_τK(E)的元素满足:

式中: i=1, 2, …, n₁; j=1, 2, …, n₂; k=1, 2, …, n₃; sign(·)为符号函数; max(·, ·)表示取两元素的最大值运算。因此式(5)的闭式解为：

式中: S_μl^-1λKl(·)是关于μ_l^-1λK_l的收缩算子，下标表示μ_l^－1λ和张量K_l的各个元素相乘。

在TNN的意义下，式(4)的闭式解表达式可由张量奇异值阈值算子(tensor-singular value thresholding，t-SVT)^[12]给出，t-SVT的实质是计算输入张量的快速傅里叶(逆)变换与一系列矩阵的SVD。由于计算SVD复杂度通常较高(对于一个m×n矩阵，其SVD的时间复杂度为O(mn²)), 这导致式(4)的求解是每一次迭代中计算复杂度最高的部分。为了降低算法的计算复杂度，本文中使用随机化奇异值阈值算子(randomized tensor singular value thresholding，rt-SVT)D_τ(·)近似替代t-SVT算子，这里τ>0为常数。D_τ(·)的计算步骤在下节中给出，值得一提的是, 参考文献[28]和[29]中提到了类似的算法, 但在构造正交投影张量时与本文中的算法略有差异。

给定常数τ>0，张量$\boldsymbol{A} \in\mathbb{C}^{n_1 \times n_2 \times n_3}$, s=min(n₁, n₂), 目标秩r≥2，过采样参数p≥2，且p+r < s，指数q=0或者1，2。本文中计算rt-SVT算子D_τ(·)的步骤如下：(a)沿A的第3维做快速傅里叶变换(fast Fourier transform, FFT)得到A=FFT(A, [], 3), 并生成标准高斯张量$\boldsymbol{G} \in \mathbb{R}^{n_2 \times(p+r) \times n_3}$, 即张量G中的每一个元素服从相互独立的标准正态分布; (b)沿G的第3维做快速傅里叶变换得到G =FFT(G, [], 3)；(c)开始循环，令循环次数i=1分别进行如下计算; (d)使用瘦正交三角分解^[30]计算矩阵(A (: , : , i)×A (: , : , i)^*)^qA (: , : , i)G (: , : , i)的列正交基，并记为矩阵Q (: , : , i)，上标*为复共轭转置运算；(e)计算矩阵Q (: , : , i)^*A (: , : , i)的SVD，记作U (: , : , i), S (: , : , i), V (: , : , i)，分别为其左奇异矩阵、奇异值矩阵与右奇异矩阵；(f)计算并更新U (: , : , i)=Q (: , : , i)U (: , : , i)；(g)使用张量收缩算子计算得到X(: , : , i)=U (: , : , i)S_τ(S (: , : , i))V (: , : , i)^*(这里X为中间变量，S_τ为张量K中元素全为1的收缩算子)；(h)令i=i+1，重复步骤(d)~(g)，直到i=「(n₃+1)/2+1时，停止循环(符号「·表示向上取整运算)；(i)开始循环，令循环次数j=「(n₃+1)/2 +1；(j)计算复共轭X(: , : , j)=conj(X(: , : , n₃－j+2))(conj(·)为取复共轭函数，即对每个元素取复共轭)；(k)令j=j+1，重复步骤(j)，直到j=n₃+1时，停止循环；(l)沿X的第3维做快速傅里叶逆变换(inverse fast Fourier transform, IFFT)得到D_τ(A)=IFFT(X, [], 3)；(m)输出D_τ(A)。

上述算法使用随机化方法，与确定型t-SVT算子^[12]的不同之处在于，该算法不直接计算矩阵A (: , : , i)的SVD。而是先使用了随机投影对矩阵进行预处理，获得了尺寸更小的矩阵Q (: , : , i)^*A (: , : , i)，再对该矩阵进行SVD，从而获得了A (: , : , i)的SVD近似估计(这种使用随机化投影方法近似计算SVD的方法便是随机化SVD^[19])。虽然在该算法中涉及了矩阵的瘦正交三角分解，但相比于直接计算原始矩阵的SVD花费的计算时间，算法的整体计算复杂度得到了降低^{[20, 28]}。

根据上述算法，式(4)的近似解可表示为:

式中: D_1/μl(·)为随机化奇异值阈值算子。

给定$\boldsymbol{A} \in \mathbb{R}^{n_1 \times n_2 \times n_3}$以及大于0的常数λ, μ₀, μ_max, ρ, c和ε，迭代终止阈值η>0，最大迭代步数N_max。使用ADMM求解式(2)的步骤如下：(a)将张量T, B, Y的各元素初始化为0，根据参考文献[24]中的方法初始化A的权重张量K，即K₀；(b)开始迭代，令迭代次数l=0，进行以下计算; (c)根据2.2.1中算法计算B_l+1=D_1/μl(A－T_l－Y_l/μ_l)；(d)计算T_l+1=S_μl^-1λKl(A-Y_l/μ_l-B_l+1)；(e)计算Y_l+1= Y_l+μ_l(T_l+1+B_l+1－A)；(f)计算μ_l+1=min(ρμ_l, μ_max)，K_l+1= K₀⊙(c/(|T_l|+ε))；(g)取l=l+1，重复步骤(c)~(f)，直到满足如下条件之一终止迭代‖B_l－B_l+1‖_∞≤η, ‖T_l－T_l+1‖_∞≤η, ‖A－B_l+1－T_l+1‖_∞≤η, l>N_max(这里‖B_l－B_l+1‖_∞为张量B_l－B_l+1的无穷范数^[12])；(h)记B=B_l, T=T_l，并输出B，T。当该算法中的输入张量为由红外图像构造的时空张量时，则计算输出的张量T即为包含弱小目标的前景张量。将其按原时空张量的构造顺序还原，便得到了前景(目标)图像序列，这一目标检测过程如图 3所示。值得注意的是，红外图像序列构成的时空张量由多个大小为n₁×n₂×n₃的张量构成，因此计算时需要将这些张量依次输入本节中算法进行计算。

Figure 3.  Procedure of the proposed method in this paper

为了客观地对算法的性能进行评价，这里使用信杂比增益(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为目标局部区域灰度值的标准差，则目标的局部信杂比R_SCR定义如下

Figure 4.  Relationship between the target and local background

目标的局部信杂比反映了图像中目标的检测难易程度，通常局部信杂比越低，目标的检测难度越大。根据目标局部信杂比的定义，设R_in、R_out分别为输入图像与输出图像(经算法处理后图像)的目标局部信杂比; σ_in、σ_out分别表示输入图像与输出图像目标局部背景区域灰度的标准差，则目标的信杂比增益R_SCRG与背景抑制因子R_BSF分别定义如下^[10]：

这两个值反映了算法对目标局部背景的抑制程度，值越大，则说明算法对背景的抑制能力越强。

本实验中使用了4组不同场景的红外图像序列，这些序列均来源于公开数据集^[31]。图 5a~图 5d分别展示了这4组不同场景的红外图像。为了便于观察，图中使用箭头指向的方框标出了弱小目标的位置。根据R_SCR的定义，4组图像中目标的平均R_SCR、图像的尺寸以及目标的平均尺寸分别见表 1。观察表 1与图 4能够看出，实验中选取的图像背景复杂，且目标具有尺寸小、信杂比低的特点。

Figure 5.  Original infrared images with four different scenes

分别使用IPI算法^[7]、RIPT算法^[10]、top-hat算法^[7]、以及本文作者提出的算法处理选取的4组数据。对比算法中各参数的取值与原文一致，可参见参考文献[10]。本文中算法的参数分别选取如下：2.2.2节中参数的取值分别为n₁=n₂=64, n₃=10, $\lambda=3 /\sqrt{n_3 \cdot \max \left(n_1, n_2\right)}$，N_max=10, η=10^－3; 2.2.1节中参数的取值分别为p+r=「n₃/3, q=1。部分图像的处理结果如图 6所示。图中第1~第4行分别为IPI、top-hat、RIPT和本文中算法的处理结果，从左至右第1~第4列分别对应图 5a~图 5d的4种不同场景。图中方框标出的区域为目标的真实区域，圆圈标出的区域为虚景。根据式(13)，计算得到各组图像的目标平均局部信杂比增益与背景抑制因子，如表 2所示。表中带下划线的数字表示各行中R_SCRG与R_BSF的最大值。通过观察图 6中虚景的抑制情况，并对比表 2中的数据能够看出，本文作者提出的算法对虚景有较强的抑制能力，因此表现出较好的目标检测性能。

Figure 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的背景抑制能力最弱, 而本文中提出的算法较好地平衡了计算的效率与目标的检测性能之间的关系，因此是有效的。

为了进一步探讨本文中提出的随机化张量算法与确定型算法的差异，将2.2.2节中的rt-SVT替换为t-SVT^[12]，即直接使用确定型算法求解式(2)。在保持与本文中随机化算法参数一致的情况下，使用该确定型算法分别处理图 5a~图 5d中4组数据得到的处理结果如图 7所示。图中方框标出的区域为目标真实区域。对比图 7与图 6发现，随机化算法与确定型算法处理上述图像序列，获得的目标检测结果几乎一致。为了量化这种差异，将该确定型算法处理结果的SCRG值与BSF值进行比较，如表 4所示。比较表 4与表 2中的数据发现，随机化算法与确定型算法的背景抑制性能具有较小的差异，在图 5a、图 5d序列中，随机化算法表现出了更好的背景抑制性能。

Figure 7.  Results of the deterministic algorithm

使用与上述实验相同的计算机与软件，获得了确定型算法处理图 5a~图 5d 4组数据中单帧图像的平均用时，如表 5所示。通过对比表 3与表 5中的数据能够发现，本文中提出的随机化算法计算速度更快。上述对比实验进一步表明，使用随机化方法对确定型算法加速，不仅能够较好地继承原确定型算法的目标检测性能，而且能够有效地提升算法的计算速度。

红外弱小目标检测^[22]是具有挑战性的问题，在工程应用方面，一些基于形态学滤波的算法使用广泛，这些算法可以通过流水线处理进行加速, 通常表现出速度快、实时性高等特点。与上述方法不同，基于张量低秩稀疏分解的方法是使用优化模型求解图像问题，相比于传统的算法，优化模型的方法通常表现出更好的图像处理效果。然而基于优化的算法求解的过程需要进行多次迭代，每一步迭代通常都需要计算一系列张量(矩阵)分解，如SVD等, 通常难以通过流水线操作加速, 需要使用多(单)核处理平台如：数字信号处理芯片等来计算，这些平台提供了与线性代数相关的函数库，为算法的实现提供了有效的工具。虽然这些库已经获得了较好的优化，但从原理上来讲，计算SVD仍然是一个需要迭代的过程^[30]，其计算的复杂度与矩阵的尺寸密切相关。因此，减小参与SVD计算的张量(矩阵)尺寸是降低算法复杂度和减少硬件资源开销的关键。

本文中将随机化张量算法引入到张量优化问题的求解中，通过随机投影方法降低了计算步骤中参与分解的张量尺寸，有效地提升了基于张量低秩稀疏分解的弱小目标检测算法的计算效率。在图像张量的构造方面，本文中使用了时空张量的构造方法，使得构造的张量较好地利用了图像序列的空间与时间信息。虽然通过数值实验验证，本文作者提出的算法表现出了较好的性能，但仍然存在一些有待解决的问题，例如：不同场景的图像数据如何自适应地设置随机化算法中采样参数、构造图像时空张量时一些参数的设置等。这些问题的解决将有助于算法综合性能的进一步提升。