☰
PIV数据到MATLAB流速云图:从散点网格化到contourf绘制全流程
2026/10/5 6:01:25 网站建设 项目流程

我见过太多做PIV实验的同学,拿着处理好的数据,信心满满地敲下contour(x, y, speed),结果MATLAB直接报错。这不是手滑,而是没搞清楚contour函数的输入约束:它要的是网格化的矩阵,不是散点向量。PIV 实验做完之后拿到手的流场数据,几乎都是规则或不规则散布的速度矢量,想让这些矢量变成论文里那种红蓝渐变、层次分明的流速云图,中间隔着的就是网格化插值和对 contour 系列函数的正确理解。这篇文章把我自己从 PIV 原始数据到发表级流速云图的完整流程梳理一遍,重点围绕 MATLAB 的 contour / contourf 函数,把散点网格化、插值方法选择、配色、矢量叠加这些绕不开的环节一次性讲透。适合刚接触 PIV 数据处理的研究生,也适合想系统提升云图质量的科研人员。

1. PIV实验之后的数据长什么样?画云图前的准备工作

1.1 拿到手的数据本质是散点矢量

做PIV实验时,我们往水槽或风洞里撒入示踪粒子,用双脉冲激光片光源照亮目标切面,高速相机在相隔Δ t的两个时刻拍摄两张粒子图像,之后通过互相关算法把图像划分成一个个查询窗口(interrogation window),在每个窗口内统计粒子群的整体位移,再除以Δ t就得到该窗口中心点的速度。所以PIV的原始结果从来不是一张连续的场,而是一堆离散的速度矢量,每个矢量对应当前窗口中心位置。

以我常用的PIVlab为例,处理完一组图像后导出的就是几列数据:x、y、u、v,有的还会附带每个矢量的信噪比、残差等质量参数。用代码读进来之后,几个关键信息要第一时间确认:坐标是像素还是物理坐标,速度单位是m/s还是mm/s,数据总共有多少个有效矢量。这几个参数决定了后面所有操作,很多人出图变形或者数值看起来离谱,源头都在这里。

% 以PIVlab导出的csv为例 data = readmatrix('piv_result.csv'); x = data(:,1); % 坐标,可能是px,也可能已标定为mm y = data(:,2); u = data(:,3); % 速度分量 v = data(:,4); disp(size(data)); % 看看有多少个矢量

1.2 坏矢量剔除:别把脏数据画进云图

PIV后处理出来的矢量并不都可靠,常见问题包括:查询窗口内粒子太少导致相关峰不突出、激光反射或壁面眩光导致的错误矢量、图像局部过曝区域周围的速度异常。这些坏矢量直接画进云图,会在流场里产生刺眼的色斑,严重时会扭曲整个流动结构的判断。

所以画云图前先做一遍粗筛。我自己的习惯是用速度中值加标准差剪切。这个方法朴素,但效果很明显:

speed = sqrt(u.^2 + v.^2); % 合速度 mu = median(speed); % 中值更鲁棒,比均值好 sig = std(speed); valid = abs(speed - mu) < 3*sig & ~isnan(u) & ~isnan(v); x = x(valid); y = y(valid); u = u(valid); v = v(valid);

阈值选择要根据实验工况微调,不要试图一个参数通吃所有数据。比如湍流度高的区域,真实速度脉动本身大,3倍标准差会把一些有效矢量误删;层流区则可以把阈值收紧到2倍。我的习惯是先画一版散点矢量图,肉眼确认哪些位置明显异常,再据此调整阈值,不要盲删。

1.3 画云图前必须想明白的三件事

第一,单位统一。同一张图里既有毫米又有米,云图比例会失衡,速度矢量方向也会跟着乱。PIV软件导出的坐标常见为mm,而速度是m/s,在MATLAB里建议全部转成国际单位后处理。

第二,区域截取。实际流场中常有模型壁面、遮挡区,这些位置本身没有有效数据,画图之前先设定好物理范围,或者做掩膜,不要让插值算法在无数据区硬算。

第三,明确要画什么量。速度云图最常用的是合速度speed = sqrt(u.^2 + v.^2),因为它反映流场的整体强度分布。但如果你关心的是管道轴向速度、射流中心线速度这类方向性特征,就单画u分量,并在图和colorbar里明确标注。

2. contour函数到底怎么用?从简单等值线到平滑云图

2.1 contour的基础语法与维度约束

contour函数做的是"等值线图",给定二维网格上的一个标量场Z,画出Z值相同的点连成的线。基础语法有几种:

contour(Z) % 只给Z矩阵,自动选择等级 contour(X, Y, Z) % 给空间坐标网格和数据场 contour(Z, n) % 指定n条等值线 contour(Z, v) % 指定具体等值线值向量v [C, h] = contour(...) % 返回等值线矩阵C和对象句柄h

这里最大的坑是X、Y和Z的维度关系。当X、Y是向量时,必须满足length(X) == size(Z,2)且length(Y) == size(Z,1)。换句话说,X是横向(列方向)坐标,Y是纵向(行方向)坐标,对应关系搞反了,云图就会旋转90度。

先看一个最简示例,感受一下contour的输入数据长什么样:

[X, Y] = meshgrid(-2:0.1:2, -2:0.1:2); Z = X .* exp(-X.^2 - Y.^2); figure; contour(X, Y, Z, 20); colorbar; title('contour example');

这个例子跑通了,就理解了"网格化数据"到底意味着什么:X和Y都是41×41的矩阵,Z也是同规模的矩阵,每个元素对应一个网格点上的物理量。

2.2 等值线层数怎么定

contour(Z)默认画出的等值线数量大概在7到10条,视Z值的取值范围而定,对PIV流场来说往往太稀疏,流场内部的梯度变化根本看不出来。所以一般要显式指定层级数或层级值。

contourf(Xg, Yg, Speed, 30, 'LineStyle', 'none'); % 30层填色等值线

层数太多,云图会显得碎片化,色块之间过渡很碎;层数太少,梯度信息又不够。我的经验值是20到40层对大多数PIV流场都合适,先画一版看效果,再微调。如果你想让特定速度值(比如分离区再附着点附近的速度)在云图里恰好处于色带分界处,就用向量形式手动指定层级:

levels = 0:0.05:0.8; % 从0到0.8,每0.05一层 contourf(Xg, Yg, Speed, levels, 'LineStyle', 'none');

2.3 contourf与contour的选择逻辑

contour画的是纯线条,适合做辅助信息,比如叠在灰度背景上表示涡量等值线。而云图的主角通常是contourf,它把各等值层级用颜色填充,这就是我们常说的"填色云图"。

这里有一个细节必须提醒:contourf默认会在填充色块之间画黑色等值线,层数多时这些黑线会让图面非常杂乱,完全盖过颜色信息。所以绘制云图时,通常要设置'LineStyle', 'none':

figure; contourf(X, Y, Z, 30, 'LineStyle', 'none'); colorbar; title('contourf example');

如果确实需要等值线来辅助读图,推荐在这个基础上用hold on再叠加一次contour,画白色或灰色细线,再用[C, h] = contour(...)拿到等值线矩阵后用clabel标注数值。这样既保留了填色云图的视觉冲击力,又提供了定量读值的便利。

2.4 散点数据直接画contour的常见错误

新手最容易踩的坑,就是拿到PIV的散点数据后,幻想可以"直接画contour"。比如前面那个报错的例子:

% 错误示例 contour(x, y, speed); % speed是列向量,维度不满足要求

MATLAB会直接报错,因为Z必须是二维矩阵。就算运气好,散点数据恰好构成规则网格,如果x、y是向量而speed也是向量,依然不满足维度约束。正确理解是:contour家族的绘图函数,输入的不是"一堆点",而是"一个完整的网格场"。这个网格场怎么从散点来,就是下一部分要重点说的griddata插值。

3. 速度场网格化:从散点到网格,这一步决定云图质量

3.1 为什么必须经过griddata插值

PIV返回的散点数据要变成云图,必须先插值到规则网格上。这一步的质量直接决定云图观感:插值做得糙,云图会出现类似马赛克的色斑;插值过度平滑,又会把小尺度的漩涡结构直接抹平。

MATLAB里散点插值的核心函数是griddata。它的基本调用形式是:

Fq = griddata(x, y, v, Xg, Yg, method);

其中x、y、v是原始散点数据的坐标和物理量,Xg、Yg是目标网格(由meshgrid生成),返回值Fq是与Xg、Yg同尺寸的网格化数组。

3.2 网格间距怎么选

网格生成用meshgrid,但间距怎么定很多人没有概念。我的经验是参考PIV查询窗口的物理尺寸。比如PIV处理时查询窗口是32×32像素,经过标定换算成物理尺度是2 mm,那么网格间距取1 mm到2 mm之间比较合理。

为什么不能随便取?网格太细:一是计算量大,二是插值出来的相邻网格点强相关,云图看起来很均匀,但并没有增加真实信息,反而让人误以为分辨率提高了。网格太粗:很多流场细节直接消失,特别是小尺度的涡结构,可能两个网格之间就跨过去了。

xmin = min(x); xmax = max(x); ymin = min(y); ymax = max(y); dx = 0.002; % 单位m,根据查询窗口物理尺寸换算 [Xg, Yg] = meshgrid(xmin:dx:xmax, ymin:dx:ymax);

3.3 griddata四种插值方法的对比与选择

griddata支持四种插值方法,各有明显的性能和效果差异:

方法计算速度平滑度潜在问题推荐场景
nearest最快最差,色块状梯度信息丢失严重快速预览、数据点极密
linear快适中,线性连续高阶导数不连续默认推荐,大多数PIV数据
cubic较慢平滑可能过冲,超出物理范围数据质量高、分布均匀时
v4最慢,内存消耗大最平滑的全局插值计算量大、同样可能过冲小数据量精细展示

其中linear是默认方法,也是我平时用得最多的。cubic和v4虽然视觉上更光滑,但有一个隐患:插值结果可能超出原始数据的物理范围。比如原始数据里速度最大值明明只有0.5 m/s,cubic插值后可能出现0.6甚至更高的局部极值,这在论文里会被审稿人一眼盯上。

3.4 对u、v分别插值,还是对speed插值?

这里有一个很多人容易忽略的细节:正确做法是先分别对u和v插值,再计算合速度speed。如果直接对sqrt(u.^2 + v.^2)做插值,会丢失方向信息,而且后面叠加速度矢量箭头、计算涡量或流线时,根本拿不到插值后的u、v分量。

% 对速度分量分别插值 Ug = griddata(x, y, u, Xg, Yg, 'linear'); Vg = griddata(x, y, v, Xg, Yg, 'linear'); % 再计算合速度 Speed = sqrt(Ug.^2 + Vg.^2);

这一步非常关键。图省事直接griddata(x, y, speed, ...)的人,画云图的时候看不出差别,但后面只要想叠加矢量,就得重新回来处理。

3.5 插值后NaN空洞的处理策略

PIV数据覆盖区域往往不是规则的矩形,激光照亮区域可能是梯形或圆形,模型壁面处也可能没有有效数据。griddata插值完成后,未覆盖区域会是NaN。这时候有两种处理方式:

第一种是保留NaN,画contourf时这些区域自动留白。好处是诚实展示数据覆盖范围,适合实验报告和学术论文,读者一眼就能看出哪里测到了、哪里没测到。第二种是对内部小空洞做填充,可以用fillmissing(R2016b及以上)或File Exchange上的inpaint_nans工具。填充后云图图面完整,但要注意这些区域本质是估算值,不要在结论里过度依赖。

validMask = ~isnan(Speed); % 保留数据有效区域掩膜 % 如果确实要填充内部空洞 Ug = fillmissing(Ug, 'nearest'); Vg = fillmissing(Vg, 'nearest'); Speed = sqrt(Ug.^2 + Vg.^2);

我的建议是:原始数据覆盖本身就比较完整、只有个别缺口时,小范围填充没问题;但如果数据本身有大面积空白,就保留NaN让图面空白,不要硬填。

3.6 一套完整的云图绘制流程

把前面的步骤串起来,就是一个可以直接拿去用的完整流程:

clearvars; close all; clc; % 1. 读取PIV数据 data = readmatrix('piv_result.csv'); x = data(:,1); y = data(:,2); u = data(:,3); v = data(:,4); % 2. 坏矢量剔除 speed_raw = sqrt(u.^2 + v.^2); mu = median(speed_raw); sig = std(speed_raw); valid = abs(speed_raw - mu) < 3*sig & ~isnan(u) & ~isnan(v); x = x(valid); y = y(valid); u = u(valid); v = v(valid); % 3. 网格化插值 dx = 0.002; % 按查询窗口物理尺寸设定 [Xg, Yg] = meshgrid(min(x):dx:max(x), min(y):dx:max(y)); Ug = griddata(x, y, u, Xg, Yg, 'linear'); Vg = griddata(x, y, v, Xg, Yg, 'linear'); Speed = sqrt(Ug.^2 + Vg.^2); % 4. 绘制云图 figure('Color', 'w', 'Position', [100 100 700 500]); contourf(Xg, Yg, Speed, 30, 'LineStyle', 'none'); hold on; colorbar; colormap(parula); axis equal; caxis([0, max(Speed(:))]); % 新版本建议用 clim([0, max(Speed(:))]) xlabel('x (m)', 'FontSize', 12); ylabel('y (m)', 'FontSize', 12); title('PIV velocity magnitude contour', 'FontSize', 12);

4. 从"能出图"到"能发表":论文级云图美化方案

4.1 配色是云图的第二层信息

MATLAB老版本默认的jet彩虹色,饱和度太高,中间还有一段绿色让人眼花,而且jet在灰度打印时几乎不可读。R2014b之后默认的parula已经不错了,色彩过渡感知相对均匀。如果想让图面更符合现代期刊审美,可以用R2020b引入的turbo,它视觉上和jet类似,但避免了jet在视觉感知上的非线性问题。

colormap(turbo); % 或者 colormap(parula);

还有一个实用的做法:使用viridis这类从深紫到黄绿的感知均匀色图。MATLAB没有内置viridis,但File Exchange上很容易找到,也可以花两分钟自己构建一段颜色渐变矩阵。我的原则是:除非期刊有特殊要求,否则别用jet去"炫彩",审稿人看到彩虹色云图会本能地怀疑作者的专业度。

4.2 colorbar是云图的信息核心

colorbar随手一画不等于画好了。三个细节值得注意:

第一,colorbar要有清晰的变量名和单位,比如 "Velocity magnitude (m/s)" 或者 "u (m/s)"。第二,字体大小要与坐标轴一致,不要图很大colorbar标注却小得看不清。第三,最关键的一点:同一组工况、多时刻云图做对比时,colorbar范围必须统一。如果每张图都自动适配自身的最大值最小值,那0.2 m/s和0.4 m/s在不同图里显示成同一个颜色,读者根本没法对比。

clim([0, 0.8]); % 所有工况统一范围

我一般先画一版所有工况的自动范围,记录下来,选一个有代表性的最大值统一范围,再重新出图。这个习惯在写论文对比图时能省大量返工时间。

4.3 在云图上叠加速度矢量:quiver的正确打开方式

云图展示速度大小分布,但方向信息是缺失的。通常我会在云图上叠加速度矢量,这样既能看到强度分布又能读出流动方向。但直接把所有网格点都画上去,箭头会密成一团黑,跟云图糊在一起。最实用的做法是抽稀:

step = 4; % 每隔4个点取一个矢量 qX = Xg(1:step:end, 1:step:end); qY = Yg(1:step:end, 1:step:end); qU = Ug(1:step:end, 1:step:end); qV = Vg(1:step:end, 1:step:end); hq = quiver(qX, qY, qU, qV, 2, 'k', 'LineWidth', 0.8);

quiver的第五个参数是箭头缩放系数。经验做法:先给2或3,然后看效果调整。太小箭头像蚂蚁,太大箭头相互交叉。一个容易被忽略的问题是数据里的少数异常大值会让自动缩放被极端值主导,其他正常区域的箭头全部缩小到看不见。此时可以先对用于绘图的qU、qV做百分位截断,再去缩放。

矢量颜色和云图的关系也值得想一下。黑色或深灰色矢量在彩色云图上最清晰,白色在暗色区域容易看不清。如果你用的是深色背景色图,可以考虑把矢量设成白色并加一个黑色描边,层次会更清楚。

4.4 叠加壁面和障碍物边界的图层顺序

PIV测区常常有模型壁面、圆柱、翼型等物体轮廓。云图是背景信息,边界是结构信息,缺了边界读者会不知道流场是绕什么物体流动的。

画边界可以用patch填充多边形,或者用plot画轮廓线。重点是图层顺序:

% 先画云图 contourf(Xg, Yg, Speed, 30, 'LineStyle', 'none'); hold on; % 再画壁面边界(压在云图上方) patch(x_wall, y_wall, [0.6 0.6 0.6], ... 'EdgeColor', 'k', 'LineWidth', 1.2); % 最后画速度矢量(保证清晰可读) quiver(qX, qY, qU, qV, 2, 'k');

顺序不能乱。先画矢量再画边界,箭头会把边界挡住很丑;边界盖在云图上刚好。如果你想让矢量箭头也显示在边界之上,就把quiver放在patch之后,但这时要注意箭头不要都扎进壁面区域。

4.5 输出高分辨率图片

论文投稿要求一般至少300 dpi。旧方法用print,新版本我推荐exportgraphics,它处理字体和图面比例更稳:

% 位图输出 exportgraphics(gcf, 'piv_contour.png', 'Resolution', 300); % 矢量输出(投稿推荐) exportgraphics(gcf, 'piv_contour.pdf', 'ContentType', 'vector');

字体大小和窗口宽高比在输出前就要调好,别想着输出后再裁剪。PDF矢量输出时,当前figure窗口的宽高比直接决定图面比例。我习惯在画图前就设好figure('Position', [100 100 700 500])这类参数,而不是画完再拖窗口大小。

5. 踩坑实录:PIV云图绘制中的五个典型问题排查

5.1 坐标尺度混乱导致流场变成"哈哈镜"

我帮课题组处理过一组风洞数据,拿数据的人没说明坐标单位,我看到x峰值1500就默认是mm,直接画图。结果云图的长宽比严重失真,速度矢量方向看起来也是斜的。后来核对原始记录才发现,一部分数据是像素坐标,另一部分是物理坐标,混在一起画当然出问题。

排查方法:读入数据第一步就统一转成国际单位,在脚本开头加注释块标注每列数据的物理含义和单位。如果PIV系统有标定文件,优先用标定后的物理坐标。没有标定信息时,先画一版散点图看坐标范围,确认是毫米还是像素量级,再做换算。

5.2 NaN区域让云图出现莫名其妙的"白洞"

如果数据覆盖区域不规则,contourf画出来后边界会出现一些尖角或白斑。这时候先用validMask = ~isnan(Speed)看看,区分"本来就没有数据的位置"和"应该有数据但没测到的壁面区域"。

如果是壁面或无数据区,我的做法是用patch画一个灰色覆盖层,把云图的留白区域与真实壁面区分开。不要用纯白色填充,因为云图背景默认就是白色,读者看不出那里是壁面还是无数据区。灰色或深灰色覆盖层配合图例说明,图面信息更清楚。

5.3 云图坐标方向与相机图像方向不一致

PIV软件导出的坐标往往是图像像素坐标,y轴方向是向下的;而物理坐标我们习惯y轴向上。如果直接画,云图会上下颠倒。这个坑在PIV数据里特别典型,尤其是从PIVlab这类基于图像处理软件导出的数据,图像坐标系和笛卡尔坐标系经常被搞混。

解决办法:绘图前加一句axis xy,或者自己在数据导入时做一次y方向的翻转。判断依据很简单:如果流场里有重力方向(比如自然对流、落水实验),先确认重力在图上应该指向哪个方向,再决定是否翻转。这比你事后发现图反了再重画高效得多。

5.4 quiver箭头要么密成一团,要么稀疏得可怜

箭头密度问题主要出在抽稀间隔和缩放系数上。另一个隐藏因素是速度数据里有少数异常大值,导致quiver自动缩放被极端值主导。排查步骤:

首先看max(abs(qU(:)))和max(abs(qV(:)))是否明显大于数据主体范围。如果是,就对用于绘图的矢量做一次截断,比如把超过99.5%分位的值替换为分位值。然后再调整缩放系数。这样箭头长度分布会合理很多,流场方向信息也能正确传达。

5.5 colorbar范围被自动缩放,多工况对比图放到一起露馅

前面强调过统一clim,但实际操作中还有一个麻烦:如果max(Speed(:))在不同工况间差异很大,统一clim会导致某些图整体颜色偏深或偏浅,看起来对比度不足。我的处理方式是:先把所有工况的自动范围都打印出来,看整体分布,选一个能代表主流量级的统一范围,个别差异极大的工况单独调整后再统一。

画完第一版对比图,一定要把图放在一起看,确认同一个速度值在不同的图里颜色一致。这一步做完,云图对比才算真正科学。

5.6 插值方法带来的"过冲"数据

用cubic或v4方法插值后,速度云图可能出现超过原始数据范围的斑点。比如原始速度最高0.5 m/s,插值后云图上却出现0.65 m/s的红色区域。这不是发现了新物理现象,纯粹是插值算法的数值振荡。

处理方式三个:

  • 直接换成linear插值,最保险;
  • 对插值结果做物理范围截断,比如Ug(Ug < 0) = 0; Ug(Ug > umax) = umax;,但截断会产生色块突变,不彻底;
  • 如果必须用cubic插值出图,检查云图极值是否合理,并在图注里说明插值方法。

特别是当你下一步要计算涡量、散度这类涉及空间导数的物理量时,cubic插值引入的振荡会被求导过程进一步放大,出现更离谱的数值。这时候宁可牺牲一点视觉平滑度,也要保证物理量的可信度。


最后再分享一个小技巧:每画完一张云图,用datacursormode在云图上点几个位置,和原始散点数据里的值做交叉核对。云图是插值出来的可视化结果,不是原始数据本身。图面再漂亮,如果关键位置的速度值和实测对不上,那这张图就只是一张好看的壁纸。以我个人经验,PIV云图画到能发表的水平,70%的功夫其实不在画图参数本身,而在预处理和插值环节——数据干净了,contourf那几行代码根本不会出岔子。

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询