S变换在地质雷达隧道衬砌检测数据处理中的应用研究

朱自强, 王玉琴, 鲁光银

(中南大学 地球科学与信息物理学院，长沙 410083)

S变换在地质雷达隧道衬砌检测数据处理中的应用研究

朱自强, 王玉琴, 鲁光银

(中南大学 地球科学与信息物理学院，长沙 410083)

通过将S变换应用于地质雷达探测隧道衬砌的数据处理中，解决了由于衬砌表层密集钢筋网的影响导致雷达数据剖面受干扰严重、深部信号弱、病害异常难以识别的问题。进行了地质雷达二维有限差分数值模拟，研究了浅部密集钢筋网对深部钢拱架和脱空区造成的影响情况，并对原始剖面进行偏移成像压制了绕射波和多次波，提高了钢筋网下方弱反射信号的成像效果；最后通过对模拟数据和实际资料的S变换时频分析，准确清晰地识别出了深部钢筋和脱空区的特征。通过利用S变换时频处理方法，为受干扰严重的隧道衬砌检测数据的处理提供了一种方案。

S变换； 地质雷达； 偏移成像； 隧道衬砌检测

0 引言

隧道衬砌检测是隧道质量安全检测最重要的环节，检测结果影响着整个隧道的施工和运营。地质雷达技术作为一种无损方式的浅层地球物理勘探方法[1-2]，以其分辨率高、图像直观、工作效率高等技术优势，使地质雷达成为隧道衬砌检测的常规仪器。

关于地质雷达的数值模拟研究已经比较成熟[3-4]，在衬砌数值模拟上，肖都[5]模拟了不同雷达天线探测衬砌的效果研究；周斌[6]对雷达衬砌检测的数值模拟和工程实例进行了研究；舒志乐等[7]、黄玲等[8]研究了隧道衬砌病害的探地雷达三维正演模拟及工程验证。但是在实际工程中，施工环境复杂，干扰严重，实际数据与数值模拟还是有一定的差别，所以要注重在复杂环境下的资料处理，数据处理的准确性直接决定了整个检测结果。由于地质雷达受钢筋的影响，整个雷达剖面中存在较多干扰信号，影响正常的数据解释，因此如何在钢筋影响下准确进行异常识别就成为了关键。张家松等[9]研究了希尔伯特变换的原理并将其应用在隧道含水裂隙识别处理中，但是对含水裂隙通道的具体尺寸大小难以断定；朱自强等[10]研究了Curvelet变换处理隧道裂隙水GPR数据，但是实际资料处理比较困难，对于强干扰下目标的识别还是不够清晰与准确直观。

Stockwell等[11]提出了 S 变换时频分析方法，该方法的时频窗口具有可调性质，频率的倒数决定了窗口的尺度大小，适合用于非平稳信号的分析与处理。S 变换集中了短时傅里叶变换和小波变换的优点[12]，其时窗宽度随频率呈反向变化，即在低频段的时窗较宽, 从而获得较高的频率分辨率；而高频段的时窗较窄，故可获得很高的时间分辨[13]。S变换方法克服了其他滤波方法的缺陷，和小波变换和短时傅立叶变换相比，S变换具有较好的无损可逆性和时频分辨率[14-15]，在地震和地质雷达上具有广泛的研究价值，孙长苹[16]将S 变换用于地震信号薄层预测处理；张先武等[17]研究了将广义S变换应用在地质雷达的层位识别中；黄航如[18]对S变换在探地雷达数据处理中的应用进行了研究。

笔者提出将F-K偏移[19]与S变换结合运用于密集钢筋干扰下的地质雷达衬砌检测数据处理中，更加清晰地识别出了深部钢筋位置和异常信息，准确地判断出缺陷的范围。

1 F-K偏移成像与S变换原理

1.1 F-K偏移成像原理

设雷达记录剖面为e(x,z=0,t),其二维傅立叶变换为

E(kx,z=0,ω)

(1)

其中：

E(kx,z=0,ω)=∬E(x,z=0,t)·e-i(kx+ωt)dxdt

在频率-波数域作波场外推，得到深度Z处的波场为式(2)。

E(kx,z,ω)=E(kx,z=0,ω)e-ikzz

(2)

由于地质雷达的电磁波传播是双程走时，而波场偏移成像是基于单程传播，所以电磁波的传播速度为实际传播速度的一半。通过对E(kx,z,ω)的傅立叶反变换得到式(3)。

(3)

(4)

1.2 S变换原理

S变换是类似于Gabor变换[20]、短时傅立叶变换、连续小波变换[21]的一种方法，它的推导方法有很多种，可以由短时傅立叶变换和小波变换退推导而来，Stockwell使用的是由小波变换推导而来，其过程如下：

信号h(t)的小波变换定义为式(5)。

(5)

将S变换定义为特殊基本小波与相位因子的乘积：

S(τ,f)=ei2πfτW(τ,d)

(6)

而基本小波定义为式(7)。

(7)

其中扩张系数d是频率f的倒数。

基本小波式(7)不是严格意义上的小波，明确写出S变换的公式为式(8)。

(8)

(9)

2 数值模拟及S变换处理

2.1 模型一

利用有限差分方法模拟了地质雷达探测隧道衬砌结构的剖面，模型为双层钢筋，浅部钢筋密度大，模拟密集钢筋网，深部钢排列间距稍大。天线中心频率为1 GHz，背景介质的介电常数为6，电导率为0.01 S/m。钢筋的参数比较大，选择适中的参数为介电常数为100，电导率为1 S/m。模型大小为长2 m,深1 m，在距模型深0.25 m处有9根间距为0.2 m、半径为0.02 m的钢筋；在深0.55 m处有5根间距0.4 m、半径0.03 m的钢筋。模型示意图和雷达剖面图见图1和图2。

图1 模型一Fig.1 Model 1

图2 模型一雷达原始剖面图Fig.2 Radar original profile model 1

图3 模型一偏移处理后的剖面图Fig.3 Profile after migrationto model 1

由图2可知，浅部密集钢筋网的双曲线反射强烈，并影响了深部钢筋的反射，使得深部钢筋的双曲线反射弱，从原始图中难以识别出钢筋的数量与位置信息，需要进一步的处理。

将图2中的雷达原始图进行F-K偏移处理，得到如图3所示的剖面图，从图3中可以看出，剖面图中的双曲线反射得到了很好地校正，消除了多次波和绕射波的干扰。浅部钢筋的数量和位置都能清晰地辨别出，但是深部的钢筋信息还是很难确定。

不同介质有不同的结构特征，内部反射波的高、低频率特征明显不同，这可以作为区分不同物质界面的依据，根据反射回波的频谱信息进行识别目标，S变换能同时反映信号的频率与时间信息，这是它作为信号处理技术的优势。我们提取剖面中水平位置0.2 m处的信号(图4)，在旅行时5 ns和10 ns处具有强振幅且此两处的频率不同，根据这个特点，我们对此处信号进行S转换，得到该处的时频图(图5)。

图4 图3中水平位置0.2 m处的wiggle图Fig.4 The wiggle diagram of Figure 3 in the horizontal position 0.2 m

图5 图3中水平位置0.2 m处信号的S变换时频图Fig.5 S transform time-frequency map of Figure 3 in the horizontal position 0.2 m

图6 图3中水平位置0.4 m处信号的S变换时频图Fig.6 S transform time-frequency map of Figure 3 in the horizontal position 0.4 m

根据图5，在旅行时5 ns和10 ns处有频率为1 GHz的异常信号，由此可分辨出两个不同频率的信号，分析是浅部钢筋网和深部钢筋的信号。对比模型，我们根据时频图分辨出了浅部密集钢筋和深部钢筋的位置信息，与实际位置有很好地吻合。对比在深部没有钢筋的位置，比如在水平位置0.4 m处的时频图(图6)。由图6得知，此处深部不存在钢筋，所以在10 ns处并不见其他频率的信号，符合实际。为了更进一部说明整个剖面的情况，我们对剖面进行S变换频率切片处理(图7)。从图7中，我们能更加清晰地识别出剖面中的钢筋信息。对这个模型的处理说明了，S变换能解决密集钢筋网干扰下，深部钢筋的识别问题。

图7 图3中S变换频率切片图Fig.7 S transform frequency slices map of Figure 3

图8 模型二Fig.8 Model 2

2.2 模型二

本次模拟地质雷达对衬砌后的空洞缺陷的检测，天线装置和介质参数与模型一相同，只是在深0.70 m～-1.00 m范围内设置了一个岩石层，且初支与岩石层的接触处存在不平整面，在岩石层中设置了一个空洞，模型示意图和剖面图见图7和图8。

图9 模型二雷达原始剖面图Fig.9 Rdar original profile model 2

图10 模型二偏移处理后的剖面图Fig.10 Model 2 profile after migration

图11 图10中水平位置0.8 m处信号的S变换时频图Fig.11 S transform time-frequency map of Figure10 in the horizontal position 0.8 m

图8中的原始雷达剖面图受多次波和绕射波干扰严重，深部钢筋和空洞的信号基本在干扰下难以识别，我们对其进行偏移归位处理，得到偏移后的剖面图(图9)。图9中剖面图能得到浅部钢筋的信息，但是仍然受强浅部信号和不平整接触面的影响,导致深部的钢筋和空洞信号很弱，对于深部的信号难以识别，所以需要进行处理，我们选择对相应位置的信号进行S变换。水平位置0.8 m、0.9 m、1.0 m三处的S变换如图10～图12所示。

图12 图10中水平位置0.9 m处信号的S变换时频图Fig.12 S transform time-frequency map of Figure 10 in the horizontal position 0.9 m

图13 图10中水平位置1.0 m处信号的S变换时频图Fig.13 S transform time-frequency map of Figure 10 in the horizontal position 1.0 m

图10显示了在5 ns处的有异常区域，且频率大小范围为0.5 GHz～1.2 GHz，可知此处只是浅部钢筋网的异常，结合图8可见时频图划定的钢筋异常与实际十分符合；图11显示了在5 ns和18 ns两处有异常区域，且频率大小范围为1 GHz左右，此处是浅部钢筋网和深部空洞的异常，结合图8可见时频图划定的异常与实际符合；图12显示了在5 ns、13 ns、19 ns三处有异常区域，且频率大小范围为0.5 GHz～1.5 GHz，可知此处是浅部钢筋网、深部钢筋、深部空洞的异常，结合图8可见时频图划定的异常与实际吻合。为了更进一部说明整个剖面的情况，我们对剖面进行S变换频率切片处理，见图14。从图14中，我们能更加清晰地识别出剖面中的钢筋和空洞信息。

图14 图10中S变换频率切片图Fig.14 S transform frequency slices map of Figure10

图16 实测数据雷达偏移处理后的剖面图Fig.16 Actual radar profile after migration

2.3 实例

为了验证S变换在实际资料处理中的效果，对实际含两层钢筋的混凝土进行探测，工作雷达天线频率是1 GHz，图15是探测到的雷达原始剖面图，图16是偏移处理后的剖面图。

由图16仍然很难识别出深部钢筋的信息，于是我们将其进行S变换频率切片处理，得到图17的切片图。从图17中，我们能更加清晰地识别出浅部和深部的钢筋信号。

图17 实际数据S变换频率切片图Fig.17 S transform frequency slices map of actual data

在S变换图中，通过对每个位置的数据道进行时频分析，我们能清楚地得到该位置上的异常变化，通过时频图中的时间信息我们可以得到异常的深部位置及深部延伸，通过频率信息，我们可以得到异常的特性。上述两个模型及实例，说明了S变换在解决此类干扰下识别异常的可行性。

3 结论

笔者利用S变换对模拟隧道衬砌检测数据和实例进行了处理，表明了S变换在异常识别上的优点，能在表层密集钢筋网的强干扰下识别出深部的异常信号。通过模型模拟与实例研究得出结论：利用不同异常体反射系数的频谱特征，S 变换可以在得到反射回波的时间和频率的信息，据此可以有效地提高信号的时频分辨率，而且还可以根据实际信号具体情况，选择合适的基本小波进行计算，灵活度高、实用性强。基于 S 变换的单频曲线或单频剖面，可以真实、细致地刻画出深部异常的分布，使其具有很高的分辨率，为雷达资料高分辨处理研究工作，为下一步雷达资料处理和解释工作提供借鉴。

[1] 李大心.探地雷达方法与应用[M].北京:地质出版社,1994.

LI D X. The methods and application of GPR[M].Beijing: Geological Publishing House,1994. (In Chinese)

[2] DAVID J DANIELS. Ground-penetrating radar. The Institution of Electrical Engineers[M]. London, United Kingdom,2004.

[3] 戴前伟,冯德山.时域有限差分法在地质雷达二维正演模拟中的应用[J]．地球物理学进展,2004,19(4):898-902.

DAI Q W, FENG D S. The apply of finite difference time domain method in ground penetrating radar two-dimension forward simulation[J]. Progress in geophysics, 2004,19(4):898-902. (In Chinese)

[4] 刘四新,曾昭发.频散介质中地质雷达波传播的数值模拟[J].地球物理学报,2007,50(1):320-326.

LIU S X, ZENG Z F. Numerical simulation for ground penetrating radar wave propagation in the dispersive medium[J].Geophys,2007,50(1):320-326. (In Chinese)

[5] 肖都,李文杰,郭鹏.基于GPRMax的隧道衬砌检测数值模拟及应用[J].物探与化探,2015,39(4):855-859.

XIAO D, LI W J, GUO P. The tunnel lining detection numerical simulation based on GPRMax and its application[J].Geophysical and Geochemical Exploration,2015,39(4):855-859. (In Chinese)

[6] 赵峰,周斌,武永胜.探地雷达在隧道衬砌空洞检测的正演模拟应用研究[J].铁道建筑,2012(8):99-103.

ZHAO F, ZHOU B, WU Y S.Study on application of forward modelling for inspecting void in tunnel lining by ground probing radar[j].Railway Engineering,2012(8):99-103. (In Chinese)

[7] 舒志乐,刘新荣.隧道衬砌病害探地雷达三维正演模拟及工程验证[J].中国铁道科学，2013,34(4):46-53.

SHU Z L, LIU S X.GPR three-dimensional forward modeling of defects in tunnel lining and engineering verification[J].China Railway Science,2013,34(4):46-53. (In Chinese)

[8] 黄玲,曾昭发.钢筋混凝土缺陷的探地雷达检测模拟与成像效果[J].物探与化探，2007,31(2):181-185.

HUANG L, ZENG Z F .Modeling and imaging of GPR detection of reinforced concrete defects[J]. Geophysical and Geochemical Exploration,2007,31(2):181-185. (In Chinese)

[9] 张家松,朱自强,鲁光银,等.偏移成像与希尔伯特变换在地质雷达探测隧道裂隙流水通道中的应用[J].地球物理学进展,2012,27(6):2722-2728.

ZHANG J S, ZHU Z Q, LU G Y, et al. Application of migration and Hilbert transform on exploration of karst fractured ground water channel in the tunnel with GPR[J]. Progress in geophys,2012,27(6):2722-2728. (In Chinese)

[10] 朱自强,朱贺,鲁光银,等.基于Curvelet变换的隧道裂隙水GPR数据处理研究[J].物探化探计算技术,2014,36(5):571-576.

ZHU Z Q, ZHU H, LU G Y, et al. Processing of GPR data in tunnel fissure water based on Curvelet transform[J].Computing Techniques for Geophysical and Geochemical Exploration,2014,36(5):571-576. (In Chinese)

[11] STOCKWELL R G, MANSINHA L, LOWE R P.Localization of the Complex Spectrum: TheSTransform[J].IEEE Transactions on signal Processing,1996,44(4):998-1001.

[12] 赵淑红,朱光明.S变换时频滤波去噪方法[J].石油地球物理勘探,2007,42(4):402-406.

ZHAO S H, ZHU G M. Time-frequency filtering to denoise by S transform[J]. Oil Geophysical Prospecting,2007,42(4):402-406. (In Chinese)

[13] 操峰,陈淑珍,肖博勋.一种改进的STFT在PGPR回波信号处理中的应用[J].计算机应用研究,2005,22(6):127-129.

CAO F, CHEN S F, XIAO B X. Analysis of echo signal by an improved STFT algorithm for phased-array ground penetrating radar[J].Application Research of Computers, 2005,22(6):127-129. (In Chinese)

[14] 王晓凯,高静怀,陈文超,等.基于高维连续小波变换的地震资料不连续性检测方法研究[J].地球物理学报,2016,59(9):3394-3407.

WANG X K, GAO J H, CHEN W C, et al.Detecting method of seismic discontinuities based on high dimensional continuous wavelet transform[J].Geophys,2016,59(9):3394-3407. (In Chinese)

[15] 孙月.基于广义S 变换和阈值函数的地震信号去噪研究[D].吉林：吉林大学,2012:13-26.

SUN Y. Seismic signal denoising research based on generalized S transformation and threshold function[D].Jilin:Jilin University, 2012:12-26. (In Chinese)

[16] 孙长苹.S变换地震信号薄层预测处理技术[D].大庆：东北石油大学,2011:36-52.

SUN C P. Thin bed prediction processing technique of seisemic signal based on S transform[D].Daqing: Northeast Petroleum University,2011:36-52. (In Chinese)

[17] 张先武,高云泽,方广有.带有低通滤波的广义S变换在探地雷达层位识别中的应用[J].地球物理学报,2013,56(1):309-316.

ZHANG X W, GAO Y Z, FANG G Y. Application of generalized S transform with lowpass filtering to layer recognition of ground penetrating radar[J]. Geophys,2013,56(1):309-316. (In Chinese)

[18] 黄航如.S变换在探地雷达数据处理中的应用研究[D].广西:桂林理工大学,2012:24-37.

HUANG H R. The study of the application of S transform on ground penetrating radar data progressing[D].Guangxi: Guilin University of Technology,2012:24-37. (In Chinese)

[19] 李学聪,刘伊克,常旭,等.基于F-K偏移和反偏移的地震道插值方法研究[J].地球物理学进展,2009,24(6):2020-2029.

LI X C, LIU Y K, CHANG X, et al. Interpolation for seismic trances based on F-K migration and demigration[J].Progress Geophys,2009,24(6):2020-2029. (In Chinese)

[20] GABOR D. Theory of communication[J].Journal ofthe Institution of Electrical Engineers,1946,93(3):429-457.

[21] GOUPILLAUD P,GROSSMANN A,MORLET J.Cycle-octave and related transforms in seismic analysis[J].Geoexploration,1984,23(1):85-102.

ApplicationofStransformindataprocessingofdetectionoftunnelliningbygroundpenetratingradar

ZHU Ziqiang, WANG Yuqin, LU Guangyin

(School of Geosciences and Info-Physics,Central South University,Changsha 410083,China)

In this paper, the S transform is applied to the data processing of ground penetrating radar to detect the tunnel lining. To solve the problem that the radar data profile is seriously disturbed and the deep signal is weak and the disease abnormity is difficult to be identified due to the influence of the dense steel mesh. We carried out two-dimensional finite-difference numerical simulation of ground penetrating radar, studied the influence of shallow dense steel mesh on the deep steel arch and the void space, used the migration imaging method to deal with the original profile, repression of the diffraction wave and multiple waves, improve the deep weak reflection signal imaging results; Finally, the characteristics of the deep rebar and the void are identified clearly and accurately by the S transform time-frequency analysis. By using the S transform time-frequency processing of the simulation data and the actual data. This paper provides a effective solution to deal with the data of tunnel lining which is seriously disturbed.

S transform; ground penetrating radar; migration imaging; detection of tunnel lining

2016-11-03 改回日期： 2017-01-04

朱自强(1964-)，男，教授，博士生导师，从事地质灾害探测与监测工作，E-mail：13507319431@139.com。

王玉琴(1991-)，男，硕士，主要从事地质雷达信号处理研究，E-mail：940603995@qq.com。

1001-1749(2017)06-0755-07

P 631.4

A

10.3969/j.issn.1001-1749.2017.06.07