一种基于最小体积与优化约束条件的高光谱解混方法
【技术领域】
[0001] 本发明涉及一种高光谱数据解混方法，具体涉及一种基于几何学最小体积算法解 混的方法。
【背景技术】
[0002] 高光谱遥感能利用很窄的电磁波波段从感兴趣的物体上获取有关数据。高光谱遥 感是在电磁波谱的可见光、近红外、中红外和热红外波段范围内，获取许多非常窄的连续光 谱数据的技术，目前先进的成像光谱仪可以收集到成百上千个非常窄的光谱波段信息。
[0003] 利用高光谱解混技术能够确定高光谱图像中的端元数、端元光谱曲线以及每个像 素点的丰度值，进而能够确定某种物质在图像观察区域的分布。在解混方法中，高光谱的混 合模型分为两类，分别为线性混合模型以及非线性混合模型。线性模型假设物质成份在空 间布局上是离散混合，忽略在不同类型物质之间的多重散射量，即每个像元点的光谱幅度 仅是这个像元中存在的各物质成份光谱信号，以相对贡献量多少的线性叠加。对于高光谱 线性混合模型来说，解混方法可以分为三大类，分别为：1)基于凸面几何学方法的线性解混 模型；2)基于数据统计的方法，代表算法有ICA等，其核心思想在于假设原信号在体积上是 独立分布的；3)基于稀疏模型的解混方法，代表算法有BP、BTON等。
[0004] 基于几何学的解混算法又能够被分为两类:一类是基于纯像元假设的解混算法， 另一类是基于最小体积法的解混算法。基于纯像元假设的解混算法需要保证图像中的每个 不同的端元至少完全占据一个像素，否则会出现端元错误估计的情况。而基于最小体积的 解混算法即使在缺失纯像元的情况下也能够较高地估计出端元，但缺点是其优化目标函数 是非线性且非凸的，很多情况下，只能得到近似的局部最优解，同时其约束条件数量非常 多，使其算法效率没有基于纯像元假设的解混算法高。
【发明内容】
[0005] 为了解决基于最小体积的高光谱解混算法中，约束条件过多所导致的数据存储空 间大，运算时间长，精确性难以提高的问题，本发明提供了一种基于最小体积与优化约束条 件的高光谱解混方法，通过优化约束条件改进基于最小体积法的高光谱解混方法，得出更 准确的解混效果。
[0006] 本发明的目的是通过以下技术方案实现的：
[0007] -种基于最小体积与优化约束条件的高光谱解混方法，对高光谱图像进行线性解 混以求取端元矩阵时，考虑初值对求解端元矩阵的影响，并将求解体积最小问题的非线性 问题转化为线性规划问题，近似求解得到非线性优化的局部最小解;并考虑求解端元的非 负物理特性，从而提高求取端元的准确率。具体包括如下步骤：
[0008] 步骤一:数据加载与预处理。
[0009]步骤二:筛选图像采样点，构造优化的约束条件，寻找满足初始条件的数据。
[0010]步骤三:将非负非线性规划问题转换为线性规划问题，对线性规划问题的目标函 数结合优化后的约束条件进行求解，计算中间变量矩阵Hne3W、gn?。
[0011] 步骤四：根据变化率检测终止条件判断是否终止迭代计算，若不满足终止条件，则 返回步骤三，继续更新中间变量矩阵11_、g ne3W。
[0012] 步骤五:解出满足非负性要求的端元矩阵，并计算丰度系数，完成图像的解混。
[0013] 本发明与现有技术相比具有如下优点：
[0014] 1)通过筛选数据采样点，同时利用某一采样点必然能包含在由其他顶点构成的单 形体内的性质，有效缩减了约束条件的个数使得计算开销变小，同时解混效率有较大提高。
[0015] 2)不论是否存在纯像元，本方法都能够经过解算，得到高光谱图像的端元，进而得 到丰度系数，完成高光谱图像的解混。
【附图说明】
[0016] 图1为本发明的流程图；
[0017] 图2为数据点集、真实端元集、估计端元集在2(P-1)维上的投影显示；
[0018]图3为真实端元光谱图与估计端元光谱图在224个波段上的显示。
【具体实施方式】
[0019] 下面结合附图对本发明的技术方案作进一步的说明，但并不局限于此，凡是对本 发明技术方案进行修改或者等同替换，而不脱离本发明技术方案的精神和范围，均应涵盖 在本发明的保护范围中。
【具体实施方式】 [0020] 一:本实施方式提供了一种基于最小体积与优化约束条件的高光谱 解混方法，假设光谱图像数据FVR&Ar，L是波段数，N是图像的像素点数量。若在光谱图像 中有P个端元，首先对光谱信息进行降维处理，将数据从L维映射到p-1维;接着构造一个单 形体，该单形体能够包含所有降维后的数据，通过循环求解线性规划子问题，不断更新单形 体顶点矩阵，最终使单形体的体积收敛至稳定，最后通过单形体的顶点矩阵得到图像的端 元矩阵，进而求出图像的丰度系数，完成高光谱图像的解混任务。
[0021 ]如图1所示，共分为五个步骤，具体步骤如下：
[0022]步骤一:数据加载与预处理。
[0023] 1)加载数据矩阵Γ0 ，其中L是波段数，N是图像的像素点数量；
[0024] 2)对数据矩阵y 行降维。已知光谱图像中，端元个数为p，将数据矩阵Y降 维至p-Ι维数据矩阵f e &具体降维方法为：
[0025] f =QfU (.1) ?
[0026] 其中，U=[yi-d，y2_d，. . .，yN_d]:
.[，数据矩阵Y表示为列向量的形式： Y=[yi,y2, . . . ,Υν] ； O' ={Ql Q) ?1 ，Q= [qi(UUT)，q2(UUT)，· · ·，qP-1(UUT)]，qi(UUT)，i = 1，2, . . .，p_l表示矩阵1]1^第1大的特征值所对应的特征向量。
[0027] 步骤二:采用优化后的约束条件，寻找满足初始条件的数据。
[0028] 1)在数据矩阵f中得到p个列向量：
[0029]
[0030] 式中，ΛΑ E鈀〃1~，向量选择满足：由这p个列向量作为顶点构成的单形体体积比 其他任意P个列向量作为顶点构成的单形体体积大的条件。
[0031] 2)在数据集f中寻找数据点集，数据集中任意一点满足不在由M〇中的列向量作为 顶点构成的单形体中的条件，最后得到满足要求的数据点集r。
[0032]具体寻找方法如下，首先令：
[0033]
[0034]
[0035]
[0036] 其中，B e 股(;，-1)x(p-11，天 e ，si为s中的第i行元素。
[0037] 对所有兄e f, 1 ?? / ?? ;V (N为像素点个数）求取其对应的向量i，若向量i不满足 0_.·< ? <1的条件，则将另加入数据集合.fWwA_.中。假设最终有m(显然p < m < N，且一般情况下m <<N)个不满足上述条件的免，则最后得到的数据集合
[0038] 3)优化的线性不等式约束条件为式(6)与式(7):
[0039]
[0040]
[0041]
[0042] 其中，qERlx(P-D，e W，： g eR("-11x1，Η与g为要优化的变量矩阵。式(6) 与式(7)将不等式约束个数降为了 (p_l)X(m+l)个，不等式约束的个数被大大减小，从而大 大提高了求解的速度。
[0043] 本步骤通过求解线性规划问题的方法，求解出一个满足线性不等式约束式(6)与 式(7)关于变量矩阵Η与g的一组可行解矩阵Ho与go,作为后续步骤的初始值。
[0044] 本步骤中的<、<、2分别表示所有元素都<、<、2。
[0045] 步骤三:对线性规划问题的目标函数结合优化后的约束条件进行求解，计算中间 变量矩阵Hne3W、g n?。迭代计算中，变量矩阵Η与g的阶段性最优变量值存储在矩阵Hne3W、gne3W中。
[0046] 构造优化函数，通过代数余子式的展开，将非线性的优化问题改造成为线性规划 子问题，通过周期循环迭代的方式，更新得到新的矩阵Hne3W、gne3W，直到满足停止迭代的条件。 [0047] 1)构造非线性目标函数：
[0048]
[0049] 2)将非线性目标函数转化为线性的目标函数：
[0050] 首先将Η的行列式det(H)通过代数余子式的方式展开：
[0051]
[0052] 式中是矩阵Η的第i行第j列的元素，是将矩阵Η移除第i行元素与第j列元素 后得到的子阵。进而目标函数转化为：
[0053]
.
[0054]由于目标优化函数| det(H) |仍然是非线性的，因此将其分解为下面两个线性目标 函数：
[0055]
[0056]
[0057] 3)将线性目标函数与优化后的线性约束条件结合，构造如下双线性规划子问题：
[0058]
[0059]
[0060]
[0061] .V/./7； r；i[/]-g, >0,'^.1',[/]￡ ΚΛ,,ι；ι,1；; ,// -1；,
[0062] 其中，lu为矩阵Η第i行的行向量，gl为向量g的第i个元素。
[0063] 4)求解上述两个线性规划子问题，根据f与cf的值，选择性地更新对应的矩阵元素 h^gl，直到将矩阵Η与g中所有的元素都更新一遍。
[0064] 具体更新hggl的规则为:若求解线性规划问题式(14)得到的解为ξ、瓦，求解线 性规划问题式（15)得到的解为，若|p1 > |q$|，则将hggl更新为[与瓦，否则，将lu与 gi更新为t与gi。
[0065] 在这一步骤中，逐行更新Η与g的值:令i = l, ...，p-l，依次得到hi与gi的值，最终得 到矩阵Η与g中所有元素都更新后的新矩阵!1_与gne3W。
[0066] 步骤四：根据变化率检测终止条件判断是否终止迭代计算，若不满足终止条件，则 返回步骤三，继续更新中间变量矩阵11_，g ne3W。
[0067] 1)设置变化率检测终止条件如下：
[0068] { |det(Hnew) |-|det(Ho) | }/|det(H) | <ε (16)。
[0069] 式中，Ho代表在步骤三开始迭代之前，矩阵Η的初值。
[0070] 2)判断中间变量矩阵Hnew是否满足条件式（16)，如果满足终止条件，则令!f = Hnew， g* = gnew，Η*与g*作为矩阵Η与g的最终更新结果。否贝lj，令H=Hnew，Ho = Hnew，g = gnew，转入步骤 三的第4)步。
[0071 ]步骤五:还原端元矩阵，计算丰度系数，完成图像的解混。
[0072] 根据f与g'还原端元矩阵，并计算丰度系数，完成图像的解混，端元矩阵应满足非 负性的要求。
[0073] 1)令αρ = (Η*)-4*，[?丨,…丨]丨十（/广）丨，α = [