CN103761742A - 一种基于同质指数的高光谱遥感图像稀疏解混方法 - Google Patents

Abstract

一种基于同质指数的高光谱遥感图像稀疏解混方法，包括对高光谱遥感图像的每个像元进行同质分析提取同质指数，根据整幅高光谱遥图像中像元同质指数的值域对每个像元的同质指数进行线性拉伸，将其作为基于全变分的变量分离增广拉格朗日稀疏解混算法中空间正则项的权重，并重新定义算法中邻接像元间丰度的差商算子。本方法提高了稀疏解混的精度，使得解混丰度的空间平滑性更接近图像真实的空间平滑性分布，并且更有效的抑制了噪音对解混结果的影响。本方法在基于高光谱遥感图像的高精度地物分类以及地面目标检测和识别方面具有重要的应用价值。

Description

一种基于同质指数的高光谱遥感图像稀疏解混方法

技术领域

本发明属于遥感图像处理技术领域，涉及一种基于同质指数的高光谱遥感图像稀疏解混方法。

背景技术

高光谱遥感图像具有较高的光谱分辨率，能更为细致精确地分析地物的物质成分，因而得到了广泛应用。然而高光谱图像的空间分辨率一般较低，混合像元普遍存在，极大地阻碍了它的应用，因此混合像元的分解就成为推动其应用突破的一项关键技术。混合像元分解过程就是要从混合像元中识别不同类型的地物（端元），并求出其在混合像元中所占的比例（丰度），是高光谱遥感图像分析的核心问题之一。

基于稀疏回归的混合像元分解方法是一种半监督的解混方法[1-4]，它利用已有的端元光谱库[5]作为先验信息，将像元表示成端元光谱库中某些端元的线性组合。该方法不需要在图像中提取端元，也不要求图像中每个端元必须有纯净像元对应，在像元普遍混合程度较高时仍能取得较好的效果。一个像元含有的端元数目通常小于整幅图像所含有的端元数目，更远远小于端元光谱库中端元的数目，即像元用端元线性表示时丰度值具有稀疏性[2-7]。基于稀疏回归的解混方法在解混过程中加入了丰度的稀疏性约束，得到的结果更符合实际情况。

大部分混合像元分解的方法都基于图像的光谱分析，却忽略了图像的空间信息。然而图像在空间上有一定的连续性，像元与其邻近像元的端元以及端元对应的丰度有很强的相关性。近年来很多学者研究如何利用高光谱图像丰富的空间信息来提高混合像元分解的精度，提出了很多改进方法。一部分通过图像空间分析更准确地提取端元光谱[8-14]。Zortea和Plaza[12]通过对图像中像元与其邻域窗口内像元的光谱相似度分析提取空间信息，将与邻近像元光谱相似的像元称为同质像元，认为纯净像元最有可能是同质像元，所以在端元提取过程中收缩非同质像元以突出同质像元。Martin和Plaza[13]通过聚类或者分割方法自适应地提取光谱纯净区域，避免了窗口大小的选择对方法的影响，然后在空间同质区提取端元。他们进一步提出空间同质指数[14]，并将空间同质指数和光谱纯净指数融合在一起辅助提取端元。另外一些则是在原来基于光谱空间分析方法的目标优化函数里加入空间平滑性约束[15-18]。其中Iordache等[18]在稀疏解混模型中加入空间全变分正则项，使得解混结果既保持稀疏性又保持一定的空间平滑性。这些方法对不同的空间邻域建立了一致的空间平滑性模型；然而端元丰度的空间分布异常复杂，平滑性并不能保持一致。

下面介绍与本发明相关的一些概念：

1.线性光谱混合模型

线性光谱混合模型假设像元的光谱反射率是其中含有的所有地物的光谱反射率的线性组合[19]。线性光谱混合模型可以表示为

y=Mα+n   （1）

其中y表示一个像元的光谱反射率（简称光谱），是一个L维的列向量（L是高光谱遥感图像的光谱波段数）；M是这个像元中含有端元的光谱组成的一个L行V列矩阵，这个像元中含有V个端元，M的每一列表示一个端元的光谱；α是一个V维向量，它的每一个元素α_v表示第v个端元对应的丰度；n是L维的加性噪音。端元的丰度是指端元在像元中所占有的比重，因此丰度值要满足以下“和为一”约束条件（Abundance Sum-to-one Constraint,ASC）和“非负”约束条件(Abundance Nonnegative Constraint,ANC)：

$\left(\mathrm{ASC}\right):\underset{v=1}{\overset{V}{\mathrm{\Σ}}}{\mathrm{\α}}_v=1---\left(2\right)$

(ANC):α_v≥0    （3）

2.基于稀疏回归的解混模型

基于稀疏回归的解混模型将端元光谱库作为先验信息，端元光谱库是人工提取的地物光谱信息。这时图像解混就不需要进行端元提取，也不需要假设图像中每个端元都存在纯净像元。假设端元光谱库中有m个端元，记实数集为令

为L×m的向量空间，A的每一列是一个端元光谱，令

表示光谱库中端元在像元中对应的丰度值组成的列向量，

为m的向量空间。则混合像元y可以用光谱库中端元线性表示：

y=Ax+n   （4）

将混合像元分解问题通过稀疏回归模型表示，(5)式是其中一种等价的稀疏回归模型[2]

$\underset{x}{\min }\frac{1}{2}{{\left|\right|\mathrm{Ax}-y\left|\right|}_2}^2+\mathrm{\λ}{\left|\right|x\left|\right|}_1$ 且x≥0  (5)

其中||Ax-y||₂ ²项反映回归的Ax对观测值y的拟合程度，称之为图像的重构误差；||x||₁项反映解混丰度的稀疏性程度；λ是调整目标函数中两项的权值。由于光谱的可变性，通常丰度的“和为一”条件并不能被满足，模型中只加入了丰度的非负性条件约束。模型（5）可以通过基于变量分离和增光拉格朗日的稀疏解混方法（Sparse Unmixing via Variable Splitting andAugmented Lagrangian，SUnSAL）[2]、迭代光谱分析方法（Iterative Spectral Mixture Analysis，ISMA）[20]、正交匹配追踪方法（Orthogonal Matching Pursuit，OMP）[21]等及其变形方法求解。文献[2]比较了这些稀疏解混方法，认为从整体上来说SUnSAL方法及其变形方法的解混效果最好。

3.结合图像空间信息的稀疏解混模型

经典的稀疏解混模型和方法都是基于图像的光谱空间分析，然而地物在空间分布具有连续性，所以端元丰度的空间分布具有相应的平滑性。考虑图像的空间特征，并将提取的空间信息加入到稀疏解混模型和方法中，可以提高解混的精度，使得解混结果更符合图像真实的空间分布特征。文献[18]通过在稀疏回归的解混模型（5）中加入空间全变分（Total Variation，TV）正则项的方式加入端元丰度的空间平滑性约束，将（5）式变为：

$\underset{x}{\min }\frac{1}{2}{{\left|\right|\mathrm{AX}-Y\left|\right|}_F}^2+\mathrm{\λ}{\left|\right|X\left|\right|}_{\mathrm{1,1}}+{\mathrm{\λ}}_{\mathrm{TV}}\mathrm{TV}\left(X\right)$ 且x_s≥0   (6)

其中

表示整幅高光谱遥感图像，N是高光谱遥感图像中像元的个数，Y的每一列表示一个像元的光谱；

是端元丰度的分布图像，X的每一列表示光谱库中的端元在一个像元中对应的丰度向量；是X的Frobenius范数；

x_s是第s个像元的丰度向量；TV(X)是TV正则项，

表示像元s和其邻接像元t的丰度差，ε表示图像中水平和垂直相邻的像元组，即只考虑4邻域内像元之间的空间关系；λ_TV表示在目标函数中空间项TV(X)对应的权重。令H_h和H_v为定义在丰度分布图像X的垂直、水平的线性算子，其定义如下：

H_h(x_i,j)＝(x_i,j-x_i,j+1)   （7）

H_v(x_i,j)＝(x_i,j-x_i+1,j)   （8）其中x_i,j表示在图像空间中位于(i,j)处像元对应的丰度向量，在图像空间位于(i,j+1)的像元是位于(i,j)的像元水平方向上右边的最邻近像元，位于(i+1,j)的像元是位于(i,j)的像元垂直方向上下边的最邻近像元；在丰度分布图像的边界处采用周期边界的定义，即若丰度分布图像是由k行l列像元组成，H_h(x_i,l)＝(x_i,l-x_i,1)，H_v(x_k,j)＝(x_k,j-x_1,j)。令H表示丰度分布图像X的4邻接差商线性算子，其定义如下：

$H\left(X\right)=\left[\begin{array}{c}{H}_h\left(X\right)\\ H_v\left(X\right)\end{array}\right]---\left(9\right)$

显然TV(X)≡||H(X)||_1,1，（6）式可以等价地表示为

$\underset{x}{\min }\frac{1}{2}{{\left|\right|\mathrm{AX}-Y\left|\right|}_F}^2+\mathrm{\λ}{\left|\right|X\left|\right|}_{\mathrm{1,1}}+{\mathrm{\λ}}_{\mathrm{TV}}{\left|\right|H\left(X\right)\left|\right|}_{\mathrm{1,1}}+{\mathrm{\ι}}_{R+}\left(X\right)---\left(10\right)$

其中函数

ι_R+(x_s)是指标函数，当向量x_s的所有元素都大于或等于0时，函数值为0，否则为+∞。文献[18]基于交替方向乘子法（Alternating Direction Method ofMultipliers，ADMM）求解策略，提出一种基于变量分离增广拉格朗日和全变分的稀疏解混方法（sparse unmixing via variable splitting and augmented Lagrangian and TV，SUnSAL-TV）方法求解（10）式。

与本发明相关的现有技术有如下参考文献：

[]BIOUCAS-DIAS J M,PLAZA A,DOBIGEOON N,et al.Hyperspectral UnmixingOverview:Geometrical,Statistical,and Sparse Regression-Based Approaches[J].IEEE Journal ofSelected Topics in Applied Earth Observations and Remote Sensing，2012，5(2):354-379.

[2]IORDACHE M D,BIOUCAS-DIAS J M,PLAZA A.Sparse Unmixing of HyperspectralData[J].IEEE Transactions on Geoscience and Remote Sensing，2011，49(6):2014-2039.

[3]IORDACHE M D,BIOUCAS-DIAS J M,PLAZA A.Recent Developments in SparseHyperspectral Unmixing[C].2010IEEE International Geoscience and Remote Sensing Symposium,2010:1281-1284.

[4]IORDACHE M D.A Sparse Regression Approach to HyperspectralUnmixing[D].Technical University of Lisbon,2011.

[5IORDACHE M D,BIOUCAS-DIAS J M,PLAZA A.On the Use of Spectral Libraries toPerform Sparse Unmixing of Hyperspectral data[C].IEEE GRSS Workshop Hyperspectral ImageSignal Process.:Evolution in Remote Sensing(WHISPERS),2010,1:1–4.

[6]Yang ZY,et al.Blind Spectral Unmixing Based on Sparse Nonnegative MatrixFactorization[J].IEEE Transactions on Image Processing.2011.20(4):p.1112-1125.

[7]BRUCKSTEIN A M,ELAD,ZIBULEVSKY M.On the Uniqueness of NonnegativeSparse Solutions to Underdetermined Systems of Equations[J].IEEE Transactions on InformationTheory,2008.54(11):4813-4820.

[8]PLAZA A,MARTINEZ P,et al.,Spatial/spectral Endmember Extraction byMultidimensional Morphological Operations[J].IEEE Transactions on Geoscience and RemoteSensing，2002.40(9):p.2025-2041.

[9]PLAZA A,MARTINEZ P,et al.A Quantitative and Comparative Analysis of EndmemberExtraction Algorithms from Hyperspectral Data[J].IEEE Transactions on Geoscience and RemoteSensing，2004.42(3):p.650-663.

[0]ROGGE D M,RIVAED B,et al.Integration of Spatial–spectral Information for theImproved Extraction of Endmembers[J].Remote Sensing of Environment,2007(110):p.287-303.

[11]MEI S H,et al.Spatial Purity Based Endmember Extraction for Spectral Mixture Analysis.IEEE Transactions on Geoscience and Remote Sensing，2010.48(9):3434-3445.

[2]ZORTEA M,PLAZA A.Spatial Preprocessing for endmember extraction[J].IEEETransactions on Geoscience and Remote Sensing，2009，47(8):2679--2693.

[13]MARTIN G,PLAZA A.Region-Based Spatial Preprocessing for Endmember Extractionand Spectral Unmixing[J].IEEE Geoscience and Remote Sensing Letters，2011，8(4):745-749.

[14]MARTIN G,PLAZA A.Spatial-Spectral Preprocessing Prior to Endmember Identificationand Unmixing of Remotely Sensed Hyperspectral Data[J].IEEE Journal of Selected Topics inApplied Earth Observations and Remote Sensing，2012，5(2):380-395.

[15]JIA S,QIAN Y T.Spectral and Spatial Complexity-based Hyperspectral Unmixing[J].IEEE Transaction on Geoscience and Remote Sensing,2007,45(12):3867–3879.

[16]JIA S,QIAN Y T.Constrained Nonnegative Matrix Factorization for HyperspectralUnmixing[J].IEEE Transaction on Geoscience and Remote Sensing,2009.47(1):p.161-173.

[7]ECHES O,DOBIGEON N,TOURNERET J Y.Enhancing Hyperspectral Image UnmixingWith Spatial Correlations[J].IEEE Transaction on Geoscience and Remote Sensing,2011.49(11SIPart1):4239-4247.

[18]IORDACHE M D,BIOUCAS-DIAS J M,PLAZA A.Total Variation SpatialRegularization for Sparse Hyperspectral Unmixing[J].IEEE Transactions on Geoscience andRemote Sensing,2012,50(11):4484-4502.

[19]KESHAVA N,MUSTARD J F.Spectral Unmixing[J].IEEE Signal Processing Magazine,2002:44-57.

[20]ROGGE D M,et al.Iterative Spectral Unmixing for Optimizing Per-pixel EndmemberSets[J].IEEE Transactions on Geoscience and Remote Sensing,2006.44(12):3725-3736.

[21]PATI Y C,REZAHFAR R,KRISHNAPRASAD P.Orthogonal matching pursuit:Recursive function approximation with applications to wavelet decomposition[C].Asilomar Conf.Signals,Systems and Computing(ASSC),2003,1:1–10.

发明内容

本发明的目的在于针对现有高光谱遥感图像结合空间信息的稀疏解混模型的缺点和不足，提供一种基于同质指数的高光谱遥感图像稀疏解混方法，通过图像的空间同质分析估计图像的空间平滑性，按照提取的同质指数来调整空间全变分正则项的权重，解混结果能够较好保持图像真实的空间平滑性。

本发明所采用的技术方案是一种基于同质指数的高光谱遥感图像稀疏解混方法，包括以下步骤：

步骤a，对高光谱遥感图像的每个像元进行同质分析提取同质指数，提取方式如下，

对于高光谱遥感图像的某像元P_i,j，设选择以像元P_i,j为中心、以n为半径的不包含像元P_i,j的方形邻域窗口为Ω(P_i,j)，计算像元P_i,j与它的邻域窗口内像元P_r,c之间的光谱相似度S(P_r,c,P_i,j)，按下式将加权平均的光谱相似度作为该像元的同质指数，

$\mathrm{hi}\left(P_{i,j}\right)=\underset{P_{r,c\∈\mathrm{\Ω}\left(P_{i,j}\right)}}{\mathrm{\Σ}}\mathrm{\ω}\left(P_{r,c}\right)\×S(P_{r,c},P_{i,j});$

其中，ω(P_r,c)表示(r,c)处的空间权重参数；

步骤b，根据整幅高光谱遥图像中像元同质指数的值域对每个像元的同质指数进行线性拉伸，线性拉伸方式如下，

对于高光谱遥感图像的某像元P_i,j，按下式拉伸后的同质指数记为

其中，参数

M为预设参数；

步骤c，定义邻接差商算子如下，

其中，

像元P_i,j+1是像元P_i,j右边的邻接像元，像元P_i+1,j是像元P_i,j下边的邻接像元，x_i,j、x_i,j+1和x_i+1,j分别是像元P_i,j、像元P_i,j+1、像元P_i+1,j对应的端元丰度；且若丰度分布图像是由k行l列像元组成，

步骤d，将步骤c中定义的邻接差商算子作为空间平滑项引入到稀疏解混的目标优化函数中，然后用基于交替方向乘子法策略迭代求解，得到高光谱遥感图像中每个端元的丰度。

而且，步骤a中，按下式计算像元P_i,j与它的邻域窗口内像元P_r,c之间的光谱相似度S(P_r,c,P_i,j)，

$S(P_{r,c},P_{i,j})=\frac{<y_{r,c},y_{i,j}>}{\left|\right|y_{r,c}\left|\right|\left|\right|y_{i,j}\left|\right|}$

其中，y_r,c和y_i,j表示像元P_i,j和P_r,c的光谱。

而且，步骤a中，空间权重参数ω(P_r,c)定义如下，

$\mathrm{\&omega;}\left(P_{r,c}\right)=\frac{e^{-[{(r-i)}^2+{(c-j)}^2]}}{\underset{P_{p,q}\&Element;\mathrm{\&Omega;}\left(P_{i,j}\right)}{\mathrm{\&Sigma;}}{e}^{-[{(p-i)}^2+{(q-j)}^2]}}$

其中，P_p,q是Ω(P_i,j)中的像元，位于图像中(p,q)处。

而且，设端元光谱库中有m个端元，记实数集为令

表示整幅高光谱遥感图像，是端元丰度的分布图像，N是高光谱遥感图像中像元的个数，L是高光谱遥感图像的光谱波段数；

步骤d中，将步骤c中定义的邻接差商算子作为空间平滑项引入到稀疏解混的目标优化函数中，得到的目标优化函数如下，

其中，λ是调整目标函数中两项的权值，λ_TV表示在目标函数中空间项对应的权重；函数

ι_R+(x_s)是指标函数，x_s是第s个像元的丰度向量，当向量x_s的所有元素都大于或等于0时，函数值为0，否则为+∞。

本发明提供的技术方案的有益效果为：对高光谱遥感图像进行空间同质分析，通过提取的同质指数来度量端元丰度的空间平滑性，像元的同质指数是它与一个邻域窗口内像元之间光谱相似度的加权平均值，在一个邻域范围进行加权平均可以降低噪音的干扰，并且权重按照距离以指数衰减可以减小远离中心像元的不相似像元对同质指数的干扰；将线性拉伸后的同质指数作为权重值重新定义邻接差商算子，在图像平滑性较强的地方赋予更强的丰度空间平滑性约束，使得解混丰度的空间分布平滑性更接近图像真实的空间分布平滑性；将邻接差商算子的周期边界定义改为0边界定义减小了在边界处解混的误差；改进后的方法提高了稀疏解混的精度，并且更有效的抑制了噪音对解混结果的影响。

附图说明

图1为本发明实施例的流程图。

图2为本发明实施例的模拟数据的端元空间分布图。

具体实施方式

为了更好地理解本发明的技术方案，下面结合附图和实施例对本发明做进一步的详细说明。

本发明的实施例是基于地物光谱库对高光谱遥感图像进行稀疏解混，主要通过高光谱图像的同质分析，计算同质指数来调整基于全变分的变量分离增广拉格朗日稀疏解混时空间正则项的权重，进而实现基于端元光谱库的高光谱遥感图像的稀疏解混。具体实施时，本发明可采用计算机软件技术实现自动运行流程。

参照图1，本发明实施例的步骤如下：

步骤a：对高光谱遥感图像的每个像元进行同质分析提取同质指数。实施例对于每个像元，选择以像元为中心，以n为半径的不包含像元的方形邻域窗口，计算每个像元与它的邻域窗口内像元之间的光谱相似度，将加权平均的光谱相似度作为该像元的同质指数。

TV正则项表示了丰度分布的空间平滑性，然而这种平滑性在图像空间中并不一致。有些像元与其邻近像元的光谱很相似，叫做同质像元。这些同质像元空间位置相近、光谱也相似，那么它们的地物组成成分也很可能一样，对应的丰度也应该相近。即在同质像元组成的区域，丰度分布保持着很强的空间平滑性。另外一些像元与邻近像元相比，光谱变化较大，表明它们的地物组成和丰度都有较大的变化，这些区域丰度分布的空间平滑性较低。图像的空间同质分析通过比较像元与其邻域内像元之间的光谱相似性得到像元的空间同质性，能较好地度量图像的空间平滑性。

定义高光谱遥感图像的某像元P_i,j的同质指数hi(P_i,j)∈[0,1]是像元与它的邻域窗口内其它像元光谱相似性的加权平均值，按照（11）式计算：

$\mathrm{hi}\left(P_{i,j}\right)=\underset{P_{r,c\&Element;\mathrm{\&Omega;}\left(P_{i,j}\right)}}{\mathrm{\&Sigma;}}\mathrm{\&omega;}\left(P_{r,c}\right)\&times;S(P_{r,c},P_{i,j})---\left(11\right)$

其中P_i,j和P_r,c是分别位于图像中(i,j)和(r,c)处的像元，ω(P_r,c)表示(r,c)处的空间权重参数。S(P_r,c,P_i,j)表示像元P_i,j和P_r,c的光谱相似性，Ω(P_i,j)是以像元P_i,j为中心的不包含P_i,j的方形邻域窗口。具体实施时，本领域技术人员可根据具体情况预先设置窗口尺寸。如果窗口开的过小，包含的像元太少，噪音会对估计的结果产生严重的干扰；如果窗口开的过大，一些距离中心像元较远的像元与中心像元光谱相似性很小，会错误地低估像元的同质指数；在本发明的实施例中设窗口半径n＝2。地理学第一定律（Tobler's First Law of Geography）指出空间距离越近，像元的相关性越高，所以设计空间权重按照离窗口中心的距离衰减。这样距离中心较远的像元对应的权重较小，所以不论是否被包含在窗口内，对同质指数的影响较小。因此降低了同质指数对窗口大小的依赖。函数值按照距离衰减的函数有很多，例如以多项式衰减、以指数衰减。在距离较大的时候，指数比多项式衰减的更快。所以选择指数衰减的空间邻域权重可以减小距离较远且与中心像元光谱相似性较低的像元对同质指数错误地估计。定义空间权重函数ω如下：

$\mathrm{\&omega;}\left(P_{r,c}\right)=\frac{e^{-[{(r-i)}^2+{(c-j)}^2]}}{\underset{P_{p,q}\&Element;\mathrm{\&Omega;}\left(P_{i,j}\right)}{\mathrm{\&Sigma;}}{e}^{-[{(p-i)}^2+{(q-j)}^2]}}---\left(12\right)$

其中，P_p,q是Ω(P_i,j)中的像元，位于图像中(p,q)处。

光谱角距离是高光谱图像分析中最常用的光谱相似性测度，它对光谱的缩放变化具有不变性。本文选择光谱角距离的余弦值即光谱的正则化内积作为两个像元间的光谱相似性测度S(P_r,c,P_i,j)：

$S(P_{r,c},P_{i,j})=\frac{<y_{r,c},y_{i,j}>}{\left|\right|y_{r,c}\left|\right|\left|\right|y_{i,j}\left|\right|}---\left(13\right)$

y_r,c和y_i,j表示像元P_i,j和P_r,c的光谱。

步骤b：对同质指数进行线性拉伸。实施例在步骤a计算了每个像元的同质指数，本步骤根据整幅高光谱遥感图像中像元同质指数的值域对每个像元同质指数进行线性拉伸。

令参数

由于同质指数hi(P_i,j)∈[m,1]比较小，对其做一个线性拉伸。像元P_i,j拉伸后的同质指数记为定义

则

当hi(P_i,j)＝1时，表示像元与其邻近像元之间的光谱完全一样，这时它们的地物丰度值也应该完全一样，应该赋予非常强的TV正则项权重，线性拉伸后是M倍于同质指数最低的像元对应的TV正则项权重，具体实施时，参数M可由本领域技术人员根据具体情况预先设置。

步骤c：定义新的邻接差商算子。

实施例将步骤b中经过线性拉伸后的同质指数作为权重乘子引入到计算像元丰度4邻接差商算子中，定义新的邻接差商算子：

像元P_i,j+1是像元P_i,j右边的邻接像元，像元P_i+1,j是像元P_i,j下边的邻接像元，x_i,j、x_i,j+1和x_i+1,j分别是像元P_i,j、像元P_i,j+1、像元P_i+1,j对应的端元丰度。

在丰度分布图像的边界处差商算子

采用0边界定义，即若丰度分布图像X是由k行l列像元组成，

步骤d：基于同质分析的稀疏解混建模和求解。实施例将步骤c中定义的邻接差商算子作为空间平滑项引入到稀疏解混的目标优化函数中，然后用基于交替方向乘子法策略迭代求解。

设端元光谱库中有m个端元，记实数集为令

表示整幅高光谱遥感图像，

是端元丰度的分布图像，N是高光谱遥感图像中像元的个数，L是高光谱遥感图像的光谱波段数；将（10）式中的丰度差商算子用步骤c中改进后的差商算子代替后的目标优化函数为：

其中，λ是调整目标函数中两项的权值，λ_TV表示在目标函数中空间项对应的权重；函数x_s是第s个像元的丰度向量，ι_R+(x_s)是指标函数，当向量x_s的所有元素都大于或等于0时，函数值为0，否则为+∞。

对（18）采用现有技术中的基于交替方向乘子法策略迭代地求解，得到图像中每个端元的丰度。

综上所述，本发明提出的基于同质指数的高光谱遥感图像稀疏解混方法针对高光谱遥感图像空间连续的不一致性特点，引入了高光谱遥感图像同质分析，通过提取的同质指数度量图像的空间连续性，根据将线性拉伸后的同质指数调整目标函数中空间正则项的权重。因此解混丰度更接近真实丰度的空间分布平滑性。

以下分别通过模拟和真实的高光谱数据实验来验证本发明实施例所提供技术方案的有效性：

1.评价指标

本发明使用均方根误差(Root Mean Square Error,RMSE）和信号重构误差（Signal toReconstruction Error,SRE）两个指标衡量稀疏解混的精度。x表示真实的端元丰度，表示解混得到的端元丰度，RMSE和SRE的定义如下：

RMSE的值越小，表示丰度的估计值越接近真实值，解混的精度越高；SRE是信号的能量与误差的能量的比值，能更好的度量解混精度，与RMSE相反，SRE的值越大，解混精度越高。

2.模拟高光谱数据解混实验

因为真实的高光谱图像像元的组分和丰度信息难以获取，无法做定量分析，而仿真数据中像元的组分和丰度都是已知的。所以本发明设计了一份模拟数据，定量地分析稀疏解混方法的精度。从USGS光谱库splib06中提取498种地物光谱（每种地物光谱有224个波段）建立光谱库为了降低光谱库中端元的相干性，从A中提取子库A′，使得A′中每个端元之间的光谱角距离大于4.44°，A′∈R^224×240。仿真数据由端元光谱库A′中随机选择的7个端元依据（1）式线性混合而成。为了验证同质分析在高光谱图像稀疏解混中的有效性，构造一份50×50个空间像元、224个波段的模拟数据。如图2所示，该数据包含4个同质区，同质区之间是5个过渡区。同质区中像元含有的地物种类和对应的丰度完全相同；其中同质区1和4分别只含有一个端元1、2，是纯净像元区；同质区2和3是同质的混合像元区，同质区2含有2个端元3、4，它们对应的丰度分别为0.3和0.7，同质区3含有3个端元5、6、7，它们对应的丰度分别为0.2，0.3和0.5。过渡区是不同的同质区地物的平滑过渡，所以设计过渡区包含有所有它相邻接的同质区域内的所有端元。令像元中端元的丰度

x_T为端元对应同质区中端元的丰度，d是像元到此同质区的距离。图像中间的过渡区是像元混合最复杂的区域，称之为复杂混合区，与所有同质区相邻，像元的混合程度最高，包含所有的7个端元。以上构造的模拟数据覆盖了从纯净度最高的区域（1个端元）到混合度最高的区域（7个端元）；既有地物丰度稳定的区域（同质区），也有地物丰度变化的区域（过渡区）；既有纯净的同质区也有不同混合程度的同质区。并且由于TV正则项表示像元与水平和垂直邻接像元的丰度差值，在构造的模拟数据中，过渡区域1和2中像元只有水平变化，过渡区域3和4中像元只有垂直变化，过渡区域5中像元既有水平变化也有垂直变化。所以模拟数据很好地包括了各种复杂的地物分布类型，而且非常适合检验SUnSAL-TV方法以及本发明方法的解混效果。在高光谱图像中噪音大部分是低通的模型误差，因此通过对0均值的独立同分布高斯噪音进行低通滤波作为模拟数据的加性噪音。在本文的实验中对模拟数据分别加上20dB、30dB和40dB的加性噪音，来检验噪音对解混算法的影响。

将本发明的方法与SUnSAL方法和SUnSAL-TV方法的解混精度在RMSE和SRE这两个评价指标下进行比较，在不同的噪音水平下，所有算法中都令权重λ＝0.05，λ_TV＝5×10^-4，M＝5。

解混结果的精度如表1所示：

表1不同噪音水平下不同解混方法的精度比较

表1的实验结果表明本发明提出的基于同质指数的高光谱遥感图像稀疏解混方法在不同的图像噪音水平下比SUnSAL方法和SUnSAL-TV方法的解混精度都要高。特别地，当噪音的分贝变大时，以上方法的解混精度都会降低，但是本发明提出方法的解混精度变化相对而言要小，即本发明提出方法的抗噪性更好。

进一步将同质区和过渡区的解混精度进行比较，表2显示了噪音干扰最大（SNR=20dB）时的精度评价结果。

表2同质区域和过渡区域不同解混方法的精度比较

表2的结果显示出，不同的稀疏解混方法下同质区的解混精度都要明显高于过渡区，复杂混合区的解混精度最差（当SNR=30dB或者40dB时仍然得到相同的结论），这表明在解混时应该对同质区和过渡区区别对待。而本发明专利提出的方法首先对图像进行同质分析，同质区的同质指数较高，而过渡区的同质指数较低。据此调整空间正则项权重，对同质区赋予较大的权重，而过渡区赋予较小的权重，所以解混的精度比SUnSAL算法和SUnSAL-TV算法更高。

3.Cuprite地区的AVIRIS数据实验

本实验使用的真实数据为公开的AVIRIS赤铜矿遥感影像数据，有224个波段，被广泛的用于验证端元提取和混合像元分解算法的有效性。实验区的矿物信息已经被深入研究过，所有的端元都在USGS的光谱库splib06中，并且有USGS公开的的矿物分类信息。本实验选取了一个50×50个像元的实验区。剔除低信噪比低的波段，剩下188个波段，即Y∈R^188×2500。选择A′作为稀疏解混的端元光谱库，对应选择188个波段。本实验对真实数据同时使用SUnSAL、SUnSAL-TV和SUnSAL-HTV方法进行解混，方法中λ和λ_TV选择经验值，λ＝λ_TV＝10^-3。因为AVIRIS赤铜矿数据信噪比很高，设置较小的同质参数M，令M=2。由于真实数据没有地物的组分和丰度信息，无法像模拟数据那样进行定量地精度评价。但是仍然可以从以下2个方面验证本发明方法的有效性：

1)与该试验区矿物分类图的对比

虽然本实验采用的Cuprite赤铜矿数据1997年采集，而矿物分类图1995年生产，不能直接比较，但是分类图仍然可以作为一种定性的评价方法。可以将本实验区域在不同解混方法下的丰度空间分布图与分类图视觉上（定性）地比较。值得注意的是分类图是一个二值分布图（属于或者不属于），而丰度图是丰度在[0,1]区间的连续分布图，所以丰度图比分类图分布要广。虽然如此，我们仍然可以从视觉上得出结论：稀疏解混的结果与实际相符。

2)丰度分布的连续性

端元丰度在空间分布具有平滑性，从所得丰度空间分布图可以明显地看出：完全依赖于光谱信息稀疏解混的SUnSAL方法的丰度空间分布平滑性很差，有很多散乱的孤立点；SUnSAL-TV方法加入了空间TV正则项约束，丰度的空间分布平滑性有明显提升，但是仍有少许孤立点；SUnSAL-HTV方法得到的丰度空间分布平滑性最高。

综上所述，本发明提出的基于同质指数的高光谱遥感图像稀疏解混方法在经典的稀疏解混方法SUnSAL和结合空间信息的稀疏解混方法SUnSAL-TV的基础上进一步提高了解混的精度，并且有效地抑制了噪音干扰，是一种可行的结合空间信息和基于稀疏回归的高光谱遥感图像解混方法。

Claims (4)

1.一种基于同质指数的高光谱遥感图像稀疏解混方法，其特征在于，包括以下步骤：

步骤a，对高光谱遥感图像的每个像元进行同质分析提取同质指数，提取方式如下，

对于高光谱遥感图像的某像元P_i,j，设选择以像元P_i,j为中心、以n为半径的不包含像元P_i,j的方形邻域窗口为Ω(P_i,j)，计算像元P_i,j与它的邻域窗口内像元P_r,c之间的光谱相似度S(P_r,c,P_i,j)，按下式将加权平均的光谱相似度作为该像元的同质指数，

$\mathrm{hi}\left(P_{i,j}\right)=\underset{P_{r,c\&Element;\mathrm{\&Omega;}\left(P_{i,j}\right)}}{\mathrm{\&Sigma;}}\mathrm{\&omega;}\left(P_{r,c}\right)\&times;S(P_{r,c},P_{i,j});$

其中，ω(P_r,c)表示(r,c)处的空间权重参数；

步骤b，根据整幅高光谱遥图像中像元同质指数的值域对每个像元的同质指数进行线性拉伸，线性拉伸方式如下，

对于高光谱遥感图像的某像元P_i,j，按下式拉伸后的同质指数记为

其中，参数

M为预设参数；

步骤c，定义邻接差商算子如下，

其中，

像元P_i,j+1是像元P_i,j右边的邻接像元，像元P_i+1,j是像元P_i,j下边的邻接像元，x_i,j、x_i,j+1和x_i+1,j分别是像元P_i,j、像元P_i,j+1、像元P_i+1,j对应的端元丰度；且若丰度分布图像是由k行l列像元组成，

步骤d，将步骤c中定义的邻接差商算子作为空间平滑项引入到稀疏解混的目标优化函数中，然后用基于交替方向乘子法策略迭代求解，得到高光谱遥感图像中每个端元的丰度。

2.根据权利要求1所述的高光谱遥感图像稀疏解混方法，其特征在于：步骤a中，按下式计算像元P_i,j与它的邻域窗口内像元P_r,c之间的光谱相似度S(P_r,c,P_i,j)，

$S(P_{r,c},P_{i,j})=\frac{<y_{r,c},y_{i,j}>}{\left|\right|y_{r,c}\left|\right|\left|\right|y_{i,j}\left|\right|}$

其中，y_r,c和y_i,j表示像元P_i,j和P_r,c的光谱。

3.根据权利要求1所述的高光谱遥感图像稀疏解混方法，其特征在于：步骤a中，空间权重参数ω(P_r,c)定义如下，

$\mathrm{\&omega;}\left(P_{r,c}\right)=\frac{e^{-[{(r-i)}^2+{(c-j)}^2]}}{\underset{P_{p,q}\&Element;\mathrm{\&Omega;}\left(P_{i,j}\right)}{\mathrm{\&Sigma;}}{e}^{-[{(p-i)}^2+{(q-j)}^2]}}$

其中，P_p,q是Ω(P_i,j)中的像元，位于图像中(p,q)处。

4.根据权利要求1或2或3所述的高光谱遥感图像稀疏解混方法，其特征在于：设端元光谱库中有m个端元，记实数集为

，令

表示整幅高光谱遥感图像，

是端元丰度的分布图像，N是高光谱遥感图像中像元的个数，L是高光谱遥感图像的光谱波段数；

步骤d中，将步骤c中定义的邻接差商算子作为空间平滑项引入到稀疏解混的目标优化函数中，得到的目标优化函数如下，

其中，λ是调整目标函数中两项的权值，λ_TV表示在目标函数中空间项对应的权重；函数

x_s是第s个像元的丰度向量，ι_R+(x_s)是指标函数，当向量x_s的所有元素都大于或等于0时，函数值为0，否则为+∞。

Images