CN112712034A - 一种高光谱图像的解混方法及系统 - Google Patents
一种高光谱图像的解混方法及系统 Download PDFInfo
- Publication number
- CN112712034A CN112712034A CN202011630071.1A CN202011630071A CN112712034A CN 112712034 A CN112712034 A CN 112712034A CN 202011630071 A CN202011630071 A CN 202011630071A CN 112712034 A CN112712034 A CN 112712034A
- Authority
- CN
- China
- Prior art keywords
- hyperspectral image
- dimensional
- dimensional hyperspectral
- end member
- loss function
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 238000000034 method Methods 0.000 title claims abstract description 75
- 239000011159 matrix material Substances 0.000 claims abstract description 71
- 238000000605 extraction Methods 0.000 claims abstract description 17
- 239000003638 chemical reducing agent Substances 0.000 claims description 41
- 230000003595 spectral effect Effects 0.000 claims description 21
- 230000003190 augmentative effect Effects 0.000 claims description 20
- 230000009977 dual effect Effects 0.000 claims description 17
- 238000001228 spectrum Methods 0.000 claims description 16
- 239000013598 vector Substances 0.000 claims description 10
- 238000004458 analytical method Methods 0.000 claims description 8
- 238000012545 processing Methods 0.000 claims description 8
- 238000006243 chemical reaction Methods 0.000 claims description 6
- OAICVXFJPJFONN-UHFFFAOYSA-N Phosphorus Chemical compound [P] OAICVXFJPJFONN-UHFFFAOYSA-N 0.000 claims description 5
- 239000000126 substance Substances 0.000 claims description 4
- 230000009466 transformation Effects 0.000 claims description 4
- 238000004364 calculation method Methods 0.000 abstract description 7
- 230000008569 process Effects 0.000 abstract description 7
- 230000000694 effects Effects 0.000 abstract description 6
- 230000006870 function Effects 0.000 description 72
- 238000005457 optimization Methods 0.000 description 7
- 238000003384 imaging method Methods 0.000 description 5
- 230000009467 reduction Effects 0.000 description 5
- 238000002474 experimental method Methods 0.000 description 4
- 238000004422 calculation algorithm Methods 0.000 description 3
- 238000004088 simulation Methods 0.000 description 3
- 230000004087 circulation Effects 0.000 description 2
- 125000004122 cyclic group Chemical group 0.000 description 2
- 238000007405 data analysis Methods 0.000 description 2
- 238000000354 decomposition reaction Methods 0.000 description 2
- 238000013507 mapping Methods 0.000 description 2
- 239000000463 material Substances 0.000 description 2
- 239000000203 mixture Substances 0.000 description 2
- 241000243251 Hydra Species 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000015572 biosynthetic process Effects 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 238000010191 image analysis Methods 0.000 description 1
- QRXWMOHMRWLFEY-UHFFFAOYSA-N isoniazide Chemical compound NNC(=O)C1=CC=NC=C1 QRXWMOHMRWLFEY-UHFFFAOYSA-N 0.000 description 1
- 238000010801 machine learning Methods 0.000 description 1
- 230000014759 maintenance of location Effects 0.000 description 1
- 238000002156 mixing Methods 0.000 description 1
- 238000013139 quantization Methods 0.000 description 1
- 230000005855 radiation Effects 0.000 description 1
- 238000011946 reduction process Methods 0.000 description 1
- 238000002310 reflectometry Methods 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Chemical compound O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
- G06V20/00—Scenes; Scene-specific elements
- G06V20/10—Terrestrial scenes
- G06V20/13—Satellite images
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T5/00—Image enhancement or restoration
- G06T5/70—Denoising; Smoothing
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
- G06V10/00—Arrangements for image or video recognition or understanding
- G06V10/40—Extraction of image or video features
- G06V10/44—Local feature extraction by analysis of parts of the pattern, e.g. by detecting edges, contours, loops, corners, strokes or intersections; Connectivity analysis, e.g. of connected components
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
- G06V20/00—Scenes; Scene-specific elements
- G06V20/10—Terrestrial scenes
- G06V20/194—Terrestrial scenes using hyperspectral data, i.e. more or other wavelengths than RGB
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
- Y02A40/00—Adaptation technologies in agriculture, forestry, livestock or agroalimentary production
- Y02A40/10—Adaptation technologies in agriculture, forestry, livestock or agroalimentary production in agriculture
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Multimedia (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Astronomy & Astrophysics (AREA)
- Remote Sensing (AREA)
- Image Processing (AREA)
Abstract
本发明公开了一种高光谱图像的解混方法及系统，包括获取待解混三维高光谱图像，对待解混三维高光谱图像进行重构处理，得到二维高光谱图像；进行端元提取，建立端元矩阵；利用端元矩阵，构建损失函数；其中，损失函数用于求解二维高光谱图像的丰度，损失函数包含先验正则项；求解损失函数，输出二维高光谱图像的丰度，二维高光谱图像的丰度即为待解混三维高光谱图像的解混结果；本发明通过获得二维高光谱图像，并对二维高光谱图像进行端元提取，根据端元矩阵确定损失函数，通过在损失函数中加入先验正则项，将图像的空间相关性的先验知识加入到图像的解混中，有效提高了解混精度；避免了手动选择正则项的类型，计算过程复杂，图像解混效果差。
Description
技术领域
本发明属于图像处理和机器学习技术领域，特别涉及一种高光谱图像的解混方法及系统。
背景技术
高光谱将表示观测场景辐射属性的光谱与表示空间和几何关系的同源图像相结合，高光谱图像利用很多很窄的电磁波段从观测场景中感知有用信息，获得的数据除含有传统图像的长、宽两个维度，还包括光谱维度；即在每一个像素上，成像设备获得光线在该点上数百个甚至上千个波段的反射率。高光谱图像的光谱分辨率可以达到纳米级，能在电磁波谱的紫外、可见光、近红外和短波红外区域获取许多波段非常窄而且光谱连续的图像数据。丰富的光谱信息使得许多原来在单通道或者多通道图像下不能解决的问题得以解决，因此高光谱图像分析在诸多领域有着广泛的应用。然而，由于成像设备空间分辨率低、物质混杂及光谱串扰等原因，高光谱图像中一个像素的观测光谱可能是多种纯物质光谱的混合。因此，对图像中的像元进行分解已经成为高光谱图像广泛应用的重要前提；混合像元可以看作是一组基向量按照一定的比例组合而成；其中，基向量是“端元”，一定的比例即为“丰度”。
目前，大多数高光谱图像解混算法都是基于线性模型提出的，线性混合模型是一种简单且物理上可解释的混合模型；该模型假设端元之间的多次散射可以忽略不计，混合像元可以看作像元区域内端元光谱信号按照不同比例分数的线性叠加；传统的光谱解混方法都是基于单像素进行的，忽略了高光谱图像潜在的空间连续性；加入合适的空间和谱间先验知识可以有效的提高高光谱图像解混的精度。
现有技术中，对优化问题添加正则项是一种被广泛使用的添加图像先验信息的方法；例如：SUnSAL-TV方法通过使用TV正则项来添加高光谱图像相邻像元连续性信息，并使用L1范数来促进解混丰度的稀疏性；Graph TV使用图正则项来增加高光谱图像空间维和光谱维的相似性。现有的添加方法虽然能为图像解混提供一定的先验信息，但是还存在以下缺点：(1)缺乏灵活性，需要手动选择正则项的种类；(2)计算复杂，对于不同的正则项需要设计不同的优化问题求解方法，求解复杂，无法实现高光谱图像的快速解混；(3)对于噪声含量高的图像解混效果差，计算结果精度低。因此本发明中提出使用一种高光谱图像的解混方法及系统，解决现有技术中的不足。
发明内容
针对现有技术中存在的技术问题，本发明提供了一种高光谱图像的解混方法及系统，以解决现有的解混方法，缺乏灵活性，计算过程复杂，图像解混效果差，计算结果精度低的技术问题。
为达到上述目的，本发明采用的技术方案为：
本发明提供了一种高光谱图像的解混方法，包括以下步骤：
获取待解混三维高光谱图像Y3D，对待解混三维高光谱图像Y3D进行重构处理，得到二维高光谱图像Y；
对二维高光谱图像Y进行端元提取，建立端元矩阵E；
利用端元矩阵E，构建损失函数；其中，损失函数用于求解二维高光谱图像Y的丰度，损失函数包含先验正则项；
求解损失函数，输出二维高光谱图像Y的丰度，二维高光谱图像Y的丰度即为待解混三维高光谱图像的解混结果。
进一步的，对目标场景进行成像，得到待解混三维高光谱图像Y3D；对待解混三维高光谱图像Y3D，进行矩阵变换，得到二维高光谱图像数据矩阵Y；
其中，待解混三维高光谱图像Y3D的表达式为：
二维高光谱图像Y的表达式为：
N＝L×W
其中，B为待解混三维高光谱图像的光谱通道个数，L为待解混三维高光谱图像的空间维长；W为待解混三维高光谱图像的空间维宽，N为二维高光谱图像的像素总个数。
进一步的，建立端元矩阵E，具体为：
确定二维高光谱图像Y的端元数量R，采用向量成分分析端元提取方法，进行端元提取，得到二维高光谱图像Y的R个端元，并利用二维高光谱图像Y的R个端元，构建端元矩阵E。
进一步的，确定二维高光谱图像数据矩阵的端元数量R，具体操作如下：
采用最小误差识别方法，对二维高光谱图像Y进行最小化投影；利用最小化投影后的噪声能量和投影信号误差能量之和，在正交子空间中，确定二维高光谱图像的端元数量。
进一步的，采用基于几何理论的向量成分分析端元提取方法，进行端元提取；具体操作为：
初始化正交子空间，将二维高光谱图像在正交子空间中进行投影；
计算二维高光谱图像中每个像元在正交子空间的投影距离；
将二维高光谱图像中在正交子空间的投影距离最大的像元，作为候选端元，将候选端元加入到端元集合中；
更新正交子空间，将二维高光谱图像在更新后的正交子空间中进行投影，并计算二维高光谱图像中每个像元在更新后的正交子空间的投影距离；
将二维高光谱图像中在更新后的正交子空间的投影距离最大的像元，作为新的候选端元，并将新的候选端元加入到已有端元集合中；
循环更新正交子空间，至端元集合中，包含R个候选端元，得到端元矩阵E。
进一步的，损失函数的表达式为：
其中，Φ(*)为先验正则项，λ为正则项强度；ai为二维高光谱图像中第i个像素的丰度；yi为二维高光谱图像中第i个像素的光谱；i为二维高光谱图像的像素个数；A为丰度矩阵。
进一步的，求解损失函数，具体包括以下步骤：
引入辅助变量矩阵Z，对损失函数进行等价形式转化，得到转化后的损失函数；
获取转化后的损失函数的增广拉格朗日函数Lρ；
利用交替方向乘子法，求解转化后的损失函数的增广拉格朗日函数Lρ，得到二维高光谱图像数据矩阵的丰度，即得到待解混高光谱图像的解混结果。
进一步的，转化后的损失函数的增广拉格朗日函数Lρ的表达式为：
其中，V为对偶变量矩阵；||EA-Z||F为EA-Z矩阵的F范数；ρ为惩罚因子。
进一步的，利用交替方向乘子法，求解转化后的损失函数的增广拉格朗日函数Lρ，具体操作如下：
利用交替方向乘子法，将转化后的损失函数的增广拉格朗日函数分解为第一子问题和第二子问题；
其中，第一子问题为：
其中，ak+1,i为第i个像素在第k+1次迭代的丰度值，ρk为第k次迭代中的惩罚因子，
为第i个像素的中间变量，zk,i为第i个像素的辅助变量，uk,i为第i个像素的对偶变量；
第二子问题为：
其中，
为整幅二维高光谱图像第k次迭代中的辅助变量；
对二维高光谱图像中的每个像素，利用FCLS方法求解第一子问题，得到每个像素对应的丰度，并将所有像素的丰度拼接整合，获得整幅二维高光谱图像在第k+1次迭代的丰度Ak+1；利用降噪器求解第二子问题，获得整幅二维高光谱图像第k+1次迭代中的辅助变量；
循环迭代，求解第一子问题和第二子问题，至迭代循环次数达到最大循环次数K，输出高光谱图像的解混结果。
本发明还提供了一种高光谱图像的解混系统，包括获取模块、端元矩阵模块、函数模块及求解模块；
获取模块，获取待解混三维高光谱图像Y3D，对待解混三维高光谱图像Y3D进行重构处理，得到二维高光谱图像Y；
端元矩阵模块，用于对二维高光谱图像Y进行端元提取，建立端元矩阵E；
函数模块，用于利用端元矩阵E，构建损失函数；其中，损失函数用于求解二维高光谱图像Y的丰度，损失函数包含先验正则项；
求解模块，用于求解损失函数，输出二维高光谱图像Y的丰度，二维高光谱图像Y的丰度即为待解混三维高光谱图像的解混结果。
与现有技术相比，本发明的有益效果为：
本发明提供了一种高光谱图像的解混方法及系统，通过将三维高光谱图像进行重构处理，获得二维高光谱图像，并对二维高光谱图像进行端元提取，得到端元矩阵；根据端元矩阵确定损失函数，通过在损失函数中加入先验正则项，将图像的空间相关性的先验知识加入到图像的解混中，有效提高了解混精度；避免了手动选择正则项的类型，计算过程复杂，图像解混效果差。
进一步的，通过采用交替方向乘子法对转化后的损失函数的增广拉格朗日函数进行求解，将高光谱图像解混问题分解为两个子问题，其中第一子问题是标准的全约束最小二乘问题，能够采用FCLS方法进行求解，第二子问题是增加先验的图像降噪问题，通过采用加入先验信息的降噪器进行求解；通过采用降噪器自动加入图像的先验信息，避免了手动选择正则项；同时，通过降噪器能够对图像起到降噪作用，提高图像解混的精度；本发明优选的，采用二维图像降噪器NLM和二维图像降噪器BM3D在解混过程引入图像的空间信息，采用三维图像降噪器BM4D能够在解混过程引入图像的空间和谱间信息；通过降噪器的使用，有效提高了低信噪比图像的解混能力，计算精度较高，解混效果较好。
进一步的，采用降噪器代替求解带正则项的优化问题求解，图像降噪器在设计时会利用图像的空间相关性等先验知识，可以将先验知识加入到图像解混中，提高解混精度；不同的降噪器可以加入不同的先验知识，通过这种操作避免了手动选择正则项的类型；本发明中的降噪器不仅能够引入先验信息，同时，降噪器的降噪特性也提高了算法对噪声的鲁棒性。
附图说明
图1为本发明所述的解混方法的流程图；
图2为不同信噪比原始高光谱图像和解混后高光谱图像的映射对比图；
图3为Urban数据实验丰度估计映射对比图。
具体实施方式
为了使本发明所解决的技术问题，技术方案及有益效果更加清楚明白，以下具体实施例，对本发明进行进一步的详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。
如附图1所示，本发明所述的一种高光谱图像的解混方法，包括以下步骤：
步骤1、对目标场景进行成像，得到待解混三维高光谱图像Y3D；对待解混三维高光谱图像Y3D，进行矩阵变换，得到二维高光谱图像Y；
其中，待解混三维高光谱图像Y3D的表达式为：
二维高光谱图像Y的表达式为：
N＝L×W
其中，B为待解混三维高光谱图像的光谱通道个数，L为待解混三维高光谱图像的空间维长；W为待解混三维高光谱图像的空间维宽，N为二维高光谱图像的像素总个数。
步骤2、对二维高光谱图像进行端元提取，建立端元矩阵E。
具体的，包括以下步骤：
步骤21、采用最小误差识别方法，对二维高光谱图像Y进行最小化投影；利用最小化投影的噪声能量和投影信号误差能量之和，在正交子空间中，确定二维高光谱图像的端元数量R。
步骤22、根据二维高光谱图像Y的端元数量R，采用向量成分分析端元提取方法，进行端元提取，得到二维高光谱图像Y的R个端元，并利用二维高光谱图像Y的R个端元，构建端元矩阵E。
其中，采用基于几何理论的向量成分分析端元提取方法，进行端元提取，具体操作为：
初始化正交子空间，将二维高光谱图像在正交子空间中进行投影；
计算二维高光谱图像中每个像元在正交子空间的投影距离；
将二维高光谱图像中在正交子空间的投影距离最大的像元，作为候选端元，将候选端元加入到端元集合中；
更新正交子空间，将二维高光谱图像进行更新后的正交子空间中进行投影，并计算二维高光谱图像中每个像元在更新后的正交子空间的投影距离；
将二维高光谱图像中在更新后的正交子空间的投影距离最大的像元，作为新的候选端元，并将新的候选端元加入到已有端元集合中；
循环更新正交子空间，至端元集合中，包含R个候选端元，得到端元矩阵E。
步骤3、利用端元矩阵E，构建损失函数；其中，损失函数用于求解二维高光谱图像的丰度，损失函数包含先验正则项；
其中，损失函数的表达式为：
由于丰度表示二维高光谱图像中对应像元形成过程中每个端元所占的面积分数，丰度满足“和为一”及“非负”两个约束；即损失函数的约束为：
其中，Φ(*)为先验正则项，λ为正则项强度；ai为二维高光谱图像中第i个像素的丰度；yi为二维高光谱图像中第i个像素的光谱；i为二维高光谱图像的像素个数；A为丰度矩阵。
步骤4、求解损失函数，输出二维高光谱图像数据矩阵Y的丰度，二维高光谱图像数据矩阵的丰度即为待解混高光谱图像的解混结果；具体包括以下步骤：
步骤41、构造求解方法
引入辅助变量矩阵Z，对损失函数进行等价形式转化，得到转化后的损失函数；其中，转化后的损失函数的表达式为：
s.t.Z＝EA
步骤42、获取转化后损失函数的增广拉格朗日函数Lρ：
其中，V为对偶变量矩阵；||EA-Z||F为EA-Z矩阵的F范数；ρ为惩罚因子。
步骤43、利用交替方向乘子法方法，求解转化后的损失函数的增广拉格朗日函数Lρ；通过分解问题进行交替求解每个未知变量，得到待解混高光谱图像的解混结果。
具体包括以下步骤：
利用交替方向乘子法，将转化后的损失函数的增广拉格朗日函数分解为第一子问题和第二子问题；其中，第一子问题为增广拉格朗日函数Lρ中包含丰度a的项，第二子问题是增广拉格朗日函数Lρ中包含辅助变量Z的项。
其中，第一子问题为：
将转化后的损失函数的增广拉格朗日函数Lρ中包含丰度a的项整合并展开，获得求解丰度的优化问题；其中，ak+1,i为第i个像素在第k+1次迭代的丰度值，ρk为第k次迭代中的惩罚因子，
为第i个像素的中间变量，zk,i为第i个像素的辅助变量，uk,i为第i个像素的对偶变量；
第二子问题为：
其中，
为整幅二维高光谱图像第k次迭代中的辅助变量；整幅二维高光谱图像第k次迭代中的辅助变量
采用将整幅二维高光谱图像的所有像素的辅助变量拼接整合得到。
对二维高光谱图像中的每个像素，利用FCLS方法求解第一子问题，得到每个像素对应的丰度，并将所有像素的丰度拼接整合，获得整幅二维高光谱图像在第k+1次迭代的丰度Ak+1；利用降噪器求解第二子问题，获得整幅二维高光谱图像第k+1次迭代中的辅助变量。
本发明中，采用第一子问题求解的整幅二维高光谱图像在第k+1次迭代的丰度Ak+1，更新整幅二维高光谱图像第k次迭代中的辅助变量
用于求解第二子问题。
整幅二维高光谱图像第k次迭代中的辅助变量
的更新公式为：
其中，Uk为整幅二维高光谱图像在第k次迭代的对偶变量；
将整幅二维高光谱图像第k次迭代中的辅助变量
转换为三维矩阵Zk，以获得高光谱图像的空间和谱间信息，通过将将二维图像转换为三维，用于获取图像的空间和谱间先验知识：其中，转换后的三维矩阵Zk的表达式为：
第二子问题可以看作是一个对辅助变量Z的降噪过程，使用降噪器来求解第二问题，获得整幅二维高光谱图像第k+1次迭代的辅助变量Zk+1。
利用对偶变量更新公式，更新计算求解第一子问题需要的二维高光谱图像第k+1次迭代对偶变量Uk+1；其中，对偶变量Uk+1的更新公式为：
Uk+1＝Uk+EAk+1-Zk+1
利用惩罚因子更新公式，更新计算求解第一子问题需要的二维高光谱图像第k+1次迭代惩罚因子ρk+1，其中，惩罚因子ρk+1的更新公式为
ρk+1＝αρk
其中，α为扩大系数。
循环迭代，求解第一子问题和第二子问题，至迭代循环次数达到最大循环次数K，得到整幅二维高光谱图像在第K次迭代的丰度AK，输出高光谱图像的解混结果。
优选的，降噪器采用二维图像降噪器NLM、二维图像降噪器BM3D或三维图像降噪器BM4D。
步骤5、利用所述解混结果对待解混图像进行高光谱数据分析。
本实施例还提供了一种高光谱图像的解混系统，包括获取模块、端元矩阵模块、函数模块及求解模块；获取模块，获取待解混三维高光谱图像Y3D，对待解混三维高光谱图像Y3D进行重构处理，得到二维高光谱图像Y；端元矩阵模块，用于对二维高光谱图像Y进行端元提取，建立端元矩阵E；函数模块，用于利用端元矩阵E，构建损失函数；其中，损失函数用于求解二维高光谱图像Y的丰度，损失函数包含先验正则项；求解模块，用于求解损失函数，输出二维高光谱图像Y的丰度，二维高光谱图像Y的丰度即为待解混三维高光谱图像的解混结果。
实施例
本实施例提供了一种高光谱图像的解混方法及系统，具体包括以下步骤：
步骤1、对目标场景进行成像，得到待解混的高光谱图像
并将三维高光谱图像经过矩阵变换，重构处理后得到二维高光谱图像
N＝L×W；
其中，B为待解混三维高光谱图像的光谱通道个数，L为待解混三维高光谱图像的空间维长；W为待解混三维高光谱图像的空间维宽，N为二维高光谱图像的像素总个数。
步骤2、对二维高光谱图像Y，采用高光谱图像的最小误差识别方法，确定高光谱图像的端元个数R；本实施例通过最小化投影后的噪声能量和投影信号误差能量之和，在生成子空间中，确定二维高光谱图像的端元个数R；
步骤3、根据步骤2中得到的二维高光谱图像的端元个数R，采用基于几何理论的向量成分分析端元提取方法，对二维高光谱图像进行端元提取；其中，基于几何理论的向量成分分析端元提取方法通过找到二维高光谱图像中单形体的投影距离最大的像元，并将二维高光谱图像中在正交子空间的投影距离最大的像元作为候选的端元；更新正交子空间，在每次循环中寻找与已知端元正交的投影矩阵，通过寻找正交子空间并计算各像元在该正交子空间中的投影距离，找到投影距离最大的像元作为新获取的端元，加入已有的端元集合当中，然后再寻找正交子空间，再投影，通过多次循环来提取端元，直到找到R个端元E，即得到端元矩阵E。
步骤4、利用端元矩阵E，构建损失函数；其中，损失函数用于求解二维高光谱图像Y的丰度，损失函数包含先验正则项；
对于二维高光谱图像Y，构建求解二维高光谱图像丰度的损失函数：
其中，Φ(*)为先验正则项，λ为正则项强度；ai为二维高光谱图像中第i个像素的丰度；yi为二维高光谱图像中第i个像素的光谱；i为二维高光谱图像的像素个数；A为丰度矩阵。
通过最小化上式，估计该像素的丰度ai；由于丰度表示该像元形成过程中每个端元所占的面积分数，其满足“和为一”及“非负”，即该损失函数的约束为：
步骤5、构造步骤4中损失函数的求解方法；具体的，引入变量矩阵Z，将损失函数转化为一个等价形式，得到转化后的损失函数；
其中，转化后的损失函数的表达式为：
s.t.Z＝EA
步骤6、获取转化后损失函数的增广拉格朗日函数Lρ：
其中，V为对偶变量矩阵；||EA-Z||F为EA-Z矩阵的F范数；ρ为惩罚因子。
V是对偶变量，|| ||F表示矩阵的F范数，ρ表示惩罚因子。
步骤7、利用交替方向乘子法方法，求解转化后的损失函数的增广拉格朗日函数Lρ；通过分解问题进行交替求解每个未知变量，得到待解混高光谱图像的解混结果。
具体包括以下步骤：
利用交替方向乘子法，将转化后的损失函数的增广拉格朗日函数分解为第一子问题和第二子问题；其中，第一子问题为增广拉格朗日函数Lρ中包含丰度a的项，第二子问题是增广拉格朗日函数Lρ中包含辅助变量Z的项。
其中，第一子问题为：
将转化后的损失函数的增广拉格朗日函数Lρ中包含丰度a的项整合并展开，获得求解丰度的优化问题；其中，ak+1,i为第i个像素在第k+1次迭代的丰度值，ρk为第k次迭代中的惩罚因子，
为第i个像素的中间变量，zk,i为第i个像素的辅助变量，uk,i为第i个像素的对偶变量；
第二子问题为：
其中，
为整幅二维高光谱图像第k次迭代中的辅助变量；整幅二维高光谱图像第k次迭代中的辅助变量
采用将整幅二维高光谱图像的所有像素的辅助变量拼接整合得到。
对二维高光谱图像中的每个像素，利用FCLS方法求解第一子问题，得到每个像素对应的丰度，并将所有像素的丰度拼接整合，获得整幅二维高光谱图像在第k+1次迭代的丰度Ak+1；利用降噪器求解第二子问题，获得整幅二维高光谱图像第k+1次迭代中的辅助变量。
本发明中，采用第一子问题求解的整幅二维高光谱图像在第k+1次迭代的丰度Ak+1，更新整幅二维高光谱图像第k次迭代中的辅助变量
用于求解第二子问题。
整幅二维高光谱图像第k次迭代中的辅助变量
的更新公式为：
其中，Uk为整幅二维高光谱图像在第k次迭代的对偶变量。
将整幅二维高光谱图像第k次迭代中的辅助变量
转换为三维矩阵Zk，以获得高光谱图像的空间和谱间信息，通过将将二维图像转换为三维，用于获取图像的空间和谱间先验知识：其中，转换后的三维矩阵Zk的表达式为：
第二子问题可以看作是一个对辅助变量Z的降噪过程，使用降噪器来求解第二问题，获得整幅二维高光谱图像第k+1次迭代的辅助变量Zk+1。
利用对偶变量更新公式，更新计算求解第一子问题需要的二维高光谱图像第k+1次迭代对偶变量Uk+1；其中，对偶变量Uk+1的更新公式为：
Uk+1＝Uk+EAk+1-Zk+1
利用惩罚因子更新公式，更新计算求解第一子问题需要的二维高光谱图像第k+1次迭代惩罚因子ρk+1，其中，惩罚因子ρk+1的更新公式为
ρk+1＝αρk
其中，α为扩大系数。
循环迭代，求解第一子问题和第二子问题，至迭代循环次数达到最大循环次数K，得到整幅二维高光谱图像在第K次迭代的丰度AK，输出高光谱图像的解混结果。
优选的，降噪器采用二维图像降噪器NLM、二维图像降噪器BM3D或三维图像降噪器BM4D。
步骤8、利用所述解混结果对待解混图像进行高光谱数据分析。
本实施例中，采用二维图像降噪器NLM或二维图像降噪器BM3D时，λ＝1×10-4，ρ＝4×10-4；K＝18，α＝1.2；采用三维图像降噪器BM4D时，λ＝5×10-4，ρ＝4×10-3；K＝18，α＝1.2。
本实施例中还提供了一种高光谱图像的解混系统，包括获取模块、端元矩阵模块、函数模块及求解模块；获取模块，获取待解混三维高光谱图像Y3D，对待解混三维高光谱图像Y3D进行重构处理，得到二维高光谱图像Y；端元矩阵模块，用于对二维高光谱图像Y进行端元提取，建立端元矩阵E；函数模块，用于利用端元矩阵E，构建损失函数；其中，损失函数用于求解二维高光谱图像Y的丰度，损失函数包含先验正则项；求解模块，用于求解损失函数，输出二维高光谱图像Y的丰度，二维高光谱图像Y的丰度即为待解混三维高光谱图像的解混结果。
具体试验
采用本发明所述的解混方法及系统，对一组模拟数据和Urban数据进行试验；本试验中，采用二维图像降噪器NLM、二维图像降噪器BM3D及三维图像降噪器BM4D，来验证本发明的灵活性和有效性。
在模拟数据集实验中，从美国地质调查局光谱库中随机选取4个端元，共包含224个光谱段；使用HYDRA工具包生成包含空间信息的丰度，生成高光谱图像的空间维大小为256×256；基于线性混合模型生成模拟数据，加入信噪比为5、10、20及30dB的零均值高斯白噪声。
Urban高光谱图像的空间维为307×307，由210个光谱波段组成，光谱范围从400nm到2500nm，光谱分辨率高达10nm；去除被密集的水蒸气和大气影响的通道(1-4，76，87，101-111，136-153，198-210)后，实验中使用162个通道的光谱信息进行实验。
为了证明方法的有效性，使用3种方法进行对比：FCLS、SUnSAL-TV和SCHU。
本实验采用根值均方误差来评价高光谱图像的解混性能，得到估计模拟丰度的根值均方误差对比结果，如下表1所示。
表1本发明与现有技术估计模拟数据丰度的根值均方误差结果对比表
由表1可以得到，本实施例提出的解混算法取得了最低的RMSE结果；本发明通过将一个优化子问题转化为使用降噪器对高光谱图像进行降噪来避免求解优化问题，且可以通过降噪器学习高光谱图像的空间和谱间先验信息，通过对先验信息的使用提高了高光谱图像解混的精度。
如附图2所示，附图2给出了不同信噪比情况下，原始高光谱图像和使用二维图像降噪器NLM解混后重构的高光谱图像的第190个光谱通道的映射图；从附图2中可以看出，得益于二维图像降噪器NLM的使用，重构高光谱图像拥有较少的噪声，图像更加平滑，进一步说明了本实施例能够充分使用图像空间和谱间信息，且解混结果对噪声鲁棒。
如附图3所示，附图3给出了本实施例与现有技术Urban数据集丰度估计映射图，从附图3中可以看出，本实施例具有更好的空间结构学习能力和空间结构保持能力，本发明提出的方法很多图像边缘和细节都更突出。
本发明所述的解混方法系统，和现有方法相比，使用交替方向乘子法将高光谱图像解混问题分解为两个子问题，其中第一子问题为标准的全约束最小二乘问题，使用FCLS方法求解，另一子问题可以看作对高光谱图像降噪，采用降噪器来求解高光谱图像降噪子问题，通过降噪器可以自动学习高光谱数据的空间和谱间信息；其中，使用二维降噪器学习图像的空间结构信息，使用三维的降噪器可以同时学习图像的空间和谱间信息，此外，降噪器的使用使本发明对图像的噪声鲁棒；本发明可以灵活使用各种降噪器求解学习高光谱图像的先验知识，避免了手动选择正则项，最大程度地使用了高光谱图像的空间和谱间信息，提高了解混的精度。
上述实施例仅仅是能够实现本发明技术方案的实施方式之一，本发明所要求保护的范围并不仅仅受本实施例的限制，还包括在本发明所公开的技术范围内，任何熟悉本技术领域的技术人员所容易想到的变化、替换及其他实施方式。
Claims (10)
1.一种高光谱图像的解混方法，其特征在于，包括以下步骤：
获取待解混三维高光谱图像Y3D，对待解混三维高光谱图像Y3D进行重构处理，得到二维高光谱图像Y；
对二维高光谱图像Y进行端元提取，建立端元矩阵E；
利用端元矩阵E，构建损失函数；其中，损失函数用于求解二维高光谱图像Y的丰度，损失函数包含先验正则项；
求解损失函数，输出二维高光谱图像Y的丰度，二维高光谱图像Y的丰度即为待解混三维高光谱图像的解混结果。
2.根据权利要求1所述的一种高光谱图像的解混方法，其特征在于，对目标场景进行成像，得到待解混三维高光谱图像Y3D；对待解混三维高光谱图像Y3D，进行矩阵变换，得到二维高光谱图像数据矩阵Y；
其中，待解混三维高光谱图像Y3D的表达式为：
二维高光谱图像Y的表达式为：
N＝L×W
其中，B为待解混三维高光谱图像的光谱通道个数，L为待解混三维高光谱图像的空间维长；W为待解混三维高光谱图像的空间维宽，N为二维高光谱图像的像素总个数。
3.根据权利要求1所述的一种高光谱图像的解混方法，其特征在于，建立端元矩阵E，具体为：
确定二维高光谱图像Y的端元数量R，采用向量成分分析端元提取方法，进行端元提取，得到二维高光谱图像Y的R个端元，并利用二维高光谱图像Y的R个端元，构建端元矩阵E。
4.根据权利要求3所述的一种高光谱图像的解混方法，其特征在于，确定二维高光谱图像数据矩阵的端元数量R，具体操作如下：
采用最小误差识别方法，对二维高光谱图像Y进行最小化投影；利用最小化投影后的噪声能量和投影信号误差能量之和，在正交子空间中，确定二维高光谱图像的端元数量。
5.根据权利要求3所述的一种高光谱图像的解混方法，其特征在于，采用基于几何理论的向量成分分析端元提取方法，进行端元提取；具体操作为：
初始化正交子空间，将二维高光谱图像在正交子空间中进行投影；
计算二维高光谱图像中每个像元在正交子空间的投影距离；
将二维高光谱图像中在正交子空间的投影距离最大的像元，作为候选端元，将候选端元加入到端元集合中；
更新正交子空间，将二维高光谱图像在更新后的正交子空间中进行投影，并计算二维高光谱图像中每个像元在更新后的正交子空间的投影距离；
将二维高光谱图像中在更新后的正交子空间的投影距离最大的像元，作为新的候选端元，并将新的候选端元加入到已有端元集合中；
循环更新正交子空间，至端元集合中，包含R个候选端元，得到端元矩阵E。
6.根据权利要求1所述的一种高光谱图像的解混方法，其特征在于，损失函数的表达式为：
其中，Φ(*)为先验正则项，λ为正则项强度；ai为二维高光谱图像中第i个像素的丰度；yi为二维高光谱图像中第i个像素的光谱；i为二维高光谱图像的像素个数；A为丰度矩阵。
7.根据权利要求6所述的一种高光谱图像的解混方法，其特征在于，求解损失函数，具体包括以下步骤：
引入辅助变量矩阵Z，对损失函数进行等价形式转化，得到转化后的损失函数；
获取转化后的损失函数的增广拉格朗日函数Lρ；
利用交替方向乘子法，求解转化后的损失函数的增广拉格朗日函数Lρ，得到二维高光谱图像数据矩阵的丰度，即得到待解混高光谱图像的解混结果。
8.根据权利要求7所述的一种高光谱图像的解混方法，其特征在于，转化后的损失函数的增广拉格朗日函数Lρ的表达式为：
其中，V为对偶变量矩阵；||EA-Z||F为EA-Z矩阵的F范数；ρ为惩罚因子。
9.根据权利要求7所述的一种高光谱图像的解混方法，其特征在于，利用交替方向乘子法，求解转化后的损失函数的增广拉格朗日函数Lρ，具体操作如下：
利用交替方向乘子法，将转化后的损失函数的增广拉格朗日函数分解为第一子问题和第二子问题；
其中，第一子问题为：
其中，ak+1,i为第i个像素在第k+1次迭代的丰度值，ρk为第k次迭代中的惩罚因子，
为第i个像素的中间变量，zk,i为第i个像素的辅助变量，uk,i为第i个像素的对偶变量；
第二子问题为：
其中，
为整幅二维高光谱图像第k次迭代中的辅助变量；
对二维高光谱图像中的每个像素，利用FCLS方法求解第一子问题，得到每个像素对应的丰度，并将所有像素的丰度拼接整合，获得整幅二维高光谱图像在第k+1次迭代的丰度Ak+1；利用降噪器求解第二子问题，获得整幅二维高光谱图像第k+1次迭代中的辅助变量；
循环迭代，求解第一子问题和第二子问题，至迭代循环次数达到最大循环次数K，输出高光谱图像的解混结果。
10.一种高光谱图像的解混系统，其特征在于，包括获取模块、端元矩阵模块、函数模块及求解模块；
获取模块，获取待解混三维高光谱图像Y3D，对待解混三维高光谱图像Y3D进行重构处理，得到二维高光谱图像Y；
端元矩阵模块，用于对二维高光谱图像Y进行端元提取，建立端元矩阵E；
函数模块，用于利用端元矩阵E，构建损失函数；其中，损失函数用于求解二维高光谱图像Y的丰度，损失函数包含先验正则项；
求解模块，用于求解损失函数，输出二维高光谱图像Y的丰度，二维高光谱图像Y的丰度即为待解混三维高光谱图像的解混结果。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202011630071.1A CN112712034B (zh) | 2020-12-30 | 2020-12-30 | 一种高光谱图像的解混方法及系统 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202011630071.1A CN112712034B (zh) | 2020-12-30 | 2020-12-30 | 一种高光谱图像的解混方法及系统 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN112712034A true CN112712034A (zh) | 2021-04-27 |
| CN112712034B CN112712034B (zh) | 2024-06-21 |
Family
ID=75547756
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202011630071.1A Active CN112712034B (zh) | 2020-12-30 | 2020-12-30 | 一种高光谱图像的解混方法及系统 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN112712034B (zh) |
Cited By (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113344049A (zh) * | 2021-05-27 | 2021-09-03 | 湖州师范学院 | 一种基于Sinkhorn距离的盲高光谱解混模型的构建方法 |
| CN113378664A (zh) * | 2021-05-26 | 2021-09-10 | 辽宁工程技术大学 | 基于半非负矩阵分解的高光谱图像变化检测方法和系统 |
| CN113673300A (zh) * | 2021-06-24 | 2021-11-19 | 核工业北京地质研究院 | 一种基于非监督训练的高光谱图像智能解混方法 |
| CN115272093A (zh) * | 2022-04-22 | 2022-11-01 | 哈尔滨师范大学 | 一种基于空间结构信息约束的高光谱图像解混方法 |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104952050A (zh) * | 2015-07-07 | 2015-09-30 | 西安电子科技大学 | 基于区域分割的高光谱图像自适应解混方法 |
| US20190180429A1 (en) * | 2016-08-26 | 2019-06-13 | Nec Corporation | Image processing device, image processing method, and computer-readable recording medium |
| CN110428454A (zh) * | 2019-08-13 | 2019-11-08 | 电子科技大学中山学院 | 一种高光谱解混方法、装置、电子设备及存储介质 |
-
2020
- 2020-12-30 CN CN202011630071.1A patent/CN112712034B/zh active Active
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104952050A (zh) * | 2015-07-07 | 2015-09-30 | 西安电子科技大学 | 基于区域分割的高光谱图像自适应解混方法 |
| US20190180429A1 (en) * | 2016-08-26 | 2019-06-13 | Nec Corporation | Image processing device, image processing method, and computer-readable recording medium |
| CN110428454A (zh) * | 2019-08-13 | 2019-11-08 | 电子科技大学中山学院 | 一种高光谱解混方法、装置、电子设备及存储介质 |
Non-Patent Citations (3)
| Title |
|---|
| FENGCHAO XIONG 等: "Hyperspectral Unmixing via Total Variation Regularized Nonnegative Tensor Factorization", 《IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING》, vol. 57, no. 04, 30 April 2019 (2019-04-30), pages 2341 - 2357, XP011716758, DOI: 10.1109/TGRS.2018.2872888 * |
| SMARTWEED: "Learning Deep CNN Denoiser Prior for Image Restoration阅读笔记", Retrieved from the Internet <URL:www.cnblogs.com/smartweed/p/10444039.html> * |
| 李宏俏: "高光谱图像的光谱解混模型与算法研究", 《中国优秀硕士学位论文全文数据库 信息科技辑》, no. 2018, 15 February 2018 (2018-02-15), pages 140 - 1437 * |
Cited By (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113378664A (zh) * | 2021-05-26 | 2021-09-10 | 辽宁工程技术大学 | 基于半非负矩阵分解的高光谱图像变化检测方法和系统 |
| CN113344049A (zh) * | 2021-05-27 | 2021-09-03 | 湖州师范学院 | 一种基于Sinkhorn距离的盲高光谱解混模型的构建方法 |
| CN113344049B (zh) * | 2021-05-27 | 2022-07-26 | 湖州师范学院 | 一种基于Sinkhorn距离的盲高光谱解混模型的构建方法 |
| CN113673300A (zh) * | 2021-06-24 | 2021-11-19 | 核工业北京地质研究院 | 一种基于非监督训练的高光谱图像智能解混方法 |
| CN115272093A (zh) * | 2022-04-22 | 2022-11-01 | 哈尔滨师范大学 | 一种基于空间结构信息约束的高光谱图像解混方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN112712034B (zh) | 2024-06-21 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN112712034B (zh) | 一种高光谱图像的解混方法及系统 | |
| US9317929B2 (en) | Decomposition apparatus and method for refining composition of mixed pixels in remote sensing images | |
| CN111260576B (zh) | 一种基于去噪三维卷积自编码网络的高光谱解混算法 | |
| Zhong et al. | Non-local sparse unmixing for hyperspectral remote sensing imagery | |
| Guo et al. | Hyperspectral image unmixing using autoencoder cascade | |
| Iordache et al. | Sparse unmixing of hyperspectral data | |
| CN109241843B (zh) | 空谱联合多约束优化非负矩阵解混方法 | |
| US8139650B2 (en) | Fast noise reduction in digital images and video | |
| CN108427934B (zh) | 一种高光谱影像混合像元分解方法 | |
| CN111008975B (zh) | 一种空间人造目标线性模型的混合像元解混方法及系统 | |
| CN111680579B (zh) | 一种自适应权重多视角度量学习的遥感图像分类方法 | |
| CN114529830A (zh) | 基于混合卷积网络遥感图像时空融合方法 | |
| CN112750091A (zh) | 一种高光谱图像解混方法 | |
| CN109727280B (zh) | 一种基于正交基的高光谱图像丰度估计方法 | |
| CN105957112A (zh) | 基于快速uncls的高光谱亚像素探测方法 | |
| CN113446998B (zh) | 一种基于高光谱目标探测数据的动态解混方法 | |
| CN113421198A (zh) | 一种基于子空间的非局部低秩张量分解的高光谱图像去噪方法 | |
| Wang et al. | Hyperspectral unmixing via plug-and-play priors | |
| Liu et al. | Comparative analysis of pixel level fusion algorithms in high resolution SAR and optical image fusion | |
| Ince et al. | Simultaneous nonconvex denoising and unmixing for hyperspectral imaging | |
| CN110956221A (zh) | 基于深度递归网络下小样本极化合成孔径雷达图像分类方法 | |
| Shah et al. | Estimating sparse signals with smooth support via convex programming and block sparsity | |
| CN109447951B (zh) | 基于吉文斯旋转的高光谱图像端元提取方法 | |
| CN112819909A (zh) | 基于低分辨率先验光谱图像区域分割的自适应编码方法 | |
| CN114332624A (zh) | 一种基于非负矩阵分解的即插即用高光谱图像解混方法 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |