一种基于导向滤波的遥感图像融合方法与流程

本发明涉及遥感图像的可视化增强处理技术，尤其是针对星载多光谱图像和全色光图像的图像融合方法。

背景技术：

由于传感器自身的物理局限性以及数据传输技术的壁垒，使得当前星载传感器仅能提供光谱分辨率高但空间分辨率低的多波段多光谱图像以及空间分辨率高但光谱分辨率低的单一波段全色光图像。在轨的星载传感器，例如quickbird、geoeye-1、worldview-2等能够同时提供多波段多光谱遥感图像和单波段高空间分辨率的全色光图像。通过融合多光谱图像和全色光图像，可以在保持多光谱图像光谱信息的同时提高其空间分辨率，达到增强可视化的效果。

全色光图像空间分辨率高于多光谱图像，因此融合时首先要将多光谱图像重采样放大至与全色光图像同尺寸，然后再将放大后的多光谱图像和原始全色光图像输入至融合模型，输出得到融合的多光谱图像。在对多光谱图像的重采样处理上，目前国内外通用的做法是选用最近邻插值法、双线性插值法或者双三次插值法对多光谱图像的各波段独立地进行插值放大(参见文献ieeetransactionsongeoscienceandremotesensing,45(10):3012-3021,2007)。事实上，由上述插值法得到的多光谱图像会出现不同程度的光谱失真与空间信息失真，从而进一步影响融合图像的信息保持(参见文献ieeegeoscienceandremotesensingletters,4(1):27-31,2007；ieeejournalofselectedtopicsinsignalprocessing,5(3):446-453,2011)。

技术实现要素：

为了克服现有技术的不足，本发明提供一种基于导向滤波的图像融合方法，通过导向滤波联合全色光图像完成对多光谱图像的重采样放大，再利用intensity-hue-saturation(ihs)变换分离光谱信息，最终将全色光图像的结构信息提取并注入至多光谱图像中，有效地提高了多光谱图像的空间分辨率。

本发明解决其技术问题所采用的技术方案包括以下步骤：

第一步，通过卷积运算对全色光图像pan进行低通滤波，低通滤波器记为lpf，则输出图像pl＝pan*lpf；对输出图像pl作下2抽样处理，输出结果记为p2；

第二步，设置图像p2为导向滤波的指导图像i，任意第j个多光谱波段mj为导向滤波的输入图像p，导向滤波输出图像记为q；像素点i位于以像素点k为中心的窗口ωk中，该窗口半径r≥1，窗口大小为(2r+1)×(2r+1)像素；滤波输出图像q在像素点i的滤波输出值qi由指导图像i在该点的像素值ii以及i点所在局部窗口的窗口系数ak和bk共同决定，矩阵形式为q＝a×i+b，其中a和b是窗口系数矩阵，窗口系数窗口系数|ω|为窗口ωk中所有像素的个数，pi为输入图像p在i处的像素值，μk和σk分别是指导图像在窗口ωk中的均值和方差，是输入图像p在窗口ωk中的均值，ε是规则化参数且ε∈(0,1)；

第三步，利用最近邻插值法对窗口系数矩阵a和b重采样放大，分别得到对应的新窗口系数和

第四步，利用原始全色光图像pan和新窗口系数和计算大尺寸的输出图像该输出图像与原始全色光图像尺寸相同；

第五步，从n波段多光谱图像中任选三波段多光谱图像mt1、mt2和mt3；重复第二步至第四步，计算出对应的三幅大尺寸输出图像和计算强度分量

第六步，对所选的三波段多光谱图像mt1、mt2、mt3进行ihs变换，得到利用双三次插值法对u、v进行重采样插值，放大至与原始全色光图像pan同尺寸，记为和

第七步，对强度分量int以及和进行ihs反变换，得到原始彩色空间三波段图像；与三波段多光谱图像mt1、mt2、mt3对应的ihs反变换输出

第八步，计算空间细节d＝pan-int；利用最近邻插值法对原始多光谱图像的所有多光谱波段m1，…，mn进行重采样插值，放大至与全色光图像pan同尺寸，记为计算各波段的空间细节调制系数利用空间细节调制系数将空间细节d调制后注入三波段多光谱图像st1、st2、st3中，输出图像ft1＝st1+gt1×d、ft2＝st2+gt2×d和ft3＝st3+gt3×d即为融合结果；

第九步，返回第五步选择不同的三波段进行融合，最终遍历所有n个波段，输出n个波段的融合结果。

本发明的有益效果是：针对遥感图像融合技术在重采样多光谱图像过程中极易引入空间和光谱信息失真的问题，通过利用导向滤波器的系数矩阵将全色光图像的空间结构细节导入放大的多光谱图像中，提高了多光谱图像的空间分辨率；再利用ihs变换分离多光谱图像的光谱信息，利用ihs反变换将光谱信息注入至多光谱图像中，有助于保持原始的光谱信息；最后通过细节的提取调制与注入，在保持多光谱图像光谱信息的基础上进一步提升了空间分辨率，得到高质量的融合图像，是一种适用于高分辨率星载多光谱与全色光图像融合的有效融合方法。

附图说明

图1是本发明的原理示意图；

图2是本发明的流程图。

具体实施方式

下面结合附图和实施例对本发明进一步说明，本发明包括但不仅限于下述实施例。

本发明包括以下步骤：

假设原始多光谱图像包含n个波段(m1，m2，…,mn)，各波段图像尺寸均为r1行×c1列，原始全色光图像尺寸为r2行×c2列。

第一步、下采样全色光图像：

通过卷积运算对全色光图像(pan)进行低通滤波，低通滤波器记为lpf，输出图像记为pl，*为卷积运算符

pl＝pan*lpf(1)

对pl作下2抽样处理，即对像素点隔行取样和隔列取样，输出结果记为p2，其尺寸减小至原始全色光图像尺寸的一半(r2/2行×c2/2列)

p2＝(pl)↓2(2)

第二步、导向滤波

设置图像p2为导向滤波的指导图像i，任意第j个多光谱波段(mj)为导向滤波的输入图像p，导向滤波输出图像记为q。按照导向滤波器的定义(导向滤波器原理参见ieeetransactionsonpatternanalysisandmachineintelligence,35(6):1397-1409,2013.)，假设滤波输出图像q与指导图像i之间存在局部线性关系，即

其中像素点i位于以像素点k为中心的邻域(或窗口)ωk中，该窗口半径为r，窗口大小为(2r+1)×(2r+1)像素(通常取半径r≥1)。ak和bk为窗口系数，假设在窗口ωk内均为常数。由此可知，滤波输出图像q在像素点i的滤波输出值qi由指导图像i在该点的像素值ii以及i点所在局部窗口的窗口系数ak和bk共同决定。窗口ωk内的系数ak和bk分别由下式计算：

其中|ω|为窗口ωk中所有像素的个数，pi为输入图像p在i处的像素值，μk和σk分别是指导图像在窗口ωk中的均值和方差，是输入图像p在窗口ωk中的均值，ε是规则化参数且ε∈(0,1)。

由于包含像素点i的邻域ωk不唯一，同时存在多个，因此可采用均值策略计算输出图像q在像素点i的像素值qi，即

因为正方形窗口的对称性，故可以将公式(6)重写成

其中

将公式(7)表示成矩阵形式，可得

q＝a×i+b(10)

其中a和b是窗口系数矩阵。

第三步、窗口系数的重采样放大

利用最近邻插值法对公式(10)的窗口系数矩阵a和b重采样放大，分别得到对应的新窗口系数和

第四步、导向滤波输出

利用原始全色光图像pan和新窗口系数和计算大尺寸的输出图像

该输出图像与原始全色光图像尺寸相同，尺寸为r2行×c2列。

第五步、计算强度分量int：

从n波段多光谱图像中任选三波段(假设记为mt1、mt2、mt3)进行彩色合成，将该三波段分别作为导向滤波的输入图像，并利用第二步至第四步分别计算出对应的三幅大尺寸输出图像，分别记为和取平均计算出强度分量int

第六步、多光谱图像的intensity-hue-saturation(ihs)变换

对所选的三波段多光谱图像mt1、mt2、mt3进行ihs变换，具体如下

利用双三次插值法对u、v进行重采样插值，放大至与原始全色光图像pan同尺寸，记为和

第七步、ihs反变换

对强度分量int以及和进行ihs反变换，得到原始彩色空间三波段图像。与三波段多光谱图像mt1、mt2、mt3对应的ihs反变换输出分别记为st1、st2、st3

第八步、空间细节的提取、调制与注入

空间细节d直接由全色光图像pan和强度分量int的差值计算

d＝pan-int(18)

利用最近邻插值法对所有多光谱波段(m1，…，mn)进行重采样插值，放大至与全色光图像pan同尺寸，记为

计算各波段的空间细节调制系数gj，j＝1,…,n具体如下

利用调制系数将空间细节d调制后注入至第七步输出的三波段多光谱图像st1、st2、st3中，输出图像分别记为ft1、ft2、ft3，即

ft1＝st1+gt1×d(20)

ft2＝st2+gt2×d(21)

ft3＝st3+gt3×d(22)

图像ft1、ft2、ft3即为输出的融合结果。对于包含n个波段的多光谱图像(m1，m2，…,mn)，可以多次选择其中的三波段进行融合，最终遍历所有n个波段，输出n个波段的融合结果。

方法实施例：

采用真实worldview-2星载遥感多光谱图像和全色光图像，多光谱图像包含八个波段(m1,m2,…,m8)，全色光图像(pan)为单波段。多光谱图像的空间分辨率为2.0m，大小为200行×200列。全色光图像空间分辨率为0.5m，大小为800行×800列。实施本发明包括以下步骤：

第一步、下采样全色光图像

由于全色光图像空间分辨率是多光谱图像空间分辨率的四倍，因此全色光图像的尺寸是多光谱图像的四倍，故需要对全色光图像下采样两次。这里低通滤波器选择使用cohen-daubechies-fauveau(cdf)9/7双正交滤波器组中的低通滤波器cdf9，即lpf＝[0.026748,-0.016864,-0.078223,0.266864,0.602949,0.266846,-0.078223,-0.016864,0.026748]利用公式(1)和(2)完成第一次下采样处理，得到p2，p2的尺寸缩小至400行×400列。再利用公式(1)和(2)对p2完成第二次下采样，即

pl2＝p2*lpf

p4＝(pl2)↓2

输出得到的p4尺寸缩小至200行×200列，与多光谱图像尺寸一致。

第二步、导向滤波

从多光谱图像中选择三个波段(mj，j＝2，3，5)，设置图像p4为导向滤波的指导图像i，mj为输入图像，利用公式(3)-(10)分别对这三个波段进行导向滤波，输出分别记为qj，对应的窗口系数矩阵记为aj和bj，即

qj＝aj×i+bj,(j＝2,3,5)

这里规则化参数ε＝0.01，窗口半径r＝2，窗口大小为5×5。

第三步、窗口系数的重采样放大

利用最近邻插值法对窗口系数矩阵aj和bj重采样放大四倍，分别得到对应的新窗口系数和

第四步、导向滤波输出

利用全色光图像pan和新窗口系数和计算输出图像即

该输出图像与原始全色光图像pan尺寸相同，为800行×800列。

第五步、计算强度分量int：

按公式(12)对输出图像和取平均，计算出强度分量int，即

第六步、多光谱图像的ihs变换

按公式(13)和(14)对所选的三波段多光谱图像m2、m3、m5进行ihs变换，具体如下

利用双三次插值法对u、v进行重采样插值，放大四倍与全色光图像pan同尺寸，记为和

第七步、ihs反变换

按照公式(15)-(17)对强度分量int以及和进行ihs反变换，得到原始彩色空间三波段图像，对应三波段输出分别记为s2、s3、s5，即

第八步、空间细节提取、调制与注入

利用公式(18)提取出空间细节d，即

d＝pan-int

利用最近邻插值法对所有多光谱波段(m1，…，m8)进行重采样插值，放大四倍至与全色光图像pan同尺寸，记为尺寸为800行×800列。

利用公式(19)计算各波段的细节调制系数gj，j＝1,…,8，即

利用公式(20-22)计算输出图像f2、f3和f5即

f2＝s2+g2×d

f3＝s3+g3×d

f5＝s5+g5×d

f2、f3和f5即为输出的融合结果。