CN108664899A - 基于模型驱动和rvm回归的高光谱图像混合像元分解方法 - Google Patents

Abstract

本发明公开了一种基于模型驱动和RVM回归的高光谱图像混合像元分解方法，包括以下步骤：S1，读取高光谱图像，利用端元提取算法计算出包含该高光谱图像所有波段的端元个数；S2，提取出高光谱图像中所有端元的光谱向量；S3，根据Hapke模型，计算出每个端元的单次散射反照率，将端元的单次散射反照率进行线性混合，转换成混合像元的反照率，本发明的有益效果是：采用了非线性光谱混合模型，能够对混合像元的非线性特性进行较好的解释。与传统的线性混合像元分解方法相比，本文的方法精度更高，且RVM回归模型具有概率输出，对于混合像元中各组分丰度值的计算结果具有统计意义。

Description

基于模型驱动和RVM回归的高光谱图像混合像元分解方法

技术领域

本发明涉及一种较大范围，具体是一种基于模型驱动和RVM回归的高光谱图像混合像元分解方法。

背景技术

高光谱影像的空间分辨率较低,其像元的光谱特征通常不是单一地物的光谱特征,而是几种地物光谱特征的混合反映。混合像元的存在,是像元级影像分类和面积量测精度难以达到使用要求的主要原因。由于多光谱影像的波段数量有限,无法实现准确的光谱解混。高光谱影像的光谱分辨率从微米级提高到纳米级,使混合像元分析、分解及其模型研究成为可能。通过高光谱影像混合像元分解,可以提高影像中提取小目标的能力、地物精细分类精度和影像处理的自动化、智能化水平。高光谱混合像元分解成为国内外的研究热点。

遥感图像混合像元分解有线性和非线性两种方法。传统的线性光谱混合模型因计算简便而广泛应用。然而由于物质之间存在散射，加上大气传输过程中的混合效应，会导致真实地物的空间分布也远比线性混合模型所假设的大尺度宏观分布要复杂的多，对于真实地物的混合像元分解，如果用线性光谱混合模型来求解势必存在误差。因此需要研究非线性光谱混合模型及分解算法。

关于非线性混合像元分解问题，国内外学者提出了基于人工神经网络、支持向量机(support vector machine,SVM)等方法。人工神经网络方法容易收敛于局部最优解，而且计算结果缺乏明确的物理意义。SVM方法在混合像元分解中取得了很好的效果。然而这种方法的局限在于与数据的特征空间性质直接相关的核函数的选取，且核函数参数的选取对计算结果有很大的影响，核函数必须满足Mercer条件；预测时间长等。

2000年，模式识别领域出现了相关向量机的研究(RVM,Relevance VectorMachine)。目前RVM已应用于模式识别、高光谱图像分类等领域。杨国鹏利用RVM分类器输出的后验概率实现了高光谱图像混合像元分解，然而上述研究是一种基于数据驱动的混合像元分解，并没有考虑到光谱非线性混合的物理特性。

发明内容

本发明的目的在于提供一种基于模型驱动和RVM回归的高光谱图像混合像元分解方法，以解决上述背景技术中提出的问题。

为实现上述目的，本发明提供如下技术方案：

一种基于模型驱动和RVM回归的高光谱图像混合像元分解方法，包括以下步骤：

S1，读取高光谱图像，利用端元提取算法计算出包含该高光谱图像所有波段的端元个数；

S2，提取出高光谱图像中所有端元的光谱向量；

S3，根据Hapke模型，计算出每个端元的单次散射反照率，将端元的单次散射反照率进行线性混合，转换成混合像元的反照率；

S4，将上述步骤计算出的高光谱图像中混合像元的反射率数据分为两部分，分别为训练样本部分和测试样本部分，分别用于相关向量机模型的训练和混合像元分解实验结果的精度验证，构造相关相量机回归模型，利用训练样本对相关向量机模型进行训练，计算出相关向量机模型的参数；

S5，基于步骤S3得到的混合像元的反射率数据，将高光谱图像的像元光谱输入到步骤S4中训练好的相关向量机模型中，得到相关向量机输出的概率预测值；

S6，将相关向量机输出的概率预测值转换成高光谱图像中每个混合像元属于各地物类别的概率值，即可输出混合像元各组分的丰度值，即像元在对应区域内的每一类别地物的地表覆盖比例。

作为本发明进一步的方案：步骤S3中，利用Hapke模型生成混合像元反照率的过程如下：

(1)将提取出的端元光谱的反射率转换为单次散射反照率

式中：μ＝cos i，μ₀＝sin i，其中i为光照入射角，g为相位角，ω为均值单次散射反照率，r为二向反射率，B为后向散射的量度，P为粒子的位相函数，H为多向散射函数，

(2)对于各向同性的散射粒子以及低反照比的表面，对上式进行如下简化：

其中，

(3)基于Hapke模型，利用端元的单次散射反照率计算混合像元的反射率R，根据Hapke模型，混合像元的平均单次散射反照率是端元的单次散射反照率的线性混合，假设影像中共有n个端元，在[0,1]之间产生n个随机数f_i，并且满足f_j≥0且则混合的反照率为：

其中n为端元数，ω为均值单次散射反照率；

(4)利用上式把端元混合的反照率逆向转换为混合的反射率R，将反射率数据附加一定信噪比的正态分布的随机噪声，此时的数据可以作为模拟的非线性混合的数据。

作为本发明再进一步的方案：步骤S4中，相关向量机模型的训练过程如下：

(1)确定训练样本集，利用稀疏贝叶斯模型，建立相关向量机预测模型，

给定训练样本集为其中x_n∈R^d为训练样本向量，对于回归问题，t_n为目标向量，相关向量机的预测式如下式：

式中，t＝(t₁,t₂,...,t_n)^T为目标向量，w＝(w₀,w₁,...,w_n)^T为为权重系数，Φ＝[φ(x₁),φ(x₂),...,φ(x_M)]为样本基函数映射组成的列向量，φ_i(x)(i＝1,…,N)是定义在训练样本点上的核函数，即φ_i(x)＝k(x,x_i)；

(2)根据贝叶斯准则，求解相关向量机模型的权值向量w的后验概率，如下式：

其中，均值为μ＝σ²ΣΦ^Tt，协方差为Σ＝(σ^-2Φ^TΦ+A)^-1，A＝diag(a0,a1,…,aN)；

(3)根据下式确定相关向量机模型超参数的似然分布：

p(t|a,σ²)＝p(t|w,σ²)p(w|α)dw＝N(t|0,C)；

其中，协方差C＝σ²I+σA^-1φ^T；

(4)计算相关向量机模型的超参数α和σ²，超参数α和σ²的计算可转化为超参数后验分布p(a,σ²|t)∝p(t|a,σ²)p(a)p(σ²)关于a和σ²的最大值问题，在一致超先验分布的情况下，最大化p(a,σ²|t)的问题可以通过下式求解：

其中，C＝σ²I+ΦA^-1Φ^T；

(5)相关向量机模型的权值的求解：对于一组新的输入数据x_*，其相应输出t_*的预测值符合高斯分布：

其中，预测均值与协方差分别为：

y_*＝μ^Tφ(x_*)

与现有技术相比，本发明的有益效果是：采用了非线性光谱混合模型，能够对混合像元的非线性特性进行较好的解释。与传统的线性混合像元分解方法相比，本文的方法精度更高，且RVM回归模型具有概率输出，对于混合像元中各组分丰度值的计算结果具有统计意义。与非线性混合像元分解的SVR方法相比，利用RVM模型进行混合像元分解，其优势在于，本发明的方法具有明确的物理意义，能够有效提高解混精度及可靠性，且RVM回归模型的核函数不需要满足Mercer条件，模型需要求解的参数较少，解的稀疏性更高，训练过程更快。

附图说明

图1示出了本发明所述的流程示意图。

图2示出了本发明所述的一优选实施例的美国内华达州南部Cuprite地区的AVIRIS高光谱图像。

图3示出了本发明所述的一优选实施例的美国内华达州南部Cuprite地区的AVIRIS高光谱图像提取的端元位置分布图。

图4示出了本发明所述的一优选实施例的美国内华达州南部Cuprite地区的AVIRIS高光谱图像使用RVM回归方法的混合像元分解得到的各组分的丰度图(由左至右依次为钙铁榴石、高岭石和胶岭石)。

图5示出了本发明所述的一优选实施例的美国内华达州南部Cuprite地区的AVIRIS高光谱图像使用SVM回归方法的混合像元分解得到的各组分的丰度图(由左至右依次为钙铁榴石、高岭石和胶岭石)。

图6示出了本发明所述的一优选实施例的美国内华达州南部Cuprite地区的AVIRIS高光谱图像使用最小二乘法线性混合像元分解得到的各组分的丰度图(由左至右依次为钙铁榴石、高岭石和胶岭石)。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

请参阅图1，本发明实施例中，一种基于模型驱动和RVM回归的高光谱图像混合像元分解方法，包括以下步骤：

S1，读取高光谱图像，利用端元提取算法计算出包含该高光谱图像所有波段的端元个数；

S2，提取出高光谱图像中所有端元的光谱向量；

S3，根据Hapke模型，将每个像元的光谱值转换成每个端元的单次散射反照率，将端元的单次散射反照率进行线性混合，转换成混合像元的反照率，其步骤如下：

(1)将提取出的端元光谱的反射率转换为单次散射反照率

式中：μ＝cos i，μ₀＝sin i，其中i为光照入射角，g为相位角，ω为均值单次散射反照率，r为二向反射率，B为后向散射的量度，P为粒子的位相函数，H为多向散射函数，

(2)对于各向同性的散射粒子以及低反照比的表面，对上式进行如下简化：

其中，

(3)基于Hapke模型，利用端元的单次散射反照率计算混合像元的反射率R。根据Hapke模型，混合像元的平均单次散射反照率是端元的单次散射反照率的线性混合，假设影像中共有n个端元，在[0,1]之间产生n个随机数f_i，并且满足f_j≥0且则混合的反照率为：

其中n为端元数，ω为均值单次散射反照率；

(4)将反射率数据附加一定信噪比的正态分布的随机噪声，此时的数据可以作为模拟的非线性混合的数据。

S4，将上述步骤计算出的高光谱图像中混合像元的反射率数据分为两部分，分别为训练样本部分和测试样本部分，分别用于相关向量机模型的训练和混合像元分解实验结果的精度验证，构造相关相量机回归模型，利用训练样本对相关向量机模型进行训练，计算出相关向量机模型的参数；相关向量机模型的训练过程如下：

(1)确定训练样本集，利用稀疏贝叶斯模型，建立相关向量机预测模型，

给定训练样本集为其中x_n∈R^d为训练样本向量，对于回归问题，t_n为目标向量，采用步骤S3所述的Hapke非线性光谱混合模型构造出训练样本X_n,目标向量是随机产生的0-1之间的百分比构成的像元的丰度矩阵，相关向量机的预测式如下式：

式中，t＝(t₁,t₂,...,t_n)^T为目标向量，w＝(w₀,w₁,...,w_n)^T为为权重系数，Φ＝[φ(x₁),φ(x₂),...,φ(x_M)]为样本基函数映射组成的列向量，φ_i(x)(i＝1,…,N)是定义在训练样本点上的核函数，即φ_i(x)＝k(x,x_i)；

本发明采用Hapke非线性混合光谱模型构造的样本作为模型的输入，在相关向量机模型的参数迭代求解过程中采用非线性光谱混合模型，因此模型的训练结果具有明确的物理意义。

(2)根据贝叶斯准则，求解相关向量机模型的权值向量w的后验概率，如下式：

其中，均值为μ＝σ²ΣΦ^Tt，协方差为Σ＝(σ^-2Φ^TΦ+A)^-1，A＝diag(a0,a1,…,aN)；

(3)根据下式求解相关向量机模型超参数的似然分布：

p(t|a,σ²)＝p(t|w,σ²)p(w|α)dw＝N(t|0,C)；

其中，协方差C＝σ²I+σA^-1φ^T；

(4)计算相关向量机模型的超参数α和σ²，超参数α和σ²的计算可转化为超参数后验分布p(a,σ²|t)∝p(t|a,σ²)p(a)p(σ²)关于a和σ²的最大值问题，在一致超先验分布的情况下，最大化p(a,σ²|t)的问题可以表示为下式：

其中，C＝σ²I+ΦA^-1Φ^T；

(5)相关向量机模型的权值的求解：对于一组新的输入数据x_*，其相应输出t_*的预测值符合高斯分布：

其中，预测均值与协方差分别为：

y_*＝μ^Tφ(x_*)

S5，基于步骤S3，将高光谱图像混合像元的光谱向量转换成反射率数据，将高光谱图像的像元反射率数据输入到步骤S4中训练好的相关向量机模型中，得到相关向量机输出的概率预测值；

S6，将相关向量机输出的概率预测值转换成高光谱图像中每个混合像元属于各地物类别的概率值，即可输出混合像元各组分的丰度值，即像元在对应区域内的每一类别地物的地表覆盖比例。

为了说明本发明的有效性，特进行如下实验论证。实验采用的计算机硬件环境为Intel Core2双核CPU，1.58GHz/3.25GB内存，软件环境为Microsoft Windows XP、MatlabR2008a。本发明所述算法由MATLAB R2008a实现。

实验目标：

(1)从解的精度、核函数的稀疏程度、混合像元分解计算时间等方面的性能比较基于RVM和SVR的两种非线性混合像元分解方法；

(2)基于真实高光谱遥感图像，将RVM回归非线性分解方法、SVR非线性分解方法和最小二乘(least squares，LS)线性分解三种算法的解混精度进行比较分析。

图2示出了本发明优选实施例的高光谱图像。该图像是ENVI软件自带的美国内华达州南部Cuprite地区的AVIRIS高光谱图像，该图像波长为2.071～2.56μm，光谱分辨率10nm，共50个波段。该地区富含各种矿物质，常被作为各种遥感图像算法测试的标准数据。本发明从原始图像中截取大小为50×50的子图作为实验数据，如图2所示。

结合附图及实施方式详细说明如下：

1.模拟数据实验

1.1主要实施步骤

(1)实验采用HFC方法对数据虚拟维数进行估计，求出端元个数为3。

(2)本发明利用正交子空间投影方法对影像进行端元提取，得到影像端元分布情况如图3所示。参照http://speclab.cr.usgs.gov所提供的矿物实地分布数据，实验区图像的3种主要矿物为：(a)钙铁榴石，(b)高岭石，(c)胶岭石。

(3)通过遥感软件ENVI取得端元光谱。利用本发明所述步骤S3生成2500个模拟数据。研究表明，入射角取值为30度时实验效果较好。本发明的实验中，SNR取值为30：1，入射角为30度。

(4)将模拟数据按数量随机等分成两部分分别作为训练样本和测试样本。采用训练样本计算RVM回归模型的参数，对相关向量机和SVR进行训练。训练过程中SVR的松弛因子ε设为0.01。Gauss核函数参数取不同值时，采用测试样本进行模拟实验的精度评定，基函数数量为各分类器的基函数个数之和。

1.2模拟数据实验精度评价

采用测试样本进行模拟实验的精度评定，实验结果如(表1RVM回归方法解混性能统计表)所示：

表1 RVM回归方法解混性能统计表

常用的核函数包括线性核，多项式核与径向基核函数。已有研究表明，径向基核函数具有良好的非线性处理能力，为便于比较，实验中RVM与SVR两种算法都采用径向基核函数，其表达式如下式：

k(x,x’)＝exp(-γ||x-x'||)

其中，γ为核函数参数，它控制了核函数宽度。

通过SVR方法进行混合像元分解时，惩罚系数C需要通过交叉验证获取。本发明的实验中惩罚系数C取2.0，表2是核函数参数γ取不同值时的SVR方法解混的性能参数。

表2 SVR方法解混性能统计表

通过RVM与SVR两种方法的比较，可以得出以下结论：

(a)在RVM模型中不存在惩罚系数C，只需要确定核函数参数，节省了计算复杂度。

(b)RVM解的稀疏程度更高，在模型训练过程中，相关向量的数目远少于支持向量。RVM与SVR两种方法的训练时间相当，但是由于两种方法都是通过核函数的线性组合进行预测的，因此当测试样本较多时，RVM实时性好，适合处理大场景高光谱影像。

(c)在解混精度方面，仅当核函数参数γ等于1时，SVR方法解混的精度略高于RVM回归方法。当核函数参数γ小于等于0.8时，RVM回归方法解混的精度均高于SVR方法。总体来讲，RVM回归模型能够用较少的相关向量获得相对较高的解混精度。

2真实图像数据实验

2.1主要实施步骤

用模拟数据完成RVM与SVR的训练后，再以真实高光谱图像的反射光谱作为输入，利用训练好的RVM与SVR即可回归出影像中像元各组分的丰度值。

RVM模型中核函数取值。从表1、表2可以看出，RVM模型的训练结果中，Gauss核函数参数γ取0.8时，分解效果最好。因此真实高光谱遥感图像分解实验中，RVM的Gauss核函数参数γ取0.8。为了与RVM进行对比，SVR的Gauss核函数参数与RVM一致。

(3)为了进行对比研究，分别用RVM的非线性混合像元分解、基于SVR的非线性混合像元分解和最小二乘线性分解三种方法，对真实高光谱遥感图像进行混合像元分解。图4是实验区图像基于RVM回归的高光谱图像混合像元分解方法。基于SVR方法的的非线性混合像元分解结果如图5，图5是最小二乘法线性分解的结果。

图4、5、6分别给出了采用3种算法(RVM回归、SVR和最小二乘法)基于真实高光谱影像解混得到的端元丰度图。丰度图中灰度越亮代表该地物在混合像元中所占的比例越大，反之代表该地物所占的比例越小。从分解效果来看，基于RVM回归与SVR两种非线性分解方法的分解效果明显优于最小二乘法线性分解的结果。图4和图5中不同矿物的区分明显，而最小二乘法线性分解的丰度图中存在大量的噪声点。

2.2真实高光谱图像数据实验结果

本发明采用模拟数据建立的测试数据集进行了混合像元分解的精度评价。以各组分的RMSE和总平均RMSE作为评价指标。其中SVR的径向基函数参数取值与RVM相同，均为0.8。表3(三种解混方法模拟数据分解误差比较)给出了三种方法产生的各类别的RMSE与总平均RMSE。

表3三种解混方法模拟数据分解误差比较

从表3可以看出，RVM回归方法的分解精度高于SVR方法和最小二乘法。RVM的非线性分解、SVR非线性分解两种分解方法各类别的RMSE均不超出3％。RVM回归方法中三种端元的分解误差与总平均RMSE均低于SVR方法。而最小二乘法分解的误差则较大。这是由于训练数据与测试数据是根据Hapke非线性模型模拟的，因此利用线性分解方法去解决非线性光谱混合问题误差较大。

对于本领域技术人员而言，显然本发明不限于上述示范性实施例的细节，而且在不背离本发明的精神或基本特征的情况下，能够以其他的具体形式实现本发明。因此，无论从哪一点来看，均应将实施例看作是示范性的，而且是非限制性的，本发明的范围由所附权利要求而不是上述说明限定，因此旨在将落在权利要求的等同要件的含义和范围内的所有变化囊括在本发明内。不应将权利要求中的任何附图标记视为限制所涉及的权利要求。

此外，应当理解，虽然本说明书按照实施方式加以描述，但并非每个实施方式仅包含一个独立的技术方案，说明书的这种叙述方式仅仅是为清楚起见，本领域技术人员应当将说明书作为一个整体，各实施例中的技术方案也可以经适当组合，形成本领域技术人员可以理解的其他实施方式。

Claims (3)

1.一种基于模型驱动和RVM回归的高光谱图像混合像元分解方法，其特征在于，包括以下步骤；

S1，读取高光谱图像，利用端元提取算法计算出包含该高光谱图像所有波段的端元个数；

S2，提取出高光谱图像中所有端元的光谱向量；

S3，根据Hapke模型，计算出每个端元的单次散射反照率，将端元的单次散射反照率进行线性混合，转换成混合像元的反照率；

S4，将上述步骤计算出的高光谱图像中混合像元的反射率数据分为两部分，分别为训练样本部分和测试样本部分，分别用于相关向量机模型的训练和混合像元分解实验结果的精度验证，构造相关相量机回归模型，利用训练样本对相关向量机模型进行训练，计算出相关向量机模型的参数；

S5，基于步骤S3得到的混合像元的反射率数据，将高光谱图像的像元光谱输入到步骤S4中训练好的相关向量机模型中，得到相关向量机输出的概率预测值；

S6，将相关向量机输出的概率预测值转换成高光谱图像中每个混合像元属于各地物类别的概率值，即可得到混合像元各组分的丰度值，即像元在对应区域内的每一类别地物的地表覆盖比例。

2.根据权利要求1所述的一种基于模型驱动和RVM回归的高光谱图像混合像元分解方法，其特征在于，步骤S3中，利用Hapke模型生成混合像元反照率的过程如下：

将提取出的端元光谱的反射率转换为单次散射反照率

式中：μ＝cosi，μ₀＝sini，其中i为光照入射角，g为相位角，ω为均值单次散射反照率，r为二向反射率，B为后向散射的量度，P为粒子的位相函数，H为多向散射函数，

(2)对于各向同性的散射粒子以及低反照比的表面，对上式进行如下简化：

其中，

(3)基于Hapke模型，利用端元的单次散射反照率计算混合像元的反射率R，根据Hapke模型，混合像元的平均单次散射反照率是端元的单次散射反照率的线性混合，假设影像中共有n个端元，在[0,1]之间产生n个随机数fi，并且满足f_j≥0且则混合的反照率为：

其中n为端元数，ω为均值单次散射反照率；

(4)利用上式把端元混合的反照率逆向转换为混合的反射率R，将反射率数据附加一定信噪比的正态分布的随机噪声，此时的数据可以作为模拟的非线性混合的数据。

3.根据权利要求1所述的一种基于模型驱动和RVM回归的高光谱图像混合像元分解方法，其特征在于，

步骤S4中，相关向量机模型的训练过程如下：

(1)确定训练样本集，利用稀疏贝叶斯模型，建立相关向量机预测模型，

给定训练样本集为其中x_n∈R^d为训练样本向量，对于回归问题，t_n为目标向量，相关向量机回归模型的预测式如下式：

式中，t＝(t₁,t₂,...,t_n)^T为目标向量，w＝(w₀,w₁,...,w_n)^T为为权重系数，Φ＝[φ(x₁),φ(x₂),...,φ(x_M)]为样本基函数映射组成的列向量，φ_i(x)(i＝1,…,N)是定义在训练样本点上的核函数，即φ_i(x)＝k(x,x_i)；

(2)根据贝叶斯准则，求解相关向量机模型的权值向量w的后验概率，如下式：

其中，均值为μ＝σ²ΣΦ^Tt，协方差为Σ＝(σ^-2Φ^TΦ+A)^-1，A＝diag(a0,a1,…,aN)；

(3)根据下式求解相关向量机模型超参数的似然分布：

p(t|a,σ²)＝p(t|w,σ²)p(w|α)dw＝N(t|0,C)；

其中，协方差C＝σ²I+σA^-1φ^T；

(4)计算相关向量机模型的超参数α和σ²，超参数α和σ²的计算可转化为超参数后验分布p(a,σ²|t)∝p(t|a,σ²)p(a)p(σ²)关于a和σ²的最大值问题，在一致超先验分布的情况下，最大化p(a,σ²|t)的问题可以表示为下式：

其中，C＝σ²I+ΦA^-1Φ^T；

(5)相关向量机模型的权值的求解：对于一组新的输入数据x_*，其相应输出t_*的预测值符合高斯分布：

其中，预测均值与协方差分别为：

y_*＝μ^Tφ(x_*)