CN108009125B - 基于l0正则化的核磁共振回波数据反演方法及装置 - Google Patents
基于l0正则化的核磁共振回波数据反演方法及装置 Download PDFInfo
- Publication number
- CN108009125B CN108009125B CN201711346882.7A CN201711346882A CN108009125B CN 108009125 B CN108009125 B CN 108009125B CN 201711346882 A CN201711346882 A CN 201711346882A CN 108009125 B CN108009125 B CN 108009125B
- Authority
- CN
- China
- Prior art keywords
- magnetic resonance
- nuclear magnetic
- determining
- matrix
- dimensional
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000005481 NMR spectroscopy Methods 0.000 title claims abstract description 193
- 238000000034 method Methods 0.000 title claims abstract description 66
- 238000001228 spectrum Methods 0.000 claims abstract description 82
- 102000008297 Nuclear Matrix-Associated Proteins Human genes 0.000 claims abstract description 47
- 108010035916 Nuclear Matrix-Associated Proteins Proteins 0.000 claims abstract description 47
- 210000000299 nuclear matrix Anatomy 0.000 claims abstract description 47
- 239000011159 matrix material Substances 0.000 claims abstract description 45
- 238000005259 measurement Methods 0.000 claims abstract description 45
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 16
- 238000000655 nuclear magnetic resonance spectrum Methods 0.000 claims abstract description 15
- 238000000685 Carr-Purcell-Meiboom-Gill pulse sequence Methods 0.000 claims description 35
- 238000001208 nuclear magnetic resonance pulse sequence Methods 0.000 claims description 33
- 239000013598 vector Substances 0.000 claims description 25
- 239000000126 substance Substances 0.000 claims description 14
- 238000009792 diffusion process Methods 0.000 claims description 12
- 238000010276 construction Methods 0.000 claims description 5
- 238000010606 normalization Methods 0.000 claims description 2
- 230000000694 effects Effects 0.000 abstract description 9
- 238000012545 processing Methods 0.000 description 13
- 238000004088 simulation Methods 0.000 description 8
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 8
- 239000012530 fluid Substances 0.000 description 6
- 238000005084 2D-nuclear magnetic resonance Methods 0.000 description 4
- 238000010586 diagram Methods 0.000 description 4
- 230000015572 biosynthetic process Effects 0.000 description 3
- 230000003595 spectral effect Effects 0.000 description 3
- 238000005516 engineering process Methods 0.000 description 2
- 239000004215 Carbon black (E152) Substances 0.000 description 1
- UFHFLCQGNIYNRP-UHFFFAOYSA-N Hydrogen Chemical compound [H][H] UFHFLCQGNIYNRP-UHFFFAOYSA-N 0.000 description 1
- 238000004364 calculation method Methods 0.000 description 1
- 239000007789 gas Substances 0.000 description 1
- 238000009499 grossing Methods 0.000 description 1
- 229930195733 hydrocarbon Natural products 0.000 description 1
- 150000002430 hydrocarbons Chemical class 0.000 description 1
- 239000001257 hydrogen Substances 0.000 description 1
- 229910052739 hydrogen Inorganic materials 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 210000004940 nucleus Anatomy 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 230000000704 physical effect Effects 0.000 description 1
- 239000011148 porous material Substances 0.000 description 1
- 238000011084 recovery Methods 0.000 description 1
- 238000005070 sampling Methods 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
-
- E—FIXED CONSTRUCTIONS
- E21—EARTH OR ROCK DRILLING; MINING
- E21B—EARTH OR ROCK DRILLING; OBTAINING OIL, GAS, WATER, SOLUBLE OR MELTABLE MATERIALS OR A SLURRY OF MINERALS FROM WELLS
- E21B49/00—Testing the nature of borehole walls; Formation testing; Methods or apparatus for obtaining samples of soil or well fluids, specially adapted to earth drilling or wells
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Mathematical Physics (AREA)
- Pure & Applied Mathematics (AREA)
- Mathematical Optimization (AREA)
- Mining & Mineral Resources (AREA)
- Theoretical Computer Science (AREA)
- Data Mining & Analysis (AREA)
- Computational Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Geology (AREA)
- Mathematical Analysis (AREA)
- Fluid Mechanics (AREA)
- Environmental & Geological Engineering (AREA)
- Computing Systems (AREA)
- Algebra (AREA)
- Geochemistry & Mineralogy (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- General Engineering & Computer Science (AREA)
- General Life Sciences & Earth Sciences (AREA)
- Magnetic Resonance Imaging Apparatus (AREA)
Abstract
本发明提供一种基于L0正则化的核磁共振回波数据反演方法及装置，该方法根据核磁共振测量参数确定核矩阵；根据核矩阵，并以L0范数作为目标函数的惩罚项，构建目标函数；基于迭代算法，确定目标函数的最优解；基于弛豫信息，根据最优解，确定核磁共振回波数据反演的核磁共振谱。本发明的方法与传统核磁共振回波数据反演方法相比，该反演方法利用L0范数对反演结果进行惩罚，提高了短弛豫组分反演谱的效果，且迭代求解过程简单，算法的鲁棒性好，在不同信噪比核磁共振数据情况下，能够得到稳定的反演结果。
Description
技术领域
本发明涉及油气勘探中核磁共振测井数据处理技术，尤其涉及一种基于L0正则化的核磁共振回波数据反演方法及装置。
背景技术
核磁共振测井是一种通过测量地层流体中氢核的弛豫特性，提供反映地层储集特性的测井新技术。核磁共振测井技术可有效应用于储层物性参数的计算，储层流体性质的识别等，对识别低电阻率、低孔低渗油气层、评价储层孔隙结构等方面具有明显优势。
核磁共振测井测量的原始回波数据需要通过反演才能获得反映地层流体信息的核磁共振弛豫谱。现有技术中，核磁共振回波数据反演通常运用正则化方法。正则化的目的是防止解过拟合，提高模型的泛化能力。目前，广泛使用的正则化方法包括L2正则化方法和L1正则化方法。其中，L2正则化方法要求残差小的同时反演谱的L2范数也相对较小，从而容易造成反演谱的过度平滑。L1正则化约束解的稀疏性，可以得到分辨率更高的反演谱。L0正则化是利用非零参数的个数约束解的稀疏性，它希望求取的参数大部分是0，这样才能保证目标函数值小。相较于L1正则化和L2正则化来说，L0正则化方法得到的反演谱的稀疏性强于L2正则化方法和L1正则化方法。反演解的稀疏性越好，反演的谱更易于解释，且反演过程中核磁共振短弛豫组分丢失的信息少。因此，L0正则化方法在核磁共振回波数据反演问题中具有很大的应用前景。
但是，L0正则化问题很难求解，目前运用正交最小二乘法、正交投影寻踪法等方法求出来的解不光滑，造成基于L0正则化的核磁共振回波数据反演谱不准确。因此，亟需一种求解L0正则化的方法以获得可靠的核磁共振回波数据的反演谱。
发明内容
本发明为了解决现有技术中基于L0正则化的核磁共振回波数据反演谱不准确的技术问题，提供一种基于L0正则化的核磁共振回波数据反演方法及装置。
本发明提供一种基于L0正则化的核磁共振回波数据反演方法，包括：
根据核磁共振测量参数确定核矩阵；
根据所述核矩阵，并以L0范数作为目标函数的惩罚项，构建所述目标函数；
基于迭代算法，确定所述目标函数的最优解；
基于弛豫信息，根据所述最优解，确定核磁共振回波数据反演的核磁共振谱。
可选的，所述根据核磁共振测量参数确定核矩阵，包括：
根据核磁共振测量的脉冲序列，确定所述核矩阵；
其中，所述核磁共振测量的脉冲序列为一维核磁共振CPMG脉冲序列、二维核磁共振SR-CPMG脉冲序列或基于变回波间隔TE的二维核磁共振CPMG脉冲序列；
若所述核磁共振测量的脉冲序列为一维核磁共振CPMG脉冲序列，则所述核矩阵为
若所述核磁共振测量的脉冲序列为二维核磁共振SR-CPMG脉冲序列，则所述核矩阵为
其中
表示两个矩阵的Kronecker积；
若所述核磁共振测量的脉冲序列为基于变回波间隔TE的二维核磁共振CPMG脉冲序列，则所述核矩阵为
其中，t为进行核磁共振测量时每个回波对应的时刻；T1为纵向弛豫时间；T2为横向弛豫时间；TW为等待时间；γ为旋磁比；G为磁场梯度；TE为回波间隔；D为扩散系数。
可选的，所述根据所述核矩阵，并以L0范数作为目标函数的惩罚项，构建所述目标函数，包括：
根据所述核矩阵K，构建目标函数：
其中，||·||0表示向量的L0范数，s表示反演的谱，以||s||0作为所述目标函数的惩罚项；m表示核磁共振回波数据；其中，s和m为向量；s≥0表示向量s中的所有元素大于等于0，τ为正则化参数。
可选的，所述基于迭代算法，确定所述目标函数的最优解，包括：
获取噪声信号的标准差σ，根据所述标准差σ，对初始正则化参数τinit进行迭代更新，得到更新后的正则化参数τ；
其中，所述初始正则化参数
所述更新后的正则化参数
式中，N为核磁共振数据点的个数，g为向量s中元素的个数，snew为前一次迭代的解s；
根据预设的初始解s0，确定第一解s1，
其中，
根据所述初始解s0、所述第一解s1，确定第n次迭代的解sn+1，
sn+1＝asn-1+bsn+cΨτ(sn+KT(m-Ksn))；其中，n表示迭代次数，且n≥1；a、b、c为预设的迭代参数；
根据所述解sn+1，确定误差
根据最小误差tol，判断所述误差ζ＜tol的条件是否成立；
若成立，则确定所述目标函数的最优解为sn+1；
若不成立，则更新n＝n+1，并返回执行所述确定第n次迭代的解sn+1的步骤。
可选的，所述基于弛豫信息，根据所述最优解，确定核磁共振回波数据反演的核磁共振谱，包括：
若反演一维核磁共振T2谱时，基于横向弛豫时间T2，根据所述最优解sn+1，确定核磁共振T2谱；
若反演二维核磁共振T2-T1谱时，对所述最优解sn+1进行排序，生成与所述最优解sn+1对应的二维矩阵S1，基于横向弛豫时间T2和纵向弛豫时间T1，根据所述二维矩阵S1，确定核磁共振T2-T1谱；
若反演二维核磁共振T2-D谱时，对所述最优解sn+1进行排序，生成与所述最优解sn+1对应的二维矩阵S2，基于横向弛豫时间T2和扩散系数D，根据所述二维矩阵S2，确定核磁共振T2-D谱。
本发明还提供一种基于L0正则化的核磁共振回波数据反演装置，包括：
确定模块，用于根据核磁共振测量参数确定核矩阵；
构建模块，用于根据所述核矩阵，并以L0范数作为目标函数的惩罚项，构建所述目标函数；
所述确定模块，还用于基于迭代算法，确定所述目标函数的最优解；
所述确定模块，还用于基于弛豫信息，根据所述最优解，确定核磁共振回波数据反演的核磁共振谱。
可选的，所述确定模块包括：
核矩阵确定子模块，用于根据核磁共振测量的脉冲序列，确定所述核矩阵；其中，所述核磁共振测量的脉冲序列为一维核磁共振CPMG脉冲序列、二维核磁共振SR-CPMG脉冲序列或基于变回波间隔TE的二维核磁共振CPMG脉冲序列；
当所述核磁共振测量的脉冲序列为一维核磁共振CPMG脉冲序列时，所述核矩阵确定子模块确定所述核矩阵为
当所述核磁共振测量的脉冲序列为二维核磁共振SR-CPMG脉冲序列时，所述核矩阵确定子模块确定所述核矩阵为
当所述核磁共振测量的脉冲序列为基于变回波间隔TE的二维核磁共振CPMG脉冲序列时，所述核矩阵确定子模块确定所述核矩阵为
可选的，所述构建模块，具体用于根据所述核矩阵K，构建目标函数：
其中，||·||0表示向量的L0范数，s表示反演的谱，以||s||0作为所述目标函数的惩罚项；m表示核磁共振回波数据；其中，s和m为向量；s≥0表示向量s中的所有元素大于等于0，τ为正则化参数。
可选的，所述确定模块包括：
获取子模块，用于获取噪声信号的标准差σ；
更新子模块，用于根据所述标准差σ，对初始正则化参数τinit进行迭代更新，得到更新后的正则化参数τ；其中，所述初始正则化参数
所述更新后的正则化参数
式中，N为核磁共振数据点的个数，g为向量s中元素的个数，snew为前一次迭代的解s；
第一确定子模块，用于根据预设的初始解s0，确定第一解s1，s1＝Ψτ(s0+KT(m-Ks0))；
其中，
第二确定子模块，用于根据所述初始解s0、所述第一解s1，确定第n次迭代的解sn+1，sn+1＝asn-1+bsn+cΨτ(sn+KT(m-Ksn))；其中，n表示迭代次数，且n≥1；a、b、c为预设的迭代参数；
误差确定子模块，用于根据所述解sn+1，确定误差
判断子模块，用于根据最小误差tol，判断所述误差ζ＜tol的条件是否成立；
最优解确定子模块，用于当所述判断子模块判断所述误差ζ＜tol的条件成立时，确定所述目标函数的最优解为sn+1；还用于当所述判断子模块判断所述误差ζ＜tol的条件不成立时，更新n＝n+1，并返回所述第二确定子模块执行所述确定第n次迭代的解sn+1的步骤。
可选的，所述确定模块，包括：
第一反演子模块，用于当反演一维核磁共振T2谱时，基于横向弛豫时间T2，根据所述最优解sn+1，确定核磁共振T2谱；
第二反演子模块，用于当反演二维核磁共振T2-T1谱时，对所述最优解sn+1进行排序，生成与所述最优解sn+1对应的二维矩阵S1，基于横向弛豫时间T2和纵向弛豫时间T1，根据所述二维矩阵S1，确定核磁共振T2-T1谱；
第三反演子模块，用于当反演二维核磁共振T2-D谱时，对所述最优解sn+1进行排序，生成与所述最优解sn+1对应的二维矩阵S2，基于横向弛豫时间T2和扩散系数D，根据所述二维矩阵S2，确定核磁共振T2-D谱。本发明提供的基于L0正则化的核磁共振回波数据反演方法及装置，通过根据核磁共振测量参数确定核矩阵；根据核矩阵，并以L0范数作为目标函数的惩罚项，构建目标函数；基于迭代算法，确定目标函数的最优解；基于弛豫信息，根据最优解，确定核磁共振回波数据反演的核磁共振谱。本发明的方法与传统核磁共振回波数据反演方法相比，该反演方法利用L0范数对反演结果进行惩罚，提高了短弛豫组分反演谱的效果，且迭代求解过程简单，算法的鲁棒性好，在不同信噪比核磁共振数据情况下，能够得到稳定的反演结果。
附图说明
图1为本发明一示例性实施例示出的基于L0正则化的核磁共振回波数据反演方法的流程图；
图2为本发明另一示例性实施例示出的基于L0正则化的核磁共振回波数据反演方法的流程图；
图3是本发明构造的含有束缚水和油的T2谱模型图；
图4(a)是本发明正演的回波间隔为0.2ms，信噪比为15的回波串；
图4(b)是本发明正演的回波间隔为0.2ms，信噪比为30的回波串；
图5(a)是基于图2所示实施例的反演方法处理图4(a)中的数据得到的核磁共振T2谱；
图5(b)是基于图2所示实施例的反演方法处理图4(b)中的数据得到的核磁共振T2谱；
图6是本发明构造的含有束缚水和油的T2-T1谱模型图；
图7(a)是本发明正演的9组不同等待时间的回波串，信噪比为50；
图7(b)是本发明正演的9组不同等待时间的回波串，信噪比为100；
图8(a)是基于图2所示实施例的反演方法处理图7(a)中的数据得到的核磁共振T2-T1谱；
图8(b)是基于图2所示实施例的反演方法处理图7(b)中的数据得到的核磁共振T2-T1谱；
图9是本发明构造的含有束缚水和油的T2-D谱模型图；
图10(a)是本发明正演的9组不同回波间隔的回波串，信噪比为50；
图10(b)是本发明正演的9组不同回波间隔的回波串，信噪比为100；
图11(a)是基于图2所示实施例的反演方法处理图10(a)中的数据得到的核磁共振T2-D谱；
图11(b)是基于图2所示实施例的反演方法处理图10(b)中的数据得到的核磁共振T2-D谱；
图12为本发明一示例性实施例示出的基于L0正则化的核磁共振回波数据反演装置的结构示意图；
图13为本发明另一示例性实施例示出的基于L0正则化的核磁共振回波数据反演装置的结构示意图。
具体实施方式
为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例，对本发明实施例中的技术方案进行清楚、完整地描述。需要说明的是，在附图或说明书中，相似或相同的元件皆使用相同的附图标记。
图1为本发明一示例性实施例示出的基于L0正则化的核磁共振回波数据反演方法的流程图，如图1所示，本实施例的方法包括：
步骤101、根据核磁共振测量参数确定核矩阵。
步骤102、根据核矩阵，并以L0范数作为目标函数的惩罚项，构建目标函数。
步骤103、基于迭代算法，确定目标函数的最优解。
步骤104、基于弛豫信息，根据最优解，确定核磁共振回波数据反演的核磁共振谱。
本实施例的基于L0正则化的核磁共振回波数据反演方法，通过根据核磁共振测量参数确定核矩阵；根据核矩阵，并以L0范数作为目标函数的惩罚项，构建目标函数；基于迭代算法，确定目标函数的最优解；基于弛豫信息，根据最优解，确定核磁共振回波数据反演的核磁共振谱。本发明的方法与传统核磁共振回波数据反演方法相比，该反演方法利用L0范数对反演结果进行惩罚，提高了短弛豫组分反演谱的效果，且迭代求解过程简单，算法的鲁棒性好，在不同信噪比核磁共振数据情况下，能够得到稳定的反演结果。
图2为本发明另一示例性实施例示出的基于L0正则化的核磁共振回波数据反演方法的流程图，如图2所示，基于上一实施例，本实施例的方法包括：
步骤201、根据核磁共振测量的脉冲序列，确定核矩阵。
在本步骤中，根据核磁共振测量参数生成核矩阵K，以及利用L0范数作为目标函数的惩罚项，可以形成L0正则化问题。其中，根据核磁共振测量的脉冲序列以及设定的测量参数，核矩阵K就被确定。不同的采样脉冲序列具有不同形式的核矩阵。例如，核磁共振测量的脉冲序列可以为一维核磁共振CPMG(Carr-Purcell-Meiboom-Gill)脉冲序列，或者为二维核磁共振SR(Saturation-Recovery)-CPMG脉冲序列，或者为基于变回波间隔TE的二维核磁共振CPMG脉冲序列。
相应的，核矩阵的确定方法具体可以包括以下三种情况：
第一种情况，若核磁共振测量的脉冲序列为一维核磁共振CPMG脉冲序列，则核矩阵为：
第二种情况，若核磁共振测量的脉冲序列为二维核磁共振SR-CPMG脉冲序列，则核矩阵为
其中：
第三种情况，若核磁共振测量的脉冲序列为基于变回波间隔TE的二维核磁共振CPMG脉冲序列，则核矩阵为：
其中，上述各个公式中，t为进行核磁共振测量时每个回波对应的时刻；
表示两个矩阵的Kronecker积；T1为纵向弛豫时间；T2为横向弛豫时间；TW为等待时间；γ为旋磁比；G为磁场梯度；TE为回波间隔；D为扩散系数。
步骤202、根据核矩阵K，构建目标函数。
在本步骤中，为了求取核磁共振谱，利用L0范数作为目标函数的惩罚项，其目标函数可以写为：
公式(5)中s表示反演的谱，m表示核磁共振回波数据，s和m都为向量，s≥0表示向量s中的所有元素大于等于0，τ为正则化参数，||·||0表示向量的L0范数，以该||s||0作为该目标函数的惩罚项，即向量中非零元素的个数。K为步骤201中得到的核矩阵。
步骤203、获取噪声信号的标准差σ。
在本步骤中，结合噪声来选取最优正则化参数，首先根据噪声信号，计算得到该噪声的标准差σ。
步骤204、根据标准差σ，对初始正则化参数τinit进行迭代更新，得到更新后的正则化参数τ。
在本步骤中，定义初始的正则化参数为：
其中，N为核磁共振数据点的个数，g为向量s中元素的个数。随后每次迭代过程中的正则化参数依据公式(7)进行更新：
其中，snew为前一次迭代的解。
步骤205、根据预设的初始解s0，确定第一解s1。
在本步骤中，利用迭代算法求取L0正则化问题的最优解的步骤为：
首先定义初始解s0以及最小误差tol，并根据公式(8)计算s1，
s1＝Ψτ(s0+KT(m-Ks0)) 公式(8)
其中，
将s1中值为负的项取值为0。
步骤206、根据初始解s0、第一解s1，确定第n次迭代的解sn+1。
在本步骤中，迭代次数用n表示，每次迭代的解存入sn+1中，若sn+1中有负项，将sn+1中的负项置0。第n(n≥1)次迭代的解可表示为：
sn+1＝asn-1+bsn+cΨτ(sn+KT(m-Ksn)) 公式(9)
式中a，b和c为迭代参数，迭代参数a，b和c的值可以通过大量数值模拟实验计算得出，优选的，a＝0.04，b＝-2.96和c＝3.92。
步骤207、根据解sn+1，确定误差。
在本步骤中，基于步骤206中的解sn+1，计算误差
以判断ζ＜tol是否成立。
步骤208、根据最小误差tol，判断误差ζ＜tol的条件是否成立；若成立执行步骤209；若不成立，更新n＝n+1并返回执行步骤206，即确定下一次迭代的解sn+1的值。
步骤209、确定目标函数的最优解为sn+1，并基于弛豫信息，根据该最优解sn+1，确定核磁共振回波数据反演的核磁共振谱。
在本步骤中，基于反演的核磁共振谱的不同，可以分为以下三种情况得到反演的核磁共振谱：
第一种，若反演一维核磁共振T2谱时，基于横向弛豫时间T2，根据最优解sn+1，确定核磁共振T2谱。
该核磁共振T2谱的反演效果可以通过以下数值模拟实验进一步说明：
首先如图3所示，构造含有束缚水和油的T2谱，两个峰对应的T2值分别为10ms和200ms，然后向正演结果中加入一定的高斯白噪声，模拟得到不同信噪比的回波串数据。例如，图4(a)中回波间隔为0.2ms，信噪比为15；图4(b)中回波间隔为0.2ms，信噪比为30。图5(a)是基于本实施例的反演方法处理图4(a)中的数据得到的核磁共振T2谱，图5(b)是基于本实施例的反演方法处理图4(b)中的数据得到的核磁共振T2谱。从上述各个图中可以看出，采用上述各个步骤的反演方法反演出的核磁共振T2谱与模型谱(图3)接近，误差较小，说明本实施例提出的反演方法能够用于核磁共振回波数据的反演，且反演的准确性好。
第二种，若反演二维核磁共振T2-T1谱时，对最优解sn+1进行排序，生成与该最优解sn+1对应的二维矩阵S1，基于横向弛豫时间T2和纵向弛豫时间T1，根据二维矩阵S1，确定核磁共振T2-T1谱。
该核磁共振T2-T1谱的反演效果可以通过以下数值模拟实验进一步说明：
首先如图6所示，构造含有束缚水和油的T2-T1谱，峰值坐标分别为(T2，T1)＝(5，15)ms和(T2，T1)＝(100，300)ms，然后向正演结果中加入一定的高斯白噪声，模拟得到不同信噪比的回波串数据。图6是数值模拟实验构造的含有束缚水和油的T2-T1谱。图7(a)和图7(b)是正演的回波串数据，其中图7(a)中回波数据的信噪比为50，图7(b)中回波数据的信噪比为100。图8(a)是基于本实施例的反演方法处理图7(a)中的回波数据得到的核磁共振T2-T1谱，图8(b)是基于本实施例的反演方法处理图7(b)中的回波数据得到的核磁共振T2-T1谱。从上述各个图中可以看出，采用上述各个步骤的反演方法可以在不同信噪比条件下都能够很好地反演出真实的T2-T1谱，且能够抑制核磁共振谱中流体信号的重叠。
第三种，若反演二维核磁共振T2-D谱时，对最优解sn+1进行排序，生成与最优解sn+1对应的二维矩阵S2，基于横向弛豫时间T2和扩散系数D，根据二维矩阵S2，确定核磁共振T2-D谱。
该核磁共振T2-D谱的反演效果可以通过以下数值模拟实验进一步说明：
首先如图9所示，构造含有束缚水和油的T2-D谱，峰值坐标分别为(T2，D)＝(5ms，2.0×10-5cm2/s)和(T2，D)＝(100ms，2.5×10-6cm2/s)，然后向正演结果中加入一定的高斯白噪声，模拟得到不同信噪比的回波串数据。图9是数值模拟实验构造的含有束缚水和油的T2-D谱。图10(a)和图10(b)是正演的回波串数据，其中图10(a)中回波数据的信噪比为50，图10(b)中回波数据的信噪比为100。图11(a)是基于本实施例的反演方法处理图10(a)中的回波数据得到的核磁共振T2-D谱，图11(b)是基于本实施例的反演方法处理图10(b)中的回波数据得到的核磁共振T2-D谱。从上述各个图中可以看出，采用上述各个步骤的反演方法可以在不同信噪比条件下都能够很好地反演出真实的T2-D谱，且能够抑制核磁共振谱中不同流体信号的重叠。
本实施例的基于L0正则化的核磁共振回波数据反演方法，其使用了L0范数对反演结果进行惩罚，提高了短弛豫组分反演谱的效果，通过该实施例中各个附图所显示的反演出的核磁共振T2谱、T2-T1谱、T2-D谱可以看出，本实施例的反演方法在不同信噪比条件下都能够很好地反演出真实谱，且能够抑制核磁共振谱中不同流体信号的重叠。同时，该方法与现有的核磁共振回波数据反演方法相比，算法的鲁棒性好，迭代求解过程简单，在不同信噪比核磁共振数据情况下，均能够得到稳定的反演结果。
图12为本发明一示例性实施例示出的基于L0正则化的核磁共振回波数据反演装置的结构示意图，如图12所示，本实施例的装置包括：
确定模块121，用于根据核磁共振测量参数确定核矩阵。
构建模块122，用于根据核矩阵，并以L0范数作为目标函数的惩罚项，构建目标函数。
确定模块121，还用于基于迭代算法，确定目标函数的最优解。
确定模块121，还用于基于弛豫信息，根据最优解，确定核磁共振回波数据反演的核磁共振谱。
本实施例的装置可以用于执行图1所示方法实施例的技术方案，其实现原理和技术效果类似，此处不再赘述。
图13为本发明另一示例性实施例示出的基于L0正则化的核磁共振回波数据反演装置的结构示意图，在上述实施例的基础上，
进一步地，确定模块121包括：
核矩阵确定子模块1211，用于根据核磁共振测量的脉冲序列，确定核矩阵；其中，核磁共振测量的脉冲序列为一维核磁共振CPMG脉冲序列、二维核磁共振SR-CPMG脉冲序列或基于变回波间隔TE的二维核磁共振CPMG脉冲序列；
当核磁共振测量的脉冲序列为一维核磁共振CPMG脉冲序列时，核矩阵确定子模块1211确定核矩阵为
当核磁共振测量的脉冲序列为二维核磁共振SR-CPMG脉冲序列时，核矩阵确定子模块1211确定核矩阵为
其中
表示两个矩阵的Kronecker积；
当核磁共振测量的脉冲序列为基于变回波间隔TE的二维核磁共振CPMG脉冲序列时，核矩阵确定子模块1211确定核矩阵为
可选的，构建模块122，具体用于根据核矩阵K，构建目标函数：
其中，||·||0表示向量的L0范数，s表示反演的谱，以||s||0作为目标函数的惩罚项；m表示核磁共振回波数据；其中，s和m为向量；s≥0表示向量s中的所有元素大于等于0，τ为正则化参数。
可选的，确定模块121包括：
获取子模块1212，用于获取噪声信号的标准差σ。
更新子模块1213，用于根据标准差σ，对初始正则化参数τinit进行迭代更新，得到更新后的正则化参数τ；其中，初始正则化参数
更新后的正则化参数
式中，N为核磁共振数据点的个数，g为向量s中元素的个数，snew为前一次迭代的解s。
第一确定子模块1214，用于根据预设的初始解s0，确定第一解s1，s1＝Ψτ(s0+KT(m-Ks0))；其中，
第二确定子模块1215，用于根据初始解s0、第一解s1，确定第n次迭代的解sn+1，sn+1＝asn-1+bsn+cΨτ(sn+KT(m-Ksn))；其中，n表示迭代次数，且n≥1；a、b、c为预设的迭代参数。
误差确定子模块1216，用于根据解sn+1，确定误差
判断子模块1217，用于根据最小误差tol，判断误差ζ＜tol的条件是否成立。
最优解确定子模块1218，用于当判断子模块1217判断误差ζ＜tol的条件成立时，确定目标函数的最优解为sn+1；还用于当判断子模块1217判断误差ζ＜tol的条件不成立时，更新n＝n+1，并返回第二确定子模块1215执行确定第n次迭代的解sn+1的步骤。
可选的，确定模块121，包括：
第一反演子模块1219，用于当反演一维核磁共振T2谱时，基于横向弛豫时间T2，根据最优解sn+1，确定核磁共振T2谱；
第二反演子模块1220，用于当反演二维核磁共振T2-T1谱时，对最优解sn+1进行排序，生成与最优解sn+1对应的二维矩阵S1，基于横向弛豫时间T2和纵向弛豫时间T1，根据二维矩阵S1，确定核磁共振T2-T1谱；
第三反演子模块1221，用于当反演二维核磁共振T2-D谱时，对最优解sn+1进行排序，生成与最优解sn+1对应的二维矩阵S2，基于横向弛豫时间T2和扩散系数D，根据二维矩阵S2，确定核磁共振T2-D谱。
本实施例的装置可以用于执行图2所示方法实施例的技术方案，其实现原理和技术效果类似，此处不再赘述。
本领域普通技术人员可以理解：实现上述各方法实施例的全部或部分步骤可以通过程序指令相关的硬件来完成。前述的程序可以存储于一计算机可读取存储介质中。该程序在执行时，执行包括上述各方法实施例的步骤；而前述的存储介质包括：ROM、RAM、磁碟或者光盘等各种可以存储程序代码的介质。
最后应说明的是：以上各实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述各实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的范围。
Claims (2)
1.一种基于L0正则化的核磁共振回波数据反演方法，其特征在于，包括：
根据核磁共振测量参数确定核矩阵；
根据所述核矩阵，并以L0范数作为目标函数的惩罚项，构建所述目标函数；
基于迭代算法，确定所述目标函数的最优解；
基于弛豫信息，根据所述最优解，确定核磁共振回波数据反演的核磁共振谱；
其中，
所述根据核磁共振测量参数确定核矩阵，包括：
根据核磁共振测量的脉冲序列，确定所述核矩阵；
其中，所述核磁共振测量的脉冲序列为一维核磁共振CPMG脉冲序列、二维核磁共振SR-CPMG脉冲序列或基于变回波间隔TE的二维核磁共振CPMG脉冲序列；
若所述核磁共振测量的脉冲序列为一维核磁共振CPMG脉冲序列，则所述核矩阵为
若所述核磁共振测量的脉冲序列为二维核磁共振SR-CPMG脉冲序列，则所述核矩阵为
其中
表示两个矩阵的Kronecker积；
若所述核磁共振测量的脉冲序列为基于变回波间隔TE的二维核磁共振CPMG脉冲序列，则所述核矩阵为
其中，t为进行核磁共振测量时每个回波对应的时刻；T1为纵向弛豫时间；T2为横向弛豫时间；TW为等待时间；γ为旋磁比；G为磁场梯度；TE为回波间隔；D为扩散系数；
其中，
所述根据所述核矩阵，并以L0范数作为目标函数的惩罚项，构建所述目标函数，包括：
根据所述核矩阵K，构建目标函数：
其中，||·||0表示向量的L0范数，s表示反演的谱，以||s||0作为所述目标函数的惩罚项；m表示核磁共振回波数据；其中，s和m为向量；s≥0表示向量s中的所有元素大于等于0，τ为正则化参数；
其中，
所述基于迭代算法，确定所述目标函数的最优解，包括：
获取噪声信号的标准差σ，根据所述标准差σ，对初始正则化参数τinit进行迭代更新，得到更新后的正则化参数τ；
其中，所述初始正则化参数
所述更新后的正则化参数
式中，N为核磁共振数据点的个数，g为向量s中元素的个数，snew为前一次迭代的解；
根据预设的初始解s0，确定第一解s1，s1＝Ψτ(s0+KT(m-Ks0))；
其中，
根据所述初始解s0、所述第一解s1，确定第n次迭代的解sn+1，
sn+1＝asn-1+bsn+cΨτ(sn+KT(m-Ksn))；其中，n表示迭代次数，且n≥1；a、b、c为预设的迭代参数；
根据所述解sn+1，确定误差
根据最小误差tol，判断所述误差ζ<tol的条件是否成立；
若成立，则确定所述目标函数的最优解为sn+1；
若不成立，则更新n＝n+1，并返回执行所述确定第n次迭代的解sn+1的步骤；
其中，
所述基于弛豫信息，根据所述最优解，确定核磁共振回波数据反演的核磁共振谱，包括：
若反演一维核磁共振T2谱时，基于横向弛豫时间T2，根据所述最优解sn+1，确定核磁共振T2谱；
若反演二维核磁共振T2-T1谱时，对所述最优解sn+1进行排序，生成与所述最优解sn+1对应的二维矩阵S1，基于横向弛豫时间T2和纵向弛豫时间T1，根据所述二维矩阵S1，确定核磁共振T2-T1谱；
若反演二维核磁共振T2-D谱时，对所述最优解sn+1进行排序，生成与所述最优解sn+1对应的二维矩阵S2，基于横向弛豫时间T2和扩散系数D，根据所述二维矩阵S2，确定核磁共振T2-D谱。
2.一种基于L0正则化的核磁共振回波数据反演装置，其特征在于，包括：
确定模块，用于根据核磁共振测量参数确定核矩阵；
构建模块，用于根据所述核矩阵，并以L0范数作为目标函数的惩罚项，构建所述目标函数；
所述确定模块，还用于基于迭代算法，确定所述目标函数的最优解；
所述确定模块，还用于基于弛豫信息，根据所述最优解，确定核磁共振回波数据反演的核磁共振谱；
其中，
所述确定模块包括：
核矩阵确定子模块，用于根据核磁共振测量的脉冲序列，确定所述核矩阵；其中，所述核磁共振测量的脉冲序列为一维核磁共振CPMG脉冲序列、二维核磁共振SR-CPMG脉冲序列或基于变回波间隔TE的二维核磁共振CPMG脉冲序列；
当所述核磁共振测量的脉冲序列为一维核磁共振CPMG脉冲序列时，所述核矩阵确定子模块确定所述核矩阵为
当所述核磁共振测量的脉冲序列为二维核磁共振SR-CPMG脉冲序列时，所述核矩阵确定子模块确定所述核矩阵为
当所述核磁共振测量的脉冲序列为基于变回波间隔TE的二维核磁共振CPMG脉冲序列时，所述核矩阵确定子模块确定所述核矩阵为
其中，
所述构建模块，具体用于根据所述核矩阵K，构建目标函数：
其中，||·||0表示向量的L0范数，s表示反演的谱，以||s||0作为所述目标函数的惩罚项；m表示核磁共振回波数据；其中，s和m为向量；s≥0表示向量s中的所有元素大于等于0，τ为正则化参数；
其中，
所述确定模块包括：
获取子模块，用于获取噪声信号的标准差σ；
更新子模块，用于根据所述标准差σ，对初始正则化参数τinit进行迭代更新，得到更新后的正则化参数τ；其中，所述初始正则化参数
所述更新后的正则化参数
式中，N为核磁共振数据点的个数，g为向量s中元素的个数，snew为前一次迭代的解；
第一确定子模块，用于根据预设的初始解s0，确定第一解s1，s1＝Ψτ(s0+KT(m-Ks0))；
其中，
第二确定子模块，用于根据所述初始解s0、所述第一解s1，确定第n次迭代的解sn+1，sn+1＝asn-1+bsn+cΨτ(sn+KT(m-Ksn))；其中，n表示迭代次数，且n≥1；a、b、c为预设的迭代参数；
误差确定子模块，用于根据所述解sn+1，确定误差
判断子模块，用于根据最小误差tol，判断所述误差ζ<tol的条件是否成立；
最优解确定子模块，用于当所述判断子模块判断所述误差ζ<tol的条件成立时，确定所述目标函数的最优解为sn+1；还用于当所述判断子模块判断所述误差ζ<tol的条件不成立时，更新n＝n+1，并返回所述第二确定子模块执行所述确定第n次迭代的解sn+1的步骤；
其中，
所述确定模块，包括：
第一反演子模块，用于当反演一维核磁共振T2谱时，基于横向弛豫时间T2，根据所述最优解sn+1，确定核磁共振T2谱；
第二反演子模块，用于当反演二维核磁共振T2-T1谱时，对所述最优解sn+1进行排序，生成与所述最优解sn+1对应的二维矩阵S1，基于横向弛豫时间T2和纵向弛豫时间T1，根据所述二维矩阵S1，确定核磁共振T2-T1谱；
第三反演子模块，用于当反演二维核磁共振T2-D谱时，对所述最优解sn+1进行排序，生成与所述最优解sn+1对应的二维矩阵S2，基于横向弛豫时间T2和扩散系数D，根据所述二维矩阵S2，确定核磁共振T2-D谱。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201711346882.7A CN108009125B (zh) | 2017-12-15 | 2017-12-15 | 基于l0正则化的核磁共振回波数据反演方法及装置 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201711346882.7A CN108009125B (zh) | 2017-12-15 | 2017-12-15 | 基于l0正则化的核磁共振回波数据反演方法及装置 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN108009125A CN108009125A (zh) | 2018-05-08 |
| CN108009125B true CN108009125B (zh) | 2020-05-12 |
Family
ID=62059151
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201711346882.7A Active CN108009125B (zh) | 2017-12-15 | 2017-12-15 | 基于l0正则化的核磁共振回波数据反演方法及装置 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN108009125B (zh) |
Families Citing this family (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110618463B (zh) * | 2019-09-02 | 2020-07-28 | 中国石油大学(北京) | 一种核磁共振数据反演方法、装置、存储介质及设备 |
| CN112710688B (zh) * | 2019-10-24 | 2023-08-22 | 中国石油天然气股份有限公司 | 核磁共振纵向弛豫获取方法及系统 |
| CN113743682B (zh) * | 2021-11-03 | 2022-02-18 | 中国科学院精密测量科学与技术创新研究院 | 一种基于有监督深度神经网络的nmr弛豫时间反演方法 |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104375108A (zh) * | 2014-11-19 | 2015-02-25 | 上海理工大学 | 一种基于lsqr的低场核磁共振二维谱反演算法 |
| CN106156505A (zh) * | 2016-07-05 | 2016-11-23 | 中国科学技术大学 | 一种基于正交匹配追踪算法的核磁共振t2谱反演方法 |
| CN107102020A (zh) * | 2017-03-27 | 2017-08-29 | 北京青檬艾柯科技有限公司 | 多维核磁共振测量方法 |
Family Cites Families (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US6859033B2 (en) * | 2002-08-28 | 2005-02-22 | Schlumberger Technology Corporation | Method for magnetic resonance fluid characterization |
| CN103116148B (zh) * | 2013-01-30 | 2015-04-01 | 上海理工大学 | 一种核磁共振二维谱反演的方法 |
| CN105891249B (zh) * | 2016-03-29 | 2018-03-13 | 上海理工大学 | 时域核磁共振谱反演的方法 |
| CN106201991B (zh) * | 2016-07-07 | 2019-03-01 | 中国石油大学(北京) | 核磁共振回波数据反演方法及装置 |
| CN106199474B (zh) * | 2016-07-21 | 2018-10-12 | 上海理工大学 | 一种低场核磁共振二维谱反演算法 |
-
2017
- 2017-12-15 CN CN201711346882.7A patent/CN108009125B/zh active Active
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104375108A (zh) * | 2014-11-19 | 2015-02-25 | 上海理工大学 | 一种基于lsqr的低场核磁共振二维谱反演算法 |
| CN106156505A (zh) * | 2016-07-05 | 2016-11-23 | 中国科学技术大学 | 一种基于正交匹配追踪算法的核磁共振t2谱反演方法 |
| CN107102020A (zh) * | 2017-03-27 | 2017-08-29 | 北京青檬艾柯科技有限公司 | 多维核磁共振测量方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN108009125A (zh) | 2018-05-08 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN110431408B (zh) | 使用短rf线圈的空间灵敏度曲线的长全岩心岩样的高空间分辨率核磁共振 | |
| CN110785682B (zh) | 对井下多维测量结果的快速测量和解释 | |
| CN109270107B (zh) | 多维核磁共振测量方法 | |
| CN108009125B (zh) | 基于l0正则化的核磁共振回波数据反演方法及装置 | |
| US9222902B2 (en) | Estimations of nuclear magnetic resonance measurement distributions | |
| CN107748126B (zh) | 一种获取岩石孔隙尺寸和孔隙表面弛豫率的核磁共振方法 | |
| US10732314B2 (en) | Estimation of petrophysical and fluid properties using integral transforms in nuclear magnetic resonance | |
| CN111596366B (zh) | 一种基于地震信号优化处理的波阻抗反演方法 | |
| CN108027409A (zh) | 时域mri | |
| US10228484B2 (en) | Robust multi-dimensional inversion from wellbore NMR measurements | |
| CN111898734B (zh) | 一种基于mlp的nmr弛豫时间反演方法 | |
| Jin et al. | A new method for permeability estimation using integral transforms based on NMR echo data in tight sandstone | |
| Parker et al. | Assigning uncertainties in the inversion of NMR relaxation data | |
| Trevizan et al. | Method for predicting permeability of complex carbonate reservoirs using NMR logging measurements | |
| CN108663711B (zh) | 一种基于τ分布的贝叶斯地震反演方法 | |
| Gu et al. | A novel method for NMR data denoising based on discrete cosine transform and variable length windows | |
| US9874619B2 (en) | Methods for performing NMR measurements on porous media | |
| Kwak et al. | Predicting Carbonate Rock Properties Using NMR Data and Generalized Interpolation-Based Techniques | |
| CN110320227B (zh) | 一种二维核磁共振d-t2谱反演方法及装置 | |
| Gruber et al. | A more accurate estimate of T2 distribution from direct analysis of NMR measurements | |
| CN115356672B (zh) | 多维磁共振成像方法、系统及存储介质 | |
| CN105759233B (zh) | 一种快速化学交换饱和转移成像方法和系统 | |
| CN107727678B (zh) | 一种核磁共振弛豫高低本征模态耦合方法 | |
| CN113743682B (zh) | 一种基于有监督深度神经网络的nmr弛豫时间反演方法 | |
| CN107861918B (zh) | 基于m-稀疏算法的核磁共振回波数据反演方法与装置 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |