基于相位差譜的非平穩(wěn)地震波合成及應(yīng)用

賈宏宇 ，藍(lán)先林 ，陳 航 ，鄭史雄 ，李 晰 ，張永水

（1.西南交通大學(xué)土木工程學(xué)院，四川 成都 610031；2.北京工業(yè)大學(xué)城市與工程安全減災(zāi)省部共建教育部重點實驗室，北京 100124；3.中鐵二院工程集團有限責(zé)任公司，四川 成都 610031；4.重慶交通大學(xué)土木工程學(xué)院，重慶400074）

大跨度結(jié)構(gòu)物的多點地震激勵分析能較為真實反映地震作用下各支撐處不同的地面運動，能更加真實地考慮實際地震對長大結(jié)構(gòu)物的作用，從而得到地震作用下更為真實的結(jié)構(gòu)動力響應(yīng)[1-4].但是，有限的實測地震波根本不可能滿足抗震分析的需要，同時也不可能在擬建橋梁的各橋墩處剛好就有實測的地震波，且這些地震波正好能為擬建橋梁的抗震分析所用，這使得人工合成地震波成為天然地震波必要的補充[5-6].

對人工地震波合成，許多學(xué)者提出了各種方法，如三角級數(shù)法[7-9]、ARMA 方法[10]、小波變換法[11]等，其中三角級數(shù)法應(yīng)用最為廣泛[12].三角級數(shù)法假設(shè)相位角在[0，2π]內(nèi)均勻分布，從而使得三角級數(shù)法只能依靠包絡(luò)函數(shù)來模擬地震動的強度非平穩(wěn)性[13]，縱所周知，人為選擇強度包絡(luò)函數(shù)具有任意性，故所生成的地震波表征的強度非平穩(wěn)性也具有任意性.地震動的非平穩(wěn)性除了強度非平穩(wěn)性以外，還具有頻率非平穩(wěn)性.但是，頻率非平穩(wěn)性在三角級數(shù)法中不能得到體現(xiàn)[14].基于上述原因，提出關(guān)于相位差譜的地震波合成方法，而且這種方法能很好地考慮地震動的強度非平穩(wěn)和頻率非平穩(wěn)性[15-16].此外，天然地震動又具有空間性和多維性.就目前關(guān)于相位差譜合成地震波的研究而言，幾乎沒有相關(guān)研究考慮地震動的多維性.為了能更好地模擬地震動的真實特性，本文基于相位差譜合成了多維多點非平穩(wěn)地震波，能彌補天然地震波的不足，也為抗震規(guī)范人工波合成提供理論參考.

本文基于隨機振動理論，模擬多維多點地震動場，基于三角級數(shù)法，推導(dǎo)了基于相位差譜合成非平穩(wěn)地震波的理論，并分析了相位差譜的分布規(guī)律，最后以某剛構(gòu)橋為依托，運用合成的非平穩(wěn)地震波對剛構(gòu)橋進(jìn)行了地震響應(yīng)分析，并與反應(yīng)譜法進(jìn)行對比，驗證其正確性和精度.

1 多維多點人工地震動模擬

1.1 多維多點地震動場

為了表示結(jié)構(gòu)m個地面支撐點處的地震激勵的相互影響，用加速度表示的互功率譜矩陣為

式中：矩陣主對角元素Skk（iω）表示任意點k處的自功率譜密函數(shù)；非主對角元素Skl（iω）表示任意兩點k、l的互功率譜密度函數(shù)；pkl（iω）表示相干函數(shù)；|pkl（iω）|表示相干函數(shù)的模，反映部分相干效應(yīng)；相干函數(shù)的指數(shù)部分- ωdkl/υapp表示相干函數(shù)的幅角，相位體現(xiàn)地震空間性中的行波效應(yīng)；dkl表示支撐點k和l之間的水平距離沿波的傳播方向的投影；υapp表示視波速，在實際地震中，視波速隨頻率變化，離散性較大.

在橋梁抗震分析中，取多個確定的視波速分別計算，其結(jié)果能包絡(luò)視波速離散性的影響，所以在考慮行波效應(yīng)的抗震分析中，取多個視波速是必要的，也是合理的.

為了研究地震動的多維性，將一維功率譜密度函數(shù)矩陣（式（1））中每一個元素按式（4）擴充為三維互功率譜譜密度函數(shù)矩陣.

式中：x、y和z分別表示兩水平向分量和豎直向分量.

各分量之間有一定比例關(guān)系，根據(jù)文獻(xiàn)[17]：水平分量各個自功率譜密度相同且相關(guān)，則可得到

豎直分量的互功率譜也可根據(jù)式（6）獲得.

將式（4）、（5）和（6）代入式（1），將只考慮一維（m×m）地震激勵的互功率譜密度函數(shù)矩陣擴展到了能同時考慮3 個方向地震激勵相互影響的三維（3m× 3m）的互功率譜矩陣.

1.2 三角級數(shù)法合成多維多點地震動

三維的地震加速度互功率密度函數(shù)矩陣是Herrmite 矩陣，而且是正定矩陣.根據(jù)式（7），對三維互功率譜矩陣進(jìn)行LDLT 分解，即

式中：P為3m×r的矩陣，r為矩陣S(ω)的秩.

P*為3m× 3m維矩陣，用lm'（iω）表示x、y和z方向中任意方向在P*矩陣中非主對角的元素，那么地震動的幅值am'和相位角θm'計算方法為

在合成第m′點任意方向dx、dy和dz的平穩(wěn)地震動時，本文考慮了第m′點任意方向地震動與另外任意點各方向的相關(guān)性，那么，第m′點dx或dy、dz方向平穩(wěn)地震動表示為

式中：am′b和 θm′b分 別為點m′與b在第k個頻率時的幅值和相位角；φbk為點m′與b在第k個頻率成分中的相位變化值；t為時間.

采用三角級數(shù)法合成人工地震波時，φbk為[0，2π]區(qū)間內(nèi)均勻分布的隨機變量.運用包絡(luò)函數(shù)模擬地震動非平穩(wěn)性中的強度非平穩(wěn)，但是強度包絡(luò)函數(shù)的選擇具有很大的任意性，也不能考慮地震動的頻率非平穩(wěn)性.真實地震動具有強度和頻率雙重非平穩(wěn)性，雙重非平穩(wěn)性與相位差有著密切的關(guān)系，因此，采用相位差譜合成人工地震波更能表征地震動的雙重非平穩(wěn)性.

2 地震動相位差譜

2.1 相位差譜和脈動相位差譜

Ohsaki[15]認(rèn)為實際地震動的相位譜是隨機變量，且具有一定相關(guān)性，并用相位差譜的形式來表示其相關(guān)性.相位差可表示為

式中：Δφ(f)為相位差，其定義域為[-2π，0]，若其值大于0，則減去2π 轉(zhuǎn)換至定義域內(nèi)；φk和φk+1為相鄰兩個頻率對應(yīng)的相位角；f為頻率.

文獻(xiàn)[18]對真實強震加速度記錄的相位差譜進(jìn)行研究，回歸得到均值函數(shù)Δ(f)為

式中：M為震級；R為震源距；R0為常數(shù)，其值取15；a1（f）、a2（f）和a3（f）為與頻率有關(guān)的回歸系數(shù)，其取值參考文獻(xiàn)[18].

ε(f)的標(biāo)準(zhǔn)差σε根據(jù)式（15）計算獲得.

式中：d1、d2和d3為回歸系數(shù)，當(dāng)處于巖石地基時，d1=0.089，d2= -1.124，d3= 0.316.

圖1給出了震距R= 100 km，不同震級時相位差譜均值變化規(guī)律.

圖1 相位差譜均值趨勢曲線Fig.1 Mean curve trend chart of phase difference spectrum

2.2 相位差譜的分布律

為了研究相位差譜的分布規(guī)律，對常用的實際地震波中的El Centro 波S00E 脈動相位差譜進(jìn)行頻數(shù)分析[18]，得出隨機變量εb滿足對數(shù)正態(tài)分布規(guī)律，即隨機變量εb取對數(shù)有 εc=lnεb，那么εc服從正態(tài)分布：

則 εb滿足式（17）形式的對數(shù)正態(tài)分布：

隨機變量εc與εb的均值、均方差有如下關(guān)系：

式（16）～（18）中：參數(shù)mc、σc分別為隨機變量εc的均值和均方差；mb、σb為εb的均值和均方差.

2.3 基于相位差譜的人工地震動合成

均值mb= 2π，σb可根據(jù)式（15）回歸得到，再根據(jù)式（18）獲得εb的均值和均方差.最后根據(jù)分布規(guī)律計算得到滿足相應(yīng)條件的隨機數(shù)εb，從而就得到脈動相位差ε（f），即

再將ε（f）與式（13）得到的 Δ(f)相加，并結(jié)合式（10），就合成了非平穩(wěn)人工地震動.

3 實際算例

為了說明本文方法的正確性和準(zhǔn)確性，采用傳統(tǒng)反應(yīng)譜法來進(jìn)行驗證工作.以某高墩剛構(gòu)橋（89 m +160 m + 89 m）為對象.該橋位處6 級震區(qū)，震源距200 km.首先，選擇對應(yīng)的設(shè)計反應(yīng)譜，特征周期Tg= 0.3 s，

抗震重要系數(shù)Ci= 1.7，場地系數(shù)Cs= 0.9，阻尼調(diào)整系數(shù)Cd= 1.0，加速度峰值A(chǔ)= 0.15g，0.20g，反應(yīng)譜最大值Smax= 2.25CiCsCdA，場地類型為I 類場地，轉(zhuǎn)化為功率譜并模擬地震動場（縱向、橫向和豎向地震波傳播速度分別為300、500、400 m/s）；然后基于相位差譜合成多維多點地震波（地震動持續(xù)時間為50 s），再迭代修正合成地震波；最后得到滿足精度的人工地震波.本橋兩橋墩底處（1# 墩和2# 墩）設(shè)計反應(yīng)譜的最大值分別為0.15g和0.20g，見圖2.

3.1 非平穩(wěn)地震動合成及修正

首先將設(shè)計反應(yīng)譜視為目標(biāo)反應(yīng)譜，且轉(zhuǎn)為功率譜，構(gòu)造具有多維性和地震動空間性的互功率譜矩陣，基于式（10）和式（20）合成多維多點非平穩(wěn)地震動時程.為了提高合成精度，對合成的地震動進(jìn)行修正.在頻域內(nèi)，對合成地震動的幅值進(jìn)行迭代修正.修正方法和步驟如下：

步驟1將合成地震動時程轉(zhuǎn)化為反應(yīng)譜，并與目標(biāo)反應(yīng)譜做商，即

圖2 連續(xù)剛構(gòu)橋示意（單位：m）Fig.2 Schematic diagram of continuous rigid frame bridge(unit：m)

式中：SaT（ω，ζ）為目標(biāo)反應(yīng)譜；Sa（ω，ζ）為合成地震動轉(zhuǎn)化的反應(yīng)譜；ζ為阻尼比.

步驟2精度不滿足要求，進(jìn)行幅值修正，再次合成地震動.

根據(jù)式（18）對幅值進(jìn)行修正.

步驟3轉(zhuǎn)到第1 步進(jìn)行反應(yīng)譜對比，如滿足精度，則跳出此步驟，不滿足就進(jìn)入第2 步.

經(jīng)上述流程，人工合成橋梁墩底處初始（未迭代修正）多維多點非平穩(wěn)地震動見圖3～5.同時，本文也給出初始合成的人工地震波轉(zhuǎn)化成的反應(yīng)譜與目標(biāo)反應(yīng)譜的對比，見圖6.

由圖3～6 可知：（1）從波形來看，基于相位差譜合成的多維多點非平穩(wěn)地震動符合實際地震動非平穩(wěn)性的3 個階段：增幅階段、平穩(wěn)階段和衰減階段；（2）從幅值來看，1# 墩底合成地震動幅值比2#墩地震動幅值大，這與兩墩底的目標(biāo)反應(yīng)譜峰值對應(yīng)，可從宏觀判斷地震波合成是否正確；（3）1# 墩和2# 墩底處初步合成地震動轉(zhuǎn)化的反應(yīng)譜和目標(biāo)反應(yīng)譜基本一致，但1# 墩處幅值最大值相差約40%，2# 墩處相差約57%，且擬合反應(yīng)譜尖峰較多，不平滑，未達(dá)到目標(biāo)擬合精度15%.為了提高精度，需進(jìn)一步迭代修正.根據(jù)式（18）對初次合成的人工非平穩(wěn)地震動時程進(jìn)行幅值修正，并再次對比擬合反應(yīng)譜和目標(biāo)反應(yīng)譜，其最終修正的人工非平穩(wěn)地震動曲線見圖7，最終擬合反應(yīng)譜與目標(biāo)反應(yīng)譜的對比見圖8.基于MATLAB 采用Duhamel 積分將新生成地震波轉(zhuǎn)換為擬合反應(yīng)譜，相關(guān)參數(shù)：阻尼比為0.05；積分時間步長為0.02 和周期為0～10 s.由于篇幅有限，此處僅列兩個橋墩處x方向人工合成的任意地震波對應(yīng)反應(yīng)譜和目標(biāo)反應(yīng)譜的對比.由圖7和圖8可知，經(jīng)過迭代合成的地震動同樣具有前述非平穩(wěn)特性，且在1# 墩和2# 墩處幅值最大值相差分別變?yōu)榱思s6%和13%，且擬合反應(yīng)譜尖峰較少，曲線平滑.同時，經(jīng)過兩個墩處地震動合成發(fā)現(xiàn)，目標(biāo)反應(yīng)譜峰值越大，擬合后峰值相差也越大，如2#墩處目標(biāo)反應(yīng)譜峰值0.2g，擬合反應(yīng)譜與目標(biāo)反應(yīng)譜差值比1# 墩處大.修正后1# 墩和2# 墩處的地震動幅值縮減最大值分別為0.745 m/s2和1.580 m/s2.

圖3 x 方向非平穩(wěn)地震動（未修正）Fig.3 Non-stationary ground motion at x direction（uncorrected）

圖4 y 方向非平穩(wěn)地震動時程（未修正）Fig.4 Non-stationary ground motion at y direction（uncorrected）

圖5 z 方向非平穩(wěn)地震動（未修正）Fig.5 Non-stationary ground motion at z direction（uncorrected）

圖6 擬合反應(yīng)譜和設(shè)計反應(yīng)譜（未修正）Fig.6 Fitting response spectrum and design response spectrum（uncorrected）

圖7 x 方向的非平穩(wěn)人工地震時程（修正）Fig.7 Non-stationary ground motion at x direction（corrected）

圖8 擬合反應(yīng)譜和目標(biāo)設(shè)計反應(yīng)譜（修正）Fig.8 Non-stationary fitting response spectrum and target design response spectrum（corrected）

3.2 橋梁地震響應(yīng)分析

為了進(jìn)一步研究修正前與修正后地震波對橋梁結(jié)構(gòu)動力響應(yīng)的影響，也為驗證人工合成地震波的正確性.在圖2剛構(gòu)橋?qū)?yīng)的有限元模型中，在1#墩和2# 墩底基礎(chǔ)處分別施加合成的非平穩(wěn)人工地震動，對橋梁結(jié)構(gòu)進(jìn)行多點激勵分析.圖9是基于本文合成地震動的時程分析法和傳統(tǒng)反應(yīng)譜法的工況說明圖例.圖9中列出了5 種工況：工況1，1# 墩和2# 墩峰值加速度都為0.15g時程分析法（一致激勵）；工況2，1# 墩和2# 墩峰值加速度分別為0.15g和0.2g時程分析法（多點激勵）；工況3，1# 墩和2# 墩峰值加速度都為0.2g時程分析法（一致激勵）；工況4，1# 墩和2# 墩峰值加速度都為0.15g反應(yīng)譜法（一致激勵）；工況5，1# 墩和2# 墩峰值加速度都為0.2g反應(yīng)譜法（一致激勵）.圖中：UE 為一致激勵；PDS為相位差譜法；ME 為多點激勵；RS 為反應(yīng)譜法.

圖10～12 給出了各工況荷載作用下，剛構(gòu)橋墩底剪力和彎矩的對比情況，由圖10可知，加速度峰值越大，墩底的剪力和彎矩值也越大，初步判斷地震作用下結(jié)構(gòu)響應(yīng)規(guī)律符合實際情況.采用反應(yīng)譜法比基于相位差譜合成地震動的時程方法（本文方法）計算得到的剪力和彎矩要大，當(dāng)加速度為0.2g時，1# 墩的剪力最大相差約3 973.5 kN，2# 墩的彎矩值最大相差889 146.2 kN·m；2# 墩比1# 墩剪力值小，但2# 墩比1# 墩彎矩值大，其中采用反應(yīng)譜法計算的剪力值相差在10%以內(nèi)，彎矩值相差在17%以內(nèi)，采用本文時程分析法計算的剪力值相差在9%，彎矩相差在12%以內(nèi).這說明了本文方法和傳統(tǒng)反應(yīng)譜法計算結(jié)果基本一致，也說明了本文方法正確性.

圖9 連續(xù)剛構(gòu)計算工況圖例Fig.9 Analysis case of continuous rigid frame bridge

圖10 一致激勵與多點激勵的墩底剪力彎矩（相位差譜法）Fig.10 Shearing and moment at the bottom of the piers under UE and ME（phase difference spectrum method）

圖11 墩底處剪力Fig.11 Shearing of the piers bottom

圖12 墩底處彎矩Fig.12 Moment of the piers bottom

4 結(jié) 論

本文基于相位差譜合成了多維多點非平穩(wěn)地震動，采用頻域方法對合成地震動進(jìn)行迭代修正.并研究了人工非平穩(wěn)地震動對剛構(gòu)橋多點激勵下的響應(yīng)情況.得到如下結(jié)論：

（1）基于相位差譜合成的多維多點非平穩(wěn)地震動不用包絡(luò)函數(shù)控制地震波的波形，較好地排除了人工選擇強度包絡(luò)函數(shù)的任意性，優(yōu)于傳統(tǒng)的三角級數(shù)法；

（2）基于相位差譜合成多維多點地震動不但考慮了地震動的空間相關(guān)性，還考慮了地震動的多維性，更加符合真實地震動特征；

（3）本文方法與傳統(tǒng)反應(yīng)譜法計算結(jié)果基本一致，采用反應(yīng)譜法計算的1# 和2# 墩剪力值相差在10%以內(nèi)，彎矩值相差在17%以內(nèi)，而本文時程分析法分別在9%和12%以內(nèi).