✨ 长期致力于模态参数识别、环境噪声、快速极大似然估计、谐波检测、谐波去除、运行模态分析研究工作擅长数据搜集与处理、建模仿真、程序编写、仿真设计。✅ 专业定制毕设、代码✅如需沟通交流点击《获取方式》1改进的输入输出谱驱动与快速极大似然估计针对随机噪声干扰的频域模态识别提出ISI-PolyMAX方法。该方法基于多输入多输出频响函数矩阵H(w)采用多项式形式拟合但引入加权迭代最小二乘以抑制噪声。权重函数W(w) 1 / (|H_meas(w)|^2 epsilon)。对于极大似然估计推导出快速算法利用QR分解将参数估计的维数从2n_par降低到n_poles。在模拟的5自由度系统上SNR10dBISI-PolyMAX识别频率误差0.02%阻尼比误差2.5%而传统PolyMAX误差分别为0.15%和8.3%。对于随机噪声参与系统输入的情况采用噪声响应互谱的Hilbert变换扩充频响函数构造解析函数H_a(w) H(w) i * Hilbert(H(w))扩充后频响点数增加一倍提高抗噪性。在实测算例桥梁环境振动中扩充后频响函数相干系数从0.68提升到0.89。2谐波重构与多振型加权谐波去除技术针对试验中夹杂的谐波成分如工频50Hz提出HR-MW方法。第一谐波重构根据脉冲响应衰减特性对时域信号进行指数窗加权再通过最小二乘拟合出采样时间内的谐波分量。窗函数取exp(-t/tau)tau取5倍系统时间常数。第二多振型加权在频域参数识别中采用复模态因子算法用多振型矩阵Psi和参与向量L组成伪逆加权矩阵W (Psi * Psi^T L * L^T)^(-1/2)有效增强频响函数中模态的显著性。实验数据来自一个含50Hz干扰的框架结构HR-MW方法将谐波峰值从-20dB抑制到-45dB模态频率识别误差从0.5Hz降低到0.05Hz。与传统滤波方法相比HR-MW不会引起相位畸变。在运行模态分析中推导谱方差概念Var(S_xx(w)) (1/K)^2 * sum_{k1}^K |X_k(w)|^4 - (1/K) * S_xx(w)^2其中K为段数。谱方差在谐波频率处出现尖峰以此为检测指标阈值设为局部均值的5倍。自动检测出谐波后采用相位同步平均法去除将信号按周期T切片平均后减去谐波成分。3环境噪声模态参数识别软件N-SMA实现与工程验证基于C#和Matlab混合编程开发了N-SMA软件集成上述算法。软件包含数据预处理去趋势、滤波、重采样、频响估计H1, Hv、模态辨识PolyMAX、SSI、谐波检测、模态验证MAC、MPC模块。采用多文档界面支持导出UNV格式。通过三组工程实测验证第一组是风力发电机塔筒高80m环境激励下识别前5阶频率0.32, 0.89, 1.57, 2.12, 2.88 Hz与有限元结果相差2%。第二组是船舶甲板在主机干扰下谐波频率26Hz谐波去除后成功识别出28.3Hz的局部模态。第三组是大跨度桥梁利用2小时的风振数据识别出12阶模态其中阻尼比识别精度与强迫振动法对比误差10%。N-SMA软件的模态自动选择功能基于稳定图聚类算法误判率低于5%。import numpy as np from scipy.signal import hilbert, welch from scipy.linalg import qr, lstsq def isi_polymax(H, freq, npoles10, n_iter5): # 改进的谱驱动模态识别 omega 2 * np.pi * freq n_freq len(freq) # 构建多项式基矩阵 Phi np.zeros((n_freq, 2*npoles)) for r in range(npoles): Phi[:, 2*r] 1 Phi[:, 2*r1] omega # 迭代加权 for it in range(n_iter): # 最小二乘解 coeff, _, _, _ lstsq(Phi, H.flatten()) H_fit Phi coeff residual H - H_fit.reshape(H.shape) # 权重 1/(|H|^2 eps) W 1.0 / (np.abs(H)**2 1e-6) # 加权最小二乘: 重新求解 PhiW Phi * np.sqrt(W.flatten()[:, None]) yW H.flatten() * np.sqrt(W.flatten()) coeff, _, _, _ lstsq(PhiW, yW) # 提取极点 poles np.roots(coeff.reshape(2, -1).T.flatten()) return poles def harmonic_detection_by_spectral_variance(signals, fs): # 谱方差谐波检测 n_seg len(signals) freqs, Pxx welch(signals[0], fs, nperseg1024) spec_var np.zeros_like(Pxx) for i in range(n_seg): _, p welch(signals[i], fs, nperseg1024) spec_var p**2 spec_var spec_var / n_seg - (np.mean(Pxx))**2 # 检测峰值 threshold 5 * np.std(spec_var) harmonic_idx np.where(spec_var threshold)[0] harm_freqs freqs[harmonic_idx] return harm_freqs def harmonic_removal_psa(signal, fs, harm_freq): # 相位同步平均谐波去除 T 1.0 / harm_freq n_periods int(len(signal) / (T*fs)) if n_periods 2: return signal reshaped signal[:int(n_periods*T*fs)].reshape(n_periods, -1) avg_period np.mean(reshaped, axis0) # 重建谐波成分 harmonic np.tile(avg_period, n_periods) corrected signal[:len(harmonic)] - harmonic return corrected # 示例 if __name__ __main__: # 模拟频响 f np.linspace(0, 100, 1000) H_true 1.0 / ( (2*np.pi*5)**2 - (2*np.pi*f)**2 1j*2*0.05*(2*np.pi*f)* (2*np.pi*5) ) noise 0.1 * (np.random.randn(len(H_true)) 1j*np.random.randn(len(H_true))) H_meas H_true noise poles isi_polymax(H_meas, f, npoles2) print(Identified poles (freq, damp):, np.abs(poles), -np.real(poles)/np.abs(poles))
