基于数字表面模型辅助的高光谱图像本征分解方法
技术领域
本发明涉及一种高光谱图像本征分解方法，属于高光谱图像融合处理
技术领域
。背景技术
高光谱图像具有丰富的光谱信息，但是高维的光谱空间也造成信息的高度冗余，不利于信息的处理解译。为了更好的挖掘信息，现在已经设计出了各种特征提取的方法，但是这些信息提取方法具有本质的缺陷。当环境光照等成像条件变化时，成像得到的光谱也会随之变化，这使得得到的光谱具有很高的不确定性，这种不确定性体现在提取的特征上就会使得特征表达的信息不可靠。为了解决这个问题，研究基于物理成像模型的本征分解方法具有极大的意义。本征分解旨在研究物理成像过程，恢复反映地物本身属性的反射率，同时也反演出环境成像要素，比如光照。本征分解的难点在于，本征分解是一个欠定问题，在模型求解上存在挑战。
发明内容
本发明目的是为了解决现有高光谱图像本征分解精度低的问题，提供了一种基于数字表面模型辅助的高光谱图像本征分解方法。
本发明所述基于数字表面模型辅助的高光谱图像本征分解方法，它包括：
S1、输入高光谱图像和数字表面模型数据，计算获得几何分量；
S2、计算获得局部先验矩阵；
S3、计算获得非局部先验矩阵；
S4、根据几何分量、局部先验矩阵和非局部先验矩阵进行本征分解，输出高光谱反射率和环境光照。
优选的，S1所述计算获得几何分量的具体方法包括：
输入高光谱图像
输入数字表面模型高程数据
其中，hk＝[hk(λ1),hk(λ2),…,hk(λd)]T,k＝1,2,...,u表示每个像素的光谱特征，k＝1,2,...,u表示每个像素的索引，λ表示波长，d表示波段数，u表示高光谱图像像素的个数，z1,z2,…,zu表示各个像素对应的高程，表示域；
计算每个像素的法线，获得法线特征：
其中，表示法线在x、y、z空间坐标轴的投影；
计算获得几何分量J：J＝[J1,J2,…,Ju]ú；
其中，
c1、c2、c3、c4、c5是五个常数。
优选的，所述c1、c2、c3、c4、c5分别是：
c1＝0.429,c2＝0.512,c3＝0.743,c4＝0.886,c5＝0.248。
优选的，S2所述计算获得局部先验矩阵的具体方法包括：
遍历每个像素的索引k＝1,2,...,u，建立字典Dk＝[h1,…,hk-1,hk+1,…,hu,Id]；
其中，Id表示d维单位矩阵；
根据下式计算hk在字典Dk中的稀疏表示系数α：
minα‖α‖1subject to hk＝Dkα；
局部先验矩阵W的每个元素Wkj由下式赋值获得：
Wkj表示W的第k行第j列的元素；
αj表示第j列的稀疏表示系数，αj-1表示第j-1列的稀疏表示系数。
优选的，S3所述计算获得非局部先验矩阵的具体方法包括：
将高光谱图像每个像素的光谱特征hk与数字表面模型高程数据zk堆栈，建立为d+1维向量空间：[hk ú,zk],k＝1,2,...,u；
在d+1维向量空间中搜索每个像素的紧邻点，建立紧邻点的集合Qk；
非局部先验矩阵T的每个元素Tkj由下式赋值获得：
Tkj表示T的第k行第j列的元素。
优选的，S4所述根据几何分量、局部先验矩阵和非局部先验矩阵进行本征分解，输出高光谱反射率和环境光照的具体方法包括：
高光谱反射率
其中，rk＝[rk(λ1),rk(λ2),...,rk(λd)]T,k＝1,2,...,u表示每个像素的反射率；λ1、λ2、...、λd分别表示第1、2、…、d个图像通道对应的波长；
环境光照为L＝[L1,L2,...,L9]ú；
L1、L2、...、L9分别表示九维球谐光照的系数；
本征分解的总的代价函数为：
σr表示局部先验代价项的权重，σz表示非局部先验代价项的权重；
求偏导数，获得：
其中，中间变量G＝Iu+σr(Iu-W)(Iu-W)+σz(Iu-T)(Iu-T)；
其中，1d表示大小为d×1的全1列向量，Iu表示u×u单位矩阵；
令得到关于R和L的线性方程组，求解获得高光谱图像反射率R和环境光照L。
优选的，所述环境光照L是九维球谐光照系数。
本发明的优点：为了解决高光谱图像本征分解精度低的问题，本发明引入了DSM(Digital Surface Model，数字表面模型)数据的高程信息做指导，能够极大地提高本征分解的精度。本发明能够在DSM数据高程信息的指导下进行高精度的高光谱图像本征分解，，生成高光谱反射率和环境光照，方法简单，运算效率高，易于实现。本发明提出的高光谱图像本征分解方法建立高光谱成像的物理模型，从信号的角度出发结合高光谱图像的光谱信息和DSM数据的高程信息，消除高光谱图像由于光照导致的光谱信息退化导致的光谱不确定性。
附图说明
图1是本发明所述基于数字表面模型辅助的高光谱图像本征分解方法的原理框图；
图2是输入的高光谱图像的示意图；
图3是输入的数字表面模型数据的示意图；
图4是输出的高光谱反射率的示意图；
图5是输出的环境光照的示意图；
图6是采用本发明所述本征分解方法获得的高光谱反射率与其他方法获得的高光谱反射率的误差比较，其中，a表示采用本发明提出基于数字表面模型辅助的高光谱图像本征分解方法获得的高光谱反射率误差，b表示基于稀疏图编码的高光谱本征分解方法获得的高光谱反射率误差，c表示基于特征提取的高光谱本征分解方法获得的高光谱反射率误差，d表示基于纹理结构分离的高光谱本征分解方法获得的高光谱反射率误差。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动的前提下所获得的所有其他实施例，都属于本发明保护的范围。
需要说明的是，在不冲突的情况下，本发明中的实施例及实施例中的特征可以相互组合。
下面结合附图和具体实施例对本发明作进一步说明，但不作为本发明的限定。
具体实施方式一：下面结合图1说明本实施方式，本实施方式所述基于数字表面模型辅助的高光谱图像本征分解方法，它包括：
S1、输入高光谱图像和数字表面模型数据，计算获得几何分量；
S2、计算获得局部先验矩阵；
S3、计算获得非局部先验矩阵；
S4、根据几何分量、局部先验矩阵和非局部先验矩阵进行本征分解，输出高光谱反射率和环境光照。
具体实施方式二：本实施方式对实施方式一作进一步说明，S1所述计算获得几何分量的具体方法包括：
输入高光谱图像
输入数字表面模型高程数据
其中，hk＝[hk(λ1),hk(λ2),…,hk(λd)]T,k＝1,2,...,u表示每个像素的光谱特征，k＝1,2,...,u表示每个像素的索引，λ表示波长，d表示波段数，u表示高光谱图像像素的个数，z1,z2,…,zu表示各个像素对应的高程，表示域；
计算每个像素的法线，获得法线特征：
其中，表示法线在x、y、z空间坐标轴的投影；
计算获得几何分量J：J＝[J1,J2,…,Ju]ú；
其中，
c1、c2、c3、c4、c5是五个常数。
进一步的，所述c1、c2、c3、c4、c5分别是：
c1＝0.429,c2＝0.512,c3＝0.743,c4＝0.886,c5＝0.248。
具体实施方式三：本实施方式对实施方式二作进一步说明，S2所述计算获得局部先验矩阵的具体方法包括：
遍历每个像素的索引k＝1,2,...,u，建立字典Dk＝[h1,…,hk-1,hk+1,…,hu,Id]；
其中，Id表示d维单位矩阵；
根据下式计算hk在字典Dk中的稀疏表示系数α：
minα‖α‖1subject to hk＝Dkα；
局部先验矩阵W的每个元素Wkj由下式赋值获得：
Wkj表示W的第k行第j列的元素；
αj表示第j列的稀疏表示系数，αj-1表示第j-1列的稀疏表示系数。
具体实施方式四：本实施方式对实施方式三作进一步说明，S3所述计算获得非局部先验矩阵的具体方法包括：
将高光谱图像每个像素的光谱特征hk与数字表面模型高程数据zk堆栈，建立为d+1维向量空间：[hk ú,zk],k＝1,2,...,u；
在d+1维向量空间中搜索每个像素的紧邻点，建立紧邻点的集合Qk；
非局部先验矩阵T的每个元素Tkj由下式赋值获得：
Tkj表示T的第k行第j列的元素。
具体实施方式五：本实施方式对实施方式四作进一步说明，S4所述根据几何分量、局部先验矩阵和非局部先验矩阵进行本征分解，输出高光谱反射率和环境光照的具体方法包括：
高光谱反射率
其中，rk＝[rk(λ1),rk(λ2),...,rk(λd)]T,k＝1,2,...,u表示每个像素的反射率；λ1、λ2、...、λd分别表示第1、2、…、d个图像通道对应的波长；
环境光照为L＝[L1,L2,...,L9]ú；
L1、L2、...、L9分别表示九维球谐光照的系数；
本征分解的总的代价函数为：
σr表示局部先验代价项的权重，σz表示非局部先验代价项的权重；
求偏导数，获得：
其中，中间变量G＝Iu+σr(Iu-W)(Iu-W)+σz(Iu-T)(Iu-T)；
其中，1d表示大小为d×1的全1列向量，Iu表示u×u单位矩阵；
令得到关于R和L的线性方程组，求解获得高光谱图像反射率R和环境光照L。
进一步的，所述环境光照L是九维球谐光照系数。
本发明中，高光谱图像与DSM高程数据具有的互补性、异构性，更具体地说，高光谱图像具有丰富的光谱信息，但是空间信息是三维到二维图像的退化；DSM能获取精确的三维空间信息，但是光谱信息相对匮乏。
为了验证本发明提出的高光谱图像本征分解方法的性能，针对一组实拍的高光谱图像和DSM数据进行了验证，如图2所示，是输入的高光谱图像，如图3所示，是输入的DSM数据，图4是采用本发明提出的高光谱图像本征分解方法获得的高光谱反射率，图5是采用本发明提出的高光谱图像本征分解方法获得的环境光照，如图6所示，是采用本发明所述本征分解方法获得的高光谱反射率与其他方法获得的高光谱反射率的误差比较，图6显示出了本发明提出的高光谱图像本征分解方法的有效性，图6显示，本发明提出的高光谱图像本征分解方法恢复的高光谱反射率误差是最小的。
虽然在本文中参照了特定的实施方式来描述本发明，但是应该理解的是，这些实施例仅仅是本发明的原理和应用的示例。因此应该理解的是，可以对示例性的实施例进行许多修改，并且可以设计出其他的布置，只要不偏离所附权利要求所限定的本发明的精神和范围。应该理解的是，可以通过不同于原始权利要求所描述的方式来结合不同的从属权利要求和本文中所述的特征。还可以理解的是，结合单独实施例所描述的特征可以使用在其他所述实施例中。
