企业官网建设流程全解析

简介：面向FVCOM海洋模型使用者的MATLAB源码集，专注解决模型前处理与后处理环节的常见问题，适用于风暴潮模拟、潮汐预报、近岸水交换等研究场景。前处理部分覆盖非结构化三角网格生成、边界强迫条件导入、初始温度盐度设定，以及卫星遥感、浮标等观测数据的插值整合，最终生成模型可读的.nc或.dat输入文件。后处理部分支持读取大体积二进制或NetCDF输出结果，开展统计分析、时空插值，绘制平面图、断面图、时间序列图与流场动画，并可进行敏感性对比分析，便于验证模型表现。压缩包共7个文件，全部为.m脚本，仅23KB，轻量易修改，已有479人学习下载。脚本内含坐标转换、节点名称管理、流速裁剪、快速一维插值、断面时间序列绘图、平均态计算与极坐标标注等实用函数，用户可依据研究区域或网格分辨率快速调整参数，也可提取核心函数嵌入自己的流程中，脚本命名与函数功能对应，便于检索调用，适合海洋科学、环境工程等领域的研究生与工程师快速上手。

1. FVCOM 模拟最容易翻车的环节，其实是数据进出的那几十个脚本

做 FVCOM 的人大多有这样的经历：求解器跑了一个月，后处理又花了两周。真正让人头疼的不是模型本身，而是前处理里网格和观测数据对不上、后处理里二进制输出读不出来这类琐碎问题。这套 fvcom-matlab.rar 工具包就是在解决这件事——它不碰 FVCOM 内核，专门负责模型跑之前的数据准备和跑完之后的场提取。压缩包里的 arrow.m、node_name.m、trim_uv.m、interp1_fast.m、plot_zsec_uts.m、plot_usmean_onz.m、suptitle.m 这些源程序，覆盖了从节点文件生成、观测数据插值到断面图绘制、多子图拼接的完整链路。适合正在跑 FVCOM 的研究生和工程师，也适合刚开始接触非结构化网格海洋模型、想找一套现成 MATLAB 后处理流程来改的人。下面我把每个脚本的用途、参数含义和接线方式拆开讲。

2. 前处理不是画网格，是把观测数据“翻译”成 FVCOM 能读的场

2.1 这套工具包里的前处理脚本到底做了什么

FVCOM 用的是非结构化三角网格，节点坐标、单元连接关系、地形水深这些信息都写在 NetCDF 或 ASCII 文件里。前处理的核心任务不是“画一张好看的网格图”，而是把杂乱的外部数据映射到网格节点上。压缩包里的 node_name.m 负责把网格节点编号和经纬度、水深整理成可读文件，interp1_fast.m 负责把站点观测或再分析数据快速插值到 FVCOM 节点，trim_uv.m 则是把插值后的流速场按边界或深度范围做裁剪。它们互相配合，形成一条从原始数据到模型输入文件的流水线。

用之前先看清脚本的输入约定。这些 MATLAB 源程序没有统一的数据结构，我一般先把 FVCOM 输出文件读成结构体，再传给各脚本：

% 读取 FVCOM 网格基本信息 nc = netcdf('fvcom_grid.nc', 'nowrite'); x = nc{'x'}(:); % 节点 x 坐标（投影后） y = nc{'y'}(:); % 节点 y 坐标 h = nc{'h'}(:); % 水深，负值表示水下 tri = nc{'nv'}(:)'; % 单元节点连接关系，3 x nele

逻辑说明：FVCOM 的 NetCDF 输出中，x和y是节点坐标，h是静态水深，nv是每个三角形单元的三个节点编号。用 MATLAB 的netcdf工具箱读出来以后，tri必须转成3 x nele的格式，因为后续做插值和绘图时，patch或trisurf都按这个约定接收。这里nc{'nv'}(:)读出来的是nele x 3，转置一下最稳妥。很多脚本报“维度不匹配”，根源就是这一步方向反了。

2.2 interp1_fast.m 的插值逻辑与参数调整

interp1_fast.m是前处理里最常用的一个脚本。名字看起来是 MATLAB 内置interp1的加速版，实际做的是“站点数据到网格节点”的空间插值。它的输入一般是四个参数：源数据点的经纬度、源数据值、目标节点的经纬度，以及插值半径或搜索点数。下面是一段典型调用：

% 假设 obs_lon, obs_lat 是观测站点坐标，obs_val 是温度观测值 % target_lon, target_lat 是 FVCOM 网格节点坐标 field = interp1_fast(obs_lon, obs_lat, obs_val, target_lon, target_lat, ... 'radius', 0.05, 'min_points', 3);

逻辑说明：这段代码把每个目标节点周围指定半径内的观测点找出来做反距离加权平均。radius的单位要和经纬度一致，如果是球面坐标，0.05 度大约 5 公里；min_points是最少参与插值的点数，少于它就把该节点置为NaN，避免外延出离谱值。如果观测站点稀疏，我会把radius放大到 0.1 度，但要在后续检查isnan(field)的比例，超过 20% 就说明网格分辨率比数据密度高太多，插值出来的场其实没有意义。

这套插值逻辑比 MATLAB 自带的scatteredInterpolant好在两点：一是支持半径限制，不会用几十公里外的站点外推一个值来“填空”；二是反距离权重天然适合海洋数据，因为近岸温度、盐度局部相关性很强。代价是慢，节点多时建议先压缩目标点，只对浅水区或关注区域插值。

2.3 把海岸线和地形数据变成网格属性的常见做法

前处理里还有一个高频需求：把岸线数据变成网格的干湿边界，把地形数据修正到网格水深上。这个工具包里没有专门的地形修饰脚本，我一般用 MATLAB 的inpolygon加interp1_fast配合实现：

% land_mask = inpolygon(x, y, coastline_lon, coastline_lat); % 把岸线外的节点标记为陆地 x_land = coastline_lon; y_land = coastline_lat; in = inpolygon(x, y, x_land, y_land); h(in) = -999; % 标记为干节点 % 用地形数据替换网格水深 h_interp = interp1_fast(bathy_lon, bathy_lat, bathy_depth, x, y, ... 'radius', 0.02, 'min_points', 2); id = ~isnan(h_interp) & ~in; h(id) = h_interp(id);

逻辑说明：inpolygon返回每个节点是否落在岸线多边形内部，陆地上的节点水深被设为-999，FVCOM 会把它当作干单元不参与计算。地形数据插值用的是同一套interp1_fast，但半径取得更小，因为高分辨率地形数据本身很密，半径太大反而会把峡谷填平。这段代码的关键在掩膜顺序：先判陆地，再插值地形，最后用id把地形值只写到非陆地节点上，避免岸线外的插值假值污染网格。

3. 后处理读取：从 NetCDF 输出到断面图的关键环节

3.1 FVCOM 输出的时间维和节点/单元数据

FVCOM 的 NetCDF 输出分节点变量和单元变量两类。节点变量包括水位zeta、温度temp、盐度salt，单元变量包括u、v流速分量。两者网格不同：节点数据存在xc/yc上，单元数据存在x/y上，这是初学者最容易踩的坑。plot_zsec_uts.m这类脚本内部会分别处理这两套坐标，读数据前先确认变量名valid_min和valid_max，可以避免把单元坐标当节点坐标画图的错位问题。

读取 FVCOM 输出的标准起点是这样一段代码：

nc = netcdf('output_0001.nc', 'nowrite'); time = nc{'time'}(:); u = nc{'u'}(:); % 维度：nele, siglay, time v = nc{'v'}(:); temp = nc{'temp'}(:); % 维度：node, siglay, time siglay = nc{'siglay'}(:); % 无量纲 sigma 层，-1 到 0 x_ele = nc{'x'}(:); y_ele = nc{'y'}(:); x_node = nc{'xc'}(:); y_node = nc{'yc'}(:);

逻辑说明：siglay是负值，-1表示底层，0表示表层。后处理计算垂向平均时，不能简单把第一层和最后一层相加，必须乘以层厚度权重。默认的 sigma 层不均匀，表层加密，所以平均时我会先把siglay变换成深度。这个工具包里的plot_usmean_onz.m就是做垂向平均流的，它的输入参数里有一个zlim选项，可以只对某一层到某一层之间做平均。

3.2 用 trim_uv.m 切出边界层流速

trim_uv.m名字里的trim不是修剪数组，而是把流速矢量限制在“有效计算域”内。FVCOM 的边缘单元可能因为插值产生异常值，或者在干湿边界附近出现非物理流速，这个脚本用掩膜把无效点剔除。典型调用方式如下：

% u,v: 单元流速，mask: 逻辑数组，1 表示有效单元 [u_clean, v_clean] = trim_uv(u, v, mask); % mask 通常由水深和干湿状态生成 mask = h_ele > -10; % 只保留水深大于 10m 的区域

逻辑说明：mask可以是逻辑数组，也可以是一个数值阈值，如果传单个数，脚本会把它当成流速绝对值上限，超过阈值的点被置为NaN。实际中我建议同时用两种：物理上通过水深排除陆地和浅滩，数值上排除超过 5 m/s 的异常值。处理完之后再用plot_usmean_onz.m画垂向平均流，图面就不会出现海底地形剧烈变化处的零散箭头。

3.3 时间序列与断面绘图的实现：plot_zsec_uts.m

plot_zsec_uts.m是画“沿某条断面的时间序列图”的脚本，也就是常说的 Hovmöller 图。输入是一条断面端点的经纬度、某个变量的三维场（节点或单元）、以及时间数组。脚本先把断面上的点找出来，再沿断面插值，最后画成 x 轴为距离、y 轴为时间、颜色为变量值的等值线图。

% 断面从 A 点到 B 点，共 50 个插值点 lonA = 122.5; latA = 30.5; lonB = 123.0; latB = 31.0; npts = 50; plot_zsec_uts(u, x_ele, y_ele, siglay, time, lonA, latA, lonB, latB, npts);

逻辑说明：plot_zsec_uts内部会先算出每个单元到断面的投影距离，选取距离小于某个阈值的单元，然后在断面方向做插值。这里u是原始流速场，siglay用来计算垂向位置。如果断面经过地形剧烈变化区域，默认阈值会选不到足够单元，此时脚本会报“Not enough points”，需要把npts调小或者检查断面经纬度是否越界。我一般会在调用前加一句disp(minmax(distance_to_section))来确认断面的实际覆盖范围。

4. 从 MATLAB 环境搭建到跑通第一个后处理脚本

4.1 MATLAB 环境与路径配置

拿到fvcom-matlab.rar以后，先别急着双击跑脚本。这套工具依赖两个底层能力：NetCDF 文件读取和m_map绘图工具箱。MATLAB 自带netcdf函数，但版本较老的 R2016a 之前是nctoolbox，脚本里如果写的是nc = netcdf(...)，在 R2020 以后会出现提示改用ncread。我在本地用 MATLAB 2023b 测试时，直接把netcdf调用改成了ncread加ncinfo的组合，改动量不大。

环境配置的步骤我用过很多次，最稳的流程是：

unzip fvcom-matlab.rar -d ~/tools/fvcom-matlab cd ~/tools/fvcom-matlab matlab -nodisplay -nodesktop

% 在 MATLAB 中执行 addpath(genpath('~/tools/fvcom-matlab')); savepath;

逻辑说明：addpath(genpath(...))会递归添加所有子目录，避免脚本之间相互调用时找不到函数。savepath把路径存进pathdef.m，这样下次启动 MATLAB 不用重新添加。如果用的是 2021b 之后的版本，注意netcdf命令已经被ncread取代，脚本里遇到netcdf报未定义函数时，用which -all ncread检查工具箱是否安装。很多人在 matlab 下载安装教程里只关注主程序，忽略了Mapping Toolbox和Parallel Computing Toolbox，而这两个工具箱是后处理流程里轮子最多的部分。

4.2 运行 node_name.m 生成节点文件

node_name.m这个脚本比较轻量，作用是把网格节点编号和经纬度、水深整理成一个文本文件，方便外部工具查看。它的输入约定往往是node_name(x, y, h, filename)，运行后生成一个三列或四列的 ASCII 文件。我用它来快速检查网格范围是否覆盖研究区：

node_name(x_node, y_node, h, 'grid_nodes.txt'); % 读取生成文件的前几行检查坐标范围 fid = fopen('grid_nodes.txt'); data = textscan(fid, '%f %f %f %f', 'HeaderLines', 1); fclose(fid); fprintf('lon range: %.3f - %.3f\n', min(data{2}), max(data{2}));

逻辑说明：输出文件第一列是节点编号，第二第三列是经纬度，第四列是水深。textscan读回来以后先看经纬度范围是否符合预期，比如中国近海案例应该在 120°E 左右。如果发现坐标偏移，多半是投影坐标系没转成经纬度，FVCOM 的x和y可能是 UTM 坐标，此时要先用deg2km或m_map的m_ll2xy逆变换转换。

4.3 常见报错与排查

我在跑这套脚本时总结过一张排错表，按出现频率排序：

排查时我会用dbstop if error打断点，这样出错后能直接看工作区变量的维度和值。这套脚本是十多年前的风格，很多函数没有输入校验，出错位置往往在脚本内部而非调用处，所以看变量维度比看报错行号更有效。

5. 把后处理脚本变成自动化流程：动画与敏感性分析

5.1 用 suptitle.m 组装多子图动画

suptitle.m的作用是给整个 figure 加一个总标题。单独用意义不大，但配合时间循环画动画时，它是串联多个子图的关键。FVCOM 输出通常几十上百个时间步，手工一张张截图不现实。我一般写一个循环，每个时间步生成四个子图：平面流场、温盐断面、水位时间序列、端到端对比，然后saveas导出 PNG，最后用 MATLAB 自带writeVideo合成视频：

vidObj = VideoWriter('fvcom_animation.mp4', 'MPEG-4'); open(vidObj); for it = 1:10:nt figure('visible', 'off'); subplot(2,2,1); plot_usmean_onz(u(:,:,it), v(:,:,it), x_ele, y_ele); subplot(2,2,2); plot_zsec_uts(temp(:,:,it), x_node, y_node, siglay, time(it)); subplot(2,2,3); plot(time(1:it), zeta(node_id, 1:it)); subplot(2,2,4); quiver(x_ele, y_ele, u(:,1,it), v(:,1,it), 2); suptitle(sprintf('FVCOM snapshot t = %.1f h', time(it)/3600)); frame = getframe(gcf); writeVideo(vidObj, frame); close(gcf); end close(vidObj);

逻辑说明：这个循环里u(:,:,it)的维度是单元数 x sigma 层数，temp是节点数 x 层数，两者在子图 2 中都被plot_zsec_uts处理。第四个子图画的是表层流场，quiver的矢量密度可以用参数2控制，太大图会糊成一团。getframe在visible off的窗口下依然有效，但首次调用会慢一些，属于正常现象。

5.2 敏感性分析时的参数批量替换技巧

敏感性分析经常要在不同底部摩擦系数、不同风场强迫下反复跑模型，然后用同一套后处理脚本对比结果。这时不要把脚本里的参数写死，而是用eval或函数句柄批量生成文件名：

cases = {'cf_0.001', 'cf_0.005', 'cf_0.01'}; for i = 1:length(cases) ncfile = sprintf('output_%s.nc', cases{i}); u = ncread(ncfile, 'u'); zeta = ncread(ncfile, 'zeta'); % 计算并保存统计量 u_mean = mean(u, 3, 'omitnan'); save(sprintf('u_mean_%s.mat', cases{i}), 'u_mean', 'zeta'); end

逻辑说明：omitnan选项在 R2018b 以后才支持，旧版本要写nanmean(u, 3)。用sprintf拼接文件名是这套流程里最容易出错的地方，因为 FVCOM 输出文件通常带_0001.nc这样的序号，直接替换中间段很容易路径对不上。我会先把文件列表用dir('output_*.nc')读一遍，再用regexp提取序号，而不是硬编码编号序列。

这套工具包的价值在于：它不是一个封装好的商用后处理软件，而是让你能改、能拆、能接进自己流程的源程序。把interp1_fast.m和trim_uv.m的输入输出接口研究透，你就能在它基础上拼出自己的前处理和后处理流程，而不是每次新项目都从零开始画图。

本文还有配套的精品资源，点击获取