2019, Vol. 40
2. 吉林大学 地球探测科学与技术学院, 吉林 长春 130026
2. College of Geo-Exploration Science and Technology, Jilin University, Changchun 130026, China
GPR是一种借助于高频电磁波探测地下介质分布与性质的地球物理方法[1-3].在近地表地球物理中, 以其操作简单、分辨率高等特点得到了广泛应用, 但在实际探测中, 电磁波在有耗介质中传播, 使得高频成分随着传播深度的增加而严重衰减, 降低了对底层介质的探测能力.同时, 受复杂环境和仪器自身的影响, GPR数据中将包含大量随机噪声, 影响数据解释的精度.
现有的GPR数据高频衰减补偿方法主要有反褶积法、反Q滤波法和谱白化法.反褶积方法需要对GPR数据进行子波估计, 反Q滤波方法需要对GPR数据进行Q值估计, 谱白化方法需要对GPR数据进行开时窗处理.三种方法受处理人员经验水平的影响较大, 不同的处理人员可能给出不同的处理结果.此外, 以上三种方法对GPR数据的信噪比也有一定要求, 在补偿有效信号的同时也会使噪声信号得到增强.
与地震波特点类似, GPR高频电磁波在完全弹性介质中传播时, 不会产生高频衰减, 不同频段深部和浅部的振幅谱应该相同, 即在频率域, 电磁波不同频段的信号应该具有相似的时间能量分布关系, 深浅层应该具有相同的能量比[4-5].在GPR实际探测中, 介质的损耗改变了各频段的时间能量分布关系.由于给定频率时, 衰减比率只与传播时间有关, 所以, 求取一个时变补偿因子, 使各频段的深浅层能量比相同, 就实现了补偿GPR信号的作用.
Curvelet变换[6]是一种非自适应多尺度多方向的几何分析方法, 具有极佳的稀疏特性, 能够以极少的系数进行信号的重构.在信号的去噪、插值、稀疏反演等数据处理方面得到了广泛应用[7-9].
本文推导了Curvelet域去噪与高频补偿同步处理的理论过程.在Curvelet域, 首先, 结合自适应阈值函数进行随机噪声的压制; 其次, 根据Curvelet变换的多尺度特性, 求取时变补偿因子, 进行高频补偿; 最后, 进行Curvelet反变换, 重构GPR数据, 实现Curvelet域去噪与高频补偿同步处理的目的.利用GPR实际数据进行测试, 验证该方法的有效性与准确性.
1 数据驱动吸收补偿原理GPR的电磁波在地下介质中传播时, 由于介质的损耗, 使得高频成分随着传播深度的增加而急剧衰减, 使GPR的分辨率降低, 且传播深度越大, 分辨率越低.地下介质的反射系数序列和反射时刻可以表示为rd(k=0, 1, …, D), Td.经过地下介质滤波作用的振幅响应可以表示为
|
(1) |
式中:
|
(2) |
τ表示双程旅行时, Q(τ)表示τ时的地层Q值, Qeq(Td)表示反射时刻Td的等效Q值, A0(f, 0)表示初始时刻的地层振幅.如果仅考虑振幅的衰减作用, 则频率域电磁波记录可表示为
|
(3) |
式中: 
|
(4) |
则当相邻反射波互不干涉, 相对于初始时刻, Td时刻的振幅比率为
|
(5) |
为了消除地层反射系数的影响, 定义衰减比率η为某一频率与特定频率的振幅比率之比:
|
(6) |
式(6)表明, η(f, Td)仅与时间和频率有关.在给定频率下, η(f, Td)仅与传播时间有关.因此, 用η(f, Td)倒数加权对应频段的电磁波信号, 就能消除介质损耗的影响.
2 Curvelet域去噪与高频补偿同步处理 2.1 Curvelet变换任意函数z∈L2(R2), Curvelet变换定义为
|
(7) |
式中:C为Curvelet正变换; c为Curvelet系数; j为尺度参数; l为剪切参数; k为平移参数; φj, l, k为Curvelet基函数.Curvelet反变换可以表示为
|
(8) |
式中C-1表示Curvelet反变换.
2.2 Curvelet域自适应阈值去噪GPR数据在Curvelet域能够得到最佳的稀疏表示.有效信号通常具有一定方向, 当与Curvelet基方向大致相同时, 有效信号仅集中在较少且较大的Curvelet系数上, 而随机噪声不具有方向性, 所以对应的Curvelet系数值较小.因此, 通过阈值函数去除较小的系数, 达到压制噪声的目的, 表示为
|
(9) |
式中Thr为阈值.根据Curvelet变换多尺度、多角度的特点, 以及有效信号和随机噪声在Curvelet域的分布特点, 本文给出随尺度和角度变化的自适应阈值算法:
|
(10) |
式中:N为GPR数据总的采样点数; E为L1范数值; σ为噪声标准差; λ为任意小正数.
2.3 Curvelet域高频补偿在实际GPR探测中, 电磁波的低频成分衰减较小, 此处忽略不计.利用Curvelet变换将电磁波信号分解为频带很窄的不同频带不同方向的cj, l(f, Tk), 进而求取振幅比率：
|
(11) |
根据有效信号的分布范围, 也可以选取特定的角度, 再用振幅比率求取衰减比率：
|
(12) |
利用求得的衰减比率倒数加权到对应尺度和角度的Curvelet系数上, 再进行Curvelet反变换, 重构GPR的电磁波信号, 就达到了补偿高频成分的目的.
图 1为Curvelet域GPR数据去噪与高频补偿同步处理的流程图.首先, 输入GPR数据, 利用Curvelet变换对GPR数据进行多尺度多角度的剖分.GPR数据的有效信号在Curvelet域对应于不同尺度不同方向上较大的Curvelet系数.利用自适应阈值函数计算不同尺度不同方向上的阈值, 进而利用式(9)去除小于阈值的Curvelet系数, 达到去除随机噪声的目的.其次, 利用去除随机噪声以后的Curvelet系数进行振幅比率和衰减比率的计算, 再倒数加权到对应尺度和方向上, 补偿高频成分.最后, 通过Curvelet反变换重构GPR数据, 实现一次Curvelet变换同时进行随机噪声压制和高频补偿的目的.
|
图 1 流程图 Fig.1 Flow diagram |
为了验证该方法的有效性与准确性, 利用实际GPR数据进行测试.实际GPR数据为一路基结构探测数据, GPR天线发射频率为100 MHz, 采样频率为1 500 MHz, 道间距为0.098 m, 实际探测结果如图 2a所示.可以发现, GPR数据剖面80~150样点范围内, 地层结构较为明显, 但剖面150~600道之间同相轴的某些部分连续性较差, 可能受到随机噪声或者存在路基破裂带的影响.同时, GPR数据150样点以下, 由于介质的损耗作用, 电磁波高频严重衰减, 底层反射波能量减弱, 从剖面上基本得不到任何有效的结构信息, 对路基底层结构的判断造成较为严重的影响.同时, GPR数据剖面中存在的随机噪声影响GPR数据剖面的成像效果.
|
图 2 GPR数据剖面 Fig.2 GPR-data profile (a)—原始数据；(b)—补偿与去噪以后数据. |
利用本文所述的Curvelet域去噪与高频补偿方法, 对上述实际GPR数据进行处理.由于实际GPR数据为一路基数据, 其地下结构理论上应该呈水平层状结构分布, 所以在进行Curvelet变换高频补偿时, 可以把补偿的角度设置小一些, 本文所选补偿角度为0°~30°, 其他角度不进行补偿, 这样就能在补偿有效信号的同时, 不对干扰方向进行补偿, 也不降低补偿以后数据的信噪比.Curvelet变换时, 利用的尺度参数为5, 角度参数为32, 首先进行随机噪声的压制, 再进行衰减比率的计算, 进而倒数加权对应频率和方向的电磁波信号, 最后将去噪与高频补偿以后的Curvelet系数进行Curvelet反变换, 获得实际GPR数据去噪与高频补偿的结果, 如图 2b所示.
从GPR实际数据可以发现, 随机噪声被有效压制, 信噪比提高, 底层有效信息得到有效的补偿, 由于介质损耗而衰减的电磁波高频信息得到较好的恢复, 底层的构造信息得到了较好的反映, 同相轴能量增强.与图 2a相比, 经过Curvelet域去噪与高频补偿以后的GPR数据信噪比提高, 底层结构信息变得更加丰富, 为路基底层结构的分析与判断提供了良好的依据.
为了能更好体现Curvelet域高频补偿方法对GPR探测方法高频损失补偿的有效性, 选取GPR实际原始数据和补偿以后的第750道进行时频分析, 绘制时频谱, 观察补偿前后GPR实际数据的时频特征.本文所用时频分析方法为传统的短时傅里叶变换, 结果如图 3所示.由图 3a可知, 原始GPR实际数据的高频成分随着传播时间的增加, 被严重衰减, 频带随之变窄, 得到的底层有效信息相对较少.相比之下, 经过Curvelet域高频补偿以后的GPR实际数据的时频谱(图 3b)衰减的高频信息得到了有效补偿, 底层频率信息变得丰富, 频带得到了明显拓宽, 地下底层的有效信息得到有效表达.利用补偿以后的GPR实际数据进行简单的结构分析和裂隙勾画, 如图 4所示.根据同相轴的趋势和不连续性, 可以较为清晰地判断出裂缝的走向, 在图 4中给出较为明显的7条断裂趋势线, 为路基的病害程度判断提供了可靠的依据.
|
图 3 GPR数据时频谱 Fig.3 Time-frequency spectrum of GPR data (a)—原始数据；(b)—补偿与去噪以后数据. |
|
图 4 补偿与去噪以后GPR数据剖面结构分析图 Fig.4 Structure analysis of compensated and de-noised GPR-data profile |
通过上述GPR实际数据的验证可知, 基于Curvelet变换的GPR数据去噪与高频补偿方法, 可以有效提高GPR实际数据的信噪比和补偿高频损失, 去噪与补偿以后的GPR数据底层结构信息变得丰富, 图像变得更加清晰, 对比补偿前后GPR数据的时频谱可以发现, 底层部分的频带信息得到了较好的拓宽, 使得底部的结构信息清晰凸显出来, 证明了Curvelet域GPR数据去噪与高频补偿方法的有效性和准确性.
4 结论1) 根据GPR数据在Curvelet域的分布特点, 以及Curvelet变换多尺度多角度的信号剖分特性, 给出随尺度和角度变化的自适应阈值, 能够更好地压制随机噪声, 保留GPR数据的有效信号.
2) 利用Curvelet变换对GPR数据进行多尺度、多角度的分频处理, 根据实际需求可以选定需要补偿的角度进行补偿, 其他方向则不对其进行补偿.该方法属于完全数据驱动, 避免了人为因素引入的误差.
3) 在Curvelet域, 同步进行随机噪声的压制与高频补偿, 提高了数据处理的效率.同时, 随机噪声的压制, 信噪比的提高, 也是有效进行高频补偿的前提.通过实际路基数据的测试, 验证了该方法的有效性与适用性.
| [1] |
Wang X N, Liu S X.
Noise suppressing and direct wave arrivals removal in GPR data based on Shearlet transform[J]. Signal Processing, 2017, 132: 227–242.
DOI:10.1016/j.sigpro.2016.05.007 |
| [2] |
佘黎煌, 王培人, 张石.
基于块目标的频率步进连续波探地雷达压缩感知重建算法[J]. 东北大学学报(自然科学版), 2018, 39(3): 316–319.
( She Li-huang, Wang Pei-ren, Zhang Shi. Reconstruction algorithm of compressed sensing for stepped-frequency continuous wave ground penetrating radar based on block objects[J]. Journal of Northeastern University(Natural Science), 2018, 39(3): 316–319. ) |
| [3] |
Liu C M, Wang D L, Sun J L, et al.
Crossline-direction reconstruction of multi-component seismic data with shearlet sparsity constraint[J]. Journal of Geophysics and Engineering, 2018, 15: 1929–1942.
DOI:10.1088/1742-2140/aac097 |
| [4] |
Liu C M, Wang D L, Sun J L, et al.
Stratigraphic absorption compensation based on multiscale shearlet transform[J]. Acta Geophysica, 2018, 66(4): 575–584.
DOI:10.1007/s11600-018-0156-8 |
| [5] |
Wang D L, Sun J L, Meng D J, et al.
Absorption compensation based on curvelet transform[J]. Journal of Seismic Exploration, 2013, 22: 19–23.
|
| [6] |
Candès E, Demanet L, Donoho D, et al.
Fast discrete curvelet transforms[J]. Multiscale Modeling & Simulation, 2006, 5(3): 861–899.
|
| [7] |
程浩, 陈刚, 王恩德, 等.
基于Shearlet变换的自适应阈值地震数据去噪方法[J]. 石油学报, 2018, 39(1): 82–91.
( Cheng Hao, Chen Gang, Wang En-de, et al. Seismic data de-noising method of adaptive threshold based on shearlet transform[J]. Acta Petrolei Sinica, 2018, 39(1): 82–91. ) |
| [8] |
Cheng H, Wang D L, Feng F, et al.
Estimating primaries from passive seismic data[J]. Exploration Geophysics, 2015, 46: 184–191.
DOI:10.1071/EG14079 |
| [9] |
Wang T X, Wang D L, Sun J, et al.
Closed-loop SRME based on 3D L1-norm sparse inversion[J]. Acta Geophysics, 2017, 65: 1145–1152.
DOI:10.1007/s11600-017-0098-6 |