CN111981968B - 基于地质雷达小波变换的隧道衬砌厚度识别方法及系统 - Google Patents
基于地质雷达小波变换的隧道衬砌厚度识别方法及系统 Download PDFInfo
- Publication number
- CN111981968B CN111981968B CN202010908318.5A CN202010908318A CN111981968B CN 111981968 B CN111981968 B CN 111981968B CN 202010908318 A CN202010908318 A CN 202010908318A CN 111981968 B CN111981968 B CN 111981968B
- Authority
- CN
- China
- Prior art keywords
- wavelet
- geological radar
- radar
- tunnel lining
- geological
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 105
- 239000004567 concrete Substances 0.000 claims abstract description 84
- 230000004044 response Effects 0.000 claims abstract description 49
- 238000012545 processing Methods 0.000 claims abstract description 47
- 238000004458 analytical method Methods 0.000 claims abstract description 28
- 230000009466 transformation Effects 0.000 claims abstract description 28
- 238000001514 detection method Methods 0.000 claims abstract description 23
- 238000013461 design Methods 0.000 claims abstract description 16
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 14
- 238000013508 migration Methods 0.000 claims abstract description 13
- 230000005012 migration Effects 0.000 claims abstract description 13
- 230000006870 function Effects 0.000 claims description 142
- 238000004088 simulation Methods 0.000 claims description 34
- 239000011378 shotcrete Substances 0.000 claims description 22
- 239000011435 rock Substances 0.000 claims description 21
- 229910000831 Steel Inorganic materials 0.000 claims description 18
- 239000010959 steel Substances 0.000 claims description 18
- 238000010276 construction Methods 0.000 claims description 15
- 230000008569 process Effects 0.000 claims description 13
- 238000005070 sampling Methods 0.000 claims description 13
- 239000010410 layer Substances 0.000 claims description 11
- 230000000694 effects Effects 0.000 claims description 6
- 239000011241 protective layer Substances 0.000 claims description 6
- 238000013519 translation Methods 0.000 claims description 6
- 238000004364 calculation method Methods 0.000 claims description 5
- 238000010521 absorption reaction Methods 0.000 claims description 4
- 239000006185 dispersion Substances 0.000 claims description 4
- 238000004134 energy conservation Methods 0.000 claims description 4
- 230000005284 excitation Effects 0.000 claims description 4
- 239000011159 matrix material Substances 0.000 claims description 4
- 239000002689 soil Substances 0.000 claims description 4
- 230000002708 enhancing effect Effects 0.000 claims 1
- 229910001294 Reinforcing steel Inorganic materials 0.000 description 11
- 239000011150 reinforced concrete Substances 0.000 description 8
- 238000010586 diagram Methods 0.000 description 7
- 238000005259 measurement Methods 0.000 description 5
- 238000009412 basement excavation Methods 0.000 description 3
- 230000007547 defect Effects 0.000 description 3
- 238000005553 drilling Methods 0.000 description 3
- 238000000605 extraction Methods 0.000 description 3
- 230000005428 wave function Effects 0.000 description 3
- 239000002131 composite material Substances 0.000 description 2
- 238000001914 filtration Methods 0.000 description 2
- 230000010354 integration Effects 0.000 description 2
- 230000004807 localization Effects 0.000 description 2
- 238000009659 non-destructive testing Methods 0.000 description 2
- WHXSMMKQMYFTQS-UHFFFAOYSA-N Lithium Chemical compound [Li] WHXSMMKQMYFTQS-UHFFFAOYSA-N 0.000 description 1
- 230000003321 amplification Effects 0.000 description 1
- 230000002457 bidirectional effect Effects 0.000 description 1
- 238000005422 blasting Methods 0.000 description 1
- 239000011083 cement mortar Substances 0.000 description 1
- 230000007797 corrosion Effects 0.000 description 1
- 238000005260 corrosion Methods 0.000 description 1
- 230000001066 destructive effect Effects 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000009440 infrastructure construction Methods 0.000 description 1
- 229910052744 lithium Inorganic materials 0.000 description 1
- 238000012423 maintenance Methods 0.000 description 1
- 238000011430 maximum method Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000003199 nucleic acid amplification method Methods 0.000 description 1
- 238000009527 percussion Methods 0.000 description 1
- 230000008092 positive effect Effects 0.000 description 1
- 238000005086 pumping Methods 0.000 description 1
- 238000004445 quantitative analysis Methods 0.000 description 1
- 230000003014 reinforcing effect Effects 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 238000010183 spectrum analysis Methods 0.000 description 1
- 239000000126 substance Substances 0.000 description 1
- 238000012795 verification Methods 0.000 description 1
- 230000000007 visual effect Effects 0.000 description 1
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01B—MEASURING LENGTH, THICKNESS OR SIMILAR LINEAR DIMENSIONS; MEASURING ANGLES; MEASURING AREAS; MEASURING IRREGULARITIES OF SURFACES OR CONTOURS
- G01B7/00—Measuring arrangements characterised by the use of electric or magnetic techniques
- G01B7/02—Measuring arrangements characterised by the use of electric or magnetic techniques for measuring length, width or thickness
- G01B7/06—Measuring arrangements characterised by the use of electric or magnetic techniques for measuring length, width or thickness for measuring thickness
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S7/00—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
- G01S7/02—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00
- G01S7/41—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00 using analysis of echo signal for target characterisation; Target signature; Target cross-section
- G01S7/411—Identification of targets based on measurements of radar reflectivity
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Radar, Positioning & Navigation (AREA)
- Remote Sensing (AREA)
- Computer Networks & Wireless Communication (AREA)
- Radar Systems Or Details Thereof (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
本发明属于隧道衬砌厚度检测技术领域，公开了一种基于地质雷达小波变换的隧道衬砌厚度识别方法及系统。根据电磁波发生反射与折射均在二维平面内的特性，建立隧道衬砌不同厚度混凝土结构地电模型。利用地质雷达天线对隧道衬砌不同厚度混凝土结构地电模型进行正演模拟；根据偏移算法，对地质雷达特征响应进行偏移处理；根据地质雷达发射子波的时频特征，采用高斯函数拟合构造地质雷达小波基。根据已构造的地质雷达小波基，编写小波变换时能密度分析程序，对地质雷达正演模拟信号和实测信号进行小波变换时能密度法分析识别，获得隧道衬砌混凝土结构的厚度值。本发明分析的地质雷达信号分辨率高，结果准确，识别误差小，小于设计值的5％。
Description
技术领域
本发明属于隧道衬砌厚度检测技术领域，尤其涉及一种基于地质雷达小波变换的隧道衬砌厚度识别方法及系统。
背景技术
目前，隧道工程是交通基础设施建设的关键通道、节点以及重点控制性工程之一。截至2018年底，包含铁路隧道、公路隧道以及城市轨道交通隧道等在内的中国隧道建设共计36103km，正在施工的各类隧道建设约20000km，计划修建的隧道建设约20000km。在隧道工程的建造过程中，由于隧道围岩光面爆破效果不佳、喷射混凝土施工工艺不规范以及泵送模筑混凝土压力不足等原因，造成围岩与衬砌结构之间以及喷射混凝土、防水板与模筑混凝土之间常存在空洞与脱空等情况的发生，从而导致隧道衬砌混凝土结构厚度的不足。当衬砌混凝土结构厚度不足时，容易产生应力集中，使得结构受力不均匀，导致隧道衬砌混凝土结构沉降，甚至会造成衬砌结构发生断裂、坍塌等事故，进而影响运营期隧道工程的正常使用。钢筋混凝土保护层厚度是隧道钢筋衬砌混凝土结构耐久性与抗腐蚀性的另一个重要影响因素。若钢筋保护层实际厚度不能满足设计要求，那么钢筋在受到有害物质的侵蚀后，会严重锈蚀钢筋衬砌混凝土结构，进而威胁整个隧道衬砌结构的安全与稳定。由此可知，对隧道衬砌混凝土结构厚度与钢筋保护层厚度进行安全检测显得尤为重要。
现有常用于隧道衬砌结构厚度检测的方法有：预埋钢筋法、凿孔法、取芯法、断面仪法和地质雷达法等。预埋钢筋法是在喷射混凝土或者围岩预定位置上，采用冲击钻钻孔，然后清孔，用水泥砂浆埋置好钢筋或膨胀螺杆。这种施工方法虽然技术成熟与施工方便，但没有考虑隧道衬砌结构的成层性，即复合式衬砌结构由初期支护和二次衬砌组成，还包括中间的防水板。预埋钢筋或膨胀螺杆法会破坏防水板，造成衬砌结构渗漏水。凿孔法与取芯法都是有损检测方法，这种方法比较直观，但难免对隧道衬砌结构整体性造成损伤，甚至使隧道衬砌的防水系统遭到破坏，造成严重后果。并且该种检测方法速度慢，只能进行少量抽样检测，不能反映隧道衬砌混凝土结构的整体质量。隧道激光断面仪法需要在同一桩号前后量测几次，通过隧道衬砌断面轮廓相减法得到衬砌结构厚度，这样容易增加检测人员的工作量，并且隧道激光断面仪前后架设在同一位置的误差也较大。与以上方法相比，地质雷达因具有探测精度高、无损检测、结果直观与操作灵活等优点，已被广泛应用于隧道衬砌结构厚度无损检测中。地质雷达在现场检测过程中，采集的原始雷达信号常被复杂的现场环境所干扰，导致原始雷达信号需采用先进的信号分析方法处理后才能解译。
近年来，随着信号分析处理领域的技术突破，小波时频分析方法在地质雷达信号阈值去噪、特征提取、频谱分析、波形预测和图像识别等方面具有了良好的解译能力。但是在用地质雷达剖面法识别混凝土，尤其是钢筋混凝土隧道衬砌厚度时，由于钢筋衬砌混凝土结构中的钢筋会对地质雷达发射的电磁波产生强烈干扰，导致电磁波能量衰减快，使得深部二次衬砌混凝土与喷射混凝土、喷射混凝土与围岩之间的反射界面处电磁波反射信号弱，难以获得钢筋衬砌混凝土结构的厚度。特别是当模筑混凝土衬砌的相对介电常数小于喷射混凝土衬砌的相对介电常数、喷射混凝土衬砌的相对介电常数大于围岩的相对介电常数时，电磁波会在各界面处发生相位反相。因此，从地质雷达时间剖面上无法识别电磁波遇到界面产生反射波的瞬时位置(究竟处于波峰还是波谷)。
通过上述分析，现有技术存在的问题及缺陷为：
(1)预埋钢筋法、凿孔法与取芯法都会对隧道衬砌结构整体性(尤其是隧道防水系统)造成破坏，而且这几种检测方法速度慢，仅能对选取的少数断面进行抽样检测，不能完全保证隧道衬砌混凝土结构的施工质量。隧道激光断面仪法很难保证架设在同一个位置，前后的测量误差也较大，不仅增加检测人员的工作量，也会在隧道有限的作业空间内影响施工进度。
(2)地质雷达法在识别钢筋混凝土衬砌厚度时，隧道衬砌钢筋混凝土中的钢筋对地质雷达电磁波检测产生强烈干扰，导致电磁波能量衰减快，使得深部二次衬砌混凝土与喷射混凝土、喷射混凝土与围岩之间的反射界面处电磁波反射信号弱，采用地质雷达剖面法难以获得钢筋衬砌混凝土结构的厚度。
(3)当模筑混凝土衬砌的相对介电常数小于喷射混凝土衬砌的相对介电常数、喷射混凝土衬砌的相对介电常数大于围岩的相对介电常数时，电磁波会在各界面处发生相位反相，而从地质雷达时间剖面上无法识别电磁波遇到界面产生反射波的瞬时位置。
解决以上问题及缺陷的难度为：仅通过地质雷达剖面上的波峰以及波幅变化难以准确识别衬砌混凝土与喷射混凝土、喷射混凝土与围岩之间的界面位置，并且由于钢筋对电磁波的强干扰作用，使得深部界面处的弱信号处理的难度大大增加。同时，为精确获取电磁波在围岩等介质内的传播距离，需要预知电磁波在不同界面间相对介电常数的差异，这样就难免因主观因素而造成误判。
解决以上问题及缺陷的意义为：利用小波变换强大的时频分析以及精细化识别特征点信号的能力，不仅能解决电磁波在各界面相位反相无法准确识别的问题，而且还能抑制钢筋对电磁波信号的强干扰，对深部界面弱信号的提取与识别也有增强作用。本发明能有效提高地质雷达检测的识别精度，为隧道衬砌结构厚度的高效且精准识别提供了理论依据与可行性方案。
发明内容
针对现有技术存在的问题，本发明提供了一种基于地质雷达小波变换的隧道衬砌厚度识别方法及系统。
本发明是这样实现的，一种基于地质雷达小波变换的隧道衬砌厚度识别方法，所述基于地质雷达小波变换的隧道衬砌厚度识别方法，包括：
步骤一，假设隧道围岩为半无限空间连续的均匀介质，依据电磁波反射与折射均在二维平面内的特性，建立隧道衬砌混凝土结构地电模型，对隧道衬砌结构不同厚度在地质雷达正演模拟信号上的反射特征进行分析；
步骤二，隧道衬砌不同厚度混凝土结构的雷达响应，利用地质雷达天线对不同厚度衬砌结构的地电模型进行正演模拟，得到隧道衬砌混凝土结构正演模拟的雷达响应图像；
步骤三，根据电磁波在界面处的反射系数为负的原则，找到地质雷达信号最大的负峰值点位置，对地质雷达响应进行去除直达波处理；
步骤四，根据偏移算法，在MATLAB语言工作平台上，对隧道衬砌不同厚度混凝土结构雷达响应进行偏移处理；
步骤五，地质雷达信号的规范处理，采用去除背景信号与直流分量以及衰减补偿方法对地质雷达响应图像进行处理；
步骤六，地质雷达小波基函数的构造与添加，根据地质雷达发射子波的衰减特性，采用高斯函数曲线拟合构造地质雷达小波基，验证小波基函数满足小波容许性条件。在MATLAB语言工作平台环境下，将构造的地质雷达小波基函数添加到小波分析工具箱，以供地质雷达小波变换时选用。
步骤七，小波变换时能密度法的定量识别，根据已构造的地质雷达小波基函数，提出一种基于地质雷达小波变换的时能密度算法，编写小波变换时能密度分析程序，采用该程序对地质雷达单道信号进行小波变换时能密度法分析识别。
进一步，所述步骤一中，地电模型参数设置为：模拟区域范围7.2m×2.0m，左上角为坐标原点，水平横坐标为距离，竖向纵坐标为探测深度；
将水平距离长为7.2m划分成六等份，每份长为1.2m，厚度值从左至右依次分别为0.40m、0.35m、0.30m、0.30m、0.35m和0.40m；
钢筋布置型号均为Ф20，每份钢筋间距均为0.24m，保护层厚度为50mm；模筑混凝土与喷射混凝土衬砌的相对介电常数分别为6和20，电导率分别为0.01S/m和0.1S/m，导磁率均为1。
进一步，所述步骤二中，根据隧道衬砌结构现场检测使用的地质雷达参数，模拟所采用的雷达天线中心频率为900MHz，边界吸收条件为完全匹配层，激励源采用Ricker子波，900MHz雷达天线的采样时窗为15ns，网格的空间步长均为0.005m，采样步长为0.01m，采样道数为701道。
进一步，所述步骤三中，地质雷达发射的脉冲电磁波从空气进入岩土层中，由于两种不同介质相对介电常数的差异，将在界面上产生反射现象，其规律符合反射定律，产生的电磁波能量大小取决于反射系数
根据电磁波在界面处的反射系数为负的原则，找到地质雷达信号最大的负峰值点位置，对地质雷达响应进行去除直达波处理。
进一步，所述步骤四中，地质雷达信号偏移处理过程为：
设
为矩阵E(x,z＝0,t)的二维傅立叶变换，则
在频率波数域内，深度z处的波场可以表示为：
设E(x,z,t)(E(x,t)在z方向的延拓)为
关于kx,w的二维傅立叶逆变换
由式(3)，可得
令t＝0，得到
然后，根据频散关系式
通过式(6)，可将
也就是
插值映射为
得到下式：
根据上述偏移算法，在MATLAB语言工作平台上，对雷达响应进行偏移处理。
进一步，所述步骤六中，地质雷达小波基函数的构造过程为：
在MATLAB语言工作平台上，选用Curve Fitting工具箱中带参数的高斯函数对地质雷达发射子波进行曲线拟合，然后比较各函数模型的拟合效果；
通过拟合精度对不同参数获得的高斯函数进行综合比选，得到的高斯函数形式为：
要想知道构造的拟合高斯函数是否可以作为小波基，需验证高斯函数是否满足小波基函数容许性的条件。也就是说，满足小波基的波动性与衰减性的要求。首先，对拟合高斯函数进行积分，可得
由于高斯函数积分获得的结果，相对于拟合高斯函数f(x)的幅值来说，量值较小，可忽略。由此，可判断拟合高斯函数f(x)满足小波基的波动性要求∫Rf(t)dt≈0；然后，通过对拟合高斯函数f(x)进行积分计算，可知
(其中α＝1)，即拟合高斯函数f(x)满足小波基的衰减性要求；
由于拟合高斯函数f(x)满足小波基函数容许性的条件，就可以将拟合高斯函数f(x)添加到小波分析工具箱中，并作为一个新的小波基函数进行选用。
进一步，所述地质雷达小波基函数的添加，在小波基函数添加之前，进行如下工作：①定义小波基函数的全称(familyName)；②定义小波基函数的缩写名(familyShortName)；③定义小波基函数的类型(familyWaveType)；④新建一个MATLAB中的mat文件(fileWaveName)；
在MATLAB语言工作平台上，可使用wavemngr命令添加一个新的小波函数，具体形式如下：wavemngr(‘add’，fn，fsn，wt，file)，其中fn为小波基函数的全称；fsn为小波基函数的缩写名；wt为小波基函数的类型；file为mat文件的文件名。
进一步，所述步骤七中，小波变换时能密度法，假设地质雷达小波基函数满足容许性条件，则地质雷达小波变换是完备的并且保持能量守恒；根据Moyal定理，有下式成立
式中表明地质雷达小波变换系数平方的积分同被分析原始信号的能量成正比；因此，把|Wf(a,b)|2/Cψa2看作是地质雷达小波时间-尺度平面上的能量密度函数，即|Wf(a,b)|2ΔaΔb/Cψa2给出了以尺度因子a和平移因子b为中心的、尺度间隔为Δa、时间间隔为Δb的能量；根据小波能量密度的概念，也可以写成
∫R|f(t)|2dt＝∫RE(b)db (10)
式中：
在小波变换中，尺度因子a在某种程度上对应于实际频率f；因此，给出了地质雷达信号所有频带的能量随平移因子b的分布情况，称为小波变换时能密度函数。
进一步，所述步骤七中，通过识别小波变换时能密度曲线上的双程旅行时间Δt＝t2-t1，利用距离公式
可以得到隧道衬砌混凝土结构的厚度值L；其中：c为电磁波在真空中的传播速度(3.0×108m/ns)，t1为直达波的反射位置(t1＝0)，t2为小波变换时能密度曲线上峰值点的位置，ε为介质的相对介电常数；
采用小波变换时能密度法对隧道衬砌不同厚度的混凝土结构进行定量识别，得到隧道衬砌不同厚度混凝土结构的识别结果，并与理论设计值进行比较，结果表明采用小波变换时能密度法识别得到的相对误差值均不超过设计值的5％。
本发明另一目的在于提供一种实施所述的基于地质雷达小波变换的隧道衬砌厚度识别方法及系统，所述基于地质雷达小波变换的隧道衬砌厚度识别系统包括：
模型构建模块，假设隧道围岩为半无限空间连续的均匀介质，电磁波反射与折射均在二维平面内进行，建立隧道衬砌不同厚度的混凝土结构地电模型；
正演模拟模块，利用雷达天线对不同厚度混凝土结构地电模型进行正演模拟，得到隧道衬砌混凝土结构的正演模拟雷达响应图像；
直达波去除处理模块，根据电磁波在界面处的反射系数为负的原则，找到地质雷达信号最大的负峰值点位置，对地质雷达响应图像进行去除直达波处理；
雷达信号偏移处理模块，根据偏移算法，在MATLAB语言工作平台上，对雷达响应图像进行偏移处理；
雷达信号规范处理模块，采用去除背景信号与直流分量以及衰减补偿方法对地质雷达响应图像进行信号去噪与弱信号增强处理；
地质雷达小波基的构造与添加，根据地质雷达子波在时频域内的衰减特性，采用高斯函数构造地质雷达小波基，并将其添加到小波分析工具箱，以供地质雷达小波变换时选用。
雷达信号定量识别模块，根据已构造的地质雷达小波基，在MATLAB语言工作平台环境下，编写小波变换时能密度分析程序，采用该程序对地质雷达单道信号进行小波变换时能密度法分析识别。
结合上述的所有技术方案，本发明所具备的优点及积极效果为：本发明根据电磁波发生反射与折射均在二维平面内的特性，建立隧道衬砌不同厚度混凝土结构地电模型。利用地质雷达天线对隧道衬砌不同厚度混凝土结构地电模型进行正演模拟，得到隧道衬砌混凝土结构的正演模拟雷达响应图像，对隧道衬砌结构厚度在地质雷达正演模拟信号上的反射特征进行分析，找到电磁波上最大的负峰值点位置，进行去除直达波处理；根据偏移算法，对地质雷达特征响应进行偏移处理；采用去除背景信号与直流分量以及衰减补偿方法对地质雷达响应图像进行处理；根据地质雷达发射子波的时频特征，采用高斯函数拟合构造地质雷达小波基，并将其添加到小波分析工具箱中，以供地质雷达小波变换时选用。根据已构造的地质雷达小波基，提出一种基于地质雷达小波变换的小波变换时能密度方法，编写小波变换时能密度分析程序，采用该程序对地质雷达正演模拟信号和实测信号进行小波变换时能密度法分析识别，可以获得隧道衬砌混凝土结构的厚度值。本发明分析的地质雷达信号分辨率高，结果准确，识别误差小，小于设计值的5％，为隧道衬砌地质雷达图像特征的解译提供参考。
(1)本发明能够有效提高地质雷达检测隧道衬砌结构的探测精度，将采用时域有限差分法对隧道衬砌不同厚度混凝土(包括素混凝土与钢筋混凝土)进行地质雷达正演模拟。采用经典Db系列小波和地质雷达小波对不同厚度的隧道衬砌混凝土结构正演模拟图像进行识别，发现Db系列不同小波基分析同一地质雷达信号，识别结果不完全相同，甚至差异很大。
(2)本发明采用地质雷达小波变换与小波变换时能密度法对隧道衬砌混凝土结构地质雷达信号进行分析，图像分辨率高，识别误差相对较小。将其与隧道衬砌混凝土结构厚度理论设计值进行比较，表明采用小波变换时能密度法识别得到的隧道衬砌结构厚度的相对误差值均不超过设计值的5％。
(3)本发明使用小波变换时能密度法对隧道混凝土衬砌结构地质雷达响应特征信号进行定量分析，大大提高钢筋强干扰下衬砌混凝土结构厚度识别的准确度与可靠度，从而为解译隧道衬砌结构的地质雷达特征响应图提供参考。
附图说明
为了更清楚地说明本申请实施例的技术方案，下面将对本申请实施例中所需要使用的附图做简单的介绍，显而易见地，下面所描述的附图仅仅是本申请的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下还可以根据这些附图获得其他的附图。
图1是本发明实施例提供的基于地质雷达小波变换的隧道衬砌厚度识别方法流程图。
图2是本发明实施例提供的基于地质雷达小波变换的隧道衬砌厚度识别系统结构示意图。
图中：1、模型构建模块；2、正演模拟模块；3、直达波去除处理模块；4、雷达信号偏移处理模块；5、雷达信号规范处理模块；6、地质雷达小波基构造与添加模块；7、雷达信号定理识别模块。
图3是本发明实施例提供的地质雷达发射子波时域图。
图4是本发明实施例提供的地质雷达发射子波功率谱图。
图5是本发明实施例提供的地质雷达发射子波曲线拟合图。
图6是本发明实施例提供的地质雷达小波基的显示图。
图7是本发明实施例提供的隧道衬砌混凝土结构的地电模型图。
图8是本发明实施例提供的隧道衬砌混凝土结构的雷达响应图。
图9是本发明实施例提供的隧道衬砌混凝土结构去除直达波后的雷达响应图。
图10是本发明实施例提供的隧道混凝土衬砌结进行F-K偏移处理后的雷达响应图。
图11是本发明实施例提供的隧道混凝土衬砌结进行规范处理后的雷达响应图。
图12是本发明实施例提供的从图11所示雷达响应抽取64条单道信号图。
图13是本发明实施例提供的1条单道信号的Db2小波变换模极大值曲线。
图14是本发明实施例提供的1条单道信号的Db4小波变换模极大值曲线。
图15是本发明实施例提供的1条单道信号的Db6小波变换模极大值曲线。
图16是本发明实施例提供的1条单道信号的Db8小波变换模极大值曲线。
图17是本发明实施例提供的1条单道信号的Db10小波变换模极大值曲线。
图18是本发明实施例提供的1条单道信号的雷达小波变换模极大值曲线。
图19是本发明实施例提供的1条单道信号的小波变换时能密度曲线图。
图20是本发明实施例提供的64条单道信号的小波变换时能密度曲线图。
图21是本发明实施例提供的采用小波变换时能密度法识别的层位线图。
图22是本发明实施例提供的实测地质雷达原始剖面图。
图23是本发明实施例提供的经规范处理后地质雷达剖面图。
图24是本发明实施例提供的从图23中抽取40条单道信号图。
图25是本发明实施例提供的采用小波变换时能密度法识别的层位线图。
具体实施方式
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。
针对现有技术存在的问题，本发明提供了一种基于地质雷达小波变换的隧道衬砌厚度识别方法及系统，下面结合附图对本发明作详细的描述。
如图1所示，本发明实施例提供的基于地质雷达小波变换的隧道衬砌厚度识别方法，包括：
S101：假设隧道围岩为半无限空间连续的均匀介质，电磁波反射与折射均在二维平面内进行，建立隧道衬砌不同厚度混凝土结构地电模型，对隧道衬砌混凝土结构厚度在地质雷达正演模拟信号上的反射特征进行分析；
S102：不同厚度隧道衬砌混凝土结构的雷达响应，利用雷达天线对隧道衬砌混凝土结构地电模型进行正演模拟，得到隧道衬砌混凝土结构的正演模拟雷达响应图像。
S103：根据电磁波在界面处的反射系数为负的原则，找到地质雷达信号最大的负峰值点位置，对地质雷达响应图像进行去除直达波处理；
S104：根据F-K偏移算法，在MATLAB语言工作平台上，对地质雷达响应图像进行偏移处理。
S105：地质雷达信号的规范处理，采用去除背景信号与直流分量以及衰减补偿方法对地质雷达响应图像进行处理。
S106：地质雷达小波基的构造与添加，根据地质雷达发射子波在时频域内的衰减特性，采用高斯函数曲线拟合构造地质雷达小波基，验证地质雷达小波基函数满足小波容许性条件。在MATLAB语言工作平台环境下，将构造的地质雷达小波基函数添加到小波分析工具箱，以供地质雷达小波变换时选用。
S107：小波变换时能密度法的定量识别，根据已构造的地质雷达小波基，提出一种基于地质雷达小波变换的时能密度法，在MATLAB语言工作平台环境下，编写小波变换时能密度分析程序，采用该程序对地质雷达单道信号进行小波变换时能密度法分析识别。
本发明提供的基于地质雷达小波变换的隧道衬砌厚度识别方法，业内的普通技术人员还可以采用其他的步骤实施，图1的本发明提供的基于地质雷达小波变换的隧道衬砌厚度识别方法仅仅是一个具体实施例而已。
本发明实施例提供的S101中，地电模型参数设置为：模拟区域范围7.2m×2.0m，左上角为坐标原点，水平横坐标为距离，竖向纵坐标为探测深度；将水平距离长为7.2m划分成六等份，每份长为1.2m，厚度值从左至右依次分别为0.40m、0.35m、0.30m、0.30m(含钢筋)、0.35m(含钢筋)和0.40m(含钢筋)；钢筋布置型号均为Ф20，每份钢筋间距均为0.24m，保护层厚度为50mm；模筑混凝土与喷射混凝土(围岩)衬砌的相对介电常数分别为6和20，电导率分别为0.01S/m和0.1S/m，导磁率均为1。
本发明实施例提供的S102中，根据隧道衬砌混凝土结构现场检测使用的地质雷达参数，模拟所采用的雷达天线中心频率为900MHz，边界吸收条件为完全匹配层，激励源采用Ricker子波，900MHz雷达天线的采样时窗为15ns，网格的空间步长均为0.005m，采样步长为0.01m，采样道数为701道。
本发明实施例提供的S103中，地质雷达发射的脉冲电磁波从空气进入岩土层中，由于两种不同介质相对介电常数的差异，将在空气与二衬、二衬与初衬界面上产生反射现象，电磁波传播规律符合反射定律，产生的电磁波能量大小取决于反射系数
根据电磁波在界面处的反射系数为负的原则，找到地质雷达信号最大的负峰值点位置，对地质雷达响应图像进行去除直达波处理。
本发明实施例提供的S104中，地质雷达信号偏移处理过程为：
设
为矩阵E(x,z＝0,t)的二维傅立叶变换，则：
在频率波数域内，深度z处的波场可以表示为：
设E(x,z,t)(E(x,t)在z方向的延拓)为
关于kx,w的二维傅立叶逆变换：
由式(3)，可得：
令t＝0，得到：
然后，根据频散关系式：
通过式(6)，可将
也就是
插值映射为
得到下式：
根据上述偏移算法，在MATLAB语言工作平台上，对雷达响应图像进行偏移处理。
本发明实施例提供的S107中，通过提取小波变换时能密度曲线上的双程旅行时间Δt＝t2-t1，利用距离公式
可以得到隧道衬砌混凝土结构的厚度值L，其中：c为电磁波在真空中的传播速度(3.0×108m/ns)，t1为直达波的反射位置(t1＝0)，t2为小波变换时能密度曲线上峰值点的位置，ε为介质的相对介电常数。
将不同厚度隧道衬砌混凝土结构正演模拟识别结果与理论设计值进行比较，采用小波变换时能密度法识别得到的衬砌混凝土结构厚度的相对误差值均不超过设计值的5％。
本发明实施例提供的S106中，地质雷达小波基函数的构造过程为：
在MATLAB语言工作平台上，选用Curve Fitting工具箱中带参数的高斯函数对地质雷达发射子波进行曲线拟合，然后比较各函数模型的拟合效果；
通过拟合精度对各函数模型进行综合比选后，得到的拟合高斯函数形式为：
要想知道构造的地质雷达子波是否可以作为小波基供小波变换时选用，需要验证拟合高斯函数是否满足小波基函数容许性的条件，即波动性与衰减性的要求。首先对拟合高斯函数进行积分，可得
由于拟合高斯函数积分后的计算结果，相对高斯函数f(x)的幅值来说，量值相对较小，可忽略不计。由此，可判断拟合高斯函数f(x)满足波动性要求，即：∫Rf(t)dt≈0；然后，通过对拟合高斯函数进行积分计算，可知
(其中α＝1)，即拟合高斯函数f(x)满足衰减性要求；
综上所述，拟合高斯函数f(x)满足小波基函数容许性的条件；再将拟合高斯函数f(x)添加到小波分析工具箱中，可以作为一个新的小波基供小波变换分析时进行选用。
所述地质雷达小波基函数的添加，在小波基函数添加之前，进行如下工作：①定义小波基函数的全称(familyName)；②定义小波基函数的缩写名(familyShortName)；③定义小波基函数的类型(familyWaveType)；④新建一个MATLAB中的mat文件(fileWaveName)；
在MATLAB语言工作平台上，可使用wavemngr命令添加一个新的小波函数，具体形式如下：wavemngr(‘add’，fn，fsn，wt，file)，其中fn为小波基函数的全称；fsn为小波基函数的缩写名；wt为小波基函数的类型；file为mat文件的文件名。
本发明实施例提供的S107中，小波变换时能密度法，假设构造的地质雷达小波基满足小波基函数的容许性条件，则地质雷达信号小波变换是完备的并且保持能量守恒；根据Moyal定理，有下式成立：
式中表明地质雷达信号小波变换后，系数平方的积分同被分析原始信号的能量成正比；因此，可以把|Wf(a,b)|2/Cψa2看作是地质雷达小波时间-尺度平面上的能量密度函数，即|Wf(a,b)|2ΔaΔb/Cψa2给出了以尺度因子a和平移因子b为中心的、尺度间隔为Δa、时间间隔为Δb的能量；根据小波变换能量密度的概念，也可以写成：
∫R|f(t)|2dt＝∫RE(b)db (10)
式中：
在小波变换中，尺度因子a在某种程度上对应于实际频率f；因此，给出了地质雷达信号所有频带的能量随平移因子b的分布情况，称为小波变换时能密度函数。
如图2所示，本发明实施例提供的基于地质雷达小波变换的隧道衬砌厚度识别系统包括：
模型构建模块1，假设隧道围岩为半无限空间连续的均匀介质，电磁波反射与折射均在二维平面内进行，建立隧道衬砌不同厚度混凝土结构地电模型。
正演模拟模块2，利用雷达天线对隧道衬砌不同厚度混凝土结构地电模型进行正演模拟，得到隧道衬砌混凝土结构的正演模拟雷达响应图像。
直达波去除处理模块3，根据电磁波在界面处的反射系数为负的原则，找到地质雷达信号最大的负峰值点位置，对地质雷达响应图像进行去除直达波处理。
雷达信号偏移处理模块4，根据偏移算法，在MATLAB语言工作平台上，对雷达响应图像进行偏移处理。
雷达信号规范处理模块5，采用去除背景信号与直流分量以及衰减补偿方法对地质雷达响应图像进行处理。
地质雷达小波基构造与添加模块6，根据地质雷达子波的衰减特性，采用高斯函数构造地质雷达小波基，并将其添加到小波分析工具箱中。
雷达信号定理识别模块7，根据已构造的地质雷达小波基，在MATLAB语言工作平台环境下，编写小波变换时能密度分析程序，采用该程序对单道信号进行小波变换时能密度法分析识别。
下面结合具体实施例对本发明的技术方案作进一步的描述。
(1)地质雷达小波基函数构造思想的提出
信号或函数f(t)的连续小波变换可以定义为：
从式(1)可以看出，连续小波变换的实质是将原始信号与小波基函数做卷积，而卷积运算从某种意义上反映了两个函数或信号之间的相似程度。如果小波基函数与原始信号的相似程度越高，则表示原始信号经过小波变换后小波变换系数的时频能量越集中，越能有效描述原始信号的时频局部化特征。为了更能反映地质雷达电磁波反射信号的时频特征，针对地质雷达发射子波的特点，构造与地质雷达发射子波波形相似度高的小波基函数，提出以下假设：根据地质雷达发射子波的主要特征，通过曲线拟合方法构造与地质雷达发射子波波形时频特征相似度高的数学表达式，使构造的曲线拟合函数满足小波基函数的允许条件，则可以作为一个新的小波基应用到地质雷达信号分析与处理中去。
(2)小波基函数的容许性条件
如果函数
能够作为一个小波基函数进行使用，则函数
必须是平方可积的，且函数傅里叶变换ψ(ω)需满足：
由式(2)可以看出，小波基函数容许性条件有两个主要特点：①波动性：由ψ(ω)的连续性及Cψ＜∞，可知ψ(0)＝0，即
表明小波基函数
的波形或幅值应是上、下波动的，且程度大致相等；②衰减性：由函数
属于能量有限空间
可知，小波基函数
具有快速衰减性。
(3)地质雷达发射子波的时频特征
这里，选用雷克子波作为地质雷达发射子波，发射子波的时域波形和功率谱图分别见图3和图4。从图3和图4可以看出，雷达发射子波具有一个半周期，在时域和频域中的能量都比较集中，且衰减速度较快，这符合小波基函数具有时频局部化分析能力的特点。因此，可以由雷克子波构造出与地质雷达反射信号相似度高的小波基函数可以作为一个新的小波基函数进行应用。
(4)基于雷克子波的地质雷达小波基函数的构造
在MATLAB语言工作平台上，选用Curve Fitting工具箱中带参数的衰减型高斯函数对如图3所示的地质雷达发射子波进行曲线拟合，然后比较各函数模型的拟合效果。拟合界面如图5所示。
经综合比选，得到的衰减型高斯函数形式为：
为了验证拟合高斯函数是否满足小波基函数容许性的条件，即波动性与衰减性的要求，首先将拟合高斯函数公式(3)进行积分，可得
由于积分计算后的结果，相对拟合高斯函数f(x)的幅值来说，量值相对较小，可忽略不计，由此可判断拟合高斯函数f(x)满足波动性要求，即∫Rf(t)dt≈0。然后，通过对拟合高斯函数进行积分，可知
(其中α＝1)，即拟合高斯函数f(x)满足衰减性要求。综上分析，拟合高斯函数f(x)满足小波基函数的容许性条件。因此，可以将拟合高斯函数f(x)添加到小波分析工具箱中作为一个新的小波基函数进行选用。
(5)地质雷达小波基函数的添加
在小波基函数添加之前，需进行如下工作：①定义小波基函数的全称(familyName)；②定义小波基函数的缩写名(familyShortName)；③定义小波基函数的类型(familyWaveType)；④新建一个MATLAB中的mat文件(fileWaveName)。
在MATLAB语言工作平台上，可使用wavemngr命令添加一个新的小波函数，具体形式如下：wavemngr(‘add’，fn，fsn，wt，file)，其中fn为小波基函数的全称；fsn为小波基函数的缩写名；wt为小波基函数的类型；file为mat文件的文件名。添加完的地质雷达小波基函数显示如图6所示。
(6)小波变换时能密度法
假设地质雷达小波基函数满足式(2)的容许性条件，则说明小波变换是完备的并且保持能量守恒。根据Moyal定理，有下式成立：
式(4)表明地质雷达小波变换系数平方的积分同被分析原始信号的能量成正比。因此，可以把|Wf(a,b)|2/Cψa2看作是地质雷达小波时间-尺度平面上的能量密度函数，即|Wf(a,b)|2ΔaΔb/Cψa2给出了以尺度因子a和平移因子b为中心的、尺度间隔为Δa、时间间隔为Δb的能量。根据小波能量密度的概念，式(4)也可以写成：
∫R|f(t)|2dt＝∫RE(b)db (5)
式中：
在小波变换中，尺度因子a在某种程度上对应于实际频率f。因此，式(6)给出了地质雷达信号所有频带的能量随平移因子b的分布情况，称为小波时能密度函数。
(7)隧道衬砌厚度中的定量识别
①模拟方案设计
假设隧道围岩为半无限空间连续的均匀介质，电磁波反射与折射均在二维平面内进行。为了研究隧道衬砌结构不同厚度在地质雷达正演模拟信号上的反射特征，设计了如图7所示混凝土结构的地电模型。
地电模型参数设置为：模拟区域范围7.2m×2.0m，左上角为坐标原点，水平横坐标为距离，竖向纵坐标为探测深度；将水平距离长为7.2m划分成六等份，每份长为1.2m，厚度值从左至右依次分别为0.40m、0.35m、0.30m、0.30m(含钢筋)、0.35m(含钢筋)和0.40m(含钢筋)；钢筋布置型号均为Ф20，每份钢筋间距均为0.24m，保护层厚度为50mm；模筑混凝土与喷射混凝土(围岩)衬砌的相对介电常数分别为6和20，电导率分别为0.01S/m和0.1S/m，导磁率均为1。
②不同厚度衬砌结构的雷达响应
根据隧道衬砌结构现场检测使用的地质雷达参数，模拟所采用的雷达天线中心频率为900MHz，边界吸收条件为完全匹配层，激励源采用Ricker子波，900MHz雷达天线的采样时窗为15ns，网格的空间步长均为0.005m，采样步长为0.01m，采样道数为701道。采用900MHz雷达天线对如图7所示的地电模型进行正演模拟，可以得到隧道衬砌混凝土结构的正演模拟雷达响应图像。
③地质雷达信号直达波的去除
直达波或表面波是地质雷达探测设备的首要干扰。如果目标体距离地表比较远时，直达波与目标反射信号在时间上是错开的，可以理解为直达波对目标反射信号影响较小；如果目标体距离地表比较近时，属于浅层目标体，这时直达波与目标反射信号在时间上有可能相互叠加干扰，将严重影响目标体的反射信号。因此，必须予以去除直达波。
地质雷达发射的脉冲电磁波从空气进入岩土层中，由于两种不同介质相对介电常数的差异，将在空气与二衬、二衬与初衬界面界面上产生反射现象，电磁波传播规律符合反射定律，产生的电磁波能量大小取决于电磁波反射系数
根据电磁波在界面处的反射系数为负的原则，找到地质雷达信号最大的负峰值点位置，对图8所示的地质雷达响应图像进行去除直达波处理，结果见图9。
④地质雷达信号偏移处理
设
为矩阵E(x,z＝0,t)的二维傅立叶变换，则：
在频率波数域内，深度z处的波场可以表示为：
设E(x,z,t)(E(x,t)在z方向的延拓)为
关于kx,w的二维傅立叶逆变换：
由式(9)，可得：
令t＝0，得到：
然后，根据频散关系式：
通过式(12)，可将
也就是
插值映射为
得到下式：
根据上述偏移算法，在MATLAB语言工作平台上，对如图9所示的雷达响应进行偏移处理，结果如图10所示。
⑤地质雷达信号的规范处理
采用去除背景信号与直流分量以及衰减补偿等方法对如图10所示地质雷达响应图像进行处理，结果见图11。
由于模筑混凝土衬砌与喷射混凝土衬砌(围岩)相对介电常数的差异，通过算法处理后使得电磁波信号在分界面处存在显著差异。从图11所示的地质雷达图像上可以看到较明显的连续反射面(层位线)。
⑥小波变换模极大值法
为了进一步研究地质雷达电磁波反射信号的瞬时位置，从图11所示衬砌混凝土结构雷达响应的依次抽取64条单道信号进行分析，结果见图12。电磁波由空气进入模筑混凝土衬砌，经模筑混凝土衬砌再进入喷射混凝土衬砌(围岩)，各界面的反射波依次被接收天线所接收。当模筑混凝土衬砌的相对介电常数小于喷射混凝土衬砌的相对介电常数、喷射混凝土衬砌的相对介电常数大于围岩的相对介电常数时，电磁波会在界面处发生相位反相。而从图12所示的地质雷达时间剖面上无法识别电磁波遇到界面产生反射波的瞬时位置(究竟处于波峰还是波谷)。
下面以1条单道信号为例进行分析，采用Matlab小波分析工具箱中的Db系列小波(Db2、Db4、Db6、Db8和Db10)分别对单道信号进行小波变换模极大值法识别，结果分别见图13～图17。将前文已构造的地质雷达小波基作为小波变换模极大值法用的小波基，对单道信号进行小波变换模极大值法识别，结果见图18。通过识别如图13～图18所示小波变换时能密度曲线上的双程旅行时间Δt＝t2-t1，利用距离公式
可以得到隧道衬砌混凝土结构的厚度值L，其中：c为电磁波在真空中的传播速度(3.0×108m/ns)，t1为直达波的反射位置(t1＝0)，t2为小波变换时能密度曲线上峰值点的位置，ε为介质的相对介电常数。分别将Db2、Db4、Db6、Db8、Db10小波和地质雷达小波变换模极大值法识别的厚度值与理论设计值进行比较，如表1所示。
表1不同方法获得的隧道衬砌厚度值
从表1可以看出，采用不同Db系列小波分析同一个地质雷达单道信号，识别结果不完全相同，甚至差异较大；与经典Db系列小波识别法相比，地质雷达小波识别法获得的图像分辨率较高，识别误差较小，相对误差小于设计值的5％。
⑦小波变换时能密度法的定量识别
根据已有的地质雷达小波基，提出一种基于地质雷达小波变换的时能密度法。在MATLAB语言工作平台环境下，编写小波变换时能密度分析程序。采用该程序对单道信号进行小波变换时能密度法分析，结果见图19。同理，可以对其他63条单道信号进行小波变换时能密度法分析，结果见图20。将小波变换时能密度曲线上峰值点连成一条线，可以获得隧道衬砌不同混凝土厚度的层位线，结果见图21。
(8)本方法的验证
湖南某隧道进口里程K5+382.000，出口里程K5+520.000，全长138m，其中K5+400.000～K5+505.000为暗挖段，长105m。两侧为明洞段，长33m。隧道建筑限界净宽13.5m，设双向两机动车道，两侧各设一条非机动车道，外侧设检修道。机动车道建筑限界净高5.0m，检修道净高2.5m。隧道暗挖段开挖跨度约16.3m，高度约11.5m，最大覆盖层厚度约26.5m。暗挖段采用新奥法工艺施工，结构采用复合式衬砌。初期支护以喷射混凝土、锚杆、钢筋网、工字钢为主要支护手段，二次衬砌采用C35钢筋混凝土，混凝土厚度为60cm，混凝土抗渗等级不小于P8。
采用LTD-2100型地质雷达，配以900MHz天线。整个系统由便携式主机、收发天线、综合控制电缆、测量轮、内置12V锂电池、数据采集和处理软件等组成。主机技术参数特点：一体化设计、体积小、重量轻、功耗低，实现滤波与放大等实时处理，操作简便，具备连续、点测或测量轮控制等测量方式，实时二维图像显示，扫描速率16Hz、32Hz、64Hz和128Hz可调，记录道长度256、512、1024和2048可调，输入带宽1Hz～16kHz，动态范围-7dB～130dB。
采用地质雷达对整个隧道二次衬砌厚度进行检测，限于篇幅，这里选取一段5m进行分析，原始地质雷达剖面如图22所示。采用去直达波、增益、去背景信号、滤波和偏移处理后，得到的地质雷达剖面图见图23。图24为从图23中抽取的40条单道信号的波形图。采用小波变换时能密度法对图24中的40条单道信号进行分析与处理，得到隧道衬砌厚度识别的层位图(见图25)，并将隧道衬砌的理论设计值60cm进行比较，就以很明显地掌握隧道衬砌厚度是否满足要求。
在本发明的描述中，除非另有说明，“多个”的含义是两个或两个以上；术语“上”、“下”、“左”、“右”、“内”、“外”、“前端”、“后端”、“头部”、“尾部”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。此外，术语“第一”、“第二”、“第三”等仅用于描述目的，而不能理解为指示或暗示相对重要性。
以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，都应涵盖在本发明的保护范围之内。
Claims (4)
1.一种基于地质雷达小波变换的隧道衬砌厚度识别方法，其特征在于，所述基于地质雷达小波变换的隧道衬砌厚度识别方法包括：
隧道围岩为半无限空间连续的均匀介质，电磁波反射与折射均在二维平面内进行，建立隧道衬砌不同厚度混凝土结构地电模型，对隧道衬砌结构不同厚度在地质雷达正演模拟信号上的反射特征进行分析；
隧道衬砌混凝土结构不同厚度的雷达响应，利用雷达天线对隧道衬砌不同厚度混凝土结构地电模型进行正演模拟，得到隧道衬砌混凝土结构的正演模拟雷达响应图像；
根据电磁波在界面处的反射系数为负的原则，找到地质雷达信号最大负峰值点位置，对地质雷达响应图像进行去除直达波处理；
根据偏移算法，在MATLAB语言工作平台上，对地质雷达响应图像进行偏移处理；
地质雷达信号的规范处理，采用去除背景信号与直流分量以及衰减补偿方法对地质雷达响应图像进行去噪与弱信号增强处理；
地质雷达小波基函数的构造与添加，根据地质雷达发射子波的衰减特性，采用高斯函数拟合构造地质雷达小波基，在MATLAB语言工作平台环境下，将构造的地质雷达小波基函数添加到小波分析工具箱中，以供小波变换时选用；
地质雷达信号的定量识别，根据已构造的地质雷达小波基，提出一种基于地质雷达小波变换的时能密度方法，在MATLAB语言工作平台环境下，编写小波变换时能密度分析程序，采用该程序对地质雷达单道信号进行小波变换时能密度法分析识别；
所述地电模型参数设置为：模拟区域范围7.2m×2.0m，左上角为坐标原点，水平横坐标为距离，竖向纵坐标为探测深度；
将水平距离长为7.2m划分成六等份，每份长为1.2m，厚度值从左至右依次分别为0.40m、0.35m、0.30m、0.30m、0.35m和0.40m；
钢筋布置型号均为Ф20，每份钢筋间距均为0.24m，保护层厚度为50mm；模筑混凝土与喷射混凝土衬砌的相对介电常数分别为6和20，电导率分别为0.01S/m和0.1S/m，导磁率均为1；
所述根据隧道衬砌结构现场检测使用的地质雷达参数，模拟所采用的雷达天线中心频率为900MHz，边界吸收条件为完全匹配层，激励源采用Ricker子波，900MHz雷达天线的采样时窗为15ns，网格的空间步长均为0.005m，采样步长为0.01m，采样道数为701道；
所述地质雷达发射的脉冲电磁波从空气进入岩土层中，由于两种不同介质相对介电常数的差异，将在界面上产生反射现象，其规律符合电磁波反射定律，产生的电磁波能量大小取决于反射系数根据电磁波在界面处的反射系数为负的原则，找到地质雷达信号最大的负峰值点位置，对地质雷达响应图像进行去除直达波处理；
所述地质雷达信号偏移处理过程为：
设
为矩阵E(x,z＝0,t)的二维傅立叶变换，则：
在频率波数域内，深度z处的波场表示为：
设E(x,z,t)(E(x,t)在z方向的延拓)为
关于kx,w的二维傅立叶逆变换：
由式
可得：
令t＝0，得到：
然后，根据频散关系式：
通过式
可将
也就是
插值映射为
得到下式：
根据上述偏移算法，在MATLAB语言工作平台上，对雷达响应图像进行偏移处理；
所述通过提取小波变换时能密度曲线上的双程旅行时间Δt＝t2-t1，利用距离公式
得到隧道衬砌混凝土结构的厚度值L；其中：c为电磁波在真空中的传播速度3.0×108m/ns，t1为直达波的反射位置t1＝0，t2为小波变换时能密度曲线上峰值点的位置，ε为介质的相对介电常数；
将不同厚度隧道衬砌混凝土结构正演模拟识别结果与理论设计值进行比较，采用小波变换时能密度法识别计算得到的隧道衬砌混凝土结构厚度的相对误差值均不超过设计值的5％；
所述地质雷达小波基函数的构造过程为：
在MATLAB语言工作平台上，选用Curve Fitting工具箱中带参数的高斯函数对地质雷达发射子波进行曲线拟合，然后比较各函数模型的拟合效果；
经综合比选，拟合得到的高斯函数形式为：
f(x)＝-23.64*exp(-((x+0.001558)/3.248)2)+0.009297*exp(-((x-1.954)/2.245)2)+；24.64*exp(-((x+0.002082)/3.116)2)+
要想了解构造的拟合高斯函数是否可以作为小波分析用的小波基，需验证高斯函数是否满足小波基的容许性条件，即波动性与衰减性的要求；首先，对拟合高斯函数进行积分，得：
通过判断，拟合高斯函数f(x)满足波动性要求，即∫Rf(t)dt≈0；然后，通过对拟合高斯函数进行积分计算，可知：
dx＝9.4233*10-4，其中α＝1，则拟合高斯函数f(x)满足衰减性要求；
通过上述分析，拟合高斯函数f(x)满足小波基函数容许性的条件；将拟合的高斯函数f(x)添加到小波分析工具箱中作为一个新的小波基函数进行选用。
2.如权利要求1所述的基于地质雷达小波变换的隧道衬砌厚度识别方法，其特征在于，所述地质雷达小波基函数的添加，在小波基函数添加之前，进行如下工作：①定义小波基函数的全称；②定义小波基函数的缩写名；③定义小波基函数的类型；④新建一个MATLAB中的mat文件；
在MATLAB语言工作平台上，可使用wavemngr命令添加一个新的小波函数，具体形式如下：wavemngr，其中fn为小波基函数的全称；fsn为小波基函数的缩写名；wt为小波基函数的类型；file为mat文件的文件名。
3.如权利要求1所述的基于地质雷达小波变换的隧道衬砌厚度识别方法，其特征在于，所述小波变换时能密度法，假设地质雷达小波基函数满足容许性条件，则地质雷达小波变换是完备的并且保持能量守恒；根据Moyal定理，有下式成立：
式中表明小波变换系数平方的积分同被分析原始信号的能量成正比；把|Wf(a，b)|2/Cψa2看作是地质雷达小波变换时间-尺度平面上的能量密度函数，即|Wf(a，b)|2ΔaΔb/Cψa2给出了以尺度因子a和平移因子b为中心的、尺度间隔为Δa、时间间隔为Δb的能量；根据小波变换能量密度的概念，写成：
∫R|f(t)|2dt＝∫RE(b)db；
式中：
在小波变换中，尺度因子a在某种程度上对应于实际频率f；因此，给出了地质雷达信号所有频带的能量随平移因子b的分布情况，称为小波时能密度函数。
4.一种实施如权利要求1～3任意一项所述的基于地质雷达小波变换的隧道衬砌厚度识别方法的系统，其特征在于，所述基于地质雷达小波变换的隧道衬砌厚度识别系统包括：
模型构建模块，假设隧道围岩为半无限空间连续的均匀介质，电磁波反射与折射均在二维平面内进行，建立隧道衬砌不同厚度混凝土结构地电模型；
正演模拟模块，利用雷达天线对隧道衬砌不同厚度混凝土结构地电模型进行正演模拟，得到隧道衬砌混凝土结构的正演模拟雷达响应图像；
直达波去除处理模块，根据电磁波在界面处的反射系数为负的原则，找到地质雷达信号最大的负峰值点位置，对地质雷达响应图像进行去除直达波处理；
雷达信号偏移处理模块，根据偏移算法，在MATLAB语言工作平台上，对雷达响应图像进行偏移处理；
雷达信号规范处理模块，采用去除背景信号与直流分量以及衰减补偿方法对地质雷达响应图像进行处理；
地质雷达小波基的构造与添加，根据地质雷达子波在时频域内的衰减特性，采用高斯函数构造地质雷达小波基，并将其添加到小波分析工具箱中，以供地质雷达小波变换时选用；
雷达信号定量识别模块，根据已构造的地质雷达小波基，提出一种基于地质雷达小波变换的时能密度方法，在MATLAB语言工作平台环境下，编写小波变换时能密度分析程序，采用该程序对地质雷达单道信号进行小波变换时能密度法识别。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010908318.5A CN111981968B (zh) | 2020-09-02 | 2020-09-02 | 基于地质雷达小波变换的隧道衬砌厚度识别方法及系统 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010908318.5A CN111981968B (zh) | 2020-09-02 | 2020-09-02 | 基于地质雷达小波变换的隧道衬砌厚度识别方法及系统 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN111981968A CN111981968A (zh) | 2020-11-24 |
| CN111981968B true CN111981968B (zh) | 2022-05-10 |
Family
ID=73448697
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202010908318.5A Active CN111981968B (zh) | 2020-09-02 | 2020-09-02 | 基于地质雷达小波变换的隧道衬砌厚度识别方法及系统 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN111981968B (zh) |
Families Citing this family (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112488120A (zh) * | 2020-12-02 | 2021-03-12 | 山东大学 | 一种基于地质雷达图像预处理方法 |
| CN112684443B (zh) * | 2021-03-15 | 2021-05-28 | 中国铁建重工集团股份有限公司 | 一种衬砌结构病害检测方法、装置及介质 |
| CN113970683B (zh) * | 2021-09-27 | 2024-09-10 | 国网浙江省电力有限公司电力科学研究院 | 基于小波变换的电缆检测方法、设备及介质 |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102288625A (zh) * | 2011-07-08 | 2011-12-21 | 李典基 | 一种渠道混凝土衬砌无损检测系统 |
| CN106382898A (zh) * | 2016-08-29 | 2017-02-08 | 西安科技大学 | 一种寒区隧洞衬砌外覆冰层厚度的测量方法及系统 |
| CN109799539A (zh) * | 2019-01-11 | 2019-05-24 | 湖南城市学院 | 基于雷达小波变换时能密度的隧道地质体三维识别方法 |
-
2020
- 2020-09-02 CN CN202010908318.5A patent/CN111981968B/zh active Active
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102288625A (zh) * | 2011-07-08 | 2011-12-21 | 李典基 | 一种渠道混凝土衬砌无损检测系统 |
| CN106382898A (zh) * | 2016-08-29 | 2017-02-08 | 西安科技大学 | 一种寒区隧洞衬砌外覆冰层厚度的测量方法及系统 |
| CN109799539A (zh) * | 2019-01-11 | 2019-05-24 | 湖南城市学院 | 基于雷达小波变换时能密度的隧道地质体三维识别方法 |
Non-Patent Citations (3)
| Title |
|---|
| 基于小波变换时能密度法的隧道空洞充填物识别;张胜;《煤炭学报》;20191115;正文1-5节 * |
| 基于车载非接触式雷达的公路隧道检测技术研究;刘小龙;《中国优秀硕士学位论文全文数据库》;20200616;正文第3.6.1节 * |
| 探地雷达钢筋正演模拟及偏移成像研究;邓方进;《隧道建设》;20190912;正文0-3节 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN111981968A (zh) | 2020-11-24 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN111981968B (zh) | 基于地质雷达小波变换的隧道衬砌厚度识别方法及系统 | |
| CN101968550B (zh) | 基于阵列光纤传感器的岩层识别装置及方法 | |
| Li et al. | An overview of ahead geological prospecting in tunneling | |
| Xiang et al. | GPR evaluation of the Damaoshan highway tunnel: A case study | |
| CN108152854B (zh) | 一种基于微动功率谱密度的无损探测方法及其应用 | |
| CN102590874B (zh) | 一种山地采煤沉陷水田地表裂缝探测方法 | |
| CN109521479A (zh) | 一种地下连续墙渗漏水的非开挖检测方法 | |
| CN109799539B (zh) | 基于雷达小波变换时能密度的隧道地质体三维识别方法 | |
| CN113419294A (zh) | 一种多维度岩溶特殊地质综合探测方法 | |
| CN115291200B (zh) | 一种基于数字显示的埋深管道定位方法 | |
| Geng et al. | Identifying void defects behind Tunnel composite lining based on transient electromagnetic radar method | |
| CN108614307B (zh) | 一种利用综合物探技术的仰拱质量检测方法及系统 | |
| Zhu et al. | Research on accuracy and error analysis of coal and rock strata detection based on air-coupled GPR | |
| Wang et al. | Research on multi-frequency ultrasonic scanning detecting technology of cavity in the test borehole | |
| Liu et al. | Applied research of comprehensive advance geological prediction in Daluoshan water diversion tunnel | |
| Zhang et al. | Time-frequency analysis of GPR simulation signals for tunnel cavern fillings based on short-time Fourier transform | |
| Zhang et al. | Wavelet Packet Analysis of Ground-Penetrating Radar Simulated Signal for Tunnel Cavity Fillings. | |
| CN111538073A (zh) | 一种隧道不良地质体超前地质预报探测模拟装置及方法 | |
| Liu et al. | A system for inspecting karst voids during construction of cast-in-place pile foundations | |
| CN111927469A (zh) | 一种岩溶隧道涌水超前预报系统化施工方法 | |
| CN212364603U (zh) | 一种隧道不良地质体超前地质预报探测模拟装置 | |
| Fan et al. | Application of Seismic Channel Wave Technology on Small Structure Exploration in Coal Mine | |
| CN117784264B (zh) | 基于探地雷达数据的输电塔间地下病害定位方法及系统 | |
| CN112987082B (zh) | 一种隧道施工不良地质三维地震波超前探测方法 | |
| Klose et al. | Spatial predictions of geological rock mass properties based on in-situ interpretations of multi-dimensional seismic data |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant | ||
| EE01 | Entry into force of recordation of patent licensing contract |
Application publication date: 20201124 Assignee: Hunan Fangwei Construction Engineering Co.,Ltd. Assignor: Hunan City University Contract record no.: X2023980052843 Denomination of invention: Tunnel lining thickness identification method and system based on ground penetrating radar wavelet transform Granted publication date: 20220510 License type: Exclusive License Record date: 20231219 |
|
| EE01 | Entry into force of recordation of patent licensing contract |