CN110082832A - 一种地面磁共振和探地雷达数据联合成像方法 - Google Patents

Abstract

本发明属于地面核磁共振数据处理领域及反演解释领域，为一种地面磁共振和探地雷达数据联合成像方法。该方法包括：在同一实验区域分别采用地面磁共振方法和探地雷达方法获得两类地球物理数据，并利用Gazdag相移算法，获得探地雷达数据的初步反射结构成像结果；根据探地雷达成像结果得到的结构信息，设计初始化地面磁共振二维反演剖分网格，进行反演成像，得到初步的含水量和弛豫时间分布图像；进一步根据磁共振得到的含水量和弛豫时间，重新计算每层介质的波速，再次对探地雷达数据进行成像；反复进行多次上述探地雷达数据和地面磁共振数据联合反演成像过程，直到满足精度要求，最终获得一致的地下水成像结果。经试验，本方法能够提高浅层地下水的成像精度。

Description

一种地面磁共振和探地雷达数据联合成像方法

技术领域

本发明属于地球物理数据处理及反演解释领域，具体来讲为一种地面磁共振和探地雷达数据联合成像方法。

背景技术

在进行浅层地下水探测时，探地雷达方法能够获得含水层的结构和边界信息，而地面磁共振方法能够获得含水层的含水量和孔隙分布信息。但是两种方法都存在自身的局限性，地面磁共振方法对含水层边界的反演分辨率较低，探地雷达方法无法直接获得准确水文参数。因此，研究地面磁共振和探地雷达数据联合成像具有重要意义。

CN104297807A公开了一种下灾害水源探测磁共振成像装置及探测和成像方法，利用梯度发射机，在探测深度层面产生稳定频率梯度场。通过增加主动梯度场，实现了在地下空间采用磁共振成像对灾害水源直接成像，提高成像分辨率。但是该方法反演结果仍然存在含水结构边界模糊的问题。

CN108387896A公开了一种基于探地雷达回波数据的自动汇聚成像方法，基于余弦距离加权的反向投影算法，采用阂值分析的图像分割技术，对回波图像进行会聚成像。该方法在复杂的城市环境下，可以快速、自动地对探地雷达回波信号进行汇聚成像，满足实际和应用需求，但该方法无法直接获得含水量信息。

CN104537714A公开了一种磁共振与瞬变电磁空间约束联合反演方法，基于连分式展开数值积分等价变换法与数字滤波法相结合的正演方法，构建了联合反演目标函数与迭代方程和三维大地模型。通过引入先验信息约束矩阵、粗糙度矩阵及层厚度约束矩阵等空间约束矩阵，实现了电阻率、层厚度、含水量、弛豫时间的光滑连续约束，解决了反演非唯一性问题。但是该方法反演结果同样存在含水结构边界分辨率低的问题。

发明内容

本发明所要解决的技术问题在于提供一种地面磁共振和探地雷达数据联合成像方法，能够提高浅层地下水的成像精度。

本发明是这样实现的，

一种地面磁共振和探地雷达数据联合成像方法，该方法包括，在同一实验区域分别采用地面磁共振方法和探地雷达方法获得两类地球物理数据；

利用Gazdag相移算法，获得探地雷达数据的初步反射结构成像结果；

根据探地雷达成像结果得到的结构信息，设计初始化地面磁共振二维反演剖分网格，对地面磁共振测量数据进行反演成像，得到初步的含水量和弛豫时间分布图像；

根据磁共振得到的含水量和弛豫时间，重新计算每层介质的波速；

根据重新计算每层介质的波速再次对探地雷达数据进行成像；反复进行多次上述探地雷达数据和地面磁共振数据联合反演成像过程，直到满足精度要求。

进一步地，所述利用Gazdag相移算法，获得探地雷达数据的初步反射结构成像结果包括：

a)对探地雷达数据v_G在水平位置x和时间t方向进行二维傅里叶变换：

其中k_x为波数矢量在水平方向的分量水平位置；z₀为初始深度，为频率；通过相移操作，计算深度z₁＝z₀+Δz处的波场：

其中v(z)为地下介质的波速；

b)对波数域数据进行傅里叶逆变换，并令t＝0：

d)对k_x进行傅里叶逆变换，得到成像值：

进一步地，所述设计初始化地面磁共振二维反演剖分网格包括：根据已知的含水模型边界信息进行三角剖分，达到网格完全匹配地质构造边界。

进一步地，对地面磁共振测量数据进行反演成像，得到初步的含水量和弛豫时间分布图像包括：根据探地雷达地质结构约束的二维剖分网格，对地面磁共振测量数据v_M进行二维QT反演，二维QT反演目标函数如下：

式中D为地面磁共振测量数据v_M的权值，K_2D为二维核函数，C为一阶平滑度矩阵，λ为正则化因子，m为含水量和弛豫时间参数组合成的向量，利用高斯牛顿迭代法求解式(5)，得到含水量和弛豫时间。

进一步地，根据磁共振得到的含水量和弛豫时间，重新计算每层介质的波速包括：

根据水文参数和介电常数的关系及电磁波速度的转换关系式(6)，得到岩石颗粒的介电常数ε_s：

其中ε_w，ε_b水和空气的介电常数，w₁为含水量，φ为孔隙度；

利用岩石颗粒的介电常数ε_s，重新计算每层介质的波速v′(z)。

本发明与现有技术相比，有益效果在于：本发明提供一种地面磁共振和探地雷达数据联合成像方法，旨在基于上述两种地球物理探测结果结构相似和物性参数耦合关系，提出一种新的地面磁共振和探地雷达联合成像机制，通过结合两种方法的优势，利用互相约束的策略，实现基于结构约束和参数耦合的浅层地下水精细成像和水文信息解释。

附图说明

图1为地面磁共振和探地雷达数据联合成像方法的流程图；

图2基于探地雷达成像结构信息的磁共振反演网格二维剖分结果；

图3基于探地雷达成像结构信息的磁共振反演结果A为含水量与深度，B为弛豫时间与深度；

图4地面磁共振和探地雷达数据联合成像结果A为含水量与深度，B为弛豫时间与深度。

具体实施方式

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

参见图1所示，一种地面磁共振和探地雷达数据联合成像方法，包括在同一实验区域分别采用地面磁共振方法和探地雷达方法获得两类地球物理数据；利用Gazdag相移算法，获得探地雷达数据的初步反射结构成像结果；根据探地雷达成像结果得到的结构信息，设计初始化地面磁共振二维反演剖分网格，对地面磁共振测量数据进行反演成像，得到初步的含水量和弛豫时间分布图像；根据磁共振得到的含水量和弛豫时间，重新计算每层介质的波速；根据重新计算每层介质的波速再次对探地雷达数据进行成像；反复进行多次上述探地雷达数据和地面磁共振数据联合反演成像过程。

具体包括：步骤1：在野外实验地点，首先利用一发多收阵列线圈测量方式，得到二维地面磁共振数据v_M，在一实施例中地面磁共振发射线圈和接收线圈匝数均为1匝，发射尺寸为60m×20m，接收线圈尺寸为20m×20m，7个接收线圈基于半覆盖型铺设。其次在同一测量区域，利用共中心点测量方式，得探地雷达数据v_G，在一实施例中探地雷达测线与地面磁共振测线重合。通过Gazdag相移算法，获得探地雷达数据v_G的初步反射结构二维成像结果P₁，其反射面结果见图2黑色虚线所示。

其中Gazdag相移算法实现方式如下：

c)对探地雷达数据v_G在水平位置x和时间t方向进行二维傅里叶变换：

其中k_x为波数矢量在水平方向的分量水平位置，z₀为初始深度，为频率。

d)通过相移操作，计算深度z₁＝z₀+Δz处的波场：

其中v(z)为地下介质的波速。

f)对波数域数据进行傅里叶逆变换，并令t＝0：

d)对k_x进行傅里叶逆变换，即得到成像值：

步骤2：根据步骤1得到的二维成像结果P₁，设计初始化地面磁共振二维反演剖分网格，即在非均匀三角剖分网格的基础上，根据已知的含水模型边界信息进行三角剖分，达到网格完全匹配地质构造边界的效果，见图2所示。根据探地雷达地质结构约束的二维剖分网格，对地面磁共振测量数据v_M进行反演，获得初步的含水量和弛豫时间分布图像(w₁,t₁)。

其中，地面磁共振二维QT反演目标函数如下：

式中D为实测数据v_M的权值，K_2D为二维核函数，C为一阶平滑度矩阵，λ为正则化因子，m为含水量和弛豫时间参数组合成的向量。利用高斯牛顿迭代法求解式(5)，即得到含水量和弛豫时间，见图3所示。一实施例中正则化因子λ为100。

步骤3：根据步骤2得到的成像结果(w₁,t₁)，基于水文参数和介电常数的关系及电磁波速度的转换关系，得到岩石颗粒的介电常数ε_s：

其中ε_w，ε_b水和空气和介电常数，w₁为含水量，φ为孔隙度。

利用地面磁共振反演结果得到的岩石颗粒的介电常数ε_s，重新计算每层介质的波速v′(z)，并将v′(z)引入Gazdag相移方法中计算每个深度对应的波场(步骤1中b))，结合地下介质重新对探地雷达测量数据v_G进行成像，得到P₂。

步骤4：根据步骤3得到的成像结果P₂，对二维大地模型重新剖分，使网格与P₂得到的地质构造边界完全匹配，再次对地面磁共振数据v_M进行反演获得含水量和弛豫时间成像结果(w₂,t₂)，见图3所示。

步骤5：通过步骤3和步骤4的互补机制，反复利用地面磁共振反演结果(w_n,t_n)和探地雷达成像结果P_n，重复进行多次探地雷达数据和地面磁共振数据联合反演成像，直到测量数据和成像结果正演数据间的拟合误差小于1％，满足精度要求，在一实施例中，经过3次联合成像，最终获得一致的地下水成像结果见图4。对比图3，可以看出图4中含水层的边界信息明显，能够通过弛豫时间有效区分，因此证明了地面磁共振和探地雷达数据联合成像结果能够提供更多有效信息，有利于浅层地下水的精细探测。

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

Claims (5)

1.一种地面磁共振和探地雷达数据联合成像方法，其特征在于，该方法包括，在同一实验区域分别采用地面磁共振方法和探地雷达方法获得两类地球物理数据；

利用Gazdag相移算法，获得探地雷达数据的初步反射结构成像结果；

根据探地雷达成像结果得到的结构信息，设计初始化地面磁共振二维反演剖分网格，对地面磁共振测量数据进行反演成像，得到初步的含水量和弛豫时间分布图像；

根据磁共振得到的含水量和弛豫时间，重新计算每层介质的波速；

根据重新计算每层介质的波速再次对探地雷达数据进行成像；反复进行多次上述探地雷达数据和地面磁共振数据联合反演成像过程，直到满足精度要求。

2.按照权利要求1所述的方法，其特征在于，其特征在于，所述利用Gazdag相移算法，获得探地雷达数据的初步反射结构成像结果包括：

a)对探地雷达数据v_G在水平位置x和时间t方向进行二维傅里叶变换：

其中k_x为波数矢量在水平方向的分量水平位置；z₀为初始深度，为频率；通过相移操作，计算深度z₁＝z₀+Δz处的波场：

其中v(z)为地下介质的波速；

b)对波数域数据进行傅里叶逆变换，并令t＝0：

d)对k_x进行傅里叶逆变换，得到成像值：

3.按照权利要求1所述的方法，其特征在于，所述设计初始化地面磁共振二维反演剖分网格包括：根据已知的含水模型边界信息进行三角剖分，达到网格完全匹配地质构造边界。

4.按照权利要求1所述的方法，其特征在于，对地面磁共振测量数据进行反演成像，得到初步的含水量和弛豫时间分布图像包括：根据探地雷达地质结构约束的二维剖分网格，对地面磁共振测量数据v_M进行二维QT反演，二维QT反演目标函数如下：

式中D为地面磁共振测量数据v_M的权值，K_2D为二维核函数，C为一阶平滑度矩阵，λ为正则化因子，m为含水量和弛豫时间参数组合成的向量，利用高斯牛顿迭代法求解式(5)，得到含水量和弛豫时间。

5.按照权利要求1所述的方法，其特征在于，根据磁共振得到的含水量和弛豫时间，重新计算每层介质的波速包括：

根据水文参数和介电常数的关系及电磁波速度的转换关系式(6)，得到岩石颗粒的介电常数ε_s：

其中ε_w，ε_b水和空气的介电常数，w₁为含水量，φ为孔隙度；

利用岩石颗粒的介电常数ε_s，重新计算每层介质的波速v′(z)。