本发明涉及高光谱遥感图像数据处理领域，特别是指一种高光谱影像混合像元分解方法。
背景技术：
高光谱遥感近年来得到突飞猛进的发展，图像中不仅仅包含了全色及彩色摄影丰富的空间信息，同时相对于多光谱图像其拥有更多更精细的光谱信息，使得图像中每一个像元都对应一条平滑而完整的光谱曲线。鉴于不同地物光谱曲线属性不同，利用高光谱遥感技术能够充分挖掘不同物质的光谱特征及形态特征，这为地物的精细探测奠定了基础。然而由于地物复杂性的影响以及空间分辨率的限制，混合像元问题在高光谱遥感图像中是不可避免的，而感兴趣目标往往处于亚像元级的状态，因此需要研究混合像元分解技术，提取不同地物端元将目标光谱与背景光谱分离开，对于提高遥感定量化应用及亚像元级目标探测具有重要意义。
根据混合像元产生机理，忽略光谱间的非线性混合，假设像元内地物为线性混合，则线性模型中观测光谱向量可表述为
y＝ms(1)
式中，
现有技术中使用的高光谱混合像元分解算法没有针对高光谱影像端元光谱类间差异做进一步处理，存在高光谱影像混合像元分解精度低的问题。
技术实现要素：
本发明要解决的技术问题是提供一种高光谱影像混合像元分解方法，以解决现有技术所存在的高光谱影像混合像元分解精度低的问题。
为解决上述技术问题，本发明实施例提供一种高光谱影像混合像元分解方法，包括：
获取高光谱影像数据；
根据端元光谱变异情况，构造类内端元光谱变异矩阵；
利用构造的类内端元光谱变异矩阵对所述高光谱影像数据进行加权处理；
引入丰度稀疏约束条件并基于加权处理结果，构建目标函数，利用分层式非负矩阵分解策略对构建的目标函数进行分解，得到端元矩阵及丰度矩阵。
进一步地，在获取高光谱影像数据之后，所述方法还包括：
对所述高光谱影像数据中每个像元光谱进行平滑处理。
进一步地，所述对所述高光谱影像数据中每个像元光谱进行平滑处理包括：
利用基于时域多项式卷积计算方法对所述高光谱影像数据中每个像元光谱进行最小二乘拟合计算，其中，影像边界点则以相应像元为中心2邻域空间像元的五点三次拟合的平均值作为初始值，推导出边界数据结果。
进一步地，所述影像边界点的求解公式为：
其中，a为平滑前的光谱曲线，ai(j)为第i个像元第j个波段光谱反射率的值，b为平滑后的光谱曲线，b(j)为第j个波段反射率的值。
进一步地，所述构造的类内端元光谱变异矩阵表示为：
其中，a表示类内端元光谱变异矩阵，m表示第p类端元的光谱曲线，μj表示第p类端元下cj个训练样本光谱均值。
进一步地，利用构造的类内端元光谱变异矩阵对所述高光谱影像数据进行加权处理，得到加权后的高光谱影像数据和加权后的端元矩阵：
其中，a表示类内端元光谱变异矩阵，m表示第p类端元的光谱曲线，μj表示第p类端元下cj个训练样本光谱均值，y表示高光谱影像数据，m表示端元矩阵，s表示丰度矩阵，
进一步地，构建的目标函数表示为：
其中，j(m,s)表示目标函数，λ表示丰度稀疏约束系数。
进一步地，所述分层式非负矩阵分解策略表示为：
其中，mi和si分别是第i层分解得到的端元矩阵和丰度矩阵，l表示非负矩阵分解层数。
进一步地，在引入丰度稀疏约束条件并基于加权处理结果，构建目标函数，利用分层式非负矩阵分解策略对构建的目标函数进行分解，得到端元矩阵及丰度矩阵之后，所述方法还包括：
利用t分布，设置置信度区间；
计算端元间光谱角距离，低于置信区间内的端元视为冗余端元；
若有冗余端元，则根据预先设置的取舍冗余端元的阈值，去除冗余端元并相应减少端元数目，再执行利用分层式非负矩阵分解策略对构建的目标函数进行分解的步骤。
进一步地，取舍冗余端元的阈值和光谱角距离分别表示为：
其中，
本发明的上述技术方案的有益效果如下：
上述方案中，获取高光谱影像数据；根据端元光谱变异情况，构造类内端元光谱变异矩阵；利用构造的类内端元光谱变异矩阵对所述高光谱影像数据进行加权处理；引入丰度稀疏约束条件并基于加权处理结果，构建目标函数，利用分层式非负矩阵分解策略对构建的目标函数进行分解，得到端元矩阵及丰度矩阵，完成混合像元分解；这样，在非负矩阵分解算法中采用分层策略，并引入丰度稀疏约束，能够增大高光谱影像端元光谱类间差异，从而使得能够提取纯净端元光谱曲线，提高高光谱影像混合像元分解精度，为实现精确地物分类及目标探测奠定基础。
附图说明
图1为本发明实施例提供的高光谱影像混合像元分解方法的流程示意图；
图2(a)为本发明实施例提供的平滑前的光谱曲线示意图；
图2(b)为本发明实施例提供的平滑后的光谱曲线示意图；
图3为本发明实施例提供的随机选取的端元光谱曲线图；
图4(a)为本发明实施例提供的模拟数据端元提取结果示意图；
图4(b)为本发明实施例提供的模拟数据端元丰度估计结果示意图；
图5为本发明实施例提供的不同信噪比光谱角距离结果示意图；
图6为本发明实施例提供的不同信噪比下均方误差示意图；
图7为本发明实施例提供的真实数据端元提取结果示意图。
具体实施方式
为使本发明要解决的技术问题、技术方案和优点更加清楚，下面将结合附图及具体实施例进行详细描述。
本发明针对现有的高光谱影像混合像元分解精度低的问题，提供一种高光谱影像混合像元分解方法。
实施例一
如图1所示，本发明实施例提供的高光谱影像混合像元分解方法，包括：
s101，获取高光谱影像数据；
s102，根据端元光谱变异情况，构造类内端元光谱变异矩阵；
s103，利用构造的类内端元光谱变异矩阵对所述高光谱影像数据进行加权处理；
s104，引入丰度稀疏约束条件并基于加权处理结果，构建目标函数，利用分层式非负矩阵分解策略对构建的目标函数进行分解，得到端元矩阵及丰度矩阵。
本发明实施例所述的高光谱影像混合像元分解方法，获取高光谱影像数据；根据端元光谱变异情况，构造类内端元光谱变异矩阵；利用构造的类内端元光谱变异矩阵对所述高光谱影像数据进行加权处理；引入丰度稀疏约束条件并基于加权处理结果，构建目标函数，利用分层式非负矩阵分解策略对构建的目标函数进行分解，得到端元矩阵及丰度矩阵，完成混合像元分解；这样，在非负矩阵分解算法中采用分层策略，并引入丰度稀疏约束，能够增大高光谱影像端元光谱类间差异，从而使得能够提取纯净端元光谱曲线，提高高光谱影像混合像元分解精度，为实现精确地物分类及目标探测奠定基础。
在前述高光谱影像混合像元分解方法的具体实施方式中，进一步地，在获取高光谱影像数据之后，所述方法还包括：
对所述高光谱影像数据中每个像元光谱进行平滑处理。
本实施例中，通过对实际利用高光谱成像仪拍摄的高光谱影像数据分析，由于混合像元分解中端元提取精度受光谱高频噪声影响较大，高频噪声的影响增加了端元变异的可能性，从而导致端元提取精度差的问题。因此，在获取高光谱影像数据之后，可以对高光谱影像数据中每个像元光谱进行平滑处理，去除像元光谱的高频噪声，得到平滑后的高光谱影像数据；平滑处理的具体步骤可为：
利用基于时域多项式卷积计算方法对高光谱影像数据中每个像元光谱进行最小二乘拟合计算，其中，影像边界点则以相应像元为中心2邻域空间像元的五点三次拟合的平均值作为初始值，推导出边界数据结果，从而对高光谱影像中的像元光谱进行平滑处理。
本实施例中，五点三次拟合是指通过选取相邻5个波段的数据点，拟合出一条三次曲线，然后用三次曲线上的相应的位置的数值作为光谱平滑滤波的结果。在影像边界值选取上，可以采用本实施例所述的改进的五点三次算法。改进的五点三次算法是在原算法的基础上以相应像元为中心的2邻域空间相邻像元求取平均值，影像边界点的求解公式如式(3-1)和式(3-2)所示，同理可推算第n-1和第n个边界点，该算法能够考虑影像边界点的数据，利用数据拟合结果推导出数据结果，平滑前、后的光谱曲线效果图如图2(a)、2(b)所示。
式(3-1)和式(3-2)中，a为平滑前的光谱曲线，ai(j)为第i个像元第j个波段光谱反射率的值，b为平滑后的光谱曲线，b(j)为第j个波段反射率的值。
本实施例中，针对平滑后的高光谱影像数据，考虑端元光谱变异情况，构造类内端元光谱变异矩阵，其中，构造的类内端元光谱变异矩阵表示为：
其中，a表示类内端元光谱变异矩阵，m表示第p类端元的光谱曲线，μj表示第p类端元下cj个训练样本光谱均值。
本实施例中，利用构造的类内端元光谱变异矩阵对所述高光谱影像数据进行加权处理，得到加权后的高光谱影像数据和加权后的端元矩阵，从而得到原高光谱影像数据的线性混合模型：
其中，a表示类内端元光谱变异矩阵，m表示第p类端元的光谱曲线，μj表示第p类端元下cj个训练样本光谱均值，y表示高光谱影像数据(是以矩阵的形式读入的，因此，高光谱影像数据也可以称为高光谱矩阵)，m表示端元矩阵，s表示丰度矩阵，
本实施例中，在非负矩阵分解的基础上，引入加权后的高光谱影像数据、端元矩阵以及丰度稀疏约束条件，构建目标函数，增大端元类间差异，去除类间端元变异影响，并利用分层式非负矩阵分解策略求解目标函数的最小值，从而同时获得最优的端元矩阵及相应丰度矩阵，完成混合像元分解。
本实施例中，构建的目标函数表示为：
其中，j(m,s)表示目标函数，λ表示丰度稀疏约束系数。
在非负矩阵求解过程中，相对于传统非负矩阵分解的单层处理，本发明采用分层式非负矩阵分解策略，具体的：采用分层式将矩阵逐层分解，设定分解层数，每层分解中采用相同的损失目标函数及迭代更新法则，算法终止条件为算法更新的误差达到一定的阈值或最大迭代次数。
本实施例中，所述分层式非负矩阵分解策略表示为：
其中，mi和si分别是第i层分解得到的端元矩阵和丰度矩阵，l表示非负矩阵分解层数。
本实施例中，通过分层分解，不仅能够减小求解过程中陷入局部极小的风险，同时能够节约算法运行时间及提高算法精度及非负矩阵分解求解性能。
本实施例中，每层使用相同的迭代准则以及终止算法条件，迭代准则更新公式为：
其中，←、.*、./分别表示迭代替换、矩阵对应元素相乘与矩阵对应元素相除，l表示第l层，t表示矩阵转置，α和τ是正则化系数约束的常数，t是迭代次数。
在前述高光谱影像混合像元分解方法的具体实施方式中，进一步地，在引入丰度稀疏约束条件并基于加权处理结果，构建目标函数，利用分层式非负矩阵分解策略对构建的目标函数进行分解，得到端元矩阵及丰度矩阵之后，所述方法还包括：
利用t分布，设置置信度区间；
计算端元间光谱角距离，低于置信区间内的端元视为冗余端元；
若有冗余端元，则根据预先设置的取舍冗余端元的阈值，去除冗余端元并相应减少端元数目，再执行利用分层式非负矩阵分解策略对构建的目标函数进行分解的步骤。
本实施例中，采用基于多元回归理论的最小误差高光谱信号辨识方法估计噪声。在获得较小的均方误差的基础上，完成噪声的估计以及信号子空间的确认，其中端元数目为p。
然后，利用光谱角距离(sad)约束对重复性端元进行去除，取舍冗余端元的阈值的选取避免人为设定，本发明方法利用t分布，计算端元间光谱角距离，设置置信区间为85％，低于置信区间内的端元视为冗余端元，若有冗余端元，则根据预先设置的取舍冗余端元的阈值，去除冗余端元并相应减少端元数目，然后再执行利用分层式非负矩阵分解策略对构建的目标函数进行分解的步骤。
本实施例中，光谱角距离主要是使用目标光谱向量与参考光谱向量夹角公式进行计算。
本实施例中，取舍冗余端元的阈值和光谱角距离分别表示为：
其中，
综上，本发明实施例所述的高光谱影像混合像元分解方法，采用线性混合模型，根据高光谱影像端元的光谱特性和丰度特点，利用基于时域多项式卷积计算方法对像元光谱进行最小二乘拟合计算，实现高光谱影像中的像元光谱平滑处理，去除像元光谱高频噪声，并在传统非负矩阵分解算法模型中，利用分层式非负矩阵分解策略，引入丰度稀疏约束条件，从而增大高光谱影像端元光谱类间差异，能够有效提取端元及丰度信息，实现高光谱遥感影像数据的混合像元分解。本发明具有算法模型简单的特点，在高光谱影像高精度地物分类及目标探测方面具有重要的应用价值。且本发明实施例所述的高光谱影像混合像元分解方法具有无需目标先验知识以及纯净像元的假设，具有算法速度快的特点。
实施例二
为了更好地理解本发明实施例所述的高光谱影像混合像元分解方法，分别用实验室仿真数据以及实际高光谱数据说明具体实施方式，并与经典混合像元分解算法进行对比。
1)实验室仿真数据
仿真数据的生成方式，在本发明实施例中可分为两种：一种为通过地物光谱仪采集的光谱数据，另一种为调用现有的光谱数据库。
本实施例中，通过野外实验采集的光谱数据由于数据采集量大，不同时间、不同高度均可获得大量的变异光谱信息，由此可计算出类内端元光谱变异矩阵；对于现有光谱库中的文件，对于同一种地物光谱唯一或样本较少，这时可设置类内端元光谱变异矩阵为单位阵；然后，在利用分层式非负矩阵分解策略进行分解，得到端元矩阵及丰度矩阵。
在本次仿真数据实验中，本实施例通过预先设置端元个数，随机调用美国地质勘探局(usgs)光谱，信噪比40，波段数224，像素总数3364，丰度值服从狄利克雷分布。
为了保证图像中不包含纯像元，把每个像元中包含任意某种端元比例值高于0.8的舍弃掉。图3为从标准光谱库中随机选取5条波谱曲线并作为高光谱影像数据中的端元组分。
本实施例中，利用光谱角距离(sad)和均方误差(rmse)作为混合像元分解的评价指标，其中，均方误差是用来衡量光谱分解后生成丰度反演影像与原参考影像之间的差异。假设
εi＝xi-asi
整幅影像的均方根误差为：
其中，n为像元总个数，l为波段数，εi,j为第i个像素第j个波段余差。
图4(a)-(b)为利用本发明方法对模拟图像混合像元分解的结果，图4(a)为端元提取的结果，实线为真实端元光谱曲线，虚线为本发明方法提取的端元光谱曲线，图4(b)选取其中三个端元丰度为例，其中，第一行为真实端元光谱图，第二行为本发明估计的丰度图。由图4(a)-(b)可看出，本发明方法能够得到精确的结果。
由图5和图6可以看出，本发明方法可同时获取较精确的端元光谱曲线及丰度估计图像；其中，图5中，iea表示迭代误差分析，vca表示顶点成分分析，n-findr表示内部最大体积法，osp表示正交子空间投影，ica表示独立成分分析，atgp表示自动目标生成方法。
2)实际高光谱数据
真实数据实验采用美国aviris机载数据拍摄的cupritenevada矿物地区，图像大小250*191像元，去除噪声以及大气水汽吸收波段最终波段数目为188，波长范围400nm到2500nm。在本次实验中，设置分解层数为10，图7为本发明方法提取端元的示意图，利用光谱角距离及均方误差公式计算误差分别为0.1029和0.056，能够得到较好的效果。
需要说明的是，在本文中，诸如第一和第二等之类的关系术语仅仅用来将一个实体或者操作与另一个实体或操作区分开来，而不一定要求或者暗示这些实体或操作之间存在任何这种实际的关系或者顺序。
以上所述是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明所述原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。