基于稀疏张量和多视图特征的遥感图像融合方法与流程

本发明属于图像处理技术领域，更进一步涉及遥感图像融合技术领域中的一种基于稀疏张量和多视图特征的遥感图像融合方法。本发明用于融合同一个传感器获得的关于同一个场景的多光谱图像和全色图像融合，在保持多光谱图像光谱信息的同时提高图像的空间分辨率，可应用于地形测绘与地图更新、地物分类、遥感监测和森林资源调查、军事侦察等领域。

背景技术：

目前，现有的遥感图像融合方法最具代表性的是基于模型类的融合方法。该方法结合了稀疏表示理论使得融合后的图像取得了更高的空间分辨率，它的优点在于实现简单、速度快。但是存在的缺点是缺少了对于多光谱图像的谱间信息的考虑，同时只是基于图像灰度信息的处理视角过于单一，使得融合后的多光谱图像产生了光谱扭曲。

jiangc,zhangh,shenh,etal在其发表的论文“two-stepsparsecodingforthepan-sharpeningofremotesensingimages”(《selectedtopicsinappliedearthobservationsandremotesensing》ieeejournalof,2014,7(5):1792-1805)中提出了一种两步稀疏编码的遥感图像融合方法。该方法的步骤是，第一步稀疏编码过程是低分辨多光谱图像的图像块由对应位置的低分辨全色图像块线性表示。第二步稀疏编码是通过低分辨字典来稀疏表示残差。两步稀疏编码结束后，得到融合后的多光谱图像。该方法存在的不足之处是，由于该方法将多光谱图像分波段融合，没能充分考虑对波段间的信息充分利用，使得融合结果产生了色彩失真。

上海交通大学在其拥有专利技术“一种基于稀疏表示的遥感图像融合方法”(申请日：2013年03月29日，申请号：201310108594.3，公开号：103208102b)中公开了一种基于稀疏表示的遥感图像融合方法。首先建立多光谱图像与其亮度分量之间的线性回归模型；其次利用训练的高、低分辨率字典分别对全色图像和多光谱图像进行稀疏表示，并根据线性回归模型获得多光谱图像亮度分量稀疏表示系数；然后根据全色图像和亮度分量的稀疏表示系数提取细节成分，并在通用分量替换融合框架下注入到多光谱图像各波段的稀疏表示系数中；最后进行图像复原得到高空间分辨率的多光谱图像。该方法存在的不足是，从单一视图对多光谱和全色图像进行考虑，使得融合后的图像的空间分辨率受到制约；其次，该方法没有考虑到多光谱的谱间关系，导致最终的融合图像的光谱有扭曲。

技术实现要素：

本发明的目的是针对上述现有技术的不足，提出一种基于稀疏张量和多视图特征的遥感图像融合方法。

本发明的思路是，本发明根据待融合图像的灰度特征、纹理特征、边缘特征构成图像的多视图特征，基于非采样轮廓变换的框架对低通系数进行融合，计算多视图特征块的特征接近度并构成多模字典，在多模字典下对带通系数融合，经过逆变换后得到高分辨的多光谱图像。由于利用了低分辨多光谱图像多视图特征矩阵和全色图像多视图特征矩阵和张量基追踪方法对融合后图像空谱信息进行增强，提高了融合图像的清晰度并光谱扭曲，克服了现有图像融合技术处理视角过于单一和欠缺对于多光谱图像谱间关系考虑的缺点。

为实现上述目的，本发明包括如下步骤：

(1)输入待融合图像：

分别输入待融合的低分辨多光谱图像和全色图像；

(2)提取低通和带通系数：

(2a)从低分辨多光谱图像中提取低分辨多光谱主成分低通系数和低分辨多光谱主成分带通系数；

(2b)利用非采样轮廓波变换方法，对全色图像进行分解，得到全色图像低通系数和全色图像带通系数；

(3)提取边缘特征：

(3a)利用一二阶梯度算子和拉普拉斯算子，对低分辨多光谱主成分进行梯度变换，得到低分辨多光谱主成分边缘特征；

(3b)利用一二阶梯度算子和拉普拉斯算子，对全色图像进行梯度变换，得到全色图像边缘特征；

(4)获得多视图特征矩阵：

(4a)将低分辨多光谱主成分、低分辨多光谱主成分带通系数、低分辨多光谱边缘特征按第三维方向依次排列，组成低分辨多光谱多视图特征矩阵；

(4b)将全色图像、全色图像带通系数、全色图像边缘特征按第三维方向依次排列，组成全色图像多视图特征矩阵；

(5)获得多视图特征小块：

(5a)利用列向量化方法，对低分辨多光谱多视图特征矩阵进行列向量化，得到低分辨多光谱多视图特征小块；

(5b)利用列向量化方法，对全色图像多视图特征矩阵进行列向量化，得到全色图像多视图特征小块；

(6)获得特征接近度：

(6a)按照图像块的排列顺序，依次选取低分辨多光谱多视图特征小块中的一个小块；

(6b)按照下式，计算所选的低分辨多光谱多视图特征小块与全部的全色图像多视图特征小块的特征接近度：

其中，fij表示第i个多光谱多视图特征小块与第j个全色图像多视图特征小块的特征接近度，min表示取最小值操作，max表示取最大值操作，var表示协方差操作，t表示多光谱多视图特征小块和全色图像多视图特征小块按第二维度方向组合得到的矩阵；

(6c)判断是否选完所有的图像块，若是，则执行步骤(7)，否则，执行步骤(6a)

(7)获得因子矩阵：

利用索引映射法，将特征接近度映射与其对应的全色图像，获得因子矩阵；

(8)将低分辨多光谱带通系数按照第三维方向组合，得到低分辨多光谱带通系数张量块；

(9)将全色图像带通系数按照第三维方向组合，得到全色图像带通系数张量块；

(10)获得多模字典：

利用张量模展开方法，对低分辨多光谱带通系数进行模方向展开，获得多模字典；

(11)获得稀疏系数：

(11a)使用张量基追踪方法，获得低分辨多光谱带通系数在多模字典下的多光谱稀疏系数；

(11b)使用张量基追踪方法，获得全色图像带通系数在多模字典下的全色稀疏系数；

(12)获得融合后的带通系数：

将多光谱稀疏系数与全色稀疏系数相加，得到融合后的带通系数；

(13)获得融合后的多光谱图像主成分：

利用非采样的轮廓波变换的逆变换，将融合后的带通系数与低分辨多光谱低通系数逆变换，得到融合后的多光谱图像主成分；

(14)获得融合后的多光谱图像：

利用主成分分析逆变换方法将，将融合后的多光谱图像主成分逆变换，得到融合后的多光谱图像。

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

第一，本发明能够根据低分辨多光谱图像多视图特征矩阵和全色图像多视图特征矩阵，对低分辨多光谱图像和全色图像进行融合，克服了现有技术存在的只是基于图像灰度信息的处理视角过于单一的问题，使得本发明在融合过程中增强了对空间信息的捕捉能力，融合结果具有更高的空间分辨率。

第二，本发明利用张量基追踪方法，对低分辨多光谱图像的带通系数进行重构，克服了现有技术中欠缺对于多光谱图像的谱间关系的考虑，使得本发明能够显著的改善融合结果的光谱扭曲现象。

附图说明

图1为本发明的流程图；

图2为本发明的仿真图。

具体实施方式

下面结合附图对本发明作进一步的描述。

参考附图1，对本发明实现的步骤作进一步的详细描述。

步骤1，输入待融合图像。

分别输入待融合的低分辨多光谱图像m和全色图像p。

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

步骤2，提取低分辨多光谱主成分c1的低通系数lm与带通系数hm。

第一步，利用主成分分析方法，对低分辨多光谱图像进行主成分分析，得到低分辨多光谱的主成分c1，c1的大小为256×256。

第二步，利用非采样轮廓波变换方法，对低分辨多光谱主成分c1进行分解，得到低分辨多光谱主成分低通系数lm和低分辨多光谱主成分带通系数hm。

本发明的实施例在非采样轮廓波变换中，使用的金字塔滤波器类型为“9-7”，方向滤波器组类型为“pkva”，低分辨多光谱主成分低通系数lm大小为256×256，低分辨多光谱主成分带通系数hm大小为256×256×4。

步骤3，提取全色图像p的低通系数lp和带通系数hp。

利用轮廓波变换方法，对全色图像进行分解，得到全色图像的低通系数lp和全色图像带通系数hp。

本发明的实施例中，全色图像低通系数lp大小为256×256，全色图像带通系数hp大小为256×256×4。

步骤4，提取待融合图像边缘特征fp。

第一步，利用一二阶梯度算子，经过梯度变换，得到主成分c1的水平方向一阶梯度f1、垂直方向上一阶梯度f2、水平方向二阶梯度f3、垂直方向二阶梯度f4。

第二步，将主成分c1的水平方向一阶梯度f1、垂直方向上一阶梯度f2、水平方向二阶梯度f3、垂直方向二阶梯度f4按第三维的方向组合得到低分辨多光谱主成分边缘特征fm。

第三步，利用一二阶梯度算子，经过梯度变换，得到全色图像p的水平方向一阶梯度f5、垂直方向一阶梯度f6、水平方向二阶梯度f7、垂直方向上二阶梯度f8；

第四步，将全色图像p的水平方向一阶梯度f5、垂直方向一阶梯度f6、水平方向二阶梯度f7、垂直方向上二阶梯度f8按第三维的方向组合得到全色图像边缘信息fp。

步骤5，获得低分辨多光谱的多视图特征矩阵bm。

将低分辨多光谱主成分c1、低分辨多光谱主成分带通系数hm、低分辨多光谱边缘特征fm按第三维方向组合排列，得到低分辨多光谱多视图特征矩阵bm。

步骤6，获得全色图像多视图特征矩阵bp。

将全色图像p、全色图像带通系数hp、全色图像边缘特征fp按第三维方向组合排列，组成低分辨多光谱多视图特征矩阵bp。

步骤7，获得低分辨多光谱多视图特征小块nm。

利用列向量化方法，对低分辨多光谱多视图特征矩阵bm进行列向量化，得到低分辨多光谱多视图特征小块nm。

步骤8，获得全色图像多视图特征小块np。

利用列向量化方法，对全色图像多视图特征矩阵bp进行列向量化，得到全色图像多视图特征小块np。

步骤9，获得特征接近度。

第一步，按照图像块的排列顺序，依次选取低分辨多光谱多视图特征小块nm的一个小块。

第二步，按照下式，计算所选的低分辨多光谱多视图特征小块与全部的全色图像多视图特征小块的特征接近度：

其中，fij表示第i个多光谱多视图特征小块与第j个全色图像多视图特征小块的特征接近度，min表示取最小值操作，max表示取最大值操作，var表示协方差操作，t表示多光谱多视图特征小块和全色图像多视图特征小块按第二维度方向组合得到的矩阵。

第三步，判断是否选完所有的图像块，若是，则执行步骤10，否则，执行步骤9的第一步。

步骤10，获得因子矩阵d1。

第一步，利用列向量化方法，对全色图像进行列向量化，得到全色图像列向量。

第二步，对所有的特征接近度由高到低排序，保留前200个特征接近度的索引。

第三步，将特征接近度的索引映射到全色图像列向量。

第四步，在全色图像列向量中，选取前200个特征接近度索引对应的全色图像列向量。

第五步，将全色图像列向量组合成矩阵，得到因子矩阵d1。

步骤11，将低分辨多光谱带通系数hm按照第三维方向组合，得到低分辨多光谱带通系数张量块e1。

步骤12，将全色图像带通系数hp按照第三维方向组合，得到全色图像带通系数张量块e2。

步骤13，获得多模字典dn。

第一步，将低分辨多光谱带通系数张量块e1按照模2方向展开，得到低分辨多光谱模2矩阵d2；

第二步，将低分辨多光谱带通系数张量块e1张量块按照模3方向展开，得到低分辨多光谱模3矩阵d3；

第三步，将因子矩阵d1、低分辨多光谱模2矩阵d2、低分辨多光谱模3矩阵d3组合得到多模字典dn。

步骤14，使用张量基追踪方法，获得低分辨多光谱带通系数张量块e1在多模字典dn下的多光谱稀疏系数a1。

所述张量基追踪公式如下：

其中，min表示取最小值操作，a表示待求解的稀疏系数，m表示带通系数，d表示多模字典，ε表示重构误差，s.t表示限制条件的符号，||·||1表示一范数操作，表示二范数平方操作。

步骤15，使用张量基追踪方法，获得全色图像带通系数张量块e2在多模字典dn下的全色稀疏系数a2。

所述张量基追踪公式如下：

其中，min表示取最小值操作，a表示待求解的稀疏系数，m表示带通系数，d表示多模字典，ε表示重构误差，s.t表示限制条件的符号，||·||1表示一范数操作，表示二范数平方操作。

步骤16，获得融合后的带通系数an。

将多光谱稀疏系数a1与全色稀疏系数a2相加，得到融合后的带通系数an＝a1+a2。

步骤17，获得融合后的多光谱图像主成分c2。

利用非采样的轮廓波变换的逆变换，将融合后的带通系数hn与低分辨多光谱低通系数lm逆变换，得到融合后的多光谱图像主成分c2。

步骤18，获得融合后的多光谱图像。

利用主成分分析逆变换方法将，将融合后的多光谱图像主成分c2逆变换，得到融合后的多光谱图像。

下面结合附图2对本发明的仿真效果做进一步说明。

1.仿真条件：

本发明的仿真所使用待融合的图像由geoeye数据集提供，包括2.0m空间分辨率的低分辨多光谱图像和0.5m空间分辨率的高分辨全色图像。低分辨多光谱图像具有四个光谱带：蓝色(b)，绿色(g)，红色(r)和近红外(nir)。

仿真环境：matlabr2014aonpcwithintel(r)core(tm)/2.50g/2.5g。

2.仿真结果与分析：

图2(a)是本发明仿真实验中输入的待融合低分辨的多光谱图像，图2(b)是本发明仿真实验中输入的待融合高分辨的全色图像，图2(c)是高分辨多光谱图像，作为本发明仿真实验中的参考图像，图2(d)是采用本发明对图2(a)和图2(b)融合后获得的高分辨多光谱图像。

由图2(d)与图2(a)的比较可见，图2(a)的图像细节不清晰，如图像的道路的边缘模糊，图像中的白色小块不清楚，而图2(d)在空间分辨率、光谱信息、视觉效果上明显优于图2(a)，且更加接近于图2(c)的视觉效果，图2(d)中道路轮廓更加清晰，白色小块的空间分辨率和光谱分辨率有明显提升。

本发明则很好地利用了附图2(a)和附图2(b)中待融合图像自身的多视图特征，对融合后图像的空间分辨率进行增强。此外，对多光谱图像谱间关系的考虑，使得本发明显著的改善了融合结果的光谱扭曲现象。

综上所述，可以看出本发明能够很好的提升待融合图像的空间信息和光谱信息，克服了现有图像融合技术中视图单一、缺少对图像谱间信息考虑的缺点。