简介：面向古生态学、火灾历史重建研究者的Matlab开源工具包，基于沉积物炭屑记录识别局部火灾事件，提供诊断工具判断峰值检测的适用性与参数选择，分析工具则用于统计和图形化汇总结果。包内共35个文件，以24个m源代码为主，辅以csv数据模板、xls参数文件、md说明文档及用户指南PDF，整体仅1.54MB，下载与使用均很方便。目前已有354人浏览学习。代码模块划分清晰，涵盖炭屑数据预处理、背景值估计、峰值提取、混合模型分析、结果输出等完整环节，并配有示例数据和逐行注释，便于研究者理解算法原理并按需修改。无论是入门学习还是进阶改进，这套工具都能帮助科研人员快速复现经典炭屑分析方法，适合具备一定Matlab基础、希望系统开展炭屑定量分析的生态学与古气候研究者。

1. 沉积物炭分析还没过时：CharAnalysis 在拆哪条曲线、替谁省时间

拿到一根沉积岩芯，数完炭屑颗粒，你大概率会得到一条锯齿状的时间序列：数值忽高忽低，看起来有不少“尖峰”，但根本分不清哪些来自本地火灾、哪些只是背景飘灰。CharAnalysis 就是专门干这个的——它把木炭浓度时间序列拆成本底、峰值和噪声，再通过显著性检验输出一份「火灾事件清单」，包括事件年份、火频率、回返间隔。这套 MATLAB 工具是系统开源的，核心代码全部可见，可以按自己的数据改参数、换平滑算法，适合做古火灾重建、炭屑形态演化、以及把同一套峰值检测逻辑迁移到其他沉积物代用指标的研究者。我把它拆过几遍，也封装过批处理入口，下面按数据从原始记录到结果图件的完整链路写一遍，每一步都给出可复现的做法和参数边界。

2. 核心分析逻辑：从木炭浓度到火灾事件的拆分与检验

2.1 CHAR 序列的物理意义：为什么直接用浓度不行

沉积物炭分析的第一步不是跑软件，而是先确认你手里的数据是什么物理量。岩芯里数出来的通常是「浓度」——每立方厘米沉积物有多少个炭屑颗粒，或者每个样品里炭屑总面积是多少平方毫米。浓度受沉积速率影响极大：沉积物送得快，炭屑被稀释，浓度偏低；沉积速率慢，炭屑相对富集，浓度偏高。两个湖的岩芯放在一起直接比浓度，就像拿两份不同稀释倍数的溶液比溶质质量，没有意义。

所以必须先换算成累积速率（Charcoal Accumulation Rate，CHAR），单位是 particles/cm²/yr，含义是每年每平方厘米沉积面上落下来多少炭屑。这个量扣除沉积速率的影响之后，湖与湖之间、层位与层位之间才具备可比性，后续的峰值检测也不会被稀释效应带偏。

换算公式很直接：

CHAR = 浓度 × 沉积速率 / 样品厚度

举例：一个样品厚度 0.5 cm，沉积速率 0.2 cm/yr，测得浓度 400 particles/cm²，那么 CHAR = 400 × 0.2 / 0.5 = 160 particles/cm²/yr。

CharAnalysis 对输入数据的第一个要求就是 CHAR 序列，或者能换算出 CHAR 的原始参数。数据导入界面里有「浓度 + 厚度 + 沉积速率」的换算入口，也可以提前在外部换算好再导入。我一般倾向于在外部换算好，因为这样单位制由自己控制，少一层黑匣子。

2.2 背景与峰值的分离：平滑窗宽决定「证据尺度」

把 CHAR 序列画成折线图，几乎所有真实数据都是锯齿状的。炭屑的搬运过程很复杂：区域范围内的火灾会带来持续数年的低水平炭屑沉降，本地火灾会形成一个突然的尖峰，湖泊沉积还有再悬浮、再沉积的二次搬运。所以分析时要做一个基本假设：

观测 CHAR = 背景（Background） + 峰值（Peak） + 噪声（Noise）

背景代表长时间尺度的炭屑通量，由植被演替、气候趋势、区域火频率等慢变量驱动；峰值代表短时间的火灾事件；噪声是采样、计数和分析带来的随机误差。CharAnalysis 的核心工作，就是把背景和峰值拆开。

拆分背景用的是平滑：对 CHAR 序列做低通滤波，常见做法是高斯平滑或 LOWESS 平滑，窗宽一般取 200～500 年。窗宽大小直接决定你识别的是哪个尺度的事件。窗宽越大，背景曲线越平，短时间的小波动都会被算进残差、变成候选峰；窗宽越小，背景越贴着原始数据走，只有非常突出的异常点才可能超越阈值。

我自己的习惯是先看采样分辨率：如果采样间隔是 5～20 年，取 250 年窗宽，参与平滑的有效数据点有十几个到几十个，背景稳定且合理；如果采样间隔已经到 50 年，250 年窗宽只能覆盖 5 个点，背景曲线会变成折线，基本失去意义，这时候要么把窗宽放大到 500 年以上，要么想办法加密采样。

下面是一段用合成数据演示背景分离逻辑的 MATLAB 脚本，可以先把整体流程跑通，再切换到自己的真实数据。

rng(42); % 固定随机种子，保证结果可复现 t = 0:900; % 模拟900年，每年一个采样点 background = 0.8 + 0.001 .* t; % 缓慢上升的背景趋势 peaks = zeros(size(t)); peaks(150) = 3; % 第150年发生一次火灾事件 peaks(410) = 4; % 第410年发生一次火灾事件 peaks(700) = 2.5; % 第700年发生一次火灾事件 noise = 0.15 * randn(size(t)); % 分析噪声 char_obs = background + peaks + noise; % 观测到的CHAR序列 % 用250年窗宽的高斯平滑估计背景 bg_est = smoothdata(char_obs, 'gaussian', 250); residual = char_obs - bg_est; figure; subplot(2,1,1); plot(t, char_obs, 'k'); hold on; plot(t, bg_est, 'r', 'LineWidth', 2); legend('观测CHAR', '背景估计'); xlabel('年份'); ylabel('CHAR (particles/cm^2/yr)'); subplot(2,1,2); plot(t, residual, 'b'); xlabel('年份'); ylabel('残差'); title('背景扣除后剩余峰值信号');

说明一下各行的作用：rng(42)固定随机种子，保证每次运行生成的噪声一致，这是复现实验结果的底线；background是一条线性趋势，模拟植被或气候导致的长期变化；三个peaks位置的赋值模拟三次本地火灾；smoothdata的'gaussian'方法与 250 窗宽是 CharAnalysis 最常见的背景估计组合。残差序列里能看到三个明显峰，但同时也混着噪声——接下来的显著性检验就是要把「真峰」从「噪声假峰」里筛出来。

2.3 SNI 与显著性检验：把每一个峰从「疑似」变成「判据」

背景扣完之后，残差里到处都是峰，怎么判断哪些是真正的火灾事件？只看峰高不够，因为噪声水平也在变：有的层位保存条件好，噪声小；有的层位炭屑浓度低，相对噪声大。CharAnalysis 给出的方案是两步检验。

第一步，把残差标准化到 Z-score 空间。常见做法是假设残差近似服从正态分布，用局部噪声估计标准差，把每个残差值转换成 Z 分数，超过阈值的点才算候选峰。阈值一般取 2.33，对应 99% 置信水平的单侧检验——也就是说，这个点高到随机噪声产生它的概率不到 1%。把阈值设成 1.5 这种偏低值，结果会得到一堆假峰，这个后面避坑章会专门讲。

第二步是 SNI（Signal-to-Noise Index），这是个实用启发式指标：峰高与局部噪声水平之间的比值。CharAnalysis 一般要求 SNI ≥ 3.0，目的是过滤那些幅度很小、即使统计显著也没有生态意义的微峰。两步同时通过的峰，才会进入最终事件清单。

这套双条件过滤在实操中非常有用。单独跑 2.33 阈值，会有一些幅度不足但刚好卡在线上的点混进来；单独看 SNI，又可能放过长尾巴宽峰。两个条件叠加，能得到相对保守、更经得起同行评审的事件表。我的经验是，如果一序列检出几百个峰，先别怀疑数据，去查滤波条件是不是被改松了。

3. MATLAB 数据入库：CSV 转 CHAR 序列与预处理

3.1 原始数据长什么样：从记录表到标准两列 CSV

从岩芯记录到 CharAnalysis 输入，中间隔着一步数据规整。工具包本身带数据导入模块，但它的要求比较固定：要么识别 GRL 格式，要么读入你已经整理好的表格。我最推荐的做法是先把原始记录整理成标准两列 CSV——第一列年份或深度，第二列计数结果——再交给工具处理，省去软件猜单位的过程。

计数结果通常有三种形式：颗粒浓度（particles/cm²）、累积速率（particles/cm²/yr）、显微镜下的面积浓度（mm²/cm²）。这三种形式工具都认，但前提是你在导入时明确告诉它用的是哪种。如果年份列还有缺失值、重复值，或者中间穿插了空行，导入阶段就会出现各种奇怪的报错。下面这段脚本可以把分散在多个 CSV 文件里的原始记录统一读进来，并做格式校验：

% 扫描 raw 目录下所有 CSV，统一读成 year + count 两列 rawFiles = dir('raw/*.csv'); for i = 1:length(rawFiles) T = readtable(fullfile(rawFiles(i).folder, rawFiles(i).name)); year = T{:, 1}; % 第一列必须是年份 count = T{:, 2}; % 第二列是浓度或计数 assert(all(diff(year) > 0), '年份列必须严格递增'); assert(~any(isnan(count)), '计数列不能有NaN'); save(sprintf('prepared_%02d.mat', i), 'year', 'count'); end

这段代码做了两件重要的事：diff(year) > 0检查年份是否严格递增，很多野外记录表按样品编号排，导入后才发现年份是倒的；isnan检查计数列是否存在空值。保存成.mat而不是继续用 CSV，是因为后续分析要反复读取年份和计数值，.mat的读取速度更快，也不会出现 CSV 列类型被自动转换的问题。

3.2 从浓度到累积速率：厚度、沉积速率与单位陷阱

前面讲过，浓度必须先换算成 CHAR。这一步有四个变量参与：浓度、样品厚度、沉积速率，以及最终单位。最容易翻车的是沉积速率的单位——测年报告里给出的可能是 mm/yr，而公式要求的是 cm/yr，不换算直接代进去，CHAR 凭空差 10 倍。

% 读入已规整的两列数据 load('prepared_01.mat', 'year', 'count'); % 换算参数：采样与沉积信息 h = 0.5; % 样品厚度，单位 cm s = 0.2; % 沉积速率，单位 cm/yr（注意不是 mm/yr） % 如果count已经是浓度(particles/cm^2)，换算成累积速率 char_rate = count .* s ./ h; % 结果单位：particles/cm^2/yr

注释里写清楚单位，不是多此一举。我在几个数据集上反复核对过，炭屑分析最隐蔽的错误就是单位不写、隔一个月回来连自己都不确定当时用的是cm/yr还是mm/yr。如果沉积速率在序列内有明显变化，不要用一个常数，建议把它写成与年份对应的向量，逐点乘进去；三角洲、湖盆边缘这类沉积环境，沉积速率上下能差三四倍，常数近似会直接扭曲 CHAR 的长期趋势，连带背景估计一起出错。

还有一点容易被忽略：如果count本身已经是累积速率，就不要再乘s / h。判断方法很简单——看原始计数时的单位说明。野外记录里写的是「每立方厘米多少颗粒」就要乘；写的是「每年每平方厘米多少颗粒」就直接用。这个判断决定结果差三个数量级还是保持原样。

3.3 插值预处理与时间标尺统一

岩芯采样很少严格等间隔：有的层位加密采样、有的层位因为沉积物量不足跳过了，再加上测年给出的是深度-年龄模型，最终年份序列基本都不均匀。CharAnalysis 的主流程要求时间轴等间隔，否则背景平滑会不稳定，峰值位置也会偏移。

常见做法是把所有序列统一插值到固定步长网格，步长选 1 年、5 年或 10 年，取决于平均采样分辨率。插值方法我固定用线性，不要用 spline。spline 虽然让曲线更光滑，但在年份间隔大的地方会产生过冲振荡，峰值位置和幅度都会被污染，尤其是在序列两端。

% 等间隔插值到1年步长 ti = (ceil(min(year)) : 1 : floor(max(year)))'; ci = interp1(year, char_rate, ti, 'linear', 'extrap'); % 检查最大间隔，判断插值结果是否可靠 gap = diff(year); fprintf('最大采样间隔: %.1f 年\n', max(gap)); % 保存为后续分析的标准输入 save('char_prepared.mat', 'ti', 'ci', 'h', 's');

interp1的参数值得留意：'linear'是线性插值，'extrap'是让外插不返回 NaN，可以处理两端点附近的缺失；但如果序列两端外推距离太大，外插值会不可信，更好的做法是直接在后续分析里把边界区段的结果丢弃，这个在避坑章会再展开。max(gap)是一个快速诊断指标：如果最大采样间隔超过你选的窗宽的一半，说明这段序列对背景估计的支撑不足，要么加密采样，要么把分析区段截掉。这是我在处理浅层沉积物时最常跑的一条检查。

4. 关键参数与主流程：从背景窗宽到 Z-score 的一次完整解算

4.1 四个必须理解的参数

CharAnalysis 的界面里参数不算多，但每个都对结果有决定性影响。很多人第一次跑的时候用默认值出了一张图，就以为万事大吉，直到改了一个数字结果全变样。下面四个参数，建议动手前先理解：

窗宽的影响前面已经说过，这里重点说 Z 阈值和 SNI 的配合。Z-score 阈值负责过滤统计意义上的「不可能噪声」，SNI 负责过滤统计显著但幅度太小的「无聊峰」。两个条件必须同时满足，才能进最终事件表。如果你做的是探索性分析，可以先把 SNI 阈值放低看候选峰全貌；但如果要出文章，建议保持在 3.0 以上，否则审稿人会拿信噪比说事。

4.2 调用主流程：GUI 操作顺序与封装入口

工具的主程序是 MATLAB 图形界面，常规操作顺序很直接：启动工具后进入数据导入界面，选择上一步保存的.mat文件；指定分析的时间区间和插值步长；设置背景平滑窗宽；设置 Z 阈值和 SNI 阈值；运行分析并导出结果。

单数据集这样点没有问题，十几个岩芯挨个点就很痛苦了。我后来把整个分析流程封装成了一个可重复执行的函数，把参数全部显式传进去，输出存到固定目录。这是个人惯用做法，不是工具自带的官方调用方式，但逻辑与界面流程一致：

function [peak_years, peak_vals, bg_curve, fire_freq] = ... run_charanalysis(year, char_rate, opts) % 基于CharAnalysis核心流程的自用封装入口 % opts 字段： % .smooth_win : 背景平滑窗宽（年），默认250 % .z_thresh : Z-score显著阈值，默认2.33 % .sni_thresh : 信噪比阈值，默认3.0 % .step : 插值步长（年），默认1 arguments year (:,1) double char_rate (:,1) double opts.smooth_win double = 250 opts.z_thresh double = 2.33 opts.sni_thresh double = 3.0 opts.step double = 1 end % 等间隔插值 ti = (ceil(min(year)):opts.step:floor(max(year)))'; ci = interp1(year, char_rate, ti, 'linear'); % 背景平滑 bg = smoothdata(ci, 'gaussian', opts.smooth_win); % 残差与Z-score检验 resid = ci - bg; sigma = movstd(resid, opts.smooth_win); % 局部噪声估计 zscore = resid ./ max(sigma, eps); % 候选峰：超过Z阈值且SNI大于阈值的点 cand = (zscore > opts.z_thresh) & (resid ./ max(sigma, eps) > opts.sni_thresh); peak_idx = find(cand); peak_years = ti(peak_idx); peak_vals = resid(peak_idx); bg_curve = [ti, bg]; % 事件频率：按100年滑动窗口统计峰数 fire_freq = movsum(cand, 100) ./ 100; end

这段代码体现了整个分析链路：插值、平滑、残差、Z 检验与 SNI 检验、以及最终的事件频率统计。arguments块是 MATLAB R2019b 之后支持的参数声明语法，如果你用的版本更老，把opts改成struct传入即可。movstd做局部噪声估计时，窗宽直接复用背景平滑窗宽，这个选择在多数场景下是稳定的；如果序列太短，movstd会在两端产生较多 NaN，后续处理要截掉边界区段。

4.3 输出结果的解读：事件年份、火频与回返间隔

分析跑完后，生成的输出一般包含几类：峰值年份与峰高、背景曲线序列、事件频率曲线，以及按事件间隔统计的回返间隔。事件年份是最核心的产出——它是后续所有生态解释的基本单元。峰值高度本身也有意义，可以代表火灾的相对强度或炭屑生产量，但解释时要谨慎，因为它同时受火灾强度、风向、距湖距离影响。

火频率曲线是另一个常用输出：把每个峰当作一次火灾事件，按 100 年或 200 年滑动窗口统计事件数，就得到随时间变化的火频率曲线。这条曲线可以直接和植被变化、气候代用指标做对比，是古火灾研究里最常见的图表。回返间隔则是相邻两次事件之间的年数，可以拟合 Weibull 或对数正态分布，用于判断火灾发生是随机过程还是存在聚集性。

拿到输出后，我会先把峰值年份列成一张表，人工核对每个峰是否与同期其他代用指标有对应关系。跑完不检查就写进文章，后面返工的成本更高。

5. 避坑清单：五个高频错误与排查方法

这几个坑我在不同数据集上基本全踩过一遍。前两个最隐蔽，出问题的时候图形看起来完全正常，但结果解释不通；后三个最容易影响交付质量，属于审稿人一眼就能看出问题的类型。

5.1 年份列是倒序的：图形左右颠倒，背景贴合失败

现象：导入数据后，图形与标记点左右颠倒，背景平滑线看起来怎么也「贴」不上原始曲线，残差里出现大量一端高一端低的假峰。

原因：野外记录表经常按深度从下往上排列，年份随之递减。工具和脚本一般都假设年份升序，拿到递减序列时不会报错，而是直接当成升序处理，导致所有位置错位。

解决：导入前强制检查年份差分。我一般在数据入库脚本里加一行断言，年龄序列不升序就直接终止并提示：

assert(all(diff(year) > 0), '年份列必须严格递增，请执行 sortrows');

如果数据本身是乱的，先按年份排序再继续：

[year, idx] = sort(year); count = count(idx);

这里的sort会同时把计数列按相同索引重排，保证年份和计数值对应关系不变。这是一个 30 秒能完成的检查，却能省掉后面半天排查时间。

5.2 背景窗宽与采样分辨率不匹配：背景变成锯齿

现象：背景平滑曲线呈锯齿状，残差里大量一个点宽的「尖峰」，Peak 检出数量异常多，查看原始数据却看不出这些峰有明确事件特征。

原因：采样间隔 50 年，却把背景窗宽设成了 250 年。250 年窗在等间隔插值后只有 5 个有效支撑点，平滑窗口覆盖的独立数据太少，背景估计不稳定，任何局部波动都会突破检验阈值。

解决：窗宽至少覆盖 10 个以上采样点。采样间隔 20 年，窗宽别低于 200 年；采样间隔 50 年，窗宽建议直接拉到 500 年以上。还有一种方案是放弃等间隔插值到 1 年，改插到 5 年或 10 年步长，让窗宽内的有效点更多，但代价是峰值位置的时间精度降低。

5.3 Z 阈值设置过低：几百个「假火灾事件」

现象：同一份数据，把 Z 阈值从 2.33 降到 1.5，检出峰数量从几十个变成几百个，事件清单看起来「哪哪都是火」。

原因：1.5 对应的置信水平只有约 93%，十个随机点里就有接近一个会超过阈值，加上噪声序列在平滑后存在自相关，连续的高点还会连成一片，虚假倍率进一步放大。

解决：Z 阈值保持 2.33 或更保守的 2.5，同时打开 SNI ≥ 3.0 双条件过滤。SNI 过滤掉「统计显著但幅度很小」的峰，比单纯提高 Z 阈值更合理——它不改变峰的位置判断标准，只改变峰是否被接受。

5.4 CHAR 值与文献差三个数量级

现象：算出的 CHAR 均值在 0.001～0.01 particles/cm²/yr 量级，而同类文献普遍在 1～100 量级，结果图里所有峰都缩在坐标轴底部。

原因：单位换算错位。常见的有三种：沉积速率用了 mm/yr 却当成 cm/yr，计数结果实际已经是累积速率却又乘了一遍s / h，或者面积浓度（mm²/cm²）和颗粒浓度（particles/cm²）混用。

解决：在换算脚本里给每个变量写清楚单位并做一次手工校验。取一个已知样品，用笔头算一遍 CHAR，再和脚本输出对比。我每次换数据集都会跑一次这个检查，两分钟不到，能拦住 90% 的单位错误。

5.5 序列两端的峰不可信：边界效应导致假峰

现象：事件清单里密集出现序列最开始和最后结束的 50 年内的峰，而且与上下文事件完全不搭。把这几段去掉，事件图反而变干净了。

原因：平滑滤波在两端存在边界问题——窗口中心靠近边界时，窗口一侧没有数据，背景估计被拉向边界值，残差也随之为负或为正，造成系统性假峰。这不是数据分析错误，是滤波算法的固有边界效应。

解决：分析时给序列前后各加一段缓冲时间。一种做法是把分析区段扩展 50～100 年，跑完后再只取中间区段的结果；如果数据本身覆盖到岩芯顶端，就干脆把顶部最近 50 年的结果标记为「受边界效应影响」，不进入结论讨论。这比事后硬解释更严谨。

6. 结果验证与进阶应用：模拟数据自检与输出诊断

工具跑完不是终点。我对任何一套新数据，都会先做一遍「已知事件自检」：构造一个合成的 CHAR 序列，里面埋好三到五个已知年份的大事件、若干小事件，然后跑完整分析流程，看检测出来的峰值年份能不能和埋入事件对得上。这一步做的是「检出率」验证——工具和参数对已知事件有响应，才能对未知事件给出可信结果。

具体做法是：先用合成数据跑出事件年份，再和已知事件表比对，统计漏检率。如果漏检率超过 20%，说明当前参数组合对本数据集的信号强度不匹配，需要调整窗宽或阈值，而不是强行解释结果。这个流程花不了几分钟，但能避免把参数问题误判成生态信号。

进阶一点，可以做回返间隔分布拟合。把所有峰值年份做差分得到间隔序列，用 Weibull 分布拟合形状参数和尺度参数。形状参数明显大于 1 时，说明火灾事件在时间上存在聚集趋势，可能受气候波动调制；接近 1 时，说明火灾发生接近随机泊松过程。这个指标比简单平均回返间隔信息量大得多，可以直接接进生态建模。工具输出的峰值年份表，本身就是回返间隔分析的输入，导出后用 MATLAB 的fitdist一行即可完成拟合：

intervals = diff(peak_years); % 相邻火灾事件的间隔 pd = fitdist(intervals, 'Weibull'); % 拟合Weibull分布

拟合前先把小于 1 年的间隔去掉——同一场火灾在插值后可能被识别成两个相邻峰，这不代表真实回返间隔，属于分析伪影。

最后说一个交付习惯。有一次我用旧版设置跑浅湖沉积数据，SNI 阈值没有正确传入，输出报告里的峰数量翻了接近三倍，交付前才发现。从那以后，我每次提交分析结果前都强制走一遍「已知事件 → 检出率 → SNI 复检」三步流程，还会把参数记录写进结果文件的注释行，保证任何一步可追溯。这套习惯比任何算法微调都更能保证结论稳定。希望帮到你。

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