CN111123373A - 一种基于波场扩展重构的探地雷达全波形反演方法 - Google Patents

一种基于波场扩展重构的探地雷达全波形反演方法 Download PDF

Info

Publication number
CN111123373A
CN111123373A CN201911334761.XA CN201911334761A CN111123373A CN 111123373 A CN111123373 A CN 111123373A CN 201911334761 A CN201911334761 A CN 201911334761A CN 111123373 A CN111123373 A CN 111123373A
Authority
CN
China
Prior art keywords
inversion
model
wave field
ground penetrating
penetrating radar
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
CN201911334761.XA
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 CN201911334761.XA priority Critical patent/CN111123373A/zh
Publication of CN111123373A publication Critical patent/CN111123373A/zh
Withdrawn legal-status Critical Current

Links

Images

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V3/00Electric or magnetic prospecting or detecting; Measuring magnetic field characteristics of the earth, e.g. declination, deviation
    • G01V3/38Processing data, e.g. for analysis, for interpretation, for correction
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01NINVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N23/00Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01SRADIO 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/00Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
    • G01S7/02Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V3/00Electric or magnetic prospecting or detecting; Measuring magnetic field characteristics of the earth, e.g. declination, deviation
    • G01V3/12Electric or magnetic prospecting or detecting; Measuring magnetic field characteristics of the earth, e.g. declination, deviation operating with electromagnetic waves

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Remote Sensing (AREA)
  • General Physics & Mathematics (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • Geology (AREA)
  • Geophysics (AREA)
  • Environmental & Geological Engineering (AREA)
  • Health & Medical Sciences (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Electromagnetism (AREA)
  • Chemical & Material Sciences (AREA)
  • Analytical Chemistry (AREA)
  • Biochemistry (AREA)
  • General Health & Medical Sciences (AREA)
  • Immunology (AREA)
  • Pathology (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

本发明涉及一种基于波场扩展重构的探地雷达全波形反演方法,包括如下步骤:步骤一:建立探地雷达全波形反演的初始模型并设定反演终止精度;步骤二:构建新的扩展后的波动方程;步骤三:基于初始模型进行正演运算并抽取道集;步骤四:利用抽取的道集信息构建目标函数,其中,目标函数中的波场为扩展后的重构波场,在接下来的运算中可看作已知;步骤五:计算梯度和近似海塞矩阵,确定模型更新方向;步骤六:确定步长并更新模型;通过扩展波场使反演包含了更多的信息,进而使反演可以在更大的搜索空间内进行搜索,因此具有更大的自由度,从而不易陷入局部极小点。扩展后的波场中同时包含入射波和反射波的信息,不需计算伴随波场。

Description

一种基于波场扩展重构的探地雷达全波形反演方法
技术领域
本发明涉及一种基于波场扩展重构的探地雷达全波形反演方法。
背景技术
探地雷达是一种高效且具有较高分辨率的无损检测手段,可应用于隧道质量检测、路面密实度检测、混凝土构件非破坏性检测、地基夯实与加固检测、地灾评估和土壤地质信息评定等许多方面。探地雷达利用反射系数确定被检测体的结构及缺陷,但由于电磁波的传播过程极其复杂,人类活动对成像结果也有很大干扰,因此获得的地质雷达剖面往往不够准确,解释起来也与被检测体的真实结构存在差异,因此对雷达数据进行优化处理以期进行合理解释迫在眉睫。
得益于地震全波形反演的发展以及麦克斯韦方程组在解的形式上与声波方程的相似性,雷达数据可借鉴地震数据的处理方法,利用全波形反演进行数据的优化。Ernst(2006)实现了时间域探地雷达全波形反演,并对亚波长异常体进行了精准定位;Meles(2010)利用钻孔雷达全波形反演获得了河边较清晰的地下地层结构;胡周文(2018)利用全波形反演对混凝土进行了无损检测,冯德山(2018)等利用GPU并行实现了雷达数据的双参数反演,可见雷达数据的全波形反演是可行的。
然而受算法的限制,全波形反演严重依赖初始模型,初始模型精度不足极易导致全波形反演陷入局部极小点,进而反演失败;此外,全波形反演计算量巨大,虽然利用共轭梯度法可避免海塞矩阵的计算,但其计算仍旧很复杂。因此,如何提高探地雷达全波形反演的稳定性,使其不陷入局部极小点,并在保证反演精度的前提下简化反演过程,也是进行雷达数据全波形反演必须要考虑的关键问题。
探地雷达全波形反演技术经过二十几年的发展,目前已取得了长足的进步,但由于算法本身的缺陷,探地雷达全波形反演在应用时仍存在许多问题,针对本发明的研究内容,现提出其中两点:
第一,全波形反演易陷入局部极小点导致反演失败。使用全波形反演本是为了提高雷达图像的精度,但使用局部优化算法易导致全波形反演陷入局部极小点,反而会降低精度,因此,要想充分发挥全波形反演高精度的优势,必须克服周期跳跃现象,使反演结果收敛于全局极小点;
第二,全波形反演过程复杂,计算量大。全波形反演通过多次迭代逐渐使反演参数模型逼近真实参数模型,当其反演精度不足时,会自动计算梯度和海塞矩阵来迭代更新模型以期提高反演精度。然而每一次迭代过程均涉及较大的计算量,尽管利用共轭梯度法可以避免直接计算海塞矩阵,但由于反演数据本身体量巨大,因此近似海塞矩阵的求取依然很耗时。如能进一步简化梯度和海塞矩阵的求取过程,则反演的速度会得到进一步的提高,探地雷达的应用效果也将得到明显改善。
发明内容
为了解决这一问题,本发明提出一种基于波场扩展重构的探地雷达全波形反演方法,能够提高雷达数据反演剖面分辨率,进而提高雷达数据解释可靠性。
为解决上述技术问题,本发明所采用的技术方案为:一种基于波场扩展重构的探地雷达全波形反演方法,包括如下步骤:
步骤一:建立探地雷达全波形反演的初始模型并设定反演终止精度;
步骤二:构建新的扩展后的波动方程;
步骤三:基于初始模型进行正演运算并抽取道集;
步骤四:利用抽取的道集信息构建目标函数,其中,目标函数中的波场为扩展后的重构波场,在接下来的运算中可看作已知;
步骤五:计算梯度和近似海塞矩阵,确定模型更新方向;
步骤六:确定步长并更新模型;
步骤七:将步骤六中得到的模型作为新的初始模型,重复步骤三到六直到反演结束,获得最终的反演模型。
作为优选,所述步骤一中,以模型为例进行方法验证,真实模型已知,故初始模型是经平滑后获得的。
作为优选,所述步骤二中,扩展后的波动方程形如
Figure BDA0002330647290000031
其中Γ为拾取因子,λ为加权系数,L为正演算子,U为波场,Dcal为理论电磁波数据,S为场源,m1和m2分别为介电常数和电导率参数;解上述方程可得到扩展后的波场U,此时U为已知,未知量仅剩模型参数m1和m2,反演时两个参数分开计算,其中一个计算时另一个作为常量。
作为优选,所述步骤四中,目标函数形如
Figure BDA0002330647290000041
作为优选,所述步骤五中,介电常数梯度形如
Figure BDA0002330647290000043
电导率梯度形如
Figure BDA0002330647290000044
介电常数近似海塞矩阵形如
Figure BDA0002330647290000045
电导率近似海塞矩阵形如
Figure BDA0002330647290000046
作为优选,所述步骤六中,参数模型的更新形式为
Figure BDA0002330647290000042
式中γi为步长,mi表示介电常数和电导率的第i次迭代模型,步长的选择使用非单调线搜索方法。
本发明所达到的有益效果:本发明的基于波场扩展重构的探地雷达全波形反演方法,通过扩展波场使反演包含了更多的信息,进而使反演可以在更大的搜索空间内进行搜索,因此具有更大的自由度,从而不易陷入局部极小点。此外,扩展后的波场中同时包含入射波和反射波的信息,不需计算伴随波场,因此简化了近似海塞矩阵的求取。
附图说明
图1为本发明的方法流程示意图;
图2为本发明的模型试算初始模型、真实模型、反演模型;
图3为本发明的真实模型、反演模型的数值对比。
具体实施方式
下面结合附图对本发明作进一步描述。以下实施例仅用于更加清楚地说明本发明的技术方案,而不能以此来限制本发明的保护范围。
图1所示为本方法的流程示意图。在反演之初首先建立初始模型,设置相应的参数以及反演的终止精度;随后构建新的波动方程并在此基础上进行正演模拟以获得重构后的波场;构建新的考虑数据误差和方程误差的目标函数;对目标函数求导获得梯度和近似海塞矩阵,利用非单调线搜索获得步长;分别更新介电常数和电导率,得到更新后的参数模型。具体步骤如下:
步骤一:建立探地雷达全波形反演的初始模型并设定反演终止精度;本发明以模型为例进行方法验证,真实模型和初始模型如图二所示,其中初始模型为经真实模型平滑所得。
步骤二:类比地震全波形反演中的方程形式,传统的探地雷达全波形反演的波动方程可写作:
L(m1 m2)U=S
其中L为正演算子,U为波场,S为场源,m1和m2分别为介电常数和电导率参数。与传统全波形反演相比,本发明提供的方法对波场进行了扩展,扩展后的波动方程形如:
Figure BDA0002330647290000051
其中L、U、S、m1和m2的意义不变,Γ为拾取因子,λ为加权系数,Dcal为理论电磁波数据。对上述新方程求解则可得到扩展后的波场U,此时U为已知,未知量仅剩模型参数m1和m2,反演时两个参数分开计算,其中一个计算时另一个作为常量;
步骤三:基于初始模型进行正演运算并抽取道集。正演时采用基于交错网格的时间域有限差分法进行数值模拟,边界采用完全匹配层边界条件;
步骤四:利用抽取的道集等信息构建目标函数。新的目标函数同时考虑了数据误差和方程误差,其形式为:
Figure BDA0002330647290000061
方程中各符号意义与步骤二一致。目标函数中的重构波场U在接下来的运算中可看作已知;
步骤五:计算梯度和近似海塞矩阵,确定模型更新方向。梯度的计算为对目标函数F求一阶导数,求介电常数的梯度时可将电导率看作常量,同理,求电导率梯度时可将介电常数看作常量。
介电常数的梯度表达式为:
Figure BDA0002330647290000062
电导率的梯度表达式为:
Figure BDA0002330647290000063
求取近似海塞矩阵时,介电常数和电导率同样分开处理。
介电常数的近似海塞矩阵表达式为:
Figure BDA0002330647290000065
电导率的近似海塞矩阵表达式为:
Figure BDA0002330647290000064
步骤六:确定步长并更新模型。步长的大小利用非单调线搜索方法来确定。在反演过程中,反演的目标是使目标函数下降至误差允许范围,而并不要求目标函数严格单调下降。一味要求目标函数的单调下降会降低目标函数的收敛速度,增大计算量,因此本发明使用相对更高效的非单调线性搜索方法。
步长αi确定后,参数模型的更新形式可表示为:
Figure BDA0002330647290000071
其中,mi表示介电常数和电导率的第i次迭代模型,
Figure BDA0002330647290000072
分别为两个参数的梯度和近似海塞矩阵的通用表示,
步骤七:将步骤六中得到的模型作为新的初始模型,重复步骤三到六直到反演结束,由此可获得最终的反演模型。
为便于介绍基于波场扩展重构的探地雷达全波形反演方法,本发明使用一实施案例加以说明。
图2所示为验算所用的初始模型、真实模型和其反演结果,其中图a1为介电常数的真实模型,图a2为介电常数的初始模型,图a3为介电常数的反演模型,图b1为电导率的真实模型,图b2为电导率的初始模型,图b3为电导率的反演模型。模型被离散为30×40的网格,网格间距0.05m,实际深度为1.5m×2m。
图3所示为真实模型与反演模型的数值对比,其中图3,a中实线A表示介电常数的反演结果,实线B表示介电常数的真实模型,对比两条实线可知反演结果与真实模型数值上比较接近;图3,b中实线C表示电导率的反演结果,实线D表示电导率的真实模型,同样对比两条实线可发现反演结果与真实模型数值上比较接近。该数值对比图证明了本发明提供的反演方法的可行性。
以上仅是本发明的优选实施方式,应当指出,对于本技术领域的普通技术人员来说,在不脱离本发明技术原理的前提下,还可以做出若干改进和变形,这些改进和变形也应视为本发明的保护范围。

Claims (6)

1.一种基于波场扩展重构的探地雷达全波形反演方法,其特征在于:反演包括如下步骤:
步骤一:建立探地雷达全波形反演的初始模型并设定反演终止精度;
步骤二:构建新的扩展后的波动方程;
步骤三:基于初始模型进行正演运算并抽取道集;
步骤四:利用抽取的道集信息构建目标函数,其中,目标函数中的波场为扩展后的重构波场,在接下来的运算中可看作已知;
步骤五:计算梯度和近似海塞矩阵,确定模型更新方向;
步骤六:确定步长并更新模型;
步骤七:将步骤六中得到的模型作为新的初始模型,重复步骤三到六直到反演结束,获得最终的反演模型。
2.根据权利要求1所述的基于波场扩展重构的探地雷达全波形反演方法,其特征在于:所述步骤一中,以模型为例进行方法验证,真实模型已知,故初始模型是经平滑后获得的。
3.根据权利要求1所述的基于波场扩展重构的探地雷达全波形反演方法,其特征在于:所述步骤二中,扩展后的波动方程形如
Figure FDA0002330647280000011
其中Γ为拾取因子,λ为加权系数,L为正演算子,U为波场,Dcal为理论电磁波数据,S为场源,m1和m2分别为介电常数和电导率参数;解上述方程可得到扩展后的波场U,此时U为已知,未知量仅剩模型参数m1和m2,反演时两个参数分开计算,其中一个计算时另一个作为常量。
4.根据权利要求1所述的基于波场扩展重构的探地雷达全波形反演方法,其特征在于:所述步骤四中,目标函数形如
Figure FDA0002330647280000021
5.根据权利要求1所述的基于波场扩展重构的探地雷达全波形反演方法,其特征在于:所述步骤五中,
介电常数梯度形如
Figure FDA0002330647280000022
电导率梯度形如
Figure FDA0002330647280000023
介电常数近似海塞矩阵形如H=λ2μ2ω4diag(U)*diag(U);
电导率近似海塞矩阵形如H=-λ2μ2ω2diag(U)*diag(U)。
6.根据权利要求1所述的基于波场扩展重构的探地雷达全波形反演方法,其特征在于:所述步骤六中,参数模型的更新形式为
Figure FDA0002330647280000024
式中γi为步长,mi表示介电常数和电导率的第i次迭代模型,步长的选择使用非单调线搜索方法。
CN201911334761.XA 2019-12-23 2019-12-23 一种基于波场扩展重构的探地雷达全波形反演方法 Withdrawn CN111123373A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201911334761.XA CN111123373A (zh) 2019-12-23 2019-12-23 一种基于波场扩展重构的探地雷达全波形反演方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201911334761.XA CN111123373A (zh) 2019-12-23 2019-12-23 一种基于波场扩展重构的探地雷达全波形反演方法

Publications (1)

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

Family

ID=70502032

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201911334761.XA Withdrawn CN111123373A (zh) 2019-12-23 2019-12-23 一种基于波场扩展重构的探地雷达全波形反演方法

Country Status (1)

Country Link
CN (1) CN111123373A (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112084655A (zh) * 2020-09-08 2020-12-15 南京众诚土地规划设计咨询有限公司 一种基于非单调线搜索的探地雷达参数反演方法

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112084655A (zh) * 2020-09-08 2020-12-15 南京众诚土地规划设计咨询有限公司 一种基于非单调线搜索的探地雷达参数反演方法

Similar Documents

Publication Publication Date Title
CN113221393B (zh) 一种基于非结构有限元法的三维大地电磁各向异性反演方法
CN106405651B (zh) 一种基于测井匹配的全波形反演初始速度模型构建方法
Ji et al. Deep neural network-based permittivity inversions for ground penetrating radar data
CN106526674A (zh) 一种三维全波形反演能量加权梯度预处理方法
CN111123374A (zh) 一种基于匹配滤波的探地雷达全波形反演方法
CN111239819B (zh) 一种基于地震道属性分析的带极性直接包络反演方法
CN109459787B (zh) 基于地震槽波全波形反演的煤矿井下构造成像方法及系统
CN113552625A (zh) 一种用于常规陆域地震数据的多尺度全波形反演方法
CN111290019B (zh) 一种应用于最小二乘逆时偏移的l-bfgs初始矩阵求取方法
CN115390033A (zh) 基于探地雷达路面下空洞及水损检测方法、系统及装置
Kong et al. Novel hybrid method to predict the ground-displacement field caused by shallow tunnel excavation
CN111123373A (zh) 一种基于波场扩展重构的探地雷达全波形反演方法
CN112666612B (zh) 基于禁忌搜索的大地电磁二维反演方法
CN110187386B (zh) 一种自动快速识别地质构造的dtw地震体属性分析方法
CN112084655A (zh) 一种基于非单调线搜索的探地雷达参数反演方法
CN109725354B (zh) 各向异性速度建模方法及系统
CN116224265A (zh) 探地雷达数据反演方法、装置、计算机设备及存储介质
CN107945271B (zh) 基于地质块体追踪的三维压力场建模方法
CN113050163B (zh) 一种振幅、相位信息可调节的全波形反演方法及装置
CN113608258B (zh) 一种构建高分辨率波阻抗反演标签的自洽深度学习方法
CN110794469B (zh) 基于最小地质特征单元约束的重力反演方法
CN111208568B (zh) 一种时间域多尺度全波形反演方法及系统
Guo et al. Research on tunnel lining image target recognition method based on YOLOv3
CN113031072A (zh) 虚同相轴层间的多次波压制方法、装置及设备
CN107526102A (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
WW01 Invention patent application withdrawn after publication
WW01 Invention patent application withdrawn after publication

Application publication date: 20200508