随着我国交通基础设施的发展，我国道路建设的速度和规模早已居于世界第一[1]。然而，由于道路建设施工质量不佳、建设工期和线路选择不合理以及道路基础地质条件复杂等因素，造成近几年我国城市道路塌陷事故频发，严重威胁人民群众的生命安全和城市交通秩序[2]。城市道路塌陷的主要诱因是道路地下空洞，地下空洞形成的原因主要有以下几方面：(1)地质、地下水作用；(2)长期振动作用；(3)存在地下人防工程；(4)基建施工的影响；(5)管道施工特别是非开挖施工；(6)道路施工中碾压不实；(7)地下管道破损。由此可见，道路塌陷隐患有成因复杂、隐蔽性强、预警难和探测难等特点。采取有效、经济的手段对道路地下隐患的位置及规模进行准确探测以避免类似事故发生已迫在眉睫。

探地雷达作为一种检测速度快、无破损的物探方法，在市政工程地基病害探测中得到了广泛应用[3]。传统的探地雷达信号处理方法是通过基于傅里叶分析法的探地雷达自带软件处理得到的探地雷达剖面图，分析地下空洞的发育情况。城市道路的市政附属设备设施较多，如地下管线（如雨水、污水、热力、燃气、供水金属管线）和通电的设备或电缆（如通信/供电电缆、电子眼、路灯、交通信号灯、电子屏等），都是常见的电磁干扰源，严重干扰探地雷达的电磁信号，在探地雷达波形剖面图上形成与真实病害相似的图形，甚至会覆盖有效信号，造成对资料解译的漏判或误判，给准确识别城市道路地下隐患带来阻碍[4]。余志雄等[5]采用复信号分析技术对探地雷达信号进行处理，但是该方法会放大信号噪声，不能增强城市道路探地雷达的抗干扰能力。张胜等[6]提出一种基于雷达小波基的小波变换时能密度法，但仍存在不同尺度小波基导致信号处理结果不同的多解性问题，过于依赖人员的经验。凌同华等[7]采用希尔伯特‒黄变换方法，提取各固有模式函数分量的瞬时剖面图，结合探地雷达单道信号的瞬时频率、相位和幅值，达到对探地雷达图像的定量解释，但仍无法避免城市道路中存在的强电磁干扰。张先武等[8]使用广义S变换对实测探地雷达数据进行了层位识别，取得了较好的效果，但尚未将广义S变换运用于探地雷达道路地下空洞探测。

针对城市道路地下空洞的特点，采用正演模拟、广义S变换和功率谱分析等方法分析探地雷达对城市道路地下空洞的波场响应、时频响应和功率谱衰减特征。结合工程实例，对提出的波场和频谱特性分析方法进行系统阐述，验证所提出方法的可行性和准确性。

1、道路地下地层物性特征

以河南某机场停机坪的工程地质条件为例，道路下方地层主要有：(1)回填土，黄褐色，成分以粉土、粉砂为主；(2)强夯碎石垫层，该区域原始地貌为沟塘，清淤后采用粉土、粉砂夹少量碎石回填，再进行强夯处理，碎石大小不一；(3)回填土，褐黄色，成分以粉土、粉砂为主；(4)强夯块石垫层，该区域原始地貌为沟塘，清淤后采用块石夹粉土回填，再进行强夯处理，块石大小差别较大，土、石体积比例差别也较大，较不均匀；(5)回填土，黄褐色，成分以粉土、粉砂为主；(6)粉土，褐黄色，稍湿，密实，含锈黄、灰褐斑，偶见蜗牛碎壳及少量钙核，局部夹少量粉质黏土。道路地下常见物质材料的相对介电常数（实部）和静态电导率如表1所示[9]。金属为导体，对电磁波具有屏蔽作用，电磁波无法穿透。一般地，探地雷达信号经过金属区域会产生强能量反射界面[10]。 

表1 道路地下常见物质材料的相对介电常数（实部）和静态电导率  

探地雷达电磁波在地下地层传播过程中，当遇到不同电性介质的界面（波阻抗界面）时，会发生反射和透射，界面两侧介质的电性差异越大，反射系数也越大，反射波能量强弱取决于反射系数。反射系数的数学表达式为

式中：r为反射系数；εr1和εr2为反射界面两侧的相对介电常数。

由表1和式（1）可知，空洞（主要充填空气或水）与回填土、强夯碎石垫层、强夯块石垫层及粉土之间具有较大的电性差异，这为探地雷达探测道路地下空洞提供了基本的物理条件。

2、道路地下病害波场特征模拟

2.1 道路地下空洞发育区域模型

对于探地雷达探测道路地下病害波场特性的分析，采用时域有限差分（FDTD）法和复频移卷积完全匹配层（CFS‒PML)[11]的吸收边界条件[12]。

复频移卷积完全匹配层通过构造一种非物理的吸收媒介与FDTD的网格外部进行连接，这种吸收媒介的波阻抗与向外散射波的频率和入射角无关，能较好地吸收低频信号在边界的反射，相较于常规的完全匹配层、单轴各向异性完全匹配层等吸收边界，复频移卷积完全匹配层对隐失波、低频波等大角度掠射波有更好的吸收效果[13]。

根据实际道路地下病害发育区的结构与电性参数，建立不同病害发育区模型，如图1所示。模型Ⅰ为道路下方存在充填水和充填空气的空洞模型，模型Ⅱ为回填不密实和地下埋藏金属物的模型。为定量分析探地雷达波场特征，将模型Ⅰ中空洞发育区域设置为形状规则的圆形，分析充填不同介质空洞的探地雷达波场响应和时频分布特征。模型Ⅱ设置了回填不密实和埋藏金属物区域，分析探地雷达对回填不密实和埋藏金属物的波场响应特征。

图1 道路地下空洞发育模型（模型Ⅰ）示意图   

模型Ⅰ和模型Ⅱ采用FDTD法分别对地下充填水、充填空气的空洞以及回填不密实和地下埋藏金属物进行正演模拟，模型的宽度×深度为8 m×10m。为了模拟实际道路地下空间，设置层厚为0.2 m的混凝土层，模拟道路路面。网格步长为Δx=Δy=0.01 m，中心频率为250 MHz，时窗为250 ns，总道数为40，天线间距为1.0 m，道间距为0.2 m。为了更好地将数值模拟结果与实测信号对比，以上数值模拟参数均与探地雷达在道路空洞检测工程应用中实测的参数相一致。正演模拟所用的道路地下空洞发育区域常见介质电性参数参考表1，围岩的平均深度分别设计为：(1)回填土，深度0.2～2.0 m;(2)强夯碎石垫层，深度2.0～2.7 m;(3)回填土，深度2.7～6.0 m;(4)强夯块石垫层，深度6.0～7.9 m;(5)回填土，深度7.9～8.8 m;(6)粉土，深度8.8～10.0 m。时间‒深度转换环节中采用的回填土电磁波波速为0.08 m·ns-1。同时，为从不同角度分析空洞发育区域探地雷达响应特征，对经过空洞区域和埋藏金属物区域的单道探地雷达信号进行广义S变换，分析时频变化特征。

2.2 道路地下空洞发育区域波场特征模拟

图1为道路地下空洞发育模型示意图。模型Ⅰ中，道路地下充气空洞球心为（2.5,2.5) m，半径为1 m，即充气空洞顶界面的深度为1.5 m；道路地下充水空洞球心为（5.5,2.5) m，半径为1.0 m，即充水空洞顶界面的深度也为1.5 m。

为了验证正演结果和理论模型在深度上的对应性，此处均为时间‒深度转换后的数值模拟结果，时间‒深度转换公式为

式中：S为深度；v1为地下地层的电磁波波速；t为时间。

由图2可知，道路地下空洞发育的正演结果能较好地反映道路地下空洞的顶界面，地下充气空洞和地下充水空洞正演的反射波同相轴均呈弧形反射形态，与4节中的探地雷达工程实测结果（见图7）相似。弧形反射的顶界面深度与模型设置的顶部深度非常一致，均为1.5 m，但充气空洞的反射波能量强于充水空洞。由表1可知，回填土的相对介电常数ε1=20，空气的相对介电常数ε2=1，水的相对介电常数ε3=80，依据式（1）计算得到充气空洞与回填土界面的反射系数r1=0.634，而充水空洞与回填土界面的反射系数r2=0.309，从定量角度进一步解释了充气空洞反射波能量更强的原因。图2中，出现4个较为明显的横向连续同相轴，推断为道路地下不同介质分界面反射波形，只有2 m深度界面和理论模型中第一层回填土深度较为一致，均为2 m，其余层面位置均与理论模型层面深度有偏差。在时间‒深度转换环节采用的是回填土的电磁波波速，而强夯碎石垫层和强夯块石垫层的电磁波波速都要大于回填土的电磁波波速，通过相应地层的真实时间更短，所以数值模拟结果中出现强夯碎石垫层和强夯块石垫层与回填土的分界面深度较理论模型更深的情况。

图2 道路地下空洞发育模型（模型Ⅰ）正演成果   

2.3 道路地下回填不密实和埋藏金属物区域波场特征模拟

图3为道路地下回填不密实和埋藏金属物区域的模型示意图。模型Ⅱ中，分别设计回填不密实区域和埋藏金属物区域。回填不密实深度范围为1.8～3.2m，距离（x轴）位置为1.9～3.3 m；埋藏金属物深度范围为2.0～3.0 m，距离（x轴）位置为5.0～6.0 m。

图3 道路地下回填不密实和埋藏金属物区域发育模型（模型Ⅱ）示意图   

由图4可知，正演结果中的顶界面强反射波能较好地反映回填不密实区域和埋藏金属物区域的顶界面，回填土不密实和埋藏金属物区域的顶界面反射波同相轴均呈现出中间直线、左右边界弯曲的形态，但金属物的反射波能量更强。回填不密实内部区域波形连续性差，波形较杂乱，模型Ⅱ中回填土不密实区域的反射波同相轴形态与4节的探地雷达工程实测结果（见图7）相似。埋藏金属物区域内部未出现明显的波形，可能是金属物将向下发射的大部分雷达电磁波屏蔽而造成的。探地雷达剖面与实际地层深度偏差的原因在2.2节中已经详细描述，此处不再赘述。

图4 道路地下回填不密实和埋藏金属物区域发育模型（模型Ⅱ）正演成果  

3、时频分析

通过对模拟的单道探地雷达数据进行时频分析，可以更详细地得到频率随着时间的局部变化规律，并显示频率的振幅。采用在S变换[14]基础上加以改进的广义S变换[8]进行时间‒频率域信号处理，计算式如下所示：

式中：t、τ为时间；f为频率；βG和βL分别为时窗和频率调节参数。通过改变βG值可以改变广义S变换的时间分辨率，通过改变βL值可以改变广义S变换的频率分辨率。

3.1 空洞发育区域时频分析

对模型Ⅰ的第12道数据进行广义S变换，该道在模型Ⅰ距离（x轴）2.4 m处，经过空洞发育的中心区域，得到的时间‒频率变化如图5所示。由图5可知，采样时间25～60 ns的中心频率在0.025～0.075GHz之间，其余频率成分分布较均匀。探地雷达信号在经过空洞上层介质区域后，发生了较为明显的高频衰减，所以在时频图中只发现空洞发育区域的低频信息。中心频率在时间轴上的分布与图2中的数值模拟结果较为相符，都在25～60 ns范围内。

图5 模型Ⅰ第12道数据的广义S变换时频图   

3.2 埋藏金属物区域时频分析

对模型Ⅱ的第27道数据进行广义S变换，该道在模型Ⅱ距离（x轴）5.4 m处，经过埋藏金属物发育的中心区域，得到的时间‒频率图如图6所示。由图6可知，采样时间50～125 ns内，出现两簇高幅值能量团，分别是0.025～0.050 GHz和0.060～0.150GHz，频率分布不集中。通过图2和图4的数值模拟可知，埋藏金属物区域边界处也会出现弧形反射，如果仅分析波形图，就会对地下异常体造成误判，应进一步分析单道的时频分布特征，去伪存真。

由于回填不密实区域地质模型较为复杂，因此频率成分也较为复杂，未能显示出较明显的规律。复杂地质模型的探地雷达时频变化规律还有待进一步研究。

4、工程实例分析

为验证理论模型正演模拟及时频分析的准确性，选取河南省某机场停机坪工程探地雷达空洞检测剖面图并对其进行时频分析。

4.1 工程地质条件

工程概况：该场地建设前地势起伏较大，自然地面标高144.08～148.78 m，第(1)层杂填土分布范围及厚度不均匀，主要分布在西南部，层底标高139.85～147.60 m。场地西侧有丈八沟通过（沟宽约40 m），故场地整体都要进行深浅不一的回填，低洼处沟塘区回填深度达11 m。根据场地西侧现状站坪的地面标高，场地的地面整平标高为148.23～150.78 m。

图6 模型Ⅱ第27道数据的广义S变换时频图   

由地形和地质条件所决定，需对场地地基进行处理。采用回填石渣垫层强夯法进行地基处理，石渣垫层厚度为0.5 m，用2 000 k N·m、3 000 k N·m点夯及800 k N·m满夯；采用素土分层回填碾压，处理厚度为2.0 m。回填部分约为70 000 m2，尺寸为256m×295 m。停机坪地层及描述详见2.1节。

场地西南侧局部出现了地面塌陷，塌陷区域近似椭圆形，塌陷区东西向长约12 m，南北向长约8m，深约1.8～2.0 m。塌陷区西边缘距离T2航站楼约80 m，南边缘距离西南指廊南端约100 m。本次塌陷区域西侧约36 m处为正在施工的管沟，管沟宽约2.5 m，深约2.5～3.0 m，沟内可见强夯碎石垫层出露，沟内局部有积水。

探测仪器为加拿大Pulse EKKO PRO型探地雷达仪，采用点测方式，天线中心频率为250 MHz，天线间距为1.0 m，点距为0.2 m，采样时窗为250 ns，叠加次数为64。

4.2 地下空洞和回填不密实病害探测实例

图7为探地雷达探测机场停机坪地下空洞和回填不密实病害的典型剖面。为简化描述，将该条探地雷达测线命名为测线1。由图7可知，测线1从左至右9.2～13.4 m、深度1.0～3.0 m内，探地雷达波形呈弧形反射，振幅能量强。深度1.0 m处，可明显识别出空洞的顶部反射波同相轴；深度2.0 m处，又出现一个明显的弧形反射，推断为空洞下界面的反映。实测波形顶部反射波形态与图2中空洞的数值模拟波形相似。

由图7可知，测线1从左至右15.0～19.0 m、深度3.0～5.6 m内，波形连续性差，波形较杂乱，推断为回填不密实的病害。实测波形形态与图4中回填不密实的数值模拟波形相似。

对测线1的第56道数据进行广义S变换，该道在测线1距离（x轴）11.2 m处，经过空洞中心区域，得到的时间‒频率变化规律如图8所示。由图8可知，采样时间25～60 ns的中心频率为0.025～0.075GHz，其余频率成分分布较均匀。探地雷达信号在经过空洞上层介质区域后发生了较为明显的高频衰减，所以在时频图中只发现了空洞发育区域的低频信息。中心频率在时间轴的分布与图7的空洞发育区域在时间轴的分布较为相符，都在25～60 ns范围内。实测数据时频分布规律与数值模拟结果（见表2）相似，证实了数值模拟的可靠性。

图7 探地雷达地下空洞实测剖面图     

表2 测线1探测结果与实际结果对比  

经钻孔验证，测线1从左至右9.0～13.0 m、深度1.0～3.5 m内，发育一个小空洞，小空洞内为空气充填（见图9）；测线1从左至右15.0～18.5 m、深度3.0～5.0 m内，回填土较为松散，出现不密实的情况。由表2可知，探测结果的形状及大小与开挖揭示结果总体较相似。通过探地雷达实测结果的剖面图和时频图，结合数值模拟结果，可以深入分析地下空洞发育区域探地雷达波场特征和时频变化特征。探地雷达探测地下空洞达到了较好的效果，说明数值模拟分析可以有效指导工程实测工作。回填不密实病害区域由于介质分布较为复杂，波场特征有一定规律，但时频变化没有表现出较明显的规律，有待进一步研究。

4.3 埋藏金属物探测实例

图10为探地雷达探测机场停机坪地下埋藏金属物的典型剖面图。为简化描述，将该条探地雷达测线命名为测线2。由图10可知，测线2从左至右49.0～54.4 m、深度2.0～10.0 m内，探地雷达波形呈弧形反射，振幅能量强，并在深度轴上重复出现多条弧形反射同相轴，推断在该范围内可能埋藏一个方形金属物。反射波形态与图4中埋藏金属物的数值模拟波形相似。

图8 探地雷达测线1第56道数据时频图   

对测线2的第260道数据进行广义S变换，该道在测线2距离（x轴）52.0 m处，经过埋藏金属物中心区域，得到的时频变化规律如图11所示。采样时间50～150 ns内，出现两簇高幅值能量团，分别是0.015～0.040 GHz和0.050～0.150 GHz。频率分布不集中。实测数据时频分布规律与数值模拟结果相似，证实了数值模拟的可靠性。

由表3可知，经钻孔验证，测线2从左至右49.5～54.0 m、深度2.0～5.0 m内，埋藏一个废弃的钢管，探测结果与开挖结果总体较相似。通过分别分析实测探地雷达地下空洞和埋藏金属物的剖面图和时频图，发现了两者在波场和时频分布特征的区别。

图9 回填土层小空洞发育照片   

图1 0 探地雷达埋藏金属物实测剖面图   

图1 1 探地雷达测线2第260道数据时频图    

表3 测线2探测结果与实际结果对比  

4.4 探地雷达信号功率衰减分析

探地雷达的设计方式不同于常规雷达，探地雷达受到探测目标体近距离分布和探测环境的影响，电磁波信号在传播过程中会衰减。路径传播损耗Lp(λ)会较大地影响探地雷达最大探测深度，其值是信号波长的函数。探地雷达接收功率计算式为[13]

式中：PR为接收功率；PT为发射功率；G为天线增益；σ为探测目标体的雷达散射截面；λ为工作波长；εr为媒质的相对介电常数；R为目标体的距离。

通过自行编写的Matlab程序计算相应记录道的分析信号，从而计算其瞬时功率（即为接收功率PR在某一个时刻的值），并得到中值和平均值衰减曲线，即剖面中相应记录道的中值和平均瞬时功率。

图12为地下空洞发育区域（测线1距离（x轴）8.0～14.0 m）瞬时功率衰减曲线。由图12可知，当地下有空洞发育时，瞬时功率变化存在以下4种情况：在0～40 ns时，随着时间的变化，均值呈对数形式快速降低，原因可能是，探地雷达信号在空洞上层介质（回填土）传播相对于空洞是强损耗媒质，而空洞是无损耗媒质，导致瞬时功率降低，时间与图7中空洞顶界面时间较为一致；在40～80 ns内，随着时间的变化，均值呈对数形式快速增加；在80～110 ns内，达到峰值；在110～300 ns内，随着时间的变化，均值呈对数形式缓慢降低。

图13为地下埋藏金属物区域（测线2距离（x轴）48.0～55.0 m）的瞬时功率衰减曲线。由图13可知，当地下埋藏金属物时，瞬时功率变化存在以下5种情况：在0～50 ns内，随着时间的变化，均值呈对数形式缓慢增加，中值呈对数形式略有降低，原因可能是，探地雷达信号传播到50 ns时，遇到金属物这种强电磁干扰媒质，导致瞬时功率发生异常，时间与图10中埋藏金属物的顶界面时间较为一致；在50～75ns内，随着时间的变化，均值呈对数形式快速增加；在75～90 ns内，达到峰值；在90～125 ns内，随着时间的变化，均值呈对数形式快速降低；在125～300 ns内，均值忽高忽低，呈无规律的变化特征，原因是金属这种强电磁干扰源对探地雷达信号尾部造成不规则干扰，该现象在很多探地雷达实测波形图中都可以发现（见图10）。

图1 2 地下空洞发育区域瞬时功率衰减曲线   

图1 3 地下埋藏金属物区域瞬时功率衰减曲线   

由于回填不密实区域地质模型较为复杂，瞬时功率衰减机制也较为复杂，未能显示出较明显的规律。复杂地质模型的探地雷达功率谱变化规律还有待进一步研究。

综上可知，探地雷达在地下空洞发育区域和地下埋藏金属物区域，瞬时功率衰减特征都是不相同的。在遇到类似难以区分的异常时，可通过分析功率衰减机制，对真、假异常进行进一步筛选，去伪存真。

5、结论

(1）根据充填不同媒质空洞和埋藏金属物探地雷达剖面图中反射波的形态、振幅强弱、相应的时频变化特征和信号功率谱衰减特性，能有效地探测并区分空洞和埋藏金属物发育区域的准确位置和构造。

(2）地下空洞发育区域的中心频率集中在某个低频区间，出现高频衰减，其余频率成分分布较均匀；埋藏金属物区域出现两簇高幅值能量团，频率成分分布较分散。回填不密实区域由于其地质模型较为复杂，频率成分也较为复杂，因此未能显示出较明显的规律。复杂地质模型的探地雷达时频变化规律还有待进一步研究。

(3）数值模拟和工程实测的对比结果表明：探地雷达经过地下空洞和埋藏金属物的道有不同的波场和时频变化特征，而且实测数据和数值模拟结果相似，证明了数值模拟的可靠性。

(4）探地雷达在地下空洞发育区域和地下埋藏金属物区域的瞬时功率衰减特征各不相同。在遇到类似难以区分的异常时，可通过分析其瞬时功率衰减曲线，对真、假异常进行进一步筛选，去伪存真。

参考文献:

[1]谢昭晖,李金铭.探地雷达技术在道路路基病害探测中的应用[J].地质与勘探,2007,43(5):92.

[2]郭士礼,段建先,张建锋,等.探地雷达在城市道路塌陷隐患探测中的应用[J].地球物理学进展,2019,34(4):1609.

[3]邓世坤.探地雷达用于深圳湾软基处理工程中的质量检测[J].岩土工程学报,2000,22(1):136.

[4]胡群芳,郑泽昊,刘海,等.三维探地雷达在城市市政管线渗漏探测中的应用[J].同济大学学报(自然科学版),2020,48(7):972.

[5]余志雄,薛桂玉,周创兵.复信号分析技术及其在地质雷达数字处理中的应用[J].岩石力学与工程学报,2005,24(5):798.

[6]张胜,何文超,黎永索,等.基于小波变换时能密度法的隧道空洞充填物识别[J].煤炭学报,2019,44(11):3504.

[7]凌同华,张胜,李升冉.地质雷达隧道超前地质预报检测信号的HHT分析法[J].岩石力学与工程学报,2012,31(7):1422.

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

[9]乔尔.探地雷达理论与应用[M].北京:电子工业出版社,2011.

[10]张劲松,丛鑫,杨伯钢,等.地下管线探测雷达图特征分析[J].地球物理学进展,2019,34(3):1244.

[12]岳建华,何兵寿.超吸收边界条件在地质雷达剖面正演中的应用[J].中国矿业大学学报,1999(5):42.

[13]石明,冯德山,吴奇.CFS-PML在探地雷达FDTD法模拟中应用[J].大连理工大学学报,2015,55(6):559.

基金资助:国家重点研发计划(2017YFC1501203);长江科学院中央级公益性科研院所基本科研业务费(CKSF2019434 SL);中央级公益性科研院所基本科研业务费专项资金(CKSF2021463/YT);

文章来源:周黎明,张杨,付代光等.道路地下空洞探地雷达波场和时频特性[J].同济大学学报(自然科学版),2024,52(01):77-85.