CN111105363A - 一种含噪高光谱图像快速解混方法 - Google Patents

Abstract

本发明公开了一种含噪高光谱图像快速解混方法，包括以下步骤：采用VCA算法提取含噪声的端元光谱；对提取的端元光谱进行SSA算法去噪；利用去噪后的端元光谱根据最小二乘法进行丰度估计；本发明方法直接对原始的高光谱图像的端元进行提取，并对端元光谱进行去噪，然后利用去噪后的端元进行丰度估计，实现高光谱图像的解混，该方法避免了现有方法在第一步中所产生的误差，同时也解决了原始图像去噪过程中重要信息丢失而降低解混精度的问题，此外，本方法对端元光谱去噪，数据量小，解决了图像去噪后解混技术效率低的问题，本方法不仅可以提高解混精度，操作简单而且处理速度快。

Description

一种含噪高光谱图像快速解混方法

技术领域

本发明涉及高光谱图像解混技术领域，尤其涉及一种含噪高光谱图像快速解混方法。

背景技术

混合像元在高光谱图像上广泛存在，它影响了基于高光谱图像的地物识别精度，图像中噪声的存在对解混结果的稀疏性及稳定性有较大影响，因此，含噪高光谱图像的快速解混在应急事件中的地物精细分类发挥着重要作用，同时也对地面目标检测有着重要的意义；

噪声问题的解决国内外的研究学者提出了各种各样的有效方法，至今，去噪方法主要分为以下三类：第一是基于一维频谱域的去噪方法，该类方法将光谱影像看作是含有噪声的光谱曲线集合，对每条光谱曲线进行噪声去除，算法实现相对简单；其次是基于二维空间域的去噪方法，这类方法仅利用高光谱的空间信息，也就是高光谱的单波段影像看作不同的灰度影像，但是按影像波段去除噪声使得频谱域的信息被忽略；第三是基于三维空间域和波谱域相结合的去噪方法，该方法充分利用了高光谱遥感影像的特征，相对与前面两类算法来说，这类算法实现难度较大；

到目前为止，解混技术已经有几十年的发展历史，国内外解混的主要方法可以分为4类：几何分析方法、统计分析方法、稀疏回归分析方法以及光谱-空间联合分析方法；以往含噪高光谱图像解混是先对整幅图像去噪然后再对图像解混，这虽然能得到较好的图像去噪效果，但同时具有以下几个缺点：1)在图像去噪过程中，可能会破坏地物的空间相关性或丢失图像中部分有用信息，对解混造成不利的影响；2)高光谱图像数据量大，对图像去噪所需的时间较长，效率低；3)去噪和解混分步进行，容易造成误差累积，影响最终的解混精度。因此，在实际应用中，以往的方法会受到图像去噪速度的制约，无法实现含噪高光谱图像的快速解混，在图像去噪过程中，受各种去噪方法的影响，原始图像可能会发生部分改变并产生新的实验误差，这可能对进一步的解混问题产生干扰，导致低的解混精度对目标检测和亚像元分类等产生不利影响，因此，本发明提出一种含噪高光谱图像快速解混方法，以解决现有技术中的不足之处。

发明内容

针对上述问题，本发明提出一种含噪高光谱图像快速解混方法，本发明方法直接对原始的高光谱图像的端元进行提取，并对端元光谱进行去噪，然后利用去噪后的端元进行丰度估计，实现高光谱图像的解混，该方法避免了现有方法在第一步中所产生的误差，同时也解决了原始图像去噪过程中重要信息丢失而降低解混精度的问题；此外，本方法对端元光谱去噪，数据量小，解决了图像去噪后解混技术效率低的问题，本方法不仅可以提高解混精度，操作简单而且处理速度快。

本发明提出一种含噪高光谱图像快速解混方法，包括以下步骤：

步骤一：采用VCA算法提取含噪声的端元光谱，基于凸面几何理论在假设数据中存在纯像元的情况下提取端元，利用凸锥来对数据进行建模，在超平面上的投影为由端元作为顶点的单体，将数据投影到选择的超平面上，VCA将所有影像上的像元投影到随机方向上，并将具有最大投影的像元作为第一个端元，通过迭代地将数据投影到与由已提取端元构成的子空间正交的方向上的方法来提取其余端元，被极限投影所对应的像元作为新的端元；

步骤二：对提取的端元光谱进行SSA算法去噪，利用端元提取算法提取端元光谱Y，计算其轨迹矩阵X，然后进行奇异值分解、分组和重构；

步骤三：利用去噪后的端元光谱根据最小二乘法进行丰度估计，丰度估计表达式如公式(1)所示：

其中，E代表去噪后的端元光谱，x_i为原始图像。

进一步改进在于：所述步骤一中采用VCA算法提取含噪声的端元光谱的具体算法过程为：

假定线性混合模型表达式如公式(2)所示：

x＝Aγs+ε (2)

其中，式中γ为比例因子，表示由于地表起伏引起的光照变化。

进一步改进在于：所述公式(2)中由于丰度的物理限制，s∈Δ_p,Δ_p是单体，每个波段可以作为L维空间的坐标轴，每个像元为L维欧氏空间的一个向量，S_x＝{x∈R^l:x＝As,s∈Δ_p}也是单体；

C_p＝{x∈R^l:x＝Aγs,s∈Δ_p,γ≥0}，由于比例参数γ的存在形成一个凸锥体。

进一步改进在于：所述丰度的物理限制包括非负限制、和为一限制。

进一步改进在于：所述凸锥体C_p在适当选择的超平面上的投影为单体S_x所对应的端点。

进一步改进在于：所述步骤二中利用端元提取算法提取端元光谱Y，计算其轨迹矩阵X的具体过程为：

将一维端元光谱信息Y_T＝(y₁,…y_T)转化为其轨迹矩阵，然后根据窗口长度L计算光谱信息的轨迹矩阵X，轨迹矩阵X为L×K阶，其中K＝N-L+1，轨迹矩阵X的计算公式如公式(3)所示：

其中，X_i＝(x_i,x_i+1,…,x_i+L-1)；轨迹矩阵X的(i,j)处的元素x_ij＝x_i+j-1；所有反对角线上的元素相等。

进一步改进在于：所述步骤二中奇异值分解具体过程为：

定义矩阵C_x＝XX^T,X^T为X的转置矩阵，然后计算矩阵C_x的特征值λ_i及特征向量U_i，按照由大到小排列特征值，其特征值依次为λ₁≥…≥λ_L≥0，对应的特征向量为U₁,…,U_L，再设d＝L^*,L^*＝min{L,K},

U_i与V_i是轨迹矩阵的左右特征向量，U称为时间经验正交函数，V称为时间主成分，轨迹矩阵X可以由初等矩阵合成，合成公式如公式(4)所示：

X＝X₁+…+X_d (4)

初等矩阵为

其中

为X的奇异谱，其最大特征值对应最大特征向量，代表地物光谱，较小的特征值对应的特征向量，代表噪声。

进一步改进在于：所述步骤二中分组具体过程为：

将初等矩阵X_i的下标{1，2，…d}分成p个不相交的子集I₁,I₂,…I_p，设I＝{i₁,i₂,…i_m}，则合成矩阵X_I＝X_i1+X_i2+…+X_im，计算集合I＝I₁,I₂,…,I_p的每个合成矩阵，将公式(4)的分解式用公式(5)表达：

X＝X_I1+X_I2+…+X_Ip (5)

其中，选取集合I₁,…I_p的过程称为分组。

进一步改进在于：所述步骤二中重构具体过程为：

将公式(5)分组分解得到的每个X_Ij转换为长度为N的新光谱数据，设Y是L×K的矩阵，其元素为y_ij，其中1≤i≤L,1≤j≤K,L^*＝min{L,K},K^*＝max{L,K}且N＝L+K-1，根据对角平均公式将矩阵Y转换为y₁,…,y_N的光谱信息，对角平均公式如公式(6)所示：

对于合成的矩阵X_Ik，由公式(6)可以生成一个重构序列

则初始光谱信息Y_T＝(y₁,…y_T)分解为总和为p的重构序列，表达式为

本发明的有益效果为：本发明方法直接对原始的高光谱图像的端元进行提取，并对端元光谱进行去噪，然后利用去噪后的端元进行丰度估计，实现高光谱图像的解混，该方法避免了现有方法在第一步中所产生的误差，同时也解决了原始图像去噪过程中重要信息丢失而降低解混精度的问题；此外，本方法对端元光谱去噪，数据量小，解决了图像去噪后解混技术效率低的问题，在端元提取中选用了顶点成分分析算法(VCA),该算法在效率和精度上都表现出明显的优势，端元光谱去噪方法选用了基于一维频谱域的奇异谱分析算法(SSA)，该算法能够较好的平滑光谱曲线中的噪声且保留原始特征，基于一维频谱域的去噪方法所需时间短，可以提高效率，本方法不仅可以提高解混精度，操作简单而且处理速度快。

附图说明

图1为本发明含噪高光谱图像快速解混方法流程示意图。

图2为本发明含噪高光谱图像快速解混方法示意图。

图3为本发明VCA算法原理示意图。

图4为本发明实施例中多种解混方法进行高光谱图像解混后的解混效果对比示意图。

具体实施方式

为了加深对本发明的理解，下面将结合实施例对本发明做进一步详述，本实施例仅用于解释本发明，并不构成对本发明保护范围的限定。

根据图1、2、3、4所示，本实施例提出一种含噪高光谱图像快速解混方法，包括以下步骤：

步骤一：采用VCA算法提取含噪声的端元光谱，基于凸面几何理论在假设数据中存在纯像元的情况下提取端元，利用凸锥来对数据进行建模，在超平面上的投影为由端元作为顶点的单体，将数据投影到选择的超平面上，VCA将所有影像上的像元投影到随机方向上，并将具有最大投影的像元作为第一个端元，通过迭代地将数据投影到与由已提取端元构成的子空间正交的方向上的方法来提取其余端元，被极限投影所对应的像元作为新的端元；

假定线性混合模型表达式如公式(2)所示：

x＝Aγs+ε (2)

其中，式中γ为比例因子，表示由于地表起伏引起的光照变化；

公式(2)中由于丰度的物理限制，包括非负限制、和为一限制，s∈Δ_p,Δ_p是单体，每个波段可以作为L维空间的坐标轴，每个像元为L维欧氏空间的一个向量，S_x＝{x∈R^l:x＝As,s∈Δ_p}也是单体；

C_p＝{x∈R^l:x＝Aγs,s∈Δ_p,γ≥0}，由于比例参数γ的存在形成一个凸锥体，凸锥体C_p在适当选择的超平面上的投影为单体S_x所对应的端点，如如3所示，本实施例中S_p＝{y∈R^l:y＝x/(x^Tμ),r∈C_p}为凸锥在平面x^Tμ＝1上的投影，μ的选择需要确保没有观测向量与其正交，在确定单体S_p之后，VCA迭代地将数据投影到与之前已提取的端元构成的子空间正交的方向上，图3表示了VCA应用于单体S_p的两次迭代过程，在每一次迭代中，数据被投影到第一个方向f₁上，投影的极值对应于端元m_a，在下一次迭代中，端元m_b对应于将数据投影到f₂上的极值，f₂与m_a正交，算法迭代执行直到确定数目的端元被提取出来；

步骤二：对提取的端元光谱进行SSA算法去噪，利用端元提取算法提取端元光谱Y，计算其轨迹矩阵X，然后进行奇异值分解、分组和重构；

利用端元提取算法提取端元光谱Y，计算其轨迹矩阵X的具体过程为：

将一维端元光谱信息Y_T＝(y₁,…y_T)转化为其轨迹矩阵，然后根据窗口长度L计算光谱信息的轨迹矩阵X，轨迹矩阵X为L×K阶，其中K＝N-L+1，轨迹矩阵X的计算公式如公式(3)所示：

其中，X_i＝(x_i,x_i+1,…,x_i+L-1)；轨迹矩阵X的(i,j)处的元素x_ij＝x_i+j-1；所有反对角线上的元素相等；

奇异值分解具体过程为：

定义矩阵C_x＝XX^T,X^T为X的转置矩阵，然后计算矩阵C_x的特征值λ_i及特征向量U_i，按照由大到小排列特征值，其特征值依次为λ₁≥…≥λ_L≥0，对应的特征向量为U₁,…,U_L，再设d＝L^*,L^*＝min{L,K},

U_i与V_i是轨迹矩阵的左右特征向量，U称为时间经验正交函数，V称为时间主成分，轨迹矩阵X可以由初等矩阵合成，合成公式如公式(4)所示：

X＝X₁+…+X_d (4)

初等矩阵为

其中

为X的奇异谱，其最大特征值对应最大特征向量，代表地物光谱，较小的特征值对应的特征向量，代表噪声；

分组具体过程为：

将初等矩阵X_i的下标{1，2，…d}分成p个不相交的子集I₁,I₂,…I_p，设I＝{i₁,i₂,…i_m}，则合成矩阵X_I＝X_i1+X_i2+…+X_im，计算集合I＝I₁,I₂,…,I_p的每个合成矩阵，将公式(4)的分解式用公式(5)表达：

X＝X_I1+X_I2+…+X_Ip (5)

其中，选取集合I₁,…I_p的过程称为分组；

重构具体过程为：

将公式(5)分组分解得到的每个X_Ij转换为长度为N的新光谱数据，设Y是L×K的矩阵，其元素为y_ij，其中1≤i≤L,1≤j≤K,L^*＝min{L,K},K^*＝max{L,K}且N＝L+K-1，根据对角平均公式将矩阵Y转换为y₁,…,y_N的光谱信息，对角平均公式如公式(6)所示：

对于合成的矩阵X_Ik，由公式(6)可以生成一个重构序列

则初始光谱信息Y_T＝(y₁,…y_T)分解为总和为p的重构序列，表达式为

步骤三：利用去噪后的端元光谱根据最小二乘法进行丰度估计，丰度估计表达式如公式(1)所示：

其中，E代表去噪后的端元光谱，x_i为原始图像。

本实施例中利用本发明方法和基于对整幅图像的去噪再解混的方法以及图像不去噪的解混方法进行对比，得出丰度图对比如图4所示，精度评定结果如表1所示：

表1本方法与基于图像去噪再解混方法、图像不去噪解混方法效果对比

其中，SAD(Spectral angular distance)代表某一种方法提取端元光谱与标准端元光谱之间的光谱角距，该值越小越好；RMSE(Root Mean Square Error)代表某一种方法估计丰度与标准丰度之间的均方根误差，该值越小越好；T(Time)代表各方法的运行时间；

从图4和表1的数据可以看出，基于本发明方法的含噪高光谱图像的解混平均只需0.1s即可完成，而基于对整幅图像去噪再解混方法5s才能完成；从图像中可以看出本方法最终的解混结果与标准图像更为接近，从表1中也可以看出本方法在SAD和RMSE效果评价中均取得最小的误差；因此，本方法在时效性和精确性上均好于对比方法，即对整幅图像去噪再解混的方法。

本发明方法直接对原始的高光谱图像的端元进行提取，并对端元光谱进行去噪，然后利用去噪后的端元进行丰度估计，实现高光谱图像的解混，该方法避免了现有方法在第一步中所产生的误差，同时也解决了原始图像去噪过程中重要信息丢失而降低解混精度的问题；此外，本方法对端元光谱去噪，数据量小，解决了图像去噪后解混技术效率低的问题，在端元提取中选用了顶点成分分析算法(VCA),该算法在效率和精度上都表现出明显的优势，端元光谱去噪方法选用了基于一维频谱域的奇异谱分析算法(SSA)，该算法能够较好的平滑光谱曲线中的噪声且保留原始特征，基于一维频谱域的去噪方法所需时间短，可以提高效率，本方法不仅可以提高解混精度，操作简单而且处理速度快。

以上显示和描述了本发明的基本原理、主要特征和优点。本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实施例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。本发明要求保护范围由所附的权利要求书及其等效物界定。

Claims (9)

1.一种含噪高光谱图像快速解混方法，其特征在于，包括以下步骤：

步骤一：采用VCA算法提取含噪声的端元光谱，基于凸面几何理论在假设数据中存在纯像元的情况下提取端元，利用凸锥来对数据进行建模，在超平面上的投影为由端元作为顶点的单体，将数据投影到选择的超平面上，VCA将所有影像上的像元投影到随机方向上，并将具有最大投影的像元作为第一个端元，通过迭代地将数据投影到与由已提取端元构成的子空间正交的方向上的方法来提取其余端元，被极限投影所对应的像元作为新的端元；

步骤二：对提取的端元光谱进行SSA算法去噪，利用端元提取算法提取端元光谱Y，计算其轨迹矩阵X，然后进行奇异值分解、分组和重构；

步骤三：利用去噪后的端元光谱根据最小二乘法进行丰度估计，丰度估计表达式如公式(1)所示：

其中，E代表去噪后的端元光谱，x_i为原始图像。

2.根据权利要求1所述的一种含噪高光谱图像快速解混方法，其特征在于：所述步骤一中采用VCA算法提取含噪声的端元光谱的具体算法过程为：

假定线性混合模型表达式如公式(2)所示：

x＝Aγs+ε (2)

其中，式中γ为比例因子，表示由于地表起伏引起的光照变化。

3.根据权利要求2所述的一种含噪高光谱图像快速解混方法，其特征在于：所述公式(2)中由于丰度的物理限制，s∈Δ_p,Δ_p是单体，每个波段可以作为L维空间的坐标轴，每个像元为L维欧氏空间的一个向量，S_x＝{x∈R^l:x＝As,s∈Δ_p}也是单体；

C_p＝{x∈R^l:x＝Aγs,s∈Δ_p,γ≥0}，由于比例参数γ的存在形成一个凸锥体。

4.根据权利要求3所述的一种含噪高光谱图像快速解混方法，其特征在于：所述丰度的物理限制包括非负限制、和为一限制。

5.根据权利要求3所述的一种含噪高光谱图像快速解混方法，其特征在于：所述凸锥体C_p在适当选择的超平面上的投影为单体S_x所对应的端点。

6.根据权利要求1所述的一种含噪高光谱图像快速解混方法，其特征在于：所述步骤二中利用端元提取算法提取端元光谱Y，计算其轨迹矩阵X的具体过程为：

将一维端元光谱信息Y_T＝(y₁,…y_T)转化为其轨迹矩阵，然后根据窗口长度L计算光谱信息的轨迹矩阵X，轨迹矩阵X为L×K阶，其中K＝N-L+1，轨迹矩阵X的计算公式如公式(3)所示：

其中，X_i＝(x_i,x_i+1,…,x_i+L-1)；轨迹矩阵X的(i,j)处的元素x_ij＝x_i+j-1；所有反对角线上的元素相等。

7.根据权利要求6所述的一种含噪高光谱图像快速解混方法，其特征在于：所述步骤二中奇异值分解具体过程为：

定义矩阵C_x＝XX^T,X^T为X的转置矩阵，然后计算矩阵C_x的特征值λ_i及特征向量U_i，按照由大到小排列特征值，其特征值依次为λ₁≥…≥λ_L≥0，对应的特征向量为U₁,…,U_L，再设d＝L^*,L^*＝min{L,K},

U_i与V_i是轨迹矩阵的左右特征向量，U称为时间经验正交函数，V称为时间主成分，轨迹矩阵X可以由初等矩阵合成，合成公式如公式(4)所示：

X＝X₁+…+X_d (4)

初等矩阵为

其中rank(X_i)＝1,

为X的奇异谱，其最大特征值对应最大特征向量，代表地物光谱，较小的特征值对应的特征向量，代表噪声。

8.根据权利要求7所述的一种含噪高光谱图像快速解混方法，其特征在于：所述步骤二中分组具体过程为：

将初等矩阵X_i的下标{1，2，…d}分成p个不相交的子集I₁,I₂,…I_p，设I＝{i₁,i₂,…i_m}，则合成矩阵X_I＝X_i1+X_i2+…+X_im，计算集合I＝I₁,I₂,…,I_p的每个合成矩阵，将公式(4)的分解式用公式(5)表达：

X＝X_I1+X_I2+…+X_Ip (5)

其中，选取集合I₁,…I_p的过程称为分组。

9.根据权利要求8所述的一种含噪高光谱图像快速解混方法，其特征在于：所述步骤二中重构具体过程为：

将公式(5)分组分解得到的每个X_Ij转换为长度为N的新光谱数据，设Y是L×K的矩阵，其元素为y_ij，其中1≤i≤L,1≤j≤K,L^*＝min{L,K},K^*＝max{L,K}且N＝L+K-1，根据对角平均公式将矩阵Y转换为y₁,…,y_N的光谱信息，对角平均公式如公式(6)所示：

对于合成的矩阵X_Ik，由公式(6)可以生成一个重构序列

则初始光谱信息Y_T＝(y₁,…y_T)分解为总和为p的重构序列，表达式为

Images