遥感图像融合是图像融合领域的一个重要分支，是指利用相关的技术手段将取自同一区域中的多源遥感图像进行整合的过程，而针对多光谱与全色图像的融合是当前遥感图像处理领域的研究热点之一。多光谱图像光谱信息丰富，具有较高的光谱分辨率，但空间分辨率较低，而全色图像可以反映更多的空间细节信息，具有较高的空间分辨率。二者融合的意义在于：能够在提高多光谱图像空间分辨率的同时保持光谱信息；结合全色图像中地物的形状、结构信息以及多光谱图像中的光谱信息，提供更强的解译能力和更加可靠的结果；提高地物的分类精度和目标的检测精度。当前随着多光谱遥感技术的快速发展，多光谱图像在灾情监测评估、地质资源勘探、军事目标检测与识别等领域发挥了重要作用。

传统的多光谱与全色图像的融合方法主要包括亮度色调饱和度(IHS)法^[1]，主成分分析(PCA)法^[2]以及正交变换(GS)法^[3]等分量替代方法。这些方法能够有效提高融合图像的空间分辨率，但存在较为严重的频谱失真问题。近年来，基于多分辨率分析的方法被广泛应用于多光谱图像融合领域，频谱失真问题得到了有效解决。多分辨率分析方法符合人眼视觉从粗到细的特性，且可以有效捕获图像的细节信息。其中小波变换是较为典型的多分辨率方法，能够得到更好的多光谱图像融合效果^[4]，但是由于小波变换的方向局限性使其不能有效地表示2维图像信号，为了寻找高维函数的最优表示方法，人们相继提出了许多新的多尺度几何分析方法。其中，非下采样剪切波变换(NSST)基于剪切波变换改进而来，具有良好的图像稀疏表示性能及低计算成本的特点，更适合于多光谱与全色图像的融合^[5]。NSST的出现为多光谱图像融合问题提供了新的思路。

脉冲耦合神经网络(PCNN)是第3代人工神经网络，具有空间邻近与特征相似聚集的特点，不通过训练和学习就可以从复杂背景中提取有效信息。应用于遥感图像融合，与多尺度几何分析相结合，可以更好地发挥其潜力，得到更好的融合效果。但PCNN模型应用于多光谱图像融合依然存在许多问题，还需进一步改进。

形态学算子在许多图像处理应用中已展示了其有效性，如图像编码^[6]、医学图像融合^[7]等，到目前为止，关于遥感图像融合只有Laporterie等人^[8]所做的一些初步的研究工作，以及Bejinariu等人^[9]、Restaino等人^[10]关于形态学信号分解应用于多光谱图像融合问题的研究工作。形态学方法具有很高的计算效率，能够有效提高融合图像的空间分辨率，可以进一步挖掘其在多光谱图像融合中的潜力。

为了有效融合全色图像的空间细节信息和多光谱图像的光谱信息，得到空间细节信息丰富、光谱保真度高的融合图像，本文提出了一种形态学滤波与改进脉冲耦合神经网络(PCNN)的NSST域多光谱与全色图像融合方法。首先分别对多光谱图像和全色图像进行NSST变换，在NSST域中，高频子带系数利用改进的PCNN进行融合，充分利用PCNN的优势，很好地保留细节信息以及光谱信息，低频子带系数采用形态学滤波和高通调制框架(HPM)进行融合，将全色图像低频子带的细节信息注入到多光谱图像低频子带中得到融合后的低频子带。最后通过NSST逆变换得到融合图像。将本文方法与亮度色调饱和度(IHS)法、主成分分析(PCA)法、基于非负矩阵分解(CNMF)、基于非下采样轮廓波变换和脉冲耦合神经网络(NSCT-PCNN)以及基于非下采样剪切波变换和脉冲耦合神经网络(NSST-PCNN)5种融合方法进行比较，验证本文方法的有效性和优越性。

Guo等人^[11]于2007年通过仿射系统将几何和多尺度结合构造了剪切波，它具有各向异性、结构简单、计算效率高、无方向数和支撑基尺寸限制等优点。当维数$n$=2时，具有合成膨胀的仿射系统表达式为

$ \begin{array}{*{20}{c}} {{\psi _{AB}}\left( \psi \right) = \left\{ {{\psi _{j,k,l}}\left( x \right) = } \right.}\\ {\left. {{{\left| {\det \mathit{\boldsymbol{A}}} \right|}^{j/2}}\psi \left( {{\mathit{\boldsymbol{B}}^l}{\mathit{\boldsymbol{A}}^j}x - k} \right),j,l \in {\bf{Z}},k \in {{\bf{Z}}^2}} \right\}} \end{array} $

(1)

式中，$\psi \in {L^2}\left( {{{\bf{R}}^2}} \right)$，$L$表示可积空间；$\mathit{\boldsymbol{A}}$和$\mathit{\boldsymbol{B}}$均为二阶可逆矩阵，且$\left| {\det \mathit{\boldsymbol{B}}} \right| = 1$。

若${\psi _{AB}}\left( \psi \right)$具有Parseval框架(紧框架)，即对任意$f \in {L^2}\left( {{{\bf{R}}^2}} \right)$满足

$ \sum\limits_{j,k,l} {{{\left| {\left\langle {f,{\psi _{j,k,l}}} \right\rangle } \right|}^2}} = {\left\| f \right\|^2} $

(2)

则${\psi _{AB}}\left( \psi \right)$的元素称为合成小波。其中$\mathit{\boldsymbol{A}}$称为具有各向异性的膨胀矩阵，${\mathit{\boldsymbol{A}}^j}$与尺度变换相关联，$\mathit{\boldsymbol{B}}$称为剪切矩阵，${\mathit{\boldsymbol{B}}^l}$与保持面积不变的几何变换相关联。$\mathit{\boldsymbol{A}}$和$\mathit{\boldsymbol{B}}$的表达式分别为

$ \mathit{\boldsymbol{A}} = \left[ {\begin{array}{*{20}{c}} a&0\\ 0&{{a^{1/2}}} \end{array}} \right],\mathit{\boldsymbol{B}} = \left[ {\begin{array}{*{20}{c}} 1&s\\ 0&1 \end{array}} \right] $

此时的合成小波即为剪切波，通常取$a = 4$，$s=1$。

如图 1所示，经过频率分解，剪切波的几何性质更为直观。对于不同的尺度，${\hat \psi _{j, k, l}}$支撑在近似大小为${2^{2j}} \times {2^j}$、线方向斜率为$l\;{2^{ - j}}$的梯形对上，支撑基的线方向斜率随着尺度$j$的变化而变化，使得这种近似长条形的支撑基具有各向异性，可以有效捕获方向信息，随着$j$的变换带来的旋转过程中，梯形面积并不会随着线方向斜率的改变而变化。关于剪切波变换的详细理论参见文献[12]。

为了使剪切波变换具备平移不变特性，Easley等人^[13]提出了非下采样剪切波变换(NSST)。NSST作为小波变换在多维方向的自然扩展，具有类似于小波变换的严格的数学基础，继承了曲波变换及轮廓波变换的优点，其基函数的楔形支撑空间具有可变性，可通过剪切和膨胀自适应对图像中的几何边缘进行表示，对于2维图像具有最优的逼近特性，是图像中边缘等特征的真正的稀疏表示。

NSST的离散化过程主要分为两个步骤：多尺度分解与方向局部化。多尺度分解通过非下采样金字塔滤波器组(NSP)来实现，每一级需要对上级中采用的滤波器组按照矩阵$\mathit{\boldsymbol{D}} = \left[ {\begin{array}{*{20}{c}} 2&0\\ 0&2 \end{array}} \right]$进行上采样，每一级分解均会产生一个低频子带和一个高频子带，后续分解均在低频分量上迭代进行，经过$k$级分解，可得到1个低频子带和$k$个高频子带，这$k + 1$个子带图像具有与源图像相同的尺寸大小。

方向局部化则通过剪切滤波器(SF)来实现。对某尺度的子带图像进行$j$级分解可以得到$\sum\limits_{j = 1}^J {{2^l}j + 1} $个与源图像尺寸相同的子带图像。其中$J$为图像分解层数，${l_j}$为尺度$j$下的方向分解级数。

在图像处理中，形态学是借助集合论的语言来描述的，是以形态结构元素为基础进行图像分析的数学工具。主要思想是利用具有一定形态的结构元素从图像中提取对表达和描绘区域形状有意义的图像分量, 从而达到图像分析与识别等目的。目前，形态学已广泛应用于图像去噪、图像增强、边缘提取、图像分割等领域。

令$\mathit{\boldsymbol{F}}$表示灰度图像，$\mathit{\boldsymbol{B}}$为结构元素，其本身也是一个子图像，利用$\mathit{\boldsymbol{B}}$对$\mathit{\boldsymbol{F}}$进行膨胀，记作$\mathit{\boldsymbol{F}} \oplus \mathit{\boldsymbol{B}}$，定义为

$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{F}} \oplus \mathit{\boldsymbol{B}}\left( {x,y} \right) = \mathop {\max }\limits_{\left( {i,j} \right)} \left\{ {\left. {\mathit{\boldsymbol{F}}\left( {x - i,y - j} \right) + \mathit{\boldsymbol{B}}\left( {i,j} \right)} \right|} \right.}\\ {\left. {\left( {x - i,y - j} \right) \in {\mathit{\boldsymbol{D}}_F},\left( {i,j} \right) \in {\mathit{\boldsymbol{D}}_B}} \right\}} \end{array} $

(3)

式中，${\mathit{\boldsymbol{D}}_F}$和${\mathit{\boldsymbol{D}}_B}$分别为灰度图像$\mathit{\boldsymbol{F}}$和结构元素$\mathit{\boldsymbol{B}}$的定义域。计算过程相当于让结构元素$\mathit{\boldsymbol{B}}$关于原点的镜像${\mathit{\boldsymbol{\hat B}}}$在图像$\mathit{\boldsymbol{F}}$的所有位置上滑过，过程中需保持$\left( {x + i, y + j} \right)$始终位于灰度图像$\mathit{\boldsymbol{F}}$之内。膨胀结果$\mathit{\boldsymbol{F}} \oplus \mathit{\boldsymbol{B}}$为在其定义域内每一点$\left( {x, y} \right)$处的取值为以该点为中心，在镜像${\mathit{\boldsymbol{\hat B}}}$规定的局部邻域内$\mathit{\boldsymbol{F}}$和${\mathit{\boldsymbol{\hat B}}}$之和的最大值。该过程与卷积相似，只是用最大值运算对卷积求和进行替代，用加法运算替代卷积运算。

利用$\mathit{\boldsymbol{B}}$对$\mathit{\boldsymbol{F}}$进行腐蚀，记$\mathit{\boldsymbol{F}} \odot \mathit{\boldsymbol{B}}$，定义为

$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{F}} \odot \mathit{\boldsymbol{B}}\left( {x,y} \right) = \mathop {\min }\limits_{\left( {i,j} \right)} \left\{ {\left. {\mathit{\boldsymbol{F}}\left( {x + i,y + j} \right) - \mathit{\boldsymbol{B}}\left( {i,j} \right)} \right|} \right.}\\ {\left. {\left( {x + i,y + j} \right) \in {\mathit{\boldsymbol{D}}_F},\left( {i,j} \right) \in {\mathit{\boldsymbol{D}}_B}} \right\}} \end{array} $

(4)

计算过程相当于让结构元素$\mathit{\boldsymbol{B}}$在灰度图像$\mathit{\boldsymbol{F}}$的所有位置上滑过。过程中需保证$\left( {x + i, y + j} \right)$始终在灰度图像$\mathit{\boldsymbol{F}}$之内。腐蚀结果$\mathit{\boldsymbol{F}} \odot \mathit{\boldsymbol{B}}$在其定义域内每一点$\left( {x, y} \right)$处的取值为以该点为中心，在$\mathit{\boldsymbol{B}}$规定的局部邻域内$\mathit{\boldsymbol{F}}$和$\mathit{\boldsymbol{B}}$之差的最小值。

在结构元素$\mathit{\boldsymbol{B}}$的值均大于零的情况下，膨胀的输出结果由于局部最大值运算的作用通常比输入图像更亮，同时消除了一些能够包含于结构元素的暗细节，也一定程度上减少了其他暗细节。而腐蚀的输出图像则暗于输入图像，小于结构元素的亮部细节其亮度会被削弱。

可以在灰度腐蚀和膨胀的基础上定义灰度图像的开运算和闭运算，灰度开运算为先进行灰度腐蚀后再灰度膨胀，灰度闭运算为先进行灰度膨胀后再灰度腐蚀，下面分别给出定义：

使用结构元素$\mathit{\boldsymbol{B}}$对图像$\mathit{\boldsymbol{F}}$进行灰度开运算，记作$\mathit{\boldsymbol{F}} \circ \mathit{\boldsymbol{B}}$，可表示为

$ \mathit{\boldsymbol{F}} \circ \mathit{\boldsymbol{B}} = \left( {\mathit{\boldsymbol{F}} \odot \mathit{\boldsymbol{B}}} \right) \oplus \mathit{\boldsymbol{B}} $

(5)

使用结构元素$\mathit{\boldsymbol{B}}$对图像$\mathit{\boldsymbol{F}}$进行灰度闭运算，记作$\mathit{\boldsymbol{F}} \bullet \mathit{\boldsymbol{B}}$，可表示为

$ \mathit{\boldsymbol{F}} \bullet \mathit{\boldsymbol{B}} = \left( {\mathit{\boldsymbol{F}} \oplus \mathit{\boldsymbol{B}}} \right) \odot \mathit{\boldsymbol{B}} $

(6)

开运算常用于去除相对于结构元素$\mathit{\boldsymbol{B}}$而言较小的高灰度区域，而对于较大的亮区域影响不大，虽然首先进行灰度腐蚀会在去除图像细节的同时降低图像的整体灰度，但随后的膨胀操作又会使图像的整体亮度得以增强，因而整体的灰度可以基本保持不变；通常利用闭运算去除图像中的暗细节部分，相对地保留图像中的高灰度部分不受影响。

在图像处理中，PCNN受生物视觉皮层模型启发而产生，它是一种反馈型网络，由若干神经元相互链接而构成，每个神经元包括接收域、调制域和脉冲发生器3个部分。由于原始PCNN模型结构比较复杂，用于图像处理时所需调整的参数较多，因而采用如图 2所示的简化的PCNN模型^[14]，其表达式为

$ \left\{ \begin{array}{l} {F_{ij}}\left( n \right) = {D_{ij}}\\ {L_{ij}}\left( n \right) = {L_{ij}}\left( {n - 1} \right) \times \exp \left( { - {\alpha _L}} \right) + \\ \;\;\;\;\;\;\;\;\;\;\;\;{V_L}\sum\limits_{pq} {{\omega _{ij,pq}}{Y_{pq}}\left( {n - 1} \right)} \\ {U_{ij}}\left( n \right) = {F_{ij}}\left( n \right)\left( {1 + \beta {L_{ij}}\left( n \right)} \right)\\ {\theta _{ij}}\left( n \right) = {\theta _{ij}}\left( {n - 1} \right) \times \exp \left( { - {\alpha _\theta }} \right) + \\ \;\;\;\;\;\;\;\;\;\;\;{V_\theta }{Y_{ij}}\left( {n - 1} \right)\\ {Y_{ij}}\left( n \right) = {\mathop{\rm sgn}} \left( {{U_{ij}}\left( n \right) - {\theta _{ij}}\left( n \right)} \right) \end{array} \right. $

(7)

式中，$\left( {i, j} \right)$表示图像中的像素位置；$n$为迭代次数；${D_{ij}}$、${F_{ij}}\left( n \right)$、${L_{ij}}\left( n \right)$、${U_{ij}}\left( n \right)$、${\theta _{ij}}\left( n \right)$分别为神经元外部输入的刺激信号、反馈输入项、链接输入项、内部活动项以及动态阈值；${\alpha _L}$、${\alpha _\theta }$为时间衰减常数；${V_L}$为链接输入的放大倍数；$\beta $为链接强度；${\omega _{ij, pq}}$为神经元的突触连接权系数；${V_\theta }$表示阈值放大系数；${Y_{ij}}\left( n \right)$为二值输出，${U_{ij}}\left( n \right) \ge {\theta _{ij}}\left( n \right)$则输出1，表示神经元被激活，发生点火，${U_{ij}}\left( n \right) < {\theta _{ij}}\left( n \right)$则输出0，表示神经元未被激活，未发生点火，继续等待捕获。

对于多光谱图像与全色图像而言，影响二者融合质量的主要因素一方面是空间细节信息的融入度，另一方面是光谱信息的保持度。经过NSST变换，源图像被分解为不同的频率分量，高频分量部分存储了图像中丰富的纹理、边缘等细节信息，是图像突变特性的反映，高频系数的选择对于空间细节信息的保留具有十分重要的作用，而低频分量部分是图像的近似子图，决定着图像的轮廓，低频系数的选择对于提高融合图像的视觉效果，很好地保持原始图像的光谱信息具有非常重要的意义。针对不同频率域需要保留的显著特征，对相应的高、低频分量的融合规则进行设计。

低频子带部分作为原始图像的近似子图，继承了原始图像的整体特性且集中了绝大部分的能量，但仍然有一些边缘和细节信息留存，本文采用一种基于形态学滤波和高通调制的细节信息注入方法对低频分量进行融合。

采用半梯度算子^[15]对全色图像的低频子带图像进行细节提取，其表达式为

$ \begin{array}{*{20}{c}} {{{\mathit{\boldsymbol{\bar \psi }}}_{{\rm{HG}},B}} = 0.5\left( {{\mathit{\boldsymbol{\rho }}^ - } - {\mathit{\boldsymbol{\rho }}^ + }} \right) = }\\ {0.5\left( {\mathit{\boldsymbol{F}} - {\varepsilon _B}\left( \mathit{\boldsymbol{F}} \right)} \right) - 0.5\left( {{\delta _B}\left( \mathit{\boldsymbol{F}} \right) - \mathit{\boldsymbol{F}}} \right)} \end{array} $

(8)

式中，${\mathit{\boldsymbol{\rho }}^ - } = \mathit{\boldsymbol{F}} - {\varepsilon _B}\left( \mathit{\boldsymbol{F}} \right)$为内梯度，${\mathit{\boldsymbol{\rho }}^ + } = {\delta _B}\left( \mathit{\boldsymbol{F}} \right) - \mathit{\boldsymbol{F}}$为外梯度，腐蚀算子${\varepsilon _B}\left( \mathit{\boldsymbol{F}} \right)$和膨胀算子${\delta _B}\left( \mathit{\boldsymbol{F}} \right)$为两个基本形态学算子，${\mathit{\boldsymbol{\rho }}^ - }$和${\mathit{\boldsymbol{\rho }}^ + }$均为应用腐蚀和膨胀算子得到的残差图像。则相应的形态学滤波可表示为

$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{\psi }}_{{\rm{HG}},B}} = \mathit{\boldsymbol{F}} - {{\mathit{\boldsymbol{\bar \psi }}}_{{\rm{HG}},B}} = }\\ {\mathit{\boldsymbol{F}} - \left[ {0.5\left( {\mathit{\boldsymbol{F}} - {\varepsilon _B}\left( \mathit{\boldsymbol{F}} \right)} \right) - 0.5\left( {{\delta _B}\left( \mathit{\boldsymbol{F}} \right) - \mathit{\boldsymbol{F}}} \right)} \right] = }\\ {0.5\left( {{\varepsilon _B}\left( \mathit{\boldsymbol{F}} \right) + {\delta _B}\left( \mathit{\boldsymbol{F}} \right)} \right)} \end{array} $

(9)

即等同于膨胀与腐蚀结果相加的0.5倍和，根据膨胀与腐蚀的计算公式，${\mathit{\boldsymbol{\psi }}_{{\rm{HG, }}B}}$得到的是在邻域范围内和函数最大值与差函数最小值的0.5倍和。

通过形态学滤波得到全色图像的低分辨率低频子带图像$\mathit{\boldsymbol{P}}_k^{{\rm{low}}}$后，基于高通调制框架(HPM)^[10]进行低频分量的融合，表达为

$ {\widehat {\mathit{\boldsymbol{MS}}}_k} = {\widetilde {\mathit{\boldsymbol{MS}}}_k} + {\widetilde {\mathit{\boldsymbol{MS}}}_k}\frac{{\mathit{\boldsymbol{P}}_k^{\rm{0}} - \mathit{\boldsymbol{P}}_k^{{\rm{low}}}}}{{\mathit{\boldsymbol{P}}_k^{{\rm{low}}}}} $

(10)

式中，$\widetilde {\mathit{\boldsymbol{MS}}} = {\left\{ {{{\widetilde {\mathit{\boldsymbol{MS}}}}_k}} \right\}_{k = 1, \cdots , N}}$为多光谱图像低频分量，$\mathit{\boldsymbol{P}}_k^0$为根据多光谱$k$通道低频子带图像$\mathit{\boldsymbol{M}}{\mathit{\boldsymbol{S}}_k}$将全色图像$PAN$的低频子带图像进行直方图均衡化得到的结果，${\widehat {\mathit{\boldsymbol{MS}}}_k}$为融合后的多光谱$k$通道低频子带图像。

这一表达式与全色图像的局部对比度密切相关^[16]，根据Weber关于对比度的定义^[17]，${C_{\rm{W}}} = \frac{{L - {L_{\rm{b}}}}}{{{L_{\rm{b}}}}}$，$L$为像素亮度分布，${{L_{\rm{b}}}}$为背景亮度分布，上式可以表示为

$ {\widehat {\mathit{\boldsymbol{MS}}}_k} = {\widetilde {\mathit{\boldsymbol{MS}}}_k}\left( {1 + {C_{\rm{W}}}} \right) $

(11)

基于形态学算子的低频分量融合步骤为：

1) 将全色图像在NSST域的低频子带图像依据多光谱图像在NSST域的低频子带图像进行直方图均衡化处理，得到图像$\mathit{\boldsymbol{P}}_k^0$；

2) 设置结构元素$\mathit{\boldsymbol{B = }}\left[ {\begin{array}{*{20}{c}} 0&1&0\\ 1&1&1\\ 0&1&0 \end{array}} \right]$；

3) 根据式(9)计算半梯度算子形态学滤波结果，得到全色图像的低分辨率低频子带图像$\mathit{\boldsymbol{P}}_k^{{\rm{low}}}$；

4) 根据低频分量融合式(10)计算低频子带系数的融合结果。

本文采用简化的PCNN模型对NSST分解得到的高频子带系数进行融合。采用PCNN模型可有效提取图像中的边缘以及纹理等细节信息，点火次数的多少则对应于图像相应位置的信息丰富程度。

现有的基于PCNN模型的融合算法，大都将PCNN模型输出中的点火次数作为判决准则进行融合系数的选取，通常采用硬限幅函数，点火输出非0即1，不能很好地反映同步脉冲激发的幅度差异。因此在简化的PCNN模型基础上，对输出部分进行改进，采用一个软限幅Sigmoid函数^[18]对各子带系数迭代过程中的点火输出幅度进行计算，将输出幅度和作为系数选择的判决依据。输出幅度为

$ {T_{ij}}\left( n \right) = \frac{1}{{1 + {{\rm{e}}^{{\theta _{ij}}\left( n \right) - {U_{ij}}\left( n \right)}}}} $

(12)

式中，${T_{ij}}\left( n \right)$为图像点$\left( {i, j} \right)$的点火输出幅度，${T_{ij}}\left( n \right) \in \left[ {0, 1} \right]$。$n$次点火输出幅度的总和为

$ {Z_{ij}}\left( n \right) = {Z_{ij}}\left( {n - 1} \right) + {T_{ij}}\left( n \right) $

(13)

基于改进PCNN的高频分量融合步骤为：

1) 采用3×3像素的区域滑动窗口分别对全色及多光谱图像进行遍历，计算NSST域中的高频分量对应的${F_{ij}}$，将其进行归一化处理作为PCNN模型输入；

2) 设置最大迭代次数${N_{\max }}$及PCNN模型初始值为${L_{ij}}\left( 0 \right) = {U_{ij}}\left( 0 \right) = {\theta _{ij}}\left( 0 \right) = 0$；

3) 根据式(7)进行模型的迭代计算，得到所有的中间结果，并计算出模型点火输出幅度的总和；

4) 当迭代次数$n < {N_{\max }}$时，重复执行步骤3)，直到迭代次数满足$n = {N_{\max }}$时停止迭代；

5) 根据模型点火输出幅度计算得到决策矩阵

$ H\left( {i,j} \right) = \left\{ {\begin{array}{*{20}{c}} \begin{array}{l} 1\\ 0 \end{array}&\begin{array}{l} Z_{ij}^{\rm{M}}\left( {{N_{\max }}} \right) > Z_{ij}^{\rm{P}}\left( {{N_{\max }}} \right)\\ 其他 \end{array} \end{array}} \right. $

基于决策矩阵完成融合结果图像中高频系数的选取，即

$ {S_{\rm{F}}}\left( {i,j} \right) = \left\{ {\begin{array}{*{20}{c}} \begin{array}{l} {S_{\rm{M}}}\left( {i,j} \right)\\ {S_{\rm{P}}}\left( {i,j} \right) \end{array}&\begin{array}{l} H\left( {i,j} \right) = 1\\ H\left( {i,j} \right) = 0 \end{array} \end{array}} \right. $

式中，${S_{\rm{M}}}\left( {i, j} \right)$，${S_{\rm{P}}}\left( {i, j} \right)$分别为NSST域多光谱图像与全色图像的高频系数。

针对全色图像与多光谱图像融合的特点，根据形态学滤波、PCNN以及NSST各自优势进行融合策略的设计，提出了一种结合形态学滤波和改进PCNN的NSST域多光谱与全色图像融合方法。方法流程如图 3所示，具体的步骤如下：

1) 对全色图像与多光谱图像进行精确配准，本文算法中的待融合图像均已完成精确配准；

2) 利用NSST变换对全色图像和多光谱图像进行多尺度多方向分解，得到全色图像的低频与高频子带系数$\left\{ {S_{\rm{P}}^{0, d}, S_{\rm{P}}^{l, d}} \right\}$、多光谱图像的低频与高频子带系数$\left\{ {S_{\rm{M}}^{0, d}, S_{\rm{M}}^{l, d}} \right\}$，设定NSST分解层数为3，方向数为{4，8，8}；

3) 针对全色图像与多光谱图像经NSST分解得到的低频部分，结合形态学滤波方法和高通调制框架将全色图像低频子带的细节信息注入到多光谱图像低频子带中，得到融合后的低频子带系数$\left\{ {S_{\rm{F}}^{0, d}} \right\}$；

4) 利用改进的PCNN方法对全色图像与多光谱图像经NSST分解得到的高频部分进行融合，得到融合后的高频子带系数$\left\{ {S_{\rm{F}}^{l, d}} \right\}$；

5) 对融合后的低频系数$\left\{ {S_{\rm{F}}^{0, d}} \right\}$和高频系数$\left\{ {S_{\rm{F}}^{l, d}} \right\}$进行NSST逆变换，得到最终的融合图像。

为了对本文提出的基于形态学滤波与改进PCNN的NSST域图像融合方法进行有效性验证，本文采用基于IHS^[1]、PCA^[2]、CNMF^[19]、NSCT-PCNN^[20]、NSST-PCNN^[21]5种融合方法与本文方法进行了比较，对大量多光谱图像与全色图像进行了融合实验。实验的硬件环境为Intel Core i5-4200U，CPU主频1.60 GHz，内存4 GB，软件环境为MATLAB R2014a。

选取3组图像数据进行实验分析，所用数据集均已完成配准。如图 4所示，3组图像来源于QuickBird卫星，于2002年11月21日拍摄于印度孙德尔本斯国家森林公园，所得全色图像与多光谱图像的大小分别为1 024×1 024像素和256×256像素，分辨率分别为2.44 m、0.61 m。由于缺少用于融合性能评估的参考图像，通过对原始多光谱及全色图像进行降采样因子为4的低通滤波降采样处理得到用于融合的图像数据，实验中，图像大小变为256×256像素和64×64像素，将原始多光谱图像作为参考图像。实验参数选取如下：NSST中，多尺度分解滤波器选取“maxflat”，尺度分解层数为3，方向数分别为{4，8，8}；改进PCNN中，$n = 200$、${\alpha _L} = 1.0$、${\alpha _\theta } = 0.2$、${V_L} = 1.0$、${V_\theta } = 20$、$\beta = 3.0$。

关于多光谱与全色图像融合方法的性能评价，分为主观和客观两类。其中视觉分析是最为直接的检测手段，同时为了对图像的融合效果进行客观的定量评价，本文选取清晰度、信息熵、相关系数、空间频率和光谱扭曲度^[22] 5种客观评价指标。

1) 清晰度，即平均梯度，反映了融合图像中的纹理和细节特征的变化程度，定义为

$ \begin{array}{*{20}{c}} {AG = \frac{1}{{M \times N}}\sum\limits_{i = 1}^M {\sum\limits_{j = 1}^N {} } }\\ {\sqrt {\begin{array}{*{20}{c}} {\frac{1}{2}\left( {{{\left( {F\left( {i,j} \right) - F\left( {i + 1,j} \right)} \right)}^2} + } \right.}\\ {\left. {{{\left( {F\left( {i,j} \right) - F\left( {i,j + 1} \right)} \right)}^2}} \right)} \end{array}} } \end{array} $

(14)

2) 信息熵，图像的信息熵用来衡量图像包含信息的丰富程度，定义为

$ EN = - \sum\limits_{i = 0}^{L - 1} {{p_F}\left( i \right){{\log }_2}\left( {{p_F}\left( i \right)} \right)} $

(15)

式中，$L$表示图像的灰度级数，${p_F}\left( i \right)$为图像$\mathit{\boldsymbol{F}}$相应灰度级分布概率。

3) 相关系数，反映了融合图像与多光谱图像之间的相关程度。定义为

$ CC = \frac{{\begin{array}{*{20}{c}} {\sum\limits_{m = 1}^M {\sum\limits_{n = 1}^N {\left[ {\left( {F\left( {m,n} \right) - \bar F} \right) \times } \right.} } }\\ {\left. {\left( {MS\left( {m,n} \right) - \overline {MS} } \right)} \right]} \end{array}}}{{\sqrt {\begin{array}{l} \sum\limits_{m = 1}^M {\sum\limits_{n = 1}^N {{{\left( {F\left( {m,n} \right) - \bar F} \right)}^2} \times } } \\ \sum\limits_{m = 1}^M {\sum\limits_{n = 1}^N {{{\left( {MS\left( {m,n} \right) - \overline {MS} } \right)}^2}} } \end{array} }}} $

(16)

式中，${\bar F}$和$\overline {MS} $分别为融合图像及多光谱图像相应波段的像素平均值。

4) 空间频率，用来度量图像空间域细节信息的丰富程度，定义为

$ {F_{\rm{S}}} = \sqrt {{F_{\rm{R}}}^2 + {F_{\rm{C}}}^2} $

(17)

式中，${F_{\rm{R}}}$、${F_{\rm{C}}}$分别表示行频及列频，即

$ {F_{\rm{R}}} = \sqrt {\frac{1}{{M \times N}}\sum\limits_{i = 1}^M {\sum\limits_{j = 2}^N {{{\left[ {F\left( {i,j} \right) - F\left( {i,j - 1} \right)} \right]}^2}} } } $

$ {F_{\rm{C}}} = \sqrt {\frac{1}{{M \times N}}\sum\limits_{i = 2}^M {\sum\limits_{j = 1}^N {{{\left[ {F\left( {i,j} \right) - F\left( {i - 1,j} \right)} \right]}^2}} } } $

5) 光谱扭曲度，反映了融合图像与多光谱图像之间的光谱扭曲程度。定义为

$ SD = \frac{1}{{M \times N}}\sum\limits_{i = 1}^M {\sum\limits_{j = 1}^N {\left| {F\left( {i,j} \right) - MS\left( {i,j} \right)} \right|} } $

(18)

图 5为3组多光谱与全色图像利用上述5种融合方法及本文方法获得的融合结果。从主观视觉效果来看，基于IHS方法的融合结果具有严重的光谱失真和空间信息缺失，如图 5(a)所示。图 5(b)(c)分别为基于PCA变换和CNMF方法的融合结果，这两种方法对于融合图像的细节信息有所加强，但依然存在较为严重的光谱扭曲。图 5(d)(e)分别为基于NSCT-PCNN方法和NSST-PCNN方法得到的融合图像，这两种方法较好地保持了多光谱图像的光谱信息，图像纹理、边缘等细节信息也较为清晰。图 5(f)为本文方法得到的融合结果，视觉效果上明显优于其他5种方法，细节信息清晰且光谱保真度高。

表 1、表 2和表 3分别给出了6种图像融合方法对于3组图像数据的定量评价结果。IHS变换方法光谱扭曲严重，空间分辨率也很差。PCA和CNMF方法空间分辨率和清晰度较优，具有优于IHS方法捕获图像细节的能力，但光谱扭曲度依然较高。NSCT-PCNN方法以空间频率${F_{\rm{S}}}$作为PCNN激励在一定程度上提高了融合结果的空间分辨率，但光谱扭曲度较高。NSST-PCNN方法相比于NSCT-PCNN方法空间频率略有降低，但光谱扭曲度相对较低，具有较高的相关系数值，反映了该方法在提升融合图像的空间细节融入度的同时，对于光谱信息的保持能力也有所提高。本文方法在融合效果上优势较为明显，各项指标与其他方法相比整体上较优，其中相关系数和光谱扭曲度两项指标的评价结果表明本文方法相比于其他5种方法能够更好地保持光谱信息，清晰度和空间频率两项指标的评价结果则展示了本文方法具有优于其他对比方法的空间细节注入能力，信息熵指标虽不是最优值，但与最优值非常接近，多波段清晰度、相关系数、空间频率及光谱扭曲度4项指标的平均值也为各方法中最优值。综合来看，本文方法在主观与客观方面的总体效果要优于其他5种方法。

针对多光谱与全色图像融合方法中存在的光谱信息损失以及纹理细节不够丰富的问题，提出了一种形态学滤波和改进脉冲耦合神经网络(PCNN)的NSST域多光谱与全色图像融合方法。该方法在对原始图像进行NSST变换的基础上，研究了基于形态学滤波和高通调制框架(HPM)的低频子带融合规则，以及基于改进脉冲耦合神经网络(PCNN)的高频子带融合规则，融合结果达到了增强图像局部空间细节的目的，同时很好地保留了多光谱图像的光谱信息。仿真实验结果表明本文方法对于融合全色图像空间细节信息和多光谱图像光谱信息的有效性和优越性。

下一步的工作需要对多尺度分解过程中保持边缘的滤波进行研究，克服该方法由于空间不连续性所带来的人造纹理和灰度不均衡现象，进一步改善多光谱图像融合质量。