道路结构型地层瑞利波相速度频散曲线的完整求取

杨天春¹，何继善²，鲁光银²，朱自强²

(1. 湖南科技大学 土木工程学院，湖南 湘潭，411201；

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

摘 要：

何完整地求取道路结构型地层的瑞利波频散曲线。在以往研究工作的基础上采用附加层法，即假定原模型的最底层介质厚度为有限，再在其下部添加一个与表层介质物性相同的半无限空间，可使原复数域的求解问题转化到实数域进行，则利用瑞利波频散方程，用二分法可得原模型的多模式相速度频散曲线。在采用附加层法进行计算时，一般假定的层厚度较大，会出现高频数值溢出现象。为此，利用二分法的特点，当计算过程中数值趋于较大时，采用较小的有限值进行替代，并且保证不改变原计算结果的正负特性，避免传递矩阵法计算中高频数值的溢出。对2个假定的三层介质模型进行模拟计算，结果表明原模型瑞利导波曲线与附加层法的计算结果完全吻合，同时附加层法计算出了高频区域泄漏波模式的相速度频散曲线，从而验证了附加层法、高频溢出处理方法的正确性和可行性。

关键词：

瑞利波；频散；道路；数值溢出；附加层法；

中图分类号：P631.414          文献标志码：A         文章编号：1672-7207(2013)02-0642-07

Holonomic calculation of Rayleigh dispersion curves for pavement systems

YANG Tianchun¹, HE Jishan², LU Guangyin², ZHU Ziqiang²

(1. College of Civil Engineering, Hunan University of Science and Technology, Xiangtan 411201, China;

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

Abstract: How to obtain the whole dispersion curves of Rayleigh waves for pavement systems in real domain was discussed. A new method, namely, appending layer method based on the past research was proposed. The last layer’s thickness of original model was supposed finite, and an infinite half space was added with the original model, and the adding layer’s characteristics were the same as the surface layer. Then, a complex domain problem was translated into a real domain problem. So, multimode dispersion curves could be obtained by Rayleigh waves’ equation and bisection method. The last layer’s thickness of original model was chosen to be very big when appending layer method was used. So, numerical overflow in high frequencies would appear. The characteristics of bisection method to avoid numerical overflow in transmission matrix approach were utilized. When an intermediate calculating value was very big, a small finite value would substitute it, and the sign of original calculating results would not be altered. The Rayleigh dispersion curves of two models were calculated. According to the modeling results, the guided wave curves of original model are consistent with the appending layer methods’. And phase velocity curves of Rayleigh leaky mode waves are also calculated in high frequency by appending layer method. All of these prove that appending layer method and disposition of numerical overflow are correct and feasible.

Key words: Rayleigh wave; dispersion; pavement; numerical overflow; appending layer method

近年来，随着我国基础设施建设投入的增加，高速公路建设已成为我国目前工程建设的重点之一。如湖南省“十一五”新开工53个高速公路项目，总计4 457km；“十二五”期间，全省共安排高速公路建设项目70个，建设总规模6 334km，估计投资3 000亿元。随着高速公路建设规模的扩大，道路建设及今后的维护急需经济、高效的无损检测技术，探地雷达技术作为一种地球物理新方法已被广泛应用于高速公路建设中^[1-2]，但这种电磁方法是以介质介电常数的差异为基础，无法检测道路结构的强度等力学性质。回弹法、拔出法、射钉法等用于砼强度无损检测的方法只能测定混凝土表面1~3cm深度的强度。所以，国内外学者近年来开展了采用瑞利波法对道路结构进行检测的技术研究^[3-4]。最近，爱尔兰UCD(University College Dublin)在欧盟框架计划FP7的支持下，也开始研究瑞利面波法探测道路结构的正、反演研究。瑞利面波法具有操作简便、无损、经济等特点，目前已被广泛应用于岩土工程领域，并引起了广大科研工作者对它的深入研究^[5-6]。通过对瑞利波频散曲线的反演，可划分地层结构并获得各地层的物理力学参数。但任何反演方法都是以正演为基础的，如果不能获得地层的瑞利波正演曲线，真正的反演就无从谈起。对于层状介质而言，其瑞利波频散曲线正演计算方法主要有Thomson-Haskell法、δ矩阵法、Schwab-Knopff法和Abo-Zena法等，这些方法都是从瑞利波频散方程出发，利用矩阵递推的方法计算频散曲线^[7-12]。对于道路结构模型而言，其地层是一个上硬下软的结构，即各层的横波速度是随深度增加的。如果直接采用上述矩阵递推的方法进行计算，则只能获得单条相速度频散曲线，频散曲线不完整，这是因为在采用这些方法求解时，求根的速度大小被限定在小于v_SN(即最底层介质的横波速度)内。为此，国内外学者一般采用垂直柔性系数法来计算，但该方法需计算瑞利波各模式在地表的垂直位移大小，求得的频散曲线实际上为合并化的频散曲线，并未求得多模式频散曲线的全解，且求解是在复数域进行，因此计算复杂且速度慢^[3]。在此，本文作者将以往提出的附加层法与避开高频数值溢出的方法相结合，研究快速求取道路结构型地层完整相速度频散曲线的方法。

1  理论计算公式

对于由N层均匀各向同性的弹性固体介质组成的层状半空间，引入柱坐标系(r,θ,z)，假定z轴垂直向下，z=0为地表自由界面，第i层介质位于平面z=z_i-1和z=z_i之间。在只考虑P-SV波的情形下，根据Zhang等^[11-12]曾采用的传递矩阵理论，可知沿层状介质传播的瑞利导波频散方程为

                (1)

其中：是由传递矩阵推出的第i层介质列矩阵E的第j个元素。并且

     (2)

即式(1)可写为，式(2)中的上标T表示矩阵转置，是与第i层介质参数相关的4×4阶矩阵的第m行、第n列元素。

对于第N层半无限介质，D^(N)为：

(3)

式中：ρ为第N层介质的密度；γ，γ_P和γ_S为与第N层介质的P波速度v_P、S波速度v_S和瑞利波速度v_R相关的参数。式(2)中的矩阵E满足如下递推公式：

，(i=2，3，…，N)       (4)

其中：传递矩阵X可分解为：

            (5)

式中：；；；A₁=-2γ_P；A₂=2γ_S；A₃=-2ργγ_P；A₄=2ργ_S(γ-1)；A₅=2ργ_P(γ-1)；A₆=-2ργγ_S；A₇=-2ρ²γγ_P(γ-1)；A₈= 2ρ²γγ_S(γ-1)；a₁=1-γ_Pγ_S；b₁=1+γ_Pγ_S；a₂=ρ(γ-1-γγ_Pγ_S)；b₂=ρ(γ-1+γγ_Pγ_S)；a₃=ρ²[(γ-1)²-γ²γ_Pγ_S]；b₃=ρ²[(γ-1)²+ γ²γ_Pγ_S]；d=ργ_S；e=ργ_P。在第i层介质中，矩阵中的，，h_i为第i层介质的厚度。

在以上各式中，各介质层的水平波数k=ω/v_R，，，。根据上述各计算式，在Matlab环境下利用二分法编程计算，即可获得层状介质瑞利导波的各模式频散曲线。

2  求解方法分析

2.1  问题的提出

表1所示为道路结构模型的物理参数，表1中v_R为用该层介质组成无限半空间时的瑞利波速度，利用前面传递矩阵计算公式，即可获得其频散曲线(见图1)。在利用二分法进行计算时，频率步长为0.05Hz，追索根的速度步长为11.1mm/s。

由图1可知：模型1的瑞利导波频散曲线只有单个模式，随着频率的增大，该模式的频散曲线逐渐递增，最终趋近于最底层介质(即路基)的横波速度122m/s，在上述设定的计算步长下导波频散曲线截止于11Hz处。众所周知，瑞利波法的勘探深度随着频率的增大而逐渐变浅，实测频散曲线应该是随着频率的增大最终趋近于表层介质(即沥青混凝土层)的瑞利波速度1 464.09m/s。显然，直接应用频散公式(1)求取道路结构频散曲线存在求解不全的问题。究其原因，这是由于瑞利波除了存在导波外，还存在一种泄漏模式波，后者在传播过程中能量不断损耗，衰减系数较大，传播距离很小，对于长距离的传播，一般将泄漏模式的影响忽略。在探测道路结构这种浅层介质时，偏移距、炮检距或道间距一般很小，并用小锤作振源以激起较强的高频信号，这时需考虑泄露模式波的存在。

图1  模型1的导波频散曲线

Fig.1  Guided-wave curve of model 1

2.2  附加层法

从数学意义上看，泄漏模式波对应频散方程的复数解，解的实部表示传播速度，虚部表示衰减大小。二分法只能对实数方程进行求解，为此，Hunaidi^[3]采用求解刚度矩阵的办法来获得泄漏模式波的解，由于计算结果中多模式的存在，又通过计算地表垂直柔性系数来判定实际勘探工作中所接收到的瑞利波模式(即垂直柔性系数最大的瑞利波模式)。在国内，夏唐代等^[13]采用有限单元法，通过求解波数k的二次特征值问题获得道路结构的瑞利波频散曲线。这2种方法都是在复数域求解，求解过程费时费力，非常繁琐。为此，根据以往提出的附加层法对道路结构模型做如下变化(见图2)，即假定路基的厚度有限，并在原模型的底部添加一层与表层介质相同的半无限空间层^[14]。

就瑞利波频散曲线而言，低频反映的是深部介质特性，高频反映的是浅部介质特性。因为瑞利波的波长不同，穿透深度也不同；当频率很低时，瑞利波的影响深度达到附加层时，则附加层法在0Hz附近的计算结果必然会出现一定的误差。如图2(b)所示，当路基的厚度h₃取得足够大时，除0Hz附近的频散曲线外，图2(a)与图2(b)所示的2个模型的频散曲线应该相同。同时，在道路结构的实际检测工作中，由于目标体在近地表且厚度较小，0Hz附近的频散曲线对实测工作无太多的意义。所以，可通过求解图2(b)的瑞利波频散曲线可获得原道路结构的频散曲线。对于图2(b)所示的模型来说，由于附加层选取的参数为原结构模型中横波速度最大的结构层(即面层)，在用式(1)进行求解时，求根的速度范围就包含了所有结构层；同时，求根是在实数域进行，可加快求解速度。

表1  道路结构模型的物理参数

Table 1  Material properties of simulated pavement site

图2  原道路结构模型与修改后的道路结构模型

Fig.2  Pavement structures and modified structures

2.3  防止高频数值溢出

当采用附加层法求解道路结构的相速度频散曲线时，由于式(5)中矩阵包含P_i和Q_i项，若路基厚度h₃选取较大时，很容易出现数值溢出现象。本文作者曾用附加层法计算模型1的合并化频散曲线时，是采取随频率的增大，计算的速度范围逐渐变小，即出现溢出的区域不予计算，所以，实际上是没有将相速度频散曲线完全计算出来^[14]。

为解决瑞利波频散曲线计算中高频数字丢失问题，Knopoff^[9]将δ矩阵法应用于瑞利波频散方程的推导；1979年，Abo-Zena^[10]通过一系列四阶反对称矩阵的循环计算得到了瑞利面波的频散方程，在一定程度上避免了高频有效数字丢失。随后，许多学者也在此方面进行了研究，但他们解决问题的思路都是从数学推导上着手，改变公式形式，解决溢出的频率范围一般比较有限。

在此，本文作者延用以前提出的方法——通过计算机编程来解决频散曲线计算中的数值溢出问题^[15]。如图3所示，对于y=f(x)的函数，在用二分法判断(x₁，x₂)区间内是否有根时，是依据f(x₁)和f(x₂)的符号，而不是其大小。当f(x₁)与f(x₂)的符号相反时，则存在根，否则(x₁，x₂)区间内无根。若在计算过程中，当f(x₁)或f(x₂)快要溢出时，如果令其值为“-1”或“+1”，保持其符号不变，则不影响求根的结果。

图3  二分法示意图

Fig.3  Schematic diagram of bisect-method

前面的计算公式中有和项，当厚度h或频率f较大时，它们很容易产生溢出。通过对频散函数的计算过程进行认真分析后可知：频散函数的正负特性主是取决于含有P_i和Q_i的数值项。因此，在编程计算时，主要对2个指数项P_i和Q_i进行修正。

具体计算时，理论上求根的速度范围为(0，v_SN)，v_SN是最底层介质的横波速度；因而第N层半无限介质中D^(N)的各元素为实数，其各项元素与h和f无关，不必考虑溢出问题。

就第i层介质而言(i=1，2，…，N-1)，由于传递矩阵X中的矩阵与f有关,并且求根的速度大小v_R与层速度参数的关系比较复杂，可分3种情况讨论矩阵，并提出相应的修正方法。

(1) 当v_R＜v_Si＜v_Pi时，γ_P和γ_S均为实数，矩阵X的每一个分量均为实数。当和中的|γ_Pkh_i|＞10^2.5时，则令P_i=e^-i；若|γ_Skh_i|＞10^2.5，则令Q_i=e^-i；并同时令中的元素”1”为”0”。

(2) 当v_Si＜v_R＜v_Pi时，γ_P为实数，而γ_S变为纯虚数，可令γ_S=iγ′_S，此处γ′_S为实数。当|γ_Pkh_i|＞10^2.5时，则令P_i=e^-i；若|γ′_Skh_i|＞10^2.5，则令Q_i=e^-i；并同时令中的元素”1”为”0”。

(3) 当v_Si＜v_Pi＜v_R时，γ_P和γ_S均为纯虚数，可令γ_P=iγ′_P，γ_S=iγ′_S，此时γ′_P和γ′_S均为实数；若|γ′_Pkh_i|＞10^2.5，则令P_i=e^-i；若|γ′_Skh_i|＞10^2.5，则令Q_i=e^-i；并同时假定中的元素”1”为”0”。

通过上述修正，就能很好地解决瑞利波相速度频散曲线计算中的高频数值溢出问题。

3  数值模拟计算

根据上述求解思路，对表1中的模型在Matlab环境下进行编程计算。图4所示为导波频散曲线计算结果的对比。图4中的实线为模型1的原瑞利导波曲线，即图1中的频散曲线；图4中的“●”和“▲”分别为当h₃=100m或h₃=300m时，采用附加层法和防高频溢出方法后的计算结果。由图4可知：当h₃=100m时，在频率f＞4Hz之后，附加层法的计算结果与原模型的瑞利导波曲线几乎重合；随着h₃取值的增大，则附加层法的计算结果与原曲线的吻合程度逐渐增高，当h₃=300m时，则在f=1Hz时，二者的计算结果就几乎相同。可见，通过附加层法也可求取原道路结果的瑞利导波频散曲线。

图4  导波频散曲线计算结果的对比

Fig.4  Contrastive results of Guided-wave curves

由上面的分析可知：在实际的勘探工作中，h₃的选取与所用检波器的自然频率f₀直接相关。就模型1而言，如果检波器的自然频率f₀为4Hz，则h₃取100m即可，因为在频率f＞4Hz之后，附加层法的计算结果与原模型的瑞利导波曲线几乎重合。

为进一步验证方法的有效性，将计算的频率范围进一步扩大，图5所示为h₃=100m时，计算至60Hz时的频散曲线。为使图5中相速度曲线显示清晰一些，相邻两条相速度曲线分别用实线和虚线绘制；另外，为使基阶模式的相速度频散曲线在低频部分复合原模型的实际情况，图5中的基阶模式在频率f≤11Hz范围内的频散曲线用图1中的曲线进行了替换。

由图5可知：新的计算方法不仅能计算出原道路结构模型的导波频散曲线(即频散方程的实数解)，还求出了波速v_R＞122m/s(即原模型最底层介质的横波速度)范围内的频散曲线，也就是原方程复数解中的实部(即瑞利波速度)，求根的范围v_R＜1 570m/s，范围包含了所有层的瑞利波波速。

从图5还可看出：当频率f较大时，h₃取100m也能计算出相速度频散曲线，说明计算程序很好地避免了数值溢出现象。从理论上而言，频率f取无穷大时，仍能正确计算出泄漏模式波的相速度大小，只是随着频率的增大，各频率点对应的模式会逐渐增多，采用二分法进行搜根的范围也必须越来越小，使得程序计算的速度会越来越慢。

同理，对于表1中的道路模型2，若假定第3层介质的厚度h₃为50m，再在第3层介质底部添加一个与表层介质完全相同的半无限空间，则根据前面的频散方程即可获得如图6所示的频散曲线。图6中，基阶模式波曲线的低频段用原三层介质模型的导波曲线进行了替换，所以，图6中的频散曲线包括了模型2的导波频散曲线和泄漏模式波频散曲线，从而将模型2的面波频散特性完整反映出来。

图5  模型1修改后的频散曲线

Fig.5  Dispersion curves of modified model 1

图6  模型2修改后的频散曲线

Fig.6  Dispersion curves of modified model 2

4  结论

(1) 提出了求解道路结构型地基模型完整面波频散曲线的方法。

(2) 通过添加附加层的办法，将传递矩阵法搜索根的范围包含所有的介质层，从而将复数域的求解问题变换到实数域，相对于国外的垂直柔性系数法而言，简化了求根的过程，节省了计算机耗时，并计算出了原模型泄漏模式波的相速度大小。同时，就附加层法而言，只要原模型最底层介质的厚度选取足够大，则在低频区段，采用附加层法所得的基阶模式波频散曲线与原模型的导波频散曲线完全吻合。

(3) 采用附加层法时可能会出现高频数值溢出问题，利用二分法求根的特点，即只关心所求函数值的正、负符号，提出了在计算过程中函数值较大时，采用一个较小的有限值进行替代，避免了数值溢出，从而可求出高频范围内瑞利波相速度频散曲线。

(4) 在公路面层、基层质量检测中，由于探测的是近地表超浅层介质的物理特性，介质层的厚度很薄，所以，采用附加层法进行计算时，第3层(即路基)介质的厚度主要取决于检波器的频率大小，不必取得过大，这样可加快程序计算速度。

参考文献：

[1] LI Shucai, LI Shuchen, ZHANG Qingsong, et al. Predicting geological hazards during tunnel construction[J]. Journal of Rock Mechanics and Geotechnical Engineering, 2010, 2(3): 232-242.

[2] 杨天春, 吕绍林, 伍永贵. 地质雷达检测道路结构的理论及应用分析[J]. 中南工业大学学报: 自然科学版, 2001, 32(2): 118-121.

YANG Tianchun, L Shaolin, WU Yonggui. Analyzing theory and application of prospecting pavement structure by ground-penetrating radar(GPR)[J]. Journal of Central South University of Technology: Natural Science, 2001, 32(2): 118-121.

[3] Hunaidi O. Evolution-based genetic algorithms for analysis of non-destructive surface wave tests on pavements[J]. NDT & E International, 1998, 31(4): 273-280.

[4] 王洁, 陈汉良, 夏唐代. 列车行驶随机激励引起周围地面振动响应分析及其应用[J]. 振动与冲击, 2004, 23(1): 87-90.

WANG Jie, CHEN Hanliang, XIA Tangdai. Analysis of ground vibration under random excitation induced by passing train and its application[J]. Journal of vibration and Shock, 2004, 23(1): 87-90.

[5] 张大洲, 熊章强, 秦臻. 基于Fourier变换的瑞雷面波分离提取及实例分析[J]. 中南大学学报: 自然科学版, 2010, 41(2): 643-648.

ZHANG Dazhou, XIONG Zhangqiang, QIN Zhen. Separation and extraction of Rayleigh wave based on Fourier transform and case analysis[J]. Journal of Central South University: Science and Technology, 2010, 41(2): 643-648.

[6] 周竹生, 刘喜亮, 熊孝雨. 弹性介质中瑞雷面波有限差分正演模拟[J]. 地球物理学报, 2007, 50(2): 567-573.

ZHOU Zhusheng, LIU Xiliang. XIONG Xiaoyu. Finite- difference modeling of Rayleigh surface wave in elastic media[J]. Chinese Journal Geophysics, 2007, 50(2): 567-573.

[7] 凡友华, 肖柏勋, 刘家琦. 层状介质中轴对称柱面瑞利面波频散函数的计算[J]. 地震工程与工程振动, 2001, 21(3): 1-5.

FAN Youhua, XIAO Baixun, LIU Jiaqi. Computation of dispersion function of axis-symmetrical cylindrical Rayleigh wave in multi-layered media[J]. Earthquake Engineering and Engineering Vibration, 2001, 21(3): 1-5.

[8] MA Yanlu, WANG Rongjiang, ZHOU Huilan. A note on the equivalence of three major propagator algorithms for computational stability and efficiency[J]. Acta Seismologica Sinica, 2012, 25(1): 55-64.

[9] XU Yixian, XIA Jianghai, Richard D. Miller. Approximation to cutoffs of higher modes of Rayleigh Waves for a layered earth model[J]. Pure and Applied Geophysics, 2009, 166(3): 339-351.

[10] Abo-Zena A. Dispersion function computations for unlimited frequency values[J]. Geophys J Rastr Soc, 1979, 58: 91-105.

[11] ZHANG Bixing, YU Ming, LAN Congqing, et al. Elastic wave and excitation mechanism of surface waves in multilayered media[J]. J Acoust Soc Am, 1996, 100(6): 3527-3538.

[12] LU Laiyu, ZHANG Bixing. Inversion of Rayleigh wave using genetic algorithm in the existence of a low-velocity layer[J]. Acoustical Physics, 2006, 52(6): 701-712.

[13] 夏唐代, 胡永生, 杨顺群, 等. 道路结构瑞利波特性及动力响应分析[J]. 郑州工业大学学报, 2000, 21(1): 19-22.

XIA Tangdai, HU Yongsheng, YANG Shunqun, et al. Characteristics of R wave and calculating dynamic response in pavement structures[J]. Journal of Zhengzhou University of Technology, 2000, 21(1): 19-22.

[14] 杨天春, 易伟建, 何继善, 等. 求取道路结构型地层瑞利波频散曲线的方法[J]. 物探与化探, 2005, 29(5): 459-462.

YANG Tianchun, YI Weijian, HE Jishan, et al. A tentative method for calculation of Rayleigh dispersion curves for the foundation of pavement systems[J]. Geophysical & Geochemical Exploration, 2005, 29(5): 459-462.

[15] 杨天春, 吴燕清, 刘新华. 对瑞利波频散曲线计算中高频数值溢出的处理[J]. 煤炭学报, 2007, 32(10): 1041-1045.

YANG Tianchun, WU Yanqing, LIU Xinhua. Treating of numerical overflow in high frequencies for computing Rayleigh wave dispersion curves[J]. Journal of China Coal Society, 2007, 32(10): 1041-1045.

(编辑  杨幼平)

收稿日期：2012-02-25；修回日期：2012-05-07

基金项目：国家自然科学基金资助项目(41174061)；湖南省自然科学基金资助项目(12JJ3035)

通信作者：杨天春(1968-)，男，湖南津市人，博士后，副教授，从事地球物理与岩土工程研究；电话：13017155739；E-mail：chtyang@yahoo.cn

摘要：研究在实数域如何完整地求取道路结构型地层的瑞利波频散曲线。在以往研究工作的基础上采用附加层法，即假定原模型的最底层介质厚度为有限，再在其下部添加一个与表层介质物性相同的半无限空间，可使原复数域的求解问题转化到实数域进行，则利用瑞利波频散方程，用二分法可得原模型的多模式相速度频散曲线。在采用附加层法进行计算时，一般假定的层厚度较大，会出现高频数值溢出现象。为此，利用二分法的特点，当计算过程中数值趋于较大时，采用较小的有限值进行替代，并且保证不改变原计算结果的正负特性，避免传递矩阵法计算中高频数值的溢出。对2个假定的三层介质模型进行模拟计算，结果表明原模型瑞利导波曲线与附加层法的计算结果完全吻合，同时附加层法计算出了高频区域泄漏波模式的相速度频散曲线，从而验证了附加层法、高频溢出处理方法的正确性和可行性。