简介本资源是一套基于时域分解TDD方法提取结构模态振型的MATLAB实现代码面向机械、土木及航空航天等领域的工程技术人员、高校师生与科研初学者解决结构动态特性分析中模态参数频率、振型难以从实测时域信号中准确分离的实践难题。压缩包共4个文件含可交互运行的Example.mlx含注释与可视化、核心算法脚本TDD.m、实测梁结构振动数据beamData.mat以及说明文档README.md整体3.97MB适配MATLAB 2014–2024a多版本开箱即用。已有154人学习下载无需额外准备数据即可复现完整TDD流程——从原始时程信号输入、模态参数辨识到振型可视化输出覆盖算法原理、代码逻辑与工程解释三层内容特别适合理解模态分析底层机制及快速开展结构健康监测、振动抑制等方向的仿真实验。1. 时域分解TDD不是频域方法它用原始振动信号直接分离模态振型——适合实验模态分析中无激励力测量、传感器数量少但采样率高的场景很多做结构健康监测或机械振动分析的工程师一看到“模态振型”就默认要FFT、谱密度、峰值拾取或ERA/ITD这类频域/时域联合方法。但TDDTime Domain Decomposition完全不同它不依赖傅里叶变换也不需要已知激励力而是把一段多通道实测振动响应时间序列比如加速度计阵列采集的10秒数据通过奇异值分解SVD和时延嵌入time-delay embedding构造汉克尔矩阵再结合物理约束如模态衰减指数、频率分布先验对子空间进行旋转与解耦最终直接输出各阶模态振型向量和对应模态坐标时间历程。整个过程完全在时域完成对非平稳信号鲁棒性更强特别适合桥梁微振动、风力机塔架低频晃动、或实验室中无法布置力传感器的悬臂梁敲击试验。本篇聚焦TDD在MATLAB中的可复现实现——不调用任何第三方工具箱仅用基础矩阵运算与优化函数代码结构清晰、参数可调、结果可验证适用于R2018b及以上版本含R2023b/R2026a等主流发行版。2. 构建汉克尔矩阵与子空间识别从原始信号到模态子空间的三步核心流程TDD的理论根基在于线性时不变系统自由响应的数学表达多通道输出 $ y(t) \in \mathbb{R}^{m \times 1} $ 可近似为 $ r $ 阶模态叠加 $ y(t) \approx \sum_{i1}^r \phi_i , a_i(t) $其中 $ \phi_i \in \mathbb{R}^{m} $ 是第 $ i $ 阶模态振型待求$ a_i(t) $ 是对应模态坐标的时域响应。TDD的关键洞察是若将 $ y(t) $ 按固定时延 $ \tau $ 构造成汉克尔矩阵 $ H \in \mathbb{R}^{Lm \times (N-L1)} $其列空间将张成由所有模态振型张成的 $ r $ 维子空间。因此第一步是构造该矩阵并提取其左奇异向量作为初始子空间估计。2.1 时延嵌入与汉克尔矩阵生成控制L与τ的物理意义汉克尔矩阵的维度由两个关键参数决定嵌入长度 $ L $即行数分块数和时延步长 $ \tau $单位采样点。设原始信号为 $ Y \in \mathbb{R}^{m \times N} $其中 $ m $ 为传感器通道数$ N $ 为总采样点数。以下MATLAB代码生成标准汉克尔结构function H build_hankel_matrix(Y, L, tau) % Y: m x N 矩阵每行一个通道每列一个时间点 % L: 嵌入长度行块数建议取 L floor(N/4) ~ floor(N/2)需满足 L*tau N % tau: 时延步长采样点数通常取1连续采样或根据Nyquist准则调整 m size(Y, 1); N size(Y, 2); if L*tau N error(L*tau must be less than N: increase N or decrease L/tau); end % 计算有效列数 n_cols N - L*tau 1; H zeros(L*m, n_cols); for k 1:n_cols for l 0:L-1 col_idx k l*tau; if col_idx N H((l*m1):(l1)*m, k) Y(:, col_idx); else break; end end end end提示tau1表示最密集嵌入保留全部时序相关性tau1可降低矩阵条件数尤其当采样率远高于模态频率时如10 kHz采样测100 Hz模态tau5~10能有效抑制高频噪声干扰。L过小如5导致子空间维数不足无法分辨密集模态L过大如200则引入冗余并放大数值误差。实践中L应覆盖至少2~3个最低阶模态周期——例如最低模态频率为5 Hz、采样率为1000 Hz则单周期200点取L50~100较稳妥。2.2 SVD分解与稳定图判据如何确定真实模态阶数r对汉克尔矩阵 $ H $ 执行SVD$ H U \Sigma V^T $其中 $ U \in \mathbb{R}^{Lm \times Lm} $ 的前 $ r $ 列 $ U_r $ 即为模态子空间的初始估计。但 $ r $ 未知需通过稳定图stability diagram判断。TDD中常用两种判据奇异值衰减比计算 $ \sigma_{i}/\sigma_{i-1} $当比值突增如10表明后续奇异值主要由噪声贡献模态置信度MAC一致性对不同L值重复SVD计算同一阶次左右奇异向量的模态保证准则MAC值高MAC值0.95对应稳定模态。以下代码实现双判据联合判定function [r_est, sig_ratio, mac_table] estimate_modal_order(H, L_list, Y, max_r) % H: 当前L下的汉克尔矩阵 % L_list: 测试的L值数组如[20,40,60,80] % Y: 原始信号用于MAC计算 % max_r: 最大搜索阶数建议取 min(10, size(H,2)/2) U_all {}; for idx 1:length(L_list) H_test build_hankel_matrix(Y, L_list(idx), 1); [~, ~, U_test] svd(H_test, econ); U_all{idx} U_test(:, 1:max_r); end % 计算奇异值衰减比 [~, S, ~] svd(H, econ); sig_vals diag(S); sig_ratio ones(size(sig_vals)); for i 2:length(sig_vals) sig_ratio(i) sig_vals(i)/sig_vals(i-1); end % 计算MAC表U_all{i}(:,j) 与 U_all{k}(:,j) 的MAC mac_table zeros(max_r, length(L_list), length(L_list)); for j 1:max_r for i 1:length(L_list) for k i:length(L_list) u_i U_all{i}(:,j); u_k U_all{k}(:,j); mac_val abs(u_i * u_k)^2 / ( (u_i*u_i) * (u_k*u_k) ); mac_table(j,i,k) mac_val; mac_table(j,k,i) mac_val; end end end % 综合判定取奇异值比5 且 平均MAC0.9 的最大j r_est 1; for j 1:max_r if sig_ratio(j1) 5 mean(mac_table(j,:,:)(mac_table(j,:,:) 0)) 0.9 r_est j; else break; end end end注意max_r不宜过大否则MAC计算量剧增。实际工程中前6阶模态已覆盖绝大多数结构动力学问题。若sig_ratio在第3阶后持续2而mac_table中第4阶平均MAC仅0.7则说明第4阶可能是噪声模态或测量误差主导应截断至r3。3. 模态振型解析与物理约束嵌入从数学子空间到可解释振型向量获得子空间 $ U_r $ 后TDD的核心挑战是如何将其映射为物理意义明确的模态振型 $ \Phi [\phi_1,\dots,\phi_r] $。纯数学SVD给出的是正交基但真实振型需满足① 各阶振型间正交质量/刚度正交② 振型幅值具有相对比例关系如某传感器响应最大对应节点位移最大③ 模态坐标 $ a_i(t) $ 应呈衰减正弦形式。因此需引入物理约束进行子空间旋转。3.1 基于模态坐标时域拟合的振型缩放TDD标准做法是假设模态坐标 $ a_i(t) $ 可表示为 $ a_i(t) e^{-\zeta_i \omega_i t} \cos(\omega_i t \theta_i) $其中 $ \zeta_i $ 为阻尼比$ \omega_i $ 为固有频率。对 $ U_r $ 的每一列 $ u_i $将其与原始信号 $ Y $ 进行最小二乘投影得到初始模态坐标估计 $ \hat{a}_i(t) $再对该时间序列进行非线性拟合反推 $ \zeta_i $ 和 $ \omega_i $最后用拟合残差修正振型缩放因子。function [Phi, A, zeta, omega] extract_mode_shapes(Ur, Y, fs) % Ur: m*L x r 子空间矩阵来自SVD % Y: m x N 原始信号 % fs: 采样频率Hz m size(Y, 1); N size(Y, 2); r size(Ur, 2); % 步骤1投影得到初始模态坐标 A0 ∈ r x N A0 Ur * Y(:); % 展开为向量后投影 A0 reshape(A0, r, N); % 恢复为 r x N % 步骤2对每阶模态坐标进行衰减正弦拟合 Phi zeros(m, r); A zeros(r, N); zeta zeros(r, 1); omega zeros(r, 1); t (0:N-1)/fs; for i 1:r ai A0(i, :); % 初始猜测FFT找主频极值点估算衰减 f_fft (0:N/2)*fs/N; Yf fft(ai); [~, idx_max] max(abs(Yf(1:floor(N/2)1))); omega0 f_fft(idx_max) * 2*pi; % rad/s % 拟合模型ai(t) exp(-zeta*omega*t) * cos(omega*t theta) opts optimoptions(lsqcurvefit,Display,off,MaxFunctionEvaluations,1000); lb [0, 0.1*omega0, -pi]; ub [0.1, 10*omega0, pi]; x0 [0.01, omega0, 0]; try x_fit lsqcurvefit(damped_cosine_model, x0, t, ai, lb, ub, opts); zeta(i) x_fit(1); omega(i) x_fit(2); % 步骤3用拟合后的ai_ref重新计算振型最小二乘 ai_ref damped_cosine_model(x_fit, t); % 解 A_ref * phi_i Y_i phi_i (A_ref^T A_ref)^{-1} A_ref^T Y_i A_ref repmat(ai_ref, m, 1); % m x N Y_i Y; Phi(:,i) (A_ref * A_ref) \ (A_ref * Y_i(:)); A(i,:) ai_ref; catch % 拟合失败时回退到SVD第一列归一化 Phi(:,i) Ur(1:m,i) / norm(Ur(1:m,i)); A(i,:) ai; zeta(i) NaN; omega(i) NaN; end end end function y_fit damped_cosine_model(x, t) % x [zeta, omega, theta] y_fit exp(-x(1)*x(2)*t) .* cos(x(2)*t x(3)); end逻辑说明Ur的前m行对应第一个时延块物理上最接近原始传感器输出因此Ur(1:m,i)可视为第i阶振型的粗略估计。但直接使用会导致振型幅值无物理意义因SVD缩放任意。本方法通过拟合模态坐标时域行为将振型缩放与系统物理参数阻尼、频率绑定使Phi(:,i)的元素代表各传感器相对于参考点的相对位移幅值。例如若Phi(3,i)2.1、Phi(7,i)0.8则说明第i阶模态下3号传感器振幅约为7号的2.6倍。3.2 振型正交性校验与MAC矩阵输出提取后的振型需验证其物理合理性。TDD要求振型满足质量正交性$ \Phi^T M \Phi I $但实际中常以模态保证准则MAC衡量振型独立性$$ \text{MAC}(\phi_i, \phi_j) \frac{|\phi_i^T \phi_j|^2}{(\phi_i^T \phi_i)(\phi_j^T \phi_j)} $$MAC≈1表示两阶振型高度相关可能为虚假模态MAC0.1表示正交性良好。function mac_matrix compute_mac(Phi) % Phi: m x r 振型矩阵 r size(Phi, 2); mac_matrix zeros(r, r); for i 1:r for j 1:r num abs(Phi(:,i) * Phi(:,j))^2; den (Phi(:,i) * Phi(:,i)) * (Phi(:,j) * Phi(:,j)); mac_matrix(i,j) num / den; end end % 输出上三角部分避免重复 fprintf(MAC matrix (upper triangle):\n); disp(triu(mac_matrix, 1)); end参数说明triu(mac_matrix, 1)仅显示ij的MAC值因MAC(i,j)MAC(j,i)且对角线恒为1。若mac_matrix(2,3)0.92说明第2、3阶振型高度耦合需检查是否为密集模态未分离或考虑增加L值重跑。4. 完整TDD流程封装与典型参数配置表一键运行可复现的MATLAB脚本将前述模块整合为可直接调用的主函数输入为多通道时间序列输出为振型矩阵、模态频率、阻尼比及验证指标。以下为完整封装代码包含默认参数推荐与错误处理。function [Phi, freq_hz, zeta, mac_mat, info] tdd_modal_analysis(Y, fs, varargin) % TDD Modal Analysis: Extract mode shapes from time-domain response only % Input: % Y: m x N matrix, each row is a sensor channel % fs: sampling frequency (Hz) % varargin: optional name-value pairs: % L - Hankel embedding length (default: floor(N/3)) % tau - time delay step (default: 1) % max_r - max modal order to search (default: 8) % L_list - L values for stability diagram (default: [20,40,60]) % Output: % Phi: m x r mode shape matrix (columns are mode shapes) % freq_hz: r x 1 vector of natural frequencies (Hz) % zeta: r x 1 vector of damping ratios % mac_mat: r x r MAC matrix % info: struct with intermediate matrices and diagnostics p inputParser; addParameter(p, L, floor(size(Y,2)/3)); addParameter(p, tau, 1); addParameter(p, max_r, 8); addParameter(p, L_list, [20,40,60]); parse(p, varargin{:}); L p.Results.L; tau p.Results.tau; max_r p.Results.max_r; L_list p.Results.L_list; % Step 1: Build Hankel matrix H build_hankel_matrix(Y, L, tau); % Step 2: Estimate modal order r [r_est, ~, ~] estimate_modal_order(H, L_list, Y, max_r); if r_est 0, r_est 1; end % Step 3: SVD to get subspace [~, ~, U] svd(H, econ); Ur U(:, 1:r_est); % Step 4: Extract mode shapes with physical constraints [Phi, A, zeta_vec, omega_vec] extract_mode_shapes(Ur, Y, fs); % Step 5: Compute frequencies and MAC freq_hz omega_vec / (2*pi); mac_mat compute_mac(Phi); % Package info info struct(... Hankel_matrix, H, ... subspace_Ur, Ur, ... modal_coordinates, A, ... estimated_order, r_est, ... singular_values, diag(svd(H, econ))(1:min(20, size(H,2))) ... ); % Normalize each mode shape to unit max amplitude for plotting for i 1:size(Phi,2) Phi(:,i) Phi(:,i) / max(abs(Phi(:,i))); end end4.1 典型工况参数配置表针对不同结构类型快速选参结构类型采样率 (Hz)推荐L值推荐taumax_r关键注意事项小型金属悬臂梁500040–8014–6高频模态密集L取中值防过拟合混凝土桥梁桥面200100–2002–53–5低频主导10 Hztau3抑制交通噪声风力机塔架100150–30012–4强非平稳性优先用tau1保时序航空发动机叶片50000200–50010–206–10超高频模态tau15避免混叠提示表中L与tau需协同调整。例如桥梁工况若fs200 Hz最低模态约1.5 Hz周期667 ms ≈ 133点取L150覆盖2个周期tau3则实际时间跨度为150×3/2002.25 s足够捕获衰减过程。若tau过大如tau10则L150对应7.5 s可能混入环境变化干扰。4.2 验证案例用仿真信号测试TDD代码可靠性构造一个双自由度系统2-DOF的自由响应作为黄金标准验证代码输出是否匹配理论振型% 生成理论2-DOF响应M[1,0;0,1], K[200,-100;-100,150], C0.02*K fs 1000; T 10; N fs*T; t (0:N-1)/fs; % 理论模态phi1[0.707;0.707], phi2[-0.707;0.707], f12.15Hz, f25.42Hz y1 0.707*exp(-0.02*2*pi*2.15*t).*cos(2*pi*2.15*t) ... (-0.707)*exp(-0.02*2*pi*5.42*t).*cos(2*pi*5.42*t 0.3); y2 0.707*exp(-0.02*2*pi*2.15*t).*cos(2*pi*2.15*t) ... 0.707*exp(-0.02*2*pi*5.42*t).*cos(2*pi*5.42*t 0.3); Y_sim [y1; y2]; % 2 x N % 运行TDD [Phi_est, freq_est, zeta_est, mac_est, ~] tdd_modal_analysis(Y_sim, fs, max_r, 3); % 对比理论振型取符号一致 Phi_true [0.707, -0.707; 0.707, 0.707]; err_norm norm(Phi_est - Phi_true, fro) / norm(Phi_true, fro); fprintf(Reconstruction error (Frobenius norm): %.3f\n, err_norm); % 若 err_norm 0.05说明代码在理想条件下可靠参数说明err_norm是重建误差的Frobenius范数相对值。在无噪声理想信号下err_norm0.02属优秀加入5%白噪声后err_norm0.08仍属可用。若误差0.15需检查L是否过小或max_r是否误设。5. 振型可视化与工程解读技巧如何从Φ矩阵读出结构动态特性TDD输出的Phi矩阵本身是数学对象必须结合传感器物理布局才能转化为工程洞见。核心技巧在于振型符号不代表方向只反映相对相位振型幅值比决定节点位置多阶振型叠加揭示复杂变形模式。5.1 传感器布局映射与振型图绘制假设4个加速度计沿简支梁等距布置位置0.2L, 0.4L, 0.6L, 0.8LPhi(:,1)为第一阶振型。以下代码生成标准振型图function plot_mode_shape(Phi, positions, mode_idx, title_str) % Phi: m x r, positions: 1 x m vector of sensor locations % mode_idx: which mode to plot (1-based) m size(Phi, 1); if length(positions) ~ m error(positions length must equal number of sensors); end figure; plot(positions, Phi(:,mode_idx), -o, LineWidth, 1.5, MarkerSize, 8); xlabel(Position (m)); ylabel(Relative Amplitude); title([title_str, - Mode , num2str(mode_idx)]); grid on; % 添加零线与节点标注 yline(0, --k, Zero line); % 查找过零点节点 zero_crossings []; for i 1:m-1 if Phi(i,mode_idx)*Phi(i1,mode_idx) 0 % 线性插值找零点 x_zero positions(i) (0-Phi(i,mode_idx)) * (positions(i1)-positions(i)) / (Phi(i1,mode_idx)-Phi(i,mode_idx)); zero_crossings [zero_crossings, x_zero]; end end if ~isempty(zero_crossings) text(zero_crossings(1), 0.1*max(abs(Phi(:,mode_idx))), Node, Color,r,FontSize,10); end end % 调用示例 positions [0.2, 0.4, 0.6, 0.8]; % 单位米 plot_mode_shape(Phi, positions, 1, Cantilever Beam);逻辑说明yline(0,--k)绘制零线直观显示节点node位置zero_crossings计算传感器间过零点即实际节点所在区间。若Phi(:,1)[0.1, 0.5, -0.4, -0.2]则节点在0.4–0.6 m之间符合一阶弯曲模态特征。5.2 多阶振型能量占比分析识别主导模态结构响应常由少数几阶模态主导。计算各阶模态振型的能量贡献比$$ \text{Energy}_i \frac{|\phi_i|2^2}{\sum{j1}^r |\phi_j|_2^2} $$function energy_ratio compute_mode_energy(Phi) % Phi: m x r r size(Phi, 2); norm_sq zeros(r, 1); for i 1:r norm_sq(i) norm(Phi(:,i))^2; end energy_ratio norm_sq / sum(norm_sq); fprintf(Mode energy ratio:\n); for i 1:r fprintf(Mode %d: %.1f%%\n, i, 100*energy_ratio(i)); end end工程解读若Mode 1: 65.2%,Mode 2: 22.1%,Mode 3: 8.7%说明结构动力学行为主要由前两阶模态决定后续模态可忽略。若Mode 4: 15.3%且freq_hz(4)接近激励源频率则提示可能存在共振风险需在设计中规避。5.3 振型置信度量化MAC与相位一致性双指标仅靠MAC不足以判断振型可靠性还需检查模态坐标相位一致性。对同一阶模态不同传感器信号经振型加权后应具有一致相位function phase_consistency check_phase_consistency(Y, Phi, mode_idx, fs) % Y: m x N, Phi: m x r, mode_idx: target mode m size(Y, 1); % 加权合成Y_weighted Phi(:,mode_idx) * Y weighted_signal Phi(:,mode_idx) * Y; % 1 x N % 计算每个传感器通道与加权信号的相位差FFT phase_diffs zeros(m, 1); for i 1:m % 互谱相位 Pxy cpsd(Y(i,:), weighted_signal, [], [], [], fs); [Pxy_f, f] cpsd(Y(i,:), weighted_signal, [], [], [], fs); phase_at_peak angle(Pxy_f(find(abs(Pxy_f)max(abs(Pxy_f)),1))); phase_diffs(i) mod(phase_at_peak pi, 2*pi) - pi; % 归到[-pi,pi] end phase_consistency std(phase_diffs); % 标准差越小相位越一致 fprintf(Phase consistency (std of phase diffs): %.3f rad\n, phase_consistency); % 0.2 rad 为优秀0.5 rad 为可接受 end参数说明phase_consistency是各传感器与模态坐标加权信号的相位差标准差。值越小说明该阶振型物理意义越强——所有传感器振动确实同步按此比例叠加。若phase_consistency0.8 rad则需怀疑该阶模态是否受局部噪声污染建议检查对应传感器安装状态。本文还有配套的精品资源点击获取