CharAnalysis沉积物炭屑分析实战：从CHAR计算到峰值识别与SNI诊断

简介CharAnalysis是一套面向沉积物炭记录分析的诊断与分析工具专用于检测炭化颗粒峰值重建局部火灾历史广泛适用于生态学、古气候和火灾生态研究。其诊断模块可辅助研究者判断是否需要进行峰检测并确定最合理的参数分析模块则提供统计与图形化结果汇总便于解读火灾频率和变化趋势。资源包共35个文件以24个MATLAB源文件为核心涵盖主程序、参数配置、背景炭化率计算、峰值识别等功能另含CSV数据模板、XLS示例文件、Markdown说明和PDF用户指南整体仅1.54MB轻量而完整。已有354人学习下载。源码注释详尽鼓励使用者深入理解算法并自行修改附带的示例数据与文档可快速上手是重建古火灾历史研究的高质量开源工具。1. CharAnalysis 是什么给沉积物炭屑分析拆掉“黑匣子”一根湖相沉积岩芯从野外运回实验室一层层取样、蒸解、过筛再放到显微镜下一个人数炭屑最后落到表格里的往往是深度、体积、计数这几列数字。可就是这几列数字要回答的是过去几千年火灾频率到底怎么变。CharAnalysis 作为沉积物-炭分析的诊断和分析工具处理的正是从计数到结论之间最容易变成黑匣子的那段链条浓度换算、沉积速率、背景分解、峰值识别、信噪比诊断。它适合两类人一是做古火灾重建的研究者二是要替实验室建立标准化分析流程的前处理人员。这篇笔记从一个长期搭 MATLAB 工作流的角度把每个环节的输入格式、参数含义和最常翻车的位置一次讲透。2. 从四要素数据到 CHARCharAnalysis 的数据入口与单位校验2.1 输入数据四要素深度、年龄、体积、计数做古火灾重建数据不能只有“深度-计数”两列。CharAnalysis 这类工具之所以强调诊断是因为它每一步都会反问你这个数是怎么来的单位对不对。要得到 CHARcharcoal accumulation rate炭屑累积速率原始数据至少要有四个要素要素常见单位作用深度cm标记样品在岩芯中的位置也是沉积速率计算的骨架年龄yr BP来自年代模型需要插值到每个样品深度体积cm³显微镜计数时实际扫过的样品体积用于浓度标准化计数粒该样品中炭屑颗粒的数量也可能是面积法测得的炭屑面积四要素缺一不可。没有体积浓度就算不出来没有年龄沉积速率无从谈起没有深度年龄插值就没有锚点。实际工作中我见过最多的问题是“计数、深度都有但体积忘了记”这组数据基本就废了因为同样的计数在 1 cm³ 和 0.2 cm³ 样品里含义完全不同。标准化逻辑是两步先算浓度conc count / volume单位是粒/cm³再用年龄-深度模型算沉积速率sed Δdepth / Δage单位是 cm/yr最后CHAR conc * sed单位是粒/cm²/yr。这串公式是后面所有峰值识别的地基任何一步单位错了出来的峰值都是假的。2.2 用 MATLAB 把四要素数据装进标准结构体拿到原始表之后我习惯先用 MATLAB 脚本做一次严格校验和结构体组装而不是直接进入平滑和找峰。这样做的好处是让输入格式固定下来后续换数据、换岩芯都走同一条路。% 读取沉积物炭屑分析的原始数据表 % 约定字段名depth_cm, age_yr, volume_cm3, count_grains raw readtable(site1_clean.csv, PreserveVariableNames, true); % 校验四个必要字段是否齐全 required {depth_cm, age_yr, volume_cm3, count_grains}; miss setdiff(required, raw.Properties.VariableNames); if ~isempty(miss) error(缺少必要字段: %s, strjoin(miss, , )); end % 组装成后续处理用的结构体 d struct(); d.depth raw.depth_cm; d.age raw.age_yr; d.volume raw.volume_cm3; d.count raw.count_grains; % 单位一致性校验年龄列可能被录成了公元纪年或负值 if any(d.age 0) any(d.depth 0) error(年龄列包含非正值请检查是否使用了公元纪年); end % 计算炭屑浓度、沉积速率和 CHAR d.conc d.count ./ d.volume; % 粒/cm3 d.sed gradient(d.depth) ./ gradient(d.age); % cm/yr前提是样品等间距 d.char d.conc .* d.sed; % 粒/cm2/yr % 打印前几条确认数量级和单位 disp(head(struct2table(d), 5));代码逻辑说明readtable负责把 CSV 读成表格PreserveVariableNames防止列名被改写成depthcm这种形式。结构体d比散落变量好维护后面传给平滑、找峰函数都不用重复写参数。gradient是 MATLAB 的数值差分函数这里直接算出每个采样深度的沉积速率。参数和坑gradient假设输入是均匀间距。如果样品在深度上不是等间距采样直接用gradient会把沉积速率算偏这种情况下要用diff手动差分再按区间中心对齐。我的习惯是先用一行代码检查深度间隔% 检查深度是否等间距 if std(diff(d.depth)) 0.01 warning(深度不是等间距改用 diff 重算沉积速率); d.sed [diff(d.depth) ./ diff(d.age); NaN]; end为什么这个校验值得写死很多岩芯在感兴趣时段加密取样上部疏、下部密非均匀是常态。等间距假设一旦失效后面所有插值和平滑窗口的样本数都会出错峰值位置偏移几十厘米都是可能的。2.3 采样分辨率与沉积速率避免 CHAR 的“小样本抖动”CHAR 序列对采样分辨率极其敏感这也是 CharAnalysis 这类工具反复强调“先看沉积速率再谈峰值”的原因。沉积速率低的地方同样 0.5 cm 的取样间隔可能代表 50 年沉积速率高的地方则可能只代表 5 年。两者在同一个序列里共存时峰值识别的窗口参数不能一刀切。沉积速率 cm/yr采样间隔 cm等效时间分辨率 yr可用窗口建议0.010.550背景窗口 ≥ 500 yr0.050.510背景窗口 300-600 yr 可接受0.11.010注意峰值最小间隔要放大经验上一个可靠的火灾峰值至少要有一到两个样品落在峰值区间内。如果等效时间分辨率大于事件本身的持续时间峰会被抹平成一次普通起伏。遇到这种情况不要迷信插值加密结果因为插值不会凭空创造信息只会让曲线看起来细腻。常见做法是把 CHAR 线性插值到等时间间隔后再做背景分解比如统一到每 10 年一个点。插值本身不是问题问题在于插值步长不能小于原始分辨率的 1/2 到 1/3。把一个 50 年分辨率序列插成每年一个点产生的虚假尖峰基本无法靠后续平滑挽回。这个决策要在进入 CharAnalysis 流程之前完成并写进数据处理记录里否则最后审稿或复核时根本说不清峰值是真实事件还是插值产物。3. 诊断核心背景分解、峰值识别和 SNI 阈值怎么配合3.1 背景 CHARLOESS 平滑的窗口宽度选择CHAR 序列可以拆成两部分低频背景和高频峰值。背景代表区域性的火灾面貌和燃料变化趋势它不应该跟着单次火灾事件剧烈起伏峰值则代表局地事件。分离两者的常用工具是 LOESS 局部加权回归MATLAB 里对应smooth的rloess选项。注意这里不要用普通移动平均——移动平均会被峰值拉高导致背景曲线“穿透”事件最后峰值残差偏小。% 背景估计鲁棒 LOESS窗口按时间宽度换算成样本数 time_axis d.age; % 年龄序列 char_series d.char; % CHAR序列粒/cm2/yr win_yr 600; % 背景窗口常见取值 500-1000 yr dt median(diff(time_axis)); % 中值时间采样间隔yr win_n max(round(win_yr / dt), 5); % 至少保留5个样本 % rloess 是鲁棒局部回归能降低极端峰值的权重 background smooth(time_axis, char_series, win_n, rloess); resid char_series - background; % 残差用于峰值识别 % 第一时间出图诊断 figure; plot(time_axis, char_series, k-, LineWidth, 0.8); hold on; plot(time_axis, background, r-, LineWidth, 1.6); legend(CHAR, Background (rloess), Location, northwest); xlabel(Age (yr)); ylabel(CHAR (grains/cm2/yr));逻辑说明win_yr是核心参数它表示背景窗口的时间宽度。dt用中值采样间隔而不是平均值这样可以降低个别异常厚样品的影响。窗口换算成样本数后至少 5 个否则 LOESS 退化成逐点拟合背景会紧贴原始序列残差里几乎只剩噪声。参数说明win_yr取多少要结合序列整体长度。序列只有 2000 年时600 年窗口合理序列有 10000 年时可以放宽到 800-1000 年。判断背景是否合适的标准很简单背景曲线是否穿过了明显的峰。如果背景把峰也包进去了说明窗口太短如果背景几乎是一条直线忽略了几千年的区域趋势说明窗口太长。这个参数值得在 400、600、800、1000 几档之间各跑一遍把图叠在一起看。3.2 峰值识别全局阈值和最小间隔背景扣掉之后残差序列里正的部分是候选火灾信号。峰值识别不能只看“高”还要满足两个条件残差超过全局阈值并且峰值之间有时间间隔。全局阈值常见取正残差的 95% 或 99% 分位数这是数据自身定义的不需要人为设定绝对量纲。最小间隔常见取 75 到 100 年用来防止同一场火灾在多个相邻样品里被反复记为“峰”。% 峰值识别基于残差阈值和最小时间间隔 thr quantile(resid(resid 0), 0.95); % 95%分位数作为阈值 min_sep_yr 75; min_sep_n round(min_sep_yr / dt); % 找到所有残差从阈值下方上穿到上方的位置 cross diff([0; resid thr]) 0; candidate find(cross); % 对每个候选区间取区间内残差最大的样本作为峰 peak_idx zeros(size(candidate)); for i 1:numel(candidate) j candidate(i); while j numel(resid) resid(j 1) thr j j 1; end % 如果区间内有多个样本取最大值位置 [~, seg_shift] max(resid(candidate(i):j)); peak_idx(i) candidate(i) seg_shift - 1; end peak_idx(peak_idx 0) []; % 按最小时间间隔过滤两个峰太近时保留残差更大的 keep true(size(peak_idx)); for i 2:numel(peak_idx) if (time_axis(peak_idx(i)) - time_axis(peak_idx(i - 1))) min_sep_yr if resid(peak_idx(i)) resid(peak_idx(i - 1)) keep(i - 1) false; else keep(i) false; end end end peak_idx peak_idx(keep); peak_age time_axis(peak_idx); peak_resid resid(peak_idx); fprintf(识别到 %d 个候选峰过滤后保留 %d 个\n, ... sum(cross), numel(peak_idx));逻辑说明这版实现没有依赖 Signal Processing Toolbox 的findpeaks因为findpeaks的默认输出不完全等价于“残差上穿阈值”的定义。diff([0; resid thr]) 0是标准的上沿检测写法找到阈值被突破的起点。循环里把每个候选区间扫到末尾再用max定位区间内最高点这样即使一个峰跨三个样品也只取最突出的那个样本。参数说明min_sep_yr取 75 是古火灾研究里常见约定。对沉积速率低的序列这个间隔可以适当放大到 100 年对高沉积速率的湖泊75 年也未必足够因为一个事件可能在沉积里留下多层信号。阈值从 95% 改成 99% 会明显减少峰值数量适合那些残差尾部特别重的数据。改参数时记住一个原则宁可漏掉弱峰也不要产出大量无法解释的假峰。3.3 SNI诊断峰值质量的核心指标识别出峰不等于这些峰都能用。CharAnalysis 真正区别于普通找峰的地方在于它会给每个峰算一个信噪比指标 SNIsignal-to-noise index。SNI 计算的思路是峰值高度是信号残差序列中非峰值部分的离散程度是噪声两者的比值就是 SNI。经验阈值是 3SNI 大于 3 的峰被认为可靠小于 3 的峰在统计火灾频率时应被排除或单独标注。% 计算每个峰的SNI并筛选可靠峰 non_peak resid; skip unique([peak_idx; peak_idx - 1; peak_idx 1]); skip skip(skip 0 skip numel(resid)); non_peak(skip) []; % 把峰值附近样本剔除避免噪声估计被拉高 % 鲁棒标准差1.4826 * MAD noise 1.4826 * median(abs(non_peak - median(non_peak))); % SNI 峰值残差 / 噪声 sni_all resid(peak_idx) / noise; % 按SNI 3筛选 reliable peak_idx(sni_all 3); reliable_sni sni_all(sni_all 3); fprintf(可靠峰 %d 个平均SNI %.1f\n, ... numel(reliable), mean(reliable_sni));逻辑说明噪声估计必须剔除峰值样本否则峰值本身会把标准差拉大SNI 全被压低。剔除范围包括峰位前后各一个样本因为残差在峰周围都是高值。用中位数绝对偏差 MAD 代替标准差是鲁棒统计的常规做法系数 1.4826 让 MAD 在正态数据下等价于标准差。SNI 的含义用一张小表解释更直观SNI 区间含义处理建议 3信号弱可能是噪声抬升不纳入火灾频率统计3 - 5可识别但幅度有限保留但敏感性分析时要检查 5强峰事件信号明确正常纳入统计实际分析中SNI 的价值不只是筛峰它更是一种诊断如果一整条序列的 SNI 普遍很低说明问题不在找峰参数而在更上游的数据质量——可能是计数太少、采样太粗或者背景窗口和真实信号尺度不匹配。这时应该退回第 2 章重新检查输入而不是继续调阈值碰运气。4. 避坑CharAnalysis 落地的 5 个数据和参数陷阱4.1 峰值数量多到离谱现象一段只有 2000 年、明显低频变化的 CHAR 序列跑完识别出十几个峰平均不到 150 年就一个“火灾事件”和区域火史常识对不上。原因最常见是背景窗口太短。窗口短的时候LOESS 把短尺度波动视为背景趋势残差里混入了大量中频波动阈值只切掉最高的一小部分剩下的都成了峰。另一种可能是阈值取了 90% 分位数过于宽松。解决先把win_yr提到 600 到 800 年看背景曲线是否变平滑同时把阈值改为 95% 或 99% 分位数再把最小间隔加到 100 年。峰值数量应能减少一半以上。调整后对比原始曲线确保留下的峰确实对应肉眼可见的尖刺。4.2 晚全新世序列出现连续的 0现象数据前半段峰值正常到了最近几百年 CHAR 突然全变成 0背景和残差都是 0找峰结果为空。原因优先检查计数和体积列。如果最近段样品体积被记录成了 0 或空值conc就会变成 0 或 NaN。另一种情况是沉积速率被年代模型算成负值char出现负值或 0这是因为年龄模型在序列顶部存在倒转或过度平滑。解决逐样本打印四要素定位第一个 0 出现的位置。对照实验室原始记录确认体积是否漏填如果体积没问题就检查年龄模型在顶部的插值结果用线性插值替换平滑样条通常能消除负沉积速率。4.3 SNI 全线低于 3现象峰识别出一堆但 SNI 算出来没有一个是可靠的全部小于 3等于白做。原因噪声估计被污染。如果剔除峰值邻域的样本数太少或者残差序列本身很短MAD 会高估噪声。更根本的原因是采样分辨率不足一个真实的火灾峰在沉积记录里只有一个样品残差高度有限信号被噪音淹没。解决先扩大峰值剔除范围把前后各 2 个样本都拿掉再算噪声。SNI 还是低的话就不要试图在找峰参数上找补了老实承认数据不适合做事件级重建。这时候 CHAR 序列还能描述背景趋势但火灾频率和火返回间隔这类统计量没有实际意义。4.4 年代模型一换峰值位置和数量剧烈变化现象同一个岩芯用两种不同插值方式生成年龄列其他参数完全一致结果一个版本 5 个峰另一个版本 9 个峰峰位还发生了偏移。原因CHAR 计算里沉积速率直接来自年龄差分的倒数年龄模型对深度的斜率只要略有差异sed就会明显变化。尤其在深度-年龄曲线弯曲处线性插值和样条插值给出的沉积速率相差能到 30% 以上。解决固定一套年代模型作为分析基准不要中途混合。同时做一次敏感性分析把年龄列整体加减一个标准差范围重新跑整个流程看哪些峰在不同年代模型下都能保留这些峰才有资格进入最终结果。报告里要明确写清楚“基于哪种年代模型”否则结果不可复核。4.5 同一份输入跑两次结果不同现象排查问题时重跑脚本参数一模一样输出的峰值个数和 SNI 却对不上。原因多数情况下是随机数。背景平滑的鲁棒迭代或某些重采样步骤用了默认随机种子每次迭代初值不同结果就有波动。也有可能是表格读取时行顺序不稳定比如数据源来自多个文件合并。解决在脚本开头固定随机数种子rng(0)一行就够。如果固定后结果仍然不同检查数据读取部分是否做了排序确保每次输入行序一致。记录参数时也把 MATLAB 版本和所用工具箱版本写进去方便复现。5. 模拟数据自检与批处理把诊断工具用成可复用管线5.1 合成数据验证管线已知火灾事件能否抓回把 CharAnalysis 流程写到什么程度才算放心我的习惯是先造一段已知答案的合成 CHAR 序列设计一个低频背景叠加上几个时间位置已知的脉冲峰再混入随机噪声。然后跑完整流程看识别出的峰位置和预设位置是否对得上。这一步能在几分钟内暴露参数问题比拿真实数据反复试高效得多。% 合成数据自检生成1000年序列每10年一个点 t (0:10:1000); background 2 0.3 * sin(2 * pi * t / 300); % 低频背景趋势 % 在 t300 和 t700 处加入两个已知火灾峰 signal background; signal(30) signal(30) 8; % 300 yr 处 signal(70) signal(70) 5; % 700 yr 处 noise 0.6 * randn(size(t)); char_sim max(signal noise, 0); % 用自编流程识别峰值 d_sim.age t; d_sim.char char_sim; % 假设背景窗口300年阈值95%最小间隔75年 % ... 调用第3章的平滑和找峰代码 ... % 最后对比 peak_age 是否落在 300 和 700 附近逻辑说明合成数据里背景和峰的间隔是已知的所以能直接检验找峰逻辑有没有系统性偏差。signal(30)对应第 30 个样本因为时间步长 10 年所以就是 300 年同理signal(70)是 700 年。随机噪声用randn生成水平控制在峰高的 10% 左右模拟较理想的高分辨率沉积。参数说明模拟时建议把噪声逐步加大比如从0.2一直试到1.5观察峰在什么噪声水平下开始识别不出来。这样能摸清这套流程对数据质量的底线。真实数据投入使用前我会把这段自检脚本固化下来每次调参数都先过一遍再跑真实岩芯。5.2 把分析封装成函数批处理多根岩芯实验室里往往不是一根岩芯而是十几根一起分析。把第 2、3 章流程封装成函数输入统一表格输出每根岩芯的峰表和统计量批处理才能落地。函数化之后参数集中在结构体里而不是散落在各脚本中。function out run_charanalysis(df, param) % df: table至少包含 depth_cm, age_yr, volume_cm3, count_grains % param: struct含 win_yr, thr_quantile, min_sep_yr d struct(); d.depth df.depth_cm; d.age df.age_yr; d.volume df.volume_cm3; d.count df.count_grains; d.conc d.count ./ d.volume; d.sed gradient(d.depth) ./ gradient(d.age); d.char d.conc .* d.sed; % 背景分解 dt median(diff(d.age)); win_n max(round(param.win_yr / dt), 5); bg smooth(d.age, d.char, win_n, rloess); resid d.char - bg; % 峰值识别和SNI thr quantile(resid(resid 0), param.thr_quantile); % ... 第3章找峰和SNI代码此处省略重复部分 ... out.peak_age peak_age; out.sni reliable_sni; out.background bg; out.char d.char; end逻辑说明函数入参是统一表格和一个参数结构体返回值里同时包含了峰表和用于诊断的背景序列。这样既能输出最终事件表也能随时重新画出原始序列、背景、残差和峰位的叠加图。调用时只要写一句parfor i 1:numel(sites) run_charanalysis(data{i}, param); end就能对多根岩芯批量执行。参数说明批量处理时最容易出问题的是“一套参数打天下”。不同岩芯的沉积速率差异大win_yr可以统一但min_sep_yr和插值步长必须结合每根岩芯的时间分辨率单独检查。我在批处理里会为每个站点生成一页诊断图先扫一遍图再合并结果表不要直接相信峰表。做完整套流程后我现在的习惯是拿到任何新数据都先造一段合成序列自检再看诊断图最后才把峰值结果写进记录里。这个习惯帮我挡住过好几波假事件也省去了很多和审稿人解释的时间。希望帮到你。本文还有配套的精品资源点击获取