CN102609929A - 一种自适应独立信息遥感影像融合方法 - Google Patents

Abstract

本发明涉及一种独立信息源遥感影像融合方法，该方法通过独立成分分析建立多光谱影像的独立信息源模型，然后运用全色波段影像向量替换主体独立成分以获得融合后影像。在获得多光谱影像主体独立信息源过程中，首先制定多光谱影像独立信息源分离标准，然后运用小波变换和独立成分分析进行迭代求解，直至混合矩阵满足独立信息源分离标准，在此基础上根据混合矩阵以及主体独立成分来获得主体独立信息源。本发明可以在大幅度提高原始多光谱影像空间分辨率的同时有效地保持原有的光谱信息，且能够在二者之间取得一个很好的平衡；此外不需要设置参数以及人工的干预，因此，本发明具有保真度好、自动化程度高等特点。

Description

一种自适应独立信息遥感影像融合方法

技术领域

本发明涉及一种遥感影像处理方法，特别是一种自适应独立信息遥感影像融合方法。具体来说是涉及遥感影像全色波段和多光谱波段的融合，用于实现遥感影像全色波段和多光谱波段的自动融合。

背景技术

随着遥感技术的发展，同一地区，不同传感器的数据被大量获得。这些数据具有不同的特点，例如全色波段影像的空间分辨率高但光谱分辨率低，而多光谱影像具有丰富的光谱信息(具有多个波段)，但其空间分辨率较低。而在实际的应用过程中，往往需要遥感数据同时具有高的空间分辨率和丰富的波谱信息。因此，采用影像融合技术充分利用多源影像中信息的互补性，使融合结果在提高空间分辨率的同时，保持原始影像的波谱信息，具有重要的研究意义和应用价值。

迄今为止，人们已经提出了一系列的遥感影像融合方法，比如色调-饱和度-亮度变换法(Hue-saturation-intensity，HSV)，BT(Brovey transform，)，主成分分析法(principlecomponent analysi s，PCA)和小波变换法(wavelet)等。虽然HSV，BT以及PCA能够使融合后影像具有较高的空间分辨率，但是相比原始的多光谱影像，它们严重地扭曲了光谱信息。小波变换法可以通过多尺度分解来调整融合后影像的空间信息和光谱信息，但是该方法的一种严重不足是它不能同时使融合后的影像同时具有高的空间分辨率和丰富的波谱信息，此外，针对不同的影像，其全色波段和多光谱波段的分辨率之比不同，比如ETM+影像其全色波段和多光谱波段的空间分辨率之比是1∶2，而IKONOS影像的全色波段和多光谱波段空间分辨率之比是1∶4，因此该方法需要人工的干预来确定最优分解级数。所述已有遥感影像融合方法存在着不能同时使融合后影像具有高的空间分辨率和丰富波谱信息的不足，以及部分融合方法需要手工设置或选取参数，不能实现自动化的弊端

发明内容

本发明目的在于提供一种自适应独立信息遥感影像融合方法，它能够克服已有技术的不足，公开一种利用子带分解和独立成分分析来实现遥感影像全色波段与多光谱波段融合方法。该方法不仅能够使融合后影像同时具有高的空间分辨率和丰富的波谱信息，而且完全做到自动化处理。

本发明目的通过如下技术方案实现。包括如下步骤：

步骤I.影像预处理

影像预处理阶段包括对影像进行几何精矫正以及图像向量化。通过几何精矫正使全色波段影像和多光谱波段影像统一在同一坐标系下，并且对多光谱波段影像进行重采样，使其具有与全色波段影像相等的像元大小。对任意波段影像按行为顺序构成一个向量，分别形成多光谱影像向量组和全色波段影像向量。

步骤II.构建多光谱假彩色合成影像独立信息源模型，三个独立信息源分别代表彩色影像的主体独立成分、色彩独立成分以及空间细节独立成分信息，并制定分离这些独立信息源所应满足的标准。

步骤III.利用子带分解和独立成分分析寻找影像向量的独立信息源。

(1)对多光谱影像向量组进行独立成分分析初步寻找影像的独立信息源，独立成分分析定义如下：

x(t)＝As(t)                                (1)

其中A是混合矩阵，判断A是否满足独立信息源标准。如果满足，转入步骤IV，否则转入本步骤的第(2)步.

(2)对多光谱影像向量组或小波结点向量组进行小波分解并寻找其互信息量最大的向量组，其中互信息定义如下：

$I(x_1,x_2,...,x_n)\&ap;\frac{1}{4}\underset{\underset{n\&NotEqual;l}{0\&le;n<l\&le;N}}{\mathrm{\&Sigma;}}{\mathrm{cum}}^2(x_n,x_l)+\frac{1}{12}\underset{\underset{n\&NotEqual;l}{0\&le;n<l\&le;N}}{\mathrm{\&Sigma;}}\left({\mathrm{cum}}^2(x_n,x_n,x_l){+\mathrm{cum}}^2(x_n,x_l,x_l)\right)+$

$\frac{1}{48}\underset{\underset{n\&NotEqual;l}{0\&le;n<l\&le;N}}{\mathrm{\&Sigma;}}({\mathrm{cum}}^2(x_n,x_n,x_n,x_l)+{\mathrm{cum}}^2(x_n,x_n,x_l,x_l)+{\mathrm{cum}}^2(x_n,x_l,x_l,x_l))---\left(2\right)$

(3)对互信息量最大的子带向量组重新进行独立成分分析，生成混合矩阵A和其逆矩阵W，然后判断其混合矩阵A是否满足信息源独立性标准。如果满足标准，转入步骤IV，否则转入本步骤的第(2)步.

步骤IV.获得多光谱影像向量X和矩阵W，计算多光谱影像的独立信息源s(t)＝Wx(t)。接着根据混合矩阵A来确定独立信息源中主体独立成分(ma in body independent component，MBIC)。

步骤V.以主体独立成分为基础，对全色波段影像向量进行直方图匹配：

PAN_new＝H(PAN，MBIC)                        (3)

其中H是直方图匹配函数，PANnew是以MBIC为基础对PAN进行直方图匹配的结果。并用其来替换MBIC，生成新的独立信息源s_new(t)。

步骤VI.利用s_new(t)和混合矩阵A进行逆独立成分变换生成融合影像。

IMG_fusion＝IICA(s_new(t)，A)

式中IMGfusion是融合结果，IICA是逆独立成分变换。

本发明与现有技术相比不仅使融合后影像能够大幅度提高空间分辨率，很好地继承全色波段影像的空间细节信息，而且能够保持多光谱波段的光谱信息。使融合后影像能够在两者之间达到一个很好的平衡。并且在处理过程中不需要人工的干预，例如设置参数等，能够做到完全的自动化。主要体现在：首先，本发明建立了构建多光谱假彩色合成影像独立信息源模型，在这个模型中，多光谱影像的三个独立成分一主体独立成分、彩色独立成分、细节独立成分分别表示多光谱影像的空间主体、彩色以及空间细节信息。其次，为了解决现实中源信号不满足统计独立性原则，采用迭代的小波分析的方法来寻找具有统计独立性的原始信号的子带。再次，通过混合矩阵A来构建独立信息源标准，以致用最小的代价自适应地寻找这些独立子带。

附图说明

图1本发明的流程图

图2IKONOS数据实验结果

图3ETM+数据试验结果

具体实施方式

下面详细描述本发明的具体实施方式

1.系统流程

步骤I.影像预处理

首先对多光谱影像和全色波段影像进行几何精矫正。分别对ETM+和IKONOS影像的多光谱波段(2，3和4)和其对应全色波段进行配准，并对多光谱影像进行重采样，使其像元大小和全色波段一样，并确保它们的配准精度在0.5个像元以内。

然后对配准后任意波段影像以行为顺序连接成一个向量，多光谱波段构成向量组x(t)＝[x_r(t)，x_g(t)，x_b(t)]，全色波段构成向量PAN(t)。

步骤II.制定影像独立信息源标准

当混合矩阵中每一行绝对值的最大值处在同一列时，表明独立信息源被分离，这样三个独立成分分别代表多光谱影像的空间主体信息、彩色信息以及空间细节信息。否则，表明独立信息源没有完全分离。

步骤III.利用子带分解和独立成分分析寻找影像向量的独立信息源

(1)对多光谱向量组进行独立成分变换，以获得混合矩阵A₀和W₀。

x(t)＝A₀s(t)

A₀W₀＝I

其中x(t)是多光谱影像向量组，A₀是混合矩阵，s(t)是隐变量既独立成分。利用步骤II制定的信息源独立性分解标准来判断多光谱影像的独立信息源被分离出来，如果完全分解，进入步骤IV，否则转向本步骤的第(2)步。

(2)利用′db1′小波对原始多光谱向量组或者子节点向量组进行下一级分解，获得每个向量的两个子带，其分别为高频子带和低频子带。

$\left\{\begin{array}{c}{x}_{k+1,r}^j\left(t\right)=W\left[x_{k,r}\left(t\right)\right]\\ x_{k+1,g}^j\left(t\right)=W\left[x_{k,g}\left(t\right)\right]\\ x_{k+1,b}^j\left(t\right)=W\left[x_{k,b}\left(t\right)\right]\end{array}\right.$

式中，x_k，r(t)，x_k，g(t)和x_k，b(t)分别表示多光谱影像向量组第k级具有最小互信息量的向量组，k取大于等于0的整数，k＝0时表示原始多光谱影像向量组；W[·]表示对向量组进行一级小波分解；和

分别为对x_k，r(t)，x_k，g(t)和x_k，b(t)进行一级小波分解的子带，j的取值为1和2，分别为低频和高频子带。

(3)根据j的取值，把

和划分为和

两组向量组，分别计算每一组的互信息量。并选择最小互信息量的一组。

$x_{k+1}\left(t\right)=\mathrm{Max}\left(I\left(x_{k+1}^j\left(t\right)\right)\right)$

其中

$I\left(x_{k+1}^j\left(t\right)\right)=I(x_{k+1,r}^j\left(t\right),x_{k+1,g}^j\left(t\right),x_{k+1,b}^j\left(t\right))$

$=\frac{1}{4}\left({\mathrm{cum}}^2(x_{k+1,r}^j,S_{k+1,g}^j)+{\mathrm{cum}}^2(x_{k+1,r}^j,S_{k+1,b}^j)+{\mathrm{cum}}^2(x_{k+1,g},S_{k+1,b}^j)\right)+$

$\frac{1}{12}\left({\mathrm{cum}}^2(x_{k+1,r}^j,x_{k+1,r}^j{,x}_{k+1,g}^j)+{\mathrm{cum}}^2(x_{k+1,r}^j{,x}_{k+1,g}^j,x_{k+1,g}^j)+{\mathrm{cum}}^2(x_{k+1,r}^j,x_{k+1,r}^j,x_{k+1,b}^j\right)+$

${\mathrm{cum}}^2(x_{k+1,r}^j,x_{k+1,b}^j,x_{k+1,b}^j)+{\mathrm{cum}}^2(x_{k+1,g}^j,x_{k+1,g}^j,x_{k+1,b}^j)+{\mathrm{cum}}^2(x_{k+1,g}^j,x_{k+1,b}^j,x_{k+1,b}^j))+$

$\frac{1}{48}({\mathrm{cum}}^2(x_{k+1,r}^j\left(t\right),x_{k+1,r}^j\left(t\right),x_{k+1,r}^j\left(t\right),x_{k+1,g}^j\left(t\right))+{\mathrm{cum}}^2(x_{k+1,r}^j\left(t\right),x_{k+1,r}^j\left(t\right),x_{k+1,g}^j\left(t\right),x_{k+1,g}^j\left(t\right))+$

${\mathrm{cum}}^2(x_{k+1,r}^j\left(t\right),x_{k+1,g}^j\left(t\right),x_{k+1,g}^j\left(t\right),x_{k+1,g}^j\left(t\right))+{\mathrm{cum}}^2(x_{k+1,r}^j\left(t\right),x_{k+1,r}^j\left(t\right),x_{k+1,r}^j\left(t\right),x_{k+1,b}^j\left(t\right))+$

${\mathrm{cum}}^2(x_{k+1,r}^j\left(t\right),x_{k+1,r}^j\left(t\right),x_{k+1,b}^j\left(t\right),x_{k+1,b}^j\left(t\right))+{\mathrm{cum}}^2(x_{k+1,r}^j\left(t\right),x_{k+1,b}^j\left(t\right),x_{k+1,b}^j\left(t\right),x_{k+1,b}^j\left(t\right))+$

${\mathrm{cum}}^2(x_{k+1,g}^j\left(t\right),x_{k+1,g}^j\left(t\right),x_{k+1,g}^j\left(t\right),x_{k+1,b}^j\left(t\right))+{\mathrm{cum}}^2(x_{k+1,g}^j\left(t\right),x_{k+1,g}^j\left(t\right),x_{k+1,b}^j\left(t\right),x_{k+1,b}^j\left(t\right))+$

${\mathrm{cum}}^2(x_{k+1,g}^j\left(t\right),x_{k+1,b}^j\left(t\right),x_{k+1,b}^j\left(t\right),x_{k+1,b}^j\left(t\right)))$

式中cum是互累积量函数，cum(x₁，x₂)，cum(x₁，x₂，x₃)和cum(x₁，x₂，x₃，x₄)分别为2阶、3阶以及4阶互累积量。

(4)对互信息量最小的向量组x_k+1(t)进行独立成分分析，计算其混合矩阵A_k+1和它的逆矩阵W_k+1.

x_k+1(t)＝A_k+1s_k+1(t)

A_k+1W_k+1＝I

式中A_k+1为混合矩阵，s_k+1(t)为子带向量组x_k+1(t)的独立信息源分量。然后判断其混合矩阵A_k+1是否满足信息源独立标准。如果满足标准，转入步骤IV，否则转向本步骤的第(2)步。

步骤IV.寻找主体独立成分(MBIC)

根据矩阵W_l和原始多光谱影像向量组x(t)，计算多光谱影像向量组的独立信息源分量：

s(t)＝W_lx(t)

式中，s(t)为多光谱影像的独立信息源分量，W_l和混合矩阵A_l的逆；1取值分两种情况，如果没有进行小波分解，那么1＝0，这时A_l是原始多光谱影像的混合矩阵；否则1＝k+1，这时A_l是k+1级最小互信息量子带向量组的混合矩阵。x(t)是原始多光谱影像向量组。其次，根据混合矩阵A_l选择主体独立分量。首先对矩阵求绝对值，然后确定三行中绝对值最大的一列，该列所对应独立成分的行即为主体独立成分。

步骤V.用全色波段向量替换主体独立成分

以主体独立分量为基础，对全色波段向量进行直方图匹配。首先计算主体独立成分向量和全色波段向量的直方图均衡化函数：

$m=T_{\mathrm{PAN}}\left(r\right)={\&Integral;}_0^r{p}_r\left(w\right)\mathrm{dw}$

$n=G_{\mathrm{MBIC}}\left(z\right)={\&Integral;}_0^z{p}_z\left(w\right)\mathrm{dw}$

式中T_PAN(r)和G_MBIC(z)分别为主体独立成分和全色波段向量的直方图均衡化函数，m和n是r和z经过直方图均衡化后的输出。在此基础上，全色波段的直方图匹配函数为：

$z=G_{\mathrm{MBIC}}^{-1}\left(T_{\mathrm{PAN}}\left(r\right)\right)$

接着，用直方图匹配后的全色波段向量替换独立成分中的主体独立成分，生成新的独立分量组s_new(t)。

s_new(t)′＝[PAN_new，IC_color，IC_detail]

步骤VI.影像融合

运用混合矩阵A_l和替换后的独立成分s_new(t)，进行逆独立成分变换以获得融合影像。

IMG_fusion＝A_ls_new(t)

式中IMGfus ion是融合后影像；A_l是混合矩阵，1含义与步骤IV相同；s_new(t)是替换后新的独立成分向量组。

2.实验结果分析及本发明的优越性

为了检验本发明的性能，两种数据IKONOS和ETM+影像，其全色和多光谱波段空间分辨率比率分别为1∶4和1∶2被用来验证本发明的有效性。此外，选择基于HSV变换，PCA变换以及小波变换方法作为对比。

实验一：IKONOS影像全色波段(空间分辨率为1米)与多光谱影像(波段4，3和2，空间分辨率为4米)，经过配准后影像大小为400×400。

实验二：ETM+影像全色波段(空间分辨率为15米)与多光谱影像(波段4，3和2，空间分辨率为30米)，经过配准后影像大小为400×400。

融合影像质量评价：对于融合影像的质量评价包括主观视角和客观指标两方面，好的融合结果不仅能够保持原始多光谱影像的光谱信息，而且能够同时较好地保持全色波段的空间信息。本发明主观视觉评价为目视感受；相比RGB彩色空间，HSI彩色空间更容易描述颜色和空间信息，因此客观评价在HSI彩色空间中进行，采取的指标包括均值，标准差和相关系数。

(1)均值和标准差：反映融合后影像与原始影像能量的接近程度。其与原始影像的值越相近，表示融合后影像从原始多光谱影像得到的信息量越大。

(2)相关系数：衡量融合后影像与原始影像的相似度，相关系数分为两种，一种用来衡量融合后影像保持原始多光谱影像彩色信息的量，其通过HSI彩色空间中融合后影像和多光谱影像H和S的相关系数来衡量。H和S的相关系数越接近1，表明融合后结果保持原始多光谱影像彩色信息的程度越高；另外一种用来融合后影像保持全色波段空间信息的量，通过HSI彩色空间中融合后影像的I分量与全色波段影像的相关系数来衡量。其值越接近1，表明融合后影像保持全色波段空间信息的程度越高。

主观评价：

从视觉上看，虽然HSV和PCA法的融合结果能够较大程度地提高原始多光谱影像的空间分辨率，但是其对多光谱影像的彩色信息扭曲严重，不符合人眼的视觉感受。基于小波的融合方法随着小波分解级数的增加，融合后影像的空间分辨率呈提高态势，但是其对多光谱影像光谱信息的保持呈下降态势，由于两种影像全色和多光谱影像空间分辨率之比不同，因此对于IKONOS数据，小波分解级数为三级时方可达到较好的融合结果，而ETM+影像其小波分解级数为二级时即可达到较好的结果。相比以上方法，本发明的融合结果不仅能够使融合后影像同时具有高的空间分辨率和丰富的波谱信息，并且使两者之间能够达到很好的平衡。

客观评价：

两个实验的客观评价结果见表1-6。

表1融合结果在HSI空间中各分量的均值和标准差(IKONOS数据)

表2融合结果与多光谱影像在HSI彩色空间的相关系数(IKONOS数据)

表3融合结果在HSI空间中I分量与全色波段的相关系数(IKONOS数据)

表4融合结果在HSI空间中各分量的均值和标准差(ETM+数据)

表5融合结果与多光谱影像在HSI彩色空间的相关系数(ETM+数据)

表6融合结果在HSI空间中I分量与全色波段的相关系数(ETM+数据)

由实验一(表1-3)的实验数据可知，相比其他方法，HSV法融合结果I分量的均值和标准差均严重偏离原始的多光谱影像，因此，即使在HSI空间中其与原始多光谱影像的H和S分量具有较高的相关系数，该方法仍然严重扭曲原始多光谱影像的光谱；而其余方法均值和方差与原始多光谱影像偏离较小。相比其他方法，本发明融合结果的I分量与全色波段具有最高的相关系数，因此，该方法的融合结果能够最大程度地保持全色波段的空间信息。此外，本发明融合方法与原始多光谱影像H和S分量的相关系数大于HSV以外的其他方法，表明该方法很好地保持了原始多光谱影像的光谱信息。

由实验二(表4-6)的实验数据可知，HSV和PCA均具有验证的光谱扭曲。小波方法随着分解级数的增加，其融合结果空间分辨率增加，但是其光谱信息损失严重，小波二级分解融合结果H分量与原始多光谱影像H分量的相关系数只有0.54。而本发明融合结果与多光谱影像H和S分量的相关系数为0.98和0.78，且I分量与全色波段的相关系数为0.99，因此由实验一和实验二可知，本发明融合结果能够同时具体有高的空间分辨率和丰富的波谱信息

Claims (1)

1.一种自适应独立信息遥感影像融合方法，利用子带分解和独立成分分析来实现遥感影像全色波段和多光谱波段融合，包括如下步骤：

步骤I，影像预处理

影像预处理阶段包括对影像进行几何精矫正以及图像向量化。通过几何精矫正使全色波段影像和多光谱波段影像统一在同一坐标系下，并且对多光谱波段影像进行重采样，使其具有与全色波段影像相等的像元大小，对任意波段影像按行为顺序构成一个向量，分别形成多光谱影像向量组和全色波段影像向量；

步骤II，构建多光谱假彩色合成影像独立信息源模型，三个独立信息源分别代表彩色影像的主体独立成分、色彩独立成分以及空间细节独立成分信息，并制定分离这些独立信息源所应满足的标准；

步骤III，利用子带分解和独立成分分析寻找影像向量的独立信息源，

(1)对多光谱影像向量组进行独立成分分析初步寻找影像的独立信息源，独立成分分析定义如下：

x(t)＝As(t)                                    (1)

其中A是混合矩阵，判断A是否满足独立信息源标准，如果满足，转入步骤IV，否则转入本步骤的第(2)步；

(2)对多光谱影像向量组或小波结点向量组进行小波分解并寻找其互信息量最大的向量组，其中互信息定义如下：

$I(x_1,x_2,...,x_n)\&ap;\frac{1}{4}\underset{\underset{n\&NotEqual;l}{0\&le;n<l\&le;N}}{\mathrm{\&Sigma;}}{\mathrm{cum}}^2(x_n,x_l)+\frac{1}{12}\underset{\underset{n\&NotEqual;l}{0\&le;n<l\&le;N}}{\mathrm{\&Sigma;}}\left({\mathrm{cum}}^2(x_n,x_n,x_l){+\mathrm{cum}}^2(x_n,x_l,x_l)\right)+$

$\frac{1}{48}\underset{\underset{n\&NotEqual;l}{0\&le;n<l\&le;N}}{\mathrm{\&Sigma;}}({\mathrm{cum}}^2(x_n,x_n,x_n,x_l)+{\mathrm{cum}}^2(x_n,x_n,x_l,x_l)+{\mathrm{cum}}^2(x_n,x_l,x_l,x_l))---\left(2\right)$

(3)对互信息量最大的子带向量组重新进行独立成分分析，生成混合矩阵A和其逆矩阵W，然后判断其混合矩阵A是否满足信息源独立性标准，如果满足标准，转入步骤IV，否则转入本步骤的第(2)步；

步骤IV，获得多光谱影像向量X和矩阵W，计算多光谱影像的独立信息源s(t)＝Wx(t)，接着根据混合矩阵A来确定独立信息源中主体独立成分(ma in body independent component，MBIC)；

步骤V，以主体独立成分为基础，对全色波段影像向量进行直方图匹配：

PAN_new＝H(PAN，MBIC)                            (3)

其中H是直方图匹配函数，PANnew是以MBIC为基础对PAN进行直方图匹配的结果，并用其来替换MBIC，生成新的独立信息源s_new(t)；

步骤VI，利用s_new(t)和混合矩阵A进行逆独立成分变换生成融合影像，

IMG_fusion＝IICA(s_new(t)，A)

式中IMGfusion是融合结果，IICA是逆独立成分变换。

Images