CN113743325A - 监督与非监督的高光谱混合像元分解方法 - Google Patents
监督与非监督的高光谱混合像元分解方法 Download PDFInfo
- Publication number
- CN113743325A CN113743325A CN202111045548.4A CN202111045548A CN113743325A CN 113743325 A CN113743325 A CN 113743325A CN 202111045548 A CN202111045548 A CN 202111045548A CN 113743325 A CN113743325 A CN 113743325A
- Authority
- CN
- China
- Prior art keywords
- algorithm
- end member
- supervised
- volume
- abundance
- 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 131
- 238000000354 decomposition reaction Methods 0.000 title claims abstract description 59
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 113
- 239000013598 vector Substances 0.000 claims abstract description 37
- 238000000605 extraction Methods 0.000 claims abstract description 25
- 238000004458 analytical method Methods 0.000 claims abstract description 24
- 230000012010 growth Effects 0.000 claims abstract description 8
- 239000011159 matrix material Substances 0.000 claims description 30
- 238000001228 spectrum Methods 0.000 claims description 20
- 238000012880 independent component analysis Methods 0.000 claims description 17
- 238000004364 calculation method Methods 0.000 claims description 7
- 230000008569 process Effects 0.000 claims description 7
- 230000009467 reduction Effects 0.000 claims description 6
- 239000000284 extract Substances 0.000 claims description 4
- 238000000926 separation method Methods 0.000 claims description 4
- 230000009466 transformation Effects 0.000 claims description 4
- 230000003416 augmentation Effects 0.000 claims description 3
- 230000006872 improvement Effects 0.000 claims description 3
- 238000005457 optimization Methods 0.000 claims description 3
- 230000003190 augmentative effect Effects 0.000 claims description 2
- 230000021332 multicellular organism growth Effects 0.000 claims 1
- 238000004088 simulation Methods 0.000 abstract description 6
- 238000012795 verification Methods 0.000 abstract description 2
- 230000000694 effects Effects 0.000 description 13
- 238000002474 experimental method Methods 0.000 description 13
- 230000003595 spectral effect Effects 0.000 description 13
- 238000010586 diagram Methods 0.000 description 10
- 238000002156 mixing Methods 0.000 description 8
- 229910052500 inorganic mineral Inorganic materials 0.000 description 7
- 239000011707 mineral Substances 0.000 description 7
- 238000012360 testing method Methods 0.000 description 7
- 238000012545 processing Methods 0.000 description 6
- 238000005516 engineering process 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
- 238000011160 research Methods 0.000 description 4
- 244000198134 Agave sisalana Species 0.000 description 3
- 229910052934 alunite Inorganic materials 0.000 description 3
- 239000010424 alunite Substances 0.000 description 3
- BERDEBHAJNAUOM-UHFFFAOYSA-N copper(I) oxide Inorganic materials [Cu]O[Cu] BERDEBHAJNAUOM-UHFFFAOYSA-N 0.000 description 3
- LBJNMUFDOHXDFG-UHFFFAOYSA-N copper;hydrate Chemical compound O.[Cu].[Cu] LBJNMUFDOHXDFG-UHFFFAOYSA-N 0.000 description 3
- 238000001514 detection method Methods 0.000 description 3
- GUJOJGAPFQRJSV-UHFFFAOYSA-N dialuminum;dioxosilane;oxygen(2-);hydrate Chemical compound O.[O-2].[O-2].[O-2].[Al+3].[Al+3].O=[Si]=O.O=[Si]=O.O=[Si]=O.O=[Si]=O GUJOJGAPFQRJSV-UHFFFAOYSA-N 0.000 description 3
- 229910052935 jarosite Inorganic materials 0.000 description 3
- 229910052901 montmorillonite Inorganic materials 0.000 description 3
- 229910052903 pyrophyllite Inorganic materials 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
- 230000008859 change Effects 0.000 description 2
- 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 2
- 238000003384 imaging method Methods 0.000 description 2
- 229910052627 muscovite Inorganic materials 0.000 description 2
- 229910000273 nontronite Inorganic materials 0.000 description 2
- 239000000126 substance Substances 0.000 description 2
- 238000010521 absorption reaction Methods 0.000 description 1
- 239000000654 additive Substances 0.000 description 1
- 230000000996 additive effect Effects 0.000 description 1
- 238000000701 chemical imaging Methods 0.000 description 1
- 238000004891 communication Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- FDKCTEWMJWRPDS-UHFFFAOYSA-N dialuminum;trimagnesium;trisilicate Chemical compound [Mg+2].[Mg+2].[Mg+2].[Al+3].[Al+3].[O-][Si]([O-])([O-])[O-].[O-][Si]([O-])([O-])[O-].[O-][Si]([O-])([O-])[O-] FDKCTEWMJWRPDS-UHFFFAOYSA-N 0.000 description 1
- 229910052852 dumortierite Inorganic materials 0.000 description 1
- 238000000295 emission spectrum Methods 0.000 description 1
- 238000011156 evaluation Methods 0.000 description 1
- 230000036039 immunity Effects 0.000 description 1
- 230000007246 mechanism Effects 0.000 description 1
- 238000012544 monitoring process Methods 0.000 description 1
- 229910052832 pyrope Inorganic materials 0.000 description 1
- 238000004445 quantitative analysis Methods 0.000 description 1
- 238000007619 statistical method 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
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F18/00—Pattern recognition
- G06F18/20—Analysing
- G06F18/21—Design or setup of recognition systems or techniques; Extraction of features in feature space; Blind source separation
- G06F18/213—Feature extraction, e.g. by transforming the feature space; Summarisation; Mappings, e.g. subspace methods
- G06F18/2133—Feature extraction, e.g. by transforming the feature space; Summarisation; Mappings, e.g. subspace methods based on naturality criteria, e.g. with non-negative factorisation or negative correlation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F18/00—Pattern recognition
- G06F18/20—Analysing
- G06F18/21—Design or setup of recognition systems or techniques; Extraction of features in feature space; Blind source separation
- G06F18/213—Feature extraction, e.g. by transforming the feature space; Summarisation; Mappings, e.g. subspace methods
- G06F18/2134—Feature extraction, e.g. by transforming the feature space; Summarisation; Mappings, e.g. subspace methods based on separation criteria, e.g. independent component analysis
-
- 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)
- Computer Vision & Pattern Recognition (AREA)
- Data Mining & Analysis (AREA)
- Theoretical Computer Science (AREA)
- Bioinformatics & Computational Biology (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Artificial Intelligence (AREA)
- Evolutionary Biology (AREA)
- Evolutionary Computation (AREA)
- Physics & Mathematics (AREA)
- General Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Image Analysis (AREA)
Abstract
本发明公开了监督与非监督的高光谱混合像元分解方法，包括监督类混合像元分解算法和非监督类混合像元分解算法，所述监督类混合像元分解算法包括端元提取算法和丰度估计算法，所述端元提取算法包括纯净像元指数法、内部最大体积法、顶点成分分析法、单形体增长分析算法、顺序最大角凸锥和分裂增广拉格朗日法，所述丰度估计算法包括最小二乘法、基于端元投影向量的算法和基于单形体体积的算法。本发明提供了监督与非监督的高光谱混合像元分解方法，本发明分别介绍了监督类和非监督类的解混方法，并通过模拟数据和真实数据对典型的方法进行了实验验证，分析了各方法的优缺点和适用条件，并对实验结果进行了评价。
Description
技术领域
本发明涉及混合像元领域，尤其涉及监督与非监督的高光谱混合像元分解方法。
背景技术
遥感器所获取的地面反射或发射光谱信号是以像元为单位记录的，一个像元内仅包含一种类型，这种像元称为纯像元。然而，多数情况下一个像元内往往包含多种地表类型，这种像元就是混合像元。混合像元记录的是多种地表类型的综合光谱信息，混合像元是指在一个像元内存在有不同类型的地物，主要出现在地类的边界处。混合像元的存在是影响识别分类精度的主要因素之一，特别是对线状地类和细小地物的分类识别影响较为突出，在土地利用遥感动态监测工作中，经常遇到混合像元的难题，解决这一问题的关键在于通过一定方法找出组成混合像元的各种典型地物的比例。
在使用遥感技术对土地进行监测时，为了对土地进行计算，需要使用高光谱混合像元分解方法，而现有的高光谱混合像元分解方法，而现有的高光谱混合像元分解方法种类繁多，人们不清楚各个算法的优缺点和适用调节，无法快速的对分解方法进行选择，为此，我们提出监督与非监督的高光谱混合像元分解方法。
发明内容
基于背景技术存在的技术问题，本发明提出了监督与非监督的高光谱混合像元分解方法，以解决技术中提出的问题。
本发明提供如下技术方案：
监督与非监督的高光谱混合像元分解方法，包括监督类混合像元分解算法和非监督类混合像元分解算法，所述监督类混合像元分解算法包括端元提取算法和丰度估计算法，所述端元提取算法包括纯净像元指数法、内部最大体积法、顶点成分分析法、单形体增长分析算法、顺序最大角凸锥和分裂增广拉格朗日法，所述丰度估计算法包括最小二乘法、基于端元投影向量的算法和基于单形体体积的算法；所述非监督类混合像元分解算法包括独立成分分析法和非负矩阵分解法。
优选的，所述纯净像元指数法为单行体的形状特点决定了其在特征空间中任意直线上的投影为线段，且线段的端点必为单形体顶点在这个方向上的投影，基于这一性质，可在特征空间中随机生成若干向量，将所有像元点投影到各个向量上，直线上所有投影点中最靠外的两个便是端元点的投影，随向量方向的不断变化，记录下图像中每个像元作为极值点的次数，最后认为出现频率最高的点就是端元。
优选的，所述内部最大体积法通过寻找具有最大体积的单形体自动获取图像中的所有端元，算法首先进行MNF变换，使数据降至p-1维，图像Y中由端元矩阵M的p个端元构成的几何体的体积可表示为
其中mi为降维后第i个端元对应的p-1维列向量，根据体积的定义式(2.1)和式(2.2)，在理想状态下，由p个端元构成的单形体体积必然是所有由等数量像元构成的单形体中最大的。首先随机选择n个像元作为端元并计算体积，然后将其中一个端元用其他的像元替换并再次计算体积。如果替换后体积增大，则接受此替换，否则放弃此替换，再用新的像元替换，直到对每个端元都用所有像元替换过。在此过程中，端元单形体的体积是不断增大的，最后得到的必然是最大体积的单形体，而构成此单行体的像元就是端元。
优选的，所述顶点成分分析法通过反复寻找正交向量并计算图像矩阵在正交向量上的投影距离来逐一提取端元，单形体的若干个顶点可以张成一个子空间，而单形体在某个与这个子空间正交的向量上的投影距离的最大值点一定是单形体的顶点，可以先设法找到一个初始端元，然后每次循环都先找一个和已经找到的端元同时正交的单位向量，再将所有向量点投影到这一单位向量上，并将投影结果的最大值记为新端元，加入端元集合并开始下一次循环，直到找到设定的端元个数为止。
优选的，所述单形体增长分析算法是在N-FINDR算法基础上做了改进：与N-FINDR同时搜索所有的端元不同，单形体增长分析算法采用逐步增加单形体定点的方法搜索各个端元，在搜索过程中单形体的体积是膨胀的，其维度也逐渐增长，最后得到端元，其计算复杂度较N-FINDR要小；所述顺序最大角凸锥是对顶点成分分析法的改进，将顶点成分分析法的正交投影放宽为斜交投影，以此为代价，可以换取满足“非负”约束的端元丰度，此法可同时实现端元提取和丰度反演。
优选的，所述分裂增广拉格朗日法在使用前之前需要根据端元数目p对图像进行降维，使得L＝p，得到降维后的图像矩阵Ypxn。当L＝p时，端元矩阵M为方阵且各个端元线性无关，从而存在Q＝M-1，寻找最小体积单形体可以描述为如下的优化问题
那么提取端元M的问题就转化成了求解Q的问题。
优选的，所述最小二乘法是目前为止应用最为广泛的丰度反演方法，根据满足丰度约束条件的程度可以分为四种不同的最小二乘法：无约束最小二乘法，“和为1”约束最小二乘法，“非负”约束最小二乘法和全约束最小二乘法；
根据线性混合模型公式，不考虑误差时，可得线性混合模型公式为
当不考虑对丰度的约束时，可以求得第1.2.1节中式(1.4)的无约束解为
sUCLS＝(MTM)-1MTy (2.5)
考虑到ASC，可得出式的“和为1”约束解为
其中Ip为p阶单位矩阵；1＝[1,1,...1]T为全1(p维)列向量；
考虑到描述的ANC，可得出非负约束解SNCLS，但此解不能由简单的算子表示出来，而是利用迭代方法来获得最优解，非负约束最小二乘法的两个迭代关系式如下：
同时满足ANC和ASC的解成称为全约束最小二乘解，可将式(2.6)代入式(2.7)求解得到；
所述基于端元投影向量的算法
根据线性光谱混合模型和几何学描述，混合像元位于单形体的L维特征空间的内部，端元位于单形体的顶点，某个端元是距离其他端元构成的L-1维空间超平面最远的点，可以根据像元到超平面的距离与端元到超平面的距离之间的比值来计算像元中端元占的比例，即丰度。
所述基于单形体体积的算法
根据线性光谱混合模型和几何学描述，某个端元用其他像元替换后所得单形体的体积比原单形体体积小，可以根据用某像元替换后的单形体与原单形体的体积比计算出该像元中端元的丰度。
优选的，所述独立成分分析法一般是假设形成遥感图像的几个端元的丰度是互相独立的，通过最小化它们之间的相关信息来达到解混的目的，假定每一个端元的丰度都是一个随机信号源，这样就可以把高光谱遥感图像看成是多源混合信号，其中的光谱解混问题就转化为一个盲信号分离的问题，独立成分分析法是解决混合像元分解的一种有效方法，可以同时解决“端元提取”和“丰度估计”两个问题。
优选的，所述非负矩阵分解法其每步的迭代均在上一步的结果上使用乘法得到，由于乘法整数前后的符号是不变的，只要初始值非负就能保证结果非负；
给定一非负矩阵Y∈Rn×m和一个正整数r＜min(m,n)，NMF的目标是找到两个非负矩阵W∈Rn×r和H∈Rr×m，使其满足
Y≈WH (2.8)
目标函数有两种：第一种是Y和WH之间的欧式距离
它的下界是0，当Y＝WH时达到此下界；
第二种是Y和WH之间的“散度”，也即概率意义上的距离
同样，当Y＝WH时达到此下界；
就欧式距离情况下推导得到的乘法迭代公式如下：
W←W.*(YHT)./(WHHT) (2.11)
H←H.*(WTY)./(WTWH) (2.12)。
本发明提供了监督与非监督的高光谱混合像元分解方法，本发明分别介绍了监督类和非监督类的解混方法，并通过模拟数据和真实数据对典型的方法进行了实验验证，分析了各方法的优缺点和适用条件，并对实验结果进行了评价。
附图说明
图1为本发明四种算法在图像中有无纯净端元情况下的提取端元效果图；
图2为本发明四种算法提取混合像元中第1个端元的光谱曲线比较图；
图3为本发明四种算法在图像混合程度不同的情况下的端元提取效果图；
图4为本发明不同信噪比的图像中算法性能比较图；
图5为本发明不同像元个数下算法性能比较图；
图6为本发明Cuprite地区图像数据分布图；
图7为本发明VCA+FCLS分解丰度图；
图8为本发明MVC-NMF分解丰度图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
请参阅图1-8，本发明提供一种技术方案：
监督与非监督的高光谱混合像元分解方法，包括监督类混合像元分解算法和非监督类混合像元分解算法，所述监督类混合像元分解算法包括端元提取算法和丰度估计算法，所述端元提取算法包括纯净像元指数法、内部最大体积法、顶点成分分析法、单形体增长分析算法、顺序最大角凸锥和分裂增广拉格朗日法，所述丰度估计算法包括最小二乘法、基于端元投影向量的算法和基于单形体体积的算法；所述非监督类混合像元分解算法包括独立成分分析法和非负矩阵分解法。
进一步地，所述纯净像元指数法(PPI)为单行体的形状特点决定了其在特征空间中任意直线上的投影为线段，且线段的端点必为单形体顶点在这个方向上的投影，基于这一性质，可在特征空间中随机生成若干向量，将所有像元点投影到各个向量上，直线上所有投影点中最靠外的两个便是端元点的投影，随向量方向的不断变化，记录下图像中每个像元作为极值点的次数，最后认为出现频率最高的点就是端元。
进一步地，所述内部最大体积法(N-FINDR)通过寻找具有最大体积的单形体自动获取图像中的所有端元，算法首先进行MNF变换，使数据降至p-1维，图像Y中由端元矩阵M的p个端元构成的几何体的体积可表示为
其中mi为降维后第i个端元对应的p-1维列向量，根据体积的定义式(2.1)和式(2.2)，在理想状态下，由p个端元构成的单形体体积必然是所有由等数量像元构成的单形体中最大的。首先随机选择n个像元作为端元并计算体积，然后将其中一个端元用其他的像元替换并再次计算体积。如果替换后体积增大，则接受此替换，否则放弃此替换，再用新的像元替换，直到对每个端元都用所有像元替换过。在此过程中，端元单形体的体积是不断增大的，最后得到的必然是最大体积的单形体，而构成此单行体的像元就是端元。
进一步地，所述顶点成分分析法(VCA)通过反复寻找正交向量并计算图像矩阵在正交向量上的投影距离来逐一提取端元，单形体的若干个顶点可以张成一个子空间，而单形体在某个与这个子空间正交的向量上的投影距离的最大值点一定是单形体的顶点，可以先设法找到一个初始端元，然后每次循环都先找一个和已经找到的端元同时正交的单位向量，再将所有向量点投影到这一单位向量上，并将投影结果的最大值记为新端元，加入端元集合并开始下一次循环，直到找到设定的端元个数为止。
进一步地，所述单形体增长分析算法(SGA)是在N-FINDR算法基础上做了改进：与N-FINDR同时搜索所有的端元不同，单形体增长分析算法采用逐步增加单形体定点的方法搜索各个端元，在搜索过程中单形体的体积是膨胀的，其维度也逐渐增长，最后得到端元，其计算复杂度较N-FINDR要小；所述顺序最大角凸锥(SMACC)是对顶点成分分析法的改进，将顶点成分分析法的正交投影放宽为斜交投影，以此为代价，可以换取满足“非负”约束的端元丰度，此法可同时实现端元提取和丰度反演。
进一步地，所述分裂增广拉格朗日法的原理来源于最小体积变换，它的最小体积是通过从内向外扩张得到的，分裂增广拉格朗日法(SISIAL)在使用前之前需要根据端元数目p对图像进行降维，使得L＝p，得到降维后的图像矩阵Ypxn。当L＝p时，端元矩阵M为方阵且各个端元线性无关，从而存在Q＝M-1，寻找最小体积单形体可以描述为如下的优化问题
那么提取端元M的问题就转化成了求解Q的问题。
进一步地，所述最小二乘法是目前为止应用最为广泛的丰度反演方法，根据满足丰度约束条件的程度可以分为四种不同的最小二乘法：无约束最小二乘法，“和为1”约束最小二乘法，“非负”约束最小二乘法和全约束最小二乘法；
根据线性混合模型公式，不考虑误差时，可得线性混合模型公式为
当不考虑对丰度的约束时，可以求得第1.2.1节中式(1.4)的无约束解为
sUCLS＝(MTM)-1MTy (2.5)
考虑到ASC，可得出式的“和为1”约束解为
其中Ip为p阶单位矩阵；1＝[1,1,...1]T为全1(p维)列向量；
考虑到描述的ANC，可得出非负约束解SNCLS，但此解不能由简单的算子表示出来，而是利用迭代方法来获得最优解，非负约束最小二乘法的两个迭代关系式如下：
同时满足ANC和ASC的解成称为全约束最小二乘解，可将式(2.6)代入式(2.7)求解得到；
所述基于端元投影向量的算法
根据线性光谱混合模型和几何学描述，混合像元位于单形体的L维特征空间的内部，端元位于单形体的顶点，某个端元是距离其他端元构成的L-1维空间超平面最远的点，可以根据像元到超平面的距离与端元到超平面的距离之间的比值来计算像元中端元占的比例，即丰度。
所述基于单形体体积的算法
根据线性光谱混合模型和几何学描述，某个端元用其他像元替换后所得单形体的体积比原单形体体积小，可以根据用某像元替换后的单形体与原单形体的体积比计算出该像元中端元的丰度。
进一步地，所述独立成分分析法(ICA)一般是假设形成遥感图像的几个端元的丰度是互相独立的，通过最小化它们之间的相关信息来达到解混的目的，假定每一个端元的丰度都是一个随机信号源，这样就可以把高光谱遥感图像看成是多源混合信号，其中的光谱解混问题就转化为一个盲信号分离的问题，独立成分分析法是解决混合像元分解的一种有效方法，可以同时解决“端元提取”和“丰度估计”两个问题，基于统计学的盲信号分离BSS是信号处理领域的一种研究方法。它是指在源信号和混合方式未知的情况下，根据源信号的统计特性，仅由观测信号恢复出源信号的过程。BSS目前已在数字通信、生物医学信号处理、语音信号处理及图像信号处理等许多领域具有重要的理论研究价值和广泛的实际应用前景。近年来，BSS也被运用到高光谱混合像元分解中来。该方法克服了基于几何学方法需要存在纯像元的缺点，充分利用数据的统计特性，在端元信息(如图像中端元的个数，端元光谱等)完全未知的情况下，将混合像元分解为端元光谱及其在像元中所占的比例，又叫非监督的光谱解混。随着盲信号处理技术的兴起和发展，混合像元盲分解技术越来越受到遥感学者的关注，成为当前的一个研究热点。这类算法包括独立成分分析ICA和非负矩阵分解NMF以及一系列改进的算法；
ICA算法对于光谱解混的缺点也是非常明显的，主要包括以下几点：
(1)线性混合光谱模型中的丰度“和为1”约束不符合ICA源信号间统计独立的原则，这是影响ICA算法对光谱解混问题适用性的主要原因；
(2)端元成分有着明确的物理意义，且在遥感定量分析中对端元成分有着严格的要求，ICA中幅度(即方差)的不确定性影响端元成分以及丰度反演的精度；
(3)在实际应用中，用户希望每次运行算法都有相同的结果，并且能按照对端元信号的关心程度最先搜索出来，ICA的次序不确定性无法保证这一点。
进一步地，所述非负矩阵分解法(NMF)其每步的迭代均在上一步的结果上使用乘法得到，由于乘法整数前后的符号是不变的，只要初始值非负就能保证结果非负；
给定一非负矩阵Y∈Rn×m和一个正整数r＜min(m,n)，NMF的目标是找到两个非负矩阵W∈Rn×r和H∈Rr×m，使其满足
Y≈WH (2.8)
目标函数有两种：第一种是Y和WH之间的欧式距离
它的下界是0，当Y＝WH时达到此下界；
第二种是Y和WH之间的“散度”，也即概率意义上的距离
同样，当Y＝WH时达到此下界；
就欧式距离情况下推导得到的乘法迭代公式如下：
W←W.*(YHT)./(WHHT) (2.11)
H←H.*(WTY)./(WTWH) (2.12)；
该方法具有非负以及自动调整步长进行迭代的优点，消除了参数选择对该算法带来的影响。但由于其目标函数具有明显的非凸性，因而存在大量局部极小值；而且其解空间较大，会有很多对不唯一解，这是NMF存在的两个最大的问题；
向NMF的目标函数中加入了表征丰度稀疏性的约束项，从而使目标函数变为：
该约束项以矩阵H中所有元素的和来表征其稀疏度，并在迭代中使该数值达到最小，以使结果尽可能接近真实值。
向原始的NMF算法中加入了W的平滑性和H稀疏性约束，使得目标函数变为
试验与结果分析
实验评价指标
使用光谱误差和光谱角距离来衡量端元提取效果，计算端元光谱与观测值的近似程度，对于第k个端元，SE被定义为
SAD被定义为
其中Mk与
分别表示该端元的光谱解混结果与参考值。
通过均方根误差(Root Mean Square Error，RMSE)来衡量丰度估计的精度，RMS被定义为
其中Skj和
分别表示该端元在某一像元的丰度的光谱解混结果与参考值用到的SAD和RMSE均为所有端元的平均值。
(2)模拟数据实验
模拟数据的产生方法为：从USGS的光谱库splib06中选取498条光谱曲线(波段数为224)，并剔除其中相似性较高的曲线(光谱角距离小于10°的认为相似性较高)，在其中选取5种端元作为端元矩阵M；非零丰度系数S随机生成，且服从Dirichlet分布，满足“和为1”和“非负”约束，并且根据需要设定不同的稀疏性；根据高光谱成像机理，采用零均值高斯加性噪声。然后根据线性混合模型的原理，将M和S相乘再加上高斯白噪声就得到了实验所需的仿真数据。
首先选取四种常用的端元提取算法：基于几何学的N-FINDR、VCA、SISAL和将统计学方法与几何学方法结合起来的MVC-NMF，并对这四种算法的性能进行实验比较。
实验一、测试算法对有无纯净像元情况的适应性
本实验使用两个数据集，一个是图像中存在纯净像元，标记为PP(Pure Pixels)，另一个中不存在纯净像元，标记为NPP(Non Pure Pixels)，像元个数均为5000个，信噪比为30dB。用四种算法分别对两组数据进行了混合像元分解，端元提取效果如图1和图2所示。其中图1显示了四种算法提取出的端元和真实端元在二维散点图中的显示。图2是四种算法在两个模拟数据集中提取的第1个端元的光谱曲线。
从图1中可以看出，在存在纯净像元的情况下，四种像元解混算法均能提取出正确的端元；在图1中，四种算法提取端元的效果相对降低，出现了一定的误差，这是由于图像中不存在纯净像元所造成。但是MVC-NMF和SISAL比基于纯净端元算法(N-FINDR和VCA)提取的端元更接近真实端元，这是因为这两种算法不需要假设图像中存在纯净像元，算法适用性更广泛。从图2中提取的第1个端元与真实端元的光谱曲线相似性比较也可以得出相同的结论。
表2.1显示了四种算法在两个数据集中运行时间和端元提取精度的比较，可以看出N-FINDR和VCA运行时间很短，其次是SISAL，时间最长的是MVC-NMF，这是由于VCA和N-FINDR算法实现中都首先需要对数据进行降维，降至p-1维，这大大减少了运算量；而MVC-NMF是在NMF上加入最小体积约束的方法，NMF算法本身收敛速度较慢，计算时间长。
表2.1四种算法在两组数据中端元提取精度和运行时间比较
实验二、测试图像混合程度对算法性能的影响
当像元混合程度不同时各算法提取端元的效果。分别产生两组数据：第一组中混合程度不高，但是去掉了位于几何体顶点上的像元；第二组数据光谱混合程度非常高，像元都位于几何体内部，在各顶点或面上几乎没有像元。四种算法的混合像元端元提取效果如图3所示。
从图3两幅图像可以看出，在不存在纯净像元情况下，基于最小体积的两种算法提取端元的效果优于基于纯净端元的算法，这与实验一结论相符；但是，当混合程度特别高时，如图3(b)，四种算法都不准确，这是由于在混合程度很高的高光谱图像中，数据点(所有像元)聚集于整个凸面几何体的中间，图像的几何特性表现不明显。这种情况下，基于几何学的高光谱混合像元分解方法失效，需要进一步研究其他算法。相对而言，MVC-NMF与真实端元最为接近，这是由于MVC-NMF不仅有最小体积约束，还有重建数据与观测数据之间的最小误差约束，因此，MVC-NMF在不同性质的高光谱混合像元分解中更具有适用性。
实验三、测试算法的抗噪声干扰性能
改变所加噪声的强度，来比较各算法的抗噪声性能。实验数据的信噪比SNR从15dB升至40dB，每5dB为一档。像元个数固定为N＝1000，图4是四种算法的SAD随SNR变化而变化的趋势。
从图4可以看出：随着信噪比的增加，所有算法的效果在整体上都有所提升。在四种算法中，两种基于几何学的算法表现最好，抗噪声性能较强；MVC-NMF受噪声影响较大，直到图像信噪比较高(SNR>30dB)时，才表现出较好的分解效果。
实验四、测试算法性能和像元个数的关系
在这个试验中我们改变模拟数据的大小，以研究算法性能和像素数的关系。实验数据像元大小范围从40ⅹ40到100ⅹ100，数据信噪比固定为SNR＝30dB，图5是四种算法的SAD随像元个数N变化而变化的情况。
从图5可以看出，随着像元个数的增加，所有算法的性能均有小幅的改善。尤其是基于几何学的两种算法，这是由于更多的像元能形成更完整的单形体，使顶点的寻找更加精确。而对于其他算法，像元个数对它们的影响不大。在计算时间上，像元越多，计算时间越长。
计算时间上，取SNR＝25dB，端元个数p＝5，像元个数N＝1000时，在存在纯净像元的图像数据中计算，得知VCA运行时间最短，MVC-NMF最长，如表2.2所示。
表2.2四种算法运行时间比较
(3)真实数据实验
真实数据采用了1997拍摄的美国Cuprite地区的AVIRIS高光谱数据，原数据像素数较大，为提高运算速度，本文实际试验中选取其中250ⅹ191像素大小的图像。AVIRIS是美国陆地探测卫星系统，其成像范围覆盖了从红外到可见光的不同波长范围，它所采集的高光谱数据目前已广泛应用于高光谱混合像元分解研究中。该成像系统的主要参数如表2.3。
表2.3 AVIRIS的主要参数[36]
该地区位于美国内华达州的沙漠中，区域内主要为裸露的矿物，并且，各种矿物之间的混合现象较为普遍，很适合用来检验算法对高混合度数据的适应能力。在进行分解之前，有36个波段因为信噪比太低或为水汽吸收波段而被移除(波段号为1～2，104～113，148～167，221～224)，留下188个波段进行进一步处理。图6是Cuprite地区数据图(引自文献，其中图6是Cuprite地区的实地探测地物分布图，其中红色方框范围为参与本文实验的图像区域；图6(b)是本文实验使用的截取了部分区域的数据图第95波段的灰度图。图6所示的数据图与实地探测地物分布图相比较，可以确定这些端元各自对应的矿物。为了进一步衡量算法性能，将USGS库中对应矿物光谱作为参考光谱，并求取解混结果与它们之间的光谱角。
用两种方法对真实数据进行混合像元分解，分别是监督类的VCA+FCLS方法和非监督类的MVC-NMF算法。在进行解混前，首先需要设置图像中端元数目，通过常用的HFC端元估计法估计出截取的Cuprite地区图像中共包含9个端元，因此，端元个数设为9，其次，在VCA+FCLS方法中先用顶点成分法VCA对整幅图像进行端元提取，然后用全约束最小二乘法FCLS进行丰度估计，MVC-NMF算法可同时进行端元提取和丰度反演两个步骤。实验结果如图7和图8所示，每张图像都是某种端元的分布范围，其颜色从0到1分别表示某端元在像元中占的比例。将提取出的端元的光谱曲线和光谱库中的曲线对比可知，9种端元(矿物质)分别为：明矾石(Alunite)、蓝线石(Dumortierite)、黄钾铁矾(Jarosite)、高岭石(Kaolinite)、蒙脱石(Montmorillonite)、白云母(Muscovite)、绿脱石(Nontronite)、镁铝石(Pyrope)和叶腊石(Pyrophyllite)。
从丰度图7可以看出，VCA+FCLS可以将大部分混合像元分解开来，如明矾石、黄钾铁矾、叶腊石等能清晰地分离出来，与真实分布符合；但是其中有三种矿物质：蒙脱石、绿脱石和镁铝石却不能清晰分离；另外，图7(d)中估计高岭石的分布多于真实分布，错将其他物质当作了高岭石。这种方法是在VCA提取端元的基础上进行丰度反演，其反演精度受VCA端元提取结果影响较大，如不能提取出准确的端元，则丰度分布肯定不准确。总的来说，这种分解方法虽然快速，但分解效果欠佳。
图8是用非监督类分解方法MVC-NMF进行分解的丰度图。可以看出，MVC-NMF分解的效果优于VCA+FCLS，各矿物质均被分离出来。但是，可以看出物质间重叠的部分较多，尤其是图8(b)和图8(i)，也就是说，部分混合像元的丰度值不准确。
以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，根据本发明的技术方案及其发明构思加以等同替换或改变，都应涵盖在本发明的保护范围之内。
Claims (9)
1.监督与非监督的高光谱混合像元分解方法，包括监督类混合像元分解算法和非监督类混合像元分解算法，其特征在于：所述监督类混合像元分解算法包括端元提取算法和丰度估计算法，所述端元提取算法包括纯净像元指数法、内部最大体积法、顶点成分分析法、单形体增长分析算法、顺序最大角凸锥和分裂增广拉格朗日法，所述丰度估计算法包括最小二乘法、基于端元投影向量的算法和基于单形体体积的算法；所述非监督类混合像元分解算法包括独立成分分析法和非负矩阵分解法。
2.根据权利要求1所述的监督与非监督的高光谱混合像元分解方法，其特征在于：所述纯净像元指数法为单行体的形状特点决定了其在特征空间中任意直线上的投影为线段，且线段的端点必为单形体顶点在这个方向上的投影，基于这一性质，可在特征空间中随机生成若干向量，将所有像元点投影到各个向量上，直线上所有投影点中最靠外的两个便是端元点的投影，随向量方向的不断变化，记录下图像中每个像元作为极值点的次数，最后认为出现频率最高的点就是端元。
3.根据权利要求1所述的监督与非监督的高光谱混合像元分解方法，其特征在于：所述内部最大体积法通过寻找具有最大体积的单形体自动获取图像中的所有端元，算法首先进行MNF变换，使数据降至p-1维，图像Y中由端元矩阵M的p个端元构成的几何体的体积可表示为
其中mi为降维后第i个端元对应的p-1维列向量，根据体积的定义式(2.1)和式(2.2)，在理想状态下，由p个端元构成的单形体体积必然是所有由等数量像元构成的单形体中最大的。首先随机选择n个像元作为端元并计算体积，然后将其中一个端元用其他的像元替换并再次计算体积。如果替换后体积增大，则接受此替换，否则放弃此替换，再用新的像元替换，直到对每个端元都用所有像元替换过。在此过程中，端元单形体的体积是不断增大的，最后得到的必然是最大体积的单形体，而构成此单行体的像元就是端元。
4.根据权利要求1所述的监督与非监督的高光谱混合像元分解方法，其特征在于：所述顶点成分分析法通过反复寻找正交向量并计算图像矩阵在正交向量上的投影距离来逐一提取端元，单形体的若干个顶点可以张成一个子空间，而单形体在某个与这个子空间正交的向量上的投影距离的最大值点一定是单形体的顶点，可以先设法找到一个初始端元，然后每次循环都先找一个和已经找到的端元同时正交的单位向量，再将所有向量点投影到这一单位向量上，并将投影结果的最大值记为新端元，加入端元集合并开始下一次循环，直到找到设定的端元个数为止。
5.根据权利要求1所述的监督与非监督的高光谱混合像元分解方法，其特征在于：所述单形体增长分析算法是在N-FINDR算法基础上做了改进：与N-FINDR同时搜索所有的端元不同，单形体增长分析算法采用逐步增加单形体定点的方法搜索各个端元，在搜索过程中单形体的体积是膨胀的，其维度也逐渐增长，最后得到端元，其计算复杂度较N-FINDR要小；所述顺序最大角凸锥是对顶点成分分析法的改进，将顶点成分分析法的正交投影放宽为斜交投影，以此为代价，可以换取满足“非负”约束的端元丰度，此法可同时实现端元提取和丰度反演。
6.根据权利要求1所述的监督与非监督的高光谱混合像元分解方法，其特征在于：所述分裂增广拉格朗日法在使用前之前需要根据端元数目p对图像进行降维，使得L＝p，得到降维后的图像矩阵Ypxn。当L＝p时，端元矩阵M为方阵且各个端元线性无关，从而存在Q＝M-1，寻找最小体积单形体可以描述为如下的优化问题
那么提取端元M的问题就转化成了求解Q的问题。
7.根据权利要求1所述的监督与非监督的高光谱混合像元分解方法，其特征在于：所述最小二乘法是目前为止应用最为广泛的丰度反演方法，根据满足丰度约束条件的程度可以分为四种不同的最小二乘法：无约束最小二乘法，“和为1”约束最小二乘法，“非负”约束最小二乘法和全约束最小二乘法；
根据线性混合模型公式，不考虑误差时，可得线性混合模型公式为
当不考虑对丰度的约束时，可以求得第1.2.1节中式(1.4)的无约束解为
sUCLS＝(MTM)-1MTy (2.5)
考虑到ASC，可得出式的“和为1”约束解为
其中Ip为p阶单位矩阵；1＝[1,1,...1]T为全1(p维)列向量；
考虑到描述的ANC，可得出非负约束解SNCLS，但此解不能由简单的算子表示出来，而是利用迭代方法来获得最优解，非负约束最小二乘法的两个迭代关系式如下：
同时满足ANC和ASC的解成称为全约束最小二乘解，可将式(2.6)代入式(2.7)求解得到；
所述基于端元投影向量的算法
根据线性光谱混合模型和几何学描述，混合像元位于单形体的L维特征空间的内部，端元位于单形体的顶点，某个端元是距离其他端元构成的L-1维空间超平面最远的点，可以根据像元到超平面的距离与端元到超平面的距离之间的比值来计算像元中端元占的比例，即丰度。
所述基于单形体体积的算法
根据线性光谱混合模型和几何学描述，某个端元用其他像元替换后所得单形体的体积比原单形体体积小，可以根据用某像元替换后的单形体与原单形体的体积比计算出该像元中端元的丰度。
8.根据权利要求1所述的监督与非监督的高光谱混合像元分解方法，其特征在于：所述独立成分分析法一般是假设形成遥感图像的几个端元的丰度是互相独立的，通过最小化它们之间的相关信息来达到解混的目的，假定每一个端元的丰度都是一个随机信号源，这样就可以把高光谱遥感图像看成是多源混合信号，其中的光谱解混问题就转化为一个盲信号分离的问题，独立成分分析法是解决混合像元分解的一种有效方法，可以同时解决“端元提取”和“丰度估计”两个问题。
9.根据权利要求1所述的监督与非监督的高光谱混合像元分解方法，其特征在于：所述非负矩阵分解法其每步的迭代均在上一步的结果上使用乘法得到，由于乘法整数前后的符号是不变的，只要初始值非负就能保证结果非负；
给定一非负矩阵Y∈Rn×m和一个正整数r＜min(m,n)，NMF的目标是找到两个非负矩阵W∈Rn×r和H∈Rr×m，使其满足
Y≈WH (2.8)
目标函数有两种：第一种是Y和WH之间的欧式距离
它的下界是0，当Y＝WH时达到此下界；
第二种是Y和WH之间的“散度”，也即概率意义上的距离
同样，当Y＝WH时达到此下界；
就欧式距离情况下推导得到的乘法迭代公式如下：
W←W.*(YHT)./(WHHT) (2.11)
H←H.*(WTY)./(WTWH) (2.12)。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111045548.4A CN113743325B (zh) | 2021-09-07 | 2021-09-07 | 监督与非监督的高光谱混合像元分解方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111045548.4A CN113743325B (zh) | 2021-09-07 | 2021-09-07 | 监督与非监督的高光谱混合像元分解方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN113743325A true CN113743325A (zh) | 2021-12-03 |
| CN113743325B CN113743325B (zh) | 2024-01-12 |
Family
ID=78736831
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202111045548.4A Active CN113743325B (zh) | 2021-09-07 | 2021-09-07 | 监督与非监督的高光谱混合像元分解方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN113743325B (zh) |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20100130870A1 (en) * | 2007-04-25 | 2010-05-27 | Ivica Kopriva | Method for real time tumour visualisation and demarcation by means of photodynamic diagnosis |
| CN108470192A (zh) * | 2018-03-13 | 2018-08-31 | 广东工业大学 | 一种高光谱分类方法及装置 |
| CN109085131A (zh) * | 2018-07-12 | 2018-12-25 | 重庆邮电大学 | 基于丰度稀疏和端元正交性约束nmf的高光谱解混方案 |
| US20190392261A1 (en) * | 2016-05-04 | 2019-12-26 | Shandong University | End-member extraction method based on segmented vertex component analysis (vca) |
-
2021
- 2021-09-07 CN CN202111045548.4A patent/CN113743325B/zh active Active
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20100130870A1 (en) * | 2007-04-25 | 2010-05-27 | Ivica Kopriva | Method for real time tumour visualisation and demarcation by means of photodynamic diagnosis |
| US20190392261A1 (en) * | 2016-05-04 | 2019-12-26 | Shandong University | End-member extraction method based on segmented vertex component analysis (vca) |
| CN108470192A (zh) * | 2018-03-13 | 2018-08-31 | 广东工业大学 | 一种高光谱分类方法及装置 |
| CN109085131A (zh) * | 2018-07-12 | 2018-12-25 | 重庆邮电大学 | 基于丰度稀疏和端元正交性约束nmf的高光谱解混方案 |
Non-Patent Citations (2)
| Title |
|---|
| 卓莉;曹晶晶;王芳;陶海燕;郑;: "采用目标端元修正的高光谱混合像元盲分解", 遥感学报, no. 02 * |
| 宋晓瑞;吴玲达;孟祥利;: "利用稳健非负矩阵分解实现无监督高光谱解混", 中国图象图形学报, no. 04 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN113743325B (zh) | 2024-01-12 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Ding et al. | Semi-supervised locality preserving dense graph neural network with ARMA filters and context-aware learning for hyperspectral image classification | |
| CN110728192B (zh) | 一种基于新型特征金字塔深度网络的高分遥感图像分类方法 | |
| Yuan et al. | Factorization-based texture segmentation | |
| CN110348399B (zh) | 基于原型学习机制和多维残差网络的高光谱智能分类方法 | |
| Velasco-Forero et al. | Improving hyperspectral image classification using spatial preprocessing | |
| Zhang et al. | A pixel shape index coupled with spectral information for classification of high spatial resolution remotely sensed imagery | |
| CN110084159A (zh) | 基于联合多级空谱信息cnn的高光谱图像分类方法 | |
| CN107992891B (zh) | 基于光谱矢量分析多光谱遥感图像变化检测方法 | |
| Harder et al. | Spectraldefense: Detecting adversarial attacks on cnns in the fourier domain | |
| CN110533077B (zh) | 用于高光谱图像分类的形状自适应卷积深度神经网络方法 | |
| Mishne et al. | Graph-based supervised automatic target detection | |
| CN105913092B (zh) | 基于子空间学习的图正则高光谱图像波段选择方法 | |
| CN107895139A (zh) | 一种基于多特征融合的sar图像目标识别方法 | |
| CN111680579B (zh) | 一种自适应权重多视角度量学习的遥感图像分类方法 | |
| CN104751181A (zh) | 一种基于相对丰度的高光谱图像解混方法 | |
| Fu et al. | A novel spectral-spatial singular spectrum analysis technique for near real-time in situ feature extraction in hyperspectral imaging | |
| Ma et al. | Hyperspectral anomaly detection based on low-rank representation with data-driven projection and dictionary construction | |
| Minh et al. | Approximate log-Hilbert-Schmidt distances between covariance operators for image classification | |
| CN110310263B (zh) | 一种基于显著性分析和背景先验的sar图像居民区检测方法 | |
| CN105957112A (zh) | 基于快速uncls的高光谱亚像素探测方法 | |
| Sha et al. | Semi-supervised classification for hyperspectral images using edge-conditioned graph convolutional networks | |
| CN109543546B (zh) | 基于深度序分布回归的步态年龄估计方法 | |
| Schneiderman et al. | A histogram-based method for detection of faces and cars | |
| Zhou et al. | Alternating direction iterative nonnegative matrix factorization unmixing for multispectral and hyperspectral data fusion | |
| CN112784777A (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 |