请注意，以下合成地震波的方法和实例只是非常简单的一个示例，仅供参考作用，合成一些有前提条件的地震波或者有其他需求的地震波此方法不一定适用！
人工合成地震波一般分成几个步骤：
- 设置地震波的主频、时间步长、持续时间等
- 根据地震模型生成相关谱
- 生成平稳过程的地震波
- 根据设置的包络函数将平稳过程转化为非平稳过程
以下为生成地震波的一个实例，其中用到的函数在下文给出：
大家可以将下面的代码复制，在Matlab中新建以下3个m文件，运行试试。
输入文件
Example.m
clearvars;close all;clc;
f = linspace(0,40,2048); % frequency vector
zeta = 0.3; % bandwidth of the earthquake excitation.
sigma = 0.3; % standard deviation of the excitation.
fn =5; % dominant frequency of the earthquake excitation (Hz).
T90 = 0.3; % value of the envelop function at 90 percent of the duration.
eps = 0.4; % normalized duration time when ground motion achieves peak.
tn = 30; % duration of ground motion (seconds).
% function call
[y,t] = seismSim(sigma,fn,zeta,f,T90,eps,tn);
% y: acceleration record
% t: time
figure
plot(t,y,'b');
xlabel('time (s)')
ylabel('ground acceleration (m/s^2)')
axis tight
set(gcf,'color','w')
guessEnvelop=[0.33,0.43,50]; % guest for envelop
guessKT = [1,1,5]; % guess for spectrum
[T90,eps,tn,zeta,sigma,fn] = fitKT(t,y,guessEnvelop,guessKT,...
'dataPlot','yes');
函数文件
seismSim.m
function [y,t] = seismSim(sigma,fn,zeta,f,T90,eps,tn)
% [y,t] = seismSim(sigma,fn,zeta,f,T90,eps,tn) generate one time series
% corresponding to acceleration record from a seismometer.
% The function requires 7 inputs, and gives 2 outputs. The time series is
% generated in two steps: First, a stationary process is created based on
% the Kanai-Tajimi spectrum, then an envelope function is used to transform
% this stationary time series into a non-stationary record. For more
% information, see [1-3].
%
% References
%
% [1] Lin, Y. K., & Yong, Y. (1987). Evolutionary Kanai-Tajimi earthquake models. Journal of engineering mechanics, 113(8), 1119-1137.
% [2] Rofooei, F. R., Mobarake, A., & Ahmadi, G. (2001).
% Generation of artificial earthquake records with a nonstationary
% Kanai-Tajimi model. Engineering Structures, 23(7), 827-837.
% [3] Guo, Y., & Kareem, A. (2016).
% System identification through nonstationary data using Time-requency Blind
% Source Separation. Journal of Sound and Vibration, 371, 110-131.
%
%
% INPUTS
%
% sigma: [1 x 1 ]: standard deviation of the excitation.
% fn: [1 x 1 ]: dominant frequency of the earthquake excitation (Hz).
% zeta: [1 x 1 ]: bandwidth of the earthquake excitation.
% f: [ 1 x M ]: frequency vector for the Kanai-tajimi spectrum.
% T90: [1 x 1 ]: value at 90 percent of the duration.
% eps: [1 x 1 ]: normalized duration time when ground motion achieves peak.
% tn: [1 x 1 ]: duration of ground motion.
%
%
% OUTPUTS
%
% y: size: [ 1 x N ] : Simulated aceleration record
% t: size: [ 1 x N ] : time
%
%
% EXAMPLE:
%
% f = linspace(0,40,2048);
% zeta = 0.3;
% sigma = 0.9;
% fn =5;
% T90 = 0.3;
% eps = 0.4;
% tn = 30;
% [y,t] = seismSim(sigma,fn,zeta,f,T90,eps,tn);
% figure
% plot(t,y);axis tight
% xlabel('time(s)');
% ylabel('ground acceleration (m/s^2)')
%
% see also fitKT.m
% Author: Etienne Cheynet - modified: 23/04/2016
%% Initialisation
w = 2*pi.*f;
fs = f(end);
dt = 1/fs;
f0= median(diff(f));
Nfreq = numel(f);
t = 0:dt:dt*(Nfreq-1);
%% Generation of the spectrum S
fn = fn *2*pi; % transformation in rad;
s0 = 2*zeta*sigma.^2./(pi.*fn.*(4*zeta.^2+1));
A = fn.^4+(2*zeta*fn*w).^2;
B = (fn.^2-w.^2).^2+(2*zeta*fn.*w).^2;
S = s0.*A./B; % single sided PSD
%% Time series generation - Monte Carlo simulation
A = sqrt(2.*S.*f0);
B =cos(w'*t + 2*pi.*repmat(rand(Nfreq,1),[1,Nfreq]));
x = A*B; % stationary process
%% Envelop function E
b = -eps.*log(T90)./(1+eps.*(log(T90)-1));
c = b./eps;
a = (exp(1)./eps).^b;
E = a.*(t./tn).^b.*exp(-c.*t./tn);
%% Envelop multiplied with stationary process to get y
y = x.*E;
endfitKT.m
function [T90,eps,tn,zeta,sigma,fn] = fitKT(t,y,guessEnvelop,guessKT,varargin)
% [T90,eps,tn,zeta,sigma,fn] = fitKT(t,y,fs,guessEnvelop,guessKT,
% varargin) fit the nonstationary Kanai-Tajimi model to ground acceleration
% records.
%
%
% INPUTS
%
% y: size: [ 1 x N ] : aceleration record
% t: size: [ 1 x N ] : time vector
% guessEnvelop: [1 x 3 ]: first guess for envelop function
% guessKT: [1 x 3 ]: first guess for Kanai-Tajimi spectrum
% varargin:
% 'F3DB' - cut off frequency for the low pass filter
% 'TolFun' - Termination tolerance on the residual sum of squares.
% Defaults to 1e-8.
% 'TolX' - Termination tolerance on the estimated coefficients
% BETA. Defaults to 1e-8.
% 'dataPlot' - 'yes': shows the results of the fitting process
%
% OUTPUTS
%
% sigma: [1 x 1 ]: Fitted standard deviation of the excitation.
% fn: [1 x 1 ]: Fitted dominant frequency of the earthquake excitation (Hz).
% zeta: [1 x 1 ]: Fitted bandwidth of the earthquake excitation.
% f: [ 1 x M ]: Fitted frequency vector for the Kanai-tajimi spectrum.
% T90: [1 x 1 ]: Fitted value at 90 percent of the duration.
% eps: [1 x 1 ]: Fitted normalized duration time when ground motion achieves peak.
% tn: [1 x 1 ]: Fitted duration of ground motion.
%
%
% EXAMPLE: this requires the function seismSim.m
%
% f = linspace(0,40,2048);
% zeta = 0.3;
% sigma = 0.3;
% fn =5;
% eta = 0.3;
% eps = 0.4;
% tn = 30;
% [y,t] = seismSim(sigma,fn,zeta,f,eta,eps,tn);
% guessEnvelop=[0.33,0.43,50];
% guessKT = [1,1,5];
% [eta,eps,tn,zeta,sigma,fn] = fitKT(t,y,guessEnvelop,guessKT,...
% 'dataPlot','yes');
%
% see also seismSim.m
% Author: Etienne Cheynet - modified: 23/04/2016
%% inputParser
p = inputParser();
p.CaseSensitive = false;
p.addOptional('f3DB',0.05);
p.addOptional('tolX',1e-8);
p.addOptional('tolFun',1e-8);
p.addOptional('dataPlot','no');
p.parse(varargin{:});
tolX = p.Results.tolX ;
tolFun = p.Results.tolFun ;
f3DB = p.Results.f3DB ;
dataPlot = p.Results.dataPlot ;
% check number of input
narginchk(4,8)
%% Get envelop parameters
dt = median(diff(t));
h1=fdesign.lowpass('N,F3dB',8,f3DB,1/dt);
d1 = design(h1,'butter');
Y = filtfilt(d1.sosMatrix,d1.ScaleValues, abs(hilbert(y)));
Y = Y./max(abs(Y));
options=optimset('Display','off','TolX',tolX,'TolFun',tolFun);
coeff1 = lsqcurvefit(@(para,t) Envelop(para,t),guessEnvelop,t,Y,[0.01,0.01,0.1],[3,3,100],options);
eps = coeff1(1);
T90 = coeff1(2);
tn = coeff1(3);
%% Get stationnary perameters for the spectrum
E =Envelop(coeff1,t);
x = y./E; % there may be better solution than this one, but I don't have better idea right now.
x(1)=0;
% calculate the PSD
[PSD,freq]=pmtm(x,7/2,numel(t),1/median(diff(t)));%%
coeff2 = lsqcurvefit(@(para,t) KT(para,freq),guessKT,freq,PSD,[0.01,0.01,1],[5,5,100],options);
zeta = coeff2(1);
sigma = coeff2(2);
fn = coeff2(3);
%% dataPLot (optional)
if strcmpi(dataPlot,'yes')
spectra = KT(coeff2,freq);
figure
subplot(211)
plot(t,y./max(abs(y)),'k',t,Envelop(coeff1,t),'b',t,Y,'r')
legend('original data','Fitted envelop','measured envelop')
title([' T_{90} = ',num2str(coeff1(2),3),'; \epsilon = ',...
num2str(coeff1(1),3),'; t_{n} = ',num2str(coeff1(3),3)]);
xlabel('time (s)')
ylabel('ground acceleration (m/s^2)')
axis tight
subplot(212)
plot(freq,PSD,'b',freq,spectra,'r')
legend('Measured','Fitted envelop')
legend('original data','Fitted Kanai-Tajimi spectrum')
title([' \zeta = ',num2str(coeff2(1),3),'; \sigma = ',...
num2str(coeff2(2),3),'; f_{n} = ',num2str(coeff2(3),3)]);
xlabel('freq (Hz)')
ylabel('Acceleration spectrum (m^2/s)')
axis tight
set(gcf,'color','w')
end
%
%% NESTED FUNCTIONS
%
function E = Envelop(para,t)
eps0 = para(1);
eta0 = para(2);
tn0 = para(3);
b = -eps0.*log(eta0)./(1+eps0.*(log(eta0)-1));
c = b./eps0;
a = (exp(1)./eps0).^b;
E = a.*(t./tn0).^b.*exp(-c.*t./tn0);
end
function S = KT(para,freq)
zeta0 = para(1);
sigma0 = para(2);
omega0 = 2*pi.*para(3);
w =2*pi*freq;
s0 = 2*zeta0*sigma0.^2./(pi.*omega0.*(4*zeta0.^2+1));
A = omega0.^4+(2*zeta0*omega0*w).^2;
B = (omega0.^2-w.^2).^2+(2*zeta0*omega0.*w).^2;
S = s0.*A./B; % single sided PSD
end
end