一种核磁共振二维谱反演的方法转让专利
申请号 : CN201310035140.8
文献号 : CN103116148B
文献日 : 2015-04-01
发明人 : 周小龙 , 聂生东 , 王远军 , 张英力 , 杨培强
申请人 : 上海理工大学
摘要 :
本发明涉及一种核磁共振二维谱反演的方法，包括第一步，噪声提取与估计：使用小波变换对采集数据CPMG回波串中的噪声进行提取并估计其标准差；第二步，数据压缩：生成反演核，并利用核矩阵的秩进行截断奇异值分解与重构来完成数据压缩；第三步，数据拟合：对压缩后数据的拟合问题进行正则化处理，并用结合了非精确一维搜索的牛顿法进行正则化因子和反演谱的迭代求解，得到反演谱。大幅提高了二维反演算法的执行效率与二维谱的分辨率，同时还具有很好的鲁棒性。
权利要求 :
1.一种核磁共振二维谱反演的方法，其特征在于，包括如下步骤：
1）噪声提取与估计：使用小波变换对采集数据CPMG回波串中的噪声进行提取并估计其标准差；
2）数据压缩：对采集数据CPMG回波串，利用核矩阵的秩对核矩阵进行截断奇异值分解与重构，实现数据压缩；
3）数据拟合：对压缩后数据的拟合问题进行正则化处理，并用结合了非精确一维搜索的牛顿法进行正则化因子和反演谱的迭代求解，得到反演谱。
2.根据权利要求1所述核磁共振二维谱反演的方法，其特征在于，所述步骤1）中噪声提取与估计：在噪声提取中引入了小波变换，使用小波去噪对CPMG回波串进行去噪，得到去噪后数据；通过比较去噪后数据与去噪前数据的幅度差异，进行“噪声”提取再进行噪声标准差的估计；具体步骤为：首先，使用sym8小波基对原始数据进行4级分解，然后使用全局阈值对小波系数进行软阈值化，最后进行信号重构，得到滤波后相对光滑的CPMG数据串；然后计算原始数据与滤波后数据的差值，使用高斯白噪声的幅度模型对差值数据进行拟合，得到对噪声标准差的估计。
3.根据权利要求1所述核磁共振二维谱反演的方法，其特征在于，所述步骤3）数据拟合首先对核矩阵各元素平方后得到一个新矩阵，将新矩阵的迹(trace)除以对角元素个数的值作为初始的正则化因子；其次对于初始的正则因子，使用牛顿法对极小值问题进行求解，在牛顿方向上使用的步长由非精确一维搜索确定。
4.根据权利要求1或3所述核磁共振二维谱反演的方法，其特征在于，所述步骤3）正则化因子的迭代求解：根据离差原理计算正则化因子的理想值，如果理想值大于原来的值，那么将原值的1/2作为新的正则化因子的值；否则，直接使用理想值作为正则化因子的新值。
5.根据权利要求1所述核磁共振二维谱反演的方法，其特征在于，所述步骤3）反演谱的迭代求解的终值条件为下面条件之一：①正则化因子下一轮循环的新值已经不大于一个预设的小值；
正则化因子在连续两次迭代的过程中变化幅度小于预设阈值；
拟合残差在连续两次迭代的过程中变化幅度小于预设阈值；
拟合残差不大于步骤1）噪声提取与估计中估计标准差。
说明书 :
一种核磁共振二维谱反演的方法
技术领域
[0001] 本发明涉及一种核磁共振领域中信号处理，特别涉及一种核磁共振二维谱反演的方法。
背景技术
[0002] 核磁共振技术在能源勘探领域已经得到了广泛应用。NMR测井能够快速、无损地为储层流体识别、岩石物性评价和产能评估等常规测井工作流中的各个步骤提供精准的信息。相比电法测井、声波测井等传统方法，NMR测井法是惟一一种不使流体流动就能对流体类型进行识别的方法。基于一维测井实验的解决方案在效率和精度上都存在明显的不足，二维NMR测井技术应运而生。使用D-T2（扩散-横向弛豫）二维谱可以快速、直观地对水、油等不同成分进行区分，通过计算T1/T2（纵向弛豫/横向弛豫）这一比值还能准确地对烃的类型进行判断，二维NMR测井方法在定性和定量分析中都有着得天独厚的优势。二维NMR检测方法的实现还拓宽了常规NMR检测实验的应用范围，为食品、农业、生物材料等领域提供了更精确、可靠的解决方案。
[0003] NMR检测方法采集的数据并不能直接使用，真正需要的“谱”信息需要通过这些实验数据进行反演后才能得到。目前世界上只有极少数NMR测井设备制造商有能力生产二维测井设备，而且这些国际巨头只向各国能源公司提供仪器租赁和有偿的解释服务，并不出售设备和反演软件，二维反演算法几近处于被垄断状态。另一方面，二维实验采集的数据量远大于一维实验，传统的一维反演方法已经不能满足二维反演要求。快速Laplace Inversion方法是一种基于BRD方法和离差原理的二维反演方法，这种方法致力于将拟合误差限定在估计的数据误差水平（离差），能够从低SNR数据中反演出合理的结果。但是，如果数据SNR较高，在不添加其它限制条件的情况下，这种方法可能无法收敛。国内研究者使用最广泛的是基于截断奇异值分解（Truncated Singular Value Decomposition, TSVD）的反演方法，TSVD方法在SNR很高时能够得到令人满意的结果。但是由于TSVD对噪声非常敏感，当SNR较低时使用该方法无法得到准确的结果。最大熵法能够得到最随机（最概然）的二维谱，但是以熵(对数项)作为罚函数的方法在求解过程中运算量过大，在普通PC机上执行速度过慢。通过引入基函数可以减少反演的计算量，但是使用基函数法进行处理的前提是，已知待求谱中值的分布满足某种线形。其它如蒙特卡洛类方法等使用概率来控制进化方向的方法，在二维反演这一大数据量问题中仅仅理论可行，大规模的种群会导致进化速度过慢。
发明内容
[0004] 本发明是针对现在二维反演存在的问题，提出了一种核磁共振二维谱反演的方法，可在保证反演谱准确性的前提下，提高反演谱的分辨率、反演算法的效率和鲁棒性。
[0005] 本发明的技术方案为：一种核磁共振二维谱反演的方法，包括如下步骤：
[0006] 1）噪声提取与估计：使用小波变换对采集数据CPMG回波串中的噪声进行提取并估计其标准差；
[0007] 2）数据压缩：生成反演核，并利用核矩阵的秩进行截断奇异值分解与重构来完成数据压缩；
[0008] 3）数据拟合：对压缩后数据的拟合问题进行正则化处理，并用结合了非精确一维搜索的牛顿法进行正则化因子和反演谱的迭代求解，得到反演谱。
[0009] 所述步骤1）中使用小波去噪的方法对CPMG数据串进行去噪，通过比较去噪后数据与滤波前数据的幅度差异，进行“噪声”提取，首先，使用sym8小波基对原始数据进行4级分解，然后使用全局阈值对小波系数进行软阈值化，最后进行信号重构，得到滤波后相对光滑的CPMG数据串；然后计算原始数据与滤波后数据的差值，使用高斯白噪声的幅度模型对差值数据进行拟合，得到对噪声标准差的估计。
[0010] 所述步骤3）数据拟合首先对核心矩阵各元素平方后得到一个新矩阵，将新矩阵的迹(trace)除以对角元素个数的值作为初始的正则化因子；其次对于给定的正则因子，使用牛顿法（Newton Method）对极小值问题进行求解，在牛顿方向上使用的步长由Wolfe-Powell非精确一维搜索确定。
[0011] 所述步骤3）正则化因子更新方式为：根据离差原理计算正则化因子的理想值，如果理想值大于原来的值，那么将原值的1/2作为新的正则化因子的值；否则，直接使用理想值作为正则化因子的新值。
[0012] 所述步骤3）
[0013] 反演谱的迭代求解的终值条件为：
[0014] ①正则化因子下一轮循环的新值已经不大于一个预设的小值；
[0015] 正则化因子在连续两次迭代的过程中变化幅度小于预设阈值；
[0016] 拟合残差在连续两次迭代的过程中变化幅度小于预设阈值；
[0017] 拟合残差不大于第一步估计的噪声水平。
[0018] 本发明的有益效果在于：本发明核磁共振二维谱反演的方法，大幅提高了二维反演算法的执行效率与二维谱的分辨率，同时还具有很好的鲁棒性。
附图说明
[0019] 图1为本发明核磁共振二维谱反演的方法流程图；
[0020] 图2为本发明核磁共振二维谱反演的方法中噪声提取与估计部分的算法流程图；
[0021] 图3为本发明核磁共振二维谱反演的方法中对一条CPMG数据串的噪声进行提取的结果图；
[0022] 图4为本发明核磁共振二维谱反演的方法中数据拟合部分使用的停机准则示意图；
[0023] 图5为本发明核磁共振二维谱反演的方法为某品牌棉籽的T1-T2谱进行处理的二维谱图；
[0024] 图6为本发明核磁共振二维谱反演的方法为某饱水岩心的D-T2谱进行处理的二维谱图；
[0025] 图7为本发明核磁共振二维谱反演的方法预设真实谱图；
[0026] 图8为本发明核磁共振二维谱反演的方法对信噪比为100的仿真数据进行处理的结果图；
[0027] 图9为本发明核磁共振二维谱反演的方法对信噪比为10的仿真数据进行处理的结果图；
[0028] 图10为本发明核磁共振二维谱反演的方法对信噪比为1的仿真数据进行处理的结果图。
具体实施方式
[0029] 图1是本发明核磁共振二维谱反演的方法流程图，包括步骤：第一步，噪声提取与估计：使用小波变换对采集数据CPMG回波串中的噪声进行提取并估计其标准差；第二步，数据压缩：生成反演核，并利用核矩阵的秩进行截断奇异值分解与重构来完成数据压缩；第三步，数据拟合：对压缩后数据的拟合问题进行正则化处理，并用结合了非精确一维搜索的牛顿法进行正则化因子和反演谱的迭代求解，得到反演谱。
[0030] 1、噪声提取与估计：
[0031] 如图2是噪声提取与估计部分的算法流程图。在噪声提取中引入了小波变换，使用小波去噪的方法对CPMG数据串进行去噪，通过比较去噪后数据与滤波前数据的幅度差异，进行“噪声”提取。本发明使用的小波去噪的方式为：首先，使用sym8小波基对原始数据进行4级分解，然后使用全局阈值对小波系数进行软阈值化，最后进行信号重构。这里的“噪声”是指影响信号光滑性的成分，是出于曲线光滑性的考虑，这正是后续处理中没有直接使用滤波后数据的主要原因。为了对噪声的方差进行估计，需要对噪声幅度进行直方图统计，然后使用式（1）对直方图进行拟合。（1）式中，A和 为拟合参数， 即为高斯噪声的方差。图3为对一条CPMG数据串的噪声进行提取的结果图。
[0032] （1）
[0033] 2、数据压缩：
[0034] 对核矩阵K，总存在如下分解：
[0035] （2）
[0036] 式中， 为 的实二维矩阵； 为 的列正交矩阵，称为左奇异矩阵； 为的列正交矩阵，称为右奇异矩阵； 为 的对角矩阵，对角位置的元素值对应着矩阵 的奇异值。由矩阵SVD计算方法可以使S中的奇异值从左上角到右下角递减。在很多情况下，前10%甚至1%的奇异值的和就占了全部的奇异值之和的99%以上了。也就是说，我们也可以用前 大的奇异值来近似描述矩阵：
[0037] （3）
[0038] 进一步对二维反演问题进行分析，如果反演核能够拆分成两个独立核（如T1-T2反演核），那么正演规律可以表示为：
[0039] （4）
[0040] 式中，M为实验测量的数据， 表示测量数据与待求谱之间的已知关系。
[0041] 现对两个核分别进行奇异值截断： ，（ 表示截断位置），那么
[0042] （5）
[0043] 其中, ， 。显然，如果 ，那么成立，采样数据 将大幅减少，而在这个过程中，待求谱 的大小没有
改变。（后续步骤中的处理对象均为使用此方法进行压缩后的数据，为了便于表述，将不再使用符号“~”对压缩后数据进行区分）数据约简以损失极小的精度为代价，能够加速后续处理步骤中的各种运算，提高算法效率。
改变。（后续步骤中的处理对象均为使用此方法进行压缩后的数据，为了便于表述，将不再使用符号“~”对压缩后数据进行区分）数据约简以损失极小的精度为代价，能够加速后续处理步骤中的各种运算，提高算法效率。
[0044] 关于截断位置，传统方法一般使用SNR或者预设条件数进行判断，这些截断方式很可能导致矩阵损失一些重要信息，后续反演得到的结果也会是欠优的（suboptimal）。实际上，如果两个核的奇异值矩阵分别为：
[0045]
[0046] 那么， 的奇异值为 ，显然，要使截断后的数据保留所有原始的特征（即所有大于0的奇异值），必须保证截断位置 和 刚好等于矩阵的秩。
[0047] 3、数据拟合：
[0048] 二维反演问题实质上是一个拟合问题（如式6所示）。如果两个核没有耦合关系，那么可以直接用张量积将问题（6）转化为问题（7）的形式。此时， ,， ，vect表示将二维矩阵按列拼接为一个列向量。如果两个核耦合在一起，那么可以按一定的方式将两个核重新编纂到一个核中，也能改为问题（7）的形式。
[0049] （6）
[0050] （7）
[0051] 二维反演问题是一个不适定问题，在求解时需要使用一定的附加约束以得到惟一的稳定解。在众多的附加约束中，对解的模或光滑性进行限制符合人们的认知，因而被广泛使用。以模平滑项作为惩罚函数的标准Tikhonov正则化的方法如式（8）所示。
[0052] （8）
[0053] 式中， 是一个正数，表征了模平滑在整个最小值问题中占据的权重，称为正则因子。其余参数意义同前（未对压缩前后的数据进行区分），2的添加只是为了方便进行后续的数学运算。若设
[0054] （9）
[0055] （10）
[0056] （11）
[0057] 由KKT条件很容易将问题转化为式（12）所述形式。
[0058] （12）
[0059] 对于最小值问题（12），由于 至少半正定， 恒成立，所以二次规划(12)是一个凸规划，具有惟一解，且 可以取任意值。又由于存在
[0060] （13）
[0061] 对于给定的α，可以直接用牛顿法进行求解得到 。然而，实际问题中所处理的数据范围很广，直接使用牛顿法处理部分数据时可能会产生震荡，无法收敛。为了解决这个问题，同时加快算法的收敛速度，本发明引入了Powell Wolfe准则这一非精确一维搜索方法来确定每一次迭代在牛顿方向上行进的步长。
[0062] 上述处理步骤都需要在已知α的前提下进行，但是在实际反演问题中，α和 都是未知的。α可以通过L曲线或者GCV来确定，但是这两种方法都需要对某个预设范围内的不同值进行遍历，效率存在明显不足。而且，实际问题中，很多数据按照L曲线绘制的曲线并非严格的“L”形，因而无法对拐角进行定位。GCV得到的曲线在根部往往过于平坦，迭代的收敛速度过慢。BRD法是一种效率更快的方法，通过迭代的方式对 进行更新时，最大限度地减少了对一定范围内多余 值的计算。但是，如果对噪声的估计不准确，或者数据的信噪比过高时，BRD法会无法收敛。
[0063] 本发明从效率出发，以BRD方法为原型，提出了一种全新的迭代方式。结合前述步骤，本发明数据拟合部分的步骤如下：
[0064] （1）给定α一个偏大的初始值（如 ）；设定最小正则化因子Thresh和最小变化率TOL；
[0065] （2）根据 求取 ，使用 。如果 ，使用 进行更新；如果 ， ；
[0066] （3）若 ，或 ，或，或 不再变化 ，停机；否则，转至步骤（2）。
[0067] 图4是数据拟合部分使用的停机准则示意图，横坐标为不同的 ，纵坐标表示不同 对应的拟合残差。由于事先给定的 是一个很大的值（按照 计算的话3
通常不小于10量级），此时得到的 不可能是最优解。（如果一组数据的最佳正则因子很大，那么根据这组数据得到的反演谱已经没有多少真实的信息和实际应用价值。）如果噪声估计准确，按照式 将会得到一个变小的正则因子，通常经过10步以下的迭代之后 就不再变化，得到一个刚好使拟合残差等于噪声水平的解。如果噪声估计过小，仍然按照上一种情况的方法进行更新时， 会不断减小，并逐渐失去正则化的作用，最终产生一个欠光滑的不稳定解。此时，为了确保最终得到的是一个稳定的光滑解，本文算法设置了终值条件 ，和 ，通常，终止条件 更容易达到。如果噪声估计
过大，通过 来更新正则因子时 会来回震荡（因为很容易得到拟合残差小于数据噪声的解，二者刚好相等的概率很小），显然，终止条件 能够解决这个问题，使拟合误差在第一次小于噪声水平时就停机。
通常不小于10量级），此时得到的 不可能是最优解。（如果一组数据的最佳正则因子很大，那么根据这组数据得到的反演谱已经没有多少真实的信息和实际应用价值。）如果噪声估计准确，按照式 将会得到一个变小的正则因子，通常经过10步以下的迭代之后 就不再变化，得到一个刚好使拟合残差等于噪声水平的解。如果噪声估计过小，仍然按照上一种情况的方法进行更新时， 会不断减小，并逐渐失去正则化的作用，最终产生一个欠光滑的不稳定解。此时，为了确保最终得到的是一个稳定的光滑解，本文算法设置了终值条件 ，和 ，通常，终止条件 更容易达到。如果噪声估计
过大，通过 来更新正则因子时 会来回震荡（因为很容易得到拟合残差小于数据噪声的解，二者刚好相等的概率很小），显然，终止条件 能够解决这个问题，使拟合误差在第一次小于噪声水平时就停机。
[0068] 使用上述停机准则的迭代算法总能得到一个合适的反演谱，具有一定的通用性。图5为使用本发明所提出的方法进行处理得到的某品牌棉籽的T1-T2二维谱，图6为使用本发明所提出的方法进行处理得到的某饱水岩心的D-T2二维谱。经试验表明，采用本发明一种核磁共振二维谱反演的新方法对不同领域的实验数据进行反演时，经过不到10次的迭代之后算法就能收敛。如图7～10在对不同信噪比的仿真数据进行反演时，能够对SNR=1的数据给出清晰的二维谱。
[0069] 与现有技术相比，本发明一、利用小波滤波进行噪声提取的方法能够更准确地对噪声进行估计。使用高斯白噪声的函数模型对差值进行拟合时，相关系数大于99%，说明文献中将CPMG回波信号中的噪声假设为高斯白噪声是合理的，也从侧面证明了本发明中噪声提取部分的精确性；二、使用的数据压缩方法中奇异值的截断方式，是惟一一种在对数据进行大幅度压缩的同时，能够保证新问题的解与原问题等价的方法。其它根据预设条件数或者信噪比来进行截断的方式都可能造成新问题的最优解是原问题的欠优解。此外，本发明使用的数据压缩方式无需认为干预，可以自适应的得到截断位置；三、即使布点数在100*100左右，也能在数分钟甚至数十秒之内得到反演谱(以普通32位机为实验平台，具体时间依据具体数据而定)。布点数得到很大的提高，二维谱的分辨率也更为精细；四、对于信噪比较高的数据，算法可以得到十分准确的二维谱，可以用于定量分析。对于信噪比较低的数据，算法能够很快的收敛并得到一个可靠的二维谱，亦可用于定量分析。对于信噪比极低的数据(SNR<1)的数据，算法能更快的收敛，反演谱中可以清晰的看到主要成分的峰值分布，但是一些小比例的峰或者相隔很近的峰可能被大峰淹没，反演谱中峰的半高宽也较宽，此时可以用来进行定性分析。