如果你手头拿到一套 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的单位都不一样了。如果代码里写死变量名，换一个文件就会直接报错。所以更稳妥的做法是写一个通用读取函数，通过候选变量名依次判断文件里存在哪个变量。

下表列出了这套产品里常见变量及其含义，方便你对照检查：

2.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或老版本的caxis：

clim([-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 end

GIF 需要注意的是颜色数量限制。由于 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调小。我在实际作图时，一般是先画一帧静态图，肉眼判断箭头密度是否合适，再决定最终步长。这样做出来的动图，既有信息量，又不会因为箭头太密而失去可读性。