医学影像学在现代医学治疗和诊断等临床应用中有着至关重要的作用，而临床医生为了获取足够的信息，通常要用到多个模态的医学图像，如计算机断层扫描(CT)、正电子发射断层扫描(PET)以及磁共振成像(MR)等。然而不同的单一成像模态的医学影像都具有局限性，不能反映医生所需的所有细节信息，因此将来自不同模态的医学图像的信息整合到一幅图像中，将不同源图像中包含的互补信息相结合，能够帮助医生为医疗诊断做出更好的决策。

医学影像融合已经成为一种有效的信息融合技术。医学图像融合技术可以分为以下两大类：基于空间域的算法和基于变换域的算法^[1]。基于空间域的算法通常是通过图像块或者像素级的梯度信息去解决一个图像融合问题，但是容易导致图像细节丢失。与基于空间域的算法相比，基于变换域的算法在表征信号的局部特征方面有更多的优势，可以弥补基于空间域算法的细节表达不足的缺点，同时由于其更符合人类视觉系统(HVS)的生理机制，所以基于变换域的多尺度融合算法被认为是多模态图像融合中的主流融合方法^[2]。典型的基于变换域的方法包括基于多尺度分析的方法、基于稀疏表示的方法以及图像金字塔的方法。多尺度分析(MTA)的方法包括离散小波变换(DWT)^[3]、轮廓波变换(CT)^[4]、剪切变换(ST)^[5]等。这些方法都能很好地表示源图像的结构信息，但是它们能提取的方向信息是有限的。此外，由于这些方法缺少平移不变性，容易导致伪吉布斯效应^[6]。而非下采样轮廓波变换(NSCT)和非下采样剪切变换(NSST)虽然具有平移不变的性质，能够很好地克服伪吉布斯效应^[7-8]，但是这些方法有较高的计算复杂度，限制了它们的实际应用。

基于稀疏表示和字典学习的融合方法由于压缩感知理论^[9]的提出也被广泛应用于图像融合^[10-12]，并且它一般优于大多数传统融合方法^[13]。在基于稀疏表示的方法中，因为相邻的图像块是互相重叠的，导致每个像素的结果是多值的，在理想情况下，每个像素的多个值应该是相等的，以保持重叠图像块的一致性^[14]，但由于稀疏编码是在每个图像块上分别独立执行的，忽略了图像像素之间的相关性，导致每个像素有多个不等的值，同时大多数融合方法是通过采取聚合和平均的策略来获得每个像素的最终值^[15]，这将导致图像细节在融合中容易被平滑甚至丢失。Yin等人^[16]通过将源图像作为训练数据获得联合字典，然后使用最大加权多范数融合规则融合图像，但是其仍存在细节丢失的问题。Zong等人^[17]提出基于分类图像块的融合方法，使用方向梯度直方图(HOG)特征对图像块分类建立子字典，虽然减少了细节丢失，但仍不可避免导致一些细节被平滑。Zhang等人^[18]提出一种改进的多任务鲁棒稀疏表示结合空间上下文信息的融合方法，但该方法与大多数基于稀疏表示的方法一样，针对局部图像块编码而不是对整个图像，因此仍会导致细节保存不佳的情况。并且通常的基于稀疏编码的融合方法都只采用了一个字典来表示源图像不同的形态部分，容易导致源图像固有信息的次优表示, 不能充分地包含图像的特征信息，造成图像的信息丢失。Rong等人^[15]提出了对图像的不同成分采用不同的字典稀疏表示，但其字典是以解析的方式得到的，难以实现自适应描述。

针对上述问题，本文提出了基于卷积稀疏表示双重字典学习与自适应脉冲耦合神经网络(PCNN)的多模态医学图像融合方法(CSRDD-APCNN)。该方法将源图像分为稀疏与低秩两部分，构建了基于卷积稀疏表示的双重字典。其中，稀疏成分包含大量细节纹理，低秩成分包含轮廓亮度以及微小细节等基本信息。为了增强双重字典的区分性，分别对低秩成分与稀疏成分添加低秩与稀疏约束改进字典学习模型，并通过交替迭代求解，得到双重字典。然后，将总变分(TV)正则化(regularization)纳入图像分解模型，并通过基于傅里叶域的交替方向乘子法(ADMM)算法^[19]求解，将源图像分解为卷积稀疏与卷积低秩分量。同时为了充分保留源图像的细节信息，分别将卷积稀疏分量的改进的空间频率和(NMSF)与卷积低秩分量的改进的拉普拉斯能量和(NSML)作为自适应脉冲耦合神经网络的输入激励，通过输出的点火次数将各分量融合。最后通过组合卷积稀疏与卷积低秩分量来得到融合图像，进一步提高了融合图像的质量。

在卷积稀疏表示中，其基本思想是将源图像$\mathit{\boldsymbol{s}}$看做是一组特征响应$\mathit{\boldsymbol{x}}_{m}$与对应的字典过滤器$\mathit{\boldsymbol{d}}_{m}$的卷积之和。它的目的是实现对整个图像的稀疏表示，而不是针对重叠的图像块分别进行稀疏逼近。其相应的最为广泛使用的形式是卷积基追踪去噪(CBPDN)，定义为

$ \arg \mathop {\min }\limits_{\left| {{\mathit{\boldsymbol{x}}_m}} \right|} \frac{1}{2}\left\| {\sum\limits_m {{\mathit{\boldsymbol{d}}_m}} * {\mathit{\boldsymbol{x}}_m} - \mathit{\boldsymbol{s}}} \right\|_2^2 + \lambda \sum\limits_m {{{\left\| {{\mathit{\boldsymbol{x}}_m}} \right\|}_1}} $

(1)

式中，$\{\mathit{\boldsymbol{d}}_{m}\}$是一组卷积核也即卷积字典，$\mathit{\boldsymbol{x}}_{m}$是对应卷积核的稀疏系数，“*”表示2维的卷积操作，$λ$是正则化参数。为了对式(1)的求解问题进行优化，Bristow等人^[20]提出了基于傅里叶域交替方向乘子法(ADMM)框架。Wohlberg^[21]提出了一种联合使用空域和傅里叶域解决卷积稀疏表示(CSR)问题的算法。CBPDN的自然字典学习优化形式为

$ \begin{array}{*{20}{c}} {\mathop {\arg }\limits_{\left| {{\mathit{\boldsymbol{d}}_m}} \right|} \mathop {\min }\limits_{\left\{ {{\mathit{\boldsymbol{x}}_m},k} \right\}} \frac{1}{2}\sum\limits_k {\left\| {\sum\limits_m {{\mathit{\boldsymbol{d}}_m} * {\mathit{\boldsymbol{x}}_{m,k}} - {\mathit{\boldsymbol{s}}_k}} } \right\|_2^2} + }\\ {\lambda \sum\limits_{m,k} {{{\left\| {{\mathit{\boldsymbol{x}}_{m,k}}} \right\|}_1}} \quad {\rm{ s}}{\rm{.t}}{\rm{. }}{{\left\| {{\mathit{\boldsymbol{d}}_m}} \right\|}_2} = 1\forall m} \end{array} $

(2)

式中，$\mathit{\boldsymbol{s}}_{k}$表示第$k$幅训练图像；$\mathit{\boldsymbol{x}}_{m, k}$表示对应于第$m$个卷积核与第$k$幅图像的卷积系数。其中对卷积字典$\mathit{\boldsymbol{d}}_{m}$的规范约束是为了避免滤波器和系数产生缩放歧义。

PCNN是脉冲耦合神经元横向连接的2维网络，每个神经元对应一个特定的像素，其受到周围神经元的影响。单个神经元如图 1所示，包含3部分：接受域、调制域和脉冲发生器。PCNN模型描述为

$ \left\{ {\begin{array}{*{20}{l}} {{F_{ij}}(n) = {S_{ij}}}\\ {{L_{ij}}(n) = {{\rm{e}}^{ - {\alpha _L}}}{L_{ij}}(n - 1) + {V_L}\sum\limits_{k,l} {{W_{ijkl}}} {Y_{kl}}(n - 1)}\\ {{U_{ij}}(n) = {F_{ij}}(n)\left( {1 + \beta {L_{ij}}(n)} \right)}\\ {{\theta _{ij}}(n) = {{\rm{e}}^{ - {\alpha _\theta }}}{\theta _{ij}}(n - 1) + {V_\theta }{Y_{ij}}(n - 1)}\\ {{Y_{ij}}(n) = \left\{ {\begin{array}{*{20}{c}} \begin{array}{l} 1\\ 0 \end{array}&\begin{array}{l} {U_{ij}}(n) > {\theta _{ij}}(n)\\ 其他 \end{array} \end{array}} \right.} \end{array}} \right. $

(3)

式中，$F_{ij}$和$L_{ij}$分别为PCNN神经元的反馈输入和连接输入；$β$是链接强度，$θ_{ij}$是动态阈值，$W_{ijkl}$表示突触权重，来调整中心神经元附近每个神经元的影响程度；$S_{ij}$、$Y_{ij}$以及$U_{ij}$分别是神经元的外部刺激、输出和内部活动；$α_{L}$和$α_{θ}$是时间衰减常数；$V_{L}$和$V_{θ}$是神经元连接输入和动态阈值的振幅增益。

该PCNN模型主要反映在$F$通道处只考虑了外部图像像素的灰度值的影响。首先，在PCNN用于图像融合时，其中外部刺激$S_{ij}$输入对应像素$(i, j)$处的归一化灰度值。随后，在调制域中执行非线性调制，获得内部活动$U_{ij}$，然后与阈值$θ_{ij}$比较，若$U_{ij}$大于$θ_{ij}$，则神经元产生脉冲输出，也就是神经元点火，$Y=1$。反之，$Y=0$。此过程反复迭代，达到设定的迭代次数为止。经过点火，最后由所有神经元的点火次数构成的点火映射图作为输出。

在PCNN融合算法中，链接强度$β$在融合过程中起关键作用。在传统的基于PCNN的融合过程中，$β$通常根据经验手动设置，并且所有神经元设置相同的链接强度。而根据文献[22]，$β$的值在不同的神经元中不应该完全相同，其与相应像素的特征有关。所以根据图像特征自适应确定$β$，有助于提升融合效果。关于自适应$β$的确定方法，有空间频率、拉普拉斯能量、区域平均梯度和梯度能量等。其中，区域平均梯度能很好地反映图像的清晰程度，图像越清晰，区域平均梯度越大。

针对目前基于稀疏表示的融合方法存在的细节保存能力有限的问题，本文提出了一种基于卷积稀疏表示双重字典学习与自适应PCNN的多模态医学图像融合方法。其整体流程如图 2所示，首先通过建立双重字典学习模型，学习出具有判别性的卷积低秩子字典与卷积稀疏子字典，然后通过图像分解模型，将源图像$\mathit{\boldsymbol{A}}$、$\mathit{\boldsymbol{B}}$分别分解为卷积稀疏与卷积低秩分量。其次通过对$\mathit{\boldsymbol{A}}$、$\mathit{\boldsymbol{B}}$的卷积稀疏与卷积低秩分量采用不同的融合规则，来得到融合后的稀疏与低秩图像，最后通过组合不同分量来得到最终的融合图像。该方法通过对整幅图像进行建模，从而能够得到更好的细节保存效果，提高医学图像融合的质量。

定义$\mathit{\boldsymbol{y}}_{m}∈{\bf {R}}^{M}(m=1, 2, …, N)$为训练样本中第$m$幅图像(大小为$\sqrt M \times \sqrt M $)，在将源图像分解为卷积稀疏与低秩分量时，为了保证字典的重构能力，假设训练图像$\mathit{\boldsymbol{y}}_{m}$为两个不同分量的叠加

$ {\mathit{\boldsymbol{y}}_m} = \mathit{\boldsymbol{y}}_m^c + \mathit{\boldsymbol{y}}_m^t $

式中，$\mathit{\boldsymbol{y}}^{c}_{m}$和$\mathit{\boldsymbol{y}}^{t}_{m}$分别是第$m$幅训练图像的低秩分量和稀疏分量，低秩分量包含图像的本质低维空间，稀疏分量包含图像的边缘纹理信息。为了保证两个分量在各自所对应的字典下系数的稀疏性，采用了$L_{0}$范数的最优凸近似$L_{1}$范数作为正则项。定义$\{\mathit{\boldsymbol{x}}^{c}_{m, k}\}$和$\{\mathit{\boldsymbol{x}}^{t}_{m, k}\}$为第$m$幅训练图像的低秩系数与稀疏系数，并引入${\mathit{\boldsymbol{d}}_{c}}$和${\mathit{\boldsymbol{d}}_{t}}$表示训练出的卷积低秩与卷积稀疏子字典，其滤波器的数量分别为$K_{c}$与$K_{t}$，使得$\mathit{\boldsymbol{y}}_m^c = \sum\limits_{k = 1}^{{K_c}} {{\mathit{\boldsymbol{d}}_{c, k}}} *\mathit{\boldsymbol{x}}_{m, k}^c, \mathit{\boldsymbol{y}}_m^t = \sum\limits_{k = 1}^{{K_t}} {{\mathit{\boldsymbol{d}}_{t, k}}} *\mathit{\boldsymbol{x}}_{m, k}^t$。在图 3中给出了字典学习过程，建立字典学习模型为

$ \begin{array}{*{20}{c}} {\mathop {\min }\limits_{\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{d}}_c},{\mathit{\boldsymbol{d}}_t}}\\ {{\mathit{\boldsymbol{x}}_c},{\mathit{\boldsymbol{x}}_t}} \end{array}} \left\{ {\frac{1}{2}\sum\limits_{m = 1}^M {\left( {\left\| {{\mathit{\boldsymbol{y}}_m} - \sum\limits_{k = 1}^{{K_c}} {\mathit{\boldsymbol{d}}_k^c} * \mathit{\boldsymbol{x}}_{m,k}^c - \sum\limits_{k = 1}^{{K_t}} {\mathit{\boldsymbol{d}}_k^t} * \mathit{\boldsymbol{x}}_{m,k}^t} \right\|_2^2 + } \right.} } \right.}\\ {\left. {\left. {{\lambda _c}\sum\limits_{k = 1}^{{K_c}} {{{\left\| {\mathit{\boldsymbol{x}}_{m,k}^c} \right\|}_1}} + {\lambda _t}\sum\limits_{k = 1}^{{K_t}} {{{\left\| {\mathit{\boldsymbol{x}}_{m,k}^t} \right\|}_1}} } \right)} \right\}} \end{array} $

(4)

式中，$λ_{c}$和$λ_{t}$为平衡参数。由文献[19, 23]，通过离散傅里叶变换可将空间域的卷积转化为频域的乘积运算，因而定义频域变量$\mathit{\boldsymbol{Y}}=(\mathit{\boldsymbol{y}}_{1}, \mathit{\boldsymbol{y}}_{2}, …, \mathit{\boldsymbol{y}}_{m})$表示训练样本，$m$是训练样本的数目；字典$\mathit{\boldsymbol{D}}_{k}$，使得${\mathit{\boldsymbol{x}}_c} = \left({\begin{array}{*{20}{c}} {{\bf{x}}_{0, 0}^c}&{{\bf{x}}_{0, 1}^c}& \cdots \\ {{\bf{x}}_{1, 0}^c}&{{\bf{x}}_{1, 1}^c}& \cdots \\ \vdots & \vdots & \ddots \end{array}} \right) \in {{\bf{R}}^{K \times N}}$分别是卷积低秩与卷积稀疏子字典，其中$M$和$K$分别表示字典的大小与滤波器的数量；定义矩阵${\mathit{\boldsymbol{x}}_i} = \left({\begin{array}{*{20}{c}} {{\bf{x}}_{0, 0}^t}&{{\bf{x}}_{0, 1}^t}&\cdots \\ {{\bf{x}}_{1, 0}^t}&{{\bf{x}}_{1, 1}^t}& \cdots \\ \vdots & \vdots & \ddots \end{array}} \right) \in {{\bf{R}}^{K \times N}}$分别表示卷积低秩与卷积稀疏系数矩阵，其中$K$和$N$分别表示字典的大小与样本图像的数量。将式(4)转化为

$ \begin{array}{*{20}{c}} {\mathop {\min }\limits_{\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{D}}_c},{\mathit{\boldsymbol{D}}_t}}\\ {{\mathit{\boldsymbol{x}}_c},{\mathit{\boldsymbol{x}}_t}} \end{array}} \left\{ {\frac{1}{2}\left\| {\mathit{\boldsymbol{Y}} - {\mathit{\boldsymbol{D}}_c}{\mathit{\boldsymbol{x}}_c} - {\mathit{\boldsymbol{D}}_t}{\mathit{\boldsymbol{x}}_t}} \right\|_2^2 + } \right.}\\ {\left. {{\mathit{\boldsymbol{\lambda }}_c}{{\left\| {{\mathit{\boldsymbol{x}}_c}} \right\|}_1} + {\lambda _t}{{\left\| {{\mathit{\boldsymbol{x}}_t}} \right\|}_1}} \right\}} \end{array} $

(5)

为了增强字典的判别能力，通过对卷积稀疏与卷积低秩分量引入稀疏约束和低秩约束，对卷积低秩子字典添加低秩约束并引入辅助变量$\mathit{\boldsymbol{Y}}_{c}$，$\mathit{\boldsymbol{Y}}_{t}$，使得$\mathit{\boldsymbol{Y}}_{c}=\mathit{\boldsymbol{D}}_{c}\mathit{\boldsymbol{x}}_{c}$，$\mathit{\boldsymbol{Y}}_{t}=\mathit{\boldsymbol{D}}_{t}\mathit{\boldsymbol{x}}_{t}$，将式(5)改进为

$ \begin{array}{*{20}{c}} {\mathop {\min }\limits_{\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{D}}_c}}&{{\mathit{\boldsymbol{D}}_t}}\\ {{\mathit{\boldsymbol{x}}_c}}&{{\mathit{\boldsymbol{x}}_t}} \end{array}} \left\{ {\frac{1}{2}\left\| {\mathit{\boldsymbol{Y}} - {\mathit{\boldsymbol{Y}}_c} - {\mathit{\boldsymbol{Y}}_t}} \right\|_2^2 + \left\| {{\mathit{\boldsymbol{Y}}_c} - {\mathit{\boldsymbol{D}}_c}{\mathit{\boldsymbol{x}}_c}} \right\|_2^2 + } \right.}\\ {\left\| {{\mathit{\boldsymbol{Y}}_t} - {\mathit{\boldsymbol{D}}_t}{\mathit{\boldsymbol{x}}_t}} \right\|_2^2 + {\lambda _c}{{\left\| {{\mathit{\boldsymbol{x}}_c}} \right\|}_1} + {\lambda _t}{{\left\| {{\mathit{\boldsymbol{x}}_t}} \right\|}_1} + }\\ {\left. {{\lambda _1}{{\left\| {{\mathit{\boldsymbol{Y}}_c}} \right\|}_ * } + {\lambda _2}{{\left\| {{\mathit{\boldsymbol{Y}}_t}} \right\|}_1} + {\lambda _3}{{\left\| {{\mathit{\boldsymbol{D}}_c}} \right\|}_ * }} \right\}} \end{array} $

(6)

式中，||·||_*是核范数，表示矩阵奇异值的和。

式(6)的求解对于各变量为非凸优化问题，本文采用交替迭代的方式将式(6)分解为3个子问题求解^[24]。包括辅助变量更新、稀疏编码与字典更新3个阶段：

1) 随机初始化字典及其他变量，在更新辅助变量$\mathit{\boldsymbol{Y}}_{c}$时，保持其他变量不变，对$\mathit{\boldsymbol{Y}}_{c}$矩阵进行更新，即

$ \mathop {\min }\limits_{{\mathit{\boldsymbol{Y}}_c}} \left\{ {\left\| {\mathit{\boldsymbol{Y}} - {\mathit{\boldsymbol{Y}}_c} - {\mathit{\boldsymbol{Y}}_t}} \right\|_2^2 + \left\| {{\mathit{\boldsymbol{Y}}_c} - {\mathit{\boldsymbol{D}}_c}{\mathit{\boldsymbol{x}}_c}} \right\|_2^2 + {\lambda _1}{{\left\| {{\mathit{\boldsymbol{Y}}_c}} \right\|}_*}} \right\} $

(7)

式中，通过引入辅助变量$\hat {\bf{Y}}$，${\widetilde {\bf{Y}}_c}$，使得$\hat {\bf{Y}}=[\mathit{\boldsymbol{Y}}－\mathit{\boldsymbol{Y}}_{t}; \mathit{\boldsymbol{D}}_{c}\mathit{\boldsymbol{x}}_{c}]$，$\widetilde {\bf{Y}}_{c}=[\mathit{\boldsymbol{Y}}_{c}; \mathit{\boldsymbol{Y}}_{c}]$。将式(7)等价转化为

$ \mathop {\min }\limits_{{{\mathit{\boldsymbol{\tilde Y}}}_c}} \left\{ {\left\| {\mathit{\boldsymbol{\hat Y}} - {{\mathit{\boldsymbol{\tilde Y}}}_c}} \right\|_2^2 + {\lambda _1}{{\left\| {{{\mathit{\boldsymbol{\tilde Y}}}_c}} \right\|}_*}} \right\} $

(8)

该式可看做为核范数最小化的问题，通过奇异值阈值算法^[25]求解。更新变量$\mathit{\boldsymbol{Y}}_{t}$时，通过保持其他变量不变，以与更新$\mathit{\boldsymbol{Y}}_{c}$同样的方式可得优化问题

$ \mathop {\min }\limits_{{{\mathit{\boldsymbol{\tilde Y}}}_t}} \left\{ {\left\| {\mathit{\boldsymbol{\hat Y}} - {{\mathit{\boldsymbol{\tilde Y}}}_t}} \right\|_2^2 + {\lambda _2}{{\left\| {{{\mathit{\boldsymbol{\tilde Y}}}_t}} \right\|}_1}} \right\} $

(9)

式中，$\hat {\mathit{\boldsymbol{Y}}}=[\mathit{\boldsymbol{Y}}－\mathit{\boldsymbol{Y}}_{c}; \mathit{\boldsymbol{D}}_{t}\mathit{\boldsymbol{x}}_{t}]$，$\widetilde {\mathit{\boldsymbol{Y}}}_{c}=[\mathit{\boldsymbol{Y}}_{t}; \mathit{\boldsymbol{Y}}_{t}]$。该式通过迭代收缩算法^[26]求解，从而得到更新的辅助变量$\mathit{\boldsymbol{Y}}_{t}$。

2) 稀疏编码阶段。在更新卷积低秩系数$\mathit{\boldsymbol{x}}_{c}$时，通过保持其余变量不变，得

$ {\mathit{\boldsymbol{x}}_c} = \arg \mathop {\min }\limits_{{\mathit{\boldsymbol{x}}_c}} \left\{ {\left\| {{\mathit{\boldsymbol{Y}}_c} - {\mathit{\boldsymbol{D}}_c}{\mathit{\boldsymbol{x}}_c}} \right\|_2^2 + {\lambda _c}{{\left\| {{\mathit{\boldsymbol{x}}_c}} \right\|}_1}} \right\} $

(10)

同理在更新卷积稀疏系数$\mathit{\boldsymbol{x}}_{t}$时，保持其余变量不变，得

$ {\mathit{\boldsymbol{x}}_t} = \arg \mathop {\min }\limits_{{\mathit{\boldsymbol{x}}_t}} \left\{ {\left\| {{\mathit{\boldsymbol{Y}}_t} - {\mathit{\boldsymbol{D}}_t}{\mathit{\boldsymbol{x}}_t}} \right\|_2^2 + {\lambda _t}{{\left\| {{\mathit{\boldsymbol{x}}_t}} \right\|}_1}} \right\} $

(11)

显然式(10)(11)为凸函数，通过对字典与系数进行离散傅里叶变换，将其变换到频域，可将空间域的卷积转化为频域的乘积，再通过交替方向乘子法(ADMM)^[19]求解，从而得到更新的低秩系数与稀疏系数矩阵$\mathit{\boldsymbol{x}}_{c}$与$\mathit{\boldsymbol{x}}_{t}$。

3) 字典更新阶段。在更新卷积低秩字典$\mathit{\boldsymbol{D}}_{c}$时，保持其余变量不变，并引入辅助变量$\mathit{\boldsymbol{Z}}$，使得$\mathit{\boldsymbol{D}}_{c}=\mathit{\boldsymbol{Z}}$，得

$ \mathop {\min }\limits_{{\mathit{\boldsymbol{D}}_c}} \left\| {{\mathit{\boldsymbol{Y}}_c} - {\mathit{\boldsymbol{D}}_c}{\mathit{\boldsymbol{x}}_c}} \right\|_2^2 + {\lambda _3}{\left\| \mathit{\boldsymbol{Z}} \right\|_*}\quad {\rm{ s}}{\rm{. t}}{\rm{. }}{\mathit{\boldsymbol{D}}_c} = \mathit{\boldsymbol{Z}} $

(12)

在更新卷积稀疏字典$\mathit{\boldsymbol{D}}_{t}$时，保持其余变量不变，得

$ \mathop {\min }\limits_{{D_t}} \left\| {{\mathit{\boldsymbol{Y}}_t} - {\mathit{\boldsymbol{D}}_t}{\mathit{\boldsymbol{x}}_t}} \right\|_2^2 $

(13)

式(12)为约束优化问题，由增广拉格朗日乘子法^[2]求解，可得更新的子字典$\mathit{\boldsymbol{D}}_{c}$。式(13)也以相同的方式求解，得更新的卷积稀疏子字典$\mathit{\boldsymbol{D}}_{t}$。最后，循环执行上述3个步骤，直至达到设定的循环次数停止，输出卷积低秩与卷积稀疏子字典$\mathit{\boldsymbol{D}}_{c}$，$\mathit{\boldsymbol{D}}_{t}$。

根据2.1节学习到的卷积低秩和卷积稀疏子字典$\{\mathit{\boldsymbol{d}}_{c}\}$及$\{\mathit{\boldsymbol{d}}_{t}\}$，将源图像$\mathit{\boldsymbol{A}}∈{\bf {R}}^{M×N}$分解为卷积稀疏和卷积低秩两部分分量。建立图像分解模型

$ \begin{array}{*{20}{c}} {\mathop {\min }\limits_{{\mathit{\boldsymbol{x}}_{c,k}},{\mathit{\boldsymbol{x}}_{t,k}}} \frac{1}{2}\left\| {A - \sum\limits_{k = 1}^{{K_c}} {{\mathit{\boldsymbol{d}}_{c,k}}} * {\mathit{\boldsymbol{x}}_{c,k}} - \sum\limits_{k = 1}^{{K_t}} {{\mathit{\boldsymbol{d}}_{t,k}}} * {\mathit{\boldsymbol{x}}_{t,k}}} \right\|_2^2 + }\\ {{\lambda _1}\sum\limits_{k = 1}^{{K_c}} {{{\left\| {{\mathit{\boldsymbol{x}}_{c,k}}} \right\|}_1}} + {\lambda _2}\sum\limits_{k = 1}^{{K_t}} {{{\left\| {{\mathit{\boldsymbol{x}}_{t,k}}} \right\|}_1}} + }\\ {{\lambda _3}{f_{{\rm{TV}}}}\left( {\sum\limits_{k = 1}^{{K_c}} {{\mathit{\boldsymbol{d}}_{c,k}}} * {\mathit{\boldsymbol{x}}_{c,k}}} \right)} \end{array} $

(14)

式中, $λ_{1}$、$λ_{2}$和$λ_{3}$表示正则化参数，采用了总变分正则化(TV)的各向异性的规则，即${f_{{\rm{TV}}}}(x) = ||{\mathit{\boldsymbol{g}}_0}*\mathit{\boldsymbol{x}}|{|_1} + ||{\mathit{\boldsymbol{g}}_1}*\mathit{\boldsymbol{x}}|{|_1}$和$\mathit{\boldsymbol{g}}_{1}$分别为计算图像沿$x$轴和$y$轴方向梯度的过滤器。

式(14)通过更新卷积稀疏与卷积低秩系数交替迭代进行，分为以下两个步骤：

1) 在更新$\mathit{\boldsymbol{x}}_{c, k}$时，通过固定$\mathit{\boldsymbol{x}}_{t, k}$，将式(14)转化为

$ \begin{array}{*{20}{c}} {\mathop {\min }\limits_{{\mathit{\boldsymbol{x}}_{c,k}}} \frac{1}{2}\left\| {\mathit{\boldsymbol{A}} - \sum\limits_{k = 1}^{{K_t}} {{\mathit{\boldsymbol{d}}_{t,k}} * {\mathit{\boldsymbol{x}}_{t,k}}} - \sum\limits_{k = 1}^{{K_c}} {{\mathit{\boldsymbol{d}}_{c,k}}} } \right\|_2^2 + }\\ {{\lambda _1}\sum\limits_{k = 1}^{{K_c}} {{{\left\| {{\mathit{\boldsymbol{x}}_{c,k}}} \right\|}_1}} + {\lambda _3}{f_{{\rm{TV}}}}\left( {\sum\limits_{k = 1}^{{K_c}} {{\mathit{\boldsymbol{d}}_{c,k}} * {\mathit{\boldsymbol{x}}_{c,k}}} } \right)} \end{array} $

(15)

该式可视为带TV正则的稀疏编码问题，通过基于离散傅里叶变换的ADMM^[19]算法求解，由此得到更新的低秩系数矩阵。

2)在更新$\mathit{\boldsymbol{x}}_{t, k}$时，通过固定$\mathit{\boldsymbol{x}}_{c, k}$，转化原式(14)为

$ \begin{array}{*{20}{c}} {{{\mathit{\boldsymbol{\hat x}}}_{t,k}} = \mathop {\arg \min }\limits_{{\mathit{\boldsymbol{x}}_{t,k}}} \frac{1}{2}\left\| {\mathit{\boldsymbol{A}} - \sum\limits_{k = 1}^{{K_e}} {{\mathit{\boldsymbol{d}}_{c,k}} * {\mathit{\boldsymbol{x}}_{c,k}}} - } \right.}\\ {\left. {\sum\limits_{k = 1}^{{K_t}} {{\mathit{\boldsymbol{d}}_{t,k}} * {\mathit{\boldsymbol{x}}_{t,k}}} } \right\|_2^2 + {\lambda _2}\sum\limits_{k = 1}^{{K_t}} {{{\left\| {{\mathit{\boldsymbol{x}}_{t,k}}} \right\|}_1}} } \end{array} $

(16)

为了便于优化，引入辅助变量$\mathit{\boldsymbol{y}}_{r}$，使得${\mathit{\boldsymbol{y}}_r} = \mathit{\boldsymbol{A}} - \sum\limits_{k = 1}^{{K_c}} {{\mathit{\boldsymbol{d}}_{c, k}}} *{\mathit{\boldsymbol{x}}_{c, k}}$，将式(16)转化为

$ \begin{array}{*{20}{c}} {{{\mathit{\boldsymbol{\hat x}}}_{t,k}} = \arg \mathop {\min }\limits_{{\mathit{\boldsymbol{x}}_{t,k}}} \frac{1}{2}\left\| {{\mathit{\boldsymbol{y}}_r} - \sum\limits_{k = 1}^{{K_t}} {{\mathit{\boldsymbol{d}}_{t,k}} * {\mathit{\boldsymbol{x}}_{t,k}}} } \right\|_2^2 + }\\ {{\lambda _2}\sum\limits_{k = 1}^{{K_t}} {{{\left\| {{\mathit{\boldsymbol{x}}_{t,k}}} \right\|}_1}} } \end{array} $

(17)

在求解式(17)时，通过ADMM框架^[23]交替迭代求解，由此可得到更新的稀疏系数矩阵。最后循环执行上述两个步骤，直至达到预设的循环次数为止，输出源图像$\mathit{\boldsymbol{A}}$的低秩系数矩阵与稀疏系数矩阵，并通过${\mathit{\boldsymbol{\hat y}}_c} = \sum\limits_{k = 1}^{{K_e}} {{\mathit{\boldsymbol{d}}_{c, k}}} *{\mathit{\boldsymbol{\hat x}}_{c, k}}与{\mathit{\boldsymbol{\hat y}}_t} = \sum\limits_{k = 1}^{{K_t}} {{\mathit{\boldsymbol{d}}_{t, k}}} *{\mathit{\boldsymbol{\hat x}}_{t, k}}$卷积操作，来得到源图像$\mathit{\boldsymbol{A}}$的卷积低秩分量与卷积稀疏分量。

由2.2节可得源图像$\mathit{\boldsymbol{A}}$与$\mathit{\boldsymbol{B}}$的稀疏分量分别为$\mathit{\boldsymbol{S}}_{A}$与$\mathit{\boldsymbol{S}}_{B}$。稀疏分量包含了图像大量细节纹理信息，其将直接影响融合图像的视觉质量。由于改进的空间频率和(NMSF)能很好地反映图像的清晰度和活动水平，同时也可检测输入纹理的梯度变化，因此对于稀疏分量的融合，采用自适应PCNN算法进行纹理的提取，先求出稀疏分量改进的空间频率和，并作为PCNN的外部刺激，再以反映图像清晰度的区域平均梯度确定链接强度$β$。最后通过PCNN输出的点火次数图对稀疏分量进行融合。

1) 计算稀疏分量改进的空间频率$\mathit{\boldsymbol{NMSF}}_{A}(i, j)$和$\mathit{\boldsymbol{NMSF}}_{B}(i, j)$。设卷积稀疏分量$\mathit{\boldsymbol{S}}_{A}$和$\mathit{\boldsymbol{S}}_{B}$的像素系数为$\mathit{\boldsymbol{S}}_{A}(i, j)$和$\mathit{\boldsymbol{S}}_{B}(i, j)$，则首先通过式(18)计算各分量的对角频率$D^{A}_{f}$和$D^{B}_{f}$，然后通过式(19)求出卷积稀疏分量分别为$\mathit{\boldsymbol{NMSF}}_{A}(i, j)$和$\mathit{\boldsymbol{NMSF}}_{B}(i, j)$

$ \begin{array}{*{20}{c}} {{D_f} = \left[ {\sqrt {\frac{1}{{\left( {M - 1} \right)\left( {N - 1} \right)}}\left[ {\sum\limits_{i = 2}^M {\sum\limits_{j = 2}^N {{{\left( {\mathit{\boldsymbol{S}}\left( {i,j} \right) - \mathit{\boldsymbol{S}}\left( {i - 1,j - 1} \right)} \right)}^2}} } } \right]} + } \right.}\\ {{{\left. {\sqrt {\frac{1}{{\left( {M - 1} \right)\left( {N - 1} \right)}}\left[ {\sum\limits_{i = 2}^M {\sum\limits_{j = 2}^N {{{\left( {\mathit{\boldsymbol{S}}\left( {i - 1,j} \right) - \mathit{\boldsymbol{S}}\left( {i,j - 1} \right)} \right)}^2}} } } \right]} } \right]}^2}} \end{array} $

(18)

$ \mathit{\boldsymbol{NMSF}} = \sqrt {\frac{1}{{M(N - 1)}}\left[ {\sum\limits_{i = 1}^M {\sum\limits_{j = 2}^N {(\mathit{\boldsymbol{S}}(} } i,j - 1) - \mathit{\boldsymbol{S}}(i,j){)^2}} \right] + \frac{1}{{(M - 1)N}}\left[ {\sum\limits_{i = 2}^M {\sum\limits_{j = 1}^N {(\mathit{\boldsymbol{S}}(} } i,j) - \mathit{\boldsymbol{S}}(i - 1,j){)^2}} \right] + {D_f}} $

(19)

2) 计算区域平均梯度$\mathit{\boldsymbol{\bar g}}^{S}_{A}(i, j)$和$\mathit{\boldsymbol{\bar g}}^{S}_{B}(i, j)$。区域平均梯度能反映图像的清晰度，区域平均梯度越大，则该图像越清晰。对于稀疏分量系数$\mathit{\boldsymbol{S}}_{A}(i, j)$和$\mathit{\boldsymbol{S}}_{B}(i, j)$，区域平均梯度为

$ \begin{array}{l} \begin{array}{*{20}{l}} {\mathit{\boldsymbol{\bar g}}(i,j) = \frac{1}{9}\sum\limits_{a = - 1}^1 {\mathop \sum \limits_{b = - 1}^1 } }\\ {{{\left\{ {\frac{{\left[ {{\mathit{\boldsymbol{g}}_1}(i + a,i + b) + {\mathit{\boldsymbol{g}}_2}(i + a,i + b)} \right]}}{2}} \right\}}^{1/2}}} \end{array}\\ \begin{array}{*{20}{l}} {{\mathit{\boldsymbol{g}}_1}(i,j) = {{[\mathit{\boldsymbol{S}}(i,j) - \mathit{\boldsymbol{S}}(i, + 1,j)]}^2}}\\ {{\mathit{\boldsymbol{g}}_2}(i,j) = {{[\mathit{\boldsymbol{S}}(i,j) - \mathit{\boldsymbol{S}}(i,j + 1)]}^2}} \end{array} \end{array} $

(20)

3) 采用PCNN模型求得卷积稀疏分量的点火频率图$\mathit{\boldsymbol{T}}^{S}_{A}(i, j)$和$\mathit{\boldsymbol{T}}^{S}_{B}(i, j)$。将$\mathit{\boldsymbol{NMSF}}_{A}(i, j)$和$\mathit{\boldsymbol{NMSF}}_{B}(i, j)$作为PCNN的刺激输入，并将$\mathit{\boldsymbol{\bar g}}^{S}_{A}(i, j)$和$\mathit{\boldsymbol{\bar g}}^{S}_{B}(i, j)$作为链接强度$β$, 即

$ F_{i,j}^S[n] = \mathit{\boldsymbol{NMSF}}_{i,j}^S $

(21)

$ \beta (i,j) = {\mathit{\boldsymbol{\bar g}}^S}(i,j) $

(22)

4) 在PCNN达到迭代次数$n$时，根据输出点火频率图$\mathit{\boldsymbol{T}}^{S}_{A}(i, j)$和$\mathit{\boldsymbol{T}}^{S}_{B}(i, j)$，求得融合图像稀疏分量$\mathit{\boldsymbol{F}}_{S}$，即

$ {\mathit{\boldsymbol{F}}_S}(i,j) = \left\{ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{S}}_A}(i,j)}&{T_{i,j,A}^S\left[ {{n_{\max }}} \right] \ge T_{i,j,B}^S\left[ {{n_{\max }}} \right]}\\ {{\mathit{\boldsymbol{S}}_B}(i,j)}&{T_{i,j,A}^S\left[ {{n_{\max }}} \right] < T_{i,j,A}^S\left[ {{n_{\max }}} \right]} \end{array}} \right. $

(23)

经2.2节后，源图像$\mathit{\boldsymbol{A}}$、$\mathit{\boldsymbol{B}}$的卷积低秩分量$\mathit{\boldsymbol{C}}_{A}$和$\mathit{\boldsymbol{C}}_{B}$包含了图像的大量基本信息，如图像的亮度信息和轮廓信息，其融合至关重要。传统的平均策略不仅会造成图像的细节模糊而且会导致亮度下降；邻域能量策略虽然较好地保留了图像重要边缘，但是会模糊较低亮度的边缘信息。改进的拉普拉斯能量和(NSML)可以较好地反映图像的边缘特征信息，在一定程度上能恰当地反映图像的聚焦特性和清晰度，因此对于卷积低秩分量的融合，先求出$\mathit{\boldsymbol{C}}_{A}$和$\mathit{\boldsymbol{C}}_{B}$改进的拉普拉斯能量和，并作为PCNN的外部刺激，再以区域平均梯度确定链接强度$β$，尽可能保留微小细节。最后通过PCNN输出的点火次数图对低秩分量进行融合。

1) 计算卷积低秩分量改进的拉普拉斯能量和$\mathit{\boldsymbol{NSML}}_{A}(i, j)$和$\mathit{\boldsymbol{NSML}}_{B}(i, j)$。对于卷积低秩分量系数$\mathit{\boldsymbol{C}}_{A}(i, j)$和$\mathit{\boldsymbol{C}}_{B}(i, j)$，其由式(24)(25)求得$\mathit{\boldsymbol{NSML}}_{A}(i, j)$和$\mathit{\boldsymbol{NSML}}_{B}(i, j)$，即

$ \mathit{\boldsymbol{NSML}}(i,j) = \sum\limits_n {\sum\limits_k {\mathit{\boldsymbol{w}}(a,b)\mathit{\boldsymbol{F}}(i + a,j + b)} } $

(24)

$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{F}}(i,j) = \left| {2\mathit{\boldsymbol{C}}(i,j) - \mathit{\boldsymbol{C}}(i - 1,j) - \mathit{\boldsymbol{C}}(i + 1,j)} \right| + }\\ {\left| {2\mathit{\boldsymbol{C}}(i,j) - \mathit{\boldsymbol{C}}(i - 1,j) - \mathit{\boldsymbol{C}}(i + 1,j)} \right|} \end{array} $

(25)

式中，$\mathit{\boldsymbol{w}}(a, b)$表示滑动窗口区域对应的权重值，$a$、$b$指窗口的大小，通常可以取$3×3$，$5×5$或者$7×7$。为了更好地匹配变化不剧烈的区域对比度与体现卷积低秩分量的频率变化，本文选取窗口大小为$3×3$的矩阵，即

$ \mathit{\boldsymbol{w}}(a,b) = \left[ {\begin{array}{*{20}{l}} {1/15}&{2/15}&{1/15}\\ {2/15}&{3/15}&{2/15}\\ {1/15}&{2/15}&{1/15} \end{array}} \right] $

2) 计算$\mathit{\boldsymbol{C}}_{A}$和$\mathit{\boldsymbol{C}}_{B}$的区域平均梯度$\mathit{\boldsymbol{\bar g}}^{C}_{A}(i, j)$和$\mathit{\boldsymbol{\bar g}}^{C}_{B}(i, j)$。其计算过程如式(20)。

3) 通过PCNN模型计算$\mathit{\boldsymbol{C}}_{A}$和$\mathit{\boldsymbol{C}}_{B}$的点火频率图$\mathit{\boldsymbol{T}}^{C}_{A}(i, j)$和$\mathit{\boldsymbol{T}}^{C}_{B}(i, j)$。将$\mathit{\boldsymbol{NSML}}_{A}(i, j)$和$\mathit{\boldsymbol{NSML}}_{B}(i, j)$作为PCNN的刺激输入，并将$\mathit{\boldsymbol{\bar g}}^{C}_{A}(i, j)$和$\mathit{\boldsymbol{\bar g}}^{C}_{B}(i, j)$作为链接强度$β$，即

$ \mathit{\boldsymbol{F}}_{i,j}^C[n] = \mathit{\boldsymbol{NSML}}_{i,j}^C $

(26)

$ \beta (i,j) = {\mathit{\boldsymbol{\bar g}}^C}(i,j) $

(27)

4) 同样根据PCNN输出的点火频率图$\mathit{\boldsymbol{T}}^{C}_{A}(i, j)$和$\mathit{\boldsymbol{T}}^{C}_{B}(i, j)$，求得融合图像低秩分量$\mathit{\boldsymbol{F}}_{C}$

$ {\mathit{\boldsymbol{F}}_C}(i,j) = \left\{ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{C}}_A}(i,j)}&{\mathit{\boldsymbol{T}}_{i,j,A}^C\left[ {{n_{\max }}} \right] \ge \mathit{\boldsymbol{T}}_{i,j,B}^C\left[ {{n_{\max }}} \right]}\\ {{\mathit{\boldsymbol{C}}_B}(i,j)}&{\mathit{\boldsymbol{T}}_{i,j,A}^C\left[ {{n_{\max }}} \right] < \mathit{\boldsymbol{T}}_{i,j,A}^C\left[ {{n_{\max }}} \right]} \end{array}} \right. $

设已经配准的医学CT/MR源图像$\mathit{\boldsymbol{A}}$，$\mathit{\boldsymbol{B}}∈{\bf {R}}^{M×N}$，融合过程如图 4所示。

1) 字典学习阶段。通过2.1节中的字典学习模型，使用训练数据集$\mathit{\boldsymbol{Y}}=(\mathit{\boldsymbol{y}}_{1}, \mathit{\boldsymbol{y}}_{2}, …, \mathit{\boldsymbol{y}}_{m})$，进行交替迭代求解得卷积稀疏字典$\mathit{\boldsymbol{D}}_{t}$与卷积低秩字典$\mathit{\boldsymbol{D}}_{c}$。

2) 图像分解阶段。对源图像$\mathit{\boldsymbol{A}}$与$\mathit{\boldsymbol{B}}$分别通过交替方向乘子法(ADMM)求解其稀疏表示系数。使得源图像$\mathit{\boldsymbol{A}}$的稀疏表示系数为$\{\mathit{\boldsymbol{x}}^{A}_{c, k}, \mathit{\boldsymbol{x}}^{A}_{t, k}\}$，源图像$\mathit{\boldsymbol{B}}$的稀疏表示系数为$\{\mathit{\boldsymbol{x}}^{B}_{c, k}, \mathit{\boldsymbol{x}}^{B}_{t, k}\}$，然后通过卷积操作重构对应的分量。

3) 卷积稀疏分量融合阶段。求出卷积稀疏分量$\mathit{\boldsymbol{S}}_{A}$与$\mathit{\boldsymbol{S}}_{B}$改进的空间频率$\mathit{\boldsymbol{NMSF}}_{A}(i, j)$和$\mathit{\boldsymbol{NMSF}}_{B}(i, j)$以及区域平均梯度$\mathit{\boldsymbol{\bar g}}^{S}_{A}(i, j)$和$\mathit{\boldsymbol{\bar g}}^{S}_{B}(i, j)$，分别作为PCNN的刺激输入和链接强度$β$，计算相应的点火频率图$\mathit{\boldsymbol{T}}^{S}_{A}(i, j)$和$\mathit{\boldsymbol{T}}^{S}_{B}(i, j)$，最后通过融合规则求得融合卷积稀疏分量$\mathit{\boldsymbol{F}}_{S}$。

4) 卷积低秩分量融合阶段。求出卷积低秩分量$\mathit{\boldsymbol{C}}_{A}$与$\mathit{\boldsymbol{C}}_{B}$改进的拉普拉斯能量和$\mathit{\boldsymbol{NSML}}_{A}(i, j)$和$\mathit{\boldsymbol{NSML}}_{B}(i, j)$以及区域平均梯度$\mathit{\boldsymbol{\bar g}}^{C}_{A}(i, j)$和$\mathit{\boldsymbol{\bar g}}^{C}_{B}(i, j)$，分别作为PCNN的刺激输入和链接强度$β$，计算相应的点火频率图$\mathit{\boldsymbol{T}}^{C}_{A}(i, j)$和$\mathit{\boldsymbol{T}}^{C}_{B}(i, j)$，通过融合规则求得融合卷积低秩分量$\mathit{\boldsymbol{F}}_{C}$。

5) 通过进一步组合融合后的分量得到最终的融合图像

$ \mathit{\boldsymbol{F}} = {\mathit{\boldsymbol{F}}_c} + {\mathit{\boldsymbol{F}}_s} $

为了验证所提方法的性能，选取了3组经过配准的脑部多模态医学图像(大小均为256×256像素)作为待融合的源图像，其中源图像数据来源于哈佛大学医学院^[28]，包括脑血管疾病CT/MR图像、脑部肿瘤MR/PET图像、脑部肿瘤MR/SPECT图像。同时采取了6种对比算法进行比较，选取的对比算法包括：1)基于多尺度变换和稀疏表示的图像融合算法(MST-SR)^[29]；2)基于自适应稀疏表示的图像融合算法(ASR)^[30]；3)基于NSCT和单元链接PCNN的自适应医学图像融合算法(NSCT-PCNN)^[31]；4)基于联合补丁聚类的医学图像融合算法(JPC)^[32]；5)基于滑动窗口和稀疏表示的医学图像融合算法(SWT-SR)^[33]；6)基于非下采样轮廓波变换算法(NSCT)^[34]。将本文方法记做(CSRDD-APCNN)。

同时在训练卷积稀疏与低秩字典时，所需的训练图像如图 5所示。对于MSR-SR方法而言，其字典训练样本是从20幅自然的图像中随机采样10万个图像块，其滑动步长设置为1，滑动窗口大小设置为8 × 8像素，字典的大小设置为64 × 512像素，并设置误差$ε$=0.03和稀疏度$τ$=8，且MST-SR方法的多尺度变换选择LP，分解水平设为3。NSCT方法使用“9-7”金字塔滤波器和“c-d”方向滤波器，分解水平为$\{2^{2}, 2^{2}, 2^{3}, 2^{4}\}$。对于其余的ASR、NSCT-PCNN、JPC及SWT-SR方法，所有的参数分别根据文献[30-33]中的设置进行设定。

本文方法涉及的参数是由多次实验调参所得，在这里对于字典学习模型的式(6)而言，将加权指数$λ_{c}$和$λ_{t}$都设置为0.005，同时将低秩与稀疏分量的加权指数$λ_{1}$和$λ_{2}$分别设置为1和0.01，将$λ_{3}$设置为1。将如式(14)所示的图像分解模型的参数$λ_{1}$、$λ_{2}$与$λ_{3}$设置为0.38、3.37和0.05。在使用ADMM算法迭代求解时，将外循环迭代次数设置为200，内循环次数设置为100，在图像融合阶段所使用的PCNN的参数设置为$α_{L}=0.3$，$α_{T}=0.1$，$β=0.2$，$V_{L}=1$，$V_{T}=10$，并将其迭代次数设置为200。图 6给出了所训练出的卷积稀疏与卷积低秩字典。

本文采用了4种评价指标对实验融合图像的质量进行定量评价：标准差(SD)^[35]、标准化互信息(MI)^[36]、空间频率(SF)以及边缘信息传递因子($Q^{AB/F}$)^[37]。其中标准差SD是用来测量融合图像的对比度，标准差越大则对应的对比度越高；互信息MI可以度量源图像和融合图像在灰度分布上的相似程度，提取的信息越多，融合质量越好，在计算互信息MI时，首先分别计算融合图像$\mathit{\boldsymbol{F}}$与源图像$\mathit{\boldsymbol{A}}$和$\mathit{\boldsymbol{B}}$的之间的互信息为

$ {I_{FA}}(f,a) = \sum\limits_{f,a} {{p_{FA}}} (f,a){\log _2}\frac{{{p_{FA}}(f,a)}}{{{p_F}(f){p_A}(a)}} $

$ {I_{FB}}(f,b) = \sum\limits_{f,b} {{p_{FB}}} (f,b){\log _2}\frac{{{p_{FB}}(f,b)}}{{{p_F}(f){p_B}(b)}} $

式中，$p$为联合概率密度函数。根据上式计算的结果得出图像融合质量测量指标结果

$ M_F^{AB} = {I_{FA}}(f,a) + {I_{FB}}(f,b) $

空间频率SF用来测量图像的整体活动和清晰度水平，该指标定义为

$ {f_{SF}} = \sqrt {f_{RF}^2 + f_{CF}^2} $

式中，$f_{RF}$表示图像的行频率，$f_{CF}$表示图像的列频率

$ {f_{RF}} = \sqrt {\frac{1}{{M(N - 1)}}\sum\limits_{i = 0}^{M - 1} {\sum\limits_{j = 0}^{N - 2} {(\mathit{\boldsymbol{F}}(} } i,j + 1) - \mathit{\boldsymbol{F}}(i,j){)^2}} $

$ {f_{CF}} = \sqrt {\frac{1}{{(M - 1)N}}\sum\limits_{i = 0}^{M - 2} {\sum\limits_{j = 0}^{N - 1} {(\mathit{\boldsymbol{F}}(} } i + 1,j) - \mathit{\boldsymbol{F}}(i,j){)^2}} $

式中，$M$和$N$表示图像的行数和列数，其中空间频率$f_{SF}$的值越大，则融合效果越好; 边缘信息传递因子$Q^{AB/F}$通过使用Sobel边缘检测算子来测量从源图像传递到融合图像的边缘和方向信息，该指标的定义为

$ {Q^{AB/V}} = \frac{{\sum\limits_{i = 1}^M {\sum\limits_{j = 1}^N {\left[ {Q_{i,j}^{AF}\mathit{\boldsymbol{w}}_{i,j}^x + Q_{i,j}^{BF}\mathit{\boldsymbol{w}}_{i,j}^y} \right]} } }}{{\sum\limits_{i = 1}^M {\sum\limits_{j = 1}^N {\left[ {w_{i,j}^x + w_{i,j}^y} \right]} } }} $

式中，$Q^{AF}$与$Q^{BF}$定义为

$ Q_{i,j}^{AF} = Q_{g,i,j}^{AF}Q_{a,i,j}^{AF} $

$ Q_{i,j}^{BF} = Q_{g,i,j}^{BF}Q_{a,i,j}^{BF} $

式中，$Q^{*F}_{g}$和$Q^{*F}_{a}$表示源图像$\mathit{\boldsymbol{A}}$和$\mathit{\boldsymbol{B}}$分别在位置$(i, j)$处对应的边缘强度和方向保持值。其中$Q^{AB/F}$的动态范围为[0, 1]，越接近于1，其融合效果越好。

在图 7中给出了脑血管疾病CT/MR图像的融合结果，其中图像$\mathit{\boldsymbol{a}}$和$\mathit{\boldsymbol{b}}$分别表示CT源图像和MR源图像。从MST-SR、ASR、NSCT-PCNN、JPC、SWT-SR方法以及NSCT方法的融合结果看，其图像都存在皮质细节平滑的问题，并且MST方法在组织和病灶部位存在一定的伪影; ASR、SWT-SR方法以及

NSCT方法都存在对比度降低的问题。从表 1中可以看出，以上方法在MI指标以及$Q^{AB/F}$指标上都较低，表明其细节保留上存在缺陷，与主观视觉一致。其中MST-SR、ASR和SWT-SR方法细节保存不足的情况可能是在进行稀疏表示时，对每个像素对应的多值取平均时造成的细节平滑。本文方法在对比度以及细节保留方面都优于前面的方法，其细节结构是清晰的，有利于医生诊断和观察脑部皮质结构，从表 1量化评估上也可以看出，本文方法从源图像保留了较丰富的细节信息，融合质量较好。

在图 8中给出了脑部肿瘤疾病MR/PET图像的融合结果，其中图像$\mathit{\boldsymbol{a}}$和$\mathit{\boldsymbol{b}}$分别表示MR源图像和PET源图像。其中PET图像提供关于病变代谢的详细信息，而MR图像显示软组织信息和病灶的确切解剖位置，因此在融合过程中需要很好地保留颜色信息。从融合结果中可以看出, MST-SR、ASR、SWT-SR、NSCT、NSCT-PCNN以及JPC方法都存在颜色丢失和对比度下降的问题。同时，在ASR方法、SWT-SR方法的结果中都出现不同程度的伪影现象。而在NSCT-PCNN方法和JPC方法中出现皮质细节区域平滑的问题，同时从表 2的量化数据分析，MST-SR、ASR、SWT-SR以及NSCT方法在SD指标上偏低，也说明了其对比度下降的问题，以及从MI指标和$Q^{AB/F}$指标来看，NSCT、NSCT-PCNN方法偏低，表明其保留的边缘信息较少，这也与主观视觉观察较为一致。本文方法在颜色保留、对比度以及细节保留等方面表现的性能都较好，从表 2中也可观察各项指标都有较高的量值，融合质量较好，有利于医生对脑肿瘤疾病的诊断。

图 9中给出脑部肿瘤疾病的MR/SPECT图像的融合结果，图像$\mathit{\boldsymbol{a}}$和$\mathit{\boldsymbol{b}}$分别表示MR源图像和SPECT源图像。这里同样要较好地保留SPECT图像中的颜色信息，而从融合结果观察来看，在MST-SR、ASR、SWT-SR、NSCT-PCNN、JPC以及NSCT方法的融合图像中仍存在颜色传递失真的问题；在MST-SR、SWT-SR的融合图像中存在对比度下降以及细节模糊问题，这可能是字典信息不够丰富，不能提取充分的细节信息和医学图像特征导致的灰度不连续效应，以及MST-SR方法缺乏平移不变性和缺乏方向性等问题。同时从表 3的指标分析，MST-SR以及SWT-SR方法在SD指标上偏低，也说明其对比度下降的问题，与主观视觉观察较为一致。相比之下，本文方法从融合结果来看不会导致颜色和软组织信息失真，并且在融合结果中更好地保留了源图像中的补充信息，图像病灶部位细节清晰，对比度较高，没有伪影，融合质量较好，有利于医生进行诊断，在表 3中的数据也佐证了这一点，与主观视觉观察一致。

针对基于稀疏的融合方法存在保留细节能力有限的问题，提出了一种通过学习具有判别性的卷积低秩与卷积稀疏双重字典模型与自适应PCNN相结合的图像融合方法，通过字典学习与卷积系数编码将图像分解为低秩分量与稀疏分量，低秩分量融合规则采用改进的拉普拉斯能量和(NSML)激励脉冲耦合神经网络的方法，稀疏分量融合规则采用改进的空间频率和(NMSF)激励脉冲耦合神经网络的方法。从实验仿真可知，所提出的融合方法在客观评估和视觉质量方面明显优于对比的方法，保留了更丰富的特征和细节信息。

但本文仍存在字典学习模型及分解模型上需预设参数的问题，若参数设置不合适，会很大程度地影响学习的字典效率以及影响图像融合质量。探索如何有效地自适应确定参数的问题，是我们下一步研究的重点。