αυ〗，· · ·，αρ]，则端元矩阵 A 为：
[0074] A = Qa+d (17)。
[0075] 2)若矩阵A是非负矩阵，计算丰度系数矩阵：
[0076]
[0077] 式中lfv为行向量，且向量中每个元素值都为1，向量中元素个数为N; 1N为列向量， 且向量中每个元素值都为1，向量中元素个数为N; lp-i为列向量，且向量中每个元素值都为 1，向量中元素个数为(P_l);X[i]为第i个像素点丰度系数向量，f已由步骤一计算得出。
【具体实施方式】 [0078] 二:本实施方式中，用于解混实验的高光谱图像遥感数据集由标准 光谱库数据混合得到，混合光谱数据中端元个数P = 3,每个像素的丰度值满足和为1并且非 负的条件。数据的信噪比为40dB。数据中的端元光谱曲线在USGS光谱库中选择。数据大小为 224X5000,224为波段数，5000为像素点数。该仿真数据集中不存在纯像元，像素的丰度系 数均小于0.8。
[0079]执行步骤一:数据加载与预处理。
[0080] 加载存储原始数据的矩阵imageO e R224x5_。已知数据中端元数p = 3。
[0081] 对数据矩阵imageO进行降维处理，由于端元数为3,因此，通过式（1)得到降至2维 的数据矩阵，保存在imagel中，imagel代表降维后的数据矩阵。
[0082] 执行步骤二:寻找满足初始条件数据。
[0083]将降维后的数据矩阵imagel作为输入矩阵，根据【具体实施方式】一种步骤二第1)步 的方法，在数据矩阵imagel中选择3个列向量，将结果存入矩阵Μ〇, <ξ肢 2x3。
[0084] 对降维后的数据矩阵imagel中所有的列向量另进行如下判断：首先，根据式(3)、 式(4)和式(5)得到矩阵I，将结果存入中间变量矩阵S_0中；然后，再判断S_0中所有数值是 否都大于〇,如果不是，则将免存入数据筛选矩阵image2中；最后，得到的数据矩筛选阵 image2e ]R2xm，m < 5000，并且在很多情况下m< <5000，对于本实施方式m = 2114。
[0085] 通过式(8)得到数据矩阵q。以中间变量矩阵Η与g的元素作为变量，根据线性不等 式(6)和式(7)，得到一组可行解，将可行解存入中间变量初始化矩阵H_0与g_0中，作为后续 迭代初始值。
[0086] 执行步骤三:循环求解双线性规划子问题，选择性更新最优值。
[0087] 若是由步骤二转入步骤三，则中间变量迭代初值矩阵!1_〇48丨11与8_ 〇48丨11有!1_ origin = H-0，g_origin = g_0〇
[0088] 令中间变量矩阵H=H_origin，g = g_origin，同时采用如下方法更新数据矩阵Η与 g的第i行的所有元素：
[0089]首先采用内点对偶法求解线性规划问题式（14)，得到矩阵行向量hp_new，元素 gp_ new为此线性规划问题的最优可行解，以及在这组可行解下的目标函数值?_抑1;然后采用 同样的方法求解线性规划问题式（15)，得到最优可行解矩阵向量hq_new，元素 gq_neW，以及 在这组可行解下的目标函数值9_^1。如果满足| p_val | > | q_val |的条件，贝lj用hp_new、gp_ new分别代替数据矩阵Η与g的第i行所有元素;否则用hq_new、gq_new分别代替数据矩阵Η与 g的第i行所有元素。
[0090] 采用上述更新Η与g的第i行所有元素的方式，将中间变量矩阵Η与g从第1行到第p-1(B卩3-1 = 2)行的所有元素都替换一遍。进而完成Η与g的全局更新。
[0091] 执行步骤四：根据变化率检测终止条件判断是否终止迭代计算，若不满足终止条 件，则返回步骤三，继续更新中间变量矩阵Η与g。
[0092] 首先，设置比较阈值th = ε = 〇 · 05;其次，计算数值e = ( I det(H)卜I det(H_origin) )/ I det(H_origin) I ;最后，若θ<ε，则令中间变量终值矩阵!1_?;[1^1=!1，8_?;[1^1 = 8;否则 令中间变量迭代初值矩阵H_origin = H，g_origin = g，并转回执行步骤三。
[0093] 执行步骤五:还原端元矩阵，计算丰度系数，完成图像的解混。
[0094] 计算矩阵 alpha_p = (H_Final)-1 Xg_Final ;计算矩阵 alpha_pi = alpha_pX 1( +(H___Final) 1 ;计算矩阵alpha= [alpha_pl alpha_p];根据式（17) 计算端元矩阵A;端元矩阵A应满足非负性条件，按照式（18)计算所有像素点的丰度系数，存 入数据矩阵X中。得到丰度系数矩阵Z e瓶,完成高光谱的解混。
[0095]图2为本实施例数据点集、真实端元集、估计端元集在p-1(即3-1 = 2)维上的投影 显示，图3为真实端元光谱图与估计端元光谱图在224个波段上的显示。由图2可以看出估计 端元与真实端元的二维表示几乎重合，图3显示估计的端元曲线与真实的端元曲线也几乎 重合，这些都充分显示出了该解混算法的有效性及高精度筛选特性。
【主权项】
1. 一种基于最小体积与优化约束条件的高光谱解混方法，其特征在于所述方法步骤如 下： 步骤一、数据加载与预处理： 1) 加载数据矩阵y。E ，其中L是波段数，N是图像的像素点数量； 2) 对数据矩阵Fe化b'v进行降维，将数据矩阵Y降维至P-I维数据矩阵胶(p IW，p为 端元个数； 步骤二、筛选图像采样点，构造优化的约束条件，寻找满足初始条件的数据； 步骤=、将非负非线性规划问题转换为线性规划问题，对线性规划问题的目标函数结 合优化后的约束条件进行求解，计算中间变量矩阵Hnew、gnew ; 步骤四、根据变化率检测终止条件判断是否终止迭代计算，若不满足终止条件，则返回 步骤S，继续更新中间变量矩阵Hnew、gnew ; 步骤五、解出满足非负性要求的端元矩阵，并计算丰度系数，完成图像的解混。2. 根据权利要求1所述的基于最小体积与优化约束条件的高光谱解混方法，其特征在 于所述步骤二的具体步骤如下： 1 )在数据矩阵F中得到P个列向量： 馬=识，敬…，对 ， 式中，始。￡吸{P-心; 2) 在数据集f中寻找数据点集，数据集中任意一点满足不在由Mo中的列向量作为顶点 构成的单形体中的条件，最后得到满足要求的数据点集f; 3) 优化的线性不等式约束式： 、['，、!! 、浏古y ; q =皆L巧WuJL矿 ; 其中，qERWp-i;，// e ，g E化c'-iw，H与g为要优化的变量矩阵. 4) 求解出满足线性不等式约束式的可行解矩阵化与go,作为后续步骤的初始值。3. 根据权利要求2所述的基于最小体积与优化约束条件的高光谱解混方法，其特征在 于所述步骤2)中，具体寻找方法如下： 首先令：其中，度e吸tp-iw'-n，：来eF, 5€吸(的^1，： SG^Pxi，si为S中的第i行元素 . 对所有A € ；f，1 < i < N求取其对应的向量《，若向量互不满足O、S V, 1的条件，则将夾加 入数据集合'iLd中；假设最终有m个不满足上述条件的5?，则最后得到的数据集合 &,。过_ G 化iix*,p<m<N，且m<<N。4. 根据权利要求1所述的基于最小体积与优化约束条件的高光谱解混方法，其特征在 于所述步骤=的具体步骤如下： 1) 构造非线性目标函数： /识，max det(巧） ; 巧居 2) 将非线性目标函数转化为线性的目标函数： 首先将H的行列式det化)通过代数余子式的方式展开： 化1(//) = ?--"; 7?, d州; 式中hu是矩阵H的第i行第j列的元素，出J是将矩阵H移除第i行元素与第j列元素后得到 的子阵，进而目标函数转化为： /(//，.g) = max V(-ir'7;,,dct(//") ; A, '& .,=1 将上述目标函数分解为下面两个线性目标函数： 户'(/?,诚）二隱X 艺(-1)' i '7?,,化1'(/'/,,) ， 吟混j=i ' (/'化，公,)=min 打-1)",.々/,化1(八"） ； hi毎M 3) 将线性目标函数与优化后的线性约束条件结合，构造如下双线性规划子问题：其中，h功矩阵H第i行的行向量，gi为向量g的第i个元素； 4) 求解上述两个线性规划子问题，根据P*与q*的值，选择性地更新对应的矩阵元素 hi与 gi，直到将矩阵H与g中所有的元素都更新一遍，最终得到矩阵H与g中所有元素都更新后的 新矩阵H new^糾 ew O5. 根据权利要求1所述的基于最小体积与优化约束条件的高光谱解混方法，其特征在 于所述步骤四的具体步骤如下： 1) 设置变化率检测终止条件如下： { I det化new) I -1 det化0) I }/ I det化）I <e， 式中，化代表在步骤=开始迭代之前，矩阵H的初值； 2) 判断中间变量矩阵Hnew是否满足终止条件，如果满足终止条件，则令护=Hnew, gnew,H与g作为矩阵H与g的巧终更新结果,否贝!J,令H = Hnew,化二Hnew,g二gnew,转入步骤二。6.根据权利要求1所述的基于最小体积与优化约束条件的高光谱解混方法，其特征在 于所述步骤五的具体步骤如下： 1) 令日P=化*)-ig*，[巧，...，ap-i] = aj:, I +(/'''）I，日=[日 1，日2,...，QpL则端元矩阵A为： A = Qa+d; 2) 若矩阵A是非负矩阵，计算丰度系数矩阵为：式中，這为行向量，且向量中每个元素值都为1，向量中元素个数为N; In为列向量，且向 量中每个元素值都为1，向量中元素个数为N; Ip-I为列向量，且向量中每个元素值都为1，向 量中元素个数为(P-I)。
【专利摘要】本发明公开了一种基于最小体积与优化约束条件的高光谱解混方法，其步骤如下：步骤一：数据加载与预处理。步骤二：筛选图像采样点，构造优化的约束条件，寻找满足初始条件的数据。步骤三：将非负非线性规划问题转换为线性规划问题，对线性规划问题的目标函数结合优化后的约束条件进行求解，计算中间变量矩阵Hnew、gnew。步骤四：根据变化率检测终止条件判断是否终止迭代计算，若不满足终止条件，则返回步骤三，继续更新中间变量矩阵Hnew、gnew。步骤五：解出满足非负性要求的端元矩阵，并计算丰度系数，完成图像的解混。本方法解决基于最小体积的高光谱解混算法中，约束条件过多所导致的数据存储空间大，运算时间长，精确性难以提高的问题。
【IPC分类】G06T7/40
【公开号】CN105513097
【申请号】CN201511022789
【发明人】张淼, 王天成, 郭威, 沈毅
【申请人】哈尔滨工业大学
【公开日】2016年4月20日
【申请日】2015年12月30日