红外弱小目标识别技术在军事和民用领域中都扮演着十分重要的角色。在自然场景中识别红外小目标这一过程的特点是，目标区域通常可以表述为一个尺寸小、均匀、紧凑的区域，与周围背景相比，具有显著不连续性或者强烈对比性。在一幅红外图像中检测小目标，会受到目标像素数量少、目标和背景之间对比度较低、目标边缘暗淡、图像背景复杂以及小目标纹理信息缺失等因素的影响。

红外弱小目标检测近年来被学者们广泛研究。在早期方法中，Deshpande等人^[1]提出了最大均值滤波和最大中值滤波的方法用来检测红外图像中的弱小目标，Zeng等人^[2]提出了顶帽形态学滤波器以及在红外目标检测中的应用，这些方法都着力于设计一种滤波器来抑制背景，从而提取小目标。Gu等人^[3]提出了一种基于内核的非参数回归方法，用于红外图像中背景预测和小目标检测，Dong等人^[4]提出了一种用于红外点目标检测的均匀背景预测模型，这些算法同样致力于背景抑制工作，但是都对复杂背景噪声颇为敏感。之后，Qi等人^[5]设计了一种基于二阶方向导数和布尔图视觉理论的小目标检测算法，这种方法在检测率上有了很大提升，但是在应对复杂背景中小目标检测时会有较高虚警率。近几年，稀疏表示理论在信号处理与分析的领域中有着出色的表现，目标样本与背景样本的有效利用在目标检测上呈现出优良效果^[6]，Wang等人^[7]根据稀疏表示理论提出了一种利用背景字典分块表示原红外图像方法来检测小目标位置的算法，Wan等人^[8]根据稀疏表示理论提出了一种利用目标字典分块表示原红外图像从而检测出小目标位置的算法。

现有方法在红外图像中的小目标检测上已经有了较好结果，但是也存在着对复杂背景适应性低、检测率较低、虚警率较高等问题，而现有稀疏相关方法存在的问题在于：直接从原图上入手构建字典，没有考虑目标的特征提取，或者没有同时构建目标和背景字典，字典的表示能力较弱。针对这些缺陷，本文提出了一种利用红外图像的Facet方向导数特征与稀疏表示理论相结合来检测弱小目标的方法。

Haralick提出^[9]:一个图像中任意一个像素点某一个邻域所有像素点灰度值所形成的灰度强度表面，都可以被一个空间的二元三次多项式拟合，现在假定这个映射函数为$ f$，可以通过计算两组1维离散多项式的张量积，来构建2维离散正交多项式基组。又由于现在要得到一个三次函数，就可以忽略高于三阶的多项式。根据Haralick模型，针对一个像素所在5×5邻域进行具体地阐述，现定义一组离散1维数组$ \mathit{\boldsymbol{R}}$={-2 -1 0 1 2}，和另一个离散1维数组$ \mathit{\boldsymbol{L}}$={-2 -1 0 1 2}，这两组1维离散数组的张量积$ \mathit{\boldsymbol{R}} \times \mathit{\boldsymbol{C}}$为

$ \mathit{\boldsymbol{R}} \times \mathit{\boldsymbol{C}} = \left\{ \begin{array}{l} 1,r,c,{r^2} - 2,rc,{c^2} - 2,\\ {r^3} - 17/5r,\left( {{r^2} - 2} \right)c,\\ r\left( {{c^2} - 2} \right),{c^3} - 17/5c \end{array} \right\} $

(1)

令$ \mathit{\boldsymbol{S}}$为图像的某一像素点$ (x, y)$的5×5邻域，对于这个5×5邻域中有$ {\rm{(}}\mathit{r}{\rm{, }}\mathit{c}{\rm{)}} \in \mathit{\boldsymbol{S}}$, 其中$ (r, c)$为一个像素点坐标，$ r$和$ c$取值都是{-2, -1, 0, 1, 2}，像素点$ (x, y)$在这个5×5邻域中对应坐标为中心点(0, 0)，以此类推。那么，之前所述用于拟合一个图像中任意一个像素点5×5邻域二元三次多项式$ f(r, c)$就可以表示为

$ \begin{array}{*{20}{c}} {f\left( {r,c} \right) = {K_1} + {K_2}r + {K_3}c + {K_4}\left( {{r^2} - 2} \right) + }\\ {{K_5}rc + {K_6}\left( {{c^2} - 2} \right) + {K_7}\left( {{r^3} - 17/5r} \right) + }\\ {{K_8}\left( {{r^2} - 2} \right)c + {K_9}r\left( {{c^2} - 2} \right) + }\\ {{K_{10}}\left( {{c^3} - 17/5c} \right)} \end{array} $

(2)

式中，$ {\mathit{K}_\mathit{i}}\left( {i = 1, \cdots , 10} \right)$是该二元三次多项式的拟合系数，它是由最小二乘法表面拟合而成。这种方法在对图像灰度强度表面进行拟合的同时，能够将一些椒盐噪声去除，因此，这种方法不仅能计算导数，还具有一定程度的滤波效果。

在本文方法中需要求解图像中任意一个像素点0°、90°、45°和-45°方向的一阶导数，如图 1所示。

$ {{f'}_0}\left| {_{\left( {0,0} \right)}} \right. = {K_3} - 2{K_8} - 3.4{K_{10}} $

(3)

$ {{f'}_{90}}\left| {_{\left( {0,0} \right)}} \right. = {K_2} - 2{K_9} - 3.4{K_7} $

(4)

$ \begin{array}{*{20}{c}} {{{f'}_{45}}\left| {_{\left( {0,0} \right)}} \right. = {{f'}_0}\left| {_{\left( {0,0} \right)}} \right. \times \cos \left( {{{45}^ \circ }} \right) + }\\ {{{f'}_{90}}\left| {_{\left( {0,0} \right)}} \right. \times \sin \left( {{{45}^ \circ }} \right)} \end{array} $

(5)

$ \begin{array}{*{20}{c}} {{{f'}_{ - 45}}\left| {_{\left( {0,0} \right)}} \right. = {{f'}_0}\left| {_{\left( {0,0} \right)}} \right. \times \cos \left( { - {{45}^ \circ }} \right) + }\\ {{{f'}_{90}}\left| {_{\left( {0,0} \right)}} \right. \times \sin \left( { - {{45}^ \circ }} \right)} \end{array} $

(6)

利用式(3)—(6)求出原图像的4个方向上的方向导数特征。

在自然场景中，目标区域通常可以表述为一个小尺寸、均匀、紧凑的区域，没有明显形状尺寸和纹理信息，与周围背景相比，具有显著不连续性或者强烈对比性。在通常情况下，一幅灰度图中小目标可以看做一个2维高斯函数模型^[10]

$ I\left( {x,y} \right) = A\exp \left\{ { - \frac{1}{2}\left[ {{{\left( {\frac{x}{{{\sigma _1}}}} \right)}^2} + {{\left( {\frac{y}{{{\sigma _2}}}} \right)}^2}} \right]} \right\} $

(7)

式中，$ I(x, y)$表示$ (x, y)$坐标处灰度值，$ A$表示函数幅值，$ {\sigma _1}$和$ {\sigma _2}$表示对应标准差。理想情况下，目标在某一个方向上呈现图 2(a)中曲线分布，其一阶导数呈现如图 2(b)中曲线分布，而背景区域在理想情况下某一个方向上一阶导数曲线近似于曲线$ y$=0型分布，利用目标和背景在一阶导数上呈现出分布差异这一特点，就可以对它们加以区分。

假如现有一个从原图像中截取的尺寸为$ m \times m$的小块，对其求某一方向上一阶导数得到一个尺寸为$ m \times m$的导数块，将导数块按列合并为一个列向量$ \mathit{\boldsymbol{Y}}$。取$ M=m \times m$，假如有一个$ M \times N_{\rm t}$的目标字典$ {\mathit{\boldsymbol{D}}_{\rm{t}}}$和一个$ M \times N_{\rm b}$的背景字典$ {\mathit{\boldsymbol{D}}_{\rm{b}}}$(关于字典获得方式将在2.2.1节加以阐述)，对导数块向量$ \mathit{\boldsymbol{Y}}$进行稀疏表示求解

$ \begin{array}{*{20}{c}} {\min {{\left\| {{\mathit{\boldsymbol{X}}_t}} \right\|}_0} + {{\left\| {{\mathit{\boldsymbol{X}}_b}} \right\|}_0}}\\ {{\rm{s}}.\;{\rm{t}}.\;\;{\mathit{\boldsymbol{D}}_{\rm{t}}}{\mathit{\boldsymbol{X}}_{\rm{t}}} + {\mathit{\boldsymbol{D}}_{\rm{b}}}{\mathit{\boldsymbol{X}}_{\rm{b}}} = \mathit{\boldsymbol{Y}}} \end{array} $

(8)

式中，$ {\mathit{\boldsymbol{X}}_{\rm{t}}}$表示目标系数，$ {\mathit{\boldsymbol{X}}_{\rm{b}}}$表示背景系数，在理想情况下，原信号$ \mathit{\boldsymbol{Y}}$被表示成了稀疏信号$ {\mathit{\boldsymbol{X}}_{\rm{t}}}$和$ {\mathit{\boldsymbol{X}}_{\rm{b}}}$，如果$ \mathit{\boldsymbol{Y}}$是一个目标导数块向量，那么求解的稀疏系数中只有$ {\mathit{\boldsymbol{X}}_{\rm{t}}}$中有数值，$ {\mathit{\boldsymbol{X}}_{\rm{b}}}$中没有数值，如果$ \mathit{\boldsymbol{Y}}$是一个背景导数块向量，那么求解的稀疏系数中只有$ {\mathit{\boldsymbol{X}}_{\rm{b}}}$中有数值，$ {\mathit{\boldsymbol{X}}_{\rm{t}}}$中没有数值，于是，目标和背景在该求导方向上就得以区分。

然而，在实际红外图像中，背景是非常复杂的，实际红外图像中背景部分可能包含：“斜坡”、“角点”、“悬崖”、“山谷”等，并不会形如理想情况下在某一个方向上一阶导数近似于一个平面，这些复杂背景成分会对目标检测造成不同程度的干扰^[11]。尽管如此，由于目标可以被看做一个2维高斯函数模型，它在各个不同方向上都能近似为图 2(b)中所示曲线分布，这是目标在导数域上的独有特性^[12]，如图 3所示，而上述背景成分均不可能在各个不同方向上都呈现出近似于形如目标曲线的分布，利用这一特性，可以用Facet模型求解原图像各个方向上的一阶导数，分别求解稀疏系数，再进行多方向上融合处理，即可将目标和背景彻底区分开。

基于上述分析，本文提取多方向导数特征，通过多方向共同求解，应对实际红外图像中背景复杂的问题。本节首先介绍字典的获取与构建方式，然后详细阐述多方向求解的方法，最后说明图像重建和目标检测的方法。

字典获取分为背景字典获取和目标字典获取两个部分。

对于背景字典，在图像库中不包含目标的背景区域里，截取1 000个尺寸为16×16像素的小块^[13]，对这1 000个小块分别利用Facet模型求4个方向上的一阶导数，然后分别按列合并得到4个256维列向量，再对它们进行归一化，由此在每一个方向上都可以得到1 000个背景样本。

为了求解出最稀疏的稀疏系数，缩短求解时间，同时还要为输入信号的稀疏表示搭建出一个很好的过完备字典，本文选择了平稳有效的K-SVD算法^[14]对背景样本进行训练。K-SVD算法效率高，而且针对图像复原和图像压缩具有较好应用效果，K-SVD算法是一种通过迭代来实现字典训练的方法，在一个总循环中，主要包含稀疏编码和字典更新两个步骤。

利用K-SVD算法训练某一方向上1 000个背景样本，得到该方向上背景字典$ {\mathit{\boldsymbol{D}}_{\rm{b}}}$，其包含500个原子。

对于目标字典的搭建，首先考虑目标小块生成过程，红外图像中小目标的特点主要包含以下3个方面：1)小目标是呈高斯分布的圆形或者椭圆形的连续、均匀、紧凑的亮区域；2)针对椭圆形小目标，有时还具有一定旋转角度；3)目标所占总像素数一般是2×2到8×8^[15]。根据小目标所具有的特点一共生成325个含有小目标的小块，在这些目标块中，圆形小目标所具有的直径包括4~16的所有整数，椭圆形小目标的长轴包括5~16的所有整数，短轴包括从4到对应长轴长度减1的所有整数。此外，椭圆形小目标在0°基础上还分别旋转了45°、90°和135°。对这些含有小目标的小块求某一方向一阶导数，再按列并成256维列向量，然后对它们进行归一化，得到该方向上含有325个原子的目标字典$ {\mathit{\boldsymbol{D}}_{\rm{t}}}$，值得说明的是，每一个原子对应一种目标尺寸和旋转角度，这将用于后续图像重建。

引入多方向对小目标检测具有很好的效果^[16]，本文取0°、90°、45°和-45°方向这4个方向进行求解。现对式(8)进行改进，令大字典$ \mathit{\boldsymbol{D}}$=[$ {\mathit{\boldsymbol{D}}_{\rm{t}}}$ $ {\mathit{\boldsymbol{D}}_{\rm{b}}}$]，总稀疏系数$ \mathit{\boldsymbol{X}}{\rm{ = }}{\left[ {\mathit{\boldsymbol{X}}_{\rm{t}}^{\rm{T}}\;\mathit{\boldsymbol{X}}_{\rm{b}}^{\rm{T}}} \right]^{\rm{T}}}$，则式(8)可化为

$ \begin{array}{*{20}{c}} {\min {{\left\| \mathit{\boldsymbol{X}} \right\|}_0}}\\ {{\rm{s}}.\;{\rm{t}}.\;\;\mathit{\boldsymbol{DX}} = \mathit{\boldsymbol{Y}}} \end{array} $

(9)

在$ n$个方向上共同求解，式(9)可扩展为

$ \begin{array}{*{20}{c}} {\min \sum\limits_{i = 1}^n {{{\left\| {{\mathit{\boldsymbol{X}}_i}} \right\|}_0}} }\\ {{\rm{s}}.\;{\rm{t}}.}\\ {\left[ {{\mathit{\boldsymbol{Y}}_1}{\mathit{\boldsymbol{Y}}_2} \cdots {\mathit{\boldsymbol{Y}}_n}} \right] = \left[ {{\mathit{\boldsymbol{D}}_1}{\mathit{\boldsymbol{D}}_2} \cdots {\mathit{\boldsymbol{D}}_n}} \right] \cdot }\\ {\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{X}}_1}}&0& \cdots &0\\ 0&{{\mathit{\boldsymbol{X}}_2}}& \cdots &0\\ \vdots&\vdots &{}& \vdots \\ 0&0& \cdots &{{\mathit{\boldsymbol{X}}_n}} \end{array}} \right]} \end{array} $

(10)

令$ \underline{{\mathit{\boldsymbol{\bar{Y}}}}}\text{=}\left[ {{\mathit{\boldsymbol{Y}}}_{\text{1}}}{{\mathit{\boldsymbol{Y}}}_{2}}\cdots {{\mathit{\boldsymbol{Y}}}_{n}} \right], \underline{{\mathit{\boldsymbol{\bar{D}}}}}\text{=}\left[ {{\mathit{\boldsymbol{D}}}_{\text{1}}}{{\mathit{\boldsymbol{D}}}_{2}}\cdots {{\mathit{\boldsymbol{D}}}_{n}} \right]$

$ \underline {\mathit{\boldsymbol{\bar X}}} = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{X}}_1}}&0& \cdots &0\\ 0&{{\mathit{\boldsymbol{X}}_2}}& \cdots &0\\ \vdots&\vdots &{}& \vdots \\ 0&0& \cdots &{{\mathit{\boldsymbol{X}}_n}} \end{array}} \right] $

(11)

则式(10)可简化为

$ \min \sum\limits_{j = 1}^n {{{\left\| {\underline {\mathit{\boldsymbol{\bar X}}} {._j}} \right\|}_0}} \;\;\;\;{\rm{s}}.\;{\rm{t}}.\;\;\underline {\mathit{\boldsymbol{\bar Y}}} = \underline {\mathit{\boldsymbol{\bar D}}} \cdot \underline {\mathit{\boldsymbol{\bar X}}} $

(12)

式中，$ {{\underline{{\mathit{\boldsymbol{\bar{X}}}}}}_{\cdot \mathit{j}}}$表示$ \underline{{\mathit{\boldsymbol{\bar{X}}}}}$的第$ j$列。该式可以用OMP算法求解^[17]。

将原红外图像从上至下，从左至右，以步长为4截取16×16像素的图像小块，对它们求0°、90°、45°和-45° 4个方向上一阶导数，利用式(12)求解稀疏系数矩阵$ \underline{{\mathit{\boldsymbol{\bar{X}}}}}$，每一个方向上稀疏系数$ {{\mathit{\boldsymbol{X}}}_{\mathit{i}}}$由目标系数$ {{\mathit{\boldsymbol{X}}}_{\text{t}\mathit{i}}}$和背景系数$ {{\mathit{\boldsymbol{X}}}_{\text{b}\mathit{i}}}$组成。对每一个导数块在每一个方向上计算残差

$ {\gamma _i} = \left\| {{\mathit{\boldsymbol{Y}}_i} - {\mathit{\boldsymbol{D}}_{{\rm{b}}i}}{\mathit{\boldsymbol{X}}_{{\rm{b}}i}}} \right\|_2^2 $

(13)

利用目标系数$ {{\mathit{\boldsymbol{X}}}_{\text{t}\mathit{i}}}$重建目标

$ {\mathit{\boldsymbol{G}}_i} = \sum\limits_{k = 1}^K {{\mathit{\boldsymbol{X}}_{{\rm{t}}i}}\left( k \right) \cdot \mathit{\boldsymbol{B}}\left( k \right)} $

(14)

式中，$ K$为目标系数$ {{\mathit{\boldsymbol{X}}}_{\mathit{ti}}}$的维数，$ \mathit{\boldsymbol{B}}$为所有用于搭建目标字典$ {\mathit{\boldsymbol{D}}_{\rm{t}}}$时对应的16×16像素目标小块的集合，该集合每一项和目标系数每一项一一对应。接着利用4个方向上$ {{\gamma }_{i}}$和$ {{\mathit{\boldsymbol{G}}}_{\mathit{i}}}$得到所截取16×16像素图像小块的重建小块

$ \mathit{\boldsymbol{P}} = \frac{1}{n}\sum\limits_{i = 1}^n {{\gamma _i} \cdot {\mathit{\boldsymbol{G}}_i}} $

(15)

最后利用所有重建小块重建原图像，小块重叠部分采用均值滤波方式取值，由此得到最终检测结果图。

针对最终得到的重建图像，即检测结果图，取一个阈值$ T$将目标与背景分开

$ T = \bar I + q \cdot {\sigma _r} $

(16)

式中，$ {\bar{I}}$是检测结果图数值的均值，$ {{\sigma }_{\mathit{r}}}$是检测结果图数值的标准差，$ q$为定常数，一般取10左右。

本文算法的流程图如图 4所示。

ROC(receiver operating characteristic)曲线是衡量目标检测结果的一个常用指标，能够反映检测率$ DR$和虚警率$ FA$之间的动态关系^[17]。记总目标数为$ T_{\rm NT}$，检测到目标数为$ T_{\rm NC}$，误检像素数为$ P_{\rm NF}$，图像总像素数为$ P_{\rm NA}$，分别计算检测率和虚警率

$ DR = \frac{{{T_{{\rm{NC}}}}}}{{{T_{{\rm{NT}}}}}};FA = \frac{{{P_{{\rm{NF}}}}}}{{{P_{{\rm{NA}}}}}} $

(17)

以虚警率$ FA$为横坐标，检测率$ DR$为纵坐标，通过改变阈值大小，得到一组图像检测结果的ROC曲线，可知，当ROC曲线越快趋近于顶线$ DR$=1，曲线下面积越大，检测结果就越好。

为了验证算法的有效性，用4组红外图像序列，第1组图像分辨率为768×576像素，一共60帧，60个目标，第2组图像分辨率为768×576像素，一共60帧，60个目标，第3组图像分辨率为768×576像素，一共60帧，60个目标，第4组图像分辨率为768×576像素，一共60帧，60个目标。4组图像序列一共240张图。

本文选择经典的Max-mean和Max-median算法作^[1]为对比算法^[18]，每种算法的检测结果如图 5所示。其中，第1行是分别从4组图像序列中随机抽取的一张红外图像原图，第2行是Max-mean算法对这4张红外图像滤波后的结果图，第3行是Max-median算法对这4张红外图像滤波后的结果图，第4行是本文算法对这4张红外图像处理后的结果图。从图 6中可以更加清晰地看到每种算法的检测结果。其中，第1行是这4张红外图像的3维图，第2行是Max-mean算法对这4张红外图像滤波后的3维图，第3行是Max-median算法对这4张红外图像滤波后的3维图，第4行是本文算法对这4张红外图像处理后的3维图。可以看出，Max-mean和Max-median算法对红外图像中的边缘非常敏感，第3张和第4张红外图像在很大程度上被来自相机和其他原因的噪声污染，而Max-mean和Max-median算法对这些噪声也异常敏感，所有噪声点几乎都没有滤除，从检测结果的3维图像上可以明显地看到，本文算法的检测结果非常显著，对噪声不敏感，能够实现非常好的目标检测效果，这说明本文算法与传统算法相比具有一定优势。

分别绘制这Max-mean和Max-median两种算法以及本文算法ROC曲线，如图 7所示。从前述分析可知，当ROC曲线越快趋近于顶线$ DR$=1，曲线下面积越大，检测结果就越好。从检测结果来看，第1~3组图像序列中，本文算法的ROC曲线下面积明显较其他两种算法更大，第4组图像序列中3种算法的检测结果ROC曲线在横坐标上已经进入10^-8数量级，这说明3种算法在该组图像序列中都获得了非常好的检测效果，且3种算法的检测结果差异很小，这是由第4组图像序列的背景部分比较简单所造成的，由此可见，本文算法在处理简单背景的红外图像时，其ROC曲线近似于阶跃函数，说明其具有十分优秀的检测效果。从总体检测结果上来看，本文算法在红外小目标上检测效果相比于传统方法是有明显改善地。

本文算法通过Facet模型提取图像方向导数信息，将红外图像方向导数特征与稀疏表示有关理论相结合，通过对单一方向上小目标特性进行详细分析，并将其扩展到多个方向上呈现的特征信息，讨论了目标和背景之间存在的区别，以稀疏表示理论为媒介，获得了小目标最终检测结果。实验结果表明，本文算法具有较高检测精度和较强抗噪声干扰能力，相比于传统检测算法具有一定优势，提高了检测率，降低了虚警率。

本文算法另一大特点在于，可根据不同情况下的检测背景，训练出不同背景字典，从而得到更好的检测效果，在实际应用中具有良好针对性。