一种基于lsqr的低场核磁共振二维谱反演算法
【专利摘要】本发明涉及一种基于LSQR(LeastSquaresQRDecomposition，最小二乘QR分解)的低场核磁共振二维谱反演算法，属于核磁共振信号处理领域。本发明的优点表现为：本发明与目前国内外文献报道的其他低场核磁共振反演算法相比，其借助“L”曲线自适应选择合适的迭代次数，得到的解更能体现真实的谱分布；通过迭代修正来实现非负限制，计算精度高，鲁棒性好，在不同信噪比数据情况下，能够得到稳定的反演结果。
【专利说明】-种基于LSQR的低场核磁共振二维谱反演算法
【技术领域】
[0001] 本发明涉及核磁共振信号处理领域，涉及一种低场核磁共振二维谱反演算法，具 体地说，是一种基于LSQR的低场核磁共振二维谱反演算法。
【背景技术】
[0002] 低场核磁共振设备凭借成本相对低廉、易于小型化，以及工业集成等优势，已经在 能源、食品、农业、医学等领域得到广泛应用。低场核磁共振分析技术是这些领域中一种重 要的检测手段。核磁共振采集到的原始信号中包含有丰富但难以理解的样品结构信息，需 要通过反演才能得到提供反应直观信息的时域谱。其中，CPMG (Carr-Purcell-Meiboom -Gill)序列速度快，是低场核磁共振中最常用的序列之一。研究者常常利用CPMG原始数据 和样品横向弛豫时间、纵向弛豫时间的分布特点，进行时域谱反演的相关研究。时域谱可分 为一维谱和二维谱，一维谱可以为研究者分析样品的组成、性质等提供重要的依据。但是， 随着低场核磁共振应用的不断深入，研究者们发现实验得到的一维谱可能存在峰的交叠。 解决这一问题通常需要进行多次额外实验或者引入其它辅助手段，这无疑会增加实验的复 杂度，而且大多数辅助手段都只有在满足特定条件时才能使用。一维谱的交叠问题已经影 响到低场核磁共振分析技术应用领域的继续拓展，二维谱技术应运而生。二维谱不仅简化 了原有流程，而且从二维的角度进行分析，具有得天独厚的优势，可以为分析样品提供更多 有价值的信息。
[0003] 二维谱反演需要处理的数据量巨大，这给二维反演带来了极大地挑战，加上国外 相关企业的技术垄断，只提供仪器租赁和有偿的解释服务，国内二维反演技术的发展受到 极大的限制。目前，国内研究者大都使用基于截断奇异值分解（Truncated singular value decomposition, TSVD)的二维反演算法，并取得了很多成果。顾兆斌等人在国内首次实现 了 TSVD的二维反演算法，在信噪比（Signal to Noise Ratio，SNR)较高的情况下能够得 到可靠的反演结果，其后很多学者在此基础上进行了大量改进，如谭茂金等在2013年在论 文《A new inversion method for (T2, d) 2D NMR logging and fluid typing》提出的 TSVD与LSQR混合反演算法，周小龙等2013年发表的论文《基于TSVD的NMR二维谱反演算 法》中提出的逐步求精的TSVD迭代算法等。国外研究者则以显式增加罚函数的正则化算法 为主，对解的某些性质进行限制，以得到合理的反演结果。一般采用的是标准Tikhonov正 则化形式，惩罚项用于限制解的范数，这种方式的问题求解策略也相对成熟，如Honerkamp J等在 1990年发表的论文〈〈Tikhonovs regularization method for ill-posed problems: A comparison of different methods for the determination of the regularization parameter》中提出的标准Tikhonov求解算法，Yi-Qiao Song等在2013年发表的论文 《Magnetic resonance of porous media (MR-PM): A perspective》中提出的基于 BRD 算法求解Tikhonov正则化问题的二维反演算法等。此外还有对斜率和曲率进行平滑的 正则化形式，如Provencher S W等在1982年发表的论文《C0NTIN: A general purpose constrained regularization program for inverting noisy linear algebraic and integral equations》中提出的 CONTIN 算法，Ricard C A 等在 2012 年的著作《Parameter Estimation and Inverse Problems (Second Edition)》中提出的高阶Tikhonov正则化的 反演算法等。这些增加罚函数的正则化算法的关键在于如何选择一个合适的正则化因子， 现有正则化因子的选择算法不能适用于所有场合，因而需要一定的人为干预。
[0004] LSQR(Least Squares QR，最小二乘QR分解）是求解最小二乘问题的一种有效算 法，已经在地震层析成像等领域的反问题求解中得到了广泛应用。LSQR迭代不需要预先形 成TSVD法中的子空间，避免了由截断引起的虚假轮廓的出现，因而比TSVD法更具鲁棒性。 同时，LSQR直接对原始问题进行求解，无需选择正则化因子。但是，与大部分迭代反演算法 类似，具有实际应用价值的解往往出现在迭代的早期，之后的迭代虽然可以进一步减小拟 合误差，解的结构却会逐渐背离真实的谱分布。
【发明内容】
[0005] 本发明的目的是针对现有技术中的不足，根据CPMG原始数据，提出一种基于LSQR 的低场核磁共振二维谱反演算法，该算法能够提高反演运算的计算精度和鲁棒性，以得到 清晰、准确的二维谱。
[0006] 为实现上述目的，本发明采取的技术方案是： 一种基于LSQR的低场核磁共振二维谱反演算法，包括如下步骤： a. 读取低场核磁共振设备采集得到的CPMG原始数据文件； b. 对原始数据进行预处理操作得到反演核心矩阵I和信号幅值《 ; c. 根据核心矩阵I的线性关系预设一个LSQR迭代的最大迭代次数； d. 进行LSQR迭代，得到每次迭代后的反演结果和残差的范数； e. 根据不同迭代次数对应的反演结果的范数和残差的范数，绘出"L形"曲线； f. 计算"L形"曲线曲率，定义拐角位置为LSQR迭代的最优解位置； g. 使用迭代修正的方式进行约束，抑制负峰的产生，输出二维谱。
[0007] 步骤c中所述的最大迭代次数是通过大量仿真实验发现的，可将其设置为核心矩 阵的秩。
[0008] 步骤d中所述的反演结果是反映真实谱分布的一维数组。
[0009] 步骤e中所述的"L形"曲线形成原理为：随着迭代的进行，拟合误差会不断减小 并向预设阈值靠近，所以迭代次数越大，拟合误差越小；由于采样噪声无法避免，并且问题 的解不连续依赖于观测数据的变化，迭代次数越大，解的复杂度也越大，如果在直角坐标系 中以反演结果的范数和残差的范数为坐标绘图，就会得到一个"L形"的曲线。
[0010] 步骤g中所述的迭代修正的方式包括如下步骤： 5. 1.设置最小拟合误差A检测反演结果s中有没有负项，有负项转到步骤5. 2,无负项 转到步骤5. 3 ; 5. 2.将反演结果s中的负项置零，得到?计算残差A +，判 断Il Affill/lkllM是否成立，是转到步骤5.3,否则将误差分摊到&转到步骤5.1; 5. 3.得到符合真实谱分布的&输出二维谱。
[0011] 本发明优点在于： 本发明与目前国内外文献报道的其他低场核磁共振反演算法相比，其借助"L"曲线自 适应选择合适的迭代次数，得到的解更能体现真实的谱分布；通过迭代修正来实现非负限 制，计算精度高；鲁棒性好，在不同信噪比数据情况下，能够得到稳定的反演结果。
【专利附图】
【附图说明】
[0012] 图1是本发明的主要操作过程示意图。
[0013] 图2是本发明"L"形曲线示意图。
[0014] 图3是本发明迭代修正步骤流程图。
[0015] 图4(a)-图4(f)是本发明信噪比为200的仿真实验结果示意图：其中图4(a)是 构造高斯峰二维谱图；图4(b)是CPMG数据串；图4(c)是本发明算法的反演结果；图4(d) 是本发明算法的"L"曲线定位结果；图4 (e)是LSQR与TSVD混合算法反演结果；图4 (f)是 直接LSQR法反演结果。
[0016] 图5(a)-图5(f)是本发明信噪比为10的仿真实验结果示意图：其中图5(a)是 构造高斯峰二维谱图；图5(b)是CPMG数据串；图5(c)是本发明算法的反演结果；图5(d) 是本发明算法的"L"曲线定位结果；图5 (e)是LSQR与TSVD混合算法反演结果；图5 (f)是 直接LSQR法反演结果。
[0017] 图6 (a)-图6(d)是本发明实验案例结果示意图：其中图6(a)是低浓度溶液反演 结果；图6(b)是低浓度溶液"L"曲线定位图；图6(c)是加入高浓度溶液反演结果；图6(d) 是加入高浓度溶液"L"曲线定位图。
【具体实施方式】
[0018] 下面结合附图对本发明提供的【具体实施方式】作详细说明。
[0019] 如附图1所示，一种基于LSQR的低场核磁共振二维谱反演算法，包括如下步骤： a.读取低场核磁共振设备采集得到的CPMG原始数据文件。
[0020] 读取低场核磁共振设备采集得到的CPMG原始数据文件，提取出数据文件中包含 的采样时间和特定时刻的采样数据#等信息。
[0021] b.对原始数据进行横向弛豫时间、纵向弛豫时间布点，矩阵的拼接、张量积等预处 理操作得到反演核心矩阵I和信号幅值?。
[0022] 二维反演问题就是求解如公式（1)所示的具有两个核的Fredholm积分方程，r2 表示各回波波峰的回波时刻，^表示等待时间，表示横向弛豫时间为T12、纵向弛 豫时间为的物质的含量，#表示特定时刻的采样数据：
【权利要求】
1. 一种基于LSQR的低场核磁共振二维谱反演算法，其特征在于，包括如下步骤： a. 读取低场核磁共振设备采集得到的CPMG原始数据文件； b. 对原始数据进行预处理操作得到反演核心矩阵I和信号幅值《 ; c. 根据核心矩阵I的线性关系预设一个LSQR迭代的最大迭代次数； d. 进行LSQR迭代，得到每次迭代后的反演结果和残差的范数； e. 根据不同迭代次数对应的反演结果的范数和残差的范数，绘出"L形"曲线； f. 计算"L形"曲线曲率，定义拐角位置为LSQR迭代的最优解位置； g. 使用迭代修正的方式进行约束，抑制负峰的产生，输出二维谱。
2. 根据权利要求1所述的基于LSQR的低场核磁共振二维谱反演算法，其特征在于，步 骤c中所述的最大迭代次数是通过大量仿真实验发现的，可将其设置为核心矩阵的秩。
3. 根据权利要求1所述的基于LSQR的低场核磁共振二维谱反演算法，其特征在于，步 骤d中所述的反演结果是反映真实谱分布的一维数组。
4. 根据权利要求1所述的基于LSQR的低场核磁共振二维谱反演算法，其特征在于，步 骤e中所述的"L形"曲线形成原理为：随着迭代的进行，拟合误差会不断减小并向预设阈 值靠近，所以迭代次数越大，拟合误差越小；由于采样噪声无法避免，并且问题的解不连续 依赖于观测数据的变化，迭代次数越大，解的复杂度也越大，如果在直角坐标系中以反演结 果的范数和残差的范数为坐标绘图，就会得到一个"L形"的曲线。
5. 根据权利要求1所述的基于LSQR的低场核磁共振二维谱反演算法，其特征在于，步 骤g中所述的迭代修正的方式包括如下步骤： 5. 1.设置最小拟合误差A检测反演结果s中有没有负项，有负项转到步骤5. 2,无负项 转到步骤5. 3 ; 5. 2.将反演结果s中的负项置零，得到?计算残差A +，判 断II A? 11/II? II 是否成立，是则转到步骤5. 3,否则将误差分摊到&转到步骤5.1; 5. 3.得到符合真实谱分布的s，输出二维谱。
【文档编号】G01R33/56GK104375108SQ201410661016
【公开日】2015年2月25日 申请日期:2014年11月19日 优先权日:2014年11月19日
【发明者】苏冠群, 聂生东, 王丽嘉, 王远军, 周小龙, 赵彬, 张英力, 杨培强 申请人:上海理工大学