一种星载全色与多光谱图像融合方法与流程

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

背景技术：

由于成像传感器自身的物理局限性以及星地间数据传输技术的壁垒，目前的星载成像传感器仅能提供光谱分辨率高但空间分辨率低的多波段多光谱图像或者空间分辨率高但光谱信息贫乏的单波段全色图像。在轨的星载传感器，例如geoeye-1、worldview-2/3/4、高景1号等能够同时提供多波段多光谱遥感图像和单波段高空间分辨率的全色图像。通过融合多光谱图像和全色图像，可以在保持多光谱图像光谱信息的同时提高其空间分辨率，融合结果可服务于高精度地物分类、城市规划、地理信息系统等多个遥感应用。

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

技术实现要素：

为了克服现有技术的不足，本发明针对星载多光谱图像和全色图像提出一种图像融合方法，首先利用高斯导向滤波(gaussianguidedfiltering)联合全色图像完成对原始多光谱图像的重采样放大，通过该联合上采样过程将全色图像的空间细节信息注入至放大的多光谱图像中，再对原始多光谱图像使用超球彩色变换(hypersphericalcolortransform,hct)提取其光谱信息，在超球彩色空间中将光谱信息和全色图像的空间信息融合，经hct反变换得到具有高空间分辨率的多光谱图像。

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

第一步、对于包含n个波段m1,m2,…,mn的原始多光谱图像，对应的全色图像记为ph，通过卷积运算对全色图像ph进行低通滤波，低通滤波器记为lpf，输出图像pl＝ph*lpf；对pl的像素点进行隔行取样和隔列取样，输出结果pl＝(pl)↓2；

第二步、设置图像pl为高斯导向滤波的指导图像，记为i，任意第j个多光谱波段mj为高斯导向滤波的输入图像，记为p，高斯导向滤波输出图像记为q；假设滤波输出图像q与指导图像i之间存在局部线性关系，输出图像q在像素点i的像素值其中像素点i位于以像素点k为中心的窗口wk中，|w|为窗口wk内所有像素的个数，该窗口半径为r，窗口大小为(2r+1)×(2r+1)像素，窗口wk内的高斯权系数定义为xi为窗口内像素点i的像素值，xk为窗口内像素点k的像素值，ak和bk为窗口系数，

其中pi为输入图像p在i点的像素值，ε是规则化参数且ε∈(0,1)；

将q表示成矩阵形式，得到q＝a×i+b，其中a和b是窗口系数矩阵；

第三步、利用最近邻插值法对a和b分别做插值放大，放大至与全色图像同尺寸，得到对应的新窗口系数矩阵和

第四步、利用全色图像ph和和计算出大尺寸的输出图像

重复第二步至第四步对原始多光谱图像的各波段进行滤波放大，得到对应的滤波输出图像

第五步、利用超球彩色变换到超球彩色空间，得到一个强度分量im和n-1个角度分量θ1，θ2，…，θn-1，

利用双线性插值法对n-1个角度分量θ1,θ2,…,θn-1分别做插值放大，放大至与全色图像同尺寸，分别记为构建新强度分量

第六步、对第五步得到的n-1个角度分量和新强度分量inew做hct反变换，得到的n波段输出图像记为f1,f2,…,fn；

第七步、计算高斯导向滤波输出图像与全色图像ph的均值与方差，得到匹配后的全色图像j＝1,2,…,n，其中μp与μj分别为全色图像ph和图像的均值，σp与σj分别为全色图像ph和图像的标准差；

计算各波段的待注入空间细节

计算第j个波段的融合结果sj＝fj+dj，得到所有波段的融合结果(s1,s2,…,sn)。

本发明的有益效果是：先后包含了两个过程去提高多光谱图像的空间分辨率，首先使用高斯导向滤波联合全色图像，将全色图像的空间细节信息部分注入至重采样放大的多光谱图像中，提高了重采样多光谱图像的空间分辨率；然后利用空间细节的再注入(第八步)过程进一步提高hct反变换输出图像的空间分辨率。此外，本发明利用hct变换有效地分离出多光谱图像的原始光谱信息，经hct反变换将光谱信息融入至重采样的多光谱图像中，有利于提高融合结果的光谱保真能力。本发明与专利：一种基于导向滤波的遥感图像融合方法(申请号：2017102421174，公开号：cn107146212a)主要的区别之处在于：1)本发明使用的高斯导向滤波不同于传统的导向滤波，能够克服传统导向滤波在滤波窗口输出沿水平与垂直方向产生的偏执、不均匀现象；2)ihs(intensity-hue-saturation)变换每次只能输入三个波段的多光谱图像，而本发明使用hct变换可以突破对输入多光谱图像波段数量的限制，可以直接拓展至高光谱图像的融合处理；3)本发明不受全色图像与多光谱图像光谱响应匹配关系的影响，在对各波段注入空间细节时，不需要对待注入的空间细节进行专门的调制，普遍性和实用性更强。

附图说明

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

图2是本发明的流程图。

具体实施方式

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

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

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

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

pl＝ph*lpf(1)

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

pl＝(pl)↓2(2)

假定全色图像的空间分辨率是原始多光谱图像空间分辨率的2倍，即r2/2＝r1，c2/2＝c1，公式(2)输出的下采样全色图像pl与原始多光谱图像同尺寸。

第二步、高斯导向滤波

设置图像pl为高斯导向滤波的指导图像，记为i，任意第j个多光谱波段(mj)为高斯导向滤波的输入图像，记为p，高斯导向滤波输出图像记为q。假设滤波输出图像q与指导图像i之间存在局部线性关系，即

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

其中pi为输入图像p在i点的像素值，ε是规则化参数且ε∈(0,1)。

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

|w|为窗口wk内所有像素的个数。因为正方形窗口的对称性，故可以将公式(6)重写成

其中

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

q＝a×i+b(10)

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

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

利用最近邻插值法对公式(10)的窗口系数矩阵a和b分别做插值放大，放大至与全色图像同尺寸，得到对应的新窗口系数矩阵，分别记为和

第四步、高斯导向滤波输出

利用全色图像ph和新窗口系数和计算出大尺寸的输出图像，记为

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

利用上述第二步、第三步和第四步对原始多光谱图像(m1,m2,…,mn)的各波段进行滤波放大，对应的滤波输出图像可分别记为

第五步、原始多光谱图像的光谱信息提取

利用超球彩色变换(hypersphericalcolortransform(hct)，参见文献asprs2010annualconference,sandiego,california,pp.26-30,april2010)将原始多光谱图像(m1,m2,…,mn)变换到超球彩色空间，得到一个强度分量(记为im)和n-1个角度分量，记为θ1，θ2，…，θn-1，hct变换过程如下：

利用双线性插值法对公式(13)得到的n-1个角度分量(θ1,θ2,…,θn-1)分别做插值放大，放大至与全色图像同尺寸，分别记为按公式(14)构建新强度分量inew：

第六步、hct反变换

对第五步得到的n-1个角度分量和新强度分量inew按公式(15)做hct反变换，得到的n波段输出图像记为f1,f2,…,fn。

第七步、空间细节再注入

hct反变换输出的多光谱图像(f1,f2,…,fn)已经具有丰富的光谱信息和较高的空间分辨率，为了进一步提高空间分辨率本步骤设计了空间细节的再注入。分别计算高斯导向滤波输出图像与全色图像ph的均值与方差，利用公式(16)得到匹配后的全色图像，记为pmj，j＝1,2,…,n。

其中μp与μj分别为全色图像ph和图像的均值，σp与σj分别为全色图像ph和图像的标准差。

计算各波段的待注入空间细节dj：

各波段的空间细节再注入按公式(18)计算得到

sj＝fj+dj(18)

sj即为第j个波段的融合结果，j＝1,2,…,n。由上述计算过程，可以得到所有波段的融合结果(s1,s2,…,sn)。

实施例：采用真实worldview-2星载遥感多光谱图像和全色图像，多光谱图像包含八个波段，记为m1,m2,…,m8，全色图像为单波段，记为ph。多光谱图像的空间分辨率为2.0米，大小为200行×200列。全色图像空间分辨率为0.5米，大小为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)完成第一次下采样处理，得到pl1，pl1的尺寸缩小至400行×400列。

pl＝ph*lpf

pl1＝(pl)↓2

再利用公式(1)和(2)对pl1完成第二次下采样，即

pl2＝pl1*lpf

pl2＝(pl2)↓2

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

第二步、高斯导向滤波

设置图像pl2为高斯导向滤波的指导图像i，任意第j个多光谱波段mj为高斯导向滤波的输入图像，利用公式(3)-(10)对mj做高斯导向滤波，滤波输出记为qj，对应的窗口系数矩阵记为aj和bj，即

qj＝aj×i+bj,(j＝1,2,...,8)

这里规则化参数ε＝10^-8，窗口半径r＝2，即窗口大小为5×5。

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

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

第四步、高斯导向滤波输出

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

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

重复步骤二、三、四对所有8波段多光谱图像进行高斯导向滤波处理，依次得到对应的8波段重采样放大图像，记为

第五步、原始多光谱图像的光谱信息提取

利用公式(13)对原始8波段多光谱图像(m1,m2,…,m8)作hct变换，得到7个角度分量，记为θ1，θ2，…，θ7

利用双线性插值法对7个角度分量(θ1，θ2，…，θ7)分别做插值放大，放大四倍至与全色图像同尺寸(800行×800列)，分别记为按公式(14)构建新强度分量inew：

第六步、hct反变换

对第五步得到的7个角度分量和新强度分量inew按公式(15)做hct反变换，得到的8波段输出图像记为f1,f2,…,f8。

第七步、空间细节再注入

hct反变换输出的8波段多光谱图像(f1,f2,…,f8)已经具有丰富的光谱信息和较高的空间分辨率，为了进一步提高空间分辨率本步骤设计了空间细节的再注入。分别计算高斯导向滤波输出图像与全色图像ph的均值与方差，利用公式(16)得到匹配后的全色图像，记为pmj，j＝1,2,…,8。

其中μp与μj分别为全色图像ph和图像的均值，σp与σj分别为全色图像ph和图像的标准差。

由公式(17)计算各波段的待注入空间细节dj：

各波段的空间细节再注入按公式(18)计算得到

sj＝fj+dj

sj即为第j个波段的融合结果，j＝1,2,…,8。由上述计算过程，可以得到所有8波段图像的融合结果(s1,s2,…,s8)。