基于深度稀疏自编码的多尺度几何遥感图像融合方法与流程

本发明属于图像处理技术领域，是一种多尺度几何遥感图像融合方法，可用于目标识别、地物分类、遥感监测和森林资源调查。

背景技术：

遥感图像，包含全色图像和多光谱图像。全色图像，是指包含了0.38-0.76μm的全部可见光波段的图像，传感器获取的全色图像一般为灰度图像。多光谱图像，是指利用两个以上的传感器各自接收地物反射的电磁波信息的不同波段所得到的一种特殊图像，传感器获取的多光谱图像一般为四个波段。

多光谱图像与全色图像的融合旨在将信息优势互补起来获得一幅同时具有高空间分辨率和高光谱分辨率的融合图像。由于设备硬件条件和成像环境方面的限制，为了获得同时具有高空间分辨率和高光谱分辨率的融合图像，需要使用遥感图像融合技术来实现图像信息的融合，现有的融合算法主要集中在像素级的多分辨率融合。现有的融合算法可归结为三类：1)成分替换法；2)多尺度细节注入法；3)基于模型的复原法。影响融合结果的因素主要分为两个方面：1)因替换或注入的高分辨空间信息不匹配或信息冗余而造成融合结果的颜色失真及光谱扭曲；2)抽取的高分辨率信息与多光谱图像缺少的信息量不一致时，使得融合后图像的空间分辨率受到制约。

多尺度分析工具具有多尺度，多方向，各向异性的特性，轮廓波变换是一种有很好的空域和频域局部特性、多方向性、多分辨率、各向异性的完备变换。可以用不同频率和尺度的子带更加精准的保留图像中的轮廓和纹理细节信息。而非采样的轮廓波变换在继承了轮廓波变换的特征的基础上还消除了吉布斯现象，并且具有平移不变性，是一种超完备的多尺度变换方法。非采样的轮廓波变换由非采样金字塔滤波器组和非采样方向滤波器组两部分组成。

近年来，压缩感知和稀疏表示的理论的发展使得基于退化模型为理论的模型复原的图像融合方法得以发展。参见Shutao Li,and Bin Yang,“A new pan-sharpening method using a compressed sensing technique,”IEEE Transactions on Geoscience and Remote Sensing,vol.49,no.2,pp.738-746,February 2011。该方法中稀疏先验信息被用来重构高分辨率多光谱图像，通过高分辨率的字典原子线性组合得到高分辨的多光谱图像。该方法的局限性在于构造字典原子的真实数据集无法通过遥感卫星获得，字典原子包含的信息不能充分的表示高分辨的多光谱图像，造成融合结果的光谱扭曲。

技术实现要素：

本发明目的在于结合非采样的轮廓波变换提取细节的精准性，克服基于模型复原法中字典原子的局限性，提出一种基于深度稀疏自编码的多尺度几何遥感图像融合方法，以消除冗余，保持光谱分辨率，提高融合结果的高空间分辨率。

实现本发明目的的技术方案是：将高空间分辨率全色图像和低空间分辨率多光谱图像信息优势进行互补，消除冗余得到一幅同时具有高空间分辨率和高光谱分辨率的融合图，其步骤包括如下：

(1)分别输入低分辨率多光谱图像和高分辨率全色图像，并提取低分辨多光谱图像的第一主成分C1；

(2)通过非采样轮廓波变换，分别获得低分辨多光谱图像的第一主成分低通系数L_M与带通系数H_M，及高分辨全色图像的低通系数L_P与带通系数H_P；

(3)将低分辨多光谱图像的第一主成分低通系数L_M与全色图像的低通系数L_P融合，获得融合后的低通系数L_N；

(4)获得空间自相似字典D_S：

4a)将多光谱图像第一主成分的带通系数H_M和全色图像的带通系数H_P进行分块，得到多光谱图像第一主成分的带通系数小块M1和全色图像的带通系数小块P1；

4a)分别对每一个光谱图像第一主成分的带通系数小块M1和全色图像的带通系数小块P1进行近邻扩充，得到扩充后的第一主成分带通系数I_M和扩充后的全色图像的带通系数I_P，该I_M和I_P的中心分别为H_M和H_P；

4b)分别对扩充后的多光谱图像的第一主成分带通系数I_M和扩充后的全色图像的带通系数I_P以1为步长划窗取块,再将块列向量化作为字典原子,将扩充后的第一主成分带通系数I_M的字典原子进行组合得到多光谱图像第一主成分字典D_M，将扩充后的全色图像的带通系数I_P的字典原子进行组合得到全色字典D_P；将多光谱图像第一主成分字典D_M与全色字典D_P进行组合，获得空间自相似字典D_S；

(5)将低分辨多光谱图像的第一主成分带通系数H_M与全色图像的带通系数H_P进行融合，获得融合后的带通系数H_N；

(6)根据融合后的带通系数H_N获得多级残差稀疏自编码的带通系数

6a)计算融合后的带通系数H_N分别与低分辨多光谱图像的第一主成分带通系数H_M的残差R1和与全色图像的带通系数H_P的残差R2；

6b)将步骤6a)获得的两个残差R1和R2按照区域能量最大的规则进行融合，获得带通残差系数

6c)将带通残差系数与融合后的带通系数H_N相加，获得残差稀疏自编码带通系数

6d)重复6a)-6c)共S次，获得多级残差稀疏自编码带通系数1≤S≤4；

(7)对融合后的低通系数L_N和多级残差稀疏自编码带通系数同时进行非采样轮廓波变换的逆变换，获得融合后的第一主成分C2；

(8)获得空域多级深度稀疏自编码主成分C_S：

8a)计算融合后的第一主成分C2与低分辨多光谱图像第一主成分C1的残差R3，计算融合后的第一主成分C2与全色图像的残差R4；

8b)将8a)中获得的两个残差R3和R4按照深度稀疏自编码的融合方法得到融合后的空域残差

8c)将融合后的空域残差与融合后的第一主成分C2相加获得空域深度稀疏自编码主成分C₁；

8d)重复8a)-8c)共S次，获得空域多级深度稀疏自编码主成分C_S，1<S≤4；

(9)将空域多级深度稀疏自编码主成分C_S进行逆主成分分析变换，获得高分辨多光谱图像。

本发明与现有的技术相比具有以下优点：

第一，本发明由于对低通系数在过完备离散余弦变换字典下进行稀疏表示，保证了细节信息的完整；同时由于本发明依据图像区域局部相似性的特征，构造带通系数的空间自相似字典，并对带通系数在空间自相似字典下进行稀疏编码，克服了因抽取的高分辨率信息与多光谱图像缺少的信息量不一致所造成的空间分辨率的降低，克服了不匹配注入造成的色彩失真及光谱扭曲，提高了融合后图像的空间分辨率，保持多光谱图像的光谱信息。

第二，本发明由于将稀疏自编码推广至空间自相似字典下的多级稀疏自编码，引入了非采样轮廓波变换和空间自相似字典下多级稀疏自编码相结合的融合算法，实现了对图像纹理和轮廓的准确提取，减少了不匹配细节的注入，使得注入到多光谱图像的高频信息更加精准，改善了融合后多光谱图像空间扭曲。

仿真结果表明，本发明得到的融合后多光谱图像比现有技术的融合图像能更好的保持光谱分辨率，提高空间分辨率。

附图说明

图1为本发明的实现流程图；

图2为本发明的仿真结果图。

具体实施方法

以下结合附图为本发明做进一步详细描述：

参照图1，本发明的具体实施方式如下：

步骤1，分别输入低分辨率多光谱图像M和高分辨率全色图像P，并提取低分辨多光谱图像的第一主成分C1；

本发明实施例中输入的低分辨多光谱图像大小为64×64×4，分辨率为2m；高分辨率的全色图像大小为256×256，分辨率为0.5m；

将低分辨多光谱图像经过4倍上采样变换，使它的大小达到256×256×4；

对上采样后的多光谱图像进行主成分分析变换PCA，提取主成分分析变换后的第一主成分，定义该主成分为多光谱图像的第一主成分C1，C1的大小为256×256。

步骤2，通过非采样轮廓波变换获得低分辨多光谱图像的第一主成分低通系数L_M与带通系数H_M。

非采样轮廓波变换，是由金字塔滤波器和方向滤波器组两部分组成，在非采样轮廓波变换过程中先由金字塔滤波器进行多尺度变换，再由方向滤波器组进行多方向变换。

本实施例中使用的金字塔滤波器类型为“9-7”，方向滤波器组类型为“pkva”。对低分辨多光谱图像的第一主成分C1进行非采样轮廓波变换，其操作如下：

首先，将多光谱图像的第一主成分先经过非采样轮廓波变换的金字塔滤波器进行多尺度变换，得到多光谱图像第一主成分的低通系数L_M；

然后，将多光谱图像的第一主成分经过非采样轮廓波变换的方向滤波器组进行多方向变换，得到多光谱图像第一主成分的带通系数H_M。

步骤3，通过非采样轮廓波变换获得高分辨全色图像P的低通系数L_P与带通系数H_P。

3.1)将全色图像先经过非采样轮廓波变换的金字塔滤波器进行多尺度变换，得到全色图像的低通系数L_P；

3.2)再将全色图像经过非采样轮廓波变换的方向滤波器组进行多方向变换，得到全色图像的带通系数H_P。

步骤4，生成二维的过完备离散余弦变换字典D。

4.1)定义二维过完备离散余弦变换字典D中原子的长度l，和二维过完备离散余弦变换字典D的原子个数N，本实施例中l＝25，N＝256；

4.2)生成一维离散余弦变换字典D¹：

$D^1=\mathrm{\&alpha;}\left(k\right)cos\frac{nk\mathrm{\&pi;}}{\sqrt{N}},$

其中：n是一维离散余弦变换字典D¹中原子的长度，k是一维离散余弦变换字典D¹的原子个数，

4.3)计算二维的过完备离散余弦变换字典是张量乘积的符号。

步骤5，获得多光谱图像第一主成分的低通稀疏系数A1。

5.1)将多光谱图像第一主成分的低通系数L_M按5×5大小进行分块，得到多光谱图像第一主成分的低通系数小块；

5.2)将多光谱图像第一主成分的低通系数小块列向量化，得到多光谱图像第一主成分的低通系数列向量；

5.3)根据每一个多光谱图像第一主成分的低通系数列向量在离散余弦变换字典D下的稀疏表示：L_M＝DA1，得到多光谱图像第一主成分低通稀疏系数A1。

步骤6，获得全色图像低通稀疏系数A2。

6.1)将全色图像的低通系数L_M按5×5大小进行分块，得到全色图像的低通系数小块；

6.2)将全色图像的低通系数小块列向量化，得到全色图像的低通系数列向量；

6.3)根据每一个全色图像的低通系数在离散余弦变换字典D下的稀疏表示：L_P＝DA2，得到全色图像低通稀疏系数A2；

步骤7，获得融合后的低通系数L_N。

7.1)将步骤5和步骤6得到的两个低通稀疏系数A1和A2进行并集，得到融合后的低通稀疏系数对于这两个低通稀疏系数A1和A2相同位置元素均为非零的情况，则选取模值最大的低通系数作为融合后的低通稀疏系数

7.2)将融合后的低通稀疏系数与过完备离散余弦变换字典D相乘，获得融合后的低通系数L_N。

步骤8，获得多光谱图像第一主成分字典D_M。

8.1)将多光谱图像第一主成分的带通系数H_M按5×5大小进行分块，得到多光谱图像第一主成分的带通系数小块M1；

8.2)本实施例中使用大小为21的邻域对多光谱图像第一主成分的带通系数小块M1进行近邻扩充，获得大小为21×21的扩充后的多光谱图像第一主成分带通系数I_M，且I_M的中心为M1；

8.3)对每一个扩充后的多光谱图像第一主成分带通系数I_M用5×5大小的窗口，以1为步长，按“Z”字形划窗取块，得到多光谱带通系数小块；

8.4)将多光谱带通系数小块列向量化，每一个多光谱带通系数小块列向量化后，得到大小为25×1的多光谱带通系数列向量；

8.5)将多光谱带通系数列向量进行组合得到大小为25×289的多光谱图像第一主成分字典D_M。

步骤9，获得全色字典D_P。

9.1)将全色图像的带通系数H_P按5×5大小进行分块，得到全色图像带通系数小块P1；

9.2)本实施例中使用大小为21的邻域对全色图像带通系数小块P1进行近邻扩充，获得大小为21×21的扩充后的全色图像带通系数I_P，且I_P的中心为P1；

9.3)对每一个扩充后的全色图像带通系数I_P用5×5大小的窗口，以1为步长，按“Z”字形划窗取块，得到全色图像系数小块；

9.4)对全色图像带通系数小块列向量化，每一个全色图像带通系数小块列向量化后，得到大小为25×1的全色带通系数列向量；

9.5)将全色带通系数列向量进行组合，得到大小为25×289的全色字典D_P。

步骤10，获得空间自相似字典D_S。

本实施例中将大小为25×289的全色字典D_M和大小为25×289的全色字典D_P进行组合，获得大小为25×578的空间自相似字典D_S。

步骤11，获得多光谱图像第一主成分的带通稀疏系数A3。

11.1)将多光谱图像第一主成分的带通系数小块M1列向量化，得到多光谱图像第一主成分的带通系数列向量X_M；

11.2)根据每一个多光谱图像第一主成分的带通系数列向量X_M在对应的空间自相似字典D_S下的稀疏表示：X_M＝D_SA3，得到多光谱图像第一主成分低通稀疏系数A3。

步骤12，获得全色图像带通稀疏系数A4。

12.1)将全色图像的带通系数小块P1列向量化，得到全色图像的带通系数列向量X_P；

12.2)根据每一个全色图像的带通系数列向量在对应的空间自相似字典D_S下的稀疏表示：X_P＝D_SA4，得到全色图像低通稀疏系数A4；

步骤13，获得融合后的带通系数H_N。

13.1)将步骤11和步骤12得到的两个带通稀疏系数A3和A4进行并集，得到融合后的带通稀疏系数当于这两个带通稀疏系数A3和A4在相同位置的元素均为非零的情况时，则选取模值最大的带通系数作为融合后的带通稀疏系数

13.2)将每一个融合后的带通稀疏系数与对应的空间自相似字典D_S相乘，获得融合后的带通系数H_N。

步骤14，获得多级残差稀疏自编码的带通系数

14.1)计算融合后的带通系数H_N与低分辨多光谱图像的第一主成分带通系数H_M的残差R1＝H_M-H_N；

14.2)计算融合后的带通系数H_N与全色图像带通系数H_P的残差R2＝H_p-H_N；

14.3)将步骤14.1)和步骤14.2)得到的两个带通系数残差R1和R2按照区域能量最大的规则进行融合，获得带通残差系数

14.4)将带通残差系数与融合后的带通系数H_N相加，获得残差稀疏自编码带通系数

14.5)重复步骤14.1)-14.4)的操作共4次，获得多级残差稀疏自编码带通系数

步骤15，获得融合后的第一主成分C2。

将步骤7获得的融合后的低通系数L_N和步骤14获得的多级残差稀疏自编码带通系数，进行非采样轮廓波的逆变换过程，得到融合后的第一主成分C2。

步骤16，获得空域多级深度稀疏自编码主成分C_S。

16.1)计算融合后的第一主成分C2与低分辨多光谱图像第一主成分C1的残差R3＝C1-C2；

16.2)计算融合后的第一主成分C2与全色图像P的残差R4＝P-C2；

16.3)将步骤16.1)和16.2)中得到的两个残差R3和R4，按照深度稀疏自编码的融合方法，得到融合后的空域残差

16.4)将融合后的空域残差与融合后的第一主成分C2相加，获得空域深度稀疏自编码主成分

16.5)重复16.1)-16.4)共4次，获得空域多级深度稀疏自编码主成分C_S，1<S≤4。

步骤17，获得高分辨多光谱图像。

将步骤16中获得的空域多级深度稀疏自编码主成分C_S，经过逆主成分分析变换，获得高分辨多光谱图像。

本发明的效果可以通过以下仿真实验作进一步说明。

1.实验仿真环境为MATLAB R2012a on PC with Intel(R)core(TM)/3.10G/4G，在该仿真环境下，

2.仿真内容

仿真1，使用本发明方法，对低分辨的多光谱图像和高分辨的全色图像进行融合，得到高分辨的多光谱图像，结果如图2。其中

图2(a)是低分辨的多光谱图像，图2(b)是高分辨的全色图像，图2(c)是参考图像，这三幅图像选自GeoEye，图2(d)是采用本发明对图2(a)和图2(b)进行融合后获得的高分辨多光谱图像。

由图2(d)与图2(a)的比较可见，图2(a)的细节不清晰而且光谱也保持的不好，如道路的清晰度，房屋的色彩等，而图2(d)在空间分辨率、光谱信息、视觉效果上明显优于图2(a)好。

仿真2，为了证明本发明的效果，用本发明的方法和现有的代表性融合算法，分别对图2(a)和图2(b)的图像进行融合，并通过常用评价指标进行对比，结果如表1。

表1各算法融合结果评价指标

表1中的常用评价指标包括：相关系数CC、通用图像质量指标UIQI、均方根误差RMSE、整体图像质量指数Q4、光谱弧度SAM、相对整体维数综合误差ERGAS。

现有的代表性融合算法包括：基于主成分分析变换的融合算法PCA、基于小波变换的融合算法AWLP、基于非采样轮廓波变换的融合算法NSCT、基于压缩感知融合算法CS。

从表1数据对比表明，本发明融合结果的相关系数CC、通用图像质量指标UIQI和整体图像质量指数Q4均大于现有代表性融合算法的指标；其均方根误差RMSE、光谱弧度SAM和相对整体维数综合误差ERGAS均小于代表性融合算法的指标。

可见，本发明的融合图像优于现有代表性融合算法的结果，且与参考图像最接近，光谱扭曲程度最小，图像整体质量最优，融合结果最好。