1.College of Earth Sciences, Guilin University of Technology, Guilin 541004, China 2.Department of Earth and Space Sciences, Southern University of Science and Technology, Shenzhen 518055, China
Fund Project:Project supported by the Natural Science Foundation of Guangxi Province, China (Grant Nos. 2018GXNSFAA281028, 2016GXNSFBA380195), the National Natural Science Foundation of China (Grant No. 41604097), the Research Foundation from Guilin University of Technology, China (Grant No. 002401003503), and the Nonferrous Metal Hidden Heposit Exploration and Material Development Collaborative Innovation Center Project of Guangxi Province, China (Grant No. GXYSXTZX2017-Ⅱ-5).
Received Date:20 November 2018
Accepted Date:08 March 2019
Available Online:01 May 2019
Published Online:05 May 2019
Abstract:Nonlinear and non-stationary ground penetrating radar (GPR) data for geophysics exploration are often mixed with various complex noise sources, such as random and coherent noise. Those complex noise sources are introduced by the acquisition system and other sources of measurement uncertainty. The data sets which are contaminated by the above noises have a serious influence on the accurate extraction of weak reflected signals and the effective identification of diffracted wave hyperbolic coaxial characteristics. The ignorance of the influence of noise will cause great errors in the interpretation of GPR data and the subsequent migration imaging with full waveform. The curvelet transform not only is related to position and frequency as compared with the wavelet transform, but also is controlled by the translation angle. With such a unique advantage, curvelets are used for ground roll whitening, coherent noise denoising and separation of overlapping events. The traditional curvelet transform denoising with a hard threshold function needs to give a reasonable threshold function control coefficient according to the noise level of GPR data. As a result, an appropriate choice of a hard threshold is a basic requirement, and this presents a challenging task in curvelet denoising. To overcome these shortcomings, an self-adaptive threshold function for GPR data denoising with curve transform is proposed in this paper. For detailing the reasonable control coefficient of the threshold function, the complex block threshold function algorithm is used to analyze the distribution of peak-signal-to-noise ratio value of the noisy synthetic GPR data contaminated with random and coherent noise by using the traditional threshold function curvelet transform. Based on the high order statistical theory, the accuracy distribution of the curvelet coefficient for useful signals in scale and rotation direction are correlatively stacked and determined by using the statistical self-adaptive function. And then the threshold range of noise removal components is given by the statistical self-adaptive function. The effectiveness of the proposed denoising algorithm for the noisy synthetic contaminated with different types of noise (i.e., Gaussian random and coherent), and field GPR data is demonstrated by comparing the denoising results via curvelet transform with those from traditional thresholding function. The presented denoising results by the statistical self-adaptive function is helpful and instructive for the accurate inference and interpretation of complex GPR data. Keywords:ground penetrating radar denoising/ self-adaptive thresholding function/ statistic/ curvelet transform
建立理论模型并采用二维有限差分程序计算合成模拟数据. 矩形目标体规模为2 m × 1 m, X轴分布–1—1 m, Z轴分布为2—3 m; 电阻率为10 Ω·m、介电常数为30 F/m; 背景介质电阻率为1000 Ω·m、介电常数为3 F/m; 探地雷达观测系统频率为100 MHz, 点距为0.1 m, X轴方向设置测点数为101, 时间域采集点数为241. 文中峰值信噪比(PSNR)定义为$PSNR = 20\;{\log _{10}}\left(\dfrac{{{\rm max} \left( s \right)}}{{MSE}}\right)$, 其中MSE为数据均方差, 表示原始信号和噪声(或去噪后)信号的近似程度, max(s)为原始信号的峰值. 图1(a)所示为无噪声模拟合成探地雷达原始时间域数据, 可见在空间域设定矩形目标模型边界引起的探地雷达反射波信号分布于不同时间范围. 图1(b)和图1(c)分别为加入随机噪声(PSNR = 17 dB)、相关噪声后的探地雷达数据(PSNR = 15.51 dB). 椒盐式无规律的随机噪声分布于整个探地雷达剖面, 目标弱反射信号完全被噪声淹没; 相关噪声根据有效反射信号范围比例分布, 可见部分假反射信号. 采用(6)式阈值函数曲波变换去噪数据的PSNR值随阈值控制系数$\delta $的变化如图1(d)所示. 随着阈值控制系数$\delta $取值由小变大, 去噪数据PSNR值趋于极值后逐渐降至原噪声数据PSNR值水平, 其中含随机噪声数据去噪结果PSNR值在阈值控制系数$\delta $ = 0.15 (与原始数据噪声水平一致)附近取得PSNR极值(27.33 dB), 含相关噪声数据去噪结果对应在$\delta $ = 0.35附近取得PSNR极值(23.27 dB). 上述结果分析表明, 采用阈值函数曲波变换去噪效果的好坏取决于阈值控制系数(原始数据噪声水平)取值是否合理. 同时, 含随机噪声与相关噪声数据采用阈值函数曲波变换去噪的合理阈值控制系数规律不同, 前者阈值控制系数取噪声水平值($\delta $ = 0.10—0.20), 而后者取略高于阈值控制系数值($\delta $ = 0.3—0.4). 图 1 含噪声模拟合成探地雷达数据的曲波变换阈值去噪结果 (a)原始合成数据; (b)含随机噪声数据, PSNR = 17 dB; (c)含相关噪声数据, PSNR = 15.51 dB; (d)去噪数据PSNR值随阈值控制系数δ的变化 Figure1. Denoised results for the synthetic ground penetrating radar (GPR) data with random and coherent noise using traditional curvelet transform: (a) Original GPR data; (b) data with random noise, PSNR = 17 dB; (c) data with coherent noise, PSNR = 15.51 dB; (d) the δ value vs. PSNR value curves
针对含噪模拟数据阈值函数曲波变换去噪问题, 预先获取的准确阈值控制系数范围可有效用于含噪数据处理. 但实测数据常常被不同类型、不同强度的噪声源污染, 合理的阈值控制系数范围难以确定. 为了有效使用阈值函数曲波变换去噪算法, 需要通过特定方式人为估计阈值控制系数范围, 如L2标准方差算法. 图2所示为含噪声数据(图1(b)和图1(c))采用L2标准方差估计方法确定阈值控制系数为0.08时的曲波变换去噪结果. 相比原始含噪数据(图1(b)和图1(c)), 阈值函数曲波变换去噪清除了部分噪声成分信号, 但矩形目标模型引起的有效反射波信号仅依稀可辩, 无法有效用于全波形反演成像计算及后续推断解译过程(见图2). 图 2δ = 0.08时含噪声数据(图1(b)和图1(c))曲波变换去噪结果 (a)随机噪声数据去噪结果, PSNR = 20.23 dB, 提高3.23 dB; (b)相关噪声数据去噪结果, PSNR = 16.75 dB, 提高1.24 dB Figure2. Denoised results for the synthetic GPR data (Fig. 1. (b) and Fig. 1. (c)) with δ = 0.08 using traditional curvelet transform: (a) Result for random noise, PSNR = 20.23 dB; (b) result for coherent noise, PSNR = 16.75 dB.
23.2.统计量自适应阈值曲波变换去噪 -->
3.2.统计量自适应阈值曲波变换去噪
为解决传统曲波变换去噪需要估计合理阈值函数范围的问题, 采用统计量自适应阈值曲波变换去噪方法对图1所示含随机噪声、相关噪声数据进行处理. 去噪结果如图3所示, 基于统计量自适应阈值曲波变换去噪方法不仅可有效衰减随机噪声(PSNR值提高8.3 dB)、相关噪声(PSNR值提高6.41 dB), 同时较好地还原剖面中有效反射波信号时空域分布特征. 在噪声处理过程中, 本文提出的去噪算法无需估计阈值函数范围, 而是采用(8)式计算统计量自适应阈值权重, 达到了保真去噪的目的. 其中, 随机噪声数据去噪效果优于相关噪声数据去噪效果, 究其缘由, 低信噪比相关噪声数据(PSNR = 15.51 dB)导致噪声信号相关性统计量计算误差, 残留部分相关性较强的噪声信号. 总体上, 统计量自适应阈值曲波变换去噪结果可有效地用于探地雷达数据全波形偏移成像处理及后续解释推断过程. 图 3 含噪声数据(图1(b)和图1(c))的统计量自适应阈值曲波变换去噪结果 (a)随机噪声数据去噪结果, PSNR = 25.3 dB, 提高8.3 dB; (b)相关噪声数据去噪结果, PSNR = 21.92 dB, 提高6.41 dB Figure3. Denoised results for the synthetic GPR data (Fig. 1. (b) and Fig. 1. (c)) using curvelet transform with statistical self-adaption: (a) Result for random noise, PSNR = 25.3 dB; (b) result for coherent noise, PSNR = 21.92 dB.
进一步, 设计双矩形目标体规模为1 m × 1 m, 水平位置分别为–2至–3 m及2至3 m, 埋深范围分别为2至3 m及3至4 m; 其他物性参数及探地雷达观测系统与前述单个矩形目标模型一致. 数值模拟的双矩形目标体探地雷达时间剖面如图4(a)所示, 加入随机噪声及相关噪声数据如图4(b)和图4(c)所示. 含随机噪声数据(PSNR = 16.04 dB)及相关噪声数据(PSNR = 15.7 dB)完全淹没双矩形目标体引起的有效反射波信号, 无法用于双矩形目标体时空域分布推断解释. 图 4 双矩形目标模型含噪声数据统计量自适应阈值曲波变换去噪结果 (a)原始合成数据; (b)含随机噪声数据, PSNR = 16.04 dB; (c)含相关噪声数据, PSNR = 15.7 dB; (d)随机噪声数据去噪结果, PSNR = 23.97 dB; (e)相关噪声数据去噪结果, PSNR = 21.05 dB Figure4. Denoised results for the synthetic GPR data of complex model with random and coherent noise using curvelet transform with statistical self-adaption: (a) Original GPR data; (b) data with random noise, PSNR = 16.04 dB; (c) data with coherent noise, PSNR = 15.7 dB; (d) result for random noise, PSNR = 23.97 dB; (e) result for coherent noise, PSNR = 21.05 dB.
图4(d)和图4(e)为采用本文提出的去噪算法处理结果. 由去噪结果可见, 含随机噪声时间剖面内双矩形目标体引起复杂有效反射波信号被完全恢复, PSNR值较原始含噪数据提高7.93 dB; 含相关噪声数据处理结果基本重建目标体引起有效反射波信号时空分布, 但仍残余部分相关噪声成分, PSNR值较原始含噪数据提高6.65 dB. 抽取图4所示去噪结果第50道数据, 绘制原始无噪声数据、含噪数据及去噪数据信号局部细节对比曲线, 如图5所示. 可见, 本文提出的去噪算法有效清除了全局椒盐式随机噪声分布(图5(a))、可有效分辨相关噪声环境下微弱有效反射波信号(图5(b)), 去噪数据曲线与原始无噪声数据曲线吻合较好, 验证了本文提出的去噪算法有效性及可行性. 图 5 双矩形目标模型含噪声数据第50道统计量自适应阈值曲波变换去噪结果 (a)随机噪声数据去噪结果; (b)相关噪声数据去噪结果 Figure5. Denoised results for the synthetic GPR data (Fig. 4.) in the 50th receiver using curvelet transform with statistical self-adaption: (a) Result for random noise; (b) result for coherent noise.
数值模拟实验结果表明, 传统曲波变换去噪效果依赖于阈值控制系数选取范围是否符合含噪数据噪声水平范围, 而不同信噪比含噪数据噪声水平各异, 因此采用传统阈值函数的曲波去噪算法抗噪性能与估计阈值控制系数范围准确度具有直接关系. 对于含有弱信号成分探地雷达数据的去噪处理, 普遍应用于传统曲波变换去噪算法理论与实践研究的L2标准差估计的噪声水平导致伪异常信号问题. 本文提出的统计量自适应阈值曲波变换去噪算法可有效确定含噪数据噪声水平范围, 在单一以及多个异常信号叠加的探地雷达数据处理方面显示了较好的应用效果. 需要指出的是, 本文考虑数据噪声水平为15%, 远远高于地球物理勘探数据采集质量控制误差水平5%, 数值算例结果可见, 含随机噪声、相关噪声数据处理结果可有效地用于后续成像及解译; 而高于15%含噪数据通常视为无效数据, 其相应去噪处理分析研究不在本文讨论范围. 4.实测算例图6(a)所示为某地探地雷达实测某测线原始时间剖面图, 由于受场地环境影响, 椒盐式随机噪声遍布整个剖面数据; 部分采集道数据受局部不均匀体影响, 出现相邻观测道数据幅值畸变. 原始数据剖面浅部(100—250时间采样点)介质的反射波能量非常微弱, 依稀可见几处绕射波同相轴, 但受噪声淹没, 双曲线特征不明显; 在250—300时间采样范围内出现较强幅值的似平行反射波组能量, 同相轴断断续续并不连续; 强幅值反射波组下部(300—400时间采样点)出现较弱的零散反射波组, 难以识别强幅值反射波组下边界. 图 6 实测探地雷达时间剖面传统曲波变换去噪和本文去噪算法处理结果 (a)实测时间剖面; (b)采用L2标准方差估计阈值曲波去噪结果; (c)统计量自适应阈值曲波去噪结果 Figure6. Denoised results for field GPR data using curvelet transform with L2 standard deviation and statistical self-adaption respectively: (a) Original field GPR data; (b) result using curvelet transform with L2 standard deviation; (c) result using curvelet transform with statistical self-adaption.
利用上述传统曲波变换去噪处理方法对图6(a)进行数据处理, 通过L2标准方差估计阈值范围, 去噪获得的探地雷达时间剖面如图6(b)所示. 由去噪结果可见, 部分地下反射波组能量信号得以加强, 清除了部分随机噪声信号. 但原始数据内有效绕射波双曲线特征、错动反射波组依然无法有效识别, 去噪数据信噪比提高并不显著. 究其原因, 实测数据中随机噪声、相关噪声及其他类型噪声源掺杂, 同时, 原始数据包含幅值强度各异的反射波组, L2标准方差估计阈值范围无法有效涵盖上述噪声特点的数据去噪范围. 人工试错法确定最佳阈值范围或能有效清除噪声, 但最佳去噪效果评价缺乏科学依据. 因此, 基于传统阈值函数曲波变换去噪数据处理后的探地雷达时间剖面进行分析解释, 无法保证资料解译的准确度. 采用本文提出的统计量自适应阈值曲波变换去噪算法对图6(a)原始数据进行处理, 去噪结果如图6(c)所示. 由于随机噪声源在曲波变换系数尺度、方向上并不具备统计相关特性; 相关噪声源变换系数在尺度、方向上虽具备一定相关性, 但相对有效反射波组信号相关性仍然较小, 因而统计量自适应阈值可有效分辨反射波组能量, 去除噪声信号, 达到保真去噪的目的. 由图6(c)可见, 地下反射波能量得到显著增强, 特别是浅部多个绕射波双曲线同相轴特征的信噪比得到有效提高; 中深度多层介质的反射组同相轴连续、分界面明显; 且分布于强幅值反射波组下部的弱幅值反射波组清晰可辩. 经过统计量自适应阈值曲波变换数据处理后的探地雷达时间剖面进行分析解释, 有助于资料解译的准确度. 抽取图6所示实测探地雷达数据第200道数据, 绘制传统曲波变换去噪结果与本文提出的去噪结果局部细节对比曲线, 如图7所示. 由曲线对比图可见, 在0—100时间采样范围内地面回波、弱幅值绕射波(100—150时间采样)、强幅值绕射波(150—230时间采样)及强幅值反射波组(230—350时间采样)受到随机噪声及相关噪声污染, 在第350时间采样范围以外的弱反射波信号受到较强幅值的上述噪声源淹没. 相比传统曲波变换去噪结果, 本文提出的去噪算法可有效恢复弱幅值绕射波、强幅值绕射波及强幅值反射波组的同相轴特征, 同时可清晰分辨弱幅值反射波信号, 验证了本文提出的去噪算法在实测探地雷达时间剖面数据处理应用中的可行性及有效性. 图 7 实测探地雷达时间剖面第200道传统曲波变换去噪和本文去噪算法处理结果 Figure7. Denoised results for field GPR data (Fig. 6) in the 200th receiver using curvelet transform with L2 standard deviation and statistical self-adaption respectively.