非平穩地震動場的隨機函數模擬

阮鑫鑫 劉子心 劉章軍，2

(1. 三峽大學 土木與建筑學院， 湖北 宜昌 443002； 2. 防災減災湖北省重點實驗室(三峽大學)， 湖北 宜昌 443002)

地震災害作用具有顯著的隨機動力特性和空間變異特性，一般用時間變量和空間變量的隨機場來描述，即所謂的時空隨機場[1]．對于復雜工程結構(如大跨橋梁、大壩等)，地震災害作用的空間變異特性不可忽視[2]．因此，在復雜工程結構的非線性隨機地震反應分析中，時空隨機場的模擬至關重要．目前，蒙特卡羅隨機模擬方法在工程實際中受到越來越多學者的關注和重視，其中譜表示法由于其精度高和易于實現而被廣泛使用．40多年來，Yang[3]、Shinozuka[4]、Deodatis[5]等在譜表示的理論及應用方面開展了卓有成效的工作，使得譜表示方法臻于完善．

然而，采用傳統的譜表示方法模擬地震動時空隨機場時，往往需要高達數萬個隨機變量來表達地震動時空隨機場，這為復雜工程結構的非線性隨機地震反應分析與抗震可靠度計算帶來巨大挑戰[6]．因此，如何有效地降低譜表示方法隨機變量的數量已成為近年來的研究熱點．為此，陳建兵等[7]發展了隨機過程的隨機諧和函數表達，通過少量的隨機諧和分量即可獲取精確的目標功率譜．隨后，劉章軍等[8-10]提出了一種隨機函數降維方法，實現了僅用1-2個基本隨機變量即可精細地模擬一維單變量平穩和非平穩隨機過程以及多變量平穩隨機過程．

為此，本文將在地震動時空隨機場的原始譜表示基礎上，引入隨機函數的約束形式，實現2D-1V地震動時空隨機場的降維模擬，并結合FFT算法極大地提高模擬的計算效率，從而為地震動時空隨機場模擬提供一種更為精細有效的方法．

1 2D-1V時空隨機場的譜表示

設f0(x，t)是一個零均值、實值的2D-1V時空隨機場，在時間域上是平穩或非平穩，在空間域上是均勻或非均勻，單邊的頻率-波數譜Gf0(x，t；κ，ω)=2Sf0(x，t；κ，ω)，其中Sf0(x，t；κ，ω)為雙邊的頻率-波數譜．于是，2D-1V時空隨機場f0(x，t)則可表達為如式(1)的積分形式[4]：

sin(κx+ωt)dVx，t(κ，ω) (1)

式中，Ux，t(κ，ω)和Vx，t(κ，ω)是定義在-∞<κ<∞和0<ω<∞上的兩個實值正交過程，其增量dUx，t(κ，ω)和dVx，t(κ，ω)滿足如式(2)的條件：

(2)

顯然，當dUx，t(κ，ω)和dVx，t(κ，ω)均與空間無關時，式(1)即為非平穩均勻的2D-1V時空隨機場的譜表示．

進一步，將式(1)近似寫為如下的有限級數形式：

sin(κmx+ωnt)ΔVx，t(κm，ωn)]+

[cos(-κmx+ωnt)ΔUx，t(-κm，ωn)+

sin(-κmx+ωnt)ΔVx，t(-κm，ωn)]} (3)

其中

κm=(m-0.5)Δκ，ωn=(n-0.5)Δω(4)

式中，M和N分別為波數和頻率的截斷項數，Δκ為波數步長，Δω為頻率步長．Δκ與Δω應足夠小，以便保證式(3)能近似代替式(1)．

在式(3)中，若定義正交增量為

將式(5)代入式(3)中，得到

式(7)即為2D-1V時空隨機場的基于正交隨機變量的譜表示，也就是所謂的原始譜表示．

(8)

于是，將式(8)代入式(7)中，即可得到2D-1V時空隨機場模擬的基于隨機相位角的傳統譜表示[5]．

2 2D-1V時空隨機場的隨機函數模擬

需要指出的是，式(10)可以采用快速傅里葉變換(FFT)來提高計算效率；同時，基本隨機變量Θ1和Θ2是采用數論方法[11]來生成代表性點集，而非蒙特卡羅方法生成．

在利用式(10)模擬2D-1V時空隨機場時，首先需要給出2D-1V時空隨機場的頻率-波數譜Sf0(x，t；κ，ω)的表達式，對于非平穩均勻的2D-1V時空隨機場的頻率-波數譜應為Sf0(t；κ，ω)．

3 地震動場模型

為了模擬非平穩均勻地震動隨機場，首先需要給出相應的平穩均勻地震動隨機場的頻率-波數譜．對于平穩均勻地震動隨機場的頻率-波數譜S(κ，ω)，可由平穩地震動過程的自功率譜S(ω)與地震動場的空間相干函數的傅里葉變換之乘積得到[2]：

(12)

在本文中，地震動場的空間相干函數采用考慮行波效應的Luco-Wong模型[2]

γ(ξ，ω)=e-α2ω2ξ2e-iωξ/νapp(13)

將式(13)代入式(12)，得到頻率-波數譜S(κ，ω)的表達式為[2]

(14)

式中，參數α=0.147，視波速νapp=1 000 m/s．

對于平穩地震動加速度過程的自功率譜S(ω)，本文采用Kanai-Tajimi譜[2]：

(15)

式中，場地卓越圓頻率ωg=5π rad/s，阻尼比ζg=0.6，譜強度因子S0=44.29 cm2/s3．

為了反映地震動過程的時-頻非平穩性，本文采用如下的時-頻調制函數[9]：

其中

式中，參數b=a+0.001，c=0.005．這樣，上述時-頻調制函數中只有一個獨立參數a，本文取a=0.2s-1．

于是，非平穩均勻地震動時空隨機場的頻率-波數譜的表達式為：

Sf0(t；κ，ω)=A2(t，ω)S(κ，ω) (17)

4 數值算例

現在，在水平地面上沿x方向選擇2 048個作用點，且兩兩間距為10 m，地震動有效持時20 s．為了模擬非平穩均勻地震動時空隨機場，采用如下計算參數：截斷波數κu=π/10 rad·m-1，截斷頻率ωu=50π rad·s-1，波數離散點數M=512，頻率離散點數N=1 024．為利用FFT技術，本文使用Peng等[12]建議的方法，將非平穩隨機地震動場分解為6個強度調制的隨機地震動場之和．選擇FFT的波數和頻率變換點數為2M和2N，因此對應空間間隔Δx=10.0 m，時間間隔Δt=0.02 s，場的空間尺度L=20 480 m，地震動模擬持時T=81.92 s．雖然地震動的模擬持時大于有效持時，但這并不會增加過多的計算工作量，因為本文方法的模擬效率足夠高．在確定模擬參數后，利用數論方法選取基本隨機變量Θ1和Θ2的代表性點集，本文選取代表性點的數量nsel=377，再利用式(10)即可生成非平穩均勻地震動時空隨機場的代表性時程．

為分析方便，選擇3個位置處的地震動隨機過程的代表性時程進行分析．其中，位置1為第201點處，即x1=2 000 m；位置2為第205點處，即x2=2 040 m；位置3為第211點處，即x3=2 100 m．圖1為本文方法生成的3個位置處的地震動加速度代表性時程，可見代表性時程在時間和頻率上具有顯著的非平穩特性．此外，3個位置處的代表性時程在波形變化方面也具有一定的相似性，這反映了3個點之間的空間相干效應．

(a)位置1處;(b)位置2處;(c)位置3處圖1 非平穩地震動加速度代表性時程

圖2為本文方法生成的位置1處377條地震動加速度代表性時程的均值及標準差與目標值的比較．從圖中可知，均值和標準差的模擬值與目標值均擬合很好．對于位置2和3處的地震動加速度代表性時程的均值和標準差，也能得到類似的結果，這初步驗證了本文方法的有效性．

圖2 位置1處的377條非平穩地震動加速度代表性時程的均值和標準差

圖3為本文方法生成的377條地震動加速度代表性時程在位置1、2、3處的自相關函數及互相關函數與目標值的比較，可見3個位置處的自相關與互相關函數的模擬值均與目標值擬合一致，這進一步驗證了本文方法的有效性．

圖3 在t=5 s時3點的相關函數模擬值與目標值的比較

為了進一步說明本文方法的優越性，在模擬效率和精度方面與傳統的譜表示方法進行比較．在模擬效率方面，模擬2 048個位置點處隨機過程的全部樣本時，本文方法的耗時為673 s，傳統譜表示方法的耗時為664 s．在模擬精度方面，本文方法的均值相對誤差和標準差相對誤差分別為2.55%和2.83%，而傳統譜表示方法的均值相對誤差和標準差相對誤差分別為4.10%和2.96%．其中，均值和標準差的相對誤差定義見文獻[8]，且最終的相對誤差均取3個位置處相對誤差的平均值．由此可見，本文方法的模擬效率與傳統方法大致相當，但本文方法的模擬精度更高，體現了優越性．

5 結 論

本文建議了一種基于隨機函數的時空隨機場降維模擬方法．該方法是在基于正交隨機變量的譜表示上，通過將譜表示中的正交隨機變量均定義為兩個基本隨機變量的正交函數形式，從而實現地震動時空隨機場的降維模擬．數值算例表明，本文方法的模擬效率與傳統方法大致相當，但本文方法的精度更高．同時，由于本文方法僅需2個基本隨機變量表達地震動隨機場，從而可利用數論方法生成具有完備概率的代表性樣本集，這為應用概率密度演化理論進行結構的隨機地震反應和抗震可靠度分析奠定了基礎．