简介这是一份基于MATLAB的超声导波xy位移数据处理与时域分析实例面向结构健康监测、无损检测及相关科研和工程人员帮助理解导波在固体介质中的传播、衰减与反射现象并掌握从原始位移数据到时间节点位移图的完整分析流程。压缩包共5个文件包含3个mat数据文件、1个m脚本和1个xlsx表格体积仅30KB其中mat文件存储x、y方向及合成位移数据m脚本实现读取与绘图xlsx保存时间节点数据便于核对和二次处理。目前已有311人浏览学习。通过这一小体量案例可学习MATLAB处理超声导波信号的基本链路如数据导入、时域波形观察、位移分量对比等还可作为进一步开展傅里叶频谱分析、频散特性研究或缺陷回波识别的可扩展基础脚本。1. 拿到P2_U1_xy_displacement时先别急着跑超声导波在板类结构里传播时质点位移同时存在面内和面外两个分量。多数实验室习惯只看单个方向的A扫信号但遇到板边反射、模态转换或缺陷散射时单一方向图很难解释为什么信号后面跟着一串多余的波包。P2_U1_xy_displacement这份数据把x、y两个方向的位移分量分别存成.mat文件再用Excel时间节点表对齐物理时间目的是把每个时间帧上的二维位移场完整重建出来。对做结构健康监测和无损检测的工程师、研究生而言它的价值在于演示了一条“原始位移数据到预处理、到时域位移图、再到频域参数”的完整链条。后面会逐个拆解文件组织方式、脚本各段代码的物理含义以及最容易被忽略的参数坑。2. 先把数据拆开x.mat、y.mat和Excel时间表里存的是什么2.1 超声导波xy位移数据的物理含义在板结构中激励超声导波时某个空间点的质点位移可以被分解为沿x和y方向的分量。如果数据来自激光测振仪扫查x.mat和y.mat记录的通常是测量网格上每个点随时间变化的位移幅值维度一般是“时间帧数乘以空间点数”或者重排成“时间帧数乘Ny乘Nx”的立方体如果来自Abaqus或COMSOL这类有限元模拟节点位移导出后也常被整理成这种格式。P2_U1这个编号大概率表示第二批实验的第一个激励-接收配置U1可能是换能器编号“xy_displacement”则直接点明数据内容是两个正交方向的位移。这个背景直接影响后续绘图方式surf和pcolor这类函数需要的是二维网格脚本里要先确认数据排列顺序再reshape成矩阵否则画出来就是一团乱线。2.2 用whos和readtable先摸清维度很多人拿到压缩包后直接load就开始跑这一步其实很容易埋雷。load会把数据一股脑塞进工作区变量名被覆盖、维度对不上时排查成本很高。我一般先做两个侦察动作把文件的内部结构看清楚再动手% 查看x.mat内部有哪些变量以及每个变量的尺寸不载入工作区 infoX whos(-file, x.mat); disp(infoX); % 读取Excel时间节点表确认物理时间轴的采样点数量和单位 tt readtable(P2时间节点数据-U1.xlsx); disp(head(tt, 5)); disp(height(tt));这段代码的思路是先用whos的-file选项列出文件内的变量名、大小和类型不污染当前工作区再用readtable读取Excelhead(tt,5)看前五行的列名和数值格式确认时间单位是秒还是微秒。height(tt)返回行数这个数字需要和x.mat第一维度的帧数严格相等否则后面做时间轴对齐时会出现索引越界或错位。提示如果whos显示x.mat里是两个向量比如time和ux说明数据是“时间序列加单点位移”的A扫格式如果显示的是三维数组说明是扫查网格的面数据。这两种情况脚本分支完全不同先分清维度再谈滤波。参数说明readtable默认把Excel第一行当变量名P2时间节点数据-U1.xlsx如果第一行是说明文字读进来后列名会变成Var1、Var2这时可以用readtable(..., VariableNamingRule, preserve)保留原始表头。时间列是微秒还是毫秒决定了后面计算采样率时要不要乘10的负六次方这个细节错一位数群速度结果就会差三个数量级。2.3 维度不一致时的对齐策略实际数据里经常出现帧数对不上的情况x.mat的帧数是Ny.mat的帧数也是N但Excel时间节点表只有N-1行最后一行被手动删掉了或者xy_u1.mat存的是复数形式的合成位移实部虚部分别对应x和y分量取的时候要用real和imag拆开。常见的做法是先统一帧数如果时间节点数少一帧就在末尾补上最后一个采样间隔或者截断位移数据里多出来的那一帧以Excel为准。frames height(tt); % 以Excel时间节点数为准 if frames size(xData, 1) xData xData(1:frames, :, :); yData yData(1:frames, :, :); fprintf(截断位移数据至 %d 帧\n, frames); elseif frames size(xData, 1) error(时间节点数(%d)多于位移帧数(%d)检查源数据导出范围, frames, size(xData, 1)); end这里以物理时间表为准而不是以.mat文件为准因为时间表是实验时记录的硬件采样节拍位移数组可能因为导出时多写了几帧缓冲而比实际长。逻辑上先做比较再截断最后用error兜底报错避免静默错位。若xy_u1.mat里的复数位移需要拆分务必备份一份原始值再进行real、imag运算因为后续去趋势和滤波会修改数组原始信息丢了想再找回就得重新解压。3. 主脚本解析把位移数据变成时域位移图3.1 载入与时间轴对齐的实现P2_U1_xy_displacement.m这个脚本的核心工作是把多个来源的数据统一到同一个时间基准上。常见的做法是先把三个.mat文件和一个Excel读进来再根据Excel里的物理时间构造出等间隔的时间向量。这段逻辑写在脚本头部作用是用%注释标记每个变量来源后续维护时能一眼看出哪一段数据对应哪个实验条件。% 载入位移数据假设x.mat中变量名为uxy.mat中为uy load(x.mat); load(y.mat); load(xy_u1.mat); % 可能是复合位移或U1接收点数据 % 读取时间节点表时间单位假设为微秒 tt readtable(P2时间节点数据-U1.xlsx); t tt.Time_us; % 若是毫秒单位统一转换成微秒 if max(t) 1e3 t t * 1e3; % 从毫秒换算到微秒 end t t - t(1); % 归零起点便于后续FFT % 若xy_u1是复数合成位移拆出实部虚部作为x/y分量 if ~exist(ux, var) exist(xy_u1, var) ux real(xy_u1); uy imag(xy_u1); end逻辑说明load不加输出参数时会把.mat里的所有变量直接丢进工作区所以脚本里必须先确认变量名避免ux和uy被其他同名变量覆盖。时间归零这一步是为了让FFT的频率轴从0开始不归零也不会出错但会多出一个由初始相位带来的直流偏置。用exist检查变量是否存在是为了兼容两种数据来源如果x.mat已经给了ux就跳过拆分分支没有的话才从xy_u1的实部虚部里取。参数说明t(1)归零后采样间隔是通过diff(t)计算的理想情况下所有间隔相等但实验硬件偶尔丢帧这时要检查median(diff(t))和实际间隔的偏差偏差超过1%就说明时间轴不是均匀的需要先做插值重采样否则后面所有基于固定采样率的运算都会失真。3.2 去趋势与带通滤波的参数怎么定原始位移数据里总混着零漂和噪声。激光测振仪和压电传感器输出经常带一个固定偏置不做去趋势会直接反映在位移图上表现为整体抬升或缓慢漂移掩盖真实的波包。滤波是另一个关键步骤导波实验一般用带通滤波器把激励频率附近的频带留下来其余全部切掉。% 从时间轴反推采样率单位Hz Fs 1e6 / median(diff(t)); % t为微秒1e6是微秒到秒的换算 % 去直流分量 xData detrend(xData, constant); yData detrend(yData, constant); % 带通滤波频带按激励频率调整 fLow 20e3; % 下限20kHz fHigh 200e3; % 上限200kHz [b, a] butter(4, [fLow fHigh] / (Fs / 2), bandpass); xDataF filtfilt(b, a, xData); yDataF filtfilt(b, a, yData);逻辑说明detrend的第二个参数constant表示只去除常量偏置不去除线性趋势因为导波位移数据本身就有缓慢变化的包络直接用linear会把低频成分误伤。butter(4,...)是四阶巴特沃斯滤波器通带内最平坦适合保留导波信号的波形特征filtfilt做零相位滤波正向和反向各滤一次互相抵消相位延迟这样波包到达时间不会被滤波器提前或拖后。参数说明fLow和fHigh的选择有依据5周期汉宁窗调制的100kHz激励信号频谱能量集中在70kHz到130kHz之间取20k到200k是给模态频散留了余量。如果激励频率是50kHz频带应该缩到10kHz到100kHz太宽的频带会放进低频结构振动太窄会把导波的旁瓣削掉时域上表现为波形前后来回振荡。滤波完成后建议做一次前后对比绘制原始信号和滤波信号的叠加图确认波包位置没有移动、幅值没有被压扁。3.3 绘制位移场与A扫的完整流程时域位移图分两种画法一是看某个固定时刻的全场分布二是看某个固定点的幅值随时间变化。前者用imagesc或pcolor后者用plot。两者数据索引方式不同但脚本里通常把两种图放在同一个figure里做联动方便对照。% 提取某个时间帧的二维位移场frameIdx为时间索引 frameIdx 50; uxFrame squeeze(xDataF(frameIdx, :, :)); % 去掉长度为1的维度 figure; imagesc(xGrid, yGrid, uxFrame); % xGrid,yGrid为空间坐标网格 axis xy; axis equal tight; colorbar; colormap(parula); caxis([-maxAbs maxAbs]); % 固定色标范围便于逐帧对比 xlabel(x / mm); ylabel(y / mm); title(sprintf(u_x位移场, t %.2f us, t(frameIdx))); % 单点A扫:probeIdx为空间索引 probeIdx 120; figure; plot(t, xDataF(:, probeIdx), b-, t, yDataF(:, probeIdx), r-); legend(u_x, u_y); xlabel(时间 / us); ylabel(位移); grid on;逻辑说明squeeze的作用是把维度为1的轴去掉因为xDataF是“时间乘以空间”的二维数组取单帧后变成一维向量如果原始数据是三维网格取单帧后是Ny乘Nx矩阵同样需要squeeze处理。imagesc和pcolor的最大区别是缩放方式不同pcolor按顶点着色会丢最后一圈数据imagesc按网格中心着色更适合位移云图。caxis固定色标范围是关键技巧默认的自动色标会让每张图的颜色都按自己最大幅值缩放两帧之间看起来对比度相同实际幅值却差很多误导性很强。A扫代码里probeIdx这个索引指向某个具体空间位置在波场图里点选这个位置就能看到该点x和y两个方向位移的完整时间历程。实际实验里x方向分量对应面内位移y方向分量对面外位移更敏感两个方向波包到达时间的差值可以直接用来判断模态类型。图类型常用函数用途注意事项A扫plot(t, u)单点振幅分析、到时读取配合hilbert包络找峰值更稳B扫imagesc(t, 位置, u)沿某条线看传播过程空间轴要等距采样波场快照imagesc/pcolor某一时刻全场位移分布固定caxis再逐帧比较4. 从时间域到频域FFT、频散曲线与群速度提取4.1 为什么时域图不够用时域位移图只回答了“波什么时候到达”这个问题但导波是多模态的同一频率下A0和S0模态同时存在传播速度不同、频散程度不同两个波包在时域上可能重叠肉眼根本无法分开。把信号变到频域才能看清每个频率成分的幅值和相位。对P2_U1这类数据典型的后续分析是沿某个传播方向取一条线上的位移数据做空间-时间二维傅里叶变换得到频率-波数图从图上直接读出各模态的频散关系。这一步做扎实后面计算群速度和相速度才有依据。4.2 逐点FFT与中心频率确认对每个空间点的时间序列做FFT得到的频谱可以确认激励信号的中心频率是否在预期范围内。P2_U1_xy_displacement.m里如果包含频谱分析部分常见的写法是这样N length(t); % 时间采样点数 dt median(diff(t)) * 1e-6; % 换算成秒 Fs 1 / dt; % 采样率单位Hz % 加汉宁窗抑制频谱泄漏 win hann(N); spec fft(xDataF(:, probeIdx) .* win); spec spec(1:floor(N/2)); % 取单边谱 f (0:floor(N/2)-1) * Fs / N; % 找峰值频率作为实测中心频率 [~, imax] max(abs(spec)); fc f(imax); fprintf(实测中心频率: %.2f kHz\n, fc / 1e3);逻辑说明hann窗是处理导波信号的标准选择它的主瓣比矩形窗宽但旁瓣衰减快得多能有效防止信号截断产生的频谱泄漏把真实频率峰值淹没。乘窗之后频谱能量会摊到相邻频点所以后续找峰值时用的是加窗后的结果这个峰值频率和激励设定值之间的偏差如果超过3%就要检查采样率估算是否有误或者激励源是否失谐。参数说明频率分辨率df Fs / N这个值决定了频谱上的最小可分辨间隔。P2_U1这份数据如果N是1024、采样率是10MHzdf约等于9.8kHz对100kHz的导波来说分辨率偏粗想提高分辨率有两个办法一是补零到更长的FFT长度这只是插值不能增加真实信息二是采集更长时间的数据这才是真正提高分辨率的路径。实测中心频率对后续滤波参数有直接影响若fc明显偏离fLow到fHigh的中心说明带通滤波器参数需要重新设计。4.3 群速度计算希尔伯特包络与峰值到时群速度是导波检测里最常被提取的物理量。原理很简单两个接收位置之间的距离除以两个波包到达时间之差。难点在于“到达时间”怎么定义才可靠。导波经过频散后波形会展宽直接用幅值最大点作为到时结果受频散影响很大。更稳定的做法是用希尔伯特变换提取包络再找包络峰值对应的时刻。% 两个接收位置的一维时间信号 u1 xDataF(:, idxPos1); u2 xDataF(:, idxPos2); % 希尔伯特包络 env1 abs(hilbert(u1)); env2 abs(hilbert(u2)); % 找包络峰值对应的时间索引 [~, i1] max(env1); [~, i2] max(env2); t1 t(i1); t2 t(i2); % 群速度 距离差 / 时间差 dist abs(pos2 - pos1); % 单位mm cg dist / ((t2 - t1) * 1e-6) / 1e3; % 换算成km/s fprintf(群速度: %.3f km/s\n, cg);逻辑说明hilbert把实信号变成解析信号取绝对值得到包络这个过程相当于把信号的瞬时能量轮廓提取出来。包络峰值比原始信号峰值更稳定因为它抹掉了载波振荡带来的假峰。速度计算里两次单位换算要特别小心时间t的单位是微秒乘1e-6转成秒距离单位是毫米除以1e3从m/s换成km/s这样得到的结果和文献里常见的单位对齐。模态典型频率范围速度特征频散特性A010k-300kHz慢约2-3km/s强频散速度随频率下降S010k-300kHz快约5-6km/s弱频散高频略快SH0板平面内剪切中速约3km/s无频散理论速度恒定实际数据中A0和S0可能同时出现包络峰值法会优先锁定幅值最大的那个。如果速度算出来明显偏高或偏低不要急着改代码先画一下两个位置各自的时域波形确认是否真的对应同一个模态。必要时在FFT后先做窄带滤波只保留某个目标频率分量再对该分量做希尔伯特包络这样算出来的群速度就是该频率下的群速度可以和理论频散曲线直接对比。5. 几个排错技巧窗函数、滤波边界与二维FFT模态识别5.1 滤波器振铃带通滤波器阶数太高或带宽太窄时输出信号会出现明显的振铃现象波包前后出现幅度逐渐衰减的假振荡看起来像新波包其实是滤波器冲击响应的尾巴。检查方法是输入一个孤立的脉冲信号观察滤波输出尾部是否干净。四阶巴特沃斯是兼顾陡降和低振铃的常用折中六阶以上慎用fHigh和fLow的比值小于2时振铃概率显著增加。5.2 窗函数对频谱的影响窗函数主瓣宽度旁瓣衰减适合场景矩形最窄-13dB频谱峰值定位易泄漏汉宁中等-31dB常规导波频谱分析布莱克曼最宽-58dB弱信号附近的频谱识别频谱泄漏的本质是信号截断带来的突然跳变泄漏强的频段可能掩盖真实存在的弱模态。做完FFT后如果发现频谱基底抬高、出现均匀分布的梳状纹路说明窗函数选得不够长或者窗类型不对。导波信号本身是窄带的汉宁足够需要识别两个频率非常接近的模态时矩形窗的主瓣窄优势才值得利用。5.3 二维FFT识别模态成分把沿传播方向的一条线数据排成“时间乘以空间”的矩阵对时间和空间两个维度分别做FFT得到频率-波数图不同模态在图上落在不同的色斑带上。MATLAB里用fft2一行就能完成% uxt为时间乘空间的二维矩阵 spec2D fft2(uxt); spec2D fftshift(abs(spec2D)); % 零频移到中心 f_axis (0:size(uxt,1)-1) * Fs / size(uxt,1); k_axis (0:size(uxt,2)-1) * 2*pi / (dx * size(uxt,2)); imagesc(k_axis, f_axis, log10(spec2D 1)); % 取对数让弱模态可见 xlabel(波数 k / rad·mm^{-1}); ylabel(频率 f / kHz);代码里的dx是空间采样间隔fftshift把零频移到坐标中心后正负波数分别对应两个传播方向图上对称出现的亮斑就是同一模态向不同方向传播的结果。取对数这一步很关键直接显示线性幅值时强模态会把弱模态压没log10后弱模态的细节才能浮现出来。把理论频散曲线换算成频率-波数关系叠加在图上就能定性判断数据里主要激励起了哪些模态P2_U1这份数据如果激励频率设定不当导致多个模态混叠在这个图上一眼就能看出来。本文还有配套的精品资源点击获取