三角级数法合成地震动

参考书籍

有关于地震动合成的主要参考书籍《地震工程学》，作者为胡聿贤。

三角级数法

其过程简述如下：

$S_a(T)$ ;
Step 2：反应谱转换为功率谱；

S(\omega_k)=\frac{2\zeta}{\pi\omega_k}S_a^2(\omega_k)/\left\{-2ln\left[-\frac{\pi}{\omega_kT_d}ln{P}\right]\right\} ,\quad \omega_k=k\cdot\Delta\omega

Step 3：功率谱转换成傅里叶谱；

A_k=A(\omega_k)=[4S(\omega_k)\Delta\omega]^{1/2}

$f(t)$ 将傅里叶谱转换成时程；

a(t)=f(t)\cdot\sum_{k=N_1}^{N_2}A_ke^{i(\omega_kt+\phi_k)}

Step 5：求该时程的反应谱；
Step 6：求该反应谱与目标反应谱的比值；根据该比值调节傅里叶谱；
重复Step 4至Step 6过程，直到求解的反应谱与目标反应谱的偏差在允许范围以内。

上面公式中:

$\omega_k$ $A_k$ $k$ 个傅里叶分量的频率和振幅；
$\phi_k$ $0$ $2\pi$ 之间均匀随机分布；
$N_1\Delta\omega<2\pi/T_M$ $N_2\Delta\omega>2\pi/T_1$ ;
$S(\omega)$ $T_d$ $\zeta$ 为阻尼比；
$P$ $P\ge0.85$ 。
$\Delta\omega$ 需足够小，即保证有足够多的三角级数。

该过程用流程图表示为：

Matlab求解

$FFT$ 进行求解，每一组随机数迭代100次，最大尝试20组随机数。


87
1
function [atOutput,t] = S2TFFT(freqSeries,spectralSeries,dampRatio,scaleFactor,toleranceUpward,toleranceDownward,totalTime,risePeriod,stablePeriod,descendFactor)
2
% freqSeries:反应谱的频率向量
3
% SpectralSeries：反应谱的谱值向量
4
% dampRatio:阻尼比
5
% scaleFactor:谱值放大系数
6
% toleranceUpward:上偏差
7
% toleranceDownard：下偏差
8
% risePeriod：上升段
9
% stablePeriod：稳定段
10
% descendFactor：衰减系数
11
fmax=max(freqSeries(1),freqSeries(end));
12
fmin=min(freqSeries(1),freqSeries(end));
13
%% 采样频率
14
fsOrder=10^(floor(log10(fmax)));
15
fs=fmax/fsOrder;
16
fInt=[1,2,2.5,4,5,8];
17
for i=1:length(fInt)
18
   if fs<=fInt(i) 
19
       fs=fInt(i);
20
       break;
21
   end   
22
end
23
fs=2*fs*fsOrder;
24
if fs<50
25
    fs=50;
26
end
27
%% 合成参数
28
df=min(fmin/20,1/totalTime);
29
nf=fs/2/df;
30
nf=2^(nextpow2(nf));
31
nfft=2*nf;
32
df=fs/nfft;
33
T=1/df;
34
%% 反应谱转换功率谱
35
P=0.9; 
36
n1=floor(fmin/df);
37
n2=ceil(fmax/df);
38
freqs=n1*df:df:n2*df;
39
Sa=interp1(freqSeries,spectralSeries,freqs,'linear','extrap');
40
Aw=zeros(nf+1,1);
41
for i=n1+1:n2+1
42
    Aw(i)=(2*dampRatio/(2*pi*pi*freqs(i-n1)))*Sa(i-n1)^2/(-2*log(-1*pi*log(P)/(2*pi*freqs(i-n1)*T)));
43
    Aw(i)=sqrt(4*Aw(i)*2*pi*df);
44
end
45
%% 包络线
46
[t,envelope]=GetEnvelope(fs,totalTime,risePeriod,stablePeriod,descendFactor);
47
%% 合成
48
for randomTimes=1:20
49
    % 随机相位
50
    phase=2*pi*rand(nf+1,1);
51
    phase(end)=0;
52
    phase(1)=0;
53
    phase=nfft*exp(phase*sqrt(-1));
54
    SseriesScaled=spectralSeries*scaleFactor;
55
    maxSseriesScaled=max(SseriesScaled);
56
    succeed=false;
57
    for cicles=1:100
58
        vecf=phase.*Aw;
59
        at=[vecf;flipud(conj(vecf(2:end-1)))];
60
        %逆变换合成时程
61
        at=ifft(at);
62
        %乘以非平稳包络线
63
        atOutput=at(1:length(t)).*envelope';
64
        % 求解反应谱
65
        Sat=tts(dampRatio,freqSeries,atOutput,1/fs);
66
        %计算上下偏差
67
        biasu=(Sat-SseriesScaled)/maxSseriesScaled;
68
        biasb=(Sat)./SseriesScaled;
69
        %如果偏差过大则调整傅里叶谱值
70
        if max(biasu)>toleranceUpward || min(biasb)<1-toleranceDownward
71
            blxs=SseriesScaled./Sat;
72
            xblxs=interp1(freqSeries,blxs,freqs,'linear','extrap');
73
            Aw(n1+1:n2+1)=Aw(n1+1:n2+1).*xblxs';
74
        else
75
            succeed=true;
76
            break;
77
        end
78
    end
79
    if succeed
80
        break;
81
    end
82
end
83
if ~succeed
84
   atOutput=[];
85
   t=[];
86
end
87
end