主题061：结构动力学参数识别

1. 概述

1.1 参数识别的意义

结构动力学参数识别是指利用实测的振动响应数据，通过数学方法和算法，识别出结构的模态参数（固有频率、阻尼比、模态振型）和物理参数（质量、刚度、阻尼矩阵）的过程。参数识别是结构健康监测、模型修正和损伤识别的基础。

参数识别的主要目标：

• 模态参数识别：识别固有频率、阻尼比、模态振型
• 物理参数识别：识别质量、刚度、阻尼矩阵
• 模型修正：根据实测数据修正有限元模型
• 损伤评估：通过参数变化评估结构损伤
  
  
  
  
  
  
  
  

1.2 参数识别方法分类

1. 频域方法

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

2. 时域方法

• Ibrahim时域法（ITD）
• 特征系统实现算法（ERA）
• 随机减量法（RDT）

3. 时频域方法

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

4. 随机子空间方法

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

5. 贝叶斯方法

• 贝叶斯参数识别
• 马尔可夫链蒙特卡洛（MCMC）

1.3 参数识别的基本流程

实测数据 → 数据预处理 → 参数识别算法 → 模态参数 → 验证与修正

数据预处理：

• 滤波去噪
• 趋势项消除
• 重采样
• 归一化

模态参数：

• 固有频率 $f_n$
• 阻尼比 ${\zeta }_n$
• 模态振型 ${\varphi }_n$

2. 理论基础

2.1 频域参数识别方法

2.1.1 峰值拾取法

峰值拾取法是最简单的频域参数识别方法，通过识别频响函数（FRF）的峰值来确定模态参数。

固有频率识别：
$$f_n=f_{peak}$$

其中，$f_{peak}$为FRF幅值最大点对应的频率。

阻尼比识别（半功率带宽法）：
$${\zeta }_n=\frac{f_2-f_1}{2f_n}$$

其中，$f_1$和$f_2$为半功率点（幅值为峰值的KaTeX parse error: Expected 'EOF', got '}' at position 11: 1/\sqrt{2}}̲）。

模态振型识别：
$${\varphi }_n=\frac{H_{ij}(f_n)}{\max (|H_{ij}(f_n)|)}$$

其中，$H_{ij}(f_n)$为在固有频率处的频响函数值。

2.1.2 频域分解法（FDD）

FDD方法通过奇异值分解（SVD）识别密集模态。

功率谱密度矩阵：
$${\mathbf{S}}_{yy}(f)=\mathbf{H}(f){\mathbf{S}}_{xx}(f){\mathbf{H}}^H(f)$$

奇异值分解：
$${\mathbf{S}}_{yy}(f)=\mathbf{U}(f)\mathbf{\Sigma }(f){\mathbf{V}}^H(f)$$

模态参数识别：

• 固有频率：奇异值峰值对应的频率
• 模态振型：左奇异向量 ${\mathbf{u}}_1(f_n)$

2.2 时域参数识别方法

2.2.1 Ibrahim时域法（ITD）

ITD方法利用自由振动响应数据识别模态参数。

状态空间方程：
$$\dot{\mathbf{x}}=\mathbf{Ax}$$

其中，$\mathbf{A}$为系统矩阵。

特征值分解：
$$\mathbf{A}{\psi }_n={\lambda }_n{\psi }_n$$

模态参数：

• 固有频率：$f_n=\frac{|\text{Im}({\lambda }_n)|}{2\pi }$
• 阻尼比：${\zeta }_n=\frac{-\text{Re}({\lambda }_n)}{|{\lambda }_n|}$
• 模态振型：从特征向量提取

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

ERA方法利用脉冲响应函数（IRF）识别模态参数。

Hankel矩阵：
$$\mathbf{H}(k)=\left[\begin{array}{cccc}\mathbf{h}(k)& \mathbf{h}(k+1)& \cdots & \mathbf{h}(k+s-1)\\ \mathbf{h}(k+1)& \mathbf{h}(k+2)& \cdots & \mathbf{h}(k+s)\\ ⋮& ⋮& \ddots & ⋮\\ \mathbf{h}(k+r-1)& \mathbf{h}(k+r)& \cdots & \mathbf{h}(k+r+s-2)\end{array}\right]$$

奇异值分解：
$$\mathbf{H}(0)=\mathbf{P\Sigma }{\mathbf{Q}}^T$$

系统矩阵：
$$\mathbf{A}={\mathbf{\Sigma }}^{-1/2}{\mathbf{P}}^T\mathbf{H}(1)\mathbf{Q}{\mathbf{\Sigma }}^{-1/2}$$

2.3 随机子空间识别方法

2.3.1 协方差驱动随机子空间（Cov-SSI）

Cov-SSI方法利用输出数据的协方差矩阵识别模态参数。

协方差矩阵：
$${\mathbf{R}}_i=E[{\mathbf{y}}_{k+i}{\mathbf{y}}_k^T]$$

Toeplitz矩阵：
$${\mathbf{T}}_{1|i}=\left[\begin{array}{cccc}{\mathbf{R}}_i& {\mathbf{R}}_{i-1}& \cdots & {\mathbf{R}}_1\\ {\mathbf{R}}_{i+1}& {\mathbf{R}}_i& \cdots & {\mathbf{R}}_2\\ ⋮& ⋮& \ddots & ⋮\\ {\mathbf{R}}_{2i-1}& {\mathbf{R}}_{2i-2}& \cdots & {\mathbf{R}}_i\end{array}\right]$$

投影：
$${\mathbf{O}}_i={\mathbf{T}}_{1|i}{\mathbf{T}}_{0|i}^{-1}{\mathbf{Y}}_{0|i}$$

系统矩阵估计：

• 对观测矩阵进行SVD分解
• 提取系统矩阵$\mathbf{A}$和$\mathbf{C}$

2.3.2 稳定图

稳定图用于确定系统的真实模态阶数。

稳定准则：

• 频率稳定：$\left|\frac{f^{(i)}-f^{(i-1)}}{f^{(i)}}\right|<{ϵ}_f$
• 阻尼稳定：$\left|\frac{{\zeta }^{(i)}-{\zeta }^{(i-1)}}{{\zeta }^{(i)}}\right|<{ϵ}_{\zeta }$
• 振型稳定：$1-MAC({\varphi }^{(i)},{\varphi }^{(i-1)})<{ϵ}_{\varphi }$

2.4 贝叶斯参数识别方法

2.4.1 贝叶斯推断

贝叶斯方法将参数识别视为参数估计问题，利用贝叶斯定理更新参数的概率分布。

贝叶斯定理：
$$p(\mathbf{\theta }|\mathbf{D})=\frac{p(\mathbf{D}|\mathbf{\theta })p(\mathbf{\theta })}{p(\mathbf{D})}$$

其中：

• $p(\mathbf{\theta })$：先验分布
• $p(\mathbf{D}|\mathbf{\theta })$：似然函数
• $p(\mathbf{\theta }|\mathbf{D})$：后验分布

2.4.2 似然函数

对于模态参数识别，似然函数通常假设为高斯分布：

$$p(\mathbf{D}|\mathbf{\theta })=\prod _{k=1}^N\frac{1}{\sqrt{2\pi }\sigma }\exp \left(-\frac{(y_k-{\hat{y}}_k(\mathbf{\theta }){)}^2}{2{\sigma }^2}\right)$$

2.4.3 MCMC采样

马尔可夫链蒙特卡洛（MCMC）方法用于从后验分布中采样。

Metropolis-Hastings算法：

1. 初始化${\mathbf{\theta }}^{(0)}$
2. 对于$t=1,2,...,N$：
   • 从提议分布$q({\mathbf{\theta }}^{\ast }|{\mathbf{\theta }}^{(t-1)})$采样${\mathbf{\theta }}^{\ast }$
   • 计算接受概率：$\alpha =\min \left(1,\frac{p({\mathbf{\theta }}^{\ast }|\mathbf{D})q({\mathbf{\theta }}^{(t-1)}|{\mathbf{\theta }}^{\ast })}{p({\mathbf{\theta }}^{(t-1)}|\mathbf{D})q({\mathbf{\theta }}^{\ast }|{\mathbf{\theta }}^{(t-1)})}\right)$
   • 以概率$\alpha$接受${\mathbf{\theta }}^{\ast }$

3. 案例分析

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.signal import welch, find_peaks, csd
from scipy.linalg import svd


class FrequencyDomainIdentification:
    """频域参数识别类"""
    
    def __init__(self, fs):
        """
        初始化
        
        Parameters:
        -----------
        fs : float
            采样频率
        """
        self.fs = fs
    
    def compute_frf(self, force, response, window='hann', nperseg=None):
        """
        计算频响函数
        
        Parameters:
        -----------
        force : array
            力信号
        response : array
            响应信号
        window : str
            窗函数类型
        nperseg : int
            每段长度
        
        Returns:
        --------
        f : array
            频率
        frf : array
            频响函数
        """
        if nperseg is None:
            nperseg = min(256, len(force) // 4)
        
        # 计算自谱和互谱
        f, S_ff = welch(force, self.fs, window=window, nperseg=nperseg)
        f, S_fx = csd(force, response, self.fs, window=window, nperseg=nperseg)
        
        # 计算频响函数
        frf = S_fx / S_ff
        
        return f, frf
    
    def peak_picking(self, f, frf, n_modes=3, height=None, distance=None):
        """
        峰值拾取法识别模态参数
        
        Parameters:
        -----------
        f : array
            频率
        frf : array
            频响函数
        n_modes : int
            识别的模态数
        height : float
            峰值最小高度
        distance : int
            峰值间最小距离
        
        Returns:
        --------
        frequencies : array
            识别的固有频率
        damping_ratios : array
            识别的阻尼比
        mode_shapes : array
            识别的模态振型
        """
        # 找峰值
        magnitude = np.abs(frf)
        peaks, properties = find_peaks(magnitude, height=height, distance=distance)
        
        # 选择最大的n_modes个峰值
        if len(peaks) > n_modes:
            peak_heights = magnitude[peaks]
            top_indices = np.argsort(peak_heights)[-n_modes:]
            peaks = peaks[top_indices]
        
        frequencies = f[peaks]
        damping_ratios = []
        
        # 计算阻尼比（半功率带宽法）
        for peak in peaks:
            peak_mag = magnitude[peak]
            half_power = peak_mag / np.sqrt(2)
            
            # 找半功率点
            left_idx = peak
            right_idx = peak
            
            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
            
            f1, f2 = f[left_idx], f[right_idx]
            fn = f[peak]
            zeta = (f2 - f1) / (2 * fn)
            damping_ratios.append(zeta)
        
        # 模态振型（简化处理，使用FRF相位）
        mode_shapes = np.angle(frf[peaks])
        
        return frequencies, np.array(damping_ratios), mode_shapes, peaks
    
    def fdd_identification(self, data, n_modes=3, nperseg=None):
        """
        频域分解法识别模态参数
        
        Parameters:
        -----------
        data : array
            多通道响应数据 (n_channels, n_samples)
        n_modes : int
            识别的模态数
        nperseg : int
            每段长度
        
        Returns:
        --------
        frequencies : array
            识别的固有频率
        mode_shapes : array
            识别的模态振型
        singular_values : array
            奇异值
        f : array
            频率
        """
        n_channels = data.shape[0]
        if nperseg is None:
            nperseg = min(256, data.shape[1] // 4)
        
        # 计算功率谱密度矩阵
        f = np.fft.rfftfreq(nperseg, 1/self.fs)
        n_freqs = len(f)
        
        psd_matrix = np.zeros((n_freqs, n_channels, n_channels), dtype=complex)
        
        for i in range(n_channels):
            for j in range(n_channels):
                _, S_ij = csd(data[i], data[j], self.fs, window='hann', 
                              nperseg=nperseg, noverlap=nperseg//2)
                psd_matrix[:, i, j] = S_ij
        
        # 对每个频率进行SVD
        singular_values = np.zeros((n_freqs, n_channels))
        mode_shapes = np.zeros((n_freqs, n_channels, n_channels))
        
        for i in range(n_freqs):
            U, S, Vh = svd(psd_matrix[i].real)
            singular_values[i] = S
            mode_shapes[i] = U
        
        # 找奇异值峰值
        frequencies = []
        identified_mode_shapes = []
        
        for mode in range(min(n_modes, n_channels)):
            peaks, _ = find_peaks(singular_values[:, mode], 
                                  height=np.max(singular_values[:, mode])*0.1,
                                  distance=5)
            if len(peaks) > 0:
                # 选择最大的峰值
                main_peak = peaks[np.argmax(singular_values[peaks, mode])]
                frequencies.append(f[main_peak])
                identified_mode_shapes.append(mode_shapes[main_peak, :, mode])
        
        return np.array(frequencies), np.array(identified_mode_shapes), singular_values, f


def simulate_structure_response(fs, duration, frequencies, damping_ratios, mode_shapes, noise_level=0.05):
    """
    模拟结构振动响应
    
    Parameters:
    -----------
    fs : float
        采样频率
    duration : float
        持续时间
    frequencies : array
        固有频率
    damping_ratios : array
        阻尼比
    mode_shapes : array
        模态振型
    noise_level : float
        噪声水平
    
    Returns:
    --------
    t : array
        时间
    response : array
        响应信号
    force : array
        激励力
    """
    t = np.arange(0, duration, 1/fs)
    n_dof = mode_shapes.shape[0]
    
    # 生成随机激励
    force = np.random.randn(len(t))
    
    # 计算响应（模态叠加）
    response = np.zeros((n_dof, len(t)))
    
    for i, (fn, zeta, phi) in enumerate(zip(frequencies, damping_ratios, mode_shapes.T)):
        omega_n = 2 * np.pi * fn
        omega_d = omega_n * np.sqrt(1 - zeta**2)
        
        # 单自由度响应
        h = (1 / omega_d) * np.exp(-zeta * omega_n * t) * np.sin(omega_d * t)
        
        # 卷积计算响应
        modal_response = np.convolve(force, h, mode='same')[:len(t)]
        
        # 叠加到物理坐标
        response += np.outer(phi, modal_response)
    
    # 添加噪声
    response += noise_level * np.std(response) * np.random.randn(*response.shape)
    
    return t, response, force


def case1_frequency_domain_identification():
    """案例1：基于频域的参数识别"""
    print("="*60)
    print("主题061 - 案例1: 基于频域的参数识别")
    print("="*60)
    
    # 真实模态参数
    true_frequencies = np.array([5.0, 15.0, 30.0])  # Hz
    true_damping_ratios = np.array([0.02, 0.015, 0.01])
    true_mode_shapes = np.array([
        [1.0, 1.0, 1.0],
        [0.8, 0.0, -0.8],
        [0.5, -0.8, 0.5]
    ])
    
    print("\n真实模态参数:")
    for i, (f, zeta) in enumerate(zip(true_frequencies, true_damping_ratios)):
        print(f"  第{i+1}阶: 频率 = {f:.2f} Hz, 阻尼比 = {zeta:.3f}")
    
    # 模拟响应
    fs = 200  # 采样频率
    duration = 50  # 持续时间
    t, response, force = simulate_structure_response(
        fs, duration, true_frequencies, true_damping_ratios, true_mode_shapes)
    
    print(f"\n模拟数据:")
    print(f"  采样频率: {fs} Hz")
    print(f"  持续时间: {duration} s")
    print(f"  数据点数: {len(t)}")
    
    # 频域参数识别
    identifier = FrequencyDomainIdentification(fs)
    
    # 1. 峰值拾取法
    f, frf = identifier.compute_frf(force, response[0])
    frequencies_pp, damping_pp, modes_pp, peaks = identifier.peak_picking(
        f, frf, n_modes=3, height=0.001, distance=10)
    
    print("\n峰值拾取法识别结果:")
    for i, (f_id, zeta_id) in enumerate(zip(frequencies_pp, damping_pp)):
        print(f"  第{i+1}阶: 频率 = {f_id:.2f} Hz, 阻尼比 = {zeta_id:.4f}")
    
    # 2. FDD方法
    frequencies_fdd, modes_fdd, sv, f_svd = identifier.fdd_identification(response, n_modes=3)
    
    print("\nFDD方法识别结果:")
    for i, f_id in enumerate(frequencies_fdd):
        print(f"  第{i+1}阶: 频率 = {f_id:.2f} Hz")
    
    # 可视化
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # 1. 时域响应
    ax = axes[0, 0]
    ax.plot(t[:1000], response[0, :1000], 'b-', linewidth=0.8)
    ax.set_xlabel('时间 (s)')
    ax.set_ylabel('响应')
    ax.set_title('结构振动响应（前5秒）')
    ax.grid(True, alpha=0.3)
    
    # 2. 频响函数和峰值
    ax = axes[0, 1]
    ax.semilogy(f, np.abs(frf), 'b-', linewidth=1, label='FRF幅值')
    if len(peaks) > 0:
        ax.plot(frequencies_pp, np.abs(frf[peaks]), 'ro', markersize=8, label='识别峰值')
    for i, f_true in enumerate(true_frequencies):
        ax.axvline(x=f_true, color='green', linestyle='--', alpha=0.5, 
                   label='真实频率' if i == 0 else '')
    ax.set_xlabel('频率 (Hz)')
    ax.set_ylabel('FRF幅值')
    ax.set_title('频响函数和峰值拾取')
    ax.legend()
    ax.grid(True, alpha=0.3)
    ax.set_xlim([0, 50])
    
    # 3. 奇异值分解结果
    ax = axes[1, 0]
    for i in range(min(3, sv.shape[1])):
        ax.semilogy(f_svd, sv[:, i], linewidth=1.5, label=f'SV {i+1}')
    ax.set_xlabel('频率 (Hz)')
    ax.set_ylabel('奇异值')
    ax.set_title('FDD奇异值分解')
    ax.legend()
    ax.grid(True, alpha=0.3)
    ax.set_xlim([0, 50])
    
    # 4. 频率识别对比
    ax = axes[1, 1]
    x = np.arange(len(true_frequencies))
    width = 0.25
    ax.bar(x - width, true_frequencies, width, label='真实值', color='green', alpha=0.8)
    
    # 峰值拾取结果
    pp_freqs_plot = []
    for i in range(len(true_frequencies)):
        if i < len(frequencies_pp):
            pp_freqs_plot.append(frequencies_pp[i])
        else:
            pp_freqs_plot.append(0)
    ax.bar(x, pp_freqs_plot, width, label='峰值拾取', color='blue', alpha=0.8)
    
    # FDD结果
    if len(frequencies_fdd) == len(true_frequencies):
        ax.bar(x + width, frequencies_fdd, 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([f'第{i+1}阶' for i in range(len(true_frequencies))])
    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()
    
    # 创建动画
    fig, axes = plt.subplots(1, 2, figsize=(14, 5))
    
    # 时域响应动画
    ax1 = axes[0]
    line1, = ax1.plot([], [], 'b-', linewidth=0.8)
    ax1.set_xlim(0, 5)
    ax1.set_ylim(np.min(response[0, :1000]), np.max(response[0, :1000]))
    ax1.set_xlabel('时间 (s)')
    ax1.set_ylabel('响应')
    ax1.set_title('时域响应动画')
    ax1.grid(True, alpha=0.3)
    
    # 频响函数动画
    ax2 = axes[1]
    line2, = ax2.semilogy([], [], 'b-', linewidth=1)
    scatter2 = ax2.scatter([], [], c='red', s=50)
    ax2.set_xlim(0, 50)
    ax2.set_ylim(np.min(np.abs(frf)) * 0.5, np.max(np.abs(frf)) * 2)
    ax2.set_xlabel('频率 (Hz)')
    ax2.set_ylabel('FRF幅值')
    ax2.set_title('频响函数峰值拾取动画')
    ax2.grid(True, alpha=0.3)
    
    def animate(frame):
        # 更新时域响应
        end_idx = min((frame + 1) * 50, 1000)
        line1.set_data(t[:end_idx], response[0, :end_idx])
        
        # 更新频响函数（逐步显示）
        freq_idx = min((frame + 1) * 5, len(f))
        line2.set_data(f[:freq_idx], np.abs(frf[:freq_idx]))
        
        # 标记峰值
        peak_in_range = peaks[peaks < freq_idx]
        if len(peak_in_range) > 0:
            scatter2.set_offsets(np.column_stack([f[peak_in_range], np.abs(frf[peak_in_range])]))
        
        return [line1, line2, scatter2]
    
    anim = FuncAnimation(fig, animate, frames=20, interval=200, blit=False)
    anim.save('frequency_identification_animation.gif', writer='pillow', fps=5, dpi=100)
    print("频域参数识别动画已保存: frequency_identification_animation.gif")
    plt.close()
    
    return identifier, frequencies_pp, damping_pp


if __name__ == "__main__":
    identifier, freqs, damping = case1_frequency_domain_identification()

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


class TimeDomainIdentification:
    """时域参数识别类"""
    
    def __init__(self, fs):
        """
        初始化
        
        Parameters:
        -----------
        fs : float
            采样频率
        """
        self.fs = fs
        self.dt = 1.0 / fs
    
    def era_identification(self, irf, n_modes=3, r=None, s=None):
        """
        特征系统实现算法（ERA）识别模态参数
        
        Parameters:
        -----------
        irf : array
            脉冲响应函数 (n_outputs, n_inputs, n_samples)
        n_modes : int
            识别的模态数
        r : int
            Hankel矩阵行数
        s : int
            Hankel矩阵列数
        
        Returns:
        --------
        frequencies : array
            识别的固有频率
        damping_ratios : array
            识别的阻尼比
        mode_shapes : array
            识别的模态振型
        """
        n_outputs, n_inputs, n_samples = irf.shape
        
        if r is None:
            r = min(20, n_samples // 4)
        if s is None:
            s = min(20, n_samples // 4)
        
        # 构建Hankel矩阵
        H0 = np.zeros((r * n_outputs, s * n_inputs))
        H1 = np.zeros((r * n_outputs, s * n_inputs))
        
        for i in range(r):
            for j in range(s):
                H0[i*n_outputs:(i+1)*n_outputs, j*n_inputs:(j+1)*n_inputs] = irf[:, :, i+j]
                H1[i*n_outputs:(i+1)*n_outputs, j*n_inputs:(j+1)*n_inputs] = irf[:, :, i+j+1]
        
        # SVD分解
        U, Sigma, Vt = svd(H0, full_matrices=False)
        
        # 截断到n_modes阶
        n = min(n_modes * 2, len(Sigma))
        Un = U[:, :n]
        Sigman = np.diag(Sigma[:n])
        Vtn = Vt[:n, :]
        
        # 计算系统矩阵
        A = np.linalg.multi_dot([
            np.linalg.inv(np.sqrt(Sigman)),
            Un.T,
            H1,
            Vtn.T,
            np.linalg.inv(np.sqrt(Sigman))
        ])
        
        # 特征值分解
        eigenvalues, eigenvectors = eig(A)
        
        # 提取模态参数
        frequencies = []
        damping_ratios = []
        
        for ev in eigenvalues:
            if np.imag(ev) != 0:
                # 连续时间特征值
                lambda_c = np.log(ev) / self.dt
                
                # 固有频率
                fn = np.abs(np.imag(lambda_c)) / (2 * np.pi)
                
                # 阻尼比
                zeta = -np.real(lambda_c) / np.abs(lambda_c)
                
                if fn > 0 and fn < self.fs / 2 and zeta > 0 and zeta < 1:
                    frequencies.append(fn)
                    damping_ratios.append(zeta)
        
        # 去重并排序
        unique_modes = []
        for f, z in zip(frequencies, damping_ratios):
            is_unique = True
            for uf, uz in unique_modes:
                if abs(f - uf) < 0.5:  # 频率差小于0.5Hz认为是同一模态
                    is_unique = False
                    break
            if is_unique:
                unique_modes.append((f, z))
        
        unique_modes.sort(key=lambda x: x[0])
        
        frequencies = np.array([m[0] for m in unique_modes[:n_modes]])
        damping_ratios = np.array([m[1] for m in unique_modes[:n_modes]])
        
        # 模态振型（简化处理）
        mode_shapes = np.random.randn(n_outputs, len(frequencies))
        
        return frequencies, damping_ratios, mode_shapes
    
    def random_decrement(self, response, trigger_level=None, n_segments=100):
        """
        随机减量法提取自由衰减响应
        
        Parameters:
        -----------
        response : array
            随机响应信号
        trigger_level : float
            触发水平
        n_segments : int
            分段数
        
        Returns:
        --------
        rd_signature : array
            随机减量特征
        time : array
            时间向量
        """
        if trigger_level is None:
            trigger_level = 0.5 * np.std(response)
        
        segment_length = len(response) // n_segments
        rd_segments = []
        
        for i in range(n_segments):
            start = i * segment_length
            end = start + segment_length
            segment = response[start:end]
            
            # 找触发点（向上穿越触发水平）
            trigger_indices = []
            for j in range(1, len(segment)):
                if segment[j-1] < trigger_level and segment[j] >= trigger_level:
                    trigger_indices.append(j)
            
            # 提取衰减段
            for idx in trigger_indices:
                if idx + 100 < len(segment):
                    rd_segments.append(segment[idx:idx+100] - trigger_level)
        
        # 平均
        if len(rd_segments) > 0:
            min_len = min(len(s) for s in rd_segments)
            rd_segments = [s[:min_len] for s in rd_segments]
            rd_signature = np.mean(rd_segments, axis=0)
        else:
            rd_signature = np.zeros(100)
        
        time = np.arange(len(rd_signature)) * self.dt
        
        return rd_signature, time
    
    def itd_identification(self, free_response, n_modes=3):
        """
        Ibrahim时域法识别模态参数
        
        Parameters:
        -----------
        free_response : array
            自由振动响应 (n_dof, n_samples)
        n_modes : int
            识别的模态数
        
        Returns:
        --------
        frequencies : array
            识别的固有频率
        damping_ratios : array
            识别的阻尼比
        mode_shapes : array
            识别的模态振型
        """
        n_dof, n_samples = free_response.shape
        
        # 构建数据矩阵
        Y0 = free_response[:, :-2]
        Y1 = free_response[:, 1:-1]
        Y2 = free_response[:, 2:]
        
        # 构建扩展矩阵
        H0 = np.vstack([Y0, Y1])
        H1 = np.vstack([Y1, Y2])
        
        # 求解特征值问题
        A = np.dot(H1, np.linalg.pinv(H0))
        eigenvalues, eigenvectors = eig(A)
        
        # 提取模态参数
        frequencies = []
        damping_ratios = []
        
        for ev in eigenvalues:
            if np.imag(ev) != 0 and np.abs(ev) < 1:
                # 固有频率
                fn = np.arctan2(np.imag(ev), np.real(ev)) / (2 * np.pi * self.dt)
                if fn < 0:
                    fn = -fn
                
                # 阻尼比
                zeta = -np.log(np.abs(ev)) / (2 * np.pi * fn * self.dt)
                
                if fn > 0 and fn < self.fs / 2 and zeta > 0 and zeta < 1:
                    frequencies.append(fn)
                    damping_ratios.append(zeta)
        
        # 去重
        unique_modes = []
        for f, z in zip(frequencies, damping_ratios):
            is_unique = True
            for uf, uz in unique_modes:
                if abs(f - uf) < 0.5:
                    is_unique = False
                    break
            if is_unique:
                unique_modes.append((f, z))
        
        unique_modes.sort(key=lambda x: x[0])
        
        frequencies = np.array([m[0] for m in unique_modes[:n_modes]])
        damping_ratios = np.array([m[1] for m in unique_modes[:n_modes]])
        
        # 模态振型
        mode_shapes = eigenvectors[:n_dof, :len(frequencies)].real
        
        return frequencies, damping_ratios, mode_shapes


def simulate_impulse_response(fs, duration, frequencies, damping_ratios, mode_shapes):
    """
    模拟脉冲响应函数
    
    Parameters:
    -----------
    fs : float
        采样频率
    duration : float
        持续时间
    frequencies : array
        固有频率
    damping_ratios : array
        阻尼比
    mode_shapes : array
        模态振型
    
    Returns:
    --------
    t : array
        时间
    irf : array
        脉冲响应函数
    """
    t = np.arange(0, duration, 1/fs)
    n_dof = mode_shapes.shape[0]
    
    irf = np.zeros((n_dof, 1, len(t)))
    
    for i, (fn, zeta, phi) in enumerate(zip(frequencies, damping_ratios, mode_shapes.T)):
        omega_n = 2 * np.pi * fn
        omega_d = omega_n * np.sqrt(1 - zeta**2)
        
        # 单自由度脉冲响应
        h = (1 / omega_d) * np.exp(-zeta * omega_n * t) * np.sin(omega_d * t)
        
        # 叠加到物理坐标
        irf[:, 0, :] += np.outer(phi, h)
    
    return t, irf


def simulate_free_vibration(fs, duration, frequencies, damping_ratios, mode_shapes, initial_displacement):
    """
    模拟自由振动响应
    
    Parameters:
    -----------
    fs : float
        采样频率
    duration : float
        持续时间
    frequencies : array
        固有频率
    damping_ratios : array
        阻尼比
    mode_shapes : array
        模态振型
    initial_displacement : array
        初始位移
    
    Returns:
    --------
    t : array
        时间
    response : array
        自由振动响应
    """
    t = np.arange(0, duration, 1/fs)
    n_dof = mode_shapes.shape[0]
    
    response = np.zeros((n_dof, len(t)))
    
    # 将初始位移转换到模态坐标
    modal_disp = np.linalg.lstsq(mode_shapes, initial_displacement, rcond=None)[0]
    
    for i, (fn, zeta, phi) in enumerate(zip(frequencies, damping_ratios, mode_shapes.T)):
        omega_n = 2 * np.pi * fn
        omega_d = omega_n * np.sqrt(1 - zeta**2)
        
        # 初始条件响应
        A = modal_disp[i]
        y = A * np.exp(-zeta * omega_n * t) * np.cos(omega_d * t)
        
        # 叠加到物理坐标
        response += np.outer(phi, y)
    
    return t, response


def case2_time_domain_identification():
    """案例2：基于时域的参数识别"""
    print("="*60)
    print("主题061 - 案例2: 基于时域的参数识别")
    print("="*60)
    
    # 真实模态参数
    true_frequencies = np.array([5.0, 15.0, 30.0])  # Hz
    true_damping_ratios = np.array([0.02, 0.015, 0.01])
    true_mode_shapes = np.array([
        [1.0, 1.0, 1.0],
        [0.8, 0.0, -0.8],
        [0.5, -0.8, 0.5]
    ])
    
    print("\n真实模态参数:")
    for i, (f, zeta) in enumerate(zip(true_frequencies, true_damping_ratios)):
        print(f"  第{i+1}阶: 频率 = {f:.2f} Hz, 阻尼比 = {zeta:.3f}")
    
    fs = 200  # 采样频率
    duration = 10  # 持续时间
    
    # 1. ERA方法（基于脉冲响应）
    print("\n1. ERA方法识别:")
    t_irf, irf = simulate_impulse_response(
        fs, duration, true_frequencies, true_damping_ratios, true_mode_shapes)
    
    identifier = TimeDomainIdentification(fs)
    freq_era, damping_era, modes_era = identifier.era_identification(irf, n_modes=3)
    
    print("  ERA识别结果:")
    for i, (f, zeta) in enumerate(zip(freq_era, damping_era)):
        print(f"    第{i+1}阶: 频率 = {f:.2f} Hz, 阻尼比 = {zeta:.4f}")
    
    # 2. ITD方法（基于自由振动）
    print("\n2. ITD方法识别:")
    initial_disp = np.array([1.0, 0.5, 0.3])
    t_free, free_response = simulate_free_vibration(
        fs, duration, true_frequencies, true_damping_ratios, true_mode_shapes, initial_disp)
    
    # 添加噪声
    free_response += 0.01 * np.std(free_response) * np.random.randn(*free_response.shape)
    
    freq_itd, damping_itd, modes_itd = identifier.itd_identification(free_response, n_modes=3)
    
    print("  ITD识别结果:")
    for i, (f, zeta) in enumerate(zip(freq_itd, damping_itd)):
        print(f"    第{i+1}阶: 频率 = {f:.2f} Hz, 阻尼比 = {zeta:.4f}")
    
    # 3. 随机减量法
    print("\n3. 随机减量法:")
    # 模拟随机响应
    t_random = np.arange(0, 20, 1/fs)
    random_response = np.random.randn(len(t_random))
    for fn, zeta, phi in zip(true_frequencies, true_damping_ratios, true_mode_shapes.T):
        omega_n = 2 * np.pi * fn
        omega_d = omega_n * np.sqrt(1 - zeta**2)
        h = (1 / omega_d) * np.exp(-zeta * omega_n * t_random) * np.sin(omega_d * t_random)
        modal_response = np.convolve(random_response, h, mode='same')[:len(t_random)]
        random_response += 0.3 * modal_response
    
    rd_signature, rd_time = identifier.random_decrement(random_response, n_segments=50)
    print(f"  提取的RD特征长度: {len(rd_signature)}")
    
    # 可视化
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # 1. 脉冲响应函数
    ax = axes[0, 0]
    for i in range(irf.shape[0]):
        ax.plot(t_irf, irf[i, 0, :], label=f'自由度{i+1}', alpha=0.7)
    ax.set_xlabel('时间 (s)')
    ax.set_ylabel('响应')
    ax.set_title('脉冲响应函数')
    ax.legend()
    ax.grid(True, alpha=0.3)
    
    # 2. 自由振动响应
    ax = axes[0, 1]
    for i in range(free_response.shape[0]):
        ax.plot(t_free, free_response[i, :], label=f'自由度{i+1}', alpha=0.7)
    ax.set_xlabel('时间 (s)')
    ax.set_ylabel('响应')
    ax.set_title('自由振动响应')
    ax.legend()
    ax.grid(True, alpha=0.3)
    
    # 3. 随机减量特征
    ax = axes[1, 0]
    ax.plot(rd_time, rd_signature, 'b-', linewidth=1.5)
    ax.set_xlabel('时间 (s)')
    ax.set_ylabel('RD特征')
    ax.set_title('随机减量特征')
    ax.grid(True, alpha=0.3)
    
    # 4. 频率识别对比
    ax = axes[1, 1]
    x = np.arange(len(true_frequencies))
    width = 0.25
    
    ax.bar(x - width, true_frequencies, width, label='真实值', color='green', alpha=0.8)
    
    # ERA结果
    era_freqs_plot = []
    for i in range(len(true_frequencies)):
        if i < len(freq_era):
            era_freqs_plot.append(freq_era[i])
        else:
            era_freqs_plot.append(0)
    ax.bar(x, era_freqs_plot, width, label='ERA', color='blue', alpha=0.8)
    
    # ITD结果
    itd_freqs_plot = []
    for i in range(len(true_frequencies)):
        if i < len(freq_itd):
            itd_freqs_plot.append(freq_itd[i])
        else:
            itd_freqs_plot.append(0)
    ax.bar(x + width, itd_freqs_plot, 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([f'第{i+1}阶' for i in range(len(true_frequencies))])
    ax.legend()
    ax.grid(True, alpha=0.3, axis='y')
    
    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(irf.shape[0]):
        line, = ax1.plot([], [], label=f'自由度{i+1}', alpha=0.7)
        lines1.append(line)
    ax1.set_xlim(0, duration)
    ax1.set_ylim(np.min(irf) * 1.1, np.max(irf) * 1.1)
    ax1.set_xlabel('时间 (s)')
    ax1.set_ylabel('响应')
    ax1.set_title('脉冲响应函数动画')
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    
    # RD特征动画
    ax2 = axes[1]
    line2, = ax2.plot([], [], 'b-', linewidth=1.5)
    ax2.set_xlim(0, rd_time[-1])
    ax2.set_ylim(np.min(rd_signature) * 1.1, np.max(rd_signature) * 1.1)
    ax2.set_xlabel('时间 (s)')
    ax2.set_ylabel('RD特征')
    ax2.set_title('随机减量特征动画')
    ax2.grid(True, alpha=0.3)
    
    def animate(frame):
        end_idx = min((frame + 1) * 50, len(t_irf))
        
        # 更新脉冲响应
        for i, line in enumerate(lines1):
            line.set_data(t_irf[:end_idx], irf[i, 0, :end_idx])
        
        # 更新RD特征
        rd_end_idx = min((frame + 1) * 5, len(rd_time))
        line2.set_data(rd_time[:rd_end_idx], rd_signature[:rd_end_idx])
        
        return lines1 + [line2]
    
    anim = FuncAnimation(fig, animate, frames=40, interval=150, blit=False)
    anim.save('time_domain_animation.gif', writer='pillow', fps=6, dpi=100)
    print("时域参数识别动画已保存: time_domain_animation.gif")
    plt.close()
    
    return identifier, freq_era, damping_era


if __name__ == "__main__":
    identifier, freqs, damping = case2_time_domain_identification()

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, block_diag


class StochasticSubspaceIdentification:
    """随机子空间识别类"""
    
    def __init__(self, fs):
        """
        初始化
        
        Parameters:
        -----------
        fs : float
            采样频率
        """
        self.fs = fs
        self.dt = 1.0 / fs
    
    def cov_ssi(self, data, i, n_modes=3, max_order=20):
        """
        协方差驱动随机子空间识别
        
        Parameters:
        -----------
        data : array
            输出数据 (n_channels, n_samples)
        i : int
            块行数/列数
        n_modes : int
            识别的模态数
        max_order : int
            最大模型阶数
        
        Returns:
        --------
        frequencies : array
            识别的固有频率
        damping_ratios : array
            识别的阻尼比
        mode_shapes : array
            识别的模态振型
        stability_diagram : dict
            稳定图数据
        """
        n_channels, n_samples = data.shape
        
        # 计算协方差矩阵
        R = []
        for k in range(2 * i):
            if k < n_samples:
                R_k = np.dot(data[:, k:], data[:, :n_samples-k].T) / (n_samples - k)
            else:
                R_k = np.zeros((n_channels, n_channels))
            R.append(R_k)
        
        # 构建块Toeplitz矩阵
        T = np.zeros((i * n_channels, i * n_channels))
        for row in range(i):
            for col in range(i):
                T[row*n_channels:(row+1)*n_channels, col*n_channels:(col+1)*n_channels] = R[abs(row - col) + 1]
        
        # SVD分解
        U, S, Vt = svd(T, full_matrices=False)
        
        # 稳定图分析
        stability_diagram = {
            'orders': [],
            'frequencies': [],
            'damping_ratios': [],
            'stable': []
        }
        
        all_frequencies = []
        all_damping = []
        
        for order in range(2, max_order + 1, 2):
            # 截断到当前阶数
            n = order
            Un = U[:, :n]
            Sn = np.diag(S[:n])
            Vtn = Vt[:n, :]
            
            # 观测矩阵
            Oi = np.dot(Un, np.sqrt(Sn))
            
            # 提取系统矩阵
            C = Oi[:n_channels, :]
            
            # 移位后的观测矩阵
            Oi_shift = Oi[n_channels:, :]
            Oi_trunc = Oi[:-n_channels, :]
            
            # 系统矩阵A
            A = np.linalg.lstsq(Oi_trunc, Oi_shift, rcond=None)[0].T
            
            # 特征值分解
            eigenvalues, eigenvectors = eig(A)
            
            # 提取模态参数
            frequencies = []
            damping_ratios = []
            
            for ev in eigenvalues:
                if np.abs(ev) < 1 and np.imag(ev) != 0:
                    # 连续时间特征值
                    lambda_c = np.log(ev) / self.dt
                    
                    # 固有频率
                    fn = np.abs(np.imag(lambda_c)) / (2 * np.pi)
                    
                    # 阻尼比
                    zeta = -np.real(lambda_c) / np.abs(lambda_c)
                    
                    if fn > 0.1 and fn < self.fs / 2 and zeta > 0 and zeta < 1:
                        frequencies.append(fn)
                        damping_ratios.append(zeta)
            
            all_frequencies.append(frequencies)
            all_damping.append(damping_ratios)
            
            stability_diagram['orders'].append(order)
            stability_diagram['frequencies'].append(frequencies)
            stability_diagram['damping_ratios'].append(damping_ratios)
        
        # 判断稳定性
        stability_diagram['stable'] = self._check_stability(all_frequencies, all_damping)
        
        # 选择最稳定的模态
        if len(all_frequencies) > 0:
            # 使用最高阶的结果
            final_freqs = all_frequencies[-1][:n_modes]
            final_damping = all_damping[-1][:n_modes]
            
            # 如果不够，用0填充
            while len(final_freqs) < n_modes:
                final_freqs = np.append(final_freqs, 0)
                final_damping = np.append(final_damping, 0)
            
            frequencies = np.array(final_freqs)
            damping_ratios = np.array(final_damping)
        else:
            frequencies = np.zeros(n_modes)
            damping_ratios = np.zeros(n_modes)
        
        # 模态振型（简化）
        mode_shapes = np.random.randn(n_channels, len(frequencies))
        
        return frequencies, damping_ratios, mode_shapes, stability_diagram
    
    def _check_stability(self, all_frequencies, all_damping, 
                         freq_tol=0.01, damping_tol=0.05):
        """
        检查模态稳定性
        
        Parameters:
        -----------
        all_frequencies : list
            各阶模型的频率
        all_damping : list
            各阶模型的阻尼比
        freq_tol : float
            频率容差
        damping_tol : float
            阻尼比容差
        
        Returns:
        --------
        stable : list
            稳定性标记
        """
        stable = []
        
        for i in range(len(all_frequencies)):
            stable_i = []
            for j, (f, d) in enumerate(zip(all_frequencies[i], all_damping[i])):
                is_stable = False
                
                # 与前一阶比较
                if i > 0:
                    for f_prev, d_prev in zip(all_frequencies[i-1], all_damping[i-1]):
                        if (abs(f - f_prev) / f < freq_tol and 
                            abs(d - d_prev) / d < damping_tol):
                            is_stable = True
                            break
                
                stable_i.append(is_stable)
            
            stable.append(stable_i)
        
        return stable
    
    def data_ssi(self, data, i, n_modes=3, max_order=20):
        """
        数据驱动随机子空间识别
        
        Parameters:
        -----------
        data : array
            输出数据 (n_channels, n_samples)
        i : int
            块行数
        n_modes : int
            识别的模态数
        max_order : int
            最大模型阶数
        
        Returns:
        --------
        frequencies : array
            识别的固有频率
        damping_ratios : array
            识别的阻尼比
        mode_shapes : array
            识别的模态振型
        """
        n_channels, n_samples = data.shape
        
        # 构建Hankel矩阵
        Yp = np.zeros((i * n_channels, n_samples - 2*i + 1))
        Yf = np.zeros((i * n_channels, n_samples - 2*i + 1))
        
        for k in range(i):
            Yp[k*n_channels:(k+1)*n_channels, :] = data[:, k:k+n_samples-2*i+1]
            Yf[k*n_channels:(k+1)*n_channels, :] = data[:, i+k:i+k+n_samples-2*i+1]
        
        # LQ分解
        L = np.linalg.cholesky(np.dot(Yp, Yp.T) + 1e-10 * np.eye(Yp.shape[0]))
        Q = np.linalg.solve(L, np.dot(Yp, Yf.T))
        
        # SVD分解
        U, S, Vt = svd(Q, full_matrices=False)
        
        # 选择模型阶数
        order = min(n_modes * 2, max_order, len(S))
        Un = U[:, :order]
        Sn = np.diag(S[:order])
        
        # 观测矩阵
        Oi = np.dot(Un, np.sqrt(Sn))
        
        # 提取系统矩阵
        C = Oi[:n_channels, :]
        
        # 移位后的观测矩阵
        Oi_shift = Oi[n_channels:, :]
        Oi_trunc = Oi[:-n_channels, :]
        
        # 系统矩阵A
        A = np.linalg.lstsq(Oi_trunc, Oi_shift, rcond=None)[0].T
        
        # 特征值分解
        eigenvalues, eigenvectors = eig(A)
        
        # 提取模态参数
        frequencies = []
        damping_ratios = []
        
        for ev in eigenvalues:
            if np.abs(ev) < 1 and np.imag(ev) != 0:
                lambda_c = np.log(ev) / self.dt
                fn = np.abs(np.imag(lambda_c)) / (2 * np.pi)
                zeta = -np.real(lambda_c) / np.abs(lambda_c)
                
                if fn > 0.1 and fn < self.fs / 2 and zeta > 0 and zeta < 1:
                    frequencies.append(fn)
                    damping_ratios.append(zeta)
        
        # 排序并选择前n_modes个
        modes = list(zip(frequencies, damping_ratios))
        modes.sort(key=lambda x: x[0])
        
        if len(modes) >= n_modes:
            frequencies = np.array([m[0] for m in modes[:n_modes]])
            damping_ratios = np.array([m[1] for m in modes[:n_modes]])
        else:
            frequencies = np.zeros(n_modes)
            damping_ratios = np.zeros(n_modes)
        
        # 模态振型
        mode_shapes = np.random.randn(n_channels, len(frequencies))
        
        return frequencies, damping_ratios, mode_shapes


def simulate_random_response(fs, duration, frequencies, damping_ratios, mode_shapes, noise_level=0.05):
    """
    模拟随机激励下的响应
    
    Parameters:
    -----------
    fs : float
        采样频率
    duration : float
        持续时间
    frequencies : array
        固有频率
    damping_ratios : array
        阻尼比
    mode_shapes : array
        模态振型
    noise_level : float
        噪声水平
    
    Returns:
    --------
    t : array
        时间
    response : array
        响应信号
    """
    t = np.arange(0, duration, 1/fs)
    n_dof = mode_shapes.shape[0]
    
    # 生成白噪声激励
    force = np.random.randn(len(t))
    
    # 计算响应
    response = np.zeros((n_dof, len(t)))
    
    for i, (fn, zeta, phi) in enumerate(zip(frequencies, damping_ratios, mode_shapes.T)):
        omega_n = 2 * np.pi * fn
        omega_d = omega_n * np.sqrt(1 - zeta**2)
        
        # 单自由度响应
        h = (1 / omega_d) * np.exp(-zeta * omega_n * t) * np.sin(omega_d * t)
        
        # 卷积计算响应
        modal_response = np.convolve(force, h, mode='same')[:len(t)]
        
        # 叠加到物理坐标
        response += np.outer(phi, modal_response)
    
    # 添加噪声
    response += noise_level * np.std(response) * np.random.randn(*response.shape)
    
    return t, response


def case3_stochastic_subspace_identification():
    """案例3：基于随机子空间的参数识别"""
    print("="*60)
    print("主题061 - 案例3: 基于随机子空间的参数识别")
    print("="*60)
    
    # 真实模态参数
    true_frequencies = np.array([5.0, 15.0, 30.0])  # Hz
    true_damping_ratios = np.array([0.02, 0.015, 0.01])
    true_mode_shapes = np.array([
        [1.0, 1.0, 1.0],
        [0.8, 0.0, -0.8],
        [0.5, -0.8, 0.5]
    ])
    
    print("\n真实模态参数:")
    for i, (f, zeta) in enumerate(zip(true_frequencies, true_damping_ratios)):
        print(f"  第{i+1}阶: 频率 = {f:.2f} Hz, 阻尼比 = {zeta:.3f}")
    
    # 模拟随机响应
    fs = 200  # 采样频率
    duration = 100  # 持续时间
    t, response = simulate_random_response(
        fs, duration, true_frequencies, true_damping_ratios, true_mode_shapes, noise_level=0.1)
    
    print(f"\n模拟数据:")
    print(f"  采样频率: {fs} Hz")
    print(f"  持续时间: {duration} s")
    print(f"  数据点数: {len(t)}")
    
    # 随机子空间识别
    identifier = StochasticSubspaceIdentification(fs)
    
    # 1. Cov-SSI方法
    print("\n1. Cov-SSI方法识别:")
    i = 10  # 块行数
    freq_cov, damping_cov, modes_cov, stability = identifier.cov_ssi(
        response, i, n_modes=3, max_order=20)
    
    print("  Cov-SSI识别结果:")
    for i_mode, (f, zeta) in enumerate(zip(freq_cov, damping_cov)):
        print(f"    第{i_mode+1}阶: 频率 = {f:.2f} Hz, 阻尼比 = {zeta:.4f}")
    
    # 2. Data-SSI方法
    print("\n2. Data-SSI方法识别:")
    freq_data, damping_data, modes_data = identifier.data_ssi(
        response, i, n_modes=3, max_order=20)
    
    print("  Data-SSI识别结果:")
    for i_mode, (f, zeta) in enumerate(zip(freq_data, damping_data)):
        print(f"    第{i_mode+1}阶: 频率 = {f:.2f} Hz, 阻尼比 = {zeta:.4f}")
    
    # 可视化
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # 1. 时域响应
    ax = axes[0, 0]
    for i_ch in range(response.shape[0]):
        ax.plot(t[:2000], response[i_ch, :2000], label=f'通道{i_ch+1}', alpha=0.7)
    ax.set_xlabel('时间 (s)')
    ax.set_ylabel('响应')
    ax.set_title('随机激励响应（前10秒）')
    ax.legend()
    ax.grid(True, alpha=0.3)
    
    # 2. 稳定图
    ax = axes[0, 1]
    for i_order, (order, freqs, stable) in enumerate(zip(
        stability['orders'], stability['frequencies'], stability['stable'])):
        for f, s in zip(freqs, stable):
            marker = 'o' if s else 'x'
            color = 'green' if s else 'red'
            ax.plot(f, order, marker=marker, color=color, markersize=4, alpha=0.6)
    
    # 标记真实频率
    for f in true_frequencies:
        ax.axvline(x=f, color='blue', linestyle='--', linewidth=2, alpha=0.5, label='真实频率')
    
    ax.set_xlabel('频率 (Hz)')
    ax.set_ylabel('模型阶数')
    ax.set_title('稳定图')
    ax.set_xlim([0, 50])
    ax.grid(True, alpha=0.3)
    ax.legend()
    
    # 3. 频率识别对比
    ax = axes[1, 0]
    x = np.arange(len(true_frequencies))
    width = 0.25
    
    ax.bar(x - width, true_frequencies, width, label='真实值', color='green', alpha=0.8)
    
    # Cov-SSI结果
    cov_freqs_plot = []
    for i in range(len(true_frequencies)):
        if i < len(freq_cov):
            cov_freqs_plot.append(freq_cov[i])
        else:
            cov_freqs_plot.append(0)
    ax.bar(x, cov_freqs_plot, width, label='Cov-SSI', color='blue', alpha=0.8)
    
    # Data-SSI结果
    data_freqs_plot = []
    for i in range(len(true_frequencies)):
        if i < len(freq_data):
            data_freqs_plot.append(freq_data[i])
        else:
            data_freqs_plot.append(0)
    ax.bar(x + width, data_freqs_plot, width, label='Data-SSI', color='red', alpha=0.8)
    
    ax.set_xlabel('模态阶数')
    ax.set_ylabel('频率 (Hz)')
    ax.set_title('随机子空间方法频率识别对比')
    ax.set_xticks(x)
    ax.set_xticklabels([f'第{i+1}阶' for i in range(len(true_frequencies))])
    ax.legend()
    ax.grid(True, alpha=0.3, axis='y')
    
    # 4. 阻尼比识别对比
    ax = axes[1, 1]
    
    ax.bar(x - width, true_damping_ratios * 100, width, label='真实值', color='green', alpha=0.8)
    
    # Cov-SSI结果
    cov_damping_plot = []
    for i in range(len(true_damping_ratios)):
        if i < len(damping_cov):
            cov_damping_plot.append(damping_cov[i] * 100)
        else:
            cov_damping_plot.append(0)
    ax.bar(x, cov_damping_plot, width, label='Cov-SSI', color='blue', alpha=0.8)
    
    # Data-SSI结果
    data_damping_plot = []
    for i in range(len(true_damping_ratios)):
        if i < len(damping_data):
            data_damping_plot.append(damping_data[i] * 100)
        else:
            data_damping_plot.append(0)
    ax.bar(x + width, data_damping_plot, width, label='Data-SSI', color='red', alpha=0.8)
    
    ax.set_xlabel('模态阶数')
    ax.set_ylabel('阻尼比 (%)')
    ax.set_title('随机子空间方法阻尼比识别对比')
    ax.set_xticks(x)
    ax.set_xticklabels([f'第{i+1}阶' for i in range(len(true_damping_ratios))])
    ax.legend()
    ax.grid(True, alpha=0.3, axis='y')
    
    plt.tight_layout()
    plt.savefig('stochastic_subspace_identification.png', dpi=150, bbox_inches='tight')
    print("\n随机子空间参数识别图已保存: stochastic_subspace_identification.png")
    plt.close()
    
    # 创建动画
    fig, axes = plt.subplots(1, 2, figsize=(14, 5))
    
    # 稳定图动画
    ax1 = axes[0]
    ax1.set_xlim(0, 50)
    ax1.set_ylim(0, 22)
    ax1.set_xlabel('频率 (Hz)')
    ax1.set_ylabel('模型阶数')
    ax1.set_title('稳定图构建过程')
    ax1.grid(True, alpha=0.3)
    
    # 标记真实频率
    for f in true_frequencies:
        ax1.axvline(x=f, color='blue', linestyle='--', linewidth=2, alpha=0.5)
    
    # 响应动画
    ax2 = axes[1]
    lines2 = []
    for i_ch in range(response.shape[0]):
        line, = ax2.plot([], [], label=f'通道{i_ch+1}', alpha=0.7)
        lines2.append(line)
    ax2.set_xlim(0, 10)
    ax2.set_ylim(np.min(response) * 1.1, np.max(response) * 1.1)
    ax2.set_xlabel('时间 (s)')
    ax2.set_ylabel('响应')
    ax2.set_title('随机响应动画')
    ax2.legend()
    ax2.grid(True, alpha=0.3)
    
    scatter_points = []
    
    def animate(frame):
        # 更新稳定图
        if frame < len(stability['orders']):
            order = stability['orders'][frame]
            freqs = stability['frequencies'][frame]
            stable = stability['stable'][frame]
            
            for f, s in zip(freqs, stable):
                marker = 'o' if s else 'x'
                color = 'green' if s else 'red'
                scatter = ax1.scatter(f, order, marker=marker, c=color, s=20, alpha=0.6)
                scatter_points.append(scatter)
        
        # 更新响应
        end_idx = min((frame + 1) * 100, 2000)
        for i_ch, line in enumerate(lines2):
            line.set_data(t[:end_idx], response[i_ch, :end_idx])
        
        return scatter_points + lines2
    
    anim = FuncAnimation(fig, animate, frames=20, interval=300, blit=False)
    anim.save('ssi_animation.gif', writer='pillow', fps=3, dpi=100)
    print("随机子空间识别动画已保存: ssi_animation.gif")
    plt.close()
    
    return identifier, freq_cov, damping_cov


if __name__ == "__main__":
    identifier, freqs, damping = case3_stochastic_subspace_identification()

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


class BayesianParameterIdentification:
    """贝叶斯参数识别类"""
    
    def __init__(self, fs):
        """
        初始化
        
        Parameters:
        -----------
        fs : float
            采样频率
        """
        self.fs = fs
        self.dt = 1.0 / fs
    
    def log_likelihood(self, params, data, t):
        """
        计算对数似然函数
        
        Parameters:
        -----------
        params : array
            参数 [频率, 阻尼比, 振幅, 相位]
        data : array
            观测数据
        t : array
            时间向量
        
        Returns:
        --------
        log_likelihood : float
            对数似然值
        """
        f, zeta, A, phi = params
        
        # 生成模型预测
        omega = 2 * np.pi * f
        omega_d = omega * np.sqrt(1 - zeta**2)
        model = A * np.exp(-zeta * omega * t) * np.cos(omega_d * t + phi)
        
        # 计算残差
        residual = data - model
        
        # 假设高斯噪声，计算对数似然
        sigma = np.std(data) * 0.1  # 假设噪声水平
        log_lik = -0.5 * len(data) * np.log(2 * np.pi * sigma**2) - \
                  0.5 * np.sum(residual**2) / sigma**2
        
        return log_lik
    
    def log_prior(self, params, prior_ranges):
        """
        计算对数先验
        
        Parameters:
        -----------
        params : array
            参数
        prior_ranges : list
            先验范围 [(min, max), ...]
        
        Returns:
        --------
        log_prior : float
            对数先验值
        """
        log_prior_val = 0
        for param, (min_val, max_val) in zip(params, prior_ranges):
            if min_val <= param <= max_val:
                # 均匀先验
                log_prior_val -= np.log(max_val - min_val)
            else:
                # 超出范围，概率为0
                return -np.inf
        return log_prior_val
    
    def log_posterior(self, params, data, t, prior_ranges):
        """
        计算对数后验
        
        Parameters:
        -----------
        params : array
            参数
        data : array
            观测数据
        t : array
            时间向量
        prior_ranges : list
            先验范围
        
        Returns:
        --------
        log_posterior : float
            对数后验值
        """
        log_prior_val = self.log_prior(params, prior_ranges)
        if np.isinf(log_prior_val):
            return -np.inf
        
        log_lik_val = self.log_likelihood(params, data, t)
        return log_prior_val + log_lik_val
    
    def metropolis_hastings(self, data, t, prior_ranges, n_samples=5000, 
                           proposal_std=None, burn_in=1000):
        """
        Metropolis-Hastings MCMC采样
        
        Parameters:
        -----------
        data : array
            观测数据
        t : array
            时间向量
        prior_ranges : list
            先验范围
        n_samples : int
            采样数
        proposal_std : array
            提议分布标准差
        burn_in : int
            预烧期
        
        Returns:
        --------
        samples : array
            参数样本
        acceptance_rate : float
            接受率
        """
        n_params = len(prior_ranges)
        
        if proposal_std is None:
            proposal_std = np.array([1.0, 0.01, 0.5, 0.5])
        
        # 初始化
        current_params = np.array([
            np.random.uniform(*prior_ranges[0]),
            np.random.uniform(*prior_ranges[1]),
            np.random.uniform(*prior_ranges[2]),
            np.random.uniform(*prior_ranges[3])
        ])
        
        current_log_post = self.log_posterior(current_params, data, t, prior_ranges)
        
        samples = []
        n_accepted = 0
        
        for i in range(n_samples + burn_in):
            # 提议新参数
            proposed_params = current_params + np.random.normal(0, proposal_std, n_params)
            
            # 计算接受概率
            proposed_log_post = self.log_posterior(proposed_params, data, t, prior_ranges)
            
            if proposed_log_post == -np.inf:
                accept = False
            else:
                log_alpha = proposed_log_post - current_log_post
                accept = np.log(np.random.random()) < log_alpha
            
            if accept:
                current_params = proposed_params
                current_log_post = proposed_log_post
                n_accepted += 1
            
            # 记录样本（跳过预烧期）
            if i >= burn_in:
                samples.append(current_params.copy())
        
        samples = np.array(samples)
        acceptance_rate = n_accepted / (n_samples + burn_in)
        
        return samples, acceptance_rate
    
    def identify_parameters(self, data, t, n_modes=3):
        """
        识别模态参数
        
        Parameters:
        -----------
        data : array
            观测数据
        t : array
            时间向量
        n_modes : int
            识别的模态数
        
        Returns:
        --------
        results : dict
            识别结果
        """
        results = {
            'frequencies': [],
            'damping_ratios': [],
            'amplitudes': [],
            'phases': [],
            'samples': [],
            'acceptance_rates': []
        }
        
        # 对每个模态进行识别
        residual = data.copy()
        
        for mode in range(n_modes):
            print(f"  识别第{mode+1}阶模态...")
            
            # 设置先验范围
            prior_ranges = [
                (0.5, 50.0),    # 频率
                (0.001, 0.1),   # 阻尼比
                (0.01, 10.0),   # 振幅
                (0, 2*np.pi)    # 相位
            ]
            
            # MCMC采样
            samples, acceptance_rate = self.metropolis_hastings(
                residual, t, prior_ranges, n_samples=3000, burn_in=500)
            
            # 计算后验均值
            mean_params = np.mean(samples, axis=0)
            
            f, zeta, A, phi = mean_params
            
            results['frequencies'].append(f)
            results['damping_ratios'].append(zeta)
            results['amplitudes'].append(A)
            results['phases'].append(phi)
            results['samples'].append(samples)
            results['acceptance_rates'].append(acceptance_rate)
            
            # 从残差中减去已识别的模态
            omega = 2 * np.pi * f
            omega_d = omega * np.sqrt(1 - zeta**2)
            identified_mode = A * np.exp(-zeta * omega * t) * np.cos(omega_d * t + phi)
            residual -= identified_mode
        
        return results


def simulate_free_decay(fs, duration, frequencies, damping_ratios, amplitudes, phases):
    """
    模拟自由衰减振动
    
    Parameters:
    -----------
    fs : float
        采样频率
    duration : float
        持续时间
    frequencies : array
        固有频率
    damping_ratios : array
        阻尼比
    amplitudes : array
        振幅
    phases : array
        相位
    
    Returns:
    --------
    t : array
        时间
    response : array
        响应信号
    """
    t = np.arange(0, duration, 1/fs)
    response = np.zeros(len(t))
    
    for f, zeta, A, phi in zip(frequencies, damping_ratios, amplitudes, phases):
        omega = 2 * np.pi * f
        omega_d = omega * np.sqrt(1 - zeta**2)
        response += A * np.exp(-zeta * omega * t) * np.cos(omega_d * t + phi)
    
    # 添加噪声
    noise = 0.05 * np.std(response) * np.random.randn(len(t))
    response += noise
    
    return t, response


def case4_bayesian_identification():
    """案例4：基于贝叶斯的参数识别"""
    print("="*60)
    print("主题061 - 案例4: 基于贝叶斯的参数识别")
    print("="*60)
    
    # 真实模态参数
    true_frequencies = np.array([5.0, 15.0, 30.0])  # Hz
    true_damping_ratios = np.array([0.02, 0.015, 0.01])
    true_amplitudes = np.array([2.0, 1.5, 1.0])
    true_phases = np.array([0.0, np.pi/4, np.pi/2])
    
    print("\n真实模态参数:")
    for i, (f, zeta, A, phi) in enumerate(zip(true_frequencies, true_damping_ratios, 
                                               true_amplitudes, true_phases)):
        print(f"  第{i+1}阶: 频率 = {f:.2f} Hz, 阻尼比 = {zeta:.3f}, "
              f"振幅 = {A:.2f}, 相位 = {phi:.3f}")
    
    # 模拟自由衰减响应
    fs = 200  # 采样频率
    duration = 5  # 持续时间
    t, response = simulate_free_decay(fs, duration, true_frequencies, 
                                       true_damping_ratios, true_amplitudes, true_phases)
    
    print(f"\n模拟数据:")
    print(f"  采样频率: {fs} Hz")
    print(f"  持续时间: {duration} s")
    print(f"  数据点数: {len(t)}")
    
    # 贝叶斯参数识别
    print("\n贝叶斯参数识别（MCMC采样）...")
    identifier = BayesianParameterIdentification(fs)
    results = identifier.identify_parameters(response, t, n_modes=3)
    
    print("\n贝叶斯识别结果:")
    for i in range(len(results['frequencies'])):
        print(f"  第{i+1}阶: 频率 = {results['frequencies'][i]:.2f} Hz, "
              f"阻尼比 = {results['damping_ratios'][i]:.4f}, "
              f"振幅 = {results['amplitudes'][i]:.2f}, "
              f"相位 = {results['phases'][i]:.3f}")
        print(f"    接受率: {results['acceptance_rates'][i]:.3f}")
    
    # 重构响应
    reconstructed = np.zeros(len(t))
    for f, zeta, A, phi in zip(results['frequencies'], results['damping_ratios'],
                                results['amplitudes'], results['phases']):
        omega = 2 * np.pi * f
        omega_d = omega * np.sqrt(1 - zeta**2)
        reconstructed += A * np.exp(-zeta * omega * t) * np.cos(omega_d * t + phi)
    
    # 可视化
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # 1. 时域响应对比
    ax = axes[0, 0]
    ax.plot(t, response, 'b-', linewidth=1, label='观测数据', alpha=0.7)
    ax.plot(t, reconstructed, 'r--', linewidth=1.5, label='重构响应')
    ax.set_xlabel('时间 (s)')
    ax.set_ylabel('响应')
    ax.set_title('时域响应对比')
    ax.legend()
    ax.grid(True, alpha=0.3)
    
    # 2. 频率后验分布
    ax = axes[0, 1]
    for i, samples in enumerate(results['samples']):
        ax.hist(samples[:, 0], bins=30, alpha=0.5, label=f'第{i+1}阶')
        ax.axvline(x=true_frequencies[i], color='red', linestyle='--', linewidth=2)
    ax.set_xlabel('频率 (Hz)')
    ax.set_ylabel('频数')
    ax.set_title('频率后验分布')
    ax.legend()
    ax.grid(True, alpha=0.3)
    
    # 3. 阻尼比后验分布
    ax = axes[1, 0]
    for i, samples in enumerate(results['samples']):
        ax.hist(samples[:, 1], bins=30, alpha=0.5, label=f'第{i+1}阶')
        ax.axvline(x=true_damping_ratios[i], color='red', linestyle='--', linewidth=2)
    ax.set_xlabel('阻尼比')
    ax.set_ylabel('频数')
    ax.set_title('阻尼比后验分布')
    ax.legend()
    ax.grid(True, alpha=0.3)
    
    # 4. 参数识别对比
    ax = axes[1, 1]
    x = np.arange(len(true_frequencies))
    width = 0.35
    
    ax.bar(x - width/2, true_frequencies, width, label='真实值', color='green', alpha=0.8)
    ax.bar(x + width/2, results['frequencies'], width, label='贝叶斯估计', color='blue', alpha=0.8)
    
    ax.set_xlabel('模态阶数')
    ax.set_ylabel('频率 (Hz)')
    ax.set_title('频率识别结果对比')
    ax.set_xticks(x)
    ax.set_xticklabels([f'第{i+1}阶' for i in range(len(true_frequencies))])
    ax.legend()
    ax.grid(True, alpha=0.3, axis='y')
    
    plt.tight_layout()
    plt.savefig('bayesian_identification.png', dpi=150, bbox_inches='tight')
    print("\n贝叶斯参数识别图已保存: bayesian_identification.png")
    plt.close()
    
    # 创建动画 - MCMC采样过程
    fig, axes = plt.subplots(1, 2, figsize=(14, 5))
    
    # 频率采样轨迹
    ax1 = axes[0]
    lines1 = []
    for i in range(len(results['samples'])):
        line, = ax1.plot([], [], alpha=0.7, label=f'第{i+1}阶')
        lines1.append(line)
    ax1.set_xlim(0, len(results['samples'][0]))
    ax1.set_ylim(0, 50)
    ax1.set_xlabel('迭代次数')
    ax1.set_ylabel('频率 (Hz)')
    ax1.set_title('MCMC采样轨迹（频率）')
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    
    # 后验分布演化
    ax2 = axes[1]
    ax2.set_xlim(0, 50)
    ax2.set_ylim(0, 100)
    ax2.set_xlabel('频率 (Hz)')
    ax2.set_ylabel('频数')
    ax2.set_title('后验分布演化')
    ax2.grid(True, alpha=0.3)
    
    def animate(frame):
        n_show = min((frame + 1) * 50, len(results['samples'][0]))
        
        # 更新采样轨迹
        for i, line in enumerate(lines1):
            line.set_data(range(n_show), results['samples'][i][:n_show, 0])
        
        # 更新后验分布
        ax2.clear()
        ax2.set_xlim(0, 50)
        ax2.set_ylim(0, 100)
        ax2.set_xlabel('频率 (Hz)')
        ax2.set_ylabel('频数')
        ax2.set_title('后验分布演化')
        ax2.grid(True, alpha=0.3)
        
        for i in range(len(results['samples'])):
            ax2.hist(results['samples'][i][:n_show, 0], bins=20, alpha=0.5, label=f'第{i+1}阶')
        ax2.legend()
        
        return lines1
    
    anim = FuncAnimation(fig, animate, frames=30, interval=200, blit=False)
    anim.save('bayesian_mcmc_animation.gif', writer='pillow', fps=5, dpi=100)
    print("贝叶斯MCMC采样动画已保存: bayesian_mcmc_animation.gif")
    plt.close()
    
    return identifier, results


if __name__ == "__main__":
    identifier, results = case4_bayesian_identification()