主题063：结构动力学模态参数识别

摘要

本主题深入探讨结构动力学中的模态参数识别理论与方法。模态参数识别是结构健康监测、振动控制和有限元模型验证的核心技术。本主题系统介绍频域识别方法（峰值拾取法、频域分解法）、时域识别方法（ITD法、STD法、ERA法）、随机子空间识别方法（SSI）以及基于机器学习的智能识别方法。通过四个完整的工程案例，详细演示如何从实测振动数据中提取结构的固有频率、阻尼比和模态振型，为工程结构的动态特性分析和健康评估提供理论基础和实用工具。

关键词

模态参数识别；频域分解；随机子空间；工作模态分析；模态振型；阻尼比估计；机器学习；结构健康监测

---

1. 理论基础

1.1 模态参数识别的基本概念

模态参数识别是通过分析结构在激励下的动力响应，提取结构的模态参数（固有频率、阻尼比、模态振型）的过程。这些参数完整地描述了结构的动态特性。

1.1.1 模态参数的定义

对于N自由度线性系统，运动方程为：

$$M\ddot{x}+C\dot{x}+Kx=f(t)$$

其中，$M$、$C$、$K$分别为质量、阻尼和刚度矩阵。通过模态变换，可得到模态坐标下的解耦方程：

$${\ddot{q}}_r+2{\zeta }_r{\omega }_r{\dot{q}}_r+{\omega }_r^2{q}_r=\frac{{\varphi }_r^{T}f(t)}{m_r}$$

其中：

• ${\omega }_r$：第r阶固有频率（rad/s）
• ${\zeta }_r$：第r阶模态阻尼比
• ${\varphi }_r$：第r阶模态振型
• $m_r$：第r阶模态质量

1.1.2 频响函数与模态参数的关系

频响函数（FRF）是模态参数识别的基础。对于比例阻尼系统，频响函数可表示为：

$$H_{ij}(\omega )=\sum _{r=1}^N\frac{{\varphi }_{ir}{\varphi }_{jr}}{m_r({\omega }_r^2-{\omega }^2+2i{\zeta }_r{\omega }_r\omega )}$$

在共振频率附近，第r阶模态占主导：

$$H_{ij}(\omega )\approx \frac{{\varphi }_{ir}{\varphi }_{jr}}{m_r({\omega }_r^2-{\omega }^2+2i{\zeta }_r{\omega }_r\omega )}$$

1.1.3 模态参数识别的分类

根据识别域的不同，模态参数识别方法可分为：

1. 频域方法：基于频响函数或功率谱密度

   • 峰值拾取法（Peak Picking）
   • 频域分解法（FDD）
   • 多参考最小二乘复频域法（p-LSCF）

2. 时域方法：基于脉冲响应函数或自由衰减响应

   • 特征系统实现算法（ERA）
   • 多参考 Ibrahim 时域法（ITD）
   • 特征值实现算法（ERA）

3. 时频域方法：结合时域和频域信息

   • 小波变换法
   • Hilbert-Huang变换

4. 随机子空间方法：基于随机过程理论

   • 数据驱动随机子空间（SSI-data）
   • 协方差驱动随机子空间（SSI-cov）

5. 机器学习方法：基于数据驱动

   • 神经网络方法
   • 支持向量机
   • 深度学习

1.2 频域识别方法

1.2.1 峰值拾取法

峰值拾取法是最简单的频域识别方法，基于频响函数幅值的峰值确定固有频率。

基本原理：
在共振频率处，频响函数的幅值达到极大值：

$$|H_{ij}({\omega }_r)|=\frac{|{\varphi }_{ir}{\varphi }_{jr}|}{2m_r{\zeta }_r{\omega }_r^2}$$

半功率带宽法估计阻尼：

$${\zeta }_r=\frac{{\omega }_b-{\omega }_a}{2{\omega }_r}$$

其中${\omega }_a$和${\omega }_b$是半功率点（幅值为峰值$1/\sqrt{2}$处的频率）。

1.2.2 频域分解法（FDD）

频域分解法利用输出谱密度的奇异值分解提取模态参数。

算法步骤：

1. 计算输出谱密度矩阵：
   $$G_{yy}(\omega )=H(\omega )G_{ff}(\omega )H^H(\omega )$$

2. 对每个频率进行SVD分解：
   $$G_{yy}(\omega )=U(\omega )\Sigma (\omega )V^H(\omega )$$

3. 在共振频率处，第一奇异值对应主导模态：
   $${\sigma }_1({\omega }_r)\gg {\sigma }_2({\omega }_r)$$

4. 模态振型从奇异向量提取：
   $${\varphi }_r\approx u_1({\omega }_r)$$

1.2.3 多参考最小二乘复频域法（p-LSCF）

p-LSCF方法通过有理多项式拟合频响函数：

$$H_{ij}(s)=\frac{{\sum }_{k=0}^n{b}_k{s}^k}{{\sum }_{k=0}^n{a}_k{s}^k},s=i\omega$$

通过最小二乘优化确定系数$a_k$和$b_k$，然后求解分母多项式的根得到极点：

$$s_r=-{\zeta }_r{\omega }_r\pm i{\omega }_r\sqrt{1-{\zeta }_r^2}$$

1.3 时域识别方法

1.3.1 特征系统实现算法（ERA）

ERA方法基于系统的脉冲响应函数，通过Hankel矩阵的奇异值分解提取模态参数。

算法步骤：

1. 构建Hankel矩阵：
   $$H(k)=\left[\begin{array}{cccc}h(k)& h(k+1)& \cdots & h(k+\beta -1)\\ h(k+1)& h(k+2)& \cdots & h(k+\beta )\\ ⋮& ⋮& \ddots & ⋮\\ h(k+\alpha -1)& h(k+\alpha )& \cdots & h(k+\alpha +\beta -2)\end{array}\right]$$

2. 对$H(0)$进行SVD分解：
   $$H(0)=U\Sigma V^T$$

3. 确定系统阶数，保留前n个奇异值：
   $$H(0)\approx U_n{\Sigma }_n{V}_n^T$$

4. 构建系统矩阵：
   $$A={\Sigma }_n^{-1/2}{U}_n^{T}H(1)V_n{\Sigma }_n^{-1/2}$$

5. 对A进行特征值分解：
   $$A{\psi }_r={\lambda }_r{\psi }_r$$

6. 提取模态参数：
   $${\omega }_r=|\ln ({\lambda }_r)|/\Delta t$$
   $${\zeta }_r=-\text{Re}(\ln ({\lambda }_r))/|\ln ({\lambda }_r)|$$

1.3.2 Ibrahim时域法（ITD）

ITD方法利用自由响应数据构建特征值问题。

基本原理：

对于自由振动响应：

$$x(t)=\sum _{r=1}^N{\varphi }_r{e}^{-{\zeta }_r{\omega }_{r}t}\sin ({\omega }_{dr}t+{\theta }_r)$$

通过构造响应矩阵和延迟响应矩阵，形成广义特征值问题：

$$A\Phi =B\Phi \Lambda$$

其中$\Lambda$包含系统的特征值，从中可提取模态参数。

1.4 随机子空间识别方法

1.4.1 随机状态空间模型

对于随机激励下的线性系统，状态空间方程为：

$$x_{k+1}=A{x}_k+w_k$$
$$y_k=C{x}_k+v_k$$

其中：

• $x_k$：状态向量
• $y_k$：观测向量
• $w_k$：过程噪声
• $v_k$：测量噪声

1.4.2 协方差驱动随机子空间（SSI-cov）

算法步骤：

1. 计算输出协方差矩阵：
   $$R_i=E[y_{k+i}{y}_k^T]$$

2. 构建块Toeplitz矩阵：
   $$T_{1|i}=\left[\begin{array}{cccc}{R}_i& R_{i-1}& \cdots & R_1\\ R_{i+1}& R_i& \cdots & R_2\\ ⋮& ⋮& \ddots & ⋮\\ R_{2i-1}& R_{2i-2}& \cdots & R_i\end{array}\right]$$

3. 对Toeplitz矩阵进行SVD分解：
   $$T_{1|i}=U\Sigma V^T$$

4. 确定系统阶数，提取观测矩阵和状态矩阵：
   $$O_i=U_n{\Sigma }_n^{1/2}$$
   $${\Gamma }_i={\Sigma }_n^{1/2}{V}_n^T$$

5. 从观测矩阵提取C，从状态矩阵提取A

6. 对A进行特征值分解，提取模态参数

1.4.3 数据驱动随机子空间（SSI-data）

SSI-data方法直接使用测量数据，无需计算协方差。

算法步骤：

1. 构建数据矩阵：
   $$Y_{0|i-1}=\left[\begin{array}{cccc}{y}_0& y_1& \cdots & y_{j-1}\\ y_1& y_2& \cdots & y_j\\ ⋮& ⋮& \ddots & ⋮\\ y_{i-1}& y_i& \cdots & y_{i+j-2}\end{array}\right]$$

2. 对数据矩阵进行QR分解或SVD分解

3. 提取扩展观测矩阵

4. 计算系统矩阵A和C

5. 特征值分解提取模态参数

1.5 稳定图与模态定阶

1.5.1 稳定图的概念

稳定图是确定系统真实模态阶数的重要工具。通过逐渐增加模型阶数，观察模态参数的稳定性。

稳定准则：

• 频率稳定：$|\frac{f^{(n+1)}-f^{(n)}}{f^{(n)}}|<{ϵ}_f$
• 阻尼稳定：$|\frac{{\zeta }^{(n+1)}-{\zeta }^{(n)}}{{\zeta }^{(n)}}|<{ϵ}_{\zeta }$
• MAC稳定：$MAC({\varphi }^{(n+1)},{\varphi }^{(n)})>{ϵ}_{MAC}$

1.5.2 模态定阶方法

1. 奇异值跳跃法：观察奇异值曲线的拐点
2. 稳定图法：识别稳定的模态
3. 信息准则：AIC、MDL等统计准则
4. 交叉验证：预测误差最小化

1.6 模态验证指标

1.6.1 模态置信准则（MAC）

MAC用于评估两个模态振型的相关性：

$$MA{C}_{ij}=\frac{|{\varphi }_i^H{\varphi }_j{|}^2}{({\varphi }_i^H{\varphi }_i)({\varphi }_j^H{\varphi }_j)}$$

• $MAC=1$：完全相关
• $MAC=0$：完全无关
• $MAC>0.9$：认为同一模态

1.6.2 模态参与因子

模态参与因子表示各模态对响应的贡献：

$$MP{F}_r=\frac{{\varphi }_r^{T}M1}{m_r}$$

1.6.3 模态复杂性

模态复杂性指标评估模态振型的实模态特性：

$$MCI=1-\frac{|\sum {\varphi }_{ir}^2{|}^2}{(\sum |{\varphi }_{ir}{|}^2{)}^2}$$

• $MCI=0$：实模态
• $MCI>0.2$：复模态

1.7 机器学习方法在模态识别中的应用

1.7.1 基于神经网络的模态识别

基本原理：
利用神经网络的非线性映射能力，从振动数据直接学习模态参数。

网络结构：

• 输入层：振动信号特征（频谱、小波系数等）
• 隐藏层：多个全连接层或卷积层
• 输出层：模态参数（频率、阻尼、振型）

损失函数：
$$L=\alpha L_{freq}+\beta L_{damping}+\gamma L_{mode}$$

1.7.2 卷积神经网络（CNN）

CNN可用于从时频图像中提取模态特征。

网络架构：

输入（时频图）→ 卷积层 → 池化层 → 卷积层 → 全连接层 → 输出（模态参数）

1.7.3 自编码器与特征提取

自编码器可用于降维和特征提取：

$$\min ||x-\hat{x}|{|}^2$$

其中$\hat{x}=D(E(x))$，$E$为编码器，$D$为解码器。

---

2. 数值方法

2.1 频响函数估计

2.1.1 H1估计

$${\hat{H}}_1(f)=\frac{G_{xy}(f)}{G_{xx}(f)}$$

适用于输出噪声占主导的情况。

2.1.2 H2估计

$${\hat{H}}_2(f)=\frac{G_{yy}(f)}{G_{yx}(f)}$$

适用于输入噪声占主导的情况。

2.1.3 Hv估计（总最小二乘）

$${\hat{H}}_v(f)=\sqrt{{\hat{H}}_1(f){\hat{H}}_2(f)}$$

2.2 谱估计方法

2.2.1 周期图法

$${\hat{G}}_{xx}(f)=\frac{1}{N}|X(f){|}^2$$

2.2.2 Welch方法

将数据分段，计算平均周期图：

$${\hat{G}}_{xx}(f)=\frac{1}{M}\sum _{i=1}^M{\hat{G}}_{xx}^{(i)}(f)$$

2.2.3 多窗法（MTM）

使用多个正交窗函数减少估计方差。

2.3 信号预处理

2.3.1 去趋势

去除信号的线性趋势：

$$y(t)=x(t)-(a+bt)$$

2.3.2 滤波

带通滤波保留感兴趣的频段：

$$H(f)=\{\begin{array}{ll}1& f_1\le f\le f_2\\ 0& \text{其他}\end{array}$$

2.3.3 重采样

根据奈奎斯特准则确定采样频率：

$$f_s>2f_{max}$$

---

3. 工程案例

3.1 案例1：基于频域的模态参数识别

问题描述：
对一个简支梁进行锤击试验，利用频响函数数据，采用峰值拾取法和频域分解法识别结构的模态参数。

Python代码实现：

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
import matplotlib
matplotlib.use('Agg')
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
from scipy.linalg import eigh, svd
from scipy.signal import find_peaks, welch
from scipy.fft import fft, ifft, fftfreq

class FrequencyDomainIdentification:
    """频域模态参数识别"""
    
    def __init__(self, fs=1000):
        """
        初始化
        
        Parameters:
        -----------
        fs : float
            采样频率 (Hz)
        """
        self.fs = fs
        self.dt = 1.0 / fs
    
    def compute_frf(self, force, response, window='hann', nperseg=None):
        """
        计算频响函数 (H1估计)
        
        Parameters:
        -----------
        force : array
            力信号
        response : array
            响应信号
        window : str
            窗函数类型
        nperseg : int
            每段长度
        
        Returns:
        --------
        f : array
            频率向量
        H : array
            频响函数
        coherence : array
            相干函数
        """
        if nperseg is None:
            nperseg = min(256, len(force) // 4)
        
        # 计算互功率谱和自功率谱
        f, Gxf = welch(force, fs=self.fs, window=window, 
                       nperseg=nperseg, noverlap=nperseg//2, 
                       detrend='constant', return_onesided=True)
        
        _, Gyx = welch(response, force, fs=self.fs, window=window,
                       nperseg=nperseg, noverlap=nperseg//2,
                       detrend='constant', return_onesided=True)
        
        _, Gxx = welch(force, fs=self.fs, window=window,
                       nperseg=nperseg, noverlap=nperseg//2,
                       detrend='constant', return_onesided=True)
        
        _, Gyy = welch(response, fs=self.fs, window=window,
                       nperseg=nperseg, noverlap=nperseg//2,
                       detrend='constant', return_onesided=True)
        
        # H1估计
        H = Gyx / (Gxx + 1e-10)
        
        # 相干函数
        coherence = np.abs(Gyx)**2 / (Gxx * Gyy + 1e-10)
        
        return f, H, coherence
    
    def peak_picking(self, f, H, peak_height=None, min_distance=5):
        """
        峰值拾取法识别模态参数
        
        Parameters:
        -----------
        f : array
            频率向量
        H : array
            频响函数
        peak_height : float
            峰值高度阈值
        min_distance : int
            峰值间最小距离
        
        Returns:
        --------
        modes : list
            识别的模态参数 [(freq, damping, mode_shape), ...]
        """
        # 计算幅值
        magnitude = np.abs(H)
        
        if peak_height is None:
            peak_height = np.max(magnitude) * 0.1
        
        # 寻找峰值
        peaks, properties = find_peaks(magnitude, height=peak_height, 
                                       distance=min_distance)
        
        modes = []
        
        for peak_idx in peaks:
            freq = f[peak_idx]
            
            # 半功率带宽法估计阻尼
            peak_mag = magnitude[peak_idx]
            half_power = peak_mag / np.sqrt(2)
            
            # 寻找半功率点
            left_idx = peak_idx
            right_idx = peak_idx
            
            while left_idx > 0 and magnitude[left_idx] > half_power:
                left_idx -= 1
            
            while right_idx < len(magnitude) - 1 and magnitude[right_idx] > half_power:
                right_idx += 1
            
            # 计算阻尼比
            if left_idx < peak_idx and right_idx > peak_idx:
                f_left = f[left_idx]
                f_right = f[right_idx]
                damping = (f_right - f_left) / (2 * freq)
            else:
                damping = 0.01  # 默认值
            
            # 模态振型（简化，实际应从多测点数据提取）
            mode_shape = np.array([1.0])  # 单点测量
            
            modes.append({
                'frequency': freq,
                'damping': damping,
                'mode_shape': mode_shape,
                'peak_idx': peak_idx
            })
        
        return modes
    
    def frequency_domain_decomposition(self, f, H_matrix):
        """
        频域分解法 (FDD)
        
        Parameters:
        -----------
        f : array
            频率向量
        H_matrix : array
            多参考频响函数矩阵 (n_freq x n_dof x n_ref)
        
        Returns:
        --------
        modes : list
            识别的模态参数
        singular_values : array
            奇异值
        """
        n_freq = len(f)
        n_dof = H_matrix.shape[1]
        
        # 计算谱密度矩阵
        Gyy = np.zeros((n_freq, n_dof, n_dof), dtype=complex)
        
        for i in range(n_freq):
            H_i = H_matrix[i, :, :]
            Gyy[i, :, :] = H_i @ H_i.conj().T
        
        # SVD分解
        singular_values = np.zeros((n_freq, n_dof))
        mode_shapes = []
        
        for i in range(n_freq):
            U, S, Vh = svd(Gyy[i, :, :])
            singular_values[i, :] = S
            mode_shapes.append(U[:, 0])  # 第一奇异向量
        
        # 从奇异值峰值识别模态
        modes = []
        
        for mode_idx in range(min(3, n_dof)):  # 识别前3阶
            sv = singular_values[:, mode_idx]
            
            # 寻找峰值
            peaks, _ = find_peaks(sv, height=np.max(sv)*0.1, distance=10)
            
            for peak_idx in peaks:
                freq = f[peak_idx]
                
                # 半功率带宽法
                peak_sv = sv[peak_idx]
                half_power = peak_sv / np.sqrt(2)
                
                left_idx = peak_idx
                right_idx = peak_idx
                
                while left_idx > 0 and sv[left_idx] > half_power:
                    left_idx -= 1
                
                while right_idx < len(sv) - 1 and sv[right_idx] > half_power:
                    right_idx += 1
                
                if left_idx < peak_idx and right_idx > peak_idx:
                    damping = (f[right_idx] - f[left_idx]) / (2 * freq)
                else:
                    damping = 0.01
                
                modes.append({
                    'frequency': freq,
                    'damping': damping,
                    'mode_shape': mode_shapes[peak_idx],
                    'sv_peak': peak_sv
                })
        
        return modes, singular_values

def simulate_impact_test(n_dof=5, fs=1000, duration=2.0):
    """
    模拟锤击试验数据
    
    Parameters:
    -----------
    n_dof : int
        自由度数量（测点数）
    fs : float
        采样频率
    duration : float
        信号持续时间
    
    Returns:
    --------
    t : array
        时间向量
    force : array
        冲击力信号
    responses : array
        响应信号 (n_dof x n_samples)
    true_params : dict
        真实模态参数
    """
    np.random.seed(42)
    
    dt = 1.0 / fs
    t = np.arange(0, duration, dt)
    n_samples = len(t)
    
    # 真实模态参数
    true_freqs = np.array([10.5, 42.3, 95.8])  # Hz
    true_dampings = np.array([0.02, 0.015, 0.01])
    
    # 模态振型（简支梁）
    L = 1.0
    x_pos = np.linspace(0, L, n_dof)
    
    true_modes = np.zeros((n_dof, len(true_freqs)))
    for i, freq in enumerate(true_freqs):
        mode_num = i + 1
        true_modes[:, i] = np.sin(mode_num * np.pi * x_pos / L)
    
    # 归一化
    for i in range(len(true_freqs)):
        true_modes[:, i] /= np.max(np.abs(true_modes[:, i]))
    
    # 生成冲击力（半正弦脉冲）
    pulse_duration = 0.01  # 10ms
    force = np.zeros(n_samples)
    pulse_samples = int(pulse_duration * fs)
    
    if pulse_samples > 0:
        pulse_t = np.linspace(0, np.pi, pulse_samples)
        force[:pulse_samples] = 1000 * np.sin(pulse_t)
    
    # 计算响应（模态叠加）
    responses = np.zeros((n_dof, n_samples))
    
    for i in range(len(true_freqs)):
        omega = 2 * np.pi * true_freqs[i]
        zeta = true_dampings[i]
        omega_d = omega * np.sqrt(1 - zeta**2)
        
        # 单自由度响应
        h = np.zeros(n_samples)
        for j in range(1, n_samples):
            tau = j * dt
            h[j] = (1 / omega_d) * np.exp(-zeta * omega * tau) * np.sin(omega_d * tau)
        
        # 卷积计算响应
        for dof in range(n_dof):
            modal_force = true_modes[dof, i] * force
            resp = np.convolve(modal_force, h, mode='full')[:n_samples]
            responses[dof, :] += resp
    
    # 添加噪声
    noise_level = 0.05
    responses += noise_level * np.std(responses) * np.random.randn(n_dof, n_samples)
    
    true_params = {
        'frequencies': true_freqs,
        'dampings': true_dampings,
        'mode_shapes': true_modes
    }
    
    return t, force, responses, true_params

def case1_frequency_domain():
    """案例1：基于频域的模态参数识别"""
    print("="*60)
    print("主题063 - 案例1: 基于频域的模态参数识别")
    print("="*60)
    
    # 模拟锤击试验数据
    n_dof = 5
    fs = 1000
    t, force, responses, true_params = simulate_impact_test(n_dof, fs)
    
    print("\n真实模态参数:")
    for i, (f, zeta) in enumerate(zip(true_params['frequencies'], true_params['dampings'])):
        print(f"  第{i+1}阶: 频率 = {f:.2f} Hz, 阻尼比 = {zeta:.3f}")
    
    # 创建识别器
    identifier = FrequencyDomainIdentification(fs)
    
    # 方法1：峰值拾取法（使用第一个测点）
    print("\n方法1: 峰值拾取法")
    f, H, coherence = identifier.compute_frf(force, responses[0, :])
    
    modes_pp = identifier.peak_picking(f, H, peak_height=np.max(np.abs(H))*0.05)
    
    print(f"识别到 {len(modes_pp)} 个模态:")
    for i, mode in enumerate(modes_pp[:3]):
        print(f"  第{i+1}阶: 频率 = {mode['frequency']:.2f} Hz, "
              f"阻尼比 = {mode['damping']:.4f}")
    
    # 方法2：频域分解法
    print("\n方法2: 频域分解法 (FDD)")
    
    # 构建多参考频响函数矩阵
    H_matrix = np.zeros((len(f), n_dof, 1), dtype=complex)
    for dof in range(n_dof):
        _, H_dof, _ = identifier.compute_frf(force, responses[dof, :])
        H_matrix[:, dof, 0] = H_dof
    
    modes_fdd, singular_values = identifier.frequency_domain_decomposition(f, H_matrix)
    
    print(f"识别到 {len(modes_fdd)} 个模态:")
    for i, mode in enumerate(modes_fdd[:3]):
        print(f"  第{i+1}阶: 频率 = {mode['frequency']:.2f} Hz, "
              f"阻尼比 = {mode['damping']:.4f}")
    
    # 可视化
    fig, axes = plt.subplots(2, 3, figsize=(15, 10))
    
    # 1. 力信号和响应
    ax = axes[0, 0]
    ax.plot(t[:500], force[:500], 'b-', linewidth=1.5, label='冲击力')
    ax.set_xlabel('时间 (s)')
    ax.set_ylabel('力 (N)')
    ax.set_title('锤击力信号')
    ax.legend()
    ax.grid(True, alpha=0.3)
    
    # 2. 响应信号
    ax = axes[0, 1]
    for i in range(min(3, n_dof)):
        ax.plot(t[:1000], responses[i, :1000], linewidth=1, label=f'测点{i+1}')
    ax.set_xlabel('时间 (s)')
    ax.set_ylabel('位移 (m)')
    ax.set_title('结构响应')
    ax.legend(fontsize=8)
    ax.grid(True, alpha=0.3)
    
    # 3. 频响函数幅值
    ax = axes[0, 2]
    ax.semilogy(f, np.abs(H), 'b-', linewidth=1.5)
    
    # 标记识别的峰值
    for mode in modes_pp:
        ax.axvline(mode['frequency'], color='r', linestyle='--', alpha=0.7)
        ax.plot(mode['frequency'], np.abs(H)[mode['peak_idx']], 'ro', markersize=8)
    
    ax.set_xlabel('频率 (Hz)')
    ax.set_ylabel('|H(f)|')
    ax.set_title('频响函数 (峰值拾取)')
    ax.set_xlim(0, 150)
    ax.grid(True, alpha=0.3)
    
    # 4. 相干函数
    ax = axes[1, 0]
    ax.plot(f, coherence, 'g-', linewidth=1.5)
    ax.axhline(0.9, color='r', linestyle='--', alpha=0.7, label='阈值 0.9')
    ax.set_xlabel('频率 (Hz)')
    ax.set_ylabel('相干函数')
    ax.set_title('相干函数')
    ax.set_xlim(0, 150)
    ax.set_ylim(0, 1)
    ax.legend()
    ax.grid(True, alpha=0.3)
    
    # 5. 奇异值
    ax = axes[1, 1]
    for i in range(min(3, n_dof)):
        ax.semilogy(f, singular_values[:, i], linewidth=1.5, label=f'SV{i+1}')
    ax.set_xlabel('频率 (Hz)')
    ax.set_ylabel('奇异值')
    ax.set_title('FDD奇异值分解')
    ax.set_xlim(0, 150)
    ax.legend()
    ax.grid(True, alpha=0.3)
    
    # 6. 模态参数对比
    ax = axes[1, 2]
    
    x = np.arange(len(true_params['frequencies']))
    width = 0.25
    
    ax.bar(x - width, true_params['frequencies'], width, label='真实值', color='green', alpha=0.8)
    
    pp_freqs = [m['frequency'] for m in modes_pp[:3]]
    while len(pp_freqs) < 3:
        pp_freqs.append(0)
    ax.bar(x, pp_freqs, width, label='峰值拾取', color='blue', alpha=0.8)
    
    fdd_freqs = [m['frequency'] for m in modes_fdd[:3]]
    while len(fdd_freqs) < 3:
        fdd_freqs.append(0)
    ax.bar(x + width, fdd_freqs, width, label='FDD', color='red', alpha=0.8)
    
    ax.set_xlabel('模态阶数')
    ax.set_ylabel('频率 (Hz)')
    ax.set_title('频率识别结果对比')
    ax.set_xticks(x)
    ax.set_xticklabels(['第1阶', '第2阶', '第3阶'])
    ax.legend()
    ax.grid(True, alpha=0.3, axis='y')
    
    plt.tight_layout()
    plt.savefig('frequency_domain_identification.png', dpi=150, bbox_inches='tight')
    print("\n频域识别结果图已保存: frequency_domain_identification.png")
    plt.close()
    
    return identifier, modes_pp, modes_fdd

if __name__ == "__main__":
    identifier, modes_pp, modes_fdd = case1_frequency_domain()

结果分析：

1. 峰值拾取法：通过识别频响函数幅值的峰值确定固有频率，使用半功率带宽法估计阻尼比。该方法简单直观，但受频率分辨率限制。

2. 频域分解法（FDD）：通过奇异值分解分离各阶模态，在密集模态情况下表现更好。第一奇异值曲线清晰显示各阶模态频率。

3. 相干函数：用于评估测量质量，相干值接近1表示测量可靠。

4. 识别精度：两种方法都能较好地识别固有频率，但阻尼比估计存在一定误差，这是频域方法的共同特点。

---

3.2 案例2：基于时域的模态参数识别

问题描述：
利用结构的自由振动响应数据，采用特征系统实现算法（ERA）和Ibrahim时域法识别模态参数。

Python代码实现：

---

"""
主题063 - 案例2: 基于时域的模态参数识别
使用特征系统实现算法(ERA)和Ibrahim时域法识别结构模态参数
"""

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
import matplotlib
matplotlib.use('Agg')
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
from scipy.linalg import svd, eig
from scipy.signal import find_peaks


class TimeDomainIdentification:
    """时域模态参数识别"""
    
    def __init__(self, fs=1000, dt=None):
        """
        初始化
        
        Parameters:
        -----------
        fs : float
            采样频率 (Hz)
        dt : float
            采样时间间隔 (s)
        """
        if dt is None:
            self.fs = fs
            self.dt = 1.0 / fs
        else:
            self.dt = dt
            self.fs = 1.0 / dt
    
    def era_identification(self, impulse_response, alpha=None, beta=None, n_modes=None):
        """
        特征系统实现算法 (ERA)
        
        Parameters:
        -----------
        impulse_response : array
            脉冲响应函数 (n_outputs x n_samples) 或 (n_samples,)
        alpha : int
            Hankel矩阵行块数
        beta : int
            Hankel矩阵列块数
        n_modes : int
            要识别的模态数
        
        Returns:
        --------
        modes : list
            识别的模态参数
        singular_values : array
            奇异值
        """
        # 确保输入是2D数组
        if impulse_response.ndim == 1:
            impulse_response = impulse_response.reshape(1, -1)
        
        n_outputs, n_samples = impulse_response.shape
        
        # 设置默认参数
        if alpha is None:
            alpha = min(20, n_samples // 4)
        if beta is None:
            beta = min(20, n_samples // 4)
        
        if n_modes is None:
            n_modes = min(5, n_outputs)
        
        # 构建Hankel矩阵
        H0 = np.zeros((alpha * n_outputs, beta))
        H1 = np.zeros((alpha * n_outputs, beta))
        
        for i in range(alpha):
            for j in range(beta):
                if i + j < n_samples:
                    H0[i*n_outputs:(i+1)*n_outputs, j] = impulse_response[:, i+j]
                if i + j + 1 < n_samples:
                    H1[i*n_outputs:(i+1)*n_outputs, j] = impulse_response[:, i+j+1]
        
        # SVD分解
        U, S, Vt = svd(H0, full_matrices=False)
        
        # 确定系统阶数
        n_states = min(n_modes * 2, len(S))
        
        # 截断
        Un = U[:, :n_states]
        Sn = np.diag(S[:n_states])
        Vn = Vt[:n_states, :].T
        
        # 构建系统矩阵
        Sn_inv_sqrt = np.diag(1.0 / np.sqrt(S[:n_states] + 1e-10))
        
        A = Sn_inv_sqrt @ Un.T @ H1 @ Vn @ Sn_inv_sqrt
        
        # 观测矩阵
        C = Un[:n_outputs, :] @ np.sqrt(Sn + 1e-10 * np.eye(n_states))
        
        # 特征值分解
        eigenvalues, eigenvectors = eig(A)
        
        # 提取模态参数
        modes = []
        
        for i, ev in enumerate(eigenvalues):
            # 只考虑复共轭对的一半
            if np.imag(ev) >= 0:
                # 计算频率和阻尼
                lambda_log = np.log(ev + 1e-10)
                
                # 避免数值问题
                if np.isnan(lambda_log) or np.isinf(lambda_log):
                    continue
                
                omega_n = np.abs(lambda_log) / self.dt  # 固有频率 (rad/s)
                freq = omega_n / (2 * np.pi)  # 转换为Hz
                
                # 阻尼比
                if np.abs(lambda_log) > 1e-10:
                    zeta = -np.real(lambda_log) / np.abs(lambda_log)
                else:
                    zeta = 0.0
                
                # 模态振型
                mode_shape = C @ eigenvectors[:, i]
                
                # 只保留物理合理的模态
                if 0.1 < freq < self.fs / 2 and 0 < zeta < 0.5:
                    modes.append({
                        'frequency': freq,
                        'damping': zeta,
                        'mode_shape': mode_shape,
                        'eigenvalue': ev
                    })
        
        # 按频率排序
        modes.sort(key=lambda x: x['frequency'])
        
        return modes, S
    
    def itd_identification(self, free_response, n_modes=None, time_delay=None):
        """
        Ibrahim时域法 (ITD)
        
        Parameters:
        -----------
        free_response : array
            自由衰减响应 (n_dof x n_samples)
        n_modes : int
            要识别的模态数
        time_delay : int
            时间延迟点数
        
        Returns:
        --------
        modes : list
            识别的模态参数
        """
        if free_response.ndim == 1:
            free_response = free_response.reshape(1, -1)
        
        n_dof, n_samples = free_response.shape
        
        if n_modes is None:
            n_modes = min(n_dof, n_samples // 4)
        
        if time_delay is None:
            time_delay = max(1, n_samples // (2 * n_modes))
        
        # 构建响应矩阵
        X0 = np.zeros((n_dof, n_modes))
        X1 = np.zeros((n_dof, n_modes))
        
        for i in range(n_modes):
            if i * time_delay < n_samples:
                X0[:, i] = free_response[:, i * time_delay]
            if (i + 1) * time_delay < n_samples:
                X1[:, i] = free_response[:, (i + 1) * time_delay]
        
        # 构建增广矩阵（包含位移和速度）
        X_aug = np.vstack([X0, X1])
        Y_aug = np.vstack([X1, np.zeros_like(X1)])
        
        # 求解广义特征值问题
        try:
            eigenvalues, eigenvectors = eig(Y_aug, X_aug + 1e-10 * np.eye(X_aug.shape[0]))
        except:
            # 如果失败，使用标准特征值问题
            eigenvalues, eigenvectors = eig(X_aug.T @ Y_aug, X_aug.T @ X_aug + 1e-10 * np.eye(n_modes))
        
        # 提取模态参数
        modes = []
        
        for i, ev in enumerate(eigenvalues):
            if np.isnan(ev) or np.isinf(ev):
                continue
            
            # 计算频率和阻尼
            if np.abs(ev) > 1e-10:
                omega_d = np.angle(ev) / self.dt
                sigma = np.log(np.abs(ev)) / self.dt
                
                omega_n = np.sqrt(omega_d**2 + sigma**2)
                freq = omega_n / (2 * np.pi)
                zeta = -sigma / omega_n if omega_n > 1e-10 else 0
                
                # 模态振型
                if i < eigenvectors.shape[1]:
                    mode_shape = eigenvectors[:n_dof, i]
                    
                    # 只保留物理合理的模态
                    if 0.1 < freq < self.fs / 2 and 0 < zeta < 0.5:
                        modes.append({
                            'frequency': freq,
                            'damping': zeta,
                            'mode_shape': mode_shape,
                            'eigenvalue': ev
                        })
        
        # 按频率排序并去重
        modes.sort(key=lambda x: x['frequency'])
        
        # 去重（基于频率接近度）
        unique_modes = []
        for mode in modes:
            is_duplicate = False
            for existing in unique_modes:
                if abs(mode['frequency'] - existing['frequency']) / existing['frequency'] < 0.05:
                    is_duplicate = True
                    break
            if not is_duplicate:
                unique_modes.append(mode)
        
        return unique_modes[:n_modes]


def simulate_free_vibration(n_dof=5, fs=1000, duration=5.0):
    """
    模拟自由振动响应
    
    Parameters:
    -----------
    n_dof : int
        自由度数量
    fs : float
        采样频率
    duration : float
        信号持续时间
    
    Returns:
    --------
    t : array
        时间向量
    responses : array
        响应信号 (n_dof x n_samples)
    impulse_response : array
        脉冲响应函数
    true_params : dict
        真实模态参数
    """
    np.random.seed(42)
    
    dt = 1.0 / fs
    t = np.arange(0, duration, dt)
    n_samples = len(t)
    
    # 真实模态参数
    true_freqs = np.array([8.5, 34.2, 76.8])  # Hz
    true_dampings = np.array([0.025, 0.018, 0.012])
    
    # 模态振型
    L = 1.0
    x_pos = np.linspace(0, L, n_dof)
    
    true_modes = np.zeros((n_dof, len(true_freqs)))
    for i, freq in enumerate(true_freqs):
        mode_num = i + 1
        true_modes[:, i] = np.sin(mode_num * np.pi * x_pos / L)
        true_modes[:, i] /= np.max(np.abs(true_modes[:, i]))
    
    # 生成脉冲响应
    impulse_response = np.zeros((n_dof, n_samples))
    
    for i in range(len(true_freqs)):
        omega = 2 * np.pi * true_freqs[i]
        zeta = true_dampings[i]
        omega_d = omega * np.sqrt(1 - zeta**2)
        
        for dof in range(n_dof):
            h = (1 / omega_d) * np.exp(-zeta * omega * t) * np.sin(omega_d * t)
            impulse_response[dof, :] += true_modes[dof, i] * h
    
    # 添加噪声
    noise_level = 0.03
    impulse_response += noise_level * np.std(impulse_response) * np.random.randn(n_dof, n_samples)
    
    # 生成自由振动响应（初始位移激励）
    responses = np.zeros((n_dof, n_samples))
    
    for i in range(len(true_freqs)):
        omega = 2 * np.pi * true_freqs[i]
        zeta = true_dampings[i]
        omega_d = omega * np.sqrt(1 - zeta**2)
        
        # 初始条件
        x0 = 0.01  # 初始位移
        v0 = 0.0   # 初始速度
        
        for dof in range(n_dof):
            A = x0
            B = (v0 + zeta * omega * x0) / omega_d
            
            x = np.exp(-zeta * omega * t) * (A * np.cos(omega_d * t) + B * np.sin(omega_d * t))
            responses[dof, :] += true_modes[dof, i] * x
    
    # 添加噪声
    responses += noise_level * np.std(responses) * np.random.randn(n_dof, n_samples)
    
    true_params = {
        'frequencies': true_freqs,
        'dampings': true_dampings,
        'mode_shapes': true_modes
    }
    
    return t, responses, impulse_response, true_params


def case2_time_domain():
    """案例2：基于时域的模态参数识别"""
    print("="*60)
    print("主题063 - 案例2: 基于时域的模态参数识别")
    print("="*60)
    
    # 模拟自由振动数据
    n_dof = 5
    fs = 1000
    t, responses, impulse_response, true_params = simulate_free_vibration(n_dof, fs)
    
    print("\n真实模态参数:")
    for i, (f, zeta) in enumerate(zip(true_params['frequencies'], true_params['dampings'])):
        print(f"  第{i+1}阶: 频率 = {f:.2f} Hz, 阻尼比 = {zeta:.3f}")
    
    # 创建识别器
    identifier = TimeDomainIdentification(fs)
    
    # 方法1：ERA算法
    print("\n方法1: ERA (特征系统实现算法)")
    
    modes_era, singular_values = identifier.era_identification(
        impulse_response, alpha=15, beta=15, n_modes=3
    )
    
    print(f"识别到 {len(modes_era)} 个模态:")
    for i, mode in enumerate(modes_era[:3]):
        print(f"  第{i+1}阶: 频率 = {mode['frequency']:.2f} Hz, "
              f"阻尼比 = {mode['damping']:.4f}")
    
    # 方法2：ITD算法
    print("\n方法2: ITD (Ibrahim时域法)")
    
    modes_itd = identifier.itd_identification(responses, n_modes=3)
    
    print(f"识别到 {len(modes_itd)} 个模态:")
    for i, mode in enumerate(modes_itd[:3]):
        print(f"  第{i+1}阶: 频率 = {mode['frequency']:.2f} Hz, "
              f"阻尼比 = {mode['damping']:.4f}")
    
    # 可视化
    fig, axes = plt.subplots(2, 3, figsize=(15, 10))
    
    # 1. 脉冲响应
    ax = axes[0, 0]
    for i in range(min(3, n_dof)):
        ax.plot(t[:500], impulse_response[i, :500], linewidth=1, label=f'测点{i+1}')
    ax.set_xlabel('时间 (s)')
    ax.set_ylabel('响应')
    ax.set_title('脉冲响应函数')
    ax.legend(fontsize=8)
    ax.grid(True, alpha=0.3)
    
    # 2. 自由振动响应
    ax = axes[0, 1]
    for i in range(min(3, n_dof)):
        ax.plot(t[:1000], responses[i, :1000], linewidth=1, label=f'测点{i+1}')
    ax.set_xlabel('时间 (s)')
    ax.set_ylabel('位移 (m)')
    ax.set_title('自由振动响应')
    ax.legend(fontsize=8)
    ax.grid(True, alpha=0.3)
    
    # 3. ERA奇异值
    ax = axes[0, 2]
    ax.semilogy(singular_values, 'bo-', markersize=6)
    ax.axvline(x=6, color='r', linestyle='--', alpha=0.7, label='截断点')
    ax.set_xlabel('奇异值序号')
    ax.set_ylabel('奇异值')
    ax.set_title('ERA奇异值分解')
    ax.legend()
    ax.grid(True, alpha=0.3)
    
    # 4. 频率识别结果对比
    ax = axes[1, 0]
    
    x = np.arange(len(true_params['frequencies']))
    width = 0.25
    
    ax.bar(x - width, true_params['frequencies'], width, label='真实值', color='green', alpha=0.8)
    
    era_freqs = [m['frequency'] for m in modes_era[:3]]
    while len(era_freqs) < 3:
        era_freqs.append(0)
    ax.bar(x, era_freqs, width, label='ERA', color='blue', alpha=0.8)
    
    itd_freqs = [m['frequency'] for m in modes_itd[:3]]
    while len(itd_freqs) < 3:
        itd_freqs.append(0)
    ax.bar(x + width, itd_freqs, width, label='ITD', color='red', alpha=0.8)
    
    ax.set_xlabel('模态阶数')
    ax.set_ylabel('频率 (Hz)')
    ax.set_title('频率识别结果对比')
    ax.set_xticks(x)
    ax.set_xticklabels(['第1阶', '第2阶', '第3阶'])
    ax.legend()
    ax.grid(True, alpha=0.3, axis='y')
    
    # 5. 阻尼比识别结果对比
    ax = axes[1, 1]
    
    ax.bar(x - width, true_params['dampings'], width, label='真实值', color='green', alpha=0.8)
    
    era_dampings = [m['damping'] for m in modes_era[:3]]
    while len(era_dampings) < 3:
        era_dampings.append(0)
    ax.bar(x, era_dampings, width, label='ERA', color='blue', alpha=0.8)
    
    itd_dampings = [m['damping'] for m in modes_itd[:3]]
    while len(itd_dampings) < 3:
        itd_dampings.append(0)
    ax.bar(x + width, itd_dampings, width, label='ITD', color='red', alpha=0.8)
    
    ax.set_xlabel('模态阶数')
    ax.set_ylabel('阻尼比')
    ax.set_title('阻尼比识别结果对比')
    ax.set_xticks(x)
    ax.set_xticklabels(['第1阶', '第2阶', '第3阶'])
    ax.legend()
    ax.grid(True, alpha=0.3, axis='y')
    
    # 6. 方法特点
    ax = axes[1, 2]
    ax.text(0.5, 0.5, '时域方法特点:\n\n' +
            'ERA:\n- 基于脉冲响应\n- 数值稳定\n- 适合多输出系统\n\n' +
            'ITD:\n- 基于自由响应\n- 计算简单\n- 对噪声敏感', 
            ha='center', va='center', fontsize=10, transform=ax.transAxes)
    ax.set_xlim(0, 1)
    ax.set_ylim(0, 1)
    ax.axis('off')
    ax.set_title('方法特点对比')
    
    plt.tight_layout()
    plt.savefig('time_domain_identification.png', dpi=150, bbox_inches='tight')
    print("\n时域识别结果图已保存: time_domain_identification.png")
    plt.close()
    
    # 创建动画
    fig, axes = plt.subplots(1, 2, figsize=(14, 5))
    
    ax1 = axes[0]
    lines1 = []
    for i in range(min(3, n_dof)):
        line, = ax1.plot([], [], linewidth=1.5, label=f'测点{i+1}')
        lines1.append(line)
    ax1.set_xlim(0, t[500])
    ax1.set_ylim(np.min(impulse_response[:, :500]) * 1.1, 
                 np.max(impulse_response[:, :500]) * 1.1)
    ax1.set_xlabel('时间 (s)')
    ax1.set_ylabel('响应')
    ax1.set_title('脉冲响应动画')
    ax1.legend(fontsize=8)
    ax1.grid(True, alpha=0.3)
    
    ax2 = axes[1]
    lines2 = []
    for i in range(min(3, n_dof)):
        line, = ax2.plot([], [], linewidth=1.5, label=f'测点{i+1}')
        lines2.append(line)
    ax2.set_xlim(0, t[1000])
    ax2.set_ylim(np.min(responses[:, :1000]) * 1.1, 
                 np.max(responses[:, :1000]) * 1.1)
    ax2.set_xlabel('时间 (s)')
    ax2.set_ylabel('位移 (m)')
    ax2.set_title('自由振动响应动画')
    ax2.legend(fontsize=8)
    ax2.grid(True, alpha=0.3)
    
    def animate(frame):
        end_idx = min((frame + 1) * 10, 500)
        for i, line in enumerate(lines1):
            line.set_data(t[:end_idx], impulse_response[i, :end_idx])
        for i, line in enumerate(lines2):
            line.set_data(t[:end_idx*2], responses[i, :end_idx*2])
        return lines1 + lines2
    
    anim = FuncAnimation(fig, animate, frames=50, interval=100, blit=False)
    anim.save('time_domain_animation.gif', writer='pillow', fps=10, dpi=100)
    print("时域识别动画已保存: time_domain_animation.gif")
    plt.close()
    
    return identifier, modes_era, modes_itd


if __name__ == "__main__":
    identifier, modes_era, modes_itd = case2_time_domain()