CN107102020B - 多维核磁共振测量方法 - Google Patents

Abstract

本发明提出了一种用于材料检测的多维核磁共振数据采集及处理的分析方法。该方法适用于不同场强不同配置下的核磁共振检测仪器。该方法涉及多种核磁共振脉冲序列，通过对得到的三维数据结构体进行相应数据处理，可分析得到检测材料的孔隙度、孔隙结构、饱含流体种类、骨架磁化系数以及非均质性等重要参数。

Description

多维核磁共振测量方法

技术领域

本发明涉及核磁共振领域，尤其涉及一种多维核磁共振测量方法的实现过程及其数据处理的过程。

技术背景

核磁共振技术作为一种先进的无损探测手段，在医学、能源、材料、农林、食品、安全监测、化工等多个领域均有着极为广泛的应用。核磁共振常见的特性参数包括纵向弛豫时间T₁、横向弛豫时间T₂，另外通过核磁共振技术可以快速获得分子的扩散系数D或固体骨架的磁化系数χ。通过在脉冲序列的不同时间段内编辑相应核磁共振特性，可以在多个维度上关联这些核磁共振特性。由这些关联实验可以对孔隙介质性质进行全面的研究，得到更加丰富的信息。核磁共振弛豫特性被应用于探测孔隙介质的孔隙结构，孔隙连通性，孔隙空间结构组成；扩散系数被用于了解内部饱含流体种类以及原油的组分等。因此通过三维核磁共振技术能够在更高的维度上对孔隙介质的孔隙尺寸、内部磁场梯度及复杂流体的分子组成进行综合研究。

孔隙介质内的核磁共振弛豫或扩散响应满足多指数衰减规律。由测量结果反演得到孔隙空间弛豫或扩散分布是一个病态问题。该反演过程的解不唯一，并且测量结果中小的扰动就会对求解结果造成较大影响。通过在求解过程中引入正则化项，对反演过程进行稳定处理。三维核磁共振的数据量大，按照一维数据处理思路进行反演耗时较长。通过对核函数分别压缩后进行张量积，对得到的新的核函数进行奇异值分解，截取满足设定条件数的奇异值，同时采用截取后的正交单位矩阵对三维数据进行压缩，加快反演速度。

发明内容

本发明的目的是阐述三种用于分析材料孔隙结构和内部填充流体特性的多维核磁共振测量方法，以及相应的多维核磁共振数据处理流程。

第一种多维核磁共振测量方法，所述方法包括：

步骤1、向被测样品施加90°射频脉冲将所述宏观磁化强度矢量M₀扳转90°；

步骤2、等待时间τ过后，向被测样品施加第二个90°脉冲，将恢复一定量的磁化矢量从沿着静磁场B₀相同的方向再次扳转90°；

步骤3、在等待T_E/2的时间之后，向被测样品施加180°射频脉冲，相同的等待时间T_E/2过后，会在ACQ通道采集产生自旋回波信号；

步骤4、重复施加180°射频脉冲，会在ACQ通道采集产生重复多个自旋回波信号，得到回波串信号；

步骤5、改变等待时间τ，重复所述步骤2-4以分别采集数个不同等待时间τ下产生的回波串信号；

步骤6、根据采集到的所述回波串信号进行核磁共振数据处理。

其中，步骤5所产生的回波串信号磁化矢量矩阵为：

M(τ，nT_E，mt_s)＝∫∫∫K₁K₂K₃F(T₁，T₂，Δχ)dT₁·dT₂·dΔχ

其中，τ为等待时间，n为180°脉冲个数，m为FID采集点数，t_s为FID采集点时间间隔，F(T₁，T₂，Δx)为被测样品的三维T₁-T₂-Δχ特性矩阵，T₁为纵向弛豫时间，T₂为横向弛豫时间，Δχ为被测样品与内部填充流体的磁化系数差异，K₁，K₂，K₃为三个核函数，其具体形式为：

K₁＝1-exp(-τ/T₁)

K₂＝exp(-nT_E/T₂)

K₃＝exp(-γ·Δχ·B₀·mt_s)

其中，γ为质子的旋磁比，B₀为静磁场强度。

所述步骤6具体为，采用快速三维数据处理算法对多维核磁共振数据进行反演，通过选取合适的正则化因子，得到被测样品的三维T₁一T₂-Δχ分布信息。具体为：

步骤1：采用数学张量乘积将前两个核矩阵函数K₁和K₂耦合为一个新的核函数矩阵K₁₂：

步骤2：将测量得到的三维核磁共振数据重新表达为：

M＝K₁₂FK₃

步骤3：对上述核函数矩阵K₁₂和K₃进行SVD分解及奇异值截取，进而对采集数据进行压缩处理，对核函数矩阵进行奇异值分解可得：

K₁₂＝U₁₂·S₁₂·V′₁₂

K₃＝U₃·S₃·V′₃

其中S₁₂和S₃对角线元素值从大到小排列，且为对角矩阵，大小分别为s₁₂×s₁₂和s₂×s₃，其中s₁₂为K₁₂非零奇异值个数，s₃为K₃非零奇异值个数；U₁₂、V₁₂和U₃、V₃为正交单位阵；对对角矩阵S₁₂和S₃进行截取，使得核函数矩阵的条件数满足设定值C，即：

假设C为1000；和分别对应K₁₂和K₃最大的奇异值，即对角矩阵S₁₂和S₃的第一个对角线元素，表示K₁₂的第i个奇异值，表示K₃的第j个奇异值；

步骤4：用截取之后的奇异值分解的单位矩阵对回波串信号磁化矢量矩阵M进行压缩，降低数据内存，压缩后的磁化矢量为：

其中，为压缩后的磁化矢量矩阵，为两个核函数K₁₂和K₃经过SVD分解和奇异值截取后的残余矩阵，U′₁₂、V′₁₂分别为K₁₂经过SVD分解和奇异值截取后的矩阵U₁₂、V₁₂的转置，V′₃为K₃经过SVD分解和奇异值截取后的矩阵V₃的转置；

步骤5：数据压缩完成之后，采用Tikhonov正则化方法对数据矩阵进行反演，其中，正则化项为：

其中，α是正则化因子，与采集数据的信噪比相关，K为反演过程中的核矩阵的一般形式，||·||项代表矩阵的Frobenius范数；

步骤6：通过选取最优的正则化因子α，得到最终解f_r，求解公式如下：

其中，为压缩截取后残余矩阵的张量积矩阵，为的转置矩阵，α_opt为最优正则化因子，I是单元矩阵，为压缩后的磁化矢量矩阵的一维矩阵形式，其中，

第二种多维核磁共振测量方法，所述方法包括：

步骤1、在TRS通道向被测样品施加90°射频脉冲将所述宏观磁化强度矢量M₀扳转90°；

步骤2、在GRD通道施加一个持续时间为δ的梯度脉冲；

步骤3、随后在TRS通道施加180°射频脉冲后，在GRD通道施加第二个持续时间为δ的梯度脉冲。本步骤中的梯度脉冲起始时间与步骤2中的梯度脉冲起始相差Δ时间；

步骤4、在整个系统等待前两个射频脉冲间隔相等的时间后，向被测样品施加90°射频脉冲，将一定量的磁化矢量从垂直静磁场B₀相同的方向扳转90°；

步骤5、施加小角度射频脉冲，在ACQ通道采集自由衰减信号；

步骤6、等待时间t_a后，重复施加小角度射频脉冲，会在ACQ通道重复采集到自由衰减信号；

步骤7、改变GRD通道中的两个梯度脉冲的强度值，重复所述步骤4-6以分别采集不同脉冲强度下产生的自由衰减信号；

步骤8、根据采集到的所述自由衰减串信号进行核磁共振数据处理。

所述步骤7得到的自由衰减信号磁化矢量矩阵为：

M(g，nt_a，mt_s)＝∫∫∫K₁K₂K₃F(D，T₁，Δχ)dD·dT₁·dΔχ

其中，g为静磁场梯度值，n为180°脉冲个数，t_a为等待时间，m为FID采集点数，t_s为FID采集点时间间隔，F(D，T₁，Δχ)为被测样品的三维D-T₁-Δχ特性矩阵，T₁为纵向弛豫时间，D为被测样品自扩散系数，Δχ为被测样品与内部填充流体的磁化系数差异，K₁，K₂，K₃为三个核函数，其具体形式为：

K₁＝exp(-Dγ²g²δ²(Δ-δ/3))

K₂＝exp(-nt_a/T₁)

K₃＝exp(-γ·Δχ·B₀·mt_s)

通过改变D编辑中的脉冲梯度强度g，T₁编辑中的小角度脉冲个数n，以及Δχ编辑中的采集点数m获取三维数据。采用快速三维数据处理算法对其进行反演，通过选取合适的正则化因子，得到被测样品的三维D-T₁-Δχ分布信息。具体为：

2.1、采用数学张量乘积将前两个核矩阵函数K₁和K₂耦合为一个新的核函数矩阵K₁₂：

2.2、将测量得到的回波串信号磁化矢量矩阵重新表达为：

M＝K₁₂FK₃

2.3、对上述核函数矩阵K₁₂和K₃进行SVD分解及奇异值截取，进而对采集数据进行压缩处理，对核函数矩阵进行奇异值分解可得：

K₁₂＝U₁₂·S₁₂·V′₁₂

K₃＝U₃·S₃·V′₃

其中S₁₂和S₃对角线元素值从大到小排列，且为对角矩阵，大小分别为s₁₂×s₁₂和s₃×s₃，其中s₁₂为K₁₂非零奇异值个数，s₃为K₃非零奇异值个数；U₁₂、V₁₂和U₃、V₃为正交单位阵；对对角矩阵S₁₂和S₃进行截取，使得核函数矩阵的条件数满足设定值C，即：

假设C为1000；和分别对应K₁₂和K₃最大的奇异值，即对角矩阵S₁₂和S₃的第一个对角线元素，表示K₁₂的第i个奇异值，表示K₃的第j个奇异值；

2.4：用截取之后的奇异值分解的单位矩阵对回波串信号磁化矢量矩阵M进行压缩，降低数据内存，压缩后的磁化矢量为：

其中，为压缩后的磁化矢量矩阵，为两个核函数K₁₂和K₃经过SVD分解和奇异值截取后的残余矩阵，U′₁₂、V′₁₂分别为K₁₂经过SVD分解和奇异值截取后的矩阵U₁₂、V₁₂的转置，V′₃为K₃经过SVD分解和奇异值截取后的矩阵V₃的转置；

2.5：数据压缩完成之后，采用Tikhonov正则化方法对数据矩阵进行反演，其中，正则化项为：

其中，α是正则化因子，与采集数据的信噪比相关，K为反演过程中的核矩阵的一般形式，||·||项代表矩阵的Frobenius范数；

2.6：通过选取最优的正则化因子α，得到最终解f_r，求解公式如下：

其中，为压缩截取后残余矩阵的张量积矩阵，为的转置矩阵，α_opt为最优正则化因子，I是单元矩阵，为压缩后的磁化矢量矩阵的一维矩阵形式，其中，

第三种多维核磁共振测量方法，所述方法包括：

步骤1、向被测样品施加90°脉冲将宏观磁化强度矢量M₀扳转90°；

步骤2、等待时间τ过后，系统向被测样品施加第二个90°脉冲，将恢复一定量的磁化矢量从沿着静磁场B₀相同的方向再次扳转90°；

步骤3、在GRD通道施加一个持续时间为δ的梯度脉冲；

步骤4、随后在TRS通道施加180°射频脉冲后，在GRD通道施加第二个持续时间为δ的梯度。本步骤中的梯度脉冲起始时间与步骤2中的梯度脉冲起始相差Δ时间；

步骤5、在整个系统等待T_E/2的时间之后，向被测样品施加180°射频脉冲，相同的等待时间T_E/2过后，会在ACQ通道采集产生自旋回波信号；

步骤6、重复施加180°射频脉冲，会在ACQ通道采集产生重复多个自旋回波信号，称之为回波串信号，记录回波串信号尖峰值；

步骤7、改变GRD通道中的两个梯度脉冲的强度值，重复所述步骤4-6以分别采集不同脉冲强度下产生的回波串信号尖峰值；

步骤8、改变等待时间τ，重复所述步骤7以分别采集不同等待时间下产生的回波串信号尖峰值；

步骤9、根据采集到的所述自由衰减串信号进行核磁共振数据处理。

步骤8得到的回波串信号尖峰值磁化矢量矩阵为：

M(τ，g，nT_E)＝∫∫∫K₁K₂K₃F(T₁，D，T₂)dT₁·dD·dT₂

其中，τ为等待时间，g为静磁场梯度值，n为180°脉冲个数，T_E为回波间隔，F(T₁，D，T₂)为被测样品的三维T₁-D-T₂特性矩阵，T₁为纵向弛豫时间，T₂为横向弛豫时间，D为被测样品自扩散系数，K₁，K₂，K₃为三个核函数，其具体形式为：

其中γ为质子的旋磁比，g为脉冲梯度场梯度强度值，δ为单个脉冲梯度持续时间。通过改变T₁编辑时间τ，D编辑中的脉冲梯度强度g和T₂编辑中的180°脉冲个数n获取三维数据。采用快速三维数据处理算法对其进行反演，通过选取合适的正则化因子，得到被测样品的三维T₁-D-T₂分布信息。具体为：

2.1、采用数学张量乘积将前两个核矩阵函数K₁和K₂耦合为一个新的核函数矩阵K₁₂：

2.2、将测量得到的回波串信号磁化矢量矩阵重新表达为：

M＝K₁₂FK₃

2.3、对上述核函数矩阵K₁₂和K₃进行SVD分解及奇异值截取，进而对采集数据进行压缩处理，对核函数矩阵进行奇异值分解可得：

K₁₂＝U₁₂·S₁₂·V′₁₂

K₃＝U₃·S₃·V′₃

其中S₁₂和S₃对角线元素值从大到小排列，且为对角矩阵，大小分别为s₁₂×s₁₂和s₃×s₃，其中s₁₂为K₁₂非零奇异值个数，s₃为K₃非零奇异值个数；U₁₂、V₁₂和U₃、V₃为正交单位阵；对对角矩阵S₁₂和S₃进行截取，使得核函数矩阵的条件数满足设定值C，即：

假设C为1000；和分别对应K₁₂和K₃最大的奇异值，即对角矩阵S₁₂和S₃的第一个对角线元素，表示K₁₂的第i个奇异值，表示K₃的第j个奇异值；

2.4、用截取之后的奇异值分解的单位矩阵对回波串信号磁化矢量矩阵M进行压缩，降低数据内存，压缩后的磁化矢量为：

其中，为压缩后的磁化矢量矩阵，为两个核函数K₁₂和K₃经过SVD分解和奇异值截取后的残余矩阵，U′₁₂、V′₁₂分别为K₁₂经过SVD分解和奇异值截取后的矩阵U₁₂、V₁₂的转置，V′₃为K₃经过SVD分解和奇异值截取后的矩阵V₃的转置；

2.5、数据压缩完成之后，采用Tikhonov正则化方法对数据矩阵进行反演，其中，正则化项为：

其中，α是正则化因子，与采集数据的信噪比相关，K为反演过程中的核矩阵的一般形式，||·||项代表矩阵的Frobenius范数；

2.6、通过选取最优的正则化因子α，得到最终解f_r，求解公式如下：

其中，为压缩截取后残余矩阵的张量积矩阵，为的转置矩阵，α_opt为最优正则化因子，I是单元矩阵，为压缩后的磁化矢量矩阵的一维矩阵形式，其中，

相较于常规的一维核磁共振测量，本申请的多维核磁共振测量由于能够提供更加丰富的弛豫及扩散分布信息，因而在多个领域中得到了广泛的应用。多维核磁共振技术可用于检测孔隙介质的孔隙结构，识别不同种类的流体及确定赋存流体的饱和度信息。

附图说明

图1为本发明实施例一核磁共振测量方法的脉冲序列图；

图2为本发明实施例二核磁共振测量方法的脉冲序列图；

图3为本发明实施例三核磁共振测量方法的脉冲序列图；

图4为本发明实施例四处理多维核磁共振数据的反演流程图。

其中，TRS代表核磁共振系统的脉冲发射通道，ACQ代表核磁共振系统的信号接受通道，GRD通道代表梯度脉冲发射通道。

具体实施方式

结合说明书附图说明本发明的具体实施方式。在此需要说明的是，对于这些实施例方式的说明用于帮助理解本发明，但并不构成对本发明的限定。

首先对本发明中涉及的相关技术术语的定义及其物理意义做如下介绍。

静磁场B₀。静磁场由磁体提供，决定核磁共振信号的信噪比。被测样品置于静磁场中，自旋系统内发生能级分裂，沿着静磁场方向会产生一个宏观磁化矢量M₀。M₀由静磁场强度B₀，温度等参数决定。

射频磁场B₁与射频脉冲。射频脉冲为电磁信号，通常由线圈产生。射频脉冲产生的磁场为射频磁场。射频磁场的方向与静磁场方向垂直，实现对在静磁场中形成的磁化矢量的扳转操作，扳转角度为：θ＝γB₁t_p。其中B₁为射频磁场强度，t_p为射频脉冲的持续时间。因此可通过控制射频脉冲的幅值或持续时间达到改变扳转角的目的。核磁共振脉冲序列由不同数量和频率属性的射频脉冲按照设定时序组成。通过调整脉冲间时间间隔，脉冲角度及脉冲的频率选择性，实现对自旋系统的弛豫、扩散等测量。

磁场梯度。磁场梯度可以记录一定时间内分子沿着梯度方向平均扩散位移，从而计算出分子的自扩散系数。此方法作为一种有效的自扩散系数测量，应用于流体种类识别及样品标定等领域。脉冲磁场梯度由梯度线圈产生，通常在施加过程中考虑脉冲梯度线圈与射频线圈的涡流效应，注意屏蔽效果。

自旋回波。自旋回波为核磁共振测量最常见的一种信号。首先对被测样品施加90°脉冲，将磁化矢量M₀扳转至垂直于静磁场方向的横向平面上。由于分子的扩散及静磁场的空间非均匀性等原因，磁化矢量M₀发生散相。这一段时间如果打开信号采集通道对信号进行采集，得到自由衰减信号。经历一定时间τ之后，施加180°脉冲。散相后的磁化矢量会在同等时间τ之后实现重聚，形成一个回波信号。改回波信号称之为自旋回波信号。自旋回波在核磁共振应用中主要有以下三个方面：(1)通过施加一连串的180°脉冲，反复形成自旋回波，记录回波串信号，此脉冲序列为CPMG脉冲序列。该信号对于研究孔隙介质的横向弛豫特性极为重要，在一定条件下可以得到孔隙尺寸相关信息；(2)在梯度磁场下通过改变梯度幅值或梯度持续时间，记录自旋回波幅值的变化，可以得到流体分子的自扩散系数；(3)通过施加成对的频率编码或相位编码梯度，解析被测样品的空间自旋密度信息，实现核磁共振成像。

弛豫。自旋系统从共振状态恢复至热平衡状态的过程。该过程在不同方向上由纵向弛豫时间T₁或横向弛豫时间T₂表征。T₁又称为自旋-晶格弛豫时间，反映自旋系统与外部环境的能量交换。T₂又称为自旋-自旋弛豫时间，反映自旋系统内部能量损耗。自旋系统弛豫过程可由Bloch方程进行描述。纵向弛豫时间T₁可采用饱和恢复脉冲序列进行测量。通过改变两个脉冲之间的时间间隔τ，记录信号幅值，反映纵向磁化矢量在不同编辑时间下的演化过程：

横向弛豫时间T₂的测量由CPMG脉冲序列完成。发射自旋回波串脉冲序列，观测到自旋回波串的幅值衰减，反映纵向磁化矢量随时间的演化过程：

其中n为采集的回波个数，T_E为回波间隔。

自扩散系数D。反映分子的扩散快慢程度。由于分子的扩散过程为随机运动，在一定时间之后的扩散传播函数或扩散概率密度符合高斯分布。当分子在梯度磁场中扩散，其一定时间内信号的改变与分子的平均扩散位移有关，可通过这一规律计算得到分子的自扩散系数。通常采用脉冲磁场梯度或静磁场梯度实现扩散系数的测量。以自旋回波脉冲序列为例，利用积分计算梯度在特定时间内对磁化矢量相位的影响。磁化矢量随脉冲参数的变化规律为：

其中，γ为质子的旋磁比，2π×42.58MHz/T，g为静磁场梯度值，Δ为脉冲梯度对的间隔，δ为单个脉冲梯度持续时间。因此，通过测量磁化矢量在磁场梯度存在下的衰减率，可以得出流体分子在自由状态下的弛豫及扩散特征。

磁化系数χ。磁化系数反映材料本身对静磁场分布的影响情况。孔隙介质的FID信号反映的是孔隙内磁化矢量受弛豫和磁场非均匀性影响衰减过程，衰减率1/T₂ ^*由以下公式可得：

通过标准水样对弛豫T₂及外部静磁场非均匀性γΔB₀进行标定后，就可以得到被测样品与内部填充流体的磁化系数差异Δχ。当已知流体的磁化系数后，即可求得固体骨架材料的磁化系数。

多维核磁共振技术的脉冲序列由数个编辑窗口组成，通过改变几个窗口的编辑参数，采集得到多维核磁共振数据。多维核磁共振数据衰减规律由核矩阵决定。核函数通常为指数函数，表征磁化矢量在不同维度上的变化规律。通过对采集数据进行数据处理，可以得到多维核磁共振特性的分布函数。

实施例一

根据测量被测样品纵向弛豫时间T₁、横向弛豫时间T₂和固体骨架的磁化系数χ参数规律，设计如图1所示的核磁共振脉冲序列，具体如下：

步骤1、向被测样品施加90°脉冲将所述宏观磁化强度矢量M₀扳转90°；

步骤2、等待时间τ过后，系统向被测样品施加第二个90°脉冲，将恢复一定量的磁化矢量从沿着静磁场B₀相同的方向再次扳转90°；

步骤3、在整个系统等待T_E/2的时间之后，向被测样品施加180°射频脉冲，相同的等待时间T_E/2过后，会在ACQ通道采集产生自旋回波信号；

步骤4、重复施加180°射频脉冲，会在ACQ通道采集产生重复多个自旋回波信号，称之为回波串信号；

步骤5、改变等待时间τ，重复所述步骤2、3和4以分别采集数个不同等待时间τ下产生的回波串；

步骤6、根据采集到的所述回波串信号进行核磁共振数据处理。

通过采集信号可得以下响应公式：

M(τ，nT_E，mt_s)＝∫∫∫K₁K₂K₃F(T₁，T₂，Δχ)dT₁·dT₂·dΔχ

三个核函数K₁，K₂，K₃的具体形式为：

K₁＝1-exp(-τ/T₁)

K₂＝exp(-nT_E/T₂)

K₃＝exp(-γ·Δχ·B₀·mt_s)

通过改变T₁编辑时间τ，T₂编辑中的180°脉冲个数n，以及Δχ编辑中的采集点数m获取三维数据。采用快速三维数据处理算法对其进行反演，通过选取合适的正则化因子，得到被测样品的三维T₁-T₂-Δχ分布信息。

实施例二

在实施例一的基础上，根据测量被测样品自扩散系数D、纵向弛豫时间T₁和固体骨架的磁化系数χ参数规律，设计如图2所示的核磁共振脉冲序列，具体为：

步骤1、在TRS通道向被测样品施加90°脉冲将所述宏观磁化强度矢量M₀扳转90°；

步骤2、在GRD通道施加一个持续时间为δ的梯度脉冲；

步骤3、随后在TRS通道施加180°射频脉冲后，在GRD通道施加第二个持续时间为δ的梯度。本步骤中的梯度脉冲起始时间与步骤2中的梯度脉冲起始相差Δ时间；

步骤4、在整个系统等待前两个射频脉冲间隔相等的时间后，向被测样品施加90°射频脉冲，将一定量的磁化矢量从垂直静磁场B₀相同的方向扳转90°；

步骤5、施加小角度射频脉冲，在ACQ通道采集自由衰减信号；

步骤6、等待时间t_a后，重复施加小角度射频脉冲，会在ACQ通道重复采集到自由衰减信号；

步骤7、改变GRD通道中的两个梯度脉冲的强度值，重复所述步骤4、5和6以分别采集不同脉冲强度下产生的自由衰减信号；

步骤8、根据采集到的所述自由衰减串信号进行核磁共振数据处理。

通过采集信号可得以下响应公式：

M(g，nt_a，mt_s)＝∫∫∫K₁K₂K₃F(D，T₁，Δχ)dD·dT₁·dΔχ

三个核函数K₁，K₂，K₃的具体形式为：

K₁＝exp(-Dγ²g²δ²(Δ-δ/3))

K₂＝exp(-nt_a/T₁)

K₃＝exp(-γ·Δχ·B₀·mt_s)

通过改变D编辑中的脉冲梯度强度g，T₁编辑中的小角度脉冲个数n，以及Δχ编辑中的采集点数m获取三维数据。采用快速三维数据处理算法对其进行反演，通过选取合适的正则化因子，得到被测样品的三维D-T₁-Δχ分布信息。

实施例三

在上述实施例的基础上，根据测量被测样品纵向弛豫时间T₁、自扩散系数D和横向弛豫时间T₂参数规律，设计如图3所示的核磁共振脉冲序列，具体如下：

步骤1、向被测样品施加90°脉冲将所述宏观磁化强度矢量M₀扳转90°；

步骤2、等待时间τ过后，系统向被测样品施加第二个90°脉冲，将恢复一定量的磁化矢量从沿着静磁场B₀相同的方向再次扳转90°；

步骤3、在GRD通道施加一个持续时间为δ的梯度脉冲；

步骤4、随后在TRS通道施加180°射频脉冲后，在GRD通道施加第二个持续时间为δ的梯度。本步骤中的梯度脉冲起始时间与步骤2中的梯度脉冲起始相差Δ时间；

步骤5、在整个系统等待T_E/2的时间之后，向被测样品施加180°射频脉冲，相同的等待时间T_E/2过后，会在ACQ通道采集产生自旋回波信号；

步骤6、重复施加180°射频脉冲，会在ACQ通道采集产生重复多个自旋回波信号，称之为回波串信号，记录回波串信号的尖峰值；

步骤7、改变GRD通道中的两个梯度脉冲的强度值，重复所述步骤4-6以分别采集不同脉冲强度下产生的回波串信号峰值；

步骤8、改变等待时间τ，重复所述步骤7以分别采集不同等待时间下产生的回波串信号峰值；

步骤9、根据采集到的所述自由衰减串信号进行核磁共振数据处理。

通过采集信号可得以下响应公式：

M(τ，g，nT_E)＝∫∫∫K₁K₂K₃F(T₁，D，T₂)dT₁·dD·dT₂

三个核函数K₁，K₂，K₃的具体形式为：

其中γ为质子的旋磁比，g为脉冲梯度场梯度强度值。通过改变T₁编辑时间τ，D编辑中的脉冲梯度强度g和T₂编辑中的180°脉冲个数n获取三维数据。采用快速三维数据处理算法对其进行反演，通过选取合适的正则化因子，得到被测样品的三维T₁-D-T₂分布信息。

实施例四

在实施例一到实施例三的基础上，对上述实施例所得到的多维核磁共振数据进行处理。结合图4对提出的多维核磁共振数据的处理算法进行阐述：

步骤1：采用数学张量乘积将前两个核矩阵函数K₁和K₂耦合为一个新的核函数矩阵K₁₂：

步骤2：将测量得到的三维核磁共振数据重新表达为：

M＝K₁₂FK₃

步骤3：对得到的两个核函数矩阵进行SVD分解及奇异值截取，进而对采集数据进行压缩处理。进行数据压缩之前，先对核磁共振数据反演过程的病态程度进行分析。矩阵的病态程度与奇异值有关。数据快速衰减为零，核函数矩阵的奇异值矩阵中的对角元素也快速衰减为零。如果在反演过程中仍然考虑全部的奇异值，整个核函数矩阵的条件数会非常大，造成反演问题病态的严重程度较高。因此利用设定需要的条件数截取奇异值，降低核函数的病态程度。对核函数矩阵进行奇异值分解可得：

K₁₂＝U₁₂·S₁₂·V′₁₂

K₃＝U₃·S₃·V′₃

其中S₁₂和S₃对角线元素值从大到小排列，且为对角矩阵，大小分别为s₁₂×s₁₂和s₃×s₃。其中s₁₂为K₁₂非零奇异值个数，s₃为K₃非零奇异值个数。U₁₂，V₁₂和U₃，V₃为正交单位阵。对对角矩阵S₁₂和S₃进行截取，使得核函数矩阵的条件数满足设定值C，即：

假设C为1000；和分别对应K₁₂和K₃最大的奇异值，即对角矩阵S₁₂和S₃的第一个对角线元素，表示K₁₂的第i个奇异值，表示K₃的第j个奇异值；

步骤4：用截取之后的奇异值分解的单位矩阵对实测数据进行压缩，降低数据内存。由于采用的是单位矩阵，因此压缩之后的数据与原始数据比对，没有丢失信息。压缩后的磁化矢量为：

步骤5：数据压缩完成之后，引入正则化项对数据矩阵进行反演。为了得到稳定准确的解f，通常采用Tikhonov正则化方法，引入正则化项项来通过测量得到的核磁共振数据矩阵M得到核磁共振特性矩阵F：

其中，α是正则化因子，与采集数据的信噪比相关，||·||项代表矩阵的Frobenius范数，K代表核磁共振信号衰减的规律，也就是反演过程中的核矩阵的一般形式。引入的正则化项决定求解结果的稳定性与准确性。正则化因子选取过大，尽管求解得到的分布越稳定，但是解的准确性越差，即所谓的过平滑；正则化因子选取过小，解的求取越准确，但是解的稳定性降低，出现的伪信号越多，即欠平滑。因此，综合考虑解的真实性和解的稳定性，使用合理的正则化因子，是该方法的重点。通过非负约束步骤可以得到特定正则化因子α下的非负约束解f，得到的解通过以下公式得到求解结果与测量结果的残差分布：

χ(α)＝||M-K·f(α)||²

通常最优的正则化因子α的选取标准为

步骤6：通过得到最优正则化因子α，得到最终解f_r，求解公式如下所示：

其中，为压缩截取后残余矩阵的张量积矩阵；为压缩截取后残余矩阵的张量积矩阵的转置矩阵，I是单元矩阵。为压缩后的磁化矢量矩阵的一维矩阵形式。

本发明提出的多维核磁共振数据处理算法的优点是依照核函数的大小对数据进行处理，有效的利用了矩阵信息，同时简化了高维度数据的反演难度，保证了反演的速度。

如上所述，可较好地实现本发明。在不脱离本发明的原理和精神的情况下对这些实施例进行变化、修改、替换、整合和变型仍落入本发明的保护范围内。

Claims (3)

1.一种多维核磁共振测量方法，所述方法用于分析材料孔隙结构和内部填充流体特性，其特征在于，所述方法包括：

(1)多维核磁共振脉冲序列的设计及多维核磁共振数据采集；

(2)对多维核磁共振数据进行反演和解释；

其中，所述步骤(1)包括如下步骤：

1.1向被测样品施加90°射频脉冲将宏观磁化强度矢量M₀扳转90°；

1.2等待时间τ过后，向被测样品施加第二个90°脉冲，将恢复一定量的磁化矢量从沿着静磁场B₀相同的方向再次扳转90°；

1.3等待T_E/2的时间之后，其中，T_E为回波间隔，向被测样品施加180°射频脉冲，等待时间T_E/2过后，会在ACQ通道采集产生自旋回波信号，其中，ACQ通道代表信号接受通道；

1.4重复施加多个180°射频脉冲，会在ACQ通道采集产生重复多个自旋回波信号，得到回波串信号；

1.5改变所述等待时间τ，重复所述步骤1.2-1.4以分别采集数个不同等待时间τ下产生的回波串信号；所述步骤1.5所产生的回波串信号磁化矢量矩阵为：

M(τ,nT_E,mt_s)＝∫∫∫K₁K₂K₃F(T₁,T₂,Δχ)dT₁·dT₂·dΔχ

其中，τ为等待时间，n为180°脉冲个数，m为FID采集点数，t_s为FID采集点时间间隔，F(T₁,T₂,Δχ)为被测样品的三维T₁–T₂-Δχ特性矩阵，T₁为纵向弛豫时间，T₂为横向弛豫时间，Δχ为被测样品与内部填充流体的磁化系数差异，K₁，K₂，K₃为三个核函数，其具体形式为：

K₁＝1-exp(-τ/T₁)

K₂＝exp(-nT_E/T₂)

K₃＝exp(-γ·Δχ·B₀·mt_s)

其中，γ为质子的旋磁比，B₀为静磁场强度。

2.根据权利要求1所述的多维核磁共振测量方法，其特征在于，所述步骤(2)采用快速三维数据处理算法对多维核磁共振数据进行反演。

3.根据权利要求1所述的多维核磁共振测量方法，其特征在于，所述步骤(2)具体为：

2.1采用数学张量乘积将前两个核矩阵函数K₁和K₂耦合为一个新的核函数矩阵K₁₂：

2.2将测量得到的回波串信号磁化矢量矩阵重新表达为：

M＝K₁₂FK₃

2.3对上述核函数矩阵K₁₂和K₃进行SVD分解及奇异值截取，进而对采集数据进行压缩处理，对核函数矩阵进行奇异值分解可得：

K₁₂＝U₁₂·S₁₂·V′₁₂

K₃＝U₃·S₃·V′₃

其中S₁₂和S₃对角线元素值从大到小排列，且为对角矩阵，大小分别为s₁₂×s₁₂和s₃×s₃，其中s₁₂为K₁₂非零奇异值个数，s₃为K₃非零奇异值个数；U₁₂、V₁₂和U₃、V₃为正交单位阵；对对角矩阵S₁₂和S₃进行截取，使得核函数矩阵的条件数满足设定值C，即：

假设C为1000；和分别对应K₁₂和K₃最大的奇异值，即对角矩阵S₁₂和S₃的第一个对角线元素，表示K₁₂的第i个奇异值，表示K₃的第j个奇异值；

2.4用截取之后的奇异值分解的单位矩阵对回波串信号磁化矢量矩阵M进行压缩，降低数据内存，压缩后的磁化矢量为：

其中，为压缩后的磁化矢量矩阵，为两个核函数K₁₂和K₃经过SVD分解和奇异值截取后的残余矩阵，U′₁₂、V′₁₂分别为K₁₂经过SVD分解和奇异值截取后的矩阵U₁₂、V₁的转置，V′₃为K₃经过SVD分解和奇异值截取后的矩阵V₃的转置；

2.5数据压缩完成之后，采用Tikhonov正则化方法对数据矩阵进行反演，其中，正则化项为：

其中，α是正则化因子，与采集数据的信噪比相关，K为反演过程中的核矩阵的一般形式，||·||项代表矩阵的Frobenius范数；

2.6通过选取最优的正则化因子α，得到最终解f_r，求解公式如下：

其中，为压缩截取后残余矩阵的张量积矩阵，为的转置矩阵，α_opt为最优正则化因子，I是单元矩阵，为压缩后的磁化矢量矩阵的一维矩阵形式，