地震动反应谱计算及分析

地震动反应谱的频域计算原理与实现

在地震工程学中，地震动反应谱是描述结构动力响应最重要的工具之一。本文将介绍一种基于频域方法的反应谱计算原理，并结合 MATLAB 实现过程，展示如何高效地计算多个阻尼比下的位移、速度与加速度谱。

一、反应谱基本定义

对一个线性单自由度振动系统，其响应受地震动激励影响，反应谱定义为该系统在不同自振周期（$T$）和阻尼比（$\xi$）下的最大响应值，具体包括：

• 加速度反应谱（$S_a(T)$）
• 速度反应谱（$S_v(T)$）
• 位移反应谱（$S_d(T)$）

二、频域法的理论基础

1. 系统动力学方程

考虑质量 $m$、阻尼比 $\xi$、固有频率 $\omega_n = \frac{2\pi}{T}$ 的单自由度系统，受到地震动加速度 $\ddot{u}_g(t)$ 激励，其相对位移 $u(t)$ 满足：

$$\ddot{u}(t) + 2\xi\omega_n\dot{u}(t) + \omega_n^2 u(t) = -\ddot{u}_g(t)$$

2. 频域响应表达式

将上式进行傅里叶变换，得到系统响应的频域表示：

$$U(\omega) = H(\omega) \cdot A_g(\omega)$$

其中：

• $A_g(\omega)$：地震动加速度的频谱（由 FFT 得到）
• $H(\omega)$：系统的频率响应函数，表示对单位输入的响应：

$$H(\omega) = \frac{1}{-\omega^2 + 2i\xi\omega\omega_n + \omega_n^2}$$

三、三种响应谱的频域计算方法

1. 位移谱 $S_d(T)$

$$U(\omega) = H(\omega) \cdot A_g(\omega)$$

取反傅里叶变换 $u(t) = \mathcal{F}^{-1}[U(\omega)]$，位移谱为其最大绝对值：

$$S_d(T) = \max |u(t)|$$

2. 速度谱 $S_v(T)$

$$V(\omega) = i\omega \cdot H(\omega) \cdot A_g(\omega)$$

反变换得 $v(t)$，求最大值：

$$S_v(T) = \max |v(t)|$$

3. 加速度谱 $S_a(T)$

计算相对加速度：

$$A_{\text{rel}}(\omega) = -\omega^2 \cdot H(\omega) \cdot A_g(\omega)$$

加上原始地震动得绝对加速度：

$$a_{\text{abs}}(t) = a_{\text{rel}}(t) + \ddot{u}_g(t)$$

谱值为其最大值除以重力加速度：

$$S_a(T) = \frac{\max |a_{\text{abs}}(t)|}{g}$$

四、阻尼比与周期的多维计算

在实际工程中，通常会选取多个阻尼比（例如 0%、1%、2%、5%、10%）以及一定周期范围（如 $T \in [0.01, 10]\, \text{s}$），分别计算对应的谱曲线，形成完整的谱图族。

五、频域方法的优势

与时域积分法相比，频域方法具有以下优点：

• 效率高：通过快速傅里叶变换（FFT）加速运算，适用于大量地震波处理；
• 直观分析：频域可揭示系统与激励的滤波关系；
• 适应性强：便于扩展至不同频率带宽和采样率的地震动。

MATLAB实现

clear; clc;

% 参数设置
damping_ratios = [ 0.01, 0.02, 0.05, 0.10];  % 阻尼比
T = linspace(0.01, 10, 1000);                    % 周期范围
dt = 0.005;                                      % 地震动采样间隔（秒）
g = 9.81;                                        % 重力加速度

% 场地目录名
site_dirs = {'1类场地', '2类场地', '3类场地', '4类场地'}; %将地震波按场地类型分类，存放至各目录下，地震波文件类型为.txt。

for i = 1:length(site_dirs)
    site_path = site_dirs{i};
    files = dir(fullfile(site_path, '*.txt'));  % 该场地所有地震动

    for j = 1:length(files)
        file_path = fullfile(site_path, files(j).name);
        acc = load(file_path);  % 读取地震动加速度，单位 g
        acc = acc(:).'.*g;
        n = length(acc);

        % 初始化谱矩阵（行：不同阻尼比，列：不同T）
        Sa_matrix = zeros(length(damping_ratios), length(T));
        Sv_matrix = zeros(length(damping_ratios), length(T));
        Sd_matrix = zeros(length(damping_ratios), length(T));

        for k = 1:length(damping_ratios)
            xi = damping_ratios(k);
            Sa = zeros(1, length(T));
            Sv = zeros(1, length(T));
            Sd = zeros(1, length(T));

            for t = 1:length(T)
                wn = 2 * pi / T(t);
                Nfft = 2^nextpow2(n);
                omega = 2 * pi * (0:Nfft-1) / (Nfft * dt);
                af = fft(acc, Nfft);

                % 傅里叶频域响应函数 H(ω)：从地震动加速度 -> 相对位移
                H = 1 ./ (-(omega.^2) + 2i * xi * omega * wn + wn^2);

                % 相对位移（位移谱）
                u = real(ifft(af .* H));
                u = u(1:n);
                Sd(t) = max(abs(u));

                % 相对速度
                v = real(ifft(af .* H .* 1i .* omega));
                v = v(1:n);
                Sv(t) = max(abs(v));  

                % 相对加速度
                a_rel = real(ifft(-af .* omega.^2 .* H));
                a_rel = a_rel(1:n);

                % 绝对加速度谱
                a_abs = a_rel + acc;
                Sa(t) = max(abs(a_abs)) / g;
            end


            % 存储当前阻尼比下所有谱
            Sd_matrix(k, :) = Sd;
            Sv_matrix(k, :) = Sv;
            Sa_matrix(k, :) = Sa;
        end

        % 绘图：加速度谱
        figure('Name', ['Sa - ', files(j).name], 'NumberTitle', 'off');
        hold on; grid on;
        colors = lines(length(damping_ratios));
        for k = 1:length(damping_ratios)
            plot(T, Sa_matrix(k, :), 'DisplayName', ['ξ = ', num2str(damping_ratios(k)*100), '%'], ...
                'Color', colors(k, :), 'LineWidth', 1.5);
        end
        xlabel('周期 T (s)'); ylabel('S_a (g)');
        title(['加速度反应谱 - ', files(j).name]);
        legend show; xlim([0 10]); ylim auto;
        saveas(gcf, fullfile(site_path, [files(j).name(1:end-4), '_Sa.png']));
        close;

        % 绘图：速度谱
        figure('Name', ['Sv - ', files(j).name], 'NumberTitle', 'off');
        hold on; grid on;
        for k = 1:length(damping_ratios)
            plot(T, Sv_matrix(k, :), 'DisplayName', ['ξ = ', num2str(damping_ratios(k)*100), '%'], ...
                'Color', colors(k, :), 'LineWidth', 1.5);
        end
        xlabel('周期 T (s)'); ylabel('S_v (m/s)');
        title(['速度反应谱 - ', files(j).name]);
        legend show; xlim([0 10]); ylim auto;
        saveas(gcf, fullfile(site_path, [files(j).name(1:end-4), '_Sv.png']));
        close;

        % 绘图：位移谱
        figure('Name', ['Sd - ', files(j).name], 'NumberTitle', 'off');
        hold on; grid on;
        for k = 1:length(damping_ratios)
            plot(T, Sd_matrix(k, :), 'DisplayName', ['ξ = ', num2str(damping_ratios(k)*100), '%'], ...
                'Color', colors(k, :), 'LineWidth', 1.5);
        end
        xlabel('周期 T (s)'); ylabel('S_d (m)');
        title(['位移反应谱 - ', files(j).name]);
        legend show; xlim([0 10]); ylim auto;
        saveas(gcf, fullfile(site_path, [files(j).name(1:end-4), '_Sd.png']));
        close;
    end
end