CN111239837A - 基于mcmc的地面核磁共振信号参数提取方法 - Google Patents

Abstract

本发明属于核磁共振数据处理领域，为一种基于MCMC的地面核磁共振信号参数提取方法，该方法包括，采集地面核磁共振数据；根据录入的数据通过已有方法处理尖峰噪声与工频噪声；使用希尔伯特变换和低通滤波处理的方法得到地面磁共振信号的复包络表达式；利用最大似然估计的方法提取参数值；以似然估计参数值为中心建立各参数的高斯分布作为马尔科夫链蒙特卡洛方法的先验分布；将参数分布带入到MCMC方法中通过迭代计算得到各参数的后验分布，其中出现次数最多，权值最大的参数值即为各参数的最优参数值，本发明方法可以加性随机噪声和乘性随机噪声同时存在时，实现地面磁共振数据中参数提取，且准确性好，提高地面磁共振数据解释和反演的可信度。

Description

基于MCMC的地面核磁共振信号参数提取方法

技术领域

本发明属于核磁共振数据处理领域，具体地来讲为一种基于MCMC的地面核磁共振信号参数提取方法。

背景技术

磁共振探测技术(MRS)以其无损、定量、直接等优势，被着重应用于地下水调查的领域。但是由于获得的信号十分微弱，易受到各种环境噪声的干扰，低信噪比仍然是SNMR通常面临的最大挑战。在实际应用中，尖峰噪声、工频谐波噪声、随机噪声等会严重影响MRS信号的质量，目前已有成熟的算法解决工频噪声与尖峰噪声的干扰。但是随机噪声水平过大，导致MRS信号淹没在随机噪声中，严重影响信号的特征参数提取和反演解释，因此研究强高斯噪声中的MRS特征信号提取具有十分重要的意义。

专利CN107957566B公开了一种基于频率选择奇异谱分析的磁共振测深信号提取方法，该方法基于频率选择的奇异谱分析，通过包括嵌入、RSVD分解、根据MRS信号幅值选择对应的奇异值进行矩阵重构和对角平均化四个步骤。实现MRS信号的提取，能够解决含噪MRS信号中随机噪声、尖峰噪声和工频谐波干扰的有效滤除，但是未考虑信号中的乘性噪声。

专利CN106772646A公开了一种基于工频谐波建模和自相关的地面核磁共振信号提取方法，该发明专利通过自相关和叠加处理减小随机噪声，但是由于乘性噪声的噪声水平为指数模型，不具有加性噪声的补偿性，因此通过自相关和叠加方法无法有效抑制。

综上，可以得出，上述专利均针对地面核磁共振信号中的加性噪声，未涉及信号中的乘性噪声。

发明内容

本发明所要解决的技术问题在于提供一种基于MCMC(Markov Chain MonteCarlo，马尔科夫链蒙特卡洛)的地面核磁共振信号参数提取方法，能够在加性噪声和乘性噪声同时存在时，准确提取信号。

本发明是这样实现的，

一种基于MCMC的地面核磁共振信号参数提取方法，包括如下的步骤：

步骤1:利用地面磁共振仪器进行野外实验，获得地面磁共振数据，数据中包含尖峰噪声，工频谐波噪声和环境噪声；

步骤2:对采集得到的磁共振数据进行去噪处理，消除数据中的尖峰噪声和工频谐波噪声，再经过希尔伯特变换和低通滤波处理获得式(1)和式(2)表示的地面核磁共振信号的复包络：

其中V_R和V_I分别为MRS数据的实部和虚部，e₀为初始振幅、

为弛豫时间、df为频率偏量，t为时间，

为初始相位，ε_r是随机噪声的实部分量，ε_i是随机噪声的虚部分量；

步骤3：使用最大似然估计方法得到步骤2中e₀、

df、

四个参数的估计值；

步骤4：以步骤3得到的e₀、

df、

参数的估计值为中心，选择e₀在0～500nV、

在0～1s间、df在-2，～2Hz、

在-π～πrad范围，生成4组服从均匀分布的数据，作为MCMC模型的初始模型；

步骤5:以步骤4得到的结果作为初始模型，e₀、

df、

四个参数的固有特征作为先验信息，输入MCMC方法中进行迭代搜索，得到对应的e₀、

df、

的后验分布；

步骤6：对步骤5得到的特征参数后验信息的样本取中位数，观察出现次数最多，权值最大的数据，为地面磁共振信号中e₀、

df、

的最优估计值。

进一步地，步骤3中对地面磁共振复包络信号进行最大似然估计的过程为：

受到随机噪声的影响，测量的数据的复包络V＝[V_R,V_I]服从高斯分布，则其似然函数为式(3)

其中,v_i为V中第i个地面核磁共振数据的复包络值，n等于V长度，θ代表未知参数

μ代表含噪地面磁共振复包络数据的期望均值，σ代表含噪地面磁共振复包络数据的标准差；

对于公式(3)两侧同时取对数，则得到对数似然函数模型ln L(θ)为：

式(4)分别对4个未知参数求偏导，求解似然方程

得到最大似然估计值四个参数的最大似然估计值分别表示为：e₀,

df,

进一步地，步骤(5)中使用MCMC方法参数提取过程为：

a.初始化t＝0，选取多元高斯分布为转换函数q(θ′|θ)，设定MCMC的初始值为θ⁰，迭代终止的次数为T；

b.令t＝t+1，根据q(θ′|θ)，生成一个候选样本

并根据式(5)，计算其接受概率

其中π(θ)代表平稳分布函数模型，即为地面磁共振复包络信号；

c.根据步骤4中的先验信息，生成一个随机数u，若u＜α(θ′|θ)，则接受候选样本

否则，拒绝候选样本并令θ_t＝θ_t-1；

d.若t＝T，停止迭代，输出对应的e₀、

df、

的后验分布。否则返回步骤b重复迭代。

本发明与现有技术相比，有益效果在于：本发明利用最大似然估计值作为初始模型，待求参数的固有属性作为先验信息，使用马尔可夫链蒙特卡洛算法估计特征参数的后验分布。实现特征参数提取。该方法不仅适用于包含加性噪声的数据，尤其对于包含乘性噪声的数据，仍然能够准确提取信号特征参数，进而提高地面磁共振地下水探测的可靠性。

经实验，本发明方法对地面磁共振数据中加性随机误差和乘性随机误差削弱效率高，准确性好，提高了地面磁共振信号测量数据的可信度。

附图说明

图1为马尔科夫链蒙特卡洛提取地面核磁共振信号参数方法的流程图；

图2a为马尔科夫链蒙特卡洛模型提取的初始振幅后验分布；

图2b为马尔科夫链蒙特卡洛模型提取的弛豫时间后验分布；

图2c为马尔科夫链蒙特卡洛模型提取的初始相位后验分布；

图2d为马尔科夫链蒙特卡洛模型提取的频率偏量后验分布；

图3a为马尔科夫链蒙特卡洛模型提取参数后的磁共振信号实部曲线拟合效果；

图3b为马尔科夫链蒙特卡洛模型提取参数后的磁共振信号虚部曲线拟合效果。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

基于马尔科夫链蒙特卡洛的地面核磁共振信号参数提取方法，为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。参见图1所示。

步骤(1)，采集得到一组MRS数据，其中MRS信号的初始振幅为100nV，弛豫时间为0.2s，拉莫尔频率为2300Hz，频率偏量为1Hz，初始相位为

rad，信号中的加性噪声水平为30nV，乘性噪声噪声水平系数为l。

步骤(2):对步骤1中的磁共振数据进行去噪处理，消除数据中的尖峰噪声和工频谐波噪声，再经过希尔伯特变换和低通滤波处理获得地面核磁共振信号的复包络数据。

步骤(3):通过最大似然估计方法得到e₀、

df、

四个参数的估计值。对应分别为84nV、0.25s、0.9Hz、1.1rad。

步骤3中对地面磁共振复包络信号进行最大似然估计的过程为：

受到随机噪声的影响，测量的数据的复包络V＝[V_R,V_I]服从高斯分布，则其似然函数为式(3)

其中,v_i为V中第i个地面核磁共振数据的复包络值，n等于V长度，θ代表未知参数

μ代表含噪地面磁共振复包络数据的期望均值，σ代表含噪地面磁共振复包络数据的标准差；

对于公式(3)两侧同时取对数，则得到对数似然函数模型ln L(θ)为：

式(4)分别对4个未知参数求偏导，求解似然方程

得到最大似然估计值四个参数的最大似然估计值分别表示为：e₀,

df,

步骤(4)：利用步骤(3)得到的结果e₀、

df、

四个参数，选择合理的范围，生成4组服从均匀分布的数据，作为MCMC方法的初始先验信息模型。本实施例中4个参数的分布范围分别为设定为(0，500)nV、(0，1)s、(-2，2)Hz、(-π,π)rad。

步骤(5):将步骤(4)得到的先验信息输入MCMC方法中进行迭代搜索，得到对应的e₀、

df、

的后验分布，如图2a-d所示。步骤(5)中使用MCMC方法参数提取过程为：

a.初始化t＝0，选取多元高斯分布为转换函数q(θ′|θ)，设定MCMC的初始值为θ⁰，迭代终止的次数为T；

b.令t＝t+1，根据q(θ′|θ)，生成一个候选样本

并根据式(5)，计算其接受概率

其中π(θ)代表平稳分布函数模型，即为地面磁共振复包络信号；

c.根据步骤4中的先验信息，生成一个随机数u，若u＜α(θ′|θ)，则接受候选样本

否则，拒绝候选样本并令θ_t＝θ_t-1；

d.若t＝T，停止迭代，输出对应的e₀、

df、

的后验分布。否则返回步骤b重复迭代。

步骤(6)：对步骤(5)得到的特征参数后验信息的样本取中位数，观察出现次数最多，权值最大的数据，即为地面磁共振信号模型中e₀、

df、

的最优估计值，其结果对应为96nV、0.22s、1.0Hz、1.1rad。通过MCMC方法提取参数后的得到的包络曲线结果如图3a、3b中灰色虚线所示，进一步给出了利用传统非线性拟合方法得到的包络结果，见图3a、3b中灰色实线。通过对比可以得出利用本申请提出的基于马尔科夫链蒙特卡洛的地面核磁共振信号参数提取方法得到的包络信号与真实值基本重合，验证了方法的有效性。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (3)

1.一种基于MCMC的地面核磁共振信号参数提取方法，其特征在于，

包括：

步骤1:利用地面磁共振仪器进行野外实验，获得地面磁共振数据，数据中包含尖峰噪声，工频谐波噪声和环境噪声；

步骤2:对采集得到的磁共振数据进行去噪处理，消除数据中的尖峰噪声和工频谐波噪声，再经过希尔伯特变换和低通滤波处理获得式(1)和式(2)表示的地面核磁共振信号的复包络：

其中V_R和V_I分别为MRS数据的实部和虚部，e₀为初始振幅、

为弛豫时间、df为频率偏量，t为时间，

为初始相位，ε_r是随机噪声的实部分量，ε_i是随机噪声的虚部分量；

步骤3：使用最大似然估计方法得到步骤2中e₀、

df、

四个参数的估计值；

步骤4：以步骤3得到的e₀、

df、

参数的估计值为中心，选择e₀在0～500nV、

在0～1s间、df在-2，～2Hz、

在-π～πrad范围，生成4组服从均匀分布的数据，作为MCMC模型的初始模型；

步骤5:以步骤4得到的结果作为初始模型，e₀、

df、

四个参数的固有特征作为先验信息，输入MCMC方法中进行迭代搜索，得到对应的e₀、

df、

的后验分布；

步骤6：对步骤5得到的特征参数后验信息的样本取中位数，观察出现次数最多，权值最大的数据，为地面磁共振信号中e₀、

df、

的最优估计值。

2.按照权利要求1所述的方法，其特征在于，步骤3中对地面磁共振复包络信号进行最大似然估计的过程为：

受到随机噪声的影响，测量的数据的复包络V＝[V_R,V_I]服从高斯分布，则其似然函数为式(3)

其中,v_i为V中第i个地面核磁共振数据的复包络值，n等于V长度，θ代表未知参数

μ代表含噪地面磁共振复包络数据的期望均值，σ代表含噪地面磁共振复包络数据的标准差；

对于公式(3)两侧同时取对数，则得到对数似然函数模型lnL(θ)为：

式(4)分别对4个未知参数求偏导，求解似然方程

得到最大似然估计值四个参数的最大似然估计值分别表示为：e₀,

df,

3.按照权利要求1所述的方法，其特征在于，步骤(5)中使用MCMC方法参数提取过程为：

a.初始化t＝0，选取多元高斯分布为转换函数q(θ′|θ)，设定MCMC的初始值为θ⁰，迭代终止的次数为T；

b.令t＝t+1，根据q(θ′|θ)，生成一个候选样本

并根据式(5)，计算其接受概率

其中π(θ)代表平稳分布函数模型，即为地面磁共振复包络信号；

c.根据步骤4中的先验信息，生成一个随机数u，若u＜α(θ′|θ)，则接受候选样本

否则，拒绝候选样本并令θ_t＝θ_t-1；

d.若t＝T，停止迭代，输出对应的e₀、

df、

的后验分布。否则返回步骤b重复迭代。

Images