利用GPR天线-目标极化的瞬时属性分析方法探测LNAPL污染土壤

王 焱, 鹿 琪,2,刘 财，佘松盛，刘四新

1.吉林大学地球探测科学与技术学院，长春 130026 2.油页岩地下原位转化与钻采技术国家地方联合工程实验室，长春 130026

0 引言

轻非水相流体(light non-aqueous phase liquid, LNAPL),指不溶于水且密度小于水的液体。石油类为常见的LNAPL，它在开采、运输、加工以及使用的过程中不可避免地会发生泄漏，严重污染土壤和地下水。应用探地雷达(ground penetrating radar, GPR)探测地下石油类污染物是近十几年地球物理探测浅层污染中的热门课题。探地雷达相比较于传统的打井方法具有原位无损、操作简单、图像清晰直观、节约成本、抗干扰性强等优点[1]。应用探地雷达探测土壤污染物成为了国内外专家学者关注的焦点。Benson[2]在亚利桑那州和犹他州进行实验研究，结果表明了GPR信号和碳氢化合物污染之间存在良好相关性。Atekwana等[3]应用探地雷达与原位电阻率探测相结合，表征出LNAPL生物降解相关的分布区域，并对其进行了污染物监测和修复。Cassidy[4]通过探地雷达信号衰减分析和介电性能检测来评估LNAPL污染物。Daniels等[5]进行了GPR探测液体污染物实验。Hwang等[6]对小体积释放的重非水相流体进行了长期探地雷达监测。对于LNAPL污染物的探地雷达四维探测，Castro[7]和Bertolla等[8]分别进行了实际泄漏地点的探测。Bano等[9]对砂箱中柴油泄漏进行GPR探测并进行了时间域建模。Sauck[10]建立了LNAPL羽状体的电阻率结构模型。尤志鑫等[11]通过模拟轻非水相液体污染物污染土壤，根据污染物及水含量对石英砂介电常数的影响进行了研究。李晔等[12]研究了不同速率点源泄漏对LNAPL在土壤包气带中迁移的影响，建立了LANPL泄漏的数值计算模型。任磊[13]针对石油勘探开发中的石油类污染进行了探地雷达监测。侯晓东等[14]从探地雷达在土壤污染中的应用、污染土壤介电常数模型的演化和土壤介电常数求取三方面阐述了探地雷达技术在土壤污染中的研究。王春辉[15]通过提取污染土壤的介电常数，并与污染物的含量建立对应关系，进行了定量评估污染情况。武晓峰等[16]对LNAPL透镜体厚度及最大厚度进行分析并给出了关系式。郭秀军等[17]对不同土壤中含油污染区的电性变化进行了研究。

实测GPR数据具有高频、短波、强介质吸收和抗干扰等特点，以傅里叶变换为基础的信号处理方法往往不能达到理想效果，一些学者尝试应用时频分析方法对GPR数据进行处理。王宪楠等[18]利用Shearlet变换对GPR随机噪声进行压制。张先武等[19]将Hilbert谱分析法应用到GPR薄层勘探中，该方法对薄层分界面有较好的增强，但同时干扰信号也被增强。冯德山等[20]以Hilbert变换为基础，提取了实测GPR信号的瞬时参数信息，并分别绘制出各个瞬时参数对应的剖面图，相互对照得出了更好的解释。

极化在探地雷达中的应用方面，许多国内外学者已经做了一些工作：Tsoflias等[21]利用探地雷达中的极化成分探测地质构造中的垂直断裂；Seol等[22]通过利用探地雷达中的极化分量寻找裂缝的走向；Radzevicius等[23]也对GPR极化特性做过研究，指出圆柱体走向和线性偶极子天线长轴方向平行或垂直时，其后向散射存在差异。通常情况下利用商用雷达探测时，沿相同方向测线进行天线方位固定的测量，但当天线相对于地下目标物的方位不同时，存在天线与目标物之间的极化，对于这种极化的分析与解释有助于提高GPR对目标物的探测与辨识能力。Maurizio Lualdi等[24]对4组形状规则的目标进行天线-目标极化探测，证实了天线-目标极化特征对目标物几何形状和物性分析是十分重要的。Tsoflias等[25]提出用多极化采集测量方式检测垂直裂缝，并基于多极化GPR数据之间的相位差来确定垂直裂缝的位置和方位角。

本实验利用探地雷达对该污染土壤分别进行天线-目标极化探测。所谓的天线-目标极化，指天线布设方位与目标物方位存在一定角度，在本实验中将天线相对于目标物的方位设定为0°和90°。为了对比分析这两种天线-目标极化方式下的振幅、相位、频率等特征，对预处理和偏移后的数据进行Hilbert变换，求出瞬时振幅、瞬时相位、瞬时频率，将归一化的瞬时振幅分别与瞬时相位、瞬时频率进行叠合运算，得到新的复合属性。结合瞬时属性和复合属性，分析在探测同一目标物时天线方位对结果的影响，并对LNAPL污染土壤范围进行讨论。基于对天线-目标极化的瞬时属性的分析，旨在揭示LNAPL土壤污染分布规律。

1 方法原理

1.1 天线-目标极化分析

对于电磁散射问题，当电磁波入射到某一目标时，根据麦克斯韦方程和相应的电磁场的边界条件，在目标体上和目标体内便有电流磁流流动，这些感应电磁流产生二次辐射，这个辐射就是目标的散射场，并沿着不同的方位以不同的幅度和相位传播。当物体被电磁波照射时，能量将朝各个方向散射，能量的空间分布依赖于物体的形状、大小和结构以及入射波的频率和特征[26]。

极化在光学中称为偏振，所谓的极化是指在空间任一固定点上电磁波的电场强度矢量的空间取向随时间的变化方式，以E的矢端轨迹来描述。地下目标体可以导致反射波的极化方向与入射波的极化方向不同。若利用GPR探测LNAPL土壤污染时，布设水平和垂直方向相互垂直的测线，分别对应0°和90°方位天线测量，受目标体形状、大小、结构的影响，天线-目标极化不同，使得探测结果存在差别。由于LNAPL的密度小于水且不溶于水，可以设想一旦泄露到达地下水面后，将会在地下水面上横向扩展，漂浮在地下水面上，形成LNAPL的透镜体。由典型地下有机污染物渗漏和迁移模型可知，散落在地表的LNAPL在地下的分布主要分为3部分：受毛管力和重力共同作用下的残留相；受重力作用漂浮在潜水面上形成油池的自由相；部分溶于水中的溶解相[27-28]。由此可见，LNAPL土壤污染物形状、大小、结构并不规则，因此可以通过0°和90°方位天线测量来分析天线-目标极化特征，了解目标物的走向、形状等地下介质信息。

1.2 瞬时属性基本原理

GPR信号瞬时属性分析的基本原理是将实信号u(t)看成复信号z(t)的实部，通过Hilbert变换得到与u(t)相对应的复信号的虚部v(t),然后从复信号中提取瞬时振幅、瞬时相位、瞬时频率等特征参数，数学表达式如下：

z(t)=u(t)+iv(t)。

(1)

其中，

v(t)=hilbert[u(t)]。

(2)

则有瞬时振幅为

(3)

瞬时相位为

(4)

瞬时频率为

(5)

瞬时振幅是反射强度的度量，主要反映能量上的变化；瞬时相位描述复信号实部和虚部组成的矢量角度的变化，它是同相轴连续性的度量，无论能量强弱，异常或分层在相位图中都会体现；瞬时频率是相位随时间的变化率或者说是瞬时相位的导数[27]。

1.3 复合属性分析

瞬时属性可以将隐藏在原始GPR数据中的某些信息凸显出来。将瞬时属性的叠合、差值、乘积属性称为复合属性，它们是由瞬时属性再导出的新的属性[28]。对归一化后的瞬时属性进行和、差、积等运算，即可得到复合属性[29]。本文就是将归一化后瞬时振幅分别与瞬时相位和瞬时频率进行叠合运算，得到新的复合属性。

2 数据采集与处理

本实验在实验室砂箱中进行(图1)，采用水和石英砂混合物模拟土壤，柴油作为LNAPL污染物，砂箱长0.53 m、宽0.39 m、高0.32 m，砂箱中潜水面深度为0.24 m。测量仪器使用天线中心频率为1 600 MHz的MALA探地雷达。本次实验共布设60条测线，水平方向30条测线，垂直方向30条测线，测线间距为0.01 m，测线长度为0.30 m(图2)。图2a的模型主视图指的是从测量模型的前面向后面所看到的视图，图2b的模型俯视图指的是从测量模型的上方向下所看到的视图。布置两组相互垂直的天线，其中T1、R1为0°方位天线组合，T2、R2为90°方位天线组合，如图2所示。为了使探地雷达在行进过程中能沿测线准确测量，首先将塑料米尺固定在测线旁，尺边和待测测线的距离等于雷达宽度的一半，然后将探地雷达的侧边与尺边贴合，从起点处开始测量。由此，固定了探地雷达的行进轨迹。

a.模型主视图；b.模型俯视图。 图2 测量模型示意图Fig.2 Experiment setting

数据处理流程如图3所示。首先将与天线方位相同的30条测线测得的数据组合为一个三维数据体，0°方位天线采集的三维数据体称为A，90°方位天线采集的三维数据体称为B。其次，从A和B数据体中找出测线方向相同，天线方位垂直的对应剖面，进行包括帯通滤波、去除平均值在内的预处理，然后为了更好地反映异常的形态与位置进行三维偏移。通过Hilbert变换公式(3)(4)(5)得到0°和90°方位天线探测时的瞬时振幅(Ai、Ac)、瞬时相位(Pi、Pc)、瞬时频率(Fi、Fc)，再计算出振幅和(Ai+c)、相位差(Pi-c)、频率差(Fi-c)，计算公式如下：

Ai+c=Ai+Ac，

(6)

Pi-c=Pi-Pc，

(7)

Fi-c=Fi-Fc。

(8)

图3 数据处理流程图Fig.3 Processing flowchart for GPR data

3 实验结果分析

为了验证三维数据体A和B是否具有良好相关性，提取A和B中同一位置的原始单道振幅进行对比，结果如图4所示，可见数据体A和B具有很好的相关性。

图4 单道原始数据振幅对比图Fig.4 Comparison of the traces acquired at the same survey point by the inline and crossline survey

从三维数据体A和B中提取出水平第15条测线对应的原始剖面以及预处理后剖面，如图5、6所示。图5a、6a分别是0°和90°方位天线探测的原始剖面图，从图中可以看到：双程走时3.6 ns和6.0 ns处两条明显的水平同相轴，分别对应潜水面和砂箱底部反射；图5b和6b分别是0°和90°方位天线探

测的预处理后剖面图，图5b中双程走时2.0～4.0 ns、水平位置0.05～0.20 m可见同相轴不连续，判断为油污染异常，潜水面的水平同相轴依然可见。由于模型大小的局限性，有较强的边界反射现象。此外，由于天线相对于异常体方向不同，图6b与5b的异常显示结果存在差异，90°方位天线探测到6.0～8.0 ns较为明显的砂箱底部边界效应。

为了更好地反映油污染异常的位置与形状，对预处理后数据进行三维偏移，结果如图7所示。对于0°方位天线测量结果，图5b与图7a比较，由于原始资料潜水面处的砂箱边界反射较弱，因此潜水面在偏移后更清楚。对于90°方位天线测量结果，图6b与图7b比较，由于原始资料潜水面处的砂箱边界反射较强，干扰了对潜水面的偏移，使得偏移后的潜水面不如偏移前清楚。但是LNAPL的位置在砂箱中部，受边界影响弱，因此得到了较好的偏移结果。

接下来进行Hilbert变换，瞬时振幅如图8、9所示。其中8a和8b分别为0°和90°方位天线测得的瞬时振幅三维切片图，图中3个切片的位置分别为0.13 m(垂直距离) , 0.15 m(水平距离), 0.12 m(深度)，图8c和8d分别为0°和90°方位天线的0.15 m(水平距离)瞬时振幅切片，它们反映了能量上的变化。由于GPR对垂直于天线的目标物探测效果明显，对平行于天线的目标物探测效果不明显，而图8c中砂箱侧边界垂直于0°方位天线，因此，图8c中深度0.10 m的侧边界反射能量、0.30 m的侧边界反射与砂箱底边界反射的叠加能量都很强，0.10～0.30 m的侧边界能量存在但是相对很弱；图8d中砂箱侧边界平行于90°方位天线，图中几乎没有砂箱侧边界反射能量。为增强LNAPL油污染异常的响应，对两者求和得到瞬时振幅之和如图9所示,其中9b为0.15 m(水平距离)切片。通过对比可以看出，求和之后油污染区域能量聚集相对更加明显，0.10～0.30 m原本就弱的侧边界反射能量更弱，图9b中0.10 m处边界能量较图8c虽然也减弱，但仍大于0.10～0.30 m的侧边界反射能量。

a.原始数据；b.预处理后数据。图5 0°方位天线剖面图(测线15)Fig.5 Vertical profile detectedby 0 azimuth antenna(Line 15)

a.原始数据；b.预处理后数据。图6 90°方位天线剖面图(测线15)Fig.6 Vertical profile detected by 90 azimuth antenna(Line 15)

a.0°方位天线；b.90°方位天线。图7 偏移后剖面图(测线15)Fig.7 Vertical profile after migration(Line 15)

进行Hilbert变换后，同时还得到0°和90°方位天线测得的瞬时相位、瞬时频率；然后利用公式(7)、(8)计算求出瞬时相位差(Pi-c)、瞬时频率差(Fi-c)；最后，将归一化后的瞬时振幅和(A*i+c)、瞬时相位差(P*i-c)、瞬时频率差(F*i-c)分别编码为RGB成像方式中的红、蓝、绿3个颜色通道，这样可以综合考虑0°和90°天线方位探测的瞬时振幅、瞬时相位和瞬时频率，重构地下目标图像。红色代表0°和90°天线方位探测的瞬时振幅叠和较为突出的区域，红色区域中两种天线-目标极化方式都具有较强能量。对于横向分布不均匀介质，0°和90°天线方位探测存在相位差[30]，蓝色代表0°和90°天线方位探测的瞬时相位差值较大的区域，绿色代表0°和90°天线方位探测的频率差值较大的位置。从RGB成像结果(图10)可以看到：深度0.10～0.20 m水平居中位置处，红色突出且蓝色褪去，即瞬时振幅和大，瞬时相位差小，代表此处油污染聚集；深度在0.30 m左右层状区域红色和蓝色均连续，此处代表砂箱底界面。

a.三维显示；b.剖面显示。图9 瞬时振幅之和Fig.9 Summation of instantaneous amplitude

图10 复合属性RGB显示Fig.10 RGB image of the composite characteristics

4 结论

1)对于LANPL这类形状不规则目标物，存在天线-目标极化。通过布设不同方位天线，可以获取天线-目标极化信息，从而更准确分析LNAPL土壤污染在不同方向上的分布。

2)对于不同天线-目标极化方式的探测结果，可以通过三瞬分析更好地进行对比，其中瞬时振幅求和能更全面地观察土壤污染物的分布。

3)将0°和90°天线方位探测结果求取瞬时振幅之和、瞬时相位之差、瞬时频率之差，并将瞬时相位之差、瞬时频率之差与瞬时振幅叠合形成复合属性。复合属性能够表现出原来瞬时振幅无法表现的较深层信息，同时复合属性还使得瞬时相位、瞬时频率所涵盖的信息更好地体现。

[1] 曾昭发，刘四新，王者江，等.探地雷达方法原理及应用[M].北京：科学出版社，2006.

Zeng Zhaofa, Liu Sixin, Wang Zhejiang, et al. Method and Application of Ground-Penetrating Radar[M].Beijing:Science Press, 2006.

[2] Benson A K. Applications of Ground Penetrating Ra-dar in Assessing Some Geological Hazards: Examples of Groundwater Contamination, Faults, Cavities[J]. Journal of Applied Geophysics, 1995, 33(1/2/3): 177-193.

[3] Atekwana E A, Sauck W A, Werkema D D. Investi-gations of Geoelectrical Signatures at a Hydrocarbon Contaminated Site[J]. Journal of Applied Geophysics, 2000, 44(2): 167-180.

[4] Cassidy N J. Evaluating LNAPL Contamination Using GPR Signal Attenuation Analysis and Dielectric Property Measurements: Practical Implications for Hydrological Studies[J]. Journal of Contaminant Hydrology, 2007, 94(1): 49-75.

[5] Daniels J J, Roberts R, Vendl M. Ground Penetrating Radar for the Detection of Liquid Contaminants[J]. Journal of Applied Geophysics, 1995, 33(1/2/3): 195-207.

[6] Hwang Y K, Endres A L, Piggott S D, et al. Long-Term Ground Penetrating Radar Monitoring of a Small Volume DNAPL Release in a Natural Groundwater Flow Field[J]. Journal of Contaminant Hydrology, 2008, 97(1): 1-12.

[7]Castro D L D. 4-D Ground Penetrating Radar Moni-toring of a Hydrocarbon Leakage Site in Fortaleza (Brazil) During Its Remediation Process: A Case History[J]. Journal of Applied Geophysics, 2003, 54(1):127-144.

[8] Bertolla L, Porsani J L, Soldovieri F, et al. GPR-4D Monitoring a Controlled LNAPL Spill in a Masonry Tank at USP, Brazil[J]. Journal of Applied Geophysics, 2014, 103: 237-244.

[9] Bano M, Loeffler O, Girard J F. Ground Penetrating Radar Imaging and Time-Domain Modelling of the Infiltration of Diesel Fuel in a Sandbox Experiment[J]. Comptes Rendus Geoscience, 2009, 341(10): 846-858.

[10] Sauck W A. A Model for the Resistivity Structure of LNAPL Plumes and Their Environs in Sandy Sediments[J]. Journal of Applied Geophysics, 2000, 44(2): 151-165.

[11] 尤志鑫，冯晅，鹿琪. LNAPL污染物及水含量对石英砂介电常数的影响[J].世界地质，2015， 34(2)：551-556.

You Zhixin, Feng Xuan, Lu Qi. Influence of LNAPL Contamination and Water Content to Dielectric Constant of Quartz Sand[J]. Global Geology, 2015, 34(2): 551-556.

[12] 李晔，鹿琪，刘财. LNAPLs 迁移的数值模拟和土壤介电性质的变化分析[J].地球物理学进展，2014，29(2)：936-943.

Li Ye, Lu Qi, Liu Cai. Numerical Simulation of LNAPLs Migration and Analysis on Variation of Soil Dielectric Properties[J]. Progress in Geophysics, 2014, 29(2): 936-943.

[13] 任磊.石油勘探开发中的石油类污染及其监测分析技术[J]. 中国环境监测，2004，20(3)：44-47.

Ren Lei. The pollution Status and Determination Methods of Petroleum Contaminations During the Processes of oil Exploration[J]. Environmental Monitoring in China, 2004, 20(3): 44-47.

[14] 侯晓冬，郭秀军，贾永刚，等.基于探地雷达回波信号获取污染土壤中污染物含量的研究进展[J].地球物理学进展，2008，23(3)：962-968.

Hou Xiaodong, Guo Xiujun, Jia Yonggang, et al. Progress on Detection of Contamination Content in the Contaminated Soil Based on Ground Penetrating Radar[J]. Progress in Geophysics, 2008, 23(3): 962-968.

[15] 王春辉.探地雷达方法测量近地表含水量及污染物探测研究[D].长春：吉林大学，2007.

Wang Chunhui. Near Surface Water Content Measurement and Contamination Detection Using Ground-Penetrating Radar: A Simulation Study[D]. Changchun: Jilin University, 2007.

[16] 武晓峰，唐杰，藤间幸久.地下水中轻质有机污染物(LNAPL)透镜体研究[J].环境污染与防治，2000，22(3): 17-20.

Wu Xiaofeng, Tang Jie, Yukihisa Fujima Y. Study on Light Non-Aqueous Phase Liquid Lens in Subsurface Water[J]. Environmental Pollution & Control, 2000, 22(3): 17-20.

[17] 郭秀军，武瑞锁，贾永刚，等.不同土壤中含油污水污染区的电性变化研究及污染区探测[J].地球物理学进展，2005，20(2)：402-406.

Guo Xiujun, Wu Ruisuo, Jia Yonggang, et al. The Study of Electrical Resistivity Change of Different Saturation Soils Contaminated with Oil Sewage and the Contaminated Area Detecting[J]. Progress in Geophysics, 2005, 20(2): 402-406.

[18] 王宪楠，刘四新，程浩.Shearlet变换在GPR数据随机噪声压制中的应用[J].吉林大学学报(地球科学版)，2017，47(6)：1855-1864.

Wang Xiannan, Liu Sixin, Cheng Hao. Application of Shearlet Transform for Suppressing Random Noise in GPR Data[J]. Journal of Jilin University(Earth Science Edition), 2017, 47(6): 1855-1864.

[19] 张先武，高云泽，方广有. Hilbert谱分析在探地雷达薄层识别中的应用[J]. 地球物理学报，2013，56(8)：2790-2798.

Zhang Xianwu, Gao Yunze, Fang Guangyou. Application of Hilbert Spectrum Analysis in Ground Penetrating Radar Thin Layer Recognition[J]. Chinese Journal of Geophysics, 2013, 56(8): 2790-2798.

[20] 冯德山，戴前伟，余凯. 基于经验模态分解的低信噪比探地雷达数据处理[J]. 中南大学学报(自然科学版)，2012，43(2)：596-604.

Feng Deshan, Dai Qianwei, Yu Kai. GPR Signal Processing Under Low SNR Based on Empirical Mode Decomposition[J]. Journal of Central South University, 2012, 43(2):596-604.

[21] Tsoflias G, Van Gestel J P, Stoffa P, et al. Detection of Vertical Fractures in Geologic Formations Using the Polarization Properties of Ground-Penetrating Radar Signal[C]//69th Annual International Meeting. Houston: SEG, 1999: 567-570.

[22] Seol S J, Kim J H, Song Y, et al. Finding the Strike Direction of Fractures Using GPR[J]. Geophysical Prospecting, 2001, 49(3): 300-308.

[23]Radzevicius S J, Daniels J J. Ground Penetrating Radar Polarization and Scattering from Cylinders[J]. Journal of Applied Geophysics, 2000, 45(2): 111-125.

[24]Lualdi M, Lombardi F. Orthogonal Polarization Approach for Three Dimensional Georadar Surveys[J]. NDT & E International, 2013, 60: 87-99.

[25]Tsoflias G P, Van Gestel J P, Stoffa P L, et al. Vertical Fracture Detection by Exploiting the Polarization Properties of Ground-Penetrating Radar Signals[J]. Geophysics, 2004, 69(3): 803-810.

[26] 李双喜，曾昭发，李静，等. 探地雷达极化探测时域有限差分法模拟效果分析[J]. 工程地球物理学报，2014， 11(4)：513-521.

Li Shuangxi, Zeng Zhaofa, Li Jing, et al. The Simulation of GPR Polarization Exploration Response with Finite-Difference Time-Domain Method[J]. Chinese Journal of Engineering Geophysics, 2014, 11(4): 513-521.

[27] 陈冬. 地震多属性分析及其在储层预测中的应用研究[D].北京：中国地质大学，2008.

Chen Dong. Research on Seismic Multi-Attributes Analysis and Its Application in Reservoir Predict[D]. Beijing: China University of Geosciences, 2008.

[28] 张军华，朱焕，高荣涛，等. 地震复合属性：地震属性提取与解释新方法[J].新疆石油地质，2007，28(4)：494-496，503.

Zhang Junhua, Zhu Huan, Gao Rongtao, et al. Compound Attribute as New Method for Pickup and Interpretation of Seismic Attributes[J]. Xinjiang Petroleum Geology, 2007, 28(4):494-496,503.

[29] 乐友喜，江凡，问雪，等. 用于地震反射界面识别的瞬时相位复合属性[J]. 物探化探计算技术，2012， 34(5)：505-509.

Yue Youxi, Jiang Fan, Wen Xue, et al. Composite Attribute Derived from Instantaneous Phase for the Recognition of Reflecting Interface[J]. Computing Techniques for Geophysical & Geochemical Exploration, 2012, 34(5):505-509.

[30]Tsoflias G P, Van Gestel J P, Stoffa P L, et al. Vertical Fracture Detection by Exploiting the Polarization Properties of Ground-Penetrating Radar Signals[J]. Geophysics, 2004, 69(3): 803-810.