一种基于带约束稀疏表示的高光谱图像异常检测方法与流程

本发明属于图像处理领域，涉及一种基于带约束稀疏表示的高光谱图像异常检测方法。

背景技术：

高光谱图像目标探测是高光谱遥感应用的重要方向之一，涵盖了环境检测、城市调查、矿物填图和军事侦察等诸多领域。异常检测不需要任何目标光谱信息，具有更为广阔的应用前景。在异常检测中，异常点往往通过测试像元光谱和它局部背景光谱的显著差异来进行检测。

经典的异常检测算法有reed-xiaoli(rx)检测器、核rx检测器、支持向量数据描述(supportvectordatadescription，svdd)检测器等等。然而，当局部背景中含有目标信号的时候，传统的异常检测算法性能急剧下降，其中以rx算法最为显著。近年来，稀疏表示(sparserepresentation，sr)已经广泛应用到高光谱图像目标探测。这些稀疏表示算法的一个主要优势就是它们不需要对背景的统计分布做出明确假设。但它们通常都是在设置一定稀疏度的条件下进行求解的，而这个稀疏度很难解释且不能事先获知。并且，稀疏度和重构误差之间脆弱的平衡关系一直没有得到有效解决。和稀疏表示不同的是，联合表示(collaborativerepresentation，cr)异常检测器利用背景字典中所有的原子对测试像元进行表示。

技术实现要素：

针对现有技术不足，为解决上述局部背景受目标信号污染的问题，以及稀疏表示中稀疏度难以确定的问题，本发明所采用的具体技术方案如下。

一种基于带约束稀疏表示的高光谱图像异常检测方法，包括以下步骤：

(s1)将高光谱图像线性规范化；

(s2)对于每一个测试像元，根据双窗模型提取局部背景字典；

(s3)根据局部背景字典，求解带约束稀疏表示模型，得到模型最优解一；

(s4)根据模型最优解一，将所有异常原子从局部背景字典中删除，得到新的背景字典；

(s5)根据新的背景字典，求解带约束稀疏表示模型，得到模型最优解二；

(s6)根据所述步骤(s5)中的模型最优解二，计算像元的检测值；

(s7)遍历整个高光谱图像，对高光谱图像的每一个像元计算检测值，输出这些检测值构成的图像，即异常检测图像。

优选地，所述步骤(s1)的具体过程为：根据下式将高光谱图像的每个像元线性规范化到0～1之间：

其中，hl,j,k表示高光谱图像在第l行、第j列、第k个波段的灰度值，表示高光谱图像在线性规划化后第l行、第j列、第k个波段的灰度值，min(h)表示高光谱图像中所有像元灰度值的最小值，max(h)表示高光谱图像中所有像元灰度值的最大值。

优选地，所述步骤(s2)中提取局部背景字典的具体过程为：

设y表示一个测试像元(测试像元为线性规范化后的高光谱图像中任一个像元)，测试像元y被内窗和外窗所包围，内窗大小为win×win，外窗大小为wout×wout，win、wout表示灰度值个数，win、wout为整数；a表示局部背景字典，由内窗和外窗之间像元的光谱向量构成，称之为训练样本，其中，m表示高光谱图像的波段数量；n为所选择的训练样本个数，则局部背景字典a表示为：其中ai为内窗和外窗之间像元的光谱向量。

优选地，所述步骤(s3)求解带约束稀疏表示模型具体为：

对于每一个测试像元y，根据smo(sequentialminimaloptimization，缩写为smo)算法求解下面二次规划模型，得到最优解一α^*：

s.t.e^tα＝1

0≤αi≤c,i＝1,...,n

其中，a表示局部背景字典，α是一个稀疏向量，它的第i个元素αi表示背景字典中第i个原子的丰度，e表示单位列向量，即大小为n×1、元素全为1的列向量；参数c为常量，表示元素αi的上限，n为所选择的训练样本总个数。

优选地，所述参数c的取值为：

对于所有待检测的异常像元，参数ν取值必须需要满足以下条件：

ν≥max(η1,η2/αa)，

其中η1＝na/n，

na表示局部背景字典a中异常原子的个数，表示局部背景字典a中第i个异常原子的异常丰度，n为所选择的训练样本个数，αa为测试异常像元y中异常点的异常丰度。在实际应用中，通常根据图像场景对所有待检测异常像元的η1和η2/αa进行估计，取它们的最大值来确定参数ν的取值。

为了对模型进行更为直观的解释，对引入参数ν，将参数c表示为进行详细说明。

根据线性光谱混合模型，每一个像元的光谱am可以表示为：

其中，aa表示纯像元异常光谱，αa表示它对应的丰度；表示第h个背景成分的纯像元光谱，表示它对应的丰度；p表示背景成分类型的数量，n表示误差项。

参数ν的设置和背景字典中异常点的比例相关，异常点占比有两种定义方式：

η1＝na/n

其中na是背景字典a中异常原子的个数，表示背景字典a中第i个异常原子的异常丰度。

要想从背景字典a中删除所有的异常原子，参数ν必须需要满足以下条件：

ν≥max(η1,η2/αa)

其中，αa为测试像元y中异常点的异常丰度。在实际应用中，考虑到高光谱图像中普遍存在的噪声，以及背景光谱的起伏变化，参数ν设置为稍大于整个图像的max(η1,η2/αa)。

优选地，所述步骤(s5)中求解带约束稀疏表示模型具体过程为：

对于每一个测试像元y，根据smo(sequentialminimaloptimization，简称smo)算法求解下面二次规划模型，得到最优解二

s.t.e^tα＝1

0≤αi≤1,i＝1,...,n

其中，表示新的背景字典，α是一个稀疏向量，它的第i个元素αi表示背景字典中第i个原子的丰度，e表示一个大小为n×1、元素全为1的列向量；n为所选择的训练样本总数。

优选地，设r表示测试像元y的检测值，则所述步骤(s6)中计算像元的检测值的公式为：

与现有技术相比，本发明所具有的有益效果为：1)本发明的带约束稀疏表示模型是一个非参数模型，它不需要任何背景统计信息。2)本发明的稀疏表示模型加入了“和为一”约束和非负约束，具有实际物理意义。同时，它删除了稀疏度上限约束，不需要设置稀疏度，提高了重构精度。3)本发明的稀疏表示模型引入了参数ν，它能够在局部背景被异常点污染的情况下检测亚像元异常目标。

附图说明

图1为本发明方法流程图；

图2为同心双窗模型示意图。

图3为实施例中sandiego高光谱图像数据；

图4为实施例中urban高光谱图像数据；

图5为对sandiego高光谱图像和urban高光谱图像分别采用本发明与现有技术中方法处理的roc曲线比较图。

具体实施方式

下面，结合附图和实施例对本发明作进一步说明。

如图1所示，为本发明方法流程图，一种基于带约束稀疏表示的高光谱图像异常检测方法，具体包括下述步骤：

第一步：对高光谱图像进行线性规范化。

根据最大最小规范化方法对高光谱图像的每一个像元进行线性规范化。根据下式将高光谱图像的每个像元hl,j,k线性规范化到0～1之间：

第二步：对于每一个测试像元，根据双窗模型提取局部背景字典。详细内容见参考文献1。参考文献1：s.matteoli,m.diani,andg.corsini,“atutorialoverviewofanomalydetectioninhyperspectralimages,”ieeeaerosp.electron.syst.mag.,vol.25,no.7,pp.5–28,jul.2010。对于每一个测试像元，利用图2所示的滑动双窗模型提取其局部背景，构造局部背景字典，使背景字典具有局部自适应性。稀疏表示异常检测模型采用滑动的同心双窗模型，图2给出了从高光谱数据立方体到同心双窗模型的构建过程。设y是一个未知的测试像元，它被内窗(大小为win×win)和外窗(大小为wout×wout)所包围，a是局部背景字典。由内窗和外窗之间像元的光谱向量构成，称之为训练样本，其中，m表示高光谱图像的波段数量，表示m维向量空间；为所选择的样本个数。因此，对于每一个测试像元字典a都适应于局部统计特性。

第三步：根据背景字典求解稀疏表示模型。将测试像元用局部背景字典进行稀疏表示，并且使稀疏表示原子的丰度和为一且大于零。同时，删除稀疏度上限约束，给每个原子的丰度加入一个上限约束，建立带约束稀疏表示异常检测模型，如下所示：

s.t.e^tα＝1

0≤αi≤c,i＝1,...,n

c＝1/νn

其中参数ν的值设置为大于整个图像的max(η1,η2/αa)，s.t.表示满足约束条件的数学符号,t表示向量或矩阵转置符号。对于每一个测试像元y，根据smo算法求解以上二次规划模型，得到最优解一α^*。具体求解方法见参考文献2和参考文献3。参考文献2：r.-e.fan,p.-h.chen,andc.-j.lin,“workingsetselectionusingsecondorderinformationfortrainingsvm,”journalofmachinelearningresearch,vol.6,pp.1889–1918,2005。参考文献3：j.c.platt,fasttrainingofsupportvectormachinesusingsequentialminimaloptimization.cambridge,ma:mitpress,1999。该算法的在不同实际应用场景下的计算复杂度为ο(n)～ο(n^2.2)，求解速度非常迅速。

第四步：将所有异常原子从背景字典中删除，构造新的背景字典。根据第三步的求解结果删除所有异常原子：将第三步求解结果中丰度达到上限的原子全部删除，从而将所有异常原子删除，得到新的背景字典。实施例中，即将丰度α^*为c的原子从背景字典a中删除，构造新的背景字典

第五步：根据新的背景字典求解稀疏表示模型。求解第三步的稀疏表示模型时，首先将原子丰度的上限约束αi≤c改为αi≤1，再利用新的背景字典进行求解该模型。对于每一个测试像元y，根据smo算法再次求解以下二次规划模型，得到最优解二

s.t.e^tα＝1

0≤αi≤1,i＝1,...,n

第六步：根据稀疏表示残差计算检测值。该模型的目标函数即为稀疏表示残差的一部分，可用来计算检测值。根据下式计算每一个测试像元的检测值：

其中为第五步求解的稀疏表示模型的目标函数值。

第七步：遍历整个线性规范化后的高光谱图像，对图像的每一个像元计算检测值，输出检测值图像。

相比于传统rx方法，本方法的带约束稀疏表示模型是一个非参数模型，它不需要任何背景统计信息。相比于传统的sr方法在固定的稀疏度下进行求解，本方法加入了和为一约束和非负约束，使得稀疏表示模型具有了物理意义，并提高了算法的检测能力。同时，它不对稀疏度进行约束，不需要设置稀疏度，提高了重构精度。本方法的稀疏表示模型引入了参数ν，它能够在局部背景被异常点污染的情况下检测亚像元异常目标。

采用两幅高光谱图像进行实验，sandiego高光谱图像和它的真实目标分布情况如图3所示，其中图3(a)表示伪彩色图像(根据专利申请相关规定要求，图中仅示出了灰度图像)，图3(b)表示目标真实位置图像。urban高光谱图像和它的真实目标分布情况如图4所示，其中图4(a)表示伪彩色图像(根据专利申请相关规定要求，图中仅示出了灰度图像)，图4(b)表示目标真实位置图像。评价指标采用接收机工作特性(receiveroperatingcharacteristic,roc)曲线和roc曲线下的面积对算法的性能进行评价。roc曲线定义为：给定一个阈值，计算检测值图像的检测概率(probabilityofdetection)和虚警概率(falsealarmrate)；遍历所有的阈值，得到一系列检测概率和虚警概率；以检测概率为纵坐标，以虚警概率为横坐标绘制得到roc曲线。auc值定义为roc曲线下方和坐标轴围成的区域的面积。

利用sandiego高光谱图像测试算法在不同的内窗win和外窗wout下的检测性能，如表1所示。内窗越小，背景字典受污染的程度越高。传统方法rx和sr的auc指标变化较大，特别是在内窗较小的情况下，它们的效果非常差。而本方法的auc值则变化不大，说明本方法对不同污染程度的背景都具有良好的检测性能。也验证了算法的鲁棒性。

将本方法和传统的rx，sr，cr，svdd等算法进行比较，它们的roc曲线如图5所示，其中图5(a)为sandiego数据；图5(b)为urban数据。auc值如表2所示。本方法的检测概率在大多数虚警率下比其它算法高。而且本方法的auc值大于其它算法，特别是比传统的rx算法和svdd算法高得多。所有算法的运行时间如表3所示。本发明方法的运行时间比svdd稍长，但检测性能却好得多。而且，和现有技术中的其它算法相比，本方法的运行时间相对较短。总的来说，本方法取得了非常优越的检测性能，同时还具有较低的算法复杂度。

表1各主流算法在sandiego高光谱图像上roc曲线下的面积(auc)

表2各主流算法在两个高光谱图像上roc曲线下的面积(auc)

表3各主流算法在两个高光谱图像上的运行时间(s)

尽管已经示出和描述了本发明的实施例，对于本领域的普通技术人员而言，可以理解在不脱离本发明的原理和精神的情况下可以对这些实施例进行多种变化、修改、替换和变型，本发明的范围由所附权利要求及其等同物限定。