CN107346535A - 基于多尺度几何分析和pca变换的高光谱遥感图像融合方法 - Google Patents

Abstract

本发明涉及一种基于多尺度几何分析和PCA变换的高光谱遥感图像融合方法，其特征在于，包括以下步骤：对N个波段的高光谱图像进行预处理、配准；分别对N个波段的高光谱图像进行多尺度几何分析变换，获取每一个波段高光谱图像的多尺度几何分析系数C_j(i)和几何流G_j(i)；用PCA变换分析N个波段高光谱图像的多尺度几何分析系数和几何流，获得各自的主成分；进行多尺度几何分析逆变换，重构得到融合新图像。通过本发明的方法获取的遥感图像质量更好，也将会更有利于遥感图像的提取、分类等后续处理。

Description

基于多尺度几何分析和PCA变换的高光谱遥感图像融合方法

技术领域

本发明涉及一种高光谱遥感图像的融合方法。

背景技术

遥感是利用任何物体都有不同的电磁波反射或辐射特征，通过获取地面目标反射或辐射出来的电磁波来得到目标信息，并对所获得的信息进行处理，实现目标的定位、识别、定性或定量的描述等。遥感科学是一门综合性科学，它以空间信息科学、电子信息科学、光学、计算机科学、生物、化学等学科为基础，通过对空间对象的观测，为国民经济和科学研究等人类活动提供了丰富的海量遥感影像数据。

多尺度几何分析变换是基于边缘的一种二维小波变换，通过把图像中的几何特征定义为矢量场，而不是简单的看作普通边缘的集，能够自适应地跟踪图像的几何正则，正则表达式的约束，到达自适应的稀疏表达图像，被称之为自适应的多尺度几何分析工具。

PCA变换是通过使用Principal Componentsa Analysis选项生成互不相关的输出波段，达到隔离噪声和减少数据集的维数的方法。

高光谱图像超多的波段数目形成了海量的数据，处理高光谱图像就面临着维数过多的问题：波段选择、运算效率低等困难。

发明内容

本发明的目的是将多尺度几何分析和PCA变换结合起来融合高光谱遥感图像。

为了达到上述目的，本发明的技术方案是提供了一种基于多尺度几何分析和 PCA变换的高光谱遥感图像融合方法，其特征在于，包括以下步骤：

步骤1、对N个波段的高光谱图像进行预处理、配准；

步骤2、分别对N个波段的高光谱图像进行多尺度几何分析变换，获取每一个波段高光谱图像的多尺度几何分析系数C_j(i)和几何流G_j(i)，式中，j表示融合的源图像的个数，i表示第j幅图像中分割的块数；

步骤3、用PCA变换分析N个波段高光谱图像的多尺度几何分析系数和几何流，获得各自的主成分，分别记为CF(i)和GF(i)；

步骤4、将CF(i)作为多尺度几何分析系数，将GF(i)作为几何流，进行多尺度几何分析逆变换，重构得到融合新图像。

优选地，所述步骤2包括：

步骤2.1、对原始的高光谱图像f(x,y)进行二维正交小波变换；

步骤2.2、对步骤2.1得到的图像进行二进四叉树分割，直到最底层块的分解尺度达到预先设定的最小尺度，获得图像的多尺度几何分析块，分割时，采用借助目标函数的最小化准则对图像各高频子带最佳分割块大小进行优化，达到最优的分割；

步骤2.3、在多尺度几何分析块内选择最佳几何方向，即几何流G_j(i)，同时计算多尺度几何分析系数。

优选地，步骤2.2中所述优化包括以下步骤：

步骤2.2.1、将L初始化为4，对每个尺寸为L×L的方块S计算最佳方向和相应的Lagrange函数值，将多尺度几何分析块L_0(S)初始化为本步骤计算得到的Lagrange函数值；

步骤2.2.2、令L＝2L，对每个尺寸为L×L的方块S计算最佳几何流方向和相应的Lagrange函数值L(S)；

步骤2.2.3、对每个尺寸为L×L的方块S，分割为四个孩子，分别记为 S1,S2,S3,S4，这4个孩子作为叶节点联合在一起的Lagrange函数值为：

式中，L_0(S、1)L_0S(、2)L_S、0(3L)S_分0别(表4)示四个子块的Lagrange值，λ表示Lagrange乘子，T表示量化阈值；

步骤2.2.4、令

步骤2.2.5、若L＝2∧(j/2)，j表示分割块的大小，返回步骤2.2.2，否则输出L_0(S)作为四叉树分割的结果。

优选地，步骤2.3中，几何流G_j(i)的计算步骤包括：

步骤2.3A.1、在多尺度几何分析块内进行角度的采样；

步骤2.3A.2、构造一个与子带同样大小网格点，然后计算每个网格点在采样角度上的正交投影误差，随后将与网格点位置相对应的子带系数进行一维小波变换，即实现了曲波变换；

步骤2.3A.3、最小Lagrange系数L(f_θ,R)的采样角度即为几何流G_j(i)：

L(f_θ,R)＝||f_θ-f□_θ||∧2+λ·T∧2(R_g+R_b)

式中，f_θ表示一维小波系数，R表示比特系数，f□_θ表示曲波变换所得的一维小波系数进行量化后所得的量化系数，R_g表示对几何流编码所需的比特率数量，R_b表示量化后的多尺度几何分析系数编码所需比特系数。

优选地，步骤2.3中，在步骤2.2得到的多尺度几何分析块内分别进行多尺度几何分析化得到多尺度几何分析系数，包括以下步骤：

步骤2.3B.1、沿几何流方向对小波系数进行重排，得到一维小波系数；

步骤2.3B.2、对一维小波系数做一维小波变换，得到多尺度几何分析系数。

本发明将多尺度几何分析和PCA变换结合起来融合高光谱遥感图像。采用计算图像的熵、平均梯度、标准差的方法进行评价后，通过本发明的方法获取的遥感图像质量更好，也将会更有利于遥感图像的提取、分类等后续处理。

附图说明

图1为本发明的流程图；

图2为采用本发明的方法得到的图像；

图3为采用多尺度几何分析方法得到的图像；

图4为采用PCA方法得到的图像。

具体实施方式

为使本发明更明显易懂，下面结合具体实施例，并附图作详细说明如下。应理解，这些实施例仅用于说明本发明。

结合图1，本发明提供的一种基于多尺度几何分析和PCA变换的高光谱遥感图像融合方法包括以下步骤：

步骤1：N个波段的高光谱图像预处理、配准；

高光谱图像不仅记录了地物的图像信息，同时也记录了光谱信息；高光谱图像的每一个波段都单独成像，而这些波段的成像又可以叠加，好像是经过分类之后，每一层单独成像，每一层地物具有相同的辐射值，即光谱信息。

步骤2:分别对N个波段的图像进行多尺度几何分析变换，获取每一个波段图像的多尺度几何分析系数C_j(i)和几何流G_j(i)，式中，j表示融合的源图像的个数，i表示第j幅图像中分割的块数；

(1)首先对原图像进行二维正交小波变换。假定原图像为f(x,y)，用二维尺度方程和ψ，沿着x轴方向和y轴方向先后进行分解，得到平滑逼近和细节的部分，然后分解对角高频信息。二维正交小波分解的公式如下：

式中，f_k(i,j)是分解到第k层的低频信息，而f_k ^H(i,j)、f_k ^V(i,j)、f_k ^D(i,j) 是水平、垂直、对角方向的高频信息，f_k-1(x,y)表示k-1层的高频图像，h(2x-i) 表示水平方向高频基函数，h(2y-i)表示垂直方向高频基函数，g(2y-i)表示水平方向低频基函数，g(2x-i)表示垂直方向低频基函数。

图像二维小波变换后，就可以进行二进四叉树分割。先将各高频等分成四个小块，每一小块在下一层的分割中又被分成四个小块，依次分割下去，直到最底层的块的分解尺度达到预先设定的最小尺度Jmin，最小值一般设为8×8或8×8。用四叉树表示分解过程，四叉树的叶节点对应的就是多尺度几何分析块。二进四叉树分割，依次达到最小预先设定最小尺度。多尺度几何分析块尺度过小，影响了计算速度；尺度过大，则没有达到几何方向分析的最优目的。

对于各高频子带最佳分割块大小，采用借助目标函数的最小化准则，达到最优的分割。设最底层小方块的尺寸大小为8×8，优化过程如下：

(a)将L初始化为4，对每个尺寸为L×L的方块S计算最佳方向和相应的 Lagrange函数值，将多尺度几何分析块L_0(S)初始化为本步骤计算得到的Lagrange函数值；

(b)令L＝2L，对每个尺寸为L×L的方块S计算最佳几何流方向和相应的 Lagrange函数值L(S)；

(c)对每个尺寸为L×L的方块S，分割为四个孩子，分别记为S1,S2,S3,S4，这4个孩子作为叶节点联合在一起的Lagrange函数值为：

式中，L_0(S、1)L_0S(、2)L_S、0(3L)S_分0别(表4)示四个子块的Lagrange值，λ表示Lagrange乘子，T表示量化阈值；

(d)令

(e)若若L＝2∧(j/2)j表示分割块的大小，返回步骤(b)，否则输出L_0(S) 作为四叉树分割的结果。

图像经过二进分割之后，将在多尺度几何分析块内选择最佳几何方向，同时计算多尺度几何分析系数。多尺度几何分析块的几何流是从最小尺度子带自下而上进行计算的，计算步骤如下：

1.角度的预采样：为计算几何流，先在多尺度几何分析块内进行角度的采样，假设块尺寸大小为L×L，那么将圆周角[0,π)等角度离散为2L-1个，即角度θ可能的取值为：

2.曲波变换：小波系数重排是多尺度几何分析化的前提，重排需要一个机制，Peyré等提出根据角度投影误差由小到大作为排序的机制。首先构造一个与子带同样大小网格点(设子带大小L×L)，然后计算每个网格点在采样角度上的正交投影误差τ：

τ＝-sinθ×X(i)+cosθ×Y(i)

式中，X(i)、Y(i)为第i个网格点坐标。网格点按误差值从小到大排序成为一维数组，然后将与网格点位置相对应的子带系数进行一维小波变换后即实现了曲波变换。

3.选择最优几何流方向：计算几何流是通过最小化Lagrange系数获取的，Lagrange系数表达系数L(f_θ,R)的表达式：

L(f_θ,R)＝||f_θ-f□_θ||∧2+λ·T∧2(R_g+R_b)，

式中，f_θ表示一维小波系数，R表示比特系数，f□_θ表示曲波变换所得的一维小波系数进行量化后所得的量化系数，R_g表示对几何流编码所需的比特率数量，R_b表示量化后的多尺度几何分析系数编码所需比特系数。

量化公式Q(x)为：

式中，x表示Bandelet系数，sgn(x)表示阶跃函数，q表示自适应系数。

获得最小Lagrange系数的采样角度即为最优几何正则方向，也就是要进行融合变换的几何流G_j(i)。

然后进行多尺度几何分析化。多尺度几何分析化的对象是小波变换的各高频子带的多尺度几何分析块，也就是经过二进制分割，和四叉树优化之后的小区域；因为尺寸足够小，在每个分割区域内，图像的几何正则可以用直线逼近；然后在每个分割区域内分别进行多尺度几何分析化。多尺度几何分析化的实施过程包含了两步：

1)沿几何流方向对小波系数进行重排，得到一维小波系数；

2)对一维小波系数做一维小波变换，得到多尺度几何分析系数。

多尺度几何分析化的过程必须是可逆的，才能保证多尺度几何分析变换是可以重构的。因为多尺度几何分析化两步骤中，第二步是一维小波变换，这个过程是可逆的；因此，小波系数的重排必须是可逆的。这样才能保证从一维信号准确地恢复出二维图像小波系数；为了操作过程简便，要求算法过程简单，控制参数尽可能少。采用直线逼近几何流，使得几何流只需用一个参数来描述，即几何流方向相对于水平方的倾角。

步骤3：用PCA变换分析N个波段图像的多尺度几何分析系数和几何流，获得各自的主成分，其主成分分别记为：CF(i)和GF(i)。

设X₁,X₂,…,X_p是某实际问题设计的p个随机变量，X＝(X₁,X₂,…,X_p)^T，其均值μ为：μ＝E(X)＝E(μ₁,μ₂,…,μ_p)，式中，μ₁,μ₂,…,μ_p表示随机变量的数值。

协方差矩阵Σ为：式中，σ_ij表示协方差阵。

协方差矩阵Σ是一个p阶半正定矩阵，设l_i＝(l_1i,l_2i,…,l_pi)(i＝1,2,…p)为p个常数向量，考虑如下线性组合，

容易知道：

式中，var(Y_i)表示Y_i的方差，Cov(Y_i,Y_j)表示协方差阵。

如果希望用Y₁来代替原来p个分量X₁,X₂,…,X_p这就要求Y₁尽可能多地反映原来p个分量的信息，即要求Y₁的方差var(Y₁)尽可能大。必须对l_i加以某种限制，否则var(Y₁)可以无限大。最常用的限制是要求所有l_i具有单位长度，即：

l_i ^Tl_i,i＝1,2,…,p

假设有大小为M×N像素的n幅配准好的同一地区的遥感图像，将这n幅数字图像各自排列成列向量，p＝M×N组成矩阵如下：

求取X的特征值λ₁,λ₂,…,λ_p，特征值按从大到小顺序排列，对应特征向量：

第一主成分就可以按照下面的公式求解：

F₁＝v₁₁X₁+v₁₂X₂+…+v_1pX_p

式中，F₁表示第一主成分分量。

步骤4：利用多尺度几何分析系数CF(i)和新的几何流GF(i)，多尺度几何分析逆变换，重构得到融合新图像。

根据GF(i)(几何流的)方向，按照小波系数CF(i)的排序方式(投影误差)，实现一维小波变换的逆变换和排序，恢复原来的二维小波系数。下面第一个公式利用几何流方向重排系数，正交投影；第二个公式是二维小波逆变换，即重构图像：

式中，分别表示水平方向、垂直方向、45°方向的低频分量。

图2、图3、图4三幅融合图像的水体和道路能够清晰识别，其它地物并不清晰，但是较融合源图像还是有了很大的提高。图3和图4稍微有些朦胧，重新审视整个融合过程，问题原因在于，受融合波段较多、融合图像本身和噪声的影响，二代多尺度几何分析变换得到的最优几何流方向和多尺度几何分析系数，出现了混乱矛盾。

为了从客观的角度评价三种融合方法，计算了融合图像的灰度标准差、熵、平均梯度，结果见表如下。

灰度标准差反映了灰度相对于灰度均值的离散情况，可用来评价图像反差的大小，如果一幅图像的标准差较大，则表明图像的反差较大，图像整体色调比较丰富，可观察到的信息较多；如果标准差较小，则情况相反。本发明方法的标准差均大于另外两种方法，说明从该方法融合的图像中可以观察到更多的信息，这从主观视觉上也可验证。熵是衡量该幅图像中信息量丰富的程度的指标，熵越大包含的信息也越多，本发明方法的熵也大于其它两种方法。平均梯度更多反映了图像局部的微小细节变化和纹理特征，平均梯度越大，表明图像的层次更分明，本发明方法融合的结果平均梯度最大，优于多尺度几何分析、PCA方法。

Claims (5)

1.一种基于多尺度几何分析和PCA变换的高光谱遥感图像融合方法，其特征在于，包括以下步骤：

步骤1、对N个波段的高光谱图像进行预处理、配准；

步骤2、分别对N个波段的高光谱图像进行多尺度几何分析变换，获取每一个波段高光谱图像的多尺度几何分析系数C_j(i)和几何流G_j(i)，式中，j表示融合源图像个数，i表示图像分割块数；

步骤3、用PCA变换分析N个波段高光谱图像的多尺度几何分析系数和几何流，获得各自的主成分，分别记为CF(i)和GF(i)；

步骤4、将CF(i)作为多尺度几何分析系数，将GF(i)作为几何流，进行多尺度几何分析逆变换，重构得到融合新图像。

2.如权利要求1所述的一种基于多尺度几何分析和PCA变换的高光谱遥感图像融合方法，其特征在于，所述步骤2包括：

步骤2.1、对原始的高光谱图像f(x,y)进行二维正交小波变换；

步骤2.2、对步骤2.1得到的图像进行二进四叉树分割，直到最底层块的分解尺度达到预先设定的最小尺度，获得图像的多尺度几何分析块，分割时，采用借助目标函数的最小化准则对图像各高频子带最佳分割块大小进行优化，达到最优的分割；

步骤2.3、在多尺度几何分析块内选择最佳几何方向，即几何流G_j(i)，同时计算多尺度几何分析系数。

3.如权利要求2所述的一种基于多尺度几何分析和PCA变换的高光谱遥感图像融合方法，其特征在于，步骤2.2中所述优化包括以下步骤：

步骤2.2.1、将L初始化为4，对每个尺寸为L×L的方块S计算最佳方向和相应的Lagrange函数值，将多尺度几何分析块L_0(S)初始化为本步骤计算得到的Lagrange函数值；

步骤2.2.2、令L＝2L，对每个尺寸为L×L的方块S计算最佳几何流方向和相应的Lagrange函数值L(S)；

步骤2.2.3、对每个尺寸为L×L的方块S，分割为四个孩子，分别记为S1,S2,S3,S4，这4个孩子作为叶节点联合在一起的Lagrange函数值为：

式中，L_0(S、1)L_0S(、2)L_S、0(3L)S_分0别(表4)示四个子块的Lagrange值，λ表示Lagrange乘子，T表示量化阈值；

步骤2.2.4、令

步骤2.2.5、若L＝2∧(j/2)，j表示图像分割大小，返回步骤2.2.2，否则输出L_0(S)作为四叉树分割的结果。

4.如权利要求2所述的一种基于多尺度几何分析和PCA变换的高光谱遥感图像融合方法，其特征在于，步骤2.3中，几何流G_j(i)的计算步骤包括：

步骤2.3A.1、在多尺度几何分析块内进行角度的采样；

步骤2.3A.2、构造一个与子带同样大小网格点，然后计算每个网格点在采样角度上的正交投影误差，随后将与网格点位置相对应的子带系数进行一维小波变换，即实现了曲波变换；

步骤2.3A.3、最小Lagrange系数L(f_θ,R)的采样角度即为几何流G_j(i)：

L(f_θ,R)＝||f_θ-f□_θ|||∧2+λ·T∧2(R_g+R_b)

式中，f_θ表示一维小波系数，R表示比特系数，表示曲波变换所得的一维小波系数进行量化后所得的量化系数，R_g表示对几何流编码所需的比特率数量，R_b表示量化后的多尺度几何分析系数编码所需比特系数。

5.如权利要求2所述的一种基于多尺度几何分析和PCA变换的高光谱遥感图像融合方法，其特征在于，步骤2.3中，在步骤2.2得到的多尺度几何分析块内分别进行多尺度几何分析化得到多尺度几何分析系数，包括以下步骤：

步骤2.3B.1、沿几何流方向对小波系数进行重排，得到一维小波系数；

步骤2.3B.2、对一维小波系数做一维小波变换，得到多尺度几何分析系数。