基于小波包方法的非平穩地震動模擬和設計譜擬合

李 亞 楠，　王 國 新

( 大連理工大學 建設工程學部, 遼寧 大連　116024 )

?

基于小波包方法的非平穩地震動模擬和設計譜擬合

李 亞 楠，王 國 新*

( 大連理工大學 建設工程學部, 遼寧 大連116024 )

摘要：提出了一種基于小波包分解和重構的模擬非平穩地震動并擬合設計譜的方法．該方法模擬出工程用地震動并調整頻率分量，使其以一定的精度擬合設計譜．首先，利用5個參數生成了構成地震動加速度時程的小波包系數矩陣，并通過小波包重構獲得加速度時程，該系數矩陣模擬了實際地震記錄中頻率隨持時逐漸減小的特點．其次，將模擬的地震動利用小波包分解方法分解為具有高分辨率非重疊的小波包系數矩陣，然后根據設計譜調整頻率分量．數值算例表明，該方法模擬的不同持時的加速度時程均能吻合同一設計譜，且迭代后仍然能保留地震動的非平穩特性．另外，該方法具有較穩定和較快的收斂過程，能在有限次調整迭代中實現較高的擬合精度．

關鍵詞：地震動模擬；譜兼容加速度時程；小波包；非平穩時程

0引言

在工程抗震設計中，反應譜技術自提出以來，一直是抗震設計規范中的重要內容之一[1-3]．而隨著設計方法和計算機技術的發展，人們已經不再滿足于抗震設計的反應譜靜力分析方法，而開始考慮全部地震過程，尤其對于重要工程結構，如核電站、大壩、大型儲物裝置、大跨橋梁、超高層建筑等．輸入結構的地震動時程的反應譜必須在一定精度下與設計加速度反應譜(以下簡稱設計譜)相擬合，也就是所謂的譜兼容時程(spectrum-compatible time history)．通常情況下，擬合方法是對實際強震記錄做修正，讓其滿足設計譜．然而，由于強震記錄的數量有限且分布不均勻，不能滿足多方面要求，地震動的模擬越來越受到重視．

一條加速度時程的反應譜是唯一的且是很容易計算的，但是由給定的設計譜來尋求加速度時程這一反問題求解卻是非常困難的．在過去的三十多年里，許多學者對此問題給出了不同的解法．Scanlan等在1974年提出了著名的三角級數模型來模擬具有給定功率譜的地震動加速度時程[4]，這一方法也是工程上應用最廣泛的方法之一．隨后，就開啟了利用傅里葉譜幅值和相位迭代的研究思路．傅里葉變換是非常直接的反應譜擬合方法，多數學者把重點放在了時域幅值非平穩的模擬上[5-7]．然而該方法存在著局限性：一是沒有很好的收斂特性；二是如果傅里葉譜修正得過多，則加速度時程的非平穩特性將發生較大的改變．之后，演化功率譜方法[8-9]被提出來反映模擬地震動的頻率特性，該方法假定相位和相位差遵循一定的概率分布．從此，頻率非平穩性的模擬受到了學者的重視．Mukherjee等[10]基于小波變換提出了一種直接對實際記錄進行頻率修正來擬合設計譜的方法，從而可以保留實際記錄的大部分頻率特性．但該方法精度不高且沒有給出擬合精度的分析結果．其他方法，如基于諧小波的方法[11]和希爾伯特黃變換法[12]也在保留地震動的瞬時頻率方面取得了較為理想的效果．國內對反應譜擬合技術的研究多采用基于相位譜和相位差譜的方法[13-16]．隨后，趙鳳新等提出了多阻尼反應譜擬合的時域疊加法[17]并引入到小波變換中[18]，其在精度上取得了較好的效果．

本文提出一種基于小波包變換的反應譜擬合方法，首先模擬具有時間和頻率非平穩性的地震動；然后采用小波包分解的方法將地震動加速度時程分解為有限頻率段；最后根據設計譜，對各個頻率段進行調整并迭代，從而達到較高的擬合精度．通過對一個核電站的設計譜進行擬合，分別得到兩條不同的地震動加速度時程，并且討論該方法的精度．

1基本原理

小波方法同時在時間域和頻域對信號做分解，從而將頻譜在時間軸上展開．小波包是小波方法的延伸，它不僅對低頻部分做分解，同時也對高頻部分繼續分解，提高了頻率的分辨率．

1.1小波包分解

信號的時間和頻率關系叫時頻圖，能夠表征信號這一關系的方法有很多，如短時傅里葉變換、離散與連續小波變換、小波包變換和偽韋格納分布．然而地震動模擬不僅需要對信號進行分解，同時需要能夠方便、完整地逆變換回時域，小波包分解與重構即符合這一要求．

小波包分解與重構的定義如下：

(1)

(2)

在給定的分解層數j下，第i尺度的第k個位置的小波包系數可以用下式計算：

(3)

(4)

(5)

(6)

圖1為1979年美國Imperial Valley地震中El Centro Array #12號臺站記錄的時頻圖．加速度時程位于正上方，單位為重力加速度g(g=9.81 m/s2)；左邊為傅里葉幅值譜．中間的黑白色方塊即小波包系數，顏色的深淺代表了系數的大小，圖中色度條的數值為系數絕對值的相對值．可以看出，小波包系數矩陣中數值大的部分同時對應地震波加速度時程的時域和頻域中幅值較大的區域．這說明，如果表達出該小波包系數矩陣的總體特征，就可以通過小波包重構的方法將地震動模擬出來，并且可以修正這些系數，從而達到擬合設計譜的目的．

圖1　時頻圖

1.2地震動模擬

從式(2)可以看出，如果得到了小波包系數矩陣，就可以模擬出一條地震動時程．而如果該小波包系數矩陣擁有與實際記錄相似的特征，則通過小波包重構得到的模擬時程同樣和實際記錄相似．本文首先通過對一條高斯白噪聲進行窗口化，來模擬地震動的時域非平穩性．所用的Saragoni窗函數[19]表達如下：

w(t)=a(t/Td)be-ct/Td

(7)

b=-θlnτ/(1+θ(lnτ-1))

(8)

c=b/θ

(9)

a=(e/θ)b

(10)

式中：w(t)是Saragoni窗函數，Td是模擬時程的持時，θ是窗函數的峰值發生點占整個高斯白噪聲時長的比例，且在高斯白噪聲結束時的幅值占峰值的比例為τ；a、b和c是中間參數．圖2展示了不同θ和τ下Saragoni窗函數歸一化后的形狀．

圖2　Saragoni窗函數在不同θ和τ時的形狀

然后，對窗口化的高斯白噪聲做小波包分解．這里采用的模擬時間間隔為0.005 s，分解層數j為9，所以分解后的小波包系數頻率分辨率為0.195 3 Hz，時間分辨率為2.560 0 s．此時小波包系數在時頻圖的頻率帶上是均勻分布的，需要對每一列系數做修正． 圖3(a)展示了美國1992年Big Bear地震中Seal Beach-Office Bldg臺站記錄到的一條加速度時程的時頻圖，可以看出其頻率在時域上是不斷變化的．在地震波開始階段，主要包含了高頻分量，隨著時間的進行和地震動幅值的減弱，頻率逐漸衰減．從圖1也可以看出該特性．通過對大量強震記錄的分析，發現大多數地震動均存在該現象．因此，在對每一列小波包系數進行調整時，本文采用了不同的歸一化函數．這里采用對數正態分布密度(logarithmic normal distribution density，LNDD)函數作為基準函數：

k=1,2,…,2N/2j

(11)

式中：Lk(·)是對第k列小波包系數進行歸一化的函數，f是頻率，μwk和σwk分別是Lk(·)的均值和方差．圖3(b)顯示了采用LNDD對每列小波包系數進行歸一化后的模擬結果，可以看出和實際記錄的時頻圖是十分相似的．

(a) 時頻圖

圖3時頻圖和模擬結果

Fig.3Time-frequency diagram and simulation results

在該調整過程中，μwk和σwk由下式計算：

當k

(12)

當k≥m時，

(13)

其中μwm和σwm是對第m列小波包系數歸一化的LNDD的均值和方差，這里，m=int (2N/2j)，int(·) 是取整函數．因此，第m列實際上就是位于加速度時程持時中間的一列．在該列之前，頻率應增大；在該列之后，頻率應逐漸減小．而p是頻率變化率參數．經過多次與實際記錄對比，為了避免過高和過低的頻率，p對于模擬不同持時的地震動時程應有不同的取值：

(14)

這樣，就得到了歸一化后的小波包系數．利用式(2)進行小波包重建，就可以獲得模擬的加速度時程x(t)．

1.3設計譜擬合

(15)

(16)

可以看出，在該步中增加了頻率的分辨率．根據測不準原理[20]，小波包分解中頻率分辨率和時間分辨率不可能同時無限擴大，在增加頻率分辨率的同時，勢必要損失一定的時間分辨率．對于設計譜擬合來說，要保證頻率分辨率的精度才能達到較好的效果．雖然損失了時間分辨率，但對于長持時的地震動而言，影響較小，而對于短持時時程，可以通過補零的方法將其持時適當地延長．

i=1,2,…,Nf

(17)

其中[PSA(fi)]target是設計譜上對應的頻率點fi的譜值，[PSA(fi)]simulated是模擬時程的偽加速度譜上對應頻率點fi的譜值．

在每次迭代后，對小波包系數矩陣進行重建，得到相應的加速度時程．其相應的偽加速度譜和設計譜的誤差可以表示為

(18)

其中Tl是反應譜的周期，E(n)(Tl)是某周期點對應的相對誤差．當平均相對誤差(Emean)不再減小時，可以判斷為迭代終止．值得注意的是，每次迭代后，不需要將加速度時程進行重新分解，而是繼續用式(15)對小波包系數矩陣進行調整即可．另外，當各個周期點中最大相對誤差(Emax)小于某一限值時，可以判斷為迭代成功．

2擬合實例

在實際的設計地震動模擬和擬合中，通常要求地震動的反應譜包絡設計譜．以核電站設計為例[1-2]，在所有控制周期(T)點中，模擬地震動反應譜低于設計譜的點數不超過5個，并且模擬反應譜與設計譜之間的相對誤差應控制在10%以內．根據上述要求，本章在周期0.03～4.00 s 范圍內選取了76個周期控制點，且設計譜的峰值設置為1g．在該例中，針對設計譜模擬了兩種不同持時的地震動加速度時程，即x1(t)和x2(t)．圖4顯示了迭代后的加速度時程x1(t)和x2(t)以及相應的速度和位移時程．可以看出，加速度時程包含了時間和頻率的非平穩性．x2(t)在整個持時前半程具有相對較高的頻率，而在后半程頻率逐漸降低．由于兩個時程均吻合設計譜，其速度與位移時程也具有相似性．

圖5顯示了模擬地震動時程的反應譜(PSA)對設計譜的擬合情況．可以看出，模擬時程的反應譜較好地包絡了設計譜，與設計譜非常接近．由于設定開始迭代的條件為初始模擬出的地震動的反應譜在各個周期上與設計譜誤差在一定范圍內，圖中兩個時程迭代前，其反應譜已經接近設計譜，這是通過設定合理的μwm和σwm來達到的．另外，通過將設計譜譜值增大一定比例的方法，就能達到規范中模擬反應譜低于設計譜點數不大于5的要求．

圖4　x1(t)和x2(t)迭代后的加速度、速度和位移時程

(a) x1(t)

(b) x2(t)

圖5擬合情況

Fig.5Fitting situation

圖6為x1(t)和x2(t)迭代過程中的收斂情況．可以看出，Emean在第1次迭代時就可以達到10%以內，并且迭代結束時，可以收斂到小于2.5%．Emax在第1次迭代時，其值相對較大，且在迭代過程中出現先下降后升高的起伏現象．這是因為在迭代過程中，出現最大相對誤差的周期點

圖6　每次迭代后的平均和最大相對誤差

在不斷地變換．最終，在20次迭代步以內，Emax可以達到7%左右的誤差．

3結語

本文基于小波包方法提出了一種新的模擬非平穩地震動且擬合設計譜的方法．該方法首先通過模擬小波包系數矩陣來生成地震動，然后根據設計譜來調整小波包系數，用迭代的方法達到模擬地震動的反應譜吻合設計譜的目的．

本文提出的方法具有以下特點：

(1)本文所模擬的地震動具有實際地震記錄的特點，即具有時間和頻率非平穩性．

(2)模擬地震動的過程簡單實用，通過指定5個參數——Td、θ、τ、μwm和σwm，可以模擬出工程需要的地震動．

(3)該方法同時包含了隨機性和確定性兩方面．隨機性表現在每次迭代失敗時將重新隨機產生高斯白噪聲；確定性表現在5個初始設定的參數，將限制模擬出的地震動在有限的范圍內變動．

(4)該方法收斂過程較快，且比較穩定，擬合設計譜相對誤差較小．

參考文獻：

[1] U.S. Atomic Energy Commission. Design response spectra for seismic design of nuclear power plants:Regulatory Guide 1.60 [R]. US:U.S. Atomic Energy Commission, 1973.

[2]國家地震局. 核電廠抗震設計規范:GB 50267—97 [S]. 北京:中國計劃出版社, 1997.

China Earthquake Administration. Code for Seismic Design of Nuclear Power Plants: GB 50267—97 [S]. Beijing:China Planning Press, 1997. (in Chinese)

[3]中華人民共和國住房和城鄉建設部. 建筑抗震設計規范:GB 50011—2010 [S]. 北京:中國建筑工業出版社, 2010.

Ministry of Housing and Urban-Rural Development of the People′s Republic of China. Code for Seismic Design of Buildings:GB 50011—2010 [S]. Beijing:China Architecture & Building Press, 2010. (in Chinese)

[4]Scanlan R H, Sachs K. Earthquake time histories and response spectra [J]. Journal of the Engineering Mechanics Division, 1974, 100(EM4):635-655.

[5]Levy S, Wilkinson J P D. Generation of artificial time-histories, rich in all frequencies, from given response spectra [J]. Nuclear Engineering and Design, 1976, 38(2):241-251.

[6]Iyengar R N, Rao P N. Generation of spectrum compatible accelerograms [J]. Earthquake Engineering & Structural Dynamics, 1979, 7(3):253-263.

[7]Preumont A. Generation of spectrum compatible accelerograms for the design of nuclear power plants [J]. Earthquake Engineering & Structural Dynamics, 1984, 12(4):481-497.

[8]Giaralis A, Spanos P D. Wavelet-based response spectrum compatible synthesis of accelerograms-Eurocode application (EC8) [J]. Soil Dynamics and Earthquake Engineering, 2009, 29(1):219-235.

[9]Cacciola P, Zentner I. Generation of response-spectrum-compatible artificial earthquake accelerograms with random joint time frequency distributions [J]. Probabilistic Engineering Mechanics, 2012, 28:52-58.

[10]Mukherjee S, Gupta V. Wavelet based generation of spectrum compatible time-histories [J]. Soil Dynamics and Earthquake Engineering, 2002, 22(9-12):799-804.

[11]Spanos P D, Giaralis A, LI Jie. Synthesis of accelerograms compatible with the Chinese GB 50011-2001 design spectrum via harmonic wavelets:Artificial and historic records [J]. Earthquake Engineering and Engineering Vibration, 2009, 8(2):189-206.

[12]Ni S H, Xie W C, Pandey M D. Tri-directional spectrum-compatible earthquake time-histories for nuclear energy facilities [J]. Nuclear Engineering and Design, 2011, 241(8):2732-2743.

[13]胡聿賢,何 訓. 考慮相位譜的人造地震動反應譜擬合[J]. 地震工程與工程振動, 1986, 6(2):37-51.

HU Yu-xian, HE Xun. Phase angle consideration in generating response spectrum-compatible ground motion [J]. Earthquake Engineering and Engineering Vibration, 1986, 6(2):37-51. (in Chinese)

[14]趙鳳新,胡聿賢. 地震動非平穩性與幅值譜和相位差譜的關系[J]. 地震工程與工程振動, 1994, 14(2):1-6.

ZHAO Feng-xin, HU Yu-xian. On the relationship of earthquake ground motion′s non-stationarity with its amplitude spectrum and phase differences spectrum [J]. Earthquake Engineering and Engineering Vibration, 1994, 14(2):1-6. (in Chinese)

[15]楊慶山,姜海鵬. 基于相位差譜的時-頻非平穩人造地震動的反應譜擬合[J]. 地震工程與工程振動, 2002, 22(1):32-38.

YANG Qing-shan, JIANG Hai-peng. Generation of response-spectrum-compatible ground motions based on phase-difference spectrum [J]. Earthquake Engineering and Engineering Vibration, 2002, 22(1):32-38. (in Chinese)

[16]韋 征,葉繼紅,包 偉,等. 考慮幅值譜與相位譜的人造地震動模擬與反應譜擬合[J]. 工程抗震與加固改造, 2007, 29(2):100-103, 106.

WEI Zheng, YE Ji-hong, BAO Wei,etal. Generation of response-spectrum-compatible ground motions consider the phase spectrum and the amplitude spectrum [J]. Earthquake Resistant Engineering and Retrofitting, 2007, 29(2):100-103, 106. (in Chinese)

[17]趙鳳新,張郁山. 多阻尼反應譜擬合的時域疊加法[J]. 核動力工程, 2008, 29(3):35-40.

ZHAO Feng-xin, ZHANG Yu-shan. Time-domain superposition method for fitting multi-damping response spectra [J]. Nuclear Power Engineering, 2008, 29(3):35-40. (in Chinese)

[18]張郁山,趙鳳新. 基于小波函數的地震動反應譜擬合方法[J]. 土木工程學報, 2014, 47(1):70-81.

ZHANG Yu-shan, ZHAO Feng-xin. Matching method of ground-motion response spectrum based on the wavelet function [J]. China Civil Engineering Journal, 2014, 47(1):70-81. (in Chinese)

[19]Saragoni G R, Hart G C. Simulation of artificial earthquakes [J]. Earthquake Engineering and Structural Dynamics, 1974, 2(3):249-267.

[20]Mallat S. A Wavelet Tour of Signal Processing [M]. 2nd ed. New York:Academic Press, 1999.

Simulation of non-stationary ground motions and fitting of design spectrum based on wavelet packet method

LIYa-nan,WANGGuo-xin*

( Faculty of Infrastructure Engineering, Dalian University of Technology, Dalian 116024, China )

Abstract：A method to simulate the non-stationary ground motions and fit the design spectrum is proposed based on wavelet packet decomposition and reconstruction. The method simulates ground motions for engineering and adjusts its frequency components so that the response spectrum will fit the design spectrum. Firstly, the wavelet packet coefficient matrix for ground motion acceleration time history is generated using 5 parameters,which characterizes the frequency attenuation with duration in real recordings. The acceleration time history can be obtained by wavelet packet reconstruction. Then, by wavelet packet decomposition the simulated ground motion is decomposed into a desired number of wavelet packet coefficient matrices with high resolution and non-overlapping frequency contents, and then the frequency component is adjusted for matching the response spectrum of the simulated accelerogram with specified design spectrum. Numerical example demonstrates that the acceleration time histories with different durations simulated by the proposed method all fit the design spectrum. The accelerograms keep the non-stationarities of ground motion even after iteration. Moreover, the method is stable and fast for convergence, and can achieve high fitting accuracy in limited iteration times.

Key words：ground motion simulation; spectrum-compatible acceleration time history; wavelet packet; non-stationary time history

中圖分類號：P315.9

文獻標識碼：A

doi:10.7511/dllgxb201603007

作者簡介:李亞楠(1987-)，男，博士生，E-mail:liyananren@163.com；王國新*(1961-)，男，博士，教授，E-mail:gxwang@dlut.edu.cn.

基金項目:國家自然科學基金資助項目(51378092，51121005).

收稿日期：2015-09-21；修回日期: 2015-11-19．

文章編號：1000-8608(2016)03-0263-07