从实验采集信号属性及提取信息的不同，核磁共振波谱技术可以大致分为傅里叶核磁共振（Fourier NMR）与拉普拉斯核磁共振（Laplace NMR）两大类^[1].Fourier NMR可获取化学位移信息，是解析分子结构的强大工具^[2]；而Laplace NMR可提供扩散系数或弛豫时间的信息，从中人们可分析分子的化学结构、物理环境，得到分子动力学及构象、相互作用等信息^[3,4].例如，扩散排序谱（DOSY，diffusion ordered spectroscopy）或者扩散NMR（Diffusion NMR）^[5⇓⇓-8]即属于Laplace NMR，利用该技术可对混合物进行虚拟分离，从而解析出样品中不同的分子成分^[9]，或利用分子扩散系数分析其物理环境、分子大小、聚合程度等等^[10].而测量横向或纵向弛豫时间的实验也同样属于Laplace NMR，在材料化学领域，研究者常通过对弛豫参数的分析获取物质成分或多孔材料的孔径大小分布信息^[11⇓-13]；在结构生物学领域，弛豫参数的测量可有助于研究者们了解蛋白质等大生物分子的构象、相互作用及动力学信息^[14].

在Fourier NMR中所采集的自由感应衰减（FID，free induction decay）信号是具有振荡与衰减的复指数信号，一个极简化的信号模型如(1)式所示：

式中t是演化时间， $Y_{\text {Fourier }}(t)$是实验采集的FID信号，σ是与弛豫时间相关的信号衰减速率，而Ω代表化学位移对应的振荡频率.从(1)式中可以看出，化学位移对应的频率Ω与演化时间t构成了一个傅里叶变换对.

而在Laplace NMR中所采集的是仅有衰减的无振荡指数信号，一个极简化的信号模型如(2)式所示：

式中b是实验中的演化量，在DOSY实验中b是与梯度场强度有关的变量^[15]，而对于弛豫实验而言b是演化时间变量； $Y_{\text {Laplace }}(b)$是所采集的实验信号；σ是与测量参数相关的信号衰减速率，在DOSY实验中即代表扩散系数D，而在横向弛豫实验中$\sigma=T_{2}^{-1}$.在纵向弛豫实验中，实验信号可经由一个线性变换表示为(2)式的形式.从(2)式中可以看出，Laplace NMR实验中待提取的信息σ与实验中的演化量b（也即实验信号的横轴）构成了一个拉普拉斯变换对.

傅里叶变换为酉变换，时域与频域之间的相互变换均有唯一的结果.而与之不同的是，拉普拉斯变换并不具有良好的反变换形式，即，实验采集信号$Y_{\text {Laplace }}(b)$无法通过解析表达式变换为待分析参数σ的分布.从指数衰减信号反推出衰减参数σ的问题是一个典型的逆问题，一般被称为拉普拉斯反演（Laplace inversion）或者拉普拉斯逆变换（ILT，inverse Laplace transform），在下文中统一使用拉普拉斯反演或以简称ILT代表.拉普拉斯反演问题不具有傅里叶变换的良好性质，从数学角度上来说该问题是不适定的，即可能无解或不存在唯一解，这也便导致不同的处理算法可能得到不一样的结果.为了深刻理解Laplace NMR信号处理中面临的挑战，我们将首先讨论其不适定性.

Laplace NMR数据反演问题存在着较强的不适定性^[16]，即可能存在无穷多个解，或者解对于输入数据的微小变化是高度敏感的.在数学中，一个适定问题应该有唯一的解，并且在输入数据的扰动下解仍然保持稳定，而不适定性问题需要额外的信息或约束条件来确保唯一性和稳定性.

在Laplace NMR数据处理时，有两大类处理方案，一类涉及参数定量问题，而另一类则是参数定性分析.无论是定量或定性分析问题，均面临不适定性的挑战，下面我们将分别举例说明.

尽管扩散过程和弛豫过程具有不同的物理本质，但在核磁共振实验中获取的实验信号具有相似的数学描述，即信号将经由指数衰减趋于平衡.在扩散系数或弛豫时间参数分布已知的情况下，其对应实验中可采集到的衰减信号可以通过一个前向的物理模型进行仿真.上一节(2)式是一个简化的单指数模型，而一个更一般的表达形式为：

(3)式中，衰减参数σ由扩散系数或弛豫时间决定（比如在横向弛豫实验中$\sigma=T_{2}^{-1}$，而在扩散实验中σ代表扩散系数D），$S(\sigma)$代表着衰减参数的分布，而b是实验中的演化量（在弛豫实验中b代表回波间隔时间，在扩散实验中b是与梯度场强的平方成正比的量）.尽管该前向物理模型是确定的，但在实际情况下，当我们采集了$Y(b)$并需要反推衰减参数的分布$S(\sigma)$时，则不具有显式的解析表达式.不仅如此，不同的$S(\sigma)$分布可能对应着很相似的$Y(b)$，这会使得逆问题$G: Y(b) \rightarrow S(\sigma)$的求解变得病态和不适定^[1].

在图1中给出一个例子.图1(a)中展示了三种$S(\sigma)$的分布情况，蓝线代表一种宽分布，比如在一些不均匀的混合分散体系（polydisperse）中常出现的扩散系数宽分布情况；红线代表单指数分布模型，如均匀且单一的分散体系（monodisperse）中的扩散系数分布情况；绿线代表双指数分布模型，如均匀且具有两个不同分子成分的混合物样品所对应的扩散系数分布情况.尽管这三种$S(\sigma)$的分布属性不同，其对应的实验采集数据$Y(b)$却非常相近，在图1(b)中所示的线性坐标下，三条衰减信号基本重合.若将其纵坐标以对数形式显示，图1(c)中的宽分布对应的衰减曲线仍旧与双指数衰减曲线基本重合.尽管图1(c)中，单指数分布与其余两种情况相比存在轻微差异，但这些差异与实际不理想实验条件下的噪声等影响相比仍旧微不足道.

因此，从图1(b)或(c)所示的衰减信号反演而准确获得图1(a)中的衰减参数分布是极具挑战的.换言之，利用实验采集的衰减信号对参数分布进行准确定量分析十分困难，因为实验中的一些小扰动可能会导致参数分布计算结果产生显著变化.因此，在分析参数分布时需要引入一些关于分布的先验知识（prior knowledge），从而将所推算的参数分布结果限定在一定的范围内.

在扩散NMR的数据处理中，很多情况下并不需要得到精准的定量结果，而仅需要将不同成分所对应的信号分离开来，这类任务通常被称为定性分析.上一小节我们提到扩散系数的定量估计可能存在多解，而定性分析是否唯一呢?

在扩散NMR（或DOSY）的数据处理中，实验采集的二维信号有一个维度代表演化时间，另一个维度代表与梯度强度相关的演化量b：

其中，γ是旋磁比，δ是梯度脉冲宽度，g为梯度场强， $\Delta^{\prime}$是校正后的扩散时间^[15].对采集数据沿着演化时间维度进行一维傅里叶变换后，时间维变换为频率维，便得到了一组衰减序列：$\left\{Y_{i}(b)\right\}_{i}$，其中i代表频率维中的第i个谱峰.而DOSY数据定性分析的目标便是把$\left\{Y_{i}(b)\right\}_{i}$依据不同的衰减速率（或衰减趋势）进行分组，从而得到类别数（对应于扩散成分的个数），以及将每一个谱峰归类于一个对应的扩散成分.

在图2(a)中给出一个仿真信号，假设该混合物中的两个成分在频率维（即化学位移维度）无重叠，则在化学位移维度中每个谱峰都对应着一个实验中采集的单指数衰减序列，而不同成分对应不同的指数衰减速率.在图2(c)中展示了所有谱峰处提取出的衰减序列$\left\{Y_{i}(b)\right\}_{i}$，其中不同颜色代表着不同成分对应的衰减序列，很明显地可以看出两个成分对应的衰减趋势有所不同.若实验获取的衰减如图2(c)所示，我们可以很容易地根据衰减趋势的异同将{Y_i(b)}_i分为两组，并对相应的谱峰进行归类，从而反推出图2(a)的DOSY谱结果.

然而，DOSY实验信号处理最大的挑战是在谱峰重叠情况下的准确定性分析.在图2(b)中即展示了这种情况，两个成分（分别具有扩散系数D₁与D₂）在化学位移的维度存在若干重叠的谱峰（如图2(b)中黄色高亮区域）.而这些重叠的谱峰对应的衰减趋势将是一种混合态，如图2(d)中的蓝线所示，将处于黑线所代表的第1类衰减（扩散系数为D₁）与红线所代表的第2类衰减（扩散系数为D₂）之间.而这些重叠信号的加入，将使衰减序列$\left\{Y_{i}(b)\right\}_{i}$的类别划分情况变得模糊，如图2(d)所示，使各谱峰的归属及混合物成分的定性分析也变得困难.

图2中显示的只是仿真信号的结果，而在实际实验条件下，因为噪声、相位扭曲等的影响，衰减序列的分组会变得更加困难.在图3中展示了两个实际的混合物样品所采集的DOSY实验数据.这两个混合物样品都各自含有3个不同的成分，因此理论上其衰减序列$\left\{Y_{i}(b)\right\}_{i}$应该包含3组不同的衰减趋势.然而，从图3中很难进行这样的分组，无论是成分个数，抑或是各谱峰的分组情况，都可能存在着很多种不同的估计结果.

因此，无论是定量分析还是定性分析，Laplace NMR的数据处理均涉及不适定的逆问题求解，这意味着可能存在多个解.不同的算法可以通过不同的约束来限制解的空间，从而得到符合预期的结果.在经典的Laplace NMR信号处理算法中，将采用人为给定的信号模型或正则化约束项，这些条件是依据理想信号的先验信息而设计的，在下一节中我们将对这类经典算法进行介绍.

根据解的构造形式的不同，经典的Laplace NMR信号处理算法可以分成两大类：一是参数拟合算法，二是网格化重建算法，另外近年来也出现了二者的结合算法，下面我们将分别进行介绍.

参数拟合算法需要事先给定理想的信号模型，不仅包含衰减函数（例如：指数衰减），还需要包含具体的成分数量.比如，若假设分散体系样本是均匀且单一的，则采集到的扩散NMR将具有单指数衰减的形式，如(2)式所示.而如果分散体系不单一，我们需要预先给定明确的信号物理模型，如：

其中P代表预先给定的成分数量， $\sigma_{k}$与$S_{K}$分别代表每个成分对应的衰减因子及其强度.在拟合算法中需要计算确定每个衰减成分的衰减因子以及强度：$\left\{\sigma_{k}, S_{k}\right\}_{k=1}^{P}$，以使$\hat{Y}(b)$与实验采集信号$Y(b)$之间误差足够小，即：

其中$\boldsymbol{b}=\left[b_{1}, b_{2}, \cdots, b_{N}\right]^{T}$是所有的实验演化量b值排成的向量，假设N是实验采集的b值点数. $\boldsymbol{Y}=\left[Y\left(b_{1}\right), Y\left(b_{2}\right), \cdots, Y\left(b_{N}\right)\right]^{T}$是相对应的实验采集序列;$\boldsymbol{C}=\left[\exp \left(-\sigma_{1} \boldsymbol{b}\right), \exp \left(-\sigma_{2} \boldsymbol{b}\right), \cdots, \exp \left(-\sigma_{P} \boldsymbol{b}\right)\right] \in \mathbb{R}^{N \times P}$是由各衰减成分组成的矩阵, $\boldsymbol{S}=\left[S_{1}, S_{2}, \cdots, S_{P}\right]^{T} \in \mathbb{R}^{P \times 1}$则为各衰减成分的强度.

依据处理序列的并行或串行形式，参数拟合算法又分为单变量拟合（Univariate Fitting）与多变量拟合（Multivariate Fitting）.在单变量拟合中^[17⇓-19]，(6)式中的$\boldsymbol{Y} \in \mathbb{R}^{N \times 1}$是一个N点序列，在处理DOSY数据时，每个化学位移维度上的谱峰都对应着一个N点序列，利用单变量拟合方法，我们将对各个化学位移处的衰减序列进行分别处理，例如将第i个化学位移处对应的衰减序列$Y_{i}$代入(6)式以获得$\left\{\sigma_{k}^{(i)}, S_{k}^{(i)}\right\}_{k=1}^{P}$.而由于不同序列对应的衰减因子是独立计算的，结果受噪声影响较大，在DOSY谱的处理中往往存在着同一个成分的不同谱峰处扩散系数不对齐的现象.另外，单变量拟合方法还存在优化算法收敛性的问题.当模型较为简单时，比如P=1时的单指数信号模型，求解(6)式可采用经典的优化算法，比如简单型法（simplex method），也被称为Nelder-Mead算法（如Matlab中的fminsearch函数）^[20]，Levenberg-Marquardt算法（在Matlab中被集成于lsqcurvefit函数中）^[21]，非负最小二乘法（如Matlab中的lsqnonneg函数）^[22].然而，若模型比较复杂，比如P > 1，则该问题的病态性常使得经典优化算法难以收敛.具体而言，在P > 1时(6)式将出现多个局部极值，优化问题变得更加复杂；另外(6)式中多个参数$\sigma_{k}$对应的衰减向量所构成的矩阵（即(6)式中的矩阵C）容易产生列不满秩的情况，比如，若$\sigma_{k} \gg 1$则列向量$\mathrm{e}^{-\sigma_{k} b}$的所有元素将趋于0，而若估计出的不同参数值特别相似时，如$\sigma_{i} \approx \sigma_{j}$，则相对应的列之间也将有较大的相关性，从而构成病态的矩阵C，进而影响对S矩阵的求解.

而多变量拟合方法则对所有序列进行并行处理，代表性算法包括：直接指数曲线解析算法（DECRA，direct exponential curve resolution algorithm）^{[23⇓⇓-26]}，NMR成分分析方法（CORE，component resolved NMR）^[27⇓-29]，成分分析加速算法（SCORE，speedy component resolution）^[30]，多参数曲线解析（MCR，multivariate curve resolution）^[31]，谱图成分优化解析算法（OUTSCORE，optimized unmixing of true spectra for component resolution）^[32]．与单变量拟合方法中依次处理各个衰减$Y_{i}$不同的是，多变量拟合方法将同时处理$\boldsymbol{Y}=\left[Y_{1}, Y_{2}, \cdots, Y_{M}\right] \in \mathbb{R}^{N \times M}$，其中M代表DOSY数据中所有待处理的衰减序列个数，也即化学位移维度的点数．例如，在CORE与SCORE中处理的目标函数为^[30]：

其中$\boldsymbol{S}_{k}=\left[S_{k}^{(1)}, S_{k}^{(2)}, \cdots, S_{k}^{(M)}\right] \in \mathbb{R}^{1 \times M}$即代表着第k个衰减成分对应的谱图.通过将所有或多个衰减序列同时考虑在优化模型中，多变量拟合方法可以对噪声更为鲁棒，从而提高同一成分谱峰的对齐性.(7)式也可以写成矩阵形式：

其中$\tilde{\boldsymbol{S}}=\left[S_{1}^{T}, S_{2}^{T}, \cdots, S_{P}^{T}\right]^{T} \in \mathbb{R}^{P \times M}, \boldsymbol{C}=\left[\exp \left(-\sigma_{1} \boldsymbol{b}\right), \exp \left(-\sigma_{2} \boldsymbol{b}\right), \cdots, \exp \left(-\sigma_{P} \boldsymbol{b}\right)\right] \in \mathbb{R}^{N \times P}$.与(6)式的矩阵化形式不同之处在于在多变量拟合算法中$\tilde{\boldsymbol{S}}$矩阵是代表着M个化学位移处的重建谱，而(6)式的矩阵化形式中S向量仅代表某一个化学位移处的分析结果.

然而，多参数拟合方法在实际应用时会因为预设信号模型(5)式与实际信号的差别而产生错误的参数估计结果.例如，若设定的成分数P不等于实际样品中的成分数，则算法结果会与实际有较大差别，甚至无法收敛.举个具体例子，选择两种标准混合物样品GSP与QGC的DOSY数据（即图3所示的两组数据），采用SCORE并选择不同的P值进行多参数拟合，算法收敛需要迭代的步数及计算时间如表1所示.这两个样品中均包含3种不同的成分，故理想的P值应该设置为3，在P = 3时SCORE计算GSP数据只需要不到100步即可收敛，而在P > 3时收敛过程将显著变慢，并且计算过程中采用非负最小二乘法常出现无法达到计算精度要求的报错提醒.

同时，我们在图4中展示了在P = 5与P = 8时SCORE对GSP数据的计算结果.P = 5时，计算中出现两个具有较大扩散系数的无效成分，如图4(a4)、(a5)所示；而P = 8时出现了3个无效成分，如图4(b1)、(b7)、(b8)所示，并且把实际中的3个成分拆分成了5个成分，出现了过拟合.无效成分的存在将使拟合方程中的系数矩阵更加病态，比如在扩散系数值D特别大的时候，序列$e^{-Db}$趋于全0向量，即(8)式中的系数矩阵C将出现近似全0的行向量.另外，当计算出的不同扩散系数之间差别减小时，对应的系数矩阵不同列之间的相关度将增大，也会加剧矩阵的病态性，使计算结果不稳定.这种病态性也是当P值选择与实际信号不一致时，计算过程易报错的原因.

与参数拟合方法不同，网格化重建算法不需要给定成分个数，该类算法把参数维划分成网格，并在由网格点构成的信号空间中进行实验数据的表征.人们习惯将这样的结果称为“重建”谱图，因为它等同于对处理信号进行了一次包含编解码的重新表征过程.从(3)式中可以看出，待求解的衰减参数分布$S(\sigma)$经由(3)式的拉普拉斯变换得到了实验中获取的衰减序列$Y(b)$，因此若希望由$Y(b)$计算得到$S(\sigma)$，则需要进行拉普拉斯反演.正因如此，这类网格化重建算法也常被称为拉普拉斯反演（或ILT）.

在将(3)式中的$\sigma $轴网格化后，我们将得到类似于(5)式的衰减信号模型，而与拟合算法不同的是，在该网格化模型中， $\sigma _{k}$是给定的网格点而非需要估计的参数.网格化重建方法所处理的优化目标可以表示为：

其中$\boldsymbol{K}=\left[\exp \left(-\sigma_{1} \boldsymbol{b}\right), \exp \left(-\sigma_{2} \boldsymbol{b}\right), \cdots, \exp \left(-\sigma_{P} \boldsymbol{b}\right)\right] \in \mathbb{R}^{N \times P} $被称为拉普拉斯变换核，代表着将拉普拉斯算子离散化形成的矩阵，而其中$\mathrm{e}^{-\sigma_{k} b} $代表第k个网格点上对应的衰减序列，P代表总的网格点数.待重建的谱图矩阵$\boldsymbol{S} \in \mathbb{R}^{P \times M} $即M个化学位移点上对应的衰减序列分解成P个单指数衰减成分后的各成分强度.值得注意的是，拉普拉斯反演算法是不适定的，求解过程中需要使用一定的先验信息去约束所期望得到的结果的属性，以规范化解空间.这便是(9)式中加入的正则化项$R{S}$的作用，其中λ表示正则化参数，用于调整保真项与正则化项之间的平衡性.

不同的正则化约束会带来不一样的谱图结果，常见的约束项有非负性约束、稀疏约束、光滑性约束等.非负性约束可保证Laplace重建谱图具有明确的物理意义，即(8)式中的S矩阵所有元素均需要大于等于0，或(5)式中每个衰减成分的强度值均需大于等于0.稀疏性约束对应于S的l₁范数正则化项（即$\boldsymbol{R}(\boldsymbol{S})=\|\boldsymbol{S}\|_{1} $），所谓的稀疏，即指非零元素少，稀疏性约束将促使谱图更加简洁，谱峰变窄（即谱图分辨率提高），杂峰减少.而光滑性约束可采用与S相关的l₂范数正则化项实现（如：$\boldsymbol{R}(\boldsymbol{S})=\|\mathrm{L} \boldsymbol{S}\|_{2}^{2} $），也可称为Tikhnov正则化，它强调谱图相邻两点间的连续性，但根据L算子的不同取法重建谱图将有不一样的效果.而在求解这类约束优化问题时，常采用约束正则化（CONTIN，constrained regularization）算法^[18,33]、快速迭代阈值算法（FISTA，fast iterative thresholding algorithms）^[34,35]、交替方向乘子法（ADMM，alternative direction multiplier method）算法^[36⇓-38]、截断牛顿内点法（TNIPM，truncated Newton interior point method）^[39,40]，等等.

具体而言，英国曼彻斯特大学的Morris等人在其重建算法中采用谱图的l₂范数正则化项^[18]，法国德特拉斯堡大学的Delsuc等用最大熵约束法^[41]，英国萨塞克斯大学的Day等人用Tikhnov正则化方法平滑重建结果和增强算法稳定性^[42]，这几种方法都可归为基于l₂范数的光滑性约束；波兰华沙大学的Urbańczyk等人用稀疏性（即l₁范数）与非负性约束保证重建结果的有效性，该方法被称为多指数衰减的迭代阈值算法（ITAMed，iterative thresholding algorithm for multiexponential decay）^[35].本课题组在ITAMed算法基础上改进了优化算法，采用TNIPM算法获得了更高效且效果更好的重建谱图^[40].

上述这类方法均属于单变量算法，即在处理类似DOSY的多衰减序列时，各衰减序列的ILT重建过程实质上是相互独立的.因此，这类算法重建的谱图对齐性较差，也就意味着同一分子的不同谱峰可能显示多个扩散系数值，而理想情况下同一分子的所有谱峰应该显示同一扩散系数值.因此，重建谱的对齐性不佳将直接影响混合物分子的组分分析.

近些年来，出现了一些多变量的网格化重建方法，相比于单变量方法而言，多变量的方法中各衰减序列的ILT重建并不是独立的，需要联合优化以实现完整的谱图重建，而一般这一过程也是由特殊的正则化项所控制的.比如，中国科学院精密测量科学与技术创新研究院的刘买利和张许课题组用低秩性约束使重建结果具有规整的对齐性，该方法被称为Laplace NMR的联合反演算法（SILT，simultaneous inversion of Laplace transform）^[43].在这一方法中，所采用的正则化项为谱图矩阵S的核范数：$\boldsymbol{R}(\boldsymbol{S})=\|\boldsymbol{S}\|_{*} $，借此对S的秩进行约束.当重建谱图的谱峰具有良好对齐性时，其对应的数据矩阵S具有各行、各列间的高度相关性，而该相关性必然导致矩阵低秩，因此约束S的低秩性能够促使重建谱图谱峰对齐.

针对于SILT的重建结果在扩散维分辨率低的问题，本课题组开发了同时约束低秩性与稀疏性的DOSY谱重建算法，得到了兼具高分辨率和规整对齐性的谱图.该算法被称为结合低秩和稀疏的Laplace反演（LRSpILT，low-rank and sparse inverse Laplace transform）^[44]，其目标函数中加入了两个正则化项：

其中$\|\boldsymbol{S}\|_{*, \mathrm{~A}} $代表着加权核范数，用以约束重建谱图的对齐性；$\|\boldsymbol{S}\|_{1, \mathrm{~A}} $是加权l₁范数，用以约束S的稀疏性.如图5所示，单指数拟合算法（图5(a)）得到的谱图中同一个成分的不同谱峰可能出现无法对齐的现象，SILT（图5(b)）重建谱图对齐性较好，但扩散系数维度分辨率低，LRSpILT则在SILT的基础上提高了谱图分辨率.

上文中提到，网格化重建方法获取的是参数的分布，但可能出现谱峰过宽、分辨率低的问题，这是因为网格化重建对成分个数没有限制，导致谱图无法进一步简化.而拟合方法可直接获取参数值，但需要指定信号模型，容易受实测数据中的噪声等不理想因素影响，使重建结果出现相同成分的谱峰不对齐、或伪峰太多等问题，这本质上是由于拟合方法的目标函数中对重建谱图属性的约束不够.

针对于以上两大类方法的问题，本课题组于2022年提出一种混合式的算法，称为协同多指数拟合（CoMeF，Coordinated Multi-exponential Fitting）^[45].在这一方法中，我们在多参数拟合模型中加入了正则化项，以限制拟合模型的复杂度，避免数据过拟合，如(11)式所示：

如拟合算法一样需要给定成分数P，但与(7)式不同的是，在(11)式中多了一个正则化项，为了约束重建谱图的稀疏性以减少伪峰的出现，我们采用的是加权l₁范数.而这种带正则化项的多参数拟合模型比单纯只含保真项的目标函数复杂，且具有非凸特性，若采用前文提及的经典优化算法均无法得到收敛且有效的结果，因此我们设计了一种基于神经网络的优化算法^[46]对其进行求解，最终得到了相比于传统算法而言更优的拟合结果.

如图6所示，传统重建方法[图6(a)]只能提供一张重建谱，若需要扩散系数的确切值则要通过进一步的谱峰提取；传统的多参数拟合方法[图6(b)]可以得到具体的扩散系数值，但会因为强信号的拖尾而产生扩散系统的错误估计；而我们提出的方法[图6(c)]可以直接提供更为准确的扩散系数估计结果.并且值得一提的是，在(11)式中加入的正则化项具有改善问题病态性的作用.当设定的成分数P大于真实数据中的成分个数时，经典算法如SCORE往往很难收敛，如表1所示，这便是因为参数过多导致问题不适定性，比如目标函数的梯度变化过大，或出现多个极值点，等等.加入正则化改善后，算法可更快收敛，并且重建谱图更加简洁，伪峰较少.

该算法中虽然引入了神经网络，但也仅仅将其作为优化器，算法的关键在于设计目标函数，其中需要加入人为定义的正则化项.因此，尽管使用了神经网络，该算法仍应归类于经典算法.

第2部分我们介绍了Laplace NMR数据处理的经典方法，从中可以看出整体上经典方法有两大问题：首先，预设信息是很重要的，在参数拟合算法中需要用户给出预设的信号模型（包括成分个数），而在网格化重建方法中需要设定先验信息，并需要根据不同的数据调整正则化参数.在这个过程中将有大量的人工干预，对于相同的数据，不同用户可能会给出不同的重建结果.其次，经典方法十分依赖优化算法的设计，涉及大量数学公式推导以及收敛条件的摸索，算法的适用性与鲁棒性均受到复杂数据条件下的挑战.

本文中所提经典算法大多数均有开源工具可供调用，比如由英国曼彻斯特大学开发的开源工具包——GNAT（General NMR Analysis Toolbox）^[47]和DOSYToolbox^[48,49]——集成了很多经典Laplace NMR数据处理算法.而本课题组的代码也是开源的，可从相应文章的附加链接中下载，或联系作者获取.

近年来，机器学习与人工智能技术在信号处理领域有了不少成功的应用，也为Laplace反演问题提供了全新的解决思路.在磁共振领域，研究者们基于深度学习开发出很多比传统方法而言更高效准确的信号处理方法.比如，厦门大学陈忠和屈小波课题组应用深度学习实现了非均匀采样核磁共振谱的高效重建^[50]，重建二维谱的时间是经典低秩算法重建时间的4%~8%，而重建三维谱的时间是经典压缩感知算法的12%~22%.中国科学院深圳先进技术研究院梁栋课题组使用自编码机对磁共振图像的先验信息进行学习，并结合传统的图像重建方法实现了磁共振图像欠采样数据的高保真度重建^[51].尽管深度学习在磁共振领域已经得到广泛应用，但在Laplace NMR的数据处理问题上基于深度学习的方法相对较少.需要注意的是，Laplace NMR与一般的基于傅里叶变换的NMR技术有着显著不同.首先，实验测量的数据并不要求等间隔采样，不同实验参数条件下的数据条件（如梯度序列或回波间隔时间等）较难以统一.此外，衰减速率的计算与频率参数的计算具有本质的不同，衰减参数计算的分辨率与实际采样条件的联系并无法用类似傅里叶变换的理论进行分析.因此，能较好处理传统NMR谱图的深度学习网络框架并不一定能适用于Laplace NMR.

目前已有一些工作将数据驱动深度学习方法成功应用于Laplace NMR领域^[52⇓-54].与经典方法不同，深度学习方法无需事先设定数据对应的信号模型，也无需人为给定先验信息，而是通过神经网络自动从训练数据集中学习输入到输出之间的映射关系，以及理想谱图的先验属性.Laplace NMR谱图重建是典型的逆问题，其前向的物理模型是已知的，即，若已知弛豫时间或扩散系数的分布情况，则可以通过物理模型仿真得到实验中可能采集到的Laplace NMR数据.在这种情况下，我们可以通过前向物理模型构建训练数据集，并训练网络完成从Laplace NMR数据到弛豫时间或扩散系数分布情况的“逆向推理”.这一方案也是深度学习用于很多科学领域逆问题求解的一种常用范例，该过程可以由示意图7代表.

在下文中我们挑选其中两种方法作为代表进行介绍.在2022年，中国石油大学Luo等人采用深度学习方法实现了横向弛豫实验数据的分析^[53]，在该工作中尝试了用全连接层搭建起的自编码机、和加入注意力机制的多尺度卷积神经网络（ATT-CNN，attention multi-scale convolutional neural network）等几种网络结构以实现T₂分布的计算.其中在ATT-CNN网络结构中，作者采用了一个全连接层构成的去噪模块对输入的含噪序列进行去噪；接着用一个注意力层提取序列中的能量变化过程，并依此对序列加权，以使网络更关注有效的信号特征；随之使用卷积层提取出不同尺度的特征，并最终通过全连接层转换为T₂的分布曲线.该方法在低信噪比的条件下能取得比经典方法更鲁棒的效果（图8）.其中，BRD代表Butler-Reeds-Dawsons方法^[55]，是类似于CONTIN的加了Tikhnov正则化约束的经典优化算法.从中可以看到，即使是在低信噪比条件下，深度学习方法仍可以得到较接近于理想分布的结果，相比于经典算法的优势很明显.

而在2023年，本课题组也提出一种深度学习方法，称为数据驱动指数曲线翻译器（DRECT，data-driven exponential curve translator）以处理Laplace NMR问题^[54]，既包括弛豫参数的计算，又包括DOSY谱重建.在处理这两类问题时，网络结构和训练集数据有所区别.在处理DOSY谱重建问题时，我们采用了简化的Transformer^[56]结构，该网络在自然语言处理与计算视觉等领域均有十分广泛的应用.而应用于Laplace NMR数据处理问题时，我们仅保留了前端的全连接层，以及由多头注意力机制（multi-head attention）与前向传播（feed forward）组成的编码模块.在图9(a)中展示了网络的结构，输入序列$\left\{Y_{i}(b)\right\}_{i} $经过全连接层将被转换到特征域，随后编码模块将对这个特征数据进行类似于去噪的处理.经过了若干个编码模块后，我们便可得到每个序列$Y_{i}(b) $对应的扩散系数分布，所有分布序列组成的矩阵即构成DOSY重建谱图.值得一提

的是，在Transformer中最引人注目的是其自注意力单元，这种自注意力机制在DOSY谱重建问题中起到了非常关键的作用.自注意力单元将对$Y_{i}(b) $对应的特征序列之间的相关性进行计算，并使网络更关注于其中相关的部分，而减弱其中不相关的部分.经过这样的“重编码”后，不同谱峰对应的$Y_{i}(b) $对应的特征域表征将会趋于简化，如图9(b)所示，其实现的效果类似于去噪.而在计算弛豫参数分布的时候，输入的仅是单个序列，框架中的自注意力机制不再有效，因此在实现弛豫实验数据处理时，我们所构建的网络将不包含该自注意力单元.

在图10中展示了深度学习方法与经典方法对一个称为M6的混合物的DOSY谱图重建结果，图中的红色箭头标示出重建谱中的伪峰或漏峰.该混合物中包含甲醇（methanol）、乙醇（ethanol）、丁醇（butanol）、苏氨酸（threonine）、赖氨酸（lysine）、以及蔗糖（sucrose）六种成分.从图10中可以看出，相比于网格化重建方法[图10(a)]，深度学习网络结果[图10(c)]具有与参数拟合算法类似的高分辨率[图10(b)]，且计算结果准确率较高.

而在表2中将经典算法与深度学习方法的使用过程进行了对比.以LRSpILT^[44]与CoMeF^[45]为例的经典方法需要经过繁琐的参数调整，针对于不同样品所设置的正则化参数均有很大的不同.尤其值得注意的是，在类似于CoMeF的参数拟合算法中需要给定组分数量，但在样本成分较复杂、不同成分扩散系数较为接近等情况下，往往需要设定一个大于理想组分数的成分数量才可以较好地分离不同组分.例如在表2中CoMeF算法的成分数设定为10时才能清楚分离出6个组分的信号，这无疑为调参过程带来很大的挑战.而深度学习方法则无需这种调参过程，对于不同样品可采用相同的预训练网络进行分析.另外，经典算法的单次运行时间大概在数分钟，而深度学习方法仅需数秒即可解决问题.

除了节省参数设置及单次运行时间之外，深度学习方法相比于经典算法还有个很重要的优势，即在噪声条件下可以更加鲁棒，从而增强Laplace NMR谱图重建的可靠性，如图8结果所示.经典的Laplace NMR谱图重建算法依赖于模型假设与人为设置的正则化项，无法灵活适应于不同类型、不同程度的噪声条件.相较而言，深度学习方法可以通过在训练集中加入各种复杂噪声以模拟实验数据中的干扰，并通过大规模的数据训练去捕捉信号与噪声之间的复杂关系，从而使得训练好的网络能够自适应地处理复杂噪声条件下的数据.

值得注意的是，本章所提两种基于深度学习的Laplace NMR数据重建方法得到的谱图具有很不一样的属性，在ATT-CNN^[53]中生成的结果类似于CONTIN等经典的网格化重建方法，可以展示出T₂的分布曲线，而在DRECT^[54]中生成的结果更类似于参数拟合方法，谱峰较窄，并不能展示出衰减参数的分布情况.这种区别的主要根源在于训练集中数据的生成方式，在DRECT中衰减序列的生成用的是简化的信号模型（(5)式，且采用P ≤ 3），并采用理想扩散系数为中心的窄高斯型谱峰组合而成理想DOSY谱图，故输入-输出对中只包含了衰减参数数值信息，并不包括其分布信息；而在ATT-CNN一文中则采用了相对较宽的高斯型谱峰组合成理想DOSY谱图，并且用(3)式的离散化形式计算对应的衰减序列，因此T₂的分布信息是包含在输入-输出对中的.因此，深度学习算法的功能极大依赖于训练数据集的构建，用户需要根据所关心的问题（如：样品参数的分布范围、或者混合物成分的定性归属）而设计不一样的训练数据集让网络进行有针对性的学习.

本文针对Laplace NMR的谱图处理问题进行了系统性的介绍.从该信号处理问题本身的不适定性开始，讨论了Laplace衰减参数精确定量的困难，以及扩散NMR中成分定性分析的模糊性.接着介绍了经典的两大类方法——参数拟合方法与网格化重建方法，拟合算法可以直接计算得到参数值，但其性能极大程度依赖于事先设定好的信号模型及成分个数；网格化重建方法无需给定成分个数，但需要将理想谱图的先验信息建模成正则化约束项加入目标函数中，不同的正则化项、不同的正则化参数设置将得到不同的重建结果.而为了减少在谱图重建过程中的人为干预，可采用数据驱动深度学习方法来解决Laplace NMR数据处理问题，其性能极大程度取决于训练集的构建，而我们所期待获取的重建谱图先验属性也应该在训练集中予以表现.

目前已出现的深度学习Laplace NMR数据处理方法并不是太多，未来我们期待在这一领域能出现更鲁棒、更智能化的方法，可针对用户的不同需求提供相应的训练集构建方式和网络模型.另外，本文仅讨论了一维Laplace NMR问题（DOSY虽然是二维数据，但处理时仅有一维需要经过Laplace反演），二维甚至更高维度的Laplace NMR^[57⇓-59]数据处理问题将具有更严重的病态性，其理论解释、不确定度分析、高精度方法设计均面临着更大的挑战.需要特别注意的是，对于Laplace NMR等磁共振定量参数分析问题而言，仅关注算法性能并不够，不确定度分析也是至关重要的.目前深度学习的不确定度量化方法已在计算机视觉^[60]以及医学成像^[61]中有所应用，将不确定度分析引入Laplace NMR谱图分析中，不仅可以量化重建谱图的可信度，还能更全面地了解模型的局限性和潜在风险，这对于用户对数据的后续分析及决策具有重大意义.

无