技术领域：本发明涉及一种地球物理探测方法及其数据的反演解释，尤其是以二维的测量方式同时实现基于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和搜索步长Δλ，利用线性搜索的方法自动选取最优正则化因子区间[λn-1,λn+1]，计算λn＝λ0-(n-1)Δλ对应的总体目标函数Φ和数据目标函数Φd，为保证数据的拟合误差逐阶减小，定义数据拟合目标项系数β通常选取[0,1]区间内较大的值；B、在区间[λn-1,λn+1]内判断的符号，若大于零，则重新搜索最优正则化因子区间；若小于零，则确定[λn-1,λn+1]为满足条件的单谷区间。确定含根区间后，求区间的中点λn＝1/(2(λn-1+λn+1))，计算λn对应的总体目标函数Φ和数据目标函数Φd的值，直到λn满足则λn为本次迭代的最优正则化参数。有益效果：本发明公开的一种基于B样条插值的二维地面核磁共振反演方法，用铺设较少的线圈，只采集几个测点的信号，结合B样条插值的反演方法，仍然可以达到用阵列式线圈半覆盖所得到的解释结果。适用于非层状不均匀水体探测，采用B样条插值的方法对核磁共振信号的初始振幅进行插值计算，扩大了数据矩阵，增加了数据间的相关性，通过二分法搜索最优正则化因子，提高了计算速度，使用高斯牛顿迭代法求解反演目标函数，获得了高精度平滑反演结果。解决了二维核磁共振探测布线时间长，探测效率低，反演结果不平滑、精度不高的难题，可应用于裂隙水和岩溶水等复杂条件地下水的高效率、高分辨率、高精度探测，极大地提高了探测效率，节省了野外测量时间。附图说明：图1为二维地面核磁共振探测线圈布设位置关系图图2为三维矢量磁场分解示意图图3为均匀B样条插值方法流程图图4为利用二分法搜索最优正则化因子流程图图5为基于B样条插值的二维地面核磁共振反演方法流程图1发射线圈，2第一接收线圈，3第二接收线圈，4第三接收线圈，5第四接收线圈，xeast正东方向，ysouth正南方向，zdown垂直向下方向。具体实施方式：下面结合附图和实施例对本发明作进一步的详细说明：如图1所示，该方法通过发射线圈1向地下发射频率为当地拉莫尔频率的交变电流，产生感应磁场。在感应磁场的作用下，地下水中氢质子产生能级跃迁，大量氢质子跃迁到高能级上。当撤去供电电流，这些高能级氢质子便逐渐回到低能级状态，释放出大量的具有拉莫尔频率的能量子，在第一接收线圈2、第二接收线圈3、第三接收线圈4和第四接收线圈5中感应出核磁共振信号。如图2所示，由于地下任意一点的感应磁场可沿x、y、z三个方向进行分解，氢质子核自旋过程中，感应磁场中只有垂直于地磁场的分量起实际作用，因此，设地磁场方向沿x轴，利用感应磁场的y和z分量计算垂直于地磁场方向的分量，利用该感应磁场的垂直分量计算灵敏度核函数。对地下空间的核磁共振响应灵敏度核函数进行计算：核磁共振响应灵敏度核函数K(q,r)的计算公式：式中，ωL为激发频率，M0为氢原子核的宏观磁化强度，γp为旋磁比，q为脉冲矩序列，I0为单位电流强度，为发射场顺时针分量，为接收场逆时针分量，ξT和ξR为发射和接收相位。求解二维核函数，可通过在直角坐标内对上述关于位置r的三维核函数沿y方向进行积分，后化简为关于深度z和一个剖面x的二维表达式：式中，K2D为关于深度z和一个剖面x的二维核函数，q为激发脉冲矩。如图3所示，采用均匀B样条插值方法对各线圈内的初始振幅E0进行插值计算，首先得到探测线圈的核磁共振信号E0，由n(n＝2)次B样条插值表达式计算待插值点数据，若待插值点数据误差大于设定误差，计算n+1次B样条待插值点数据，再对待插值点数据误差和设定误差进行比较，直至插值点数据误差小于设定误差，输出插值后的探测线圈之间区域初始振幅E0的近似值假设n+1个控制点的坐标n次B样条曲线段的参数表达式为式中，Fi,n(t)为n次B样条基函数，其形式为式中，如图4所示，B样条插值后的数据矩阵增大，需要通过二分法搜索反演目标函数的最优正则化因子，提高目标函数的求解速度。首先给定初始值λ0和搜索步长Δλ，利用线性搜索的方法自动选取最优正则化因子区间[λn-1,λn+1]，计算λn＝λ0-(n-1)Δλ对应的总体目标函数Φ和数据目标函数Φd，为保证数据的拟合误差逐阶减小，定义数据拟合目标项系数β通常选取[0,1]区间内较大的值。在区间[λn-1,λn+1]内判断的符号，若大于零，则重新搜索最优正则化因子区间；若小于零，则确定[λn-1,λn+1]为满足条件的单谷区间。确定含根区间后，求区间的中点λn＝1/(2(λn-1+λn+1))，计算λn所对应的总体目标函数Φ和数据目标函数Φd的值，直到λn满足则λn为本次迭代的最优正则化参数；如图5所示，为基于B样条插值的二维地面核磁共振反演方法流程图一种基于B样条插值的二维地面核磁共振反演方法，包括以下步骤：a、如图1所示，在水平地面沿东西方向铺设长方形发射线圈1，在发射线圈1内部等距离铺设第一接收线圈2、第二接收线圈3、第三接收线圈4和第四接收线圈5，采用商业化的有限元软件COMSOL，建立三维可视化模型，令线圈中通入单位电流脉冲，利用自适应网格剖分技术精确计算地下任意一点的三维矢量磁场，即式中，B为三维矢量磁场；Bx为沿x方向的矢量磁场；By为沿y方向的矢量磁场；Bz为沿z方向的矢量磁场；和为直角坐标系的方向向量；b、如图2所示，令地磁场方向为x轴，水平面上与地磁场方向垂直为y轴，垂直地面向下为z轴，计算感应磁场垂直于地磁场方向的分量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过程，直到反演数据误差小于设定误差，获得高精度的反演结果，输出并快速成像，基于B样条的二维地面核磁共振反演方法流程如图5所示。实施例1基于B样条插值的二维地面核磁共振反演方法是一种适用于非层状不均匀水体的高效率高分辨率探测，既能节省测量时间，又能可靠获得高精度平滑的反演结果的反演方法，该方法通过发射线圈1向地下发射频率为当地拉莫尔频率的交变电流，产生感应磁场。在感应磁场的作用下，地下水中氢质子产生能级跃迁，大量氢质子跃迁到高能级上。当撤去供电电流，这些高能级氢质子便逐渐回到低能级状态，释放出大量的具有拉莫尔频率的能量子，在第一接收线圈2、第二接收线圈3、第三接收线圈4和第四接收线圈5中感应出核磁共振信号。令地磁场方向为x轴，水平面上与地磁场方向垂直为y轴，垂直地面向下为z轴，计算感应磁场的垂直分量获得灵敏度核函数，借助该灵敏度核函数得到模型数据，使用均匀B样条插值方法对各线圈内的初始振幅进行插值计算，得到探测线圈之间区域初始振幅的近似值，再通过二分法搜索最优正则化因子，对目标函数进行高斯牛顿迭代法求解，得到高精度平滑的二维含水量反演结果。一种基于B样条插值的二维地面核磁共振反演方法，包括下列顺序和步骤：a、使用常规尺寸线圈探测模式，在水平地面沿东西方向铺设长方形发射线圈1，长100米，宽25米，在发射线圈1内每隔20米铺设第一接收线圈2、第二接收线圈3、第三接收线圈4和第四接收线圈5，接收线圈为边长10米的正方形线圈，采用商业化的有限元软件COMSOL，建立三维可视化模型，通过自适应网格剖分技术精确计算地下任意一点的三维矢量磁场，即式中，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过程，直到反演数据误差小于设定误差，获得高精度的反演结果，输出并快速成像。实施例2基于B样条插值的二维地面核磁共振反演方法是一种适用于非层状不均匀水体的高效率高分辨率探测，既能节省测量时间，又能可靠获得高精度平滑的反演结果的反演方法，该方法通过发射线圈1向地下发射频率为当地拉莫尔频率的交变电流，产生感应磁场。在感应磁场的作用下，地下水中氢质子产生能级跃迁，大量氢质子跃迁到高能级上。当撤去供电电流，这些高能级氢质子便逐渐回到低能级状态，释放出大量的具有拉莫尔频率的能量子，在第一接收线圈2、第二接收线圈3、第三接收线圈4和第四接收线圈5中感应出核磁共振信号。令地磁场方向为x轴，水平面上与地磁场方向垂直为y轴，垂直地面向下为z轴，计算感应磁场的垂直分量获得灵敏度核函数，借助该灵敏度核函数得到模型数据，使用均匀B样条插值方法对各线圈内的初始振幅进行插值计算，得到探测线圈之间区域初始振幅的近似值，再通过二分法搜索最优正则化因子，对目标函数进行高斯牛顿迭代法求解，得到高精度平滑的二维含水量反演结果。一种基于B样条插值的二维地面核磁共振反演方法，包括下列顺序和步骤：a、使用小尺寸线圈探测模式，在水平地面沿东西方向铺设长方形发射线圈1，长20米，宽5米，在发射线圈1内每隔5米铺设第一接收线圈2、第二接收线圈3、第三接收线圈4和第四接收线圈5，接收线圈为边长1米的正方形线圈，采用商业化的有限元软件COMSOL，建立三维可视化模型，通过自适应网格剖分技术精确计算地下任意一点的三维矢量磁场，即式中，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过程，直到反演数据误差小于设定误差，获得高精度的反演结果，输出并快速成像。