CN109035154A - 一种考虑光谱变异性的高光谱图像非线性解混方法 - Google Patents
一种考虑光谱变异性的高光谱图像非线性解混方法 Download PDFInfo
- Publication number
- CN109035154A CN109035154A CN201810579080.9A CN201810579080A CN109035154A CN 109035154 A CN109035154 A CN 109035154A CN 201810579080 A CN201810579080 A CN 201810579080A CN 109035154 A CN109035154 A CN 109035154A
- Authority
- CN
- China
- Prior art keywords
- matrix
- spectrum
- abundance
- coefficient
- variable
- 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
- 238000001228 spectrum Methods 0.000 title claims abstract description 73
- 238000000034 method Methods 0.000 title claims abstract description 60
- 238000002156 mixing Methods 0.000 title claims abstract description 23
- 230000003595 spectral effect Effects 0.000 claims abstract description 17
- VMXUWOKSQNHOCA-UKTHLTGXSA-N ranitidine Chemical compound [O-][N+](=O)\C=C(/NC)NCCSCC1=CC=C(CN(C)C)O1 VMXUWOKSQNHOCA-UKTHLTGXSA-N 0.000 claims abstract description 11
- 238000009499 grossing Methods 0.000 claims abstract description 7
- 239000011159 matrix material Substances 0.000 claims description 65
- 239000013598 vector Substances 0.000 claims description 12
- 238000013507 mapping Methods 0.000 claims description 8
- 230000009977 dual effect Effects 0.000 claims description 6
- 230000003416 augmentation Effects 0.000 claims description 4
- 230000035772 mutation Effects 0.000 claims description 3
- 230000003247 decreasing effect Effects 0.000 claims description 2
- 238000005457 optimization Methods 0.000 claims description 2
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 claims 1
- 238000012545 processing Methods 0.000 abstract description 2
- 230000006870 function Effects 0.000 description 20
- 230000000694 effects Effects 0.000 description 11
- 239000000126 substance Substances 0.000 description 10
- BERDEBHAJNAUOM-UHFFFAOYSA-N copper(I) oxide Inorganic materials [Cu]O[Cu] BERDEBHAJNAUOM-UHFFFAOYSA-N 0.000 description 7
- LBJNMUFDOHXDFG-UHFFFAOYSA-N copper;hydrate Chemical compound O.[Cu].[Cu] LBJNMUFDOHXDFG-UHFFFAOYSA-N 0.000 description 7
- 238000002310 reflectometry Methods 0.000 description 7
- 239000000203 mixture Substances 0.000 description 6
- 238000002474 experimental method Methods 0.000 description 5
- 238000012360 testing method Methods 0.000 description 5
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 5
- YGANSGVIUGARFR-UHFFFAOYSA-N dipotassium dioxosilane oxo(oxoalumanyloxy)alumane oxygen(2-) Chemical compound [O--].[K+].[K+].O=[Si]=O.O=[Al]O[Al]=O YGANSGVIUGARFR-UHFFFAOYSA-N 0.000 description 4
- 238000005516 engineering process Methods 0.000 description 4
- 238000000605 extraction Methods 0.000 description 4
- NLYAJNPCOHFWQQ-UHFFFAOYSA-N kaolin Chemical compound O.O.O=[Al]O[Si](=O)O[Si](=O)O[Al]=O NLYAJNPCOHFWQQ-UHFFFAOYSA-N 0.000 description 4
- 229910052622 kaolinite Inorganic materials 0.000 description 4
- 229910052627 muscovite Inorganic materials 0.000 description 4
- 229910052934 alunite Inorganic materials 0.000 description 3
- 239000010424 alunite Substances 0.000 description 3
- 238000004458 analytical method Methods 0.000 description 3
- 239000002245 particle Substances 0.000 description 3
- 238000011160 research Methods 0.000 description 3
- KPZTWMNLAFDTGF-UHFFFAOYSA-D trialuminum;potassium;hexahydroxide;disulfate Chemical compound [OH-].[OH-].[OH-].[OH-].[OH-].[OH-].[Al+3].[Al+3].[Al+3].[K+].[O-]S([O-])(=O)=O.[O-]S([O-])(=O)=O KPZTWMNLAFDTGF-UHFFFAOYSA-D 0.000 description 3
- 102000008297 Nuclear Matrix-Associated Proteins Human genes 0.000 description 2
- 108010035916 Nuclear Matrix-Associated Proteins Proteins 0.000 description 2
- 238000010521 absorption reaction Methods 0.000 description 2
- 229910052833 almandine Inorganic materials 0.000 description 2
- 238000009472 formulation Methods 0.000 description 2
- 238000003384 imaging method Methods 0.000 description 2
- 229910000273 nontronite Inorganic materials 0.000 description 2
- 210000000299 nuclear matrix Anatomy 0.000 description 2
- 230000005855 radiation Effects 0.000 description 2
- 239000002689 soil Substances 0.000 description 2
- 238000012546 transfer Methods 0.000 description 2
- 241000183024 Populus tremula Species 0.000 description 1
- 230000009471 action Effects 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 230000000052 comparative effect Effects 0.000 description 1
- 238000000354 decomposition reaction Methods 0.000 description 1
- 230000007423 decrease Effects 0.000 description 1
- 238000009795 derivation Methods 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 230000018109 developmental process Effects 0.000 description 1
- 235000013399 edible fruits Nutrition 0.000 description 1
- 230000007613 environmental effect Effects 0.000 description 1
- 238000011156 evaluation Methods 0.000 description 1
- 239000004744 fabric Substances 0.000 description 1
- 238000005286 illumination Methods 0.000 description 1
- 238000010191 image analysis Methods 0.000 description 1
- 238000003709 image segmentation Methods 0.000 description 1
- 238000003331 infrared imaging Methods 0.000 description 1
- 229910052500 inorganic mineral Inorganic materials 0.000 description 1
- 230000003993 interaction Effects 0.000 description 1
- 238000012804 iterative process Methods 0.000 description 1
- 239000000463 material Substances 0.000 description 1
- 230000007246 mechanism Effects 0.000 description 1
- 239000011707 mineral Substances 0.000 description 1
- 238000005065 mining Methods 0.000 description 1
- 238000012544 monitoring process Methods 0.000 description 1
- 230000008569 process Effects 0.000 description 1
- 239000004576 sand Substances 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
- 238000004611 spectroscopical analysis Methods 0.000 description 1
- 239000004575 stone Substances 0.000 description 1
Classifications
-
- 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
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T5/00—Image enhancement or restoration
- G06T5/50—Image enhancement or restoration using two or more images, e.g. averaging or subtraction
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/10—Image acquisition modality
- G06T2207/10032—Satellite or aerial image; Remote sensing
- G06T2207/10036—Multispectral image; Hyperspectral image
Landscapes
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- Theoretical Computer Science (AREA)
- Image Processing (AREA)
- Investigating Or Analysing Materials By Optical Means (AREA)
Abstract
本发明属于遥感图像处理技术领域，具体为一种考虑光谱变异性的高光谱图像非线性解混方法。本发明首先利用核方法将原始数据映射到高维特征空间，在高维空间中考虑光谱变异系数进行线性解混；同时，依据地物分布的空间连续性，对丰度和变异系数添加局部平滑约束，使得二者具有空间上的平滑性。本方法在Hapke和GBM两种非线性混合模型中存在光谱变异性时，能进行有效的无监督非线性光谱解混。本发明能克服不同非线性混合场景中存在的光谱变异性问题，提高光谱解混的精度，在实际应用中具有重要的意义。
Description
技术领域
本发明属于遥感图像处理技术领域，具体涉及一种高光谱图像非线性解混方法。
背景技术
遥感技术是本世纪六十年代发展起来的新兴综合技术，与空间、电子光学、计算机、地理学等科学技术紧密相关，是研究地球资源环境的最有力的技术手段之一。高光谱遥感是将成像技术与光谱技术相结合的多维信息获取技术，其图像通常包含数百个波段，具有较高的光谱分辨率，可以在获取地物空间分布的同时得到丰富的光谱信息，被广泛应用在军事侦察、环境监测和地质勘探等诸多领域[1]。但是，由于传感器较低的空间分辨率以及地面物质构成的复杂性，导致混合像元大量存在于高光谱图像中，严重阻碍了高光谱图像的高精度应用。光谱解混可以将高光谱图像中的混合像元分解为纯物质光谱(即端元)和端元在混合像元中所占的比例(即丰度)，从而将图像分析深入到亚像元级，推动了定量遥感的发展。
线性混合模型(Linear Mixing Model,LMM)是在以往的研究中使用比较多的一种模型[2]。LMM模型简单、物理意义明确，它基于光子只与地面场景中的一种物质发生作用的假设，因而，像元光谱被表示为由端元光谱以一定的比例系数线性组合而成。在LMM模型中，所有像元由一组共同的端元光谱表示，但是在实际情况中，由于光照、采集条件和物质的固有性质等因素的影响，端元光谱在像元中会发生改变，产生光谱变异性问题[3]。而忽略光谱变异性的存在会给解混结果带来较大的误差，为此学者针对LMM中的光谱变异性问题开展了许多研究。但是LMM的假设在真实场景中并不能总被满足，特别是对于非线性混合效应显著的复杂场景情况[4]。因此，在开展非线性光谱解混研究的同时，考虑场景中的光谱变异性，在实际应用中具有重要的意义。
通常，非线性混合场景中主要关注两种非线性混合模型。第一种是Hapke等模型，其主要针对沙地、矿物等地区存在的紧密混合现象。Hapke模型[5]依据辐射传输理论，把双向反射率表示为与场景相关的粒子密度和大小以及单次散射反照率等物理参数的非线性函数，从而描述光线在微观尺度上与不同粒子之间的多次散射。第二种是广义双线性模型(Generalized Bilinear Model,GBM)[6]，一般被用来描述植被覆盖、建筑林立的地区存在的多层次混合现象。GBM在LMM的基础上添加了两两端元间存在的二次散射效应并忽略更高次散射的影响。尽管这两种非线性混合模型在混合机理和表达形式上都有很大的不同，但是，核方法可以利用合理的核函数形式，将低维空间的非线性混合的原始数据投射到高维特征空间中，从而将低维空间的非线性解混问题转变为高维空间中的线性解混问题[7]。因为它不需要考虑具体的非线性形式，所以可适用于不同的非线性模型。Zhu等人提出了Bi-objective NMF算法[8]，同时考虑了原始空间和核方法映射后的高维特征空间中的数据重构误差，并调整两者间的权重因子实现解混。文献[9]中提出了一种基于丰度约束核非负矩阵分解ASSKNMF的解混算法，该算法用核方法对数据进行非线性映射的同时，也根据地物分布特性在丰度上添加了稀疏和平滑约束来进行非线性解混。但是，以上所述算法都仅关注非线性解混问题，而没有考虑到非线性场景中存在的光谱变异性，这样仍会使实际的解混结果出现较大的偏差。
Hapke混合模型
对于高光谱数据X＝[x1,x2,…,xN]∈RL×N，其每列xn∈RL×1(n＝1,2,...N)都对应一个具有L个波段的像元向量，共存在N个这样的像元。以A＝[a1,a2,...,aP]∈RL×P表示端元矩阵(P是端元数)，spn是端元ap在像元xn中的丰度，en是模型误差。LMM假设下的像元是端元按丰度的线性组合，为满足物理意义，丰度需满足非负与“和为一”的约束条件[2]：
Hapke模型是依据辐射传输理论，同时假设粒子为球形且各向散射同性，且表面反照比低，光谱反射率x可以表示为：
其中，μ0,μ分别是光线入射角和出射角的余弦值，ω是单次散射反照率，H是多项散射函数，其表达式如下：
根据式(2)和(3)，可以在已知反射率x的情况下推导出单次散射反照率ω，R-1表示反函数：
根据Hapke模型[5]，像元的反射率x由(4)非线性映射为单次散射反照率ω后，便能被端元的单次散射反照率线性表示：
其中ωp是端元ap的单次散射反照率，sp,n是对应的丰度系数。对端元的单次散射反照率按丰度系数进行线性组合后，再经由非线性映射R(ω)变为高光谱数据中的反射率。
GBM混合模型
GBM在LMM的基础上增加了两两端元间的二次散射效应作为非线性项：
其中，ap⊙aq＝[ap,1aq,1,ap,2aq,2,...ap,Laq,L]T为端元之间的Hadamard乘积，用于描述光线在两种端元物质间发生的二次散射非线性交互作用，参数ζp,q,n用于控制非线性程度，使得GBM模型可以更加灵活地描述非线性混合效应。
核方法
考虑非线性映射：
将原始数据映射到高维空间中，使得数据在高维空间中变得线性可分。通常情况下，非线性映射φ的形式很复杂且难以求得，但可通过核函数的方法来避免这个问题[10]。两个向量在高维空间中的内积能用核函数的形式表示：
一般来说，满足Mercer定理或正定性的函数就可以作为核函数，本文采用RBF函数[10]：
其中σ是核参数。
发明内容
本发明的目的在于提出一种适用于不同非线性模型的考虑光谱变异性的高光谱非线性解混方法。
本发明提出的考虑光谱变异性的高光谱图像非线性解混算法，具体步骤为：首先利用核方法将原始数据映射到高维特征空间，在高维空间中考虑光谱变异系数进行线性解混。与此同时，依据地物分布的空间连续性，对丰度和变异系数添加局部平滑约束，使得二者具有空间上的平滑性。所提出的方法在Hapke和GBM两种非线性混合模型中存在光谱变异性时，能进行有效的无监督非线性光谱解混。具体内容介绍如下：
针对所提及的Hapke和GBM两类非线性混合模型，可将原始高光谱数据先用核方法映射到高维空间，然后在高维空间中使基础端元矩阵乘上对应的比例系数向量来描述光谱变异性，最后同时求解端元、丰度和光谱变异系数进行无监督的线性解混。这样既能不需要考虑具体的非线性混合形式，也可以顾及到光谱变异性的影响。因此，在核空间中，像元相应的表达形式为：
φ(xn)＝φ(A)diag(zn)sn＝φ(A)(zn⊙sn) (10)。
其中，φ(A)＝[φ(a1),φ(a2),...,φ(aP)]，结合实际地物的分布特性加入平滑约束后可以得到目标函数：
这里的第一项考虑了光谱变异性的重构误差；和ψ(Z)用于反映物质空间分布的连续性，是对丰度和变异系数添加的局部平滑约束，α和β分别是相应的惩罚系数；τC(S)是丰度满足“非负”与“和为一”约束的指示函数。为了方便目标函数的求解，我们引入了一个新的变量B，得到新的目标函数：
其中，约束参数λ衡量了所引入变量bn与zn⊙sn间的差异，像元光谱差异wn,m可以更好的体现中心像素和邻域像素在丰度和变异系数上的相似性。针对上述目标函数，对其中的端元矩阵A，变量矩阵B，丰度矩阵S和变异系数矩阵Z逐一交替优化进行求解。
算法初始化：端元矩阵A和丰度矩阵S分别采用VCA和FCLS进行初始化，变异系数矩阵Z每个元素均初始取1，变量矩阵B用A和S初始值的乘积初始化。
(1)优化端元矩阵A
保留式(12)中与端元矩阵A相关的项，将式(12)重新写成：
根据前面介绍的核方法，f(A)的核矩阵形式为：
其中，KXX＝φ(X)Tφ(X)，KAA＝φ(A)Tφ(A)，KXA＝φ(X)Tφ(A)是对应的核矩阵。由于核函数的存在，这里采用投影梯度法(Projected Gradient,PG)[11]对端元矩阵A进行更新：
其中，max(·,0)用于保证非负性，t是用投影梯度法更新时的迭代次数，是优化函数对端元矩阵A的偏导，这里我们采取逐个端元求导：
尺度因子ηt的计算采用回溯Armoji搜索法[12]，使得在该过程中，目标函数值都能满足一定的单调下降：
f(At+1)-f(At)≤γηtvec(▽Af(At))Tvec(At+1-At) (17)
vec(·)是将矩阵转为向量，γ是下降程度，通常可设为0.01。
(2)优化变量矩阵B
保留式(12)中与变量矩阵B相关的项，将式(12)重新写成：
对矩阵B逐个像素进行求解，表达式如下：
(3)优化丰度矩阵S
保留式(12)中与丰度矩阵S相关的项，将式(12)重新写成：
因为难以对式(20)直接逐像素求丰度，所以这里采用了ADMM的方法[13]进行求解，并引入两个变量g1和g2：
对应的增广拉格朗日表达式为：
其中，和分别是g1和g2的列向量，v1和v2是尺度化的对偶变量，和是它们对应的列向量，ρ是惩罚参数。我们逐一对式(22)中的变量进行求解，可以得到：
其中，L＝D-W是拉普拉斯矩阵，D是对角矩阵，对角元素满足
(4)优化变异系数矩阵Z
保留式(12)中与丰度矩阵Z相关的项，将式(12)重新写成：
其中，是g3的列向量，对应的增广拉格朗日表达式是：
其中，v3是尺度化的对偶变量，是它的列向量，θ是惩罚参数，上述三个变量的求解表达式如下：
根据上述内容，本发明采用的算法中采用的具体步骤归纳如下：
输入：高光谱数据X∈RL×N；
输出：端元矩阵A∈RL×P、丰度矩阵S∈RP×N、变异系数矩阵Z∈RP×N。
步骤1优化端元矩阵A：
1.1：按式(16)对端元矩阵求偏导；
1.2：按式(15)用投影梯度法对端元矩阵进行优化；
1.3：根据Armoji搜索法对尺度因子进行调整，使式(17)可以被满足。
步骤2优化变量矩阵B：
2.1：按式(19)对变量矩阵B进行逐像素求解。
步骤3优化丰度矩阵S
3.1：按式(23)对丰度矩阵S进行逐像素求解；
3.2：按式(24)对引入的变量g1进行求解；
3.3：为了满足非负和“和为1”，按式(24)对变量g2进行更新。
3.4：按式(26)更新对偶变量v1和v2。
步骤4优化变异系数矩阵Z
3.1：按式(29)对变异系数矩阵S进行逐像素求解；
3.2：按式(30)对引入的变量g3进行求解；
3.3：按式(31更新对偶变量v3。
步骤5：满足以下其中一个条件时跳出循环，输出结果，否则返回步骤1：
条件1：迭代次数达到200
条件2：每次迭代与前一次迭代之间的目标函数值相对变化小于10-2。
本发明实施方式中，一些参数取如下固定值：
λ取为0.01；
σ取为50；
ρ取为0.1；
α取为0.001；
θ取为0.1；
β取为0.001；
最大迭代次数为200。
迭代过程中，通过判断每次迭代与前一次迭代之间的目标函数值相对变化是否小于10-2来判断是否收敛。
本发明的有益效果在于：通过核函数将原始高光谱图像数据隐式地映射到高维特征空间中，从而在该空间中结合光谱变异性进行线性解混；与此同时，依据实际地物的分布特性，添加丰度和光谱变异系数的局部平滑约束，可以克服不同非线性混合场景中存在的光谱变异性问题。在基于高光谱遥感图像的高精度解混以及地面目标的检测和识别方面具有重要的应用价值。
仿真和真实高光谱遥感图像的实验表明，与相关算法相比，本发明方法中采取的算法在有效地解释非线性混合效应的同时，较好地克服了端元光谱变异性问题，对实际应用具有重要的意义。
附图说明
图1USGS光谱库中5条端元光谱。
图2端元的丰度图(上面)和变异系数图(下面)，从左到右依次为：(a)绿脱石Nontronite,(b)高岭石Kaolinite,(c)铁铝榴石Almandine,(d)山杨Aspen,(e)古铜辉石Bronzite。
图3真实高光谱遥感图像。其中，(a)AVIRIS Cuprite及子区域(b)AVIRIS Moffett及子区域。
图4Cuprite数据各算法估计的丰度图和UNSUSC-SV估计的变异系数图。其中，从左到右每列为(a)Bi-objective NMF，(b)ASSKNMF，(c)UNSU-SV和(d)UNSUSC-SV算法所获得的丰度图；(e)UNSUSC-SV算法所获得的变异系数图；从上到下每行分别为白云母，明矾石和高岭石三种物质。
图5Moffett数据各算法估计的丰度图和UNSUSC-SV估计的变异系数图。其中，从左到右每列为(a)Bi-objective NMF，(b)ASSKNMF，(c)UNSU-SV和(d)UNSUSC-SV算法所获得的丰度图；(e)UNSUSC-SV算法所获得的变异系数图；从上到下每行分别为植被，水体和土壤三种物质。
具体实施方式
下面分别用模拟数据和实际遥感图像数据为例说明本发明的具体的实施方式。
采用的考虑光谱变异性的高光谱遥感图像非线性解混算法用UNSUSC-SV表示。
1、模拟数据实验
在本节将UNSUSC-SV算法与双目标非负矩阵分解算法Bi-objective NMF[8]，基于丰度约束核非负矩阵分解的非线性解混算法ASSKNMF[9]，进行解混性能的比较，同时将本发明提出的不加丰度和光谱变异系数平滑约束的方法记为UNSU-SV，以考察添加平滑约束的作用。利用端元的光谱角距离SAD(Spectral Angle Distance)，丰度的均方根误差RMSE(S)(Root Mean Square Error)，端元矩阵的均方根误差RMSE(An)和重构误差RE(Reconstruction Error)比较解混结果：
其中，ap,sn和xn是求得的端元光谱，丰度值和像元光谱，ap′,sn′和xn′是真实的端元光谱，丰度值和像元光谱。式(35)度量了基础端元矩阵的求解精度，但因为光谱变异性的存在，每个像元的端元矩阵都不同，所以进一步用式(35)的RMSE(An)来度量求得的每个像元的端元矩阵An和真实的端元矩阵An′之间的差异。
从美国地质调查局(USGS)光谱库中选取五种不同地物的光谱作为端元生成模拟数据，该光谱库中地物光谱的波长在0.38μm-2.5μm之间，分辨率10nm，共有224个波段，图1显示了这5种端元的光谱曲线。分别针对Hapke和GBM两种非线性混合模型生成数据的方式如下：
(1)首先生成丰度，将36×36大小的图像分割为互不重叠的大小为6×6的正方形块，并随机地分配不同的端元给这些块；用7×7的低通滤波器对地物的丰度进行滤波；变异系数矩阵的生成方式和丰度矩阵类似，并在1.5～0.5的范围内随机生成，图2显示了5个端元对应的丰度图和变异系数图；
(2)生成Hapke模型数据时，先将端元光谱反射率按照式(4)转为单次散射反照率，然后与生成的变异系数相乘，以添加光谱变异性。接着，与1)中产生的丰度按LMM的形式混合得到单次散射反照率，再按式(2)转换为混合光谱反射率数据。最后，向生成的模拟数据中添加不同信噪比的高斯白噪声；
(3)生成GBM模型数据时，将端元矩阵乘以每个像元的变异系数，然后将产生的丰度带入式(6)逐像素地产生双线性混合数据，参数ζ在[0,1]中随机取值。生成模拟数据后，再添加不同信噪比的高斯白噪声。
通过改变端元数目和噪声并比较不同端元数目下运行时间来全面地评价算法的性能。每个实验都在相同条件下独立运行10次。
实验1端元数目和信噪比的影响分析：在该实验中比较了各算法在不同的端元数目和信噪比下的解混精度。图像大小是36×36，端元数目取3、4、5。在端元数目固定的情况下，用端元光谱分别生成Hapke和GBM模型的混合数据，再添加高斯白噪声，信噪比取50dB、40dB、30dB三种情况。表1～3分别显示了端元数目为3、4、5的情况下各算法的解混结果，最优结果由粗体标出。对表中的数据进行综合分析，可以得出以下结论：
(1)端元提取方面：UNSU-SV和本发明提出的方法UNSUSC-SV考虑了光谱变异性的存在，SAD比没有考虑光谱变异性的Bi-objective NMF和ASSKNMF算法要低，同时还可以逐像素求得每个像元的端元矩阵，很好地凸显了光谱变异性；UNSUSC-SV在UNSU-SV的基础上，对丰度和变异系数添加了局部平滑约束，端元提取精度有所提高。Bi-objective NMF在解混中同时权衡了原始低维空间和高维特征空间中的数据重构误差，但是效果不理想，ASSKNMF通过添加丰度的稀疏和平滑约束，结果略优于Bi-objective NMF；
(2)丰度估计方面：本发明所提出的方法UNSUSC-SV既考虑了光谱变异性，又添加了局部平滑约束，RMSE是最低的，只有其他三种方法的十分之一左右。UNSU-SV只考虑了光谱变异性，没有对丰度添加任何约束，所以结果比添加了稀疏和平滑约束的ASSKNMF方法要差；Bi-objective NMF既没有考虑光谱变异性，也没有添加约束，丰度估计结果是最不好的；
(3)重构误差方面：UNSU-SV和本发明提出的方法UNSUSC-SV同时考虑了非线性解混和光谱变异性，得到的RE基本在1e-4数量级，Bi-objective NMF和ASSKNMF算法只进行了非线性解混，没有考虑光谱变异性，得出的RE是它们的80至20倍左右，效果差很多；而相对UNSU-SV的RE来说，UNSUSC-SV还要更小一点。
端元数目为3时各算法结果比较见表1。端元数目为4时各算法结果比较见表2。端元数目为5时各算法结果比较见表3。
实验2运算时间分析：在该实验中比较了各算法在不同端元数目下的运算时间。图像大小是36×36，信噪比是40dB，端元数目是3、4、5。表4给出了端元数目分别为3、4、5情况下各算法的运行时间。从表中的数据可以看出：1)各算法在两种非线性模型下的运算时间几乎相同，相差基本都在1～2s。2)端元数目的增加会使得各算法的运算时间有所上升，Bi-objective NMF算法受影响最大，其余算法只是略微影响。3)Bi-objective NMF的运算时间最长，ASSKNMF和UNSU-SV时间最短，而本发明所提出的方法UNSUSC-SV的运算时间只比ASSKNMF和UNSU-SV多1～5s，但是解混效果是各算法中最好的。不同端元数目下算法运行时间比较(s)见表4。
2、真实遥感图像实验
本节首先采用了机载可见光/红外成像光谱仪(Airborne Visible/InfraredImaging Spectrometer,AVIRIS)获取于1997年6月19日的美国内华达州Cuprite矿区高光谱图像来分析和评价各个算法的丰度估计结果。该图像大小为512×614(如图3(a)所示)，在0.4-2.5μm的波长区间内有224个波段，光谱分辨率为10nm，空间分辨率为20m。解混前剔除噪声影响较大及水汽吸收波段(1-2,104-113,148-167,221-224)后剩余188个波段。另外，还将AVIRIS收集的Moffett地区高光谱图像用于实验(如图3(b)所示)，该地区波长范围和空间、光谱分辨率和Cuprite地区相同，剔除信噪比低和水吸收波段后剩余188个波段。通过端元光谱角距离SAD和重构误差RE来评价各算法解混性能，并展示了各算法求解的丰度图。
从第一幅Cuprite图像中选取一块大小为50×50的子区域作为测试数据，如图3(a)所示。该地区主要有白云母(Muscovite)、明矾石(Muscovite)和高岭石(Kaolinite)三种地物，表5中列出了各算法对该区域解混得到的SAD以及RE值，UNSUSC-SV的结果依然好于其他三种方法，Alunite的光谱角距离比UNSU-SV略高。图4描述了各方法得到的对应物质丰度图以及UNSUSC-SV求得的变异系数图。从第二幅Moffett图像中截取大小为50×50的子区域作为测试数据，如图3(b)所示。根据以往研究可知，该区域主要包含了植被，水体和土壤三种物质，表6中列出了各算法对该子区域的端元提取结果和重构误差。本发明提出的方法UNSUSC-SV的解混结果要优于其它比较方法，SAD和RE的平均值都是最低的，水体的提取结果比ASSKNMF略差，RE还是在1e-4左右，比Bi-objective NMF和ASSKNMF算法低一个数量级，图5展示了四种方法求解的端元对应的丰度图以及UNSUSC-SV求得的变异系数图。
从以上实际数据的解混结果可以看出，本发明提出的方法UNSUSC-SV在实际的非线性数据中，解混效果在各方法中是最好的，同时从UNSUSC-SV求得的两个地区的光谱变异数。图可以看出，物质在自身分布多的地方相对来说越容易发生光谱变异。不同算法在Cuprite数据上的光谱角距离和重构误差比较见表5。不同算法在Moffett数据上的光谱角距离和重构误差比较见表6。
综上可知，对于模拟和实际高光谱数据来说，本发明提出的算法相对于其他算法而言，能够克服在不同的非线性混合场景中存在的光谱变异性问题，改善高光谱图像的解混精度。
表1端元数目为3时各算法结果比较
表2端元数目为4时各算法结果比较
表3端元数目为5时各算法结果比较
表4不同端元数目下算法运行时间比较(s)
表5不同算法在Cuprite数据上的光谱角距离和重构误差比较
表6不同算法在Moffett数据上的光谱角距离和重构误差比较
参考文献：
[1]Tong Q,Xue Y,Zhang L.Progress in hyperspectral remote sensingscience and technology in china over the past three decades[J].IEEEJ.Sel.Topics Appl.Earth Observ.Remote Sens.,2014,7(1):70–91.
[2]Bioucas-Dias J M,Plaza A,Dobigeon N,et al.Hyperspectral unmixingoverview:Geometrical,statistical,and sparse regression-based approaches[J].IEEE J.Sel.Topics Appl.Earth Observ.Remote Sens.,2012,5(2):354–379.
[3]Somers B,Asner G P,Tits L,et al.Endmember variability in spectralmixture analysis:A review[J].Remote Sensing of Environment,2011,115(7):1603-1616.
[4]YANG Bin,WANG Bin.Review of nonlinear unmixing for hyperspectralremote sensing imagery[J].J.Infrared Millim.Waves.(杨斌,王斌.高光谱遥感图像非线性解混综述.红外与毫米波学报)，2017,36(2):173–185.
[5]Hapke B W.Bidirectional reflectance spectroscopy.I.Theory[J].J.Geophys.Res.,1981,86:3039–3054.
[6]Halimi A,Altmann Y,Dobigeon N,et al.Nonlinear unmixing ofhyperspectral images using a generalized bilinear model[J].IEEETrans.Geosci.Remote Sens.,2011,49(11):4153–4162.
[7]Ammanouil R,Ferrari A,Richard C,et al.Nonlinear unmixing ofhyperspectral data with vector-valued kernel functions[J].IEEE Trans.ImageProcess.,2017,26(1):340-354.
[8]Zhu F,Honeine P.Biobjective nonnegative matrix factorization:Linear versus kernel-based models[J].IEEE Trans.Geosci.Remote Sens.,2016,54(7):4012-4022.
[9]智通祥，杨斌，王斌.基于丰度约束核非负矩阵分解的高光谱图像非线性解混.复旦学报自然科学版(已录用).
[10]Smola A J,B.Learning with kernels[C].Cambridge:MITPress,2008.
[11]Lin C J.Projected gradient methods for nonnegative matrixfactorization[J].Neural computation,2007,19(10):2756-2779.
[12]Huck A,Guillaume M,Blanc-Talon J.Minimum dispersion constrainednonnegative matrix factorization to unmix hyperspectral data[J].IEEETrans.Geosci.Remote Sens.,2010,48(6):2590-2602.
[13]Boyd S,Parikh N,Chu E,et al.Distributed optimization andstatistical learning via the alternating direction method of multipliers[J].Foundations andin Machine Learning,2011,3(1):1-122.。
Claims (9)
1.一种考虑光谱变异性的高光谱图像非线性解混方法，其特征在于，通过核函数将原始高光谱图像数据隐式地映射到高维特征空间中，从而在该空间中结合光谱变异性进行线性解混；同时依据实际地物的分布特性，在丰度和光谱变异系数上添加局部平滑约束；其中：
针对Hapke和GBM两类非线性混合模型中存在光谱变异性时，高光谱遥感图像中的每个像元xn∈RL×1(n＝1,2,...N)表示为：
φ(xn)＝φ(A)diag(zn)sn＝φ(A)(zn⊙sn) (1)
这里，φ(·)是非线性映射，A＝[a1,a2,…,aP]∈RL×P表示端元矩阵，P是端元数，和分别是像元xn的光谱变异系数和丰度；再结合地物分布特性加入对应的平滑约束，得到目标函数：
这里引入了新变量B＝{b1，，，bn，，，bN}，约束参数λ，λ衡量所引入变量bn与zn⊙sn之间的差异；将像元光谱差异wn,m作为权重，使和ψ(Z)较好地反应中心像素丰度和变异系数与邻域像素对应取值间的相似性；α和β分别是相应的惩罚系数；τC(S)是丰度，满足“非负”与“和为一”约束的指示函数；针对目标函数的变量，采取逐一交替优化的方法；
算法初始化：端元矩阵A和丰度矩阵S分别采用VCA和FCLS进行初始化，变异系数矩阵Z每个元素均初始取1，变量矩阵B用A和S初始值的乘积初始化；
执行以下步骤循环：
(1)优化端元矩阵A
步骤1.1：保留式(2)中与端元矩阵A相关的项，将式(2)重新写成：
步骤1.2：由于核函数的存在，采用投影梯度法对端元矩阵A进行更新：
At+1＝max(A-ηt▽Af(At),0) (4)
其中，t是用投影梯度法更新时的迭代次数，尺度因子ηt的计算采用Armoji搜索法，使得目标函数值能满足一定的单调下降；
(2)优化变量矩阵B
步骤2.1：保留式(2)中与变量矩阵B相关的项，将式(2)重新写成：
步骤2.2：对矩阵B逐个像素进行求解：
其中，Qp,q＝κ(ap,aq),cp,n＝κ(ap,xn)，IP×P是一个单位矩阵；κ是核函数，表达式是：
(3)优化丰度矩阵S
步骤3.1：保留式(2)中与丰度矩阵S相关的项，将式(2)重新写成：
对于式(8)，采用ADMM的方法进行求解，引入两个变量g1和g2，对应的增广拉格朗日表达式是：
其中，和分别是g1和g2的列向量，v1和v2是尺度化的对偶变量，和是它们对应的列向量，ρ是惩罚参数；
步骤3.2：对式(9)中的变量进行逐一求解：
(3.2a)
(3.2b)g1＝ρ(S+v1)(αL+ρIN×N)-1；
其中，L＝D-W是一个拉普拉斯矩阵，D是一个对角矩阵，对角元素满足
(3.2c)
(3.2d)
(4)优化变异系数矩阵Z
步骤4.1：变异系数矩阵Z的求解过程和丰度矩阵S类似，引入变量g3，对应的增广拉格朗日表达式是：
其中，是g3的列向量，v3是尺度化的对偶变量，是它的列向量，θ是惩罚参数；
步骤4.2：对式(10)中的变量进行逐一求解：
(4.2a)
(4.2b)g3＝θ(Z+v3)(βL+θIN×N)-1；
(4.2c)
(4.2d)如达到最大迭代次数或满足收敛精度则输出矩阵A,B,S,Z，否则回到步骤(1)。
2.根据权利要求1所述的解混方法，其特征在于，步骤2.2中的λ取为0.01。
3.根据权利要求1所述的解混方法，其特征在于，步骤2.2中的σ取为50。
4.根据权利要求1所述的解混方法，其特征在于，步骤(3.2a)中的ρ取为0.1。
5.根据权利要求1所述的解混方法，其特征在于，步骤(3.2b)中的α取为0.001。
6.根据权利要求1所述的解混方法，其特征在于，步骤(4.2a)中的θ取为0.1。
7.根据权利要求1所述的解混方法，其特征在于，步骤(3.2b)中的β取为0.001。
8.根据权利要求1所述的解混方法，其特征在于，步骤(4.2d)中的最大迭代次数为200。
9.根据权利要求1所述的解混方法，其特征在于，步骤(4.2d)中，通过判断每次迭代与前一次迭代之间的目标函数值相对变化是否小于10-2来判断是否收敛。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810579080.9A CN109035154B (zh) | 2018-06-07 | 2018-06-07 | 一种考虑光谱变异性的高光谱图像非线性解混方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810579080.9A CN109035154B (zh) | 2018-06-07 | 2018-06-07 | 一种考虑光谱变异性的高光谱图像非线性解混方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN109035154A true CN109035154A (zh) | 2018-12-18 |
| CN109035154B CN109035154B (zh) | 2021-08-20 |
Family
ID=64612063
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201810579080.9A Active CN109035154B (zh) | 2018-06-07 | 2018-06-07 | 一种考虑光谱变异性的高光谱图像非线性解混方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN109035154B (zh) |
Cited By (9)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109800690A (zh) * | 2019-01-04 | 2019-05-24 | 北京科技大学 | 一种非线性高光谱影像混合像元分解方法及装置 |
| CN109829946A (zh) * | 2019-01-18 | 2019-05-31 | 北京理工大学 | 基于快速混合迭代的map-tv高光谱亚像元定位方法 |
| CN109871774A (zh) * | 2019-01-22 | 2019-06-11 | 中国科学院南京土壤研究所 | 一种基于局域相近像元的混合像元分解方法 |
| CN110070004A (zh) * | 2019-04-02 | 2019-07-30 | 杭州电子科技大学 | 一种应用于深度学习的近地高光谱数据扩展方法 |
| CN111008975A (zh) * | 2019-12-02 | 2020-04-14 | 北京航空航天大学 | 一种空间人造目标线性模型的混合像元解混方法及系统 |
| CN112396029A (zh) * | 2020-12-03 | 2021-02-23 | 宁波大学 | 一种聚类分割与耦合端元提取协同的高光谱滨海湿地亚像元变化检测方法 |
| CN112819769A (zh) * | 2021-01-26 | 2021-05-18 | 复旦大学 | 基于核函数和联合字典的非线性高光谱图像异常探测算法 |
| CN113624691A (zh) * | 2020-05-07 | 2021-11-09 | 南京航空航天大学 | 一种基于空间-光谱相关性的光谱图像超分辨率映射方法 |
| CN114596482A (zh) * | 2022-02-10 | 2022-06-07 | 复旦大学 | 一种基于扩展多线性混合模型的高光谱图像非线性解混方法 |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106056524A (zh) * | 2016-05-25 | 2016-10-26 | 天津商业大学 | 基于微分搜索的高光谱图像非线性解混方法 |
| CN106778530A (zh) * | 2016-11-28 | 2017-05-31 | 复旦大学 | 一种基于双线性混合模型的高光谱图像非线性解混方法 |
| US20170191929A1 (en) * | 2015-12-31 | 2017-07-06 | Abb, Inc. | Spectral modeling for complex absorption spectrum interpretation |
-
2018
- 2018-06-07 CN CN201810579080.9A patent/CN109035154B/zh active Active
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20170191929A1 (en) * | 2015-12-31 | 2017-07-06 | Abb, Inc. | Spectral modeling for complex absorption spectrum interpretation |
| CN106056524A (zh) * | 2016-05-25 | 2016-10-26 | 天津商业大学 | 基于微分搜索的高光谱图像非线性解混方法 |
| CN106778530A (zh) * | 2016-11-28 | 2017-05-31 | 复旦大学 | 一种基于双线性混合模型的高光谱图像非线性解混方法 |
Non-Patent Citations (4)
| Title |
|---|
| DRUMETZ L ET AL.: ""Blind hyperspectral unmixing using an Extended Linear Mixing Model to address spectral variability"", 《IEEE TRANS.IMAGE PROCESS.》 * |
| HALIMI A ET AL.: ""Nonlinear unmixing of hyperspectral images using a generalized bilinear model"", 《IEEE TRANS． GEOSCI． ＲEMOTE SENS．》 * |
| HEYLEN Ｒ ET AL.: ""A review of nonlinear hyperspectral unmixing methodsHeylen Ｒ"", 《IEEE J． SEL． TOPICS APPL》 * |
| ZHU FEI ET AL.: ""Biobjective Nonnegative Matrix Factorization: Linear Versus Kernel-Based Models "", 《IEEE TRANS． GEOSCI． ＲEMOTE SENS．》 * |
Cited By (15)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109800690B (zh) * | 2019-01-04 | 2020-10-13 | 北京科技大学 | 一种非线性高光谱影像混合像元分解方法及装置 |
| CN109800690A (zh) * | 2019-01-04 | 2019-05-24 | 北京科技大学 | 一种非线性高光谱影像混合像元分解方法及装置 |
| CN109829946A (zh) * | 2019-01-18 | 2019-05-31 | 北京理工大学 | 基于快速混合迭代的map-tv高光谱亚像元定位方法 |
| CN109871774A (zh) * | 2019-01-22 | 2019-06-11 | 中国科学院南京土壤研究所 | 一种基于局域相近像元的混合像元分解方法 |
| CN109871774B (zh) * | 2019-01-22 | 2020-12-18 | 中国科学院南京土壤研究所 | 一种基于局域相近像元的混合像元分解方法 |
| CN110070004A (zh) * | 2019-04-02 | 2019-07-30 | 杭州电子科技大学 | 一种应用于深度学习的近地高光谱数据扩展方法 |
| CN110070004B (zh) * | 2019-04-02 | 2021-07-20 | 杭州电子科技大学 | 一种应用于深度学习的近地高光谱数据扩展方法 |
| CN111008975A (zh) * | 2019-12-02 | 2020-04-14 | 北京航空航天大学 | 一种空间人造目标线性模型的混合像元解混方法及系统 |
| CN111008975B (zh) * | 2019-12-02 | 2022-09-09 | 北京航空航天大学 | 一种空间人造目标线性模型的混合像元解混方法及系统 |
| CN113624691A (zh) * | 2020-05-07 | 2021-11-09 | 南京航空航天大学 | 一种基于空间-光谱相关性的光谱图像超分辨率映射方法 |
| CN112396029A (zh) * | 2020-12-03 | 2021-02-23 | 宁波大学 | 一种聚类分割与耦合端元提取协同的高光谱滨海湿地亚像元变化检测方法 |
| CN112396029B (zh) * | 2020-12-03 | 2022-02-18 | 宁波大学 | 一种聚类分割与耦合端元提取协同的高光谱滨海湿地亚像元变化检测方法 |
| CN112819769A (zh) * | 2021-01-26 | 2021-05-18 | 复旦大学 | 基于核函数和联合字典的非线性高光谱图像异常探测算法 |
| CN114596482A (zh) * | 2022-02-10 | 2022-06-07 | 复旦大学 | 一种基于扩展多线性混合模型的高光谱图像非线性解混方法 |
| CN114596482B (zh) * | 2022-02-10 | 2023-04-07 | 复旦大学 | 一种基于扩展多线性混合模型的高光谱图像非线性解混方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN109035154B (zh) | 2021-08-20 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN109035154A (zh) | 一种考虑光谱变异性的高光谱图像非线性解混方法 | |
| Landgrebe | Information extraction principles and methods for multispectral and hyperspectral image data | |
| CN104952050B (zh) | 基于区域分割的高光谱图像自适应解混方法 | |
| Longepe et al. | Assessment of ALOS PALSAR 50 m orthorectified FBD data for regional land cover classification by support vector machines | |
| Chen et al. | Spectral–spatial preprocessing using multihypothesis prediction for noise-robust hyperspectral image classification | |
| CN104978573B (zh) | 一种应用于高光谱图像处理的非负矩阵分解方法 | |
| CN103413292B (zh) | 基于约束最小二乘的高光谱图像非线性丰度估计方法 | |
| CN103824084A (zh) | 基于sdit和svm的极化sar图像分类方法 | |
| Yokoya et al. | Coupled non-negative matrix factorization (CNMF) for hyperspectral and multispectral data fusion: Application to pasture classification | |
| CN106228130B (zh) | 基于模糊自编码网络的遥感影像云检测方法 | |
| CN107274387A (zh) | 基于进化多目标优化的高光谱遥感影像的端元提取方法 | |
| Sun et al. | Weighted Collaborative Sparse and L 1/2 Low-Rank Regularizations With Superpixel Segmentation for Hyperspectral Unmixing | |
| CN103927551A (zh) | 基于超像素相关矩阵的极化sar图像半监督分类方法 | |
| Fertig et al. | Assimilating non-local observations with a local ensemble Kalman filter | |
| CN106778530A (zh) | 一种基于双线性混合模型的高光谱图像非线性解混方法 | |
| Moudden et al. | Blind component separation in wavelet space: Application to CMB analysis | |
| Huang et al. | Extending RAPID model to simulate forest microwave backscattering | |
| Meganem et al. | Physical modelling and non-linear unmixing method for urban hyperspectral images | |
| Yang et al. | Combined nonlocal spatial information and spatial group sparsity in NMF for hyperspectral unmixing | |
| Acito et al. | CWV-Net: A deep neural network for atmospheric column water vapor retrieval from hyperspectral VNIR data | |
| Sun et al. | Blind unmixing of hyperspectral images based on L 1 norm and tucker tensor decomposition | |
| Schott et al. | A subpixel target detection technique based on the invariance approach | |
| Li et al. | A multispectral remote sensing data spectral unmixing algorithm based on variational Bayesian ICA | |
| CN114596482B (zh) | 一种基于扩展多线性混合模型的高光谱图像非线性解混方法 | |
| Zhang et al. | Sparse hyperspectral unmixing using spectral library adaptive adjustment |
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 |