基于稀疏表示的高光谱混合像元分解方法与流程

1.本发明涉及混合像元领域，尤其涉及基于稀疏表示的高光谱混合像元分解方法。


背景技术：

2.遥感器所获取的地面反射或发射光谱信号是以像元为单位记录的，一个像元内仅包含一种类型，这种像元称为纯像元。然而，多数情况下一个像元内往往包含多种地表类型，这种像元就是混合像元。混合像元记录的是多种地表类型的综合光谱信息，混合像元是指在一个像元内存在有不同类型的地物，主要出现在地类的边界处。混合像元的存在是影响识别分类精度的主要因素之一，特别是对线状地类和细小地物的分类识别影响较为突出，在土地利用遥感动态监测工作中，经常遇到混合像元的难题，解决这一问题的关键在于通过一定方法找出组成混合像元的各种典型地物的比例。
3.监督类的方法大多基于高光谱图像的几何特性，旨在图像的特征空间中寻找端元像元组成的凸面几何体的顶点，以此作为端元。然而，如果成像光谱仪的空间分辨率不够高，很可能导致图像中没有纯净像元存在，纯净像元的图像中运用几何学方法提取端元可能会导致端元不准确，提取出的仍然是混合了其他地物的“非纯净端元”，利用这些不纯净的端元进行后续的丰度反演步骤，势必会导致进一步的误差为此，我们提出基于稀疏表示的高光谱混合像元分解方法。


技术实现要素：

4.基于背景技术存在的技术问题，本发明提出了基于稀疏表示的高光谱混合像元分解方法，以解决背景技术中提出的问题。
5.本发明提供如下技术方案：
6.基于稀疏表示的高光谱混合像元分解方法，该方法为：
7.a、稀疏表示：信号的线性表示和高光谱线性混合模型的原理和表达都是一致的，可表示为：
8.y＝ψa
ꢀꢀ
(1.1)
9.其中a为l维列向量，若a中的非零元素个数k＜＜l，即a是k
‑
稀疏信号，则称信号y是具有k
‑
稀疏性可压缩信号，ψ通常也被称为字典，当ψ列数大于行数的长方阵时，此式的线性方程组(1.1)对应这一个欠定系统，a存在无穷多解，稀疏表示的目的，就是在这些所有可能的解中寻找一个尽可能稀疏的解，即该解的0
‑
范数||a||0最小，该过程可表示为下面的优化问题：
[0010][0011]
考虑到噪声，式(1.2)可以改写成为：
[0012]
[0013]
然而，式(1.2)和(1.3)的优化问题都是np问题，难以求解，因此，大家一般用基于1
‑
范数最小化问题来近似表示：
[0014][0015][0016]
其变化形式为：
[0017][0018]
其中第一项衡量解的拟合程度，第二项衡量解的稀疏程度。正则参数λ在求解过程中调整这两项所占的比重。由此，np难题就成了具有唯一最优解且易解的凸优化问题；
[0019]
b、su
‑
admm算法：(1)、admm原理
[0020]
一个由两个函数组成的非限制性优化问题表示为：
[0021][0022]
想要变量分离，引入一个新的变量u作为函数f2的因变量，令u＝gx，就得到一个限制优化问题：
[0023]
min f1(x) f2(u)s.t.gx＝u
ꢀꢀ
(1.8)
[0024]
这样就可以用增广拉格朗日乘数方法求解问题(1.7)，加入一个与x无关的常数项d：
[0025][0026]
d
k 1
＝d
k
‑
(gx
k 1
‑
u
k 1
)
ꢀꢀ
(1.10)
[0027]
设k＝0，admm的迭代步骤为：
[0028][0029][0030]
d
k 1
←
d
k
‑
(gx
k 1
‑
u
k 1
)
ꢀꢀ
(1.13)
[0031]
b、su
‑
admm算法描述及流程
[0032]
假设一幅高光谱图像中每个像元的观测值为列向量y，端元矩阵用光谱库a表示，y中每个端元的丰度向量为x。根据admm原理，可以将稀疏表示的高光谱混合像元分解问题(1.6)中的两项当作admm的两个函数：代表解的拟合程度的第一项当作f1(x)；代表解的稀疏程度的第二项当作f2(gx)，即f2(x)＝λ||x||1，g≡i。用admm解决高光谱混合像元分解问题，还要加上对丰度的“非负”约束和“和为1”约束。引入一个新的变量u作为函数f2的因变量，令u＝gx，则高光谱混合像元分解问题可表示为admm形式：
[0033][0034]
其中约束条件为：s.t.1
t
x＝1；u≥0；x＝u。参照式(1.11)
‑
(1.13)的迭代公式即可求解，改变步长k
←
k 1，直到达到一定的精度时停止迭代，其中参数λ的设置取决于互相实际包含端元数目与光谱库中端元数目之间的比例，如果比例很小，则丰度矩阵将具有很大的稀疏性，λ的值应该大一些，才能加大丰度稀疏性的约束力度；如果比例很大，则λ的值应该小一些。总的来说，λ的值应该在0
‑
0.5之间。
[0035]
本发明提供了基于稀疏表示的高光谱混合像元分解方法，本发明采用su
‑
admm算法将已有的光谱库作为参考端元矩阵，利用丰度的稀疏性，运用稀疏表示的方法进行高光谱混合像元分解，端元是光谱库里的子集，所以算法提取出的端元准确度较高，不会存在提取出不存在的物质光谱；丰度分布也较准确、符合实际。
附图说明
[0036]
图1为本发明信号重建误差随端元个数变化情况图；
[0037]
图2为本发明真实端元的光谱曲线图；
[0038]
图3为本发明su
‑
admm分解丰度图；
[0039]
图4为本发明从光谱库中选取出的各矿物质光谱曲线图。
具体实施方式
[0040]
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0041]
请参阅图1
‑
4，本发明提供一种技术方案：
[0042]
基于稀疏表示的高光谱混合像元分解方法，该方法为：
[0043]
a、稀疏表示：信号的线性表示和高光谱线性混合模型的原理和表达都是一致的，可表示为：
[0044]
y＝ψa
ꢀꢀ
(1.1)
[0045]
其中a为l维列向量，若a中的非零元素个数k＜＜l，即a是k
‑
稀疏信号，则称信号y是具有k
‑
稀疏性可压缩信号，ψ通常也被称为字典，当ψ列数大于行数的长方阵时，此式的线性方程组(1.1)对应这一个欠定系统，a存在无穷多解，稀疏表示的目的，就是在这些所有可能的解中寻找一个尽可能稀疏的解，即该解的0
‑
范数||a||0最小，该过程可表示为下面的优化问题：
[0046][0047]
考虑到噪声，式(1.2)可以改写成为：
[0048][0049]
然而，式(1.2)和(1.3)的优化问题都是np问题，难以求解，因此，大家一般用基于
1
‑
范数最小化问题来近似表示：
[0050][0051][0052]
其变化形式为：
[0053][0054]
其中第一项衡量解的拟合程度，第二项衡量解的稀疏程度。正则参数λ在求解过程中调整这两项所占的比重。由此，np难题就成了具有唯一最优解且易解的凸优化问题；
[0055]
b、su
‑
admm算法：(1)、admm原理
[0056]
一个由两个函数组成的非限制性优化问题表示为：
[0057][0058]
想要变量分离，引入一个新的变量u作为函数f2的因变量，令u＝gx，就得到一个限制优化问题：
[0059]
min f1(x) f2(u)s.t.gx＝u
ꢀꢀ
(1.8)
[0060]
这样就可以用增广拉格朗日乘数方法求解问题(1.7)，加入一个与x无关的常数项d：
[0061][0062]
d
k 1
＝d
k
‑
(gx
k 1
‑
u
k 1
)
ꢀꢀ
(1.10)
[0063]
设k＝0，admm的迭代步骤为：
[0064][0065][0066]
d
k 1
←
d
k
‑
(gx
k 1
‑
u
k 1
)
ꢀꢀ
(1.13)
[0067]
b、su
‑
admm算法描述及流程
[0068]
假设一幅高光谱图像中每个像元的观测值为列向量y，端元矩阵用光谱库a表示，y中每个端元的丰度向量为x。根据admm原理，可以将稀疏表示的高光谱混合像元分解问题(1.6)中的两项当作admm的两个函数：代表解的拟合程度的第一项当作f1(x)；代表解的稀疏程度的第二项当作f2(gx)，即f2(x)＝λ||x||1，g≡i。用admm解决高光谱混合像元分解问题，还要加上对丰度的“非负”约束和“和为1”约束。引入一个新的变量u作为函数f2的因变量，令u＝gx，则高光谱混合像元分解问题可表示为admm形式：
[0069]
[0070]
其中约束条件为：s.t.1
t
x＝1；u≥0；x＝u。参照式(1.11)
‑
(1.13)的迭代公式即可求解，改变步长k
←
k 1，直到达到一定的精度时停止迭代，其中参数λ的设置取决于互相实际包含端元数目与光谱库中端元数目之间的比例，如果比例很小，则丰度矩阵将具有很大的稀疏性，λ的值应该大一些，才能加大丰度稀疏性的约束力度；如果比例很大，则λ的值应该小一些。总的来说，λ的值应该在0
‑
0.5之间。
[0071]
su
‑
admm的算法流程图如表5.1所示。
[0072]
表5.1su
‑
admm算法流程
[0073][0074]
实验结果与分析
[0075]
模拟数据实验
[0076]
实验一、解混效果随端元个数变化情况
[0077]
首先将usgs光谱库中光谱曲线进行挑选处理，计算并将其中任意两条光谱曲线的光谱角度距离小于0.05rad(大概等同于角度小于3
°
)的曲线去除，剩下差异较大的342条光谱曲线，波段数为224，将此新的光谱库记为a。模拟数据的产生方式同前面一样，设定像元混合程度很高，超出了基于几何学算法的适用范畴。像元个数n＝1000，在使用算法进行像元解混时，将a做端元矩阵，估计出丰度矩阵构成模拟图像的端元子集的选取分为两种：第一种是光谱曲线差异较大的，并且它们中两两曲线最小光谱角距离≧7
°
，标记为θ
min
(a)≥7
°
；第二种是光谱曲线差异相对小的，两两曲线最小光谱角距离≤4
°
，标记为θ
min
(a)≤4
°
。
[0078]
稀疏回归后的结果，本文用信号重建误差(signal
‑
to
‑
reconstructionerror，sre)衡量
[0079][0080]
其中x是每个像元的实际丰度值，是使用su
‑
admm算法估计出的每个像元的丰度向量，<
·
>代表所有样本(n＝1000个像元)的平均值。图1是信号重建误差随着模拟数据中
实际端元个数变化而变化的情况。
[0081]
图1(a)中设定的噪声为0，对丰度矩阵的约束为0的情况下丰度估计情况，此时算法退化为约束的最小二乘法(cls)。可以看出，随着端元子集中端元个数的增加，信号重建误差均降低，而端元曲线差异较大时(θ
min
(a)≥7
°
)，信号重建误差也相对较高。图2(b)中，加入了一定的噪声，信噪比snr＝25db，由于实际端元个数≤10，相对于光谱库来说具有很大的稀疏性，因此，取稀疏约束程度为λ＝5
×
10
‑2。可以看出，加入约束的算法的解混效果优于不加约束的，即在图像中像元混合程度较高的情况下，基于稀疏表示的su
‑
admm方法解混效果优于cls。
[0082]
实验二、验证su
‑
admm算法求解丰度值的解混精度
[0083]
模拟图像的丰度矩阵s已知，将s与该算法计算出的丰度矩阵比较，计算rmse，可以验证此算法丰度估计的精度。将a当作端元矩阵，求解a所对应的丰度矩阵，找出丰度值最大的p条曲线，即为此算法求得的端元曲线将提取的端元结果与真实的端元m相比较，计算光谱角距离sad，从而验证此算法能否正确地找出图像的端元。本次试验取端元个数p＝5，n＝1000，snr＝25db，λ＝5
×
10
‑2。其中5条真实端元光谱曲线如图2所示。
[0084]
此次试验经过稀疏回归分解后，求出所有光谱曲线对应的5个最大的丰度值为：0.2207、0.1888、0.1421、0.1403、0.0987，把它们对应的光谱曲线作为端元，发现与构成模拟图像的5条曲线完全相同。由此可看出，su
‑
admm算法能够正确地找到图像的端元。但是，实验结果中还存在一些其他的光谱曲线对应的丰度值不为零，这主要是因为光谱库中存在相似的曲线。因此，在端元数目未知的情况下，需要设定丰度值阈值，才能确定合理的端元个数。另一方面，需要对初始光谱库进行挑选处理，防止过多相似曲线的存在。
[0085]
根据su
‑
admm估计出的丰度值与真实的丰度值计算均方根误差rmse，并用同样的数据采用非负约束的最小二乘法和全约束最小二乘法进行丰度估计，与su
‑
admm方法进行比较，结果如表5.2所示。
[0086]
表5.2三种算法进行丰度估计的精度比较
[0087][0088][0089]
从表5.2中可以看出，基于稀疏回归的分解算法进行丰度估计精度虽然比fcls差，但是比cls要好得多；其中两种最小二乘法是以真实端元作为端元矩阵进行丰度反演的，需要端元的先验信息，而su
‑
admm是以光谱库为端元矩阵，不知道具体的端元的情况下进行分解，因此具有一定的优势。
[0090]
真实数据实验
[0091]
数据采用了1997年拍摄的美国cuprite地区的aviris高光谱数据数据，端元矩阵(字典)选取方法为：获取了usgs光谱库数据，数据包括了498条端元光谱，波段数为420，数值为地物的反射率。将其重采样，得到224波段的498条光谱。将此矩阵作为端元矩阵m，大小为224
ⅹ
498。
[0092]
通过实验将su
‑
admm的解混效果与另外两种解混方法进行比较，其他两种方法分
别为监督类混合像元分解方法vca fcls和非监督类解混方法mvc
‑
nmf。vca fcls和mvc
‑
nmf分解后各个矿物质的丰度图。图3是su
‑
admm方法分解的丰度图。
[0093]
从图3可以看出，各矿物质被清晰地分离出来，而且相互之间重叠的部分较少，与分布图对比，发现丰度反演很准确，明显优于其他两种方法。可见，在光谱库可用的前提下，基于稀疏回归的混合像元分解方法很适用于高光谱图像。图4是经过分解后确定的端元子集中端元的光谱曲线。
[0094]
几种解混方法之间的关系
[0095]
对比几种混合像元分解方法，发现cls、fcls、mvc
‑
nmf与su
‑
admm有很大的相似性。根据1.3中su
‑
admm算法描述可知，su
‑
admm算法的原始目标函数为
[0096][0097]
约束条件为x≥0。当式(1.15)中λ＝0时，问题就变成了“非负”约束的最小二乘法，记为pcls；在此基础上再加入“和为1”约束，即1
t
x＝1，就成了全约束最小二乘法的求解问题，记为pfcls。而第二章中介绍的mvc
‑
nmf的求解问题，也是在fcls的基础上在目标函数加入了对端元矩阵a的约束，变成了min||y
‑
ax||
f
/2 λf(a)，根据高光谱图像的特征，式中a也要满足非负约束，记为pmvc
‑
nmf；而第三章中本文作者提出的vsc
‑
nmf算法，则是在fcls的目标函数中对两个矩阵a和x都进行了不同方式的约束，将第3.2.2节式(3.6)中端元矩阵m、丰度矩阵s分别表示为a和x，vsc
‑
nmf目标函数就变成
[0098][0099]
s.t.:a≥0,x≥0
ꢀꢀꢀ
(1.16)
[0100]
记为pvsc
‑
nmf。将上述各解混算法之间的联系及比较总结为表5.3。
[0101]
表5.3几种混合像元分解算法比较
[0102][0103]
以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，根据本发明的技术方案及其
发明构思加以等同替换或改变，都应涵盖在本发明的保护范围之内。