磁共振成像(magnetic resonance imaging，MRI)以非接触、非电离、无需造影剂的优点，成为了医学辅助诊疗领域的一项关键技术。传统MRI设备利用单线圈获得全采样的测量数据，采样时间长，不利于实际应用，因此探索高效采样方式十分必要。近年来，结合并行采样和欠采样这2种快速采样策略可构建并行压缩感知磁共振成像(parallel compressive sensing magnetic resonance imaging，p-CSMRI)^[1]。p-CSMRI对应的算法需解决的问题是如何通过采集到的多线圈部分(欠采样)k域数据恢复原始信号。该恢复问题是一个病态问题，通常需要利用图像先验来限制解空间。如何利用先验来恢复高精度磁共振图像成为一大挑战。

传统基于模型的p-CSMRI算法首先通过正则化模型表征原始图像的先验信息，例如梯度域^[2]和小波变换域^[3]的稀疏性，然后通过合适的求解器求解相应的优化模型。如Block等^[4]提出了从多线圈径向欠采样数据中重建原始磁共振(magnetic resonance，MR)图像的算法，由惩罚线圈边缘轮廓信息和整幅图像的梯度稀疏性来引入先验知识。在梯度域和小波变换域利用全局稀疏变换法引入图像的全局稀疏先验，通常是固定的且非最优的。此外，上述先验忽略了图像的局部先验信息。为解决这些问题，Weller等^[5]提出了将自适应字典学习与L1-SPIRiT^[6]相结合的快速MRI算法，该算法利用了图像在自适应字典下的稀疏性，获得了比全局稀疏变换法更好的重建结果。另一种常用的先验是非局部相似性^[7-8]，如Zhang等^[9]提出了STDLR-SPIRiT算法，利用k域空间双向低秩特性和多线圈数据之间的相关性，通过无奇异值分解的数值算法解决优化问题。在这些方法中，基于全局稀疏变换的方法重建速度快但重建质量较差，特别是在低采样率下，重建图像存在严重伪影。相比利用全局先验的方法，基于自适应字典下稀疏性或非局部相似性的方法重建的图像质量更高，但重建速度慢。

随着人工智能的兴起，深度学习方法在图像处理任务中的应用引起广泛关注^[10]。深度学习方法通过深度神经网络(deep neural network，DNN)从大量数据中学习有效的特征或映射以刻画图像的潜在分布^[11]。这种方法可以充分利用数据集的分布信息，避免过度依赖人工经验设计的先验。将DNN用于p-CSMRI后研究出的诸多算法解决了传统方法重建速度慢、效果差的问题，如Wang等^[12]设计并训练了卷积神经网络以学习零填充(zero filled，ZF)重建和真实图像之间的映射关系，该数据驱动的端到端DNN能够克服基于传统先验的p-CSMRI算法的缺点，实现较好的重建质量，但其模型架构是不可解释的。架构不可解释不利于算法的收敛性分析，而收敛性会影响网络的稳定性。因此，研究者提出将迭代算法展开成模型驱动的DNN^[13-16]，如Wang等^[17]提出的PARCEL算法是一个展开的并行网络架构，从增强的多线圈部分采样k域数据中对比学习网络2个分支中的参数，并且2个网络训练时的损失函数相同；Yang等^[18]基于半二次分裂算法提出MD-DAN网络，该网络通过学习相应的近端算子来学习交叉对比度先验。研究表明，利用事先训练好的模型驱动DNN能够完成快速成像且取得较好的重建效果。虽然上述模型驱动的DNN在p-CSMRI任务上取得了一定的效果，但仍存在以下不足：1) 现有模型驱动的并行MRI网络利用深度神经网络学习了隐式的正则化项，但其先验网络架构缺乏模型可解释性；2) 由于网络架构的复杂性，现有模型驱动的p-CSMRI网络没有显式地给出对应迭代算法的收敛性分析。

针对上述问题，受Chan等^[19]提出的利用有界去噪器分析算法收敛性的启发，本文基于对偶标架提出了一种模型架构可解释且可证明有界的深度去噪器，并将其作为先验模块融入到模型驱动的p-CSMRI网络中，构建了收敛性可分析的模型驱动p-CSMRI算法。首先，本文基于对偶标架并结合深度阈值网络，构建满足有界条件的深度去噪器；其次，构建了基于对偶标架的p-CSMRI优化模型，并将迭代求解步骤展开成了一个可端到端有监督训练的深度神经网络；最后，利用构建的有界深度去噪器作为先验模块，提出了收敛性可分析的模型驱动p-CSMRI成像算法。通过利用有界去噪器作为先验网络构建可证明收敛的模型驱动p-CSMRI算法，网络的稳定性得到保证，同时实现了低采样率下的高质量图像重建。

1 基于可训练对偶标架的去噪器

有界去噪器的定义如下所示。

定义1  对于任意输入x∈$\mathbb{R}^N$，具有控制去噪器强度参数σ的去噪器D_σ(·)：$\mathbb{R}^N$→$\mathbb{R}^N$，如果满足

$ \left\|D_{\sigma}(\boldsymbol{x})-\boldsymbol{x}\right\|_{2}^{2} / N \leqslant \sigma^{2} C . $

则该去噪器为有界去噪器。其中C是独立于N和σ的通用常数。

有界Gauss去噪器^[19]是指任意一幅图像与该图像通过该Gauss去噪器后的差值能量有界且该边界与噪声方差呈线性关系。在噪声强度逐渐降低的前提下^[20]，有界Gauss去噪器对即插即用成像算法的收敛性保证至关重要^{[19, 21]}。然而Gauss去噪器的网络结构通常使用深度神经网络并根据经验设计，引入了很多深度学习的技巧，比如残差连接、批归一化等，属于深度去噪器。因此很难显式地证明设计的先验网络满足有界条件。为解决这些问题，本文基于对偶标架表示模型提出了一个可证明有界的深度去噪器，并首次尝试分析基于学习的模型驱动p-CSMRI算法的收敛性。

在成像算法中，标架和字典通常被用于稀疏表示的基。其中，紧标架能够表示整幅图像且稀疏编码速度较快^[22]。紧标架学习所对应的优化问题可描述为

$ \left\{\begin{array}{l} \min\limits_{\boldsymbol{W}, \boldsymbol{a}}\|\boldsymbol{a}-\boldsymbol{W} \boldsymbol{x}\|_{2}^{2}+\lambda\|\boldsymbol{a}\|_{1}\\ \text { s. t. } \boldsymbol{W}^{\mathrm{T}} \boldsymbol{W}=\boldsymbol{I} . \end{array}\right. $

(1)

其中：x∈$\mathbb{R}^N$表示原始图像，W∈$\mathbb{R}^{M \times N}$表示紧标架，a表示稀疏系数，I表示N×N的单位矩阵，λ表示正则化因子。对式(1)可以采用交替优化紧标架和稀疏系数的方式进行求解^[23]。但是上述传统交替优化的求解方式非常耗时，并且很难结合深度学习的优势进行图像表示，因此传统紧标架学习速度慢且紧标架表示能力有待提升。针对这些问题，本文将单紧标架扩展到对偶标架^[24]，其去噪或滤波对应的函数可表示为

$ D_{\varepsilon}(\bullet ; \sigma)=\mathit{\boldsymbol{ \boldsymbol{\varPhi} }} T_{\varepsilon(\sigma)}[\mathit{\boldsymbol{ \boldsymbol{\varPsi}}}(\bullet)] . $

其中：Ψ∈$\mathbb{R}^{M \times N}$表示解析标架，Φ∈$\mathbb{R}^{N \times M}$表示综合标架，二者满足对偶性质ΦΨ=I；T_ε(σ)(·)表示软阈值算子，soft_ε(σ)(·)=sign(·)max(|·|－ε(σ), 0)。深度去噪器滤波的过程可以看成求解包含隐式正则项的优化问题，可以写成以下邻近算子的形式：

$ D_{\varepsilon}(\bullet ; \sigma)=\hat{d}=\underset{d}{\operatorname{argmin}}\left\{\frac{1}{2}\|d-\bullet\|_{2}^{2}+\lambda R_{\text {dual }}(d)\right\}. $

其中R_dual(d)表示可学习的隐式正则化项。

图像像素在空域具有相关性，并且在优化过程中估计图像包含的噪声通常是空间变化的，因此阈值也应该是空间变化且与噪声相关的。本文设置阈值与噪声标准差σ呈线性相关^[25]。

$ \varepsilon_{l}=\sigma c_{l} . $

其中：ε_l是待处理特征图的第l个通道的阈值，c_l是相应的比例常数。不同于Isogawa等^[25]提出的利用反向传播算法来学习比例常数的思想，本文构建深度阈值网络从输入图像中确定这些常数，该方式能够有效提取空间变化的c_l。阈值向量e可以表示为

$ \boldsymbol{e}=\left[\varepsilon_{1}, \varepsilon_{2}, \cdots, \varepsilon_{q}\right]^{\mathrm{T}}=\boldsymbol{c} \otimes \boldsymbol{m}. $

其中：⊗表示逐元素乘积；c表示由c_l组成的比例常数向量；q表示通道的个数；m表示一个由噪声标准差σ扩展为与c相同大小的矩阵，称为噪声水平强度图。为了避免任意性，c_l的大小被限制为[c_min, c_max]。深度阈值网络的具体架构如图 1所示，该网络主要由卷积层(convolutional layer，Conv)、线性整流函数(rectified linear unit，ReLU)和拼接(Concat)操作构成。为了进一步有效提取常数，本文引入了稠密残差连接和通道注意力机制(对应图 1中的SE^[26]模块)。相关描述在2.2节中结合算法作详细展开。

本文将有界去噪器满足的条件(即有界且该界为噪声强度的线性函数)定义为有界条件。满足有界条件的深度去噪器有利于收敛性保证。引理1表明构建的深度去噪器满足有界条件。

引理1  对于任意输入x∈$\mathbb{R}^N$和独立于M(此处M表示阈值向量e的维度)的常数L=τ₀c_max²(τ₀是一个很小的正值)，所提出的基于对偶标架的可训练去噪器满足

$ \left\|\boldsymbol{x}-\mathit{\boldsymbol{ \boldsymbol{\varPhi} }} T_{\varepsilon(\sigma)}(\mathit{\boldsymbol{ \boldsymbol{\varPsi}}} \boldsymbol{x})\right\|_{2}^{2} \leqslant M \sigma^{2} L. $

本文假设通过所提出的先验网络学习到的解析标架和综合标架是满足严格对偶性质的。在实验中，上述2个标架可能不满足严格对偶条件，而是满足松弛的对偶条件。基于此种情况，本文在2.3节的推论1进行了算法的收敛性分析，并给出了相应证明。在此，假设学习到的解析标架和综合标架满足严格的对偶约束，即有ΦΨ=I。因此，

$ \begin{gather*} \left\|\boldsymbol{x}-\mathit{\boldsymbol{ \boldsymbol{\varPhi}}} T_{\varepsilon(\sigma)}(\mathit{\boldsymbol{ \boldsymbol{\varPsi}}} \boldsymbol{x})\right\|_{2}^{2}= \\ \left\|\mathit{\boldsymbol{ \boldsymbol{\varPhi}}} \mathit{\boldsymbol{ \boldsymbol{\varPsi}}} \boldsymbol{x}-\mathit{\boldsymbol{ \boldsymbol{\varPhi}}} T_{\varepsilon(\sigma)}(\mathit{\boldsymbol{ \boldsymbol{\varPsi}}} \boldsymbol{x})\right\|_{2}^{2} \leqslant \\ \|\mathit{\boldsymbol{ \boldsymbol{\varPhi}}}\|_{2}^{2}\left\|\mathit{\boldsymbol{ \boldsymbol{\varPsi}}} \boldsymbol{x}-T_{\varepsilon(\sigma)}(\mathit{\boldsymbol{ \boldsymbol{\varPsi}}} \boldsymbol{x})\right\|_{2}^{2} \end{gather*} $

(2)

由于||Φ||₂²=λ_max(Φ^TΦ)≤τ₀(此处λ_max(·)表示最大特征值)，式(2)的上界可以进一步被定义为

$ \left\|\boldsymbol{x}-\mathit{\boldsymbol{ \boldsymbol{\varPhi}}} T_{\varepsilon(\sigma)}(\mathit{\boldsymbol{ \boldsymbol{\varPsi}}} \boldsymbol{x})\right\|_{2}^{2} \leqslant \tau_{0}\left\|\mathit{\boldsymbol{ \boldsymbol{\varPsi}}} \boldsymbol{x}-T_{\varepsilon(\sigma)}(\mathit{\boldsymbol{ \boldsymbol{\varPsi}}} \boldsymbol{x})\right\|_{2}^{2}. $

(3)

式(3)可以被改写为

$ \left\|\boldsymbol{x}-\mathit{\boldsymbol{ \boldsymbol{\varPhi}}} T_{\varepsilon(\sigma)}(\mathit{\boldsymbol{ \boldsymbol{\varPsi}}} \boldsymbol{x})\right\|_{2}^{2} \leqslant \tau_{0}\|\mathit{\boldsymbol{ \boldsymbol{\varPsi}}} \boldsymbol{x}-T(\mathit{\boldsymbol{ \boldsymbol{\varPsi}}} \boldsymbol{x}, \boldsymbol{e})\|_{2}^{2}. $

(4)

其中：阈值向量e的元素被用于收缩Ψx。令(Ψx)_m和ε_m分别代表Ψx和e的第m个元素，ε_m为非负值，则软阈值算子T[(Ψx)_m, ε_m]可以定义为

$ \begin{array}{*{20}{c}} T\left[(\mathit{\boldsymbol{ \boldsymbol{\varPsi}}} \boldsymbol{x})_{m}, \boldsymbol{\varepsilon}_{m}\right]= \\ \begin{cases}(\mathit{\boldsymbol{ \boldsymbol{\varPsi}}} \boldsymbol{x})_{m}+\varepsilon_{m}, & (\mathit{\boldsymbol{ \boldsymbol{\varPsi}}} \boldsymbol{x})_{m}<-\varepsilon_{m} ; \\ 0, & \left|(\mathit{\boldsymbol{ \boldsymbol{\varPsi}}} \boldsymbol{x})_{m}\right| \leqslant \varepsilon_{m} ; \\ (\mathit{\boldsymbol{ \boldsymbol{\varPsi}}} \boldsymbol{x})_{m}-\varepsilon_{m}, & (\mathit{\boldsymbol{ \boldsymbol{\varPsi}}} \boldsymbol{x})_{m}>\varepsilon_{m} .\end{cases} \end{array} $

收缩过程T(Ψx, e)会发生以下4种情况之一：

C1：Ψx的任意元素满足(Ψx)_m < －ε_m，则T(Ψx, e)=Ψx+e。因此有

$ \|\mathit{\boldsymbol{ \boldsymbol{\varPsi}}} \boldsymbol{x}-T(\mathit{\boldsymbol{ \boldsymbol{\varPsi}}} \boldsymbol{x}, \boldsymbol{e})\|_{2}^{2}=\|\mathit{\boldsymbol{ \boldsymbol{\varPsi}}} \boldsymbol{x}-\mathit{\boldsymbol{ \boldsymbol{\varPsi}}} \boldsymbol{x}-\boldsymbol{e}\|_{2}^{2}=\|\boldsymbol{e}\|_{2}^{2}. $

(5)

C2：Ψx的任意元素满足|(Ψx)_m|≤ε_m，则T(Ψx, e)=0。因此有

$ \|\mathit{\boldsymbol{ \boldsymbol{\varPsi}}} \boldsymbol{x}-T(\mathit{\boldsymbol{ \boldsymbol{\varPsi}}} \boldsymbol{x}, \boldsymbol{e})\|_{2}^{2}=\|\mathit{\boldsymbol{ \boldsymbol{\varPsi}}} \boldsymbol{x}-\mathbf{0}\|_{2}^{2} \leqslant\|\boldsymbol{e}\|_{2}^{2}. $

(6)

C3：Ψx的任意元素满足(Ψx)_m>ε_m，则T(Ψx, e)=Ψx－e。因此有

$ \|\mathit{\boldsymbol{ \boldsymbol{\varPsi}}} \boldsymbol{x}-T(\mathit{\boldsymbol{ \boldsymbol{\varPsi}}} \boldsymbol{x}, \boldsymbol{e})\|_{2}^{2}=\|\mathit{\boldsymbol{ \boldsymbol{\varPsi}}} \boldsymbol{x}-\mathit{\boldsymbol{ \boldsymbol{\varPsi}}} \boldsymbol{x}+\boldsymbol{e}\|_{2}^{2}=\|\boldsymbol{e}\|_{2}^{2}. $

(7)

C4：上述3种情况的任意2种，或C1—C3均发生。

因此当C4情况发生时，||Ψx－T(Ψx, e)||₂²的上界是C1—C3情况下该值的最大上界。基于式(5)—(7)可得

$ \|\mathit{\boldsymbol{ \boldsymbol{\varPsi}}} \boldsymbol{x}-T(\mathit{\boldsymbol{ \boldsymbol{\varPsi}}} \boldsymbol{x}, \boldsymbol{e})\|_{2}^{2} \leqslant\|\boldsymbol{e}\|_{2}^{2}. $

(8)

根据阈值向量的定义e=c⊗m，其中比例常数向量c的每一个元素的范围均为c_l∈[c_min, c_max]。令ε_max=σc_max表示e的最大元素，则有

$ \|\boldsymbol{e}\|_{2}^{2} \leqslant M \varepsilon_{\max }^{2} \leqslant M c_{\max }^{2} \sigma^{2}. $

(9)

根据式(4)、(8)和(9)可得

$ \left\|\boldsymbol{x}-\mathit{\boldsymbol{ \boldsymbol{\varPhi}}} T_{\varepsilon(\sigma)}(\mathit{\boldsymbol{ \boldsymbol{\varPsi}}} \boldsymbol{x})\right\|_{2}^{2} \leqslant \tau_{0} M \sigma^{2} c_{\max }^{2} \leqslant M \sigma^{2} L. $

其中L=τ₀c_max²。

2 基于对偶标架的模型驱动p-CSMRI网络 2.1 优化问题

p-CSMRI的采样模型可以描述为

$ \boldsymbol{y}=\sum\limits_{i=1}^{n_{\mathrm{c}}} \boldsymbol{M F S}_{i} \boldsymbol{x}+\boldsymbol{n}. $

式中，x∈$\mathbb{C}^N$表示原始的复值MR图像，y∈$\mathbb{C}^M$表示部分采样k域数据(M < N)，n∈$\mathbb{C}^M$表示Gauss噪声，M∈$\mathbb{R}^{M \times N}$表示部分采样矩阵，F∈$\mathbb{C}^{N \times N}$表示Fourier变换矩阵，S_i∈$\mathbb{C}^{N \times N}$表示第i个线圈灵敏度，n_c表示线圈总数。

根据最大后验概率准则和Bayes理论，从测量数据y中重建原始图像x的优化模型可以表示为

$ \min\limits_{\boldsymbol{x}} \frac{1}{2} \sum\limits_{i=1}^{n_{\mathrm{c}}}\left\|\boldsymbol{M F S}_{i} \boldsymbol{x}-\boldsymbol{y}_{i}\right\|_{2}^{2}+\lambda R_{\text {dual }}(\boldsymbol{x}). $

(10)

式中，y_i∈$\mathbb{C}^M$表示第i个接收线圈的部分采样k域数据。第1项是数据保真项，该项保证了重建图像与测量数据的一致性；第2项是正则化项，本文通过DNN从大量的数据中学习图像先验分布，进而对图像施加约束；λ用以控制这2项的比重。

2.2 模型求解

为求解式(10)定义的优化问题，引入辅助变量z∈$\mathbb{C}^N$和{m_i∈$\mathbb{C}^N$}_i=1^n_c，将式(10)改写为

$ \left\{\begin{array}{l} \min\limits_{z, \boldsymbol{m}_i, \boldsymbol{x}} \frac{1}{2} \sum\limits_{i=1}^{n_{\mathrm{c}}}\left\|\boldsymbol{M F m}_i-\boldsymbol{y}_i\right\|_2^2+\lambda R_{\text {dual }}(\boldsymbol{z}) \\ \text { s. t. } \boldsymbol{m}_i=\boldsymbol{S}_i \boldsymbol{x}, \boldsymbol{z}=\boldsymbol{x} . \end{array}\right. $

(11)

利用半二次分裂算法求解式(11)定义的优化问题，可将式(11)转化为

$ \begin{array}{*{20}{c}} & \min\limits_{z, \boldsymbol{m}_{i}, \boldsymbol{x}} \frac{1}{2} \sum\limits_{i=1}^{n_{\mathrm{c}}}\left\|\boldsymbol{M F} \boldsymbol{m}_{i}-\boldsymbol{y}_{i}\right\|_{2}^{2}+\lambda R_{\text {dual }}(\boldsymbol{z})+ \\ & \frac{\alpha}{2} \sum\limits_{i=1}^{n_{\mathrm{c}}}\left\|\boldsymbol{m}_{i}-\boldsymbol{S}_{i} \boldsymbol{x}\right\|_{2}^{2}+\frac{\beta}{2}\|\boldsymbol{z}-\boldsymbol{x}\|_{2}^{2} . \end{array} $

(12)

其中α和β为惩罚参数。

式(12)通过交替优化z、m、x进行求解，在第t次迭代中，交替迭代式(13)—(15)。

$ \boldsymbol{z}^{(t)}=\underset{z}{\operatorname{argmin}}\left\{\frac{\beta}{2}\left\|\boldsymbol{z}-\boldsymbol{x}^{(t-1)}\right\|_{2}^{2}+\lambda R_{\text {dual }}(\boldsymbol{z})\right\} . $

(13)

$ \begin{array}{*{20}{c}} \boldsymbol{m}_{i}^{(t)}= & \underset{\boldsymbol{m}_{i}}{\operatorname{argmin}}\left\{\frac{1}{2} \sum\limits_{i=1}^{n_{{\rm{c}}}}\left\|\boldsymbol{M F m}_{i}-\boldsymbol{y}_{i}\right\|_{2}^{2}+\right. \\ & \left.\frac{\alpha}{2} \sum\limits_{i=1}^{n_{{\rm{c}}}}\left\|\boldsymbol{m}_{i}-\boldsymbol{S}_{i} \boldsymbol{x}^{(t-1)}\right\|_{2}^{2}\right\} .\\ & \boldsymbol{x}^{(t)}=\underset{\boldsymbol{x}}{\operatorname{argmin}}\left\{\frac{\alpha}{2} \sum\limits_{i=1}^{n_{c}} \| \boldsymbol{m}_{i}^{(t)}-\right. \end{array} $

(14)

$ \left.\boldsymbol{S}_{i} \boldsymbol{x}\left\|_{2}^{2}+\frac{\beta}{2}\right\| \boldsymbol{z}^{(t)}-\boldsymbol{x} \|_{2}^{2}\right\} . $

(15)

具体步骤如下：

1) 固定m_i和x，更新辅助变量z

式(13)可通过满足有界条件的对偶标架去噪器进行近似求解。

$ \boldsymbol{z}^{(t)}=\mathit{\boldsymbol{ \boldsymbol{\varPhi}}} T_{\varepsilon(\sigma)}\left[\mathit{\boldsymbol{ \boldsymbol{\varPsi}}}\left(\boldsymbol{x}^{(t-1)}\right)\right]. $

(16)

由于估计图像包含的噪声随着迭代次数增加而降低，因此输入到去噪器的噪声强度也应逐渐降低。本文给定σ的初始值，并按照σ^(t)=ξσ^(t－1) (0 < ξ < 1)进行更新。为了满足这样的变化趋势，在实验中给定σ的初始值σ_init和最终值σ_final，使其随着展开次数的增长呈指数级衰减。此外，由于本文最后将迭代算法展开成网络，因此将求解辅助变量z所对应的网络称为先验模块。

2) 固定z和x，更新辅助变量m_i

式(14)的目标函数包括2个关于m_i的简单二次项，可对其求导并令导数等于0，对应的闭式解满足

$ \left(\boldsymbol{F}^{H} \boldsymbol{M}^{H} \boldsymbol{M} \boldsymbol{F}+\alpha \boldsymbol{I}\right) \boldsymbol{m}_{i}^{(t)}=\boldsymbol{F}^{H} \boldsymbol{M}^{H} \boldsymbol{y}_{i}+\alpha \boldsymbol{S}_{i} \boldsymbol{x}^{(t-1)}. $

其中F作为一个酉矩阵，满足F^HF=FF^H=I。因此上述方程可以写为

$ \left(\boldsymbol{M}^{H} \boldsymbol{M}+\alpha \boldsymbol{I}\right) \boldsymbol{F} \boldsymbol{m}_{i}^{(t)}=\boldsymbol{M}^{H} \boldsymbol{y}_{i}+\alpha \boldsymbol{F} \boldsymbol{S}_{i} \boldsymbol{x}^{(t-1)} . $

即

$ \boldsymbol{F m}_{i}^{(t)}=\left(\boldsymbol{M}^{H} \boldsymbol{M}+\alpha \boldsymbol{I}\right)^{-1}\left(\boldsymbol{M}^{H} \boldsymbol{y}_{i}+\alpha \boldsymbol{F} \boldsymbol{S}_{i} \boldsymbol{x}^{(t-1)}\right) . $

(17)

其中：M^HM表示N×N的对角矩阵，其对角线上的元素为0和1，值为1的位置表示k域中对应变换系数被采样，值为0的位置表示未被采样；M^Hy_i表示ZF重建的k域测量数据。

令(k_u, k_v)表示k域的二维索引值，Ω表示被采样到的像素点索引值的集合，则在k域中m_i的更新可表示为

$ \boldsymbol{F m}_{i}^{(t)}\left(k_{u}, k_{v}\right)=\left\{\begin{array}{l} P^{(t)}\left(k_{u}, k_{v}\right), \left(k_{u}, k_{v}\right) \notin \mathit{\boldsymbol{ \boldsymbol{\varOmega}}} ; \\ \frac{P^{(0)}\left(k_{u}, k_{v}\right)+\alpha P^{(t)}\left(k_{u}, k_{v}\right)}{1+\alpha}, \\ \left(k_{u}, k_{v}\right) \in \mathit{\boldsymbol{ \boldsymbol{\varOmega}}} . \end{array}\right. $

其中P⁽⁰⁾=M^Hy_i，P^(t)=FS_ix^(t－1)。对式(17)进行Fourier反变换，可以得到m_i^(t)。

3) 固定z，m_i，更新图像x

式(15)同样也包含2个关于x的二次项，对代价函数求导并令导数等于0得到其闭式解为

$ \boldsymbol{x}^{(t)}=\left(\beta \boldsymbol{I}+\sum\limits_{i=1}^{n_{{\rm{c}}}} \alpha \boldsymbol{S}_{i}^{H} \boldsymbol{S}_{i}\right)^{-1}\left(\beta \boldsymbol{z}^{(t)}+\alpha \sum\limits_{i=1}^{n_{{\rm{c}}}} \boldsymbol{S}_{i}^{H} \boldsymbol{m}_{i}^{(t)}\right) $

(18)

通过交替更新图像x、辅助变量m_i和z这3个步骤求解式(12)以优化模型。为获得更高的重建质量，本文将上述迭代过程展开成网络，每次迭代对应深度展开网络的一个阶段(Stage)，且每一个阶段都包含上述每一步的更新过程，网络总的展开次数即总的阶段数记为K。将所提算法记为CSMRI-Dual算法，其具体实现步骤如算法1所示。

算法1  CSMRI-Dual算法.

1) 输入：测量数据y，线圈灵敏度S_i，掩膜M，噪声标准差的初始值σ_init和最终值σ_final

2) FOR t=1∶K

3) 根据式(16)更新辅助变量z

4) 根据式(17)更新辅助变量m_i

5) 根据式(18)更新图像x

6) 结束

7) 输出：重建图像$\hat{\mathit{\boldsymbol{x}}}$

算法具体的数据流程图如图 1所示。图中展示了第t个阶段的具体架构，包含辅助变量m_i更新模块、辅助变量z更新模块(先验模块)和图像x更新模块。

2.3 CSMRI-Dual算法的收敛性分析

收敛的p-CSMRI算法有益于算法稳定且获得较高的重建质量。基于先验模块的有界性，本文分析当K→∞时，即迭代次数趋于无限次的时候，算法是否收敛。定理1表明所提出的CSMRI-Dual算法是收敛的。

定理1  对于p-CSMRI的采样模型y=$\sum\limits_{i=1}^{n_ {\rm{c}}}$MFS_ix+n，给定{y，$\sum\limits_{i=1}^{n_{\mathrm{c}}}$MFS_i}，x通过CSMRI-Dual算法进行迭代求解，那么x^(t)收敛到一个不动点。也就是说，存在x^*，当t→∞时，使得||x^(t)－x^*||₂→0。

证明：对于所有t，令Θ表示x^(t)的域。在域Θ中，定义一个距离函数d：Θ×Θ→$\mathbb{R}$，满足

$ d\left(\boldsymbol{x}^{(t)}, \boldsymbol{x}^{(j)}\right)=\frac{1}{\sqrt{N}}\left\|\boldsymbol{x}^{(t)}-\boldsymbol{x}^{(j)}\right\|_{2}. $

Θ⊆$\mathbb{R}^N$和$\mathbb{R}^N$是一个完备的度量空间。因此，只要证明{x^(t)}_t=1^∞在Θ中是一个具有距离函数d的Cauchy序列，则x^(t)是收敛的。

根据上述定义，可构建

$ \begin{array}{*{20}{c}} {\left\|\boldsymbol{x}^{(t)}-\boldsymbol{x}^{(t-1)}\right\|_{2}=\left\|\boldsymbol{x}^{(t)}-\boldsymbol{z}^{(t)}+\boldsymbol{z}^{(t)}-\boldsymbol{x}^{(t-1)}\right\|_{2} \leqslant}\\ {\left\|\boldsymbol{x}^{(t)}-\boldsymbol{z}^{(t)}\right\|_{2}+\left\|\boldsymbol{z}^{(t)}-\boldsymbol{x}^{(t-1)}\right\|_{2}}. \end{array} $

(19)

由于z^(t)=D_ε[x^(t－1); σ]，且根据引理1中可训练对偶标架去噪器的边界，可得

$ \begin{gather*} \left\|\boldsymbol{z}^{(t)}-\boldsymbol{x}^{(t-1)}\right\|_{2}=\| D_{\varepsilon}\left[\boldsymbol{x}^{(t-1)} ; \sigma\right]- \\ \boldsymbol{x}^{(t-1)} \|_{2} \leqslant \sqrt{M L} \sigma^{(t-1)} . \end{gather*} $

(20)

将式(20)代入式(19)，可得

$ \begin{gather*} \left\|\boldsymbol{x}^{(t)}-\boldsymbol{x}^{(t-1)}\right\|_{2} \leqslant\left\|\boldsymbol{x}^{(t)}-\boldsymbol{z}^{(t)}\right\|_{2}+ \\ \left\|\boldsymbol{z}^{(t)}-\boldsymbol{x}^{(t-1)}\right\|_{2} \leqslant\left\|\boldsymbol{x}^{(t)}-\boldsymbol{z}^{(t)}\right\|_{2}+\sqrt{M L} \sigma^{(t-1)} . \end{gather*} $

(21)

令

$ f(\boldsymbol{x})=\frac{\alpha}{2} \sum\limits_{i=1}^{n_{{\rm{c}}}}\left\|\boldsymbol{m}_{i}^{(t)}-\boldsymbol{S}_{i} \boldsymbol{x}\right\|_{2}^{2} . $

(22)

则式(15)可改写为

$ \boldsymbol{x}^{(t)}=\underset{\boldsymbol{x}}{\operatorname{argmin}}\left\{f(\boldsymbol{x})+\frac{\beta}{2}\left\|\boldsymbol{z}^{(t)}-\boldsymbol{x}\right\|_{2}^{2}\right\}. $

(23)

对式(23)定义的代价函数求导并令导数等于0，可以得到

$ \nabla f(\boldsymbol{x})+\beta\left[\boldsymbol{x}-\boldsymbol{z}^{(t)}\right]=0 . $

即

$ \boldsymbol{x}-\boldsymbol{z}^{(t)}=-\frac{\nabla f(\boldsymbol{x})}{\beta}. $

对于第t次迭代，可得

$ \boldsymbol{x}^{(t)}-\boldsymbol{z}^{(t)}=-\frac{\nabla f\left(\boldsymbol{x}^{(t)}\right)}{\beta}. $

即

$ \left\|\boldsymbol{x}^{(t)}-\boldsymbol{z}^{(t)}\right\|_{2}=\frac{\left\|\nabla f\left(\boldsymbol{x}^{(t)}\right)\right\|_{2}}{\beta} . $

(24)

由式(22)可以求得，f(x^(t))在点x^(t)的梯度▽f(x^(t))为

$ \nabla f\left(\boldsymbol{x}^{(t)}\right)=\alpha \sum\limits_{i=1}^{n_{\mathrm{c}}} \boldsymbol{S}_{i}^{H}\left[\boldsymbol{S}_{i} \boldsymbol{x}^{(t)}-\boldsymbol{m}_{i}^{(t)}\right]. $

因此，

$ \begin{gather*} \left\|\nabla f\left(\boldsymbol{x}^{(t)}\right)\right\|_{2}=\left\|\alpha \sum\limits_{i=1}^{n_{{\rm{c}}}} \boldsymbol{S}_{i}^{H}\left[\boldsymbol{S}_{i} \boldsymbol{x}^{(t)}-\boldsymbol{m}_{i}^{(t)}\right]\right\|_{2} \leqslant \\ \alpha \sum\limits_{i=1}^{n_{{\rm{c}}}}\left\|\boldsymbol{S}_{i}^{H}\right\|_{2}\left\|\boldsymbol{S}_{i} \boldsymbol{x}^{(t)}-\boldsymbol{m}_{i}^{(t)}\right\|_{2} \leqslant \\ \alpha \sum\limits_{i=1}^{n_{{\rm{c}}}}\left\|\boldsymbol{S}_{i}^{H}\right\|_{2}\left\|\boldsymbol{S}_{i}\right\|_{2}\left\|\boldsymbol{x}^{(t)}\right\|_{2}+ \\ \alpha \sum\limits_{i=1}^{n_{{\rm{c}}}}\left\|\boldsymbol{S}_{i}^{H}\right\|_{2}\left\|\boldsymbol{m}_{i}^{(t)}\right\|_{2} . \end{gather*} $

(25)

将式(17)代入式(25)，可得

$ \begin{array}{*{20}{c}} & \left\|\nabla f\left(\boldsymbol{x}^{(t)}\right)\right\|_{2} \leqslant \alpha \sum\limits_{i=1}^{n_{\mathrm{c}}}\left\|\boldsymbol{S}_{i}^{H}\right\|_{2}\left\|\boldsymbol{S}_{i}\right\|_{2}\left\|\boldsymbol{x}^{(t)}\right\|_{2}+ \\ & \alpha \sum\limits_{i=1}^{n_{{\rm{c}}}}\left\|\boldsymbol{S}_{i}^{H}\right\|_{2} \| \boldsymbol{F}^{H}\left[\left(\boldsymbol{M}^{H} \boldsymbol{M}+\alpha I\right)^{-1} \cdot\right. \\ & \left.\left(\boldsymbol{M}^{H} \boldsymbol{y}_{i}+\alpha \boldsymbol{F} \boldsymbol{S}_{i} \boldsymbol{x}^{(t-1)}\right)\right] \|_{2} \leqslant \\ & \alpha \sum\limits_{i=1}^{n_{\mathrm{c}}}\left\|\boldsymbol{S}_{i}^{H}\right\|_{2}\left\|\boldsymbol{S}_{i}\right\|_{2}\left\|\boldsymbol{x}^{(t)}\right\|_{2}+ \\ & \alpha\left\|\left(\boldsymbol{M}^{H} \boldsymbol{M}+\alpha \boldsymbol{I}\right)^{-1}\right\|_{2} \sum\limits_{i=1}^{n_{\mathrm{c}}}\left\|\boldsymbol{S}_{i}^{H}\right\|_{2} \\ & \left\|\boldsymbol{F}^{H}\left[\boldsymbol{M}^{H} \boldsymbol{y}_{i}+\alpha \boldsymbol{F} \boldsymbol{S}_{i} \boldsymbol{x}^{(t-1)}\right]\right\|_{2} . \end{array} $

(26)

其中，M^HM是一个N×N的对角矩阵，其对角线元素为0或1，所以M^HM+αI是一个N×N的对角矩阵，其对角线元素为α或α+1。由此可推导出

$ \begin{array}{*{20}{c}} {\left\|\left(\boldsymbol{M}^{H} \boldsymbol{M}+\alpha \boldsymbol{I}\right)^{-1}\right\|_{2}=}\\ {\sqrt{\lambda_{\max }\left\{\left[\left(\boldsymbol{M}^{H} \boldsymbol{M}+\alpha \boldsymbol{I}\right)^{-1}\right]^{H}\left[\left(\boldsymbol{M}^{H} \boldsymbol{M}+\alpha \boldsymbol{I}\right)^{-1}\right]\right\}}=}\\ {\sqrt{\left(\frac{1}{\alpha}\right)^{2}}=\frac{1}{\alpha} .} \end{array} $

(27)

式(27)代入式(26)可得

$ \begin{array}{*{20}{c}} \left\|\nabla f\left(\boldsymbol{x}^{(t)}\right)\right\|_{2} \leqslant \alpha \sum\limits_{i=1}^{n_{{\rm{c}}}}\left\|\boldsymbol{S}_{i}^{H}\right\|_{2}\left\|\boldsymbol{S}_{i}\right\|_{2}\left\|\boldsymbol{x}^{(t)}\right\|_{2}+ \\ \alpha\left\|\left(\boldsymbol{M}^{H} \boldsymbol{M}+\alpha I\right)^{-1}\right\|_{2} \sum\limits_{i=1}^{n_{\mathrm{c}}}\left\|\boldsymbol{S}_{i}^{H}\right\|_{2} \cdot \\ \left\|\boldsymbol{F}^{H}\left[\boldsymbol{M}^{H} \boldsymbol{y}_{i}+\alpha \boldsymbol{F} \boldsymbol{S}_{i} \boldsymbol{x}^{(t-1)}\right]\right\|_{2} \leqslant \\ \alpha \sum\limits_{i=1}^{n_{\mathrm{c}}}\left\|\boldsymbol{S}_{i}^{H}\right\|_{2}\left\|\boldsymbol{S}_{i}\right\|_{2}\left\|\boldsymbol{x}^{(t)}\right\|_{2}+\\ \sum\limits_{i=1}^{n_{{\rm{c}}}}\left\|\boldsymbol{S}_{i}^{H}\right\|_{2}\left\|\boldsymbol{F}^{H} \boldsymbol{M}^{H} \boldsymbol{y}_{i}+\alpha \boldsymbol{S}_{i} \boldsymbol{x}^{(t-1)}\right\|_{2} \leqslant \\ \alpha \sum\limits_{i=1}^{n_{{\rm{c}}}}\left\|\boldsymbol{S}_{i}^{H}\right\|_{2}\left\|\boldsymbol{S}_{i}\right\|_{2}\left\|\boldsymbol{x}^{(t)}\right\|_{2}+ \\ \sum\limits_{i=1}^{n_{{\rm{c}}}}\left\|\boldsymbol{S}_{i}^{H}\right\|_{2}\left\|\boldsymbol{F}^{H} \boldsymbol{M}^{H}\right\|_{2}\left\|\boldsymbol{y}_{i}\right\|_{2}+ \\ \alpha \sum\limits_{i=1}^{n_{{\rm{c}}}}\left\|\boldsymbol{S}_{i}^{H}\right\|_{2}\left\|\boldsymbol{S}_{i}\right\|_{2}\left\|\boldsymbol{x}^{(t-1)}\right\|_{2}. \end{array} $

(28)

其中x^(t)与x^(t－1)分别表示第t次与第t—1次迭代过程中的重建图像，其大小被限制为

$ \|\boldsymbol{x}\|_{2} \leqslant C_{0}. $

(29)

考虑到采样模型y_i=MFS_ix+n中||n||₂≤κ，且||MFS_ix||₂≤||FS_ix||₂=||S_ix||₂≤||S_i||₂||x||₂。因此，

$ \begin{array}{*{20}{c}} \left\|\boldsymbol{y}_{i}\right\|_{2}=\left\|\boldsymbol{M F S}_{i} \boldsymbol{x}+\boldsymbol{n}\right\|_{2} \leqslant\left\|\boldsymbol{M F S}_{i} \boldsymbol{x}\right\|_{2}+ \\ \|\boldsymbol{n}\|_{2} \leqslant\left\|\boldsymbol{S}_{i}\right\|_{2}\|\boldsymbol{x}\|_{2}+ \\ \|\boldsymbol{n}\|_{2} \leqslant C_{0}\left\|\boldsymbol{S}_{i}\right\|_{2}+\kappa . \end{array} $

(30)

由于||F^HM^H||₂=$\sqrt{\lambda_{\max }\left(\boldsymbol{M F} \boldsymbol{F}^H \boldsymbol{M}^H\right)}$，其中λ_max(MFF^HM^H)代表矩阵MFF^HM^H的最大特征值，M∈$\mathbb{R}^{M \times N}$表示采样矩阵，其元素为1的位置表示采样点的对应索引。因此MFF^HM^H=MM^H=I，其中I∈$\mathbb{R}^{M \times M}$是一个单位矩阵。因此，可以得到λ_max(MFF^HM^H)=1。所以有

$ \left\|\boldsymbol{F}^{H} \boldsymbol{M}^{H}\right\|_{2}=\sqrt{\lambda_{\max }\left(\boldsymbol{M F} \boldsymbol{F}^{H} \boldsymbol{M}^{H}\right)}=1 . $

(31)

将式(29)、(30)和(31)代入式(28)中，可得

$ \begin{array}{*{20}{c}} \left\|\nabla f\left(\boldsymbol{x}^{(t)}\right)\right\|_{2} \leqslant \alpha \sum\limits_{i=1}^{n_{{\rm{c}}}}\left\|\boldsymbol{S}_{i}^{H}\right\|_{2}\left\|\boldsymbol{S}_{i}\right\|_{2}\left\|\boldsymbol{x}^{(t)}\right\|_{2}+ \\ \sum\limits_{i=1}^{n_{{\rm{c}}}}\left\|\boldsymbol{S}_{i}^{H}\right\|_{2}\left\|\boldsymbol{F}^{H} \boldsymbol{M}^{H}\right\|_{2}\left\|\boldsymbol{y}_{i}\right\|_{2}+\\ {\alpha \sum\limits_{i=1}^{n_{{\rm{c}}}}\left\|\boldsymbol{S}_{i}^{H}\right\|_{2}\left\|\boldsymbol{S}_{i}\right\|_{2}\left\|\boldsymbol{x}^{(t-1)}\right\|_{2} \leqslant C_{1}+C_{2}+C_{3}.} \end{array} $

其中C₁=αC₀$\sum\limits_{i=1}^{n_{\mathrm{c}}}$||S_i^H||₂||S_i||₂，C₂=$\sum\limits_{i=1}^{n_{\mathrm{c}}}$||S_i^H||₂(C₀||S_i||₂+κ)，C₃=αC₀$\sum\limits_{i=1}^{n_{\mathrm{c}}}$||S_i^H||₂||S_i||₂。即||f(x^(t))||₂存在一个上界。根据设定的参数更新准则β^(t)=γβ^(t－1) (γ>1)，式(24)可以写为

$ \begin{gather*} \left\|\boldsymbol{x}^{(t)}-\boldsymbol{z}^{(t)}\right\|_{2}=\frac{\left\|\nabla f\left(\boldsymbol{x}^{(t)}\right)\right\|_{2}}{\beta} \leqslant \\ \frac{C_{1}+C_{2}+C_{3}}{\gamma^{t-1} \beta^{(0)}}=\frac{C_{4}}{\gamma^{t-1}} . \end{gather*} $

(32)

其中C₄=$\frac{C_1+C_2+C_3}{\beta^{(0)}}$。

根据输入噪声标准差的更新准则σ^(t)=ξσ^(t－1)(0 < ξ < 1)并将式(32)代入式(21)，可得

$ \begin{gather*} \left\|\boldsymbol{x}^{(t)}-\boldsymbol{x}^{(t-1)}\right\|_{2} \leqslant\left\|\boldsymbol{x}^{(t)}-\boldsymbol{z}^{(t)}\right\|_{2}+\sqrt{M L} \sigma^{(t-1)} \leqslant \\ \frac{C_{4}}{\gamma^{t-1}}+\sqrt{M L} \xi^{t-1} \sigma^{(0)} \leqslant \frac{C_{4}}{\gamma^{t-1}}+C_{5} \xi^{t-1} . \end{gather*} $

(33)

其中C₅=$\sqrt{M L} \sigma^{(0)}$。

综上所述，对于任意t和常数C₄、C₅，{x^(t)}_t=1^∞序列总满足

$ d\left(\boldsymbol{x}^{(t)}, \boldsymbol{x}^{(t-1)}\right) \leqslant \frac{C_{4}}{\gamma^{t-1}}+C_{5} \xi^{t-1}(\gamma>1, 0<\xi<1). $

因此，对于任意有限的j和t(j>t)，总满足

$ \begin{array}{*{20}{c}} {d\left(\boldsymbol{x}^{(j)}, \boldsymbol{x}^{(t)}\right) \leqslant \sum\limits_{n=t+1}^{j} \frac{C_{4}}{\gamma^{n-1}}+\sum\limits_{n=t+1}^{j} C_{5} \xi^{n-1}=}\\ {C_{4} \frac{\gamma^{j-t}-1}{\gamma^{j-1}(\gamma-1)}+C_{5} \xi^{t} \frac{1-\xi^{j-t}}{1-\xi} .} \end{array} $

而

$ \begin{aligned} & \lim _{\substack{j \rightarrow \infty \\ t \rightarrow \infty \\ j>t}} C_{4} \frac{\gamma^{j-t}-1}{\gamma^{j-1}(\gamma-1)}+C_{5} \xi^{t} \frac{1-\xi^{j-t}}{1-\xi}= \\ & \lim _{\substack{j \rightarrow \infty \\ t>\infty \\ j>t}} \frac{1}{\gamma^{t-1}(\gamma-1)}-\frac{1}{\gamma^{j-1}(\gamma-1)}+ \\ & \xi^{t}-\xi^{j}=0(\gamma>1, 0<\xi<1) . \end{aligned} $

由于j>t且γ>1，0 < ξ < 1, 所以t→∞时，d(x^(j), x^(t))→0。因此{x^(t)}_t=1^∞是一个Cauchy序列。又因为在完备空间$\mathbb{R}^N$上的Cauchy序列是收敛的，所以一定存在x^*，使得

$ d\left(\boldsymbol{x}^{(t)}, \boldsymbol{x}^{*}\right) \rightarrow 0 . $

因此可得||x^(t)－x^*||₂→0。

推论1  在定理1中，本文分析了当通过先验网络学习到的标架是满足严格的对偶性质时，算法会收敛到一个不动点。然而，通过损失函数对标架进行的对偶约束，经验上并不是严格对偶。在近似对偶约束的条件下，对偶约束项不是0，但有一个上界，即||ΦΨ－I||_F²≤τ₁ (此处τ₁是一个很小的正值)，此时随着迭代次数增加，即t→∞，有||x^(t)－x^(t－1)||₂→C₆。

证明：由于||ΦΨ－I||_F²≤τ₁，进而可得||ΦΨx－x||₂²≤τ₂，则有下式：

$ \begin{gathered} \left\|\boldsymbol{x}-\mathit{\boldsymbol{ \boldsymbol{\varPhi}}} T_{\varepsilon(\sigma)}(\mathit{\boldsymbol{ \boldsymbol{\varPsi}}} \boldsymbol{x})\right\|_{2}^{2}= \\ \left\|\boldsymbol{x}-\mathit{\boldsymbol{ \boldsymbol{\varPhi}}} \mathit{\boldsymbol{ \boldsymbol{\varPsi}}} \boldsymbol{x}+\mathit{\boldsymbol{ \boldsymbol{\varPhi}}} \mathit{\boldsymbol{ \boldsymbol{\varPsi}}} \boldsymbol{x}-\mathit{\boldsymbol{ \boldsymbol{\varPhi}}} T_{\varepsilon(\sigma)}(\mathit{\boldsymbol{ \boldsymbol{\varPsi}}} \boldsymbol{x})\right\|_{2}^{2} \leqslant \\ \|\boldsymbol{x}-\mathit{\boldsymbol{ \boldsymbol{\varPhi}}} \mathit{\boldsymbol{ \boldsymbol{\varPsi}}} \boldsymbol{x}\|_{2}^{2}+\left\|\mathit{\boldsymbol{ \boldsymbol{\varPhi}}} \mathit{\boldsymbol{ \boldsymbol{\varPsi}}} \boldsymbol{x}-\mathit{\boldsymbol{ \boldsymbol{\varPhi}}} T_{\varepsilon(\sigma)}(\mathit{\boldsymbol{ \boldsymbol{\varPsi}}} \boldsymbol{x})\right\|_{2}^{2} \leqslant \\ \tau_{2}+\|\mathit{\boldsymbol{ \boldsymbol{\varPhi}}}\|_{2}^{2}\left\|\mathit{\boldsymbol{ \boldsymbol{\varPsi}}} \boldsymbol{x}-T_{\varepsilon(\sigma)}(\mathit{\boldsymbol{ \boldsymbol{\varPsi}}} \boldsymbol{x})\right\|_{2}^{2} . \end{gathered} $

根据引理1的证明，可证明下式成立：

$ \left\|\boldsymbol{x}-\mathit{\boldsymbol{ \boldsymbol{\varPhi}}} T_{\varepsilon(\sigma)}(\mathit{\boldsymbol{ \boldsymbol{\varPsi}}} \boldsymbol{x})\right\|_{2}^{2} \leqslant \tau_{2}+M \sigma^{2} L . $

在近似对偶约束的条件下，将式(33)改写，可得

$ \left\|\boldsymbol{x}^{(t)}-\boldsymbol{x}^{(t-1)}\right\|_{2} \leqslant \frac{C_{4}}{\gamma^{t-1}}+\left(\sqrt{M L} \sigma^{(t-1)}+\sqrt{\tau_{2}}\right) . $

(34)

根据输入噪声标准差的更新准则σ^(t)=ξσ^(t－1) (0 < ξ < 1)，式(34)可以改写为

$ \begin{gathered} \left\|\boldsymbol{x}^{(t)}-\boldsymbol{x}^{(t-1)}\right\|_{2} \leqslant \frac{C_{4}}{\gamma^{t-1}}+\left(\sqrt{M L} \xi^{t-1} \sigma^{(0)}+\sqrt{\tau_{2}}\right) \leqslant \\ \frac{C_{4}}{\gamma^{t-1}}+C_{5} \xi^{t-1}+C_{6} . \end{gathered} $

其中，C₆=$\sqrt{\tau_2}$。随着迭代次数增加，即t→∞，最终可以得到||x^(t)－x^(t－1)||₂→C₆。

3 实验结果与分析 3.1 参数设置

本文从纽约大学的膝盖磁共振数据集(http://mridata.org)中分别随机选取了500、100、100张不同的图片作为训练、验证和测试数据集。原始数据由15个通道的k域数据组成，尺寸大小为320×320。本文通过ESPIRiT^[27]算法，从全采样k空间中心的大小为30×30的数据块中预先计算线圈灵敏度图。采用2倍、4倍和6倍加速因子(acceleration factor，AF)对原始k空间数据进行了Déscartes采样，对于不同加速因子采用相同的全采样k空间中心进行估计。训练的轮数(Epoch)设置为100，使用ADAM作为优化器，学习率初始值设置为0.000 5，每经过50个Epoch的优化，学习率衰减10%。损失函数定义为

$ \text { Loss }=\left\|\boldsymbol{x}_{\mathrm{gt}}-\hat{\boldsymbol{x}}\right\|_{2}^{2}+\frac{\eta}{K G} \sum\limits_{t=1}^{K}\|\mathit{\boldsymbol{ \boldsymbol{\varPhi}}} \mathit{\boldsymbol{ \boldsymbol{\varPsi}}}-\boldsymbol{I}\|_{\mathrm{F}}^{2} . $

式中，$\hat{\mathit{\boldsymbol{x}}}$表示CSMRI-Dual算法的输出图像；x_gt表示参考图像，即线圈灵敏度加权的全采样图像$\sum\limits_i^{n_ {\rm{c}}}$S_i^HF^－1f_i；f_i为第i个线圈的全采样k域数据；η是正则化参数；G表示标架滤波器的大小，具体为121；K设置为8次。其他参数设置为：α=1，β=0.5，λ=0.1，η=1×10^－4，σ_init=70，σ_final=30。所有实验均在Intel(R) Core(TM) i7-9700KF CPU、主频3.60 GHz、内存64 GB、NVIDIA GTX 2080ti GPU的平台上进行。

3.2 不同倍数加速因子下多种p-CSMRI算法的比较

针对100幅测试图像在不同倍数(2、4、6)AF情况下进行图像重建，将本文提出的CSMRI-Dual算法与Modl^[13]、VN^[15]和VS-Net^[14]算法在模型参数量、平均每张图片的推断时间上进行比较。如表 1所示，在参数量方面，VN算法采用梯度下降方法结合变分模型进行图像重建，算法简单，参数量最少；CSMRI-Dual算法在提升算法性能与模型稳定性的前提下具有较高的模型参数量，这是因为深度阈值网络和对偶紧标架采用非共享的方式，如何减小这些参数量且保证算法的重建性能有待进一步研究。在推断时间方面，CSMRI-Dual算法用时少于Modl算法与VN算法，略多于VS-Net算法，计算复杂度适中；Modl算法用时最长，这是因为其展开层数为10，迭代次数最多。

为了对算法重建性能进行客观评估，本文采用了多种指标评价，包括：平均峰值信噪比(peak signal-to-noise ratio，PSNR)，结构相似性指数(structure similarity index measure，SSIM)和归一化均方误差(normalized mean squared error，NMSE)，结果如表 2所示。由表可知，通过Modl算法重建的图像的SSIM在所有采样率的情况下都比ZF重建法要低，但是其PSNR和NMSE远高于ZF重建法，这是因为Modl算法首次将深度展开技术应用于p-CSMRI算法中，网络架构较简单。VN算法使用变分模型结合梯度下降算法，能够获得优于Modl算法的重建质量，但仍略逊于VS-Net算法。VS-Net算法采用交替优化的方式更新辅助变量和重建图像，能够获得优于VN算法和VS-Net算法的重建质量，但参数量在上述3种对比算法中最大。CSMRI-Dual算法由于使用了基于对偶标架的先验网络且可证明收敛，能够取得最优的重建质量，在4倍AF情况下，其图像重建的平均PSNR比Modl算法、VN算法和VS-Net算法分别有1.70、1.45和0.46 dB的提升。

为进一步评估本文提出的CSMRI-Dual算法的重建效果，图 2和图 3分别展示了在4倍和6倍AF时，不同算法重建图像的区域放大图。由图可知，ZF重建图像非常模糊，并且图像过于平滑，缺少细节信息；Modl算法恢复出的重建图像同样丢失了很多的细节信息；VN算法和VS-Net算法重建图像类似，比Modl算法恢复出的重建图像有明显提升，但是细节信息依旧有部分丢失；CSMRI-Dual算法重建图像在4倍AF下视觉效果较好，有效消除了噪声，获得了更多细节，且保留了丰富的纹理信息；在6倍AF下，重建细节更接近原始图像，重建图像较为清晰。由此可见，CSMRI-Dual算法表现出了良好的重建效果。

4 结论

本文基于可训练对偶标架提出了可证明收敛的模型驱动并压缩感知磁共振成像算法。针对现有的模型驱动p-CSMRI网络，其先验网络架构不可解释，算法收敛性难以分析的问题，本文构建了基于可训练对偶标架的深度有界去噪器，并提出了在噪声强度逐渐下降的条件下可证明收敛的并行压缩感知磁共振成像算法框架。本文将满足有界条件的去噪器作为先验模块引入到所构建的模型驱动并行压缩感知磁共振成像网络中，在理论上证明了所提出的深度去噪器满足有界条件，分析了提出模型驱动网络的收敛性，保证了网络的稳定性。仿真实验表明，提出的算法能够取得高质量的MR图像重建。