技术领域
[0001] 本
发明涉及核磁共振领域，特别涉及一套核磁共振弛豫高低本征模态耦合新方法的原理和实际实现过程。技术背景
[0002] 目前各种物理方法上的多维实验，如光学，
中子散射，
超声波及核磁共振等领域，通过在不同维度上将对应的物理量进行耦合关联，可以得到对检测样品更为详细的信息。核磁共振技术本身与分子位移密切相关。因而对于孔隙介质来说，通过分析不同时域阶段核磁共振多维
信号，可以得到孔隙介质的孔隙结构等重要信息。尽管如此，多维技术通常会遇到特殊的技术问题而导致得到的图谱的幅度为负数。这样的技术问题主要是由于不同时域间非
正交特征解造成的。
[0003] 核磁共振技术已经演化为一种重要的研究分子动
力学手段。核磁共振技术与分子热运动扩散本身关系密切，因此当用这样的一种技术来研究孔隙介质中饱含
流体分子的扩散特性，进而反推出孔隙介质的孔隙特性成为一种趋势。另外，在求解过程中也可以采用特殊的边界条件来描述流体分子与固体孔隙骨架间的相互碰撞关系。尽管如此，如果多维核磁共振技术中不同时域阶段对应的本征函数为非正交函数时，整体核磁共振响应会变得极为复杂，有可能在最终的响应图谱中造成负值。
发明内容
[0004] 为了解决这一问题，在本发明中提出一种核磁共振弛豫高低本征模态耦合新方法。所述方法包括，
数据采集步骤，数据提取步骤和数据反演步骤。
[0005] 优选的所述数据采集步骤，具体包括，
[0006] 步骤1.1、在TRS通道中向被测样品施加90°射频脉冲将宏观磁化强度矢量M0扳转90°；所述被测样品内的孔隙空间内饱含流体分子；
[0007] 步骤1.2、等待时间te后，在TRS通道中向被测样品施加第二个90°射频脉冲，将恢复一定量的所述磁化强度矢量M0再次扳转90°；
[0008] 步骤1.3、等待时间tdiff，所述流体分子在所述孔隙空间内发生
热力学扩散和/或在所述孔隙空间表面发生碰撞；
[0009] 步骤1.4、之后在TRS通道中向所述被测样品施加第三个90°射频脉冲，将所述流体分子扩散后演化的所述磁化强度矢量M0再次扳转90°至横向磁化矢量平面上；
[0010] 步骤1.5、继续等待时间te后，在TRS通道中系统向所述被测样品施加第四个90°射频脉冲，将恢复一定量的所述磁化强度矢量M0再次扳转90°至纵向磁化矢量方向上；
[0011] 步骤1.6、随后在TRS通道施加小
角度α射频脉冲，开始操作扳转纵向方向上的磁化强度矢量M0，在ACQ通道采集自由衰减信号；
[0012] 步骤1.7、在GRD通道上施加一个梯度脉冲，所述梯度脉冲用于
加速消除残留于横向磁化矢量平面上的磁化强度矢量；从所述步骤1.6中在TRS通道施加小角度α射频脉冲始至在GRD通道上施加梯度脉冲止，整个时序
片段持续时间为ta；
[0013] 步骤1.8、重复所述时序片段N次，共采集得到N个自由衰减信号；
[0014] 步骤1.9、按照对数布点规律，系统性延长步骤1.3中的等待时间tdiffL次并重复所述步骤1.8，共采集得到L*N个自由衰减信号，得到自由衰减串信号E(tdiff,Nta)。
[0015] 优选的，所述数据提取步骤，具体包括，
[0016] 步骤2.1、绘制所述自由衰减串信号E(tdiff,Nta)的首行一维数据和首列一维数据；
[0017] 步骤2.2、在对数坐标轴下，将所述首行一维数据和所述首列一维数据的收敛值进行做差运算，得到比例因子I0,0；
[0018] 步骤2.3、将所述首行一维数据进行拉普拉斯逆变换得到特征值分布函数，通过上述分布函数对首列一维数据进行重构得到新的向量，利用所述比例因子I0,0，将所述新的向量乘以所述比例因子I0,0后与所述首列一维数据进行做差，得到新的首列一维数据；
[0019] 步骤2.4、重复上述步骤2.3，将每行一维数据经处理后乘以比例因子I0,0与每列一维数据进行做差，得到的各列一维数据组合构成新的矩阵，所述新的矩阵即为最低本征模态与一阶高模态的耦合信号M(tdiff,Nta)。
[0020] 优选的，所述数据反演步骤，具体包括，
[0021] 步骤3.1、将所述耦合信号M(tdiff,Nta)重新表达为：
[0022] M＝K1FK2
[0023] 其中，K1,K2为核函数，具体为：
[0024]
[0025]
[0026] 其中， 和 为最低T1本征模态和一阶T1高本征模态下特征函数对应的特征值；
[0027] 步骤3.2、对得到的所述核函数K1,K2进行SVD分解及奇异值截取，对核函数矩阵进行奇异值分解可得：
[0028] K1＝U1·S1·V′1
[0029] K2＝U2·S2·V′2
[0030] 其中，S1和S2对角线元素值从大到小排列，且为对角矩阵，大小分别为s1×s1和s2×s2；其中s1为K1非零奇异值个数，s2为K2非零奇异值个数；U1，V1和U2，V2为正交单位矩阵，V′1为V1的转置矩阵，V′2为V2的转置矩阵；对对角矩阵S1和S2进行截取，使得核函数矩阵K1,K2的条件数满足设定值C，即：
[0031]
[0032] 假设C为1000，σ1max和σ2max分别对应K1和K2的最大奇异值，即对角矩阵S1和S2的第一个对角线元素，σ1(i)表示K1的第i个奇异值，σ2(j)表示K2的第j个奇异值；
[0033] 步骤3.3、用截取之后的奇异值分解的正交单位矩阵U1，V1和U2，V2对所述耦合信号M(tdiff,Nta)进行压缩，降低数据内存；由于采用的是单位矩阵，因此压缩之后的数据与原始数据比对，没有丢失信息；压缩后的耦合信号为
[0034] 步骤3.4、
数据压缩完成之后，采用Tikhonov正则化方法对数据矩阵 进行反演，引入正则化项来求解，所述正则化项为：
[0035]
[0036] 其中，p是正则化因子，与采集数据的
信噪比相关，||·||项代表矩阵的Frobenius范数；
[0037] 步骤3.5、获取Tikhonov正则化方法对应的最优正则化因子，包括：
[0038] 通过非负约束步骤可以得到正则化因子p下的非负约束解f，所述非负约束解f通过以下公式得到其与测量结果M的残差分布：
[0039] χ(p)＝‖M-Kf(p)‖2
[0040] 最优的平滑因子p的选取标准公式为:
[0041]
[0042] 不断更新
迭代正则化因子p，直至满足上述选取标准公式，此时的正则化因子p即为最优正则化因子popt；
[0043] 步骤3.6、通过得到的最优正则化因子popt，得到最终解fr，求解公式如下所示：
[0044]
[0045] 其中， 表示压缩截取后残余矩阵的张量积矩阵。
[0046] 该方法首次将同一弛豫类型的不同本征模态关联在二维技术中。由于同一弛豫类型的不同本征函数在天然上为正交关系，因此可避免在二维数据采集及
数据处理中负图谱幅度的出现，为后期在孔隙介质应用中获取孔隙特性信息提供更为准确的途径。
[0047] 通过在求解过程中引入正则化项，对反演过程进行稳定处理。并通过对核函数分别压缩后进行张量积，对得到的新的核函数进行奇异值分解，截取满足设定条件数的奇异值，同时采用截取后的正交单位矩阵对多维数据进行压缩，加快反演速度。
附图说明
[0048] 图1为本发明核磁共振方法的数据采集脉冲序列图；
[0049] 图2为本发明核磁共振方法的数据提取处理
流程图；
[0050] 图3为本发明核磁共振方法的数据反演
算法处理流程图；
[0051] 其中，TRS代表核磁共振系统的脉冲发射通道，GRD代表梯度脉冲发射通道，ACQ代表核磁共振系统的信号接受通道。
具体实施方式
[0052] 结合
说明书附图说明本发明的具体实施方式。
[0053] 本发明涉及的一种核磁共振弛豫高低本征模态耦合新方法均以核磁共振理论为
基础。因此，下面首先对核磁共振的基本原理进行介绍。
[0054] 静
磁场B0。静磁场由磁体提供，决定核磁共振信号的信噪比。被测样品置于静磁场中，自旋系统内发生能级分裂，沿着静磁场方向会产生一个宏观磁化矢量M0。M0由静磁场强度B0，
温度等参数决定。
[0055] 射频磁场B1与射频脉冲。射频脉冲为电磁信号，通常由线圈产生。射频脉冲产生的磁场为射频磁场。射频磁场的方向与静磁场方向垂直，实现对在静磁场中形成的磁化矢量的扳转操作，扳转角度为：θ＝γB1tp。其中γ为质子的磁旋比，B1为射频磁场强度，tp为射频脉冲的持续时间。因此可通过控制射频脉冲的幅值或持续时间达到改变扳转角的目的。
[0056]
磁场梯度。磁场梯度可以记录一定时间内分子沿着梯度方向平均扩散位移，从而计算出分子的自扩散系数。此方法作为一种有效的自扩散系数测量，应用于流体种类识别及样品标定等领域。
[0057] 自由衰减信号。自由衰减信号为核磁共振测量最常见的一种信号。首先对被测样品施加90°射频脉冲，将磁化矢量M0扳转至垂直于静磁场方向的横向平面上。由于分子的扩散及静磁场的空间非均匀性等原因，磁化矢量M0发生散相。这一段时间如果打开信号采集通道对信号进行采集，得到自由衰减信号。
[0058] 弛豫。自旋系统从共振状态恢复至热平衡状态的过程。该过程在不同方向上由纵向弛豫时间T1或横向弛豫时间T2表征。T1又称为自旋-晶格弛豫时间，反映自旋系统与外部环境的
能量交换。T2又称为自旋-自旋弛豫时间，反映自旋系统内部能量损耗。自旋系统弛豫过程可由Bloch方程进行描述。纵向弛豫时间T1可采用饱和恢复脉冲序列进行测量。通过改变两个脉冲之间的时间间隔τ，记录信号幅值，反映纵向磁化矢量在不同编辑时间下的演化过程：
[0059]
[0060] 横向弛豫时间T2的测量由CPMG脉冲序列完成。发射自旋回波串脉冲序列，观测到自旋回波串的幅值衰减，反映纵向磁化矢量随时间的演化过程：
[0061]
[0062] 其中M0为初始时刻的磁化强度，n为采集的回波个数，TE为回波间隔。
[0063] 弛豫模态。孔隙介质内流体分子的扩散机制与自由流体分子不同，除了与分子固有扩散特性有关外，还与分子与孔隙表面的相互作用有关。由Bloch-Torrey方程可知，单位体积孔隙介质中磁化矢量变化可由下式表述：
[0064]
[0065] 其中，m(r,t)为时间t时孔隙空间内任意点的磁化矢量，r为孔隙内部空间的坐标，D为流体分子的自扩散系数，ν为流体的核磁共振体弛豫率。采用分离变量法可对上公式进行求解，得到具体解形式为：
[0066]
[0067] 其中， 为求解结果中的特征函数，An为特征函数 对应幅度，τn为特征函数 对应的特征值，n为一系列特征函数/特征解的序号。An的具体表达形式如下所示：
[0068]
[0069] 最终，通过对以上m(r,t)进行一维积分，得到测量的磁化矢量随时间变化规律为：
[0070]
[0071] 其中，m0为此单位体积孔隙空间的初始时刻的磁化强度，In为求解结果中的不同本征模态对应的相对幅度，τn为特征函数 对应的特征值。序号n＝0情况下对应的是低本征模态，n>0情况下对应的是高本征模态。
[0072] 对于多维核磁共振弛豫实验来说，在不同时序片段中对弛豫本征模态进行编辑，得到的磁化矢量衰减规律为：
[0073]
[0074] 其中，Pα,Pβ代表不同时序片段的持续时间；Im,n代表不同模态前后耦合关联的相对幅度值， 和 为前后片段对应的特征值。对于多维核磁共振弛豫测量技术来说，正是由于在不同时序片段中涉及的弛豫本征模态相互之间为非严格正交关系，Im,n幅度求解过程中 项无法保证大等于
0，因此会在图谱中出现一些负值结果。
[0076] 本发明提出了一种能够耦合单一弛豫的高低本征模态方法。鉴于单一弛豫的不同模态间能够保持正交关系，因此上述由于非正交问题而引起的图谱负值问题可以得到根本解决。具体实现的测量方法如图1脉冲序列所示。
[0077] 步骤1、在TRS通道中向被测样品施加90°射频脉冲将所述宏观磁化强度矢量M0扳转90°；所述被测样品的空隙空间内饱含流体分子；
[0078] 步骤2、等待时间te后，系统在TRS通道中向被测样品施加第二个90°脉冲，将恢复一定量的磁化矢量再次扳转90°；
[0079] 其中，所述时间te的选取与使用的测量仪器磁场强度有关，在几百微秒到几毫秒之间，优选为200微秒到8毫秒之间。
[0080] 步骤3、整个自旋质子系统等待tdiff时间，分子在孔隙空间内发生热力学扩散并有可能在孔隙表面发生碰撞；所述时间tdiff的长度为几毫秒到几秒之间，优选为2毫秒到5秒之间。
[0081] 步骤4、之后在TRS通道中向被测样品施加第三个90°射频脉冲，将分子扩散后演化的磁化矢量再次扳转90°至横向磁化矢量平面上；
[0082] 步骤5、继续等待时间te后，在TRS通道中系统向被测样品施加第四个90°脉冲，将恢复一定量的磁化矢量再次扳转90°至纵向磁化矢量方向上；
[0083] 步骤6、随后在TRS通道施加小角度α射频脉冲，将整个质子自旋系统中的纵向方向上的磁化矢量M0扳转α度，随后在ACQ通道采集自由衰减信号；其中，角度α小于15度；
[0084] 步骤7、在GRD通道上施加一个梯度脉冲。该梯度脉冲用于加速消除残留于横向平面上的磁化矢量，从TRS通道施加小角度α射频脉冲至GRD通道梯度脉冲结束的持续时序时间为ta；
[0085] 步骤8、重复TRS通道小角度α射频脉冲开始至GRD通道梯度脉冲结束这一时序片段N次，共采集得到自由衰减信号N个；
[0086] 步骤9、按照对数布点规律系统性的延长步骤3中的自旋质子系统等待时间tdiffL次并重复步骤8，整体实验中共采集得到自由衰减信号L*N个；
[0087] 步骤10、对采集到的L*N个自由衰减串信号进行核磁共振数据处理和分析，得到弛豫高低模态耦合信息。
[0088] 将图1所示脉冲序列按照循环时序分为前后两个片段。在循环前的4个90°射频脉冲期间，自旋质子经历的是纵向弛豫时间T1本征模态。受到te时间段的特殊激发，本处的磁化矢量将主要由最低本征模态和一阶高本征模态贡献。这一部分的信号演化规律为：
[0089]
[0090] 其中，I0和I1代表最低T1本征模态和一阶T1高本征模态的相对幅度值， 和 为两种模态对应的特征值。
[0091] 在循环后的一系列小角度射频脉冲施加期间，由于磁化矢量记录的为纵向方向上的磁化矢量，因此自旋质子经历的是仍然是纵向弛豫时间T1本征模态。本处自旋质子未受到任何特殊激发效应，因此磁化矢量将主要由最低本征模态贡献。这一部分的信号演化规律在
相位循环前为：
[0092]
[0093]
[0094] 在相位循环后为： 其中α为小角度射频脉冲磁化矢量扳转角度。
[0095] 最终本二维实验的磁化矢量演化规律可被描述为：
[0096]
[0097] 根据不同时序片段内T1本征模态的贡献情况，可将上公式简化为：
[0098]
[0099] 实施例二
[0100] 在实施例一的基础上，对核磁共振测量得到自由衰减串信号E(tdiff,Nta)进行提取处理，具体的步骤为：
[0101] 步骤2.1、绘制所述自由衰减串信号E(tdiff,Nta)的首行一维数据和首列一维数据；
[0102] 步骤2.2、在对数坐标轴下，将所述首行一维数据和所述首列一维数据的收敛值进行做差运算，得到比例因子I0,0；
[0103] 步骤2.3、将所述首行一维数据进行拉普拉斯逆变换得到特征值分布函数，通过上述分布函数对首列一维数据进行重构得到新的向量，利用所述比例因子I0,0，将所述新的向量乘以所述比例因子I0,0后与所述首列一维数据进行做差，得到新的首列一维数据；
[0104] 步骤2.4、重复上述步骤2.3，将每行一维数据经处理后乘以比例因子I0,0与每列一维数据进行做差，得到的各列一维数据组合构成新的矩阵，所述新的矩阵即为最低本征模态与一阶高模态的耦合信号M(tdiff,Nta)。
[0105] 经过上述处理，可将上述公式内前一项最低本征模态自耦合项进行消除，最终得到最低本征模态与一阶高模态的耦合信号M(tdiff,Nta)：
[0106]
[0107] 对于复杂的自旋质子系统来讲，上述公式为：
[0108]
[0109] 其中，核函数K1,K2的具体形式为：
[0110]
[0111]
[0112] 实施例三
[0113] 在上述实施例的基础上，对提取得到的最低本征模态与一阶高模态的耦合信号M(tdiff,Nta)，进行数据处理，具体步骤为：
[0114] 步骤1：采用数学表述将上述处理得到的数据M(tdiff,Nta)重新表达为：
[0115] M＝K1FK2
[0116] 步骤2：对得到的两个核函数矩阵进行SVD分解及奇异值截取，进而对采集数据进行压缩处理。进行数据压缩之前，先对核磁共振数据反演过程的病态程度进行分析。矩阵的病态程度与奇异值有关。数据快速衰减为零，核函数矩阵的奇异值矩阵中的对角元素也快速衰减为零。如果在反演过程中仍然考虑全部的奇异值，整个核函数矩阵的条件数会非常大，造成反演问题病态的严重程度较高。因此利用设定需要的条件数截取奇异值，降低核函数的病态程度。对核函数矩阵进行奇异值分解可得：
[0117] K1＝U1·S1·V′1
[0118] K2＝U2·S2·V′2
[0119] 其中S1和S2对角线元素值从大到小排列，且为对角矩阵，大小分别为s1×s1和s2×s2。其中s1为K1非零奇异值个数，s2为K2非零奇异值个数。U1，V1和U2，V2为正交单位阵。对对角矩阵S1和S2进行截取，使得核函数矩阵的条件数满足设定值C，即：
[0120]
[0121] 通常设定C为1000。σ1max和σ2max分别对应K1和K2最大的奇异值，即对角矩阵S1和S2的第一个对角线元素。σ1(i)表示K1的第i个奇异值，σ2(j)表示K2的第j个奇异值。
[0122] 步骤3：用截取之后的奇异值分解的单位矩阵对实测数据进行压缩，降低数据内存。由于采用的是单位矩阵，因此压缩之后的数据与原始数据比对，没有丢失信息。压缩后的磁化矢量为
[0123] 步骤4：数据压缩完成之后，引入正则化项对数据矩阵进行反演。为了得到稳定准确的解f，通常采用Tikhonov正则化方法，引入平滑项来求解该问题：
[0124]
[0125] 其中，p是正则化因子，与采集数据的信噪比相关，||·||项代表矩阵的Frobenius范数。引入的正则化项决定求解结果的
稳定性与准确性。正则化因子选取过大，尽管求解得到的分布越稳定，但是解的准确性越差，即所谓的过平滑；正则化因子选取过小，解的求取越准确，但是解的稳定性降低，出现的伪信号越多，即欠平滑。因此，综合考虑解的真实性和解的稳定性，使用合理的正则化因子，是该方法的重点。通过非负约束步骤可以得到特定正则化因子p下的非负约束解f，得到的解通过以下公式得到求解结果与测量结果的残差分布：
[0126] χ(p)＝‖M-Kf(p)‖2
[0127] 通常最优的平滑因子p的选取标准为:
[0128]
[0129] 步骤5：通过得到最优平滑因子p，得到最终解fr，fr为最终通过输入最优正则化因子popt，得到的最终非负解。求解公式如下所示：
[0130]
[0131] 其中， 为压缩截取后残余矩阵的张量积矩阵。
[0132] 如上所述，可较好地实现本发明。在不脱离本发明的原理和精神的情况下对这些实施例进行变化、
修改、替换、整合和变型仍落入本发明的保护范围内。
