CN113156427A - 探地雷达数据的反演方法和装置 - Google Patents

探地雷达数据的反演方法和装置 Download PDF

Info

Publication number
CN113156427A
CN113156427A CN202110502385.1A CN202110502385A CN113156427A CN 113156427 A CN113156427 A CN 113156427A CN 202110502385 A CN202110502385 A CN 202110502385A CN 113156427 A CN113156427 A CN 113156427A
Authority
CN
China
Prior art keywords
layer
data
wave
amplitude
interface
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
CN202110502385.1A
Other languages
English (en)
Inventor
刘四新
薛美琪
鹿琪
李宏卿
邓莉
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Jilin University
Original Assignee
Jilin University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Jilin University filed Critical Jilin University
Priority to CN202110502385.1A priority Critical patent/CN113156427A/zh
Publication of CN113156427A publication Critical patent/CN113156427A/zh
Pending legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
    • G01S13/88—Radar or analogous systems specially adapted for specific applications
    • G01S13/885—Radar or analogous systems specially adapted for specific applications for ground probing

Landscapes

  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Electromagnetism (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • General Physics & Mathematics (AREA)
  • Geophysics And Detection Of Objects (AREA)
  • Radar Systems Or Details Thereof (AREA)

Abstract

本申请属于地质探测技术领域，具体涉及一种探地雷达数据的反演方法和装置。探地雷达数据的反演方法包括：获取探地雷达采用剖面法对待测区域进行探测扫描得到的实际测量数据；对实际测量数据进行数据预处理，得到待反演数据；从待反演数据中提取各个反射界面的反射波振幅和反射波的双程旅行时；将偏移距离、入射波振幅、预先确定的电磁波在第一层介质中的传播速度、各反射界面的反射波振幅和反射波的双程旅行时作为输入数据，迭代计算各层的厚度和反射波传播速度；基于各层的速度得到各层介质的相对介电常数。通过本申请的方法能够快速地采集探地雷达数据并实现准确、高效地探地雷达数据反演。

Description

探地雷达数据的反演方法和装置
技术领域
本申请属于地质探测技术领域，具体涉及一种探地雷达数据的反演方法和装置。
背景技术
探地雷达发射天线发射高频电磁波，电磁波在向下传播时，如果遇到具有电性差异的物体，将会发生反射、折射和散射，接收天线接收到反射信号并放大为接收信号，即A-Scan。随着天线在侧线上移动，多个测点的A-Scan构成雷达剖面图并通过计算机系统实时显示在屏幕上，即B-Scan。通过处理分析雷达剖面上的振幅、波形、双程旅行时等，可以判断测线下方目标的位置、几何形态、地下结构分层情况与各层介质的电性参数等重要信息。
通过对层状结构的反分析,可以得到每层结构的相对介电常数和层厚，满足工业应用中对层状构造分析(如路面厚度检测和覆土层探测等)的要求。然而，在目前的反演理论中，对分层结构的电磁波传播和介质特征模型的反演不够严谨，极大地限制了探地雷达无损技术的发展。因此，探地雷达的信号分析通常依赖于简化公式或根据相关解释经验调整参数。现有的研究成果大多只能保证结构层单层厚度的准确检测，对于多层结构层厚度和参数反演的研究以及其它道路指标的检测结果并不令人满意。
另一方面，用于探地雷达成像的速度场通常是由单个共中心道集分析得到的。在多偏移距GPR数据上进行的速度分析结果通常不非常准确，主要因为探地雷达的速度通常随深度而降低，这与地震速度不同。因此，经典速度分析中的一些假设，如将均方根速度同化为叠加速度，在应用于探地雷达数据时会导致结果不准确。且大多数商用探地雷达系统都只配备了单接收天线，因此对多偏移数据集的采集要求非常高、操作复杂。
综上，如何实现快速地采集数据并进行准确、高效地探地雷达数据反演，成为现有技术中亟待解决的问题。
发明内容
(一)要解决的技术问题
鉴于现有技术的上述缺点、不足，本申请提供一种探地雷达数据的反演方法和装置。
(二)技术方案
为达到上述目的，本申请采用如下技术方案：
第一方面，本申请提供一种探地雷达数据的反演方法，该方法包括：
S10、获取探地雷达采用剖面法对待测区域进行探测扫描得到的实际测量数据；
S20、对所述实际测量数据进行数据预处理，得到待反演数据；
S30、从所述待反演数据中提取各个反射界面的反射波振幅和反射波的双程旅行时；
S40、将所述探地雷达发射天线和接收天线间的偏移距离、所述探地雷达的入射波振幅、预先确定的电磁波在第一层介质中的传播速度、各反射界面的反射波振幅和反射波的双程旅行时作为输入数据，迭代计算各层的厚度和反射波传播速度；
S50、基于各层速度计算得到各层介质的相对介电常数。
可选地，所述待测区域为近水平层状地下结构，层数小于等于4，所述探地雷达对所述待测区域进行探测扫描时的探测扫描波长的确定方法，包括：
S01、获取预设扫描波长，将所述预设扫描波长作为初始扫描波长；
S02、对所述待测区域以所述初始扫描波长进行雷达探测，通过步骤S10-S40确定地层各层的厚度；
S03、对所述初始扫描波长进行递减，将得到的波长作为预设扫描波长；
S04、迭代执行步骤S01-03，直至满足预设的终止条件，并将最后一次迭代时的初始扫描波长作为探测扫描波长；预设的终止条件包括：
Figure BDA0003056910700000031
其中，λ表示探测扫描的波长，lmin表示地层的厚度最小值。
可选地，步骤S20包括：
对所述探地雷达测量数据进行直流滤波以去除数据中的直流分量，得到测量第一数据；
对所述测量第一数据进行零时刻校正，得到测量第二数据；
对所述测量第二数据进行反增益处理以恢复所述探测雷达反射波的振幅，得到测量第三数据；
对所述测量第三数据进行球面扩散振幅补偿得到测量第四数据；
对所述测量第四数据进行带通滤波，过滤掉信号中的高频噪声信号和低频振荡信号，得到待反演数据。
可选地，选取雷达直达波信号的负峰值点时刻作为零时刻点，对所述测量第一数据进行零时刻校正，得到测量第二数据。
可选地，当选取直达波信号的负峰值点时刻作为零时刻点时，对反演得到地层深度进行修正，修正公式为：
d'＝d+v×(T3-T2)
其中，d'表示修正后的地层深度，d表示反演得到的地层深度，v表示电磁波在每层介质中传播的平均速度，T3表示直达波信号的负峰值点时刻，T2表示直达波信号起跳点与负峰值点对应时刻的中点时刻。
可选地，对所述测量第三数据进行球面扩散振幅补偿得到测量第四数据的计算公式为：
A0＝Avteβt
其中，A0表示测量第四数据，为探地雷达发射天线发射出的没有经过扩散与介质吸收影响的真实振幅值；A表示测量第三数据；v表示电磁波在每层介质中传播的平均速度，t表示接收的反射振幅的双程旅行时，β表示每层介质的吸收系数。
可选地，对所述测量第四数据进行带通滤波时，采用正弦函数的平方作为镶边函数的滤波器，所述滤波器为：
Figure BDA0003056910700000041
其中，H(f)表示待反演数据，f1表示低通截止频率，f2表示低通频率，f3表示高通频率，f4表示高通截止频率。
可选地，迭代计算各层的厚度和反射波传播速度，包括：
S41、基于所述偏移距离、电磁波在第一层介质中的传播速度和第一层介质中的双程旅行旅行时，计算得到第一层介质的厚度与电磁波在第一层的入射角；
Figure BDA0003056910700000042
Figure BDA0003056910700000043
其中，x表示偏移距离，v1表示电磁波在第一层介质中的传播速度，TWT1表示第一介质中的双程旅行时，h1表示第一层介质的厚度，θ1表示电磁波在第一层介质中的入射角；
S42、基于所述入射波振幅和第一层反射波振幅，计算出第一个反射界面的反射系数，计算公式为：
Figure BDA0003056910700000044
其中，Ai1表示入射波振幅，As1表示第一个反射界面的反射波振幅，R1表示第一个反射界面的反射系数；
S43、令n＝2；
S44、基于第n-1个界面的反射系数、第n-1层的反射波传播速度、第n-1个界面的电磁波的入射角，由Snell方程计算得到第n层的反射波传播速度；
S45、基于前n-1层的厚度、前n层的反射波传播速度和第n层界面反射波的双程旅行时，根据雷达波反射路线与介质层厚度、所述偏移距离的几何关系计算得到第n层的厚度；
S46、基于所述偏移距离、前n层的反射波传播速度、前n层的厚度，计算得到得到沿第n次反射波路径各界面的电磁波的入射角，计算公式为：
Figure BDA0003056910700000051
其中，k表示当前界面的层数，k＝1,2,…,n；i表示界面层数；
S47、基于沿第n次反射波路径各界面的电磁波的入射角，计算前n-1层界面的透射系数，计算公式为：
Figure BDA0003056910700000052
其中，θk+1、θk分别表示电磁波在第k层和第k+1层介质中的入射角；
S48、基于提取的各反射界面的反射波振幅、前n-1层界面的透射系数，计算得到第n个界面的反射系数，计算公式为：
Figure BDA0003056910700000053
其中，Asn表示从所述待反演数据中提取的第n个反射界面的反射波振幅，i表示界面的层数，n表示界面的总层数，Ti表示第i层界面的透射系数；
S49、令n＝n+1，迭代执行步骤S44-S48，直至n＝4。
可选地，步骤S50包括：
S51、基于各层的速度计算得到各层介质的相对介电常数，计算公式为：
Figure BDA0003056910700000061
其中，c表示真空中电磁波的传播速度，v表示介质的电磁波传播速度，εr表示相对介电常数；
S52、采用平滑滑动平均对计算得到的相对介电常数进行处理，过滤掉异常值。
第二方面，本申请提供一种探地雷达数据的反演装置，该装置包括：
测量数据获取模块，用于获取探地雷达采用剖面法对待测区域进行探测扫描得到的实际测量数据；
数据预处理模块，用于对所述实际测量数据进行数据预处理，得到待反演数据；
反射波振幅提取模块，用于从所述待反演数据中提取各个反射界面的反射波振幅和反射波的双程旅行时；
迭代计算模块，用于将所述探地雷达发射天线和接收天线间的偏移距离、所述探地雷达的入射波振幅、预先确定的电磁波在第一层介质中的传播速度、各反射界面的反射波振幅和反射波的双程旅行时作为输入数据，迭代计算各层的厚度和反射波传播速度；
相对介电常数计算模块，用于基于各层的速度计算得到各层介质的相对介电常数。
(三)有益效果
本申请的有益效果是：本申请提出了一种探地雷达数据的反演方法包括：采用剖面法对待测区域进行探测扫描，得到实际测量数据；对实际测量数据进行数据预处理，得到待反演数据；从待反演数据中提取各个反射界面的反射波振幅和反射波的双程旅行时；将偏移距离、入射波振幅、预先确定的电磁波在第一层介质中的传播速度、各反射界面的反射波振幅和双程旅行时作为输入数据，迭代计算各层的厚度和反射波传播速度；基于各层的速度得到各层介质的相对介电常数。通过本申请的方法采集探地雷达数据简单、快速，可用于网格或直线上的大范围测量；并可实现准确、高效地探地雷达数据反演，极大地推动其在工程领域中的应用。
附图说明
本申请借助于以下附图进行描述：
图1为本申请一个实施例中的探地雷达数据的反演方法流程示意图；
图2为本申请另一个实施例中的某铁路原始共偏移雷达剖面示例图；
图3为本申请另一个实施例中的探地雷达数据处理流程示意图；
图4为本申请另一个实施例中的某铁路的零时刻校正后的共偏移雷达剖面示例图；
图5为本申请另一个实施例中的仪器增益系数及折线示例图；
图6为本申请另一个实施例中的振幅恢复后探地雷达剖面示例图；
图7为本申请另一个实施例中的第三次反射传播路径入射和反射路线图；
图8为本申请另一个实施例中的某铁路路基反演结果示例图；
图9为本申请另一个实施例中的平滑滑动平均后的反演结果示例图；
图10为本申请又一个实施例中的探地雷达数据的反演装置架构示意图。
具体实施方式
为了更好的解释本发明，以便于理解，下面结合附图，通过具体实施方式，对本发明作详细描述。
探地雷达(Ground Penetrating Radar，GPR)以探测地下不同介质的电磁性质差异为前提，利用高频短脉冲电磁波来确定地下介质分布规律。由于它具有操作灵活、非破坏性、高分辨率、费用低等优势，近年来被广泛应用于道路病害检测、地下管线探测、隧道衬砌质量检测、地质灾害、水文地质、考古勘探等领域。
目前探地雷达反演方法种类繁多，但通常会通过增加反演算法的时间复杂度以来降低反演结果与真实数据的误差，在反演的过程中增加了计算的迭代次数，导致运算时间大幅增加。为了解决现有技术中存在的问题，本申请提出了一种准确且高效的反演方法。
需要说明的是，本申请中的反演方法的假设条件，即：假设传播信号是平面电磁波，每道附近的地下介质为水平层状，每一层是均匀、非磁性、非导电、各向同性的。从动力学的角度考虑了反射和透射系数，消除了几何扩散的影响；忽略了天线耦合、固有衰减和散射效应。常规的探地雷达测量采用横向电(TE)配置。因此，本发明只考虑TE模式
实施例一
图1为本申请一个实施例中的探地雷达数据的反演方法流程示意图，如图1所示，该方法包括：
S10、获取探地雷达采用剖面法对待测区域进行探测扫描得到的实际测量数据；
S20、对实际测量数据进行数据预处理，得到待反演数据；
S30、从待反演数据中提取各个反射界面的反射波振幅和反射波的双程旅行时；
S40、将探地雷达发射天线和接收天线间的偏移距离、探地雷达的入射波振幅、预先确定的电磁波在第一层介质中的传播速度、各反射界面的反射波振幅和反射波的双程旅行时作为输入数据，迭代计算各层的厚度和反射波传播速度；
S50、基于各层的速度计算得到各层介质的相对介电常数。
通过本申请的方法采集探地雷达数据简单、快速，可用于网格或直线上的大范围测量；并可实现准确、高效地探地雷达数据反演，可以极大的推动其在工程领域中的应用。以下对本实施例中各个步骤展开进行具体说明。
本实施例中，待测区域为近水平层状地下结构，层数小于等于4，探地雷达对待测区域进行探测扫描时的探测扫描波长的确定方法，包括：
S01、获取预设扫描波长，将预设扫描波长作为初始扫描波长；
S02、对待测区域以初始扫描波长进行雷达探测，通过步骤S10-S40确定地层各层的厚度；
S03、对初始扫描波长进行递减，将得到的波长作为预设扫描波长；
S04、迭代执行步骤S01-03，直至满足预设的终止条件，并将最后一次迭代时的初始扫描波长作为探测扫描波长；预设的终止条件包括：
Figure BDA0003056910700000091
其中，λ表示探测扫描的波长，lmin表示地层的厚度最小值。
基于确定的探测扫描波长，探地雷达采用剖面法对待测区域进行探测扫描，得到待测区域的实际测量数据。剖面法测量时，发射天线和接收天线之间间距固定、二者同时沿侧线移动，从而得到雷达剖面。剖面图的横轴表示天线的位置，纵轴是反射波的双程走时。由于剖面法只需要发射和接收通道，系统设计非常简单，因此目前大多数商用雷达都使用这种观测方法。
需要说明的是，本实施例的反演方法适用于地层为4层以内的近水平层状地下结构，如果层数较多时，对于较深的地层来说，反演计算出的层速度比实际值偏大，层厚也比实际值偏大，且计算值与理论值之间的误差随着深度的增加而不断增加，反演结果虽能看出地下结构电性参数的变化大致趋势与层深，但反演结果与真实地下介质的电性结构会有差异，如果仅仅只是想了解近地表的几层结构，那么无论地下层数多或少，都可以使用该反演方法。
本实施例中，步骤S20包括以下步骤：
S21、对探地雷达测量数据进行直流滤波以去除数据中的直流分量，得到测量第一数据；
S22、选取雷达直达波信号的负峰值点时刻作为零时刻点，对测量第一数据进行零时刻校正，得到测量第二数据；
S23、对测量第二数据进行反增益处理以恢复探测雷达反射波的振幅，得到测量第三数据；
S24、对测量第三数据进行球面扩散振幅补偿得到测量第四数据；
对测量第三数据进行球面扩散振幅补偿得到测量第四数据的计算公式为：
A0＝Avteβt
其中，A0表示测量第四数据，为探地雷达发射天线发射出的没有经过扩散与介质吸收影响的真实振幅值；A表示测量第三数据；v表示电磁波在介质中传播的平均速度，t表示接收的反射振幅的双程旅行时，β表示介质的吸收系数。
S25、对测量第四数据进行带通滤波，过滤掉信号中的高频噪声信号和低频振荡信号，得到待反演数据；
对测量第四数据进行带通滤波时，采用正弦函数的平方作为镶边函数的滤波器，滤波器为：
Figure BDA0003056910700000101
其中，H(f)表示待反演数据，f1表示低通截止频率，f2表示低通频率，f3表示高通频率，f4表示高通截止频率。
需要说明的是，当选取直达波信号的负峰值点时刻作为零时刻点时，需对反演得到地层深度进行修正，修正公式为：
d'＝d+v×(T3-T2)
其中，d'表示修正后的地层深度，d表示反演得到的地层深度，v表示电磁波在每层介质中传播的平均速度，T3表示直达波信号的负峰值点时刻，T2表示直达波信号起跳点与负峰值点对应时刻的中点时刻。
以负相位峰值为时间起跳点进行反演，再以T3-T2为时间校正量对反演结果进行修正，可以得到更精准的地层厚度、深度以及速度信息。
本实施例中，迭代计算各层的厚度和反射波传播速度，包括：
S41、基于偏移距离、电磁波在第一层介质中的传播速度和第一层介质中的双程旅行旅行时，计算得到第一层介质的厚度与电磁波在第一层的入射角；
Figure BDA0003056910700000111
Figure BDA0003056910700000112
其中，x表示偏移距离，v1表示电磁波在第一层介质中的传播速度，TWT1表示第一介质中的双程旅行时，h1表示第一层介质的厚度，θ1表示电磁波在第一层介质中的入射角；
S42、基于入射波振幅和第一层反射波振幅，计算出第一个反射界面的反射系数，计算公式为：
Figure BDA0003056910700000113
其中，Ai1表示入射波振幅，As1表示第一个反射界面的反射波振幅，R1表示第一个反射界面的反射系数；
S43、令n＝2；
S44、基于第n-1个界面的反射系数、第n-1层的反射波传播速度、第n-1个界面的电磁波的入射角，由Snell方程计算得到第n层的反射波传播速度；
S45、基于前n-1层的厚度、前n层的反射波传播速度和第n层界面反射波的双程旅行时，根据雷达波反射路线与介质层厚度、偏移距离的几何关系计算得到第n层的厚度；
S46、基于偏移距离、前n层的反射波传播速度、前n层的厚度，计算得到得到沿第n次反射波路径各界面的电磁波的入射角，计算公式为：
Figure BDA0003056910700000121
其中，k表示当前界面的层数，k＝1,2,…,n；i表示界面层数。
S47、基于沿第n次反射波路径各界面的电磁波的入射角，计算前n-1层界面的透射系数，计算公式为：
Figure BDA0003056910700000122
其中，θk+1、θk分别表示电磁波在第k层和第k+1层介质中的入射角；
S48、基于提取的各反射界面的反射波振幅、前n-1层界面的透射系数，计算得到第n个界面的反射系数，计算公式为：
Figure BDA0003056910700000123
其中，Asn表示从所述待反演数据中提取的第n个反射界面的反射波振幅，i表示界面的层数，n表示界面的总层数，Ti表示第i层界面的透射系数；
S49、令n＝n+1，迭代执行步骤S44-S48，直至n＝4。
本实施例中，步骤S50包括：
S51、基于各层的速度计算得到各层介质的相对介电常数，计算公式为：
Figure BDA0003056910700000124
其中，c表示真空中电磁波的传播速度，v表示介质的电磁波传播速度，εr表示相对介电常数；
S52、采用平滑滑动平均对计算得到的相对介电常数进行处理，过滤掉异常值。
以下对相对介电常数的计算公式进行说明。
探地雷达的传播参数包括速度、振幅和衰减程度。介质对电磁波的影响不仅仅体现在速度的快慢，还体现在速度与频率的关系上。介质对电磁波衰减系数的影响也会影响探地雷达的探测深度。
介质的电磁波传播速度为：
Figure BDA0003056910700000131
其中，σ、ε、μ分别表示介质的电导率、介电常数和磁导率，α、ω分别表示向心加速度、角速度。
大多数岩石介质是非磁性和非导电的，当频率很高时，它们通常满足
Figure BDA0003056910700000132
因此可得：
Figure BDA0003056910700000133
本实施例利用探地雷达数据反演方法估计层状模型的各层层厚与速度，采用递推方法计算了各层反射界面处的反射系数，进而求出各层介质的厚度和速度，进一步得到各层介质的相对介电常数以完成反演，反演结果准确，反演过程高效、简便。
实施例二
本申请第二方面通过另一实施例提供了一种探地雷达数据的反演方法，本实施例中的探地雷达数据是采用SIR-20型雷达对中国某段铁路进行共偏移距测量，天线中心频率为400MHz，时窗范围60ns，采样点数512，而测量得到数据。图2为本申请另一个实施例中的某铁路的原始共偏移雷达剖面图，如图2所示，其中横轴为雷达水平方向的移动距离，纵轴为电磁波旅行时间。
图3为本申请另一个实施例中的探地雷达数据处理流程示意图，以下结合图3对数据处理的各个步骤进行展开说明。
S1、读入探地雷达数据。
读入的探地雷达数据可以是.dzt或.mat格式的数据。需要说明的是，这里仅仅是对数据格式进行示例性地说明，并不对数据格式进行限定。
S2、数据预处理。
实测数据中会出现杂波干扰，为了使各个界面的反射振幅能够正确反映出反射界面两端介质的电磁性质，在对实测数据反演之前，需要做一系列的数据处理操作。
图4为本申请另一个实施例中的某铁路的零时刻校正后的共偏移雷达剖面图，其中横轴为雷达水平方向的移动距离，纵轴为电磁波旅行时间。首先，对实测数据进行DC滤波以去除直流，再进行零时刻校正得到正确的双程旅行时间，处理后的剖面图如图4所示。
在振幅拾取前，必须对探地雷达数据采用合适的处理流程。由于本文所用的实测数据在采集过程中被仪器增益过，图4可看出16ns附近反射层的反射振幅由于增益作用甚至大于直达波的振幅，这明显是不合常理的。所以要对实测雷达数据进行反增益处理以恢复原始的反射波振幅。
图5为本申请另一个实施例中的仪器增益系数及折线图，图5中(a)为仪器增益系数，图5中(b)为根据增益系数绘制成的折线图，横轴为时间，纵轴为增益系数。利用数据处理软件读取.dzt格式实测数据的头文件，得到其采集数据时设置的增益系数，增益方式为五点线性增益，如图5中(a)所示。将增益系数绘制成折线图，如图5中(b)所示，图中横坐标为时间，纵坐标为增益系数。
普通的增益处理会破坏反射界面的振幅与直达波振幅的相对关系，可能会导致反演结果不准确，不建议使用增益处理。因此考虑到使用振幅恢复的方式来得到10ns以下的反射层信息。
由于反演方法基于各个界面的反射振幅，而探地雷达接收到的反射波振幅由于波前扩散和介质对电磁波的吸收作用，在时间轴上会逐渐衰减，回波幅度逐渐减弱，导致难以分辨深部的反射界面。为了使深部反射界面的反射振幅更清晰、并且让反射波振幅与反射界面上下介质的电磁性质有关，需要进行振幅恢复处理。
假设地下介质均匀，则距离发射天线r处的电磁波的振幅A为:
Figure BDA0003056910700000141
其中，A0是探地雷达发射天线发射出的没有经过扩散与介质吸收影响的真实振幅值；
Figure BDA0003056910700000151
是雷达波的波前扩散因子；β是介质的吸收系数；t为接收的反射振幅的双程走时。逆推上式得到反射波的原始振幅值A0:
A0＝Areβt
反射波实际的旅行路径r＝vt，其中v是电磁波在介质中传播的平均速度,将其带入上式可以得到:
A0＝Avteβt
反射波的真实振幅可由上式恢复。对于具有可忽略的内在衰减的介质(β≈0)，应用恒定速度的扩散校正是一种可接受的振幅衰减补偿。
图6为本申请另一个实施例中的振幅恢复后探地雷达剖面图，其中横轴为雷达水平方向的移动距离，纵轴为电磁波旅行时间。如图6所示，深部的振幅明显增强，同时各种干扰也被加强了，图像比较杂乱，因此还需进行滤波处理。
带通滤波可以过滤掉信号中的高频噪声和低频振荡。由于大部分噪声的频率是在雷达天线信号频带之外的，因此带通滤波可以有效去除大多数相干和随机噪声成分。带通滤波器采用实施例1中的滤波器，带通滤波的下限截止频率通常为天线主频率的一半，上限截止频率通常为天线主频率的2倍。因此对于400MHz天线来说，截止频率应为200～800MHz。
S3、拾取反射界面的振幅及走时。
图7为本申请另一个实施例中的第三次反射传播路径入射和反射路线图；S代表发射天线，T代表接收天线，x是天线之间的偏移距，hi是第i层介质层的厚度,vi是电磁波在第i层介质的速度，θi是电磁波在第i层介质的入射角度，Δxi传播路径的水平投影。以下结合图7对本申请中的反演方法原理进行说明。在通过实施例一中步骤S41计算得到第一层介质的厚度h1、第一层介质中的入射角θ1、第一个反射界面的反射系数R1后，利用Snell方程求得第二层介质中的雷达波速度v2：
Figure BDA0003056910700000152
其中，v1表示电磁波在第一层介质中的传播速度，θ2表示电磁波在第二层介质中的入射角，可以通过重新排列的TE模式的菲涅尔方程得到：
Figure BDA0003056910700000161
由于在第n个循环,第n层界面的反射波旅行时间TWTn可表示为：
Figure BDA0003056910700000162
Figure BDA0003056910700000163
其中，hi和vi分别是第i层介质的层厚和电磁波速度。
如果已知前n-1层的厚度、前n层的速度和第n层界面反射波的双程旅行时，则第n层的厚度hn为下列三次方程的唯一正解：
Figure BDA0003056910700000164
其中：
Figure BDA0003056910700000165
Figure BDA0003056910700000166
Figure BDA0003056910700000167
Figure BDA0003056910700000168
第k个界面的入射角和第k-1个界面的透射角在水平层状介质中是相等的。该角θk与第k层传播路径的水平投影Δxk有关：
Figure BDA0003056910700000169
由于θk比较小，再结合Snell方程，可以得到：
Figure BDA00030569107000001610
根据几何条件：
Figure BDA00030569107000001611
得到沿第n次反射波路径各界面的入射角θk，即：
Figure BDA0003056910700000171
然后利用TE模式的菲涅尔方程，计算前n-1层界面的反射系数Rk和透射系数Tk：
Figure BDA0003056910700000172
Tk＝1+Rk(k＝1,2,…,n-1)
考虑到传播路径的对称性，第n界面的入射振幅Ain和反射振幅Arn分别为：
Figure BDA0003056910700000173
Figure BDA0003056910700000174
那么，第n个界面的反射系数Rn为：
Figure BDA0003056910700000175
第n+1层的速度vn+1由Snell方程给出：
Figure BDA0003056910700000176
Figure BDA0003056910700000177
将该方法应用于所有需要解释的探地雷达数据道，就可以得到地下各层的层厚和电磁波速度。
S4、执行反演得到反演结果图。
图8为本申请另一个实施例中的某铁路路基反演结果示例图，图9为本申请另一个实施例中的平滑滑动平均后的反演结果示例图，其中横轴为雷达水平方向的移动距离，纵轴为深度。图8中可以清晰的看出两层反射界面的位置，但由于反射振幅的不连续，导致某些道的反射振幅选取不恰当，进而引起反演图中介质介电常数横向变化较大，为了突出剖面的主要速度变化并过滤掉异常值，采用平滑滑动平均对反演结果进行处理，处理结果如图9所示。
实施例三
本申请第二方面通过实施例三提供了一种探地雷达数据的反演装置，该装置可执行上述的探地雷达数据的反演方法。图10为本申请又一个实施例中的探地雷达数据的反演装置架构示意图，如图10所示，该装置包括：
测量数据获取模块101，用于获取探地雷达采用剖面法对待测区域进行探测扫描得到的实际测量数据；
数据预处理模块102，用于对实际测量数据进行数据预处理，得到待反演数据；
反射波振幅提取模块103，用于从待反演数据中提取各个反射界面的反射波振幅和反射波的双程旅行时；
迭代计算模块104，用于将探地雷达发射天线和接收天线间的偏移距离、探地雷达的入射波振幅、预先确定的电磁波在第一层介质中的传播速度、各反射界面的反射波振幅和反射波的双程旅行时作为输入数据，迭代计算各层的厚度和反射波传播速度；
相对介电常数计算模块105，用于基于各层的速度计算得到各层介质的相对介电常数。
上述探地雷达数据的反演装置可执行本申请实施例所提供的探地雷达数据的反演方法，具备执行方法相应的功能模块和有益效果。至于其中各个功能模块所执行的处理方法，例如，测量数据获取模块101、数据预处理模块102、反射波振幅提取模块103、迭代计算模块104、相对介电常数计算模块105可参照上述方法实施例中的描述，此处不再进行赘述。
显然，本领域的技术人员可以对本发明进行各种修改和变型而不脱离本发明的精神和范围。这样，倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之内，则本发明也应该包含这些修改和变型在内。

Claims (10)

1.一种探地雷达数据的反演方法，其特征在于，该方法包括：
S10、获取探地雷达采用剖面法对待测区域进行探测扫描得到的实际测量数据；
S20、对所述实际测量数据进行数据预处理，得到待反演数据；
S30、从所述待反演数据中提取各个反射界面的反射波振幅和反射波的双程旅行时；
S40、将所述探地雷达发射天线和接收天线间的偏移距离、所述探地雷达的入射波振幅、预先确定的电磁波在第一层介质中的传播速度、各反射界面的反射波振幅和反射波的双程旅行时作为输入数据，迭代计算各层的厚度和反射波传播速度；
S50、基于各层速度计算得到各层介质的相对介电常数。
2.根据权利要求1所述的方法，其特征在于，所述待测区域为近水平层状地下结构，层数小于等于4，所述探地雷达对所述待测区域进行探测扫描时的探测扫描波长的确定方法，包括：
S01、获取预设扫描波长，将所述预设扫描波长作为初始扫描波长；
S02、对所述待测区域以所述初始扫描波长进行雷达探测，通过步骤S10-S40确定地层各层的厚度；
S03、对所述初始扫描波长进行递减，将得到的波长作为预设扫描波长；
S04、迭代执行步骤S01-03，直至满足预设的终止条件，并将最后一次迭代时的初始扫描波长作为探测扫描波长；预设的终止条件包括：
Figure FDA0003056910690000011
其中，λ表示探测扫描的波长，lmin表示地层的厚度最小值。
3.根据权利要求1所述的方法，其特征在于，步骤S20包括：
对所述探地雷达测量数据进行直流滤波以去除数据中的直流分量，得到测量第一数据；
对所述测量第一数据进行零时刻校正，得到测量第二数据；
对所述测量第二数据进行反增益处理以恢复所述探测雷达反射波的振幅，得到测量第三数据；
对所述测量第三数据进行球面扩散振幅补偿得到测量第四数据；
对所述测量第四数据进行带通滤波，过滤掉信号中的高频噪声信号和低频振荡信号，得到待反演数据。
4.根据权利要求3所述的方法，其特征在于，选取雷达直达波信号的负峰值点时刻作为零时刻点，对所述测量第一数据进行零时刻校正，得到测量第二数据。
5.根据权利要求4所述的方法，其特征在于，当选取直达波信号的负峰值点时刻作为零时刻点时，对反演得到地层深度进行修正，修正公式为：
d'＝d+v×(T3-T2)
其中，d'表示修正后的地层深度，d表示反演得到的地层深度，v表示电磁波在每层介质中传播的平均速度，T3表示直达波信号的负峰值点时刻，T2表示直达波信号起跳点与负峰值点对应时刻的中点时刻。
6.根据权利要求3所述的方法，其特征在于，对所述测量第三数据进行球面扩散振幅补偿得到测量第四数据的计算公式为：
A0＝Avteβt
其中，A0表示测量第四数据，为探地雷达发射天线发射出的没有经过扩散与介质吸收影响的真实振幅值；A表示测量第三数据；v表示电磁波在每层介质中传播的平均速度，t表示接收的反射振幅的双程旅行时，β表示每层介质的吸收系数。
7.根据权利要求3所述的方法，其特征在于，对所述测量第四数据进行带通滤波时，采用正弦函数的平方作为镶边函数的滤波器，所述滤波器为：
Figure FDA0003056910690000031
其中，H(f)表示待反演数据，f1表示低通截止频率，f2表示低通频率，f3表示高通频率，f4表示高通截止频率。
8.根据权利要求3所述的方法，其特征在于，迭代计算各层的厚度和反射波传播速度，包括：
S41、基于所述偏移距离、电磁波在第一层介质中的传播速度和第一层介质中的双程旅行旅行时，计算得到第一层介质的厚度与电磁波在第一层的入射角；
Figure FDA0003056910690000032
Figure FDA0003056910690000033
其中，x表示偏移距离，v1表示电磁波在第一层介质中的传播速度，TWT1表示第一介质中的双程旅行时，h1表示第一层介质的厚度，θ1表示电磁波在第一层介质中的入射角；
S42、基于所述入射波振幅和第一层反射波振幅，计算出第一个反射界面的反射系数，计算公式为：
Figure FDA0003056910690000034
其中，Ai1表示入射波振幅，As1表示第一个反射界面的反射波振幅，R1表示第一个反射界面的反射系数；
S43、令n＝2；
S44、基于第n-1个界面的反射系数、第n-1层的反射波传播速度、第n-1个界面的电磁波的入射角，由Snell方程计算得到第n层的反射波传播速度；
S45、基于前n-1层的厚度、前n层的反射波传播速度和第n层界面反射波的双程旅行时，根据雷达波反射路线与介质层厚度、所述偏移距离的几何关系计算得到第n层的厚度；
S46、基于所述偏移距离、前n层的反射波传播速度、前n层的厚度，计算得到得到沿第n次反射波路径各界面的电磁波的入射角，计算公式为：
Figure FDA0003056910690000041
其中，k表示当前界面的层数，k＝1,2,…,n；i表示界面层数；
S47、基于沿第n次反射波路径各界面的电磁波的入射角，计算前n-1层界面的透射系数，计算公式为：
Figure FDA0003056910690000042
其中，θk+1、θk分别表示电磁波在第k层和第k+1层介质中的入射角；
S48、基于提取的各反射界面的反射波振幅、前n-1层界面的透射系数，计算得到第n个界面的反射系数，计算公式为：
Figure FDA0003056910690000043
其中，Asn表示从所述待反演数据中提取的第n个反射界面的反射波振幅，i表示界面的层数，n表示界面的总层数，Ti表示第i层界面的透射系数；
S49、令n＝n+1，迭代执行步骤S44-S48，直至n＝4。
9.根据权利要求3所述的方法，其特征在于，步骤S50包括：
S51、基于各层的速度计算得到各层介质的相对介电常数，计算公式为：
Figure FDA0003056910690000044
其中，c表示真空中电磁波的传播速度，v表示介质的电磁波传播速度，εr表示相对介电常数；
S52、采用平滑滑动平均对计算得到的相对介电常数进行处理，过滤掉异常值。
10.一种探地雷达数据的反演装置，其特征在于，该装置包括：
测量数据获取模块，用于获取探地雷达采用剖面法对待测区域进行探测扫描得到的实际测量数据；
数据预处理模块，用于对所述实际测量数据进行数据预处理，得到待反演数据；
反射波振幅提取模块，用于从所述待反演数据中提取各个反射界面的反射波振幅和反射波的双程旅行时；
迭代计算模块，用于将所述探地雷达发射天线和接收天线间的偏移距离、所述探地雷达的入射波振幅、预先确定的电磁波在第一层介质中的传播速度、各反射界面的反射波振幅和反射波的双程旅行时作为输入数据，迭代计算各层的厚度和反射波传播速度；
相对介电常数计算模块，用于基于各层的速度计算得到各层介质的相对介电常数。
CN202110502385.1A 2021-05-08 2021-05-08 探地雷达数据的反演方法和装置 Pending CN113156427A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110502385.1A CN113156427A (zh) 2021-05-08 2021-05-08 探地雷达数据的反演方法和装置

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110502385.1A CN113156427A (zh) 2021-05-08 2021-05-08 探地雷达数据的反演方法和装置

Publications (1)

Publication Number Publication Date
CN113156427A true CN113156427A (zh) 2021-07-23

Family

ID=76874118

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110502385.1A Pending CN113156427A (zh) 2021-05-08 2021-05-08 探地雷达数据的反演方法和装置

Country Status (1)

Country Link
CN (1) CN113156427A (zh)

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114111553A (zh) * 2021-11-25 2022-03-01 安徽理工大学 一种快速获取新增耕地重构土体填充层厚度的方法
CN115421132A (zh) * 2022-11-04 2022-12-02 北京锐达仪表有限公司 一种具有安装角度误差自修正功能的3d扫描雷达
CN116840807A (zh) * 2023-09-01 2023-10-03 中国科学院空天信息创新研究院 一种基于探地雷达系统的全波反演介电常数估计方法
CN117890904A (zh) * 2024-03-14 2024-04-16 中南大学 一种探地雷达衍射属性场的提取方法、存储介质及设备

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106405061A (zh) * 2016-09-22 2017-02-15 北京林业大学 一种基于雷达波的木质体内部异常无损探测系统
CN109633762A (zh) * 2019-01-07 2019-04-16 吉林大学 基于正弦函数的相关性约束条件联合反演重磁数据的方法

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106405061A (zh) * 2016-09-22 2017-02-15 北京林业大学 一种基于雷达波的木质体内部异常无损探测系统
CN109633762A (zh) * 2019-01-07 2019-04-16 吉林大学 基于正弦函数的相关性约束条件联合反演重磁数据的方法

Non-Patent Citations (5)

* Cited by examiner, † Cited by third party
Title
MEIQI XUE: "A Co-offset GPR Data Inversion Based on Ray Theory", 《EARTH AND ENVIRONMENTAL SCIENCE660(2021) 012047》, pages 1 - 6 *
李成方: "振幅处理技术在探地雷达资料处理中的应用", 《物探化探计算技术》, pages 223 - 227 *
王超: "两种物探方法在岩溶塌陷中的综合应用效果研究", 《中国优秀硕士学位论文全文数据库基础科学辑》, pages 37 - 38 *
蔡佳琪: "利用探测雷达探测铁路路基含水率", 《中国优秀硕士学位论文全文数据库 工程科技Ⅱ辑》, pages 37 *
袁明德: "浅析探地雷达的分辨率", 《物探与化探》, pages 28 - 32 *

Cited By (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114111553A (zh) * 2021-11-25 2022-03-01 安徽理工大学 一种快速获取新增耕地重构土体填充层厚度的方法
CN114111553B (zh) * 2021-11-25 2023-12-08 安徽理工大学 一种快速获取新增耕地重构土体填充层厚度的方法
CN115421132A (zh) * 2022-11-04 2022-12-02 北京锐达仪表有限公司 一种具有安装角度误差自修正功能的3d扫描雷达
CN115421132B (zh) * 2022-11-04 2023-02-21 北京锐达仪表有限公司 一种具有安装角度误差自修正功能的3d扫描雷达
CN116840807A (zh) * 2023-09-01 2023-10-03 中国科学院空天信息创新研究院 一种基于探地雷达系统的全波反演介电常数估计方法
CN116840807B (zh) * 2023-09-01 2023-11-10 中国科学院空天信息创新研究院 一种基于探地雷达系统的全波反演介电常数估计方法
CN117890904A (zh) * 2024-03-14 2024-04-16 中南大学 一种探地雷达衍射属性场的提取方法、存储介质及设备

Similar Documents

Publication Publication Date Title
CN113156427A (zh) 探地雷达数据的反演方法和装置
CA2712439C (en) Characterizing spatial variablility of surface waves in seismic processing
CN102590862B (zh) 补偿吸收衰减的叠前时间偏移方法
WO2009124446A1 (zh) 三维小面元电磁连续阵列数据采集方法
CN107479042B (zh) 一种表层岩溶带空间蓄水能力的估算方法
CN109884709B (zh) 一种基于面波旅行时层析的转换波静校正方法
CN111103621A (zh) 一种主动源共成像点叠加多道面波分析方法
CN102073064B (zh) 一种利用相位信息提高速度谱分辨率的方法
CN109541690B (zh) 一种浅层介质结构面松散程度评价方法
CN106526678B (zh) 一种反射声波测井的波场分离方法及装置
CN103424777A (zh) 一种提高地震成像分辨率的方法
CN114235962B (zh) 一种面向各向异性结构的超声导波成像方法及系统
CN106291542A (zh) 一种隧道三维成像方法
CN114460649A (zh) 一种深海近底拖曳式多道地震接收阵列形态重建方法
CN112666554A (zh) 一种沥青路面雷达振幅特征裂缝宽度识别方法
CN111142165A (zh) 一种利用探地雷达获取含水层的水位信息的方法
CN104714251A (zh) 用于同相轴自动拾取的倾斜叠加峰值振幅处边缘检测法
CN104793237A (zh) 一种获得宽频可控震源扫描信号的方法和装置
CN113391283A (zh) 基于探地雷达的土壤分层信息识别方法和装置
CN108508093A (zh) 一种工件缺陷高度的检测方法及系统
CN114415234B (zh) 基于主动源面波频散和h/v确定浅地表横波速度的方法
CN109541689B (zh) 一种基于反射波能量特征的介质密实度评价方法
CN102798888A (zh) 一种利用非零井源距数据计算纵横波速度比的方法
CN111352153A (zh) 一种基于瞬时相位互相关加权的微地震干涉定位方法
CN106443674A (zh) 一种基于衍射和成像与最小熵技术的探地雷达波速估计方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination