[发明专利]一种基于随机投影的高光谱图像稀疏解混方法

摘要

一种基于随机投影的高光谱图像稀疏解混方法，它有四大步骤：一、计算机在MATLAB??R2008b环境下读取数据；二、计算机对高光谱图像数据和高光谱库数据进行随机投影；三、构建稀疏解混的目标函数，使用分裂Bregman算法优化目标函数求极值，直至达到收敛停止条件。四、设定合适的阈值处理丰度矩阵，获得最终的丰度图和端元。本发明利用了高光谱数据库来选择端元，克服了以往算法所求出的端元与标准高数据库中的纯物质光谱无法严密对应的缺点；并使用随机投影技术对原始数据进行降维，从而达到了节省内存和减少计算量的效果。本发明实现了对高光谱图像快速定量分析，它在高光谱遥感图像分析领域里具有实用价值和广阔的应用前景。

主权项

一种基于随机投影的高光谱图像稀疏解混方法，其特征在于：该方法具体步骤如下：步骤一：用计算机读取数据：计算机在MATLAB R2008b环境下读取高光谱图像数据，数据来源于成像光谱仪采集到的遥感图像，得到数据立方体，高光谱图像数据去除被水汽吸收的波段和信噪比较低的波段，将高光谱图像数据逐像素点排列，得到原始的高光谱图像矩阵V；设高光谱图像有m个波段，共有n个像素点，V＝[v1，v2，...，vn]，vi是高光谱图像第i个像素点的光谱列向量，是一个m维向量；读取现有的高光谱数据库数据，选择光谱库中纯物质光谱数据构建光谱库矩阵S；设纯物质光谱的数量为q，S＝[s1，s2，...，sq]，si是光谱库中第i个纯物质的光谱列向量，同样是一个m维向量；步骤二：高光谱数据随机投影：用计算机读取数据后，对高光谱图像数据和高光谱库数据进行随机投影，就是用一个所包含的元素符合正态分布的矩阵Rd与矩阵V，S分别相乘；设矩阵Rd的大小为d×m，d＜m，得到降维后的数据RdV和RdS，其矩阵大小分别为d×n，d×q，代替原始的高光谱图像数据V和高光谱数据库S参与下述步骤三的计算，其中d是随机投影后光谱波段的个数；步骤三：高光谱图像稀疏解混：数据进行随机投影后，就对降维后的高光谱图像数据进行光谱解混；设高光谱数据库中含有足够丰富的纯物质光谱特征，那么高光谱图像中含有的端元只占极少一部分，也就是说高光谱图像中每个像素点的光谱曲线是由高光谱数据中少数的纯物质的光谱曲线线性组合构成，这体现了端元的稀疏性表达；端元在高光谱数据库中是稀疏的，他们所对应的丰度也具有稀疏性，即丰度矩阵是稀疏的，以丰度的稀疏性作为正则化项，符合实际的物理意义；在这里，采用线性混合模型，即将一条检测得到的光谱曲线，分解为光谱库中纯物质光谱线性组合的形式，其系数即为相应丰度；设丰度矩阵为H，其大小为q×n，其所有元素满足非负性，线性混合模型如下表示：V＝SH对原始高光谱图像稀疏解混的目标函数为： min H | | H | | 1 s.t.V＝SH，H≥0||H||1表示H包含的所有元素绝对值的和，因为丰度具有非负性，||H||1越小，丰度矩阵H含有的非零元素值将越少，其表现出的稀疏性就较强；将上式转变为无约束优化问题： min H 1 2 | | V - SH | | 2 2 + λ | | H | | 1 在上式目标函数中，l2项是拟合项，体现的是解混的结果与原高光谱图像的拟合程度；l1即稀疏性正则化项；λ是正则化参数，调节着拟合项与正则化项的权重，其大小根据实际高光谱数据设置；因为V，S数据量很大，在迭代求解的过程中占用内存较大，并且计算时间较长；采用随机投影算法，能够有效保持上式中的l2距离项，实现数据量压缩，随机投影后的目标函数为： min H 1 2 | | R d V - R d SH | | 2 2 + λ | | H | | 1 针对上述目标函数，采用分裂Bregman算法求解上述函数极值，从而得到该优化模型下的最优解——最优丰度矩阵H；引入辅助变量D，并且D＝H，上式目标函数转化为： min H , d 1 2 | | R d V - R d SH | | 2 2 + λ | | D | | 1 s . t . D = H 转化无约束形式为： min H , d λ | | D | | 1 + μ 2 | | R d V - R d SH | | 2 2 + 1 2 | | D - H | | 2 2 其中μ是正则化参数，其大小根据实际高光谱数据设置。，引入误差反馈变量B，上述优化问题转变为迭代求解过程，如下： H k + 1 = arg min { μ 2 | | R d V - R d SH | | 2 2 + 1 2 | | D k - H - B k | | 2 2 } D k + 1 = arg min { λ | | D | | 1 + 1 2 | | D - H k + 1 - B k | | 2 2 } 具体求解步骤如下：(1)初始化丰度矩阵H以及辅助变量矩阵B，D，Vd，H0＝B0＝D0＝0；Vd0＝RdV(2)迭代变量H，D，B，Vd，交替迭代过程如下： H k + 1 = [ μ ( R d S ) T ( R d S ) + I ] - 1 [ μ ( R d S ) T V d k - B k + D k ] D k + 1 = max ( H k + 1 + B k - λ , 0 ) B k + 1 = B k + H k + 1 - D k + 1 V d k + 1 = V d k + R d V - R d SH k + 1 这里k为迭代次数，Vdk+1的迭代过程可以加速收敛；(3)停止迭代条件：当丰度矩阵H相邻两次迭代的相对误差不大时，迭代停止；迭代停止条件为： | | H k + 1 - H k | | 2 2 | | H k | | 2 2 < tol 这里，取tol＝10 6，若迭代步数少于10次就能收敛，使用迭代次数控制迭代过程，设迭代次数为1000步；步骤四：获取丰度图及真正的端元：在获得丰度矩阵后，设置合适的阈值，丰度矩阵中小于阈值的元素设为零，不小于阈值的元素不做处理；找出丰度矩阵中含有非零元素的行标号，丰度矩阵的行标号与光谱库的列标号对应；取出光谱矩阵中列标号所对应的列，即为真正的端元光谱；取出丰度矩阵中行标号所对应的行，即为端元所对应的丰度图；在丰度图中颜色较亮的区域表示端元所占的比例较大，较暗的区域对应的比例小，颜色为黑色的区域表示不含有该端元。