CN101739563A - 一种基于信源可变的约束单象素ica的混合像元分解方法 - Google Patents

Abstract

本发明为一种信源可变、基于约束单象素ICA的混合像元分解方法。ICA是一种常用的盲源分离方法，但由于混合像元分解问题存在其他盲源分离问题所不具有的约束条件，ICA用于混合像元分解问题并不可以取得好的分解效果，同时，混合像元分解领域里存在两大难题：如何解决地物类型大于遥感图像通道数的混合像元分解问题，如何解决线性模型中因地物类型始终固定所带来的精度问题。为有效利用ICA模型解决混合像元分解问题和有效解决混合像元分解的上述两个问题，本发明对传统ICA进行改进，提出了一种新的混合像元分解的方法-信源可变的约束单象素ICA方法。试验证明本发明的方法可取得优于其他分解方法的结果。

Description

一种基于信源可变的约束单象素ICA的混合像元分解方法

技术领域：

本发明属于遥感图像混合像元分解技术领域，具体涉及一种信源可变的、基于约束单象素ICA的混合像元分解方法。

背景技术：

传感器空间分辨率的限制，和地表地物复杂多样的影响，使得混合像元普遍存在于各类遥感图像中，对于地面地物分布比较复杂的区域尤其如此(参考对比文件1)。如果不考虑混合像元的影响而将其当作纯像元来考虑，势必会带来分类误差，导致精度下降，不能反映真实的地物覆盖情况。

在现有混合像元模型中线性模型是迄今为止最受欢迎、使用最多的一种模型。但线形模型中，混合像元分解的地物类型的种类始终固定，而实际中大多数的影像区域只是由全部地物类型中的部分类型组成。因此，在混合像元分解过程中，待分解像元的地物类型如果可以随着象素点的变化而变化，会相应的提高混合像元的分解精度。关于混合像元分解过程中分解类型随象素点的变化而变化的思想已引起了较大的关注(参考对比文件2)。同时，遥感图像混合像元分解中的绝大部分理论方法都要求遥感图像的通道数目大于地物类型的种类(参考对比文件3)，而某些遥感图像(如合成孔径雷达图像)并没有那么多的通道数目(参考对比文件4)，于是，地物种类大于遥感图像通道数的混合像元分解问题也已成为混合像元分解问题的一个瓶颈问题。

独立成分分析(Independent Component Analysis，ICA)是近年发展起来的一种多道信号处理方法，在核磁共振成像、光谱遥感图像分类、图像特征提取、遥感图像变化检测、数字水印、信号检测等领域也已得到较好的应用。但现有的ICA的目标函数不能满足图像混合像元分解中所要求的一些约束条件(即分解结果之和为1，分解结果中每一类地物的比例应在[0，1]范围内，等等)，ICA还未在图像混合像元分解问题中取得较好应用。

本发明在混合像元分解的线形模型及ICA相关理论的基础上，提出一种新的混合像元分解方法-信源可变的约束单象素ICA分解方法。本发明中的方法动态调整地物种类，在提高了分类精度降低分类误差的同时也实现了地物种类多于通道数目的混合像元的分解。

对比文件1：C.-I Chang，H.Ren，C.-C.Chang，J.O.Jensen and F.DˇAmico，Estimation of subpixel target size forremotely sensed imagery[J]，IEEE Trans.on Geoscience and Remote Sensing，2004，42(6)：1309-1320

对比文件2：Plaza A，Martinex P，Perez R.A quantitative and comparative analysis of endmember extracionalgorithms from hyperspectral data[J].IEEE Transaction on Geoscience and Remote Sensing，2004，42(3)：650-663

对比文件3：ZHANG Yan，SHAO Mei-Zhen.Decomposing mixed pixels by means of RBF neural network[J].Journal of Remote Sensing，2004，7；285-289.

对比文件4：Thomas M.Lillesand，Ralph W.Kiefer.Remote Sensing And Image Interpretation[M]，New York：JohnWiley & Sons，Inc，2003

发明内容：

本发明的目的在于解决混合像元分解领域里的两大难题：地物类型大于遥感图像通道数的混合像元分解问题，线性模型中因地物类型始终固定所带来的精度问题。为解决上述两个问题，本发明对传统ICA模型进行改进，提出了一种新的混合像元分解的方法-信源可变的约束单象素ICA方法。动态调整待分解像元的地物类型，利用混合像元与纯净像元的相似性，判断尚待分解的混合像元中包含的地物类型，结合本发明中提出的新的ICA--带约束条件的单象素ICA进行混合像元分解。

其中动态调整待分解像元的地物类型的特征在于：采用交叉相关光谱匹配技术，计算像元光谱和参考光谱之间的响应值，以此来判断两光谱之间的相似程度，选择与参考光谱相似程度最高的类型作为当前待分解像元的地物类型。像元光谱和参考光谱之间的响应值函数如下式：

$r=\frac{n(\mathrm{\Σ}{R}_r{R}_t-\mathrm{\Σ}{R}_r\mathrm{\Σ}{R}_t)}{\sqrt{[\mathrm{n\Σ}{R_r}^2-{\left(\mathrm{\Σ}{R}_r\right)}^2\stackrel{-}{]}[\mathrm{n\Σ}{R_t}^2-{\left(\mathrm{\Σ}{R}_r\right)}^2]}}$

其中带约束条件的单象素ICA的特征在于：把负熵和误差函数统一在一个框架内，形成了带有约束条件的组合目标函数，同时在随机梯度方法的基础上改进得到双梯度下降方法。组合目标函数的具体构造如下式：

$\underset{\mathrm{\λ},A}{\max }L(S,A)=J\left(S\right)-{\mathrm{\λ}}^T(\mathrm{AS}-X)-({\mathrm{\λ}}_0-1)$

$(\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{s}_i-1)-\underset{i=1}{\overset{m}{\mathrm{\Σ}}}[s_i^{2h}+{(1-s_i)}^{2h}]$

式中，J(S)是负熵，

是误差函数，误差函数保证了混合像元的分解结果满足同一象素点的不同地物类型的分解结果之和为1，同一象素点的不同地物类型的分解结果均大于0小于1，整个式子保证了目标函数满足混合像元分解的约束条件

其中动态调整分解像元的地物类型和利用带约束条件的单象素ICA进行分解的联系在于：在利用带约束条件的单象素ICA进行分解的过程中，待分解像元的地物类型一直在动态调整中

其中带约束条件的单象素ICA中，约束条件的特征在于：利用误差函数约束单象素ICA的分解结果，以保证同一象素点的不同地物类型的分解结果之和为1，同一象素点的不同地物类型的分解结果均大于0小于1.误差函数的具体构造如下式：

E(S，A)＝λ^T(AS-X)+(λ₀-1)

$(\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{s}_i-1)-\underset{i=1}{\overset{m}{\mathrm{\Σ}}}[s_i^{2h}+{(1-s_i)}^{2h}]$

式中第一项是负熵，第二项是误差函数，误差函数保证了混合像元的分解结果满足同一象素点的不同地物类型的分解结果之和为1，同一象素点的不同地物类型的分解结果均均为大于0小于1的实数，从而保证了目标函数满足混合像元分解的约束条件。

附图说明：

图1是从北京市ASAR遥感图像(VV，VH)中截取的512*512大小的图像。

图2是用PCA对图1进行混合像元分解的结果，从左向右依次为：人工建筑物/植被、道路/空地。

图3是用改进的ICA(本发明提出方法的子方法)对图1进行混合像元分解的结果，从左向右依次为：人工建筑物/植被、道路/空地。

图4是用本发明方法对图1进行混合像元分解的结果，从左向右依次为：人工建筑物、道路/空地、植被。

图5是从北京市IK遥感图像中截取的北京师范大学地区的图像，图中包含此遥感图像的1通道，2通道，4通道。

图6是用遥感处理软件ERDAS中的分解方法对图5进行混合像元分解的结果，从左向右依次为：道路/空地、人工建筑物、植被、水体。

图7是用本发明方法对图5进行混合像元分解的结果，从左向右依次为：道路/空地、人工建筑物、植被、水体。

图8是图6，图7分解精度的平均误差比较图。

图9是本发明的整体流程图。

具体实施方式：

本发明提出的信源可变的约束单象素ICA方法，具体步骤如下：为动态调整待分解像元的地物类型(即信源)，首先利用混合像元与纯净像元的相似性，判断出尚待分解的混合像元中包含的地物类型，而后结合本发明中提出的新的ICA方法：带约束条件的单象素ICA方法进行分解。

上述方法中，所述“动态调整待分解像元的地物类型”如下：

传统的混合像元分解对每个像元采用统一的标准，即针对每个像元的信源数都是不变的。但实际上不同像元的信源可能不同。根据混合像元的光谱响应特性，采用交叉相关光谱匹配技术，计算像元光谱和参考光谱之间的响应值，以此来判断两光谱之间的相似程度，从而保证当前待分解像元所选择的地物类型是与参考光谱相似程度最高的类型。而响应值小于某一较小值的辅助像元可忽略不计，不列入此像元分解的地物类型范围内。

上文中，“响应值”应按照如下公式计算求得：

$r=\frac{\mathrm{\Σ}(R_r-\stackrel{\‾}{R_r})(R_t-\stackrel{\‾}{R_t})}{\sqrt{\left[\mathrm{\Σ}{(R_r-\stackrel{\‾}{R_r})}^2\right]\left[\mathrm{\Σ}{(R_t-\stackrel{\‾}{R_t})}^2\right]}}$

式中，R_r，R_t分别为参考光谱和像元光谱，r是要计算求解的响应值。由协方差的性质可知，上式等同于：

$r=\frac{n(\mathrm{\Σ}{R}_r{R}_t-\mathrm{\Σ}{R}_r\mathrm{\Σ}{R}_t)}{\sqrt{[\mathrm{n\Σ}{R_r}^2-{\left(\mathrm{\Σ}{R}_r\right)}^2\stackrel{-}{]}[\mathrm{n\Σ}{R_t}^2-{\left(\mathrm{\Σ}{R}_r\right)}^2]}}$

式中，n为极化方式数目。

上述方法中，所述“带约束条件的单象素ICA方法”如下：

在中心化和白化之后，本发明中的方法以负熵为基础设计目标函数，得到带约束条件的单象素ICA的目标函数如下：

$\underset{\mathrm{\λ},A}{\max }L(S,A)=J\left(S\right)-{\mathrm{\λ}}^T(\mathrm{AS}-X)-({\mathrm{\λ}}_0-1)(\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{s}_i-1)-\underset{i=1}{\overset{m}{\mathrm{\Σ}}}[s_i^{2h}+{(1-s_i)}^{2h}]$

式中，X为已知的观测信号量，也是待分解象素所观测的象素值.S是原信号向量，即尚待分解象素中地物类型所占的比例，也是我们分解的目标，S，A均为未知，上式的目标即是通过迭代求解出使得L(S，A)取最大值S，A。J(S)是S的负熵，负熵描述了非高斯信号的熵，其定义为J(S)＝H(S_guass)-H(S)，其中H表示熵，S_guass是具有与S相同方差的高斯随机向量。在X具有零均值和单位方差的前提下，可用一种基于最大熵原理的估计方法来估计负熵：J(S)∝[E{G(S)}-E{G(S_guass)}]²，其中G可取

本申请所述方法采用的就是这个负熵求解公式。

是误差函数，这一误差函数中，λ是拉格朗日乘子向量，λ₀是自由能向量，S是原信号向量，也即是要待分解象素中地物类型所占的比例，也即是我们分解的目标，A是混合矩阵，s_i是S中的各个分量，也即是S＝(s₁，s₂，...，s_n)^T，h是一个足够大的量，它保证s_i ^2h和(1-s_i)^2h是一个近似为0的值.误差函数中

保证了混合像元的分解结果满足同一象素点的不同地物类型的分解结果之和为1的条件，

保证了同一象素点的不同地物类型的分解结果满足分解结果均大于0小于1的条件，λ^T(AS-X)保证了混合矩阵与分解结果的乘积AS与X的吻合程度.整个式子保证了目标函数满足混合像元分解的约束条件.

上述目标函数把负熵和误差函数统一在一个框架内，形成带有约束条件的组合目标函数。确定目标函数后，本发明在随机梯度方法的基础上改进得到双梯度下降方法。采用如下的双梯度下降方式：

式中ΔA是每次迭代时混合矩阵A的变化量，Δλ是每次迭代时....λ的变化量，

是L(S，A)对A求偏导，

是L(S，A)对λ求偏导.a是一个系数，经过试验验证a＝0.1时迭代效果最好。

本方法的具体迭代步骤如下：

1在遥感影像中选取N个不同的纯净像元，以其光谱值作为参考光谱，利用上述方法中“动态调整待分解像元的地物类型”部分所述方法选择带分解像元的地物类型。

2用本发明中提出的新的混合像元分解方法：带约束条件的单象素ICA方法进行分解，分解过程如下：

2.1中心化、白化，等同于传统的ICA算法中的中心化、白化步骤。

2.2按本发明所述方法确立目标函数。

2.3按本发明所述迭代方法进行迭代，直至分类精度足够精确为止.

仿真结果

具体的仿真结果如下：

采用北京市ASAR遥感图像(VV，VH)作试验，从北京地区ENVISAT-ASAR图像(VV、VH)中截取512*512大小的图像作为数据源。从图像和实际调查知图中主要包含3类地物：道路/空地、人工建筑物和植被。

原图及实验结果如图1-图4所示：

图1的两幅图中，类似于标记1处的地方均为人工建筑物/植被，类似于标记2处的地方均为道路/空地，在原始图像中，人工建筑和植被混合在一起。

试验结果如图2-图4。图2的两幅图是PCA的分解结果。图2的两幅图中左侧的是人工建筑物/植被图，右侧的是道路/空地图。图3的两幅图是约束单像素ICA的分解结果。图3的两幅图中左侧的是人工建筑物/植被图，右侧的是道路/空地图。通过图2和图3的标记1的比较可见：图2中结果的人工建筑物/植被的轮廓已变得模糊不好辨认，而图3中的人工建筑物/植被的轮廓就比较清晰。比较图2与图3的标记2可见：图2中结果中，建筑物中间的细微的道路已完全模糊变得不可辨认，而图3中的建筑物中间的细微的道路还比较清晰。但图2图3所示的分解方法均只可以分解出两种地物类型。图4所示的本申请中的方法可以分解出3类地物类型，图4中最左侧的一幅是人工建筑物图，中间一幅是道路/空地图，最右边一幅是植被图。通过图2图3和图4的标记1的比较可见：图4中的人工建筑物的轮廓明显要比图2图3的清晰。比较图2图3和图4的标记2可见：图4中细微的道路要比图2图3中的细微道路均清晰。而且如图4中标记3所示，本申请的方法可以分解出多于原图数据的分解类型，在这一试验实例中，即为：植被图。上述结果说明：PCA、约束单像素ICA方法只能解决地物类型小于遥感图像通道数的混合像元分解问题，将其建筑物和植被分为一类，而本申请提出的方法不仅有着优于PCA、ICA的目视效果，还能解决地物类型大于遥感图像通道数的混合像元分解问题。

采用均方根误差RMSE和系统误差SE两个指标来衡量分解结果：

${\mathrm{RMSE}}_i=\sqrt{\frac{\underset{i=1,j=1}{\overset{i=m,j=N}{\mathrm{\Σ}}}{(\widehat{s_i\left(j\right)}-s_i\left(j\right))}^2}{N}},$ ${\mathrm{SE}}_i=\frac{\underset{i=1,j=1}{\overset{i=m,j=N}{\mathrm{\Σ}}}|\widehat{s_i}-s_i\left(j\right)|}{N}$

RMSE＝∑RMSE_i/∑i，SE＝∑SE_i/∑i

式子中，s_i(j)代表真实的地物比例，是本申请中算法分解中的地物比例，N是图像中像素的总数目。

表1

本发明的均方根误差和系统误差均最小，上述数据定量的说明了本文发明方法的优越性。

进一步采用北京市IK遥感图像作试验，从中截取北京师范大学地区的图像作为数据源。从实际调查知图中主要包含4类地物：道路/空地、人工建筑物、植被、水体。

原图及实验结果如图5-图7所示：

图5是从北京市IK遥感图像中截取的北京师范大学地区的图像，图中包含此遥感图像的1通道，2通道，4通道。图5的类似于标记1处的地方是道路/空地较多的地方，类似于标记2处的地方是人工建筑，类似于标记3处的地方是植被/水体，在原始图像图5中，植被和水体混合在一起。

试验结果如图6-图7。图6的四幅图是用遥感处理软件ERDAS中的分解方法对图5进行混合像元分解后的结果。图6的四幅图从左向右依次为：道路/空地、人工建筑物、植被、水体。图7的四幅图是本申请提出的方法的分解结果。图7的四幅图中从左向右依次是道路/空地、人工建筑物、植被、水体。通过图6和图7的标记1的比较可见：图7中结果的细微轮廓保存的完好，而图6结果的细微轮廓已经变得不好辨认。比较图6与图7中的标记1的、标记2可见：图7中的水体比图6中的水体清晰明朗。这一试验实例也说明着本申请中所提算法的较好的分解目视效果。

为比较分解结果，从本发明方法的处理结果中分别随机选取25、50、75、100、125、150、175、200个象素点的分解结果。同时在ENVIEW中读出此像素对应的地理坐标，再通过实地考查得到这一地理坐标下的真实地物。采用如下方法对分解精度作定量衡量。计算真实地物的成分含量比例，设为p₁，p₂，p₃，真实地物成分含量的比例和必为1。分解精度为：

$k=1-\frac{1}{3}*(\frac{|p_1-s_1|}{p_1}+\frac{|p_2-s_2|}{p_2}+\frac{|p_3-s_3|}{p_3})$

式中用

计算平均误差。

附图8所示的数据说明本申请提出的方法可以分解出比通道数目更多地物种类，同时具有较高的分解精度与和较好的目视结果。

本发明是一种信源可变的、基于约束单象素ICA的混合像元分解方法。按照本发明的方法，可以在分解出的地物类型大于遥感图像通道数的前提下，达到很好的混合像元分解效果。

本发明所述的方法不限于具体实施方式中所述的实施例，本领域技术人员根据本发明的技术方案得出其他的实施例，只要是通过信源可变的、基于约束单象素ICA来实现混合像元的分解，包括实现相应功能的装置，也应当同样属于本发明的创新范围.

Claims (5)

1.一种基于信源可变的约束单象素ICA的混合像元分解方法，其特征在于动态调整待分解像元的地物类型，利用混合像元与纯净像元的相似性，判断尚待分解的混合像元中包含的地物类型，采用带约束条件的单象素ICA进行混合像元分解。

2.如权利要求1所述的基于信源可变的约束单象素ICA的混合像元分解方法，其中动态调整待分解像元的地物类型的特征在于：采用交叉相关光谱匹配技术，计算像元光谱和参考光谱之间的响应值，以此来判断两光谱之间的相似程度，选择与参考光谱相似程度最高的类型作为当前待分解像元的地物类型。像元光谱和参考光谱之间的响应值函数如下式：

$r=\frac{n(\mathrm{\Σ}{R}_r{R}_t-\mathrm{\Σ}{R}_r\mathrm{\Σ}{R}_t)}{\sqrt{[\mathrm{n\Σ}{R_r}^2-{\left(\mathrm{\Σ}{R}_r\right)}^2\stackrel{\‾}{]}[\mathrm{n\Σ}{R_t}^2-{\left(\mathrm{\Σ}{R}_r\right)}^2]}}$

式中，R_r，R_t分别为参考光谱和像元光谱，r即要计算求解的响应值。由协方差的性质可知，上式等同于：

$r=\frac{n(\mathrm{\Σ}{R}_r{R}_t-\mathrm{\Σ}{R}_r\mathrm{\Σ}{R}_t)}{\sqrt{[\mathrm{n\Σ}{R_r}^2-{\left(\mathrm{\Σ}{R}_r\right)}^2\stackrel{\‾}{]}[\mathrm{n\Σ}{R_t}^2-{\left(\mathrm{\Σ}{R}_r\right)}^2]}}$

式中，n为极化方式数目。

3.如权利要求1所述的基于信源可变的约束单象素ICA的混合像元分解方法，其中带约束条件的单象素ICA的特征在于：把负熵和误差函数统一在一个框架内，形成了带有约束条件的组合目标函数，同时在随机梯度方法的基础上改进得到双梯度下降方法。组合目标函数的具体构造如下式：

$\underset{\mathrm{\λ},A}{\max }L(S,A)=J\left(S\right)-{\mathrm{\λ}}^T(\mathrm{AS}-X)-({\mathrm{\λ}}_0-1)(\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{s}_i-1)-\underset{i=1}{\overset{m}{\mathrm{\Σ}}}[s_i^{2h}+{(1-s_i)}^{2h}]$

式中，X是观测信号量，即待分解象素所观测的象素值，是已知量.S是原信号向量，也即待分解象素中地物类型所占的比例，是分解的目标，S，A均为未知，上式的目标即是通过迭代求解出使得L(S，A)取最大值的S，A.J(S)是S的负熵，负熵描述了非高斯信号的熵，其定义为J(S)＝H(S_guass)-H(S)，其中H表示熵，S_guass是具有与S相同方差的高斯随机向量。在X具有零均值和单位方差的前提下，可用一种基于最大熵原理的估计方法来估计负熵：J(S)∝[E{G(S)}-E{G(S_guass)}]²，其中G可取

本申请所述方法采用的就是这个负熵求解公式。

4.如权利要求1所述的基于信源可变的约束单象素ICA的混合像元分解方法，其中动态调整分解像元的地物类型和利用带约束条件的单象素ICA进行分解的联系在于：在利用带约束条件的单象素ICA进行分解的过程中，待分解像元的地物类型一直在动态调整中。

5.如权利要求1所述的基于信源可变的约束单象素ICA的混合像元分解方法，其中带约束条件的单象素ICA中，约束条件的特征在于：利用误差函数约束单象素ICA的分解结果，以保证同一象素点的不同地物类型的分解结果之和为1，同一象素点的不同地物类型的分解结果均大于0小于1.误差函数的具体构造如下式：

E(S，A)＝λ^T(AS-X)+(λ₀-1)

$(\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{s}_i-1)+\underset{i=1}{\overset{m}{\mathrm{\Σ}}}[s_i^{2h}+{(1-s_i)}^{2h}]$

上述误差函数保证了混合像元的分解结果满足同一象素点的不同地物类型的分解结果之和为1，同一象素点的不同地物类型的分解结果均为大于0小于1的实数，从而保证了目标函数满足混合像元分解的约束条件。

这一误差函数中，λ是拉格朗日乘子向量，λ₀是自由能向量，S是原信号向量，即待分解象素中地物类型所占的比例，即我们分解的目标，A是混合矩阵，s_i是S中的各个分量，即S＝(s₁，s₂，...，s_n)^T，h是一个足够大的量，它保证s_i ^2h和(1-s_i)^2h是一个近似为0的值.误差函数中

保证了混合像元的分解结果满足同一象素点的不同地物类型的分解结果之和为1的条件，

保证了同一象素点的不同地物类型的分解结果满足分解结果均大于0小于1的条件，λ^T(AS-X)保证了混合矩阵与分解结果的乘积AS与X的吻合程度.

整个式子保证了目标函数满足混合像元分解的约束条件.

Images