技术领域本发明属于遥感图像处理技术领域，涉及一种高光谱图像混合像元分解算法，尤其涉及一种基于双重丰度约束非负矩阵分解的高光谱图像混合像元分解算法。
背景技术：
高光谱图像由于包含多个波段，含有丰富的光谱信息，可以完整地反映不同地物的诊断性光谱特征。上百个波段的光谱信息为地物识别提供了有力的识别信息。然而由于空间分辨率的限制和地物的复杂性，混合像元在高光谱图像中普遍存在。大量的混合像元阻碍了对像元所属物质类别的准确识别和解译。混合像元分解是解决混合像元问题的关键技术。其中非负矩阵分解由于其得到的矩阵的非负性，被应用到混合像元分解中。然而非负矩阵分解方法由于其非凸性，求得的结果仍然存在一些限制：1)非负矩阵分解是已知一个非负矩阵，将其分解为两个非负矩阵的乘积。这是一个非凸优化问题，存在很多局部最优解；2)为了缓解局部最优问题，一些方法对非负矩阵分解加入丰度平滑约束项进行混合像元分解。假设相似的像元具有相似的丰度，并且空间越相邻的像元越相似，利用空间邻域像元到当前像元的空间距离远近来衡量像元之间的相似性。然而，这一衡量方式在复杂的真实影像上不一定合适；3)一些算法利用高光谱图像的光谱信息来计算当前像元和空间邻域像元之间的相似性。但是对空间邻域内的所有像元都进行平滑约束，没有考虑到位于不同物质类别边界的像元和异常像元的存在而导致的不平滑现象，因此会导致一些不正确的约束。丰度平滑特征是高光谱图像的一个重要特征，但是现有方法没能充分考虑当前像元的邻域环境，可能造成一些不恰当的约束。因此，需要一种能恰当地表达高光谱图像的丰度平滑信息的方法。
技术实现要素：
为了解决上述技术问题，本发明主要提供了一种能去除不恰当的丰度平滑约束、保留相似像元的丰度平滑约束的基于双重丰度约束非负矩阵分解的高光谱图像混合像元分解算法。本发明所采用的技术方案是：一种高光谱图像混合像元分解算法，其特征在于，包括以下步骤：步骤1：估计高光谱图像的端元数，初始化端元矩阵和丰度矩阵；步骤2：计算每个像元和它的空间邻域像元的相似性；步骤3：将得到的所有像元和它们的邻域像元之间的相似性进行从小到大的排序，将排于后的部分像元对不纳入平滑约束；步骤4：计算每个像元和其邻域像元的平滑约束权值；步骤5：对丰度矩阵进行和为1约束；步骤6：进行端元矩阵的更新；步骤7：计算丰度的平滑约束梯度矩阵和分离约束梯度矩阵，进行丰度矩阵的更新；步骤8：判断迭代是否收敛，若是，则结束迭代；否则转步骤6步。作为优选，步骤1中是利用端元数目估计方法估计高光谱图像的端元数，模拟实验使用最大化光谱信息散度方法初始化端元矩阵；真实实验使用顶点成分分析方法初始化端元矩阵。利用图像数据和初始化的端元矩阵求得的无约束最小二乘结果作为初始化的丰度矩阵，并将丰度矩阵中小于0的值赋值为0。作为优选，步骤2中是利用高光谱图像的光谱特征，计算当前像元和空间邻域像元之间的光谱角来衡量它们的相似性，邻域窗口大小为5×5。作为优选，步骤3中是通过经验值将排于后40％的像元对不纳入平滑约束。作为优选，步骤4中是将相似性值输入热核公式进行计算，得到每个像元和其邻域像元的平滑约束权值。作为优选，热核参数通过计算当前像元和它的所有邻域像元的相似性平均值来自动获取。作为优选，步骤6中是用梯度下降法进行端元矩阵的更新。作为优选，步骤7中是用梯度下降法进行丰度矩阵的更新。作为优选，步骤8中所述判断迭代是否收敛，其判断准则为：若影像的重构残差小于设定的阈值或者达到最大迭代次数，则迭代收敛。本发明相对于现有技术，其有益效果为：(1)利用高光谱图像的光谱信息来衡量每个像元对所对应的丰度之间的平滑程度，使光谱越相似的像元对所对应的丰度之间有更高的平滑权值。这种方法比利用空间距离决定平滑权值更接近真实情况；(2)将与当前像元不相似的空间邻域像元不纳入平滑约束，以此来避免由于边界点和异常点的存在带来的不平滑现象；(3)引入丰度分离约束，来避免产生过平滑的丰度估计结果。丰度分离约束和平滑约束之间的有效权衡可以使得混合像元分解的结果更稳定；(4)实现了更加准确地描述高光谱图像的丰度特征，以这些特征对非负矩阵分解进行约束，可以得到更好的丰度估计结果。附图说明图1：为本发明实施例的丰度平滑约束原理示意图；图2：为本发明实施例的丰度平滑约束原理中邻域内有边界点的平滑约束情况示意图；图3：为本发明实施例的丰度平滑约束原理中邻域内有异常点的平滑约束情况示意图；图4：为本发明实施例的丰度分离约束原理示意图。具体实施方式为了便于本领域普通技术人员理解和实施本发明，下面结合附图及实施例对本发明作进一步的详细描述，应当理解，此处所描述的实施示例仅用于说明和解释本发明，并不用于限定本发明。本发明提供的一种高光谱图像混合像元分解算法，其特征在于，包括以下步骤：步骤1：利用端元数目估计方法估计高光谱图像的端元数，模拟实验使用最大化光谱信息散度方法初始化端元矩阵；真实实验使用顶点成分分析方法初始化端元矩阵。利用图像数据和初始化的端元矩阵求得的无约束最小二乘结果作为初始化的丰度矩阵，并将丰度矩阵中小于0的值赋值为0。步骤2：利用高光谱图像的光谱特征，计算当前像元和空间邻域像元之间的光谱角来衡量它们的相似性，邻域窗口大小为5×5；步骤3：将得到的所有像元和它们的邻域像元之间的相似性进行从小到大的排序，通过经验值将排于后40％的像元对不纳入平滑约束；步骤4：将相似性值输入热核公式进行计算，得到每个像元和其邻域像元的平滑约束权值；考虑到人工设置统一的热核参数不能较好地表达数据的结构特征，热核参数通过计算当前像元和它的所有邻域像元的相似性平均值来自动获取。步骤5：对丰度矩阵进行和为1约束；步骤6：用梯度下降法进行端元矩阵的更新；步骤7：计算丰度的平滑约束梯度矩阵和分离约束梯度矩阵，用梯度下降法进行丰度矩阵的更新；请见图1-3。图1中可以看出对于异常点和边界点来说，它们与有的邻域像元并不相似，如果对它们加入丰度平滑约束，会导致不正确的结果。图2和图3分别表示邻域内有边界点和异常点的平滑约束情况。圆圈表示邻域像元和中心像元并不属于同一类，所以对它们的平滑约束权值赋为0；星号表示邻域像元和中心像元属于同一类，但是每个像元对的相似性不一样，所以为相似性较高的像元对赋予较大的平滑约束权值(用较大的星号表示)。请见图4，可以看出，属于不同类别的像元所对应的丰度在特征空间中有不同的分布领域，所以可以认为，不同的丰度向量之间是有分离性质的，即弱相关性。因此加上丰度的分离约束来避免过平滑约束现象。将上述两个约束项设计完毕后，即可将它们加入到非负矩阵分解的基本模型里作为正则项，利用梯度下降法来进行目标函数的求取。步骤8：判断迭代是否收敛，若是，则结束迭代；否则转步骤6步。判断迭代是否收敛，其判断准则为：若影像的重构残差小于设定的阈值或者达到最大迭代次数，则迭代收敛。以上是本发明涉及的高光谱图像混合像元分解的实现步骤。通过丰度平滑约束和丰度分离约束的引入，可以充分利用高光谱图像丰度的特征来进行混合像元分解。在具体实现的时候还有以下注意事项：丰度根据其物理性质需要满足大于等于0与和为1的条件。非负矩阵分解本身的非负性满足了第一个条件，第二个条件在求解过程中是通过对端元矩阵和高光谱数据矩阵加入一个全1行来达到的。应当理解的是，本说明书未详细阐述的部分均属于现有技术。应当理解的是，上述针对较佳实施例的描述较为详细，并不能因此而认为是对本发明专利保护范围的限制，本领域的普通技术人员在本发明的启示下，在不脱离本发明权利要求所保护的范围情况下，还可以做出替换或变形，均落入本发明的保护范围之内，本发明的请求保护范围应以所附权利要求为准。