简介面向古生态学、火灾历史重建研究者的Matlab开源工具包基于沉积物炭屑记录识别局部火灾事件提供诊断工具判断峰值检测的适用性与参数选择分析工具则用于统计和图形化汇总结果。包内共35个文件以24个m源代码为主辅以csv数据模板、xls参数文件、md说明文档及用户指南PDF整体仅1.54MB下载与使用均很方便。目前已有354人浏览学习。代码模块划分清晰涵盖炭屑数据预处理、背景值估计、峰值提取、混合模型分析、结果输出等完整环节并配有示例数据和逐行注释便于研究者理解算法原理并按需修改。无论是入门学习还是进阶改进这套工具都能帮助科研人员快速复现经典炭屑分析方法适合具备一定Matlab基础、希望系统开展炭屑定量分析的生态学与古气候研究者。1. 沉积物炭分析还没过时CharAnalysis 在拆哪条曲线、替谁省时间拿到一根沉积岩芯数完炭屑颗粒你大概率会得到一条锯齿状的时间序列数值忽高忽低看起来有不少“尖峰”但根本分不清哪些来自本地火灾、哪些只是背景飘灰。CharAnalysis 就是专门干这个的——它把木炭浓度时间序列拆成本底、峰值和噪声再通过显著性检验输出一份「火灾事件清单」包括事件年份、火频率、回返间隔。这套 MATLAB 工具是系统开源的核心代码全部可见可以按自己的数据改参数、换平滑算法适合做古火灾重建、炭屑形态演化、以及把同一套峰值检测逻辑迁移到其他沉积物代用指标的研究者。我把它拆过几遍也封装过批处理入口下面按数据从原始记录到结果图件的完整链路写一遍每一步都给出可复现的做法和参数边界。2. 核心分析逻辑从木炭浓度到火灾事件的拆分与检验2.1 CHAR 序列的物理意义为什么直接用浓度不行沉积物炭分析的第一步不是跑软件而是先确认你手里的数据是什么物理量。岩芯里数出来的通常是「浓度」——每立方厘米沉积物有多少个炭屑颗粒或者每个样品里炭屑总面积是多少平方毫米。浓度受沉积速率影响极大沉积物送得快炭屑被稀释浓度偏低沉积速率慢炭屑相对富集浓度偏高。两个湖的岩芯放在一起直接比浓度就像拿两份不同稀释倍数的溶液比溶质质量没有意义。所以必须先换算成累积速率Charcoal Accumulation RateCHAR单位是 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 平滑窗宽一般取 200500 年。窗宽大小直接决定你识别的是哪个尺度的事件。窗宽越大背景曲线越平短时间的小波动都会被算进残差、变成候选峰窗宽越小背景越贴着原始数据走只有非常突出的异常点才可能超越阈值。我自己的习惯是先看采样分辨率如果采样间隔是 520 年取 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 这种偏低值结果会得到一堆假峰这个后面避坑章会专门讲。第二步是 SNISignal-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 的界面里参数不算多但每个都对结果有决定性影响。很多人第一次跑的时候用默认值出了一张图就以为万事大吉直到改了一个数字结果全变样。下面四个参数建议动手前先理解参数默认/常用值作用调整逻辑背景平滑窗宽250 年决定背景与峰值的分离尺度采样间隔大则调大序列短则调小Z-score 阈值2.33残差显著性的判定门槛噪声大时适当提高到 2.5SNI 阈值3.0过滤低幅度峰的启发式判据需要保守结果时提高探索时降低插值步长1 或 5 年统一时间标尺的网格密度按平均采样间隔的 1/31/5 取窗宽的影响前面已经说过这里重点说 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.0010.01 particles/cm²/yr 量级而同类文献普遍在 1100 量级结果图里所有峰都缩在坐标轴底部。原因单位换算错位。常见的有三种沉积速率用了 mm/yr 却当成 cm/yr计数结果实际已经是累积速率却又乘了一遍s / h或者面积浓度mm²/cm²和颗粒浓度particles/cm²混用。解决在换算脚本里给每个变量写清楚单位并做一次手工校验。取一个已知样品用笔头算一遍 CHAR再和脚本输出对比。我每次换数据集都会跑一次这个检查两分钟不到能拦住 90% 的单位错误。5.5 序列两端的峰不可信边界效应导致假峰现象事件清单里密集出现序列最开始和最后结束的 50 年内的峰而且与上下文事件完全不搭。把这几段去掉事件图反而变干净了。原因平滑滤波在两端存在边界问题——窗口中心靠近边界时窗口一侧没有数据背景估计被拉向边界值残差也随之为负或为正造成系统性假峰。这不是数据分析错误是滤波算法的固有边界效应。解决分析时给序列前后各加一段缓冲时间。一种做法是把分析区段扩展 50100 年跑完后再只取中间区段的结果如果数据本身覆盖到岩芯顶端就干脆把顶部最近 50 年的结果标记为「受边界效应影响」不进入结论讨论。这比事后硬解释更严谨。6. 结果验证与进阶应用模拟数据自检与输出诊断工具跑完不是终点。我对任何一套新数据都会先做一遍「已知事件自检」构造一个合成的 CHAR 序列里面埋好三到五个已知年份的大事件、若干小事件然后跑完整分析流程看检测出来的峰值年份能不能和埋入事件对得上。这一步做的是「检出率」验证——工具和参数对已知事件有响应才能对未知事件给出可信结果。具体做法是先用合成数据跑出事件年份再和已知事件表比对统计漏检率。如果漏检率超过 20%说明当前参数组合对本数据集的信号强度不匹配需要调整窗宽或阈值而不是强行解释结果。这个流程花不了几分钟但能避免把参数问题误判成生态信号。进阶一点可以做回返间隔分布拟合。把所有峰值年份做差分得到间隔序列用 Weibull 分布拟合形状参数和尺度参数。形状参数明显大于 1 时说明火灾事件在时间上存在聚集趋势可能受气候波动调制接近 1 时说明火灾发生接近随机泊松过程。这个指标比简单平均回返间隔信息量大得多可以直接接进生态建模。工具输出的峰值年份表本身就是回返间隔分析的输入导出后用 MATLAB 的fitdist一行即可完成拟合intervals diff(peak_years); % 相邻火灾事件的间隔 pd fitdist(intervals, Weibull); % 拟合Weibull分布拟合前先把小于 1 年的间隔去掉——同一场火灾在插值后可能被识别成两个相邻峰这不代表真实回返间隔属于分析伪影。最后说一个交付习惯。有一次我用旧版设置跑浅湖沉积数据SNI 阈值没有正确传入输出报告里的峰数量翻了接近三倍交付前才发现。从那以后我每次提交分析结果前都强制走一遍「已知事件 → 检出率 → SNI 复检」三步流程还会把参数记录写进结果文件的注释行保证任何一步可追溯。这套习惯比任何算法微调都更能保证结论稳定。希望帮到你。本文还有配套的精品资源点击获取