兩種計算程序對弱非線性硬場地地震響應的計算比較*

于嘯波，孫 銳，袁曉銘，陳龍偉

(中國地震局工程力學研究所地震工程與工程振動重點實驗室，黑龍江哈爾濱150080)

0 引言

土層地震動計算在地震工程中具有重要意義。常用的地震反應分析方法是一維等效線性化方法(齊文浩，2004)，LSSRLI-1與SHAKE2000是這種方法的代表程序 (廖振鵬，1989):前者代表了20世紀90年代的國際水平，是我國目前工程地震安全性評價工作規定使用的分析程序;后者最初由Schnable等 (1972)編寫，并陸續推出了91、2000等改進版本，代表著目前國際上的先進水平，這兩個程序都能在一定工況下得到可信的結果。按照我國抗震規范美國La Cienega臺陣場地屬于Ⅱ類場地，該臺陣記錄到了1999年M7.1地震從地下252 m處到地表的地震動，地下記錄輸入到SHAKE91中，計算得到的地表地震動與實際記錄較吻合 (Haddadi，Shakal，2006)。將位于Ⅱ類場地響嘡臺陣3號測井的地震記錄輸入到LSSRLI-1與SHAKE91中，地表反應都與記錄相近，這表明輸入地震動強度較低時，LSSRLI-1與SHAKE91可以基本滿足工程需要 (齊文浩，2004)。以上結果都是基于計算與實測的對比，從一個側面證明了兩程序的可靠性，但存在實測記錄少、工況數量不足的局限。

鑒于我國井下臺陣較少，一些學者在沒有強震記錄的情況下，也進行了一些兩程序間的對比分析，李亞東 (2011)對LSSRLI-1與SHAKE91進行比較發現，LSSRLI-1結果偏大，特別是在土層剛度小強震作用下更明顯。另外，LSSRLI-1在計算中還可能出現計算不收斂情況 (齊文浩，2008)。但是，目前對于SHAKE與LSSRLI-1的對比分析，還缺乏系統性，沒有進行深入的對比及分析差異原因。

場地條件對于等效線性化計算分析會有較大的影響 (景立平等，2005)。實際場地有不同分類方式，除去軟、硬夾層等特殊情況 (李恒等，2014)，大致可分為軟硬兩種類型 (巖土工程勘察規范，GB50021—2001)。在我國抗震規范中，場地分類與巖土勘察規范的意義有所不同，不僅考慮巖土體的剪切波速，還考慮了覆蓋層的厚度。Ⅰ、Ⅱ類場地土層具有相對高的剪切波速，覆蓋層厚度相對小，是抗震有利的場地。以往的研究結果顯示，在較硬場地上，等效線性化方法往往能夠有較好的表現 (齊文浩，2004;Haddadi，Shakal，2006)，但由于對比的工況較少，仍需進一步研究。理論上來講，等效線性化方法更適用于土的動力非線性較弱時，即土的模量衰減與阻尼比的增長曲線相對平緩，而硬場地往往具有這樣的動力特性。本文選擇較硬的Ⅰ、Ⅱ類場地為研究對象，并將這些工況簡稱為弱非線性硬場地。

本文以構造的2個Ⅰ類和4個Ⅱ類場地土層剖面為研究對象，輸入3種不同頻率類型的地震波，采用兩種等效線性化分析程序 LSSRLI-1與SHAKE計算土層地震動，以期掌握兩種典型土層地震反應分析方法的異同性。

1 計算理論與計算模型

1.1 水平成層土中波動方程的解

LSSRLI-1與SHAKE2000都是一維土層反應分析程序，以一維剪切波在水平成層介質中的地震反應穩態解為基礎編制的。

首先，底層輸入頻域地震位移為ENexp(ikNz)的地震波時，第n層地震波位移的頻域穩態一般解為

式中，kn=(1-idn)·ω/cn，dn為第n層的阻尼比;cn=(1+idn)，μn為第 n層剪切模量。En、Fn為波幅系數，需要用傳遞矩陣進行求解:

式中，αn=ρncn/ρn+1cn+1，為相鄰兩層的復波阻抗比值。(2)式右邊左乘的矩陣可以根據每一層的土層信息計算出來。根據 (2)式由E1、F1推出各層的En、Fn。但E1、F1開始是未知的，需要根據已知的波幅系數EN來得到，關系如下

式中，EN為已知的輸入地震波波幅系數;eN通過(2)式中的傳遞矩陣得到。由 (2)、(3)式可遞推出從1層到N層的波幅系數E、F。這樣可以通過 (1)式得到各層的位移，再進行傅式逆變換就得到了時域位移，進而可得速度、加速度。但是因為實際的土層是非線性系統，因此疊加原理失效，不能進行傅式變換，所以需要引入等效線性化的假設，讓系統成為可以使用傅式變換的線性系統 (廖振鵬，2003)。

1.2 計算模型

根據建筑抗震設計規范 (GB50011—2010)中場地分類的條款，構造了2個Ⅰ類、4個Ⅱ類場地的土層剖面。Ⅰ類場地模型因為比較薄，所以在程序中僅分為4個計算層，前3層為土層，第4層為基巖，土的物理力學特征見表1。Ⅱ類場地模型的計算層厚度一般為每層2～3 m，土層物理力學特征見表2。

表1 Ⅰ類場地土層的物理力學性質Tab.1 Physical and mechanical properties of class-I sites

表2 Ⅱ類場地土的物理力學性質Tab.2 Physical and mechanical properties of class-II sites

考慮到輸入地震波對兩程序計算結果 (PGA、反應譜和剪應變)差異程度有影響，輸入的地震波需要具有一定的代表性。本文選擇5%阻尼比反應譜中長周期成分占優、短周期成分占優以及成分比較均勻的地震波，確定輸入地震波為El-Centro地震波、AKTH19地震波和KSRH地震波(圖1)。AKTH19與KSRH均為日本Kik-net臺站代碼，輸入地震波信息見表3。從反應譜來看 (圖2)，KSRH地震波短周期成分豐富，AKTH19地震波長周期成分豐富，El-Centro地震波介于之間，以此代表頻率成分不均勻的地震波。將地震波都調幅至0.1 g、0.2 g與0.4 g進行輸入，分別代表Ⅶ、Ⅷ和Ⅸ度區。這樣，Ⅰ類場地的兩個剖面、3條地震波和3個幅值，共有18個組合工況，同樣Ⅱ類場地共有36個組合工況。

本文使用的非線性曲線是陳紅娟 (2009)提出的粘土的模量衰減曲線上限與阻尼比增長曲線的下限的組合，代表著土為弱非線性的情況，數值如表4所示。

表3 輸入地震波基本信息Tab.3 Basic information of the inputting waves

表4 土層非線性曲線Tab.4 Non-linear curves of profiles

圖1 輸入El Centro(a)、AKTH19(b)和KSRH(c)不同地震波時的波形Fig.1 Waveform of the inputting El Centro(a)，AKTH19(b)and KSRH(c)earthquakes

圖2 輸入El Centro(a)、AKTH19(b)和KSRH(c)不同地震波時的加速度反應譜 (5%阻尼比)Fig.2 Acceleration response spectrum of the inputting El Centro(a)，AKTH19(b)and KSRH(c)earthquakes(5%damping ratio)

2 結果比較

本文主要比較LSSRLI-1與SHAKE兩個程序計算出的地表加速度峰值、加速度反應譜和最大剪應變。相對差以SHAKE2000結果作為基準進行計算。

2.1 加速度峰值

結果表明，兩個程序計算出的Ⅰ類場地的地表加速度峰值PGA的相對差變化范圍是 (0～1.4)%，平均值為0.6%，差異很小。從圖3可見，加速度峰值的相對差隨埋深變化，而且靠近地表的時候相對差較大;剖面Ⅰ-2的相對差都比較小，沿深度變化不大。

兩個程序計算出的Ⅱ類場地的地表加速度峰值PGA的相對差變化范圍是 (0.03～30.7)%，平均值為6.6%，可以看出，與Ⅰ類場地相比，差異程度的變化范圍變大。由圖4可見，不同深度相對差有所不同，而且最大的相對差并不一定出現在地表處。

不同輸入地震波對PGA結果的影響見圖5。由圖可見，輸入3種地震波的計算結果具有很相近的分布，相對差都集中在10%以內。這表明，在弱非線性硬場地條件下，對于不同地震波，LSSRLI-1與SHAKE2000的大部分PGA計算結果是相近的。

2.2 反應譜

為了對比反應譜間的差異，得到兩個程序計算的地表加速度反應譜的譜值比為

其中，SS是SHAKE2000的反應譜結果 (g)，SL是LSSRLI-1的反應譜結果 (g)。

圖3 Ⅰ類場地加速度峰值相對差隨埋深的變化Fig.3 Relative difference of acceleration peak value varies with depth in class-I sites

圖4 Ⅱ類場地加速度峰值相對差隨埋深的變化Fig.4 Relative deviation of acceleration peak value varies with depth in class-II sites

圖5 不同輸入地震波下Ⅰ、Ⅱ類場地加速度峰值相對差統計直方圖Fig.5 Relative deviation histogram of acceleration peak in class-Ⅰ and class-Ⅱsites under different inputting waves

本文進一步采用平均譜值比來表示兩種程序計算的加速度反應譜的平均差異，按下式計算:

計算表明，Ⅰ類場地平均譜值比均值是0.3%，變化范圍是 (0.1～0.6)%;Ⅱ類場地的平均譜值比均值是2.3%，變化范圍是 (0.2～12.3)%。

圖6a是Ⅰ類場地計算結果。從各工況結果可見，在周期T=0.03 s附近，RS=0.6(相當于相對差為45.1%)，兩程序在這個周期附近的反應譜結果差異較大;0.1 s之后，反應譜結果差異很小。

圖6b是Ⅱ類場地計算結果。從各工況結果的曲線看，反應譜結果的相對差在3 s以前存在較大的波動，與Ⅰ類場地相比，反應譜的較大差異從短周期段發展到了長周期段。

譜比值反映的是反應譜間的相對差，不同輸入地震波下Ⅰ、Ⅱ類場地反應譜相對差統計結果見圖7。由圖可見，輸入3種地震波計算的結果具有相近的分布，而且不超過20%。這表明，在弱非線性硬場地條件下，兩個程序在輸入不同地震波時，反應譜計算結果是相近的。

2.3 最大剪應變

LSSRLI-1與SHAKE2000給出了每一計算層內的剪應變峰值，在不同埋深處，剪應變峰值也不同，筆者對剪應變峰值隨埋深變化曲線上的最大值進行比較。結果表明，Ⅰ類場地剪應變相對差變化范圍是 (0.7～9.0)%，均值為4.4%;Ⅱ類場地的剪應變相對差變化范圍是 (0.7～303.0)%，均值為46.4%。

圖6 Ⅰ類 (a)、Ⅱ類 (b)場地反應譜譜值比匯總Fig.6 Spectral ratio of response spectrum of class-Ⅰ(a)and class-Ⅱsites(b)

圖7 不同輸入地震波下Ⅰ、Ⅱ類場地反應譜相對差統計直方圖Fig.7 Relative deviation histogram of response spectrum in class-Ⅰ and class-Ⅱsites under different inputting waves

圖8為最大剪應變相對差的頻率分布直方圖，從圖中可以看出兩個程序計算出的最大剪應變在Ⅰ類場地上比較相近，相對差在10%內波動 (圖8a);而Ⅱ類場地的相對差變化范圍則較大，主要集中在50%以內，但最大可達300%(圖8b)。

按照不同的輸入地震波，將剪應變峰值相對差進行統計，結果見圖9。由圖可見，3種地震波造成的剪應變峰值差異程度具有相近的分布，但剪應變峰值相對差大于100%的情況均來自于KSRH地震波輸入。

2.4 剪應變相對差與PGA、反應譜相對差的相關性

經過對比分析發現，剪應變、反應譜和PGA相對差間有一定的關聯，其中一個較大時，其余兩個也較大。圖10a、b給出了一工況中剪應變差異與相應的反應譜差異，從圖中可以看出，SHAKE2000計算得到的最大剪應變為17.7×10-4，而LSSRLI-1為17.6×10-4，相對差為0.5%。SHAKE2000計算的PGA為0.28 g，LSSRLI-1也為0.28 g，結果相同。反應譜最大相對差為16%。圖10c、d為另一工況的計算結果，由圖可知，SHAKE2000計算得到的最大剪應變為83.1×10-4，而LSSRLI-1的為144×10-4，相對差為72.6%;SHAKE2000計算的PGA為0.81 g，而LSSRLI-1的為0.61 g，兩者之間的相對差為24%。反應譜最大相對差為59%。由此可見，當剪應變差異越大，則PGA和反應譜的差異越大，三者之間存在一定的關聯性。

2.5 差異原因初步分析

兩個程序的基本原理相同，輸入的參數相同，計算結果的差異極有可能是程序內部某些參數的計算方法或算法不同造成的。在迭代過程中，會變化的參數是根據等效剪應變來取值的土層模量比、阻尼比。等效剪應變計算都是剪應變峰值乘以折減系數，計算中該系數均為0.65，即等效剪應變計算方法相同。但從前文第2.3節可以看出，兩程序的剪應變峰值是不同的，所以兩程序迭代計算得到的等效剪應變也是不同的。因此，迭代中剪應變的不同計算方法是導致不同計算結果的原因。此外，迭代中模量比、阻尼比的插值方法也可能不同。有研究表明，土層地震反應結果對模量比與阻尼比很敏感(劉紅帥等，2005;孫銳等，2009)，這也能導致兩程序的計算結果有差異。

圖8 Ⅰ類 (a)、Ⅱ類 (b)場地最大剪應變相對差頻率分布直方圖Fig.8 Frequency distribution histogram of maximum shear strain deviation in class-Ⅰsite(a)and class-Ⅱsite(b)

圖9 不同輸入地震波下Ⅰ、Ⅱ類場地剪應變峰值相對差統計直方圖Fig.9 Relative deviation histogram of shear strain peak in class-Ⅰ and class-Ⅱsites under different inputting waves

為了表明兩程序的計算結果差異來自內部計算方法不同，筆者給出一個算例。采用圖10c、d的工況，即在剖面Ⅱ-1輸入0.4 g的El-Centro地震波。在這個工況中，SHAKE2000與LSSRLI-1的計算結果經過若干次迭代已經收斂，現稱圖10c、d中的結果為“收斂結果”。在兩程序已經收斂的結果上，用SHAKE2000再將其各自迭代一次，稱這次結果為“再迭代結果”。如果SHAKE2000與LSSRLI-1內部計算方法相同，那么LSSRLI-1的再迭代結果應該與LSSRLI-1的收斂結果相同或相近。實施這次“額外”迭代的具體做法如下:將圖10c中SHAKE2000與LSSRLI-1的各層剪應變計算結果提取出來，結合前文給出的土非線性曲線 (等效剪應變系數取0.65)，可以插值獲得對應兩個程序中各土層的剪切模量與阻尼比，見表5。將這些模量與阻尼比當作初始參數都輸入到SHAKE2000中進行1次迭代計算，LSSRLI-1的再迭代計算流程見圖11。因為只有1次迭代，所以剪應變計算結果只由各程序的剪應變計算方法決定，排除了多次迭代帶來的其他差異累積，同樣的還有反應譜的計算結果。剪應變、反應譜收斂、再迭代結果見圖12。

圖10 剖面Ⅱ-1輸入El-Centro地震波的剪應變及反應譜Fig.10 Shear strain and response spectrum of inputting El-Centro wave in profileⅡ-1

圖11 對LSSRLI-1收斂結果實施再迭代流程圖Fig.11 Flow chart for re-iterate on convergence results of LSSRLI-1

圖12 中的LS代表輸入表5中LSSRLI-1收斂結果的G、D列，再用SHAKE2000進行1次計算的再迭代結果，第一個L代表使用LSSRLI-1的收斂結果，第二個S代表用SHAKE2000進行計算;L代表圖10c、d中LSSRLI-1收斂結果;SS、S具有相似的意義，代表SHAKE2000相關結果。dS與dS1代表剪應變與反應譜中 SS與 S的差距，即SHAKE2000再迭代結果與其收斂結果的差距;dL與dL1代表各圖中LS與L的差距，即LSSRLI-1再迭代結果與其收斂結果的差距;dLS與dLS1代表各圖中L與S的差距，即圖10c、d中SHAKE2000與LSSRLI-1收斂結果的差距。

表5 從各層的剪應變計算對應的剪切模量G、阻尼比DTab.5 Shear modulus G and damping ratio D calculated by shear strain in every layers

圖12 不同迭代次數下剪應變結果 (a)和反應譜結果 (b)(剖面Ⅱ-1輸入El Centro地震波)Fig.12 Shear strain(a)and response spectrum results(b)under different iteration number(inputting El Centro wave in profileⅡ-1)

在圖12a中SS與S很接近，即dS很小，這說明SS僅僅是在S本已收斂的結果基礎上又迭代了一次;而LS與L有大的偏差，與dS相比，dL很大，這說明SHAKE2000程序與LSSRLI-1程序具有不同的剪應變計算方法。而從圖12b中可以看到，SS與S的反應譜很接近，dS1很小，而LS與L之間的距離則略大，dL1＞dS1。SS與S的反應譜也僅僅是多迭代一次而已，所以很接近;而LS與L反應譜有較小的差異，可能是由輸入參數上存在的誤差導致，也可能是反應譜計算方法存在不同，這些需要進一步研究。

從以上的對比可知，當剪應變不參與計算，即不依據剪應變來插值求得下一次計算所需各層土的剪切模量和阻尼比時，兩程序計算結果并無太大差異，如圖12b中L與LS兩條線已非常相近。但剪應變參與到計算當中時，兩方法所得到的地表反應就相差較大，如圖12b中S與L兩條反應譜。因而，可以初步推斷，造成兩程序間差異的主要原因是剪應變計算方法不同，當然，算法等方面也會有一些差異，但所帶來的影響并不大。

3 結論

本文討論了土為具有弱非線性的硬場地 (Ⅰ類、Ⅱ類場地)下采用LSSRLI-1與SHAKE2000計算土層地震動的異同性，主要結論為:

(1)Ⅰ類場地兩個程序計算出的PGA十分接近，相對差變化范圍不超過1.4%，平均值為0.6%;Ⅱ類場地兩程序間平均相對差為6.6%，可以忽略，但在某些情況下較大，達到30%。

(2)Ⅰ類場地兩個程序計算出的反應譜差別很小，譜值比均值是0.3%，變化范圍不超過0.6%;Ⅱ類場地的反應譜存在差異，平均譜值比均值是2.3%，變化范圍是0.2%～12.3%。

(3)兩個程序計算出的Ⅰ類場地反應譜中，差異出現在0.1 s前，0.1 s后差異很小;Ⅱ類場地反應譜的差異則從小于0.1 s延伸到0.5 s之后，且差異程度增大。

(4)就兩個程序計算出的最大剪應變而言，Ⅰ類場地相對差變化范圍是 (0.7～9.0)%，均值為4.4%;而Ⅱ類場地的相對差則有明顯增長，相對差變化范圍為 (0.7～300)%，均值為46%。

(5)PGA相對差和反應譜相對差與剪應變相對差之間存在一定相關性。

(6)LSSRLI-1與 SHAKE2000計算結果不同，一個主要原因在于剪應變計算方法的不同，其他的原因需要進一步對程序研究。

陳紅娟.2009.土動力非線性的變異性及其對地震動影響的概率分析［D］.哈爾濱:中國地震局工程力學研究所.

景立平，卓旭煬，王祥建.2005.復雜介質對地震波傳播的影響［J］.巖土工程學報，27(4):393-397.

李恒，張靜波，吳建超.2014.軟夾層和硬夾層對地表地震動特性的影響［J］.地震工程學報，36(3):441-444.

李亞東.2011.區域性土的動力特性及對地下結構抗震分析的影響［D］.廣州:廣州大學.

廖振鵬.2003.工程波動理論導論(第二版)［M］.北京:科學出版社.

廖振鵬.1989.地震小區劃理論與實踐［M］.北京:地震出版社.

劉紅帥，薄景山，吳兆營，等.2005.土體參數對地表加速度峰值和反應譜的影響［J］.地震研究，28(2):167-171.

齊文浩.2004.土層地震反應分析方法的比較研究［D］.哈爾濱:中國地震局工程力學研究所.

齊文浩.2008.土層非線性地震反應分析方法研究［D］.哈爾濱:中國地震局工程力學研究所.

孫銳，袁曉銘，劉曉鍵.2009.動剪切模量比與剪切波速對地震動影響及等量關系研究［J］.巖土工程學報，31(8):1267-1274.

Haddadi H.，Shakal A..2006.Analysis of some strong ground motion records obtained at CSMIP downhole arrays［C］//France:Third International Symposium on the effect of Surface Geology on Seismic Motion，Grenoble.

Schnabel P.B.，Lysmer J.，Seed H.B..1972.SHAKE:A computer program for earthquake response analysis of horizontally layered sites［R］//Report NO.EERC 72-12.Berkely:College of Engineering，University of California.

GB50011—2010，建筑抗震設計規范［S］.

GB50021—2001，巖土工程勘察規范［S］.