CN103593852A - 基于同质图斑的高光谱影像异常探测方法 - Google Patents

Abstract

本发明公开了一种基于同质图斑的高光谱异常探测方法，包括步骤：S1、采用光谱降维法对原始高光谱数据X进行变换，得到变换后的特征向量Y；S2、采用多光谱影像分割方法对降维后的影像Y进行分割，得到标记影像L；S3、搜索标记影像L中标记号相同的像素，形成像素集，每一个像素集对应着一个同质斑块；S4、采用加权法计算每个同质斑块的均值及方差；S5、以同质斑块为单位，依次计算每个像素的异质度；S6、确定异质度的阈值，依次判断每个像素是否为异常点。本发明克服了背景信息统计时背景种类复杂、维数多以及目标参与计算等缺点，提高了背景信息计算精度及正态分布函数与背景分布模型的符合度，提高了在低虚警率条件下检测概率。

Description

基于同质图斑的高光谱影像异常探测方法

技术领域

本发明涉及工业自动化检测技术领域，特别涉及一种基于同质图斑的高光谱影像异常探测方法。

背景技术

高光谱遥感以其纳米级的光谱分辨率和在可见光到短波红外甚至热红外波段的超多波段成像特点，引起广大从事遥感和应用人士的高度重视。其应用领域也越来越广泛，涉及全球环境、土地利用、资源调查等许多领域。特别是在目标探测方面，高光谱尤其重要。由于受地面采样距离的限制,在很多情况下感兴趣的目标都以亚像元的形式存在，亚像元是高光谱影像上目标存在的主要方式。此时,利用传统的遥感手段很难对这些小目标进行有效探测,高光谱遥感以其波段多、光谱分辨率高、信息量大的特点为亚像元方式存在目标的检测提供了可能。

目前存在多种高光谱影像目标检测方法，依据先验知识的有无大致可分为四种：（1）目标已知、背景未知；（2）目标已知、背景已知；（3）目标未知、背景已知；（4）目标未知、背景未知。在很多实际情况下，研究者总没有足够的先验知识来表征目标类别及背景的信息。因此，只能通过对原始数据进行统计分析，并依据目标在特征空间的特性进行目标探测，这一类方法称之为异常检测。这一类方法首先统计背景的统计信息，如均值、协方差矩阵、相关矩阵等，然后设计一个滤波器，并利用滤波器探测目标。这一类的算法有RX算法及其一系列改进算法、低概率目标探测算法（LPTD）、均衡目标探测算法（UTD）等。

综上所述，这些方法在探测目标时最重要的是计算背景的统计信息，且一般用正态分布模型来对背景数据进行模拟。但是一幅光谱影像上地物种类繁多，背景数据有时并不能很好的符合正态分布模型。且因为目标的存在，更有可能使得背景模型存在较大的误差。如此一来，基于正态分布模型的这些异常探测方法所获得的结果就会存在虚警率高、探测效率低等缺陷。

发明内容

（一）要解决的技术问题

本发明要解决的技术问题是：如何提供一种基于同质图斑的高光谱影像异常提取方法，解决背景信息统计时背景种类复杂、目标参与背景特征计算、虚警率高以及探测效率低的问题。

（二）技术方案

为解决上述问题，本发明提供一种基于同质图斑的高光谱异常探测方法，包括：步骤S1、采用光谱降维法对原始高光谱数据X进行变换，得到变换后的特征向量Y；

S2、采用多光谱影像分割方法对降维后的影像Y进行分割，得到标记影像L；

S3、搜索标记影像L中标记号相同的像素，形成像素集，每一个像素集对应着一个同质斑块；

S4、采用加权法计算每个同质斑块的均值及方差；

S5、以同质斑块为单位，依次计算每个像素的异质度；

S6、确定异质度的阈值，依次判断每个像素是否为异常点。

优选地，在步骤S1中变换方法具体为：

Y=W_PCAX    (1)

式中:W_PCA为PCA的投影矩阵，其由光谱数据协方差矩阵的特征值按大小排列后，前m个特征值所对应的列特征向量组成的矩阵。

优选地，在步骤S2中，采用基于背景梯度分布模型的分水岭影像分割方法对多光谱影像Y进行分割，得到标记图像L。

优选地，所述步骤S2具体包括：

其中，影像的梯度影像IG具体为：

$\mathrm{IG}=\underset{i}{\max }\left\{{\mathrm{IG}}_i\right\}---\left(2\right)$

式中，

为每一个波段的梯度影像，其中

水平方向的梯度算子及竖直方向梯度算子分别为 ${\mathrm{GD}}_x=-\mathrm{xexp}\left(\frac{x^2+y^2}{2{\mathrm{\σ}}^2}\right)/\left({\mathrm{\π\σ}}^2\right)$ 及 ${\mathrm{GD}}_y=-y\·\exp \left(\frac{x^2+y^2}{2\·{\mathrm{\σ}}^2}\right)/(\mathrm{\π}\·{\mathrm{\σ}}^2),$ σ为高斯函数的方差，x、y为影像的行列坐标值；

${\mathrm{IG}}_B=\mathrm{IG}\⊗\mathrm{LD}---\left(4\right)$

式中，为高斯低通滤波算子；

IM=1，if(IG>IG_B)    (5)

其中，IM为标记影像(当梯度值大于梯度趋势影像值时为1，否则为0)，IG为梯度影像，IG_B为梯度趋势影像，采用分水岭影像分割方法进行影像分割，得到最终的标记影像L。

优选地，步骤S3中，采用种子生长法搜索出每个标记对应的像素集L_i，i为标记号。

优选地，在步骤S4，采用加权法计算每个标记对应集合的均值向量μ_i及方差矩阵Σ_i，如式（6）

$\begin{array}{c}{\mathrm{\μ}}_i=\underset{j\∈L_i}{\mathrm{\Σ}}{w}_{\mathrm{ij}}{X}_{L_{\mathrm{ij}}}\\ {\mathrm{\Σ}}_i=\underset{j\∈L_i}{\mathrm{\Σ}}{w}_{\mathrm{ij}}{(X_{L_{\mathrm{ij}}}-{\mathrm{\μ}}_i)}^T(X_{L_{\mathrm{ij}}}-{\mathrm{\μ}}_i)\end{array}---\left(6\right)$

式中，w_ij为同质斑块i第j像素对应的权值，该值定义为该点对应特征向量与该区域均值向量的欧式距离的倒数,

为第i个斑块第j个像素的特征矢量，L_i为第i个标记对应斑块的像素集合。

优选地，在步骤S5中，依次计算同质块中每个像素的异质度，计算公式如式（7）

$D\left(x\right)={(x-u_i)}^T{(\frac{M}{M+1}{\mathrm{\Σ}}_i+\frac{1}{M+1}{(x-u_i)}^T(x-u_i))}^{-1}(x-u_i)---\left(7\right)$

其中x为像素的特征矢量，M为第i标记对应斑块包含像素的个数，μ_i及Σ_i为标记i对应斑块的均值向量及协方差矩阵。

优选地，在步骤S6中，确定异质度阈值T，对于每一个像素若其异质度大于T，则认为该点为异常点，否则为背景点。

优选地，所述阈值确定具体包括：统计异质度图像的直方图，从而获得异质度值的概率分布函数以及概率分布函数，利用预先确定的虚警率确定阈值T。

（三）有益效果

本发明利用图像分割技术对降维后主要包含背景信息的光谱数据分割，形成同质图斑，并以此图斑为基础，计算背景统计信息以及计算每个像素的异质度，进而检测出异常点。克服了背景信息统计时背景种类复杂、维数多以及目标参与计算等缺点，提高了背景信息计算精度及正态分布函数与背景分布模型的符合度，提高了在低虚警率条件下检测概率。

附图说明

图1为依照本发明实施例的基于同质图斑的高光谱异常探测方法流程示意图；

图2为依照本发明实施例的基于梯度分布多光谱影像分割方法流程示意图。

具体实施方式

下面结合附图和实施例，对本发明的具体实施方式作进一步详细描述。以下实施例用于说明本发明，但不用来限制本发明的范围。

本发明所提供的基于同质图斑的高光谱异常检测方法的流程图，如图1-2所示，主要包括以下步骤：

步骤101对高光谱影像进行光谱降维。

首先计算将3维的高光谱数据变为2维，其中一维为光谱维，另外一维为像素的序号维，每个像素的序号定义为Id_ij=i×col+j,i、j分别为像素在二维影像上的行与列，col则为影像的列数。记高光谱影像X变为二维后为X2，其维数为p×n、p为波段数，n为像素个数。

接着计算Y2的协方差矩阵，计算公式为(均值向量为

)

$\mathrm{Cov}=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{({X2}_i-u)}^T({X2}_i-u)$

对协方差矩阵Cov进行奇异值分解，有Cov=UΣV^T。取U前m个列向量，得矩阵U_m=U(:,1:m),因此主成分变换矩阵

，这里取m为6。

然后对矩阵X2进行变换，有Y2=W_PCAX2,将2维矩阵X2变换为三维矩阵Y,即完成了对高光谱数据的降维处理。

步骤102对多光谱影像进行分割。

首先采用高斯梯度算子计算影像的梯度影像，每一像素对应有多个波段的梯度值，取最大的梯度值作为该像素的梯度，如式

$\mathrm{IG}=\underset{i}{\max }\left\{{\mathrm{IG}}_i\right\}$

式中，

为每一个波段的梯度影像，其中

水平方向的梯度算子及竖直方向梯度算子分别为 ${\mathrm{GD}}_x=-\mathrm{xexp}\left(\frac{x^2+y^2}{2{\mathrm{\σ}}^2}\right)/\left({\mathrm{\π\σ}}^2\right)$ 及 ${\mathrm{GD}}_y=-y\·\exp \left(\frac{x^2+y^2}{2\·{\mathrm{\σ}}^2}\right)/(\mathrm{\π}\·{\mathrm{\σ}}^2),$ σ为高斯函数的方差，x、y为影像的行列坐标值。这里取方差σ为1.0。

接着采用高斯低通滤波算子对梯度影像进行滤波得到梯度趋势影像 ${\mathrm{IG}}_B=\mathrm{IG}\⊗\mathrm{LD}.$ 其中， $\mathrm{LD}=\exp \left(\frac{x^2+y^2}{2\·{{\mathrm{\σ}}_b}^2}\right)$ 为高斯低通滤波算子。为了得到梯度分布的整体趋势，方差σ_b取值一般较大，这里取其值为10.0。

然后将梯度影像IG与梯度趋势影像进行比较。若某一位置处，趋势影像IG_B的值小于梯度影像值时，则认为该点为候选标记点（记为1，否则为0），形成标记影像IM。

IM=1，if(IG>IG_B)

处理完所有点后，对候选标记点进行后处理，采用种子生长法搜索所有候选标记块，如果每块点数大于阈值，则保留此区域块，否则将该块内所有点的值在影像IM均设置为0。

最后以IM为标记影像，IG为梯度影像，采用分水岭影像分割方法进行影像分割，得到最终的标记影像L。

步骤103搜索标记影像中标记号相同的像素，形成同质斑块。

首先对标记影像L上的像素因此进行扫描，判断其是否已被标记为已处理像素，如果已被处理，则扫描下一个像素。否则，以该点为种子点。

接着，以上一步得到的种子点采用种子生长方法进行跟踪，跟踪所有标记号与种子点相同的像素点。

然后记录所有跟踪到的像素，形成同像素集（同质图斑），并将这些像素标记为已被处理。

继续对标记图像进行处理，直到所有像素都已被处理。

步骤104采用加权法计算每个同质图斑的均值及方差。

对于每一个同质图斑，采用如下步骤计算该图斑的均值及方差。

首先计算图斑内所有点的均值，如式

然后计算图斑内每个像元到均值矢量的距离， $d_{L_{\mathrm{ij}}}={(X_{L_{\mathrm{ij}}}-{\tilde{u}}_i)}^T(X_{L_{\mathrm{ij}}}-{\tilde{u}}_i).$

接着计算每个像元的权值，

最后按下式计算该图斑的均值及方差

${\mathrm{\μ}}_i=\underset{j\∈L_i}{\mathrm{\Σ}}{w}_{\mathrm{ij}}{X}_{L_{\mathrm{ij}}}$

${\mathrm{\Σ}}_i=\underset{j\∈L_i}{\mathrm{\Σ}}{w}_{\mathrm{ij}}{(X_{L_{\mathrm{ij}}}-{\mathrm{\μ}}_i)}^T(X_{L_{\mathrm{ij}}}-{\mathrm{\μ}}_i)$

步骤105以同质斑块为单位，依次计算每个像素的异质度。

对图像上每一个像素，其对应的特征向量为x，首先查找其所在的图斑，获取该像素对应背景的均值及方差参数u_i、Σ_i。

然后利用下列公式计算该像素的异质度

$D\left(x\right)={(x-u_i)}^T{(\frac{M}{M+1}{\mathrm{\Σ}}_i+\frac{1}{M+1}{(x-u_i)}^T(x-u_i))}^{-1}(x-u_i)$

步骤106确定异质度的阈值，依次判断每个像素是否为异常点。

异质度的阈值由预先给定的虚警率及异质度图像的直方图确定。

统计异质度图像的直方图。

统计直方图时，首先确定异质度的最大最小值，然后将最大最小值确定的范围等分为k个间隔，k值事先确定，这里设置为500。并将记录每个范围存在像素个数的数组nHist所有元素设置为0。

接着依次扫描影像上的每一个像素，计算该像素所在范围，然后将该范围对应数组nHist的元素加一。

其次统计异质度值的分布函数，也就是从nHist的第1个元素开始，将其一个元素的值加到后一个元素上，即nHist[i]=nHist[i-1],直到所有元素被处理完。

再次计算背景像素的个数，N_B=(1-α)×M×N，其中α为虚警水平，M×N为影像的尺寸。

然后从第0个元素开始，依次将nHist中每个元素与N_B进行比较，如果N_B大于该元素，则继续比较下一个元素；否则，停止比较，并将该元素对应值的范围中值作为判断像素点是否为异常点的阈值T。

确定阈值T后，就可以判断每个像素是否为异常点。将图像上每个像素点的异质度与T比较，如果小于T，则为背景点，否则为异常点。

通过结合附图对本发明具体实施例的描述，本发明的其它方面及特征对本领域的技术人员而言是显而易见的。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明技术原理的前提下，还可以做出若干改进和替换，这些改进和替换也应视为本发明的保护范围。

Claims (9)

1.一种基于同质图斑的高光谱异常探测方法，其特征在于，包括步骤：

S1、采用光谱降维法对原始高光谱数据X进行变换，得到变换后的特征向量Y；

S2、采用多光谱影像分割方法对降维后的影像Y进行分割，得到标记影像L；

S3、搜索标记影像L中标记号相同的像素，形成像素集，每一个像素集对应着一个同质斑块；

S4、采用加权法计算每个同质斑块的均值及方差；

S5、以同质斑块为单位，依次计算每个像素的异质度；

S6、确定异质度的阈值，依次判断每个像素是否为异常点。

2.如权利要求1所述的方法，其特征在于，在步骤S1中变换方法具体为：

Y=W_PCAX    (1)

式中:W_PCA为PCA的投影矩阵，其由光谱数据协方差矩阵的特征值按大小排列后，前m个特征值所对应的列特征向量组成的矩阵。

3.如权利要求1所述的方法，其特征在于，在步骤S2中，采用基于背景梯度分布模型的分水岭影像分割方法对多光谱影像Y进行分割，得到标记图像L。

4.如权利要求3所述的方法，其特征在于，所述步骤S2具体包括：

其中，影像的梯度影像IG具体为：

$\mathrm{IG}=\underset{i}{\max }\left\{{\mathrm{IG}}_i\right\}---\left(2\right)$

式中，

为每一个波段的梯度影像，其中

水平方向的梯度算子及竖直方向梯度算子分别为 ${\mathrm{GD}}_x=-\mathrm{xexp}\left(\frac{x^2+y^2}{2{\mathrm{\σ}}^2}\right)/\left({\mathrm{\π\σ}}^2\right)$ 及 ${\mathrm{GD}}_y=-y\·\exp \left(\frac{x^2+y^2}{2\·{\mathrm{\σ}}^2}\right)/(\mathrm{\π}\·{\mathrm{\σ}}^2),$ σ为高斯函数的方差，x、y为影像的行列坐标值；

${\mathrm{IG}}_B=\mathrm{IG}\⊗\mathrm{LD}---\left(4\right)$

式中，

为高斯低通滤波算子；

标记影像定义为

IM=1，if(IG>IG_B)    (5)

其中，IM为标记影像，当梯度值大于梯度趋势影像值时为1，否则为0，IG为梯度影像，IG_B为梯度趋势影像，采用分水岭影像分割方法进行影像分割，得到最终的标记影像L。

5.如权利要求1所述的方法，其特征在于，步骤S3中，采用种子生长法搜索出每个标记对应的像素集L_i，i为标记号。

6.如权利要求1所述的方法，其特征在于，在步骤S4，采用加权法计算每个标记对应斑块的均值μ_i及方差矩阵Σ_i，如式（6）

$\begin{array}{c}{\mathrm{\μ}}_i=\underset{j\∈L_i}{\mathrm{\Σ}}{w}_{\mathrm{ij}}{X}_{L_{\mathrm{ij}}}\\ {\mathrm{\Σ}}_i=\underset{j\∈L_i}{\mathrm{\Σ}}{w}_{\mathrm{ij}}{(X_{L_{\mathrm{ij}}}-{\mathrm{\μ}}_i)}^T(X_{L_{\mathrm{ij}}}-{\mathrm{\μ}}_i)\end{array}---\left(6\right)$

式中，w_ij为同质斑块i第j像素对应的权值，该值定义为该点对应特征向量与该区域均值向量的欧式距离的倒数,

为第i个斑块第j个像素的特征矢量，L_i为第i个标记对应斑块的像素集合。

7.如权利要求1所述的方法，其特征在于，在步骤S5中，依次计算同质块中每个像素的异质度，计算公式如式（7）

$D\left(x\right)={(x-u_i)}^T{(\frac{M}{M+1}{\mathrm{\Σ}}_i+\frac{1}{M+1}{(x-u_i)}^T(x-u_i))}^{-1}(x-u_i)---\left(7\right)$

其中x为像素的特征矢量，M为第i标记对应斑块包含像素的个数，μ_i及Σ_i为标记i对应斑块的均值向量及协方差矩阵。

8.如权利要求1所述的方法，其特征在于，在步骤S6中，确定异质度阈值T，对于每一个像素若其异质度大于T，则认为该点为异常点，否则为背景点。

9.如权利要求8所述的方法，其特征在于，所述阈值确定具体包括：统计异质度图像的直方图，从而获得异质度值的概率分布函数以及概率分布函数，利用预先确定的虚警率确定阈值T。

Images