一种利用滚动导向滤波的遥感图像融合方法与流程

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

背景技术：

随着星载成像传感器的快速发展，许多卫星可以提供高分辨率的遥感图像。一些在轨的高分辨率商业卫星，例如ikonos、quickbird、worldview-2能同时捕捉全色光图像和多光谱图像。由于物理约束，全色光图像具有非常高的空间分辨率，但是只具有单一光谱通道，与此相反，多光谱图像具有较低的空间分辨率但更高的光谱分辨率。通过融合多光谱图像和全色光图像，可以得到一幅同时具有高空间分辨率和光谱信息丰富的融合图像，其在保持光谱特性的同时，增强了图像的可视化效果。

在多光谱和全色光图像融合方法中，空间细节模型主要思想是将有用的空间细节信息从多光谱图像中提取出来，然后用预先设定的注入规则加在多光谱图像中。在该类模型中边缘保持滤波器因其边缘感知能力成为研究热点，近几年在研究中使用的该类滤波器有非局部均匀滤波、加权最小二乘滤波、导向滤波(参见文献ieeecomputersocietyconferenceoncomputervisionandpatternrecognition,2(7):60-65,2005；acmtransactionsongraphics,27(3):15-19,2008；ieeetransactionsonpatternanalysisandmachineintelligence,35(6):1397-1409,2013)。然而边缘保持滤波器很难分离结构与细节，得到的融合结果光谱保真度与空间信息都有一定的欠缺。

技术实现要素：

为了克服现有技术的不足，本发明提供一种基于滚动导向滤波的遥感图像融合方法，不同于边缘保持滤波器，滚动导向滤波器能够提取不同尺度结构而保留其他的信息。本发明通过滚动导向滤波将全色光和多光谱图像的空间信息准确的提取出来，再使用gramschmidt(gs)变换将空间信息注入到多光谱图像中。

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

第一步，对于原始多光谱图像包含的n个波段图像ms1、ms2、…、msn，通过双三次插值将各波段图像上采样到和全色光图像pan尺寸相同的图像rms1、rms2、…、rmsn；

第二步，将rms1、rms2、…、rmsn波段图像对应像素点的灰度值取平均作为低分辨率全色光图像pl对应像素点灰度值，即将pl与rms1、rms2、…、rmsn波段图像作为gs正变换的输入图像得到变换后的波段图像gs1，gs2，…，gsn+1；

第三步，令亮度分量图像i等于pl，分别计算图像i与pan所有像素点灰度值的均值与标准差，将图像i与pan做直方图匹配得到式中μi与μp分别为亮度分量图像i和全色光图像pan所有像素点灰度值的均值，σi与σp分别为亮度分量图像i和全色光图像pan的所有像素点灰度值的标准差；

第四步，构建滚动导向滤波器，并将滚动导向滤波算法记为rgf(·)，分两个阶段进行：

第一阶段对输入图像做高斯滤波，该阶段输入图像为q，输出图像记为g，其中σs为q的所有像素点的灰度值标准偏差，n(p)为以p为中心、w×w窗口大小的像素集，q为n(p)中不同于p的像素点，||p-q||表示两个像素点的欧式距离；

第二阶段采用导向滤波算法gf(·)，该阶段输入图像为q，与第一阶段输入图像相同，图像j作为迭代更新的指导图像，其中第一阶段的输出g记为j¹，第t次迭代的结果j^t+¹＝gf(j^t，q)；

第五步，将亮度分量图像i与直方图匹配全色光图像pan'分别作为滚动导向滤波器的输入图像，计算出对应的输出图像，即li＝i-rgf(i)和lp＝pan'-rgf(pan')，再计算i与pan'的细节图像d1＝i-li与d2＝pan'-lp；

第六步，构造能量显著映射sj＝|box(dj)|,j＝1,2，其中box(·)为r×r大小的平均滤波器，|·|为取绝对值操作，dj为第j幅细节图像；再利用能量显著映射构造加权映射其中sj(m)表示第j幅图像中第m个像素的显著映射值；

将dj作为导向滤波器的指导图像，wj作为输入图像，计算出输出图像oj＝gf(dj,wj)，最后计算细节图像

第七步，将gs变换后得到的第一波段图像gs1作为滚动导向滤波的输入图像，滚动导向滤波器的输出图像记为lgs1，即lgs1＝rgf(gs1)。将细节图像d注入lgs1作为新的gs变换的第一波段图像gs1new＝d+lgs1；

第八步，对新的变换波段gs1new，gs2，…，gsn+1进行gs反变换，得到n个波段融合结果图像f1，f2，…，fn。

本发明的有益效果是：针对多光谱与全色光图像融合过程中由于细节提取与注入模型不适导致融合结果产生光谱分辨率与空间分辨率失真问题，采用高斯滤波与导向滤波相结合的滚动导向滤波器，该滤波器最主要的优势是其尺度感知性能，其可将输入图像分解为两层：一层是包含大尺度结构的基础层，另一层是包含小尺度纹理的细节层。在发明中利用滚动导向滤波器将多光谱与全色光图像不同尺度的结构分离，从而可以提取出所需的空间信息，并结合gramschmidt(gs)变换将细节注入到多光谱图像中得到融合结果。融合结果不仅提高了原始多光谱图像的空间分辨率，也保持了光谱信息，是一种适合于高分辨率星载多光谱与全色光图像融合的有效方法。

附图说明

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

图2是本发明的流程图。

具体实施方式

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

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

第一步、上采样多光谱图像：

假设原始多光谱图像包含n个波段(ms1，ms2，…，msn)，通过双三次插值将各波段图像上采样到和全色光图像尺寸相同的图像(rms1，rms2，…，rmsn)。

第二步、gs正变换：

将n个波段(rms1，rms2，…，rmsn)对应像素点的灰度值取平均设为低分辨率全色光图像pl对应像素点灰度值，将pl与这n个波段(rms1，rms2，…，rmsn)作为gs正变换(参见文献u.s.patent6011875,2000)的输入图像得到变换后的波段(gs1，gs2，…，gsn+1)。

第三步、直方图匹配：

令亮度分量图像i等于pl，分别计算图像i与pan所有像素点灰度值的均值与标准差，将图像i与pan做直方图匹配得到pan'：

i＝pl(2)

式中μi与μp分别为亮度分量图像i和全色光图像pan的所有像素点灰度值的均值，σi与σp分别为亮度分量图像i和全色光图像pan的所有像素点灰度值的标准差。

第四步、滚动导向滤波：

滚动导向滤波器原理参见2014europeanconferenceoncomputervision(eccv2014),8691:815-830,2014。

设图像q为滚动导向滤波的输入图像，并将滚动导向滤波算法记为rgf(·)，该滤波器分两个阶段进行：

第一阶段使用高斯滤波器对输入图像做高斯滤波，该阶段输出图像记为g，具体公式如下：

其中σs为q的所有像素点的灰度值标准偏差，n(p)为以p为中心、w×w窗口大小的像素集，q为n(p)中不同于p的像素点，||p-q||表示两个像素点的欧式距离。

第二阶段使用联合滤波器进行迭代操作，联合滤波器可采用联合双边滤波器，导向滤波器等。在该发明中采用导向滤波器，该阶段输入图像为q，与第一阶段输入图像相同，图像j作为迭代更新的指导图像，其中第一阶段的输出g记为j¹，j^t+1为第t次迭代的结果。将gf(·)记为导向滤波算法，则该阶段公式如下：

j^t+1＝gf(j^t，q)(5)

第五步、细节的提取：

将多光谱图像的亮度分量图像i与直方图匹配全色光图像pan'分别作为滚动导向滤波器的输入图像，并利用第四步分别计算出对应的输出图像，分别记为li和lp，再计算i与pan'的细节图像d1与d2：

li＝i-rgf(i)(6)

lp＝pan'-rgf(pan')(7)

d1＝i-li(8)

d2＝pan'-lp(9)

第六步、细节图像融合：

将细节图像d1、d2通过低通滤波器并取绝对值，构造能量显著映射：

sj＝|box(dj)|,j＝1,2(10)

其中box(·)为r×r大小的平均滤波器，dj为第j幅细节图像，|·|为取绝对值操作，得到的显著映射提供了图像中较好的细节特征。再利用能量显著映射构造加权映射：

其中sj(m)表示第j幅图像中第m个像素的显著映射值。

将dj作为导向滤波器的指导图像，wj作为输入图像，计算出输出图像oj，最后计算细节图像d：

oj＝gf(dj,wj),j＝1,2(12)

第七步、细节注入：

将gs变换后得到的第一波段gs1作为滚动导向滤波的输入图像，利用第四步计算出输出图像记为lgs1，将细节图像d注入lgs1作为新的gs变换的第一波段，并记为gs1new：

lgs1＝rgf(gs1)(14)

gs1new＝d+lgs1(15)

第八步、gs反变换：

对新的变换波段(gs1new，gs2，…，gsn+1)进行gs反变换，得到n个波段融合结果图像(f1，f2，…，fn)。

方法实施例：

采用真实ikonos星载遥感多光谱和全色光图像，多光谱图像包含四个波段(ms1-ms4)，全色光图像(pan)为单波段。多光谱图像的空间分辨率为4.0m，大小为100行×100列。全色光图像空间分辨率为1.0m，大小为400行×400列。实施本发明包括以下步骤：

第一步、上采样多光谱图像：

由于全色光图像空间分辨率是多光谱图像空间分辨率的4倍，因此通过双三次插值将4个波段(ms1，ms2，ms3，ms4)上采样到和全色光图像尺寸相同的图像(rms1，rms2，rms3，rms4)。

第二步、gs正变换：

将4个波段(rms1，rms2，rms3，rms4)对应像素点的灰度值取平均设为低分辨率全色光图像pl对应像素点灰度值，将pl与这4个波段作为gs正变换的输入图像得到变换后的波段(gs1，gs2，…，gs5)。

pl＝(rms1+rms2+rms3+rms4)/4

第三步、直方图匹配：

利用公式(2)得到亮度分量图像i，再利用公式(3)将图像i与pan做直方图匹配得到pan'。

第四步、滚动导向滤波：

设置图像i与pan'分别为滚动导向滤波的输入图像，利用公式(4)和(5)对这两幅图像进行滚动导向滤波，其中σs、w取值均为2，t取值为4。根据公式(6)和(7)输出图像分别记为li与lp。

第五步、细节的提取：

利用公式(8)和(9)计算i与pan'的细节图像记为d1与d2。

第六步、细节图像融合：

将细节图像d1与d2通过3×3的平均滤波器并取绝对值，构造能量显著映射：

s1＝|box(d1)|

s2＝|box(d2)|

再利用s1与s2构造加权映射w1与w2。

将d1与d2分别作为导向滤波器的指导图像，相应的加权映射w1与w2分别作为输入图像，利用导向滤波算法gf(·)计算出输出图像o1与o2。

o1＝gf(d1,w1)

o2＝gf(d2,w2)

最后根据公式(13)计算融合的细节图像d。

第七步、细节注入：

将gs正变换后得到的第一波段gs1作为滚动导向滤波的输入图像，利用公式(4)和(5)计算出lgs1，其中σs、w取值均为2，t取值为4。根据公式(15)将细节图像d注入lgs1作为新的gs变换的第一波段gs1new，即：

gs1new＝d+lgs1

第八步、gs反变换：

对新的变换波段(gs1new，gs2，…，gs5)进行gs反变换，得到4个波段融合结果图像(f1，f2，f3，f4)。