模态参数识别中的特征实现算法-CSDN博客

本文链接：https://blog.csdn.net/cici15874/article/details/147956214

模态参数识别中的特征实现算法（ERA）是一种基于状态空间系统识别的方法，适用于从响应信号中提取固有频率、阻尼比和模态振型。

ERA算法步骤

数据预处理

去除信号均值，消除直流分量。
截断信号长度，确保足够的数据点。

构建Hankel矩阵
使用时域响应数据构造Hankel矩阵，形式为：

H=y(t1)y(t2)⋮y(tn)y(t2)y(t3)⋮y(tn+1)⋯⋯⋱⋯y(tm)y(tm+1)⋮y(tn+m−1)

其中，n为行数，m为列数，通常满足 n+m−1=N（总采样点数）。

奇异值分解（SVD）
对Hankel矩阵进行SVD分解：

H=UΣV∗

选择有效奇异值（对应真实模态），舍弃噪声主导的小奇异值。

计算留数矩阵
利用有效奇异向量构建状态空间矩阵 A，并通过留数矩阵 R 提取极点：

R=i=1∑rσiuivi∗

其中，σi、ui、vi 分别为第i个奇异值及其对应的左右奇异向量。

提取极点（自然频率和阻尼比）
求解留数矩阵的特征方程：

det(sI−A)=0

解得极点 s=−ζωn±jωd，其中 ζ 为阻尼比，ωn 为固有频率，ωd 为阻尼固有频率。

计算模态振型
通过留数矩阵中的元素计算振型系数：

ϕ=UrRVr∗

Python代码

import numpy as np
from scipy.linalg import hankel, svd

def era_modal_parameter_identification(response, fs, max_modes=10):
    """
    ERA算法实现模态参数识别
    :param response: 响应信号 (1D numpy array)
    :param fs: 采样频率 (Hz)
    :param max_modes: 最大提取模态数
    :return: 固有频率 (rad/s), 阻尼比, 模态振型
    """
    # 1. 数据预处理
    response = response - np.mean(response)
    N = len(response)

    # 2. 构建Hankel矩阵
    L = min(max_modes * 2, N // 2)  # 行数
    H = hankel(response[:L], response[N-L:])

    # 3. SVD分解
    U, S, Vh = svd(H)

    # 4. 选择有效奇异值（根据能量阈值或固定数量）
    r = min(max_modes, np.sum(S > 1e-6))  # 自动选择有效秩
    U_r = U[:, :r]
    S_r = np.diag(S[:r])
    V_r = Vh[:r, :]

    # 5. 计算留数矩阵
    R = U_r @ S_r @ V_r

    # 6. 提取极点（自然频率和阻尼比）
    poles = np.linalg.eigvals(-R)  # 极点s = -sigma（假设R为状态矩阵）
    freq = np.angle(poles) * fs / (2 * np.pi)  # 固有频率 (Hz)
    zeta = -np.real(poles) / np.abs(poles)     # 阻尼比

    # 7. 计算模态振型（示例简化，实际需结合留数矩阵）
    modes = U_r @ R  # 振型近似

    return freq, zeta, modes

# 示例使用
if __name__ == "__main__":
    # 生成模拟响应信号（示例）
    t = np.linspace(0, 10, 1000)
    fs = 100  # Hz
    # 单自由度系统响应（示例）
    omega_n = 2 * np.pi * 5  # 5 Hz
    zeta = 0.05
    response = np.exp(-zeta * omega_n * t) * np.sin(omega_n * np.sqrt(1 - zeta**2) * t)

    # 调用ERA
    freq, zeta_est, modes = era_modal_parameter_identification(response, fs)

    print("Estimated Natural Frequencies (Hz):", freq)
    print("Estimated Damping Ratios:", zeta_est)