PIV实验做完拿着一堆互相关算出来的速度矢量很多人第一反应是直接出图。但矢量图密密麻麻的箭头或者只用颜色显示速度大小的散点图其实都不够直观。真正能拿出去展示、能用来分析流场结构的往往是流速云图。而MATLAB里画云图的主力函数就是contour尤其是contourf。这篇东西就把我从PIV数据处理到最终出云图的完整经验整理出来包括contour函数的各种细节坑适合刚接触PIV数据后处理、被一堆矢量文件搞得头疼的同学参考。1. 从PIV实验到流速云图为什么云图是最后一步的关键1.1 PIV实验数据到底长什么样PIVParticle Image Velocimetry粒子图像测速的核心原理不复杂在流场里撒示踪粒子用激光片光照明用高速相机连续拍两张或一组照片然后对图像划分成许多查询窗口对每个窗口做互相关计算找到粒子团在Δt时间内的位移除以时间间隔就得到速度矢量。但真正拿到手的数据和教材上画的简图差距很大。一套常规二维PIV实验做下来你会得到这样的东西一个文件夹里塞满了几百上千对图片经过Davis、PIVlab、OpenPIV这类软件计算后导出的是若干组数据文件。常见的导出格式是xyuv四列有时候还带一个表示矢量质量的参数比如信噪比、相关性峰值。这里的x和y是每个查询窗口中心的空间坐标u和v是对应位置的流向速度和法向速度。也就是说你手里的数据本质上是分布在离散网格上的矢量场不是一张现成的云图。横竖网格点数量通常取决于查询窗口大小和重叠率。比如相机分辨率是2048×2048像素查询窗口是64×64步长是32像素那算出来大概是63×63个矢量点。这些点在空间上虽然规则但因为实验过程中可能存在标定系数变化、坐标方向设定等问题直接拿去画图经常出各种幺蛾子。1.2 云图和矢量图怎么选别一上来就裸画箭头有些同学拿到数据后第一件事就是quiverxyuv画出一大片箭头。箭头图适合看方向比如旋涡的旋转方向、剪切层的流动方向确实一目了然。但箭头图有两个很明显的问题。第一个问题是信息密度太低。63×63个矢量点画出来的箭头密密麻麻屏幕上全是黑色的线条和箭头尾巴根本分不清哪里速度快哪里速度慢。你得把图放大一个个箭头去比长度非常费眼睛。第二个问题是无法直观反映速度大小的空间分布。流速云图等值线图恰恰解决了这个需求横纵坐标代表空间位置颜色或等值线代表速度大小整个流场的高低速区域一眼就能看出来哪里是射流核心区、哪里是回流区、哪里是边界层低速带全部清清楚楚。所以我的习惯是先用云图看整体分布确定关键区域后再叠加箭头图看方向两张图配合使用。这也是为什么contour函数在PIV后处理里几乎是必用工具。2. contour函数的核心细节与参数玩法2.1 contour与contourf一字之差效果天差地别MATLAB里画等值线的函数有一对孪生兄弟contour和contourf。很多新手第一次用直接contourxyU画出来的图只有一圈一圈的黑线没有颜色填充看起来像地形图而不是云图。contour函数默认只画等值线线条本身不带颜色填充。如果你想要那种颜色连续变化的云图效果两个办法用contourfxyUlevelsf是fill的意思会在相邻等值线之间填充颜色。用contourxyUFillonLineColornone这是老版本MATLAB的常用写法效果和contourf类似。新版MATLABR2014b之后里contourf的行为也有一些调整。现在contourfxyU默认画出来的是填充的云图但注意不同版本对levels参数的解释略有差异。我建议明确写levels不要靠默认值。另外还有一个常见误区contourf默认会在等值线之间画出黑色的边界线这在云图上会显得很杂乱。通常我们不需要这些黑线所以我画云图时基本固定搭配contourf(xq, yq, U_mag, 50, LineColor, none);如果你还在用特别老的MATLAB版本R2014a之前可能不支持LineColor参数那就用contourf(xq, yq, U_mag, 50, LineStyle, none);或者直接[c,h] contourf(xq, yq, U_mag, 50); set(h, EdgeColor, none);这三种写法在功能上是等价的选一种自己习惯的就行。我的建议是记住第一种够用了。2.2 levels参数等值线数量和控制列表levels这个参数是contour函数里最容易让人困惑的一个点。它有两种完全不同的含义取决于你传进去的是单个数字还是向量。如果你传的是单个数字比如contourfxyU30意思是让MATLAB自动生成30条等值线。这个数字不一定是精确的等值线数量更准确地说是颜色的分层数。MATLAB会在数据的最小值和最大值之间自动划分30个颜色级别。如果你传的是向量比如contourfxyU[0.10.020.5]意思是让MATLAB在指定的数值处画等值线。这个等值线列表可以是不等间距的完全由你控制。实际做PIV云图时我强烈建议用向量方式不要用数字。原因很简单如果你用contourfxyU30颜色范围是由当前数据min和max决定的。当你对比不同实验工况时每个图的颜色范围都不一样看起来就无法直接对比。统一用clim把色标范围固定之后再用数字levels也可以但最好还是用向量把等值线间距设成一个固定步长图与图之间才有可比性。举例说明假设你要对比三个不同流速下的流场速度范围分别是0~2 m/s、0~4 m/s和0~6 m/s。如果都用30条等值线第一张图的颜色从0到2渐变第三张图从0到6渐变视觉上完全没法对比。正确的做法是设定统一的levels向量levels 0:0.1:6; contourf(xq, yq, U1, levels, LineColor, none);这样三张图的颜色级别完全对应放一起看才有意义。2.3 颜色映射与色标范围云图好看的关键级别和颜色范围的设定直接决定了云图的信息表达。这里有三个东西要分开理解colormap、clim老版本叫caxis、colorbar。colormap是颜色映射表MATLAB里最常用的包括parulaMATLAB默认的颜色映射蓝到黄渐变颜色区分度好对色盲相对友好做科研出图推荐。jet经典的彩虹色蓝→青→绿→黄→红视觉冲击强但存在两个问题黄色区域亮度过高容易造成视觉假象且两端深蓝色和深红色对于色盲人群不友好。现在很多期刊已经明确不推荐用jet。turbo新版MATLAB提供的jet替代方案比jet更平滑色彩分布更均匀。hot黑→红→黄→白适合热场。cmocean系列专门为海洋和大气科学设计的颜色映射比如cmoceanbalance、cmoceanthermal如果装了Climate Data Toolbox可以直接用效果非常专业。我个人的习惯是内部处理用jet方便观察细节出正式图用parula或turbo投期刊再按期刊要求调整。在代码里用colormap命令设置colormap(parula); % 或者 colormap(turbo); % 或者自定义 cmap cmocean(thermal); colormap(cmap);clim控制的是颜色映射的范围也就是数据的最小值和最大值对应到色标的最小值和最大值。很多时候数据里有异常的大值或小值如果不限制clim云图的颜色范围会被这些异常值拉得很宽导致主体区域颜色没有区分度。这时候就得手动固定clim([0 5]); % 新版本写法 % 老版本用 caxis([0 5]);固定clim之后数据超出范围的部分会显示成最红或最蓝的颜色这样反而能把异常区域标出来方便你回到原始数据排查。2.4 从离散矢量到连续场网格化的关键一步这一步是PIV云图最容易翻车的环节也是我一开始踩坑最多的地方。PIV软件导出的xy坐标虽然来自规则网格但由于标定、畸变矫正等原因这些点往往不是严格等间距的矩形网格。而contour函数对输入数据格式是有要求的它需要的是矩阵形式的数据也就是xq、yq各自是网格坐标矩阵U_mag和它们维度一致。很多同学直接做contourxyU_mag发现报错说维度不对或者画出来的图乱七八糟就是因为x、y、U_mag都是向量而不是矩阵。解决办法就是对原始离散点做网格插值。MATLAB里常用的插值函数有两个griddata老牌函数用法是[xqyq] meshgrid(xiyi)Uq griddataxyU_magxqyqlinear。简单直接缺点是复杂数据时速度较慢。scatteredInterpolant推荐使用速度快且内存效率高。用法F scatteredInterpolant(x, y, U_mag, natural, linear); Uq F(xq, yq);scatteredInterpolant支持三种插值方法linear线性插值速度快但会在数据点周围产生折线效果、natural自然邻域插值效果好光滑适合PIV数据、nearest最近邻插值速度快但结果块状感强。对于PIV数据我推荐用natural它能比较好地保留流场细节又不会像linear那样出现明显的插值痕迹。如果你用的是griddata那对应的方法是v4效果也还可以但数据量大的时候非常慢。还有一点要注意插值的目标网格范围不能超出原始数据的凸包范围否则插值结果是NaN云图上会出现空洞。这个我在后文的常见问题部分再细说。3. 完整实操从原始矢量场到一张能进论文的流速云图3.1 准备数据与坏矢量处理假设你手里已经有一个从PIV软件导出的数据文件比如是CSV格式包含x、y、u、v四列。第一步自然是读进来data readmatrix(piv_data.csv); x data(:,1); y data(:,2); u data(:,3); v data(:,4);但这里有一个大坑PIV互相关计算出来的原始矢量场几乎必定含有坏矢量。所谓坏矢量就是明显错误的计算值可能是由于粒子图像质量差、窗口内粒子数量少、跨帧丢失等原因造成的。这些坏矢量的典型特征是速度大小异常大、方向突变和周围矢量完全不一致。如果不处理坏矢量直接插值后果非常严重。因为插值会把坏矢量的错误值扩散到周围区域云图上会出现一块假的低速或高速区整个流场看起来就像长了疮一样。所以读入数据后必须先做坏矢量检测和剔除。我常用的方法有三个第一个是速度阈值法。根据实验条件估计一个合理的速度上限。比如水管实验流速估计最大不超过2 m/s那么abs(sqrt(u.^2v.^2)) 3就可以判断为坏矢量。speed sqrt(u.^2 v.^2); max_valid 3; % 根据实验情况调整 valid speed max_valid;第二个是邻域中值法。把每个矢量点与周围8个邻近点的中值比较如果偏差超过中值的两到三倍判定为坏矢量。这个在PIVlab里叫standard deviation filter或者median filter。第三个是最原始的人工检查。把矢量场画出来用鼠标一个个把明显的错误箭头挑出来。这个办法适用于矢量点很少的情况数据量大了不现实。在实际处理中我一般先用阈值法粗筛再对剩余的矢量做一次邻域中值滤波双保险。处理完之后被剔除的位置暂时是NaN后面插值时把这些点排除掉u_clean u(valid); v_clean v(valid); x_clean x(valid); y_clean y(valid);3.2 插值、速度合成与平滑坏矢量剔除后就可以开始插值了。推荐用scatteredInterpolant因为它对NaN点的处理比较自然。先把有效点构造成插值对象F_u scatteredInterpolant(x_clean, y_clean, u_clean, natural, linear); F_v scatteredInterpolant(x_clean, y_clean, v_clean, natural, linear);然后定义插值网格。注意网格的定义范围要尽量贴近原始数据的有效区域不要随便包一个巨大的范围否则边缘会有大片NaNxi linspace(min(x_clean), max(x_clean), 100); yi linspace(min(y_clean), max(y_clean), 100); [xq, yq] meshgrid(xi, yi); uq F_u(xq, yq); vq F_v(xq, yq);插值完成之后计算速度大小U_mag sqrt(uq.^2 vq.^2);这里还可以顺带计算涡量用MATLAB的curl函数[curlz, cav] curl(xq, yq, uq, vq);涡量云图在某些场合比速度云图更能反映流场结构比如边界层分离、旋涡脱落涡量云图一看就明白。插值之后还有一个细节值得注意要不要做平滑。PIV原始数据本身带有一定的噪声插值后这些噪声会被保留下来云图上表现为颜色过渡不自然、有颗粒感。我一般会做一次二维平滑最简单的是移动平均但要注意不能过度平滑导致小尺度结构被抹掉。用imgaussfilt做高斯平滑是最省事的U_mag_smooth imgaussfilt(U_mag, 1.5);不过要记住平滑是一种修图手段不要为了好看而失真。论文里如果用平滑最好在方法部分写清楚。3.3 绘制云图与叠加矢量数据准备好了画图本身其实很快figure; contourf(xq, yq, U_mag, 50, LineColor, none); hold on; quiver(xq(1:5:end,1:5:end), yq(1:5:end,1:5:end), ... uq(1:5:end,1:5:end), vq(1:5:end,1:5:end), 2, k); hold off;quiver叠加箭头时注意两个技巧第一个技巧是抽稀。云图本身已经足够密了箭头不需要每格都画否则黑压压一片。用15end每隔5个点画一个箭头即可具体间隔看网格密度调整。第二个技巧是颜色。箭头用黑色或者白色比较合适在彩色云图上识别度最高。如果用红色蓝色反而会和云图颜色混淆。箭头密度和长短也要调可以用quiver的缩放参数比如上面的2表示箭头自动缩放2倍调一次到看起来疏密合适为止。如果不想用quiverMATLAB的streamline函数可以做流线图在某些场景下比箭头更优雅。但对PIV数据来说streamline需要在规则网格上工作我们的插值网格正好满足要求h streamline(xq, yq, uq, vq, startx, starty); set(h, Color, k, LineWidth, 0.5);画出来是流线叠加云图效果非常专业。startx和starty是流线起点的坐标可以按需设置几个。3.4 坐标、色标与出图设置云图画出来后还有一系列细节需要处理。很多PIV数据在采集时X方向是流向水平Y方向垂直向下因为相机图像坐标习惯从左上角开始。但流场分析的惯例是Y轴向上所以画出来的图经常是上下颠倒的。解决办法很简单set(gca, YDir, reverse);或者反过来如果你希望图像和实际流场对应用set(gca, YDir, normal)。这个方向问题没有绝对标准关键是确保你知道自己的Y轴方向代表什么标注清楚。接下来设置色标colormap(parula); clim([min_val max_val]); colorbar; ylabel(colorbar, Velocity (m/s));这里重点强调统一范围的问题。你如果要出一组不同工况的图务必用同一个clim否则放在一起没有可比性。我通常在代码前部定义一个变量v_range [0 2]; % 统一速度范围 % 每张图都写 clim(v_range);坐标轴单位换算也要留意。PIV软件输出速度的单位取决于标定设置常见的单位有m/s、mm/s、pixel/frame。如果是pixel/frame还需要乘以标定系数和时间间隔才能换算成实际速度。换算公式是U_mag_ms U_mag * calib_mm_per_pixel / dt_s / 1000;calib_mm_per_pixel是每个像素对应的毫米数dt_s是两帧之间的时间间隔秒除以1000是把毫米换算成米。这一步很容易被忽略导致云图上的数值量级完全不对。坐标轴的标注文字也不要偷懒。xlabel、ylabel写清楚物理量和单位标题里写上实验工况编号。这个习惯在写论文时能省很多事。最后是出图导出高清图片有两个办法% 老办法 print(gcf, result.png, -dpng, -r300); % 新办法 exportgraphics(gcf, result.png, Resolution, 300);exportgraphics是R2020a之后引入的更好地保留了图的原始样式推荐使用。300 dpi是期刊投稿的最低标准如果要放大打印建议600 dpi。4. 常见问题与排查技巧实录4.1 云图出现空洞或锯齿状边界这个问题的根本原因是插值范围超出了原始数据的凸包范围。contourf对NaN数据很敏感如果你的插值网格上有大片NaN云图就会出现空洞。排查思路先用isnan检查插值结果中有多少NaNsum(isnan(U_mag(:)))。检查插值网格范围是否超出原始数据范围。如果确实需要在整个矩形区域显示云图可以考虑用scatteredInterpolant的linear方法并配合外推但是要谨慎外推出来的数据并不代表真实流场。解决方法是把网格范围缩小到数据覆盖范围内。还有一个办法是生成网格后把数据范围以外的点直接置为NaN再用contourf绘制得到有边界的云图。4.2 云图颜色乱七八糟和速度对不上这种情况十有八九是clim没设置。特别是数据里有几个坏矢量没清干净的情况下速度最大值可能是一个几十甚至几百的错误值导致整张图颜色全部被这个异常值拉走真实流场区域全是同一个颜色。处理思路很简单设置合理的clim比如clim([0 prctile(U_mag(:)98)])用98分位数作为上界能保留大部分有效信号又不会被尾部极值干扰。4.3 坐标轴方向反了图上下颠倒这个很常见。PIV软件导出的Y方向可能是从上往下也可能从下往上取决于你的标定设置。如果你发现云图的上下和实验流场的上下对不上直接set(gcaYDirreverse)翻转即可。如果在某些情况下坐标轴的tick标签也需要翻转可能需要手动设置YTickLabel。一个更隐蔽的问题是有些同学把图像坐标系的y和PIV物理坐标系的y搞混了。图像坐标系是左上角为原点y向下物理坐标系是左下角为原点或者根据实际布置定义y向上。画图前先确认数据里的y到底是哪种避免出现瀑布往下流但云图上显示往上的笑话。4.4 性能优化与批量出图试验数据一多比如要画100个工况的云图如果每个都手动调参数工作量巨大。我的习惯是写一个批处理脚本把所有数据文件放在一个文件夹里循环处理自动出图。files dir(*.csv); for k 1:length(files) data readmatrix(files(k).name); % ... 处理、插值、画图 ... exportgraphics(gcf, [files(k).name(1:end-4) _cloud.png], ... Resolution, 300); end批量出图要注意内存管理。每次循环结束用close all关掉图形窗口否则积攒几十个figure会把内存耗尽。另外如果数据量特别大比如几百MB的PIV数据建议用matfile或者只读入所需列避免MATLAB卡死。我试过一次性读入5GB的CSV文件MATLAB直接内存不足。后来改用readmatrix读取一列处理完再释放就流畅多了。最后再分享一个提高出图速度的小技巧在做云图调试阶段先用低分辨率网格比如50×50快速画图调整好表现参数后再用高分辨率200×200出正式图。不要一上来就用超高分辨率否则每调整一次参数都要等半天非常浪费时间。