如果你手头拿到一套 AVISO 这类卫星高度计网格化产品想在 MATLAB 里画一张“SLA u”的动图那这篇文章应该能帮你跳过不少弯路。这里的 SLA 指海表面高度异常u 是地转流的纬向分量。我按“读取 NetCDF → 掩膜与清洗 → 由 SLA 算地转流 u → 单帧绘图 → GIF/MP4 动图”的顺序把整条链路完整走了一遍并把每个环节里最容易踩的坑单独拎出来讲。整个流程是通用的目标海域换成任何区域都适用无论你是做物理海洋研究、渔业资源分析还是纯粹想用卫星高度计数据做一个可视化作品都可以直接参考这份流程。需要先说明的是我不打算只丢几个代码片段。因为这种数据处理任务绝大多数问题不是出在最后那几行画图而是出在数据读取、变量名匹配、掩膜处理、方向符号这些看起来琐碎的地方。所以下面每一节我都会把“为什么这么做”讲清楚再把可以直接复制的 MATLAB 代码给出来。1. 先搞清楚“SLA u”这三个术语到底代表什么1.1 这套高度计产品能直接拿到哪些变量很多第一次接触 AVISO 高度计数据的同学打开 NetCDF 文件后会看到十几个变量容易一下子看晕。除了最常见的sla还有adt、mdt、ugos、vgos、err以及经纬度、时间、质量标记等。这里 SLA 指的是海表面高度异常是相对于多年平均海面或平均动力地形的高度偏差单位一般是米。它直接反映了涡旋、流轴、锋面等中尺度信号的空间形态。ugos和vgos则是地转流的纬向分量和经向分量单位一般是米每秒。如果文件里直接提供这两个变量就不用再从 SLA 梯度去推流速了直接读取即可。需要提醒的是不同版本的产品速度变量的命名习惯不一样。有的用ugos、vgos有的直接叫u、v甚至还有u_geo、v_geo这种自定义命名。所以写代码前第一件事不是急着读数据而是先用ncdisp看一眼结构。1.2 沿轨数据和网格化数据怎么选如果你的目标是做动图我强烈建议直接选网格化产品。卫星高度计的原始沿轨数据虽然在沿轨方向分辨率很高但需要自己做轨道插值、跨轨道的空间合成、不同卫星之间的偏差校正处理成本非常高不适合快速出图。网格化产品已经把多颗卫星的高度计观测融合成了规则经纬度网格空间分辨率常见的有 0.25° 和 0.125°时间分辨率则有日平均、周平均、月平均。对动图而言日平均通常是最合适的因为周平均和月平均会把中尺度信号的瞬变过程抹掉不少。如果你的目标是看季节变化或年际变化周平均、月平均反而更平滑渲染出来的动画也更稳定看你的具体需求。1.3 什么时候该自己由 SLA 算 u什么时候直接读 u这里有一个很容易纠结的问题既然产品里可能直接给u为什么还有很多人选择自己算原因主要有两类。一类是你的数据源只提供了 SLA没有速度变量另一类是你希望用一套完全统一的计算流程来处理多套数据避免不同版本速度变量的处理差异。但自己算也有条件限制。地转近似在赤道附近是不成立的因为科里奥利力趋近于零用 SLA 梯度算出来的流速会大得离谱所以通常需要在纬度 ±5° 以内做遮挡或插值。如果产品里直接有ugos、vgos而且你的研究区域不在高纬那直接读取是最省事的方案。但如果文件只给了 SLA或者你想自己控制每一步处理细节那么从 SLA 梯度推导 u/v 也完全可行下面会给出完整代码。2. 数据预处理是决定动图质量的地基2.1 文件结构检查与变量名差异拿到 NetCDF 文件后先运行一次ncdisp把文件里的变量名、维度、单位看清楚。这一步能避免后面很多“变变量名不存在”的报错。dataFile sla_u_sample.nc; ncdisp(dataFile);这是我实际遇到过的情况同一套高度计网格化产品2010 年的文件里速度变量叫ugos、vgos换到 2018 年的文件后变量名变成了u、v连time的单位都不一样了。如果代码里写死变量名换一个文件就会直接报错。所以更稳妥的做法是写一个通用读取函数通过候选变量名依次判断文件里存在哪个变量。下表列出了这套产品里常见变量及其含义方便你对照检查变量名含义典型单位sla海表面高度异常madt绝对动力地形mmdt平均动力地形mugos地转流纬向分量m/svgos地转流经向分量m/serr数据误差估计mtime时间轴hours/days since 1950-01-012.2 陆地掩膜和无效值不处理会画出一片假海洋这是最容易踩的坑。很多网格化产品在陆地区域的 SLA 值不是 NaN而是 0如果不做掩膜画出来的图会出现一块从海岸线一直延伸到内陆的红色或蓝色区域而且在海陆交界处产生非常大的伪梯度。用这个伪梯度去算地转流岸线附近会出现一条速度极大的箭头带直接毁掉整张图。处理方式有三种一是优先读取文件自带的陆地掩膜变量例如landmask二是如果文件没有掩膜变量就根据 SLA 第一帧的缺失值来构造掩膜三是如果 SLA 在陆地上是 0则需要准备一份独立的陆地掩膜数据把所有陆地格点设为 NaN。代码实现如下sla ncread(dataFile, sla); % 维度可能是 (lon, lat, time) lon ncread(dataFile, longitude); lat ncread(dataFile, latitude); % 如果文件里有 landmask 变量优先使用 try landmask ncread(dataFile, landmask); landmask squeeze(landmask); catch landmask ~isnan(sla(:,:,1)); end % 将陆地或缺失值统一设为 NaN sla(repmat(~landmask, 1, 1, size(sla,3))) NaN; % 剔除极端异常值 sla(abs(sla) 5) NaN;这里要注意repmat在数据量很大时非常耗内存。如果文件是全球 0.25° 日平均直接把 sla 的所有时间帧都加载进来再乘一个逻辑掩膜矩阵内存很容易爆掉。更好的做法是在后面的循环里逐帧处理而不是一次性把三维数组全部改造完。2.3 空间平滑和异常值剔除高度计网格化产品本身已经做过不少处理理论上是比较干净的。但如果你发现单帧图上有很多孤立的椒盐噪声点那就有必要做一次轻度空间平滑。MATLAB 里用imgaussfilt或者medfilt2都可以我自己更常用imgaussfilt因为它对边界更友好。比较稳妥的处理方式是逐帧进行高斯平滑slaTmp squeeze(sla(:,:,t)); slaTmp imgaussfilt(slaTmp, 1.2, Padding, replicate); sla(:,:,t) slaTmp;不过要记住空间平滑是把双刃剑。sigma 取值过大中尺度涡旋的边界会被磨得过于平滑真实的小尺度信号就被抹掉了。一般 sigma 取 1 到 1.5 个格距就够如果你想保留更多细节甚至可以不做空间平滑只做时间方向的滑动平均。时间方向的平滑也很常见。如果是逐日动画但原始日数据里含有较多轨道采样噪声可以先沿时间方向做 3 天或 5 天滑动平均。MATLAB 里直接用movmean就可以例如slaSm movmean(sla, 5, 3)这里的第三维就是时间维。这样得到的动画会稳定很多代价是时间上的一些瞬时变化会被减弱。3. 从 SLA 到地转流 u公式、单位与边界限制3.1 地转近似下 SLA 与 u 的物理关系地转流是物理海洋里最基础的概念在大尺度运动中水平压力梯度力与科里奥利力基本平衡于是海流沿着等压线流动。对海表面高度异常 SLA 来说水平压力梯度主要体现在海面斜率上所以地转流的纬向分量 u 和经向分量 v 可以由 SLA 的空间梯度直接计算[ u -\frac{g}{f} \frac{\partial SLA}{\partial y}, \quad v \frac{g}{f} \frac{\partial SLA}{\partial x} ]其中 (g \approx 9.81\ \mathrm{m/s^2})科里奥利参数 (f 2\omega \sin(\phi))地球自转角速度 (\omega 7.2921 \times 10^{-5}\ \mathrm{rad/s})(\phi) 是纬度。理解这个公式的最好方式是用天气图类比。天气图上气流沿等压线流动北半球低压逆时针旋转、高压顺时针旋转。海洋里的地转流也一样SLA 的等值线就相当于海面的“等压线”SLA 梯度越大的地方流速越快。所以当你在 SLA 图上看到密集的等值线时不用算也知道那一带流很强。下面给一个具体的量级估算。假设 SLA 在纬度方向 25 km 内变化了 0.1 m纬度 30° 处的科里奥利参数约为 (3.65 \times 10^{-5}\ \mathrm{s^{-1}})。那么[ \frac{\partial SLA}{\partial y} \frac{0.1}{25000} 4 \times 10^{-6} ][ u \approx -\frac{9.81}{3.65 \times 10^{-5}} \times 4 \times 10^{-6} \approx -1.07\ \mathrm{m/s} ]这个量级在中尺度涡旋里是合理的说明公式的系数没有明显量级问题。3.2 中央差分实现与纬度限制在 MATLAB 里gradient函数用起来很方便但它对坐标轴的顺序和方向很敏感而且如果经度方向的格距随纬度变化gradient会让dx的处理变得不够精确。所以我更建议直接写一个手动中央差分的函数把经纬度距离换算清楚。假设你已经在某个时间帧上拿到了二维 SLA 矩阵sla2d并且它的维度是(lat, lon)经度向量和纬度向量都是递增的。那么可以这样写function [u, v] geo_velocity_from_sla(sla2d, lon, lat, g, omega) nlat length(lat); nlon length(lon); dlat_deg abs(lat(2) - lat(1)); dlon_deg abs(lon(2) - lon(1)); dy dlat_deg * 111320.0; % 纬度方向距离单位m dSLA_dy zeros(size(sla2d)); dSLA_dx zeros(size(sla2d)); % 纬度方向中央差分 dSLA_dy(2:end-1, :) (sla2d(3:end, :) - sla2d(1:end-2, :)) / (2 * dy); % 经度方向中央差分注意经度方向格距随纬度变化 for i 2:nlat-1 dx dlon_deg * 111320.0 * cosd(lat(i)); dSLA_dx(i, 2:end-1) (sla2d(i, 3:end) - sla2d(i, 1:end-2)) / (2 * dx); end f 2 * omega * sind(lat); fMat repmat(f(:), 1, nlon); u -g ./ fMat .* dSLA_dy; v g ./ fMat .* dSLA_dx; % 赤道附近地转近似失效 u(abs(lat) 5, :) NaN; v(abs(lat) 5, :) NaN; % 边界处梯度不可靠直接置 NaN u(:, 1) NaN; u(:, end) NaN; v(1, :) NaN; v(end, :) NaN; end这个函数的代码量虽然比一行gradient多一点但每一步的物理含义都很清楚也便于自己检查和调试。需要注意的是边界处理。水平边界的 SLA 梯度是没法用中央差分算的所以函数直接把这些位置设成了 NaN。如果你后面画图时发现边界处有空白这不是 bug而是特意为之。海陆边界附近的梯度本来就不适合用于地转流计算尤其是岸线附近SLA 插值误差很大箭头方向很容易“放飞自我”。如果目标区域包含赤道那么 ±5° 范围内会出现一条空白带这是地转近似本身的限制。你可以在图上明确标注“赤道附近地转关系不适用”别让读者误以为数据缺失。4. 单帧静态图先画好一帧再谈动画4.1 填色图与矢量图的排版要点做动画之前一定要先把单帧静态图调试好。因为动画里的每一帧其实就是一帧静态图如果单帧里的坐标、颜色、箭头都没有对齐做出来的动图只会更乱。首先要做的是区域切片。如果文件是全球网格建议先把目标区域的经纬度范围截取出来再进行后续计算。比如某研究区域非常粗糙地设定为经度 105°E 到 130°E、纬度 10°N 到 28°N那么在 MATLAB 里可以这样切片lonLim [105 130]; latLim [10 28]; idxLon lon lonLim(1) lon lonLim(2); idxLat lat latLim(1) lat latLim(2); lonS lon(idxLon); latS lat(idxLat); slaS squeeze(sla(idxLon, idxLat, t)); % 维度变成 (lon, lat) uS squeeze(u(idxLon, idxLat, t)); vS squeeze(v(idxLon, idxLat, t));之后画图时SLA 矩阵需要从(lon, lat)转置成(lat, lon)否则pcolor的横纵坐标会对不上。这也是一个很典型的坑画出来的图可能整体错位海岸线怎么都对不上最后发现是维度顺序错了。单帧画图的代码大致如下figure(Color,white,Position,[100 100 1000 500]); pcolor(lonS, latS, slaS); shading interp; hold on; % 箭头降采样避免满屏箭头 qstep 3; quiver(lonS(1:qstep:end), latS(1:qstep:end), ... uS(1:qstep:end, 1:qstep:end), ... vS(1:qstep:end, 1:qstep:end), ... 1.5, k, LineWidth, 0.8); load coastlines; plot(coastlon, coastlat, k); clim([-0.3 0.3]); colorbar; xlim(lonLim); ylim(latLim); axis equal;这里需要留意quiver的输入。uS和vS如果维度是(lon, lat)那么传输给quiver前同样需要转置。并且由于箭头密度取决于矩阵尺寸降采样步长qstep要按实际网格大小调整。如果网格很细不降采样的话箭头会叠在一起视觉上完全糊掉。4.2 海岸线、色标和颜色轴的固定策略海岸线叠加看似简单但却是验证数据坐标是否正确的关键步骤。如果你画出来的 SLA 色块明显不在海里而是整体偏移到陆地上说明在读取或转置时维度搞错了。MATLAB 基础版提供了一份默认的海岸线数据直接load coastlines就能用然后用plot叠加相当方便。色标选择上我习惯使用蓝白红三段式发散色标让零值附近显示为浅色正异常偏红负异常偏蓝。MATLAB 自带的是像parula或jet这种连续色带但 SLA 图用发散色标更符合海洋学惯例。你可以用下面的方式构造一个简单的蓝白红色标cmap [ ... linspace(0, 1, 128), linspace(0.2, 1, 128), linspace(1, 1, 128); ... linspace(1, 1, 128), linspace(1, 0.2, 128), linspace(1, 0, 128)]; colormap(cmap);不过这个手工色标只是一个简化版本真正好看的蓝白红还需要在后面调一调。重点是静态图里颜色轴可以自动变化但动画图里一定要固定。因为动画是逐帧播放的如果每一帧都自动调整色标范围观众会看到颜色不停“呼吸”几秒钟后就会产生视觉疲劳而且会误读信号强度。固定颜色轴的方法就是用clim或老版本的caxisclim([-0.3 0.3]);这个范围不是随便定的。建议先统计所有时段 SLA 的 1% 和 99% 分位数再根据分位数取一个对称范围。不要用全局最大值否则大多数帧的颜色都会被压缩到很浅的范围内显著信号反而不突出。5. 动图制作GIF 和 MP4 两条路的完整代码5.1 基于 imwrite 的 GIF 循环动图最常见的形式是 GIF。MATLAB 里生成 GIF 的标准思路是每一帧用getframe抓取图像再把图像数据转换成索引图像最后用imwrite写出。下面是一个完整的循环gifFile sla_u_anim.gif; firstFrame true; figure(Color,white,Position,[100 100 1000 500]); for t 1:ntime % 更新当前帧数据 slaS squeeze(sla(idxLon, idxLat, t)); uS squeeze(u(idxLon, idxLat, t)); vS squeeze(v(idxLon, idxLat, t)); % 重新绘制当前帧 cla; pcolor(lonS, latS, slaS); shading interp; hold on; quiver(...); % 和单帧图相同 load coastlines; plot(coastlon, coastlat, k); clim([-0.3 0.3]); title(datestr(timeVec(t))); % 抓取并写入 frame getframe(gcf); [A, map] rgb2ind(frame2im(frame), 256); if firstFrame imwrite(A, map, gifFile, gif, LoopCount, Inf, DelayTime, 0.2); firstFrame false; else imwrite(A, map, gifFile, gif, WriteMode, append, DelayTime, 0.2); end endGIF 需要注意的是颜色数量限制。由于 GIF 最多支持 256 种颜色连续色标在压缩成索引图像时很可能出现色带尤其当同时显示 SLA 填色图、黑色箭头、海岸线、色标时颜色条会非常紧张。如果你只做汇报用的简短动画GIF 是够用的但如果要做长时间序列动画GIF 文件会变得非常大而且播放卡顿。5.2 基于 VideoWriter 的高清 MP4 写法如果希望动画更清晰、更流畅MP4 是更好的选择。MATLAB 的VideoWriter可以直接输出 MP4 格式代码量比 GIF 还少videoFile sla_u_anim.mp4; v VideoWriter(videoFile, MPEG-4); v.FrameRate 5; open(v); figure(Color,white,Position,[100 100 1000 500]); for t 1:ntime % 更新与绘制当前帧和上面一样 cla; pcolor(lonS, latS, slaS); shading interp; hold on; quiver(...); load coastlines; plot(coastlon, coastlat, k); clim([-0.3 0.3]); title(datestr(timeVec(t))); writeVideo(v, getframe(gcf)); end close(v);这里如果手动指定FrameRate 5表示每秒播放 5 帧。假设你的数据是日平均那么 5 帧每秒意味着实际时间轴上每秒钟播放 5 天视觉节奏比较合适。如果你的数据量很大还可以进一步提高帧率比如设成 10但帧率太高反而会造成“信息过密”的感觉。VideoWriter 在旧版 MATLAB 里可能不支持 MPEG-4只支持 Motion JPEG AVI。如果遇到这种情况可以把文件名后缀改成.avi编码改成Motion JPEG AVI。不建议用未压缩的 RGB AVI文件体积会爆得离谱。5.3 帧间平滑与时间游标即使数据本身是逐日的原始动画播起来也可能非常“跳”。这是因为高度计融合产品里仍有不少短时间的噪声SLA 的逐日变化未必是平滑的。在做动画之前我通常建议先对 SLA 和 u/v 做时间方向的 3 到 5 天滑动平均这样动图的每一帧之间过渡会更自然观众更容易抓住中尺度信号的演变过程。时间方向平滑的代码非常简单slaSm movmean(sla, 5, 3);但要注意如果你用了movmean那么在绘图时所有帧的时间标签也要相应调整最好标注滑动窗口的中心日期而不是窗口的第一天或最后一天。否则观众看到的时间轴和实际信号的中心时间会产生偏差。动画里的时间游标也很关键。在每次更新标题时建议显示具体的时间点例如datestr(timeVec(t))并把它放在图形右上角或标题处。因为动图一旦循环播放如果没有时间标签观众根本不知道每一帧代表哪一天。时间轴的制作需要先把原始时间单位解析出来。常见格式是days since 1950-01-01你可以用 MATLAB 的datetime直接转换baseTime datetime(1950, 1, 1); timeVec baseTime days(time);如果单位是秒钟则改成baseTime seconds(time)。这一步很小但能避免很多动画“时间轴乱跳”的错觉。6. 我踩过的坑和最终效果优化6.1 变量名版本坑同一产品不同年份数据结构不同最开始我处理数据时直接写死了变量名ugos和vgos。结果从 2018 年之后的产品文件开始速度变量名变成了u和v一运行就是“变量不存在”的报错。后来我再也不敢偷懒每次处理前都会先ncdisp并把变量名做成可配置的候选列表。可以写一个简单的读取函数自动从文件里找候选变量名function varData ncread_alt(dataFile, candidateNames) for i 1:length(candidateNames) if ~isempty(ncdisp(dataFile, candidateNames{i})) varData ncread(dataFile, candidateNames{i}); return; end end error(No candidate variable found: %s, strjoin(candidateNames, , )); end这个函数虽然简单但可以大幅减少换数据源时改代码的次数。ncdisp判断变量是否存在时你还可以用ncinfo来检查变量列表效果一样。6.2 纬度递减时梯度方向搞反另一个让我印象深刻的坑是纬度数组排列方向。正常情况下纬度应该是从南到北递增也就是低纬度在前、高纬度在后。但有的数据文件偏偏把纬度排成了从北到南递减导致 SLA 矩阵的纬度方向反了。如果此时直接用中央差分公式算 u/v画出来的箭头方向会和真实流向完全相反。排查这个问题其实不难。先打印出lat(1)和lat(end)确认方向。如果确实递减要把纬度翻转到递增同时把 SLA 矩阵沿纬度维翻转。假设lat对应sla2d的第一维行那就执行lat flipud(lat); sla2d flipud(sla2d);关键问题是翻转到底作用在哪一维取决于你自己定义的矩阵维度。建议在自己封装的计算函数里加一个断言如果lat(1) lat(end)就自动翻转。这样就不会再因为这个问题画反了。6.3 自动色标导致画面闪跳与误读我第一次做动图时每一帧都没有固定色标范围结果播放时颜色不停闪跳一个本来很强的涡旋在颜色变化上反而不明显。原因很简单自动色标会让每一帧都按照当前帧的最大值来分配颜色导致同一数值在不同帧里显示完全不同的颜色。后来我养成两个习惯。第一在开始做动画前把所有帧的 SLA 按百分位数统计一遍确定一个固定的clim范围第二在动画每帧的标题或色标旁注明单位是多少让读者不用猜。如果你发现有少数极端值把颜色轴拉得太宽可以对这些极端帧不做饱和处理而是直接固定颜色轴让超出范围的值在图中表现为饱和红或饱和蓝这种做法是合理且有意的。6.4 内存不足与批量循环的优化思路全球 0.25° 日平均数据一天大概是 1440×720 个格点一年就是 1440×720×365 个单精度数换算下来约 1.5 GB。如果你一次性读取五年甚至十年的数据内存直接爆炸。我早期犯过这个错误读数据时兴致勃勃等 MATLAB 解算完所有文件内存已经红了绘图时卡到几乎无法操作。正确做法是逐帧读取。用ncread的start和count参数每次只读取一个时间片。例如% 假设索引数组 idxLon、idxLat 已经确定 startIdx [idxLon(1), idxLat(1), t]; countIdx [length(idxLon), length(idxLat), 1]; slaT ncread(dataFile, sla, startIdx, countIdx);这样每次循环只读一帧目标区域的数据内存占用可以控制在几十 MB 以内。如果还嫌慢可以先用ncread读取全部经纬度、时间再在循环里逐帧读取变量。把计算和绘图分离开也是一个更好的模式先做一次数据处理循环把目标区域的结果保存成一个小的.mat文件后续动图绘制直接加载.mat会轻松很多。最后再分享一个小技巧箭头降采样步长不一定固定。目标区域网格更密时qstep可以调大网格较粗时qstep调小。我在实际作图时一般是先画一帧静态图肉眼判断箭头密度是否合适再决定最终步长。这样做出来的动图既有信息量又不会因为箭头太密而失去可读性。