基于非均匀输入参数网格的井中雷达最小二乘反演方法

摘要

本发明公开了一种基于非均匀输入参数网格的井中雷达最小二乘反演方法，其特征在于，包括如下步骤：S1、通过偏移成像确定目标形状；S2、通过k‑means算法进行初始模型网格划分；S3、将网格划分完毕的初始模型输入反演代码进行反演计算；S4、将模型参数更新矩阵进行迭代，迭代终止输出反演结果。本发明将脉冲信号作为未知的反演方法，借助于正演模拟生成脉冲信号，通过搭建非均匀模型网格进行了计算量的优化，通过并行计算方案加快计算速度，获得仿真模型的目标反演，提高了计算速度以及计算精度，优化反演收敛情况。

权利要求

1.基于非均匀输入参数网格的井中雷达最小二乘反演方法，其特征在于，包括如下步骤：S1、通过偏移成像确定目标形状：包括以下步骤：S11、对采集获取的雷达B扫数据进行信号预处理；

S12、对预处理完毕的B扫数据进行偏移成像，获得包含目标形状的初始模型成像结果图像数组；

S2、通过k‑means算法进行初始模型网格划分；

S3、将网格划分完毕的初始模型输入反演代码进行反演计算，反演计算具体包括如下步骤：0

S31、模型初始电磁参数m为初始模型介电常数矩阵，Δm为模型电磁参数改变量，目标函数近似为：将求解反演目标函数的问题转化为求解 函数的反演问题，根据最优化原理，目标是将T T函数达到极小，其对模型电磁参数改变量偏导为0，反演方程组简写如下：J·JΔm＝J b；

0 0 T

J为E(m)在m处的雅可比矩阵，即E(m)对于模型矩阵m的偏导矩阵，b为 J 表示矩阵的转置；

步骤S32、通过差分扰动法求雅可比矩阵J的解，公式如下：T T

步骤S33、对雅可比矩阵J·J加上阻尼系数，改善矩阵的病态特性，将J·J记为A，阻尼‑系数也即 反演方程组 k的解为Δm＝(A+λI)1

k；

S4、将模型参数更新矩阵进行迭代，迭代终止输出反演结果。

2.根据权利要求1所述的基于非均匀输入参数网格的井中雷达最小二乘反演方法，其特征在于：所述步骤S11中的信号预处理，具体是指通过平均滤波法去除直达波和直流偏移，通过带通滤波算法去除带外噪声干扰。

3.根据权利要求1所述的基于非均匀输入参数网格的井中雷达最小二乘反演方法，其特征在于：所述步骤S12中的偏移方法为F‑K偏移，具体是总结每个惠更斯二次能源产生的能量，并将其映射到发电点，通过衍射求和将每个衍射双曲线收缩到其顶点。

4.根据权利要求1所述的基于非均匀输入参数网格的井中雷达最小二乘反演方法，其特征在于：所述步骤S2中的通过k‑means算法进行初始模型网格划分是指将步骤S12生成的成像结果图像数组进行非均匀网格划分，得到初始模型的介电常数，具体步骤包括：对网格划分之后的矩阵中所有观测值分别划分到k集群，每个观测值属于与最近的平均相近的簇作为原型的集群，使用集群的平均值代替远观测值，使用欧几里得距离用作度量，方差用作聚类散布的量度；将大面积的均匀背景进行合并，从而将细节部分凸显出来，将数量庞大的差分网格计算转化为非均匀的网格计算。

5.根据权利要求1所述的基于非均匀输入参数网格的井中雷达最小二乘反演方法，其特征在于：所述的步骤S4中迭代所使用的并行计算矩阵方法为多线程计算方法。

说明书

基于非均匀输入参数网格的井中雷达最小二乘反演方法

技术领域

[0001] 本发明涉及井中雷达波形反演成像技术领域，具体涉及一种基于非均匀输入参数网格的井中雷达最小二乘反演方法。

背景技术

[0002] 井中雷达将收发天线、发射源、接收采集以及通信等雷达部件置于井中，从而进行井中探测。井中雷达的雷达波发射频率通常工作在较低的雷达频率(10MHz～1GHz)范围，采用高频脉冲信号，从而达到宽频带、高分辨率的目标。井中雷达不同于常规的地表探地雷达，由于其特殊的结构形式，可以下放的钻孔中，从而对钻孔周围地质环境进行探测，因此井中雷达有其独特的优势。首先，相较于地表探地雷达，井中雷达能够放置于钻孔内，探测范围不仅取决于垂直于井壁的径向深度，同时还包括井的轴向深度，因此探测的方式不同，并且探测范围更大；其次，雷达置于井中，由于地下电磁环境十分干净，受到的地面电磁干扰更小，从而能够更精确的反应地下异常目标。

[0003] 井中雷达数据可视化解释的一个重要方案就是雷达的波形反演，反演是实现资料从定性分析到定量解释的可行手段，其原理主要是以采集的井中雷达数据为基础，通过分析回波数据中波形的振幅、相位、频率等信息来确定或者定位目标区域介电参数分布情况，同雷达成像相似，雷达波形反演也需要对雷达信号的预处理，但是对于雷达的波形反演，更需要保证雷达信号能够到达一个非常高信噪比程度，相较于成像来说，对雷达信号质量的要求更加苛刻。

[0004] 在过去的几年中，由于计算能力的不断增强，对雷达数据进行全波反演已经成为解决地球物理反演问题的新解决方案，探地雷达全波反演是基于传统目标函数的时域全波形反演，利用井中雷达波形的全部信息进行反演。利用全部信息反演,可以获得高分辨率的反演结果，因此亟需一种能解决传统全波形反演非线性问题的方法。

发明内容

[0005] 为解决现有技术中存在的问题，本发明提供了一种基于非均匀输入参数网格的井中雷达最小二乘反演方法，将脉冲信号作为未知的反演方法，借助于正演模拟生成脉冲信号，通过搭建非均匀模型网格进行了计算量的优化，通过并行计算方案加快计算速度，获得仿真模型的目标反演，解决了上述背景技术中提到的问题。

[0006] 为实现上述目的，本发明提供如下技术方案：基于非均匀输入参数网格的井中雷达最小二乘反演方法，包括如下步骤：

[0007] S1、通过偏移成像确定目标形状；

[0008] S2、通过k‑means算法进行初始模型网格划分；

[0009] S3、将网格划分完毕的初始模型输入反演代码进行反演计算；

[0010] S4、将模型参数更新矩阵进行迭代，迭代终止输出反演结果。

[0011] 优选的，所述步骤S1通过偏移成像确定目标形状具体包括以下步骤：

[0012] S11、对采集获取的雷达B扫数据进行信号预处理；

[0013] S12、对预处理完毕的B扫数据进行偏移成像，获得包含目标形状的初始模型成像结果图像数组。

[0014] 优选的，所述步骤S11中的信号预处理，具体是指通过平均滤波法去除直达波和直流偏移，通过带通滤波算法去除带外噪声干扰。

[0015] 优选的，所述步骤S12中的偏移方法为F‑K偏移，具体是总结每个惠更斯二次能源产生的能量，并将其映射到发电点，通过衍射求和将每个衍射双曲线收缩到其顶点。

[0016] 优选的，所述步骤S2中的通过k‑means算法进行初始模型网格划分是指将步骤S12生成的成像结果图像数组进行非均匀网格划分，得到初始模型的介电常数，具体步骤包括：对网格划分之后的矩阵中所有观测值分别划分到k集群，每个观测值属于与最近的平均相近的簇作为原型的集群，使用集群的平均值代替远观测值，使用欧几里得距离用作度量，方差用作聚类散布的量度；将大面积的均匀背景进行合并，从而将细节部分凸显出来，将数量庞大的差分网格计算转化为非均匀的网格计算。

[0017] 优选的，所述步骤S3中的反演计算具体包括如下步骤：

[0018] S31、模型初始电磁参数m0为初始模型介电常数矩阵，Δm为模型电磁参数改变量，目标函数近似为：

[0019]

[0020] 将求解反演目标函数的问题转化为求解 函数的反演问题，根据最优化原理，目标T是将 函数达到极小，其对模型电磁参数改变量偏导为0，反演方程组简写如下：J·J Δm＝T 0 0
Jb；J为E(m)在m 处的雅可比矩阵，即E(m)对于模型矩阵m 的偏导矩阵，b为T
J表示矩阵的转置；

[0021] 步骤S32、通过差分扰动法求雅可比矩阵J的解，公式如下：

[0022]

[0023] 步骤S33、对雅可比矩阵J·JT加上阻尼系数，改善矩阵的病态特性，将J·JT记为A，阻尼系数也即 反演方程组 的解为Δm＝‑1
(A+λI) k；

[0024] 优选的，所述步骤S4将模型参数更新矩阵，参数更新为：mk＝mk‑1+Δm，输入到步骤0
S3中，作为新的m 进行迭代，直到迭代终止条件 当ηk<0.01时迭代停
止输出反演结果。

[0025] 优选的，所述的步骤S4中迭代所使用的并行计算矩阵方法为多线程计算方法。

[0026] 本发明的有益效果是：

[0027] 1)本发明适用于井中雷达数据反演，能够对井中雷达数据进行地层的电磁参数的计算，通过成像及网格划分，其反演输入参数矩阵为非均匀网格，对于小目标、稀疏矩阵特征的模型具有更好的，能够有效提高计算速度以及计算精度，本发明使用阻尼因子修正反演方程组，优化了雅可比矩阵的病态特征，优化反演收敛情况。

[0028] 2)本发明所提出的井中雷达全波反演方法是基于最小二乘方法，结合成像算法对反演初始模型进行构建，通过奇异值分解法计算、优化了反演方程组，使用经典的k‑means算法进行网格分割，最终形成了全波反演的整个算法框架。

附图说明

[0029] 图1为本发明方法流程示意图；

[0030] 图2为本发明成像效果示意图，图2(a)为B扫数据信号预处理后结果图，图2(b)为偏移成像结果图；

[0031] 图3为本发明模型与网格划分结果示意图，图3(a)为仿真模型示意图，图3(b)为k‑means聚类后的非均匀网格划分结果示意图；

[0032] 图4为本发明反演结果示意图；图4(a)为模型结果二维示意图，图4(b)为模型反演结果示意图。

具体实施方式

[0033] 下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

[0034] 实施例

[0035] 请参阅图1‑4，本发明提供一种技术方案：一种基于非均匀输入参数网格的探地雷达最小二乘反演方法，分为初始模型搭建与反演两个部分。初始模型搭建部分是通过井中雷达到模型的回波数据，通过分析回波数据进行初始模型的搭建。反演部分是将初始模型作为井中雷达反演的迭代入口，对井中雷达模型进行反演迭代，最终生成反演结果，包含目标的电磁参数矩阵，本发明方法流程如图1所示，包括下列步骤：

[0036] 步骤1：通过偏移成像确定目标形状

[0037] (1)对采集获取的雷达B扫数据进行信号预处理，去除直达波、直流偏移、噪声等干扰，如图2(a)所示。

[0038] (2)对预处理完毕的B扫数据进行偏移成像，获得包含目标形状的初始模型成像结果图像数组，偏移成像结果如图2(b)所示。偏移方法为F‑K偏移，即总结每个惠更斯二次能源产生的能量，并将其映射到发电点，通过衍射求和将每个衍射双曲线收缩到其顶点。

[0039] 步骤2：通过k‑means算法进行反演初始模型网格划分

[0040] 将步骤1(2)生成的成像结果图像数组进行非均匀网格划分，得到初始模型的介电常数，具体步骤包括：对网格划分之后的矩阵中所有观测值分别划分到k集群，每个观测值属于与最近的平均相近的簇作为原型的集群，使用集群的平均值代替远观测值，使用欧几里得距离用作度量，方差用作聚类散布的量度；将大面积的均匀背景进行合并，从而将细节部分凸显出来，将数量庞大的差分网格计算转化为非均匀的网格计算，如图3所示，图3(a)为模型示意图，图3(b)为k‑means聚类后的非均匀网格划分结果。

[0041] 步骤3：将网格划分完毕的初始模型输入反演代码进行反演计算

[0042] 1)模型初始电磁参数m0为初始模型介电常数矩阵，包括介电常数以及电导率、磁导率等，Δm为模型电磁参数改变量，也即迭代步长，目标函数近似为：

[0043]

[0044] 将求解反演目标函数的问题转化为求解 函数的反演问题，根据最优化原理，目标T是将 函数达到极小，其对模型电磁参数改变量偏导为0，反演方程组简写如下：J·J Δm＝T 0 0
Jb；J为E(m)在m处的雅可比矩阵，即E(m)对于模型矩阵m的偏导矩阵，b为
T
J表示矩阵的转置。

[0045] 2)通过差分扰动法求雅可比矩阵J的解，公式如下：

[0046]

[0047] 3)对雅可比矩阵J·JT加上阻尼系数，改善矩阵的病态特性，将J·JT记为A，阻尼系‑数也即 反演方程组 的解为Δm＝(A+λI)
1
k。

[0048] 步骤4：将模型参数更新矩阵进行迭代，迭代终止输出反演结果。

[0049] 将模型参数更新矩阵，参数更新为：mk＝mk‑1+Δm，输入到步骤3中，作为新的m0进行迭代，直到迭代终止条件 当ηk<0.01时迭代停止输出反演结果。其中，由于雅可比行列式的计算设计到大量的并行计算，使用线程池达到并发执行雅可比行列式计算，为多线程计算方法，反演结果如图4所示，图4(a)为模型二维示意图，目标为白色部分，介电常数为1左右，其他部分为地层介电常数为7左右，图4(b)为提取trace number＝
15处的模型，点划线为初始模型，点线为真实模型，实线为反演计算结果。

[0050] 实施例中，仿真原理为FDTD，仿真软件为GprMax。仿真模型如图3(a)所示，其中地层相对介电常数为εr＝7，电导率为σ＝1e‑4S/m，目标模型为三角形目标建立二维FDTD仿真模型，目标填充物为空气(介电常数为εr＝1，电导率σ＝1e‑4S/m)，距离井眼距离1.3m，发射天线和接收天线间距0.2m，发射波形为180MHz零阶高斯脉冲信号，时窗长度为6e‑8s，采样率为8.5GHz。

[0051] 成像结果如图2(b)所示，能够较直观地发现目标模型位置、大小等，通过网格分割，分割结果如图3(b)，能够看到，目标模型被分割成为非均匀的网格，以此作为初始模型，作为反演入口进行模型优化。

[0052] 经过反演计算获得的模型结果二维图如图4(a)所示，通过反演计算，能够有效计算出地层介电常数为7左右，目标模型节点常数为1左右，与实际模型接近。截取模型第15行介电常数分部情况，如图4(b)所示。其中，点划线为初始模型，目标区域相对介电常数介于2至17，地层部分介电常数为3；点线为真实模型，反演计算结果为实线线条。

[0053] 本发明所提出的井中雷达全波反演方法是基于最小二乘方法，结合成像算法对反演初始模型进行构建，通过奇异值分解法计算、优化了反演方程组，使用经典的k‑means算法进行网格分割，最终形成了全波反演的整个算法框架。本发明适用于井中雷达数据反演，能够对井中雷达数据进行地层的电磁参数的计算，通过成像及网格划分，其反演输入参数矩阵为非均匀网格，对于小目标、稀疏矩阵特征的模型具有更好的，能够有效提高计算速度以及计算精度，本发明使用阻尼因子修正反演方程组，优化了雅可比矩阵的病态特征，优化反演收敛情况。

[0054] 尽管参照前述实施例对本发明进行了详细的说明，对于本领域的技术人员来说，其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换,凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。