一种基于数据空间正交基的遥感图像混合像元分解方法

专利名称：一种基于数据空间正交基的遥感图像混合像元分解方法
技术领域：
本发明属于遥感图像处理技术领域，具体涉及一种基于数据空间正交基的遥感图像混合像元分解方法。
背景技术：
遥感是本世纪六十年代发展起来的新兴综合技术，与空间、电子光学、计算机、地理学等科学技术紧密相关，是研究地球资源环境的最有力的技术手段之一。近年来，随着成像技术的进步，多波段遥感图像在越来越多的领域得到了广泛应用。由于成像系统空间分辨率的限制和地表的复杂多样，所获得的遥感图像中的一个像元往往包含着多种地物类型，这就形成了混合像元。如何从混合像元广泛存在的多波段遥感图像中准确的提取端元信号，并有效的对混合像元进行分解，已成为了遥感图像定量分析的一个重要研究课题[1]。
单形体几何学方法是目前遥感图像混合像元分解问题中十分重要的一类方法，该类方法运用遥感图像数据的几何特性进行端元的提取。近年来，一些基于该方法的遥感图像端元提取算法被提出并得到了广泛研究和应用，它们中具有代表性的有N-FINDR[2]、SGA[3]等，它们的一些主要缺点为算法最终结果易受初始值选择影响；计算复杂度高，运算十分耗时；需要另外的一些方法来帮助确定端元个数以及完成丰度解混工作。
针对以上问题，在遥感图像混合像元分解的研究中，如何快速准确的确定端元个数，提取出端元并获得各端元地物的分布情况已成为目前研究的热点。
下面介绍与本发明相关的一些概念1、线性光谱混合模型近年的研究中，线性光谱混合模型被广泛的应用于遥感图像中的混合像元分解问题，该模型假设图像中的每个像元都为各个端元像元通过线性混合得到。设X为多通道遥感图像中单一像元的多光谱或高光谱矢量，A为由各类纯地物信号(端元)的多光谱或高光谱矢量所组成的反射特性矩阵，S为该像元中各类地物所占的百分比(即丰度)，N为模型的误差，则依此模型有如下关系式X＝AS+N.
若遥感图像有n个通道，其中有m类地物类型，则式中X为n×1的向量，A为n×m的矩阵，S为m×1的向量，N为n×1的向量。对于高光谱遥感图像，一般有n＞＞m。
同时，基于混合像元分解问题的实际物理意义，S应满足如下两个约束条件1)混合像元中各成分的比例si之和应该等于1，即Σi=1msi=1.]]>2)分解所得各成分的比例si应该在
的范围内，即0≤si≤1，(i＝1，2，...，m).
2、高维空间中关于单形体的相关定义和定理定义1定义n维空间中由原点和k个点α1，α2，...，αk构成的k维凸平行多面体有向体积为k重向量α1α2...αk，由原点到这k个点形成的k个矢量α1，α2，...，αk也即为该k维凸平行多面体的支撑棱。
定义2n维空间中由原点和k个点α1，α2，...，αk构成的k维凸单形体体积为V(E)=1k!E,]]>其中E为k重向量α1α2...αk的模|α1α2...αk|，它通过Gram行列式计算得到E=α1·α1α1·α2...α1·αk............αk·α1αk·α2...αk·αk12]]>定理1Gram行列式的正交向量计算法假设有一组线性无关向量组α1，α2，...，αk，将其正交化为向量组β1，β2，...，βk，则有 其中
β1＝α1β2=α2-α2·β1β1·β1β1]]>......
βk=αk-αk·βk-1βk-1·βk-1βk-1-...αk·β1β1·β1β1]]>发明内容本发明的目的在于提出一种基于数据空间正交基的遥感图像混合像元分解方法。该发明可以快速有效的从多通道遥感图像中提取出各个端元信号，并能确定出合适的端元个数以及进行高精度的丰度解混工作。
本发明提出的遥感图像混合像元分解方法，具体内容如下顺序的搜索端元像元，在每一步搜索过程中，图像内所有像元中作为新的顶点与已得到的端元像元构成的单形体具有最大的体积的像元为该步得到的新端元像元；其中，将基于行列式的单形体体积计算等价于一组正交基的模的乘积运算，并且在将端元形成的支撑棱正交化为数据空间的正交基运算中引入了递推概念；由所得到的正交基模的单调下降性质，自行确定合适的端元个数；由所得到的正交基两两正交性质，可以同时获得各端元成分的丰度解混结果。该解混方法具有良好的抗噪声特性。具体计算步骤如下1、端元提取(顺序的搜索端元像元)。
对于一幅由线性光谱混合模型描述的n波段遥感图像，在无噪声环境下，其所有像元在n维空间中正好构成了一个m-1维的单形体，而端元则位于这个单形体的顶点上；该m个端元通过顺序的搜索得到，首先，初始化选取合适的第一端元，然后，在后续的每一步搜索过程中，图像内所有像元中作为新的顶点与已得到的端元像元构成的单形体具有最大体积的像元为该步得到的新端元像元。
对上述端元提取方法进一步说明如下由背景技术中的定理1可获得如下结论对于一个由k条支撑棱撑起的凸单形体体积计算，可以先将这k条支撑棱所表示的矢量正交化为一组正交向量，然后直接计算这组向量模的乘积即可得到该单形体的体积。由于该定理所计算的是k个点和原点形成的单形体体积，而实际计算的是由k个点独立形成的单形体体积，因此实际应用时先要初始化选取第一个端元矢量e0作为上面定义定理中的原点。在已提取k个端元时，将同时获得k-1条支撑棱α1，α2，...，αk以及由其正交化得到的k-1个正交向量β1，β2，...，βk(这里是k-1而不是k，因为第一个端元实际上等于做为原点处理)，在第k+1步搜索时，图像中所有像元中与已提取到的k个端元构成的单形体具有最大体积的被取为第k+1端元，由定理1，这里的单形体体积计算被等价于一组正交向量模乘积计算，对于图像中不同的像元构成的不同单形体体积计算，其前k-1个正交向量β1，β2，...，βk都是相同的已知量，因此对于图像中的每一像元p，其形成的新支撑棱为αk＝p-e0，将αk正交化为新的正交向量βk，则具有最大模的βk所对应的像元即为提取的第k+1端元。具体而言，本发明将基于行列式的单形体体积计算等价于一组正交基的模的乘积运算，并且在将端元形成的支撑棱正交化为数据空间的正交基运算中引入了递推概念，具体步骤如下n维空间中对于一个k维凸单形体，如果其支撑棱为α1，α2，...，αk，则其体积V(E)=1k!α1·α1α2·α2...α1·αk............αk·α1αk·α2...αk·αk12=1k!|β1|·|β2|···|βk|,]]>其中β1，β2，...，βk由α1，α2，...，αk通过下述式子得到β1＝α1β2=α2-α2·β1β1·β1β1]]>......
βk=αk-αk·βk-1βk-1·βk-1βk-1-...αk·β1β1·β1β1.]]>对于βk的更新计算，每个去相关性的减法运算以及αk在不同的搜索步骤k中都是不变的，在每次搜索步骤完成时将其保留以方便下一搜索过程的计算，分别由下式中γki和li实现已提取端元为e0，e1，...ek，其形成的支撑棱和正交基分别为α1，α2，...，αk，β1，β2，...，βk，以及图像中每个像元pi的γki，取ek+1=argmaxpi(|γki(pi)|),]]>并取βk+1=argmaxγki(|γki|),]]>对于图像中每个像元pi，更新γk+1i=γki-Ii·βk+1βk+1·βk+1βk+1.]]>2、端元个数确定已往的端元提取方法需要预先确定出参数m，比较常用的方法是PCA、MNF[5]等特征值分析以及近年来Chang等人提出的一种基于假设检验的VD方法[1][6]。本发明由于采用的是顺序的搜索过程，因此并不需要预先确定端元个数。在端元提取的同时，本发明还会获得由这些端元张成的子空间的一组正交基，理论上可以证明随着端元的递推提取，这组基的模是单调下降的。对于一个无噪声的理想环境，其将会最终降为0；对于一实际环境，其会下降到一个变化平稳的小量。因此，可以在端元提取的同时观察正交基模的变化曲线，当其下降到变化平稳且为一小量的部分时即可终止算法运行。
3、丰度解混不同于传统的基于最小二乘或数值优化方法[7]，本发明将每一混合像元所形成支撑棱向k个正交基分别投影，然后通过求解一具有唯一解的线性方程组，获得k个端元各自对应的丰度，第1个端元对应的丰度由1减去其余k个端元对应的丰度和所得到。设待分解混合像元为p，提取的端元矢量为e0，e1，...，ek，其在p中分别对应的丰度为p0，p1，p2，...，pk，得到支撑棱为α1，α2，...，αk，正交基为β1，β2，...，βk，则所求解的线性方程组为βk·βk00...0αk·βk-1βk-1·βk-10...0...............αk·β1αk-1·β1αk-2·β1...β1·β1pkpk-1...p1=(p-e0)·βk(p-e0)·βk-1...(p-e0)·β1]]>p0=1-Σi=1kpi]]>本发明的优点本发明为一种基于数据空间正交基的遥感图像混合像元分解方法。其优点在与分步搜索端元方法的使用确保了可以获得固定的端元提取结果；递推正交基的计算显著提升了方法的计算速度；可以自行确定合适的端元个数；可以同时获得准确的丰度解混结果。本发明在基于多光谱和高光谱遥感图像的高精度的地物分类以及地面目标的检测和识别方面具有重要意义。


图1为本发明与两种传统的基于单形体的端元提取方法N-FINDR，SGA计算时间复杂度比较。
图2模拟遥感图像示意图。其中，(a)25个包含测试像元的方格示意图，(b)往(a)中添加背景像元以及信噪比(SNR)为20db高斯白噪声形成的模拟遥感图像。
图3中，(a)端元提取结果，(b)正交基模单调下降曲线。
图4中，(a)为美国Cuprite地区AVIRIS遥感图象，(b)为从该数据中提取到的端元位置示意图。
图5中，(a)为得到的正交基模单调下降曲线，(b)为作为比较的MNF变换得到的特征值曲线。
图6为得到的五种典型矿物丰度解混结果，其中(a)Buddingtointe (b)Muscovite (c)Calcite (d)Alunite (e)Kaolinite。
图7中(a)为美国Indiana地区AVIRIS遥感图象，(b)为从该数据中提取到的端元位置示意图。
图8中，(a)为得到的正交基模单调下降曲线，(b)为作为比较的MNF变换得到的特征值曲线。
图9中为得到的七种典型端元丰度解混结果，其中(a)谷物 (b)树林 (c)铁塔 (d)干草堆 (e)大豆 (f)公路 (g)植被。
具体实施例方式
1.端元提取1)初始化对于图像中所有像元pi(上标i表示所有不同的N个像元)a)选取所有像元矢量中具有最大模的像元，令其为e0b)选取所有像元矢量中离e0最远的像元，令其为e1，计算α1＝e1-e0，β1＝α1c)对于所有像元pi，计算li=pi-e0,γ1i=Ii-Ii·β1β1·β1β1,(i=1,...,N)]]>2)当k＞1时，已提取端元为e0，e1，...，ek，其形成的支撑棱和正交基为α1，α2，...，αk，β1，β2，...，βk，以及图像中每个像元pi的γki，搜索满足下式的pi，令其为ek+1ek+1=argmaxpi(|γki(pi)|)]]>并取βk+1=argmaxγki(|γki|)---(9)]]>3)对于图像中每个像元pi，更新γk+1i=γki-Ii·βk+1βk+1·βk+1βk+1---(10)]]>
返回2)2.端元个数确定由得到的β1，β2，...，βk绘出其模下降曲线，观察正交基模的变化曲线，当其下降到变化平稳且为一小量的部分时可终止上述算法的运行。
3.丰度解混当端元全部提取完毕后，设待分解混合像元为p，所提取的端元矢量为e0，e1，...，ek，其在p中分别对应的丰度为p0，p1，p2，...，pk，得到支撑棱为α1，α2，...，αk，正交基为β1，β2，...，βk，求解线性方程组βk·βk00...0αk·βk-1βk-1·βk-10...0...............αk·β1αk-1·β1αk-2·β1...β1·β1pkpk-1...p1=(p-e0)·βk(p-e0)·βk-1...(p-e0)·β1]]>计算p0=1-Σi=1kpi.]]>由上式左边可以看出该线性方程组的系数矩阵为一下三角矩阵，因此其求解的速度是十分快速的，同时其右边的非齐次项是由混合像元得到的支撑棱分别往k个正交基的方向投影得到，因此对于噪声中正交于端元矢量张成子空间的成分，其投影将为0，也就是说该解混方法可以自动去除噪声中正交于数据空间的部分。
下面，我们分别以模拟和实际遥感图像数据为例说明具体的实施方式1.模拟遥感图像数据图2(a)为实验所用模拟图像示意图，其大小为200×200，所有像元均由Alunite(A)、Buddingtointe(B)、Calcite(C)、Kaolinite(K)和Muscovite(M)共5种矿物以不同丰度混合而成。在这幅图中有5×5共25个不同大小的方格做为测试像元，其中第一列的5个方格各为4×4像元大小，每个方格对应一种矿物的纯像元；第二列的5个方格各为2×2像元大小，同样分别对应五种矿物的纯像元；第三列方格大小为2×2像元，第四、第五列的方格均为一个像元大小，这三列方格的像元均为混合像元，混合的丰度如表1和表2所示。除了这25个方格中像元外的所有像元均为背景像元(BKG)，它们由五种矿物均匀混合而成，也即BKG＝20％A+20％B+20％C+20％K+20％M。最后，模拟的遥感图象中被加入了信噪比(SNR)为20db的高斯白噪声形成最终实验所使用的模拟遥感图像如图2(b)所示。
表1.第3列方格中所包含的像元混合丰度 表2.第4及第5列方格中所包含的像元混合丰度


将本发明方法应用于图2所示的模拟遥感图像，获得的端元提取结果以及正交基模值下降曲线分别如图3(a)和(b)所示，从图3(a)可以看到，所发明方法准确的提取出了图像中包含的所有5种端元，由图3(b)所显示的正交基模值下降曲线可以看到，第5个正交基的模值已下降到约等于零，因此可以确定原模拟图像中的端元个数为5个。
表3.所提议方法丰度解混结果

同时我们将表2中带*的像元p5，111做为测试像元对所发明的丰度解混方法进行了验证，该像元被人为添加了信噪比为6分贝的正交于数据空间的噪声成分。实验结果如表3所示。表3中的每行表示5种不同矿物各自的丰度。三列中的左列为丰度真实值，中间列为所发明方法解混结果，右列为最小二乘(LS)法的解混结果。从表中列出的结果可以看出所发明方法方法可以自动去除噪声成分中正交于数据子空间的部分，分解精度明显优于最小二乘方法。
2.Cuprite地区的AVIRIS数据实验使用ENVI软件自带的Cuprite地区的AVIRIS数据如图4(a)，图像大小为400×300，波长范围是1.99-2.48μm，光谱分辨率为10nm，共有172-221波段间的50个波段数据。该地区位于美国内华达州南部，其地表基本无植物覆盖，多为裸露矿物。实地勘测表明该地区广泛分布的矿物主要为下列五种Alunite(A)，Buddingtointe(B)，Calcite(C)，Kaolinite(K)，Muscovite(M)，文献[8]提供了其实地勘测的分布图。
图5(a)为本发明方法得到的正交基模值下降曲线，图5(b)为对该数据做MNF变换得到的特征值图。MNF变换通过计算大于1的特征值个数来确定源个数，由图5(a)和(b)均可知，该地区端元个数约为11-13个左右。将本发明的端元提取方法算法用于该实际数据提取出13个端元，将结果与U.S.Geological Survey(USGS)光谱库进行对比可知，所提取的端元中第3、第4、第6、第7、第10端元分别对应Kaolinite(K)，Alunite(A)，Calcite(C)，Muscovite(M)，Buddingtointe(B)五种典型矿物，其在原图像中的位置如图4(b)所示，再应用所发明的丰度解混方法对该数据进行丰度解混，获得的五种典型矿物各自的分布图如图6所示，解混结果与实地勘测结果有很好的吻合。
表(4)中给出了分别使用最佳N-FINDR，SGA以及我们提出的基于数据空间正交基的端元提取方法对该数据进行端元提取所需要的计算时间对比。计算环境为CPUIntel(R)Pentium(R)M Processor 1.60GHz；Memory1GBytes；OSWindows XP Matlab 7.0.
表4.N-FINDR，SGA以及所提议算法各自运算时间对比

3.Indiana地区的AVIRIS数据该部分实验中我们使用的是成像于1992年7月的一幅AVIRIS高光谱遥感数据。该数据包含0.4-2.5μm共220个波段的数据，光谱分辨率为10nm，空间分辨率为17m，大小为145×145(共21025pixels)。该地区地表的主要覆盖类型有各种农作物(包括玉米、大豆、小麦等)、植被(包括草地、树林等)、以及各种人工建筑(高速公路、铁塔、房屋等)。取第70、86、136波段分别作为R、G、B分量合成伪彩色图如图7(a)所示。该数据由美国Purdue大学提供网上下载[9]，同时，该研究组也给出一份该地区实地勘测结果可供参考[10]，它在将不同土壤开垦情况下的同一作物看成不同类型，忽视土壤、部分植被等背景以及一些小目标的情况下将该成像区域划分为16类。在我们的实验分析前，该数据的第1-4，78-82，103-115，148-166以及211-220波段由于水吸收波段或很低的信噪比已被事先舍弃，因此，剩下的总共169个波段数据被用于方法验证工作。
图8(a)为所提议算法得到的正交基内积值下降曲线，图8(b)为对该数据做MNF变换得到的特征值图，由图8可知该地区可能的端元个数在12-14个左右。所提议算法提取的7个典型端元位置如图7(b)所示，其中各字母代表的地物成分为a谷物 b树林 c铁塔d干草堆 e大豆 f公路 g植被。丰度解混得到的各端元在该区域对应的分布情况如图9所示。将图9的解混结果与实地调查的情况[10]比较，可以看出，解混结果与实地调查结果非常吻合[10]。
表2中同样给出了对该数据分别使用最佳N-FINDR，SGA以及我们提出的基于数据空间正交基的端元提取方法进行端元提取所需要的计算时间对比。
表2.N-FINDR，SGA以及所提议算法各自运算时间对比

参考文献[1]C.-I Chang，Hyperspectral ImagingTechniques for Spectral Detection and Classification.New YorkPlenum，2003. M.E.Winter，“N-FINDRAn algorithm for fast autonomous spectral end-memberdetermination in hyperspectral data，”in Proc.SPIE Conf.Imaging Spectrometry V，1999，pp.266-275. Chein-I Chang，Chao-Cheng Wu，Wei-min Liu，and Yen-Chieh Ouyang，“A New GrowingMethod for Simplex-Based Endmember Extraction Algorithm”，IEEE Transactions onGeoscience and Remote Sensing，Vol.44，No.10，October 2006[4]沈文选，“单形论导引三角形的高维推广研究”，湖南师范大学出版社，2000[5]Andrew A.Green，Mark Berman，Paul Switzer，and Maurice D.Craig，“A Transformationfor Ordering Multispectral Data in Terms of Image Quality with Implications for NoiseRemoval”，IEEE Transactions on Geoscience and Remote Sensing，Vol.26，No.1，January1988[6]C.-I Chang and Q.Du，“Estimation of number of spectrally distinct signal sources inhyperspectral imagery，”IEEE Transactions on Geoscience and Remote Sensing，Vol.42，No.3，March 2004[7]Daniel C.Heinz，and Chein-I Chang，“Fully Constrained Least Squares Linear SpectralMixture Analysis Method for Material Quantification in Hyperspectral Imagery”，IEEETransactions on Geoscience and Remote Sensing，Vol.39，No.3，March 2001 Http://speclab.cr.usgs.gov/cuprite.html[9]Http://cobweb.ecn.purdue.edu/～biehl/Multispec/documentation.html[10]D.Landgrebe，“Multispectral data analysisA signal theory perspective，”School of Electr.Comput.Eng.，Purdue Univ.，West Lafayette，IN，1998.
权利要求
1.一种基于数据空间正交基的遥感图像混合像元分解方法，其特征在于顺序的搜索端元像元，在每一步搜索过程中，图像内所有像元中作为新的顶点与已得到的端元像元构成的单形体具有最大的体积的像元为该步得到的新端元像元；其中，将基于行列式的单形体体积计算等价于一组正交基的模的乘积运算，并且在将端元形成的支撑棱正交化为数据空间的正交基运算中引入了递推概念；由所得到的正交基模的单调下降性质，自行确定合适的端元个数；由所得到的正交基两两正交性质，可以同时获得各端元成分的丰度解混结果。
2.根据权利要求1所述的基于数据空间正交基的遥感图像混合像元分解方法，其特征在于所述顺序的搜索端元像元的步骤如下对于一幅由线性光谱混合模型描述的n波段遥感图像，在无噪声环境下，其所有像元在n维空间中正好构成了一个m-1维的单形体，而端元则位于这个单形体的顶点上；该m个端元通过顺序的搜索得到，首先，初始化选取合适的第一端元，然后，在后续的每一步搜索过程中，图像内所有像元中作为新的顶点与已得到的端元像元构成的单形体具有最大体积的像元为该步得到的新端元像元。
3.根据权利要求2所述的基于数据空间正交基的遥感图像混合像元分解方法，其特征在于将基于行列式的单形体体积计算等价于一组正交基的模的乘积运算，并且在将端元形成的支撑棱正交化为数据空间的正交基运算中引入了递推概念，具体步骤如下n维空间中对于一个k维凸单形体，如果其支撑棱为α1，α2，...，αk，则其体积V(E)=1k!α1·α1α1·α2···α1·αk············αk·α1αk·α2···αk·αk12=1k!|β1|·|β2|···|βk|,]]>其中β1，β2，...，βk由α1，α2，...，αk通过下述式子得到β1＝α1β2=α2-α2·β1β1·β1β1]]>……βk=αk-αk·βk-1βk-1·βk-1βk-1-···αk·β1β1·β1β1]]>对于βk的更新计算，每个去相关性的减法运算以及αk在不同的搜索步骤k中都是不变的，在每次搜索步骤完成时将其保留以方便下一搜索过程的计算，分别由下式中γki和li实现已提取端元为e0，e1，...ek，其形成的支撑棱和正交基分别为α1，α2，...，αk，β1，β2，...，βk，以及图像中每个像元pi的γki，取ek+1=argmaxpi(|γki(pi)|),]]>并取βk+1=argmaxγki(|γki|),]]>对于图像中每个像元pi，更新γk+1i=γki-Ii·βk+1βk+1·βk+1βk+1.]]>
4.根据权利要求3所述的基于数据空间正交基的遥感图像混合像元分解方法，其特征在于在端元提取的同时观察正交基模的变化曲线，当其下降到变化平稳且为一小量的部分时即终止算法运行，得到合适的端元个数。
5.根据权利要求4所述的基于数据空间正交基的遥感图像混合像元分解方法，其特征在于所述同时得到各个端元成分的丰度解混结果的方法如下通过求解一具有唯一解的线性方程组，获得k个端元各自对应的丰度，第1个端元对应的丰度由1减去其余k个端元对应的丰度和所得到；设待分解混合像元为p，提取的端元矢量为e0，e1，...ek，其在p中分别对应的丰度为p0，p1，p2，...，pk，得到支撑棱为α1，α2，...，αk，正交基为β1，β2，...，βk，则所求解的线性方程组为βk·βk00···0αk·βk-1βk-1·βk-10···0···············αk·β1αk-1·β1αk-2·β1···β1·β1pkpk-1···p1=(p-e0)·βk(p-e0)·βk-1···(p-e0)·β1]]>p0=1-Σi=1kpr.]]>
全文摘要
本发明属于遥感图像处理技术领域，具体为一种基于数据空间正交基的遥感图像混合像元分解方法。在由数据集形成的具有最大体积的单形体中，该方法通过递推寻找该单形体的一个新顶点来确定一个新的端元。同时，在每一个端元的提取过程中，将基于行列式的单形体体积计算等价于一组正交基模乘积的计算，从而可显著提高方法的计算效率并确保本方法总能获得相同的端元提取结果。此外，本发明不仅能快速有效的完成端元提取工作，还可以同时完成端元个数确定与丰度解混两项工作，其性能优于传统的基于单形体的遥感图像混合像元分解方法。新方法在基于多光谱和高光谱遥感图像的高精度的地物分类以及地面目标的检测和识别方面具有特别重要的应用价值。
文档编号G06T7/00GK101030299SQ200710038629
公开日2007年9月5日 申请日期2007年3月29日 优先权日2007年3月29日
发明者陶雪涛, 王斌, 张立明 申请人:复旦大学