CN109581516B - 曲波域统计量自适应阈值探地雷达数据去噪方法及*** - Google Patents

Abstract

本发明属于数据去噪技术领域，公开了一种曲波域统计量自适应阈值探地雷达数据去噪方法及***，方法包括：引入块状复数域阈值函数算法，分析传统阈值函数曲波变换去噪效果随阈值函数控制系数变化规律，用于后续曲波域统计量自适应阈值对比；利用高阶统计量理论，对曲波变换系数在尺度、方向上进行相关性叠加，通过相关性统计量自适应确定有效信号在曲波变换系数分布尺度、旋转方向；确定出清除噪声成分阈值范围，构建统计量自适应阈值函数曲波变换去噪算法。本发明对包含随机噪声、相关噪声合成探地雷达数据及实测探地雷达数据处理结果与现有技术对比，成果对复杂探地雷达数据精确推断解译具有指导意义。

Description

曲波域统计量自适应阈值探地雷达数据去噪方法及***

技术领域

本发明属于数据去噪技术领域，尤其涉及一种曲波域统计量自适应阈值探地雷达数据去噪方法及***。

背景技术

目前，业内常用的现有技术是这样的：

探地雷达采用天线发射不同频率电磁波，利用地下介质电磁性质差异，根据接收回波的幅值、波形等运动学及动力学特征推断目标介质空间、物性分布，其广泛应用于地质灾害监测、工程与环境地质勘察、水文地质勘查及军事侦查等领域。然而，随着探地雷达勘探环境越来越复杂、目标探测越来越精细(张先武等,2014；阮超等,2015)，如何在强能量的各种复杂干扰噪声湮没环境中提取微弱的目标信号是一项艰巨的任务(白大为等，2017；Zheng et al.,2017)。

常规探地雷达数据去噪技术多基于简单的最优化或正交的变换算法，其形成的是固定滤波窗口，导致时频域重叠的微弱有效信号失真，具有特定噪声处理局限性，如中值滤波法(王磊等，2009)、S变换去噪(高静怀等,2004；Jeng et al.,2009；张先武等,2013a)及F-K滤波方法等(张先武等,2013b；Wang and Zhang,2015)。由于各种噪声源掺杂，探地雷达数据常为非线性、非平稳信号序列(Wang and Liu,2017)。基于傅里叶变换理论研发的信号处理技术可有效解决平稳信号分析处理问题，非自适应时频窗口局限性使其无法解决探地雷达复杂数据的去噪问题。小波变换用于探地雷达数据去噪的前提是假设有效信号与噪声信号的频谱分离，通过采用具有伸缩和平移特性基函数在特定频段的稀疏表示，设置噪声频段阈值以达到保真去噪的目的。如王超和沈斐敏(2015)利用高频阈值函数小波变换去噪理论开展探地雷达弱信号提取研究，为推断异常位置提供依据。改进阈值的提升小波变换用于混凝土脱空探地雷达数据去噪(许后磊和储冬冬，2010；Ghozzi et al.,2017)，提高了解释推断精度；可协调方向阈值函数的小波变换用于考古学和岩土勘察探地雷达数据去噪及杂波压制(Tzanis,2013、2015)，改善了数据信噪比等。因而，有效的阈值函数可有效提升去噪效果，实现保真去噪的目的。

为突破小波变换采用矩形时频窗口去噪的限制，曲波变换采用小波域伸缩及平移特性，并增加一个方向参量而具有更好的方向识别能力(Candès and Donoho,2004；Candèset al.,2006)。曲波变换自提出便迅速广泛应用于地球物理数据处理领域，特别是针对地震数据噪声压制、多次波分离等领域，研究成果不断涌现(张华和陈小宏,2013；董烈乾等,2016；齐少华等,2016；张华等,2017)。其中，Neelamani等(2008)将曲波变换应用于三维地震数据，通过设定阈值函数有效清除随机噪声、相关噪声，提高了地震资料的信噪比；张金良等(2013)快速曲波变换域的SAR图像去噪算法，使用均值阈值滤波器使得图像的可视化和解译变得精确；Terrasse等(2017)等通过人工挑选探地雷达数据曲波变换域系数，明确需要去除杂乱回波及噪声源的尺度、方位信息，加以硬阈值函数实现去噪处理。可见，曲波变换应用于探地雷达数据噪声压制的前提条件是获取较为准确的噪声阀值参数及其所属曲波系数在尺度、方位上分布范围。

但实测复杂探地雷达数据常包含随机噪声、相关噪声及其他未知噪声类型，所需阀值参数及所属曲波系数在尺度、方位上分布范围亦有所不同。如朱自强等(2014)在角度窗函数下选择变换域中设置方向因子为零、同时估计噪声方差确定阈值函数实现去除表层直达波、噪声源，在强度较高直达波背景下提取弱化了钢筋层和裂隙水层的反射信号。Bao等(2014)认为背景噪声主要能量集中在方向90°区域附近，而随机噪声会相对均匀地分布在整个曲波域，设计二维滤波器估计噪声分布。Tzanis(2017)通过数值模拟人为建立不同裂隙构造对应探地雷达数据曲波变换域尺度、方位分布范围，继而实现特定方向发育裂隙结构探地雷达反射波信号提取。因而，如何有效确定目标体引起的探地雷达数据曲波变换域尺度、方位上分布范围是实现高效去噪、精确偏移成像目的的关键(雷林林等，2015)。当前，曲波变换配置研发了诸多计算阀值的方法，如分段线性滤波方法、L2标准差法、曲波正变换法、对角实数及复数阀值法等，需根据实际数据噪声类型选用，总体上无法满足自适应保真去噪的目的。因而，曲波域自适应阈值的构建是学者关注的研究热点。

综上所述，现有技术存在的问题是：

(1)非线性、非平稳探地雷达数据常掺杂各种复杂噪声源，对精确提取弱反射波信号、识别绕射波双曲线同相轴特征有严重影响，忽略噪声的影响会给探地雷达探测数据全波形偏移成像及后续解译造成很大误差。传统阈值函数曲波变换去噪需要根据探地雷达数据噪声水平给定合理阈值函数控制系数。

(2)曲波变换应用于探地雷达数据噪声压制的前提条件是获取较为准确的噪声阀值参数及其所属曲波系数在尺度、方位上分布范围。但实测复杂探地雷达数据常包含随机噪声、相关噪声及其他未知噪声类型，所需阀值参数及所属曲波系数在尺度、方位上分布范围亦有所不同。

(3)当前，曲波变换配置研发了诸多计算阀值的方法，如分段线性滤波方法、L2标准差法、曲波正变换法、对角实数及复数阀值法等，需根据实际数据噪声类型选用，总体上无法满足自适应保真去噪的目的。

解决上述技术问题的难度和意义：

探地雷达勘探环境越来越复杂、目标探测越来越精细需求，使得在强能量的各种复杂干扰噪声湮没环境中提取微弱的目标信号成为一项极具挑战的课题。如何突破人工试错法选取阈值函数去噪的瓶颈、构建自适应阈值函数的去噪算法是当前研究的重点。通过引入曲波变换与高阶相关统计分析的优势，本发明提出在曲波变换多尺度、多方向系数高阶相关统计量分析理论基础上构建自适应阈值去噪算法，为探地雷达数据自适应阈值保真去噪提供技术手段。

发明内容

针对现有技术存在的问题，本发明提供了一种曲波域统计量自适应阈值探地雷达数据去噪方法及***。本发明结合曲波变换与高阶相关统计分析优势，本发明提出在曲波变换多尺度、多方向分析理论基础上采用高阶相关统计量构建自适应阈值算法，探索探地雷达数据曲波域自适应保真去噪方法技术。

本发明是这样实现的，一种曲波域统计量自适应阈值探地雷达数据去噪方法，所述曲波域统计量自适应阈值探地雷达数据去噪方法包括：

引入块状复数域阈值函数算法，通过选取多种不同块状复数域阈值，分析传统阈值函数曲波变换去噪效果随阈值函数控制系数变化规律，明确可有效用于噪声去除的块状复数域阈值，用于后续曲波域统计量自适应阈值对比；

曲波域统计量即利用高阶相关统计量理论，对曲波变换系数在尺度、方向上进行相关性叠加，通过相关性统计量自适应确定有效信号在曲波变换系数分布尺度、旋转方向，用于后续自适应阈值范围构建；

基于随机噪声、相关噪声成分与有效成分分布的不同尺度、旋转方向，确定出有效信号分布尺度和旋转方向、清除噪声成分阈值范围，构建统计量自适应阈值函数曲波变换去噪算法。

进一步，通过相关性统计量自适应确定有效信号在曲波变换系数分布尺度、旋转方向；包括：

对探地雷达数据做二维的曲波正变换获取不同尺度、不同旋转角度的系数分布；对每一道探地雷达数据的在不同尺度、不同角度下计算高阶相关统计量、根据所述高阶相关统计量分布明确阈值分布权重，实现噪声滤除；最后采用二维的曲波反变换重建去噪后的探地雷达数据；其中，未偏移的三阶相关函数表示为：

归一化高阶相关统计量为：

式中，i为探地雷达数据接收点，t为时间采样点，q为平移因子，取值为1，

为第i道探地雷达数据曲波域第k个尺度、第θ个方向正变换系数；

为第i道探地雷达数据曲波域第k个尺度、第θ个方向二阶自相关函数值；

为第i+1道探地雷达数据曲波域第k个尺度、第θ个方向二阶自相关函数值；

所述通过高阶相关统计量计算自适应确定有效信号在曲波变换系数分布尺度、旋转方向具体包括：

步骤一：对原始探地雷达数据做曲波正变换获取所有尺度j和旋转角度θ系数分布，并提取大尺度和小尺度的曲波系数

和

步骤二：采用曲波变换数值分析确定参与高阶相关统计量计算的尺度范围；

步骤三：For j＝1,…,Jdo；

For i＝1,…,M-1do；每次相关计算使用到前一个观测点数据，因而只有M-1个观测点数据参与计算。

I按公式

计算相关系数

并按公式

归一化

按公式

进行探地雷达数据叠加

对

重复过程I，获取

的叠加高阶相关统计量；

End；

i＝i+1；θ＝θ+1；j＝j+1；

End；

步骤四：对经过曲波域高阶相关统计量叠加计算

和原有大尺度曲波系数

整合为完整的曲波系数，并采用曲波反变换，获取最终噪声压制的探地雷达数据。

本发明的另一目的在于提供一种实现所述曲波域统计量自适应阈值探地雷达数据去噪方法的计算机程序。

本发明的另一目的在于提供一种实现所述曲波域统计量自适应阈值探地雷达数据去噪方法的信息数据处理终端。

本发明的另一目的在于提供一种计算机可读存储介质，包括指令，当其在计算机上运行时，使得计算机执行所述的曲波域统计量自适应阈值探地雷达数据去噪方法。

本发明的另一目的在于提供一种曲波域统计量自适应阈值探地雷达数据去噪控制***包括：

复数域阈值函数分析单元，用于引入块状复数域阈值函数算法，分析传统阈值函数曲波变换去噪效果随阈值函数控制系数变化规律；

曲波变换系数分布尺度单元，用于利用高阶相关统计量理论，对曲波变换系数在尺度、方向上进行相关性叠加，通过相关性统计量自适应确定有效信号在曲波变换系数分布尺度、旋转方向；

统计量自适应阈值函数曲波变换去噪算法构建单元，用于确定出清除噪声成分阈值范围，构建统计量计算自适应阈值函数曲波变换去噪算法。

本发明的另一目的在于提供一种搭载所述曲波域统计量自适应阈值探地雷达数据去噪控制***的地质灾害监测设备。

本发明的另一目的在于提供一种搭载所述曲波域统计量自适应阈值探地雷达数据去噪控制***的工程与环境地质勘察设备。

本发明的另一目的在于提供一种搭载所述曲波域统计量自适应阈值探地雷达数据去噪控制***的水文地质勘查设备。

综上所述，本发明的优点及积极效果为：

本发明引入块状复数域阈值函数算法，分析传统阈值函数曲波变换去噪效果随阈值函数控制系数变化规律；利用高阶相关统计量理论，对曲波变换系数在尺度、方向上进行相关性叠加，通过相关性统计量计算自适应确定有效信号在曲波变换系数分布尺度、旋转方向，由此确定清除噪声成分阈值范围，构建统计量自适应阈值函数曲波变换去噪算法。通过传统阈值函数曲波变换去噪与本发明提出去噪算法对包含随机噪声、相关噪声合成探地雷达数据及实测探地雷达数据处理结果对比，检验了本发明算法的有效性及可行性。研究成果对复杂探地雷达数据精确推断解译具有指导意义。

本发明基于传统曲波变换探地雷达数据去噪需估计阈值函数问题，对探地雷达数据做二维的曲波正变换获取不同尺度、不同旋转角度的系数分布；对每一道探地雷达数据在不同尺度、不同角度下计算高阶相关统计量、根据高阶相关统计量分布明确阈值分布权重，开展了高阶相关统计量自适应函数算法研究，实现不需要估计曲波变换阀值函数及其所属系数范围，而是通过统计理论建立目标体探地雷达有效信息所属的尺度及旋转角度范围，从而构建了统计量自适应阈值函数曲波变换去噪算法。

基于块状复数域阈值函数理论，设计简单矩形及埋深不同双矩形模型，含随机噪声、相关噪声合成探地雷达数据传统阈值函数曲波变换去噪分析表明：块状复数域阈值范围需给定阈值函数控制系数，去噪效果的优劣取决于与含噪数据匹配的阈值函数控制系数；理论合成含噪数据可准确估计阈值函数、实现高效去噪，但实测含噪数据难以通过估计阈值函数去噪实现高效保真去噪目的。针对复杂噪声源环境下采集的探地雷达实测数据，开展微弱绕射波双曲线同相轴特征提取，中深部强幅值、弱幅值似平行非连续性反射波组能量恢复及弱幅值零散反射波组分析处理，获得相应去噪结果。将相应结果与传统阈值函数曲波变换去噪结果对比，分析了高阶统计量算法在复杂噪声源弱信号提取方面具有高效保真弱反射波信号、有效去除复杂噪声、恢复强反射波组同相轴等优势，有助于探地雷达探测资料进行可靠、准确的解译。

附图说明

图1是本发明实施例提供的曲波域统计量自适应阈值探地雷达数据去噪方法流程图。

图2是本发明实施例提供的含噪声模拟合成探地雷达数据曲波变换阈值去噪示意图。

图2中：(a)、原始合成数据；(b)、含随机噪声数据，PSNR＝17dB；(c)、含相关噪声数据，PSNR＝15.51dB；(d)、去噪数据PSNR值随阈值控制系数δ取值变化曲线。

图3是本发明实施例提供的含噪声数据(图2(b)、(c))阈值控制系数为0.08曲波变换去噪结果图。

图3中：(a)、随机噪声数据去噪结果，PSNR＝20.23dB，提高3.23dB；(b)相关噪声数据去噪结果，PSNR＝16.75dB，提高1.24dB。

图4是本发明实施例提供的含噪声数据(图2(b)、(c))统计量自适应阈值曲波变换去噪结果图。

图4中：(a)、随机噪声数据去噪结果，PSNR＝25.3dB，提高8.3dB；(b)、相关噪声数据去噪结果，PSNR＝21.92dB，提高6.41dB。

图5是本发明实施例提供的双矩形目标模型含噪声数据统计量自适应阈值曲波变换去噪结果图。

图5中：(a)、原始合成数据；(b)、含随机噪声数据，PSNR＝16.04dB；(c)、含相关噪声数据，PSNR＝15.7dB；(d)、随机噪声数据去噪结果，PSNR＝23.97dB；(e)、相关噪声数据去噪结果，PSNR＝21.05d。

图6是本发明实施例提供的双矩形目标模型含噪声数据第50道统计量自适应阈值曲波变换去噪结果图。

图6中：(a)、随机噪声数据去噪结果；(b)、相关噪声数据去噪结果。

图7是本发明实施例提供的实测探地雷达时间剖面传统曲波变换去噪和本发明去噪算法处理结果图。

图7中：(a)、实测时间剖面；(b)、采用L2标准方差估计阈值曲波去噪结果；(c)、统计量自适应阈值曲波去噪结果。

图8是本发明实施例提供的实测探地雷达时间剖面第200道传统曲波变换去噪和本发明去噪算法处理结果图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

当前，曲波变换配置研发了诸多计算阀值的方法，如分段线性滤波方法、L2标准差法、曲波正变换法、对角实数及复数阀值法等，需根据实际数据噪声类型选用，总体上无法满足自适应保真去噪的目的。

如图1，本发明实施例提供的曲波域统计量自适应阈值探地雷达数据去噪方法，包括：

S101:引入块状复数域阈值函数算法，分析传统阈值函数曲波变换去噪效果随阈值函数控制系数变化规律；用于后续曲波域统计量自适应阈值对比；

S102:曲波域统计量即利用高阶相关统计量理论，对曲波变换系数在尺度、方向上进行相关性叠加，通过相关性统计量自适应确定有效信号在曲波变换系数分布尺度、旋转方向；

S103:确定出清除噪声成分阈值范围，构建统计量自适应阈值函数曲波变换去噪算法。

本发明实施例提供的曲波域统计量自适应阈值探地雷达数据去噪控制***包括：

复数域阈值函数分析单元，用于引入块状复数域阈值函数算法，分析传统阈值函数曲波变换去噪效果随阈值函数控制系数变化规律；

曲波变换系数分布尺度单元，用于利用高阶相关统计量理论，对曲波变换系数在尺度、方向上进行相关性叠加，通过相关性统计量自适应确定有效信号在曲波变换系数分布尺度、旋转方向；

统计量自适应阈值函数曲波变换去噪算法构建单元，用于确定出清除噪声成分阈值范围，构建统计量自适应阈值函数曲波变换去噪算法。

下面结合具体分析对本发明的应用作进一步描述。

1方法原理

1.1曲波变换

本发明拟采用第二代截断离散曲波变换算法(Candès and Donoho,2004；Candèset al.,2006)开展去噪算法研究。曲波变换变量包含频率、尺度及方位(角度)，其变换表达式为：

则频率窗口U、尺度窗口W和角度窗口V可表述为：

U_j(r,θ)＝2^-3/4W(2^-jr)V(2[j/2]·θ/2π) (2)

式中，j为尺度，m为旋转方向，r和t分别为空间及时间域参数。角度窗口V为环形域，并满足|x|∈[2^j,2^j+1]及极坐标定义θ_j,m＝2πm·2^-[j/2]。对尺度j、旋转角度θ_j,m及空间位置

对f(x)∈L²(R²)，曲波系数定义为：

曲波域噪声压制主要思路为：1)对探地雷达数据进行二维快速傅里叶正变换获取系数分布，2)按公式(2)-(4)配置频率、尺度和角度窗口，3)根据给定阀值函数截断去噪、或将噪声所属系数范围设置为零实现噪声压制，4)对3)部分处理剩余曲波系数利用二维快速傅里叶反变换算法获取处理数据结果。曲波域噪声压制技术的关键是获取较为准确的阀值函数及其所属系数范围。对于合成探地雷达数据，可根据加入的噪声类型、强度给出准确的阀值函数及确定其所属系数范围；对于实测探地雷达数据，常规操作是根据实际情况，人为选取合适的估计阀值函数及所属系数范围。

第二代曲波变换算法(Candès et al,2004；2006)开展含噪数据处理对比研究，其中，曲波变换涉及的阈值窗口在所有尺度设置为

而δ为阈值分布范围控制系数，可人为设定；E_c为曲波变换系数L2范数。δ值越大，过滤窗口越大，残余信息越少；反之，δ值越小，噪声及有效信息量越大。关于数据过滤判断范围的估计，有多种方法供参考，如计算数据L2范数、块状复数阈值估计等(Saha et al,2015)。对于复杂数据处理，块状复数阈值函数更具有优势，因而本发明引入块状复数域阈值用于曲波变换去噪过程。其表达式如式(6)所示：

式中，Ψ为复数域阈值函数，用于估计判断滤除范围，具体推导请参考文献(Sahaet al,2015)，Γ₁为相邻道曲波系数权重值，Γ₂为曲波系数开平方归一化和权重值。

1.2曲波域统计量自适应阈值去噪

对探地雷达数据做二维的曲波正变换获取不同尺度、不同旋转角度的系数分布；对每一道探地雷达数据的在不同尺度、不同角度下计算高阶相关统计量、根据所述高阶相关统计量分布明确阈值分布权重，实现噪声滤除；最后采用二维的曲波反变换重建去噪后的探地雷达数据；其中，未偏移的三阶相关函数表示为：

归一化高阶相关统计量为：

式中，i为探地雷达数据接收点，t为时间采样点，q为平移因子，取值为1，

为第i道探地雷达数据曲波域第k个尺度、第θ个方向正变换系数；

为第i道探地雷达数据曲波域第k个尺度、第θ个方向二阶自相关函数值；

为第i+1道探地雷达数据曲波域第k个尺度、第θ个方向二阶自相关函数值；

所述通过高阶相关统计量计算自适应确定有效信号在曲波变换系数分布尺度、旋转方向具体包括：

步骤一：对原始探地雷达数据做曲波正变换获取所有尺度j和旋转角度θ系数分布，并提取大尺度和小尺度的曲波系数

和

步骤二：采用曲波变换数值分析确定参与高阶相关统计量计算的尺度范围；

步骤三：For j＝1,…,Jdo；

For i＝1,…,M-1do；每次相关计算使用到前一个观测点数据，因而只有M-1个观测点数据参与计算。

I按公式

计算相关系数

并按公式

归一化

按公式

进行探地雷达数据叠加

对

重复过程I，获取

的叠加高阶相关统计量；

End；

i＝i+1；θ＝θ+1；j＝j+1；

End；

步骤四：对经过曲波域高阶相关统计量叠加计算

和原有大尺度曲波系数

整合为完整的曲波系数，并采用曲波反变换，获取最终噪声压制的探地雷达数据。

2模拟试验

2.1传统曲波变换去噪

建立理论模型并采用二维有限差分程序计算合成模拟数据。矩形目标体规模大小为2×1m，X轴分布-1～1m，Z轴分布为2～3m；电阻率为10Ω.m、介电常数为30F/m；背景介质电阻率为1000Ω.m、介电常数为3F/m；探地雷达观测***频率为100MHz，点距为0.1m，X轴方向设置测点数为101，时间域采集点数为241。文中峰值信噪比定义为

其中MSE为数据均方差，表示原始信号和噪声(或去噪后)信号的近似程度，max(s)为原始信号的峰值。

图2(a)所示为无噪声模拟合成探地雷达原始时间域数据，可见空间域设定目标矩形模型边界引起的探地雷达反射波信号分布于不同时间域。图2(b)和图2(c)分别为加入随机噪声(PSNR＝17dB)、相关噪声探地雷达数据(PSNR＝15.51dB)，随机噪声椒盐式无规律分布于整个探地雷达剖面，弱目标反射信号完全被噪声淹没；相关噪声根据有效反射信号范围比例分布，可见部分假反射信号。采用公式(6)阈值函数曲波变换去噪数据的PSNR值随阈值控制系数δ取值变化曲线如图2(d)所示。随着阈值控制系数δ取值由小变大，去噪数据PSNR值趋于极值后逐渐降至原噪声数据PSNR值水平，其中随机噪声数据去噪结果PSNR值在阈值控制系数δ＝0.15(与原始数据噪声水平一致)附近取得PSNR极值(27.33dB)，相关噪声对应在δ＝0.35附近取得PSNR极值(23.27dB)。可见，采用阈值函数曲波变换去噪效果好坏取决于阈值控制系数(原始数据噪声水平)取值是否合理。同时，随机噪声与相关噪声数据采用阈值函数曲波变换去噪的合理阈值控制系数规律不同，前者阈值控制系数取噪声水平值(δ＝0.10～0.20)，而后者取略高于阈值控制系数值(δ＝0.3～0.4)。

针对含噪模拟数据阈值函数曲波变换去噪，预先获取的准确阈值控制系数范围可有效用于含噪数据处理。但实测数据常常被不同类型、不同强度的噪声源污染，合理的阈值控制系数范围难以确定，而为了有效使用阈值函数曲波变换去噪算法，需要通过特定方式人为估计阈值控制系数范围，如L2标准方差算法。图3(a)-(b)所示为含噪声数据(图2(b)、(c))采用估计方法确定阈值控制系数为0.08时曲波变换去噪结果。相比较原始含噪数据(图2(b)、(c))，阈值函数曲波变换去噪清除了部分噪声成分信号，但目标矩形模型引起的有效反射波信号仅依稀可辩，无法有效用于全波形反演成像计算及后续推断解译过程。

2.2统计量自适应阈值曲波变换去噪

为避免传统曲波变换去噪需要估计合理阈值函数范围，引入统计量自适应阈值算法曲波变换对图2(a)-(d)所示随机噪声及相关噪声数据去噪处理，如图4(a)-(b)所示。从去噪结果可见，基于统计量自适应阈值曲波变换去噪方法不仅可有效衰减随机噪声(PSNR值提高8.3dB)、相关噪声(PSNR值提高6.41dB)，还可以较好地还原剖面中有效反射波信号空间及时间域分布。完整噪声衰减过程中，无需估计阈值函数范围，通过采用公式(8)计算统计量自适应阈值权重，达到了保真去噪的目的。其中，随机噪声数据去噪效果优于相关噪声数据去噪效果，究其缘由，相关噪声数据信噪比低(PSNR＝15.51dB)导致噪声信号相关性统计量计算误差，残留部分相关性较强的噪声信号。总体上，统计量自适应阈值算法曲波变换去噪结果可有效用于数据全波形偏移成像处理及后续解释推断过程。

设计双矩形目标体规模大小1×1m，水平位置为-2～-3m及2～3m，埋深范围为2～3m及3～4m；其他物性参数及探地雷达观测***与前述模型一致。数值模拟双矩形目标体探地雷达时间剖面如图5(a)所示，加入随机噪声及相关噪声数据如图5(b)、(c)所示。含随机噪声数据(PSNR＝16.04dB)及相关噪声数据(PSNR＝15.7dB)完全淹没目标体有效反射波信号，无法用于目标体分布推断解释。

图5(d)、(e)为采用本发明提出去噪算法处理结果，随机噪声时间剖面内双矩形目标体引起复杂有效反射波信号被完全恢复，PSNR值较原始含噪数据提高7.93dB；相关噪声数据基本重建目标体有效反射波信号时空分布，但仍残余部分相关噪声成分，PSNR值较原始含噪数据提高6.65dB。图5第50道原始无噪声数据、含噪数据及去噪数据信号细节对比曲线如图5(a)-(e)所示，本发明提出的去噪算法有效清除了全局椒盐式随机噪声分布(图6(a))、准确分辨相关噪声环境下微弱有效反射波信号(图6(b))，去噪数据曲线与原始无噪声数据曲线吻合较好，验证了本发明提出的去噪算法有效性及可行性。

下面结合实测算例对本发明的应用作进一步描述。

图7(a)所示为某地探地雷达实测测线原始时间剖面图，可见受场地环境影响，椒盐式随机噪声遍布整个剖面数据；个别采集道受局部不均匀体影响，出现相邻观测道幅值畸变。原始数据剖面浅部(100～250时间采样点)介质的反射波能量非常微弱，依稀可见几处绕射波同相轴，但受噪声淹没，双曲线特征不明显；250～300时间采样范围中深部出现较强幅值似平行反射波组能量，同相轴断断续续并不连续；强反射波组下(300～400时间采样点)较弱零散反射波组，难以识别强反射波组能量下边界。

利用上述传统曲波变换去噪处理方法对图7(a)进行数据处理，通过L2标准方差估计阈值范围，去噪获得的雷达剖面如图7(b)所示。由7(b)可见，部分地下反射波组信号能量得以加强，随机噪声源一定程度上被清除；但原始数据内有效绕射波双曲线特征、断续反射波组依然无法有效识别，去噪数据信噪比提高有限。究其原因，实测数据中随机噪声、相关噪声及其他类型噪声源掺杂，同时，原始数据包含幅值强度各异反射波组，估计阈值范围无法有效覆盖上述特点数据去噪范围。人工试错法确定最佳阈值范围或能有效清除噪声，但最佳去噪效果评价缺乏科学依据。由此，仅对经过传统阈值函数曲波变换去噪数据处理后的雷达剖面进行分析解释，无法保证资料解译的准确度。

采用本发明提出的统计量自适应阈值曲波变换去噪算法对图7(a)原始数据进行处理，去噪结果如图7(c)所示。由于随机噪声源在曲波变换系数尺度、方向上并不具备统计相关特性；相关噪声源变换系数在尺度、方向上虽具备一定相关性，但相对有效反射波组信号相关性仍然较小，因而统计量自适应阈值可有效分辨反射波组能量，去除噪声信号，达到保真去噪的目的。由图7(c)可见，地下反射波能量得到明显加强，特别是浅部多个绕射波双曲线同相轴特征的信噪比得到有效提高，中深度多层介质的反射组同相轴连续、分界面明显，且来自于强反射波组下弱幅值的反射波能清晰可见。经过统计量自适应阈值曲波变换数据处理后的雷达剖面进行分析解释，有助于资料解译的准确度。

图7(a)-(c)实测探地雷达数据第200道传统曲波变换去噪、本发明提出的去噪结果细节对比曲线如图8所示。由曲线对比可见，0～100时间采样范围内地面回波、弱幅值绕射波(100～150)、强幅值绕射波(150～230)及强反射波组(230～350)受到弱幅值的随机噪声及相关噪声污染，大于350时间采样范围弱反射波信号受到较强幅值的上述噪声源淹没。相对较传统曲波变换去噪结果，本发明提出的去噪算法可有效恢复弱幅值绕射波(100～150)、强幅值绕射波(150～230)及强反射波组(230～350)同相轴特征，同时可清晰分辨弱反射波信号，验证了本发明提出去噪算法的可行性及有效性。

下面结合效果对本发明的应用作进一步描述。

基于传统曲波变换探地雷达数据去噪需估计阈值函数问题，对探地雷达数据做二维的曲波正变换获取不同尺度、不同旋转角度的系数分布；对每一道探地雷达数据在不同尺度、不同角度下计算高阶相关统计量、根据统计量分布明确阈值分布权重，开展了统计量自适应函数算法研究，实现不需要估计曲波变换阀值函数及其所属系数范围，而是通过统计理论建立目标体探地雷达有效信息所属的尺度及旋转角度范围，从而构建了统计量自适应阈值函数曲波变换去噪算法。

基于块状复数域阈值函数理论，设计简单矩形及埋深不同双矩形模型，含随机噪声、相关噪声合成探地雷达数据传统阈值函数曲波变换去噪分析表明：块状复数域阈值范围需给定阈值函数控制系数，去噪效果的优劣取决于与含噪数据匹配的阈值函数控制系数；理论合成含噪数据可准确估计阈值函数、实现高效去噪，但实测含噪数据难以通过估计阈值函数去噪实现高效保真去噪目的。

针对复杂噪声源环境下采集的探地雷达实测数据，开展微弱绕射波双曲线同相轴特征提取，中深部强幅值、弱幅值似平行非连续性反射波组能量恢复及弱幅值零散反射波组分析处理，获得相应去噪结果。将相应结果与传统阈值函数曲波变换去噪结果对比，分析了高阶统计量算法在复杂噪声源弱信号提取方面具有高效保真弱反射波信号、有效去除复杂噪声、恢复强反射波组同相轴等优势，有助于探地雷达探测资料进行可靠、准确的解译。

在上述实施例中，可以全部或部分地通过软件、硬件、固件或者其任意组合来实现。当使用全部或部分地以计算机程序产品的形式实现，所述计算机程序产品包括一个或多个计算机指令。在计算机上加载或执行所述计算机程序指令时，全部或部分地产生按照本发明实施例所述的流程或功能。所述计算机可以是通用计算机、专用计算机、计算机网络、或者其他可编程装置。所述计算机指令可以存储在计算机可读存储介质中，或者从一个计算机可读存储介质向另一个计算机可读存储介质传输，例如，所述计算机指令可以从一个网站站点、计算机、服务器或数据中心通过有线(例如同轴电缆、光纤、数字用户线(DSL)或无线(例如红外、无线、微波等)方式向另一个网站站点、计算机、服务器或数据中心进行传输)。所述计算机可读取存储介质可以是计算机能够存取的任何可用介质或者是包含一个或多个可用介质集成的服务器、数据中心等数据存储设备。所述可用介质可以是磁性介质，(例如，软盘、硬盘、磁带)、光介质(例如，DVD)、或者半导体介质(例如固态硬盘SolidState Disk(SSD))等。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (8)

1.一种曲波域统计量自适应阈值探地雷达数据去噪方法，其特征在于，所述曲波域统计量自适应阈值探地雷达数据去噪方法包括：

引入块状复数域阈值函数算法，分析传统阈值函数曲波变换去噪效果随阈值函数控制系数变化规律，用于后续曲波域统计量自适应阈值对比；

曲波域统计量即利用高阶相关统计量理论，对曲波变换系数在尺度、方向上进行相关性叠加，通过高阶相关统计量计算自适应确定有效信号在曲波变换系数分布尺度、旋转方向；

确定出清除噪声成分阈值范围，构建高阶相关统计量自适应阈值函数曲波变换去噪算法；

通过高阶相关统计量计算自适应确定有效信号在曲波变换系数分布尺度、旋转方向；包括：

对探地雷达数据做二维的曲波正变换获取不同尺度、不同旋转角度的系数分布；对每一道探地雷达数据的在不同尺度、不同角度下计算高阶相关统计量、根据所述高阶相关统计量分布明确阈值分布权重，实现噪声滤除；最后采用二维的曲波反变换重建去噪后的探地雷达数据；其中，未偏移的三阶相关函数表示为：

归一化高阶相关统计量为：

式中，i为探地雷达数据接收点，t为时间采样点，q为平移因子，取值为1，F_i ^k,θ为第i道探地雷达数据曲波域第k个尺度、第θ个方向正变换系数；

为第i道探地雷达数据曲波域第k个尺度、第θ个方向二阶自相关函数值；

为第i+1道探地雷达数据曲波域第k个尺度、第θ个方向二阶自相关函数值；

所述通过高阶相关统计量计算自适应确定有效信号在曲波变换系数分布尺度、旋转方向具体包括：

步骤一：对原始探地雷达数据做曲波正变换获取所有尺度j和旋转角度θ系数分布，并提取大尺度和小尺度的曲波系数

和

步骤二：采用曲波变换数值分析确定参与高阶相关统计量计算的尺度范围；

步骤三：For j＝1,…,J do；

For i＝1,…,M-1 do；每次相关计算使用到前一个观测点数据，因而只有M-1个观测点数据参与计算；

过程I：按公式

计算获得相关系数

按公式

归一化相关系数

得到结果为

按公式

下标索引i的变化进行探地雷达数据叠加获得结果

对

重复过程I，获取

的叠加高阶相关统计量；

End；

i＝i+1；θ＝θ+1；j＝j+1；

End；

步骤四：对经过曲波域高阶相关统计量叠加计算

和原有大尺度曲波系数

整合为完整的曲波系数，并采用曲波反变换，获取最终噪声压制的探地雷达数据。

2.一种实现权利要求1所述曲波域统计量自适应阈值探地雷达数据去噪方法的计算机程序。

3.一种实现权利要求1所述曲波域统计量自适应阈值探地雷达数据去噪方法的信息数据处理终端。

4.一种计算机可读存储介质，包括指令，当其在计算机上运行时，使得计算机执行如权利要求1所述的曲波域统计量自适应阈值探地雷达数据去噪方法。

5.一种实现权利要求1所述曲波域统计量自适应阈值探地雷达数据去噪方法的曲波域统计量自适应阈值探地雷达数据去噪控制***，其特征在于，所述曲波域统计量自适应阈值探地雷达数据去噪控制***包括：

复数域阈值函数分析单元，用于引入块状复数域阈值函数算法，分析传统阈值函数曲波变换去噪效果随阈值函数控制系数变化规律，用于后续曲波域统计量自适应阈值对比；

曲波变换系数分布尺度单元，用于利用高阶统计量理论，对曲波变换系数在尺度、方向上进行相关性叠加，通过高阶相关统计量计算自适应确定有效信号在曲波变换系数分布尺度、旋转方向；

统计量自适应阈值函数曲波变换去噪算法构建单元，用于确定出清除噪声成分阈值范围，构建统计量自适应阈值函数曲波变换去噪算法。

6.一种搭载权利要求5所述曲波域统计量自适应阈值探地雷达数据去噪控制***的地质灾害监测设备。

7.一种搭载权利要求5所述曲波域统计量自适应阈值探地雷达数据去噪控制***的工程与环境地质勘察设备。

8.一种搭载权利要求5所述曲波域统计量自适应阈值探地雷达数据去噪控制***的水文地质勘查设备。

Images