CN111046800A - 一种基于低秩与稀疏分解的高光谱图像异常目标检测方法 - Google Patents

Abstract

本发明属于遥感影像处理技术领域，公开了一种基于低秩与稀疏分解的高光谱图像异常目标检测方法，包括：S1，将高光谱图像转换成二维矩阵；S2，利用图像中背景的低秩特性和目标的稀疏特性，对所述二维矩阵进行低秩与稀疏矩阵分解，得到低秩矩阵和稀疏矩阵；S3，针对所述低秩背景矩阵，计算其均值，然后计算全局协方差矩阵；S4，利用双滑动窗遍历所述高光谱图像，得到检测结果。本发明体现了简单、快速和高精度的特点，对采用高光谱技术在进行质量检测、环境监测以及军事侦察等方面具有很好的实用价值。

Description

一种基于低秩与稀疏分解的高光谱图像异常目标检测方法

技术领域

本发明方法涉及图像异常检测技术领域，具体地，涉及一种基于低秩与稀疏分解的高光谱图像异常目标检测方法。

背景技术

高光谱图像异常目标检测是高光谱图像处理领域的热门和关键技术，其实质是一种二分类技术，即将高光谱图像分为背景和与背景有差异的异常两部分。它能够在不需要任何先验信息的情况下有效检测出高光谱图像中的异常目标。由于其无应用条件限制的特点，因此广受关注，在环境监测、矿物探测、食品质量检测和军事侦察方面有巨大的应用价值和广阔的应用前景。

目前，高光谱图像异常目标检测领域最经典的算法是Reed和Xiaoli Yu提出的RX算法。RX算法是在高光谱图像背景的概率密度函数服从高斯分布的假设基础上推导而来的，通过计算高光谱图像中每一像元与背景分布的马氏距离判断该像元是否为异常目标。由于缺乏先验信息，背景的分布特征(均值和协方差矩阵)通常由整幅图像作为数据集计算而得。然而，由于实际场景的复杂性和多样性，高光谱图像的背景往往不能完全符合高斯分布的假设，因此会产生较大的误差，导致检测精度不高。虽然RX算法的改进算法有很多，如SSRX、UTD和KRX等，但精度提升不高。

基于表示的方法被越来越多的应用到高光谱异常目标检测领域，基于协同表示的高光谱图像异常目标检测算法CRD即是其中的典型代表。CRD的核心思想是：如果一个像元属于背景，那么它就能够由其周围像元的线性组合近似表示，否则，即可判定该像元为异常。计算待测像元与其周围像元的欧式距离作为权重，利用待测像元与其周围像元线性组合差值的二范数以及权重与线性组合系数乘积的二范数构建目标函数，通过最小化目标函数求解确定线性组合系数，并以之重建高光谱图像，计算待测像元的原始高光谱图像数据与重建后的高光谱图像数据之差，其二范数即可作为检测结果。然而，当背景复杂或干扰大时，容易将背景误检为异常，且当待检测像元周围有异常时，算法的输出结果有可能会变得很小，导致检测精度下降。

近年来，基于高光谱图像背景低秩特性的方法开始用于检测异常，其中张玉香等人提出的LSMAD算法取得了较好的检测结果。LSMAD首先将高光谱图像转换为二维矩阵，然后将其分解为低秩矩阵、稀疏矩阵和噪声矩阵三个部分。其中低秩矩阵可被视为背景矩阵，根据低秩矩阵计算均值和协方差矩阵，而后进行基于马氏距离的检测。此方法中，低秩和稀疏矩阵的分解通过GoDec算法实现。此方法需要设置的参数较多，且在检测精度和运行时间上仍有改进的空间。

发明内容

针对现有技术的问题，本发明提出了一种基于低秩与稀疏分解的高光谱图像异常目标检测方法，可用于快速有效地检测高光谱图像中的异常目标。所述方法包括：S1，将高光谱图像转换成二维矩阵；S2，利用图像中背景的低秩特性和目标的稀疏特性，对所述二维矩阵进行低秩与稀疏矩阵分解，得到低秩矩阵和稀疏矩阵；S3，针对所述低秩背景矩阵，计算其均值，然后计算全局协方差矩阵；S4，利用双滑动窗遍历所述高光谱图像，得到检测结果。

进一步，在S4中，对每一个待测像元，以其为中心建立双滑动窗，以外窗和内窗之间的像元作为数据集求均值，并结合步骤S3所得的协方差矩阵求其与待测像元的马氏距离作为检测结果。

进一步，当检测完所述高光谱图像中的每一个像元后，将每一个像元的检测结果作为灰度值，合成灰度图像，选择设定的阈值进行二值分割，所得结果为最终的异常检测结果。

进一步，在S2中，采用SSGoDec算法进行低秩与稀疏矩阵分解。

进一步，步骤S2包括：S21，根据高光谱图像数据选择合适的背景秩上界参数r，最大迭代次数Iter和误差阈值tol；S22，对所述低秩矩阵和所述稀疏矩阵进行初始化；S23，对所述低秩矩阵进行循环迭代；S24，当满足迭代结束条件时，将最后一次迭代中的结果作为最终结果输出，得到低秩矩阵。

进一步，在S23中，当所述二维矩阵和所述低秩矩阵和所述稀疏矩阵之间的噪声小于所述误差阈值或达到最大迭代次数，则结束循环迭代。

进一步，在S1中，所述高光谱图像转换成N×p的二维矩阵后由下式表示：

X＝B+S+N

其中，

是低秩矩阵，可以视为低秩背景，

为稀疏矩阵，可以表示稀疏异常部分，

则代表图像中存在的噪声。

进一步，在S22中，初始化包括：令迭代次数t＝0，B⁽⁰⁾＝X，S⁽⁰⁾＝0，A₁＝randn(p,r)，其中上标表示迭代次数，randn是一个随机函数，A₁是使用此随机函数生成的p×r的随机矩阵。

进一步，在S23中，通过如下方式进行循环迭代，

迭代次数t＝t+1，计算：

Y₁＝(X-S^(t-1))A₁

A₂＝Y₁

Y₂＝(X-S^(t-1))^TA₂

判断

是否成立，如果成立，令

并返回S21，计算本次迭代中的低秩矩阵和稀疏矩阵，其中λ＝8，

进一步，在S3中，通过如下方式计算全局协方差矩阵：

均值：

协方差矩阵：

本发明方法的有益效果是，本发明可用于从高光谱图像中检测出与背景有差异的异常点目标和小区域目标。本发明采用新的低秩与稀疏分解方法，提高了获取低秩背景矩阵的运算效率以及精度。而且，该方法参数较少且设置简单。同时，利用双滑动窗策略，使用全局协方差矩阵计算马氏距离，简化了计算流程，因此可以快速地获取高精度的异常检测结果。本发明操作简单、方法可靠、适用性好、实用性强，所得实验结果精度高，应用前景广阔。

附图说明

为了更容易理解本发明，将通过参照附图中示出的具体实施方式更详细地描述本发明。这些附图只描绘了本发明的典型实施方式，不应认为对本发明保护范围的限制。

图1是本发明方法的一个实施方式的流程图。

图2是本发明方法的一个实施方式的流程图。

图3是双滑动窗示意图。

图4是PaviaC高光谱图像，其中(a)为PaviaC高光谱数据第70波段的灰度图，(b)为异常点的位置。

图5是本发明的方法和对比方法的异常检测结果，其中(a)为RX检测结果，图(b)为LSMAD检测结果，(c)为本发明的方法的检测结果。

图6是本发明的方法和对比方法的ROC曲线。

图7显示了三种方法的AUC值。

具体实施方式

下面参照附图描述本发明的实施方式，其中相同的部件用相同的附图标记表示。在不冲突的情况下，下述的实施例及实施例中的技术特征可以相互组合。

图1-2显示了本发明方法的流程图。本发明的方法包括以下步骤：将高光谱图像数据转换为二维矩阵；进行低秩与稀疏矩阵分解；计算协方差矩阵；利用滑动窗遍历整幅图像，得到检测结果。下面进行详细说明。

参照图1，在S1中，将高光谱图像数据转换为二维矩阵。

对大小为m×n×p的高光谱图像X，将其转换为二维N×p的矩阵X。其中，m表示高光谱图像空间维中行的总数，n表示高光谱图像空间维中列的总数，p表示高光谱图像中所有光谱波段的总数，N表示高光谱图像中总的像元数。在转换所得的矩阵中，每一行代表高光谱图像中的每一个像元，列中的元素即为该像元在每一波段中的反射值。按照一定的顺序排列，例如：对于一个高光谱图像，先提取中的第一行第一列的像元，将其转换成列向量，而后再提取第二行第一列的像元将其转换成列向量，并将其置于第一行第一列的像元列向量之后，作为整个矩阵的第二列，以此类推，当转换完第一列的像元之后，再对第二列的像元进行转换，转换所得列向量依次置于前一个列向量的后面，作为矩阵对的一列，直到完成图像中所有像元的转换，即可得到一个N×p的二维矩阵。

再次参照图1，在S2中，对S1中得到的所述二维矩阵进行低秩与稀疏矩阵分解。

采用SSGoDec算法对S1中得到的高光谱图像二维矩阵进行分解，将其分解为低秩矩阵、稀疏矩阵和噪声矩阵：

X＝B+S+N

其中，

是低秩矩阵，可以视为低秩背景，

为稀疏矩阵，可以表示稀疏异常部分，

则代表图像中存在的噪声。

具体地，所述分解过程可以包括如下步骤：

S21，设置参数：根据高光谱图像数据选择合适的背景秩上界参数r，最大迭代次数Iter和误差阈值tol。其中秩上界参数可利用已有算法进行选择，迭代次数和误差阈值按算法的默认值。

S22，对所述低秩矩阵B和所述稀疏矩阵S进行初始化。

令迭代次数t＝0，B⁽⁰⁾＝X，S⁽⁰⁾＝0，A₁＝randn(p,r)

上标表示迭代次数，randn是一个随机函数，A₁是使用此随机函数生成的p×r的随机矩阵。

S23，对所述低秩矩阵B进行循环迭代，当所述二维矩阵和所述低秩矩阵B和所述稀疏矩阵之间的噪声小于所述误差阈值或达到最大迭代次数，则结束循环迭代。

通过推导这些公式，使求解所得矩阵B的秩小于r，且矩阵S具有稀疏特性。低秩表示的是图像中的背景通常仅由少数几种地物构成，而稀疏则表征的是图像中目标所占像素很少，在位置上分布稀疏。

对于第t次迭代，判断

是否成立，或t是否大于于最大迭代次数Iter，若二者有一个成立，则结束循环迭代，否则，进入循环。

迭代次数t＝t+1，计算：

Y₁＝(X-S^(t-1))A₁

A₂＝Y₁

Y₂＝(X-S^(t-1))^TA₂

判断

是否成立，如果成立，令

并返回S21。

计算本次迭代中的低秩矩阵和稀疏矩阵，其中λ＝8，

该公式(1)基于BRP(bilateral random projection)的低秩近似理论所得，其目的是当S确定时，求得满足秩小于r且

最小的矩阵B。

该公式(2)基于增广拉格朗日乘子法推导而得，利用软阈值法求解，其目的是当B确定时，求得

最小的矩阵S，其中‖S‖₁确保了S具有稀疏性。

上述两个公式(1)和(2)即是如下两个子问题的求解结果。

S24，输出结果：当满足迭代结束条件时，将最后一次迭代中的结果作为最终结果输出，得到低秩矩阵B，作为低秩背景矩阵。

再次参照图1，在S3中，计算协方差矩阵：针对步骤S2所得的低秩背景矩阵B，首先计算均值，而后计算全局协方差矩阵：

再次参照图1，在S4中，利用双滑动窗遍历所述高光谱图像，得到检测结果。

对每一个待测像元，以其为中心建立双滑动窗，以外窗和内窗之间的像元作为数据集求均值，并结合步骤S3所得的协方差矩阵求其与待测像元的马氏距离，计算所得马氏距离即为检测结果。

对像元

马氏距离的计算公式为：

其中，μ_local为以该像元为中心的滑动窗中内外窗之间像元的均值。

当滑动窗遍历整幅图像，即检测完图像中的每一个像元后，将每一个像元的检测结果作为灰度值，合成灰度图像，选择合适的阈值进行二值分割，所得结果即为最终的异常检测图。

下面参照图3-7，以ROSIS传感器在意大利北部Pavia城市获取的高光谱图像为待检测图像进行本发明的方法的说明。此传感器的空间分辨率为1.3m，图像大小为108×120个像素，包含430-860nm的102个波段。

首先将此高光谱图像转换为12960×102的二维矩阵，矩阵的每一行代表图像中的一个像元。设置低秩与稀疏矩阵分解所需要的参数，r＝3，Iter＝200，tol＝0.001。

根据分解所得低秩背景矩阵，计算全局协方差矩阵。而后利用双滑动窗进行检测，双滑动窗的构造如图3所示。内外窗的尺寸分别为7×7和9×9，利用外窗与内窗间的像元计算均值，再求对应的马氏距离。

下面结合仿真实验对本发明的效果做进一步的说明。

仿真条件：仿真实验是在主频2.60GHz的Intel(R)Core(TM)i7CPU、内存为8GB的硬件环境和MATLAB R2018a的软件环境下进行的。

仿真内容与结果分析：该仿真实验是采用本发明方法和现有技术的RX和LSMAD高光谱异常检测方法进行对比仿真验证。

参照图4，图4(a)为来自于ROSIS高光谱遥感卫星采集的真实高光谱图像，为PaviaC高光谱数据第70波段的灰度图。图4中的(b)中的白色亮点表示真实高光谱图像中的异常点。

图5显示了本发明的方法和对比方法的异常检测结果，图5中的(a)为RX检测结果，图5中的(b)为LSMAD检测结果，图5中的(c)为本发明的方法的检测结果。比较图5中的三幅检测图可以看出：三种检测方法虽然都有很好的检测效果，但是与本发明相比，另外两种常用高光谱异常检测方法将多个背景像元误检为异常，使得虚警率上升。因此本发明比现有技术的RX方法和LSMAD方法检测效果更好。

为了更精确的比较三种检测方法的优劣，选择ROC曲线和AUC值作为评价标准。ROC的横坐标为虚警率，采用对数坐标轴，纵坐标为检测率。显而易见，ROC曲线的位置越靠近坐标空间的左上，代表其对应方法的检测性能越好。图6显示了本发明的方法和对比方法的ROC曲线，由图6可见，本发明方法的ROC曲线比另外两种对比方法的ROC曲线都要靠近左上方，说明本发明的方法在三种方法中效果是最好的。为定量比较三种方法的优劣，采用ROC曲线与横坐标之间的面积AUC值进行衡量。

从图7可见，本发明与现有技术的RX和LSMAD高光谱异常检测方法相比，使用本发明方法得到的异常检测效果更好，精度更高。

以上所述的实施例，只是本发明较优选的具体实施方式，本领域的技术人员在本发明技术方案范围内进行的通常变化和替换都应包含在本发明的保护范围内。

Claims (10)

1.一种基于低秩与稀疏分解的高光谱图像异常目标检测方法，其特征在于，包括：

S1，将高光谱图像转换成二维矩阵；

S2，利用图像中背景的低秩特性和目标的稀疏特性，对所述二维矩阵进行低秩与稀疏矩阵分解，得到低秩矩阵和稀疏矩阵；

S3，针对所述低秩背景矩阵，计算其均值，然后计算全局协方差矩阵；

S4，利用双滑动窗遍历所述高光谱图像，得到检测结果。

2.根据权利要求1所述的基于低秩与稀疏分解的高光谱图像异常目标检测方法，其特征在于，

在S4中，对每一个待测像元，以其为中心建立双滑动窗，以外窗和内窗之间的像元作为数据集求均值，并结合步骤S3所得的协方差矩阵求其与待测像元的马氏距离作为检测结果。

3.根据权利要求2所述的基于低秩与稀疏分解的高光谱图像异常目标检测方法，其特征在于，

当检测完所述高光谱图像中的每一个像元后，将每一个像元的检测结果作为灰度值，合成灰度图像，选择设定的阈值进行二值分割，所得结果为最终的异常检测结果。

4.根据权利要求1所述的基于低秩与稀疏分解的高光谱图像异常目标检测方法，其特征在于，

在S2中，采用SSGoDec算法进行低秩与稀疏矩阵分解。

5.根据权利要求4所述的基于低秩与稀疏分解的高光谱图像异常目标检测方法，其特征在于，步骤S2包括：

S21，根据高光谱图像数据选择合适的背景秩上界参数r，最大迭代次数Iter和误差阈值tol；

S22，对所述低秩矩阵和所述稀疏矩阵进行初始化；

S23，对所述低秩矩阵进行循环迭代；

S24，当满足迭代结束条件时，将最后一次迭代中的结果作为最终结果输出，得到低秩矩阵。

6.根据权利要求5所述的基于低秩与稀疏分解的高光谱图像异常目标检测方法，其特征在于，

在S23中，当所述二维矩阵和所述低秩矩阵和所述稀疏矩阵之间的噪声小于所述误差阈值或达到最大迭代次数，则结束循环迭代。

7.根据权利要求5所述的基于低秩与稀疏分解的高光谱图像异常目标检测方法，其特征在于，

在S1中，所述高光谱图像转换成N×p的二维矩阵后由下式表示：

X＝B+S+N

其中，

是低秩矩阵，可以视为低秩背景，

为稀疏矩阵，可以表示稀疏异常部分，

则代表图像中存在的噪声。

8.根据权利要求7所述的基于低秩与稀疏分解的高光谱图像异常目标检测方法，其特征在于，在S22中，初始化包括：

令迭代次数t＝0，B⁽⁰⁾＝X，S⁽⁰⁾＝0，A₁＝randn(p，r)

上标表示迭代次数，randn是一个随机函数，A₁是使用此随机函数生成的p×r的随机矩阵。

9.根据权利要求8所述的基于低秩与稀疏分解的高光谱图像异常目标检测方法，其特征在于，在S23中，通过如下方式进行循环迭代，

迭代次数t＝t+1，计算：

Y₁＝(X-S^(t-1))A₁

A₂＝Y₁

Y₂＝(X-S^(t-1))^TA₂

判断

是否成立，如果成立，令

并返回S21，

计算本次迭代中的低秩矩阵和稀疏矩阵，其中λ＝8，

10.根据权利要求6所述的基于低秩与稀疏分解的高光谱图像异常目标检测方法，其特征在于，在S3中，通过如下方式计算全局协方差矩阵：

均值：

协方差矩阵：

Images