程浩¹, 王德利², 王恩德¹, 付建飞¹
1. 东北大学 深部金属矿山开采教育部重点实验室, 辽宁 沈阳 110819;
2. 吉林大学 地球探测科学与技术学院, 吉林 长春 130026
收稿日期：2018-05-15
基金项目：国家自然科学基金资助项目(41804103)；国家重点研发计划项目(2016YFC0801603，2017YFC1503101)；中央高校基本科研业务费专项资金资助项目(N160103001)。
作者简介：程浩(1988-), 男, 辽宁沈阳人, 东北大学师资博士后;
王恩德(1957-), 男，辽宁盖州人，东北大学教授，博士生导师。

摘要：提出GPR数据Curvelet域随机噪声压制与高频补偿同步处理方法.首先, 将GPR数据变换到Curvelet域, 结合其多角度、多尺度的稀疏性, 给出随尺度和角度变化的自适应阈值函数进行随机噪声的压制; 其次, 根据电磁波在完全弹性介质中的传播规律, 结合Curvelet的多尺度多角度特性, 求取时变补偿因子, 倒数加权对应的尺度、角度, 补偿高频损失; 最后, 进行Curvelet反变换, 获得随机噪声压制与高频补偿以后的GPR数据.该方法属于完全数据驱动, 克服了传统方法人为因素的影响.
关键词：GPR数据Curvelet变换随机噪声高频衰减去噪方法吸收补偿
Noise Suppression and High-Frequency Compensation of GRP Data in Curvelet Domain
CHENG Hao¹, WANG De-li², WANG En-de¹, FU Jian-fei¹
1. Key Laboratory of Ministry of Education on Safe Mining of Deep Metal Mines, Northeastern University, Shenyang 110819, China;
2. College of Geo-Exploration Science and Technology, Jilin University, Changchun 130026, China
Corresponding author: CHENG Hao, E-mail: chenghao@mail.neu.edu.cn
Abstract: This paper proposed a synchronous processing method using random noise suppression and high-frequency compensation for GPR data in curvelet domain. First, the GPR data were transformed into curvelet domain in order to use its multi-angle and multi-scale sparsity. The random noise thus could be suppressed by the adaptive threshold function changed with angle and scale. Second, according to the propagation laws of electromagnetic waves in elastic medium and the curvelet multi-scale feature, the time-varying compensation factor could be calculated so as to inverse weighting the corresponding scale and angle for compensation of the high-frequency absorption. Finally, inverse curvelet transformation was conducted thus the GPR data after random noise suppression and high-frequency compensation could be obtained. This method is totally a data-driven method thus can overcome the influence from artificial factors in traditional method.
Key words: GPR datacurvelet transformrandom noisehigh-frequency attenuationde-noising methodabsorption-compensation
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的分辨率降低, 且传播深度越大, 分辨率越低.地下介质的反射系数序列和反射时刻可以表示为r_d(k=0, 1, …, D), T_d.经过地下介质滤波作用的振幅响应可以表示为

(1)

式中:

(2)

τ表示双程旅行时, Q(τ)表示τ时的地层Q值, Q_eq(T_d)表示反射时刻T_d的等效Q值, A₀(f, 0)表示初始时刻的地层振幅.如果仅考虑振幅的衰减作用, 则频率域电磁波记录可表示为

(3)

式中: 为电磁波记录的傅里叶变换; 为子波ω(t)的傅里叶变换.则在T_d时刻的振幅谱为

(4)

则当相邻反射波互不干涉, 相对于初始时刻, T_d时刻的振幅比率为

(5)

为了消除地层反射系数的影响, 定义衰减比率η为某一频率与特定频率的振幅比率之比:

(6)

式(6)表明, η(f, T_d)仅与时间和频率有关.在给定频率下, η(f, T_d)仅与传播时间有关.因此, 用η(f, T_d)倒数加权对应频段的电磁波信号, 就能消除介质损耗的影响.
2 Curvelet域去噪与高频补偿同步处理2.1 Curvelet变换任意函数z∈L²(R²), 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变换将电磁波信号分解为频带很窄的不同频带不同方向的c_{j, l}(f, T_k), 进而求取振幅比率：

(11)

根据有效信号的分布范围, 也可以选取特定的角度, 再用振幅比率求取衰减比率：

(12)

利用求得的衰减比率倒数加权到对应尺度和角度的Curvelet系数上, 再进行Curvelet反变换, 重构GPR的电磁波信号, 就达到了补偿高频成分的目的.
图 1为Curvelet域GPR数据去噪与高频补偿同步处理的流程图.首先, 输入GPR数据, 利用Curvelet变换对GPR数据进行多尺度多角度的剖分.GPR数据的有效信号在Curvelet域对应于不同尺度不同方向上较大的Curvelet系数.利用自适应阈值函数计算不同尺度不同方向上的阈值, 进而利用式(9)去除小于阈值的Curvelet系数, 达到去除随机噪声的目的.其次, 利用去除随机噪声以后的Curvelet系数进行振幅比率和衰减比率的计算, 再倒数加权到对应尺度和方向上, 补偿高频成分.最后, 通过Curvelet反变换重构GPR数据, 实现一次Curvelet变换同时进行随机噪声压制和高频补偿的目的.
图 1(Fig. 1)

图 1 流程图Fig.1 Flow diagram

3 实际数据应用为了验证该方法的有效性与准确性, 利用实际GPR数据进行测试.实际GPR数据为一路基结构探测数据, GPR天线发射频率为100 MHz, 采样频率为1 500 MHz, 道间距为0.098 m, 实际探测结果如图 2a所示.可以发现, GPR数据剖面80~150样点范围内, 地层结构较为明显, 但剖面150~600道之间同相轴的某些部分连续性较差, 可能受到随机噪声或者存在路基破裂带的影响.同时, GPR数据150样点以下, 由于介质的损耗作用, 电磁波高频严重衰减, 底层反射波能量减弱, 从剖面上基本得不到任何有效的结构信息, 对路基底层结构的判断造成较为严重的影响.同时, GPR数据剖面中存在的随机噪声影响GPR数据剖面的成像效果.
图 2(Fig. 2)

图 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(Fig. 3)

图 3 GPR数据时频谱Fig.3 Time-frequency spectrum of GPR data (a)—原始数据；(b)—补偿与去噪以后数据.

图 4(Fig. 4)

图 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