一种基于加权平均曲率滤波分解的遥感图像融合方法

1.本发明涉及多光谱和全色遥感图像融合技术领域，尤其是涉及一种基于加权平均曲率滤波分解的遥感图像融合方法。


背景技术：

2.随着遥感技术的快速发展，可以获得很多具有不同分辨率的遥感图像。多光谱(multispectral，ms)图像包含丰富的光谱信息，在灾害监测、海洋研究、目标识别等领域得到广泛的应用。全色(panchromatic，pan)图像具有高空间分辨率和低光谱分辨率。单一传感器获取的遥感图像存在信息不完全的缺点。因此，ms和pan图像融合技术对于获得高空间和光谱分辨率的融合图像至关重要。ms和pan图像融合又叫全色锐化。全色锐化技术主要分为两类：基于分量替换(component substitution，cs)的算法和基于多尺度分解(multi-scale decomposition，msd)的算法。对于基于cs的算法，ms图像被变换到新的投影空间，分解为光谱成分和空间成分，然后用全色图像替换空间分量。最后，进行相应的逆变换重建融合后的图像。基于cs的经典算法包括强度色调饱和(intensity hue saturation，ihs)变换、主成分分析、gram-schmidt(gs)等。其可以有效提升空间分辨率，但可能会出现一定程度的光谱失真。常见的msd算法包括：小波变换、曲线变换、轮廓波变换和剪切波变换(non-subsampled shearlet transform，nsst)。对于基于msd的算法，源图像通过多尺度分解方法分解为不同尺度的子图像。在不同尺度的子图像上根据不同的融合规则对子图像进行融合。最后通过逆变换重建融合后的图像。基于msd的算法可以获得更好的光谱分辨率，但空间分辨率低于cs算法。
3.近年来，边缘保持滤波(edge-preserving filter，epf)在图像处理中得到了广泛的应用。常见的epf算法包括双边滤波、引导图像滤波和高斯曲率滤波。epf可以在平滑图像的同时保留图像的边缘信息，充分利用图像的空间信息。基于这一优越特性，epf经常被用作图像融合中的一种图像分解方法。另一种重要的算法是基于脉冲耦合神经网络(pulse coupled neural network,pcnn)。pcnn经过多次迭代可以更准确地提取图像特征，更有利于图像融合。


技术实现要素：

4.本发明的目的是提供一种基于加权平均曲率滤波分解的遥感图像融合方法，解决了融合图像空间失真和光谱失真问题，通过与十种典型的融合方法进行比较，在主观视觉效果和客观评价方面都取得较好的结果。
5.为实现上述目的，本发明提供了一种基于加权平均曲率滤波分解的遥感图像融合方法，包括以下步骤：
6.s1、计算多光谱强度分量：采用基于局部空间频率的自适应加权平均方法，用于将多光谱图像各波段融合生成强度分量i，加权局部空间频率利用水平、垂直以及主副对角线方向的八邻域内的多个像素点，依据距离中心像素的欧氏距离赋予不同的权重，加权之后
再参与计算，根据加权局部空间频率设计自适应加权平均的系数ωi，具体如下所示：
[0007][0008][0009]
其中n是多光谱图像的波段，wlsfi表示多光谱图像第i个波段在(x，y)位置处的加权局部空间频率值，msi(x,y)表示多光谱图像第i个波段在(x，y)位置的像素值，ωi(x,y)表示多光谱图像第i个波段在(x，y)位置的加权系数，i(x,y)表示强度分量i在(x，y)位置处的像素值；
[0010]
s2、光谱估计：将i作为初始α通道，根据下列公式计算前景图像f和背景图像b，后续利用f、b进行图像重构：
[0011][0012]
其中i是第i颜色通道，α
ix
和α
iy
值是光谱前景fk、光谱背景bk和α的水平和垂直导数；
[0013]
s3、多尺度分解：利用加权平均曲率滤波(wmcf)和高斯滤波(gf)对全色图像进行多尺度分解，得到大尺度图像l、小尺度图像s和基本图像b，通过这种分解方法，将输入图像分解成三个不同尺度的子图像，分解过程为利用wmcf和gf对输入图像m进行处理，分别获得滤波图像iw和ig，通过m和iw的差值得到小尺度图像，大尺度图像是差分结果：
[0014]
ig＝gf(m,σ,ρ).
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(4)
[0015]iw
＝wmcf(m,p).
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(5)
[0016]
s＝m-iw.
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(6)
[0017]
l＝i
w-ig.
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(7)
[0018]
b＝ig.
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(8)；
[0019]
s4、分量融合：利用图像的局部能量e来动态调节连接系数β的大小，根据图像的特征自适应变化保留图像边缘的特征信息，自适应连接系数β描述为：
[0020][0021]
其中，e(x,y)表示(x,y)处的局部能量，β
xy
表示(x,y)处的连接系数，局部能量的计算具体描述如下：
[0022][0023]
其中c(x,y)表示图像在(x，y)处的像素值，w为加权系数矩阵，具体细节如下：
[0024][0025]
将加权局部空间频率值作为脉冲耦合神经网络的输入，利用图像的局部能量e来
动态调节连接系数β的大小，根据图像的特征自适应变化，当达到最大迭代次数时，迭代停止，得到强度分量i的点火次数总和tb和基本图像b的点火次数总和，根据点火次数总和取大原则进行融合，如下式所示：
[0026][0027]
其中fa(x,y)代表融合图像fa在(x,y)位置的像素值，i(x,y)代表强度分量i在(x,y)位置的像素值，b(x,y)代表基本图像b在(x,y)位置的像素值，最终，通过加权局部空间频率引导的脉冲耦合神经网络融合基本图像b和强度分量i得到融合图像fa；
[0028]
s5、图像重构：通过s、l和fa的线性组合重建融合后的图像fb作为最终的α通道，根据如下公式，组合α、f和b进行重构，得到最终的融合结果：
[0029]ii
＝α
ifi
(1-αi)biꢀꢀꢀꢀꢀꢀꢀꢀ
(13)。
[0030]
优选的，所述步骤s2根据图像抠图模型，通过线性合成模型将输入图像区分为前景图像f和背景图像b，即第i个像素的颜色是对应的前景颜色和背景颜色的线性组合：
[0031]ii
＝α
ifi
(1-αi)biꢀꢀꢀꢀꢀꢀꢀꢀ
(14)
[0032]
其中fi是第i个像素的前景颜色，bi是第i个像素的背景颜色，其中，α是前景图像f的不透明度，根据图像抠图模型的局部线性假设，如果α包含该图像的大部分边缘信息，则该图像的对应的光谱前景和背景色在空间上是平滑的，在确定输入图像和α的同时，通过求解以下函数估计前景图像f和背景图像b：
[0033][0034]
其中i是第i颜色通道，α
ix
和α
iy
值是光谱前景fj、光谱背景bj和α的水平和垂直导数。
[0035]
优选的，所述步骤s4中脉冲耦合神经网络模型具体为：
[0036]
标准的pcnn模型由以下公式描述：
[0037][0038][0039]uij
[n]＝f
ij
[n](1 βl
ij
[n]).
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(18)
[0040][0041][0042]
t
ij
[n]用于阈值调节，u
ij
[n]用于产生调制信号与阈值进行比较产生脉冲信号y
ij
[n]，共有9个参数，9个参数分属四类：时间衰减常数、放大系数、连接系数以及权值矩阵，时间衰减常数为αf,α
l
,α
θ
，这三个参数分别来自输入域f、连接输入域l以及动态门限t，αf,α
l
控制着f，l通道信号的衰减速度，α
θ
控制着阈值衰减的速度，放大系数为vf、v
l
、ve，这三个参数分别是输入域f、连接输入域l以及动态门限t的放大系数，连接系数β来自调制域u起到调节神经元之间相互影响的作用，同时该参数还影响着中心神经元的点火周期，权值矩阵为w、
m，这两个参数分别来自输入域f、连接输入域l，矩阵中的元素表示神经元之间的连接强弱关系；
[0043]
接收域f接收外部刺激s、自身的刺激f
ij
[n-1]，l接收周围神经元的刺激，调制域针对接收域的信号f和l进行调制得到内部活动项u，将脉冲产生域神经元对应位置的内部活动项u与动态阈值t
ij
进行比较，若u
ij
＞t
ij
，则脉冲发生器启动神经元点火输出y
ij
＝1，否则神经元不点火，输出y
ij
＝0，得到最后结果y，上述过程不断重复，直到满足条件为止。
[0044]
优选的，所述步骤s1和s4中加权局部空间频率具体为：
[0045]
传统的空间频率描述水平和垂直两个方向的空间信息，传统的空间频率计算如下式：
[0046][0047]
其中，rf、cf分别是行频率、列频率，rf、cf包含水平、垂直两个方向；
[0048]
对空间频率的计算进行改进，在原有基础上增加对角方向的计算，此外，基于欧几里得度量，确定行、列频率的加权系数为改进后的空间频率称为加权局部空间频率(wlsf)，用于反映各像素邻域范围内的活跃程度，wlsf同时包含行、列和对角方向，对八个方向的空间频率进行加权求和，加权局部空间频率反映各像素邻域范围内的活跃程度，加权局部空间频率的计算公式如下：
[0049][0050]
其中，rf、cf和df分别是行频率、列频率和对角频率，rf、cf包含水平、垂直两个方向，df包含正副对角线方向，具体定义如下：
[0051]
rf＝rf1 rf2,cf＝cf1 cf2ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(23)
[0052]
df＝df1 df2 df3 df4ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(24)
[0053][0054][0055][0056][0057][0058]
[0059][0060][0061]
其中，(2m 1)(2n 1)代表局部区域的大小，f(x,y)代表(x,y)位置处的像素值。
[0062]
优选的，所述步骤s3中加权平均曲率具体为：
[0063]
图片u的加权平均曲率(wmc)定义为：
[0064][0065]
对于二维图像即n＝2，该等式变为：
[0066][0067]
其中δ是各向同性拉普拉斯算子，u
x
和uy是x和y方向的偏导数，u
xy
是相应的二阶偏导数；
[0068]
对于3
×
3窗口，考虑8个法线方向和半窗口，8种情况产生相应的内核：
[0069][0070]
根据8个内核计算出八个距离：
[0071][0072]
其中*表示卷积操作，加权平均曲率的离散形式表示为：
[0073]hw
≈dk,k＝argmin{|di|；i＝1,2,...,8}.
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(37)
[0074]
wmcf的简化过程表示为：
[0075]iout
＝wmc(i
in
).
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(38)
[0076]
其中i
in
和i
out
是输入和滤波后的图像，wmc(
·
)表示wmcf操作。
[0077]
因此，本发明采用上述一种基于加权平均曲率滤波分解的遥感图像融合方法，具备以下有益效果：
[0078]
(1)引入图像抠图模型有效保持多光谱遥感图像的光谱分辨率；
[0079]
(2)将加权平均曲率滤波应用于遥感图像融合领域，具有包括尺度和对比度鲁棒
性、快速计算和边缘保持的优点；
[0080]
(3)将加权局部空间频率作为脉冲耦合神经网络的输入，利用图像的局部能量动态调节连接系数的大小，根据图像的特征自适应变化，通过加权局部空间频率引导的脉冲耦合神经网络有效融合子图像；
[0081]
(4)本发明提出的遥感图像融合方法在大多数情况下都取得了良好的性能。
[0082]
下面通过附图和实施例，对本发明的技术方案做进一步的详细描述。
附图说明
[0083]
图1为基于wmcf和gf的分解方法流程图；
[0084]
图2为本发明方法流程图；
[0085]
图3为pcnn模型的基本结构图；
[0086]
图4为pcnn神经元的连接模型；
[0087]
图5为δu,h(u)和hwu之间的关系图；
[0088]
图6为三组ms和pan图像；
[0089]
图7为十种不同方法在第一组图像上的融合结果：(1)bt；(2)gsa；(3)gf；(4)ihs；(5)ct；(6)pca；(7)pracs；(8)nsst-sr；(9)bdsd-pc；(10)wt；(11)本发明的方法；(12)参考的多光谱图像；
[0090]
图8为十种不同方法在第二组图像上的融合结果：(1)bt；(2)gsa；(3)gf；(4)ihs；(5)ct；(6)pca；(7)pracs；(8)nsst-sr；(9)bdsd-pc；10)wt；(11)本发明的方法；(12)参考的多光谱图像；
[0091]
图9为十种不同方法在第三组图像上的融合结果：(1)bt；(2)gsa；(3)gf；(4)ihs；(5)ct；(6)pca；(7)pracs；(8)nsst-sr；(9)bdsd-pc；(10)wt；(11)本发明的方法；(12)参考的多光谱图像。
具体实施方式
[0092]
本发明提供了一种基于加权平均曲率滤波分解的遥感图像融合方法，具体的融合过程包括如下步骤：
[0093]
(1)计算多光谱强度分量：
[0094]
如果采用简单的平均规则，会丢失原图像的一些细节、纹理等重要信息。所以，有必要根据图像中不同区域包含信息的重要程度进行加权，然后再进行融合。加权系数的准确选择决定了融合图像的质量，本发明采用基于加权局部空间频率(weighted local spatial frequency，wlsf)的自适应加权平均方法，用于将多光谱图像各波段融合生成强度分量i。
[0095]
图像的邻近像素之间存在较大的相关性，任何单个像素都无法准确表达局部的区域特征。加权局部空间频率可以充分利用局部区域内包括水平、垂直以及主副对角线方向的八邻域内的多个像素点。依据距离中心像素的欧氏距离赋予不同的权重，加权之后再参与计算。加权局部空间频率可以作为空间域中细节、纹理等重要信息的评价指标。加权局部空间频率值较大的像素被认为是权重较大的信息，如细节、纹理等，在融合过程中赋予较大的权重。因此，本发明根据加权局部空间频率设计自适应加权平均的系数ωi，具体如下所
示。
[0096][0097][0098]
其中n是多光谱图像的波段，wlsfi表示多光谱图像第i个波段在(x，y)位置处的加权局部空间频率值。msi(x,y)表示多光谱图像第i个波段在(x，y)位置的像素值。ωi(x,y)表示多光谱图像第i个波段在(x，y)位置的加权系数。i(x,y)表示强度分量i在(x，y)位置处的像素值。
[0099]
(2)光谱估计：
[0100]
将i作为初始α通道，根据如下公式计算前景图像f和背景图像b，方便后续利用f、b进行图像重构。
[0101][0102]
其中i是第i颜色通道。α
ix
和α
iy
值是光谱前景fk、光谱背景bk和α的水平和垂直导数。
[0103]
f和b包含丰富的光谱信息，但不包含空间信息，后续步骤的主要目的是通过融合从全色图像中获得空间细节信息。
[0104]
(3)多尺度分解
[0105]
加权平均曲率滤波(weighted mean curvature filter，wmcf)可以在平滑图像的同时保留边缘信息，因此可以用于图像分解。此外，高斯滤波(gaussian filter，gf)是一种广泛使用的图像平滑算子(表示为gf(
·
))。在本发明中，多尺度分解是基于wmcf和gf的。利用wmcf和gf对pan图像进行多尺度分解，得到大尺度图像l、小尺度图像s和基本图像b。通过这种分解方法，可以将输入图像分解成三个不同尺度的子图像。分解过程为利用wmcf和gf对输入图像m进行处理，分别获得滤波图像iw和ig。其次，通过m和iw的差值得到小尺度图像(s)。大尺度图像(l)也是差分结果。基于wmcf和gf的图像分解方法流程图如图1所示。
[0106]
ig＝gf(m,σ,ρ).
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(4)
[0107]iw
＝wmcf(m,p).
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(5)
[0108]
s＝m-iw.
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(6)
[0109]
l＝i
w-ig.
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(7)
[0110]
b＝ig.
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(8)。
[0111]
(4)分量融合：
[0112]
在图像融合中，pcnn模型起着重要的作用，但它有一个缺点：对于一幅图像，一个像素对应一个神经元。一旦单独使用pcnn，就很难找到边缘、纹理等细节信息。在这种情况下，融合结果可能会发生局部模糊。连接系数β表示局部区域内像素值的变化情况，传统的融合方法大多选取固定的连接系数。但是，每个像素点之间的连接强度各不相同，所以固定的β不能改善融合图像的质量。局部能量是评价图像细节信息的重要指标，所以本发明利用
图像的局部能量e来动态调节连接系数β的大小，根据图像的特征自适应变化，以便充分保留图像边缘的特征信息。自适应连接系数β描述为：
[0113][0114]
其中，e(x,y)表示(x,y)处的局部能量。β
xy
表示(x,y)处的连接系数。局部能量越大，空间细节就越清晰。根据如上公式可看出，此时连接系数β的值相应也越大，从而会促进pcnn点火。局部能量的计算具体描述如下：
[0115][0116]
其中c(x,y)表示图像在(x，y)处的像素值。w为加权系数矩阵，具体细节如下：
[0117][0118]
将加权局部空间频率值作为脉冲耦合神经网络的输入。利用图像的局部能量e来动态调节连接系数β的大小，根据图像的特征自适应变化。当达到最大迭代次数时，迭代停止，得到强度分量i的点火次数总和ti和基本图像b的点火次数总和tb。根据点火次数总和取大原则进行融合，如下式所示。
[0119][0120]
其中fa(x,y)代表融合图像fa在(x,y)位置的像素值，i(x,y)代表强度分量i在(x,y)位置的像素值，b(x,y)代表基本图像b在(x,y)位置的像素值。
[0121]
最终，通过加权局部空间频率引导的脉冲耦合神经网络融合基本图像b和强度分量i得到融合图像fa。
[0122]
(5)图像重构：
[0123]
通过s、l和fa的线性组合重建融合后的图像fb作为最终的α通道。根据如下公式，组合α、f和b进行重构，得到最终的融合结果。
[0124]ii
＝α
ifi
(1-αi)biꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(13)
[0125]
其中fi是第i个像素的前景颜色，bi是第i个像素的背景颜色，ii第i个像素的颜色。
[0126]
经过上面步骤之后，便能得到最终的融合结果，图2为基于加权平均曲率滤波分解的遥感图像融合方法即本发明的流程图。
[0127]
补充说明：
[0128]
1.图像抠图模型：
[0129]
根据图像抠图模型，通过线性合成模型将输入图像区分为前景图像f和背景图像b，即第i个像素的颜色是对应的前景颜色和背景颜色的线性组合。
[0130]ii
＝α
ifi
(1-αi)biꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(14)
[0131]
其中fi是第i个像素的前景颜色，bi是第i个像素的背景颜色，其中，α是前景图像f的不透明度。根据图像抠图模型的局部线性假设，如果α包含该图像的大部分边缘信息，则该图像的对应的光谱前景和背景色在空间上将是平滑的。在确定输入图像和α的同时，通过
求解以下函数估计前景图像f和背景图像b。
[0132][0133]
其中i是第i颜色通道。α
ix
和α
iy
值是光谱前景fj、光谱背景bj和α的水平和垂直导数。
[0134]
2.脉冲耦合神经网络：
[0135]
pcnn是johnson等人研究哺乳动物视觉皮层工作机理，以eckhorn模型和rybak模型为基础改进和优化提出。单个神经元由输入部分、连接部分和脉冲发生器三个部分组成。pcnn模型的基本结构如图3所示。与人工神经网络不同，pcnn不需要任何训练过程。在pcnn模型中，图像像素和神经元存在一对一的对应关系。pcnn神经元的连接模型如图4所示。
[0136]
标准的pcnn模型由以下的公式描述：
[0137][0138][0139]uij
[n]＝f
ij
[n](1 βl
ij
[n]).
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(18)
[0140][0141][0142]
t
ij
[n]用于阈值调节，u
ij
[n]用于产生调制信号与阈值进行比较产生脉冲信号y
ij
[n]。共有9个参数，9个参数分属四类：时间衰减常数、放大系数、连接系数以及权值矩阵等。
[0143]
时间衰减常数：αf,α
l
,α
θ
，这三个参数分别来自输入域f、连接输入域l以及动态门限t。αf,α
l
控制着f，l通道信号的衰减速度，α
θ
控制着阈值衰减的速度。
[0144]
放大系数：vf、v
l
、ve。这三个参数分别是输入域f、连接输入域l以及动态门限t的放大系数。
[0145]
连接系数β：该参数来自调制域u主要起到调节神经元之间相互影响的作用，同时该参数还影响着中心神经元的点火周期。
[0146]
权值矩阵w，m：这两个参数分别来自输入域f、连接输入域l。矩阵中的元素表示神经元之间的连接强弱关系。
[0147]
工作原理：
[0148]
首先，接收域f接收外部刺激s、自身的刺激f
ij
[n-1]，l接收周围神经元的刺激。调制域针对接收域的信号f和l进行调制得到内部活动项u；将脉冲产生域神经元对应位置的内部活动项u与动态阈值t
ij
进行比较，若u
ij
＞t
ij
，则脉冲发生器启动神经元点火输出y
ij
＝1，否则神经元不点火，输出y
ij
＝0，得到最后结果y。上述过程不断重复，直到满足条件为止。
[0149]
3.加权局部空间频率：
[0150]
加权局部空间频率可以反映各像素邻域范围内的活跃程度，空间频率值越大，局部区域内各像素点的活跃度越高。但是，传统的空间频率只描述水平和垂直两个方向的空间信息而缺乏对角信息，导致融合结果可能丢失一些纹理、细节信息。然而，这些信息对融
合来说至关重要。传统的空间频率计算如下式。
[0151][0152]
其中，rf、cf分别是行频率、列频率。rf、cf包含水平、垂直两个方向。
[0153]
本发明对空间频率的计算进行改进，在原有基础上增加对角方向的计算。此外，基于欧几里得度量，确定行、列频率的加权系数为改进后的空间频率称为加权局部空间频率(weighted local spatial frequency，wlsf)。wlsf同时包含行、列和对角方向，对八个方向的空间频率进行加权求和。加权局部空间频率可以反映各像素邻域范围内的活跃程度，空间频率值越大，局部区域内各像素点的活跃度越高。加权局部空间频率的计算公式如下：
[0154][0155]
其中，rf、cf和df分别是行频率、列频率和对角频率。rf、cf包含水平、垂直两个方向，df包含正副对角线方向，具体定义如下：
[0156]
rf＝rf1 rf2,cf＝cf1 cf2ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(23)
[0157]
df＝df1 df2 df3 df4ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(24)
[0158][0159][0160][0161][0162][0163][0164][0165][0166]
其中，(2m 1)(2n 1)代表局部区域的大小。f(x,y)代表(x,y)位置处的像素值。
[0167]
4.加权平均曲率滤波：
[0168]
图片u的加权平均曲率(weighted mean curvature，wmc)定义为：
[0169][0170]
对于二维图像即n＝2，该等式变为：
[0171][0172]
其中δ是各向同性拉普拉斯算子，u
x
和uy是x和y方向的偏导数，u
xy
是相应的二阶偏导数。这种wmc定义既可以解释为按平均曲率加权的梯度，也可以解释为按梯度加权的平均曲率。wmc与梯度和mc的关系如图3所示。(33)和(34)是对wmc的连续形式的定义，但是实际图像是由离散的像素点组成的。因此，有必要对离散形式进行定义。
[0173]
对于3
×
3窗口，将考虑8个可能的法线方向和半窗口。所有8种情况都会产生相应的内核：
[0174][0175]
根据八个内核可以计算出八个距离：
[0176][0177]
其中*表示卷积操作。离散形式可以表示为：
[0178]hw
≈dk,k＝argmin{|di|；i＝1,2,...,8}.
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(37)
[0179]
wmc的简化过程可以表示为：
[0180]iout
＝wmc(i
in
).
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(38)
[0181]
其中i
in
和i
out
是输入和滤波后的图像，wmc(
·
)表示wmcf操作。
[0182]
以下通过实施例对本发明的技术方案作进一步说明。
[0183]
实施例
[0184]
本发明使用了一个包含27组图像的数据集进行实验，这些图像是由landsat 7etm 拍摄的。数据集的图像大小为400
×
400像素，在6个波段(红、绿、蓝、近红外、中红外和短波红外)工作。多光谱图像的空间分辨率为30m，全色图像的空间分辨率为15m。因此，由于数据集中没有高分辨率的ms图像作为参考图像，所以我们先将原始的多光谱图像进行上采样，
得到像素大小为400
×
400的多光谱图像。接着，将像素大小为400
×
400的多光谱图像和全色图像进行下采样，得到像素大小为200
×
200的多光谱图像和全色图像作为实验图像。随机选择三对不同场景的图像进行对照实验。最后，将原始ms图像作为参考图像，与各方法融合后的图像进行比较。图6展示了三组ms和pan图像，随后将用于实验分析。
[0185]
本发明将十种具有代表性的融合方法:(1)bt(brovey transform-based method)，(2)gsa(gram schmidt adaptive-based method)，(3)gf(guided filter-based method)，(4)ihs(intensity-hue-saturation-based method)，(5)ct(curvelet transform-based method)，(6)pca(principal component analysis-based method)，(7)pracs(partial replacement adaptive component substitution-based method)，(8)nsst-sr(nsst and sparse representation-based method)，(9)bdsd-pc(band dependent spatial-detail with physical constrains-basedmethod)，(10)wt(wavelettransform-based method)和本发明的融合方法进行比较。
[0186]
通常可以通过主观评价和客观评价来衡量一种遥感图像融合方法的性能。在进行主观评价时，通常会考虑目标的清晰度以及融合图像与原始多光谱图像光谱的接近程度。然而，仅凭主观评价难以对融合质量进行准确比较。为了对图像融合方法进行定量评价，采用若干指标来评价不同融合方法的性能。在实验中，使用六个众所周知的客观评价指标，并详细介绍如下。
[0187]
(1)相关系数(correlation coefficient，cc)计算参考图像和融合结果之间的相关性。该值越大，表示融合结果越接近参考图像，理想值是1。
[0188]
(2)光谱角制图(spectralangle mapper，sam)反映融合图像和参考图像之间的光谱失真。sam值越小，融合图像中的光谱失真越小。
[0189]
(3)光谱信息散度(spectral information divergence，sid)评测光谱之间的差异，理想值为0。
[0190]
(4)无参考质量评价指标(quality with no reference，qnr)可以在没有参考图像的情况下评估融合图像的质量，由光谱失真指数d
λ
、空间失真指数ds和全局qnr值三部分组成。对于全局qnr，值越大则融合效果越好，理想值为1。
[0191]
(5)d
λ
是qnr的亚度量，可以测量光谱失真。其值越小，融合效果越好，理想值为0。
[0192]
(6)ds是qnr的亚度量，可以测量空间畸变。其值越小，融合效果越好，理想值为0。
[0193]
图7、8、9中的(1)-(12)分别给出第一组图像、第二组图像、第三组图像对分别经bt、gsa、gf、ihs、ct、pca、pracs、nsst-sr、bdsd-pc、wt方法与本发明方法得到的融合结果，并与参考的多光谱图像作对比。为了更加直观地比较融合结果之间的差异，对融合结果进行局部放大，放大图像置于右下角。表1、2、3分别列出三组图像融合结果的六种包含光谱和空间质量评价的客观指标数值。对于所有的融合质量评价指标，最好的结果以红色粗体显示。
[0194]
表1第一组图像融合结果的客观评价
[0195][0196]
表2第二组图像融合结果的客观评价
[0197][0198][0199]
表3第三组图像融合结果的客观评价
[0200][0201]
从图7、图8、图9中可以看出，(1)bt、(2)gsa、(3)gf、(6)pca方法的融合图像在整体区域出现光谱失真，和参考多光谱图像相比光谱差异较大，在局部放大区域比较明显。(10)wt方法的融合图像空间细节失真较为严重。(3)gf方法的融合图像空间细节和光谱均发生一定程度的失真。其余方法保持了较好的光谱特性，但是局部放大部分空间细节比较模糊。从主观视觉角度分析，本发明方法的空间细节更加清晰，在改善空间细节的同时具有很好的光谱保持特性。
[0202]
从表1、2、3中可以看出，本发明方法在六个评价指标中的cc、sid、ds、qnr上表现最好，尤其是qnr值，在所有方法中遥遥领先，比较接近最优值1。从客观评价指标角度分析，本发明方法具有良好的空间细节保持特性和光谱保持特性，整体效果较好。
[0203]
综合本方法在第一、二、三组图像融合结果的客观评价结果以及视觉效果，本发明的方法在六个评价指标中的cc、sid、ds、qnr上表现较好，空间细节较为清晰且光谱分辨率较高。该方法具有良好的空间细节保持特性和光谱保持特性，能够在保留多光谱图像更多光谱信息的同时，从全色图像中获得更多的空间细节，充分验证了本发明方法的有效性。
[0204]
因此，本发明采用上述一种基于加权平均曲率滤波分解的遥感图像融合方法，解决了融合图像空间失真和光谱失真问题，通过与十种典型的融合方法进行比较，在主观视觉效果和客观评价方面都取得较好的结果。
[0205]
最后应说明的是：以上实施例仅用以说明本发明的技术方案而非对其进行限制，尽管参照较佳实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对本发明的技术方案进行修改或者等同替换，而这些修改或者等同替换亦不能使修改后的技术方案脱离本发明技术方案的精神和范围。