基于Curvelet变换的隧道裂隙水GPR数据处理研究

朱自强，朱 贺*，鲁光银，王凡，谭 洁

（1．中南大学地球科学与信息物理学院，长沙 410083；2．中国水电顾问集团贵阳勘测设计研究院，贵阳 550081）

基于Curvelet变换的隧道裂隙水GPR数据处理研究

朱自强1，朱 贺1*，鲁光银1，王凡2，谭 洁1

（1．中南大学地球科学与信息物理学院，长沙 410083；2．中国水电顾问集团贵阳勘测设计研究院，贵阳 550081）

应用地质雷达探测隧道内部裂隙水的分布状况时，受隧道钢筋网的影响以及噪声干扰，导致探测信号中弱有效信号难以识别。Curvelet变换可以对二维信号从频率、角度和空间位置实现有效反射波和干扰波的分离及降噪处理，应用Curvelet变换对含钢筋网干扰的裂隙水模型的探地雷达模拟数据进行波场分离和降噪处理，在压制直达波、去除钢筋层反射信号及多次反射的基础上，可以有效地提取出含裂隙水所引起的反射波信号。将该方法应用到汝郴高速某隧道实测探地雷达数据中，较好地去除了钢筋层反射信号的影响，准确地判断出了裂隙水的空间分布位置，证明了所提出方法的有效性。

Curvelet变换；探地雷达；隧道裂隙水；波场分离；降噪

0 引言

探地雷达法［1－2］是近些年来迅速发展起来的一种高分辨率、高效率的勘探方法，目前在隧道勘测中，特别是探测岩溶裂隙水方面应用较多［3－6］。由于隧道周围裂隙水被密集的钢筋层覆盖，用探地雷达对隧道裂隙水探测时会受到钢筋层反射的干扰，因此，需要对雷达数据进行处理，去除干扰和噪声，从而提取有效的异常信号。刘斌等［7］将复信号技术应用到探地雷达对地下裂隙水的探测中，对雷达数据中的异常分别从振幅、频率和相位等各方面进行综合分析，提高了预报精度和准确性，但没能从根本上克服探地雷达法抗干扰能力弱的缺点。刘四新等［8］用钻孔探地雷达对岩溶裂隙进行勘探，勘探效果较好，但具有探测上方位不确定性和成本较高的缺点。邹海林等［9］基于小波变换对探地雷达信号进行处理，李才明［10］基于小波能谱分析对岩溶区探地雷达进行目标识别。小波变换很好地解决了时频同步问题，而且对一维分段平稳信号的表达具有最佳性能，但是对于描述二维图像信息却存在诸多缺点［11－13］。

Candes等人［14－15］在1999年提出了第一代Curvelet变换，并于2002年提出了实现更简单、运行速度更快、更便于直观理解的第二代Curvelet变换算法。在此之后，他们又于2005年提出了两种基于第二代Curvelet 变换理论的快速离散实现方法［16］，即USFFT（unequally-spaced fast Fourier transforms）和Wrapping两种快速离散算法。2008年以来，该方法在地球物理资料处理尤其是地震信号处理中得到了迅速地研究和应用［17－20］。在国内，仝中飞等人［21］提出了迭代阈值法压制地震信号中的随机噪声，取得了很好的效果。潘雪辉［22］提出了一种用Curvelet变换压制地震资料中线性噪声的方法，为地震信号中线性噪声的压制问题提供了较好地解决方法。

Curvelet变换建立在小波变换的基础上，增加了一个方位参数，解决了小波变化在处理二维信号时的不足，而探地雷达检测数据主要为描述一定深度范围内剖面的二维信号，这让作者萌生思路，利用Curvelet变换对探地雷达数据进行处理研究。在对Curvelet变换方法研究基础上，将其应用到隧道裂隙水的探地雷达检测资料处理中。本研究首先模拟了在含有覆盖钢筋层和裂隙水的混凝土中的探地雷达波场，并通过Curvelet变换对模拟波场进行去除直达波、提取水层信息、降噪等一系列数据处理流程，提取出有效的裂隙水异常信息；随后将该方法应用到汝郴高速某隧道实测地质雷达数据，获取了隧道内裂隙水异常的空间分布，通过钻孔注浆验证，证明了此方法的有效性。

1 Curvelet变换基本原理

小波变换在信号处理中的应用得到了很大的发展，遗憾的是，由一维小波所生成的可分离小波只具有有限的方向，因此在处理图像边缘时的效果不是很好。Curvelet变换除了与小波变换一样具有尺度和位移参数，还增加了一个方位参数，因此具有更好的方位识别能力，可以精确表达出图像中边缘的方向信息。

1．1 连续Curvelet变换

Curvelet变换是利用基函数与信号（或函数）的内积形式来实现信号（或函数）的稀疏表示，表示为：

其中：φj，k，l表示Curvelet函数；j、k、l分别表示尺度、位置和方向参数。

Curvelet变换在频域内采用窗函数U来实现的。分别定义径向窗函数W（r）、r∈（1／2，2）和角度窗函数V（t），t∈［－1，1］，二者满足式（2）与式（3）。

对于每一个j≥j0，在频域中定义窗函数Uj如式（4）。

综上所述，我们可以得到Curvelet变换公式：

图1 Curvelet变换示意图Fig．1 Curvelet transform schematic diagram

1．2 离散Curvelet变换

在直角坐标系下f［t1，t2］，0≤t1，t2＜n为输入，Curvelet变换离散形式为：

1．3 离散Curvelet变换快速实现方法

本研究采用基于USFFT算法的快速离散Cur-velet变换方法，实现过程如下：

1）对于给定的直角坐标系下二维函数f［t1，t2］，0≤t1，t2＜ω进行2DFFT，得到其二维频率域表达式：

2）在频率域，对每一对（i，j），对函数＾f［n1，n2］重新采样，得到函数采样值：

其中：Pj＝｛（n1，n2）：n1，0≤n1≤n1，0＋L1，j，n2，0≤n2≤n2，0＋L2，j｝；L1，j表示窗函数支撑区间的长度；L2，j表示支撑区间的宽度。

1．4 Curvelet变换的非线性逼近能力

Curvelet变换之所以比小波变换更适合表示二维信号，其主要原因在于Curvelet变换比小波在二维平面具有更好的稀疏性。假设f∈L2［0，1］为一条曲线的光滑部分，分别用小波变换和Curvelet变换对其分解，再分别选择最佳的M个系数对函数进行重构，得到三种变换下的的最佳逼近它们的逼近误差分别为：

从重构误差式中可以看出，Curvelet变换的逼近率要优于小波变换。Curvelet变换之所在表示二维曲线时有更好的逼近率，主要在于它在小波变换的基础上增加了方位参数，使得它具有多方向性、多尺度性、波动性和各向异性的特点，这也使得用更少的Curvelet系数就能表示二维信号的主要特征。

图2 小波和Curvelet对二维曲线逼近示意图Fig．2 Schematic of two－dimensional curve approximation with wavelet and Curvelet

2 模拟数据处理分析

针对表层存在钢筋干扰下的含裂隙水混凝土模型进行数值模拟，模型大小为5．2 m×1．28 m，混凝土介电常数ε0为9，电导率σ0为0．01 s／m。距离表层0．3 m处设置一排钢筋，半径0．01 m。钢筋层下方设置一个水平和一个倾斜含裂隙水，横截面大小均为0．8 m×0．08 m，介电常数ε1为81，电导率σ1为0．05 s／m，模拟天线频率为400 MHz。模型剖面和雷达模拟剖面如图3所示。

图3 雷达数值模拟Fig．3 GPR numerical simulation

在图3中，表层直达波强度较高，弱化了下面钢筋层和裂隙水层的反射信号。将总波场变换到Curvelet域后，然后在角度窗函数下选择变换域中θ＝0（θ为方向因子）附近的Curvelet系数进行重构，得到直达波的波场信号，进而分离得到的直达波波场信号和剩余波场信号（图4）。

图4 波场分离Fig．4 Wave field separation

当直达波去除后，下面钢筋层和裂隙水层的反射信息加强，但是二次反射和干扰信息的存在干扰了异常信号，通过Curvelet变换的降噪功能可以去掉二次反射及其他干扰信号。基于Curvelet变换的特点，前人提出了各种关于降噪的阈值方法［23－24］，本研究采取一种自适应阈值方法［25］，即随着尺度的改变，每一层阈值可以相应地改变。阈值函数表达式如下：其中 sigma为噪声的方差估计值；N为Curvelet域系数矩阵大小；s为Curvelet域尺度分解下的层数。

图5 波场降噪Fig．5 Wavefield denoising

图5为sigma系数设为700时降噪的效果图，降噪后，钢筋和裂隙水的异常信息已非常明显。原始二维信号在经过Curvelet变换后会得到的一系列Curvelet系数矩阵，将这些系数矩阵在MATLAB中生成时间－频率图像，通过对比原始图像异常位置可以识别出包含有效信号和干扰信号的矩阵。如果我们只需要裂隙水的异常信息，只需对这些离散化的系数矩阵进行人工识别，并挑选出含裂隙水信息的系数矩阵，进行Curvelet反变换，得到只含裂隙水异常信息的二维数据。图6为应用Curvelet变换波场分离方法得到的裂隙水异常信息数据。

图6 提取裂隙水波场信息Fig．6 Wave signal extraction of fissure water

处理结果显示，用Curvelet变换对雷达模拟数据的处理效果良好，波场分离和降噪都能够基本达到目的。两个方向有差异的裂隙水的异常信息在处理过程中表现出不同的处理效果，这也说明了Curvelet变换在处理图像信息时带有很强的方向选择性。由于水平的裂隙水在方向上与上面钢筋层保持平行，用Curvelet变换离散化图像信息时，水平裂隙水与钢筋层的信息在角度窗函数选择下没能够完全分离，只能够通过尺度来分离部分信息。这样倾斜的裂隙水异常信息才会显得比水平裂隙更明显。

3 实测数据处理分析

在Curvelet变换方法研究和模拟处理的基础上，对汝郴高速公路某隧道内采集所得探地雷达数据进行处理。探测目的是找出隧道边墙和底部存在的裂隙水，确定这些裂隙水的位置及相关信息，为工程治理工作提供较好的参考信息。探测中使用的天线为400 MHz，采样点数设为512，参考设计方案，钢筋间距约为50 cm。

图7（a）为截取采集数据中某一异常位置信息，直达波和钢筋产生的干扰信息强度较强，同时数据中含有一些噪声干扰。首先我们要将有效异常信息和钢筋层干扰信息及直达波信息分离开，提取含有有效异常信息的数据。图7（b）为分离出的钢筋层信息及一些干扰信息，图7（c）为提取得到的含有裂隙水层异常信息的波场。由于波场分离更多在于角度窗下的选择，图7（b）中钢筋层干扰和直达波与图7（c）中有效信息被基本分离开，但仍然有一些因为与裂隙水异常信息方向一致而未被分离的遗留信息，且存在一些噪声。遗留下的干扰信息强度较弱，可以通过Curvelet降噪功能，应用自适应尺度阈值法，通过尺度窗的选择进行降噪，取sigma系数为48 000得到降噪后的异常图7（d）。图7（d）异常信息清晰，进而可以判断裂隙水的空间位置。

根据探测和处理结果，对发现的异常位置进行钻孔验证如图8所示。验证结果支持了该方法对于裂隙水空间位置的判断，并在该异常位置进行了注浆治理，总注浆量为10 T左右水泥和38 T左右化学浆。对于本次探测中发现的其他几个异常位置，用文中所述方法进行分析并打钻注浆，在判断得到的异常较大位置都能注入不同量的浆液，验证了该处理方法的准确性与可行性。

4 结论

1）基于Curvelet变换，从波场分离的角度来提高探地雷达异常信号的辨识度，较好地解决了直达波和钢筋层强反射对下层异常信息的干扰问题。

2）角度参数的增加使得Curvelet变换在二维信号的分离上有了很强的优势，但是在干扰信息与有效信息方向性一致的情况下，分离效果并不理想。

图7 工程应用实例Fig．7 Examples of application engineering

图8 现场钻孔照片Fig．8 Photos of drilling

3）阈值函数的选取直接关系到Curvelet降噪的效果，采取适应不同尺度的自适应阈值函数，基本能够达到消除噪声的目的，但阈值函数thresh没有从角度去适应系数矩阵，找到同时适应尺度和角度的阈值函数将是下一步急需解决的问题。

［1］ 李大心．探地雷达方法与应用［M］．北京：地质出版社，1994．

［2］ 杨天春，吕绍林，王齐仁．探地雷达检测道路厚度结构的应用现状及进展［J］．物探与化探，2003，27（1）：79－82．

［3］ 杨龙江，阎满存．探地雷达在隐伏岩溶区工程地质勘探中的应用［J］．物探化探计算技术，1999，21（2）：172－176．

［4］ 葛双成，邵长云．岩溶勘察中的探地雷达技术及应用［J］．地球物理学进展，2005，20（2）：476－481．

［5］ 李孟娟，李川．探地雷达检测隧道衬砌厚度的研究［J］．物探化探计算技术，2008，30（3）：231－234．

［6］ 李仁海，杨磊，许新刚，等．地质雷达探测技术在岩溶地形勘察中的应用［J］．物探化探计算技术，2009，31（5）：442－446．

［7］ 刘斌，李术才，李树忱，等．复信号分析技术在地质雷达预报岩溶裂隙水中的应用研究［J］．岩土力学，2009，30（7）：2191－2196．

［8］ 刘四新，曾昭发，徐波．利用钻孔雷达探测地下含水裂缝［J］．地球物理学进展，2006，21（2）：620－624．

［9］ 邹海林，宁书年，林捷．小波理论在探地雷达信号处理中的应用［J］．地球物理学进展，2004，19（2）：268－275．

［10］ 李才明，王良书，徐鸣洁，等．基于小波能谱分析的岩溶区探地雷达目标识别［J］．地球物理学报，2006，49（5）：1499－1504．

［11］WALKER J S，CHEN Y．Image denoising using treebased wavelet subband correlations and shrinkage［J］．Optical Engineering，2000，39（11）：2900－2908．

［12］杨凤娟，罗省贤．小波变换在探地雷达检测钢筋中的应用［J］．物探化探计算技术，2009，31（4）：354－360．

［13］黄敏，朱德兵，郭政学，等．连续小波变换在探地雷达信号分析中的应用研究［J］．物探化探计算技术，2012，34（5）：593－598．

［14］CANDES E J，DONOHO D L．Curvelets：A surprisingly effective nonadaptive representation for objects with edges［R］．DTIC Document，2000．

［15］STARCK J，CANDèS E J，DONOHO D L．The curvelet transform for image denoising［J］．Image Processing，IEEE Transactions on．2002，11（6）：670－684．

［16］CANDES E，DEMANET L，DONOHO D，et al．Fast discrete curvelet transforms［J］．Multiscale Modeling＆Simulation，2006，5（3）：861－899．

［17］Neelamani R，Baumstein A I，Gillard D G，et al．Coherent and random noise attenuation using the curvelet transform［J］．The Leading Edge．2008，27（2）：240－248．

［18］姜宇东，杨勤勇，何柯，等．基于曲波变换的地面微地震资料去噪方法研究［J］．石油物探，2012，51（6）：620－624．

［19］吴爱弟，赵秀玲．基于曲波变换的地震信号去噪方法［J］．中国石油大学学报：自然科学版，2010，34（003）：30－33．

［20］张云强，张培林，王国德，等．基于曲波变换和色度模型的彩色图像去噪［J］．中国图象图形学报，2012，17（012）：1472－1477．

［21］仝中飞，王德利，刘冰．基于Curvelet变换阈值法的地震数据去噪方法［J］．吉林大学学报：地球科学版，2008：48－52．

［22］潘雪辉．曲波变换在地震信号去噪中的应用［D］．西安：西安理工大学，2009．

［23］HERRMANN F J，HENNENFENT G．Non-parametric seismic data recovery with curvelet frames［J］．Geophysical Journal International，2008，173（1）：233－248．

［24］HERRMANN F J，MOGHADDAM P，STOLK C C．Sparsity-and continuity-promoting seismic image recovery with curvelet frames［J］．Applied and Computational Harmonic Analysis，2008，24（2）：150－173．

［25］刘静．地震图像噪声压制在Curvelet域的阈值选取与噪声定位［D］．成都：西南交通大学，2011．

Processing of GPR data in tunnel fissure water based on Curvelet transform

ZHU Zi-qiang1，ZHU He1*，LU Guang-yin1，WANG Fan2，TAN Jie1
（1．School of Geosciences and Info-Physics，Central South University，Changsha 410083，China；2．HYDROChina Guiyang Engineering Corporation，Guiyang 550081，China）

When we detect the distribution of fissure water in tunnel with ground penetrating radar，it is hard to identify the weak effective signal because of the influence of mesh reinforcement and noise．Curvelet transform can achieve effective wave field separation and denoising of two－dimensional signal from the frequency，angle and spatial location．The Curvelet transform is applied to wave field separation and signal denoise processing by using of GPR simulated data from fissured water which is disturbed seriously by the existence of mesh reinforcement．After the direct signal suppression，removal of the reflected signal and multiple reflections from mesh reinforcement，we can effectively extract the reflected wave signal which is caused by the water bearing structure crack．Our method is applied to real data measured in a tunnel of the Ruchen highway．The results show that after removal of the reflected signal from mesh reinforcement we can accurately determine the location of fissure water，which proves the effectiveness of the proposed method．

Curvelet transform；ground-penetrating radar（GPR）；fissure water in tunnel；wavefield separation；denoising

P 631．4

A

10．3969／j．issn．1001-1749．2014．05．10

1001-1749（2014）05-0571-06

2014-03-33 改回日期：2014-06-22

国家自然科学基金项目（41174061）；中南大学自由探索计划（2011QNZT011）

朱自强（1964-），男，教授，博士生导师，主要从事地质灾害探测与监测工作，E-mail：13507319431＠139．com。

*通讯作者：朱贺（1988-），男，硕士，主要从事地球物理信号处理研究，E-mail：zhuhe24＠sina．com。