简介：本资源是一套基于时域分解（TDD）方法提取结构模态振型的MATLAB实现代码，面向机械、土木及航空航天等领域的工程技术人员、高校师生与科研初学者，解决结构动态特性分析中模态参数（频率、振型）难以从实测时域信号中准确分离的实践难题。压缩包共4个文件：含可交互运行的Example.mlx（含注释与可视化）、核心算法脚本TDD.m、实测梁结构振动数据beamData.mat，以及说明文档README.md，整体3.97MB，适配MATLAB 2014–2024a多版本，开箱即用。已有154人学习下载，无需额外准备数据即可复现完整TDD流程——从原始时程信号输入、模态参数辨识到振型可视化输出，覆盖算法原理、代码逻辑与工程解释三层内容，特别适合理解模态分析底层机制及快速开展结构健康监测、振动抑制等方向的仿真实验。
1. 时域分解（TDD）不是频域方法，它用原始振动信号直接分离模态振型——适合实验模态分析中无激励力测量、传感器数量少但采样率高的场景
很多做结构健康监测或机械振动分析的工程师一看到“模态振型”就默认要FFT、谱密度、峰值拾取或ERA/ITD这类频域/时域联合方法。但TDD（Time 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_{i=1}^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-L+1)} $，其列空间将张成由所有模态振型张成的 $ 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*m+1):(l+1)*m, k) = Y(:, col_idx); else break; end end end end提示：
tau=1表示最密集嵌入，保留全部时序相关性；tau>1可降低矩阵条件数，尤其当采样率远高于模态频率时（如10 kHz采样测100 Hz模态），tau=5~10能有效抑制高频噪声干扰。L过小（如<5）导致子空间维数不足，无法分辨密集模态；L过大（如>200）则引入冗余并放大数值误差。实践中，L应覆盖至少2~3个最低阶模态周期——例如最低模态频率为5 Hz、采样率为1000 Hz，则单周期200点，取L=50~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 且 平均MAC>0.9 的最大j r_est = 1; for j = 1:max_r if sig_ratio(j+1) < 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阶可能是噪声模态或测量误差主导，应截断至r=3。
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表示两阶振型高度相关（可能为虚假模态），MAC<0.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)仅显示i<j的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值 | 推荐tau | max_r | 关键注意事项 |
|---|---|---|---|---|---|
| 小型金属悬臂梁 | 5000 | 40–80 | 1 | 4–6 | 高频模态密集，L取中值防过拟合 |
| 混凝土桥梁桥面 | 200 | 100–200 | 2–5 | 3–5 | 低频主导（<10 Hz），tau=3抑制交通噪声 |
| 风力机塔架 | 100 | 150–300 | 1 | 2–4 | 强非平稳性，优先用tau=1保时序 |
| 航空发动机叶片 | 50000 | 200–500 | 10–20 | 6–10 | 超高频模态，tau=15避免混叠 |
提示：表中
L与tau需协同调整。例如桥梁工况若fs=200 Hz，最低模态约1.5 Hz（周期667 ms ≈ 133点），取L=150覆盖2个周期，tau=3则实际时间跨度为150×3/200=2.25 s，足够捕获衰减过程。若tau过大（如tau=10），则L=150对应7.5 s，可能混入环境变化干扰。
4.2 验证案例：用仿真信号测试TDD代码可靠性
构造一个双自由度系统（2-DOF）的自由响应作为黄金标准，验证代码输出是否匹配理论振型：
% 生成理论2-DOF响应：M=[1,0;0,1], K=[200,-100;-100,150], C=0.02*K fs = 1000; T = 10; N = fs*T; t = (0:N-1)/fs; % 理论模态：phi1=[0.707;0.707], phi2=[-0.707;0.707], f1=2.15Hz, f2=5.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_norm<0.02属优秀；加入5%白噪声后err_norm<0.08仍属可用。若误差>0.15，需检查L是否过小或max_r是否误设。
5. 振型可视化与工程解读技巧：如何从Φ矩阵读出结构动态特性
TDD输出的Phi矩阵本身是数学对象，必须结合传感器物理布局才能转化为工程洞见。核心技巧在于：振型符号不代表方向，只反映相对相位；振型幅值比决定节点位置；多阶振型叠加揭示复杂变形模式。
5.1 传感器布局映射与振型图绘制
假设4个加速度计沿简支梁等距布置（位置：0.2L, 0.4L, 0.6L, 0.8L），Phi(:,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(i+1,mode_idx) < 0 % 线性插值找零点 x_zero = positions(i) + (0-Phi(i,mode_idx)) * (positions(i+1)-positions(i)) / (Phi(i+1,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{j=1}^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_consistency=0.8 rad，则需怀疑该阶模态是否受局部噪声污染，建议检查对应传感器安装状态。
本文还有配套的精品资源，点击获取
)
