CN104751181B

CN104751181B - 一种基于相对丰度的高光谱图像解混方法

Info

Publication number: CN104751181B
Application number: CN201510155418.4A
Authority: CN
Inventors: 刘治; 聂明钰; 肖晓燕; 张伟; 唐波; 孙育霖
Original assignee: Shandong University
Current assignee: Shandong University
Priority date: 2015-04-02
Filing date: 2015-04-02
Publication date: 2018-01-12
Anticipated expiration: 2035-04-02
Also published as: CN104751181A

Abstract

本发明公开了一种基于相对丰度的高光谱图像解混方法，该方法包括以下步骤：对高光谱数据进行小波分析去噪，然后通过迭代，优化目标函数，获取端元矩阵及相对端元丰度矩阵。本发明算法简单，计算量小，提高了高光谱混合像元解混的算法实用性。具体事例表明，该方法可以较好的完成高光谱混合像元的解混。

Description

一种基于相对丰度的高光谱图像解混方法

技术领域

本发明属于图像处理技术领域，涉及基于高光谱图像解混方法，尤其涉及一种基于相对丰度的高光谱图像解混方法。

背景技术

高光谱图像是同时反映物体反射光谱信息和空间信息的三维数据图像。利用图像分光仪及传感器在紫外线，可见光，近红外，红外获取相对连续的一系列窄波段下的光谱特征。在三维图像数据中，两维代表空间信息，另一维代表一个像元的光谱信息。高光谱图像由其快速无损等特性，常用于遥感、食品、医药、石油化工等行业。

由于传感器的空间分辨率限制以及自然界地物的复杂多样性，一些像元中常常不止含有一种物质，这样的像元被称为混合像元。由于混合像元的的存在，使得传统的分类方法不适于地物的精确细分和鉴别。为了提高遥感分类的准确性，必须解决混合像元的问题，因此混合像元解混在地物分类和识别之前成为关键的一步。通常情况，假设高光谱图像满足线性混合模型(LMM)，该模型是指图像中的像元是由组成图像的基本物质(端元)在不同波段下的光谱特性线性组合而成。因此混合像元解混，就是结合高光谱自身的信息，利用图像处理方法，获得组成像元的端元，和端元在对应像元中所占的比例，又称丰度。端元的丰度需要满足非负(ANC)和和为1的限制(ASC)。

传统的混合像元的解混一般可以分为两个步骤：端元提取和丰度反演。当像元中含有纯净像元时，即该像元中只含有一种物质，从几何角度出发，常用解混方法有纯像元指数(PPI)，N-FINDR,迭代误差分析(IEA)，点成分分析(VCA)和自动端元提取(AEE)等。当像元中无纯净像元时，采用带有最小体积约束的方法，如最小体积单形体分析(MVSA)，通过分裂增强拉格朗日单形体识别(SISAL)及基于最小闭合单形体体积的凸分析(MVES)。在所有的端元提取出之后，常常利用全限制最小二乘法(FCLS)对相应的端元进行丰度反演。除了上述的传统混合像元解混，盲信号分解算法也常常用于高光谱的解混中。在盲信号分解中，同时分解出图像中的端元及其对应的丰度值。常用的方法有独立成分分析(ICA)，光谱和空间复杂度分析，非负矩阵分解(NMF)及其相关算法等。如果解混之前可以用光谱辐射计从地面获取光谱数据库，那么解混就变成了从这个大的光谱数据库中优化出一个子集，然后求解出子集中端元对应的丰度。常用稀疏回归算法来实现此优化。

上述的各种解混方法都各有利弊，如传统算法中有些端元提取必须要求像元中含有纯净端元；带有最小体积约束的算法过于复杂；盲信号分解算法用于解混时，优化过程计算量很大，耗时长；稀疏回归算法必须要求提前获取光谱库。因此，针对不同的高光谱图像数据，需要合理的选择像元解混方法。

发明内容

本发明的目的就是为了解决上述问题，提供一种基于相对丰度的高光谱图像解混方法，它具有用理论相对简单，耗时短的算法实现混合像元的解混，同时获取像元中的端元及其相对丰度的优点。

为了实现上述目的，本发明采用如下技术方案：

一种基于相对丰度的高光谱图像解混方法，包括如下步骤：

步骤(1)：输入高光谱图像数据矩阵及端元个数K，所述高光谱图像数据矩阵是B×N的矩阵，其中B为高光谱图像包含的总波段数，N为图像中像元总个数；并利用值在0，1间的随机数初始化相对丰度矩阵U，使其满足约束条件；所述相对丰度矩阵U是K×N的矩阵，其中K表示端元个数，N表示图像中像元总个数；

步骤(2)：利用小波变换对高光谱图像数据矩阵在不同波段方向进行降噪处理，得到去噪后的高光谱数据矩阵X；

步骤(3)：用步骤(2)降噪后的高光谱数据矩阵X以及步骤(1)中初始化的相对丰度矩阵U，计算端元矩阵M；

步骤(4)：利用步骤(1)得到的初始化相对丰度矩阵U和步骤(3)中计算的端元矩阵M，计算目标函数J；

步骤(5)：如果步骤(4)所计算的目标函数J小于预先设定阈值或者J的改变量小于设定阈值，就进入步骤(6)；否则，就计算相对丰度矩阵U，并且返回步骤(3)；

步骤(6)：算法收敛，输出高光谱图像解混后的端元矩阵M和相对丰度矩阵U，算法结束。

所述步骤(1)的约束条件为：

式中，K为端元个数；u_ij为第j个像元中第i个端元的相对丰度。

所述步骤(4)计算计算目标函数J的计算公式为：

式中，N为像元总数；K为端元个数；x_j为第j个像元；m_i为M中第i个端元；u_ij为第j个像元中第i个端元的相对丰度；d为像元到端元的欧式距离；l∈[1,∞)是一个加权指数，通过最优化J，使得端元与像元之间的非相似性指数的价值函数达到最小。

所述步骤(3)计算端元矩阵M的计算公式为：

其中u_ij为第j个像元中第i个端元的相对丰度，j取值范围为1～N，i取值范围为1～K；m_i为M中第i个端元；K为端元个数，N为高光谱数据X的像元总数；x_j为高光谱数据X的第j个像元；l∈[1,∞)是一个加权指数。

所述步骤(5)计算相对丰度矩阵U的计算公式为：

其中，d_ij为高光谱数据X的第j个像元到端元矩阵M的第i个端元的欧式距离，i取值范围为1～K，j取值范围为1～N；u_ij为第j个像元中第i个端元的相对丰度；l∈[1,∞)是一个加权指数；d_tj为高光谱数据X的第j个像元到端元矩阵M的第t个端元的欧式距离，t的取值范围为1～K。l一般取值为2。

本发明的有益效果：

1.本发明不需要提前假设高光谱图像数据中含有纯净端元，不需要提前获取所有端元的光谱库，只需要设定端元个数，然后可直接对高光谱数据进行解混。

2.本发明具有耗时短，速度快，算法简单的优点，解决了一般解混方法算法复杂，计算速度慢的缺点。

3.本发明可以同时获取高光谱图像的端元矩阵和相对丰度矩阵，不需要分步进行。

附图说明

图1为本发明的高光谱解混流程图；

图2(a)原始图像及选取作为端元的像元；

图2(b)端元的光谱曲线；

图3(a)-图3(f)使用基于相对丰度解混方法提取的端元(虚线)光谱线与从原始图像中提取的端元(实线)的比较图；

图4(a)-图4(f)使用基于相对丰度解混方法解混结果相对丰度显示(0为白色，1为黑色)。

具体实施方式

下面结合附图与实施例对本发明作进一步说明。

如图2(a)所示，图为后续实验所用数据中的一个波段下的高光谱图像，该图像来源于HYDICE Urban高光谱数据，维数为307×307。图像数据中包含六种物质：路(Road)，泥土(Dirt)，树(Trees)，草(Grass)，屋顶(Roof)，金属(Metal)。

如图2(b)所示，为从后续实验所用高光谱数据中，提取出的各种不同物质在不同波段下的光谱特性，提取位置为图2(a)中序号标注处。

该种解混方法是在已知端元个数的情况下，利用目标函数和迭代算法求出高光谱图像中的端元和相对丰度。

之所以称为相对丰度而不是丰度，是因为该方法的解混结果是不同端元在不同像元中的相对含量，并不是实际含量，即解混结果U并不满足X＝MW，即U不是W，其中X为混合像元矩阵，M为端元矩阵，W为丰度矩阵。该种解混结果中的U只是代表像元中几种端元的相对含量大小，因此称作相对丰度矩阵。

为优化J，构造如下新的目标函数，可使上式达到最小值的必要条件：

这里，j＝1,...,N，是等式的N个约束式的拉格朗日乘子。

对所有输入参量求导，使式(1)达到最小的必要条件为：

从式(4)看出，若像元j与端元i的相似性越大，d越小，则u_ij越大，即像元j含有端元i的含量越高，相对丰度越大，符合高光谱图像解混的物理意义。

由上述两个必要条件，基于相对丰度的高光谱解混方法是一个简单的迭代过程。在批处理方式运行时，该解混方法采用下列步骤确定端元矩阵M和相对丰度矩阵U：

步骤1用值在0，1间的随机数初始化相对丰度矩阵U，使其满足式(1)中的约束条件。

步骤2用式(3)计算K个端元，i＝1,...,K。

步骤3根据式(1)计算目标函数。如果它小于某个设定的阈值，或它相对上次价值函数值的改变量小于某个阈值，则算法停止。

步骤4用式(4)计算新的U阵。返回步骤2。

当算法收敛时，就得到了各个端元和各个端元在像元中对应的相对丰度，从而完成了解混过程。

在进行解混之前，首先利用小波分析去除噪声。小波分析是一种时间-频率的变换，它具有多分辨率分析的特性，可以在时域和频域中提取表征信号局部特征的信息。

如图1所示，本发明的的高光谱解混具体包括以下步骤：

1.利用小波分析，对原始数据进行去噪。

2.构建目标函数

其中，N为像元总数；K为端元个数；x_j为第j为像元；m_i为M中第i个端元；u_ij为第j个像元中第i个端元的相对丰度；d为像元到端元的欧式距离；l∈[1,∞)是一个加权指数，在实验中取值l＝2。

3.通过迭代最小化(1)式J，其中d²(x_j,m_i)计算所需的m_i由(3)式替换，u_ij由(4)式替换。

如果J小于某一个阈值或者两次循环改变量小于某个阈值，则停止迭代。试验中迭代次数常常取值50-500，改变量阈值常常取值10^-3-10^-5。

具体实施例：

实验数据来源于HYDICE Urban高光谱数据，该数据包含210个光谱波段，维数为307×307×210。图像数据中包含六种物质：路，泥土，树，草，屋顶，金属。在实验中，去掉水吸收影响的波段，该数据剩下178个波段。利用小波分析去掉噪声，然后利用基于相对丰度的高光谱图像解混方法对数据进行解混。试验中迭代次数取值500，该变量阈值取值10^-5，端元个数为6。实验结果如图3(a)、图3(b)、图3(c)、图3(d)、图3(e)、图3(f)、图4(a)、图4(b)、图4(c)、图4(d)、图4(e)、图4(f)所示。

图3(a)、图3(b)、图3(c)、图3(d)、图3(e)、图3(f)分别为泥土(Dirt)、金属(Metal)、路(Road)、草(Grass)、树(Trees)、屋顶(Roof)的光谱曲线比较图，虚线为使用基于相对丰度解混方法提取的端元光谱线，实线为从原始图像中提取的端元光谱线。由比较结果可以看出，除了屋顶的解混提取与原图像提取光谱曲线有较大的差异之外，其他物质的解混提取端元光谱曲线与原图像提取的光谱曲线十分接近，说明基于相对丰度解混方法解混结果比较好。

图4(a)、图4(b)、图4(c)、图4(d)、图4(e)、图4(f)分别为使用基于相对丰度解混方法解混得到的泥土(Dirt)、金属(Metal)、路(Road)、草(Grass)、树(Trees)、屋顶(Roof)的相对丰度图。其中，0为白色，1为黑色，颜色越重，说明该区域中该种物质相对含量越多。例如，图4(a)为泥土(Dirt)在该图像中的相对丰度图，颜色越重的区域说明泥土的相对含量越多。

上述虽然结合附图对本发明的具体实施方式进行了描述，但并非对本发明保护范围的限制，所属领域技术人员应该明白，在本发明的技术方案的基础上，本领域技术人员不需要付出创造性劳动即可做出的各种修改或变形仍在本发明的保护范围以内。

Claims

1.一种基于相对丰度的高光谱图像解混方法，其特征是，包括如下步骤：

所述步骤(1)的约束条件为：

<mfenced open = "" close = ""> <mtable> <mtr> <mtd> <mrow> <mi>s</mi> <mo>.</mo> <mi>t</mi> <mo>.</mo> </mrow> </mtd> <mtd> <mrow> <munderover> <mo>&Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>K</mi> </munderover> <mrow> <mo>(</mo> <msub> <mi>u</mi> <mrow> <mi>i</mi> <mi>j</mi> </mrow> </msub> <mo>)</mo> </mrow> <mo>=</mo> <mn>1</mn> </mrow> </mtd> <mtd> <mrow> <mo>&ForAll;</mo> <mi>j</mi> <mo>=</mo> <mn>1</mn> <mo>,</mo> <mo>...</mo> <mo>,</mo> <mi>N</mi> </mrow> </mtd> </mtr> </mtable> </mfenced>

式中，K为端元个数；u_ij为第j个像元中第i个端元的相对丰度；

所述步骤(3)计算端元矩阵M的计算公式为：

<mrow> <msub> <mi>m</mi> <mi>i</mi> </msub> <mo>=</mo> <mfrac> <mrow> <munderover> <mi>&Sigma;</mi> <mrow> <mi>j</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>N</mi> </munderover> <msup> <msub> <mi>u</mi> <mrow> <mi>i</mi> <mi>j</mi> </mrow> </msub> <mi>l</mi> </msup> <msub> <mi>x</mi> <mi>j</mi> </msub> </mrow> <mrow> <munderover> <mi>&Sigma;</mi> <mrow> <mi>j</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>N</mi> </munderover> <msup> <msub> <mi>u</mi> <mrow> <mi>i</mi> <mi>j</mi> </mrow> </msub> <mi>l</mi> </msup> </mrow> </mfrac> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>3</mn> <mo>)</mo> </mrow> </mrow>

其中u_ij为第j个像元中第i个端元的相对丰度，j取值范围为1～N，i取值范围为1～K；m_i为M中第i个端元；K为端元个数，N为高光谱数据X的像元总数；x_j为高光谱数据X的第j个像元；l∈[1,∞)是一个加权指数；

所述步骤(4)计算目标函数J的计算公式为：

<mrow> <mi>min</mi> <mi> </mi> <mi>J</mi> <mrow> <mo>(</mo> <mi>U</mi> <mo>,</mo> <mi>M</mi> <mo>)</mo> </mrow> <mo>=</mo> <munderover> <mo>&Sigma;</mo> <mrow> <mi>j</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>N</mi> </munderover> <munderover> <mo>&Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>K</mi> </munderover> <msup> <mrow> <mo>(</mo> <msub> <mi>u</mi> <mrow> <mi>i</mi> <mi>j</mi> </mrow> </msub> <mo>)</mo> </mrow> <mi>l</mi> </msup> <msup> <mi>d</mi> <mn>2</mn> </msup> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>j</mi> </msub> <mo>,</mo> <msub> <mi>m</mi> <mi>i</mi> </msub> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow>

式中，N为像元总数；K为端元个数；x_j为第j个像元；m_i为M中第i个端元；u_ij为第j个像元中第i个端元的相对丰度；d为像元到端元的欧式距离；l∈[1,∞)是一个加权指数，通过最优化J，使得端元与像元之间的非相似性指数的价值函数达到最小；

所述步骤(5)计算相对丰度矩阵U的计算公式为：

<mrow> <msub> <mi>u</mi> <mrow> <mi>i</mi> <mi>j</mi> </mrow> </msub> <mo>=</mo> <mfrac> <mn>1</mn> <mrow> <munderover> <mo>&Sigma;</mo> <mrow> <mi>t</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>K</mi> </munderover> <msup> <mrow> <mo>(</mo> <mfrac> <msub> <mi>d</mi> <mrow> <mi>i</mi> <mi>j</mi> </mrow> </msub> <msub> <mi>d</mi> <mrow> <mi>t</mi> <mi>j</mi> </mrow> </msub> </mfrac> <mo>)</mo> </mrow> <mrow> <mn>2</mn> <mo>/</mo> <mrow> <mo>(</mo> <mi>l</mi> <mo>-</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow> </msup> </mrow> </mfrac> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>4</mn> <mo>)</mo> </mrow> </mrow>

其中，d_ij为高光谱数据X的第j个像元到端元矩阵M的第i个端元的欧式距离，i取值范围为1～K，j取值范围为1～N；u_ij为第j个像元中第i个端元的相对丰度；l∈[1,∞)是一个加权指数；d_tj为高光谱数据X的第j个像元到端元矩阵M的第t个端元的欧式距离，t的取值范围为1～K；