基于自适应局部邻域约束的小波域NMF图像解混方法

基于自适应局部邻域约束的小波域nmf图像解混方法
技术领域e
1.本发明涉及高光谱解混的技术领域，尤其涉及一种基于自适应局部邻域约束的小波域nmf图像解混方法。


背景技术：

2.近年来，高光谱遥感成像技术逐渐成为学者们研究的热点。高光谱遥感图像的处理与分析仍是制约高光谱遥感应用的重要瓶颈，而高光谱遥感图像解混的研究则是伴随高光谱成像技术发展的新型数据处理问题，是高光谱数据处理的重要任务之一。
3.高光谱遥感影像信息空间分辨率的限制以及复杂地物的影响，图像中出现大量混合像元，即一个像元同时包含几个地物的光谱信息。若一个像元只包含一种地物，则该像元称之为端元。高光谱解混主要包括端元提取和丰度估计，即确定组成混合像元的地物种类以及确定各个种类所占的比例。
4.现有的高光谱遥感影像丰度估计的方法主要是采用全约束的最小二乘方法，但是此方法在处理含有噪声的数据时，丰度估计的误差较大，解混的精度会降低。或者是采用非负矩阵分解(non-negativematrix factorization,nmf)方法。为避免nmf算法陷入局部最优解，学者们提出了各类针对nmf的约束，其中利用空间信息的局部邻域约束在近期取得了不错的解混效果，然而该算法是在时域中进行，但由于存在病态矩阵问题，其解混的精度往往会受到影响。


技术实现要素：

5.本发明提供一种基于自适应局部邻域约束的小波域nmf图像解混方法，用以解决现有的时域中nmf约束引入的病态矩阵而导致解混精度不高的问题。
6.本发明提供一种基于自适应局部邻域约束的小波域nmf图像解混方法，包括：
7.步骤一、获取高光谱遥感图像；
8.步骤二、获取高光谱数据矩阵x、端元矩阵a和丰度矩阵s；
9.步骤三、分别对所述高光谱数据矩阵和所述端元数据矩阵进行d层双正交小波包变换，分别得到高光谱数据矩阵对应的小波包树、得到端元数据矩阵对应的小波包树；
10.步骤四、获取所述端元数据矩阵对应的小波包树每个节点处的系数矩阵，计算每个所述系数矩阵的条件数，找到最低条件数的节点，记为min cond；
11.步骤五、将所述最低条件数的节点min cond对应的高光谱数据矩阵的系数矩阵、该节点对应的端元数据矩阵的系数矩阵，以及丰度矩阵，作为基于自适应局部邻域加权约束的nmf模型输入，迭代更新丰度矩阵和端元矩阵，直到所述模型收敛，输出收敛时迭代更新后的丰度矩阵和端元矩阵；
12.步骤六、将与mincond节点处在同一层的其他节点的系数矩阵，以及所述迭代更新后的端元矩阵作为重构函数的输入，输出得到重构后的端元矩阵；
13.其中在所述步骤五中，基于自适应局部邻域加权约束的nmf模型为
[0014][0015][0016]
其中λ是一个可调参数用于平衡近似误差和约束项，w
ij
为给定像元xi的局部邻域像元xj对像元xi权重的贡献，xi表示高光谱数据矩阵的第i列，xj表示高光谱数据矩阵的第j列，n(i)表示定像元xi的局部邻域像元的列集合，j∈n(i)；si表示的丰度矩阵s的第i列，sj表示的丰度矩阵s的第j列，p是端元数量，n是像元的总数；
[0017]
其中在所述步骤五中，迭代更新丰度矩阵s和端元矩阵a，
[0018]a←a·
*(xs
t
)./(ass
t
)
[0019][0020]
其中
·
*表示矩阵对应元素相乘，./表示矩阵对应元素相除，矩阵χi表示xi的局部邻域权重之和，sw的第i列表示为xf和af分别是x和a的增广矩阵。
[0021]
与现有的技术比较，本发明的有益效果是：
[0022]
提出了一种自适应局部邻域约束的nmf模型，根据丰度矩阵的特点可以自适应确定给定像元的局部邻域,计算出的权重充分地利用了给定像元和邻域内像元的空间信息，使得解混精度更高。
[0023]
进而，不同于传统的时域nmf的解混图像方法，将自适应局部邻域约束的nmf模型与双正交小波包变换相结合，获得基于自适应局部邻域约束的小波域nmf图像解混方法。采用双正交小波作为小波基对数据进行表示,小波基线性相位和较短支撑集的特性，可以更好地对数据进行压缩和消减噪声对于解混结果的影响，进一步克服了在解混过程中出现的病态矩阵问题，进一步提升了解混精度效果；
[0024]
与此同时，计算选取最小条件数节点进行后续计算和重构，有效的降低了数据量处理，以及图像信息得到了较为完整的保留。
附图说明
[0025]
为了更清楚地说明本发明或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作一简单地介绍，显而易见地，下面描述中的附图是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
[0026]
图1为本发明提供的一种基于自适应局部邻域约束的小波域nmf图像解混方法的流程示意图；
[0027]
图2为本发明提供的局部区域划分图；
具体实施方式
[0028]
为使本发明的目的、技术方案和优点更加清楚，下面将结合本发明中的附图，对本发明中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0029]
现有的手势识别分类技术中普遍存在多维手势特征提取不够有效便捷，手势分类的准确率不高的问题。下面结合图1描述本发明的基于小波域的高光谱遥感图像自适应局部邻域约束的nmf解混方法。图1为本发明的基于自适应局部邻域约束的小波域nmf图像解混方法的流程示意图，如图1所示，该方法包括：
[0030]
步骤一、获取高光谱遥感图像；
[0031]
具体地，获取大小为n
×
m的高光谱遥感图像x。
[0032]
步骤二、获取高光谱数据矩阵、端元矩阵和丰度矩阵；
[0033]
具体地，使用vca算法和fcls算法对高光谱遥感图像初始化，得到端元矩阵a_init和丰度矩阵s_init。初始化的目的是在开始使用nmf算法解混之前,先将端元矩阵a和丰度矩阵s限制在一定的范围内，这样在进行迭代时,可以加快算法收敛，减少运行时间，提升算法的精度。
[0034]
nmf是在矩阵所有元素非负的条件下对一大型矩阵实现非负分解，得到两个非负低秩矩阵，其数学模型为:
[0035]
v＝wh
ꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(1)
[0036]
式中：v为原始数据；w为基矩阵；h为系数矩阵.且有w∈rn×r，h∈rr×m，其中m，n，r均为代表矩阵大小的常数，v，w，h均为非负矩阵。
[0037]
如果不考虑噪声的影响，那么线性混合模型和nmf分解模型形式一样.将遥感图像x表示nmf中的矩阵v，端元光谱矩阵a表示nmf分解矩阵中的w，则丰度矩阵s可以用nmf中的h表示,同时对h也加上“和为一”约束。这样利用nmf进行遥感图像解混可以同时求取出混合像元中存在的端元和丰度。
[0038]
x＝as+e
ꢀꢀ
(2)
[0039]
步骤三、分别对所述高光谱数据矩阵和所述端元数据矩阵进行d层双正交小波包变换，分别得到高光谱数据矩阵对应的小波包树、得到端元数据矩阵对应的小波包树；
[0040]
高光谱图像数据在一个合适的变换域内是可压缩的，如傅里叶、小波等。傅里叶变换具有良好的频率局部化特性，但空间局部化较差。为了在空间局部化和频率局部化之间达成折衷，小波变换是一种方便的多分辨率分析工具。此外，小波基比傅里叶某具有更好的压缩性。另一方面，小波变换更适合于数据压缩相关应用，小波变换通过更好地约束解空间，有助于解决盲谱分解的不适定性问题。因此，我们可以使用小波基并利用其紧凑的表示能力来表示高光谱数据。此外，为了解决对称性和信号重构的不相容性，引入了biorthohonal(双正交小波)，双正交小波解决了线性相位和正交性要求的矛盾。
[0041]
小波包分析能够为信号提供一种更精细的分析方法，它将频带进行多层次划分，对多分辨率分析中没有细分的高频部分进一步分解，并能够根据被分析信号的特征，自适应地选择相应频段，使之与信号频谱相匹配，从而提高了时频分辨率。
[0042]
小波包变换是在小波变换的基础上发展而来的，它克服了小波变换在高频部分频
域分辨率较低的缺点。任给一个时域信号x(t)，经j层小波包分解后可以写成：
[0043][0044]
其中：
[0045][0046]
小波包系数为：
[0047][0048]
式中：ψ
j,k,i
(t)为具有尺度指标j、位置指标k和频率指标i的小波包。因为ψ
j,k,i
(t)为一组标准正交基，故当m≠n时，有：
[0049][0050]
信号x(t)的总能量为：
[0051][0052]
由小波包的正交性可得：
[0053][0054]
其中：
[0055][0056]
式中：为第i频段内的小波包节点能量，而信号的小波包节点总能量等于各频段的信号节点能量之和。
[0057]
双正交小波包变换使用的小波包基具有正交性，能够保留高光谱数据的线性结构，将高光谱数据移至小波域中，高光谱数据x的小波包变换是一个内积，定义为：
[0058]
xw:＝《x,ψ
a,b
》
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(10)
[0059]
式中，x是高光谱数据，xw是wpt变换后的高光谱数据，ψ
a,b
为小波包基，《》为内积算子。由式(10)可以得到迭代规则在小波域中的变换：
[0060][0061]
式中，xw表示大小为k
×
n的高光谱数据矩阵，ew表示大小为k
×
n的高斯噪声矩阵,对于e个端元成员，小波域内的端元成员矩阵aw大小为k
×
e，si表示大小为e
×
n的丰度矩阵。其中，k是的小波包分解后的数据维度，具体的数值大小由分解的节点和小波包基的选择来确定。
[0062]
将小波域中的lmm模型应用到nmf算法中，实现小波变换和nmf算法的结合。通常欧氏距离的平方表达比较直观且应用得更加广泛，在此基础上，小波域中的nmf的最优化问题可以表述为：
[0063][0064]
式中，s表示丰度矩阵。
[0065]
具体地，在本实施例中，分别对x和a_init进行d层双正交小波包变换，得到x的小波包树和a_init的小波包树。小包分解的层数d决定了光谱波段压缩的程度，d的取值越大,则光谱波段压缩的程度越高。但是在实际应用中,分解层数每增加一层，光谱波段会被压缩为原来的1/2左右(具体取值受小波基的选取的影响)，故在最d进行取值时，要考虑原始数据波段的大小。如果d的值过小，则压缩不充分，易造成数据冗余问题，无法良好的解决病态矩阵问题；如果d的值过大,则造成光谱数据压缩过于严重，无法以低维数据有效表达原始高维数据，造成解混精度的下降。
[0066]
步骤四、获取所述端元数据矩阵对应的小波包树每个节点处的系数矩阵，计算每个所述系数矩阵的条件数，找到最低条件数的节点，记为min cond；
[0067]
对于线性模型x＝as而言，矩阵a的条件数用于衡量矩阵乘法或逆的输出对输入误差的敏感性，条件数越大表明敏感性越差。对于一般的矩阵a，存在a的逆矩阵a-1
(即a为非奇异矩阵)，则其条件数为：
[0068]
cond(a)＝‖a‖
·
‖a-1
‖
ꢀꢀꢀꢀꢀꢀꢀꢀ
(13)
[0069]
式中，cond(a)表示矩阵a的条件数，‖
·
‖表示矩阵的范数。根据定义来看，任意一种矩阵范数都可以用来定义条件数。则对应的2-范数的条件数为：
[0070]
cond2(a)＝‖a‖2·
‖a-1
‖2ꢀꢀ
(14)
[0071]
又已证明||a||2＝σ
max
，||a-1
||2＝1/σ
min
，故2-范数的条件数可以通过矩阵的奇异值分解计算得到：
[0072][0073]
式中，cond2(a)表示矩阵a的二范数条件数，σ
max
表示矩阵的最大奇异值，σ
min
表示矩阵的最小奇异值。
[0074]
在同质区域像元之间都存在一定的相似性，但是考虑到像元间相似性时，考虑到算法的计算量，不便将中心像元与其他所有像元进行相似性判断。所以要借助局部邻域的思想解决此问题，在考虑像元xi相似性时，仅考虑xi的局部邻域n(i)内的像元，xi表示高光谱数据矩阵y∈r
l
×n的第i列。
[0075]
局部邻域n(i)的选取方式有多种，常见的局部邻域结构有方形与圆形结构，这些结构操作起来比较容易，在处理过渡区域、边界区域时不能很好地利用空间信息。结合丰度数据的特点，该算法提出了一种可以充分利用空间信息的局部邻域确定算法。
[0076]
该方法依据的原理为是将整幅高光谱图像分为多个不同的相似区域，每个相似区域都比较接近于一种地物，即一种端元。图2为高光谱图像一个局部区域的划分图。假设在区域1的地物比较接近于端元k，那么端元k在1区域中的占比就比较大，而其他端元在1区域中的丰度值就比较小。因此，通过比较不同丰度图在同一位置的丰度值即可判断出与该位置像元相近的端元。随后，利用相似区域内对应端元的丰度值比较接近的点确定局部邻域。
[0077]
首先，确定每一个像元的空间坐标。假设高光谱图像排列在一个r
×
c的网格上，并有r
×
c＝n，网格上的每一个位置都对应一个l维的像素向量，即一个像元。将高光谱数据矩阵y∈r
l
×n的n个像元按照从左到右，从上到下的顺序排列到网格当中。这样，就可以得到像元yi和它空间位置(p,q)的关系
[0078]
i＝(p-1)*r+q
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(16)
[0079]
i＝(q-1)*c+p
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(17)
[0080]
其次，改变丰度数据的表示方式，以获得像元间的相似性。让丰度矩阵s的每行，按照其像元的空间位置整理成对应的矩阵。以其中任意一个端元的丰度sq∈r1×n为例，设其被整理为sq∈rr×c，其中r和c分别表示矩阵的行数和列数。p个端元对应的p行丰度分别被整理成中的元素表示为元素值的大小等于丰度值。
[0081]
然后，在确定空间坐标为(i,j)的给定像元的局部邻域前，需要先判断该像元接近于哪种端元。通过比较空间坐标为(i,j)的各丰度矩阵对应的丰度值的大小，找到最大的即可判断该像元比较接近于端元k。
[0082]
最后，根据相似区域对应丰度值比较接近的特点，在端元k对应的丰度矩阵中，以为中心,3
×
3窗口内找到所有丰度值在范围内的位置作为局部邻域n(i)，τ是一个可控的参数。遍历每一个像元即可得到每一个像元的局部邻域区域。
[0083]
为了确定局部邻域像元xj(j∈n(i))对给定像元xj的贡献，要考虑两个方面的因素：空间距离和局部邻域像元的相似性。因为局部邻域内的像元对应的丰度相当接近，所以可以通过丰度的相似性来判断像元的相似性。
[0084]
基于上述的分析，局部邻域像元xj对xi权重的贡献可以定义为：
[0085][0086]
其中，α表示像元之间的空间距离,β表示像元之间的相似性,和β的取值范围都在(0,1]，都是同一个量级。假设第i个像元和第j个像元的空间坐标分别是(r,s)，(k,l)，我们简单定义α和β为：
[0087]
α＝|r-k|+|s-l|
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(19)
[0088]
β＝《si·
sj》
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(20)
[0089]
其中|
·
|表示绝对值，《
·
》表示内积，α表示像元xi和xj之间的曼哈顿距离，β表示丰度si和sj的相似性，si表示的丰度矩阵s的第i列。当距离比较近，丰度相似度比较高时，第j个像元对第i个像元的影响就越大。显然，w
ij
可以衡量第j个像元在空间距离和相似性方面对第i个像元的贡献。
[0090]
对于给定像元xi，局部邻域内的全部像元对它的影响可以表示为χi,χi可以通过计算每个邻域像元xj(j∈n(i))对xi贡献的和来表式:
[0091]
χi＝∑
j∈n(i)wij
,i＝1,2,
…
,n
ꢀꢀꢀꢀꢀꢀꢀ
(21)
[0092]
具体地，在本实施例中，在a的小波包树中可以得到(2
d+1-2)个节点，计算每个节点系数矩阵的2-范数条件数，找到最低条件数的节点，记为min cond。
[0093]
步骤五、将所述最低条件数的节点min cond对应的高光谱数据矩阵的系数矩阵、该节点对应的端元数据矩阵的系数矩阵，以及丰度矩阵，作为基于自适应局部邻域加权约束的nmf模型输入，迭代更新丰度矩阵和端元矩阵，直到所述模型收敛，输出收敛时迭代更新后的丰度矩阵和端元矩阵；
[0094]
因为nmf模型的目标函数的非凸性，把nmf算法直接用于高光谱图像解混会有局部
极小值问题。要更好的对端元a和丰度s进行估计，就得在nmf的模型中加入些许辅助约束。根据给定像元的丰度与局部邻域内像元的丰度相似的特点，将局部邻域丰度的先验信息作为约束项，并对其进行加权操作，充分地利用了高光谱的空间信息和局部邻域像元的相似性，提高了解混的精度。
[0095]
基于自适应局部邻域加权约束的nmf模型为
[0096][0097]
其中λ是一个可调参数用于平衡近似误差和约束项，w
ij
为给定像元xi的局部邻域像元xj对像元xi权重的贡献，xi表示高光谱数据矩阵的第i列，xj表示高光谱数据矩阵的第j列，n(i)表示定像元xi的局部邻域像元的列集合，j∈n(i)；si表示的丰度矩阵s的第i列，sj表示的丰度矩阵s的第j列，p是端元数量，n是像元的总数；
[0098]
经过推导，得到端元a和s的更新规则：
[0099]a←a·
*(xs
t
)./(ass
t
)
ꢀꢀꢀꢀꢀꢀꢀ
(23)
[0100][0101]
其中.*表示矩阵对应元素相乘，./表示矩阵对应元素相除。矩阵χi表示xi的局部邻域权重之和，sw的第i列表示为
[0102]
xf和af是增广矩阵，如式(25)：
[0103][0104]
式中δ是一个正常数，用于调整丰度和为1,δ的约束越大，s的每一列的和越接近1。
[0105]
具体地，在本实施例中，将x的小波包树中的min cond节点处的系数矩阵、a_init的小波包树中的min cond节点处的系数矩阵以及初始化的丰度矩阵s_init，作为基于自适应局部邻域约束的nmf解混算法的输入，迭代求解得到丰度矩阵s和小波域中的端元矩阵aw。
[0106]
步骤六、将与mincond节点处在同一层的其他节点的系数矩阵，以及所述迭代更新后的端元矩阵作为重构函数的输入，输出得到重构后的端元矩阵，即使用条件数最低一层的小波包树节点对端元进行重构。
[0107]
本发明实施例的关键点在于：1、提出一种自适应局部邻域约束的nmf模型，根据丰度矩阵的特点可以自适应确定给定像元的局部邻域,计算出的权重充分地利用了给定像元和邻域内像元的空间信息，使得解混精度更高；2、其次，不同于传统的时域nmf的解混图像方法，将自适应局部邻域约束的nmf模型与双正交小波包变换相结合，获得小波域的高光谱遥感图像自适应局部邻域约束的nmf解混方法，采用双正交小波作为小波基对数据进行表示,小波基线性相位和较短支撑集的特性，可以更好地对数据进行压缩和消减噪声对于解
混结果的影响，进一步克服了在解混过程中出现的病态矩阵问题，进一步提升了解混精度效果；3、并未选取所有的节点输入自适应局部邻域约束的nmf模型，以及参与重构端元矩阵，而是选取最小条件数节点。该节点的计算选取应用，有效的降低了数据量处理，以及图像信息得到了较为完整的保留，即兼顾了图像信息的完整性，计算量小且高效。
[0108]
最后应说明的是：以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的精神和范围。