CN101908211A - 一种基于变分方法的高光谱图像融合方法 - Google Patents

Abstract

本发明提供一种基于变分方法的高光谱图像融合方法，它是一种基于变分方法的高光谱图像和高分辨率可见光图像的融合方法，有三大步骤：步骤一：计算机在MATLAB7.1环境下读取高光谱图像数据；步骤二：建立基于变分方法的泛函表达式，选择合适的权重系数；步骤三：利用分裂Bregman迭代算法对上述泛函求解极值函数，得到融合结果。本发明克服了现有技术不宜直接移植到高光谱图像融合、光谱保持性不好的缺点，取得了较好的融合结果，在高光谱图像融合技术领域里具有广阔的应用前景。

Description

一种基于变分方法的高光谱图像融合方法

(一)技术领域：

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

(二)背景技术：

高光谱遥感是高光谱分辨率遥感(Hyperspectral Remote Sensing)的简称，它是在电磁波谱的可见光和红外波段范围内，获取许多光谱连续的影像数据的技术。利用成像光谱仪，可以采集到高光谱数据，所得高光谱数据包含了丰富的空间和光谱信息，即所谓“图谱合一”，其突出特点在于其极高的光谱分辨率(一般在10nm左右)，精细的光谱分辨率反映了地物光谱的细微特征，使得在光谱域内进行遥感定量分析成为可能。与传统遥感技术相比，高光谱遥感技术可以以较窄的波段区间、较多的波段数量提供遥感信息，在军事、地质、环境研究等领域中得到了广泛的应用。但是，由于高光谱数据一般是来自于空基或者天基的遥感平台，其空间分辨率一般在几米、几十米甚至上百米不等，空间分辨率较低。另一方面，利用其它的遥感手段，我们可以得到一些空间分辨率很高的图像，比如可见光图像。因此利用图像融合技术，在保持上述两种源图像优点的基础上，将两种图像结合，得到光谱分辨率和空间分辨率都比较高的结果图像，具有较高的价值。

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

(三)发明内容：

1、目的：本发明的目的是提供一种基于变分方法的高光谱图像融合方法，该方法利用光谱约束和图像细节约束来进行高光谱图像融合。

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

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

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

本发明中的高光谱数据为100谱段*256像素*200像素大小的三维数据立方体，其中100表示谱段数，256*200表示每一谱段图像大小，以下用H_n表示第n个谱段的高光谱源图像，用u_n表示第n个谱段融合的高光谱图像，用u表示一幅普通的图像；高分辨率数据大小为256像素*200像素，以下用M表示。

步骤二：建立泛函表达式。基于变分方法处理图像的基本框架，首先建立高光谱图像融合的基本泛函表达式，它包括三项：

第一项为图像的细节保持项E_a，该项用以实现对高光谱数据几何形状的保持。形式如下：

$E_a=\underset{n=1}{\overset{N}{\mathrm{\Σ}}}\left[\mathrm{\γ}\underset{\mathrm{\Ω}}{\∫}\right|\▿u_n|\mathrm{d\Ω}+\mathrm{\η}\underset{\mathrm{\Ω}}{\∫}\mathrm{div}\left(\mathrm{\θ}\right)*u_n\mathrm{d\Ω}]$

式中，δ，η为预先设定的常系数，u_n表示所求第n个谱段的融合图像，θ表示图像中每个像素点处的单位法向量。其定义如下：

表示对u求梯度，图像的水平和垂直方向分别用x和y表示，定义如下：

$\▿u=\left(\begin{array}{c}{u}_x\\ u_y\end{array}\right)$

div(θ)表示对θ求散度，具体求解时采用下式进行，

div(θ)＝θ_x+θ_y

以上两式中，u_x，u_y，θ_x，θ_y分别表示u对x，y的偏导数，θ对x，y的偏导数。对于图像而言，偏导数的计算可以用差分代替。

第二项为基于与源图相似性度量的保真项E_b，用以实现对高光谱数据灰度的保持，形式如下：

$E_b=\mathrm{\υ}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}\underset{\mathrm{\Ω}}{\∫}[{(u_n-H_n)}^2+\mathrm{\ρ}{(u_n-M)}^2]\mathrm{d\Ω}$

式中，υ，ρ为常系数，ρ用以平衡源数据中高光谱数据与高分辨率图像数据在融合结果中所占的不同比重。M表示高分辨率图像，H_n表示第n个谱段源高光谱数据。

第三项为光谱保持项E_c，该式用以实现对高光谱数据的光谱约束。形式如下：

$E_c=\mathrm{\&mu;}\underset{i,j=1i<j}{\overset{N}{\mathrm{\&Sigma;}}}\underset{\mathrm{\&Omega;}}{\&Integral;}\left[{(u_i*H_j-u_j*H_i)}^2\right]\mathrm{d\&Omega;}$

式中，μ为预先设定常数，u_i，u_j表示融合后第i，j谱段高光谱数据的图像，H_i，H_j表示原始高光谱数据第i，j两个谱段的图像。

将以上三项直接加和，建立的泛函表达式为：

$E=E_a+E_b+E_c$

$=\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}\left[\mathrm{\&gamma;}\underset{\mathrm{\&Omega;}}{\&Integral;}\right|\&dtri;u_n|\mathrm{d\&Omega;}+\mathrm{\&eta;}\underset{\mathrm{\&Omega;}}{\&Integral;}\mathrm{div}\left(\mathrm{\&theta;}\right)*u_n\mathrm{d\&Omega;}]+\mathrm{\&upsi;}\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}\underset{\mathrm{\&Omega;}}{\&Integral;}[{(u_n-H_n)}^2+\mathrm{\&rho;}{(u_n-M)}^2]\mathrm{d\&Omega;}$

$+\mathrm{\&mu;}\underset{i,j=1,i<j}{\overset{N}{\mathrm{\&Sigma;}}}\underset{\mathrm{\&Omega;}}{\&Integral;}\left[{(u_i*H_j-u_j*H_i)}^2\right]\mathrm{d\&Omega;}$

步骤三：对此能量泛函求极值函数，则得到的极值函数即为所求融合图像，本发明利用分裂Bregman迭代算法对此泛函求解极值函数。具体步骤如下：

(1)n＝1，开始循环

(2)参量初始化：

$u^0=H_i,{\left(d_n^0\right)}_x={\left(d_n^0\right)}_y={\left(b_n^0\right)}_x={\left(b_n^0\right)}_y=0$

(3)计算第k+1次迭代得到的u(即融合图像结果)：

$(2\mathrm{\&upsi;}+2\mathrm{\&upsi;}*\mathrm{\&rho;}+2\mathrm{\&mu;}\underset{j=1,j\&NotEqual;n}{\overset{N}{\mathrm{\&Sigma;}}}{\left(H_j\right)}^2-\mathrm{\&lambda;\&Delta;})u_n^{k+1}=2\mathrm{\&upsi;}*H_n+2\mathrm{\&upsi;}*\mathrm{\&rho;}*M-\mathrm{\&eta;}\mathrm{div}\left(\mathrm{\&theta;}\right)+2\mathrm{\&mu;}*H_n(\underset{j=1,j\&NotEqual;n}{\overset{N}{\mathrm{\&Sigma;}}}{u}_j^{k^{*}}*H_n)$

$+\mathrm{\&lambda;}\mathrm{div}(d_n^k-b_n^k)$

式中，λ为选定常数，参量d和b为迭代中产生的变量。这样取值，若j＜n，则取

若j＞n，则取

(4)进行参量d和b的更新：

${\left(d_n^{k+1}\right)}_x=\max (s^k-\frac{1}{\mathrm{\&lambda;}},0)*({\&dtri;}_x{u}_n^{k+1}+{\left(b_n^k\right)}_x)/s^k$

${\left(d_n^{k+1}\right)}_y=\max (s^k-\frac{1}{\mathrm{\&lambda;}},0)*({\&dtri;}_y{u}_n^{k+1}+{\left(b_n^k\right)}_y)/s^k$

${\left(b_n^{k+1}\right)}_x={\left(b_n^k\right)}_x+({\&dtri;}_x{u}_n^{k+1}-{\left(d_n^{k+1}\right)}_x)$

${\left(b_n^{k+1}\right)}_y={\left(b_n^k\right)}_y+({\&dtri;}_y{u}_n^{k+1}-{\left(d_n^{k+1}\right)}_y)$

其中 $s^k=\sqrt{{|{\&dtri;}_x{u}_n^{k+1}+{\left(b_n^k\right)}_x|}^2+{|{\&dtri;}_y{u}_n^{k+1}+{\left(b_n^k\right)}_y|}^2}$

(5)迭代停止条件。当融合结果量u在相邻两次迭代的过程中变化不大时停止迭代，本发明中停止条件为：

||u^k+1-u^k||＞tol

其中u^k+1表示当前迭代结果，u^k表示前一次迭代结果。本发明中，取tol＝0.4。得到第n个谱段的融合结果，存入新的数据立方体中。

(6)n＝n+1，继续第(2)步，直到n＞N，得到融合的高光谱数据立方体，至此完成对高光谱图像与高分辨率可见光图像的融合工作。

3、优点及功效：本发明的优点是：针对常见图像融合方法不宜直接移植到高光谱图像融合，而已有基于变分方法的融合手段建模复杂，求解繁琐，对图像配准要求高的特点，通过引入基于与源图相似性度量量来进行改进，通过对高光谱数据的光谱保持从光谱角度进行约束，以相似性度量对图像质量进行约束，最终得到了图像和光谱质量都比较好的融合结果。

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

(四)附图说明：

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

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

图2(b)本发明实验高光谱数据(50，128)点处融合前后的光谱曲线，细虚线为融合前光谱曲线，粗实线为融合后光谱曲线

(五)具体实施方式：

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

本发明在MATLAB7.1语言环境下实现。计算机读取高光谱遥感图像数据和配准好的高分辨率图像数据，首先基于与源图像的相似性度量建立正确的能量泛函表达式，选择合适的权重系数，然后对泛函表达式求解极值函数，极值函数即为所求的高光谱融合图像。在求解过程中，本发明延续已有方法中的分裂Bregman迭代算法。

本发明一种基于变分方法的高光谱图像融合方法，该方法的流程见图1所示，计算机配置采用：Intel dual-core处理器，主频1.6GHz，内存2GB。该融合方法包括以下步骤：

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

步骤二、建立基于变分方法的能量泛函。建立基于与源图像相似性度量的能量泛函，选择合适的权重系数。本发明实验中选择的参数如下：γ＝1，η＝0.5，v＝2，ρ＝4，μ＝2，λ＝1。

步骤三、对此泛函求解极值函数。利用分裂Bregman迭代求解。

(1)n＝1，开始迭代

(2)参量初始化：

$u^0=H_i,{\left(d_n^0\right)}_x={\left(d_n^0\right)}_y={\left(b_n^0\right)}_x={\left(b_n^0\right)}_y=0$

(3)计算第k+1次迭代得到的u(即融合图像结果)：

$(2\mathrm{\&upsi;}+2\mathrm{\&upsi;}*\mathrm{\&rho;}+2\mathrm{\&mu;}\underset{j=1,j\&NotEqual;n}{\overset{N}{\mathrm{\&Sigma;}}}{\left(H_j\right)}^2-\mathrm{\&lambda;\&Delta;})u_n^{k+1}=2\mathrm{\&upsi;}*H_n+2\mathrm{\&upsi;}*\mathrm{\&rho;}*M-\mathrm{\&eta;}\mathrm{div}\left(\mathrm{\&theta;}\right)$

$+2\mathrm{\&mu;}*H_n\left(\underset{j=1,j\&NotEqual;n}{\overset{N}{\mathrm{\&Sigma;}}}{u}_j^{k^{*}}{*H}_n\right)+\mathrm{\&lambda;}\mathrm{div}(d_n^k-b_n^k)$

(4)进行参量d和b的更新：

${\left(d_n^{k+1}\right)}_x=\max (s^k-\frac{1}{\mathrm{\&lambda;}},0)*({\&dtri;}_x{u}_n^{k+1}+{\left(b_n^k\right)}_x)/s^k$

${\left(d_n^{k+1}\right)}_y=\max (s^k-\frac{1}{\mathrm{\&lambda;}},0)*({\&dtri;}_y{u}_n^{k+1}+{\left(b_n^k\right)}_y)/s^k$

${\left(b_n^{k+1}\right)}_x={\left(b_n^k\right)}_x+({\&dtri;}_x{u}_n^{k+1}-{\left(d_n^{k+1}\right)}_x)$

${\left(b_n^{k+1}\right)}_y={\left(b_n^k\right)}_y+({\&dtri;}_y{u}_n^{k+1}-{\left(d_n^{k+1}\right)}_y)$

其中 $s^k=\sqrt{{|{\&dtri;}_x{u}_n^{k+1}+{\left(b_n^k\right)}_x|}^2+{|{\&dtri;}_y{u}_n^{k+1}+{\left(b_n^k\right)}_y|}^2}$

(5)迭代停止条件。当融合结果量u在相邻两次迭代的过程中变化不大时停止迭代，其停止条件为：

||u^k+1-u^k||＞tol

其中u^k+1表示当前迭代结果，u^k表示前一次迭代结果。本发明中，取tol＝0.4。得到第n个谱段的融合结果，存入新的数据立方体中。

(6)n＝n+1，继续第(2)步，直到n＞N，得到融合的高光谱数据，输出高光谱高分辨率数据立方体，至此已完成高光谱图像和高分辨率图像的融合工作，融合结果可以用于后续的目标检测等其它图像处理工作。

有益效果：

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

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

Claims (1)

1.一种基于变分方法的高光谱图像融合方法，其特征在于：将高光谱图像融合的问题转化为对泛函求解极值函数的优化问题，同时考虑了光谱约束和图像的细节约束，该方法具体步骤如下：

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

其中的高光谱数据为100谱段*256像素*200像素大小的三维数据立方体，100表示谱段数，256*200表示每一谱段图像大小，以下用H_n表示第n个谱段的高光谱源图像，用u_n表示第n个谱段融合的高光谱图像，用u表示一幅普通的图像；高分辨率数据大小为256像素*200像素，以下用M表示；

步骤二：建立泛函表达式；基于变分方法处理图像的基本框架，首先建立高光谱图像融合的基本泛函表达式，它包括三项：

第一项为图像的细节保持项E_a，该项用以实现对高光谱数据几何形状的保持，形式如下：

$E_a=\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}\left[\mathrm{\&gamma;}\underset{\mathrm{\&Omega;}}{\&Integral;}\right|\&dtri;u_n|\mathrm{d\&Omega;}+\mathrm{\&eta;}\underset{\mathrm{\&Omega;}}{\&Integral;}\mathrm{div}\left(\mathrm{\&theta;}\right)*u_n\mathrm{d\&Omega;}]$

式中，δ，η为预先设定的常系数，u_n表示所求第η个谱段的融合图像，θ表示图像中每个像素点处的单位法向量。其定义如下：

表示对u求梯度，图像的水平和垂直方向分别用x和y表示，定义如下：

$\&dtri;u=\left(\begin{array}{c}{u}_x\\ u_y\end{array}\right)$

div(θ)表示对θ求散度，具体求解时采用下式进行，

div(θ)＝θ_x+θ_y

以上两式中，u_x，u_y，θ_x，θ_y分别表示u对x，y的偏导数，θ对x，y的偏导数。对于图像而言，偏导数的计算可以用差分代替。

第二项为基于与源图相似性度量的保真项E_b，用以实现对高光谱数据灰度的保持，形式如下：

$E_b=\mathrm{\&upsi;}\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}\underset{\mathrm{\&Omega;}}{\&Integral;}[{(u_n-H_n)}^2+\mathrm{\&rho;}{(u_n-M)}^2]\mathrm{d\&Omega;}$

式中，υ，ρ为常系数，ρ用以平衡源数据中高光谱数据与高分辨率图像数据在融合结果中所占的不同比重。M表示高分辨率图像，H_n表示第n个谱段源高光谱数据。

第三项为光谱保持项E_c，该式用以实现对高光谱数据的光谱约束。形式如下：

$E_c=\mathrm{\&mu;}\underset{i,j=1i<j}{\overset{N}{\mathrm{\&Sigma;}}}\underset{\mathrm{\&Omega;}}{\&Integral;}\left[{(u_i*H_j-u_j*H_i)}^2\right]\mathrm{d\&Omega;}$

式中，μ为预先设定常数，u_i，u_j表示融合后第i，j谱段高光谱数据的图像，H_i，H_j表示原始高光谱数据第i，j两个谱段的图像。

将以上三项直接加和，建立的泛函表达式为：

$E=E_a+E_b+E_c$

$=\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}\left[\mathrm{\&gamma;}\underset{\mathrm{\&Omega;}}{\&Integral;}\right|\&dtri;u_n|\mathrm{d\&Omega;}+\mathrm{\&eta;}\underset{\mathrm{\&Omega;}}{\&Integral;}\mathrm{div}\left(\mathrm{\&theta;}\right)*u_n\mathrm{d\&Omega;}]+\mathrm{\&upsi;}\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}\underset{\mathrm{\&Omega;}}{\&Integral;}[{(u_n-H_n)}^2+\mathrm{\&rho;}{(u_n-M)}^2]\mathrm{d\&Omega;}$

$+\mathrm{\&mu;}\underset{i,j=1,i<j}{\overset{N}{\mathrm{\&Sigma;}}}\underset{\mathrm{\&Omega;}}{\&Integral;}\left[{(u_i*H_j-u_j*H_i)}^2\right]\mathrm{d\&Omega;}$

步骤三：对此能量泛函求极值函数，则得到的极值函数即为所求融合图像，本发明利用分裂Bregman迭代算法对此泛函求解极值函数。具体步骤如下：

(1)n＝1，开始循环

(2)参量初始化：

$u^0=H_i,{\left(d_n^0\right)}_x={\left(d_n^0\right)}_y={\left(b_n^0\right)}_x={\left(b_n^0\right)}_y=0$

(3)计算第k+1次迭代得到的u(即融合图像结果)：

$(2\mathrm{\&upsi;}+2\mathrm{\&upsi;}*\mathrm{\&rho;}+2\mathrm{\&mu;}\underset{j=1,j\&NotEqual;n}{\overset{N}{\mathrm{\&Sigma;}}}{\left(H_j\right)}^2-\mathrm{\&lambda;\&Delta;})u_n^{k+1}=2\mathrm{\&upsi;}*H_n+2\mathrm{\&upsi;}*\mathrm{\&rho;}*M-\mathrm{\&eta;}\mathrm{div}\left(\mathrm{\&theta;}\right)$

$+2\mathrm{\&mu;}*H_n\left(\underset{j=1,j\&NotEqual;n}{\overset{N}{\mathrm{\&Sigma;}}}{u}_j^{k^{*}}{*H}_n\right)+\mathrm{\&lambda;}\mathrm{div}(d_n^k-b_n^k)$

式中，λ为选定常数，参量d和b为迭代中产生的变量。

这样取值，若j＜n，则取

若j＞n，则取

(4)进行参量d和b的更新：

${\left(d_n^{k+1}\right)}_x=\max (s^k-\frac{1}{\mathrm{\&lambda;}},0)*({\&dtri;}_x{u}_n^{k+1}+{\left(b_n^k\right)}_x)/s^k$

${\left(d_n^{k+1}\right)}_y=\max (s^k-\frac{1}{\mathrm{\&lambda;}},0)*({\&dtri;}_y{u}_n^{k+1}+{\left(b_n^k\right)}_y)/s^k$

${\left(b_n^{k+1}\right)}_x={\left(b_n^k\right)}_x+({\&dtri;}_x{u}_n^{k+1}-{\left(d_n^{k+1}\right)}_x)$

${\left(b_n^{k+1}\right)}_y={\left(b_n^k\right)}_y+({\&dtri;}_y{u}_n^{k+1}-{\left(d_n^{k+1}\right)}_y)$

其中 $s^k=\sqrt{{|{\&dtri;}_x{u}_n^{k+1}+{\left(b_n^k\right)}_x|}^2+{|{\&dtri;}_y{u}_n^{k+1}+{\left(b_n^k\right)}_y|}^2}$

(5)迭代停止条件。当融合结果量u在相邻两次迭代的过程中变化不大时停止迭代，本发明中停止条件为：

||u^k+1-u^k‖＞tol

其中u^k+1表示当前迭代结果，u^k表示前一次迭代结果。本发明中，取tol＝0.4。得到第n个谱段的融合结果，存入新的数据立方体中。

(6)n＝n+1，继续第(2)步，直到n＞N，得到融合的高光谱数据立方体，至此完成对高光谱图像与高分辨率可见光图像的融合工作。

Images