基于速度估计与伪影抑制的探地雷达高精度逆时偏移方法技术

本发明专利技术公开了一种基于速度估计与伪影抑制的探地雷达高精度逆时偏移方法，包括建立正演模型，获取时间域正演记录；加载正演记录，确定介电常数模型序列；计算反传电磁波场；利用互相关成像条件进行逆时偏移成像；计算清晰度评价函数值；输出清晰度评价函数最大值处索引及其对应的偏移成像剖面；构建TV去噪函数，求解TV去噪函数的最优化结果；得到最终的雷达偏移剖面。本发明专利技术通过自动聚焦技术，自动获取最适当的偏移介电常数模型，以减少逆时偏移过程中人为可视化检验与修正的次数；同时成功抑制了基于互相关成像条件的逆时偏移成像剖面中的伪影，重建轮廓，提高边缘锐利度，大大提升偏移后雷达剖面的分辨率与信噪比，从而利于地质解释。解释。解释。

【技术实现步骤摘要】
基于速度估计与伪影抑制的探地雷达高精度逆时偏移方法


[0001]本专利技术属于地球物理领域，具体涉及一种基于速度估计与伪影抑制的探地雷达高精度逆时偏移方法。

技术介绍

[0002]探地雷达(Ground Penetrating Radar，简称GPR)是一种将高频电磁波传输到地下各种介质中并记录反射信号，以此呈现地下介质结构分布情况的地球物理勘探方法；探地雷达具有快速、分辨率高、探测无损等优点，已普遍应用于建筑物结构调查、工程质量检测与勘察、矿产资源勘查、考古测量等领域。由于地下介质的分布及其组成复杂多样，呈现出不规则的随机分布(例如电磁波遇到物性分界面产生多次波、绕射波)，且在GPR测区中充斥着电磁干扰，直接导致GPR的接收信号夹杂着各种杂波和噪声，使得GPR原始剖面无法准确体现地下结构的真实情况。
[0003]为了获取精准的地下复杂介质的结构分布信息，实现反射波正确归位、绕射波准确收敛到其真实空间位置，偏移技术被提出，与有限差分偏移、Kirchhoff积分偏移、F
‑
K偏移等单程方法相比，基于双程波动方程的逆时偏移(Reverse Time Migration，简称RTM)有效利用了全波场信息，无需分离上行波和下行波，克服了单程波偏移方法受横向速度变化和倾角限制影响较大的不足，具有精度高、相位准确、可沿任意方向传播的优点，更加适用于复杂地质构造的成像。
[0004]目前逆时偏移技术在探地雷达领域中已逐渐成熟，该技术依赖于给定的初始速度模型，且偏移剖面中存在杂波干扰和伪像。为此，众多地球物理学者开展了一系列的深入研究，提出了介电常数法、金属反射法、共中心点道集法、Hough变换法、多点拟合法和自动聚焦等速度估计方法，有效提高了速度估计模型的准确性。但是，逆时偏移在工程实例中受制于需要先验信息、计算效率低、内存消耗高等影响，与大规模的实践推广仍存在一定的距离。

技术实现思路

[0005]本专利技术的目的在于提供一种基于速度估计与伪影抑制的探地雷达高精度逆时偏移方法，该方法能够自动获取偏移模型、有效抑制偏移成像的伪影，从而提高雷达剖面精度，利于GPR资料解释。
[0006]本专利技术提供的这种基于速度估计与伪影抑制的探地雷达高精度逆时偏移方法，包括如下步骤：
[0007]S1.建立正演模型，获取时间域正演记录；
[0008]S2.加载正演记录，确定偏移介电常数模型序列；
[0009]S3.计算反传电磁波场；
[0010]S4.利用互相关成像条件进行逆时偏移成像；
[0011]S5.计算清晰度评价函数值；
[0012]S6.重复步骤S3～S5直至遍历介电常数模型序列结束，得到最终的清晰度评价函数曲线；
[0013]S7.输出清晰度评价函数最大值处索引及其对应的偏移成像剖面；
[0014]S8.构建TV去噪函数，加载步骤S7中原始的偏移成像剖面；
[0015]S9.引入Bregman分裂方法，求解TV去噪函数的最优化结果；
[0016]S10.计算更新参数，完成TV去噪，得到最终的雷达偏移剖面。
[0017]所述的步骤S2，具体包括，由电磁场和电磁波的理论，二维情况下，探地雷达波的传播规律遵循麦克斯韦方程，具体为：
[0018][0019]其中，表示接收点处磁场强度在x轴方向上的分量；表示接收点处磁场强度在z轴方向上的分量；表示接收点处电场强度在y轴上的分量；ε表示介质的介电常数；μ表示介质的磁导率；σ表示介质的电导率；t表示时间；
[0020]将接收信号作为边界条件从最大时刻进行逆时外推，进行逆时外推的电磁波场和正向传播的电磁波场采用如下公式计算：
[0021][0022][0023]其中，x表示接收点位置；T表示时间窗口；表示逆时外推的电场强度；表示正向传播的电场强度；表示逆时外推的磁场强度；表示正向传播的磁场强度；t表示时间。
[0024]所述的步骤S3，具体包括，以Yee网格离散空间，采用时域有限差分方法进行数值模拟，将接收信号作为边界条件从最大时刻进行逆时外推，按照时间步长交替更新电场和磁场，电磁波场进行逆时外推的更新迭代，具体方程为：
[0025][0026][0027][0028][0029]其中，n表示当前时间步；(i,j)为空间网格点的索引；表示逆时外推的磁场强度在x方向上的分量，n
‑
1/2表示前一时间步；表示逆时外推的磁场强度在x方向上的分量，n+1/2表示下一时间步；m表示更新方程的左端的磁场的分量或电场的分量在网格节点处的空间位置；CQ(m)表示更新方程的左端磁场的分量或电场的分量在网络节点处空间位置的第一电性参数，CQ(m)＝Δt/μ；表示逆时外推的电场强度在y方向上的分量；表示逆时外推的电场强度在y方向上的分量；表示逆时外推的磁场强度在z方向上的分量；表示逆时外推的磁场强度在z方向上的分量；CB(m)表示更新方程的左端磁场的分量或电场的分量在网络节点处空间位置的第二电性参数，CB(m)＝2Δt/(2ε+σ
·
Δt)；CA(m)表示更新方程的左端磁场的分量或电场的分量在网络节点处空间位置的第三电性参数，同时CA(m)定义为：CA(m)＝(2ε
‑
σ
·
Δt)/(2ε+σ
·
Δt)；表示逆时外推的磁场强度在x方向上的分量；Δt为时间变化量；ε为介质的介电常数；μ为介质的磁导率；σ为介质的电导率。
[0030]所述的步骤S4，互相关成像条件具体为，基于最大相干性原理，利用正传波场到达某一个反射界面时与接收点的信号记录逆时延拓得到的反传波场具有最大相干性的特性，将所有激发点和接收点的波场叠加，从而获得地下空间分布的成像剖面，具体表达为：
[0031][0032]其中，T
max
表示最大时刻；表示地下空间的成像结果；表示源点电磁波场；表示接收点反传电磁波场。
[0033]所述的步骤S5，自动聚焦技术计算清晰度评价函数值，具体包括如下聚焦函数计算清晰度评价值：
[0034]A1.计算图像集中度函数COT；
[0035]A2.计算平均强度函数AIT；
[0036]A3.计算图像对比度函数ICBT。
[0037]所述的步骤A1，图像集中度函数COT为：
[0038][0039]其中，M表示采集道数；N表示每道的采样点数；(i,j)为空间网格点的索引；表示像素点(x
i
,z
j
)处的灰度强度值。
[0040]所述的步骤A2，平均强度函数AIT为：
[0041][0042]其中，M表示采集道数；N表示每道的采样点数；(i,j)为空间网格点的索引；表示像素点(x
i
,z
j
)处的灰度强度值。
[0043]所述的步骤A3，图像对比度函数ICBT为：
[0044][0045]其中，M表示采集道数；N表示每道的本文档来自技高网...

【技术保护点】

【技术特征摘要】
1.一种基于速度估计与伪影抑制的探地雷达高精度逆时偏移方法，其特征在于包括如下步骤：S1.建立正演模型，获取时间域正演记录；S2.加载正演记录，确定介电常数模型序列；S3.计算反传电磁波场；S4.利用互相关成像条件进行逆时偏移成像；S5.计算清晰度评价函数值；S6.重复步骤S3～S5直至遍历介电常数模型序列结束，得到最终的清晰度评价函数曲线；S7.输出清晰度评价函数最大值处索引及其对应的偏移成像剖面；S8.构建TV去噪函数，加载步骤S7中原始的偏移成像剖面；S9.引入Bregman分裂方法，求解TV去噪函数的最优化结果；S10.计算更新参数，完成TV去噪，得到最终的雷达偏移剖面。2.根据权利要求1所述的基于速度估计与伪影抑制的探地雷达高精度逆时偏移方法，其特征在于所述的步骤S2，具体包括，由电磁场和电磁波的理论，二维情况下，探地雷达波的传播规律遵循麦克斯韦方程，具体为：其中，表示接收点处磁场强度在x轴方向上的分量；表示接收点处磁场强度在z轴方向上的分量；表示接收点处电场强度在y轴上的分量；ε表示介质的介电常数；μ表示介质的磁导率；σ表示介质的电导率；t表示时间；将接收信号作为边界条件从最大时刻进行逆时外推，进行逆时外推的电磁波场和正向传播的电磁波场采用如下公式计算：传播的电磁波场采用如下公式计算：其中，x表示接收点位置；T表示时间窗口；表示逆时外推的电场强度；表示正向传播的电场强度；表示逆时外推的磁场强度；表示正向传播的磁场强度；t表示时间。3.根据权利要求2所述的基于速度估计与伪影抑制的探地雷达高精度逆时偏移方法，
其特征在于所述的步骤S3，具体包括，以Yee网格离散空间，采用时域有限差分方法进行数值模拟，将接收信号作为边界条件从最大时刻进行逆时外推，按照时间步长交替更新电场和磁场，电磁波场进行逆时外推的更新迭代，具体方程为：为：为：其中，n表示当前时间步；(i,j)为空间网格点的索引；表示逆时外推的磁场强度在x方向上的分量，n
‑
1/2表示前一时间步；表示逆时外推的磁场强度在x方向上的分量，n+1/2表示下一时间步；m表示更新方程的左端的磁场的分量或电场的分量在网格节点处的空间位置；CQ(m)表示更新方程的左端磁场的分量或电场的分量在网络节点处空间位置的第一电性参数，CQ(m)＝Δt/μ；表示逆时外推的电场强度在y方向上的分量；表示逆时外推的电场强度在y方向上的分量；表示逆时外推的磁场强度在z方向上的分量；表示逆时外推的磁场强度在z方向上的分量；CB(m)表示更新方程的左端磁场的分量或电场的分量在网络节点处空间位置的第二电性参数，CB(m)＝2Δt/(2ε+σ
·
Δt)；CA(m)表示更新方程的左端磁场的分量或电场的分量在网络节点处空间位置的第三电性参数，同时CA(m)定义为：CA(m)＝(2ε
‑
σ
·
Δt)/(2ε+σ
·
Δt)；表示逆时外推的磁场强度在x方向上的分量；Δt为时间变化量；ε为介质的介电常数；μ为介质的磁导率；σ为介质的电导率。4.根据权利要求3所述的基于速度估计与伪影抑制的探地雷达高精度逆时偏移方法，其特征在于所述的步骤S4，互相关成像条件具体为，基于最大相干性原理，利用正传波场到达某一个反射界面时与接收点的信号记录逆时延拓得到的反传波场具有最大相干性的特性，将所有激发点和接收点的波场叠加，从而获得地下空间分布的成像剖面，具体表达为：其中，T
max
表示最大时刻；表示地下空间的成像结果；表示源点电磁波场；表示接收点反传电磁波场。
5.根据权利要求4所述的基于速度估计与伪影抑制的探地雷达高精度逆时偏移方法，其特征在于所述的步骤S5，具体为采用自动聚焦技术计算清晰度评价函数值，包括如下聚焦函数计算清晰度评价值：A1.计算图像集中度函数COT；A2.计算平均强度函数AIT；A3.计算图像对比度函数ICBT。6.根据权利要求5所述的基于速度估计与伪影抑制的探地雷达高精度逆时偏移方法，其特征在于所述的步骤A1，图像集中度函数COT为：其中，M表示采集道数；N表示每道的采样点数；(i,j)为空间网格点的索引；表示像素点(x
i
,z
j
)处的灰度强度值。7.根据权利要求6所述的基于速度估计与伪影抑制的探地雷达高精度逆时偏移方法，其特征在于所述的步骤A2，平均强度函数AIT为：其中，M表示采集道数；N表示每道的采样点数；(i,j)为空间网格点的索引；表示像素点(x
i
,z
j
)处的灰度强度值。8.根据权利要求7所述的基于速度估计与伪影抑制的探地雷达高精度逆时偏移方法，其特征在于所述的步骤A3，图像对比度函数ICBT为：其中，M表示采集道数；N表示每道的采样点数；(i,j)为空间网格点的索引；表示像素点(x
i
,z
j
)处的灰度强度值。9.根据权利要求8所述的基于速度估计与伪影抑制的探地雷达高精度逆时偏移方法，其特征在于所述的步骤S8，具体包括，采用L...

【专利技术属性】
技术研发人员：李广场，冯德山，李婷，魏名地，王珣，
申请(专利权)人：中南大学，
类型：发明
国别省市：