CN103473561A - 一种基于高斯分布的自适应高光谱分类方法 - Google Patents

Abstract

本发明公开了一种基于高斯分布的自适应高光谱分类方法。该方法根据光谱曲线的相似度，即目标与背景都符合高斯分布，将最小关联窗口与中心像元的相似度量通过高斯函数拟合获得数据分布特征，进而获得自适应阈值；然后利用像元最小关联窗口合并相邻相似像元为像块完成降元，实现高光谱图像分类。本发明方法不但实现了自适应，而且聚类精度高、边缘辨识度好及鲁棒性强。

Description

一种基于高斯分布的自适应高光谱分类方法

技术领域

本发明属于高光谱图像处理领域，特别是一种基于高斯分布的自适应高光谱分类方法。

背景技术

对于高光谱图像的目标分类，国内外的研究者一直致力于寻求高执行效率、低虚警概率、高分类概率的分类方法。高光谱分类最常见、最具代表性的是K-MEANS和ISODATA两种方法，这两种方法假设类内像元满足多元高斯分布，以高斯混合模型建模，采用欧拉距离作为基本测度手段。由于上述两种方法较简单、使用普遍,已经被集成在某些遥感图像处理软件中，而且这类算法只考虑了光谱信息忽略了图像空间特性，毗邻像元具有更高相似可能性——图像空间一致性。基于此特征的分类将更符合人类对图像的解析。

Reed和Yu两位学者利用一幅图像的局部灰度统计特性可以用高斯分布描述的性质，将一幅图像分成若干子块，认为如果背景灰度值在各个子块上近似符合高斯分布，则那些不符合高斯分布的像素点便可能是需要获得的目标像素点。以此为基础，他们提出了著名的RX算法，RX算法根据背景样本点的光谱向量均值和协方差矩阵，计算检测点光谱向量和背景样本的马氏距离，实现高光谱分类。但RX算法要求预先知道目标的形状和大小信息，而通常情况下这些信息可能是无法预先获得的。

发明内容

本发明提出一种基于高斯分布的自适应高光谱分类方法。该方法聚类精度高、边缘辨识度好及鲁棒性强，在对高光谱图像分类时，执行效率高，虚警概率低。

为了解决上述技术问题，本发明提供一种基于高斯分布的自适应高光谱分类方法，具体包括以下步骤：

步骤一：对高光谱图像二维灰度直方图的最佳一维投影进行高斯拟合得到最佳阈值；

步骤二：对高光谱图像矩阵依次完成降元操作，所述降元操作是指：对中心像元判断与其直接关联的8个像元的相似性,将具有相似性的像元划入同一像块，像元具有相似性的判断标准是该像元余弦夹角小于最佳阈值；

步骤三：利用一元回归线性模型在步骤二所获得的像块中选择用于合并的像块，并合并像块实现图像分类。

本发明与现有技术相比，其显著优点在于：

本发明方法采用了基于二维灰度直方图的最佳一维投影代替传统的一维直方图进行高斯拟合，这种一维投影既具有二维直方图的分割精度，又有一维直方图的分割速度，因而大大提高了图像分类的精度和速度；本发明方法充分利用统计直方图的峰谷特性，自动确定分割类数，并对各峰区的截断样本分布数据进行深入分析，利用部分样本点的最小偏度所对应的灰度值作为估计均值，计算各样本点分布数据与估计均值之比获得估计方差样本，取其截尾均值作为估计方差，消除了背景噪声的影响，避免了目标数据间的数据干扰，使拟合值更加精确，从而最终得到最佳分割阈值；通过降元操作和像块合并的方法，实现高光谱分类。

附图说明

图1是本发明方法流程图。

图2是本发明所选取的用于实验的两种场景的高光谱图像，其中图2(a)是场景一的高光谱图像，图2(b)是场景二的高光谱图像。

图3是图2中两种场景的高光谱图像经过本发明方法步骤一得到的高斯拟合曲线图，其中图3(a)是图2(a)的高斯拟合曲线图，图3(b)是图2(b)的高斯拟合曲线图。

图4是分别使用本发明方以及K-MEANS算法和ISODATA算法对图2场景的分类结果图，其中图4(a)是图2(a)使用本发明方法的分类结果；图4(b)是图2(a)使用K-MEANS算法的分类结果；图4(c)是图2(a)使用ISODATA算法的分类结果；图4(d)是图2(b)使用本发明方法的分类结果；图4(e)是图2(b)使用K-MEANS算法的分类结果；图4(f)是图2(b)使用ISODATA算法的分类结果。

具体实施方式

如图1所示，本发明方法包括以下步骤：

步骤一：对高光谱图像二维灰度直方图的最佳一维投影进行高斯拟合得到最佳阈值。

进一步，所述步骤一的一种具体实施方式为：

1.1在高光谱图像的二维灰度函数f(x,y)坐标中任取一像素点(x,y)，计算其局部灰度平均值g(x,y),计算方法如式（1）所示，

$g(x,y)=\frac{1}{n^2}\underset{i=-l}{\overset{l}{\mathrm{\Σ}}}\underset{j=-l}{\overset{l}{\mathrm{\Σ}}}f(x-i,y-j),n=2l+1---\left(1\right)$

式（1）中，f(x,y)是一幅M×N的数字图像表示成的二维灰度函数，M×N为图像大小，M为行，N为列；x为横坐标，y为纵坐标；(-l,l)为不包括(x,y)点的领域中像素坐标的集合，即去心领域；n表示集合(-l,l)内像素的总数；

1.2计算图像f(x,y)的邻域均值灰度直方图H(k),计算方法如式（2）所示，

$H\left(k\right)=\underset{x=1}{\overset{M}{\mathrm{\Σ}}}\underset{y=1}{\overset{N}{\mathrm{\Σ}}}g(x,y)\mathrm{\δ}(k-g(x,y))---\left(2\right)$

式（2）中，k为灰度值，δ(k-g(x,y))为冲激函数；

1.3将邻域均值灰度直方图H(k)的包络看作一条曲线，通过计算曲线的极值点获得直方图H(k)的峰谷位置和个数；谷点对应于曲线的极小值点，极小值点需同时满足

峰点对应于曲线的极大值点，极大值点需同时满足

$\frac{{\&PartialD;}^{2}H\left(k\right)}{{\&PartialD;k}^2}<0;$

1.4计算邻域均值灰度直方图H(k)的一阶差分H'(k)，计算方法如式（3）所示，

H'(k)＝H(k+1)-H(k)(k＝0,1,...,L-1)  （3）

式（3）中，L是图像的灰度级；

1.5计算H'(k)的过零点，由正到负为峰点，由负到正为谷点,定义相邻谷点间的区域为一个峰区；

1.6计算峰值点i的均值估计μ_i，具体方法为：计算峰区中以各灰度级为中心的样本偏度，取其绝对值最小的偏度所对应的灰度级作为该峰区的均值估计μ_i，从而保证其分布的对称性最佳；

1.7计算峰值点i的标准差估计σ_i和正态拟合函数p_i(μ_i)，具体方法为：

设拟合分布的概率密度函数为

其中A_i为常量系数；取x=μ_i,峰值点i的正态拟合函数

经进一步推导得出

其中k为灰度值；再进一步推导解出标准差估计

由于噪声的干扰，由估计均值两侧的样本可以得到多个σ_i值，但这些σ_i满足聚集态，可以通过聚类分析得到σ_i的最优拟合值；

1.8计算最佳阈值T_i，具体方法为：

设相邻两峰区正态拟合函数为p_i(x)和p_i+1(x)，最佳阈值T_i应满足p_i(T_i)＝p_i+1(T_i)，经进一步推导得出式（4），

$\left\{\begin{array}{c}{\mathrm{AT}}_i^2+{\mathrm{BT}}_i+C=0\\ A={\mathrm{\&sigma;}}_i^2-{\mathrm{\&sigma;}}_{i+1}^2\\ B=2({\mathrm{\&mu;}}_i{\mathrm{\&sigma;}}_{i+1}^2-{\mathrm{\&mu;}}_{i+1}{\mathrm{\&sigma;}}_i^2)\\ C={\mathrm{\&sigma;}}_i^2{\mathrm{\&mu;}}_{i+1}^2-{\mathrm{\&sigma;}}_i^2{\mathrm{\&mu;}}_i^2+{2\mathrm{\&sigma;}}_i^2{2\mathrm{\&sigma;}}_{i+1}^2\ln (p_{i+1}\left({\mathrm{\&mu;}}_{i+1}\right)/p_i\left({\mathrm{\&mu;}}_i\right))\end{array}---\left(4\right)\right.$

式（4）中，μ_i+1是样本峰值点i+1点的均值估计，σ_i+1是样本峰值点i+1点的标准差估计，p_i+1(μ_i+1)是样本峰值点i+1的正态拟合函数，

将此前计算得到的μ_i，μ_i+1，σ_i，σ_i+1以及p(μ_i)，p(μ_i+1)代入式（4）即可进一步解出最佳阈值T_i。

步骤二：对高光谱图像矩阵依次完成降元操作，所述降元操作是指：对中心像元判断与其直接关联的8个像元的相似性,将具有相似性的像元划入同一像块，像元具有相似性的判断标准是该像元余弦夹角小于最佳阈值。

进一步，所述步骤二的一种具体实施方式可以为：

2.1将高光谱图像生成光谱向量Z矩阵，Z矩阵大小与高光谱图像大小相同，Z矩阵中每个坐标点对应一个w维向量，w维向量由高光谱图像中每个坐标点对应的灰度值组成；

2.2利用像元最小关联窗口模型遍历Z矩阵，计算最小关联窗口的中心像元与其直接关联的8个像元的余弦夹角SAM_(a+i,b+j)，

所述最小关联窗口模型如式（5）所示，

Z_mod＝Z(x_a+i,y_b+j){-1≤i≤1,-1≤j≤1}  (5)

式（5）中，(x_a,y_b)为中心像元的坐标；当变量i,j都不为0时，(x_a+i,y_b+j)表示与最小关联窗口的中心像元(x_a,y_b)关联的8个像元的坐标；

所述余弦夹角SAM_(a+i,b+j)的计算方法如式（6）所示，

${\mathrm{SAM}}_{(a+i,b+j)}={\cos }^{-1}\frac{Z(x_{a+i},y_{b+j})\&CenterDot;Z(x_a,y_b)}{\left|Z\right(x_{a+i},y_{b+j}\left)\right|\left|Z(x_a,y_b)\right|},\{-1\&le;i\&le;1,-1\&le;j\&le;1\}---\left(6\right)$

式（6）中，Z(x_a+i,y_b+j)表示关联像元光谱向量，Ζ(x_a,y_b)表示中心像元光谱向量；

2.3比较余弦夹角SAM_(a+i,b+j)与最佳阈值T_i的大小，将余弦夹角SAM_(a+i,b+j)小于最佳阈值T_i的像元划入同一像块。

降元具有继承性，如果最小关联窗口模型形成的像块中有前面已经划定的像元，则此像块并入前面已划定的像块，依此类推，完成整个降元过程。

步骤三：利用一元回归线性模型在步骤二所获得的像块中选择用于合并的像块，并合并像块实现图像分类。

进一步，建立所述一元回归线性模型的具体方法为：

3.1在降元操作中,设定阈值T_i为零时，像块内的所有像元光谱向量严格满足线性关系，用式（7）表示，

${\stackrel{\&RightArrow;}{V}}_i=k_i\&CenterDot;\stackrel{\&RightArrow;}{V}---\left(7\right)$

式（7）中，

是像块第i个像元的光谱向量，i取正整数；k_i是标量系数；

是此像块的基准向量；

3.2设定阈值T_i不为零时，像块内像元光谱曲线形状不完全相同，光谱向量不满足严格的线性关系,存在一定误差，此时像元光谱向量如式（8）所示，

${\stackrel{\&RightArrow;}{V}}_i=k_i\&CenterDot;\stackrel{\&RightArrow;}{V}+{\stackrel{\&RightArrow;}{\mathrm{\&epsiv;}}}_i---\left(8\right)$

式（8）中， ${\stackrel{\&RightArrow;}{V}}_i=\left[\begin{array}{c}{v}_{\mathrm{il}}\\ .\\ .\\ .\\ v_{\mathrm{ip}}\end{array}\right],\stackrel{\&RightArrow;}{V}=\left[\begin{array}{c}{v}_1\\ .\\ .\\ .\\ v_p\end{array}\right],{\stackrel{\&RightArrow;}{\mathrm{\&epsiv;}}}_i=\left[\begin{array}{c}{\mathrm{\&epsiv;}}_{\mathrm{il}}\\ .\\ .\\ .\\ {\mathrm{\&epsiv;}}_{\mathrm{ip}}\end{array}\right],$

式（8）与式（7）相比，多了向量余项

本发明称之为余向量，余向量描述的是像块中各个像元光谱向量与基准向量的偏差，余向量直接影响像块数据的线性显著性，若余向量不大时，即像块中各个像元光谱向量与基准向量的偏差不大时，像元间是一种显著的线性关系，此时式（9）成立，

$\frac{{\stackrel{\&RightArrow;}{V}}_i}{{\stackrel{\&RightArrow;}{V}}_j}=\frac{k_i\stackrel{\&RightArrow;}{V}+{\stackrel{\&RightArrow;}{\mathrm{\&epsiv;}}}_i}{k_j\stackrel{\&RightArrow;}{V}+{\stackrel{\&RightArrow;}{\mathrm{\&epsiv;}}}_j}\&ap;\frac{k_i}{k_j}---\left(9\right)$

式（9）中，

是像块第j个像元的光谱向量，k_j是标量系数；

3.3根据式（8），可以推导出像元光谱向量

内部两个光谱向量v_ir和v_is满足式（10）所示关系，

$\frac{v_{\mathrm{ir}}}{v_{\mathrm{is}}}=\frac{k_i\&CenterDot;v_r+{\stackrel{\&RightArrow;}{\mathrm{\&epsiv;}}}_{\mathrm{ir}}}{k_i\&CenterDot;v_s+{\stackrel{\&RightArrow;}{\mathrm{\&epsiv;}}}_{\mathrm{is}}}---\left(10\right)$

式（10）中，r,s表示任意两个谱段r和s的值；

在余向量

不大的情况下，式（10）可以进一步表达为式（11），

$\frac{v_{\mathrm{ir}}}{v_{\mathrm{is}}}=\frac{k_i\&CenterDot;v_r+{\stackrel{\&RightArrow;}{\mathrm{\&epsiv;}}}_{\mathrm{ir}}}{k_i\&CenterDot;v_s+{\stackrel{\&RightArrow;}{\mathrm{\&epsiv;}}}_{\mathrm{is}}}\&ap;\frac{v_r}{v_s}=P(r,s)---\left(11\right)$

式（11）中，P(r,s)是一个常数，表示v_ir和v_is满足线性关系，可以看出各光谱向量内部谱段间也满足显著线性关系，且与像元光谱向量

内部两个谱段基准向量v_r和v_s有关；

3.4根据式（9）中得出的各像元光谱向量的线性关系和式（11）中得出的像元内部谱段光谱向量的线性关系，可以建立如式（12）所示的一元回归线性模型来同时描述这两种线性关系，

$v_{i(r+1)}={\mathrm{\&beta;}}_0+{\mathrm{\&beta;}}_1{v}_{\mathrm{ir}}+\stackrel{\&RightArrow;}{\mathrm{\&epsiv;}}---\left(12\right)$

式（12）中，r和r+1是相邻的两个谱段；v_i(r+1)是第r+1谱段的光谱向量,v_ir是第r谱段的光谱向量；i表示像块中的第i个像元，β₀为常数项，β₁为回归系数；

将式(12)给定显著性水平α＝0.05，使用“F检验”判断该像块的线性显著性，在一个基准向量估计过程中，当不满足线性显著性的谱段超过5%时，则丢弃此像块；当满足线性显著性的谱段超过5%时，则选择此像块用于合并。

本发明方法的效果可以通过以下实验结果进一步说明：

从图4看出，本发明方法与K-MEANS算法和ISODATA算法相比，本发明方法对边缘更敏感，结果中各类的边缘和轮廓更清晰、准确，因而更易于辨识；本发明方法不仅分类准确率更高，而且分类准确率大小跨度最小，因而更具鲁棒性、更稳定；再者有更低的分类噪声，错误分类像元少，提高了分类准确性。

Claims (4)

1.一种基于高斯分布的自适应高光谱分类方法，其特征在于，包括以下步骤：

步骤一：对高光谱图像二维灰度直方图的最佳一维投影进行高斯拟合得到最佳阈值；

步骤二：对高光谱图像矩阵依次完成降元操作，所述降元操作是指：对中心像元判断与其直接关联的8个像元的相似性,将具有相似性的像元划入同一像块，像元具有相似性的判断标准是该像元余弦夹角小于最佳阈值；

步骤三：利用一元回归线性模型在步骤二所获得的像块中选择用于合并的像块，并合并像块实现图像分类。

2.如权利要求1所述基于高斯分布的自适应高光谱分类方法，其特征在于，所述步骤一具体为：

2.1在高光谱图像的二维灰度函数f(x,y)坐标中任取一像素点(x,y)，计算其局部灰度平均值g(x,y),计算方法如式（1）所示，

$g(x,y)=\frac{1}{n^2}\underset{i=-l}{\overset{l}{\mathrm{\&Sigma;}}}\underset{j=-l}{\overset{l}{\mathrm{\&Sigma;}}}f(x-i,y-j),n=2l+1---\left(1\right)$

式（1）中，f(x,y)是一幅M×N的数字图像表示成的二维灰度函数，M×N为图像大小，M为行，N为列；x为横坐标，y为纵坐标；(-l,l)为不包括(x,y)点的领域中像素坐标的集合，即去心领域；n表示集合(-l,l)内像素的总数；

2.2计算图像f(x,y)的邻域均值灰度直方图H(k),计算方法如式（2）所示，

$H\left(k\right)=\underset{x=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{y=1}{\overset{N}{\mathrm{\&Sigma;}}}g(x,y)\mathrm{\&delta;}(k-g(x,y))---\left(2\right)$

式（2）中，k为灰度值，δ(k-g(x,y))为冲激函数；

2.3将邻域均值灰度直方图H(k)的包络看作一条曲线，通过计算曲线的极值点获得直方图H(k)的峰谷位置和个数；谷点对应于曲线的极小值点，极小值点需同时满足

峰点对应于曲线的极大值点，极大值点需同时满足

$\frac{{\&PartialD;}^{2}H\left(k\right)}{{\&PartialD;k}^2}<0;$

2.4计算邻域均值灰度直方图H(k)的一阶差分H'(k)，计算方法如式（3）所示，

H'(k)＝H(k+1)-H(k)(k＝0,1,...,L-1)   （3）

式（3）中，L是图像的灰度级；

2.5计算H'(k)的过零点，由正到负为峰点，由负到正为谷点,定义相邻谷点间的区域为一个峰区；

2.6计算峰值点i的均值估计μ_i，具体方法为：计算峰区中以各灰度级为中心的样本偏度，取其绝对值最小的偏度所对应的灰度级作为该峰区的均值估计μ_i；

2.7计算峰值点i的标准差估计σ_i和正态拟合函数p_i(μ_i)，具体方法为：

设拟合分布的概率密度函数为

其中A_i为常量系数；取x=μ_i,峰值点i的正态拟合函数经进一步推导得出其中k为灰度值；再进一步推导解出标准差估计

2.8计算最佳阈值T_i，具体方法为：

设相邻两峰区正态拟合函数为p_i(x)和p_i+1(x)，最佳阈值T_i应满足p_i(T_i)＝p_i+1(T_i)，经进一步推导得出式（4），再进一步可以由式（4）解出最佳阈值T_i，

$\left\{\begin{array}{c}{\mathrm{AT}}_i^2+{\mathrm{BT}}_i+C=0\\ A={\mathrm{\&sigma;}}_i^2-{\mathrm{\&sigma;}}_{i+1}^2\\ B=2({\mathrm{\&mu;}}_i{\mathrm{\&sigma;}}_{i+1}^2-{\mathrm{\&mu;}}_{i+1}{\mathrm{\&sigma;}}_i^2)\\ C={\mathrm{\&sigma;}}_i^2{\mathrm{\&mu;}}_{i+1}^2-{\mathrm{\&sigma;}}_i^2{\mathrm{\&mu;}}_i^2+{2\mathrm{\&sigma;}}_i^2{2\mathrm{\&sigma;}}_{i+1}^2\ln (p_{i+1}\left({\mathrm{\&mu;}}_{i+1}\right)/p_i\left({\mathrm{\&mu;}}_i\right))\end{array}---\left(4\right)\right.$

式（4）中，μ_i+1是样本峰值点i+1点的均值估计，σ_i+1是样本峰值点i+1点的标准差估计，p_i+1(μ_i+1)是样本峰值点i+1的正态拟合函数。

3.如权利要求1所述基于高斯分布的自适应高光谱分类方法，其特征在于，所述步骤二具体为：

3.1将高光谱图像生成光谱向量Z矩阵，Z矩阵大小与高光谱图像大小相同，Z矩阵中每个坐标点对应一个w维向量，w维向量由高光谱图像中每个坐标点对应的灰度值组成；

3.2利用像元最小关联窗口模型遍历Z矩阵，计算最小关联窗口的中心像元与其直接关联的8个像元的余弦夹角SAM_(a+i,b+j)，

所述最小关联窗口模型如式（5）所示，

Z_mod＝Z(x_a+i,y_b+j){-1≤i≤1,-1≤j≤1}    (5)

式（5）中，(x_a,y_b)为中心像元的坐标；当变量i,j都不为0时，(x_a+i,y_b+j)表示与最小关联窗口的中心像元(x_a,y_b)关联的8个像元的坐标；

所述余弦夹角SAM_(a+i,b+j)的计算方法如式（6）所示，

${\mathrm{SAM}}_{(a+i,b+j)}={\cos }^{-1}\frac{Z(x_{a+i},y_{b+j})\&CenterDot;Z(x_a,y_b)}{\left|Z\right(x_{a+i},y_{b+j}\left)\right|\left|Z(x_a,y_b)\right|},\{-1\&le;i\&le;1,-1\&le;j\&le;1\}---\left(6\right)$

式（6）中，Z(x_a+i,y_b+j)表示关联像元光谱向量，Ζ(x_a,y_b)表示中心像元光谱向量；

3.3比较余弦夹角SAM_(a+i,b+j)与最佳阈值T_i的大小，将余弦夹角SAM_(a+i,b+j)小于最佳阈值T_i的像元划入同一像块。

4.如权利要求1所述基于高斯分布的自适应高光谱分类方法，其特征在于，步骤三所述一元回归线性模型如式（7）所示，

$v_{i(r+1)}={\mathrm{\&beta;}}_0+{\mathrm{\&beta;}}_1{v}_{\mathrm{ir}}\stackrel{\&RightArrow;}{\mathrm{\&epsiv;}}---\left(7\right)$

式（7）中，r和r+1是相邻的两个谱段；v_i(r+1)是第r+1谱段的光谱向量,v_ir是第r谱段的光谱向量；i表示像块中的第i个像元，β₀为常数项，β₁为回归系数；

所述选择合并像块的方法具体为:将式(7)给定显著性水平α＝0.05，使用“F检验”判断该像块的线性显著性，在一个基准向量估计过程中，当不满足线性显著性的谱段超过5%时，则丢弃此像块；当满足线性显著性的谱段超过5%时，则选择此像块用于合并。

Images