[0057] 以下通过特定的具体实施例说明本发明的实施方式，本领域技术人员可由本说明书所揭露的内容轻易地了解本发明的其他优点与功效。本发明还可以通过另外不同的具体实施方式加以实施或应用，本说明书中的各项细节也可以基于不同观点与应用，在没有背离本发明的精神下进行各种修饰或改变。需说明的是，在不冲突的情况下，以下实施例及实施例中的特征可以相互组合。
[0058] 由于传统基于非负矩阵分解的高光谱解混方法容易陷入局部最小值，本发明提供了一种基于约束非负矩阵分解的高光谱解混方法，将丰度稀疏、平滑约束、端元平滑约束整合到经典NMF方法的目标函数中，经典NMF方法的目标函数为：
[0059]
[0060] 具体地，由于高光谱图像中每个像素通常是由较少的端元光谱混合而成，而不是全部端元混合而成，因此，丰度的列是稀疏的。为了考虑丰度列的稀疏性，传统方法采用L0，L1范数，L0虽然可以获得稀疏解，但存在NP难问题。对于L1范数，有时不满足丰度和为1约束。本发明采用加权稀疏正则化来约束丰度矩阵的稀疏性，获得比L1更稀疏且满足丰度和为1的解。利用全变差(TV)正则化来促进丰度图的平滑特性，采用马尔可夫随机场模型的自适应势函数来保证端元光谱的平滑性，以此提高解混的精度，本发明可以同时获取端元和丰度，相比其他传统方法，减少求解丰度时对端元提取精度的依赖。本发明方法较好的实现了高光谱图像解混，在端元提取和丰度估计方面都取得了较好的结果。
[0061] 参照图1，本实施例提供了一种基于约束非负矩阵分解的高光谱解混方法，本实施例提供了一种基于约束非负矩阵分解的高光谱解混方法，针对具有N个像素，每个像素为BB×N维光谱向量的高光谱数据矩阵Y∈R ，其中表示光谱波段数，N表示高光谱图像中所有像素个数，端元数量为K，高光谱图像的空间尺寸大小为m×n。
[0062] 方法包括步骤：
[0063] S1、将丰度矩阵的稀疏性约束、丰度图的平滑性约束以及端元的平滑性约束整合到NMF方法的目标函数中，以得到所需的整体目标函数；
[0064] S2、初始化端元矩阵A、丰度矩阵S、辅助矩阵L和权重矩阵W；
[0065] S3、对端元矩阵A、丰度矩阵S、辅助矩阵L和权重矩阵W进行迭代更新求解，并判断是否达到终止条件，若未达到则利用更新后的结果继续进行迭代更新，若达到则停止迭代并输出最终的端元矩阵以及丰度矩阵。
[0066] 具体的，
[0067] (一)：
[0068] 1、步骤S1中，对于丰度矩阵的稀疏性约束，采用的加权稀疏约束定义为：
[0069]
[0070] 其中，ω∈RK表示K维的权值向量，用于控制丰度向量s的稀疏性，K表示端元个数；B B×K
y∈R表示高光谱图像中某一像素的B维光谱向量，B表示光谱波段数，A∈R 表示整个高光谱图像中的K个B维端元光谱向量组成的端元光谱矩阵，每一列为一个端元的B维光谱向量，K
s∈R 为K维丰度向量，表示K个端元在同一像素中各占的比例，该向量中元素值总和为1，R表示实数域，||·||1表示向量的1范数。s.t.后面的内容表示约束条件。
[0071] 通过证明，在合适的权值向量下，加权L1正则化器可以比普通L1正则化器获得更稀疏的解，本实施例利用迭代重加权算法解决加权L1最小化问题，其中根据当前丰度矩阵计算下一次迭代的权重，即：
[0072]
[0073] 其中 表示第k次迭代的丰度矩阵，i,j为矩阵S的第i行第j列元素，表示第i个端元在第j个像素中所占的比例，eps为一很小的正数以避免除数为0， 表示第k+1次迭代的权重矩阵，因此将加权L1正则化引入到NMF模型中，构成RSNMF模型。
[0074] 2、由于全变差(TV)范数促使相邻区域具有相似的值，因此对于丰度图的平滑性约束，采用全变差正则化来促进丰度图的平滑性，对于一幅大小为m×n的单波段图像y，将其表示成矩阵形式，其TV范数可以定义为
[0075]
[0076] 其中，yi,j表示y的第i行第j列元素值，|·|表示绝对值运算；
[0077] 高光谱丰度矩阵的全变差范数可以表示为K个单波段图像的TV范数之和，其定义为：
[0078]
[0079] 其中，S∈RK×N表示K行N列的丰度矩阵，K表示端元个数，N表示高光谱图像中的像素j数，其中Si,j为第i行第j列元素，表示第i个端元在第j个像素中所占的比例，S 表示丰度矩j K×N
阵的第j行，F将S这一行的N个元素转化成m×n的矩阵，即N＝m×n，S∈R 的K行经过F变换K×N
之后可以表示为K个m×n的矩阵，然后分别求TV范数，HTV表示丰度矩阵S∈R 的TV范数，即K个TV范数之和，TV表示大小为m×n的矩阵的TV范数。
[0080] 将TV正则化整合到RSNMF模型中，构成TV‑RSNMF模型，其目标函数为：
[0081]
[0082]
[0083] 其中，λ表示丰度稀疏正则化参数，λ越大，表示结果越稀疏。τ表示丰度平滑正则化参数，τ越大，丰度图越平滑，W为用于控制丰度矩阵S的权值矩阵，T表示矩阵的转置。
[0084] 显然，TV‑RSNMF模型的目标函数对于A,S来说是非凸的，为了有效的求解该问题，文中引入辅助变量L，则目标函数转化为：
[0085]
[0086]
[0087] 其中可以将S看作L的噪声版本，将S＝L约束变成无约束，可以得到：
[0088]
[0089]
[0090] μ控制L和S的相似程度，μ越大，两者相似度越高。
[0091] 3、对于端元的平滑性约束，端元光谱a为端元矩阵A的某一列，即某一种地物的光谱曲线。g(a‑aN)反映光谱特征a的平滑性，其第i项g(ai‑aNi)定义为：
[0092]
[0093] 其中Ni＝{i‑1,i+1}，g(·)为平滑函数，ai表示a中第i个元素。
[0094] 为了表征光谱数据的平滑性，本文平滑函数g(·)采用马尔科夫随机场(MRF)模型的自适应势函数，其定义为：
[0095]
[0096] 其中常数1保证结果是非负的，γ控制模型的平滑度。
[0097] 因此端元光谱平滑约束可以表示为：
[0098] L(A)＝β，
[0099] 其中AN中各位置的元素为A中对应位置元素的邻域，且A中某一列a的第i个元素对应的邻域为i的上下两个元素，即Ni＝{i‑1,i+1}，其中，某一列的第一个和最后一个元素的邻域分别为该列第2个元素和倒数第二个元素，其余均为上下两个元素，β表示端元光谱平滑的正则参数，用于控制端元光谱平滑约束的影响，β越大，端元光谱越平滑，<·>表示矩阵元素的总和，g(·)为平滑函数。
[0100] 综上所述，本方案的目标函数为：
[0101]
[0102]
[0103] 其中，Y∈RB×N表示具有N个像素，每个像素为B维光谱向量的高光谱数据矩阵，B表示光谱波段数，N表示像素数，λ表示丰度稀疏正则化参数，λ越大，表示结果越稀疏，μ表示惩罚参数，用于控制矩阵L和S的相似程度，τ表示丰度平滑正则化参数，τ越大，丰度图越平滑，β表示端元光谱平滑正则化参数，T表示矩阵的转置。
[0104] (二)：
[0105] 步骤S2中，首先利用传统的端元提取算法顶点成分分析方法(VCA)提取端元矩阵，在端元提取过程中，由于VCA向随机方向上投影，为获得更精确的初始化，本实施例迭代运行30次VCA，选择具有体积的最大单形体对应的顶点作为解混过程端元的初始化。接着利用约束最小二乘反演丰度矩阵S，然后根据丰度矩阵S初始化辅助矩阵L和权重矩阵W。
[0106] A、利用传统的端元提取算法顶点成分分析方法(VCA)提取端元矩阵，重复30次，选择具有体积的最大单形体对应的顶点作为解混过程端元的初始化；
[0107] B、根据A中获取的端元矩阵，利用约束最小二乘法反演丰度矩阵，作为解混过程丰度矩阵的初始化；
[0108] C、利用丰度矩阵S初始化辅助矩阵；
[0109] D、利用公式 初始化权重矩阵W。
[0110] 其中eps表示很小的正数，避免除数为0，Si,j表示矩阵S的第i行第j列元素，K表示第K次迭代的结果。
[0111] (三)：
[0112] 步骤S3中，具体包括以下步骤：
[0113] S3.1、根据权重矩阵W的迭代公式对W进行更新；
[0114] S3.2、根据端元矩阵A的迭代公式对A进行更新；
[0115] S3.3、将步骤S3.1以及S3.2中更新后的结果用于丰度矩阵S的更新迭代中，以对S进行更新；
[0116] S3.4、根据辅助矩阵L的迭代公式对L进行单独求解，以对L进行更新；
[0117] S3.5、判断是否达到终止条件，若未达到则根据更新后的结果重复步骤S3.1‑S3.4，若达到则停止迭代并输出最终的端元矩阵以及丰度矩阵；
[0118] 其中：
[0119] 1、权重矩阵W的更新规则：
[0120] 权重矩阵W的迭代公式为：
[0121]
[0122] 其中 表示第k次迭代的丰度矩阵，i,j为矩阵S的第i行第j列元素，表示第i个端元在第j个像素中所占的比例，eps为一很小的正数以避免除数为0， 表示第k+1次迭代的权重矩阵；
[0123] 2、端元矩阵A的更新规则：
[0124] 将整体目标函数对A求梯度，然后结合KKT条件，即可得到端元矩阵A的迭代公式：
[0125]
[0126] 其中，每次迭代中(k
γ控制平滑程度，Ni＝{i‑1,i+1}表示端元光谱矩阵A中任意一列a的第i个元素的邻域,A+1)
表示第k+1次迭代的端元光谱矩阵；
[0127] 3、对于丰度矩阵S的更新规则：
[0128] 为了满足丰度和为1约束(ASC)，在每次迭代更新S之前需要计算高光谱数据矩阵Y(k)和端元矩阵A 的增广矩阵，如下式所示：
[0129]
[0130] 其中δ控制ASC约束的影响，1表示全1的行向量，将整体目标函数对S求梯度，结合KKT条件可得丰度矩阵S的迭代公式为：
[0131]
[0132] 其中S(k+1)表示第k+1次迭代的丰度矩阵。
[0133] 4、对于辅助矩阵L的更新迭代公式，转化为求解下式右侧部分
[0134]
[0135] 其中 表示先求Frobenius范数再进行平方运算，K表示端元个数，F将行向量的N个元素表示成m×n的矩阵，利用FGP算法对辅助矩阵L的更新迭代公式右侧部分求解K次，每(k) j次求解得到L的一行，经过K次求解，即完成一次L矩阵的更新，(L ) 表示L矩阵第k次迭代的第j行。
[0136] (四):
[0137] 步骤S3中，终止条件为连续10次误差小于预设值，或者达到预先设定的迭代次数。
[0138] (五)：
[0139] 以下通过具体试验数据解释本发明方法的优越性：
[0140] 该试验使用的真实数据集是Jasper Ridge数据集。数据集原始数据中有512×614个像元，每个像元记录从0.38到2.5μm范围的224个通道上，光谱分辨率高达9.46nm。本文使用其中大小为100×100的子图像进行实验。第一像素从原始图像中的第(105,269)像素开始。由于受到密集水蒸气及大气的影响，高光谱解混实验之前通常需要预处理，去除受影响的波段。因此在实验之前首先移除波段1‑3、108‑112、154‑166和220‑224，保留198个波段用于实验。该数据集中有“树”、“土壤”，“水”和“道路”4个端元，Jasper Ridge数据立方体和4种地物的参考光谱特征如图2所示。
[0141] 为了定量评价本文所提算法的解混性能，本文对端元估计值和参考光谱值之间的SAD值以及估计丰度和真实丰度之间的RMSE进行比较，每组实验重复5次，取结果的平均值进行制表，结果如表1、表2所示。表中最优结果用粗体标注，次优结果用下划线标注。从表中可以得出，本文所提算法获得的各端元SAD及RMSE值都小于TV‑RSNMF方法，且改进算法估计出的大部分端元与参考端元的SAD值为最优，且平均SAD比其他算法都低，这表明改进算法在提取端元方面优于其他算法，这体现引入端元约束的优越性。
[0142] 表1 Jasper Ridge数据各种算法的SAD对比
[0143]
[0144] 表2 Jasper Ridge数据各种算法的RMSE对比
[0145]
[0146]
[0147] 考虑到光谱库中的光谱是在理想条件下获得的，而估计的端元是从受大气及其他环境影响的条件下获得的，所以难免会存在偏差，图3‑6为本试验的结果图，包括4种光谱库中获得的端元光谱和估计的光谱曲线，以及估计的丰度图。从图中可以看出，估计的端元光谱与光谱库中的光谱虽然反射率存在一定差异，但整体趋势大体一致，这是由于真实场景中存在光谱变异性导致的。
[0148] 实验结果表明，整合端元和丰度的约束方法与其他方法相比，无论在端元提取还是在丰度估计精度方面都有一定的提升。
[0149] 以上所述的实施例仅仅是对本发明的优选实施方式进行描述，并非对本发明的范围进行限定，在不脱离本发明设计精神的前提下，本领域普通技术人员对本发明的技术方案作出的各种变形和改进，均应落入本发明的保护范围内。
