CN109613617B - 基于磁共振响应信号参数提取的地下水探测方法及系统 - Google Patents

Abstract

本发明公开了一种基于磁共振响应信号参数提取的地下水探测方法及系统。所述探测方法包括：获取实测磁共振响应观测数据；建立所述观测数据中磁共振响应信号向量的极大似然函数，所述磁共振响应信号向量包括平均衰减时间和拉莫尔频率；利用粒子群优化算法，获取所述极大似然函数的最优解，得到最优的平均衰减时间；获取磁共振响应信号的初始振幅；根据所述最优的平均衰减时间和所述初始振幅反演得到地下水信息。本发明提供的探测方法提高了磁共振响应信号特征参数获取的精度。

Description

基于磁共振响应信号参数提取的地下水探测方法及系统

技术领域

本发明涉及地下水探测领域，特别涉及基于磁共振响应信号参数提取的地下水探测方法及系统。

背景技术

在众多的地下水探测技术中，核磁共振探测地下水是一种新的地球物理探测方法，近十几年得到飞速发展。与过去的地下水探测办法相比，利用核磁共振技术找水技术有如下几点优势：一是与间接探测方法相比，只要有水且水层深度在探测的范围内，基于核磁共振技术的探测方法都能直接探测出结果，找水的效率更高，速度明显更快；二是经过设计利用核磁共振找水仪能获得更多的地下水的有关信息，比如含水层深度，含水量的多少，地下含水层的孔隙度等地质信息参数；三是与传统探测技术相比，经济性更强，勘测找水的全过程只需要很短的时间，如果选择钻孔的勘测方案，不仅仅浪费十倍以上的时间，也需要花费近十倍的人力物力。

核磁共振找水仪的基本原理是通过发射线圈向地下水激发能量，使地下水中的氢质子核外电子发生能级跃迁，然后再用接收线圈接收核外电子由高能级向低能级回迁时释放的能量，从而获取磁共振响应信号，由该磁共振响应信号即可反演得到地下水的相关信息。

然而，核磁共振找水仪采集得到的观测数据中除了有磁共振响应信号之外，还有种类多样的人为噪声和环境噪声，产生原因多且十分复杂。噪声干扰无处不在，和需要的磁共振响应信号共存，当噪声比较大的时候，甚至信号会淹没在噪声中。当信号的信噪比很低时，会给信号的提取和信号信息特征的分析带来很大的麻烦。

目前已经有一些信号处理方法被用于噪声背景下核磁共振找水仪磁共振响应信号的参数提取。如叠加法、数字滤波、小波、自适应滤波、自相关、高阶累积量、独立成分分析、经验模态分解等，这一类方法都是采用先消噪来提高信噪比，再进行曲线拟合来求取磁共振响应信号参数的策略。然而，在消噪过程中磁共振响应信号的能量必然有所损失，后续的曲线拟合也会带来提取误差，因此这一类提取方法过程复杂且精度不高，进而影响地下水探测的精度。

发明内容

本发明的目的是提供一种基于磁共振响应信号参数提取的地下水探测方法及系统，以提升地下水探测的精度。

为实现上述目的，本发明提供了如下方案：

本发明提供一种基于磁共振响应信号参数提取的地下水探测方法，所述探测方法包括如下步骤：

获取实测磁共振响应观测数据；

建立所述观测数据中磁共振响应信号向量的极大似然函数，所述磁共振响应信号向量包括平均衰减时间和拉莫尔频率；

利用粒子群优化算法，获取所述极大似然函数的最优解，得到最优的平均衰减时间和最优的拉莫尔频率；

获取磁共振响应信号的初始振幅；

根据所述最优的平均衰减时间和所述初始振幅反演得到地下水信息。

可选的，所述建立所述观测数据中磁共振响应信号向量的极大似然函数，所述磁共振响应信号向量包括平均衰减时间和拉莫尔频率，具体包括：

基于所述观测数据组成的磁共振响应信号观测向量 y＝[y(0)，y(1)，…，y(N-1)]^T，获得高斯白噪声向量Λ＝y-s(T₂，f)E₀，其中， y(0)，y(1)和y(N-1)分别表示第0个、第1个和第N-1个采样点的观测数据，N 为采样点数；s(T₂，f)表示磁共振响应信号向量；

E₀表示初始振幅，T₂表示平均衰减时间，f表示实测地点地磁场的拉莫尔频率；

建立所述高斯白噪声向量Λ的概率密度函数P，

其中，σ²为高斯白噪声的方差；

获取所述概率密度函数P的对数L， L＝-N lnπσ²[y-s(T₂，f)E₀]^T[y-s(T₂，f)E₀]；

根据所述对数L取最大值时，所述磁共振响应信号向量获得极大似然估计，建立所述磁共振响应信号向量的极大似然函数，

可选的，所述利用粒子群优化算法，获取所述极大似然函数的最优解，得到最优的平均衰减时间，具体包括：

构造平均衰减时间和拉莫尔频率的粒子群，确定所述粒子群的每个粒子的初始位置向量和初始速度向量，得到初始粒子群，将每个粒子的初始位置向量设置为每个粒子的初始的个体当前最优解，利用适应度函数，计算初始粒子群中每个粒子的适应度函数值，并将所有粒子的适应度函数值的最大值所对应的粒子的位置向量设置为初始的全局最优解；

更新所述粒子群；

利用适应度函数，计算更新后的粒子群中每个粒子的适应度函数值；

判断第i个粒子的适应度函数值是否大于第i个粒子的个体当前最优解所对应的粒子的适应度函数值，得到第一判断结果；

若所述第一判断结果表示第i个粒子的适应度函数值大于第i个粒子的个体当前最优解所对应的粒子的适应度函数值，则将第i个粒子的位置向量设置为第i个粒子的个体当前最优解；分别令i＝1，2，…M，确定粒子群中每个粒子的个体当前最优解，其中，M表示粒子群中粒子的个数；

确定所有粒子的适应度函数值的最大值，判断所述适应度函数值的最大值是否大于所述全局最优解所对应的粒子的适应度函数值，得到第二判断结果；

若所述适应度函数值的最大值大于所述全局最优解所对应的粒子的适应度函数值，则将所述适应度函数值的最大值所对应的粒子的位置向量设置为全局最优解；

判断迭代次数是否小于最大迭代次数，得到第三判断结果；

若所述第三判断结果表示迭代次数小于最大迭代次数，则使迭代次数增加1，返回步骤“更新所述粒子群”；

若所述第三判断结果表示迭代次数不小于最大迭代次数，则将所述全局最优解输出为最优解，得到最优的平均衰减时间。

可选的，所述构造所述磁共振响应信号向量的粒子群，具体包括：

利用公式

构造第i个粒子的位置向量；其中，X_i ^k表示第 k次迭代第i个粒子的位置向量，x_i1 ^k表示第k次迭代第i个粒子的平均衰减时间，x_i2 ^k表示第k次迭代第i个粒子的拉莫尔频率；

分别令i＝1，2，…，M，构造每个粒子的位置向量，得到第k次迭代粒子群的位置矩阵X^k＝(X₁ ^k，X₂ ^k，…，X_M ^k)^T；

利用公式V_i ^k＝(v_i1 ^k，v_i2 ^k)，构造第i个粒子的速度向量；其中，V_i ^k表示第k 次迭代第i个粒子的速度向量，v_i1 ^k表示第k次迭代第i个粒子的平均衰减时间的速度，v_i2 ^k表示第k次迭代第i个粒子的拉莫尔频率的速度；

分别令i＝1，2，…M，构造每个粒子的速度向量得到第k次迭代粒子群的速度矩阵V^k＝(V₁ ^k，V₂ ^k，…，V_M ^k)^T。

可选的，所述确定所述粒子群的每个粒子的初始位置向量和初始速度向量，具体包括：

利用公式x_i1 ⁰＝T_l+rand_i(0，1)(T_h-T_l)和x_i2 ⁰＝f_l+rand_i(0，1)(f_h-f_l)，确定第i个粒子的初始位置向量，x_i1 ⁰表示第i个粒子的初始的平均衰减时间， x_i2 ⁰表示第i个粒子的初始的拉莫尔频率，rand_i(0，1)表示利用随机数函数获取的第i个粒子的随机参数，T_l表示平均衰减时间的搜索范围的下限，T_h表示平均衰减时间的搜索范围的上限，f_l表示拉莫尔频率的搜索范围的下限， f_h表示拉莫尔频率的搜索范围的上限；

利用公式v_id ⁰＝v_min+rand_i(0，1)(v_max-v_min)，确定第i个粒子的初始速度向量，其中，d表示搜索维度，d＝1，2，当d＝1时，v_id ⁰表示第i个粒子的平均衰减时间的初始速度，当d＝2，v_id ⁰表示第i个粒子的拉莫尔频率的初始速度，其中，v_max表示速度上限，v_min表示速度下限。

可选的，所述适应度函数为：fit(X_i ^k)＝y^Ts(X_i ^k)[s^T(X_i ^k)s(X_i ^k)]^-1s^T(X_i ^k)y，

其中，fit(X_i ^k)表示第k次迭代第i个粒子的适应度函数值，y表示磁共振响应信号观测向量，s(X_i ^k)表示第k次迭代第i个粒子的磁共振响应信号向量。

可选的，所述更新粒子群，具体包括：

利用公式

更新每个粒子的速度向量；

利用公式

更新每个粒子的位置向量；

其中，d表示搜索维度，d＝1，2，当d＝1时，v_id ^k+1和x_id ^k+1分别表示第k+1 次迭代过程中，第i个粒子的平均衰减时间的速度和位置，v_id ^k和x_id ^k分别表示第k次迭代过程中，第i个粒子的平均衰减时间的速度和位置，

表示第 k次迭代第i个粒子的平均衰减时间的当前最优解，

表示第k次迭代平均衰减时间的全局最优解；当d＝2时，v_id ^k+1和x_id ^k+1分别表示第k+1次迭代过程中，第i个粒子的拉莫尔频率的速度和位置，v_id ^k和x_id ^k分别表示第k次迭代过程中，第i个粒子的拉莫尔频率的速度和位置，

表示第k次迭代第i 个粒子的拉莫尔频率的当前最优解，

表示第k次迭代拉莫尔频率的全局最优解；i＝1，2，…M，ω表示惯性权重系数，c₁和c₂分别表示第一加速常数和第二加速常数，c₁，c₂∈[0，2]，r₁表示第一随机数，r₂表示第二随机数，r₁，r₂∈[0，1]。

一种基于磁共振响应信号参数提取的地下水探测系统，所述探测系统包括：

信号获取模块，用于获取实测磁共振响应观测数据；

极大似然函数建立模块，用于建立所述观测数据中磁共振响应信号向量的极大似然函数，所述磁共振响应信号向量包括平均衰减时间和拉莫尔频率；

最优解获取模块，用于利用粒子群优化算法，获取所述极大似然函数的最优解，得到最优的平均衰减时间；

初始振幅获取模块，用于获取磁共振响应信号的初始振幅；

地下水信息获取模块，用于根据所述最优的平均衰减时间和所述初始振幅反演得到地下水信息。

可选的，所述极大似然函数建立模块，具体包括：

噪声向量获取子模块，用于基于所述观测数据组成的磁共振响应信号观测向量y＝[y(0)，y(1)，…，y(N-1)]^T，获得高斯白噪声向量Λ＝y-s(T₂，f)E₀，其中， y(0)，y(1)和y(N-1)分别表示第0个、第1个和第N-1个采样点的观测数据，N 为采样点数；s(T₂，f)表示磁共振响应信号向量；

E₀表示初始振幅，T₂表示平均衰减时间，f表示实测地点地磁场的拉莫尔频率；

概率密度函数建立子模块，用于建立所述高斯白噪声向量Λ的概率密度函数P，

其中，σ²为高斯白噪声的方差；

对数获取子模块，用于获取所述概率密度函数P的对数L， L＝-N lnπσ²[y-s(T，f)E₀]^T[y-s(T，f)E₀]；

极大似然函数建立子模块，用于根据所述对数L取最大值时，所述磁共振响应信号向量获得极大似然估计，建立所述磁共振响应信号向量的极大似然函数，

可选的，所述最优解获取模块，具体包括：

初始化子模块，用于构造平均衰减时间和拉莫尔频率的粒子群，确定所述粒子群的每个粒子的初始位置向量和初始速度向量，得到初始粒子群，将每个粒子的初始位置向量设置为每个粒子的初始的个体当前最优解，利用适应度函数，计算初始粒子群中每个粒子的适应度函数值，并将所有粒子的适应度函数值的最大值所对应的粒子的位置向量设置为初始的全局最优解；

粒子群更新子模块，用于更新所述粒子群；

适应度函数值计算子模块，用于利用适应度函数，计算更新后的粒子群中每个粒子的适应度函数值；

第一判断子模块，用于判断第i个粒子的适应度函数值是否大于第i个粒子的个体当前最优解所对应的粒子的适应度函数值，得到第一判断结果；

第一判断结果处理子模块，用于若所述第一判断结果表示第i个粒子的适应度函数值大于第i个粒子的个体当前最优解所对应的粒子的适应度函数值，则将第i个粒子的位置向量设置为第i个粒子的个体当前最优解；分别令i＝1，2，…M，确定粒子群中每个粒子的个体当前最优解，其中，M表示粒子群中粒子的个数；

第二判断子模块，用于确定所有粒子的适应度函数值的最大值，判断所述适应度函数值的最大值是否大于所述全局最优解所对应的粒子的适应度函数值，得到第二判断结果；

第二判断结果处理子模块，用于若所述适应度函数值的最大值大于所述全局最优解所对应的粒子的适应度函数值，则将所述适应度函数值的最大值所对应的粒子的位置向量设置为全局最优解；

第三判断子模块，用于判断迭代次数是否小于最大迭代次数，得到第三判断结果；

第三判断结果处理子模块，用于若所述第三判断结果表示迭代次数小于最大迭代次数，则使迭代次数增加1，返回步骤“更新所述粒子群”；若所述第三判断结果表示迭代次数不小于最大迭代次数，则将所述全局最优解输出为最优解，得到最优的平均衰减时间。

根据本发明提供的具体实施例，本发明公开了以下技术效果：

本发明公开了一种基于磁共振响应信号参数提取的地下水探测方法及系统。本发明提供的探测方法，基于极大似然方法，并利用粒子群算法获取磁共振响应信号参数(平均衰减时间和拉莫尔频率)，并获取磁共振响应信号的初始振幅，利用获取的平均衰减时间和初始振幅反演得到地下水信息，不需要先消噪再拟合过程，克服了消噪过程造成磁共振响应信号的能量损失，拟合过程造成误差的技术缺陷，提高了磁共振响应信号获取的精度。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。

图1为本发明提供的一种基于磁共振响应信号参数提取的地下水探测方法的流程图；

图2为本发明提供的利用粒子群优化算法，获取极大似然函数最优解的流程图；

图3为本发明提供的一种基于磁共振响应信号参数提取的地下水探测系统的结构意图。

具体实施方式

本发明的目的是提供一种基于磁共振响应信号参数提取的地下水探测方法及系统，以提升地下水探测的精度。

为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对发明作进一步详细的说明。

核磁共振找水仪的磁共振响应信号的实数信号模型如下：

其中，n为采样时刻，E₀为初始振幅，初始振幅的大小和地下含水量成正比，也包含了地下含水层深度、厚度、单位体积含水量等信息。T₂为平均衰减时间(也称驰豫时间)，平均衰减时间反映了地下含水层平均孔隙度的信息。ω₀＝2πf(f为拉莫尔频率)为地磁场的角频率。由实测磁共振响应观测数据中提取关键参数，应用反演软件进行反演就能够获得地下水含量、深度、含水层孔隙度等信息。然而，核磁共振找水仪采集得到的观测数据中除了有磁共振响应信号之外，还有种类多样的人为噪声和环境噪声，产生原因多且十分复杂。基于此，本发明提供一种基于磁共振响应信号参数提取的地下水探测方法，以提取观测数据中的磁共振响应信号关键特征参数，提高地下水探测的精度，核磁共振找水仪获取的磁共振响应信号s(n)中最关键的两个参数是平均衰减时间T₂和初始振幅E₀。

实施例1

本发明实施例1提供一种基于磁共振响应信号参数提取的地下水探测方法。

如图1所示，本发明提供一种基于磁共振响应信号参数提取的地下水探测方法，所述探测方法包括如下步骤：步骤101，获取实测磁共振响应观测数据；步骤102，建立所述观测数据中磁共振响应信号向量的极大似然函数，所述磁共振响应信号向量包括平均衰减时间和拉莫尔频率；步骤103，利用粒子群优化算法，获取所述极大似然函数的最优解，得到最优的平均衰减时间；步骤104，获取磁共振响应信号的初始振幅；步骤105，根据所述最优的平均衰减时间和所述初始振幅反演得到地下水信息。

实施例2

本发明实施例2提供一种基于磁共振响应信号参数提取的地下水探测方法的一个优选的实施方式，本发明的基于磁共振响应信号参数提取的地下水探测方法，不限于本发明实施例2限定的实施方式。

步骤101所述获取实测磁共振响应观测数据，具体包括：利用核磁共振找水仪获取实测核磁共振响应观测数据，在核磁共振找水仪获取的实测核磁共振响应观测数据y(n)中，一般会附加有方差为σ²的高斯白噪声η(n)，则实测核磁共振响应观测数据y(n)可以表示为：

y(n)＝s(n)+η(n)

其中n＝0，1，2...，(N-1)，N为采样点数；s(n)为不含噪声的真实磁共振响应信号；E₀为初始振幅；T₂为平均衰减时间，也称弛豫时间，平均衰减时间反映了地下含水层平均孔隙度的信息；f为实测地点地磁场的拉莫尔频率。

步骤102所述建立所述观测数据中磁共振响应信号向量的极大似然函数，所述磁共振响应信号向量包括平均衰减时间和拉莫尔频率，具体包括：基于所述观测数据组成的磁共振响应信号观测向量y＝[y(0)，y(1)，…，y(N-1)]^T，获得高斯白噪声向量Λ＝y-s(T₂，f)E₀，其中，y(0)，y(1)和y(N-1)分别表示第0 个、第1个和第N-1个采样点的观测数据，N为采样点数；s(T₂，f)表示磁共振响应信号向量；

E₀表示初始振幅，T₂表示平均衰减时间，f表示实测地点地磁场的拉莫尔频率；建立所述高斯白噪声向量Λ的概率密度函数P，其中，σ²为高斯白噪声的方差；获取所述概率密度函数P的对数L， L＝-N lnπσ²[y-s(T₂，f)E₀]^T[y-s(T₂，f)E₀]；根据所述对数L取最大值时，所述磁共振响应信号向量获得极大似然估计，建立所述磁共振响应信号向量的极大似然函数，

粒子群优化算法(Particle Swarm Optimization，PSO)是于1995年由心理学家Eberhart博士和电气工程师Kennedy博士在人工生命的探究中受到启迪而提出的一种基于群体智能的随机寻优算法，其灵感来源于对鸟类群体觅食行为的探究，基本思想是利用鸟群中个体的相互协作和信息共享体制，经过以随机寻优为主的群体智能，使粒子群整体的移动从随机分布到秩序井然，并最终到达最优位置获得最优解。

本发明将全局模式标准PSO算法应用于对磁共振响应信号极大似然函数的优化。具体的，步骤103所述利用粒子群优化算法，获取所述极大似然函数的最优解，得到最优的平均衰减时间，如图2所示，具体包括：

初始化阶段包括：

加速常数(第一加速常数和第二加速常数)、种群数量M、最大迭代次数等参数初始化，具体的，设粒子种群的种群规模是M，最大迭代次数T_max，因为一个磁共振响应信号对应两个参数(平均衰减时间和拉莫尔频率)，所以搜索维数D＝2，归一化平均衰减时间的搜索范围为[T_l，T_h]，其中下限 T_l＝0.001、上限T_h＝0.1。归一化拉莫尔频率的搜索范围为[f_l，f_h]，其中，下限f_l＝0、上限f_h＝1；速度上限设为v_max＝1，速度下限设为v_min＝-1。

粒子群位置、速度进行初始化，具体的，对于该优化问题，每一个解决方案向量(即粒子位置)都由需要初始化，利用公式 x_i1 ⁰＝T_l+rand_i(0，1)(T_h-T_l)和x_i2 ⁰＝f_l+rand_i(0，1)(f_h-f_l)，确定第i个粒子的初始位置向量，x_i1 ⁰表示第i个粒子的初始的平均衰减时间，x_i2 ⁰表示第 i个粒子的初始的拉莫尔频率，rand_i(0，1)表示利用随机数函数获取的第i 个粒子的随机参数，T_l表示平均衰减时间的搜索范围的下限，T_h表示平均衰减时间的搜索范围的上限，f_l表示拉莫尔频率的搜索范围的下限，f_h表示拉莫尔频率的搜索范围的上限；利用公式v_id ⁰＝v_min+rand_i(0，1)(v_max-v_min)，确定第i个粒子的初始速度向量，其中，d表示搜索维度，d＝1，2，当d＝1 时，v_id ⁰表示第i个粒子的平均衰减时间的初始速度，当d＝2，v_id ⁰表示第i 个粒子的拉莫尔频率的初始速度，其中，v_max表示速度上限，v_min表示速度下限。

全局最优解和个体当前最优解初始化，具体的，构造平均衰减时间和拉莫尔频率的粒子群，进一步的，利用公式X_i ^k＝(x_i1 ^k，x_i2 ^k)，构造第i个粒子的位置向量；其中，X_i ^k表示第k次迭代第i个粒子的位置向量，x_i1 ^k表示第k 次迭代第i个粒子的平均衰减时间，x_i2 ^k表示第k次迭代第i个粒子的拉莫尔频率；分别令i＝1，2，…，M，构造每个粒子的位置向量，得到第k次迭代粒子群的位置矩阵X^k＝(X₁ ^k，X₂ ^k，…，X_M ^k)^T；利用公式V_i ^k＝(v_i1 ^k，v_i2 ^k)，构造第 i个粒子的速度向量；其中，V_i ^k表示第k次迭代第i个粒子的速度向量，v_i1 ^k表示第k次迭代第i个粒子的平均衰减时间的速度，v_i2 ^k表示第k次迭代第i 个粒子的拉莫尔频率的速度；分别令i＝1，2，…M，构造每个粒子的速度向量得到第k次迭代粒子群的速度矩阵V^k＝(V₁ ^k，V₂ ^k，…，V_M ^k)^T。将每个粒子的初始位置向量设置为每个粒子的初始的个体当前最优解，利用适应度函数，计算初始粒子群中每个粒子的适应度函数值，并将所有粒子的适应度函数值的最大值所对应的粒子的位置向量设置为初始的全局最优解；所述适应度函数为：fit(X_i ^k)＝y^Ts(X_i ^k)[s^T(X_i ^k)s(X_i ^k)]^-1s^T(X_i ^k)y，其中，fit(X_i ^k)表示第k次迭代第i个粒子的适应度函数值，y表示磁共振响应信号观测向量，s(X_i ^k)表示第 k次迭代第i个粒子的磁共振响应信号向量。

更新阶段包括：

更新所述粒子群，包括由粒子速度、位置更新公式对粒子进行速度、位置更新，具体的，利用公式

更新每个粒子的速度向量；

利用公式

更新每个粒子的位置向量；

其中，d表示搜索维度，d＝1，2，当d＝1时，v_id ^k+1和x_id ^k+1分别表示第k+1 次迭代过程中，第i个粒子的平均衰减时间的速度和位置，v_id ^k和x_id ^k分别表示第k次迭代过程中，第i个粒子的平均衰减时间的速度和位置，

表示第 k次迭代第i个粒子的平均衰减时间的当前最优解，

表示第k次迭代平均衰减时间的全局最优解；当d＝2时，v_id ^k+1和x_id ^k+1分别表示第k+1次迭代过程中，第i个粒子的拉莫尔频率的速度和位置，v_id ^k和x_id ^k分别表示第k次迭代过程中，第i个粒子的拉莫尔频率的速度和位置，

表示第k次迭代第i 个粒子的拉莫尔频率的当前最优解，

表示第k次迭代拉莫尔频率的全局最优解；i＝1，2，…M，ω表示惯性权重系数，c₁和c₂分别表示第一加速常数和第二加速常数，c₁，c₂∈[0，2]，r₁表示第一随机数，r₂表示第二随机数，r₁，r₂∈[0，1]。

利用适应度函数，计算更新后的粒子群中每个粒子的适应度函数值；

判断第i个粒子的适应度函数值是否大于第i个粒子的个体当前最优解所对应的粒子的适应度函数值，得到第一判断结果；

若所述第一判断结果表示第i个粒子的适应度函数值大于第i个粒子的个体当前最优解所对应的粒子的适应度函数值，则将第i个粒子的位置向量设置为第i个粒子的个体当前最优解；分别令i＝1，2，…M，确定粒子群中每个粒子的个体当前最优解；

确定所有粒子的适应度函数值的最大值，判断所述适应度函数值的最大值是否大于所述全局最优解所对应的粒子的适应度函数值，得到第二判断结果；

若所述适应度函数值的最大值大于所述全局最优解所对应的粒子的适应度函数值，则将所述适应度函数值的最大值所对应的粒子的位置向量设置为全局最优解；

判断迭代次数是否小于最大迭代次数，得到第三判断结果；

若所述第三判断结果表示迭代次数小于最大迭代次数，则使迭代次数增加1，返回步骤“更新所述粒子群”；

若所述第三判断结果表示迭代次数不小于最大迭代次数，则将所述全局最优解输出为最优解，得到最优的平均衰减时间。

步骤104所述获取磁共振响应信号的初始振幅，具体包括：本发明采用采用快速傅氏变换(FFT，Fast Fourier Transformation)计算初始振幅E₀。

本发明利用得到的最优的平均衰减时间

将实测磁共振响应信号数据 y(n)＝s(n)+η(n)转换成正弦信号加噪声的形式，具体为：

此时，实测磁共振响应信号数据被转换为正弦信号加噪声的形式。

利用快速傅氏变换计算正弦信号幅值的原理，对x(n)进行快速傅氏变换，可以得到初始振幅。

步骤105所述的地下水信息包括含水层深度、含水量的多少、地下含水层的孔隙度等地质信息参数。

实施例3

本发明实施例3提供一种基于磁共振响应信号参数提取的地下水探测系统。

如图3所示，所述探测系统包括：

信号获取模块301，用于获取实测磁共振响应观测数据；

极大似然函数建立模块302，用于建立所述观测数据中磁共振响应信号向量的极大似然函数，所述磁共振响应信号向量包括平均衰减时间和拉莫尔频率；

最优解获取模块303，用于利用粒子群优化算法，获取所述极大似然函数的最优解，得到最优的平均衰减时间；

初始振幅获取模块304，用于获取磁共振响应信号的初始振幅；

地下水信息获取模块305，用于根据所述最优的平均衰减时间和所述初始振幅反演得到地下水信息。

本发明实施例4提供一种基于磁共振响应信号参数提取的地下水探测系统的一个优选的实施方式，本发明的基于磁共振响应信号参数提取的地下水探测系统，不限于本发明实施例4限定的实施方式。

所述极大似然函数建立模块302，具体包括：噪声向量获取子模块，用于基于所述观测数据组成的磁共振响应信号观测向量 y＝[y(0)，y(1)，…，y(N-1)]^T，获得高斯白噪声向量Λ＝y-s(T₂，f)E₀，其中， y(0)，y(1)和y(N-1)分别表示第0个、第1个和第N-1个采样点的观测数据，N 为采样点数；s(T₂，f)表示磁共振响应信号向量；

E₀表示初始振幅，T₂表示平均衰减时间，f表示实测地点地磁场的拉莫尔频率；

概率密度函数建立子模块，用于建立所述高斯白噪声向量Λ的概率密度函数P，

其中，σ²为高斯白噪声的方差；

对数获取子模块，用于获取所述概率密度函数P的对数L， L＝-N lnπσ²[y-s(T，f)E₀]^T[y-s(T，f)E₀]；

极大似然函数建立子模块，用于根据所述对数L取最大值时，所述磁共振响应信号向量获得极大似然估计，建立所述磁共振响应信号向量的极大似然函数，

所述最优解获取模块303，具体包括：初始化子模块，用于构造平均衰减时间和拉莫尔频率的粒子群，确定所述粒子群的每个粒子的初始位置向量和初始速度向量，得到初始粒子群，将每个粒子的初始位置向量设置为每个粒子的初始的个体当前最优解，利用适应度函数，计算初始粒子群中每个粒子的适应度函数值，并将所有粒子的适应度函数值的最大值所对应的粒子的位置向量设置为初始的全局最优解；

粒子群更新子模块，用于更新所述粒子群；

适应度函数值计算子模块，用于利用适应度函数，计算更新后的粒子群中每个粒子的适应度函数值；

第一判断子模块，用于判断第i个粒子的适应度函数值是否大于第i个粒子的个体当前最优解所对应的粒子的适应度函数值，得到第一判断结果；

第一判断结果处理子模块，用于若所述第一判断结果表示第i个粒子的适应度函数值大于第i个粒子的个体当前最优解所对应的粒子的适应度函数值，则将第i个粒子的位置向量设置为第i个粒子的个体当前最优解；分别令i＝1，2，…M，确定粒子群中每个粒子的个体当前最优解，其中，M表示粒子群中粒子的个数；

第二判断子模块，用于确定所有粒子的适应度函数值的最大值，判断所述适应度函数值的最大值是否大于所述全局最优解所对应的粒子的适应度函数值，得到第二判断结果；

第二判断结果处理子模块，用于若所述适应度函数值的最大值大于所述全局最优解所对应的粒子的适应度函数值，则将所述适应度函数值的最大值所对应的粒子的位置向量设置为全局最优解；

第三判断子模块，用于判断迭代次数是否小于最大迭代次数，得到第三判断结果；

第三判断结果处理子模块，用于若所述第三判断结果表示迭代次数小于最大迭代次数，则使迭代次数增加1，返回步骤“更新所述粒子群”；若所述第三判断结果表示迭代次数不小于最大迭代次数，则将所述全局最优解输出为最优解，得到最优的平均衰减时间。

根据本发明提供的具体实施例，本发明公开了以下技术效果：

本发明公开了一种基于磁共振响应信号参数提取的地下水探测方法及系统。本发明提供的探测方法，基于极大似然方法，并利用粒子群算法获取磁共振响应信号参数(平均衰减时间和拉莫尔频率)，并获取磁共振响应信号的初始振幅，利用获取的平均衰减时间和初始振幅反演得到地下水信息，不需要先消噪再拟合过程，克服了消噪过程造成磁共振响应信号的能量损失，拟合过程造成误差的技术缺陷，提高了磁共振响应信号获取的精度。

本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。对于实施例公开的系统而言，由于其与实施例公开的方法相对应，所以描述的比较简单，相关之处参见方法部分说明即可。

本文中应用了具体个例对发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想，所描述的实施例仅仅是本发明的一部分实施例，而不是全部的实施例，基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

Claims (8)

1.一种基于磁共振响应信号参数提取的地下水探测方法，其特征在于，所述探测方法包括如下步骤：

获取实测磁共振响应观测数据；

建立所述观测数据中磁共振响应信号向量的极大似然函数，所述磁共振响应信号向量包括平均衰减时间和拉莫尔频率；具体包括：

基于所述观测数据组成的磁共振响应信号观测向量y＝[y(0),y(1)，…，y(N-1)]^T，获得高斯白噪声向量Λ＝y-s(T₂,f)E₀，其中，y(0),y(1)和y(N-1)分别表示第0个、第1个和第N-1个采样点的观测数据，N为采样点数；s(T₂,f)表示磁共振响应信号向量；

E₀表示初始振幅，T₂表示平均衰减时间，f表示实测地点地磁场的拉莫尔频率；

建立所述高斯白噪声向量Λ的概率密度函数P，

其中，σ²为高斯白噪声的方差；

获取所述概率密度函数P的对数L，L＝-Nlnπσ²[y-s(T₂,f)E₀]^T[y-s(T₂,f)E₀]；

根据所述对数L取最大值时，所述磁共振响应信号向量获得极大似然估计，建立所述磁共振响应信号向量的极大似然函数，

利用粒子群优化算法，获取所述极大似然函数的最优解，得到最优的平均衰减时间；

获取磁共振响应信号的初始振幅；

根据所述最优的平均衰减时间和所述初始振幅反演得到地下水信息。

2.根据权利要求1所述的一种基于磁共振响应信号参数提取的地下水探测方法，其特征在于，所述利用粒子群优化算法，获取所述极大似然函数的最优解，得到最优的平均衰减时间，具体包括：

构造平均衰减时间和拉莫尔频率的粒子群，确定所述粒子群的每个粒子的初始位置向量和初始速度向量，得到初始粒子群，将每个粒子的初始位置向量设置为每个粒子的初始的个体当前最优解，利用适应度函数，计算初始粒子群中每个粒子的适应度函数值，并将所有粒子的适应度函数值的最大值所对应的粒子的位置向量设置为初始的全局最优解；

更新所述粒子群；

利用适应度函数，计算更新后的粒子群中每个粒子的适应度函数值；

判断第i个粒子的适应度函数值是否大于第i个粒子的个体当前最优解所对应的粒子的适应度函数值，得到第一判断结果；

若所述第一判断结果表示第i个粒子的适应度函数值大于第i个粒子的个体当前最优解所对应的粒子的适应度函数值，则将第i个粒子的位置向量设置为第i个粒子的个体当前最优解；分别令i＝1,2,…M,确定粒子群中每个粒子的个体当前最优解，其中，M表示粒子群中粒子的个数；

确定所有粒子的适应度函数值的最大值，判断所述适应度函数值的最大值是否大于所述全局最优解所对应的粒子的适应度函数值，得到第二判断结果；

若所述适应度函数值的最大值大于所述全局最优解所对应的粒子的适应度函数值，则将所述适应度函数值的最大值所对应的粒子的位置向量设置为全局最优解；

判断迭代次数是否小于最大迭代次数，得到第三判断结果；

若所述第三判断结果表示迭代次数小于最大迭代次数，则使迭代次数增加1，返回步骤“更新所述粒子群”；

若所述第三判断结果表示迭代次数不小于最大迭代次数，则将所述全局最优解输出为最优解，得到最优的平均衰减时间。

3.根据权利要求2所述的一种基于磁共振响应信号参数提取的地下水探测方法，其特征在于：所述构造所述磁共振响应信号向量的粒子群，具体包括：

利用公式X_i ^k＝(x_i1 ^k,x_i2 ^k)，构造第i个粒子的位置向量；其中，X_i ^k表示第k次迭代第i个粒子的位置向量，x_i1 ^k表示第k次迭代第i个粒子的平均衰减时间，x_i2 ^k表示第k次迭代第i个粒子的拉莫尔频率；

分别令i＝1,2,…，M，构造每个粒子的位置向量，得到第k次迭代粒子群的位置矩阵X^k＝(X₁ ^k,X₂ ^k,…,X_M ^k)^T；

利用公式V_i ^k＝(v_i1 ^k,v_i2 ^k)，构造第i个粒子的速度向量；其中，V_i ^k表示第k次迭代第i个粒子的速度向量，v_i1 ^k表示第k次迭代第i个粒子的平均衰减时间的速度，v_i2 ^k表示第k次迭代第i个粒子的拉莫尔频率的速度；

分别令i＝1,2,…M，构造每个粒子的速度向量得到第k次迭代粒子群的速度矩阵V^k＝(V₁ ^k,V₂ ^k,…,V_M ^k)^T。

4.根据权利要求3所述的一种基于磁共振响应信号参数提取的地下水探测方法，其特征在于，所述确定所述粒子群的每个粒子的初始位置向量和初始速度向量，具体包括：

利用公式x_i1 ⁰＝T_l+rand_i(0,1)(T_h-T_l)和x_i2 ⁰＝f_l+rand_i(0,1)(f_h-f_l)，确定第i个粒子的初始位置向量，x_i1 ⁰表示第i个粒子的初始的平均衰减时间，x_i2 ⁰表示第i个粒子的初始的拉莫尔频率，rand_i(0,1)表示利用随机数函数获取的第i个粒子的随机参数，T_l表示平均衰减时间的搜索范围的下限，T_h表示平均衰减时间的搜索范围的上限，f_l表示拉莫尔频率的搜索范围的下限，f_h表示拉莫尔频率的搜索范围的上限；

利用公式v_id ⁰＝v_min+rand_i(0,1)(v_max-v_min)，确定第i个粒子的初始速度向量，其中，d表示搜索维度，d＝1,2，当d＝1时，v_id ⁰表示第i个粒子的平均衰减时间的初始速度，当d＝2，v_id ⁰表示第i个粒子的拉莫尔频率的初始速度，其中，v_max表示速度上限，v_min表示速度下限。

5.根据权利要求2所述的一种基于磁共振响应信号参数提取的地下水探测方法，其特征在于，所述适应度函数为：fit(X_i ^k)＝y^Ts(X_i ^k)[s^T(X_i ^k)s(X_i ^k)]^-1s^T(X_i ^k)y；

其中，fit(X_i ^k)表示第k次迭代第i个粒子的适应度函数值，y表示磁共振响应信号观测向量，s(X_i ^k)表示第k次迭代第i个粒子的磁共振响应信号向量。

6.根据权利要求2所述的一种基于磁共振响应信号参数提取的地下水探测方法，其特征在于，所述更新粒子群，具体包括：

利用公式

更新每个粒子的速度向量；

利用公式

更新每个粒子的位置向量；

其中，d表示搜索维度，d＝1,2，当d＝1时，v_id ^k+1和x_id ^k+1分别表示第k+1次迭代过程中，第i个粒子的平均衰减时间的速度和位置，v_id ^k和x_id ^k分别表示第k次迭代过程中，第i个粒子的平均衰减时间的速度和位置，

表示第k次迭代第i个粒子的平均衰减时间的当前最优解，

表示第k次迭代平均衰减时间的全局最优解；当d＝2时，v_id ^k+1和x_id ^k+1分别表示第k+1次迭代过程中，第i个粒子的拉莫尔频率的速度和位置，v_id ^k和x_id ^k分别表示第k次迭代过程中，第i个粒子的拉莫尔频率的速度和位置，

表示第k次迭代第i个粒子的拉莫尔频率的当前最优解，

表示第k次迭代拉莫尔频率的全局最优解；i＝1,2,…M，ω表示惯性权重系数，c₁和c₂分别表示第一加速常数和第二加速常数，c₁,c₂∈[0,2]，r₁表示第一随机数，r₂表示第二随机数，r₁,r₂∈[0,1]。

7.一种基于磁共振响应信号参数提取的地下水探测系统，其特征在于，所述探测系统包括：

信号获取模块，用于获取实测磁共振响应观测数据；

极大似然函数建立模块，用于建立所述观测数据中磁共振响应信号向量的极大似然函数，所述磁共振响应信号向量包括平均衰减时间和拉莫尔频率；具体包括：

噪声向量获取子模块，用于基于所述观测数据组成的磁共振响应信号观测向量y＝[y(0),y(1)，…，y(N-1)]^T，获得高斯白噪声向量Λ＝y-s(T₂,f)E₀，其中，y(0),y(1)和y(N-1)分别表示第0个、第1个和第N-1个采样点的观测数据，N为采样点数；s(T₂,f)表示磁共振响应信号向量；

E₀表示初始振幅，T₂表示平均衰减时间，f表示实测地点地磁场的拉莫尔频率；

概率密度函数建立子模块，用于建立所述高斯白噪声向量Λ的概率密度函数P，其中，σ²为高斯白噪声的方差；

对数获取子模块，用于获取所述概率密度函数P的对数L，L＝-Nlnπσ²[y-s(T,f)E₀]^T[y-s(T,f)E₀]；

极大似然函数建立子模块，用于根据所述对数L取最大值时，所述磁共振响应信号向量获得极大似然估计，建立所述磁共振响应信号向量的极大似然函数，

最优解获取模块，用于利用粒子群优化算法，获取所述极大似然函数的最优解，得到最优的平均衰减时间；

初始振幅获取模块，用于获取磁共振响应信号的初始振幅；

地下水信息获取模块，用于根据所述最优的平均衰减时间和所述初始振幅反演得到地下水信息。

8.根据权利要求7所述的一种基于磁共振响应信号参数提取的地下水探测系统，其特征在于，所述最优解获取模块，具体包括：

初始化子模块，用于构造平均衰减时间和拉莫尔频率的粒子群，确定所述粒子群的每个粒子的初始位置向量和初始速度向量，得到初始粒子群，将每个粒子的初始位置向量设置为每个粒子的初始的个体当前最优解，利用适应度函数，计算初始粒子群中每个粒子的适应度函数值，并将所有粒子的适应度函数值的最大值所对应的粒子的位置向量设置为初始的全局最优解；

粒子群更新子模块，用于更新所述粒子群；

适应度函数值计算子模块，用于利用适应度函数，计算更新后的粒子群中每个粒子的适应度函数值；

第一判断子模块，用于判断第i个粒子的适应度函数值是否大于第i个粒子的个体当前最优解所对应的粒子的适应度函数值，得到第一判断结果；

第一判断结果处理子模块，用于若所述第一判断结果表示第i个粒子的适应度函数值大于第i个粒子的个体当前最优解所对应的粒子的适应度函数值，则将第i个粒子的位置向量设置为第i个粒子的个体当前最优解；分别令i＝1,2,…M,确定粒子群中每个粒子的个体当前最优解，其中，M表示粒子群中粒子的个数；

第二判断子模块，用于确定所有粒子的适应度函数值的最大值，判断所述适应度函数值的最大值是否大于所述全局最优解所对应的粒子的适应度函数值，得到第二判断结果；

第二判断结果处理子模块，用于若所述适应度函数值的最大值大于所述全局最优解所对应的粒子的适应度函数值，则将所述适应度函数值的最大值所对应的粒子的位置向量设置为全局最优解；

第三判断子模块，用于判断迭代次数是否小于最大迭代次数，得到第三判断结果；

第三判断结果处理子模块，用于若所述第三判断结果表示迭代次数小于最大迭代次数，则使迭代次数增加1，返回步骤“更新所述粒子群”；若所述第三判断结果表示迭代次数不小于最大迭代次数，则将所述全局最优解输出为最优解，得到最优的平均衰减时间。

Images