具体实施方式

现在将参照附图来详细描述本发明的各种示例性实施例。应注意到：除非另外具体说明，否则在这些实施例中阐述的部件和步骤的相对布置、数字表达式和数值不限制本发明的范围。

以下对至少一个示例性实施例的描述实际上仅仅是说明性的，决不作为对本发明及其应用或使用的任何限制。

对于相关领域普通技术人员已知的技术、方法和设备可能不作详细讨论，但在适当情况下，所述技术、方法和设备应当被视为说明书的一部分。

在这里示出和讨论的所有例子中，任何具体值应被解释为仅仅是示例性的，而不是作为限制。因此，示例性实施例的其它例子可以具有不同的值。

应注意到：相似的标号和字母在下面的附图中表示类似项，因此，一旦某一项在一个附图中被定义，则在随后的附图中不需要对其进行进一步讨论。

本发明针对传统二维磁共振T_1ρ定量成像T_1ρ准备脉冲后只能采集一幅图像、扫描时间长的问题，提供可进行多层扫描的二维T_1ρ定量成像方案和二维高加速倍数欠采样方案，并基于信号补偿的低秩加稀疏分解模型，从高度欠采的欠采样数据中重建出高质量的T_1ρ参数加权图像，从而获得了更准确的T_1ρ参数图。

具体地，参见图1所示，本发明实施例的磁共振旋转坐标系下自旋晶格弛豫成像方法包括以下步骤：

步骤S1000，对于不同自旋锁定时间，配置为分两次采集目标图像的T_1ρ加权图像数据,获得T_1ρ参数弛豫模型。

图2是二维多层T_1ρ定量成像序列的示意图，其中图2(a)是二维多层T_1ρ定量成像序列示意，以采集两次目标图像为例，图2(b)是第一次采集的T_1ρ准备脉冲示意，图2(c)是第二次采集的T_1ρ准备脉冲示意。在图2的实施例中，分两次采集图像，第一次采集时T_1ρ准备脉冲中最后一个90°脉冲沿-x轴方向施加，而第二次采集时T_1ρ准备脉冲中最后一个90°度脉冲沿x轴方向施加，每次采集后设置恢复时间Trec(如可根据需要设置为1ms、5ms、10ms等)，使得磁化矢量恢复到稳态。

具体地，第一次采集的纵向磁化矢量表示为：

M₁(TSL)＝M₀+(M_inite^-TSL/T1ρ-M₀)e^-Trec/T1 (1)

第二次采集的纵向磁化矢量表示为：

M₂(TSL)＝M₀+(-M_inite^-TSL/T1ρ-M₀)e^-Trec/T1 (2)

其中，M_init是施加T_1ρ准备脉冲时刻前的纵向磁化矢量，M₀是平衡状态的纵向磁化矢量，T₁是纵向驰豫的时间常数，TSL表示自旋锁定时间，T_1ρ是旋转坐标系下的自旋晶格驰豫时间。

将公式(1)减公式(2)后得到：

M(TSL)＝2M_inite^-Trec/T1e^-TSL/T1ρ (3)

其中M(TSL)＝M₁(TSL)-M₂(TSL)，令A＝2M_inite^-Trec/T1，由于Trec是一固定值，因此A为一常数，公式(3)可等价为：

M(TSL)＝Ae^-TSL/T1ρ (4)

由公式(4)可知，通过采集不同自旋锁定(TSL)的T_1ρ加权图像，经过拟合后即可得到定量的T_1ρ参数值。

在本发明实施例中，在进行二维多层T_1ρ定量成像时，引入了一个恢复时间，并且分两次采集，两次采集的T_1ρ准备脉冲不同。

步骤S2000，在图像采集过程中，采用频率编码方向全采，相位编码方向变密度欠采的采样方式，获得二维多层T_1ρ成像数据。

可选地，为了加快成像速度，本发明在图像采集中采用频率编码方向全采，相位编码方向变密度欠采的采样方式。

具体地，参见图3的K_x和K_y方向及K_y和K_t方向的欠采样示意，其中K_x方向全采，而在K_y-K_t平面，与传统的基于稀疏采样的压缩感知理论欠采样方式不同，本发明实施例中K空间的中心部分采用均匀欠采样方式，而K空间中心以外区域采用变密度欠采样方式，且采样密度按离K空间中心的距离远近而降低，例如，距离近时采用较高的采样密度，而距离近时采用较低的采样密度。图3(a)是K_x和K_y方向矩阵大小256×256的欠采样示意，图3(b)是K_y和K_t方向矩阵大小256×5的欠采样示意。

根据压缩感知的理论，只要信号是稀疏的或是压缩的，经过一个非相干的测量，利用优化方法通过求解最小化问题，就可以从高度欠采的数据中精确重建出原始信号。因此，利用步骤S1000的二维多层T_1ρ成像序列和步骤S2000的欠采样方式，分两次采集不同TSL的T_1ρ加权图像数据(其中第一次采集和第二次采集的扫描参数完全一样)在提高扫描效率的同时保证了准确的参数加权图像和参数值。

在本发明实施例中，通过变密度欠采样方式在提高扫描效率的同时也能够保证重建出的图像质量。

步骤S3000，采集低分辨率的图像数据，用于重建K空间中心数据和估计多通道线圈敏感度矩阵。

在该步骤中，采集一个低分辨率的数据，用于后续重建K空间中心数据及估计多通道线圈敏感度矩阵，其中多通道线圈敏感度矩阵的估计方法，可采用现有技术的方案，在此不再赘述。

步骤S4000，对二维多层T_1ρ成像数据进行重建，获得重建的T_1ρ参数加权图像。

在此步骤中，可结合现有技术，对欠采的二维多层T_1ρ成像数据进行重建。

例如，结合现有技术的基于信号补偿的低秩加稀疏重建模型(SCOPE)，参见文献(“signal compensation for low-rank plus sparse decomposition”，Phys Med Biol2018；63(18):185009，Zhu Y,Liu Y,Ying L,Peng X等)，对欠采的二维多层T_1ρ成像数据进行重建包括以下步骤：

步骤S141，首先，利用VCC-GRAPPA方法和低分辨率图像，重建出K空间中心部分数据。

其中VCC-GRAPPA方法可参考现有文献(“Improving GRAPPA reconstructionusing joint nonlinear kernel mapped and phase conjugated virtual coils”，Physic in Medicine and Biology，2019，64，14NT01(10pp)，DOI:10.1088/1361-6560/ab274d)

步骤S142，利用重建出的K空间中心部分数据，估计多通道线圈的敏感度矩阵。

该步骤可参考文献(“sensitivity encoding for fast MRI”，Magn Reson Med1999；42(5):952-962，Pruessmann KP,Weiger M等)。

步骤S143，在步骤S141和步骤S142的基础上，对每层的T_1ρ加权图像分别进行重建，求解模型表示为：

min_{X,L,S}‖S‖₁ s.t.C(X)＝L+S,E(X)＝d,Rank(L)＝1 (5)

其中，‖·‖₁是代表l₁范数；C(·)是一个操作算子，表示对图像进行像素级的信号补偿；是要重建的图像序列，且其表示成大小为体素数×TSL数(T)的矩阵；L是用矩阵形式表示的图像的低秩部分，S表示图像和低秩部分L的残差；E是多通道线圈编码矩阵，其等于欠采傅里叶算子与线圈的敏感度矩阵的乘积；Rank(L)表示矩阵L的秩，d表示欠采的K空间数据。

基于公式(4)的T_1ρ弛豫模型，信号补偿可具体表示为将图像中的每个像素乘以一个补偿系数，补偿系数可由下式得到：

Coef＝exp(TSL_k/T_1ρ)，k＝1，2，…，T (6)

其中Coef表示补偿系数，TSL_k是第k个自旋-锁定时间，T是自旋-锁定时间(TSL)的个数。

在一个实施例中，公式(5)的求解过程包括以下步骤：

步骤S151，对VCC-GRAPPA重建的K空间中心数据，将其经过傅里叶变换转化到图像域，根据公式(4)的T_1ρ弛豫模型对图像进行拟合，估计初始的T_1ρ参数，并根据公式(6)得到补偿系数的初始值Coef⁰；

步骤S152，设定循环为i＝1,2…，在第i次迭代中，执行：

步骤S152-1根据补偿系数，对图像进行补偿，即其中U表示补偿后的图像；

步骤S152-2，初始化S＝0，设定外循环次数为J，在第j＝1，2，…，J次迭代中，执行：

a)、更新其中SVT(·)表示奇异值阈值操作算子，其定义为：

SVT_λ(M)＝UΛ_λ(Σ)V^H (7)

其中M＝UΣV^H表示奇异值分解(SVD)，U、V分别为左、右奇异值向量组成的矩阵，V^H表示V的共轭转置，Σ是由M的奇异值组成的对角矩阵，Λ_λ(Σ)表示保留Σ中最大的奇异值不变，其他全为0，在本发明实施例中，只取L的最大奇异值，使得做奇异值阈值操作后L的秩Rank(L)＝1；

b)、更新S_j：ST(·)是一个软阈值操作算子，定义为：

其中p是图像矩阵的一个元素，v是阈值。

c)、更新数据保真项：其中E^*表示E的逆操作，即等于对多通道线圈K空间数据做傅里叶逆变换后再进行线圈组合，得到图像；

d)、更新图像其中C^-1(·)表示将图像基于每个像素除以补偿系数Coefⁱ；

e)、终止内循环迭代

步骤S152-3，根据步骤S152-2得到的Xⁱ，结合公式(4)中的参数弛豫模型，更新并更新补偿系数

步骤S152-4，当算法达到迭代终止条件时(例如，迭代次数大于最大的迭代次数或相邻两次迭代之间的重建误差小于预设值时)，终止循环迭代，得到最终重建的参数加权图像X。

步骤S5000，根据重建的T_1ρ参数加权图像和T_1ρ参数弛豫模型，对图像中的所有像素进行非线性拟合，拟合出最终的T_1ρ参数图。

在上述实施例中，在重建二维T_1ρ加权图像时，先利用VCC-GRAPPA重建K空间中心部分，去除由于K空间中心欠采样而造成的卷摺伪影，然后基于信号补偿的低秩加稀疏分解模型，迭代重建每一层的T_1ρ加权图像。其中基于信号补偿的低秩加稀疏分解重建过程是：首先基于T_1ρ参数弛豫模型采用信号补偿的方法来增强数据在TSL方向的低秩性。然后，将图像序列按参数方向排成一个空间-参数矩阵，其中矩阵的每一列表示某一TSL时刻采集到的磁共振图像，然后将空间-参数矩阵(这个矩阵实际上代表着图像)分解成低秩部分(low-rankcomponent，L)和稀疏部分(sparse component，S)。接下来，结合(L+S)重建模型，并对L做奇异值阈值操作，对S做软阈值操作，得到迭代更新的L和S，通过对更新后的L和S求和，可得到更新的空间-参数矩阵(即为图像)。在图像重建的迭代过程中，每次迭代都会根据新重建出的T_1ρ加权图像和T_1ρ参数弛豫模型，更新T_1ρ参数图，并将更新后的T_1ρ参数图用于下一次迭代中的信号补偿，如此反复迭代，直到算法达到迭代终止条件，停止重建。最后，利用T_1ρ参数弛豫模型对重建的参数加权图像进行拟合，得到最终的T_1ρ参数图。

综上所述，与现有的二维T_1ρ定量成像技术相比，本发明可实现多层T_1ρ定量成像，且设计的变密度欠采样方式能极大地加快数据扫描速度，减少T_1ρ定量成像时间，在图像重建时，本发明能够精确地从高度欠采的数据中重建出参数加权图像，通过定量计算进一步提高了图像信噪比。

本发明可以是系统、方法和/或计算机程序产品。计算机程序产品可以包括计算机可读存储介质，其上载有用于使处理器实现本发明的各个方面的计算机可读程序指令。

计算机可读存储介质可以是可以保持和存储由指令执行设备使用的指令的有形设备。计算机可读存储介质例如可以是――但不限于――电存储设备、磁存储设备、光存储设备、电磁存储设备、半导体存储设备或者上述的任意合适的组合。计算机可读存储介质的更具体的例子(非穷举的列表)包括：便携式计算机盘、硬盘、随机存取存储器(RAM)、只读存储器(ROM)、可擦式可编程只读存储器(EPROM或闪存)、静态随机存取存储器(SRAM)、便携式压缩盘只读存储器(CD-ROM)、数字多功能盘(DVD)、记忆棒、软盘、机械编码设备、例如其上存储有指令的打孔卡或凹槽内凸起结构、以及上述的任意合适的组合。这里所使用的计算机可读存储介质不被解释为瞬时信号本身，诸如无线电波或者其他自由传播的电磁波、通过波导或其他传输媒介传播的电磁波(例如，通过光纤电缆的光脉冲)、或者通过电线传输的电信号。

这里所描述的计算机可读程序指令可以从计算机可读存储介质下载到各个计算/处理设备，或者通过网络、例如因特网、局域网、广域网和/或无线网下载到外部计算机或外部存储设备。网络可以包括铜传输电缆、光纤传输、无线传输、路由器、防火墙、交换机、网关计算机和/或边缘服务器。每个计算/处理设备中的网络适配卡或者网络接口从网络接收计算机可读程序指令，并转发该计算机可读程序指令，以供存储在各个计算/处理设备中的计算机可读存储介质中。

用于执行本发明操作的计算机程序指令可以是汇编指令、指令集架构(ISA)指令、机器指令、机器相关指令、微代码、固件指令、状态设置数据、或者以一种或多种编程语言的任意组合编写的源代码或目标代码，所述编程语言包括面向对象的编程语言—诸如Smalltalk、C++等，以及常规的过程式编程语言—诸如“C”语言或类似的编程语言。计算机可读程序指令可以完全地在用户计算机上执行、部分地在用户计算机上执行、作为一个独立的软件包执行、部分在用户计算机上部分在远程计算机上执行、或者完全在远程计算机或服务器上执行。在涉及远程计算机的情形中，远程计算机可以通过任意种类的网络—包括局域网(LAN)或广域网(WAN)—连接到用户计算机，或者，可以连接到外部计算机(例如利用因特网服务提供商来通过因特网连接)。在一些实施例中，通过利用计算机可读程序指令的状态信息来个性化定制电子电路，例如可编程逻辑电路、现场可编程门阵列(FPGA)或可编程逻辑阵列(PLA)，该电子电路可以执行计算机可读程序指令，从而实现本发明的各个方面。

这里参照根据本发明实施例的方法、装置(系统)和计算机程序产品的流程图和/或框图描述了本发明的各个方面。应当理解，流程图和/或框图的每个方框以及流程图和/或框图中各方框的组合，都可以由计算机可读程序指令实现。

这些计算机可读程序指令可以提供给通用计算机、专用计算机或其它可编程数据处理装置的处理器，从而生产出一种机器，使得这些指令在通过计算机或其它可编程数据处理装置的处理器执行时，产生了实现流程图和/或框图中的一个或多个方框中规定的功能/动作的装置。也可以把这些计算机可读程序指令存储在计算机可读存储介质中，这些指令使得计算机、可编程数据处理装置和/或其他设备以特定方式工作，从而，存储有指令的计算机可读介质则包括一个制造品，其包括实现流程图和/或框图中的一个或多个方框中规定的功能/动作的各个方面的指令。

也可以把计算机可读程序指令加载到计算机、其它可编程数据处理装置、或其它设备上，使得在计算机、其它可编程数据处理装置或其它设备上执行一系列操作步骤，以产生计算机实现的过程，从而使得在计算机、其它可编程数据处理装置、或其它设备上执行的指令实现流程图和/或框图中的一个或多个方框中规定的功能/动作。

附图中的流程图和框图显示了根据本发明的多个实施例的系统、方法和计算机程序产品的可能实现的体系架构、功能和操作。在这点上，流程图或框图中的每个方框可以代表一个模块、程序段或指令的一部分，所述模块、程序段或指令的一部分包含一个或多个用于实现规定的逻辑功能的可执行指令。在有些作为替换的实现中，方框中所标注的功能也可以以不同于附图中所标注的顺序发生。例如，两个连续的方框实际上可以基本并行地执行，它们有时也可以按相反的顺序执行，这依所涉及的功能而定。也要注意的是，框图和/或流程图中的每个方框、以及框图和/或流程图中的方框的组合，可以用执行规定的功能或动作的专用的基于硬件的系统来实现，或者可以用专用硬件与计算机指令的组合来实现。对于本领域技术人员来说公知的是，通过硬件方式实现、通过软件方式实现以及通过软件和硬件结合的方式实现都是等价的。

以上已经描述了本发明的各实施例，上述说明是示例性的，并非穷尽性的，并且也不限于所披露的各实施例。在不偏离所说明的各实施例的范围和精神的情况下，对于本技术领域的普通技术人员来说许多修改和变更都是显而易见的。本文中所用术语的选择，旨在最好地解释各实施例的原理、实际应用或对市场中的技术改进，或者使本技术领域的其它普通技术人员能理解本文披露的各实施例。本发明的范围由所附权利要求来限定。