【MATLAB】海洋热浪 m_mhw 工具箱

——从零开始在 MATLAB 里检测与分析海洋热浪 (MHW) / 海洋冷潮 (MCS)

Contents

1 | 准备工作：5 步搞定 m_mhw 的运行环境

目标：让第一次接触 m_mhw 的你，也能在 10 分钟内跑出第一个 Marine Heatwave 检测脚本。

步骤 ①　获取源码

最快方式：点击 GitHub ➜ https://github.com/ZijieZhaoMMHW/m_mhw1.0 → “Code ▼” → “Download ZIP”。
熟悉 Git 的同学：
```
git clone https://github.com/ZijieZhaoMMHW/m_mhw1.0.git
```
下载完就是一个普通文件夹，没有安装程序，也不用管理员权限。

步骤 ②　将工具箱加入 MATLAB 路径

在 MATLAB 命令行输入（把路径替换成你的真实位置）：

addpath(genpath('D:/toolboxes/m_mhw1.0'));
savepath          % 保存到 permanent path，下次启动仍然可用

这样 m_mhw 的所有函数 (detect.m, event_line.m 等) 就能被任何脚本直接调用。

步骤 ③　检查依赖

必需：MATLAB 2019b 及以上；Statistics and Machine Learning Toolbox（用于百分位计算、回归等）。
推荐：可视化地图时安装 m_map。若只做数值分析，可不安装。
平台：Windows / macOS / Linux 均已验证可用。

步骤 ④　准备示例数据

工具箱附带一份 1982–2015 东塔斯马尼亚日均 SST（example_sst.mat），新手练习足够。若想用自己数据：

确保时间分辨率为每日，单位为 °C（开尔文需减 273.15）。
用 NaN 标记陆地或缺测。
把时序转成 MATLAB 的 datenum（见后文示例脚本）。

步骤 ⑤　可选的并行加速

detectc 在高分辨率网格上可用多核并行，大幅缩短计算时间。首次使用可这样启动默认本地池：

parpool               % 或 parpool('local')

随后 detectc 会自动在 parfor 循环里调用多核。对 ≤0.25° 的区域网格数据，速度可提 2–5 倍。

小贴士

路径忘了 savepath？下次开 MATLAB 会找不到函数；重新加路径即可。

没装 “Statistics Toolbox”？detect 会提示 internal.stats.parseArgs 缺失，需要勾选安装。

数据太大内存吃紧？可以先裁剪研究区或转换为 single 精度再运行。

2 | 理论预热：3 个关键词搞懂 “极端海温事件”

先厘清概念，再上手编程；理解阈值与强度分级，你运行 m_mhw 时才知道结果在说什么。

2.1 什么是 Marine Heatwave（MHW）？

想像你家附近的海面温度，每年这个时节大概在 24 °C 左右。MHW 的判定就像是给这条“每年今日温度”排了一个成绩：把过去 30 年同一日的温度从低到高排队，位于 90 % 位置的那一格，就是“平年上限”。如果今年这片海域的温度，连续至少 5 天 都“踩线”甚至超线，就正式进入“海洋热浪”状态。

关键字 90 % 分位 + 连续 5 天
工具箱默认窗口平滑气候态，因此不会被单点噪声误判

2.2 Marine Cold Spell（MCS）——海洋冷潮

冷潮是热浪的“镜像”：温度跌破 10 % 分位且持续 5 天以上。对生态系统而言，大幅降温也会造成珊瑚低温白化、浮游动物被“冻”得迟钝等影响。调用 m_mhw 时，只需把参数 'Event','MCS','Threshold',0.1，剩下逻辑一样。

2.3 强度分级（Category I–IV）——判断火力大小

m_mhw 里每个事件自动按 Hobday 等人（2018） 的方案打分：

分级	热浪强度 vs. 阈值差	日常译法
I	1–2 倍	“轻度”
II	2–3 倍	“中度”
III	3–4 倍	“重度”
IV	≥ 4 倍	“极端”

代码实现很直观：在 detect.m 里，将 实际距气候态的正温差 除以 阈值与气候态之差 (mca)，向下取整；结果超过 4 的一律写成 4。

mhwcategory = floor(manom ./ mca);   % manom = 实测 - 气候态
mhwcategory(mhwcategory > 4) = 4;    % 截顶

2.4 为什么要关心这些指标？

生态学：热浪可导致海草场、珊瑚群落瞬间“高烧”，冷潮则可能让暖水物种被迫迁移。
渔业管理：了解强度与持续时间，有助于制定休渔期或养殖补救措施。
气候属性：长时间序列的 MHW/MCS 趋势能反映海洋变暖或变冷的速率。

读懂了 90 % / 10 % 分位与 Category I–IV，你就掌握了 m_mhw 最核心的判断逻辑。接下来，我们将在 第 3 章 开始动手准备数据并运行检测脚本。

3 | 数据标准化与 QA/QC：让原始 SST “干净”起来

为什么要做这一步？
Marine Heatwave 检测对时间轴和阈值计算极为敏感。如果输入数据漏了一天、温度单位混成了 K、或者把陆地当海洋，后面跑出的结果都会“离谱”。本节带你一步步把原始 NOAA OISST v2.1 整理成 m_mhw 接受的格式。

3.1 统一时间轴 ➜ 日分辨率

m_mhw 默认一格就是一天，所以要保证每天都有数据，哪怕是 NaN 也要占位。
下面示例把 time 变量（以 1800-01-01 为零点的天数）转换成 MATLAB 的 datenum：

sst  = ncread('sst.day.mean.1982-2023.nc','sst');      % (lon,lat,time)
rawt = ncread('sst.day.mean.1982-2023.nc','time');     % days since 1800-01-01
time = datenum(1800,1,1) + double(rawt);               % 转为 datenum

3.2 闰年检查：有没有 2 月 29 日？

气候态要按历日对齐，少一天／多一天都会让阈值错位。快速扫一眼：

hasFeb29 = any(day(datetime(time),'dayofyear') == 60);   % 60 = Feb-29
fprintf('Leap day present? %d\n',hasFeb29);

1：闰日存在，OK
：说明数据不含 2-29，m_mhw 会自动插值，但最好确认原始数据是否被裁切

3.3 单位校正：全部用 °C

有些卫星或模式数据默认是 K，只需减去 273.15：

if max(sst(:)) > 200    % 粗略判断 >200 说明大概率是 K
    sst = sst - 273.15;
end

3.4 陆地掩膜 & 缺测标记

m_mhw 把 NaN 视为“此处不检测”，所以：

先找一帧代表性图层（如首日），把本身就是 NaN 的像素当作陆地掩膜；
在所有时间维度上把对应像素置 NaN。

landmask = isnan(sst(:,:,1));                     % 假设首日已标好陆地
sst(repmat(landmask,1,1,size(sst,3))) = NaN;      % 3-D 扩展掩膜

小贴士：如果数据把陆地填成 0 或 9999，需要先改成 NaN，否则会被误认为“极冷”或“极热”！

3.5 （可选）网格重采样 / 缺口插值

分辨率太高跑不动？
用 MATLAB griddedInterpolant 重采样到较粗网格。
偶发缺测点？
可以先用 inpaint_nans（File Exchange）或简单空间插值，再交给 m_mhw。

3.6 准备就绪 → 保存中间文件

整理好的数据最好立即存为 .mat，后面脚本可直接 load，省得每次都解 netCDF：

save clean_sst.mat sst time               % 压缩可用 -v7.3

完成以上 QA/QC，你的 sst 和 time 已经满足：

日尺度连续（闰日 OK）
单位统一 °C
陆地 / 缺测 = NaN

现在可以进入 第 4 章：调用 detect / detectc 开始检测 MHW / MCS。

4 | 核心实战：跑出你的第一批 MHW / MCS 事件

4.1 选择时间窗

气候态区段（用来算阈值）——常用 30 年：1983-01-01 → 2012-12-31
检测区段（真正要研究的时长）——这里举例整段：1982-01-01 → 2023-12-31

cli_start = datenum(1983,1,1);
cli_end   = datenum(2012,12,31);
det_start = datenum(1982,1,1);
det_end   = datenum(2023,12,31);

4.2 认识 `detect` / `detectc` 的关键开关

下表改成逐条说明，不再使用矩阵式表格。

'Event' — 事件类型
- 'MHW'（默认）检测热浪
- 'MCS' 检测冷潮
'Threshold' — 分位阈值
- 0.9 → 90 %（热浪常用）
- 0.1 → 10 %（冷潮常用）
- 想更严格？热浪可设 0.95
'WindowHalfWidth' — 日气候态滑窗半宽
- 默认 5 ⇒ 实际窗口 11 天
- 若关注快速爆发，可收窄到 3；做月尺度研究可放宽
'smoothPercentileWidth' — 对每日阈值再平滑
- 默认 31 天；做季节或年分析时推荐调到 61
'minDuration' — 事件最短持续天数（默认 5）
'maxGap' — 合并事件时允许的最大缺口（默认 2 天）
'percentile' — 分位算法
- 'matlab'（默认）：速度快，细节略与 numpy 不同
- 'python'：完全对齐 python / R 版算法

4.3 检测脚本：热浪 & 冷潮一次搞定

%% 1) Marine Heatwaves
[MHW, mclim, m90, mhw_ts] = detectc( ...
    sst, time, cli_start, cli_end, det_start, det_end, ...
    'Event', 'MHW',  'Threshold', 0.9, ...
    'minDuration', 5, 'maxGap', 2);

%% 2) Marine Cold Spells
[MCS, mclim10, p10, mcs_ts] = detectc( ...
    sst, time, cli_start, cli_end, det_start, det_end, ...
    'Event', 'MCS',  'Threshold', 0.1, ...
    'minDuration', 5, 'maxGap', 2);

save result_MHW_MCS.mat MHW MCS mclim mclim10 m90 p10 mhw_ts mcs_ts

为什么用 detectc？
它把每个网格的事件表存在 cell 里，明显节省内存、还能配合 parfor 加速。若你更喜欢完整的 table 格式，把 detectc 换成 detect 即可。

脚本执行完会得到四个核心输出：

MHW / MCS — 事件清单
mclim / mclim10 — 366 天气候态
m90 / p10 — 阈值曲线（已自动对 2-29 插值再平滑）
mhw_ts / mcs_ts — 逐日强度场
- 事件日 ⇒ 温度距气候态的正差（°C）
- 非事件 ⇒
- 陆地 / 缺测 ⇒ NaN

4.4 读懂事件清单（以 `MHW{i,j}` 为例）

运行

head(MHW{35,18})   % 查看某一网格的事件列表

你将看到一个 N × 10 数值矩阵，10 列分别是：

mhw_onset — 事件开始日期（YYYYMMDD）
mhw_end — 结束日期
mhw_dur — 持续天数
int_max — 峰值强度 (°C)
int_mean — 平均强度 (°C)
int_var — 强度方差 (°C²)
int_cum — 累积强度 (°C × day)
xloc — 网格列号
yloc — 网格行号
category — Hobday I–IV（>4 已截顶为 4）

小贴士

想快速算“年均 MHW 天数”或“年最大强度” ➜ 直接跳到下一章 mean_and_trend。

如果你的数据分辨率特别高，而且只关心区域平均，可以先把 mhw_ts 空间求均后再做后续分析，速度更快、文件更小。

到这里，你已经完成 阈值计算 ➜ 事件检测 ➜ 结果保存 的闭环。接下来就可以进入“画图＋统计”的阶段，看看你的海域究竟经历了怎样的热浪和冷潮！

5 | 结果可视化与快速验证

本节展示两种最常用的“眼见为实”方法：

单点时间剖面——确认算法有没有把热浪画在正确的地方
空间统计地图——从大局看年均发生天数以及它们在变快还是变慢

5.1 看一条站点“心电图”——`event_line`

目的：先挑一个你关心的像素（例如近岸浮标或礁区），把原始温度、气候态、阈值和检测到的事件叠在一张图里，肉眼检查逻辑是否合理。

loc = [35 18];                     % x=35, y=18 —— 先用网格编号
event_line(sst, MHW, mclim, m90, ...
           loc, 1982, ...          % 数据起始年份
           [2019 1 1], [2021 12 31], ...  % 想看的时间段
           'LineWidth',1.5, ...
           'Color',[1 0.3 0.3], ... % 热浪填充色
           'Alpha',0.4);           % 透明度
title('Near-shore SST & MHW (2019-2021)');
ylabel('SST (°C)');

读图小贴士

灰橙曲线 → 平年气候态
绿色曲线 → 90 % 阈值
蓝色曲线 → 实际温度
红色半透明块 → 被认定为 MHW 的区段

如果发现蓝线明显超过绿色却没有被填红，说明 minDuration 或阈值设置太苛刻；反之，填红但蓝线几乎没超，也许阈值算错了。

5.2 把热浪“涂”在地图上——`mean_and_trend`

场景：科研文章、报告或 PPT 往往要回答两个问题——

“哪里更容易起热浪？” → 年均发生天数

“近年来有没有变得更频繁？” → 趋势斜率

[meanDays, ~, trendDays, pDays] = mean_and_trend( ...
        MHW, mhw_ts, 1982, 'Metric', 'Days');

年均 MHW 天数

figure
pcolor(meanDays'); shading flat
colorbar; caxis([0 100])
title('1982–2023  年均 Marine Heatwave 天数');

蓝绿色区域说明一年里只有零星几天热浪；橙红色说明接近全年都在高温状态。

趋势——一年增加/减少多少天？

figure
trendDays(pDays==0) = NaN;   % 过滤掉不显著的格点 (p ≥ 0.05)
pcolor(trendDays'); shading flat
colorbar
title('MHW 天数趋势  (天 · 年^{-1}，p<0.05)');

正值（暖色） = 每年热浪天数在增加
负值（冷色） = 正在减少或出现冷潮倾向
白色空洞 = 变化不显著

mean_and_trend 内部流程：

把 mhw_ts 切成逐年栅格指标序列
逐格调用 regress 做线性拟合 → 斜率 (trendDays) 与 p 值 (pDays)

这样，既能看到“长期平均”又能量化“变化速度”，一张图清晰回答“哪里在变、更变得多快”。

现在你拥有 逐点剖面 与 全域统计 两把可视化利器，可以轻松核对检测结果、截屏做报告，甚至直接投稿配图。接下来，你可以尝试：

换成 'Metric','Duration' 画“热浪持续天数”地图
或调用下章的合成分析，把热浪与 ENSO / PDO 位相对比

6 | 事件-合成分析：`composites` 一把梭

目标场景：

对照型——比较 “事件年” vs “平常年” 的平均态差异

前后型——查看某些日期集合（如 MHW 峰值日）对任何 3-D 变量的“平均贡献”

composites 函数把“你指定的日子”作为权重，自动在时间轴上对齐并求平均；闰年和缺测天数无需额外处理。

6.1 搭积木：准备索引 & 数据

%% ① 读取 ENSO 指标（ONI：Oceanic Niño Index）
oni = load('oni.txt');     % 两列：YYYYMM  和  值
is_enso = oni(:,2) > 0.5;  % 例：ONI > 0.5 定义为 El Niño 月

%% ② 把“事件月份” 转成 datenum 日期
enso_dates = datenum( ...
      fix(oni(is_enso,1)/100), ...   % 年
      mod(oni(is_enso,1),100), ...   % 月
      15);                           % 取月中 15 号代表整月

技巧：enso_dates 可以是任何长度的日期向量——单独某几天、某些时间段起止、甚至另一段 time 长度的布尔索引。

6.2 （可选）剔除长期趋势

合成分析通常关注 年际异常，而不是全球变暖信号。自己写一个简单三维去趋势函数，例如：

function out = detrend3D(dat)
    sz = size(dat);
    out = nan(sz);
    for i = 1:sz(1)
        for j = 1:sz(2)
            tmp = squeeze(dat(i,j,:));
            t   = (1:length(tmp))';
            if all(isnan(tmp)), continue, end
            p   = polyfit(t, tmp, 1);    % 线性趋势
            out(i,j,:) = tmp - polyval(p,t);
        end
    end
end

然后：

sst_detrend = detrend3D(sst);   % 把多年线性趋势先扣掉

6.3 一行搞出 El Niño 合成场

c_enso = composites(sst_detrend, time, enso_dates);

背后逻辑（源码精简版）

对输入 time 建立一个 event indicator：落在 enso_dates 中的天记 1，其余 0。
计算每个时刻 “事件权重 = 1 ÷ 事件数”，再扩展到与 d 同形状。
逐日加权求和（nansum），最终得到与原空间维度一致的矩阵。
- 闰年多出的 2-29 天，若不在 enso_dates，权重即为 0——自动忽略
- 若某网格该日为 NaN（缺测/陆地），nansum 自动跳过

6.4 可视化与快速解读

figure
pcolor(c_enso'); shading flat; colorbar
title('\DeltaSST during El Niño (°C)');

正值（红-黄）：El Niño 时比“背景态”偏暖
负值（蓝-绿）：偏冷
若想评估显著性，可用 ttest 或自旋置换 (spatial permutation) 做网格级检验。

6.5 延伸玩法

思路	做法示例
MHW 峰值日合成风速	把 `mhw_ts > 0` 的日期提出来当 `index`，对风速场做 `composites`
冷潮期间表层盐度异常	先跑 `'Event','MCS'`，再对 `salinity` 做合成
前后对比	自己构建 “事件中心 ±15 天” 索引，分两次调用 `composites` 并相减

完成这一节，你已经会用 composites 从“一串日期” 快速萃取任何 3-D 变量的典型异常。接下来可进入性能调优或深入统计，发掘极端事件背后的物理机制。

7 | 性能调优：大文件、高分辨率也能飞

下面按使用场景，把最实用的加速技巧写成清单，方便你对号入座。每一点都是作者在 ¼° 全球日度数据或 0.05° 区域高分辨率实验中踩坑总结出来的，亲测有效。

7.1 全球 ¼° 网格（1440 × 720 × 15 000 天）

首选 detectc + 并行
```
parpool              % 启动默认核数
[MHW,~,~,~] = detectc( ... );   % 内部自动 parfor
```
在 16-core 工作站上，全球跑 40 多年可以从 10 小时压到 2–3 小时。

按年份分块

for yr = 1982:2023
   idx = time >= datenum(yr,1,1) & time <= datenum(yr,12,31);
   detectc(sst(:,:,idx), time(idx), ...);   % 逐年写临时文件
end

这样单次内存负担小，且便于分布式任务调度（Slurm／PBS）。

7.2 机器内存吃紧

先裁剪感兴趣区：用经纬度索引直接 sst = sst(lonIdx, latIdx, :);。

改用 single：

sst = single(sst);   % 体感 40 % 内存节省，精度足够

按块写 .mat -v7.3：再大也能边读边算。

7.3 多次试参 / 对比实验

第一次跑 detectc 时，加

'SaveClim', true   % <— 自定义开关，或自行把 mclim/m90 保存

之后把保存好的 mclim、m90 传入可选参数 'ClimTemp','ClimTime'，直接跳过阈值重算，比原来快一半以上。

7.4 I/O 加速

表格改 Parquet/Feather：
```
parquetwrite('MHW.parquet', struct2table(MHW_flat))
```
读取速度比 .mat 提速 3–10 倍，且可直接被 Python/Pandas、R 调用。
分区写 NetCDF：对长期生产流程，用多文件 NetCDF-4 分年存储，方便 xarray 或 CDO 后处理。

7.5 进阶小技巧

GPU：gpuArray + 自写卷积核平滑，可把阈值计算再提速 3×。
内存映射：memmapfile 让多个进程共享同一份只读 SST 大文件，避免重复读盘。
Tall Array / Datastore：如果你只想做统计，不必一次把完整时间序列搬进内存，可以考虑 MATLAB 的大数据工具箱（但实测写自定义块循环通常更快）。

一句话总结：
并行 + 分块是王道；再配合 single 精度、预缓存阈值，就能让高分辨率全球数据在普通工作站“跑得动、跑得快”。

8 | 进阶玩法：把 m_mhw 用到飞起

下面列了 4 个常被忽视、却能显著提升研究深度或跨语言协作的小技巧。每条都配了动手示例，让你秒懂怎么改代码。

8.1 和 Python / R 完全对齐 —— `percentile = 'python'`

m_mhw 默认使用 MATLAB 的 quantile 算法；如果你的合作者用的是 heatwaveR 或 marineHeatWaves，同一数据可能出现“事件数多 / 少几例”的偏差。解决办法极简单：

[MHW,~,~,~] = detectc( ...
      sst, time, cli_start, cli_end, det_start, det_end, ...
      'percentile', 'python');      % 这里切换算法

幕后原理：改成 numpy-style 插值规则，边界处理与 R 一致；速度几乎不变。

8.2 定义属于你的 “Category V”

想在图表上突出史诗级事件？

打开 detect.m（或 detectc.m）

找到这一行（大约 270 行）

mhwcategory = floor(manom ./ mca);
mhwcategory(mhwcategory > 4) = 4;

改成：

mhwcategory = floor(manom ./ mca);
mhwcategory(mhwcategory > 5) = 5;   % 保留 5 级

现在，大于 5 × 阈值差的日子会被标为 Category V。记得在后续绘图或统计里同步调色盘和图例说明。

8.3 写成 netCDF，让 QGIS / xarray 秒读

MATLAB .mat 文件跨平台并不友好；用 NetCDF-4 既能压缩，又可被 Python、R、Panoply、QGIS 等直接打开。核心 API 只有两行：

nccreate('MHW_daily.nc','mhw_ts', ...
         'Dimensions',{'lon',xN,'lat',yN,'time',tN}, ...
         'DeflateLevel',5,'Format','netcdf4');
ncwrite('MHW_daily.nc','mhw_ts', mhw_ts);

DeflateLevel = 1–9 控制压缩比，数值越高越小但更耗时。
如果文件太大，推荐分年写，然后用 ncrcat 拼接或直接在 xarray 里 lazy-load 目录。

8.4 把“强度方差”玩出花 —— EOF on `int_var`

int_var（事件内温度方差）常被忽略，但它能揭示“热浪内部波动强烈程度”。快速流程示例：

%% 1) 取每年每格的 int_var 平均
[~, annualVar, ~, ~] = mean_and_trend( ...
        MHW, mhw_ts, 1982, 'Metric', 'MeanInt');   % 先拿平均强度
[~,~,nYears] = size(annualVar);
varField = reshape(annualVar, [], nYears)';        % (year × pixel)

%% 2) 去均值、EOF 分解
varField = detrend(varField, 'constant');          % 年际中心化
[EOFs,~,latent] = pca(varField, 'NumComponents',3);

%% 3) 画出前两个主成分
figure; scatter(1:nYears, EOFs(:,1)); title('PC1 of Int-Var');

解释思路：

PC1 可能代表“背景变暖导致整体波动增强”
PC2 或许对应 ENSO 位相调制的区域性差异

加上时间序列相关性（如与 ONI / PDO），即可写出一段机制讨论。

小结

同步算法、自定义分级、跨平台输出、再到 新指标深挖——四招带你从“跑通”走向“玩转”。
任何时候记住：m_mhw 是 MATLAB 脚本，可读可改；大胆打开源码，你会发现更多潜力。

下一步？也许你能把 MHW 分类结果送进机器学习模型，预测下一次极端热浪的来袭 —— 探索永无止境，祝你科研顺利！

9 | 常见错误 × 调试清单（无表格版）

下面把最容易踩坑的 4 类报错／异常按 “症状 → 原因 → 解决方案” 的顺序写成清单，方便你在半夜 Debug 时快速对号入座。

9.1 “全图一个事件都没有！”

可能原因
1. 阈值设得太离谱：Threshold 0.99（过高）或 0.5（过低）。
2. 温度单位混用：数据仍是 K，实际值高到 300 + 。
3. minDuration 太大：例如硬性设成 15 天，结果很多 5–10 天的小热浪被过滤。

快速自检

max(sst(:))      % >200 基本就是 K，需要减 273.15

一键修复

sst = sst - 273.15;          % 把 K 改 °C

9.2 “Feb 29 报错 / 时间向量不连续”

症状：index exceeds array bounds 或 datenum 相关错误，常在 2-29 附近出现。
原因：时间轴缺了闰年那一天，或被别人裁切成 365×N 天。

排查

missingLeap = ~any(day(datetime(time),'dayofyear')==60);
fprintf("Leap day missing? %d\n", missingLeap);

解决
- 若确实缺 2-29，可手动插一行 NaN，或用气候值插入。
- 确保 time 向量与 sst 第三维长度一致。

9.3 “内存溢出 / MATLAB 崩溃”

典型场景：全球 ¼°、40 年数据一次性读进来。
招式
1. 换 detectc 比 detect 少一半内存，占用视情况可减到 -70 %。
2. 裁剪区域 只留研究海域。
3. 改 single sst = single(sst); 立省 40 % 内存。
4. 按块循环 每次只跑一年或一个经向带，然后合并结果。

9.4 “结果跟 heatwaveR / marineHeatWaves 对不上”

常见差异：事件数量少几例、强度略偏高。
三件事先确认
1. 'percentile','python' 是否打开？
2. WindowHalfWidth 与对方脚本一致吗？
3. 平滑窗口 smoothPercentileWidth 是否同为 31 天？
调和示例
```
[MHW,~,~,~] = detectc(sst,time,...,'percentile','python');
```
这样 m_mhw 就与 numpy 分位实现保持一致，绝大多数差异会消失。

最后的小技巧
如果仍然有疑惑，第一步永远是画图：

用 event_line 检查蓝线（实测）和绿线（阈值）的相对位置；

或者随机抽一个 “莫名其妙未被检测出” 的时间段，看看是不是因为缺测 (NaN) 或闰日错位导致。

10 | 参考文献 & 资源（扩展版）

想把 m_mhw 玩得更深、更广，下面这份“加厚版”资源单可以充当你的 阅读地图 + 工具导航。全部文献都附 DOI；数据和代码给出官方入口，方便你随时跳转。

10.1 理论基石 & 方法学

Hobday, A. J. et al. (2016). A hierarchical approach to defining marine heatwaves. Progress in Oceanography, 141, 227-238. DOI: 10.1016/j.pocean.2015.12.014
Schlegel, R. W., Oliver, E. C. J., Wernberg, T., Smit, A. J. (2017). Nearshore and offshore co-occurrences of marine heatwaves and cold-spells. Progress in Oceanography, 151, 189-205. DOI: 10.1016/j.pocean.2017.01.004
Schlegel, R. W. & Smit, A. J. (2018). heatwaveR: A central algorithm for the detection of heatwaves and cold-spells. Journal of Open Source Software, 3 (27), 821. DOI: 10.21105/joss.00821
Zhao, Z. & Marin, M. (2019). A MATLAB toolbox to detect and analyze marine heatwaves. Journal of Open Source Software, 4 (33), 1124. DOI: 10.21105/joss.01124 — m_mhw 官方论文
Oliver, E. C. J. et al. (2018). Marine heatwaves off eastern Tasmania: Trends, interannual variability, and predictability. Progress in Oceanography, 161, 116-130. DOI: 10.1016/j.pocean.2018.02.007
Holbrook, N. J. et al. (2020). A global assessment of marine heatwaves and their drivers. Nature Communications, 11, 2674. DOI: 10.1038/s41467-020-15799-1 — 综述型权威汇编
Sen Gupta, A. et al. (2020). Drivers and impacts of the largest recorded marine heatwave on the Great Barrier Reef. Nature Communications, 11, 1792. DOI: 10.1038/s41467-020-15418-5

10.2 工具箱 & 代码仓库

语言	名称	说明 / 链接
MATLAB	m_mhw	https://github.com/ZijieZhaoMMHW/m_mhw1.0 — 本教程主角，含多套示例
Python	marineHeatWaves	https://github.com/ecjoliver/marineHeatWaves — Oliver 维护，API 与 m_mhw 类似
R	heatwaveR	https://github.com/robwschlegel/heatwaveR — R 社区热浪/冷潮“标配”
Julia	HeatwaveDetection.jl	https://github.com/JuliaClimate/HeatwaveDetection.jl — 对接高性能气候生态工作流
可视化	m_map（MATLAB）、cmocean（色盘）	制作专业海图与色阶 https://matlab-mapping-tools.github.io/m_map/; https://matplotlib.org/cmocean/

10.3 开放数据集（SST & 相关场）

NOAA OISST v2.1 — 全球 0.25° 日度，1982-至今
FTP: ftp://ftp.cdc.noaa.gov/Datasets/noaa.oisst.v2.highres
GHRSST (L4 MUR, OSTIA 等) — 1 km–9 km 多源融合 SST，即时性强
ERA5 Reanalysis — 0.25°，提供海表温度（sst）、风速、辐射通量等
CMIP6 — 多模式模拟，用于未来情景热浪预估
Argo — 全球剖面浮标，可验证热浪深度结构

Tip：在 m_mhw 前期练习可直接使用 OISST；做机理研究可叠加 ERA5 风、气压场或 Argo 剖面。

10.4 气候指数 & 辅助数据库

ONI / Niño3.4 SST — https://origin.cpc.ncep.noaa.gov/products/analysis_monitoring/ensostuff/ONI_v5.php
PDO — https://www.ncdc.noaa.gov/teleconnections/pdo/
IOD (Dipole Mode Index) — https://www.jamstec.go.jp/aplinfo/sintexf/e/iod/dipole%20mode%20index.html
SAM, NAO, AO — 见 NOAA CPC 站点
这些指数可与 m_mhw 输出交互检验：composites、相关/回归、事件同步性分析。

10.5 参考配色 & Shapefile 资源

cmocean — 专为海洋科学设计的无失真色阶
GSHHG / Natural Earth — 全球海岸线、国界、海盆分区 Shapefile
ETOPO1 — 1 arc-min 地形/海底地形，做背景 DEM

10.6 引用示例（LaTeX / BibTeX）

@article{Hobday2016MHW,
  author  = {Hobday, A. J. and Alexander, L. V. and Perkins, S. E. et al.},
  title   = {A hierarchical approach to defining marine heatwaves},
  journal = {Progress in Oceanography},
  year    = {2016},
  volume  = {141},
  pages   = {227--238},
  doi     = {10.1016/j.pocean.2015.12.014}
}

复制即用；对中文版投稿请改为 GB/T 7714 格式。

本教程基于 m_mhw 函数注释与官方 README（二次整理）编写

2025 年 4 月
一	二	三	四	五	六	日
	1	2	3	4	5	6
7	8	9	10	11	12	13
14	15	16	17	18	19	20
21	22	23	24	25	26	27
28	29	30