简介面向海洋动力学与灾害预警研究的MATLAB源码集锦聚焦海啸传播中的浅水波模型涵盖从地震或滑坡触发、远场传播到近岸淹没的模拟场景适合地球科学领域研究人员、防灾减灾从业者及海洋相关专业学生用于理解模型原理与动手实践。压缩包共40个文件体积仅1.55MB轻量易部署主要包含19个MATLAB脚本、13个Fortran源文件、3个控制文件及可执行程序和COMCOT用户手册PDF构成从模型配置、数值计算到结果输出的完整工具链。已有659人学习下载。通过修改地形、初始扰动等参数可直观观察不同条件下的海啸生成、传播与衰减过程帮助读者将理论方程转化为实际计算。Fortran底层代码适合钻研数值离散方法MATLAB脚本便于数据前处理与可视化随包手册则提供模型使用与算法参考是教学演示和科研入门的实用资料。1. 海啸传播模型从浅水方程到可跑的 MATLAB 源码海啸在深水区几乎是“隐形”的波高不过一米波长却有上百公里船上的乘客根本感觉不到。真正让它具备破坏力的是近岸浅水效应——波速从每秒 200 米量级骤降到每秒十几米能量被压缩进越来越小的水体里波高因此成倍增长。这个物理过程的数学载体是浅水方程而“用 MATLAB 实现一套海啸传播模型源码”这件事的难点从来不在方程本身而在格式稳定性、地形数据预处理和边界处理这三道坎上。这篇文章不打算贴一份“完整工程”而是把一套常见做法讲透用守恒型浅水方程配合有限体积格式在 MATLAB 里搭出可运行的传播模型骨架。适合正在做课程设计、毕业设计或前期科研验证的工程师也适合想评估特定海域海啸风险的从业者。你拿到手的不只是代码而是每一步改参数时知道自己在改什么的底气。2. 控制方程与数值格式为什么越洋传播必须用守恒形式2.1 从三维 NS 方程到二维浅水方程的简化逻辑海啸的垂向尺度与水平尺度之比通常在 10⁻³ 量级垂向加速度远小于重力加速度因此压力分布近似静压。对这一条件做垂向积分就把三维 Navier-Stokes 方程化为二维浅水方程。常见的写法是水深积分形式的守恒方程组∂U/∂t ∂F/∂x ∂G/∂y S其中守恒变量和通量分别为U [h; hu; hv] F [hu; hu² gh²/2; huv] G [hv; huv; hv² gh²/2] S [0; -gh ∂zb/∂x; -gh ∂zb/∂y]h 为总水深u、v 为沿水深平均的水平流速g 为重力加速度zb 为海底高程。“守恒形式”这四个字容易被初学者跳过但它是能否稳定模拟越洋传播的分水岭。非守恒形式直接对速度 u 离散在地形剧烈变化时质量不守恒波前会“凭空”衰减或抬升守恒形式则让水量、动量在离散层面严格保持这也是工业级海啸模型普遍采用有限体积法的根本原因。用 shallow water 模型模拟海啸时波速 c sqrt(g·h)4000 米深海对应约 198 m/s即 712 km/h这个量级决定了时间步长必须精细到秒级甚至亚秒级。2.2 有限体积 HLL 格式的离散骨架有限体积法把计算域划分为网格单元更新的是单元平均值通量通过相邻单元的“黎曼问题近似解”计算。HLLHarten-Lax-van Leer格式是海啸模拟中最常用的近似黎曼求解器它只用左右两端的波速估计值构造通量不需要精确解稳定性和耗散特性都很均衡。在 MATLAB 里HLL 通量函数可以写成下面这段。这里用逐段代码说明逻辑也是后续所有算例的核心function flux hll_flux(UL, UR, g) % UL, UR: 左右单元守恒变量 [h; hu; hv] % 计算左右单元的水深、流速与波速 hL UL(1); uL UL(2) / max(hL, 1e-6); vL UL(3) / max(hL, 1e-6); hR UR(1); uR UR(2) / max(hR, 1e-6); vR UR(3) / max(hR, 1e-6); cL sqrt(g * hL); cR sqrt(g * hR); % HLL 波速估计取左右及中间波速的极值组合 SL min(uL - cL, uR - cR); SR max(uL cL, uR cR); % 物理通量 FL [hL*uL; hL*uL^2 0.5*g*hL^2; hL*uL*vL]; FR [hR*uR; hR*uR^2 0.5*g*hR^2; hR*uR*vR]; % HLL 平均通量 if SL 0 flux FL; elseif SR 0 flux FR; elseif (SL 0) (SR 0) flux (SR*FL - SL*FR SL*SR*(UR - UL)) / (SR - SL); else flux (FL FR) / 2; % 退化保护理论上不会触发 end end这段代码里最需要理解的是波速估计段。SL 取左右两侧“左行波速”的最小值SR 取“右行波速”的最大值目的是包裹真实黎曼解的所有可能波速。如果 SL 0说明整个波系在向右运动直接用左物理通量SR 0 时反过来。中间分支是真正的 HLL 平均通量它的权重与波速差成反比波速差越大数值耗散越大这也是海啸模型在间断处不产生非物理振荡的原因。2.3 干湿边界与水深下限的约定海啸模拟必然涉及“水陆交界”水深趋近于零的位置流速仍要除 h直接除零会导致 NaN 扩散到整个计算域。常见做法是设置一个干湿阈值 h_dry当单元水深小于该值时不参与通量计算或者强制流速为零。下表是几个关键参数的经验取值我一般在不同场景下会按这个范围调整参数含义常见取值调整依据h_dry干单元判定阈值1e-3 ~ 1e-2 m数值实验中出现负水深则调大CFL时间步长安全系数0.3 ~ 0.6地形越陡取值越小g重力加速度9.81 m/s²固定不变地形平滑核去除高频地形噪声3~5 个网格宽网格越粗核越大提示干湿阈值不能一味调大否则会把真实浅水区的涌波“切掉”导致近岸波高被系统性低估。3. 用 MATLAB 把 ETOPO 海床数据喂成模型的网格地形3.1 netCDF 读取与经纬度坐标的本地化真实海啸模拟必须用真实地形。NOAA 的 ETOPO1 和 ETOPO 2022、GEBCO 系列都是常用数据源分辨率从 30 角秒到 15 角秒不等。MATLAB 从 R2019a 起 ncread 函数已经非常成熟不再需要额外工具箱如果你还在用 R2014a 之前的版本则需要通过 netcdf 底层的 netcdf.open 接口操作体验会差不少。% 读取 ETOPO2022 的 netCDF 文件 fname ETOPO_2022_v1_15s_N90W180_surface.nc; ncdisp(fname); % 先看变量名和维度避免猜错 lon ncread(fname, lon); lat ncread(fname, lat); z ncread(fname, z); % 海拔单位米海床为负值 % 截取目标海域例如日本东部外海 lon_lim [135, 150]; lat_lim [28, 45]; idx_lon find(lon lon_lim(1) lon lon_lim(2)); idx_lat find(lat lat_lim(1) lat lat_lim(2)); lon_sub lon(idx_lon); lat_sub lat(idx_lat); z_sub z(idx_lat, idx_lon); % 注意 lat 是第一维很多人在这一步栽跟头ncread 返回的数组维度顺序与常识相反。ETOPO 系列通常先存纬度再存经度所以 z_sub 的第一维是 lat、第二维是 lon。如果不确认先 ncdisp 看变量声明里的 dims 顺序再决定索引方式避免地形沿经纬度“躺倒”。3.2 从经纬度到等距网格为什么要做投影转换浅水方程本身是笛卡尔坐标下的形式直接用经纬度喂进方程会引入球面曲率误差。常见的处理有三种一是小区域直接用等经纬度近似dx 取该纬度上的实际公里数二是使用 UTM 投影将经纬度转为米制坐标三是球坐标浅水方程在方程里加入曲率项。对越洋海啸模拟我一般用等距圆柱投影加局部修正选定参考纬度 lat0dx R·Δλ·cos(lat0)dy R·Δφ其中 R 取 6371000 米。这个近似在中纬度误差在 1% 以内对海啸波传播这类大尺度现象完全可接受代码量比完整球坐标方程少一个数量级。3.3 地形重采样与平滑的坑原始 ETOPO 数据分辨率 15 角秒约 450 米而模型网格往往取 1~5 公里。直接用 interp2 或 griddedInterpolant 做插值会在地形断裂带产生振铃尤其在海沟边缘一个负水深异常就可能让 HLL 通量计算出现负水深。我通常的做法是先用 imgaussfilt 对地形做高斯平滑再做降采样插值% 高斯平滑核宽度根据地形的尖锐程度调整 z_smoothed imgaussfilt(z_sub, 2); % 构造插值器linear 即可避免 spline 的过冲 F griddedInterpolant({lat_sub, lon_sub}, z_smoothed, linear, none);注意 griddedInterpolant 的 extrapolation 参数必须设为 none让越界查询返回 NaN再做一次 fillmissing 把边缘补上。地形数组里绝对不能有 NaNHLL 通量一旦遇到 NaN整个计算域的波动会在几十个时间步内全部污染。4. 完整跑通一次越洋海啸模拟初始波、CFL 步长与可视化4.1 初始水面场的两种生成方式海啸的初始条件本质上是海底地震造成的垂向位移映射为水面位移。严格做法是使用 Okada 位错模型计算断层面滑动引起的地表形变再假设海面瞬时响应于海底位移。源码层面常见做法是先用简化高斯型波面验证程序再替换为 Okada 输出。% 高斯型初始波面设震源位于 (x0, y0)特征宽度 sigma [X, Y] meshgrid(x, y); x0 500e3; y0 500e3; % 震源位置米 sigma 80e3; % 初始扰动半径米 A 1.5; % 初始波高米 eta0 A * exp(-((X-x0).^2 (Y-y0).^2) / (2*sigma^2)); h0 eta0 - z_grid; % z_grid 为海底高程注意符号真实地震的初始波高通常只有 0.5~2 米但因为波长极大传播数百公里后波高几乎不衰减只有进入浅水区后才急剧放大——这正是海啸与普通风浪的根本区别。用高斯型初始波面验证时建议把网格设成 100×100 以上能清晰看到圆形波前的扩散过程。4.2 主循环中的时间步长控制与干湿处理时间步长由 CFL 条件决定dt CFL · dx / max(sqrt(2·g·h_max))。h_max 取计算域内最大水深在深海区 6000 米水深对应的波速约 242 m/s5 公里网格下 dt 约 6 秒。实际海啸模拟经常采用自适应步长每步重新计算全局最大波速避免固定步长在地形浅水区浪费算力。dt CFL * dx / sqrt(2 * g * max(H(:))); H_new H; % 对每个单元计算 x 方向通量差分 for j 2 : ny-1 for i 2 : nx-1 % 左右单元守恒量 UL U(:, j, i-1); UR U(:, j, i); flux_x hll_flux(UL, UR, g); % 同理计算 y 方向通量 flux_y H_new(j, i) H(j, i) - dt/dx * ... (flux_x(1) - flux_left) - dt/dy * (flux_y(1) - flux_bottom); end end实际工程中几乎不会手动写双层循环因为 MATLAB 的矩阵运算并行化远快于任何手工循环优化。我是先用小网格的循环版本验证正确性再改用矢量化实现跑正式算例。4.3 用 surface 和视频输出观察波前传播可视化阶段只需一个循环帧动画即可。s surf(x/1000, y/1000, H) 并把 z 方向压缩能直观看到波前推移。建议把 ZLim 固定在 [-100, 100] 米避免近岸波峰拉爆色标。输出视频用 VideoWriter 逐帧写入帧率 15 即可波速与帧率比例合适时能看到明显的传播过程。5. 三个让模型离开课堂的进阶技巧辐射边界、GPU 加速与结果校验5.1 Sommerfeld 辐射边界让波“走出去”而不是“弹回来”固定边界在海啸模型里必须警惕波到达边界会被完全反射反射波与真实波叠加后污染整个场。常见的做法是采用一维 Sommerfeld 辐射条件把法向通量按波的传播方向外推。实现上每步对四个边界的流量做一次松弛修正代码量不大效果立竿见影。% 右边界辐射条件∂η/∂t c·∂η/∂x 0 c_bnd sqrt(g * max(H(end,:), h_dry)); eta_bnd_new eta(end,:) - c_bnd * dt/dx * (eta(end,:) - eta(end-1,:)); eta(end,:) 0.5 * eta(end,:) 0.5 * eta_bnd_new;5.2 gpuArray 加速与算力权衡海啸模型是显式格式空间离散的每个网格互不依赖天然适合 GPU 并行。MATLAB 里将地形数组和守恒变量转为 gpuArray主循环的矢量化计算会自动走 GPU 路径。但要注意GPU 显存有限5 公里网格、1000×1000 区域单精度浮点需要约 16 GB中端卡已经紧张。没人规定必须用 GPU多核 CPU 的 parfor 在中小规模算例上常常更省心。想复用现有脚本做参数扫描时我会把主循环拆成只接受地形和初始场两个入参的函数外层用 parfor 跑多组震源位置。这也让 AI 工具帮上忙——比如把 hll_flux 的内核逻辑交给代码生成工具去重构时输入输出明确、无全局变量才能保证改完能跑通。在 MATLAB 这类闭源环境中让代理像执行 Python 脚本那样操作任务前提永远是脚本本身模块化且可断言验证。5.3 用 DART 浮标数据校验传播结果计算模型跑完只是第一步可信度要靠实测数据说话。NOAA 的 DART 浮标链在太平洋部署了多个实时站点每 15 秒记录海底压力变化换算出的水位数据与模型结果对比是检验传播模型的标准动作。我通常的做法是对比首个波峰到达时间与首个波峰振幅到达时间误差在 5% 以内、振幅误差在 30% 以内算基本合格。误差超出时优先怀疑地形分辨率不足和初始波面参数标定问题而不是先去调格式参数。模型中保留一个“虚拟浮标”功能在网格中预设几个观测点坐标主循环每完成一步就把该点的 h 值记录下来模拟结束后与实测数据画在一起对比。这个功能建议从第一天就开始写否则等模型跑通再加观测代码往往要重跑一遍完整耗时模拟浪费的时间远多于写函数的时间。本文还有配套的精品资源点击获取 MATLAB海啸传播模型实战：浅水方程与有限体积法源码解析
简介面向海洋动力学与灾害预警研究的MATLAB源码集锦聚焦海啸传播中的浅水波模型涵盖从地震或滑坡触发、远场传播到近岸淹没的模拟场景适合地球科学领域研究人员、防灾减灾从业者及海洋相关专业学生用于理解模型原理与动手实践。压缩包共40个文件体积仅1.55MB轻量易部署主要包含19个MATLAB脚本、13个Fortran源文件、3个控制文件及可执行程序和COMCOT用户手册PDF构成从模型配置、数值计算到结果输出的完整工具链。已有659人学习下载。通过修改地形、初始扰动等参数可直观观察不同条件下的海啸生成、传播与衰减过程帮助读者将理论方程转化为实际计算。Fortran底层代码适合钻研数值离散方法MATLAB脚本便于数据前处理与可视化随包手册则提供模型使用与算法参考是教学演示和科研入门的实用资料。1. 海啸传播模型从浅水方程到可跑的 MATLAB 源码海啸在深水区几乎是“隐形”的波高不过一米波长却有上百公里船上的乘客根本感觉不到。真正让它具备破坏力的是近岸浅水效应——波速从每秒 200 米量级骤降到每秒十几米能量被压缩进越来越小的水体里波高因此成倍增长。这个物理过程的数学载体是浅水方程而“用 MATLAB 实现一套海啸传播模型源码”这件事的难点从来不在方程本身而在格式稳定性、地形数据预处理和边界处理这三道坎上。这篇文章不打算贴一份“完整工程”而是把一套常见做法讲透用守恒型浅水方程配合有限体积格式在 MATLAB 里搭出可运行的传播模型骨架。适合正在做课程设计、毕业设计或前期科研验证的工程师也适合想评估特定海域海啸风险的从业者。你拿到手的不只是代码而是每一步改参数时知道自己在改什么的底气。2. 控制方程与数值格式为什么越洋传播必须用守恒形式2.1 从三维 NS 方程到二维浅水方程的简化逻辑海啸的垂向尺度与水平尺度之比通常在 10⁻³ 量级垂向加速度远小于重力加速度因此压力分布近似静压。对这一条件做垂向积分就把三维 Navier-Stokes 方程化为二维浅水方程。常见的写法是水深积分形式的守恒方程组∂U/∂t ∂F/∂x ∂G/∂y S其中守恒变量和通量分别为U [h; hu; hv] F [hu; hu² gh²/2; huv] G [hv; huv; hv² gh²/2] S [0; -gh ∂zb/∂x; -gh ∂zb/∂y]h 为总水深u、v 为沿水深平均的水平流速g 为重力加速度zb 为海底高程。“守恒形式”这四个字容易被初学者跳过但它是能否稳定模拟越洋传播的分水岭。非守恒形式直接对速度 u 离散在地形剧烈变化时质量不守恒波前会“凭空”衰减或抬升守恒形式则让水量、动量在离散层面严格保持这也是工业级海啸模型普遍采用有限体积法的根本原因。用 shallow water 模型模拟海啸时波速 c sqrt(g·h)4000 米深海对应约 198 m/s即 712 km/h这个量级决定了时间步长必须精细到秒级甚至亚秒级。2.2 有限体积 HLL 格式的离散骨架有限体积法把计算域划分为网格单元更新的是单元平均值通量通过相邻单元的“黎曼问题近似解”计算。HLLHarten-Lax-van Leer格式是海啸模拟中最常用的近似黎曼求解器它只用左右两端的波速估计值构造通量不需要精确解稳定性和耗散特性都很均衡。在 MATLAB 里HLL 通量函数可以写成下面这段。这里用逐段代码说明逻辑也是后续所有算例的核心function flux hll_flux(UL, UR, g) % UL, UR: 左右单元守恒变量 [h; hu; hv] % 计算左右单元的水深、流速与波速 hL UL(1); uL UL(2) / max(hL, 1e-6); vL UL(3) / max(hL, 1e-6); hR UR(1); uR UR(2) / max(hR, 1e-6); vR UR(3) / max(hR, 1e-6); cL sqrt(g * hL); cR sqrt(g * hR); % HLL 波速估计取左右及中间波速的极值组合 SL min(uL - cL, uR - cR); SR max(uL cL, uR cR); % 物理通量 FL [hL*uL; hL*uL^2 0.5*g*hL^2; hL*uL*vL]; FR [hR*uR; hR*uR^2 0.5*g*hR^2; hR*uR*vR]; % HLL 平均通量 if SL 0 flux FL; elseif SR 0 flux FR; elseif (SL 0) (SR 0) flux (SR*FL - SL*FR SL*SR*(UR - UL)) / (SR - SL); else flux (FL FR) / 2; % 退化保护理论上不会触发 end end这段代码里最需要理解的是波速估计段。SL 取左右两侧“左行波速”的最小值SR 取“右行波速”的最大值目的是包裹真实黎曼解的所有可能波速。如果 SL 0说明整个波系在向右运动直接用左物理通量SR 0 时反过来。中间分支是真正的 HLL 平均通量它的权重与波速差成反比波速差越大数值耗散越大这也是海啸模型在间断处不产生非物理振荡的原因。2.3 干湿边界与水深下限的约定海啸模拟必然涉及“水陆交界”水深趋近于零的位置流速仍要除 h直接除零会导致 NaN 扩散到整个计算域。常见做法是设置一个干湿阈值 h_dry当单元水深小于该值时不参与通量计算或者强制流速为零。下表是几个关键参数的经验取值我一般在不同场景下会按这个范围调整参数含义常见取值调整依据h_dry干单元判定阈值1e-3 ~ 1e-2 m数值实验中出现负水深则调大CFL时间步长安全系数0.3 ~ 0.6地形越陡取值越小g重力加速度9.81 m/s²固定不变地形平滑核去除高频地形噪声3~5 个网格宽网格越粗核越大提示干湿阈值不能一味调大否则会把真实浅水区的涌波“切掉”导致近岸波高被系统性低估。3. 用 MATLAB 把 ETOPO 海床数据喂成模型的网格地形3.1 netCDF 读取与经纬度坐标的本地化真实海啸模拟必须用真实地形。NOAA 的 ETOPO1 和 ETOPO 2022、GEBCO 系列都是常用数据源分辨率从 30 角秒到 15 角秒不等。MATLAB 从 R2019a 起 ncread 函数已经非常成熟不再需要额外工具箱如果你还在用 R2014a 之前的版本则需要通过 netcdf 底层的 netcdf.open 接口操作体验会差不少。% 读取 ETOPO2022 的 netCDF 文件 fname ETOPO_2022_v1_15s_N90W180_surface.nc; ncdisp(fname); % 先看变量名和维度避免猜错 lon ncread(fname, lon); lat ncread(fname, lat); z ncread(fname, z); % 海拔单位米海床为负值 % 截取目标海域例如日本东部外海 lon_lim [135, 150]; lat_lim [28, 45]; idx_lon find(lon lon_lim(1) lon lon_lim(2)); idx_lat find(lat lat_lim(1) lat lat_lim(2)); lon_sub lon(idx_lon); lat_sub lat(idx_lat); z_sub z(idx_lat, idx_lon); % 注意 lat 是第一维很多人在这一步栽跟头ncread 返回的数组维度顺序与常识相反。ETOPO 系列通常先存纬度再存经度所以 z_sub 的第一维是 lat、第二维是 lon。如果不确认先 ncdisp 看变量声明里的 dims 顺序再决定索引方式避免地形沿经纬度“躺倒”。3.2 从经纬度到等距网格为什么要做投影转换浅水方程本身是笛卡尔坐标下的形式直接用经纬度喂进方程会引入球面曲率误差。常见的处理有三种一是小区域直接用等经纬度近似dx 取该纬度上的实际公里数二是使用 UTM 投影将经纬度转为米制坐标三是球坐标浅水方程在方程里加入曲率项。对越洋海啸模拟我一般用等距圆柱投影加局部修正选定参考纬度 lat0dx R·Δλ·cos(lat0)dy R·Δφ其中 R 取 6371000 米。这个近似在中纬度误差在 1% 以内对海啸波传播这类大尺度现象完全可接受代码量比完整球坐标方程少一个数量级。3.3 地形重采样与平滑的坑原始 ETOPO 数据分辨率 15 角秒约 450 米而模型网格往往取 1~5 公里。直接用 interp2 或 griddedInterpolant 做插值会在地形断裂带产生振铃尤其在海沟边缘一个负水深异常就可能让 HLL 通量计算出现负水深。我通常的做法是先用 imgaussfilt 对地形做高斯平滑再做降采样插值% 高斯平滑核宽度根据地形的尖锐程度调整 z_smoothed imgaussfilt(z_sub, 2); % 构造插值器linear 即可避免 spline 的过冲 F griddedInterpolant({lat_sub, lon_sub}, z_smoothed, linear, none);注意 griddedInterpolant 的 extrapolation 参数必须设为 none让越界查询返回 NaN再做一次 fillmissing 把边缘补上。地形数组里绝对不能有 NaNHLL 通量一旦遇到 NaN整个计算域的波动会在几十个时间步内全部污染。4. 完整跑通一次越洋海啸模拟初始波、CFL 步长与可视化4.1 初始水面场的两种生成方式海啸的初始条件本质上是海底地震造成的垂向位移映射为水面位移。严格做法是使用 Okada 位错模型计算断层面滑动引起的地表形变再假设海面瞬时响应于海底位移。源码层面常见做法是先用简化高斯型波面验证程序再替换为 Okada 输出。% 高斯型初始波面设震源位于 (x0, y0)特征宽度 sigma [X, Y] meshgrid(x, y); x0 500e3; y0 500e3; % 震源位置米 sigma 80e3; % 初始扰动半径米 A 1.5; % 初始波高米 eta0 A * exp(-((X-x0).^2 (Y-y0).^2) / (2*sigma^2)); h0 eta0 - z_grid; % z_grid 为海底高程注意符号真实地震的初始波高通常只有 0.5~2 米但因为波长极大传播数百公里后波高几乎不衰减只有进入浅水区后才急剧放大——这正是海啸与普通风浪的根本区别。用高斯型初始波面验证时建议把网格设成 100×100 以上能清晰看到圆形波前的扩散过程。4.2 主循环中的时间步长控制与干湿处理时间步长由 CFL 条件决定dt CFL · dx / max(sqrt(2·g·h_max))。h_max 取计算域内最大水深在深海区 6000 米水深对应的波速约 242 m/s5 公里网格下 dt 约 6 秒。实际海啸模拟经常采用自适应步长每步重新计算全局最大波速避免固定步长在地形浅水区浪费算力。dt CFL * dx / sqrt(2 * g * max(H(:))); H_new H; % 对每个单元计算 x 方向通量差分 for j 2 : ny-1 for i 2 : nx-1 % 左右单元守恒量 UL U(:, j, i-1); UR U(:, j, i); flux_x hll_flux(UL, UR, g); % 同理计算 y 方向通量 flux_y H_new(j, i) H(j, i) - dt/dx * ... (flux_x(1) - flux_left) - dt/dy * (flux_y(1) - flux_bottom); end end实际工程中几乎不会手动写双层循环因为 MATLAB 的矩阵运算并行化远快于任何手工循环优化。我是先用小网格的循环版本验证正确性再改用矢量化实现跑正式算例。4.3 用 surface 和视频输出观察波前传播可视化阶段只需一个循环帧动画即可。s surf(x/1000, y/1000, H) 并把 z 方向压缩能直观看到波前推移。建议把 ZLim 固定在 [-100, 100] 米避免近岸波峰拉爆色标。输出视频用 VideoWriter 逐帧写入帧率 15 即可波速与帧率比例合适时能看到明显的传播过程。5. 三个让模型离开课堂的进阶技巧辐射边界、GPU 加速与结果校验5.1 Sommerfeld 辐射边界让波“走出去”而不是“弹回来”固定边界在海啸模型里必须警惕波到达边界会被完全反射反射波与真实波叠加后污染整个场。常见的做法是采用一维 Sommerfeld 辐射条件把法向通量按波的传播方向外推。实现上每步对四个边界的流量做一次松弛修正代码量不大效果立竿见影。% 右边界辐射条件∂η/∂t c·∂η/∂x 0 c_bnd sqrt(g * max(H(end,:), h_dry)); eta_bnd_new eta(end,:) - c_bnd * dt/dx * (eta(end,:) - eta(end-1,:)); eta(end,:) 0.5 * eta(end,:) 0.5 * eta_bnd_new;5.2 gpuArray 加速与算力权衡海啸模型是显式格式空间离散的每个网格互不依赖天然适合 GPU 并行。MATLAB 里将地形数组和守恒变量转为 gpuArray主循环的矢量化计算会自动走 GPU 路径。但要注意GPU 显存有限5 公里网格、1000×1000 区域单精度浮点需要约 16 GB中端卡已经紧张。没人规定必须用 GPU多核 CPU 的 parfor 在中小规模算例上常常更省心。想复用现有脚本做参数扫描时我会把主循环拆成只接受地形和初始场两个入参的函数外层用 parfor 跑多组震源位置。这也让 AI 工具帮上忙——比如把 hll_flux 的内核逻辑交给代码生成工具去重构时输入输出明确、无全局变量才能保证改完能跑通。在 MATLAB 这类闭源环境中让代理像执行 Python 脚本那样操作任务前提永远是脚本本身模块化且可断言验证。5.3 用 DART 浮标数据校验传播结果计算模型跑完只是第一步可信度要靠实测数据说话。NOAA 的 DART 浮标链在太平洋部署了多个实时站点每 15 秒记录海底压力变化换算出的水位数据与模型结果对比是检验传播模型的标准动作。我通常的做法是对比首个波峰到达时间与首个波峰振幅到达时间误差在 5% 以内、振幅误差在 30% 以内算基本合格。误差超出时优先怀疑地形分辨率不足和初始波面参数标定问题而不是先去调格式参数。模型中保留一个“虚拟浮标”功能在网格中预设几个观测点坐标主循环每完成一步就把该点的 h 值记录下来模拟结束后与实测数据画在一起对比。这个功能建议从第一天就开始写否则等模型跑通再加观测代码往往要重跑一遍完整耗时模拟浪费的时间远多于写函数的时间。本文还有配套的精品资源点击获取