基于选择性核主成份分析的高光谱图像异常点的检测方法

摘要

基于选择性核主成分分析的高光谱图像异常点的检测方法，它涉及的是高光谱图像分析检测的技术领域。它是为了解决现有高光谱图像检测技术存在不能高效的对异常点的进行特征提取，而产生较多虚警，及在有严重的背景干扰下无法对异常点有效检测的问题。它的步骤为：对数据进行归一化，并执行核主成分分析；在主分量内构造滑动窗；计算滑动窗内像素的三阶矩和四阶矩并与设定值比较；遍历主分量后，记值；所有主分量得到处理；选最大主分量；用RX算子对所选主分量进行异常点检测，输出检测结果。本发明能高效的对高光谱图像中异常点目标的特征进行提取和选择，而降低虚警率，并实现在有严重的背景干扰情况下正常对异常点进行检测。

说明书

技术领域

本发明涉及的是高光谱图像分析检测的技术领域。

背景技术

异常信号也称为突变信号，信号中的异常点及不规则的突变部分往往携带着非常重要的信息，它是信号的重要特征之一。在图像分析和处理中，异常点通常被定义为在较小的邻域内异于背景类别的部分，或是发生突变的部分，它是参照不同的背景模型定义的。异常点在图像中往往仅占据少量的像元，甚至是存在于一个像元之内。背景模型是通过从图像中选取大范围区域或者测试像元附近局部邻域的参考数据确定的，主要分为全局背景模型和局部背景模型。全局和局部异常性定义的主要区别在于，观测量来源于整个图像背景类范围还是来源于局部背景类范围。相应地，异常检测算子可以划分为全局异常检测算子和局部异常检测算子。

高光谱遥感技术(Hyperspectral Remote Sensing)是在20世纪80年代发展起来的一种全新的遥感技术。该技术利用星载或机载的成像光谱仪设备对地面进行成像，成像光谱仪在对目标的空间特征成像的同时，对每个空间像元经过色散形成几十个乃至几百个窄波段以进行连续的光谱覆盖，从而形成谱分辨率达到纳米数量级的遥感数据。这种数据由于光谱分辨率高，通常称为高光谱数据或高光谱图像。高光谱数据的光谱分辨率在10纳米左右，比多光谱(Multispectral)图像高出几十甚至上百倍。伴随着成像光谱技术的小断发展，局光谙数据已经被应用到了众多领域中。从民用领域广泛应用的环境监测、城市规划、农作物估产、洪涝灾害调查、国土资源调查，到军事领域的卫星侦察、目标检测识别等等。

高光谱图像的突出特点是在获得目标图像二维空间景像信息的同时，还可以获得高分辨率的一维表征其物理属性的光谱信息，即图谱合一。与多光谱图像相比，由于高光谱图像光谱分辨力的提高，使许多原先在多光谱下不能解决的问题，现在在高光谱下可以得到解决。通过处理高光谱图像中目标图像的空间特征和光谱特征，可以以较高的可信度辨别和区分地物目标。这对遥感图像军事侦察、真／假目标识别、战场态势评估等都具有重要应用意义和巨大的潜力。

具有突出特点和巨大应用潜力的高光谱图像数据、具有重要实际意义和坚实理论基础的异常检测技术和异常检测应用需求三者的结合，使得基于高光谱图像数据的异常检测逐步成为国际遥感领域研究的热点。

早在1990年，Reed和Xiaoli Yu就提出了用于遥感图像处理的多光谱异常检测算法，称为RX算法。RX算法是一个具有恒虚警率CFAR，constantfalse-alarm rate)的自适应异常检测算法。在该算法中，异常检测被描述为一个二值假设检验，分为存在目标和不存在目标两种情况，并采用多变量高斯统计模型来描述两种假设的概率分布。RX算法直接应用于高光谱图像时，主要存在三方面的问题：一是高数据维使得高斯分布不能够很好地描述地物真实分布情况；二是高光谱图像的高数据维和波段间的高阶相关性使得协方差计算困难；三是大量混合像元的存在严重影响了检测效果和异节目标的精确定位。此外，直接在原始高光谱图像上进行RX检测会大大增加检测步骤的计算量。针对RX算法中存在的部分问题，Susan等提出一种自适应聚类的异常目标检测算法，该算法利用三维高斯Markov随机场描述高光谱图像局部区域的背景统计模型，综合利用空间信息和光谱信息，并利用似然算法估计参数。但是，该方法的检测性能受限于异常目标的空间尺寸，只能检测到具有一定空间尺寸的异常目标，不能较好地检测亚像元的异常目标，并且统计模型参数的估计比较复杂。在异常检测基础上，Chang等人进一步地提出了对检测到的异常目标进行分类的思想，利用不同类目标的均值实现对异常目标的分类，并形成了一种实时的检测算法。Stellman等人也提出了一种有效的实时高光谱图像异常检测算法。Kraut等研究了一种更为一般的子空间异常检测模型，并构造了相应的广义检测似然比。在RX算法基础上，Kwon等人进一步提出了一个非线性的RX算法，该方法能够获得优于RX算法的检测性能。Meng等人利用全色图像和多光谱图像融合得到高空间分辨率的多光谱图像，然后作为图像源，对国际上现有的两种主要异常检测算法(高斯马尔可夫随机场方法和RX算法)的性能进行了比较分析，说明了这两种方法各自的不足。Bernhardt在他的文章里，研究了一种不同于传统多变量高斯分布的统计模型，称之为重尾(fat-tailed)分布统计模型，用它来描述异常检测和混合像素解译中的自然背景分布。  国内的高光谱图像处理研究大多集中在波段选择、特征提取和分类方面，在高光谱图像的目标检测、异常检测方面的研究相对较少。国内主要有哈尔滨工业大学、国防科技大学、北京航空航天大学和西北工业大学等单位开展了高光谱图像异常目标检测研究，提出了基于多波段图像融合、基于高阶统计量的异常检测算法等。

现有方法大多采用主成分变换和图像融合等来解决高数据维问题。这些方法在降低数据维的同时，能够在一定程度上获得较好的检测效果。但是，这些方法并不能获取对于异常检测来说最有效的信息，因此也不能获得令人满意的效果。此外，国外国内方法都未能有效地解决混合像元给检测带来的背景干扰问题。

综合国内外的现有技术来看，目前，迫切需要解决的问题(现有方法未能解决)主要包括两方面：

1、异常点特征提取和选择。高数据维数是高光谱图像异常检测中的一个难点。高数据维使得已有的检测方法对背景协方差的估计变得非常困难，所估计的协方差往往是异常的，使得后续检测容易产生较多虚警，且计算量较大。因此，对于高光谱图像异常检测来说，降低数据维数、挖掘异常点目标信息是需要重点解决的问题。有效的特征提取(或变换)方法能够降低原始数据维数，并挖掘出数据中包含的异常点目标信息；合理的选择包含异常目标的变换特征能够有效地提高检测性能。

2、严重的背景干扰。低空间分辨率等因素造成的混合像元给检测带来了严重的背景干扰，不同类别背景和异常目标的混合使得简单的局部高斯统计模型不能较好地描述实际分布情况，从而降低了已有检测算法的性能。因此，有效地去除背景干扰对降低检测方法的虚警来说是非常重要的。

发明内容

本发明的目的是为了解决现有高光谱图像检测技术存在不能高效的对异常点的特征进行提取和选择，而会产生较多虚警，且计算量较大，及在有严重的背景干扰下无法对异常点检测的问题。进而提出了一种基于选择性核主成份分析的高光谱图像异常点的检测方法。

本发明的高光谱图像异常点的检测方法步骤为：

步骤一、输入三维高光谱图像I(i，j，s)，其中，f＝1，2，…，P表示一个波段图像的行，j＝1，2，…，Q表示图像的列，s=1，2，…，N表示高光谱数据的波段数目，I(i，j，s)的大小为P×Q×N，按照波段图像方式，I(i，j，s)可以表示为，其中：

上述公式中，为第波段(i，j)点对应的灰度值；

步骤二、求高光谱图像I(i，j，s)的灰度最大值$I_{\max }=\underset{i,j,s}{\max }\left(I(i,j,s)\right)$，即P×Q×N个点所对应的灰度的最大值，将高光谱图像I(i，j，s)的每个点灰度值除以该最大值I_max，得到归一化的高光谱图像I(i，j，s)＝(i，j，s)/I_max；

步骤三、将归一化的高光谱图像I(i，j，s)转换成二维矩阵，即，对于某一个波段图像I_s，s＝1，2，3，…，N逐个波段图像按照由左至右、由上至下的顺序，将每一行像素值顺序放到一个行向量中为：

写成列向量形式，为：

这样得到的数据按列方向放置于二维矩阵中，得到大小为T×N的数据，其中T=P×Q，$y_s=I_{s\prime }^T$， s＝1，2，3，…，N；可以看作是由N个列向量y_s(s＝1，2，3，…，N)构成的样本矩阵；

步骤四、输入样本矩阵Y，计算核函数矩阵K，(K)_i，j＝(φ(y_i)φ(y_j))＝k(y_i，y_j)表示核函数矩阵K的每个元素，k(y_i，y_j)代表输入样本y_i，y_j(i，j＝1，2，…，l)之间的的核函数值，具体核函数可以选择为高斯核函数(K)_ij＝k(y_i，y_j)＝exp{(y_i-y_j)²/2σ²}(σ为常数)；

步骤五、对核函数矩阵进行归一化，即K_c＝K-1_N·K-K·1_N+1_N·K·1_N，

上述公式中，1_N是一个N×N的矩阵，该矩阵的元素均为1/N，N是样本个数；

步骤六、对归一化核函数矩阵K_c进行特征值分解，即求解如下特征值方程：

λ_kβ^k＝Kβ^k，

$({\beta }^k={({\beta }_1^k,{\beta }_2^k,···,{\beta }_l^k)}^T);$

步骤七、根据特征值分解得到的特征向量构成投影变换矩阵，即：

${\beta }^k{=({\beta }_1^k,{\beta }_2^k,···,{\beta }_l^k)}^T;$

步骤八、根据下式，将输入的二维形式的高光谱数据Y＝[y₁y₂…y_N]投影到特征向量${\beta }^k{=({\beta }_1^k,{\beta }_2^k,···,{\beta }_l^k)}^T$上，得到变换后的数据z；

$Z=(V^k·\phi \left(y\right))=\underset{i=1}{\overset{N}{\Sigma }}{\beta }_i^k\left(\phi \left(y_i\right)\phi \left(y\right)\right)=\underset{i=1}{\overset{N}{\Sigma }}{\beta }_i^{k}k(y_i,y),$

将z写成矩阵形式为z＝[z₁z₂…z_s…z_N]，上述公式中，

步骤九、按照步骤三的逆过程，将投影后的高光谱数据转换成图像波段形式，得到变换后的与原始输入高光谱图像具有相同数目的非线性变换主分量，即将转换成：

代表第s个变换后的非线性主分量，s＝1，2，…，N；

步骤十、令t＝0，f表示计数变量；

步骤十一、对第s个非线性主分量(s＝1，2，…，N)，构造同检测窗口大小一致的滑动窗w，滑动窗大小设定为，n₁×n₂；设滑动窗w内对应的像素灰度值逐行排列放置到一个行向量中，得到

步骤十二、按照如下公式计算滑动窗内像素的三阶矩(扭曲度)和四阶矩(鞘度)的值：

三阶矩：$\widehat{{\gamma }_3}=\underset{i=1}{\overset{D}{\Sigma }}{(h_i-\hat{m})}^3/(D-1){\hat{\sigma }}^3,$

四阶矩：$\widehat{{\gamma }_4}=3+\underset{i=1}{\overset{D}{\Sigma }}{(h_i-\hat{m})}^4/(D-1){\hat{\sigma }}^4,$

上述公式中，表示三阶矩，表示四阶矩，表示的均值，表示的方差；

步骤十三、将计算得到的三阶矩和四阶矩的值同阈值进行比较，三阶矩的阈值为T₃=τ₃θ₃，四阶矩的阈值为T₄＝τ₄θ₄；θ₃和θ₄和乱通常取值为10，τ₃和τ₄分别取值为和，这样$T_3=10\sqrt{6/(n_1\times n_2)}$，$T_4=20\sqrt{6/(n_1\times n_2)};$

步骤十四、当三阶矩和四阶矩同时大于各自的给定阈值，令计数变量加1(t＝t+1)，否则计数变量保持不变；

步骤十五、判断当前非线性主分量是否处理完毕：判断结果为是，则返回运行步骤十一；判断结果为否，则输出该计数变量，的值；

步骤十六、令波段序号s加1(s＝s+1)，判断s是否小于等于N：判断结果为是，则返回运行步骤十；判断结果为否，则表示所有非线性主分量都得到处理，运行下一步；

步骤十七、寻找存储的计数变量t最大值所对应的非线性主分量标号s，选择该标号对应的非线性主分量作为检测的输入；

步骤十八、利用如下的RX算子对选择的非线性主分量进行异常点判决：

$\mathrm{RX}\left(x\right)={(x-\hat{\mu })}^T{\left({\hat{\Gamma }}_x\right)}^1(x-\hat{\mu })\left(\begin{array}{cc}{H}_1& \\ >& \\ <& \eta \\ H_0& ,\end{array}\right)$

上述公式中，和分别是检测窗口内像素的均值和方差，x表示检测窗口的中心像素值，即待判决的像素，η是检测阈值，是根据所需要的检测概率自动确定的；

步骤十九、用RX算子在所选择的非线性主分量上逐像素点进行计算，最终得到高光谱图像异常点检测结果图像。

本发明能高效的对高光谱图像中异常点目标的特征进行提取和选择，而显著的降低了虚警率，并能有效的去除背景干扰，实现在有严重的背景干扰情况下正常对异常点进行检测。它还能应用在乳腺钙化点检测等医疗诊断的技术领域中。

附图说明

图1是原始高光谱图像的第六波段图，图2是异常点目标的真实分布图，图3是RX方法检测的结果图，图4是K-RX方法的结果图，图5是本发明方法检测的结果图，图6是RX方法、K-RX方法、发明方法检测结果的ROC比较图。

具体实施方式

具体实施方式一：本实施方式的检测方法步骤为：

步骤一、输入三维高光谱图像I(i，j，s)，其中，i＝1，2，…，P表示一个波段图像的行，j＝1，2，…，Q表示图像的列，s=1，2，…，N表示高光谱数据的波段数目，I(i，j，s)的大小为P×Q×N，按照波段图像方式，I(i，j，s)可以表示为，其中：

上述公式中，为第s波段扫(i，j)点对应的灰度值；

步骤二、求高光谱图像I(i，j，s)的灰度最大值$I_{\max }=\underset{i,j,s}{\max }\left(I(i,j,s)\right)$，即P×Q×N个点所对应的灰度的最大值，将高光谱图像I(i，j，s)的每个点灰度值除以该最大值I_max，得到归—化的高光谱图像I(i，j，s)＝I(i，j，s)/I_max；

步骤三、将归一化的高光谱图像I(i，j，s)转换成二维矩阵，即，对于某一个波段图像I_s，s＝1，2，3，…，N逐个波段图像按照由左至右、由上至下的顺序，将每一行像素值顺序放到一个行向量中为：

写成列向量形式，为：

这样得到的数据按列方向放置于二维矩阵中，得到大小为T×N的数据，其中T=P×Q，$y_s=I_{s\prime }^T$，s=l，2，3，…，N；可以看作是由N个列向量y_s(s＝1，2，3，…，N)构成的样本矩阵；

步骤四、输入样本矩阵尸，计算核函数矩阵$K,{\left(K\right)}_{i,j}=\left(\phi \left(y_i\right)\phi \left(y_i\right)\right)=k\left(y_i{y}_j\right)$表示核函数矩阵K的每个元素，k(y_i，y_j)代表输入样本y_i，y_J(i，j＝1，2，…，l)之间的的核函数值，具体核函数可以选择为高斯核函数(K)_ij＝k(y_i，y_j)＝exp{(y_i-y_j)²/2σ²}(σ为常数)；

步骤五、对核函数矩阵进行归一化，即K_c＝K-1_N·K-K·1_N+1_N·K·1_N，

上述公式中，1_N是一个N×N的矩阵，该矩阵的元素均为1/N，N是样本个数；

步骤六、对归一化核函数矩阵Kc进行特征值分解，即求解如下特征值方程：

${\lambda }_k{\beta }^k=K{\beta }^k,({\beta }^k={({\beta }_1^k,{\beta }_2^k,···,{\beta }_l^k)}^T);$

步骤七、根据特征值分解得到的特征向量构成投影变换矩阵，即：

${\beta }^k={({\beta }_1^k,{\beta }_2^k,···,{\beta }_l^k)}^T;$

步骤八、根据下式，将输入的二维形式的高光谱数据投影到特征向量${\beta }^k{=({\beta }_1^k,{\beta }_2^k,···,{\beta }_l^k)}^T$上，得到变换后的数据z；

$Z=(V^k·\phi \left(y\right))=\underset{i=1}{\overset{N}{\Sigma }}{\beta }_i^k\left(\phi \left(y_i\right)\phi \left(y\right)\right)=\underset{i=1}{\overset{N}{\Sigma }}{\beta }_i^{k}k(y_i,y),$

将z写成矩阵形式为，上述公式中，

步骤九、按照步骤三的逆过程，将投影后的高光谱数据转换成图像波段形式，得到变换后的与原始输入高光谱图像具有相同数目的非线性变换主分量，即将转换成：

，代表第s个变换后的非线性主分量，s＝1，2，…，N；步骤十、令t=0，t表示计数变量；

步骤十一、对第s个非线性主分量(s＝1，2，…，N)，构造同检测窗口大小一致的滑动窗w，滑动窗大小设定为n₁×n₂；设滑动窗w内对应的像素灰度值逐行排列放置到一个行向量中，得到

步骤十二、按照如下公式计算滑动窗内像素的三阶矩(，扭曲度)和四阶矩(鞘度)的值：

三阶矩：$\widehat{{\gamma }_3}=\underset{i=1}{\overset{D}{\Sigma }}{(h_i-\hat{m})}^3/(D-1){\hat{\sigma }}^3,$

四阶矩：$\widehat{{\gamma }_4}=3+\underset{i=1}{\overset{D}{\Sigma }}{(h_i-\hat{m})}^4/(D-1){\hat{\sigma }}^4,$

上述公式中，表示三阶矩，表示四阶矩，表示的均值，表示的方差；

步骤十三、将计算得到的三阶矩和四阶矩的值同阈值进行比较，三阶矩的阈值为T₃＝τ₃θ₃，四阶矩的阈值为T₄＝τ₄θ₄；θ₃和θ₄通常取值为10，τ₃和τ₄分别取值为和，这样$T_3=10\sqrt{6/(n_1\times n_2)}$，$T_4=20\sqrt{6/(n_1\times n_2)};$

步骤十四、当三阶矩和四阶矩同时大于各自的给定阈值，令计数变量加I(t＝t+1)，否则计数变量保持不变；

步骤十五、判断当前非线性主分量是否处理完毕：判断结果为是，则返回运行步骤十一；判断结果为否，则输出该计数变量，的值；

步骤十六、令波段序号s加1(s＝s+1)，判断s是否小于等于N：判断结果为是，则返回运行步骤十；判断结果为否，则表示所有非线性主分量都得到处理，运行下一步；

步骤十七、寻找存储的计数变量t最大值所对应的非线性主分量标号s，选择该标号对应的非线性主分量作为检测的输入：

步骤十八、利用如下的RX算子对选择的非线性主分量进行异常点判决：

$\mathrm{RX}\left(x\right)={(x-\hat{\mu })}^T{\left({\hat{\Gamma }}_x\right)}^1(x-\hat{\mu })\left(\begin{array}{cc}{H}_1& \\ >& \\ <& \eta \\ H_0& ,\end{array}\right)$

上述公式中，和分别是检测窗口内像素的均值和方差，x表示检测窗口的中心像素值，即待判决的像素，η是检测阈值，是根据所需要的检测概率自动确定的；

步骤十九、用RX算子在所选择的非线性主分量上逐像素点进行计算，最终得到高光谱图像异常点检测结果图像。

举例说明：将本发明方法同已有技术进行对比描述。其主要进行对比的是两种技术：一是国外学者提出的经典方法，即传统RX检测算子；二是，近年来提出的非线性RX检测算子，标记为K-RX，截至目前，该方法仍处于国际领先水平。结合图1、图2、图3、图4、图5说明。

为了进一步定量地说明本方法相对于已有技术的优势(结合图6说明)。图6中，传统检测算子标记为RX，国际新方法标记为K-RX，本发明提出的方法为基于选择性核主成份分析的高光谱图像异常点的检测方法，标记为SKPCA-RX(Selective Kernel Principal Component-RX)。

我们利用接收机操作特性曲线(即ROC，该指标是专门用于评价目标检测性能的)进行评价。ROC描述了检测概率Pd和虚警率Pf之间的变化关系，能够提供检测性能定量的描述。基于地面真实分布情况，可以获取真实异常点目标在图像场中的坐标，当检测出的疑似异常点属于异常点目标坐标范围，则判定为目标，否则判定为虚警。在整个图像场内，针对某方法的检测结果执行目标和虚警的辨识，这样可以定义最终的检测概率和虚警概率：

P_d＝N_d/N_T和P_f＝N_f /N_A

上式中，N_d是检测到的像素数目，N_T是图像场内真实目标像素数目，N_f是虚警像素个数，N_A是图像像素总数。图6中，横轴是虚警率，纵轴是检测概率，由图6可以看出，本发明方法具有良好的检测性能，大大超过了RX方法和K-RX方法。

异常点常携带重要的信息，且异常检测不需要已知目标的先验信息。因此，多维信号的异常检测，特别是基于高光谱图像的异常点检测，在许多研究领域有着广泛的应用，并起着非常重要的作用，例如遥感领域的环境变化检测，生物医学领域的癌症诊断、乳腺X光片微钙化点自动检测，军事领域的战场侦察，异常目标、伪装目标、隐蔽目标的检测等等。本发明方法将会推动高光谱图像异常点检测的广泛应用，进而将在国民经济领域和国防领域取得重大的经济效益和社会效益。