2. 地球信息探测仪器教育部重点实验室(吉林大学), 长春 130026;
3. 重庆大学土木工程学院, 重庆 400000;
4. 吉林省勘查地球物理研究院, 长春 130062
2. Key Laboratory of Geophysical Exploration Equipment, Ministry of Education(Jilin University), Changchun 130026, China;
3. College of Civil Engineering, Chongqing University, Chongqing 400000, China;
4. Jilin Geophysics Prospecting Institute, Changchun 130062, China
0 引言
在探地雷达(ground-penetrating radar, GPR)研究中，在不考虑磁导率等因素影响的前提下，阻抗参数可以通过经验公式估计介电常数、层速度等其他物性参数。利用阻抗信息可以获得不同深度的介电常数估计结果，对探地雷达参数反演估计有重要的研究意义[1]。RoyLindseth[2]早在20世纪70年代提出了阻抗反演技术并应用于地震勘探。探地雷达与地震方法有相同的物理机制，但探地雷达阻抗反演的研究起步较晚。国外学者[3]2012年首次将阻抗反演与探地雷达技术结合，从探地雷达反射数据中利用阻抗反演估计得到高分辨率的含水量分布；Li等[4-5]基于混合非均匀随机土壤模型开展了有限带宽阻抗参数反演，并通过实测数据验证得到了较好的效果；2018年，刘钰等[6]提出了结合GPR数据和钻孔资料来增强考古目标的成像和特征显示，改善了以往钻孔资料仅约束钻孔周围的情况；佘松盛等[7-8]提出利用探地雷达阻抗反演等方法来探测轻非水相污染物(简称LNAPL)。
在上述研究中，阻抗反演都需要利用钻孔介电常数资料作为约束项，约束项的精度直接影响最后的反演结果。此外，在横向变化较大的区域，钻孔资料往往只能提供单点信息，无法反映区域的地下介电常数信息。在缺少钻孔资料的实测数据中，如何获得不同深度的介电常数信息作为阻抗反演约束项是该方法研究的重要内容之一。
速度分析是估计地下速度结构分布的重要方法，通常采用多次覆盖的资料来进行速度分析以求取介电常数(速度)参数[9]。在探地雷达应用中，通过多偏移测量方式得到近似多次覆盖的数据来估计地下介质的介电常数[10-11]，进而通过经验公式转换得到含水量分布[12-13]。常规速度谱分析需要对不同深度的能量反射点进行拾取从而估计得到不同深度的速度结构。对于复杂地下结构以及海量探地雷达数据，常规人工拾取以及基于能量最大值的拾取方法在计算效率和精度方面存在一定的不足。K-means方法[14]是一种常用的非监督机器学习方法，是对数据集中在某些方面相似的数据成员进行分类组织的过程，简单、高效的应用特点使其在聚类算法中受到了广泛的关注。
本文提出了基于K-means方法速度谱估计作为阻抗反演约束项，开展探地雷达阻抗参数反演。首先介绍了阻抗反演以及基于K-means方法开展速度谱估计的基本原理，通过典型模型介绍阻抗反演方法的实现流程；在此基础上，通过复杂土壤模型测试对比了基于速度谱估计和常规钻孔约束阻抗反演方法的效果，并测试噪声适应性；最后，根据美国密歇根州前沃特史密斯空军基地(Wurtsmith AFB, in Oscoda, MI)的LNAPL区域多偏移距实测雷达数据[15-16]开展阻抗反演测试，验证本文方法在探地雷达实测数据中的应用效果。
1 基于速度分析的探地雷达阻抗反演原理 1.1 探地雷达阻抗反演原理通常对于无黏土的沙砾层，可以假设高频电磁波在非磁性和低损耗的介质中传播，此时电磁波阻抗与相对介电常数和速度的关系式分别为[17]：
(1) 
(2) 式中：Z为阻抗；Z0≈377 Ω，为真空中的阻抗；εr为相对介电常数；v为地下介质对应的电磁波传播速度；c=3×108 m/s，为真空中的电磁波传播速度。对于实际雷达数据：
(3) 式中：D(t) 为雷达记录；R(t) 为反射系数；W(t) 为雷达发射源子波；δ(t) 为噪声；⊗表示卷积操作；t为时间。
由于数据采集时受到噪声等因素的影响，需要对探地雷达数据进行预处理，如零时校正和滤波等处理。根据电磁波的衰减特性，对预处理雷达记录进行增益。本文采用指数增益方法恢复雷达信号的振幅能量。在获得增益补偿探地雷达数据后，利用反褶积方法求取雷达剖面的反射系数。由式(3)来求取反射系数是一个病态的问题[18]，本文采用稀疏反褶积方法求取反射系数[19-22]。利用求取出来的反射系数进行阻抗递推：
(4) 式中：Zn为第n层的阻抗；Rn为第n层介质与第n+1层介质界面处的反射系数。如果假设反射系数都小于0.3，式(4)可近似为
(5) 获得首层(地面)的阻抗，理论上就可以由式(4)、(5)来递推出任意层的阻抗，但反褶积方法难以获得准确的反射系数，导致递推出的阻抗幅度变化超出误差范围；因此，还需要引入约束项来约束阻抗结果。本文采用有限带宽阻抗(BLIMP)反演，其主要步骤如下[23]。
1) 由公式(4)(5)递推得到波阻抗剖面Za1。
2) 导入速度约束，若是单道约束则应扩充为与Za1水平相匹配的约束剖面Zw，将Za1时域上与Zw匹配。
3) 计算Zw的线性趋势值Zb，并从Zw减去Zb，结果称为Za2。
4) 分别对Za1和Za2做傅里叶变换到频率域。
5) 对频率域中的Za2做低通滤波，并选择合适阈值对Za1进行增益，在频率域将两者结果相加。
6) 对上述结果做傅里叶反变换，并将其与Zb`相加。
7) 根据速度约束对最终结果进行深时转换。
综上所述，图 1为探地雷达阻抗反演的流程图。
|
| 图 1 探地雷达阻抗反演流程图 Fig. 1 Workflow of GPR impedance inversion |
|
|
当地下介质为水平层状时，在宽角法测量方式下,
(6) 式中：Δti为第i个反射波正常时差；t0为零偏移距时的双程走时；xi为第i个收发天线间距；vrms为均方根叠加速度；N为接收天线数目。在探地雷达记录中，反射波双程走时满足
(7) 本文通过宽角法测量方式采集的雷达数据研究提取速度层约束的流程如图 2所示。
|
| 图 2 基于速度分析提取阻抗约束流程图 Fig. 2 Flow chart of extracting impedance constraint based on velocity analysis |
|
|
1) 首先对采集到的雷达数据进行预处理操作，包括去除直达波、滤波和零时校正等处理。
2) 对预处理完的数据进行线性增益，恢复振幅能量。
3) 进行速度分析。速度分析的流程是通过给定一个时间范围[t0min, t0max]和速度范围[vmin, vmax]，给定一个时间步长Δt和速度步长Δv，从t0min开始，对每一个速度按照式(7)进行扫描，然后对下一个t0min + Δt重复上述步骤，得到一个关于t0和vrms的网格，即速度谱。在速度谱中，与真实地下层状相对应的叠加速度处可以得到一个能量团，其中心位置处即真实叠加速度与双程旅行时。
4) 引入提取叠加速度和双程旅行时的K-means方法[14]。其基本流程如下：①设立k个簇族；②将k点放置于代表每个簇初始质心的空间中；③将每个点分配给最接近质心的聚类；④分配所有点后，重新计算每个聚类的中心；⑤重复③与④，直到收敛为止。此方法应对自动批处理海量探地雷达数据具有较好的效果(图 3)。
5) 得到各个层位的叠加速度和双程旅行时，利用Dix公式求取层速度[24]：
(8) 式中：vn为第n层速度；t0, n为到第n层底部的零偏移距双程走时；vrms, n为第n层底部的叠加速度。
6) 利用式(2)将层速度转化为介电常数约束项。
2 探地雷达阻抗反演模型测试 2.1 探地雷达阻抗反演流程图 4介绍了层状介质模型基于速度谱分析的探地雷达阻抗反演流程。首先建立如图 4a所示的四层层状介质模型，相对介电常数依次为5，10，16，25。采用时域有限差分法(FDTD)开展探地雷达正演计算，其天线为主频300 MHz的Ricker子波，分别获得共偏移距和某固定发射源的多偏移距探地雷达剖面。其中，共偏移距的参数设置为移动步长0.05 m，时间窗口90 ns，采样时间间隔0.045 2 ns。经过预处理和增益后的共偏移距剖面如图 4b所示。图 4c为采用稀疏反褶积估计得到的反射系数。多偏移距的发射天线固定在水平0.1 m的位置，接收天线间隔0.1 m布设。得到的多偏移距雷达剖面经过预处理和增益后如图 4d所示。采用速度谱分析和K-means方法自动提取层速度，并转化成阻抗约束项(图 4e)。基于反射系数与阻抗约束进行有限带宽阻抗反演可以得到图 4f。对比阻抗估计的介电常数剖面与原始模型，可以看到无论从空间分布上还是数值分布上基于速度分析的阻抗反演方法有较好的效果(图 4g)。
|
| a. 真实介电常数剖面；b. 共偏移距剖面；c. 反射系数；d. 多偏移距剖面；e. 多偏移距雷达速度谱；f. 反演阻抗；g. 反演介电常数。 图 4 速度谱约束探地雷达阻抗反演流程 Fig. 4 Workflow of GPR Impedance inversion constrained by velocity spectrum |
|
|
设计如图 5a所示的复杂土壤模型开展本文提出的速度谱约束和常规钻孔约束阻抗反演方法测试，该复杂模型尺寸为7.2 m × 9.0 m，相对介电常数呈随机分布，局部范围变化较大。共偏移方式采集参数如下：天线主频100 MHz，移动步长0.1 m，时间窗口280 ns，采样时间间隔0.039 9 ns。多偏移距采集参数如下：发射天线固定在水平0.1 m的位置，接收天线间距为0.2 m，天线主频、时窗大小和采样时间同共偏移距采集方式相同。对探地雷达数据进行加噪处理以测试阻抗反演方法在无噪声和强噪声数据条件下的反演效果差异(图 5b)。基于加噪前后速度分析提取出的约束和真实单道介电常数约束(取自初始模型水平0.5 m处)在介电常数范围内对比如图 5c。从图 5c中可以看出，速度约束能反映出地下介电常数层位的大致信息，但是细节部分却有所缺失，这是因为速度约束的分辨率较低，应对小范围内的参数变化难以得到好的效果，噪声在提取约束时有一定影响，但仍保留大致层状结构特征。不加噪声和加噪声阻抗反演介电常数结果分别如图 5d和图 5e，图 5f为常规钻孔参数约束探地雷达阻抗反演结果。对比图 5d、e、f可以看出，在复杂模型条件下，速度谱约束反演结果在细节方面相比于钻孔约束较差，这主要是因为我们采用单个多偏移距剖面所得速度谱估计难以刻画横向较大的介质变化情况。对比图 5d和图 5e反演结果表明，噪声对于反演结果的影响可以忽略。
|
| a. 真实介电常数剖面；b. 加噪前后单道数据(上面是加噪前，下面是加噪后)；c. 真实单道介电常数与加噪前后速度谱分析计算所得介电常数对比；d. 不加噪速度谱约束阻抗反演结果；e. 加噪速度谱约束阻抗反演结果；f. 钻孔资料约束阻抗反演结果。 图 5 土壤模型速度约束阻抗反演测试 Fig. 5 Soil model GPR impedance inversion with velocity constraint test |
|
|
图 6为沿垂向方向拾取单道原始介电常数(取自初始模型水平0.5 m处)和加噪前后反演介电常数结果(取自反演剖面水平0.5 m处)对比。反演结果在变化趋势和数值大小上与真实介电常数都吻合较好。但随着深度的增加，存在一定的估计误差，这主要是由阻抗反演本身反演精度所决定的。综合来看，在噪声因素干扰条件下，本文提出的方法仍然具有稳定的反演结果。
|
| 图 6 真实介电常数与加噪前后反演的介电常数单道对比 Fig. 6 Comparison between actual permittivity, estimated permittivity and estimated permittivity after noise addition |
|
|
实测探地雷达数据来自于美国密歇根州奥斯科达现已退役的沃特史密斯空军基地(Wurtsmith AFB, in Oscoda, MI)，该地区在历史上曾有很长时间作为消防培训场所，进行过许多消防训练演习，燃烧了数千加仑的喷气燃料和溶剂(1加仑约等于3.785 4 L)，存在大量非轻水相污染[15]。表 1为采集参数表，采集了202次多偏移距探地雷达数据，单个多偏移距剖面如图 7a。对其进行提取速度约束操作如图 7b，综合考虑图像中的有效信息，选择对前100 ns旅行时间内的数据进行处理，提取每个多偏移距的零偏移距数据组成共偏移距剖面进行预处理、增益等操作求取反射系数。对于202个多偏移距剖面进行提取速度约束的流程操作后，也相应地截取前100 ns的约束信息，得到了基于速度分析的阻抗约束剖面图 7c，然后进行有限带宽阻抗反演并转化成介电常数得到了图 7d。
|
| a. 单个多偏移距雷达剖面；b. 速度谱分析结果；c. 速度谱约束阻抗剖面；d. 介电常数阻抗反演结果。 图 7 实际探地雷达速度约束阻抗反演 Fig. 7 Real GPR velocity constrained impedance inversion |
|
|
为了证明反演结果的可靠性，采用同野外相同的采集参数对反演得到的介电常数模型进行正演模拟。首先将介电常数模型加一层空气层，对其正演结果如图 8a，与原始野外记录图 8b对比，反射信号的空间位置分布基本吻合，验证了基于速度谱约束的介电常数反演结果，其能很好地估计地下介质的分布情况。
|
| a. 正演结果；b. 实测数据。 图 8 反演介电常数模型所得雷达剖面与原始实测探地雷达剖面对比 Fig. 8 Comparison between the radar profile obtained by inversion of dielectric constant model and the original measured GPR profile |
|
|
结合前人[15]对该地区的解释结果分析，该地区浅地表处大致分为3层(图 9)：0~1 m处为混凝土地面层；1~4 m为沙土层；4 m以下为地下水层，整体沿测线逐渐下沉，推测在4 m以下主要为非轻水相污染。采用速度谱约束阻抗反演可以很好地勾画估计地下有机污染物的分布和扩散区域。
|
| 图 9 地下结构分层解释剖面 Fig. 9 Interpretation section of underground structure |
|
|
1) 本文提出了基于K-means方法的速度估计作为约束项开展探地雷达阻抗参数反演。
2) 通过理论模型和实测数据测试验证了采用速度谱估计结果可以代替钻孔资料作为阻抗反演的约束项。
3) 最终的介电常数估计结果与常规钻孔资料约束结果吻合。速度分析约束的探地雷达阻抗反演有很大的应用潜力和研究价值。
4) K-means方法可以在复杂速度谱结构中获得高精度的估计约束，为速度谱估计提供自动、高精度的拾取结果。根据处理结果对比常规探地雷达阻抗反演可以归纳几个优点：①依靠探地雷达自身系统提取约束项，真正实现了无损探测。②可以多次采集约束信息，拓宽整体反演结果分辨率，提高精度。③约束效果在大致层位上可以和测井资料相媲美，契合阻抗反演中约束的本质。同时，我们也需要指出阻抗反演方法需要局部介电常数分布作为约束项，在横向介质变化较大且缺少更多钻孔约束时，阻抗反演方法存在一定的误差。
| [1] |
Hansen T B, Johansen P M. Inversion Scheme for Ground-Penetrating Radar that Takes into Account the Planar Air-Soil Interface[J]. IEEE Transactions on Geoscience and Remote Sensing, 2000, 38(1): 496-506. DOI:10.1109/36.823944 |
| [2] |
张进铎. 地震解释技术现状及发展趋势[J]. 地球物理学进展, 2006, 37(2): 578-587. Zhang Jinduo. CurrentSituation and Development Trend of Seismic Interpretation Technology[J]. Progress in Geophysics, 2006, 37(2): 578-587. |
| [3] |
Schmelzbach C, Tronicke J, Dietrich P. High-Resolution Water Content Estimation from Surface-Based Ground-Penetrating Radar Reflection Data by Impedance Inversion[J]. Water Resources Research, 2012, 48(8): 64-69. |
| [4] |
Li J, Zeng Z, Chen L, et al. Estimation of Mixed Soil Water Content by Impedance Inversion of GPR Data[C]//Proceedings of the 15th International Conference on Ground Penetrating Radar. Brussels: IEEE, 2014: 149-154.
|
| [5] |
Zeng Z F, Chen X, Li J, et al. Recursive Impedance Inversion of Ground-Penetrating Radar Data in Stochastic Media[J]. Applied Geophysics, 2015, 12(4): 615-625. DOI:10.1007/s11770-015-0514-0 |
| [6] |
Liu Y, Shi Z, Wang B, et al. GPR Impedance Inversion for Imaging and Characterization of Buried Archaeological Remains: A Case Study at Mudu City Cite in Suzhou, China[J]. Journal of Applied Geophysics, 2018, 148: 226-233. DOI:10.1016/j.jappgeo.2017.12.002 |
| [7] |
佘松盛. 探地雷达方法对近地表LNAPL污染物的探测研究[D]. 长春: 吉林大学, 2019. She Songsheng. Study on Near-Surface the LNAPL Pollution by GPR Method[D]. Changchun: Jilin University, 2019. |
| [8] |
王焱, 鹿琪, 刘财, 等. 利用GPR天线-目标极化的瞬时属性分析方法探测LNAPL污染土壤[J]. 吉林大学学报(地球科学版), 2018, 48(2): 491-500. Wang Yan, Lu Qi, Liu Cai, et al. Using GPR Antenna-Target Polarization Instantaneous Attribute Analysis Method to Detect LNAPL Contaminated Soil[J]. Journal of Jilin University (Earth Science Edition), 2018, 48(2): 491-500. |
| [9] |
王春辉. 探地雷达方法测量近地表含水量及污染物探测研究[D]. 长春: 吉林大学, 2007. Wang Chunhui. Near Surface Water Content Measurement and Contamination Detection Using Ground-Penetrating Radar: A Simulation Study[D]. Changchun: Jilin University, 2007. |
| [10] |
Forte E, Pipan M. Review of Multi-Offset GPR Applications: Data Acquisition, Processing and Analysis[J]. Signal Processing, 2017, 132: 210-220. DOI:10.1016/j.sigpro.2016.04.011 |
| [11] |
Parsekian A D. Inverse Methods to Improve Accuracy of Water Content Estimates from Multi-Offset GPR Parsekian: Improved Accuracy of Water Content from GPR[J]. Journal of Environmental and Engineering Geophysics, 2018, 23(3): 349-361. DOI:10.2113/JEEG23.3.349 |
| [12] |
董泽君, 鹿琪, 冯晅, 等. 探地雷达测量土壤含水量的应用研究[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. |
| [13] |
蔡佳琪. 利用探测雷达探测铁路路基含水率[D]. 长春: 吉林大学, 2017. Cai Jiaqi. Detection Railway Subgrade Moisture Content by GPR[D]. Changchun: Jilin University, 2017. |
| [14] |
Chen Yuqing, Gerard S. Automatic Semblance Picking by a Bottom-Up Clustering Method[C]//SEG 2018 Workshop: SEG Maximizing Asset Value Through Artificial Intelligence and Machine Learning. Beijing: SEG, 2018: 27-30.
|
| [15] |
Bradford J H. GPR Offset Dependent Reflectivity Analysis for Characterization of a High-Conductivity LNAPL Plume[C]//16th EEGS Symposium on the Application of Geophysics to Engineering and Environmental Problems. San Abtonio: SEG, 2003: 238-252.
|
| [16] |
Bermejo J L, Sauck W A, Atekwana E A. Geophysical Discovery of a New LNAPL Plume at the Former Wurtsmith AFB, Oscoda, Michigan[J]. Groundwater Monitoring & Remediation, 1997, 17(4): 131-137. |
| [17] |
Annan A P. Ground-Penetrating Radar[M]//Near-Surface Geophysics. Tulsa: Society of Exploration Geophysicists, 2005: 357-438.
|
| [18] |
刘钰. 探地雷达数据波阻抗反演方法及其应用研究[D]. 杭州: 浙江大学, 2018. Liu Yu. The Study of Ground Penetrating Radar Impendance Inversion Method and Its Application[D]. Hangzhou: Zhejiang University, 2018. |
| [19] |
刘钰, 石战结, 王帮兵, 等. 探地雷达子波确定性稀疏脉冲反褶积技术[J]. 浙江大学学报(工学版), 2018, 52(9): 1828-1836. Liu Yu, Shi Zhanjie, Wang Bangbing, et al. Deterministic-Wavelet Sparse Spike Deconvolution Technique for Ground Penetrating[J]. Journal of Zhejiang University (Engineering Edition), 2018, 52(9): 1828-1836. |
| [20] |
康治梁, 张雪冰. 基于L1/2正则化理论的地震稀疏反褶积[J]. 石油物探, 2019, 58(6): 855-863. Kang Zhiliang, Zhang Xuebing. Seismic Sparse Deconvolution Based on L1/2 Regularization[J]. Geophysical Prospecting for Petroleum, 2019, 58(6): 855-863. |
| [21] |
王铁兴, 王德利, 孙婧, 等. 基于三维稀疏反演的混合震源数据分离与一次波估计[J]. 吉林大学学报(地球科学版), 2020, 50(3): 895-904. Wang Tiexing, Wang Deli, Sun Jing, et al. Separation and Primary Estimation of Blended Data by 3D Sparse Inversion[J]. Journal of Jilin University (Earth Science Edition), 2020, 50(3): 895-904. |
| [22] |
Li G F, Qin D H, Peng G X, et al. Experimental Analysis and Application of Sparsity Constrained Deconvolution[J]. Applied Geophysics, 2013, 10(2): 191-200. DOI:10.1007/s11770-013-0377-1 |
| [23] |
李静. 随机等效介质探地雷达探测技术和参数反演[D]. 长春: 吉林大学, 2014. Li Jing. Ground Penetrating Radar Detection and Parameter Inversion in Stochastic Effective Medium[D]. Changchun: Jilin University, 2014. |
| [24] |
Dix C H. Seismic Velocities from Surface Measurements[J]. Geophysics, 1955, 20(1): 68-86. DOI:10.1190/1.1438126 |