CN102542547A - 一种基于光谱约束的高光谱图像融合方法 - Google Patents

Abstract

本发明一种基于光谱约束的高光谱图像融合方法，它是一种基于非负矩阵分解方法的高光谱图像和高空间分辨率可见光图像的融合方法。该方法共有三大步骤，步骤一：计算机在MATLAB7.1环境下读取高光谱图像数据；步骤二：建立基于非负矩阵分解方法的融合模型；步骤三：利用梯度下降法对上述融合模型求解极值函数，并得到融合结果。本发明克服了现有技术不宜直接移植到高光谱图像融合、光谱保持性不好的缺点，取得了较好的融合结果，在高光谱图像融合技术领域里具有广阔的应用前景。

Description

一种基于光谱约束的高光谱图像融合方法

(一)技术领域：

本发明涉及一种基于光谱约束的高光谱图像融合方法，属于高光谱遥感图像融合技术领域。

(二)背景技术：

高光谱遥感是高光谱分辨率遥感(Hyperspectral Remote Sensing)的简称，是目前遥感技术的前沿，它是一项集空间科学、计算机技术、探测器和精密光学机械、微弱信号检测、高速信息处理技术于一体的综合性技术，利用很多很窄的电磁波波段从感兴趣的物体中获取有关数据。在成像过程中，高光谱遥感技术利用成像光谱仪以纳米级的光谱分辨率，以几十或几百个波段同时对地表地物成像，能够获得地物的连续光谱信息，实现了地物空间信息、辐射信息、光谱信息的同步获取，大大提高了对地面目标解析的可靠程度。高光谱数据的基本特点是光谱分辨率高、数据量大、信息量丰富，用上百个连续且窄的光谱波段描述一个像元，能够产生一条完整而连续的光谱曲线。高光谱图像的光谱分辨率高，但空间分辨率较低，而可见光图像具有较高的空间分辨率，两者互补，采用一定的融合算法处理，可以获取高质量融合图像，既保持图像的光谱特性又能获取较高的空间分辨率。

目前存在的图像融合方法，大多为在像素级别针对多光谱图像的融合方法。其中经典的方法有直接加权、主成分分析(Principal Component Analysis，PCA)、小波变换等方法。它们在多光谱图像融合领域被研究的比较充分，但是在高光谱图像融合领域，相关的研究还比较少，多是将现有的方法直接移植，由于移植过程中一般不对数据的光谱做针对性的研究，结果不是很理想。本发明针对上述情况，提出一种基于光谱约束的方法以实现高光谱图像融合，融合效果较好。

(三)发明内容：

1、目的：本发明的目的是提供一种基于光谱约束的高光谱图像融合方法，该方法克服了现有技术的不足，它是基于非负矩阵分解算法，并引入光谱约束来进行高光谱图像融合。

2、技术方案：本发明是通过以下技术方案实现的：

本发明一种基于光谱约束的高光谱图像融合方法，它包括如下步骤：

步骤一：计算机读取数据。首先计算机在MATLAB 7.1下读取成像光谱仪所得到的高光谱图像数据(本发明中高光谱数据来源于中国HJ-1A卫星搭载的成像光谱仪拍摄的青海省托素湖遥感图像，含有100个谱段，空间分辨率为100m)，再读取已经手工配准过的高分辨率图像数据(本发明数据来源于Google-Map提供的高分辨率影像)。

本发明中的高光谱图像数据为100谱段*256像素*200像素大小的三维数据立方体，其中100表示谱段数，256*200表示每一谱段图像大小。以下用V表示融合前的图像集，V_f表示融合后的的图像集，他们的大小都是N*100。其中，N＝row*col，row和col分别表示图像的行数和列数，V的第i行包含这一像素点的光谱信息。

步骤二：建立目标函数。首先建立高光谱图像融合的目标函数，它包括两项：

第一项为高光谱图像进行非负矩阵分解的目标函数，形式如下：

$F(W,H)=\frac{1}{2}{\left|\right|V-\mathrm{WH}\left|\right|}^2$

s.t.W≥O，H≥O

式中，V是将100个波段的高光谱图像作为列向量组成的图像矩阵；W是矩阵V进行非负矩阵分解得到的特征基；H是权重系数矩阵。

第二项为光谱约束项S(W₀，H)，该式用于实现对高光谱数据的光谱约束。形式如下：

$S(W_0,H)=\mathrm{tr}(W_0{\mathrm{HH}}^T{W}_0^T\·{\mathrm{VV}}^T)-\mathrm{tr}(W_0{\mathrm{HV}}^T\·W_0{\mathrm{HV}}^T)$

其中，W₀＝W+P；P的每一列均为高分辨率图像，和W的尺寸相同；tr(*)表示矩阵*的迹；·表示点乘。

将光谱约束项S(W₀，H)加入第一项，得到最终的目标函数。形式如下：

$\min G(W,H)=F(W,H)+\mathrm{\βS}(W_0,H)$

$=\frac{1}{2}{\left|\right|V-\mathrm{WH}\left|\right|}^2+\mathrm{\β}\{\mathrm{tr}(W_0{\mathrm{HH}}^T{W}_0^T\·{\mathrm{VV}}^T)-\mathrm{tr}(W_0{\mathrm{HV}}^T\·W_0{\mathrm{HV}}^T)\}$

s.t.W≥O，H≥O

式中，β为约束项权重，用于协调融合图像的光谱分辨率和空间分辨率。

步骤三：采用梯度下降法求目标函数的极值，得到此时的W和H，于是融合后的图像为

V_f＝(W+P)H

其中，步骤三中所述的采用梯度下降法求目标函数的极值，其具体的求解步骤如下：

(1)矩阵W和H初始化：

将融合前的图像矩阵V进行VCA初始化，得到矩阵W和H的初值，其中W的列数取5。

(2)n＝1，开始循环。

(3)计算第k+1次迭代得到的特征基矩阵W和权重系数矩阵H：

$W_0^{(k+1)}=W^{\left(k\right)}+P$

$V_0^{(k+1)}=W_0^{(k+1)}{H}^{\left(k\right)}$

以下分别用B和D表示目标函数G(W，H)对W和H的偏导数，其中(*)^T表示矩阵*的转置矩阵；diag(*)表示矩阵*的对角矩阵；max(*，0)表示取矩阵*中元素和0中的最大值，用于保证矩阵*中的元素非负。

$B_1^{(k+1)}={\left[W_0^{(k+1)}\right]}^T\mathrm{diag}\left({\mathrm{VV}}^T\right)W_0^{(k+1)}{H}^{\left(k\right)}$

$B_2^{(k+1)}={\left[W_0^{(k+1)}\right]}^T\mathrm{diag}\left(V_0^{(k+1)}{V}^T\right)V$

$B^{(k+1)}={\left[W^{\left(k\right)}\right]}^T(W^{\left(k\right)}{H}^{\left(k\right)}-V)+2\mathrm{\β}{B}_1^{(k+1)}-2\mathrm{\β}{B}_2^{(k+1)}$

$B_0^{(k+1)}=\max (B^{(k+1)},0)$

$D_1^{(k+1)}=\mathrm{diag}\left({\mathrm{VV}}^T\right)W_0^{(k+1)}{H}^{\left(k\right)}{\left[H^{\left(k\right)}\right]}^T$

$D_2^{(k+1)}=\mathrm{diag}\left(V_0^{(k+1)}{V}^T\right)V{\left[H^{\left(k\right)}\right]}^T$

$D^{(k+1)}=(W^{\left(k\right)}{H}^{\left(k\right)}-V){\left[H^{\left(k\right)}\right]}^T+2\mathrm{\β}{D}_1^{(k+1)}-2\mathrm{\β}{D}_2^{(k+1)}$

$D_0^{(k+1)}=\max (D^{(k+1)},0)$

W和H的更新方式如下，其中α为梯度下降法学习率。

$W^{(k+1)}=\max (W^{\left(k\right)}-\mathrm{\α}{B}_0^{(k+1)},0)$

$H^{(k+1)}=\max (H^{\left(k\right)}-\mathrm{\α}{D}_0^{(k+1)},0)$

(4)迭代停止条件。当在相邻两次迭代的过程中变化不大时停止迭代，本发明中停止迭代的条件是：

$\frac{\left|\right|W^{\left(k\right)}{H}^{\left(k\right)}-W^{(k+1)}{H}^{(k+1)}\left|\right|}{\left|\right|W^{\left(k\right)}{H}^{\left(k\right)}\left|\right|}<\mathrm{tol}$

其中W^(k+1)H^(k+1)表示当前迭代结果，W^(k)H^(k)表示前一次迭代结果。本发明中，取tol＝10^-3。当停止迭代的条件满足时转去第(6)步，否则转第(5)步。

(5)n＝n+1，转第(3)步，直到n＞N，其中N为迭代次数，取N＝1000。

(6)计算融合后的图像V_f＝(W+P)H，并将获得的高光谱图像数据V_f存入新的数据立方体中，至此完成对高光谱图像与高分辨率可见光图像的融合工作。

3、优点及功效。本发明的优点是：针对常见图像融合方法不宜直接移植到高光谱图像融合，而现有的基于非负矩阵分解的融合手段对于光谱曲线的保持效果不佳，本发明中通过添加光谱约束对高光谱数据的光谱特性进行约束，利用融合前后的光谱曲线的相似度进行约束，衡量的标准是融合前后两条光谱曲线的光谱角度，最终得到了图像和光谱质量都较好的融合结果。

通过随机选择的实际两像素点融合前后的光谱曲线可看出，使用本方法进行高光谱图像融合，光谱保持性非常好，事实上，几乎所有的像素点都保持了融合前光谱曲线的基本特征，利用光谱角度进行光谱约束，得到的融合图像与源高光谱图像各像素点中的光谱角度余弦平均值为0.9938，同时融合图像的视觉质量也比较好，因而利用发明的添加光谱约束项的非负矩阵分解方法处理高光谱图像融合问题效果比较好。

(四)附图说明：

图1本发明所述方法融合流程框图

图2(a)本发明实验高光谱数据(129，100)点处融合前后的光谱曲线示意图

其中，细虚线为融合前光谱曲线，粗实线为融合后光谱曲线

图2(b)本发明实验高光谱数据(200，200)点处融合前后的光谱曲线示意图

其中，细虚线为融合前光谱曲线，粗实线为融合后光谱曲线

(五)具体实施方式：

为了更好地理解本发明的技术方案，以下结合附图对本发明的实施方式作进一步描述：

本发明在MATLAB7.1语言环境下实现。计算机读取高光谱遥感图像数据和配准好的高分辨率图像数据，首先基于非负矩阵分解建立初步的目标函数，然后在此基础上添加光谱约束项，得到最终的融合模型，对此融合模型求极值函数得到基图像和权值矩阵，用高分辨率图像乘以标准图像与高分辨率图像加权后的图像，所得结果即为所求的高光谱融合图像。

本发明是一种基于光谱约束的高光谱图像融合方法，该方法的流程见图1所示，计算机配置采用：Intel dual-core处理器，主频1.6GHz，内存2GB，操作系统为Windows Vista Home Basic，编程环境为Matlab7.1。该方法包括以下步骤：

步骤一、用计算机读取数据。计算机在MATLAB7.1环境下读取高光谱图像数据和高分辨率图像数据。

步骤二、建立基于光谱约束的非负矩阵分解方法的融合模型。本发明的实验中选择的参数如下：约束项权重β＝10^(-2)，梯度下降法学习率α＝10^(-5)。

步骤三、对所得融合模型求解极值函数，得到融合结果。利用梯度下降法迭代求解，以下是具体实施过程。

(1)矩阵W和H初始化：

将融合前的图像矩阵V进行VCA初始化，得到矩阵W和H的初值，其中W的列数取5。

(2)n＝1，开始循环。

(3)计算第k+1次迭代得到的特征基矩阵W和权重系数矩阵H：

$W_0^{(k+1)}=W^{\left(k\right)}+P$

$V_0^{(k+1)}=W_0^{(k+1)}{H}^{\left(k\right)}$

以下分别用B和D表示目标函数G(W，H)对W和H的偏导数，其中(*)^T表示矩阵*的转置矩阵；diag(*)表示矩阵*的对角矩阵；max(*，0)表示取矩阵*中元素和0中的最大值，用于保证矩阵*中的元素非负。

$B_1^{(k+1)}={\left[W_0^{(k+1)}\right]}^T\mathrm{diag}\left({\mathrm{VV}}^T\right)W_0^{(k+1)}{H}^{\left(k\right)}$

$B_2^{(k+1)}={\left[W_0^{(k+1)}\right]}^T\mathrm{diag}\left(V_0^{(k+1)}{V}^T\right)V$

$B^{(k+1)}={\left[W^{\left(k\right)}\right]}^T(W^{\left(k\right)}{H}^{\left(k\right)}-V)+2\mathrm{\&beta;}{B}_1^{(k+1)}-2\mathrm{\&beta;}{B}_2^{(k+1)}$

$B_0^{(k+1)}=\max (B^{(k+1)},0)$

$D_1^{(k+1)}=\mathrm{diag}\left({\mathrm{VV}}^T\right)W_0^{(k+1)}{H}^{\left(k\right)}{\left[H^{\left(k\right)}\right]}^T$

$D_2^{(k+1)}=\mathrm{diag}\left(V_0^{(k+1)}{V}^T\right)V{\left[H^{\left(k\right)}\right]}^T$

$D^{(k+1)}=(W^{\left(k\right)}{H}^{\left(k\right)}-V){\left[H^{\left(k\right)}\right]}^T+2\mathrm{\&beta;}{D}_1^{(k+1)}-2\mathrm{\&beta;}{D}_2^{(k+1)}$

$D_0^{(k+1)}=\max (D^{(k+1)},0)$

W和H的更新方式如下，其中α为梯度下降法学习率。

$W^{(k+1)}=\max (W^{\left(k\right)}-\mathrm{\&alpha;}{B}_0^{(k+1)},0)$

$H^{(k+1)}=\max (H^{\left(k\right)}-\mathrm{\&alpha;}{D}_0^{(k+1)},0)$

(4)迭代停止条件。当在相邻两次迭代的过程中变化不大时停止迭代，本发明中停止迭代的条件是：

$\frac{\left|\right|W^{\left(k\right)}{H}^{\left(k\right)}-W^{(k+1)}{H}^{(k+1)}\left|\right|}{\left|\right|W^{\left(k\right)}{H}^{\left(k\right)}\left|\right|}<\mathrm{tol}$

其中W^(k+1)H^(k+1)表示当前迭代结果，W^(k)H^(k)表示前一次迭代结果。本发明中，取tol＝10^-3。当停止迭代的条件满足时转去第(6)步，否则转第(5)步。

(5)n＝n+1，转第(3)步，直到n＞N，其中N为迭代次数，取N＝1000。

(6)计算融合后的图像V_f＝(W+P)H，并将获得的高光谱图像数据V_f存入新的数据立方体中，至此完成对高光谱图像与高分辨率可见光图像的融合工作，融合结果可以直接用于后续的目标检测等其它图像处理工作。

为了验证本发明的有效性，使用上述方法进行实验，取得了较好的融合效果。本发明中实验所用数据为来自于中国HJ-1A星载成像光谱仪得到的三维数据立方体，分析图2(a)、图2(b)可见，利用本发明的方法，得到了比较理想的融合结果，光谱保持性非常理想，基本上融合图像像素点的融合前后光谱曲线一致。本发明在保持高光谱图像数据的光谱分辨率的基础上提高了数据的空间分辨率。

从实验结果来看，本发明中的方法很好的解决了高光谱图像和高空间分辨率可见光图像的融合问题，特别是在光谱保持性方面，因而此方法可以应用到高光谱图像融合中，具有广阔的应用前景和价值。

Claims (2)

1.一种基于光谱约束的高光谱图像融合方法，其特征在于：该方法具体步骤如下：

步骤一：计算机读取数据；首先计算机在MATLAB 7.1下读取成像光谱仪所得到的高光谱图像数据，再读取已经手工配准过的高分辨率图像数据；

高光谱图像数据为100谱段*256像素*200像素大小的三维数据立方体，其中100表示谱段数，256*200表示每一谱段图像大小；以下用V表示融合前的图像集，V_f表示融合后的的图像集，他们的大小都是N*100；其中，N＝row*col，row和col分别表示图像的行数和列数，V的第i行包含这一像素点的光谱信息；

步骤二：建立目标函数；首先建立高光谱图像融合的目标函数，它包括两项：

第一项为高光谱图像进行非负矩阵分解的目标函数，形式如下：

$F(W,H)=\frac{1}{2}{\left|\right|V-\mathrm{WH}\left|\right|}^2$

s.t.W≥O，H≥O

式中，V是将100个波段的高光谱图像作为列向量组成的图像矩阵；W是矩阵V进行非负矩阵分解得到的特征基；H是权重系数矩阵；

第二项为光谱约束项S(W₀，H)，该式用于实现对高光谱数据的光谱约束；形式如下：

$S(W_0,H)=\mathrm{tr}(W_0{\mathrm{HH}}^T{W}_0^T\&CenterDot;{\mathrm{VV}}^T)-\mathrm{tr}(W_0{\mathrm{HV}}^T\&CenterDot;W_0{\mathrm{HV}}^T)$

其中，W₀＝W+P；P的每一列均为高分辨率图像，和W的尺寸相同；tr(*)表示矩阵*的迹；·表示点乘；

将光谱约束项S(W₀，H)加入第一项，得到最终的目标函数；形式如下：

$\min G(W,H)=F(W,H)+\mathrm{\&beta;S}(W_0,H)$

$=\frac{1}{2}{\left|\right|V-\mathrm{WH}\left|\right|}^2+\mathrm{\&beta;}\{\mathrm{tr}(W_0{\mathrm{HH}}^T{W}_0^T\&CenterDot;{\mathrm{VV}}^T)-\mathrm{tr}(W_0{\mathrm{HV}}^T\&CenterDot;W_0{\mathrm{HV}}^T)\}$

s.t.W≥O，H≥O

式中，β为约束项权重，用于协调融合图像的光谱分辨率和空间分辨率；

步骤三：采用梯度下降法求目标函数的极值，得到此时的W和H，于是融合后的图像为

V_f＝(W+P)H  。

2.根据权利要求1所述的一种基于光谱约束的高光谱图像融合方法，其特征在于：步骤三中所述的采用梯度下降法求目标函数的极值，其具体的求解步骤如下：

(1)矩阵W和H初始化：

将融合前的图像矩阵V进行VCA初始化，得到矩阵W和H的初值，其中W的列数取5；

(2)n＝1，开始循环；

(3)计算第k+1次迭代得到的特征基矩阵W和权重系数矩阵H：

$W_0^{(k+1)}=W^{\left(k\right)}+P$

$V_0^{(k+1)}=W_0^{(k+1)}{H}^{\left(k\right)}$

以下分别用B和D表示目标函数G(W，H)对W和H的偏导数，其中(*)^τ表示矩阵*的转置矩阵；diag(*)表示矩阵*的对角矩阵；max(*，0)表示取矩阵*中元素和0中的最大值，用于保证矩阵*中的元素非负；

$B_1^{(k+1)}={\left[W_0^{(k+1)}\right]}^T\mathrm{diag}\left({\mathrm{VV}}^T\right)W_0^{(k+1)}{H}^{\left(k\right)}$

$B_2^{(k+1)}={\left[W_0^{(k+1)}\right]}^T\mathrm{diag}\left(V_0^{(k+1)}{V}^T\right)V$

$B^{(k+1)}={\left[W^{\left(k\right)}\right]}^T(W^{\left(k\right)}{H}^{\left(k\right)}-V)+2\mathrm{\&beta;}{B}_1^{(k+1)}-2\mathrm{\&beta;}{B}_2^{(k+1)}$

$B_0^{(k+1)}=\max (B^{(k+1)},0)$

$D_1^{(k+1)}=\mathrm{diag}\left({\mathrm{VV}}^T\right)W_0^{(k+1)}{H}^{\left(k\right)}{\left[H^{\left(k\right)}\right]}^T$

$D_2^{(k+1)}=\mathrm{diag}\left(V_0^{(k+1)}{V}^T\right)V{\left[H^{\left(k\right)}\right]}^T$

$D^{(k+1)}=(W^{\left(k\right)}{H}^{\left(k\right)}-V){\left[H^{\left(k\right)}\right]}^T+2\mathrm{\&beta;}{D}_1^{(k+1)}-2\mathrm{\&beta;}{D}_2^{(k+1)}$

$D_0^{(k+1)}=\max (D^{(k+1)},0)$

W和H的更新方式如下，其中α为梯度下降法学习率；

$W^{(k+1)}=\max (W^{\left(k\right)}-\mathrm{\&alpha;}{B}_0^{(k+1)},0)$

$H^{(k+1)}=\max (H^{\left(k\right)}-\mathrm{\&alpha;}{D}_0^{(k+1)},0)$

(4)迭代停止条件：当在相邻两次迭代的过程中变化不大时停止迭代，停止迭代的条件是：

$\frac{\left|\right|W^{\left(k\right)}{H}^{\left(k\right)}-W^{(k+1)}{H}^{(k+1)}\left|\right|}{\left|\right|W^{\left(k\right)}{H}^{\left(k\right)}\left|\right|}<\mathrm{tol}$

其中，W^(k+1)H^(k+1)表示当前迭代结果，W^(k)H^(k)表示前一次迭代结果；这里，取tol＝10^-3；当停止迭代的条件满足时转去第(6)步，否则转第(5)步；

(5)n＝n+1，转第(3)步，直到n＞N，其中N为迭代次数，取N＝1000；

(6)计算融合后的图像V_f＝(W+P)H，并将获得的高光谱图像数据V_f存入新的数据立方体中，至此完成对高光谱图像与高分辨率可见光图像的融合工作。

Images