基于加权核范数正则化的单颗粒冷冻电镜联合优化重构

Joint-optimization reconstruction of single-particle cryo-EM based on weighted nuclear norm regularization

摘要

单颗粒冷冻电镜（cryo-Electron Microscopy，cryo-EM）三维图像重构利用大量未知方向的二维投影图像来重建生物大分子的三维结构。现有的重构方法多将投影方向估计与三维图像重构分开处理，容易出现信息丢失及过平滑问题。本文提出了一种基于三维加权核范数的联合优化重构算法。在分子密度图更新阶段，算法引入了加权核范数正则化方法，采用交替方向乘子法进行求解。在投影方向更新阶段，算法采用块坐标下降法直接对连续空间的方向参数进行迭代更新，避免了传统方法对投影空间的离散化网格搜索操作。两种合成数据集上的仿真实验结果表明。本文算法具有更高重构精度与鲁棒性。

Abstract

3D reconstruction of single-particle cryo-electron microscopy （cryo-EM） aims to recover the 3D structure of biological molecules from a large number of 2D projection images.In most existing reconstruction methods， the estimation of projection angles and the 3D reconstruction are addressed separately， which may lead to information loss and over-smoothing issues.In this paper， a joint optimization reconstruction algorithm based on the 3D weighted nuclear norm regularization is proposed.In the algorithm， the density map is updated by using the alternating-direction method of multipliers， while the orientations are updated through a block-coordinate gradient descent， thereby avoiding the discretized grid search operation in traditional methods.Two simuation examples are used to verify the accuracy and robustness of the proposed algorithm.

Graphical abstract

关键词

单颗粒冷冻电镜 / 加权核范数 / 联合优化 / 三维重构

Key words

single-particle cryo-EM / weighted nuclear norm minimization / joint optimization / three-dimensional reconstruction

蛋白质等生物大分子参与几乎所有生命活动，其结构研究对于理解其功能至关重要。作为一种革命性结构生物学研究方法，冷冻电镜^［1］（cryo-Electron Microscopy， cryo-EM）能够解析近原子分辨率的三维生物大分子的结构。

冷冻电镜的主要操作流程如下。

1） 样品制备：样品溶液中的生物大分子颗粒以随机的运动方向分布在溶液中，经过快速冷冻使其保持天然构象。

2） 显微成像：将冷冻后的样品在电子显微镜下进行显微成像，此时溶液中随机分布的生物大分子以不同运动方向被捕获，形成包含数十万乃至数百万个二维粒子投影的原始数据集。

3） 图像预处理：由于原始数据集易受噪声和点扩散效应（即成像系统对点源的模糊效应）的影响，进行重构前需要对图像进行分类、对齐及再平均等操作，提升信噪比。

4） 三维图像重构：通过各种图像处理方法估算投影角度，运用直接傅里叶反演方法、反投影算法及迭代算法等方法重建分子的三维结构。

单颗粒分析（Single-Particle Analysis， SPA）^［2］是用冷冻电镜重建生物大分子三维结构的核心技术之一。通过采集大量单颗粒二维投影图像，SPA结合统计学方法来估算单颗粒的三维取向和重构结构，无需结晶即可得到生物大分子的结构。在SPA重建任务中，投影方向是未知的，易受到噪声与点扩散效应的影响。

目前，在主流的SPA重建方法中，变分模型得到广泛应用。总变差（Total Variance，TV）模型和低秩模型是其中两种常见的方法。TV模型通过最小化图像的梯度总和来保留图像的边缘信息^［3］。鉴于高维信号一般具有低秩（矩阵）结构，该方法可以有效地从不完全观测数据中恢复其原始结构。在低秩特征模型中，核范数最小化^［4］（Nuclear Norm Minimization，NNM）方法的应用较为广泛。NNM方法的一个显著优势在于，在一定的数据保真度下，方法是非凸低秩矩阵分解的最紧凸松弛^［5］。文献［6］已经证明，大多数低秩矩阵都可以通过考虑NNM问题得到完全恢复。此外，文献［7］还证明，基于NNM方法的低秩矩阵逼近问题可以通过对观测矩阵的奇异值进行软阈值运算来进行求解。但是，NNM方法存在一个关键缺陷：该方法平等地对待每一个奇异值，在软阈值运算过程中对每个奇异值以相同的阈值进行削减，可能忽略奇异值的先验信息。

一般而言，较大的奇异值对应图像的主要特征，如结构、纹理、边沿等，较小的奇异值则对应噪声等冗余信息。在进行奇异值削减时，合理做法应该是有倾向地削减噪声部分对应的奇异值，同时尽可能保留图像主要特征对应的奇异值。加权核范数最小化（Weighted Nuclear Norm Minimization，WNNM）方法^［8］正是这样一种方法。该方法根据奇异值所表征的物理特性（如信号能量、结构显著性及噪声强度，等）为不同的奇异值分配不同的权重，以增强关键成分、抑制噪声或冗余信息。Li等^［9］提出了一种基于WNNM的非凸模型，用于解决相位恢复问题。Zhao等^［10］提出了一种基于WNNM的自适应加权策略，用于恢复低秩矩阵和张量完备化。Huang等^［11］提出了一种基于WNNM的冷冻电镜三维重建方法。以上研究都证明WNNM方法是鲁棒的。

目前，对于单颗粒冷冻电镜三维图像重构问题，研究者已经提出一些稳健的迭代算法。这些方法均将信号的先验信息纳入考虑，并将三维重构问题表述为正则化逆问题的迭代求解问题^［12］。Barnett等^［13］给出了一种联合重建算法，其中的密度图参数通过梯度下降法进行迭代更新，投影方向通过对三维旋转群SO(\mathrm{3})空间进行穷举搜索进行估算。Zhang等^［14］提出了一种基于矩量子空间的联合重建方法。Zehni等^［15］提出了一种基于TV正则方法的联合优化算法，用于重建分子的三维结构与估算投影方向。值得注意的是，TV正则方法虽然对边缘具有良好的保持能力，但对低对比度区域纹理信息的恢复保持较差，可能造成边缘伪影。此外，TV正则方法还对噪声点较为敏感，容易将噪声误认为边界，从而导致阶梯效应，造成重建质量下降。另一方面，基于核范数加权思想，WNNM方法可以对图像中的异构结构赋予不同的重要性权重，从而在低秩近似方面表现出优越性能。此外，基于WNNM方法的重构方法比基于TV方法的重构方法可能更加稳健。

为了实现有效提取三维结构的细节（特别是边缘和结构特征）并解决现有方法中去噪与结构保持之间的权衡难题，本文提出了一种基于WNNM方法的联合优化算法，以联合重建分子的三维结构并估算每张投影图像对应的投影方向。本文用WNNM方法进行联合优化，并将WNNM模型扩展到三维情形，进而提出单颗粒冷冻电镜的三维重构联合优化方法，恢复三维密度图和相应的投影方向。其中，联合优化问题的解决方法主要通过交替方向乘子法（Alternating Direction Method of Multiplier，ADMM）和块坐标下降法，交替地更新分子密度图和投影方向。

1 模型

假设分子密度图为\mathit{c}\in {\mathbb{R}}^{N\times N\times N}。对于P张投影图像\mathit{g}=\{{\mathit{g}}_p\in {\mathbb{R}}^M{\}}_{p=\mathrm{1}}^P，每张投影图像所对应的投影角度为\mathbf{\Theta }=\{{\mathit{\theta }}_p{\in \mathrm{\Theta }\}}_{p=\mathrm{1}}^P和平移量\mathbf{\Gamma }=\{{\mathit{t}}_p\in {\mathbb{R}}^{\mathrm{2}}{\}}_{p=\mathrm{1}}^P。那么，单颗粒冷冻电镜的全局成像模型^［16］为：

\mathit{g}=\mathit{H}(\mathbf{\Theta },\mathbf{\Gamma })\mathit{c}+\mathit{\epsilon }（1）

其中，\mathit{g}={\left[\begin{array}{c}{\mathit{g}}_{\mathrm{1}},{\mathit{g}}_{\mathrm{2}},\cdots ,{\mathit{g}}_P\end{array}\right]}^{\mathrm{T}}，全局成像算子\begin{array}{c}\mathit{H}(\mathbf{\Theta },\mathbf{\Gamma })={\left[\mathit{H}({\mathit{\theta }}_{\mathrm{1}},{\mathit{t}}_{\mathrm{1}}),\mathit{H}({\mathit{\theta }}_{\mathrm{2}},{\mathit{t}}_{\mathrm{2}}),\cdots ,\mathit{H}({\mathit{\theta }}_P,{\mathit{t}}_P)\right]}^{\mathrm{T}},\end{array}\mathit{H}({\mathit{\theta }}_i,{\mathit{t}}_i)为给定投影角度{\mathit{\theta }}_i和平移量{\mathit{t}}_i的非线性成像算子，\mathit{\epsilon }={\left[{\epsilon }_{\mathrm{1}},{\epsilon }_{\mathrm{2}},\cdots ,{\epsilon }_P\right]}^{\mathrm{T}}为噪声项。

在全局成像模型（1）中，核心的未知量主要包含两类，一类是几何参数(\mathbf{\Theta },\mathbf{\Gamma })，另一类是待恢复的三维密度图参数\mathit{c}。与传统反问题（仅从已知的投影方向信息重建三维结构）相比，此类问题的求解面临以下挑战。

1） 参数维度爆炸：从基于已知几何参数(\mathbf{\Theta },\mathbf{\Gamma })去估计密度图参数\mathit{c}的传统方式扩展至同时求解未知几何参数(\mathbf{\Theta },\mathbf{\Gamma })和密度图参数\mathit{c}，解空间的维度随投影图片数量的增加而线性增长，进而引发严重的维度灾难。

2） 非线性强耦合：在成像算子\mathit{H}(\mathbf{\Theta },\mathbf{\Gamma })中，几何参数与三维结构通过投影变换深度耦合，形成非凸非线性的约束关系，解耦求解困难。

3） 数据稀疏性与噪声问题：投影图像通常具有较高的噪声水平，且每个投影的几何视角与取样是稀疏的，给三维重建带来巨大挑战，往往需要借助正则化技术来抑制噪声并有效提取有用信息、保留全局结构。

2 联合角度细化与重构算法

本文提出一种基于加权核范数正则化的联合优化算法，用于交替估计投影角度和恢复三维结构。

2.1 正则化模型

在冷冻电镜的全局成像模型中，解反问题通常面临不适定性问题。这意味着在给定的观测数据下解可能不是唯一的，或者解对噪声非常敏感。

解决该问题的一种有效途径是正则化方法。具体而言，加权核范数正则化方法通过对解的复杂度进行控制来增强低秩结构的建模能力，避免过拟合并适应噪声和数据的稀疏性，提高解的稳定性和准确性。

本文将问题表述为以下的三维加权核范数正则最小二乘模型：

(\widehat{\mathit{c}},\widehat{\mathbf{\Theta }},\widehat{\mathbf{\Gamma }})\in \underset{\mathit{c},\mathbf{\Theta },\mathbf{\Gamma }}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}\left\{\mathrm{𝒥}(\mathit{c},\mathbf{\Theta },\mathbf{\Gamma })\right\}（2）

\mathrm{𝒥}(\mathit{c},\mathbf{\Theta },\mathbf{\Gamma })=\frac{\mathrm{1}}{\mathrm{2}}{‖\mathit{g}-\mathit{H}(\mathbf{\Theta },\mathbf{\Gamma })\mathit{c}‖}^{\mathrm{2}}+\lambda {‖\mathit{c}‖}_{\mathit{w},\mathrm{*}}（3）

其中{‖\cdot ‖}_{\mathit{w},\mathrm{*}}={\sum }_t{w}_t{\sigma }_t\left({\mathit{c}}_i^j\right)是加权核范数，\mathit{w}是权重向量，{\sigma }_t\left({\mathit{c}}_i^j\right)是\mathit{c}{}_i{}^j的第t个奇异值，{\mathit{c}}_i\in {\mathbb{R}}^{N\times N}是第i个中心切片，\mathit{c}{}_i{}^j是切片{\mathit{c}}_i按某种划分方式得到的第j个子区域，

w_t=\frac{\alpha \sqrt[]{n}}{{\sigma }_i({\mathit{c}}_i^j)+{\epsilon }^{\text{'}}}

是对应块\mathit{c}{}_i{}^j的权重分量，\alpha \mathrm{0}是常数，n是密度图\mathit{c}中相似块的数量，{\epsilon }^{\text{'}}是一个很小的常数，以避免分母为0。

针对以上优化问题，本文使用ADMM方法将全局问题分割成一系列局部子问题。算法1给出了联合优化流程。

算法1 基于WNNM的联合优化算法

输入： 密度图扩展系数的初始值{\mathit{c}}^{\mathrm{0}}，投影角度初始值{\mathbf{\Theta }}^{\mathrm{0}}，平移量初始值{\mathbf{\Gamma }}^{\mathrm{0}}；直到收敛：

1） 更新密度图系数\mathit{c}，

{\mathit{c}}^{k+\mathrm{1}}=\mathrm{A}\mathrm{D}\mathrm{M}\mathrm{M}\mathrm{W}\mathrm{N}(\mathrm{𝒥}(\cdot ,{\mathbf{\Theta }}^k,{\mathbf{\Gamma }}^k),{\mathit{c}}^k)，

2） 更新潜在变量\mathbf{\Theta },\mathbf{\Gamma }，

({\mathbf{\Theta }}^{k+\mathrm{1}},{\mathbf{\Gamma }}^{k+\mathrm{1}})=\mathrm{G}\mathrm{D}(\mathrm{𝒥}({\mathit{c}}^{k+\mathrm{1}},\cdot ,\cdot ),{\mathbf{\Theta }}^k,{\mathbf{\Gamma }}^k)，

3） k=k+\mathrm{1}；

输出： 估算的密度图扩展系数{\mathit{c}}^k，投影角度{\mathbf{\Theta }}^k，平移量{\mathbf{\Gamma }}^k。

2.2 密度图系数更新

对于给定的投影角度\mathbf{\Theta }和平移量\mathbf{\Gamma }，求解密度图扩展系数\mathit{c}等价于求解如下重构问题：

\widehat{\mathit{c}}=\mathrm{a}\mathrm{r}\mathrm{g}\underset{\mathit{c}}{\mathrm{m}\mathrm{i}\mathrm{n}}\mathit{𝒥}(\mathit{c},\mathbf{\Theta },\mathbf{\Gamma })（4）

本文引入辅助变量\mathit{u}，将无约束问题转化为如下的有约束问题：

\widehat{\mathit{c}}=\mathrm{a}\mathrm{r}\mathrm{g}\underset{\mathit{c}}{\mathrm{m}\mathrm{i}\mathrm{n}}\frac{\mathrm{1}}{\mathrm{2}}{‖\mathit{g}-\mathbf{H}(\mathbf{\Theta },\mathbf{\Gamma })\mathit{c}‖}_{\mathrm{2}}^{\mathrm{2}}+\lambda {‖\mathit{u}‖}_{\mathit{w},\mathrm{*}},\mathrm{s}.\mathrm{t}.\mathit{u}=\mathit{c}（5）

为求解问题（5），首先写出问题对应的增广拉格朗日函数，

L_{\mathbf{\Theta },\mathbf{\Gamma }}(\mathit{c},\mathit{u},\tilde{\mathit{u}})=\frac{\mathrm{1}}{\mathrm{2}}{‖\mathit{g}-\mathit{H}(\mathbf{\Theta },\mathbf{\Gamma })\mathit{c}‖}_{\mathrm{2}}^{\mathrm{2}}+

\lambda {‖\mathit{u}‖}_{\mathit{w},\mathrm{*}}+{\tilde{\mathit{u}}}^{\mathrm{{\rm T}}}(\mathit{u}-\mathit{c})+\frac{\rho }{\mathrm{2}}{‖\mathit{u}-\mathit{c}‖}_{\mathrm{2}}^{\mathrm{2}}（6）

其中，\tilde{\mathit{u}}对应约束\mathit{u}=\mathit{c}的拉格朗日乘子，\rho为惩罚参数，用于调节对偶更新步骤的步长，控制算法的收敛速度和稳定性。相应于问题（6）的局部子问题为：

{\mathit{u}}^{k+\mathrm{1}}=\mathrm{a}\mathrm{r}\mathrm{g}\underset{\mathit{u}}{\mathrm{m}\mathrm{i}\mathrm{n}}{L}_{\mathbf{\Theta },\mathbf{\Gamma }}({\mathit{c}}^k,\mathit{u},{\tilde{\mathit{u}}}^k)（7）

{\mathit{c}}^{k+\mathrm{1}}=\mathrm{a}\mathrm{r}\mathrm{g}\underset{\mathit{c}}{\mathrm{m}\mathrm{i}\mathrm{n}}{L}_{\mathbf{\Theta },\mathbf{\Gamma }}(\mathit{c},{\mathit{u}}^{k+\mathrm{1}},{\tilde{\mathit{u}}}^k)（8）

对于子问题（7），求解辅助变量\mathit{u}时可得到如下优化问题：

\begin{array}{l}{\mathit{u}}^{k+\mathrm{1}}=\mathrm{a}\mathrm{r}\mathrm{g}\underset{\mathit{u}}{\mathrm{m}\mathrm{i}\mathrm{n}}\lambda {‖\mathit{u}‖}_{\mathit{w},\mathrm{*}}+\\ ({\tilde{\mathit{u}}}^k{)}^{\mathrm{{\rm T}}}({\mathit{c}}^{k+\mathrm{1}}-\mathit{u})+\frac{\rho }{\mathrm{2}}{‖{\mathit{c}}^k-\mathit{u}‖}_{\mathrm{2}}^{\mathrm{2}}\end{array}（9）

以上优化问题等价于求解以下问题：

{\mathit{u}}^{k+\mathrm{1}}=\mathrm{a}\mathrm{r}\mathrm{g}\underset{\mathit{u}}{\mathrm{m}\mathrm{i}\mathrm{n}}{‖\mathit{u}‖}_{\mathit{w},\mathrm{*}}+\frac{\rho }{\mathrm{2}\lambda }{‖({\mathit{c}}^k+{\tilde{\mathit{u}}}^k/\rho )-\mathit{u}‖}_{\mathrm{2}}^{\mathrm{2}}（10）

\mathit{u}的求解由以下迭代公式给出^［8］：

{\mathit{u}}^{k+\mathrm{1}}=\mathit{U}{\mathrm{𝒟}}_{\tilde{\mathit{w}}}(\mathrm{\Sigma }){\mathit{V}}^{\mathrm{T}}（11）

其中{\mathit{c}}^k+{\tilde{\mathit{u}}}^k/\rho =\mathit{U}\mathrm{\Sigma }{\mathit{V}}^{\mathrm{T}}是{\mathit{c}}^k+{\tilde{\mathit{u}}}^k/\rho的奇异值分解，{\mathrm{𝒟}}_{\tilde{\mathit{w}}}(\mathrm{\Sigma })为软阈值算子，\tilde{\mathit{w}}=\frac{\mathrm{2}\lambda }{\rho }\mathit{w}。

对于子问题（8），本文考虑以下目标函数：

\begin{array}{l}{L}_{\mathbf{\Theta },\mathbf{\Gamma }}(\mathit{c},{\mathit{u}}^{k+\mathrm{1}},{\tilde{\mathit{u}}}^k)=\\ \frac{\mathrm{1}}{\mathrm{2}}{‖\mathit{g}-\mathit{H}(\mathbf{\Theta },\mathbf{\Gamma })\mathit{c}‖}_{\mathrm{2}}^{\mathrm{2}}+\frac{\rho }{\mathrm{2}}{‖\mathit{c}-{\mathit{u}}^{k+\mathrm{1}}+{\tilde{\mathit{u}}}^k/\rho ‖}_{\mathrm{2}}^{\mathrm{2}}\end{array}（12）

注意到以上目标函数是关于变量\mathit{c}的凸二次函数，最小化该目标函数等价于求解以下问题：

\begin{array}{l}({\mathit{H}}^{\mathrm{T}}\mathit{H}(\mathbf{\Theta },\mathbf{\Gamma })+\rho \mathit{I})\mathit{c}=\\ {\mathit{H}}^{\mathrm{T}}(\mathbf{\Theta },\mathbf{\Gamma })\mathit{g}+\rho ({\mathit{u}}^{k+\mathrm{1}}-{\tilde{\mathit{u}}}^k/\rho )\end{array}（13）

以上问题的解由下面的迭代公式给出：

{\mathit{c}}^{k+\mathrm{1}}=({\mathit{H}}^{\mathrm{T}}\mathit{H}{(\mathbf{\Theta },\mathbf{\Gamma })+\rho \mathit{I})}^{-\mathrm{1}}[{\mathit{H}}^{\mathrm{T}}(\mathbf{\Theta },\mathbf{\Gamma })\mathit{g}+

\rho ({\mathit{u}}^{k+\mathrm{1}}-{\tilde{\mathit{u}}}^k/\rho )]（14）

本文用共轭梯度法求解以上问题，并在计算{\mathit{H}}^{\mathrm{T}}\mathit{H}(\mathbf{\Theta },\mathbf{\Gamma })+\rho \mathit{I}和{\mathit{H}}^{\mathrm{T}}(\mathbf{\Theta },\mathbf{\Gamma })\mathit{g}时使用快速表达式进行求解^［16］。此外，对于拉格朗日乘子\tilde{\mathit{u}}，本文使用对偶上升迭代进行更新，即

{\tilde{\mathit{u}}}^{k+\mathrm{1}}={\tilde{\mathit{u}}}^k+\rho ({\mathit{u}}^{k+\mathrm{1}}-{\mathit{c}}^{k+\mathrm{1}})（15）

最终，交替方向乘子法在3个步骤之间交替进行，详见算法2。使用加权核范数求解时，可以通过奇异值分解，用软阈值法求解辅助变量\mathit{u}，然后利用共轭梯度法和{\mathit{H}}^{\mathrm{T}}\mathit{H}(\mathbf{\Theta },\mathbf{\Gamma })的快速表达式迭代来求解密度图的扩展系数\mathit{c}，最后再对\tilde{\mathit{u}}做梯度上升迭代更新。

算法2 ADMM算法

输入： 投影图像\mathit{g}，投影角度初始值{\mathbf{\Theta }}^{\mathrm{0}}，平移量初始值{\mathbf{\Gamma }}^{\mathrm{0}}，密度图系数{\mathit{c}}^{\mathrm{0}}，参数\lambda、\rho和\alpha，迭代次数K_{\mathrm{A}\mathrm{D}\mathrm{M}\mathrm{M}}；

初始化：{\mathit{u}}^{\mathrm{0}}={\mathit{c}}^{\mathrm{0}},{\tilde{\mathit{u}}}^{\mathrm{0}}={\mathit{u}}^{\mathrm{0}}k=\mathrm{0}，

For k{K}_{\mathrm{A}\mathrm{D}\mathrm{M}\mathrm{M}}

1） {\mathit{u}}^{k+\mathrm{1}}=\mathit{U}{\mathrm{𝒟}}_{\tilde{\mathit{w}}}(\mathrm{\Sigma }){\mathit{V}}^{\mathrm{T}}，

2） \mathit{b}={\mathit{H}}^{\mathrm{T}}(\mathbf{\Theta },\mathbf{\Gamma })\mathit{g}+\rho ({\mathit{u}}^{k+\mathrm{1}}-{\tilde{\mathit{u}}}^k/\rho )，

3） {\mathit{c}}^{k+\mathrm{1}}=({\mathit{H}}^{\mathrm{T}}\mathit{H}{(\mathbf{\Theta },\mathbf{\Gamma })+\rho \mathit{I})}^{-\mathrm{1}}\mathit{b}，

4） {\tilde{\mathit{u}}}^{k+\mathrm{1}}={\tilde{\mathit{u}}}^k+\rho ({\mathit{u}}^{k+\mathrm{1}}-{\mathit{c}}^{k+\mathrm{1}})，

5） k=k+\mathrm{1}；

输出： {\mathit{c}}^{K_{\mathrm{A}\mathrm{D}\mathrm{M}\mathrm{M}}}。

2.3 投影方向估计

投影方向求解即以下的优化问题：

(\widehat{\mathbf{\Theta }},\widehat{\mathbf{\Gamma }})\in \mathrm{a}\mathrm{r}\mathrm{g}\underset{\mathbf{\Theta },\mathbf{\Gamma }}{\mathrm{m}\mathrm{i}\mathrm{n}}\left\{\left.\mathrm{𝒥}(\widehat{\mathit{c}},\mathbf{\Theta },\mathbf{\Gamma })\right\}\right.（16）

\mathrm{𝒥}(\widehat{\mathit{c}},\mathbf{\Theta },\mathbf{\Gamma })=\frac{\mathrm{1}}{\mathrm{2}}{‖\mathit{g}-\mathit{H}(\mathbf{\Theta },\mathrm{\Gamma })\widehat{\mathit{c}}‖}_{\mathrm{2}}^{\mathrm{2}}+\lambda {‖\widehat{\mathit{c}}‖}_{\mathit{w},\mathrm{*}}（17）

由于WNNM正则化是与投影方向无关的量，且

\begin{array}{l}\frac{\mathrm{1}}{\mathrm{2}}{‖\mathit{g}-\mathit{H}(\mathbf{\Theta },\mathbf{\Gamma })\widehat{\mathit{c}}‖}_{\mathrm{2}}^{\mathrm{2}}=\\ \frac{\mathrm{1}}{\mathrm{2}}{\displaystyle \sum _{p=\mathrm{1}}^P}{‖{\mathit{g}}_p-\mathit{H}({\mathit{\theta }}_p,{\mathit{t}}_p)\widehat{\mathit{c}}‖}_{\mathrm{2}}^{\mathrm{2}}\end{array}（18）

本文将优化问题（16）（17）重写为

({\widehat{\mathit{\theta }}}_p,{\widehat{\mathit{t}}}_p)\in \mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}{\mathrm{𝒥}}_p(\mathit{\theta },\mathit{t})（19）

目标函数

{\mathrm{𝒥}}_p(\mathit{\theta },\mathit{t})=\frac{\mathrm{1}}{\mathrm{2}}{‖{\mathit{g}}_p-\mathit{H}(\mathit{\theta },\mathit{t})\widehat{\mathit{c}}‖}_{\mathrm{2}}^{\mathrm{2}}（20）

目标函数{\mathrm{𝒥}}_p(\mathit{\theta },\mathit{t})关于\mathit{\theta }和\mathit{t}是可微的，其梯度的显式表达式由文献［15］中定理给出。算法3给出了使用块坐标梯度下降法估算投影角度与平移量的算法框架。

算法3 GD算法

输入： 投影角度初始值{\mathbf{\Theta }}^{\mathrm{0}}，平移量初始值{\mathbf{\Gamma }}^{\mathrm{0}}，密度图系数\mathit{c}，初始步长{\alpha }_{\mathit{\theta }}^{\mathrm{0}}\mathrm{0}和{\alpha }_{\mathit{t}}^{\mathrm{0}}\mathrm{0}，迭代次数K_{\mathrm{G}\mathrm{D}}，\eta \in (\mathrm{0,1})；

初始化：{\mathit{u}}^{\mathrm{0}}={\mathit{c}}^{\mathrm{0}},{\tilde{\mathit{u}}}^{\mathrm{0}}={\mathit{u}}^{\mathrm{0}}，k=\mathrm{0}，

当k{K}_{\mathrm{G}\mathrm{D}}时，

1） 更新投影角度{\mathit{\theta }}_p，设置初始步长{\alpha }_{\mathit{\theta }}={\alpha }_{\mathit{\theta }}^{\mathrm{0}}，梯度更新，即

{\mathit{\theta }}_p^{k+\mathrm{1}}={\mathit{\theta }}_p^k-{\mathit{\alpha }}_{\mathit{\theta }}{\nabla }_{\mathit{\theta }}{\mathit{𝒥}}_p({\mathit{\theta }}_p^k,{\mathit{t}}_p^k)，

回溯线搜索，即若损失函数未减小，则缩小步长并重新梯度更新；

2） 更新平移量{\mathit{t}}_p，设定初始步长{\alpha }_{\mathit{t}}={\alpha }_{\mathit{t}}^{\mathrm{0}}，梯度更新，即

{\mathit{t}}_p^{k+\mathrm{1}}={\mathit{t}}_p^k-{\alpha }_{\mathit{t}}{\nabla }_{\mathit{t}}{\mathit{𝒥}}_p({\mathit{\theta }}_p^{k+\mathrm{1}},{\mathit{t}}_p^k)，

回溯线搜索，即若损失函数未减小，则缩小步长并重新梯度更新；

3） k=k+\mathrm{1}；

输出：{\mathbf{\Theta }}^{K_{\mathrm{G}\mathrm{D}}}，{\mathbf{\Gamma }}^{K_{\mathrm{G}\mathrm{D}}}。

3 模拟分析

3.1 数据集

本文在2个合成数据集上测试优化算法的有效性。第1个数据集是霍利迪交叉复合体（Holliday Junction Complex，HJC）。该数据集是高度受控的\lambda噬菌体病毒DNA从大肠杆菌宿主染色体切除过程的关键中间体，其结构揭示了高度受控的DNA交易机制。第2个数据集是Human Patched 1（PTCH1）蛋白质。该蛋白质是Hedgehog信号通路的重要组成部分。对每个数据集，本文使用Chimer-X软件对得到的电子密度图及蛋白质数据库中的相关原子模型生成合成数据。各数据集的具体信息参见表1。

基于冷冻电镜成像模型，生成P张投影图像，具体实施步骤如下。

1） 投影方向选取：使用P个等距离分布的采样点对投影方向进行采样，得到\{{\theta }_{\mathrm{1},p},{\theta }_{\mathrm{2},p}{\}}_{p=\mathrm{1}}^P，其中平面内旋转角是在区间[\mathrm{0,2}\mathrm{\pi })内均匀采样P个点。

2） 平移量选取：将投影图像的中心点放在水平或垂直方向上随机移动，最大平移量为m_t个像素。

3） 噪声投影图像：给定投影方向根据冷冻电镜成像模型进行前向投影，得到无噪声投影图像{\mathit{g}}^{\mathrm{c}\mathrm{l}\mathrm{e}\mathrm{a}\mathrm{n}}，给这些投影图像添加均值为0，方差为{\sigma }^{\mathrm{2}}的高斯噪声，得到最终实验所需的噪声投影图片，所有投影图片的平均信噪比（Signal-to-Noise，SNR）的计算公式为：

\mathrm{S}\mathrm{N}{\mathrm{R}}_{\mathrm{d}\mathrm{a}\mathrm{t}\mathrm{a}}=\frac{\mathrm{1}}{P}{\displaystyle \sum _{p=\mathrm{1}}^P}\frac{{‖{\mathit{g}}_p^{clean}‖}_{\mathrm{2}}^{\mathrm{2}}}{{\sigma }^{\mathrm{2}}}（21）

3.2 参数设置

本文进行如下的参数初始化设置。

1） 密度图初始化{\mathit{c}}^{\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{t}}：使用 RELION 软件所生成的\mathrm{3}\mathrm{D}Ab⁃initio密度图作为初始近似密度图。

2） 投影角度初始化{\mathit{\theta }}^{\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{t}}：在真实角度基础上添加扰动，即

{\mathit{\theta }}_p^{\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{t}}={\mathit{\theta }}_p^{\mathrm{t}\mathrm{r}\mathrm{u}\mathrm{e}}+{\mathit{\epsilon }}_{\mathit{\theta },p}（22）

其中，{\mathit{\epsilon }}_{\mathit{\theta },p}~U(-e_{\mathit{\theta }},e_{\mathit{\theta }})，p\in \{\mathrm{1},\cdots ,P\}。

3） 面内平移量初始化{\mathit{t}}^{\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{t}}：设置

{\mathit{t}}_p^{\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{t}}=\mathrm{0}， p\in \{\mathrm{1},\cdots ,P\}。

4） 迭代参数设置：将算法1中的迭代次数设定为10，用于更新求解密度图参数和几何参数的迭代次数分别设定为

K_{\mathrm{A}\mathrm{D}\mathrm{M}\mathrm{M}}=\mathrm{200}，K_{\mathrm{G}\mathrm{D}}=\mathrm{3}。

3.3 评价方法

三维重构算法完成重构后，通常用傅里叶壳相关（Fourier Shell Correlation，FSC）曲线来评价重构效果。给定重构结果{\mathit{V}}_{\mathrm{r}\mathrm{e}\mathrm{c}}和真实结构{\mathit{V}}_{\mathrm{g}\mathrm{t}}，FSC的定义为

\mathrm{F}\mathrm{S}\mathrm{C}(\mathit{r})=\frac{{\sum }_{r_i\in \mathit{r}}{\widehat{\mathit{V}}}_{\mathrm{r}\mathrm{e}\mathrm{c}}(r_i){\widehat{\mathit{V}}}_{\mathrm{g}\mathrm{t}}(r_i{)}^{\mathrm{*}}}{\sqrt[]{{\sum }_{r_i\in \mathit{r}}|{\widehat{\mathit{V}}}_{\mathrm{r}\mathrm{e}\mathrm{c}}(r_i){|}^{\mathrm{2}}{\sum }_{r_i\in \mathit{r}}|{\widehat{\mathit{V}}}_{\mathrm{g}\mathrm{t}}(r_i){|}^{\mathrm{2}}}}（23）

其中，

\mathit{r}=\{(x_i,y_i,z_i):\hspace{1em}|x_i^{\mathrm{2}}+y_i^{\mathrm{2}}+z_i^{\mathrm{2}}-r|\le {\epsilon }_r\}，

{\epsilon }_r\mathrm{0}表示离散傅里叶域中所有位于球壳内的点的集合，该球壳以原点为中心，内径为r-{\epsilon }_r，外径为r+{\epsilon }_r。因此，傅里叶壳相关系数旨在计算傅里叶域中密度图对应球壳之间的相关性。

在仿真实验中，选择\mathrm{F}\mathrm{S}\mathrm{C}(\mathit{r})=\mathrm{0.143}对应的空间频率来评价重建质量。图1和图2分别示出了不同算法对两个数据集的重建结果，其中在每个子图中上排为分子的三维结构，下排为三维结构的中心切片密度图。

为评价算法估算投影方向(\mathit{\theta },\mathit{t})的准确性，本文计算了真实投影方向与估算投影方向的均方误差（Mean Square Error，MSE），并与以下3种算法的估计结果进行比较。

1） 非精确投影方向重构法：直接使用未细化初始角度代入本文算法重建密度图。

2） 精确投影方向重构法：将真实投影方向代入本文算法重构密度图。本文将该方法作为重构标准，以帮助量化优化过程的改进效果。

3） 基于TV正则的联合重构方法（Joint-TV） ^［15］：与该方法的比较能够验证加权核范数正则在抑制噪声、克服TV正则带来的边缘平滑和阶梯效应等方面的效果。

3.4 结果与分析

3.4.1 密度图重构比较

正如预期那样，当使用未细化的\mathrm{3}\mathrm{D}投影方向时，重建是失败的（参见图1c和图2c），表明细化投影方向对于算法的重建效果至关重要。同时，如果已知真实的\mathrm{3}\mathrm{D}投影方向，重建也会成功。

此外，可以看到，基于加权核范数正则与基于TV正则的联合优化方法重构得到的密度图与直接利用真实投影方向再进行重构的三维密度图非常相似，表明本文算法能够恰当地重建三维密度图。

最后可以看到，Joint-WNNM重构得到的FSC曲线更接近根据真实投影方向重构的FSC曲线，说明基于WNNM的联合优化方法的结果比基于TV正则的联合优化方法好，Joint-WNNM的重建结果更接近真实三维结构（见图3）。进一步，WNNM联合优化算法能够揭示分子更多的细节，而TV正则因依赖梯度信息使得出现边缘过度平滑现象和阶梯效应，Joint-WNNM算法在重构边缘细节上能更好弥补TV正则的不足。

3.4.2 投影方向细化比较

表2列出了使用本文算法对两个样本进行估算得到的投影方向与真实投影方向之间的MSE值。可以看到，本文算法的估算性能优于基于TV正则的优化算法。以HJC数据集为例，本文算法的整体平均误差比基于TV正则的优化算法下降约13%，投影角度估计误差平均降低21%。

4 结论

本文将三维加权核范数应用于单颗粒冷冻电镜三维重构任务，设计了一种联合优化算法。所得方法能够更好地抑制噪声的影响，克服TV正则固有的阶梯效应。该算法采取交替步骤来迭代更新密度图与投影方向：采用ADMM求解密度图，同时利用半坐标梯度下降法求解投影方向，以极大减少计算成本。仿真实验表明，本文算法的估计效果明显优于基于TV正则的优化算法，重建质量更好。

参考文献

[1] Yip K M， Fischer N Paknia E， et al.Atomic-resolution protein structure determination by cryo-EM ［J］.Nature， 2020， 587（7832）： 157-161.

[2] Nakane T， Kotecha A， Sente A， et al.Single-particle cryo-EM at atomic resolution ［J］.Nature， 2020， 587（7832）： 152-156.

[3] Dai Q， Fang F， Li J， et al.Edge-guided composition network for image stitching ［J］.Pattern Recogn， 2021， 118： 108019.

[4] Recht B， Fazel M， Parrilo P A.Guaranteed minimum-rank solutions of linear matrix equations via nuclear norm minimization ［J］.SIAM Rev， 2010， 52（3）： 471-501.

[5] Jiang T X， Ng M K， Zhao X L， et al.Framelet representation of tensor nuclear norm for third-order tensor completion ［J］.IEEE Trans Image Process， 2020， 29： 7233-7244.

[6] Candes E， Recht B.Exact matrix completion via convex optimization ［J］.Commun ACM， 2012， 55（6）： 111-119.

[7] Cai J F， Candès E J， Shen Z.A singular value thresholding algorithm for matrix completion ［J］.SIAM J Optim， 2010， 20（4）： 1956-1982.

[8] Gu S， Zhang L， Zuo W， et al.Weighted nuclear norm minimization with application to image denoising ［C］//Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition.Piscataway： IEEE， 2014： 2862-2869.

[9] Li Z， Yan M， Zeng T， et al.Phase retrieval from incomplete data via weighted nuclear norm minimization ［J］.Pattern Recogn， 2022， 125： 108537.

[10] Zhao Q， Lin Y， Wang F， et al.Adaptive weighting function for weighted nuclear norm based matrix/tensor completion ［J］.Int J Mach Learn Cybern， 2024， 15（2）： 697-718.

[11] Huang C， Wu T， Li J， et al.Single-particle reconstruction in cryo-EM based on three-dimensional weighted nuclear norm minimization ［J］.Pattern Recogn， 2023， 143： 109736.

[12] Sorzano C O S， Vargas J， Otón J， et al.A survey of the use of iterative reconstruction algorithms in electron microscopy ［J］.BioMed Res Int， 2017， 2017（1）： 6482567.

[13] Barnett A， Greengard L， Pataki A， et al.Rapid solution of the cryo-EM reconstruction problem by frequency marching ［J］.SIAM J Imaging Sci， 2017， 10（3）： 1170-1195.

[14] Zhang A， Mickelin O， Kileel J， et al.Moment-based metrics for molecules computable from cryo-EM images ［PP/OL］.V1.arXiv （2024-01-26） ［2026-04-01］.

[15] Zehni M， Donati L， Soubies E， et al.Joint angular refinement and reconstruction for single-particle cryo-EM ［J］.IEEE Trans Image Process， 2020， 29： 6151-6163.

[16] Donati L， Nilchian M， Sorzano C O S， et al.Fast multiscale reconstruction for Cryo-EM ［J］.J Struct Biol， 2018， 204（3）： 543-554.

基金资助

科技部重点研发计划项目(2020YFA0714003)

四川大学青年教师科技创新能力提升项目(2024SCOQJTX002)