一种基于B样条插值的二维地面核磁共振反演方法技术

本发明专利技术涉及一种基于B样条插值的二维地面核磁共振反演方法，用铺设较少的线圈，只采集几个测点的信号，结合B样条插值的反演方法，仍然可以达到用阵列式线圈半覆盖所得到的解释结果。适用于非层状不均匀水体探测，采用B样条插值的方法对核磁共振信号的初始振幅进行插值计算，扩大了数据矩阵，增加了数据间的相关性，通过二分法搜索最优正则化因子，提高了计算速度，使用高斯牛顿迭代法求解反演目标函数，获得了高精度平滑反演结果。解决了二维核磁共振探测布线时间长，探测效率低，反演结果不平滑、精度不高的难题，可应用于裂隙水和岩溶水等复杂条件地下水的高效率、高分辨率、高精度探测，极大地提高了探测效率，节省了野外测量时间。

【技术实现步骤摘要】

：本专利技术涉及一种地球物理探测方法及其数据的反演解释，尤其是以二维的测量方式同时实现基于B样条插值的地面核磁共振反演方法。
技术介绍
：地面核磁共振(SurfaceNuclearMagneticresonance,简称SNMR)是近年国际上发展起来的一种新的地球物理直接探测地下水的方法，这种在地面直接探测地下介质中氢核丰度的技术，不仅可以用于缺水地区的地下水资源勘查与评价，还可以在地下水引起的堤坝渗漏、矿井突水、隧道涌水、滑坡等地质灾害水源的探测预警中发挥独特的作用。然而，随着SNMR技术的广泛应用，针对裂隙水和岩溶水等复杂条件地下水的探测需求逐渐增多，对于这种非层状不均匀水体的探测，SNMR技术需要从一维垂向测深转向二维地下水成像，有必要探索和发展适合于复杂水体的高效率高分辨率的探测方法以及高精度的二维反演方法。CN201410252327.8公开了一种“基于和声搜索算法的地面核磁共振反演方法”，其先建立了约束R-TLS反演模型，在背景区域电阻率值分布序列和测量信号的初始振幅序列都存在误差的条件下，提高反演精度，后提出了IHS改进的和声搜索算法以求解通过公式推导将该模型转化为受条件约束的非线性优化的问题，消除了反演模型中对含水层最大划分层数的限制，适用于现有层状地质模型中SNMR反演算法对背景区域电阻率的分布存在估计误差的情况。专利CN201410252243.4公开了一种“地面核磁共振二维反演方法”，其用拉直>变换方法对二维的正演模型降维，将其抽象为矩阵方程求解模型，并用最小二乘奇异值分解(LS-SVD)与改进的随机梯度下降法(ISGD)相结合的方法进行反演求解，采用LS-SVD求取矩阵方程的粗略解，在该粗略解的基础上，用ISGD求取其精细解。蒋川东在GeophysicalJournalInternational[2015，200(2)，824-836]上发表的论文“Magneticresonancetomographyusingelongatedtransmitterandin-loopreceiverarraysfortime-efficient2-Dimagingofsubsurfaceaquiferstructures”，首先设计阵列线圈的二维测量模式，即一次测量中，采用一个较大线圈作为发射线圈，多个小线圈同时采集磁共振信号，并给出了基于分辨率半径分析的阵列线圈组合模式的优化设计方法，最后通过二维初始振幅反演(IVI)获取地下含水量信息，对比分析了正方形以及长方形阵列线圈的边对边型和半覆盖型测量模式，得出了当仪器系统的采集通道足够多时，长方形阵列线圈半覆盖测量模式是实现快速高分辨率二维地下水探测的最优测量模式。上述专利技术的基于和声搜索算法只能解决一维地面核磁共振反演问题，横向分辨率低，反演结果不稳定；LS-SVD与ISGD相结合的二维反演方法能对简单的二维含水模型进行反演成像，然而，针对实际探测过程中的裂隙、岩溶等情况复杂含水体，该方法难以保证反演精度；长方形阵列线圈半覆盖测量模式与正方形阵列线圈相比，提高了探测分辨率，但受仪器系统采集通道个数的限制，需要进行多次测量，存在耗时耗力、测量效率低、反演结果精度不高、不平滑等缺点，此外，为提高水平分辨率，采用短边长的线圈对已知长度的测线进行探测及反演成像时，长方形阵列线圈半覆盖模式的测量效率将严重降低。
技术实现思路
：本专利技术的目的在于针对上述现有技术的不足，提供一种适用于非层状不均匀水体的高效率高分辨率探测，既能节省测量时间，又能得到高精度平滑的反演结果，对裂隙水和岩溶水等复杂条件的一种基于B样条插值的二维地面核磁共振反演方法。本专利技术的目的是通过以下技术方案实现的：一种基于B样条插值的二维地面核磁共振反演方法，包括以下步骤：a、采用商业化的有限元软件COMSOL，建立三维可视化模型，在水平地面沿东西方向铺设长方形发射线圈1，在发射线圈1内部等距离铺设第一接收线圈2、第二接收线圈3、第三接收线圈4和第四接收线圈5，通过自适应网格剖分技术精确计算地下任意一点的三维矢量磁场：式中，B为三维矢量磁场；Bx为沿x方向的矢量磁场；By为沿y方向的矢量磁场；Bz为沿z方向的矢量磁场；和为直角坐标系的方向向量；b、计算感应磁场垂直于地磁场方向的分量：式中，B⊥为感应磁场B垂直于地磁场方向的分量；为地磁场的方向向量；c、计算灵敏度二维核函数K(q,r)的表达式；d、提取第一接收线圈2、第二接收线圈3、第三接收线圈4和第四接收线圈5探测的地面核磁共振信号的初始振幅E0(q)1,2,3,4；e、采用均匀B样条插值方法对各线圈探测的初始振幅E0进行插值计算，得到探测线圈之间区域初始振幅E0的近似值扩大数据矩阵，增加数据相关性；f、由B样条插值后新的初始振幅数据矩阵对应的探测线圈位置关系，再次计算地下空间位置的灵敏度核函数KB；g、由数据目标函数Φd和模型目标函数Φm建立反演总体目标函数Φ，表示为：Φ＝Φd+λΦm式中，λ为正则化因子，数据目标函数表示为：式中，Dε为数据的权值，由信号或噪声的不确定度计算得到，模型目标函数表示为：式中，C是平滑度矩阵；h、利用线性搜索的方法自动选取最优正则化因子区间，再利用二分法搜索最优正则化因子，首先给定正则化因子的初始值λ0和搜索步长Δλ，搜索λ的取值区间[λn-1,λn+1]使它成为总体目标函数Φ的单谷区间，确定含根区间后，利用二分法把区间一分为二，逐步减小搜索区间，直至满足式中，λn为本次迭代的最优正则化参数；i、将需要求解的含水量w表示成迭代格式：wk+1＝wk+ηkΔwk式中，k为当前迭代次数，ηk为搜索步长，Δwk为模型增量。将反演总体目标函数Φ对w求偏导得式中，T为矩阵的转置；用高斯牛顿迭代法求解，得到模型增量Δwk，进而求得总体目标函数Φ的最小值对应的含水量w；j、重复上述i过程，直到反演数据误差小于设定误差，获得高精度的反演结果，输出并快速成像。步骤e为：首先得到探测线圈的核磁共振信号E0，由n(n＝2)次B样条插值表达式计算待插值点数据，若待插值点数据误差大于设定误差，计算n+1次B样条待插值点数据；再对待插值点数据误差和设定误差进行比较，直至插值点数据误差小于设定误差，输出插值后的探测线圈之间区域初始振幅E0的近似值假设n+1个控制点的坐标n次B样条曲线段的参数表达式为：式中，Fi,n(t)为n次B样条基函数，其形式为：式中，步骤h包括：A、给定初始值λ0和搜索本文档来自技高网...

【技术保护点】
一种基于B样条插值的二维地面核磁共振反演方法，其特征在于，包括以下步骤：a、采用商业化的有限元软件COMSOL，建立三维可视化模型，在水平地面沿东西方向铺设长方形发射线圈1，在发射线圈1内部等距离铺设第一接收线圈2、第二接收线圈3、第三接收线圈4和第四接收线圈5，通过自适应网格剖分技术精确计算地下任意一点的三维矢量磁场：B=Bxe^x+Bye^y+Bze^z]]>式中，B为三维矢量磁场；Bx为沿x方向的矢量磁场；By为沿y方向的矢量磁场；Bz为沿z方向的矢量磁场；和为直角坐标系的方向向量；b、计算感应磁场垂直于地磁场方向的分量：B⊥=B-(b^0·B)b^0]]>式中，B⊥为感应磁场B垂直于地磁场方向的分量；为地磁场的方向向量；c、计算灵敏度二维核函数K(q,r)的表达式；d、提取第一接收线圈2、第二接收线圈3、第三接收线圈4和第四接收线圈5探测的地面核磁共振信号的初始振幅E0(q)1,2,3,4；e、采用均匀B样条插值方法对各线圈探测的初始振幅E0进行插值计算，得到探测线圈之间区域初始振幅E0的近似值扩大数据矩阵，增加数据相关性；f、由B样条插值后新的初始振幅数据矩阵对应的探测线圈位置关系，再次计算地下空间位置的灵敏度核函数KB；g、由数据目标函数Φd和模型目标函数Φm建立反演总体目标函数Φ，表示为：Φ＝Φd+λΦm式中，λ为正则化因子，数据目标函数表示为：Φd=||Dϵ(E0B-KBw)||22]]>式中，Dε为数据的权值，由信号或噪声的不确定度计算得到，模型目标函数表示为：Φm=||Cw||22]]>式中，C是平滑度矩阵；h、利用线性搜索的方法自动选取最优正则化因子区间，再利用二分法搜索最优正则化因子，首先给定正则化因子的初始值λ0和搜索步长Δλ，搜索λ的取值区间[λn‑1,λn+1]使它成为总体目标函数Φ的单谷区间，确定含根区间后，利用二分法把区间一分为二，逐步减小搜索区间，直至满足Φλn<Φλn-1<Φλn+1]]>式中，λn为本次迭代的最优正则化参数；i、将需要求解的含水量w表示成迭代格式：wk+1＝wk+ηkΔwk式中，k为当前迭代次数，ηk为搜索步长，Δwk为模型增量。将反演总体目标函数Φ对w求偏导得(KBTDϵTDϵKB+λCTC)Δwk=KBTDϵTDϵ(V-KBwk)-λCTCwk]]>式中，T为矩阵的转置；用高斯牛顿迭代法求解，得到模型增量Δwk，进而求得总体目标函数Φ的最小值对应的含水量w；j、重复上述i过程，直到反演数据误差小于设定误差，获得高精度的反演结果，输出并快速成像。...

【技术特征摘要】
1.一种基于B样条插值的二维地面核磁共振反演方法，其特征在于，包括以下步骤：
a、采用商业化的有限元软件COMSOL，建立三维可视化模型，在水平地面沿东西方向铺
设长方形发射线圈1，在发射线圈1内部等距离铺设第一接收线圈2、第二接收线圈3、第三接
收线圈4和第四接收线圈5，通过自适应网格剖分技术精确计算地下任意一点的三维矢量磁
场：
B = B x e ^ x + B y e ^ y + B z e ^ z ]]>式中，B为三维矢量磁场；Bx为沿x方向的矢量磁场；By为沿y方向的矢量磁场；Bz为沿z方
向的矢量磁场；和为直角坐标系的方向向量；
b、计算感应磁场垂直于地磁场方向的分量：
B ⊥ = B - ( b ^ 0 · B ) b ^ 0 ]]>式中，B⊥为感应磁场B垂直于地磁场方向的分量；为地磁场的方向向量；
c、计算灵敏度二维核函数K(q,r)的表达式；
d、提取第一接收线圈2、第二接收线圈3、第三接收线圈4和第四接收线圈5探测的地面
核磁共振信号的初始振幅E0(q)1,2,3,4；
e、采用均匀B样条插值方法对各线圈探测的初始振幅E0进行插值计算，得到探测线圈之
间区域初始振幅E0的近似值扩大数据矩阵，增加数据相关性；
f、由B样条插值后新的初始振幅数据矩阵对应的探测线圈位置关系，再次计算地下空
间位置的灵敏度核函数KB；
g、由数据目标函数Φd和模型目标函数Φm建立反演总体目标函数Φ，表示为：
Φ＝Φd+λΦm式中，λ为正则化因子，数据目标函数表示为：
Φ d = | | D ϵ ( E 0 B - K B w ) | | 2 2 ]]>式中，Dε为数据的权值，由信号或噪声的不确定度计算得到，模型目标函数表示为：
Φ m = | | C w | | 2 2 ]]>式中，C是平滑度矩阵；
h、利用线性搜索的方法自动选取最优正则化因子区间，再利用二分法搜索最优正则化
因子，首先给定正则化因子的初始值λ0和搜索步长Δλ，搜索λ的取值区间[λn-1,λn+1]使它成
为总体目标函数Φ的单谷区间，确定含根区间后，利用二分法把区间一分为二，逐步减小搜
索区间，直至满足
Φ λ n < Φ λ n - 1 < Φ λ n + 1 ]]>式中，λn为本次迭代的最优正则化参数；
i、将需要求解的含水量w表示成迭代格式：
wk+1＝wk+ηkΔwk式中，k为当前迭代次数，ηk为搜索步长，Δwk为模型增量。将反演总体目标函数Φ对w求
偏导得
( K B T D ϵ T D ϵ K B + λC T C ) Δw k = K B T D ϵ T D ϵ ( V - K B w k ) - λC T ...

【专利技术属性】
技术研发人员：林婷婷，冷舒喆，张扬，蒋川东，赵静，万玲，
申请(专利权)人：吉林大学，
类型：发明
国别省市：吉林;22