CN111123374A - 一种基于匹配滤波的探地雷达全波形反演方法 - Google Patents

一种基于匹配滤波的探地雷达全波形反演方法 Download PDF

Info

Publication number
CN111123374A
CN111123374A CN201911359612.9A CN201911359612A CN111123374A CN 111123374 A CN111123374 A CN 111123374A CN 201911359612 A CN201911359612 A CN 201911359612A CN 111123374 A CN111123374 A CN 111123374A
Authority
CN
China
Prior art keywords
data
full
ground penetrating
model
waveform inversion
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.)
Withdrawn
Application number
CN201911359612.9A
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.)
Changzhou Institute of Technology
Original Assignee
Changzhou Institute of Technology
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 Changzhou Institute of Technology filed Critical Changzhou Institute of Technology
Priority to CN201911359612.9A priority Critical patent/CN111123374A/zh
Publication of CN111123374A publication Critical patent/CN111123374A/zh
Withdrawn 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
    • G01S7/00—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
    • G01S7/02—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00
    • 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
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V3/00—Electric or magnetic prospecting or detecting; Measuring magnetic field characteristics of the earth, e.g. declination, deviation
    • G01V3/12—Electric or magnetic prospecting or detecting; Measuring magnetic field characteristics of the earth, e.g. declination, deviation operating with electromagnetic waves
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V3/00—Electric or magnetic prospecting or detecting; Measuring magnetic field characteristics of the earth, e.g. declination, deviation
    • G01V3/38—Processing data, e.g. for analysis, for interpretation, for correction

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Electromagnetism (AREA)
  • Geology (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • Geophysics (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Environmental & Geological Engineering (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

本发明涉及一种对探地雷达全波形反演技术的改进，具体为一种基于匹配滤波的探地雷达全波形反演方法,包括如下步骤：步骤一：建立探地雷达全波形反演的真实模型；步骤二：基于真实模型进行正演运算并抽取道集，获得观测雷达数据，该数据相当于实际工程中采集到的雷达数据；步骤三：建立探地雷达全波形反演的初始模型并设定反演终止精度；步骤四：基于初始模型进行正演运算并抽取道集，获得理论雷达数据；通过合理设置滤波算子并对理论雷达数据进行改造，降低了理论数据与实际采集数据之间的差异，从而降低了周期跳跃现象发生的可能性，使反演不易陷入局部极小点，进而保障探地雷达全波形反演的稳定性和精度。

Description

一种基于匹配滤波的探地雷达全波形反演方法
技术领域
本发明涉及一种对探地雷达全波形反演技术的改进，具体为一种基于匹配滤波的探地雷达全波形反演方法。
背景技术
探地雷达是一种高效且具有较高分辨率的无损检测手段，可应用于隧道质量检测、路面密实度检测、混凝土构件非破坏性检测、地基夯实与加固检测、地灾评估和土壤地质信息评定等许多方面。探地雷达利用反射系数确定被检测体的结构及缺陷，但由于电磁波的传播过程极其复杂，人类活动对成像结果也有很大干扰，因此获得的地质雷达剖面往往不够准确，解释起来也与被检测体的真实结构存在差异，因此对雷达数据进行优化处理以期进行合理解释迫在眉睫。
得益于地震全波形反演的发展以及麦克斯韦方程组在解释形式上与声波方程的相似性，雷达数据可借鉴地震数据的处理方法，利用全波形反演进行数据的优化。Ernst(2006)实现了时间域探地雷达全波形反演，并对亚波长异常体进行了精准定位；Meles(2010)利用钻孔雷达全波形反演获得了河边较清晰的地下地层结构；胡周文(2018) 利用全波形反演对混凝土进行了无损检测，冯德山(2018)等利用 GPU并行实现了雷达数据的双参数反演，可见雷达数据的全波形反演是可行的。
尽管全波形反演可以较好的提高探地雷达的分辨率，但反演本身也有一定的局限性。全波形反演使用局部优化算法进行模型优化，其弊端是当反演的初始模型精度不足时，全波形反演会由于初始波形和理论波形间较大的差异而陷入局部极小点，从而导致反演失败，因此利用全波形反演进行雷达数据成像时必须要考虑的一个重要问题就是如何阻止反演陷入局部极小，进而保障反演的稳定性。
探地雷达全波形反演技术经过二十几年的发展，目前已取得了长足的进步，但由于算法本身的缺陷，该方法在应用时仍存在许多问题，其中的一个问题便是周期跳跃。周期跳跃是由于理论波形和实际波形之间差异较大造成的，而周期跳跃一旦发生便会导致反演陷入局部极小点，进而导致探地雷达全波形反演无法收敛。因此，要想充分发挥全波形反演高精度的优势，使雷达数据反演剖面上各点正确归位，就必须克服周期跳跃现象，确保反演结果正确收敛于全局极小点。
发明内容
为了解决这一问题，本发明提出一种能够提高雷达数据反演剖面的分辨率和雷达数据的数据解释可靠性的基于匹配滤波的探地雷达全波形反演方法。
为解决上述技术问题，本发明所采用的技术方案为：一种基于匹配滤波的探地雷达全波形反演方法，包括如下步骤：
步骤一：建立探地雷达全波形反演的真实模型；
步骤二：基于真实模型进行正演运算并抽取道集，获得观测雷达数据，该数据相当于实际工程中采集到的雷达数据；
步骤三：建立探地雷达全波形反演的初始模型并设定反演终止精度；
步骤四：基于初始模型进行正演运算并抽取道集，获得理论雷达数据；
步骤五：构建改造雷达波形的滤波算子；
步骤六：利用滤波算子对理论雷达数据进行改造；
步骤七：构建目标函数；
步骤八：分别求取介电常数和电导率的梯度；
步骤九：确定步长并更新模型；
步骤十：将步骤九中得到的模型作为新的初始模型，重复步骤四到九直到反演结束，获得最终的反演模型。
作为优选，所述步骤一中，真实模型为部分marmousi模型，并根据介电常数和电导率的数值变化范围进行赋值。
作为优选，所述步骤二中，正演运算采用基于交错网格的时间域有限差分法进行数值模拟，边界采用完全匹配层边界条件。
作为优选，所述步骤三中，初始模型是经平滑处理后获得。
作为优选，所述步骤五中，滤波算子的构建基于最小二乘法，将滤波算子作用于模拟数据Dcal，并使作用后的数据与基于真实模型的观测数据Dobs之间差异最小，即可得到理想的滤波算子f。
作为优选，所述步骤七中，目标函数用二范数形式表示，形如
Figure BDA0002336832670000031
其中Dcalf表示滤波后的模拟数据。
作为优选，所述步骤八中，分别对目标函数求取关于介电常数和电导率的一阶导数，其梯度通式可表示为
Figure BDA0002336832670000032
其中，m＝(m1，m2),m1和m2分别为介电常数和电导率，T表示转置，A为稀疏算子，Δd为滤波后的模拟数据和观测数据之差。
作为优选，所述步骤九中，所述步长的选择使用非单调线搜索方法。
本发明所达到的有益效果：本发明的基于匹配滤波的探地雷达全波形反演方法通过合理设置滤波算子并对理论雷达数据进行改造，降低了理论数据与实际采集数据之间的差异，从而降低了周期跳跃现象发生的可能性，使反演不易陷入局部极小点，进而保障探地雷达全波形反演的稳定性和精度。
附图说明
图1为本发明的方法流程示意图；
图2为本发明的介电常数的真实模型、初始模型、常规反演模型和本发明反演模型；
图3为本发明的电导率的真实模型、初始模型、常规反演模型和本发明反演模型。
具体实施方式
下面结合附图对本发明作进一步描述。以下实施例仅用于更加清楚地说明本发明的技术方案，而不能以此来限制本发明的保护范围。
图1所示为本方法的流程示意图。在反演之初首先建立用于模拟实际采集数据的真实模型，并在此基础上进行正演模拟得到观测数据；接下来建立初始模型，设置相应的参数以及反演的终止精度；随后基于初始模型进行正演模拟以获得理论数据；接下来构建滤波算子并利用滤波算子对理论雷达数据进行改造；下一步构建目标函数并在此基础上对目标函数求一阶导数，分别求取介电常数和电导率的梯度；再下一步，利用非单调线搜索获得步长，分别更新介电常数和电导率，得到更新后的参数模型。具体步骤如下：
步骤一：建立探地雷达全波形反演的真实模型。本发明以模型为例进行方法验证，真实模型为部分marmousi模型，并根据介电常数和电导率的数值变化范围进行赋值，介电常数真实模型和电导率真实模型分别如图2(a)和图3(a)所示；
步骤二：基于真实模型进行正演运算并抽取道集，获得观测雷达数据，该数据相当于实际工程中采集到的雷达数据。正演时采用基于交错网格的时间域有限差分法进行数值模拟，边界采用完全匹配层边界条件。类比地震全波形反演中正演部分的方程形式，探地雷达全波形反演的波动方程可写作：A(m1 m2)U＝S
其中A为正演算子，U为波场，S为场源,m1和m2分别为介电常数和电导率参数。当m1和m2均表示真实模型时，由波场U抽取道集得到的即是观测雷达数据；
步骤三：建立探地雷达全波形反演的初始模型并设定反演终止精度。初始模型为真实模型经平滑所得，介电常数初始模型和电导率初始模型分别如图2(b)和图3(b)所示；
步骤四：基于初始模型进行正演运算并抽取道集，获得理论雷达数据。步骤四与步骤二相似，正演方程仍为A(m1 m2)U＝S
其中A、U、S、m1和m2的意义同步骤四。当m1和m2均表示初始模型时，由波场U抽取道集得到的即是理论雷达数据；
步骤五：构建改造雷达波形的滤波算子。滤波算子的构建基于最小二乘法。通过将滤波算子作用于模拟数据Dcal，并使作用后的数据与基于真实模型的观测数据Dobs之间差异最小，即可得到理想的滤波算子f。通过最小化下式C＝‖Dobs-Dcalf‖2
即可得到最佳滤波算子f；
步骤六：利用滤波算子对理论雷达数据进行改造；
步骤七：构建目标函数。目标函数用二范数形式表示，目标函数表示为
Figure BDA0002336832670000061
其中Dcalf表示滤波后的模拟数据，也就是改造后的雷达数据；
步骤八：分别求取介电常数和电导率的梯度。其梯度通式可表示为：
Figure BDA0002336832670000062
其中，m＝(m1，m2),m1和m2分别为介电常数和电导率，T表示转置，A为稀疏算子，Δd为滤波后的模拟数据和观测数据之差。
本发明中两个参数顺序反演，即反演其中一个时，另一参数作为常量，梯度计算时也如此，计算介电常数梯度时，将电导率固定为常量，计算电导率梯度时，将介电常数看作常量；
步骤九：确定步长并更新模型。步长的大小利用非单调线搜索方法来确定。在反演过程中，反演的目标是使目标函数下降至误差允许范围，而并不要求目标函数严格单调下降。一味要求目标函数的单调下降会降低目标函数的收敛速度，增大计算量，因此本发明使用相对更高效的非单调线性搜索方法。
步长γi确定后，参数模型的更新形式可表示为：
Figure BDA0002336832670000071
其中，mi表示介电常数和电导率的第i次迭代模型，
Figure BDA0002336832670000072
表示模型为mi时计算出的介电常数和电导率的梯度；
步骤十：将步骤九中得到的模型作为新的初始模型，重复步骤四到九直到反演结束，由此可获得最终的反演模型。
为便于介绍基于匹配滤波的探地雷达全波形反演方法，本发明使用一实施案例加以说明。
图2和图3所示为验算该方法所用的真实模型、初始模型和其反演结果，其中图2a为介电常数的真实模型，图2b为介电常数的初始模型，图2c为用常规方法反演出的介电常数模型，图2d为用本发明中的方法反演出的介电常数模型。图3a为电导率的真实模型，图3b为电导率的初始模型，图3c为用常规方法反演出的电导率模型，图3d为用本发明中的方法反演出的电导率模型。模型实际深度为1.5m×2m，被离散为30×40的网格，网格间距0.05m。
通过对波形进行改造避免了周期跳跃现象的出现，从而提高了探地雷达全波形反演的稳定性和精度。根据上述介绍可知全波形反演受其算法影响易陷入局部极小点导致反演失败，针对这一问题，本发明提出了基于匹配滤波的探地雷达全波形反演方法。与传统方法相比，本发明中的方法通过合理设置滤波算子并对理论雷达数据进行改造，降低了理论数据与实际采集数据之间的差异，从而降低了周期跳跃现象发生的可能性，使反演不易陷入局部极小点，进而保障了探地雷达全波形反演的稳定性和精度。
以上仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明技术原理的前提下，还可以做出若干改进和变形，这些改进和变形也应视为本发明的保护范围。

Claims (8)

1.一种基于匹配滤波的探地雷达全波形反演方法，其特征在于：包括如下步骤：步骤一：建立探地雷达全波形反演的真实模型；
步骤二：基于真实模型进行正演运算并抽取道集，获得观测雷达数据，该数据相当于实际工程中采集到的雷达数据；
步骤三：建立探地雷达全波形反演的初始模型并设定反演终止精度；
步骤四：基于初始模型进行正演运算并抽取道集，获得理论雷达数据；
步骤五：构建改造雷达波形的滤波算子；
步骤六：利用滤波算子对理论雷达数据进行改造；
步骤七：构建目标函数；
步骤八：分别求取介电常数和电导率的梯度；
步骤九：确定步长并更新模型；
步骤十：将步骤九中得到的模型作为新的初始模型，重复步骤四到九直到反演结束，获得最终的反演模型。
2.根据权利要求1所述的基于匹配滤波的探地雷达全波形反演方法，其特征在于：所述步骤一中，真实模型为部分marmousi模型，并根据介电常数和电导率的数值变化范围进行赋值。
3.根据权利要求1所述的基于匹配滤波的探地雷达全波形反演方法，其特征在于：所述步骤二中，正演运算采用基于交错网格的时间域有限差分法进行数值模拟，边界采用完全匹配层边界条件。
4.根据权利要求1所述的基于匹配滤波的探地雷达全波形反演方法，其特征在于：所述步骤三中，初始模型是经平滑处理后获得。
5.根据权利要求1所述的基于匹配滤波的探地雷达全波形反演方法，其特征在于：所述步骤五中，滤波算子的构建基于最小二乘法，将滤波算子作用于模拟数据Dcal，并使作用后的数据与基于真实模型的观测数据Dobs之间差异最小，即可得到理想的滤波算子f。
6.根据权利要求1所述的基于匹配滤波的探地雷达全波形反演方法，其特征在于：所述步骤七中，目标函数用二范数形式表示，形如
Figure FDA0002336832660000021
其中Dcalf表示滤波后的模拟数据。
7.根据权利要求1所述的基于匹配滤波的探地雷达全波形反演方法，其特征在于：所述步骤八中，分别对目标函数求取关于介电常数和电导率的一阶导数，其梯度通式可表示为
Figure FDA0002336832660000022
其中，m＝(m1，m2),m1和m2分别为介电常数和电导率，T表示转置，A为稀疏算子，Δd为滤波后的模拟数据和观测数据之差。
8.根据权利要求1所述的基于匹配滤波的探地雷达全波形反演方法，其特征在于：所述步骤九中，所述步长的选择使用非单调线搜索方法。
CN201911359612.9A 2019-12-25 2019-12-25 一种基于匹配滤波的探地雷达全波形反演方法 Withdrawn CN111123374A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201911359612.9A CN111123374A (zh) 2019-12-25 2019-12-25 一种基于匹配滤波的探地雷达全波形反演方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201911359612.9A CN111123374A (zh) 2019-12-25 2019-12-25 一种基于匹配滤波的探地雷达全波形反演方法

Publications (1)

Publication Number Publication Date
CN111123374A true CN111123374A (zh) 2020-05-08

Family

ID=70502467

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201911359612.9A Withdrawn CN111123374A (zh) 2019-12-25 2019-12-25 一种基于匹配滤波的探地雷达全波形反演方法

Country Status (1)

Country Link
CN (1) CN111123374A (zh)

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111722287A (zh) * 2020-06-19 2020-09-29 南京大学 一种基于pda策略的震相特征识别波形反演方法
CN112084655A (zh) * 2020-09-08 2020-12-15 南京众诚土地规划设计咨询有限公司 一种基于非单调线搜索的探地雷达参数反演方法
CN113376629A (zh) * 2021-05-13 2021-09-10 电子科技大学 基于非均匀输入参数网格的井中雷达最小二乘反演方法
CN113447536A (zh) * 2021-06-24 2021-09-28 山东大学 一种混凝土介电常数反演与病害识别方法及系统

Cited By (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111722287A (zh) * 2020-06-19 2020-09-29 南京大学 一种基于pda策略的震相特征识别波形反演方法
CN111722287B (zh) * 2020-06-19 2021-10-08 南京大学 一种基于渐进数据同化方法的震相特征识别波形反演方法
CN112084655A (zh) * 2020-09-08 2020-12-15 南京众诚土地规划设计咨询有限公司 一种基于非单调线搜索的探地雷达参数反演方法
CN113376629A (zh) * 2021-05-13 2021-09-10 电子科技大学 基于非均匀输入参数网格的井中雷达最小二乘反演方法
CN113447536A (zh) * 2021-06-24 2021-09-28 山东大学 一种混凝土介电常数反演与病害识别方法及系统
CN113447536B (zh) * 2021-06-24 2022-09-30 山东大学 一种混凝土介电常数反演与病害识别方法及系统

Similar Documents

Publication Publication Date Title
CN111123374A (zh) 一种基于匹配滤波的探地雷达全波形反演方法
CN110779795B (zh) 裂缝性储层地质力学建模网格单元大小确定方法
CN107765302B (zh) 不依赖震源子波的时间域单频波形走时反演方法
CN104077451A (zh) 一种用于深厚软土地铁基坑土体参数反演分析的方法
CN110837669A (zh) 基于多源异构数据融合的滑坡不确定模型动态构建方法
CN109459787B (zh) 基于地震槽波全波形反演的煤矿井下构造成像方法及系统
CN108733952B (zh) 一种基于序贯模拟的土壤含水量空间变异性三维表征方法
CN103105622B (zh) 基于数据库技术的同型波时差定位方法
Wang et al. Determination of discontinuity persistent ratio by Monte-Carlo simulation and dynamic programming
CN115390033A (zh) 基于探地雷达路面下空洞及水损检测方法、系统及装置
CN112489208A (zh) 基于蚂蚁算法的裂缝片提取方法和三维地质模型构建方法
CN112084655A (zh) 一种基于非单调线搜索的探地雷达参数反演方法
CN117092702A (zh) 孔-隧激发极化探水结构的施工方法及反演探水方法
CN114236624B (zh) 基于电磁法估算压裂改造空间体积的方法和系统
CN111123373A (zh) 一种基于波场扩展重构的探地雷达全波形反演方法
CN111680380A (zh) 基于地质力学特征空间展布的全三维压裂设计方法
CN105319594A (zh) 一种基于最小二乘参数反演的傅里叶域地震数据重构方法
CN114492206A (zh) 一种基于破碎岩体节理模型确定隧道开挖进尺量的计算方法
CN110794469B (zh) 基于最小地质特征单元约束的重力反演方法
CN113031076A (zh) 一种在电阻率法超前探测中考虑围岩各向异性影响的方法
CN112285702A (zh) 一种基于lasso算法的探地雷达成像方法
CN109960776A (zh) 一种用于水力走时和水力信号衰减反演计算的改进算法
JP7495769B1 (ja) 高温地熱田の貯留層構造孔隙度の定量的特徴付けの方法及びシステム
RU2815465C1 (ru) Система полудетерминированного моделирования трещиноватости на основе беспорядочной матрицы
Li et al. Stability evaluation model for loess deposits based on PCA-PNN

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
WW01 Invention patent application withdrawn after publication
WW01 Invention patent application withdrawn after publication

Application publication date: 20200508