Rayleigh波頻散曲線兩種計算方法的對比

秦波,曹艷梅,夏禾

(北京交通大學土木建筑工程學院,北京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.