秦波,曹艳梅,夏禾

(北京交通大学土木建筑工程学院,北京100044)

1887年,Rayleigh首先发现了Rayleigh波的存在并揭示了其在弹性半空间介质中的传播特性。20世纪50年代初,人们又发现了Rayleigh波的频散特性,随之开始了利用天然地震记录中的Rayleigh波探测地球内部结构的研究。目前关于瑞雷波的研究主要集中在两大方面：一是瑞雷波的正演,即利用层介质参数,通过土体的传递矩阵来计算瑞雷波的频散特性(瑞雷波的相速度与频率之间的关系)[1-4];二是瑞雷波的反演,即利用现有的仪器测定瑞雷波的波速,提取瑞雷波的频散曲线,对地基土的参数等进行反演解释[5-6]。由于瑞雷波的频散关系是一个隐式关系,正演中关于理论频散曲线的推导不易求解,对此国内外的研究人员提出了一些用在正演中的研究方法,如直角坐标系下的Thomson-Haskell传递矩阵法[7],Knopoff算法[8]、Schwab算法[9],柱面坐标系下的δ矩阵算法和快速标量算法[1-3]以及薄层刚度法[4]等。本文主要对目前采用比较多的两种方法—快速标量法和薄层刚度法进行分析比较,通过算例分析给出每种方法的优缺点以及适用范围,为更好地进行Rayleigh波的正演和反演等研究提供合理的研究方法。

1 快速标量法

快速标量法是由张碧星[10]、凡友华[1-3]等在无量纲实数传递矩阵算法和快速矢量传递算法的基础上改进提出的。为降低有效数字的损失以及快速矢量传递算法的计算量,用标量的传递形式重新组织传递过程。

按下式定义五个标量x1,x2,x3,x4,x5

式中Vi—第i层介质参数。

式中 vSi、vPi—第 i层介质的横波波速和纵波波速;ν—Rayleigh波波速。

按照各层的介质参数从第n层到第1层对5个标量依次进行传递,每次传递的程序为

式中k和h分别表示Rayleigh波沿水平自由界面传播方向上的波数和层厚度,a=cos(γPkh),b= cos(γSkh),d=sin(γSkh)/γS,e=sin(γPkh)/γP,r =,s=,g=1-γ=v2/(2)。

得到第一层介质对应的5个标量,令此时的行x3为x3(1),则频散方程为

可以看出,快速标量算法通过将矩阵转化为标量运算,可以大大提高频散曲线的计算速度。但其不能进行复数运算,因此不能更好地考虑介质阻尼对模态的影响。

2 薄层刚度法

薄层刚度法是Lysmer和Wass在研究Love波时,将分层细化成很多薄层。Kausel和Roesset将此方法推广至分层介质中P-SV波场研究,将层刚度矩阵中的三角函数(若含有复数,则为超越函数)进行了化简,三角函数或超越函数运算最容易导致数字损失,从这一点讲,薄层刚度法并不比传递矩阵优越。然而当三角函数值很小时(若为复数,则要求复数的模很小),三角函数则可用简单的代数表示。

根据泰勒级数,无论 x是实数还是复数,当 x趋近于x0时

当x很小时,正余弦三角函数分别近似为

薄层刚度法层刚度矩阵中的三角函数一般表示为krαd和krβd,其中rα,rβ由下式给定求出

式中k—Rayleigh波沿水平自由界面传播方向上的波数;VRm、VPm,d—第m层Rayleigh波相速度、纵波波速和薄层的厚度。

由于rα,rβ的值有限,要使krαd,krβd很小,kd =2π d/λ必须很小,即波长相对分层厚度很大时,分层刚度矩阵可以表示为

式中 λ、μ—地基土的 Lamé常数;ρ—土体密度。

上式表明,将层厚划分为很多薄层,层厚相对波长很小时,这些薄层的刚度矩阵就可以用代数形式表示。将位移矢量以层排列方式转化为以水平向、竖直向排列方式,即将U=[u0,w0,u1,w1,…un-2,wn-2,un-1,wn-1]T改写为

式中u、w—分层界面处的水平位移分量和竖向位移分量;0—自由界面;i—第i层界面(i=1,2,3,…,n-1)。

薄层刚度矩阵相应改写为

将薄层刚度矩阵集成得到总刚度矩阵

从上式可以看出,二次特征值方程的求解可以转化为求解方程的一次特征值问题,求解方程的特征值就可以得到相应的Rayleigh波波数。

和直角坐标系下的算法(Thomson-Haskell的传递矩阵法、Knopoff改进算法)以及柱面坐标系下的算法(δ矩阵算法,快速标量法)方法不同,薄层刚度法采取了化简的形式,将三角函数进行了化简,提高了数字计算速度。用泰勒级数代替三角函数,并且根据每层的层刚度矩阵来组集总刚度矩阵,得到关于刚度矩阵的行列式,通过求解该行列式的特征值,得到波长值,他们和Rayleigh波的相速度是一一对应的关系,这种算法的优点就是求解速度快,方程为显式,关系对应明确。

虽然三角函数或超越函数容易导致数字损失,但当它用代数表示后,利用刚度矩阵可以直接求出矩阵的特征值,速度快精度高,并且可以通过复数求解来考虑阻尼效应。

3 算例

分别采用快速标量法和薄层刚度法计算分析对分层土介质模型Rayleigh波的频散性。

3.1 两层土介质模型

以中软场地土为例,地基土参数见表1。土层假设为两层,层厚度依次为2m,半无限空间(近似以100m计),计算频率范围400Hz。

分别采用快速标量法和薄层刚度法计算其频散曲线如图1所示。

表1 两层土介质模型参数Tab.1 Parameters of two-layer soil model

图1中列出了0至350Hz频率范围内的前六阶模式对应的频散曲线。从图中可以看出,对于Rayleigh波的前三阶模态,两种方法的计算结果除了在上一阶模态的截止频率附近计算差异较大外,在随后的频率范围内,二者并无差异。

3.2 三层土介质模型

第二种地基土模型选用中硬场地土,土层假设为三层,层厚度依次为2m,5m,半无限空间(近似以100m计),计算频率范围400Hz。

分别采用快速标量法和薄层刚度法计算其频散曲线如图2所示。

图2中列出了0至350Hz范围内三层土介质模型的频散曲线。可以观察出,在50Hz内,尤其对于基阶模式,快速标量法图形有一些不规则,但高阶模态和薄层刚度法差异不大。在高频区,各阶模态两种算法计算的频散曲线并无差别。

比较图1和图2,可以看出：

1)Rayleigh波相速度范围趋于层剪切波速最小最大值之间,两类方法计算的频散曲线结果比较接近,总体趋势一致,差异不是很大。

2)在50Hz之后,频率越大,两类方法计算结果越接近,尤其是基阶模态,两者的差异性最小。

3)相速度越小,两类方法计算结果越接近。

表2 三层土介质模型参数Tab.2 Parameters of three-layer soil model

4)采用快速标量法计算的Rayleigh波相速度比薄层刚度法计算的值大,但是随着频率的增大,二者的差异逐渐减小。

5)在频率较小时,尤其是在50Hz内,当相速度趋于层剪切波速最大值时,快速标量法的曲线不规则,并出现多条曲线。相反,薄层刚度法则趋于层剪切波速最大值。因此,在零频率附近薄层刚度法的计算效果优于快速标量法。

6)对于瑞雷波相速度最小值,两类方法差异性不大,均接近层剪切波速最小值。

4 结论

1)快速标量法计算频散曲线的速度比较快,并且它避免了一些不必要的虚数运算,大大提高了计算的稳定性。

2)与快速标量法相比,薄层刚度矩阵法可做复数运算,可以更好地考虑介质阻尼对模态传播特性影响,此外薄层刚度矩阵法可避免传递矩阵连乘导致的数值不稳定,但该方法要求所细化的薄层厚度要远小于波长,因此计算速度会受到一定的限制。

3)经计算发现在低频阶段薄层刚度法计算效果优于快速标量法,在高频区域二者差别不大。

4)与快速标量法相比,薄层刚度法进行简化后,可以利用特征值求出频率波速对应关系,计算简便。

[1]凡友华,刘家琦.层状介质中瑞雷面波的频散研究[J].哈尔滨工业大学学报,2001,33(5)：577-581.

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

[3]凡友华,肖柏勋,刘家琦.计算层状介质中轴对称柱面瑞利面波频散函数的δ矩阵法[J].物探与化探,2001, 25(2)：109-115.

[4]柴华友.弹性介质中表面波理论及其在岩土工程中应用[M].北京：科学出版社,2008.

[5]王赞文.瑞雷面波法正反演及应用研究[D].西安：长安大学,2005.

[6]吴燕青,杨天春.瑞利波频散曲线的反演[J].煤炭学报,2008,33(10)：1097-1101.

[7]HASKELL N A.The dispersion of surface waveson multilayered media[J].Bulletin of the Seismological Society of America,1953(43)：17-34.

[8]KNOPOFF L.A matrix method for elastic wave problem[J]. Bull Seism Soc Am,1964,(54)：431-438.

[9]SCHWAB F,KNOPOFF L.Surface wave dispersion computations[J].Bulletin of the Seismological Society of America, 1970,60(2)：321-344.

[10]张碧星,喻明,熊伟,兰从庆.层状介质中的声波场及面波研究[J].声学学报,1997,22(3)：230-241.