基于共中心点道集约束的探地雷达波阻抗反演

戴前伟，宁晓斌，张 彬

基于共中心点道集约束的探地雷达波阻抗反演

戴前伟1,2，宁晓斌1，张 彬1,2

(1. 中南大学 地球科学与信息物理学院，湖南 长沙 410083； 2. 中南大学 有色金属成矿预测与地质环境监测教育部重点实验室，湖南 长沙 410083)

探地雷达(GPR)波阻抗反演是一种准确获取地下介质本征参数的有效方法，该方法依赖于测井资料提供的低频信息，而在GPR实际应用中，钻孔资料很少。为此，提出利用共中心点(CMP)速度分析为波阻抗反演提供大尺度纵向约束，实现在CMP速度分析结果的约束框架下，精细重构介质的介电参数信息。首先，以层状模型为算例，验证了CMP速度分析结果作为波阻抗反演的初始模型约束的可行性；在此基础上，开展了2个随机介质模型的波阻抗反演测试，反演结果的整体结构与模型接近，细微结构得到了较好的重构，与理论值的平相对误差为8.73%。结果表明，该方法在随机介质模型的探地雷达的波阻抗反演中更高效和经济，并且成像结果中包含着丰富的细节信息，在土壤介质其他物理参数估计中具有可行性和适用性。

探地雷达；波阻抗反演；共中心点；速度分析；参数估计

准确探测土壤的介电常数对水文地质调查、环境领域、工程领域、考古研究具有重要意义，例如，地下水含量、污染物运移扩散、管线渗漏、古墓挖掘等[1-4]。地球物理勘探是现代探测地下相关参数的一种先进科学技术，常规的钻孔取心和地球物理测井方法能够提供详细地质信息，但探测范围仅限井眼附近的有限区域[5]。探地雷达(Ground Penetrating Radar，简称GPR)是通过用高频电磁波来确定地下埋藏物内部或结构的无损探测技术，电磁波在介质中传播时，其传播路径、波场强度和波形随所通过介质的电性及几何形态而变化[6]。在探地雷达探测时，地下介质的电性参数影响着电磁波的传播，从而控制了探地雷达的回波响应。因此，获取这些介质的电性参数信息对探地雷达的研究具有重要价值。共中心点(Common Midpoint，简称CMP)速度分析是从地面探地雷达多偏移距数据中估计目标大尺度介电常数(速度分布)的常用方法，但其成像分辨率有限，难以刻画目标结构细节信息[7-9]。近年，利用走时、振幅和相位等雷达波信息的全波形反演(Full Waveform Inversion，FWI)能够获取高分辨率的地下物性参数分布[10-12]，但FWI计算花费的时间太多[13]。新兴的神经网络、遗传算法、粒子群算法等完全非线性算法，也因其需要在较庞大的解空间中寻找最优解，由此导致大量的计算量，而不适用于大范围探测数据的解释[14-16]。

在地震勘探中，狭义的地震反演是指从有限频带宽度的地震数据中恢复出宽带的波阻抗参数模型。R. O. Lindseth[17]开发的阻抗反演方法，其表示地震资料缺少低频数据，必须通过其他方法补充低频数据才可行；D. A. Cooke等[18]详细阐述了地震资料的广义线性反演方法，拉开了地震阻抗反演的序幕；随后，周竹生等[19]利用地质、地震和测井资料进行联合反演，克服了单一线性反演的缺陷；R. J. Ferguson等[20]基于声波测井技术提出有限带宽阻抗反演，该方法结合了测井资料纵向的强约束和地震横向分辨率高的优点，实现了全宽带资料的反演；C. Schmelzbach等[21]通过地震阻抗反演流程，实现了探地雷达高分辨率地下含水量的估计，并在数值模拟及实测数据中都验证了其方法的可行性；李静[22]利用阻抗反演方法实现了三维探地雷达复杂随机层状介质含水量等目标本征属性参数的提取；Li Jing等[23]利用探地雷达波阻抗反演估算了月球风化层的相对介电常数，得到了地下详细的结构，为了解月球地下结构提供了可靠的方法；刘钰[24]采用基于模型的宽带约束探地雷达阻抗反演方法进行古墓探测，经过苏州木渎古城的探测工作被验证是一种有效的方法，但其表示基于模型的反演需要建立较好的初始模型。初始模型的低频信息是否准确关乎重建探测目标本征属性参数的成败，地震勘探常常通过测井资料和叠后速度分析建立准确的初始速度模型。由于探地雷达在工程项目中钻孔资料极少，共中心点速度分析成为了探地雷达效仿地震勘探获取二维速度场，建立初始模型的常用方法。S. Busch等[25]利用GPR共中心点道集进行速度分析，将速度分析结果定义为全波形反演的初始模型，得到了可靠的反演结果。张彬[26]指出共中心点道集能够较方便地利用叠加速度谱求取不同深度的平均波速，能够为偏移处理提供初始速度模型，但在地下横行非均匀介质中效果不佳。J. H. Bradford等[27]提出利用偏移后的共中心点数据做反射层析成像估算雷达波速方法，提高了解译近地表水分布的能力。

探地雷达波阻抗反演是了解浅地表层状介质物性参数的有效方法，但是探地雷达工程项目中稀缺的钻孔资料无法为波阻抗反演建立准确的初始模型，补充所需的低频信息。为此，笔者开展基于共中心点道集约束的地面探地雷达波阻抗反演研究，依据采集的地面雷达数据，通过CMP速度分析技术获取地层的速度分布，为波阻抗反演建立初始模型，补充低频信息，获取地下目标介质的中尺度本征属性参数。

1 方法原理

1.1 速度分析原理

地下介质速度分布图是GPR数据解释的关键，其中，获取地下介质速度分布通过时深转换可以定位异常体和反射界面的位置，也可以对GPR剖面做偏移处理提高探测精度[28]。CMP速度分析的前提是假设地层水平，介质均匀且各向同性，则当偏移距不大时，反射信号满足双曲线时距方程：

当偏移距远小于地下目标体埋深时，叠加速度近似等于均方根速度，根据速度谱拾取叠加速度后可以使用迪克斯(dix)公式[31]求取层速度，表达式为：

得到各层层速度后，结合各层的双程走时便可计算得到各层的层厚度，第层层厚度计算公式为：

1.2 探地雷达波阻抗反演原理

因此，阻抗反演可以通过获得的相对介电常数转化成其他参数，如采用Topp，CRIM公式估计含水量等[33]。

基于上述假设，自激自收的地面探地雷达反射系数可以表示为[34]：

全频带阻抗值可以根据以下公式获得：

一维有限带宽阻抗方法反演流程步骤[35]可以表示为：

①雷达反射波能量在地下传播时会随着时间、几何扩散、吸收和传输损耗等因素而衰减，所以在计算反射系数前需要对振幅进行真振幅恢复。如若子波振幅未知，则可以利用先验信息通过式(8)估算振幅补偿因子；

②计算反射系数前需要通过预测反褶积和偏移等预处理技术提升GPR剖面的分辨率；

③通过稀疏脉冲反褶积从步骤b雷达数据中计算反射系数；

2 算例分析

为了验证基于共中心点道集约束的探地雷达波阻抗反演的效果，设置了3组试验。试验一为共中心点速度分析作为波阻抗反演的初始模型有效性试验；试验二利用一维模型正演模拟高信噪比信号，试验基于共中心点道集约束的GPR波阻抗反演对高信噪比资料的反演效果，并与测井约束的波阻抗反演对比；试验三以二维随机介质为模型，试验基于共中心点道集约束的GPR阻抗反演在随机介质中的反演效果。

2.1 速度分析

为了验证CMP速度分析程序的有效性，建立了如图1a所示的深1.5 m，水平距离2 m的三层水平均匀模型。其中，第一层介质的埋深为0.5 m，其相对介电常数为10；第二层介质的埋深为0.9 m，其相对介电常数为15；第三层介质的相对介电常数为20；假设三层介质的电导率均为0.001 S/m。

图1 层状模型正演速度分析

表1 速度分析结果的相对误差分析

2.2 一维随机介质模型

为了验证CMP资料约束下的探地雷达阻抗反演方法的有效性，在大尺度相对介电常数背景下加入随机噪声建立如图2a所示的一维随机介质模型，模型的深度为2.6 m，分为五层，其相对介电常数分别为2.9、7、9、11、13，网格数为260，空间步长为0.01 m，假设所有地层的电导率为0.001 S/m。

采用时域有限差分算法对该模型进行正演模拟，采用中心频率为500 MHz的雷克子波作为激励源，采样率为0.02 ns，时窗长度为52 ns，采用共偏移距的方式记录雷达数据。图2b为正演模拟预处理后的雷达波形数据。图2c为稀疏脉冲反褶积迭代20次计算真振幅数据得到的反射系数。图2d为CMP速度分析初始模型、常规测井初始模型和理论模型相对介电常数的对比；图2e黑线为理论模型相对介电常数，蓝线为基于CMP为初始模型的阻抗反演结果，红线为基于测井低频信息为初始模型的阻抗反演结果。分别将两种方法的相对介电常数反演结果与理论模型对比，其中，测井约束的阻抗反演估计相对介电常数结果的相对误差平均值为6.81 %，CMP约束的阻抗反演估计相对介电常数结果的相对误差平均值为8.73 %。从分析结果看出，虽然测井约束的阻抗反演结果比速度分析约束下的阻抗反演结果与理论模型更为接近，但速度分析约束下的阻抗反演的细节信息与模型有着相近的趋势，其结果能够为地下介质细节的刻画提供较好的分辨率；而且，运用该技术反演，采集成本更经济，处理快捷。通过算例分析，本文提出的基于共中心点道集约束的阻抗反演方法能够高效、准确地获取层状随机介质中尺度参数。

(a) 模型相对介电常数；(b) 预处理后仿真数据；(c) 反射系数；(d) 两种方法的初始模型与理论模型；(e) 两种方法的反演估计结果与理论模型

2.3 二维随机介质模型

图3 三层随机介质模型及有限差分正演结果

Fig.3 Three layer random soil medium model and FDTD forward modeling result

图4 三层随机介质模型阻抗反演结果

Fig.4 Impedance inversion results of a three layer random medium model

3 结论

a. 层状模型测试结果显示，与实际测试模型的参数相比，共中心点道集速度分析的上下层速度的相对误差分别为5.21 %和6.41 %，上下层厚度的相对误差分别为4.04 %和3.30 %，表明共中心点道集速度分析能较准确地获取地下大尺度的速度结构信息，证明了CMP速度分析作为GPR波阻抗反演初始模型的策略可行性。

b.通过随机介质模型的GPR波阻抗反演测试，获取了随机介质模型的相对介电常数信息，结果与测井约束的阻抗效果一致，为随机介质中GPR波阻抗反演提供了一种经济高效的大尺度约束方法。

c.复杂随机介质模型测试结果显示，相对介电常数的相对误差值为8.73 %，随机介质细微结构得到了较好的重构，证明了基于共中心点道集速度分析的GPR波阻抗反演方法的精确性，在复杂的土壤随机介质物理参数估计和反演中极具可行性和适用性。

d．在实际工作中，地层结构复杂多变，如何利用CMP速度分析建立准确的二维速度初始模型，以期提高反演效果，有待研究。

致谢：感谢吉林大学地球探测科学与技术学院李静副教授在阻抗反演理论及计算上给予的帮助和讨论；同时感谢两位匿名审稿专家提出的宝贵意见。在此一并致谢！

请听作者语音介绍创新技术成果等信息，欢迎与作者进行交流

[1] GREAVES J R，LESMES P D，LEE M J，et al. Velocity variations and water content estimated from multi offset，ground penetrating radar[J]. Geophysics，1996，61(3)：683–695.

[2] 王仙丽. 地下油类污染区探地雷达(GPR)探测能力研究[D]. 青岛：中国海洋大学，2011. WANG Xianli. Research on detectability of ground penetrating radar for oil contaminated site[D]. Qingdao：Ocean University of China，2011.

[3] 金鑫. 管线渗漏异常探地雷达数据的电场分量成像分析[J]. 煤田地质与勘探，2018，46(2)：159–163. JIN Xin. Study on electric field component imaging of leakage in pipeline using GPR[J]. Coal Geology & Exploration，2018，46(2)：159–163.

[4] 吴秋霜，王齐仁，皮海康. 水泥混凝土路面脱空的探地雷达图像特征分析[J]. 煤田地质与勘探，2018，46(4)：181–185. WU Qiushuang，WANG Qiren，PI Haikang. Analysis on the image features of ground penetrating radar for cavity of concrete pavement[J]. Coal Geology & Exploration，2018，46(4)：181–185.

[5] KOBR M，MAREŠ S，PAILLET F. Geophysical well logging：Borehole geophysics for hydrogeological studies：Principles and applications in hydrogeophysics[M]. Berlin：Springer，2005.

[6] 李大心. 探地雷达方法与应用[M]. 北京：地质出版社，1994. LI Daxin. Method and application of Ground Penetrating Radar[M]. Beijing：Geological Publishing House，1994.

[7] 刘四新，蔡佳琪，傅磊，等. 利用探地雷达精确探测铁路路基含水率[J]. 地球物理学进展，2017，32(2)：878–884. LIU Sixin，CAI Jiaqi，FU Lei，et al. Accurate detection of moisture content of subgrade by GPR[J]. Progress in Geophysics，2017，32(2)：878–884.

[8] 董泽君，鹿琪，冯晅，等. 探地雷达测量土壤含水量的应用研究[J]. 地球物理学进展，2017，32(5)：2207–2213. DONG Zejun，LU Qi，FENG Xuan，et al. Estimation of soil water content using ground penetrating radar[J]. Progress in Geophysics，2017，32(5)：2207–2213.

[9] 王洪华，戴前伟. 探地雷达有限元正演及介电参数反演[M]. 长沙：中南大学出版社，2016. WANG Honghua，DAI Qianwei. Ground Penetrating Radar finite element numerical simulation and dielectric parameter inversion[M]. Changsha：Central South University Publishing House，2016.

[10] 林朋，彭苏萍，卢勇旭，等. 基于共轭梯度法的全波形反演[J]. 煤田地质与勘探，2017，45(1)：131–136. LIN Peng，PENG Suping，LU Yongxu，et al. Full waveform inversion based on the conjugate gradient method[J]. Coal Geology & Exploration，2017，45(1)：131–136.

[11] 孟旭，刘四新，吴俊军，等. 时间域跨孔雷达全波形反演及实际应用[J]. 世界地质，2016，35(1)：256–263. MENG Xu，LIU Sixin，WU Junjun，et al. Full waveform inversion of time-domain cross-hole radar and its field application[J]. Global Geology，2016，35(1)：256–263.

[12] 冯德山，王珣. 基于GPU并行的时间域全波形优化共轭梯度法快速GPR双参数反演[J]. 地球物理学报，2018，61(11)：4647–4659. FENG Deshan，WANG Xun. Fast ground penetrating radar double-parameter inversion based on GPU-parallel by time-domain full waveform optimization conjugate gradient method[J]. Chinese Journal of Geophysics，2018，61(11)：4647–4659.

[13] MELES G A，GREENHALGH S A，GREEN A G，et al. GPR full waveform sensitivity and resolution analysis using an FDTD adjoint method[J]. IEEE Transactions on Geoscience & Remote Sensing，2012，50(5)：1881–1896.

[14] 刘敦文，徐国元，黄仁东，等. 一种基于神经网络的探地雷达信号解释研究[J]. 地球物理学进展，2004，19(1)：179–182. LIU Dunwen，XU Guoyuan，HUANG Rendong，et al. Study on a signal interpretation of GPR based on BP neural network[J]. Progress in Geophysics，2004，19(1)：179–182.

[15] 袁克阔. 粒子群算法改进及内变量本构模型参数反演[J]. 煤田地质与勘探，2017，45(2)：112–117. YUAN Kekuo. Improved particle swarm optimization and parameter inversion in internal variable constitutive model[J]. Coal Geology & Exploration，2017，45(2)：112–117.

[16] 王升，陈洪松，付智勇，等. 基于探地雷达的典型喀斯特坡地土层厚度估测[J]. 土壤学报，2015，52(5)：1024–1030. WANG Sheng，CHEN Hongsong，FU Zhiyong，et al. Estimation of thickness of soil layer on typical karst hillslopes using a ground penetrating radar[J]. Acta Pedologica Sinica，2015，52(5)：1024–1030.

[17] LINDSETH R O. Synthetic sonic logs：A process for stratigraphic interpretation[J]. Geophysics，1979，44(1)：3–26.

[18] COOKE D A，SCHNEIDER W A. Generalized linear inversion of reflection seismic data[J]. Geophysics，1983，48(6)：665–676.

[19] 周竹生，周熙襄. 宽带约束反演方法[J]. 石油地球物理勘探，1993，28(5)：523–536. ZHOU Zhusheng，ZHOU Xixiang. Band-constrained inversion[J]. Oil Geophysical Prospecting，1993，28(5)：523–536.

[20] FERGUSON R J，MARGRAVE G F. A simple algorithm for band-limited impedance inversion[R]. CREWES Research Report，1996，8(21)：1–10.

[21] SCHMELZBACH C，TRONICKE J，DIETRICH P. High-resolution water content estimation from ground-penetrating radar reflection data by impedance inversion[J]. Water Resources Research，2012，48(8)：8505.

[22] 李静. 随机等效介质探地雷达探测技术和参数反演[D]. 长春：吉林大学，2014. LI Jing. Ground penetrating radar detection and parameter in stochastic effective medium[D]. Changchun：Jilin University，2014.

[23] LI Jing，ZENG Zhaofa，LIU Cai，et al. A study on lunar regolith quantitative random model and lunar penetrating radar parameter inversion[J]. IEEE Geoscience and Remote Sensing Letters，2017，14(11)：1953–1957.

[24] 刘钰. 探地雷达数据波阻抗反演方法及其应用研究[D]. 杭州：浙江大学，2018. LIU Yu. The study of ground penetrating radar impedance inversion method and its application[D]. Hangzhou：Zhejiang University，2018.

[25] BUSCH S，KRUK J V D，BIKOWSKI J，et al. Quantitative permittivity and conductivity estimation using full waveform inversion of on ground GPR data[J]. Geophysics，2012，77(6)：79–91.

[26] 张彬. 探地雷达中的逆时偏移及速度估计[D]. 长沙：中南大学，2010.ZHANG Bin. Reverse time migration and velocity estimation of ground penetrating radar[D]. Changsha：Central South University，2010.

[27] BRADFORD J H. Measuring water content heterogeneity using multifold GPR with reflection tomography[J]. Vadose Zone Journal，2008，7(1)：184–193.

[28] 薛桂霞，邓世坤，刘秀娟. 逆时偏移在探地雷达信号处理中的应用[J]. 煤田地质与勘探，2004，32(1)：55–57. XUE Guixia，DENG Shikun，LIU Xiujuan. An application of reverse-time migration in the ground-penetrating radar data processing[J]. Coal Geology & Exploration，2004，32(1)：55–57.

[29] PARSEKIAN A D，SLATER L，NTARLAGIANNIS D，et al. Uncertainty in peat volume and soil carbon estimated using ground-penetrating radar and probing[J]. Soil Science Society of America Journal，2012，76(5)：1911–1918.

[30] 高妍. 几种速度分析方法的对比分析与应用[J]. 工程地球物理学报，2014，11(4)：436–440. GAO Yan. Comparison of several velocity methods and application[J]. Chinese Journal of Engineering Geophysics，2014，11(4)：436–440.

[31] TOLDI J. Velocity analysis without picking[J]. Geophysics，1989，54(2)：191–199.

[32] ANNAN A P. Ground Penetrating Radar，in near surface geophysics[M]. Tulsa：Society of Exploration Geophysicists，2005.

[33] HUSMAN J A，SPERL C，BOUTEN W，et al. Soil water content measurements at different scales：Accuracy of time domain reflectometry and ground-penetrating radar[J]. Journal of Hydrology，2001，245(1/2/3/4)：48–58.

[34] BERTEUSSEN K A，URSIN B. Approximate computation of the acoustic impedance from seismic data[J]. Geophysics，1983，48(10)：1351–1358.

[35] ZENG Zhaofa，CHEN Xiong，LI Jing，et al. Recursive impedance inversion of ground-penetrating radar data in stochastic media[J]. Applied Geophysics，2015，12(4)：615–625.

[36] 冯德山，戴前伟，何继善. 探地雷达的正演模拟及有限差分波动方程偏移处理[J]. 中南大学学报(自然科学版)，2006，37(2)：361–365. FENG Deshan，DAI Qianwei，HE Jishan. Forward simulation of ground penetrating radar and its finite difference method wave equation migration processing[J]. Journal of Central South University(Science and Technology)，2006，37(2)：361–365.

[37] 冯德山，戴前伟. 探地雷达时域多分辨法(MRTD)三维正演模拟[J]. 地球物理学进展，2008，23(5)：1621–1625.FENG Deshan，DAI Qianwei. Application of the multi- resolution time domain method in three dimensional forward simulation of ground penetrating radar[J]. Progress in Geophysics，2008，23(5)：1621–1625.

Common midpoint gather constraint-based impedance inversion of ground penetrating radar

DAI Qianwei1,2, NING Xiaobin1, ZHANG Bin1,2

(1. School of Geosciences and Info-Physics, Central South University, Changsha 410083, China; 2. Key Laboratory of Metallogenic Prediction of Nonferrous Metal and Geological Environment Monitoring, Ministry of Education, Central South University, Changsha 410083, China)

GPR impedance inversion is an effective method to obtain accurately the intrinsic parameters of the subsurface medium. This method relies on low frequency information provided by logging data, but it is rarely accompanied by drilling in practical applications. To solve the problem, the common midpoint velocity analysis is employed to provide more low frequency component information for impedance inversion of GPR to finely reconstruct the dielectric parameter of subsurface in the framework of velocity-constrained inversion technique. Firstly, a layer model is specifically set up as an example to verify the feasibility of the developed velocity- constrained inversion to regulate the initial model. Then, the impedance inversion test is performed on two random media models, result of mean relative error from the true model is 8.73%, which show that the overall structure is more consistent with the real model, and more importantly, the microstructure is also finely depicted. The proposed method is more efficient and economical in impedance inversion of ground penetrating radar with random medium model, with more detailed information contained in imaging results, the method is very feasible and applicable in the estimation of other physical parameters in soil investigation.

ground penetrating radar(GPR); impedance inversion; common midpoint(CMP); velocity analysis; parameter estimation

P631

A

10.3969/j.issn.1001-1986.2020.03.030

1001-1986(2020)03-0211-08

2019-10-11；

2020-01-27

国家自然科学基金项目(41704128，41874148)；湖南省自然科学基金项目(2018JJ3636)

National Natural Science Foundation of China(41704128，41874148)；Hunan Provincial Natural Science Foundation Project，China(2018JJ3636)

戴前伟，1968年生，男，湖南涟源人，博士，教授，从事电磁法理论及工程地球物理勘探的研究. E-mail：qwdai@csu.edu.cn

宁晓斌，1994年生，男，广西玉林人，硕士，从事探地雷达方面研究工作. E-mail：631118892@qq.com

戴前伟，宁晓斌，张彬. 基于共中心点道集约束的探地雷达波阻抗反演[J]. 煤田地质与勘探，2020，48(3)：211–218.

DAI Qianwei，NING Xiaobin，ZHANG Bin.Common midpoint gather constraint-based impedance inversion of ground penetrating radar[J]. Coal Geology & Exploration，2020，48(3)：211–218.

(责任编辑 聂爱兰)