基于高分辨率Radon变换的瑞利波反演在役路基动模量

杨 博， 龙友明， 刘境奇， 汤跃文

(1.重庆交通大学 土木工程学院, 重庆 400074; 2. 浙江大学 地球科学院,杭州 310027; 3. 重庆交通大学 交通土建工程材料国家地方联合工程实验室, 重庆 400074)

动态回弹模量(以下简称: 模量)作为表征公路路基在车辆循环动荷载作用下抗变形能力的主要性能指标，其优劣程度直接影响道路承载能力和使用寿命[1-3]。然而，在役路基受上覆既有路面结构层的遮蔽，传统方法须破损开挖路面直至路基顶面才能进行测试，耗时、费力且操作繁琐。因此，探索精细化与智能化的在役路基模量无损感知理论、技术及评价方法，是未来公路施工质量控制和长期性能监测的发展要求和有待解决的重要科学问题之一。

瑞利波(Rayleigh wave, 以下简称: R波)勘探的理论研究从地球内部地壳和上地幔的构造研究(大尺度勘探)到石油地震勘探(中尺度勘探)再到岩土地基结构物勘探(小尺度勘探)以及物体内部特性变化的检查(微尺度调查)均有涉及且应用效果良好[4-6]。而在勘探过程中，运用合理的数学分析与信号处理方法提取精确、可靠的R波频散曲线是决定勘探成败的关键因素。R波常规提取方法主要有互功率谱法、相位移、f-k变换和τ-p变换方法等[7]，Rosyidi等[8]应用互功率谱法对路面结构层模量进行了反演，杨博等[9-10]应用互功率谱法和相位移法对刚柔式复合路面及新建路基的模量进行了反演。然而，受面层和基层等路面结构覆盖的在役路基工作区厚度往往约为1 m左右，并且，根据交通荷载受力特征及稳定性要求，路表从上至下各结构层的强度和刚度一般按逐层递减进行设计[11]，属于逆频散浅层结构[12]。这就对R波提取方法及其在频域的分辨精度提出了更高的要求。高分辨率线性Radon变换算法(linear radon transform, LRT)由Radon在1917年首次提出，属于一种利用预加权矩阵的共轭梯度算法[13]，Luo等[14]研究表明其能够有效压制结构模型中波场噪声的影响，利于频散能量的稀疏解的计算，从而显著提高R波频散特征的分辨率，现已逐渐被引入地震信号的分析与处理。为此，本文通过理论分析、数值仿真和对比试验，基于LRT算法探索R波反演在役路基动模量的理论与方法。

1 在役路基结构中R波理论频散特征

在役路基结构通常可视为轴对称条件下的N层弹性半空间，如图1所示。图1中：下标m为结构层的序号，m= 1, 2, 3, …,N；VPm为第m层结构材料中的纵波(P波)波速；VSm为第m层结构材料中的横波(S波)波速；ρm为第m层结构材料的密度。在此结构中产生的R波为P波与S波耦合而成一种沿自由表面传播的波导，相应位移矢量可表示为

(1)

式中：φ为P波位移势函数；ψ为S波位移势函数；k为波数，物理意义为单位距离内波个数，k=2πf/VR=2π/λR，其中，f为频率，VR为瑞雷波相速度，λR为波长；ez为z方向的空间单位向量。

图1 弹性层状体系Fig.1 Elastic layered system

根据动力Navier方程引入Helmholtz势分解，通过Fourier变换，则图1中由R波引起各结构层的位移与应力分量在ω-k域的解答[15]为

(2)

位移、应力与势函数的关系可按式(2)进一步用矩阵的形式表示为

Sm=Mm·Φm

(3)

(4)

同时，由于第N层介质无限深处不存在上行波，则有BN=0和DN=0。通过结合以上应力边界条件、层间接触条件和无限深处约束条件，通过系数矩阵Mm逐层传递，可建立表面应力位移、应力矢量S1(z0)和zN-1处势函数矢量ΦN(zN-1)的等式关系为

ΦN(zN-1)=K·S1(z0)

(5)

式中：K为系数传递矩阵；其余符号同前。

式(5)中若保证表面位移有非零解，则须满足

(6)

式(6)则为轴对称柱坐标下R波在弹性层状半空间结构中传播的理论频散方程。

根据材料的模量E与VS及VP之间的解析关系

(7)

式中：μ为泊松比；其余符号同前。则式(6)中R波理论频散方程可进一步表示为

VR=G*(f,E1,ρ1,h1,E2,ρ2,h2,…,E0,ρ0)

(8)

式中：E0和ρ0为最下方第N层路基的模量与密度；其余符号同前。

式(8)即为R波在在役路基结构中传播的理论频散特征显示方程。则当各层的E,ρ和h给定，按式(8)可计算得到一条R波相速度VR随频率f变化而变化的理论频散特征曲线。为此，本文根据JTG D50—2017《公路沥青路面设计规范》中规定两种基准路基路面结构，具体参数见表1。其中，根据JTG D50—2017《公路沥青路面设计规范》要求，基层须按其材料的模量乘以0.5的调整系数作为其结构层模量。据此，利用式(8)计算得到R波在基准路基路面结构中的理论频散曲线，如图2所示。

表1 基准路基路面结构的力学与几何参数参数Tab.1 Mechanical and geometric parameters in reference structure of subgrade and pavement

图2 基准路基路面结构R波理论频散曲线Fig.2 Theoretical dispersion characteristics of Rayleigh wave in reference structure of subgrade and pavement

计算结果表明，两种不同基层类型的基准路基路面结构的R波理论频散曲线在f<25 Hz的低频带，其相应R波相速度VR收敛于139.7 m/s，接近路基的VS(140.9 m/s)，而在f>7 000 Hz的高频带，其相应R波相速度VR收敛于1 050.4 m/s，接近于沥青混凝土面层的VS(1 142.9 m/s)，呈现良好的极限特性。f在25～7 000 Hz过渡带，两种不同基层类型基准路基路面结构的R波理论频散特征区别明显，此过渡频带相应无机结合料稳定类基层结构的VR比松散粒料类基层结构VR平均大142.8 m/s，最大相差321.0 m/s。对于无机结合料稳定类基层的沥青路面结构而言，无机结合料基层厚度是路面结构验算考虑的重要因素，因此，分别改变基准路基路面结构中路基模量E0与基层厚度h2，相应R波理论频散曲线如图3和图4所示。

图3 VR随基层厚度h变化的频散特征Fig.3 Dispersion characteristics of VR varying with base thicknesses

图4 VR随路基模量E0变化的频散特征曲线Fig.4 Dispersion characteristics of VR varying with modulus of subgrade

图3显示：当沥青混凝土面层E1和路基E0分别取8 000 MPa和100 MPa保持不变，而无机结合料基层厚度h2从0变化至60 cm过程中，各理论频散曲线在f<25 Hz的低频带和f>7 000 Hz的高频带相应VR均收敛于E0为100 MPa和E1为8 000 MPa对应横波波速VS0和VS1，而f在25～7 000 Hz的过渡频带向低频移动，最大平移幅值约2 400 Hz。图4显示：当路基模量E0分别取40 MPa，100 MPa，150 MPa和200 MPa而其余参数不变，R波频散曲线在f>7 000 Hz的高频带，VR均收敛于1 050.4 m/s，与VS1接近；而在f<25 Hz的低频带，VR分别收敛于88.3 m/s，139.7 m/s，170.0 m/s和194.4 m/s，与路基模量E0取40 MPa，100 MPa，150 MPa和200 MPa对应的VS0为89.1 m/s，140.9 m/s，172.5 m/s和199.2 m/s接近。以上表明R波的理论频散特征能直观反映在役路基结构层模量与厚度的差异。

2 基于LRT算法的R波频散曲线提取

2.1 路基路面弹性波场响应的数值仿真

考虑当前高等级公路普遍设置无机结合料稳定类底基层，本文在基准路基路面三层结构基础上，将在役路基结构模型增至四层，结合JTG D50—2017《公路沥青路面设计规范》中材料参数取值范围，假定模型的参数如表2所示。

表2 模型参数Tab.2 Parameters of the model

针对表2模型，本文按时间二阶和空间四阶的交错网格有限差分方法[16]进行弹性波场响应的数值仿真。其中，将水平径向和结构总厚度设置成6.0 m，四周和底面采用采用吸收边界条件[17-18]，以模拟水平无限远和路基无限深半空间工况，模型空间网格大小取为0.005 m，按文献[19]中所提出的两个单载正弦波叠加的半正弦波冲击荷载作为激励源函数，具体为

(9)

式中：s(t)为冲击荷载归一化半正弦函数；f0为激励主频；τ为延迟时间或不应期；h(t)为heaviside阶跃函数。

高阶交错网格有限差分计算过程中，为了使高频带能全面包含厚度最薄层信息，最短波长须小于最薄层厚度，则有

(10)

式中：λmin为最短波长；VSmin为所有结构层中VS的最小值；hmin为所有结构层中最薄层的厚度。

按表2中具体参数，按式(10)则f0不得小于705 Hz，同时，基于线弹性材料的模量不受加载频率变化的基本假定和图2、图3中基准路基路面结构中R波理论频散特征，本研究将上述四层在役路基结构模型的激励主频f0取为2 600 Hz，激励延迟τ取0.000 6 s。

图5 R波测试技术指标Fig.5 Testing technical indicators of Rayleigh wave

图6 波场仿真记录A(r, t)Fig.6 Wave field simulation record A(r, t)

R波测试在役路基结构过程中的技术指标见图5，为了避免近场效应，最小偏移距为[20]

(11)

式中：D为最小偏移距；ZN-1为最下层顶距自由表面的距离，μ0为最下层结构材料的泊松比。

结合表2参数可算得D为0.72 m。据此，本研究将偏移距取0.8 m。同时，根据奈奎斯特采样定理，要使R波在提取过程中在波数域不产生混叠，相邻信号接收点的间距即道间距d和最大波数kmax须满足

(12)

式中：kmax为最大波数；λmin为最小波长；d为道间距。

由式(12)可知，d必须小于最小波长λmin的一半，根据半波长理论[21]可知，R波最浅探测深度为λmin的一半，d不得超过最浅探测深度，而工程中作为顶层的沥青混凝土面层通常进一步分为上、中、下面结构分层，每分层通常厚度约为0.05 m，据此，d可取为0.05 m。这样，根据高阶交错网格有限差分计算结果，以道间距0.05 m为间隔，取距激励源水平径向0.8～1.95 m的24道路表加速度垂直分量时程，通过归一化后定义为波场仿真记录A(r,t)，如图6所示。

2.2 LRT提取方法及对比分析

将波场仿真记录A(r,t)按时域傅里叶变换得到A(r,f)，相应A(r,f)的LRT变换对为[22]

(13)

式中：p为慢度，表示为时间的倒数；m(p,f)为波场加速度在频域-慢度域的变换结果。

同时，式(13)可进一步用矩阵表示为

m=LTA(r,f)

(14)

式中：m为m(p,f)矩阵；A为A(r,f)矩阵；LT为LRT正变换算子矩阵的转置。

式(14)意味着若能将m确定，通过将p域取倒变为v域，即可确定A(r,t)在f-v的振幅谱，以此作为R波频散曲线的提取方法。然而，若直接采用式(14)计算，由于LT的分辨率很低，且容易造成波场能量延伸和假频现象。为此，可按式(15)构建目标函数[23]，以设计一个反演系统来求解m，按极值条件求导可将式(15)进一步改写为式(16)

(15)

(16)

(17)

式中：J为目标函数；Wm为模型权重矩阵；λ为阻尼因子；Wd为误差权重的对角矩阵，按式(17)确定。

联立式(16)和式(17)通过反复迭代Wd以不断优化m矩阵，从而获得高分辨的f-v振幅谱，以此提取R波频散曲线。将图6中波场仿真记录A(r,t)分别通过LRT算法和相位移法提取R波频散曲线，并与按式(8)计算得到的R波理论频散曲线进行对比，结果如图7(a)～图7(e)所示。

图7 R波频散曲线提取方法比对Fig.7 Comparison of Rayleigh wave dispersion curve for different extraction methods

由图7(a)～图7(f)结果可见，在采用同样道数加速度信号的条件下，LRT方法得到f-v域的振幅谱能量更为集中，辨析度更高。通过式(18)计算误差发现，采用24道信号相位移法、24道信号按式(14)算法、24道信号LRT算法迭代2次、12道信号LRT算法迭代5次和12道信号LRT算法迭代3次提取频散曲线与理论频散曲线的相对误差ε分别为：3.22%，5.98%，1.86%，2.19%和3.01%，见图7(g)，意味着与利用24道速度信号通过相位移法提取R波频散曲线的精度相比，仅用12道加速度信号通过LRT算法迭代3次提取频散曲线就能达到更高的精度，这样富余的12道加速度信号可以按滚动排列测试[24]可提取得到13条中站处的频散曲线，见图8，从而在已有测试条件不变的基础上，大幅提高工程中R波测试精度与效率。

(18)

图8 R波滚动排列测试Fig.8 Roll-along acquisition format

3 路基模量反演与对比测试

3.1 模量反演分析

运用高分辨率LRT算法迭代3次提取得到的R波频散曲线，以式(19)为目标函数，将各层密度ρ和厚度h作为已知条件，结合R波理论频散特征方程式(8)，采用非线性阻尼最小二乘Levenberg-Marquardt算法(以下简称: L-M算法)反演各结构层模量E，为了验证反演方法的有效性，本文将各结构层模量的初值取为模型值的一半，相应结果如图9(a)所示。

(19)

式中：Wobj为目标函数；其余符号同前。

结果表明，反演通过9次迭代，Wobj由329.2 m/s降至49.3 m/s，相应路基路面各层E的反演结果分别为：9 788.4 MPa，7 541.1 MPa，2 875.8 MPa和109.7 MPa，与表2中模型各层E的相对误差均在9.7%以下。通过将以上各层E反演结果作为已知约束条件，将路基顶面以下1.0 m范围内视为路基工作区，按0.1 m厚度细化分为10层，进一步反演路基工作区内10个分层的E，如图9(b)所示，从而提高R波感知路基E沿深度方向的分辨率。同时，结合图8沿水平径向滚动排列测试提取得到的多个中站处频散曲线反演模量的结果，通过kriging插值[25]可实现在役路基E的二维反演成像。

图9 基于R波频散特征的模量反演Fig.9 Modulus back-calculation base on dispersion characteristics of Rayleigh wave

3.2 对比测试

为了验证基于R波LRT提取方法以反演在役路基模量的可行性与合理性，本研究在重庆交通大学交通土建工程材料国家地方联合实验室的道路结构足尺综合试验路槽内开展现场对比测试。试槽长60 m×宽5 m×深2 m，路基填筑压实厚度为1.4 m。其中，以K0+030～K0+060为对比测试段，该段包括三种路面结构，如图10所示，路基填料为低液限黏土，土性参数见表3。路基压实施工过程中，距路基顶面向下0～0.3 m范围的压实控制标准为96%，0.3～1.0 m范围的压实控制标准为94%，1.0 m以下按的压实控制标准为90%。对比测试步骤为：①在路基顶面沿K0+030～K0+060中线的行车方向，每隔1.0 m布置测点，点对点依次进行贝克曼梁测试弯沉l、承载板测试静模量Eb、落锤式弯沉仪(falling weight deflectmeter,FWD)测试模量EFWD；②待路面各层结构摊铺碾压成型后，在路表通过DZQ-3A级联式工程地震仪进行R波测试。其中，检波器共计12道，按图11(a)使测点置于第6道和第7道检波器中心，另外，根据材料的异同，路面结构可简化为三层模型(见图10)，在综合考虑路基的条件下，可按四层模型反演测点处的路基模量ER；③在K0+032～K0+038长6 m的路表进行R波滚动排列测试，每次整体向前水平推进0.3 m，共滚动测试21次，将1.4 m厚路基细分为0.1 m厚的14个分层，以此基于LRT进行路基模量反演成像。为了保证R波测试能反映足够深度且具备良好的分辨率，本试验将偏移距和道间距分别定为0.6 m和0.1 m，R波激励采用主频f0为2 600 Hz的半正弦压波，激励幅值参考标准轴载取0.7 MPa，相应测试与反演的结果如图11(b)所示。

图10 路面结构的材料与厚度Fig.10 Materials and thicknesses of pavement structures

表3 土性基本物理参数Tab.3 Physical parameters of soil properties

图11 对比测试及结果Fig.11 Comparison tests and corresponding results

图11(b)中结果显示，R波反演路基模量ER、FWD测试路基模量和承载板测试路基模量Eb的均值分别为130.3 MPa，101.2 MPa和46.2 MPa，ER分别为EFWD和Eb的1.29倍和2.82倍。同时，ER与EFWD和Eb成正比，而与贝克曼梁弯沉l成反比。为此，对所有测试结果进行回归分析，得到关系曲线如图12所示，回归结果见式(20) ～式 (22)。

(20)

(21)

(22)

式中：n为测试样本个数；R为相关系数。

从式(17)～式(19)可以看出，ER与其他测试结果均呈良好的幂函数关系，所有相关系数均超过0.80。表明R波反演路基模量值与其他各测试指标具有良好的相关性。

通过R波在K0+032～K0+038长6 m的路表滚动排列测试后，通过LRT提取21根中站对应的频散曲线，将其按四层结构反演的路面各层模量并作为已知约束条件，将1.4 m厚路基细分为0.1 m厚的14个分层，进一步反演21个中站处的路基VS剖面，然后按kriging插值进行二维成像，如图13(a)所示，最后结合式(7)，即可得到该段路基模量E的反演成像结果，如图13(b)所示。结果显示：K0+032～K0+038段路基在压实控制标准为96%的0～0.30 m深度范围内，VS从表面220 m/s沿深度向下降至200 m/s左右，相应E从表面200 MPa沿深度向下降至160 MPa左右；在压实控制标准为93%的0.30 m～1.0 m深度范围内，VS由200 m/s沿深度方向降至170 m/s左右，相应E从160 MPa沿深度方向降至120 MPa左右；在压实控制标准为90%的1.0 m以下深度范围内，VS从170 m/s沿深度方向降至140 m/s左右，相应E从120 MPa沿深度方向降至80 MPa左右。由此可见，E反演成像结果能较好地反映路基的压实工况。然而，为了保证反演计算效率，本研究反演成像过程中路基各分层的ρ按土样的ρdmax取，在后续研究中，拟将各分层的密度视为反演变量，进一步提高E的反演成像精度。以上研究表明，基于高分辨率LRT变换的R波测试在役路基模量进行反演和成像是合理、可行的。

图12 相关性分析Fig.12 Correlation analysis

图13 反演成像结果Fig.13 Results of inversion imaging

4 结 论

(1) 基于R波在弹性层状介质结构中的理论频散方程，揭示了R波在在役路基结构中的理论频散特征。特征表明：在f>7 000 Hz的高频带，其相应R波相速度VR收敛于沥青混凝土面层的VS1；在f<25 Hz低频带，其相应R波相速度VR收敛于路基的VS0；当基层厚度h2从0变化至60 cm过程中，VR由高频向低频移动，平移幅值约为1 903 Hz。

(2) 基于高阶交错网格有限差分方法对在役路基结构模型的R波响应进行了数值仿真。通过高分辨率LRT和相位移法提取R波频散曲线，并进行误差对比分析。结果显示，与利用24道振动信号通过相位移法提取频散曲线与理论频散曲线之间相对误差3.22%相比，仅利用12道振动信号通过高分辨率LRT算法迭代3次提取频散曲线与理论频散曲线的相对误差即能降至3.01%，具有较高的精度。

(3) 基于高分辨率LRT提取R波频散曲线，运用L-M算法反演了在役路基结构模型的各层模量。结果显示，各结构模量反演结果的相对误差均在9.7%以下。据此，将反演结果作为已知约束条件，进一步将路基工作区细化分层，联合kriging插值方法建立了路基工作区模量的反演成像方法。

(4) 通过足尺试验路槽内对比测试发现：R波测试路基模量ER与其他测试方法相应结果均呈良好的幂函数关系，所有相关系数均超过0.80；通过R波滚动排列测试，对路基模量E进行了二维反演成像，成像结果与路基压实工况具有较好地一致性。