CN108197650A - 局部相似性保持的高光谱图像极限学习机聚类方法 - Google Patents

Abstract

本发明公开了一种局部相似性保持的高光谱图像极限学习机聚类方法，其步骤为：组织高光谱像元矩阵；计算隐层神经元的线性随机响应；计算隐层神经元的非线性激活值；隐层特征数据三维重组；空间引导滤波；滤波后的隐层特征数据二维重组；构造局部相似性保持正则项及优化模型；计算局部相似性保持投影特征，并进行K‑means聚类得到最终的聚类标签。本发明在传统极限学习机的基础上，通过引导滤波综合局部邻域的高光谱图像空间信息，并充分利用高光谱的光谱局部相似性，通过模型优化计算具有局部保持性的投影，提取空谱联合信息，提高了聚类精度，降低了计算复杂度，可广泛应用于国土资源、矿产调查和精准农业领域的高光谱无监督分类。

Description

局部相似性保持的高光谱图像极限学习机聚类方法

技术领域

本发明属于遥感图像处理技术领域，具体涉及一种局部相似性保持的高光谱图像极限学习机聚类方法。

背景技术

高光谱图像由于其光谱相关性及丰富的空间信息而被广泛应用于军事监测、精细农业和矿物监测等领域，其中，高光谱图像聚类是最重要的研究内容之一。高光谱图像聚类的基本原理是在图像聚类的基础上，结合高光谱图像的光谱特性对目标图像进行无监督分类。其理论依据是相同的像元间具有相同或相似的光谱特征，反之，不同的像元对应的光谱特征不同。

目前，已经有许多针对高光谱图像的聚类算法被提出，其中最为有效的包括K均值聚类、谱聚类、拉普拉斯特征映射、稀疏子空间聚类和低秩子空间聚类，但由于高光谱图像波段数多及光谱分辨率高等特性给上述聚类算法带来计算复杂度的难题。极限学习机凭借随机分配输入权重和求解最小二乘，即可有效获得网络输出权重的解析解的快速学习过程被广泛应用到了遥感图像识别中。2014年，Huang等人在流形正则化的基础上提出半监督和无监督的极限学习机方法[Huang G,Song S,Gupta J N,et al.Semi-supervised andunsupervised extreme learning machines[J].IEEE Transactions on Cybernetics,2014,44(12):2405]，取得了良好的聚类效果。然而，上述方法仅仅利用了高光谱的逐像素的光谱信息，没有有效联合空间-光谱信息，聚类精度较低，而且当数据存在噪声时算法性能下降。

发明内容

本发明的目的在于针对地物覆盖分析、精准农业和矿物调查等领域的高光谱无监督分类问题，提出一种无需人工标记分类样本、对噪声鲁棒的局部相似性保持的高光谱图像极限学习机聚类方法。

实现本发明目的的技术解决方案为：一种局部相似性保持的高光谱图像极限学习机聚类方法，其特征在于，该方法包括如下步骤：

步骤1，组织高光谱像元矩阵；

步骤2，计算隐层神经元的线性随机响应；

步骤3，计算隐层神经元的非线性激活值；

步骤4，隐层特征数据三维重组；

步骤5，空间引导滤波；

步骤6，滤波后的隐层特征数据二维重组；

步骤7，构造局部相似性保持正则项及优化模型；

步骤8，计算局部相似性保持投影特征，并进行K-means聚类得到最终的聚类标签。

进一步的，步骤1具体为：

输入一幅高光谱图像X₀∈R^D×W×H，其中，D表示高光谱图像的波段数，W和H分别表示图像空间维的宽度和高度；

将原始高光谱数据X₀按行逐像素排列形成矩阵X＝[x₁,x₂,...,x_N]∈R^D×N作为模型的输入，其中，N＝W×H表示高光谱像元的个数，x_i∈R^D表示一个高光谱像元。

进一步的，步骤2具体为：

步骤2-1，计算第j个隐层神经元对第i个高光谱像元的线性随机响应：

其中，w_j＝[w_j1,w_j2,...,w_jD]∈R^D是随机生成的权重向量，b_j＝[b_j1,b_j2,...,b_jD]∈R^D是随机偏置，D表示高光谱像元的波段数，M表示隐层神经元的个数，1≤i≤N；

步骤2-2，构造线性随机投影向量，即对于高光谱像元x_i，将M个隐层神经元拼装，构成对应的线性随机投影行向量o_i＝[o_i1,o_i2,...,o_iM]。

进一步的，步骤3具体为：

(1)对隐层线性随机投影行向量进行非线性映射，可得

h_i＝[g(o_i1)…g(o_iM)]^T

其中，h_i∈R^M是对应于x_i的隐层输出，g(·)是激活函数；

(2)构造整体输入数据X对应的M×N大小的隐层输出矩阵：

H＝[h₁,h₂,...,h_N]∈R^M×N

进一步的，步骤4中隐层特征数据三维重组具体为：

将二维的隐层输出矩阵H∈R^M×N按列取矩阵的元素转化成立方体数据T∈R^M×W×H，其中，W和H分别表示图像空间维的宽度和高度，M表示隐层输出特征维度；

记p＝(p₁,p₂)表示空间上的任意点，其中，1≤p₁≤W，1≤p₂≤H；T(p)∈R^M表示立方体数据在该位置p的隐层输出特征向量。

进一步的，步骤5中空间引导滤波具体为：

(1)构造与位置p的像元相邻的像元集合其中邻域像元

(2)对任意位置p处，计算引导滤波权重，其计算式为：

w_p,q＝w_p,q-1D(q,q-1)R(p,q)

其中，D(q,q-1)表示像元q和像元q-1之间的相似性权重，R(p,q)表示像元p和像元q之前的相似性权重，具体定义如下：

其中，σ_d是D(q,q-1)的高斯核方差，σ_r是R(p,q)的高斯核方差；

(3)对立方体数据进行空间引导滤波，其计算公式为：

其中，Z(p)是归一化因子，T为原始的数据，为空间引导滤波后的数据。

进一步的，步骤6中滤波后的隐层特征数据二维重组具体为：

将空间引导滤波后所得的逐像元重排为二维的特征矩阵其中，将作为新的融合空间信息的隐层输出矩阵。

进一步的，步骤7具体为：

(1)构造光谱相似性连接图：

结合基于图的局部相似性保持方法，原始光谱数据根据任意两个光谱像元的相关性可以建模为包含N个顶点的图，每条边对应两个顶点的相似性权重，其中，权重的定义如下：

其中，s_i,j表示像元x_i和x_j的相似度，σ是带宽参数，σ≥0，表示x_i的n个最相似近邻的集合，n＝5；

(2)建立光谱相似性保持正则项：

上式等价于：

其中，Tr(·)表示矩阵的迹，L＝D-S是图拉普拉斯矩阵，D是对角矩阵，d_ii＝∑_js_ij，A是要求解的局部相似性保持投影矩阵。

(3)建立局部光谱相似性保持的最小代价函数模型：

其中，λ是正则化参数，I_C是C×C的单位矩阵；

(4)模型求解：

利用拉格朗日乘子法，上述模型的解可表述为：

其中，v₁,v₂,...,v_C+1是求解的前C+1个最小特征值γ₁,γ₂,...,γ_C+1对应的特征向量，γ₁≤γ₂≤...≤γ_C+1，是标准化的特征向量，其中，C表示输出特征的维度，I_M表示M×M的单位矩阵；

进一步的，步骤8具体为：计算局部相似性保持投影特征，并进行K-means聚类得到最终的聚类标签；

最终的输出特征为将的每一行看作一个像元的局部相似性保持投影特征向量，用K均值的方法将N个像元聚成k类，得到最终类别标签向量y。

本发明与现有技术相比，其显著优点为：通过引导滤波综合局部邻域的高光谱图像空间信息，并充分利用高光谱的光谱局部相似性，通过模型优化计算具有局部保持性的投影，提取空谱联合信息，与传统的子空间聚类方法相比，提高了聚类精度，增强了噪声的鲁棒性，降低了计算复杂度；本发明可广泛应用于国土资源、矿产调查和精准农业领域的高光谱无监督分类。

下面结合附图对本发明作进一步详细描述。

附图说明

图1是本发明局部相似性保持的高光谱图像极限学习机聚类方法流程图。

图2是引导滤波权重的定义图。

图3是引导滤波权重的计算流程图。

图4(a)为Salinas-A数据集的真实地物分布图。

图4(b)为Salinas-A数据集采用K-means方法的聚类效果图。

图4(c)为Salinas-A数据集采用Fuzzy C均值聚类方法的聚类效果图。

图4(d)为Salinas-A数据集采用谱聚类方法的聚类效果图。

图4(e)为Salinas-A数据集采用拉普拉斯特征映射聚类方法的聚类效果图。

图4(f)为Salinas-A数据集采用堆栈式自编码机聚类方法的聚类效果图。

图4(g)为Salinas-A数据集采用稀疏子空间聚类方法的聚类效果图。

图4(h)为Salinas-A数据集采用低秩子空间聚类方法的聚类效果图。

图4(i)为Salinas-A数据集采用无监督极限学习机聚类方法的聚类效果图。

图4(j)为Salinas-A数据集采用本发明所提的局部相似性保持的极限学习机聚类方法的聚类效果图。

具体实施方式

结合图1，一种局部相似性保持的高光谱图像极限学习机聚类方法，包括如下步骤：

步骤1：重组高光谱像元矩阵：输入一幅高光谱图像X₀∈R^D×W×H，其中，D表示高光谱图像的波段数，W和H分别表示图像空间维的宽度和高度；将原始高光谱数据X₀逐像素排列形成矩阵X＝[x₁,x₂,...,x_N]∈R^D×N作为模型的输入，其中，N＝W×H表示高光谱像元的个数，x_i∈R^D表示一个高光谱像元。

步骤2：计算隐层神经元的线性随机响应：对任意高光谱像元x_i，1≤i≤N，计算其线性随机投影行向量o_i∈R^M。

步骤3：计算隐层神经元的非线性激活值：对隐层线性输出数据o_i进行非线性映射，可得：h_i＝[g(o_i1)…g(o_iM)]^T，其中，h_i∈R^M是对应于x_i的隐层输出，g(·)是激活函数。整体输入数据X对应的隐层输出矩阵：H＝[h₁,h₂,...,h_N]∈R^M×N。

步骤4：隐层特征数据三维重组：将二维的隐层输出矩阵H∈R^M×N按列取矩阵的元素转化成立方体数据T∈R^M×W×H。记p＝(p₁,p₂)表示空间上的任意点，其中，1≤p₁≤W，1≤p₂≤H。T(p)∈R^M表示立方体数据在该位置p的隐层输出特征向量。

步骤5：空间引导滤波：对立方体数据T做空间双边滤波处理，得到新的立方体特征数据

步骤6：滤波后的隐层特征数据二维重组：将空间引导滤波后所得的逐像元重排为二维的特征矩阵其中，将作为新的融合空间信息的隐层输出矩阵。

步骤7：构造局部相似性保持正则项及优化模型：结合基于图的局部相似性保持方法，建立光谱相似性保持正则项：通过求解局部光谱相似性保持的最小代价函数模型得到局部相似性保持投影A。

步骤8：计算局部相似性保持投影特征：计算局部相似性保持投影特征将的每一行看作一个像元，用K均值的方法将N个像元聚成k类，得到最终类别标签向量y。

下面结合实施例具体说明本发明的内容。

实施例

结合图1，下面详细说明本发明的实施过程，步骤如下：

步骤1：重组高光谱像元矩阵：输入一幅高光谱图像X₀∈R^D×W×H，以图4(a)所示图像波段数D＝204，图像宽度W＝86，图像高度H＝83的Selina_A数据集为实验用例。

将原始高光谱数据X₀逐像素排列形成矩阵X＝[x₁,x₂,...,x_N]∈R^D×N作为模型的输入，其中，N＝W×H表示高光谱像元的个数，x_i∈R^D表示一个高光谱像元。

步骤2：计算隐层神经元的线性随机响应：具体过程如下：

2.1计算第j个隐层神经元对第i个高光谱像元的线性随机响应：

其中，1≤i≤N，w_j＝[w_j1,w_j2,...,w_jD]∈R^D是随机生成的权重向量，b_j＝[b_j1,b_j2,...,b_jD]∈R^D是随机偏置，D表示高光谱像元的波段数，M表示隐层神经元的个数，本发明中M＝2000。

2.2构造线性随机投影向量，即，对于高光谱像元x_i，将M个隐层神经元拼装，构成对应的线性随机投影行向量o_i＝[o_i1,o_i2,...,o_iM]。

步骤3：计算隐层神经元的非线性激活值：对隐层线性输出数据o_i进行非线性映射，可得：h_i＝[g(o_i1)…g(o_iM)]^T，其中，h_i∈R^M是对应于x_i的隐层输出，g(·)是激活函数。整体输入数据X对应的隐层输出矩阵：H＝[h₁,h₂,...,h_N]∈R^M×N。

步骤4：隐层特征数据三维重组：将二维的隐层输出矩阵H∈R^M×N按列取矩阵的元素转化成立方体数据T∈R^M×W×H，其中，本发明实验中W＝86和H＝83分别表示图像空间维的宽度和高度，M＝2000表示隐层输出特征维度。记p＝(p₁,p₂)表示空间上的任意点，其中，1≤p₁≤W，1≤p₂≤H。T(p)∈R^M表示立方体数据在该位置p的隐层输出特征向量。

步骤5：空间引导滤波：对立方体数据T做空间双边滤波处理，得到新的立方体特征数据引导滤波的具体过程如下：

5.1构造与位置p的像元相邻的像元集合其中邻域像元

5.2对任意位置p处，计算引导滤波权重，其计算式为：

w_p,q＝w_p,q-1D(q,q-1)R(p,q)

其中，D(q,q-1)表示像元q和像元q-1之间的相似性权重，R(p,q)表示像元p和像元q之前的相似性权重，具体定义如下：

其中，σ_d是D(q,q-1)的高斯核方差，σ_r是R(p,q)的高斯核方差。引导滤波权重的定义图解如图2所示。

5.3对立方体数据进行空间引导滤波，其计算公式为：

其中，Z(p)是归一化因子。T为原始的数据，为空间引导滤波后的数据。引导滤波的计算流程如图3所示。

步骤6：滤波后的隐层特征数据二维重组：将空间引导滤波后所得的逐像元重排为二维的特征矩阵其中，将作为新的融合空间信息的隐层输出矩阵。

步骤7：构造局部相似性保持正则项及优化模型：结合基于图的局部相似性保持方法，建立光谱相似性保持正则项，并通过求解局部光谱相似性保持的最小代价函数模型得到局部相似性保持投影A，具体过程如下:

7.1构造光谱相似性连接图：

结合基于图的局部相似性保持方法，原始光谱数据根据任意两个光谱像元的相关性可以建模为包含N个顶点的图，每条边对应两个顶点的相似性权重，其中，权重的定义如下：

其中，s_i,j表示像元x_i和x_j的相似度，σ是带宽参数(σ≥0)，表示x_i的n个最相似近邻的集合，n＝5。

7.2建立光谱相似性保持正则项：

上式等价于：

其中，Tr(·)表示矩阵的迹，L＝D-S是图拉普拉斯矩阵，D是对角矩阵，d_ii＝∑_js_ij，A是要求解的输出权重矩阵。

7.3建立局部光谱相似性保持的最小代价函数模型：

其中，λ是正则化参数，I_C是C×C的单位矩阵；

7.4模型求解：

利用拉格朗日乘子法，上述模型的解可表述为：

其中，v₁,v₂,...,v_C+1是求解的前C+1个最小特征值γ₁,γ₂,...,γ_C+1(γ₁≤γ₂≤...≤γ_C+1)对应的特征向量，i＝2,...,C+1是标准化的特征向量，其中参数C表示输出特征的维度，在本发明的实验中C＝6，I_M表示M×M的单位矩阵；

步骤8：计算局部相似性保持投影特征：计算局部相似性保持投影特征将的每一行看作一个像元，用K均值的方法将N个像元聚成k类，得到最终类别标签向量y。

本发明的效果可通过以下仿真实验进一步说明：

(1)仿真条件

仿真实验采用三组真实高光谱数据：Salinas-A数据集、Pavia University数据集和Indian Pines数据集。Salinas-A数据集是由加利福尼亚州Salinas Valley的AVIRIS传感器收集的Salinas图像的子集，去除20个吸水带(108-112,154-167,224)，共包含204个波段，图像的大小为86×83。Pavia University数据集是由帕维亚的ROSIS传感器采集，共包含115个波段，图像大小为610×340，在去除噪声波段之后，选择剩下的103个波段作为研究对象。考虑到计算复杂度问题，本发明切割了一个大小为200×100的子图。Indian Pines数据集为机载可见红外成像光谱仪(AVIRIS)在美国印第安纳州Indian Pines实验区采集的高光谱遥感图像。该图像共包含220个波段，空间分辨率为20m，图像大小为145×145。去除20个水汽吸收和低信噪比波段后(波段号为104-108,150-163,220)，选择剩下的200个波段作为研究对象。该地区共包含16种已知地物共10366个样本。有关数据集的详细信息如表1所示。仿真实验均在Windows 7操作系统下采用MATLAB R2012a完成。

表1三组数据集信息表

本发明采用的评价指标是聚类精度的评价方法(ACC，Calculation method ofclustering accuracy)及标准化互信息(NMI，Normalized Mutual Information)。

(2)仿真内容

本发明采用真实高光谱数据集检验算法的聚类性能。为测试本发明算法的性能，将提出的局部相似性保持的高光谱图像极限学习机聚类方法(S-USELM)与目前国际上流行的聚类算法对比。对比方法包括：K-means，Fuzzy C-均值(FCM)，堆栈式自编码算法(SAE)，谱聚类(SC)，拉普拉斯特征映射(LE)，稀疏子空间聚类(SSC),低秩子空间聚类(LRSC)，无监督极限学习机算法(USELM)，本发明所提出的融合空间信息的无监督极限学习机算法(S-USELM)。

(3)仿真实验结果分析

表2为三组高光谱数据在不同聚类算法下聚类精度和标准化互信息聚类精度的对比结果。

表2三组数据集的聚类精度及互信息精度

由表2可以看出，在Selina_A数据集中，USELM凭借局部相似性保持的特性，与传统的子空间聚类方法相比，聚类精度提高了5.31％，标准化互信息聚类精度提高了10.32％，本发明所提的融合了空间信息的局部相似性保持的极限学习机聚类方法在USELM的基础上提高了1.4％，标准化互信息聚类精度提高了1.22％，聚类效果最显著。

表3训练时间对比

由表3可以看出，USELM的计算时间远远小于SSC和LRSC算法的计算时间，本发明所提出的S-USELM算法虽然加入了空间滤波后计算时间增加，但是与传统的聚类算法相比仍有很强的竞争性。图4(a)～图4(j)为Salinas-A数据集在不同聚类算法下的聚类效果图，由图4(j)可以看出本发明所提出的融合了空间信息的局部相似性保持的极限学习机聚类方法与图4(i)所示的USELM方法相比提出了聚类噪声，聚类效果最显著。

Claims (9)

1.一种局部相似性保持的高光谱图像极限学习机聚类方法，其特征在于，该方法包括如下步骤：

步骤1，组织高光谱像元矩阵；

步骤2，计算隐层神经元的线性随机响应；

步骤3，计算隐层神经元的非线性激活值；

步骤4，隐层特征数据三维重组；

步骤5，空间引导滤波；

步骤6，滤波后的隐层特征数据二维重组；

步骤7，构造局部相似性保持正则项及优化模型；

步骤8，计算局部相似性保持投影特征，并进行K-means聚类得到最终的聚类标签。

2.根据权利要求1所述的局部相似性保持的高光谱图像极限学习机聚类方法，其特征在于，步骤1具体为：

输入一幅高光谱图像X₀∈R^D×W×H，其中，D表示高光谱图像的波段数，W和H分别表示图像空间维的宽度和高度；

将原始高光谱数据X₀按行逐像素排列形成矩阵X＝[x₁,x₂,...,x_N]∈R^D×N作为模型的输入，其中，N＝W×H表示高光谱像元的个数，x_i∈R^D表示一个高光谱像元。

3.根据权利要求1所述的局部相似性保持的高光谱图像极限学习机聚类方法，其特征在于，步骤2具体为：

步骤2-1，计算第j个隐层神经元对第i个高光谱像元的线性随机响应：

其中，w_j＝[w_j1,w_j2,...,w_jD]∈R^D是随机生成的权重向量，b_j＝[b_j1,b_j2,...,b_jD]∈R^D是随机偏置，D表示高光谱像元的波段数，M表示隐层神经元的个数，1≤i≤N；

步骤2-2，构造线性随机投影向量，即对于高光谱像元x_i，将M个隐层神经元拼装，构成对应的线性随机投影行向量o_i＝[o_i1,o_i2,...,o_iM]。

4.根据权利要求1所述的局部相似性保持的高光谱图像极限学习机聚类方法，其特征在于，步骤3具体为：

(1)对隐层线性随机投影行向量进行非线性映射，可得

h_i＝[g(o_i1)…g(o_iM)]^T

其中，h_i∈R^M是对应于x_i的隐层输出，g(·)是激活函数；

(2)构造整体输入数据X对应的M×N大小的隐层输出矩阵：

H＝[h₁,h₂,...,h_N]∈R^M×N。

5.根据权利要求1所述的局部相似性保持的高光谱图像极限学习机聚类方法，其特征在于，步骤4中隐层特征数据三维重组具体为：

将二维的隐层输出矩阵H∈R^M×N按列取矩阵的元素转化成立方体数据T∈R^M×W×H，其中，W和H分别表示图像空间维的宽度和高度，M表示隐层输出特征维度；

记p＝(p₁,p₂)表示空间上的任意点，其中，1≤p₁≤W，1≤p₂≤H；T(p)∈R^M表示立方体数据在该位置p的隐层输出特征向量。

6.根据权利要求1所述的局部相似性保持的高光谱图像极限学习机聚类方法，其特征在于，步骤5中空间引导滤波具体为：

(1)构造与位置p的像元相邻的像元集合其中邻域像元

(2)对任意位置p处，计算引导滤波权重，其计算式为：

w_p,q＝w_p,q-1D(q,q-1)R(p,q)

其中，D(q,q-1)表示像元q和像元q-1之间的相似性权重，R(p,q)表示像元p和像元q之前的相似性权重，具体定义如下：

其中，σ_d是D(q,q-1)的高斯核方差，σ_r是R(p,q)的高斯核方差；

(3)对立方体数据进行空间引导滤波，其计算公式为：

其中，Z(p)是归一化因子，T为原始的数据，为空间引导滤波后的数据。

7.根据权利要求1所述的局部相似性保持的高光谱图像极限学习机聚类方法，其特征在于，步骤6中滤波后的隐层特征数据二维重组具体为：

将空间引导滤波后所得的逐像元重排为二维的特征矩阵其中，将作为新的融合空间信息的隐层输出矩阵。

8.根据权利要求1所述的局部相似性保持的高光谱图像极限学习机聚类方法，其特征在于，步骤7具体为：

(1)构造光谱相似性连接图：

结合基于图的局部相似性保持方法，原始光谱数据根据任意两个光谱像元的相关性可以建模为包含N个顶点的图，每条边对应两个顶点的相似性权重，其中，权重的定义如下：

其中，s_i,j表示像元x_i和x_j的相似度，σ是带宽参数，σ≥0，表示x_i的n个最相似近邻的集合，n＝5；

(2)建立光谱相似性保持正则项：

上式等价于：

其中，Tr(·)表示矩阵的迹，L＝D-S是图拉普拉斯矩阵，D是对角矩阵，d_ii＝∑_js_ij，A是要求解的局部相似性保持投影矩阵。

(3)建立局部光谱相似性保持的最小代价函数模型：

其中，λ是正则化参数，I_C是C×C的单位矩阵；

(4)模型求解：

利用拉格朗日乘子法，上述模型的解可表述为：

其中，v₁,v₂,...,v_C+1是求解的前C+1个最小特征值γ₁,γ₂,…,γ_C+1对应的特征向量，γ₁≤γ₂≤...≤γ_C+1，是标准化的特征向量，其中，C表示输出特征的维度，I_M表示M×M的单位矩阵。

9.根据权利要求1所述的局部相似性保持的高光谱图像极限学习机聚类方法，其特征在于，步骤8具体为：计算局部相似性保持投影特征，并进行K-means聚类得到最终的聚类标签；

最终的输出特征为将的每一行看作一个像元的局部相似性保持投影特征向量，用K均值的方法将N个像元聚成k类，得到最终类别标签向量y。