基于rank-order局部邻域约束的NMF遥感图像解混方法及系统

基于rank-order局部邻域约束的nmf遥感图像解混方法及系统
技术领域
1.本发明涉及遥感图像处理领域，尤其涉及一种基于rank-order局部邻域约束的nmf遥感图像解混方法及系统。


背景技术：

2.近年来，高光谱遥感成像技术逐渐成为学者们研究的热点。高光谱遥感图像的处理与分析仍是制约高光谱遥感应用的重要瓶颈，而高光谱遥感图像解混的研究则是伴随高光谱成像技术发展的新型数据处理问题，是高光谱数据处理的重要任务之一。
3.目前基于空间信息约束的nmf算法，在计算空间距离时，得到的是中心点与其邻域点间的绝对距离，即局部邻域权重仅与绝对距离有关，没有考虑局部邻域像元间的相互作用。当中心像元与其局部邻域像元所属类型不同时，若继续使用绝对距离进行权重计算，那么计算出的权重与同类像元间的权重相等，这显然是不合理的。
4.上述内容仅用于辅助理解本发明的技术方案，并不代表承认上述内容是现有技术。


技术实现要素：

5.为解决上述技术问题，本发明提供一种基于rank-order局部邻域约束的nmf遥感图像解混方法，包括：
6.s1：获取高光谱数据，通过vca-fcls对高光谱数据进行初始化，获得初始端元矩阵和初始丰度矩阵，将初始端元矩阵作为当前端元矩阵，将初始丰度矩阵作为当前丰度矩阵；
7.s2：通过基于rank-order距离的自适应局部邻域算法对当前端元矩阵和当前丰度矩阵进行更新，获得更新后的端元矩阵和更新后的丰度矩阵；
8.s3：若达到最大迭代次数或者误差小于规定的阈值，则输出更新后的端元矩阵和更新后的丰度矩阵；否则将更新后的端元矩阵作为当前端元矩阵，将更新后的丰度矩阵作为当前丰度矩阵，返回步骤s2。
9.优选的，步骤s2具体为：
10.s21：定义xi为高光谱数据x的第i列，xj为x的第j列，将xi设置为中心像元，将xj设置为局部邻域像元，计算获得中心像元与局部邻域像元间的不对称距离d(xi,xj)，以及局部邻域像元与中心像元间的不对称距离d(xj,xi)；
11.s22：通过d(xi,xj)和d(xj,xi)计算获得中心像元xi与局部邻域像元xj间的rank-order距离υ；
12.s23：计算获得局部邻域像元xj对中心像元xi的距离权重p
ij
，以及局部邻域像元xj对中心像元xi的总权重χi；
13.s24：通过总权重χi计算获得权重矩阵w，通过距离权重p
ij
计算获得变量矩阵sw；
14.s25：对当前端元矩阵进行更新，获得第一增广矩阵xf、第二增广矩阵af和更新后的
端元矩阵；
15.s26：通过权重矩阵w、变量矩阵sw、第一增广矩阵xf和第二增广矩阵af对当前丰度矩阵进行更新，获得更新后的丰度矩阵。
16.优选的，步骤s23中距离权重p
ij
的计算公式为：
[0017][0018]
其中，υ表示rank-order距离，σ表示曲线的曲率调节参数，β表示目标点与局部邻域间丰度的相似性；
[0019]
总权重χi的计算公式为：
[0020][0021]
其中，n(i)表示高光谱数据第i列像元的集合，n表示高光谱数据的像元总数。
[0022]
优选的，步骤s24具体为：
[0023]
权重矩阵w的计算公式为：
[0024][0025]
其中，χi为局部邻域像元xj对中心像元xi的总权重；p表示端元数量；n表示高光谱数据的像元总数；
[0026]
sw的i列表示为的i列表示为
[0027]
sj表示丰度矩阵s的第j行，n(i)表示高光谱数据第i列像元的集合。
[0028]
优选的，步骤s25具体为：
[0029]
当前端元矩阵的更新计算公式为：
[0030]
a'
←
a.*(xs
t
)./(ass
t
)
[0031]
其中，a'表示更新后的端元矩阵，a表示当前端元矩阵，.*表示矩阵对应元素相乘，./表示矩阵元素对应相除，(.)
t
表示矩阵的转置，s表示丰度矩阵；
[0032][0033]
其中，δ是一个正常数，用于调整丰度和为1的约束，δ越大，s的每一列的和越接近1。
[0034]
优选的，步骤s26中更新后的丰度矩阵的计算公式为：
[0035][0036]
其中，s'表示更新后的丰度矩阵，s表示当前丰度矩阵.*表示矩阵对应元素相乘，./表示矩阵元素对应相除，λ表示可调参数。
[0037]
优选的，步骤s3中阈值的计算公式为：
[0038]
[0039][0040]
其中，j(a,s)表示阈值，s表示丰度矩阵，a表示当前端元矩阵，s表示丰度矩阵，si表示丰度矩阵s的第i列，sj表示丰度矩阵s的第j行，λ表示可调参数，n(i)表示高光谱数据第i列像元的集合，n表示高光谱数据的像元总数。
[0041]
一种基于rank-order局部邻域约束的nmf遥感图像解混系统，包括：
[0042]
初始化模块，用于获取高光谱数据，通过vca-fcls对高光谱数据进行初始化，获得初始端元矩阵和初始丰度矩阵，将初始端元矩阵作为当前端元矩阵，将初始丰度矩阵作为当前丰度矩阵；
[0043]
更新模块，用于通过基于rank-order距离的自适应局部邻域算法对当前端元矩阵和当前丰度矩阵进行更新，获得更新后的端元矩阵和更新后的丰度矩阵；
[0044]
输出模块，若达到最大迭代次数或者误差小于规定的阈值，则输出更新后的端元矩阵和更新后的丰度矩阵；否则将更新后的端元矩阵作为当前端元矩阵，将更新后的丰度矩阵作为当前丰度矩阵，返回更新模块。
[0045]
本发明具有以下有益效果：
[0046]
使用rank-order算法计算局部邻域间距离时，可以考虑局部邻域点之间的相互作用，使得同类像元拥有的权重较高，而不同类像元拥有的权重较低，优化了之前算法中的得到权重都相同的不足，同时还具有一定的聚类效果，能够有效的减少在计算空间距离加权时所产生的误差，对解混精度的提升有积极作用。
附图说明
[0047]
图1为局部邻域选取范围；
[0048]
图2为高光谱数据图像简易图；
[0049]
图3为样本点xi和xj共享邻近点信息；
[0050]
图4为本发明方法流程图；
[0051]
图5为本发明方法与其他算法解混后丰度图对比；
[0052]
图6为本发明方法端元提取与各算法光谱曲线对比；
[0053]
本发明目的的实现、功能特点及优点将结合实施例，参照附图做进一步说明。
具体实施方式
[0054]
应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。
[0055]
(1)基于局部邻域约束的nmf算法：
[0056]
高光谱数据在同质区域的像元之间都会存在一定的相似性，但在考虑像元间的相似性时，考虑到算法的计算量，不便将中心像元与其他所有像元进行相似性判断，所以需要借助局部邻域的思想解决此问题，在考虑像元xi的相似性时，仅考虑xi的局部邻域n(i)内的像元，xi表示高光谱数据矩阵x∈r
l
×n的第i列。局部邻域n(i)的选取方式有很多种，常见的局部邻域结构有方形或圆形结构，这些结构操作起来较为简单，但是结构固定，在处理过渡区域、边界区域时不能很好地利用空间信息。结合丰度数据的特点，liu等人提出了一种可以更加充分的利用空间信息的局部邻域确定算法。
[0057]
参考图1为局部邻域的选取范围，首先，需要确定每一个像元在空间上的坐标位
置，方便后续的求解；假设高光谱图像的大小是r
×
c，并有r
×
c＝n，得到像元xi和它空间位置(p,q)的关系：
[0058]
i＝(p-1)*r q
ꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(1)
[0059]
i＝(q-1)*c p
ꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(2)
[0060]
其次，将丰度矩阵s的每一行，按照空间位置分割成每一个小矩阵；假设分割出的一个矩阵为sr∈r1×n，在将这个矩阵展开成像元的形式，如其中r和c分别表示矩阵的行数和列数；以为中心，在3*3窗口内先找到局部邻域集合，如果局部邻域中的丰度值没有在的范围内，那么就可以初步判定这些像元是差异过大的异类像元，本文为了剔除丰度值差异过大的像元，设置了一个可以进行自我调节的参数τ；
[0061]
局部邻域像元xj对xi的权重可以表示为：
[0062][0063]
其中，α表示像元之间的距离也就是像元间的空间距离，在进行空间距离计算时，考虑到计算量较大，故采用较为简单的曼哈顿距离；β表示丰度之间的相似性，内积可以表示向量间的投影关系，所以可以很好的表示两个向量间的相似性关系；假设第i个像元和第j个像元的空间坐标分别表示为(n,m)和(v,b)，在本文中α和β的定义为：
[0064]
α＝|n-m| |v-b|
ꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(4)
[0065]
β＝《si·
sj》
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(5)
[0066]
其中，α表示像元xi和xj之间的曼哈顿距离，|
·
|表示绝对值，《
·
》表示内积，β表示丰度si和sj的相似性,si表示的丰度矩阵s的第i列。每个邻域像元xj(j∈n(i))对xi贡献的和χi表示为：
[0067]
χi＝∑
j∈n(i)wij
,i＝1,2,
…
,n
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(6)
[0068]
(2)rank-order距离算法：
[0069]
rank-order距离算法在2011年被提出；rank-order算法与传统的欧氏距离、曼哈顿距离等绝对距离的算法相比，拥有聚类的功能，在局部邻域间计算距离具有良好的效果；
[0070]
假设当前拥有一张高光谱数据的图像如图2所示；在图像中主要包含树木、土壤、水等地物；局部邻域取中心点3x3的区域，局部邻域点在使用绝对距离对中心点进行加权时，从图中可以发现，中心点与局部邻域点的距离相近，但不属于同一种地物，这会造成加权时的误差。如果使用rank-order距离进行计算，则会在之前步骤的基础上，再以局部邻域点为新的中心点，与其局部邻域点再求一次绝对距离，并将结果集进行排序，计算出原中心点作为局部邻域的邻域点时，其绝对距离是否依然相近；若两次计算的绝对距离都相对较近，则局部邻域点对中心的权重较高；距离较远，则权重较低；从图2可以直观的看到，中心点树木对周围的土地的绝对距离较为相近；而当局部邻域点作为中心点时，此时树木则成为了一个距离较远的点，故在加权时将带来较小的权重，以减少误差。由此可以看出，rank-order距离具有一定的聚类效果，适用于利用局部邻域约束的高光谱图像解混中；
[0071]
rank-order距离的主要原理是利用邻域间的共享信息，求出两点间的相对距离。rank-order距离计算步骤主要有三步：
[0072]
计算样本点xi和其局部邻域点的欧氏距离，并维持一个距离从小到大的顺序表；
[0073]
计算每xi与xj之的非对称距离d(xi,xj)，如式(7)：
[0074][0075]
由式(8)可知，d(xi,xj)表示的是给定的目标点xi的局部邻域点在xj的局部邻域点顺序表中所在位置下标之和；
[0076]
如图3所示，样本点xi和xj之间的非对称rank-order距离为：
[0077][0078]
计算对称的rank-order距离为：
[0079][0080]
参照图4，本发明将基于局部邻域约束的nmf算法和rank-order距离算法进行结合，提供一种基于rank-order局部邻域约束的nmf遥感图像解混方法，将同类像元赋予较高的权重，不同类的像元赋予较低的权重，在考虑空间距离的同时又具有一定的聚类效果；
[0081]
包括：
[0082]
s1：获取高光谱数据，通过vca-fcls对高光谱数据进行初始化，获得初始端元矩阵和初始丰度矩阵，将初始端元矩阵作为当前端元矩阵，将初始丰度矩阵作为当前丰度矩阵；
[0083]
s2：通过基于rank-order距离的自适应局部邻域算法对当前端元矩阵和当前丰度矩阵进行更新，获得更新后的端元矩阵和更新后的丰度矩阵；
[0084]
s3：若达到最大迭代次数或者误差小于规定的阈值，则输出更新后的端元矩阵和更新后的丰度矩阵；否则将更新后的端元矩阵作为当前端元矩阵，将更新后的丰度矩阵作为当前丰度矩阵，返回步骤s2。
[0085]
本实施例中，步骤s2具体为：
[0086]
s21：定义xi为高光谱数据x的第i列，xj为x的第j列，将xi设置为中心像元，将xj设置为局部邻域像元，计算获得中心像元与局部邻域像元间的不对称距离d(xi,xj)，以及局部邻域像元与中心像元间的不对称距离d(xj,xi)；
[0087]
具体的，在通过基于局部邻域约束的nmf算法确定了局部邻域之后，在计算权重时有所不同；假设中心像元为xi，其某一局部邻域像元为xj，计算xj对xi的距离权重，应当先计算xj与xi之间的不对称距离，即d(xi,xj)与d(xj,xi)，为计算不对称距离，首先，以xi为中心分别与其局部邻域像元计算欧氏距离，并将结果存入一个集合中，即l(xi)，同理计算出l(xj)；接下来将两集合中的距离从小到大依次进行排序，并分别求出各自集合中的点，在另一集合中所在的位置序号，在这里用与表示；利用式(10)可以分别求得中心像元与局部邻域像元间的不对称距离d(xi,xj)与d(xj,xi)：
[0088]
[0089][0090]
s22：通过d(xi,xj)和d(xj,xi)计算获得中心像元xi与局部邻域像元xj间的rank-order距离υ；
[0091]
具体的，利用式(11)即可求出xi与xj间的对称距离，即rank-order距离；在本发明中将求得的rank-order距离定义为υ；经过分析υ的取值范围在[0, ∞)；考虑其空间上的流行结构，需要选择一种非线性模型对距离权重进行加权；考虑到高斯核函数在υ为0时，权重为1，在υ为 ∞时，权重趋近于0，可以很好的符合对权重分配的预期；
[0092]
s23：计算获得局部邻域像元xj对中心像元xi的距离权重p
ij
，以及局部邻域像元xj对中心像元xi的总权重χi；
[0093]
s24：通过总权重χi计算获得权重矩阵w，通过距离权重p
ij
计算获得变量矩阵sw；
[0094]
s25：对当前端元矩阵进行更新，获得第一增广矩阵xf、第二增广矩阵af和更新后的端元矩阵；
[0095]
s26：通过权重矩阵w、变量矩阵sw、第一增广矩阵xf和第二增广矩阵af对当前丰度矩阵进行更新，获得更新后的丰度矩阵。
[0096]
本实施例中，步骤s23中距离权重p
ij
的计算公式为：
[0097][0098]
其中，υ表示rank-order距离，σ表示曲线的曲率调节参数，β表示目标点与局部邻域间丰度的相似性；
[0099]
总权重χi的计算公式为：
[0100]
χi＝∑
j∈n(i)
p
ij
,i＝1,2,
…
,n
ꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(13)
[0101]
其中，n(i)表示高光谱数据第i列像元的集合，n表示高光谱数据的像元总数。
[0102]
具体的，通过调节σ还可以调整曲线的曲率，在处理不同数据时更加灵活，故采用高斯核函数来对距离进行加权，经过高斯核函数处理后，最终的权重在0到1之间；采用混合丰度间投影的方式来判断目标点与局部邻域间丰度的相似性，这里用β表示，其中β的取值范围在0到1之间，经过高斯核函数加权的距离υ也在0到1之间，二者属于同一个数量级。
[0103]
本实施例中，步骤s24具体为：
[0104]
权重矩阵w的计算公式为：
[0105][0106]
其中，χi为局部邻域像元xj对中心像元xi的总权重；p表示端元数量；n表示高光谱数据的像元总数；
[0107]
sw的i列表示为的i列表示为
[0108]
sj表示丰度矩阵s的第j行，n(i)表示高光谱数据第i列像元的集合。
[0109]
本实施例中，步骤s25具体为：
[0110]
当前端元矩阵的更新计算公式为：
[0111]
a'
←
a.*(xs
t
)./(ass
t
)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(15)
[0112]
其中，a'表示更新后的端元矩阵，a表示当前端元矩阵，.*表示矩阵对应元素相乘，./表示矩阵元素对应相除，(.)
t
表示矩阵的转置，s表示丰度矩阵；
[0113][0114]
其中，δ是一个正常数，用于调整丰度和为1的约束，δ越大，s的每一列的和越接近1；在实验中δ通常设置为15。
[0115]
本实施例中，步骤s26中更新后的丰度矩阵的计算公式为：
[0116][0117]
其中，s'表示更新后的丰度矩阵，s表示当前丰度矩阵.*表示矩阵对应元素相乘，./表示矩阵元素对应相除，λ表示可调参数。
[0118]
本实施例中，步骤s3中阈值的计算公式为：
[0119][0120][0121]
其中，j(a,s)表示阈值，s表示丰度矩阵，a表示当前端元矩阵，s表示丰度矩阵，si表示丰度矩阵s的第i列，sj表示丰度矩阵s的第j行，λ表示可调参数，n(i)表示高光谱数据第i列像元的集合，n表示高光谱数据的像元总数。
[0122]
具体的，目标像元的丰度与其局部邻域像元丰度相似性越高，且对应的rank-order距离越小，对应的权重就越大，反之，丰度相似性越低且对应的rank-order距离越大，对应的权重就越小。
[0123]
本发明的实验效果：
[0124]
实验中的算法统一都经过了vca-fcls算法来进行了预处理，并采用解混中常用的评价指标，均方根误差(rmse)和光谱角距离(sad)来对比各算法的解混效果。具体定义如下：
[0125]
光谱角距离sad：衡量原始端元光谱特征a与其估计值之间的差异,公式如下：
[0126][0127]
均方根误差rmse：衡量标准丰度值s
ij
与解混后估计的丰度值之间的差异，公式如下：
[0128][0129]
模拟数据：
[0130]
在模拟数据实验中设置的一些关键参数，λ为0.1、δ为0.8、τ为15、端元p为4、最大迭代次数为300次。
[0131]
为了证明在模拟数据上基于rank-order局部邻域约束的nmf算法的优越性，本文用基于rank-order局部邻域约束的nmf算法对比了利用空间信息的gnm和alnwc-nmf算法，各算法在不同端元的sad与rmse的值如表1与表2所示：
[0132]
表1模拟数据下算法sad值的对比
[0133][0134]
表2模拟数据下算法rmse值的对比
[0135][0136]
在上述结果对比中可以发现，rlnwc-nmf在sad的评估标准下，误差均小于与其对比的gnmf和alnwc-nmf算法，在rmse评估标准下，改进后得算法，仅端元4的误差略高于alnwc-nmf算法，其余端元均有精度上的提升。由此可以证明，在模拟数据下rlnwc-nmf算法具有很好的解混效果，且较改进前的算法有着较大的提升。
[0137]
真实数据：
[0138]
在此实验中，参数设置如下，λ为0.1、τ为0.8、δ为15、端元p为12、最大迭代次数t为1000。
[0139]
图5为本发明方法与其他算法解混后丰度图对比，下面将对比真实数据下算法的解混效果，由于原矿图数据较大，故本次选取数据集的一部分来进行实验，在选取的数据中主要有10种地物，其分别为铋矿(buddingtonite)、蒙脱石(montmorillonite)、高岭土(kaolin)、明矾(alunite)、玉髓(chacedony)、高岭石(kaolinite)、白云母(muscovite)、镁铝榴石(pyrope)、蓝线石(dumortierite)，选取前六种地物解混后的丰度图如图5所示；
[0140]
解混后的端元与标准光谱库的对比如图6所示，(a)镁铝榴石,(b)高岭石#2,(c)蓝线石,(d)高岭石#1,(e)白云母,(f)蒙脱石；其中红色实线为nmf解混算法得到的结果，蓝色实线为标准光谱库结果；
[0141]
通过图6也可以看到本发明方法的光谱曲线与标准光谱库的曲线总体趋于拟合，且解混后所提取的各端元的丰度图较为清晰，具有较好的解混效果；下面通过表3，将本发
明方法与gnmf算法和alnwc-nmf算法作比较，来证明改进后算法在解混精度方面的优势。
[0142]
表3各算法解混出的端元sad值对比
[0143][0144]
通过表3中各种算法sad的对比，可以发现，本发明方法在10种地物中，有8种地物的sad值优于其他算法，其中，玉髓(chacedony)的解混精度提升最大。总体上，本发明方法的解混精度优于gnmf算法、alnwc-nmf算法。
[0145]
本发明提供一种基于rank-order局部邻域约束的nmf遥感图像解混系统，包括：
[0146]
初始化模块，用于获取高光谱数据，通过vca-fcls对高光谱数据进行初始化，获得初始端元矩阵和初始丰度矩阵，将初始端元矩阵作为当前端元矩阵，将初始丰度矩阵作为当前丰度矩阵；
[0147]
更新模块，用于通过基于rank-order距离的自适应局部邻域算法对当前端元矩阵和当前丰度矩阵进行更新，获得更新后的端元矩阵和更新后的丰度矩阵；
[0148]
输出模块，若达到最大迭代次数或者误差小于规定的阈值，则输出更新后的端元矩阵和更新后的丰度矩阵；否则将更新后的端元矩阵作为当前端元矩阵，将更新后的丰度矩阵作为当前丰度矩阵，返回更新模块。
[0149]
需要说明的是，在本文中，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、物品或者系统不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、物品或者系统所固有的要素。在没有更多限制的情况下，由语句“包括一个
……”
限定的要素，并不排除在包括该要素的过程、方法、物品或者系统中还存在另外的相同要素。
[0150]
上述本发明实施例序号仅仅为了描述，不代表实施例的优劣。在列举了若干装置的单元权利要求中，这些装置中的若干个可以是通过同一个硬件项来具体体现。词语第一、第二、以及第三等的使用不表示任何顺序，可将这些词语解释为标识。
[0151]
以上仅为本发明的优选实施例，并非因此限制本发明的专利范围，凡是利用本发明说明书及附图内容所作的等效结构或等效流程变换，或直接或间接运用在其他相关的技术领域，均同理包括在本发明的专利保护范围内。