### 基于 Python 的模态参数识别程序设计（适用于环境激励） 以下是针对仅包含环境激励输出数据的情况，基于 Python 编写的模态参数识别程序。此方案主要采用自然激励技术（Natural Excitation Technique, NExT）结合奇异值分解（SVD）和特征系统实现算法（ERA）完成模态参数的提取。 #### 1. 导入库文件 为了支持数据分析与处理功能，需加载必要的科学计算库： ```python import numpy as np from scipy.linalg import hankel from scipy.signal import correlate import matplotlib.pyplot as plt ``` #### 2. 准备环境激励下的响应数据 假定已获得一组由多个传感器采集的时间序列响应数据 `responses`，这些数据反映了结构在环境激励作用下的动态行为。 ```python # 模拟环境激励下的一组响应数据 def simulate_responses(num_sensors=3, duration=10, sampling_rate=500): t = np.linspace(0, duration, int(sampling_rate * duration), endpoint=False) responses = [] for i in range(num_sensors): freq = 5 + i * 2 # 不同传感器可能捕捉到不同频率成分 noise_level = 0.1 * (i + 1) # 随机噪声水平逐渐增加 signal = np.sin(2 * np.pi * freq * t) + noise_level * np.random.randn(len(t)) responses.append(signal) return t, np.array(responses).T time_series, responses = simulate_responses() ``` #### 3. 计算互相关函数矩阵 通过计算各通道之间的互相关函数，构建互相关函数矩阵 $ R(\tau) $，这是 NExT 方法的关键步骤之一。 ```python def compute_cross_correlation_matrix(data): num_channels = data.shape[1] cross_corr_matrix = [] for channel_i in range(num_channels): row = [] for channel_j in range(num_channels): corr_func = correlate(data[:, channel_i], data[:, channel_j], mode='full') normalized_corr = corr_func / max(corr_func) # 归一化处理 row.append(normalized_corr[len(normalized_corr)//2:]) # 只保留正延迟部分 cross_corr_matrix.append(row) return np.array(cross_corr_matrix) cross_corr_matrix = compute_cross_correlation_matrix(responses) ``` #### 4. 构建 Hankel 矩阵并进行 SVD 分解 Hankel 矩阵能有效地捕获时间序列中的动态特性，随后对其进行奇异值分解以分离不同的振动模式。 ```python def build_hankel_and_decompose(matrix, window_size=None): if not window_size: window_size = matrix.shape[-1] // 2 # 默认窗口大小设为一半长度 hankel_matrices = [ hankel(channel[:window_size], channel[window_size:]) for channel in matrix.flatten() ] svd_results = [np.linalg.svd(hm, full_matrices=False) for hm in hankel_matrices] singular_values = np.array([sv[1] for sv in svd_results]) modes = sum((sv > 1e-6).sum() for sv in singular_values) # 统计有效奇异值数量作为估计的模态数 return modes, svd_results modes_estimated, svd_results = build_hankel_and_decompose(cross_corr_matrix) print(f"Estimated number of modes: {modes_estimated}") ``` #### 5. 应用 ERA 提取模态参数 最后一步是利用特征系统实现算法（ERA），从降维后的状态空间模型中提取固有频率、阻尼比和其他模态属性。 ```python def estimate_modal_parameters(svd_results, dt=0.002): U_all = np.concatenate([res[0][:, :int(modes_estimated)] for res in svd_results], axis=-1) Sigma_all = np.diagflat(np.concatenate([res[1][:int(modes_estimated)] for res in svd_results])) Vh_all = np.concatenate([res[2].conj().T[:, :int(modes_estimated)] for res in svd_results], axis=-1) A_id = U_all @ Sigma_all @ Vh_all.T.conj() # Identified system matrix A eigenvalues, eigenvectors = np.linalg.eig(A_id) frequencies = np.log(eigenvalues).imag / (2 * np.pi * dt) damping_ratios = -np.log(abs(eigenvalues)) / (2 * np.sqrt(-np.log(abs(eigenvalues)).real ** 2)) valid_indices = (frequencies > 0) & (~np.isnan(damping_ratios)) & (damping_ratios >= 0) modal_frequencies = frequencies[valid_indices] modal_dampings = damping_ratios[valid_indices] return modal_frequencies, modal_dampings modal_frequencies, modal_dampings = estimate_modal_parameters(svd_results, dt=1/500) ``` #### 6. 展示结果 将提取得到的模态参数可视化以便直观理解。 ```python plt.figure(figsize=(10, 6)) plt.scatter(modal_frequencies, modal_dampings, c='blue', marker='o', s=80, edgecolors='k', alpha=0.7) plt.title("Modal Frequencies vs Damping Ratios", fontsize=16) plt.xlabel("Frequency (Hz)", fontsize=14) plt.ylabel("Damping Ratio ($\zeta$)", fontsize=14) plt.grid(True, which="both", ls="--", linewidth=0.5) plt.xticks(fontsize=12) plt.yticks(fontsize=12) plt.show() ``` --- ### 总结说明 以上流程展示了一种完全基于环境激励输出数据的模态参数识别方法[^1]。它无需显式的输入信号即可工作，非常适合无法获取精确输入条件的实际工程应用场景。其中涉及的技术要点包括但不限于互相关函数计算、Hankel 矩阵构造及其奇异值分解(SVD)，以及最终借助 ERA 进行模态参数辨识的过程。