本发明涉及一种基于循环相关的磁共振响应信号的参数提取方法及系统，用于核磁共振找水仪的磁共振响应信号中参数的提取。通过循环相关处理方法，抑制了磁共振响应信号中的高斯白噪声成分和工频谐波噪声成分，从而提高了磁共振响应信号参数的提取精度。然后对得到的待解算循环相关函数进行参数估计，提高了对地下水的勘探精度，为核磁共振测深技术在地下水勘探的应用提供了技术支持。
1.一种基于循环相关的磁共振响应信号参数提取方法，其特征在于，所述方法包括：接收磁共振响应数据；
利用循环相关方法对所述磁共振响应数据进行处理，得到待解算循环相关函数，具体包括：对所述磁共振响应数据进行循环相关处理，得到所述磁共振响应数据的第一循环相关函数，所述磁共振响应数据的第一循环相关函数 为： 其中，x(n)为磁共振响应数据，n为采样时刻，N为采样点数，τ为移位个数，α为循环频率，n＝1,
2,…,N；选取预设循环频率带入所述第一循环相关函数中，得到待解算循环相关函数，当循环频率α取不同值时，有：
其中，为初始相位，E0为初始振幅，T2为平均衰减时间；
取α＝2fL为预设循环频率，将预设循环频率代入以上公式中，就得到了仅与s(n)s (n+τ)子项相关的循环函数：
其中fL为地磁场的拉摩尔频率，在采集磁共振响应信号前由地磁仪测量得到；
利用参数估计方法对所述待解算循环相关函数进行求解，得到子空间旋转算符Φ，所述子空间旋转算符Φ计算得到磁共振响应信号的平均衰减时间，完成所述磁共振响应信号参数的提取。
2.根据权利要求1所述的基于循环相关的磁共振响应信号参数提取方法，其特征在于，所述预设循环频率为仅与磁共振响应信号相关的循环频率。
3.根据权利要求1所述的基于循环相关的磁共振响应信号参数提取方法，其特征在于，所述利用参数估计方法对所述待解算循环相关函数进行求解时，具体采用旋转不变参数估计方法对所述待解算循环相关函数进行求解。
4.根据权利要求3所述的基于循环相关的磁共振响应信号参数提取方法，其特征在于，所述采用旋转不变参数估计方法对所述待解算循环相关函数进行求解具体包括：α
将所述待解算循环相关函数带入所述循环相关矩阵R 和所述辅助循环相关矩阵 中，并构造矩阵束
对所述矩阵束 进行广义特征值分解，得到子空间旋转算符Φ；
根据所述子空间旋转算符Φ计算得到平均衰减时间T2。
5.根据权利要求4所述的基于循环相关的磁共振响应信号参数提取方法，其特征在于，所述根据所述子空间旋转算符Φ计算得到平均衰减时间T2具体包括：根据公式 计算平均衰减时间T2；其中，Re函数表示取复数的实部。
6.一种基于循环相关的磁共振响应信号参数提取系统，其特征在于，所述系统包括：接收模块，用于接收磁共振响应数据；
循环相关函数计算模块，用于利用循环相关方法对所述磁共振响应数据进行处理，得到待解算循环相关函数，具体包括：对所述磁共振响应数据进行循环相关处理，得到所述磁共振响应数据的第一循环相关函数，所述磁共振响应数据的第一循环相关函数 为：其中，x(n)为磁共振响应数据，n为采样时刻，N为采样点数，τ为移位个数，α为循环频率，n＝1,2,…,N；选取预设循环频率带入所述第一循环相关函数中，得到待解算循环相关函数，当循环频率α取不同值时，有：其中， 为初始相位，τ为移位个数，E0为初始振幅，n为采样时刻，T2为平均衰减时间；
取α＝2fL为预设循环频率，将预设循环频率代入以上公式中，就得到了仅与s(n)s (n+τ)子项相关的循环函数：
其中fL为地磁场的拉摩尔频率，在采集磁共振响应信号前由地磁仪测量得到；
参数估计模块，用于利用参数估计方法对所述待解算循环相关函数进行求解，得到子空间旋转算符Φ，所述子空间旋转算符Φ计算得到所述磁共振响应信号的平均衰减时间，完成磁共振响应信号参数的提取。
7.根据权利要求6所述的基于循环相关的磁共振响应信号参数提取系统，其特征在于，所述参数估计模块包括旋转不变参数估计单元，所述旋转不变参数估计单元用于具体采用旋转不变参数估计方法对所述待解算循环相关函数进行求解。
8.根据权利要求7所述的基于循环相关的磁共振响应信号参数提取系统，其特征在于，所述旋转不变参数估计单元具体包括：α
矩阵构造子单元，用于构造循环相关矩阵R和辅助循环相关矩阵α
矩阵束构造子单元，用于将所述待解算循环相关函数带入所述循环相关矩阵R和所述辅助循环相关矩阵 中，并构造矩阵束特征值分解子单元，用于对所述矩阵束 进行广义特征值分解，得到子空间旋转算符Φ；
参数求解子单元，用于根据所述子空间旋转算符Φ计算得到平均衰减时间T2。
一种基于循环相关的磁共振响应信号参数提取方法及系统
技术领域
[0001] 本发明涉及核磁共振技术领域，特别是涉及一种基于循环相关的磁共振响应信号参数提取方法及系统。
背景技术
[0002] 核磁共振测深技术是目前世界上唯一直接有效探测地下水的地球物理勘探方法。核磁共振探测技术的原理是地下水中的氢核通过人工感应电磁场形成核磁矩。与过去的地
下水探测办法相比，利用核磁共振技术找水有如下几点优势：一是与间接探测方法相比，只
要有水且水层深度在探测的范围内，基于核磁共振技术的探测方法都能直接探测出结果，
找水的效率更高，速度明显更快；二是经过设计，利用核磁共振找水仪能获得更多的地下水
的有关信息，比如含水层深度，含水量的多少，地下含水层的孔隙度等等地质信息参数；三
是与传统探测技术相比，经济性更强，勘测找水的全过程只需要很短的时间，如果选择钻孔
的勘测方案，不仅仅浪费十倍以上的时间，也需要花费近十倍的人力物力。
[0003] 核磁共振找水仪的基本原理是通过发射线圈向地下水激发能量，使地下水中的氢质子核外电子发生能级跃迁，然后再用接收线圈接收核外电子由高能级向低能级回迁时释
放的能量，从而获取核磁共振响应信号，在核磁共振响应信号中提取关键参数，应用反演软
件即可反演得到地下水的相关信息。
[0004] 然而，核磁共振找水仪采集得到的观测数据中除了有磁共振响应信号之外，还有种类多样的人为噪声和环境噪声，产生原因多且十分复杂。噪声干扰无处不在，和我们需要
的磁共振响应信号共存，当噪声比较大的时候，甚至信号会淹没在噪声中，给信号的提取和
信号特征的分析带来很大的困难。
[0005] 基于此，本发明提供一种能够有效滤除各种噪声的高精度磁共振响应信号参数提取方法及系统。
发明内容
[0006] 本发明的目的是提供一种能够有效滤除各种噪声的磁共振响应信号参数提取方法及系统，基于循环相关方法，能够极大地提升磁共振响应数据的信噪比，继而有效提升磁
共振响应信号的参数提取精度。
[0007] 为实现上述目的，本发明提供了如下方案：
[0008] 一种基于循环相关的磁共振响应信号参数提取方法，所述方法包括：
[0009] 接收磁共振响应数据；
[0010] 利用循环相关方法对所述磁共振响应数据进行处理，得到待解算循环相关函数；
[0011] 利用参数估计方法对所述待解算循环相关函数进行求解，得到磁共振响应信号的平均衰减时间，实现磁共振响应信号参数的提取。
[0012] 可选的，所述利用循环相关方法对磁共振响应数据进行处理，得到待解算循环相关函数具体包括：
[0013] 对所述磁共振响应信号进行循环相关处理，得到所述磁共振响应信号的第一循环相关函数；
[0014] 选取预设循环频率带入所述第一循环相关函数中，得到待解算循环相关函数。
[0015] 可选的，所述磁共振响应信号的第一循环相关函数 为：
[0016]
[0017] 其中，x(n)为磁共振响应数据，n为采样时刻，N为采样点数，τ为移位个数，α为循环频率，n＝1,2,…,N。
[0018] 可选的，所述预设循环频率为仅与磁共振响应信号相关的循环频率。
[0019] 可选的，所述利用参数估计方法对所述待解算循环相关函数进行求解时，具体采用旋转不变参数估计方法对所述待解算循环相关函数进行求解。
[0020] 可选的，所述采用旋转不变参数估计方法对所述待解算循环相关函数进行求解具体包括：
α
[0021] 构造循环相关矩阵R和辅助循环相关矩阵
[0022] 将所述待解算循环相关函数带入所述循环相关矩阵Rα和所述辅助循环相关矩阵中，并构造矩阵束
[0023] 对所述矩阵束 进行广义特征值分解，得到子空间旋转算符Φ；
[0024] 根据所述子空间旋转算符Φ计算得到平均衰减时间T2。
[0025] 可选的，所述根据所述子空间旋转算符Φ计算得到平均衰减时间T2具体包括：
[0026] 根据公式 计算平均衰减时间T2，其中，Re函数表示取复数的实部。
[0027] 本发明还提供了一种基于循环相关的磁共振响应信号参数提取系统，所述系统包括：
[0028] 接收模块，用于接收磁共振响应信号；
[0029] 循环相关函数计算模块，用于利用循环相关方法对所述磁共振响应信号进行处理，得到待解算循环相关函数；
[0030] 参数估计模块，用于利用参数估计方法对所述待解算循环相关函数进行求解，得到所述磁共振响应信号的平均衰减时间，完成磁共振响应信号参数的提取。
[0031] 可选的，所述参数估计模块包括旋转不变参数估计单元，所述旋转不变参数估计单元用于具体采用旋转不变参数估计方法对所述待解算循环相关函数进行求解。
[0032] 可选的，所述旋转不变参数估计单元具体包括：α
[0033] 矩阵构造子单元，用于构造循环相关矩阵R和辅助循环相关矩阵
[0034] 矩阵束构造子单元，用于将所述待解算循环相关函数带入所述循环相关矩阵Rα和所述辅助循环相关矩阵 中，并构造矩阵束
[0035] 特征值分解子单元，用于对所述矩阵束 进行广义特征值分解，得到子空间旋转算符Φ；
[0036] 参数求解子单元，用于根据所述子空间旋转算符Φ计算得到平均衰减时间T2。
[0037] 根据本发明提供的具体实施例，本发明公开了以下技术效果：
[0038] 本发明提供了一种磁共振响应信号参数提取方法及系统，通过采用循环相关方法，设置合适的循环频率，有效滤除了不同频率干扰，最大程度的保留了信号信息，极大地
提升了磁共振响应信号的信噪比，从而提高了磁共振响应信号参数提取的精确度。
附图说明
[0039] 为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施
例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获
得其他的附图。
[0040] 图1为本发明实施例提供的一种基于循环相关的磁共振响应信号参数提取方法流程图；
[0041] 图2为本发明实施例提供的一种基于循环相关的磁共振响应信号参数提取系统框图；
[0042] 图3为本发明实施例提供的旋转不变参数估计单元结构示意图。
具体实施方式
[0043] 下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于
本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他
实施例，都属于本发明保护的范围。
[0044] 本发明的目的是提供一种磁共振响应信号参数提取方法及系统，能够有效滤除各种噪声，提高对磁共振响应信号参数提取的精度。
[0045] 采用本发明的磁共振响应信号参数提取方法及系统能够解决现有技术中利用磁共振响应信号探寻地下水的关键难点问题，即磁共振响应信号极其微弱(通常只有几十纳
伏)，而工频谐波噪声以及白噪声强度较高，严重影响磁共振响应信号重要参数提取精度的
问题。
[0046] 为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。
[0047] 实施例1
[0048] 如图1所示，本实施例提供一种基于循环相关的磁共振响应信号参数提取方法，具体包括：
[0049] 步骤101：接收磁共振响应数据；
[0050] 步骤102：利用循环相关方法对所述磁共振响应数据进行处理，得到待解算循环相关函数；
[0051] 步骤103：利用参数估计方法对所述待解算循环相关函数进行求解，得到磁共振响应信号的平均衰减时间，完成磁共振响应信号参数的提取。
[0052] 本实施例中通过循环相关方法抑制了磁共振响应数据中的高斯白噪声和工频谐波噪声成分，最大程度的保留了磁共振响应信号信息，从而提高了核磁共振找水仪的测量
精度。
[0053] 为了更清楚的展示本实施例提供的参数提取方法的优越性，下面对所述方法进行进一步的解释说明。
[0054] 首先，将核磁共振找水仪磁共振响应信号的复数信号模型表示如下：
[0055]
[0056] 其中，s(n)表示磁共振响应信号，E0为初始振幅，T2为平均衰减时间(也称驰豫时间)，n为采样时刻， 为初始相位，ω0＝2πfL(fL为拉摩尔频率)为地磁场的角频率。
[0057] 磁共振响应信号中的关键参数与地下水的含量、深度、含水层孔隙度等密切相关。如，初始振幅的大小和地下含水量成正比，也包含了地下含水层深度、厚度、单位体积含水
量等信息；平均衰减时间反应了地下含水层平均孔隙度的信息；初始相位反映了含水层的
导电率，通常，导电率为一个常数。因此，通过对磁共振响应信号参数进行提取，并通过软件
反演就能获取地下水的相关信息。
[0058] 然而，在实际应用过程中，核磁共振找水仪获取的实测磁共振响应观测信号数据中，往往会附加有白噪声和工频谐波噪声。因此，接收到的核磁共振找水仪的磁共振响应数
据x(n)可以表示为：
[0059] x(n)＝s(n)+η(n)+ε(n) (2)
[0060] 其中，η(n)为高斯白噪声，ε(n)为工频谐波噪声。工频谐波噪声主要来源于变压器、电线、工厂及居民用电。离散的工频谐波噪声表达形式为 其中
fc为工频谐波的基频，通常为50Hz或60Hz，Ab和φb分别为工频谐波分量的幅值和相位，B为
工频谐波分量个数。
[0061] 在接收到磁共振响应数据x(n)后，可以对x(n)进行预处理，滤除x(n)中明显易去除的噪声，提高精确度的同时为后续步骤的运算处理降低复杂度。
[0062] 然后利用循环相关方法对磁共振响应数据x(n)进行处理，抑制x(n)中的高斯白噪声和工频谐波噪声，具体包括：
[0063] 对磁共振响应数据x(n)进行循环相关处理，得到x(n)的第一循环相关函数
[0064]
[0065] 其中，α为循环相关函数的循环频率，τ为移位个数。
[0066] 然后将包含不含噪声的磁共振响应信号s(n)、白噪声η(n)和工频谐波噪声ε(n)三部分的磁共振响应数据x(n)带入到(3)式中：
[0067]
[0068] 由(4)式可以看出，具有相关性的子式仅有s(n)s*(n+τ)、η(n)η*(n+τ)和ε(n)ε*(n+τ)，根据循环相关函数性质，可知其余项的循环相关函数均为0。由此对第一循环相关函数
可以变换为：
[0069]
[0070] 其中，s(n)s*(n+τ)子项对应的循环相关函数为：
[0071]
[0072] 由公式(6)可以看出，当循环频率α取不同值时，有：
[0073]
[0074] 取α＝2fL为预设循环频率，将预设循环频率带入所述公式(7)中，就得到了仅与s*
(n)s(n+τ)子项相关的循环函数：
[0075]
[0076] 其中，fL为地磁场的拉摩尔频率，在采集磁共振响应信号前由地磁仪测量得到。
[0077] 若工频谐波频率和白噪声频率均与磁共振响应信号频率不同，则取循环频率为α* *
＝2fL时，公式(5)中的η(n)η(n+τ)和ε(n)ε(n+τ)子项所对应的循环相关函数均为0。因此，
当循环频率为α＝2fL时，公式(5)中的所述第一循环相关函数为：
[0078]
[0079] 公式(9)即为待解算循环相关函数。
[0080] 本实施例中，具体选择预设循环频率为2fL。但本领域的技术人员应当理解，预设循环频率的选择是与磁共振响应信号s(n)的特性密切相关的，本实施例中将2fL作为预设
循环频率仅仅是为了更加清楚的对本发明的技术方案进行说明，并不应理解为对本发明的
具体限制，实际上，选择任何能够实现噪声滤除的预设循环频率都将落入本发明的保护范
围内。
[0081] 本实施例中通过结合磁共振响应信号的自身特性，对磁共振响应数据进行划分来选择预设循环频率，将循环相关函数中与循环频率无关的白噪声和工频谐波噪声部分进行
剔除，从而仅保留了与循环频率相关的磁共振响应信号信息，为进一步的参数估计提供了
准确、简洁的数据基础，从而提高了磁共振响应信号参数提取的精度，进而提高了核磁共振
找水仪对地下水的检测精度。
[0082] 在得到待解算循环相关函数后，利用参数估计方法对待解算循环相关函数进行求解。为了更清楚的对求解过程进行解释，本实施例以具体采用旋转不变参数估计方法对所
述待解算循环相关函数进行求解为例进行详细说明。
[0083] ESPRIT方法借助旋转不变技术估计信号参数(Estimation of Signal Parameters via Rotational Invariance Techniques)，采用旋转不变参数估计方法对所
述待解算循环相关函数进行求解具体包括：
[0084] 首先，构造一个p×p(p>1)维的循环相关矩阵Rα和辅助循环相关矩阵 然后将τ＝‑p+1,‑p+2,…,0,…,p‑2,p‑1带入矩阵中，得到：
[0085]
[0086]
[0087] 由于Rα和 之间具有子空间旋转不变关系，定义Γ为与矩阵束 相对应的广义特征值矩阵，则矩阵Γ可以表示为：
[0088]
[0089] 其中， 称为旋转算符；
[0090] 将待解算循环相关函数式(9)代入矩阵Rα和 中，计算矩阵对 的广义特征值分解，所得的位于单位圆上的广义特征值z即为子空间旋转算符Φ，剩余的p‑1个广义
特征值恒等于零。
[0091] 由此，根据公式 计算平均衰减时间T2，其中，Re函数表示取复数的实部。
[0092] 当然，采用旋转不变信号参数估计方法对待解算循环相关函数进行求解来得到磁共振响应信号的平均衰减时间仅是本实施例所提供的一个具体实施方式，并不是对本发明
保护范围的进一步限制。其他参数估计方法如将循环相关函数的结构特征作为先验信息，
构建冗余字典来描述循环频率结构来对循环相关函数进行估计等也落入本发明的保护范
围内，只要能够实现对待解算循环相关函数进行求解估计即可。
[0093] 实施例2
[0094] 本实施例提供了一种基于循环相关的磁共振响应信号参数提取系统，如图2所示，所述系统包括：
[0095] 接收模块M1，用于接收磁共振响应数据；
[0096] 循环相关函数计算模块M2，用于利用循环相关方法对所述磁共振响应数据进行处理，得到待解算循环相关函数；
[0097] 参数估计模块M3，用于利用参数估计方法对所述待解算循环相关函数进行求解，得到所述磁共振响应信号的平均衰减时间，完成磁共振响应信号参数的提取。
[0098] 为了更好的实现功能，所述参数估计模块M2包括旋转不变参数估计单元，所述旋转不变参数估计单元用于具体采用旋转不变参数估计方法对所述待解算循环相关函数进
行求解。
[0099] 如图3所示，所述旋转不变参数估计单元具体包括：α
[0100] 矩阵构造子单元M31，用于构造循环相关矩阵R和辅助循环相关矩阵
[0101] 矩阵束构造子单元M32，用于将所述待解算循环相关函数带入所述循环相关矩阵Rα
和所述辅助循环相关矩阵 中，并构造矩阵束
[0102] 特征值分解子单元M33，用于对所述矩阵束 进行广义特征值分解，得到子空间旋转算符Φ；
[0103] 参数求解子单元M34，用于根据所述子空间旋转算符Φ计算得到平均衰减时间T2。
[0104] 本说明书中每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。对于实施例公开的系统而言，由于其与实施例公开的方法
相对应，所以描述的比较简单，相关之处参见方法部分说明即可。
[0105] 本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据
本发明的思想，在具体实施方式及应用范围上均会有改变之处。综上所述，本说明书内容不
应理解为对本发明的限制。
