沉积物炭分析诊断工具CharAnalysis：原理、参数设置与排错指南

简介面向从事古生态与火灾历史重建的研究者CharAnalysis是一款用于沉积物-炭记录的诊断与分析的开源MATLAB工具。它通过峰值检测辅助判断局部火灾事件并提供参数敏感性测试、背景炭水平估算、置信区间等诊断手段帮助判断预处理与峰阈值设定的合理性。资源包共35个文件其中24个m文件构成本次分析的核心算法包含数据预处理、平滑、峰识别、结果汇总等模块5个csv及2个xls文件为示例数据与参数模板2个md与1个pdf分别提供更新说明和详细使用指南另有1个html便于快速浏览整体压缩包仅1.54MB。目前已有354人学习下载。借助良好注释的完整源码与配套文档使用者既能按默认流程输出统计图表也可深入理解算法逻辑并修改参数适合科研新手与进阶用户灵活应用。1. 沉积物炭分析的诊断工具CharAnalysis 到底解决什么问题你用重力采样器取回岩芯按厘米间距挑出木炭碎片Excel 里画出来是一根毛刺很多的曲线。你不能确定哪几个波峰是局地火灾哪些只是沉积噪声。CharAnalysis 是一套基于 MATLAB 的沉积物炭分析诊断工具能把木炭浓度序列拆成背景与峰值识别局地火灾事件并输出 SNI 信噪比指数告诉你这套数据值不值得信。这篇笔记会围绕它的原理、数据整理、参数设置和排错展开适合手里有岩芯炭数据、正准备做火干扰重建的研究者也适合想给现有成果找一套可对照分析流程的团队。2. 从木炭序列到火灾信号CharAnalysis 的核心原理与选型理由2.1 木炭浓度序列为什么不能直接设阈值很多初学者拿到炭浓度曲线后先找一个“看起来像背景”的基线然后超过两倍就认为是一次火灾。这种做法在数据量大、样品分辨率好、沉积速率稳定时勉强能用但大多数岩芯数据并不具备这些条件。木炭碎片的沉积是一个随机计数过程序列的方差往往与均值相关浓度高的深度段波动也更剧烈。直接设固定阈值高浓度段会被切碎低浓度段的局地火灾又会被漏掉。更麻烦的是浓度单位带来的陷阱。同样一批炭碎片如果沉积速率慢单位时间沉积的层位更薄炭碎片在单位体积内会被“浓缩”浓度数值虚高。只有把浓度换算成炭通量也就是每平方厘米每年沉积多少颗粒才能让数据在时间轴上可比。CharAnalysis 的处理起点就是围绕通量展开的先结合年龄模型把深度换算成时间再把浓度与沉积速率相乘得到通量序列。之后的信号模型可以写成观测到的炭通量等于一个低频背景组分加上高频峰值组分再加上噪声。低频背景刻画的是该地区长期、大尺度的炭背景输入比如风运或者坡面搬运造成的平均炭通量高频峰值代表采样点在某个时段内发生了局部燃烧事件炭输入显著超过背景。检测出峰值之后再结合年龄模型把深度换算成时间就能得到事件清单、事件频率和每次事件的火旋回长度。这里有一点常被忽略背景组分必须用局部拟合而不是全局求平均。因为湖泊或者泥炭剖面的沉积速率会在不同深度发生变化气候阶段的变化也会让背景水平缓慢漂移。直接用整段岩芯的平均值当背景深层和浅层的信号都会被扭曲。所以 CharAnalysis 会把整个序列划分成多个局部窗口在每个窗口内独立估计背景值窗口与窗口之间通过重叠和平滑衔接起来。2.2 核心诊断SNI 指数与背景拟合SNI 是这套工具最常被引用的一个输出全称是信噪比指数用来度量峰值检测是否可靠。它的思路并不复杂把序列分离完背景之后取残余的“信号”强度与“噪声”强度作比较。如果噪声水平过高说明样品间距太大、计数太少或者沉积速率过低峰值辨识只能停留在猜测层面。常见的经验阈值是 3高于 3 说明信号可以分离低于 3 建议加密采样或者合并层位而不是强行解释检测出来的事件。背景拟合的窗口宽度会直接影响 SNI。窗口太宽背景曲线过于僵硬残差里会留下大尺度的低频波动这些波动被误判成峰值SNI 反而被拉低窗口太窄背景会追着峰值跑真实火灾信号被吸收进背景SNI 也会下降。实际使用中需要先预估研究区的火旋回长度再把背景窗口设置为火旋回的 1.5 到 2 倍。北方森林火旋回普遍在几十年到两百年窗口可以选择 300 到 500 年湿润热带或高寒区域的火旋回更长窗口要适当放宽到 800 年以上。SNI 并不是一个恒定指标它随插值分辨率、背景窗口和样品密度变化。也就是说同样一条岩芯数据参数不同SNI 可能从 2 跳到 5。因此SNI 更应该被理解为一个诊断工具而不是验证真理的标准。当你改变插值分辨率后第一时间要做的事情不是看峰多了还是少了而是看 SNI 是否越过 3 这条线。2.3 这套工具的边界与选型理由为什么选 CharAnalysis 而不是自己在 MATLAB 里写一组处理函数最现实的理由是同行认可。目前古火研究里大家普遍在同一套定义下比较火干扰频率和火旋回长度发表出去的图和多站点汇总集成都容易复现。工具的输入文件简单直接一列深度、一列炭浓度再加上一份年龄模型。输出结果也都整理成文本和图形方便直接进论文图表。但边界同样明显。这套工具只做时间序列诊断不做炭屑形态鉴定也不能还原火强度。如果你的研究目标是区分草本炭和木本炭或者判断燃烧温度那应当回到显微镜鉴定而不是拿统计工具打天下。选型上还要提前确认数据形态炭数据需要转换成通量因此年龄模型是必需品。如果只有深度和浓度没有可靠的年代框架工具跑出来的事件清单都落在不可靠的时间轴上后续讨论没有意义。另外要注意运行环境。这套工具是用 MATLAB 写的老版本迭代时间较长某些新版 MATLAB 可能会提示旧式图形接口不可用或者缺少编译器属于正常现象。常见做法是把工具目录加进 MATLAB 路径遇到编译类问题就安装一个与系统匹配的编译器包不必为了实现一次分析去重装整个环境。3. 把数据喂给 CharAnalysis输入格式准备与最小跑通流程3.1 从岩芯记录到 .txt/.csv数据整理的规范动作CharAnalysis 需要一个基础的文本表第一列是深度第二列是炭浓度单位建议是碎片数每立方厘米。除此之外还需要一个年龄模型文件常见字段是深度与中位年龄。下面这段 MATLAB 代码演示如何把一份乱序、含缺失值的原始记录整理成工具可读的输入。% core_data.csv 原始数据depth_cm, concentration, sample_volume_cm3 data readmatrix(core_data.csv); depth_raw data(:,1); conc_raw data(:,2); vol_raw data(:,3); % 按深度升序排列岩芯采样常按由表到底记录也可能反着写 [depth, idx] sort(depth_raw); conc conc_raw(idx); vol vol_raw(idx); % 浓度统一为每立方厘米颗粒数 if all(vol 0) conc conc ./ vol; end % 删除缺失深度处的空值 valid isfinite(depth) isfinite(conc) (conc 0); depth depth(valid); conc conc(valid); % 用线性插值补掉个别空洞 conc fillmissing(conc, linear); % 写出标准输入文件 output table(depth, conc, VariableNames, {depth_cm, conc_part_cm3}); writetable(output, char_input.csv);逻辑说明排序和缺失值清洗是最容易忽略的一步。岩芯样品在实验室经过袋装、标签、称量后深度记录经常出现错乱或重复直接拿去算年龄模型会产生错误的时间轴。先 sort 再 fillmissing能保证后续插值的结果在深度上严格单调。参数说明重点关注 sample_volume_cm3这个字段决定浓度是按体积还是按重量折算。有些研究者习惯先记录“每克干重含多少碎片”这时候就不能除以体积而要除以样品干重得到每克炭浓度。两种单位对应的通量公式不一样后期换算火干扰频率时需要留意。如果你只有干重和碎片数那就在代码里把 vol 换成样品质量。3.2 结合年龄模型把深度转换为时间单独一份炭浓度文件不足以启动分析。需要把深度转换成时间工具才能计算炭通量和火旋回长度。年龄模型通常由实验室提供常见格式是深度与校正后年龄的对应表。下面是整合年龄模型的一个简单过程。% age_model.csvdepth_cm, median_age_cal_yr age_data readmatrix(age_model.csv); age_depth age_data(:,1); age_median age_data(:,2); % 对炭数据每个深度插值出年龄 age_interp interp1(age_depth, age_median, depth, linear, extrap); % 计算相邻采样深度之间的沉积速率 g/(cm2*yr) 或 cm/yr sed_rate diff(depth) ./ diff(age_interp); % cm per year % 如果年龄模型与炭数据深度重叠不足检查结果 if any(isnan(age_interp)) warning(存在超出年龄模型范围的深度请检查岩芯底部。); end逻辑说明interp1 默认只做已知范围内的插值超出部分需要显式声明 extrap否则返回 NaN。底层采样深度经常超过年龄模型覆盖范围特别是岩芯底部年代缺失时这里宁可直接报警也不要让工具用 NaN 往下传播。参数说明sed_rate 的单位会影响后续通量计算建议在输入文件里固定写成 cm/yr后续乘上炭浓度得到每平方厘米每年颗粒数。如果你拿到的是每克浓度要先将炭浓度乘以干密度再转换成通量这步少做数据量级会差出几个数。3.3 运行 CharAnalysis 的最小流程常见做法是把工具目录放进 MATLAB 路径然后在命令窗口运行入口脚本。入口脚本一般不接受直接传参而是弹出一个文件选择窗口让你依次选择炭浓度文件、年龄模型文件和分析参数文件。我一般会先在当前工作目录放好 char_input.csv 和 age_model.csv然后运行入口脚本。% 加入工具目录指定你解压后的绝对路径 addpath(D:\tools\CharAnalysis); % 运行入口脚本工具会读取当前目录下的标准文件名 CharAnalysis;逻辑说明addpath 是 MATLAB 查找函数的机制不加会报“未定义函数或变量”。运行入口脚本后如果代码按文件名自动查找最常见的位置就是当前工作目录所以我会先把两个输入文件放到当前目录。参数说明如果不希望每次弹窗可以在脚本目录下寻找一个名为 parameters 的模板文件许多版本用默认的文本文件来读取插值分辨率和阈值。具体键名各版本略有差异但内容大同小异。第一次跑建议全部保持模板默认值先把流程贯通再回来看参数怎么调。3.4 第一次运行后的输出怎么读跑完之后工具会在输出目录生成一组文本和图形。典型产物包括输出文件内容用途interpolated series插值后的炭通量序列查看插值效果background背景组分检查是否平滑peak series去除背景后的峰值序列识别局地火灾fire events局地火灾事件表事件清单含年份和峰值强度SNI tableSNI 值判断数据可靠性cumulative plots累积图直观对比背景与峰值拿到文件后第一件事不是数事件个数而是看 SNI 表。如果整段序列的 SNI 中位数低于 3说明数据质量偏低这时的峰值事件清单应当谨慎使用优先回到 3.1 节检查样品是否足够密集、是否做过体积归一化。这也是诊断工具的核心用途之一先跑诊断模式再跑正式分析。4. 参数到底怎么调插值分辨率、峰值窗口和 SNI 阈值4.1 插值分辨率先看沉积速率再做选择插值分辨率指数据在时间轴上的采样间隔单位是年每次采样。这个参数决定了后续所有计算的粒度。太细样品之间全靠数学插值会产生大量不真实的峰值太粗局地火灾信息会被平均掉。合理选择依赖沉积速率取深度区间除以年龄区间得到平均沉积速率再让插值间隔接近实际采样间隔。比如一段岩芯 50 cm 代表 1000 年平均每厘米对应 20 年样品按每厘米取样那么时间分辨率大约 20 年插值分辨率可以设成 10 或 20 年。我个人的经验是宁可略微细化不要过于加粗。细化后用 SNI 去做筛选如果 SNI 过低再合并。与其一开始就用粗分辨率牺牲事件不如保留高频细节让诊断指标替你把关。不过要小心插值分辨率设成 1 年而实际样品平均间隔是 50 年产出的光滑曲线会在事件位置产生大量假峰这类结果几乎没法用。4.2 峰值窗口怎么定局地事件与区域炭库的分离尺度峰值检测是在去除背景之后的残差序列上进行的。背景窗口过短部分火灾信号融入背景峰值被“吸走”背景窗口过长背景曲线又太僵硬残差里会留下大尺度的低频波动这些波动在后续检测中被当成多个事件。这个参数一般以年为单位常见设置在几百到上千年具体取决于研究区域的火灾发生周期。调整顺序很重要先定插值分辨率再定背景窗口最后管 SNI 阈值。背景窗口的初始值可以按预估火旋回长度来如果研究区历史文献记录显示大约每 80 到 120 年发生一次局地火灾那背景窗口设成 300 到 400 年是合理的。跑完以后看残差图如果事件峰周围还有明显隆起说明窗口可能还是偏窄如果背景曲线把峰顶都削平了就需要暂时排除峰值样本后重新拟合背景。4.3 SNI 阈值与显著性什么时候能信这套结果SNI 是决定“是否继续使用这套数据”的闸门。推荐起始值是 3它大致对应一组 99% 置信界限。低于 3意味着采样密度不足以分辨局地事件。此时不是简单地调低阈值而是回到数据本身把相邻样品合并到更宽的深度间隔重新计数再跑一次。很多研究者在低沉积速率岩芯里通过合并样品把 SNI 从 2 拉到 4这是提升数据质量的标准做法。SNI 的另一个用途是筛查沉积物的受扰动层位。如果一个深度段里炭浓度明显异常可能是浊流或再堆积的结果这一段的 SNI 会显著偏低。遇到这种情况我会把该段在输入表里剔除而不是让异常段参与背景拟合。工具并不替你做沉积学判断它只负责用统计方式告诉你哪里不对劲。参数建议起始值调整方向插值分辨率实际采样间隔的中位数SNI 太高时可粗峰值过少时可细化背景窗口预估火旋回的 1.5 到 2 倍峰粘连时加宽背景穿峰时收窄峰值检测窗口与背景窗口配合局部阈值在低浓度段更灵敏SNI 阈值3低于 3 优先优化数据而非降阈值5. 沉积物炭分析避坑指南五个典型问题和排错思路5.1 深度序列混乱导致数组越界现象MATLAB 报 Index exceeds array bounds或者图形只出一半。原因原始表里深度不是单调排列排序被跳过插值后长度错位。很多实验室记录表格里会把补采样品插在中间结果深度顺序被打乱。解决先排序再清洗。用第 3.1 节的代码段处理一遍确保 depth 严格递增。如果岩芯是倒序记录的不排序也能跑但年龄模型的时间轴方向是反的后面的事件年份会全错。排完序后顺便看一眼相邻深度差值有没有异常为零的重复值。5.2 峰值检测结果为空图上找不到事件现象fire events 表是空的或者只有一行。原因插值分辨率太粗短命的火灾信号被平均到背景里或者背景窗口太窄把峰值吸进背景。解决将插值分辨率从 20 年调到 10 年同时把背景窗口从 500 年放宽到 750 年再跑一次。如果仍然为空回到原始浓度序列看看最大峰值是否超过背景一倍以上。本身数据没有强信号时工具不是魔法不能凭空造出火事件。这时候要考虑是不是岩芯取自一个从未发生过局地火的沉积环境。5.3 背景曲线跟着峰值走残差里没有东西现象背景曲线在每个峰值处都隆起残差非常平事件表里什么都没剩下。原因背景拟合窗口太短或者背景估计时把峰值点带了进去。背景拟合的目标是刻画长期趋势如果窗口内正好包含一个高两倍的峰值局部背景会被抬高。解决加长背景窗口或者开启“峰值掩膜”选项。很多版本支持先粗检峰值再从背景拟合中排除峰值样本这一步很关键。如果工具不支持也可以手动把明显的峰值样本在输入表里打个标签但不建议手动改动最终事件表。一个可操作的变通做法是先跑一次获得粗事件表然后只对非事件样本做背景拟合把这个平滑结果作为背景再跑一次。5.4 SNI 低到没法用数据质量被审稿人质疑现象输出表里 SNI 在 1.5 上下峰值却不少画出来的图很好看但审稿意见会问 SNR 是否足够。原因样品量太少低计数区的泊松噪声会淹没信号最典型的场景是沉积速率极慢的剖面每一厘米代表 100 年以上。解决合并相邻样品。把每 2 cm 或每 5 cm 的炭碎片数和体积累加重新计算浓度。这样噪声减小SNI 通常能抬升到 3 以上。代价是时间分辨率变低但至少得到一个可量化的结果。合并后如果 SNI 依然很低可以考虑放弃局地火灾重建只报告背景通量的长期变化趋势。5.5 同一根岩芯两组参数得到两份不同结果现象某研究者用默认参数跑出一个 8 次事件的火旋回长度另一个人微调参数后变成 5 次两版结果对不上论文没法写。原因参数记录不全或背景窗口和峰检测阈值不一致。这是这类工具最典型的翻车现场问题不在算法而在人。解决确立固定参数规范。我会先把插值分辨率定为样品间隔中位数背景窗口定为预估火旋回的 1.5 倍SNI 阈值固定为 3并把这份规范连同工具版本信息写进论文的补充材料。任何人复现这个结果只需要同样的输入表和同样的参数文件。参数每次改动都在文件名后加后缀比如 site A_v1_res10_win500一劳永逸。6. 用残差分析和交叉验证收尾怎样判定你的火干扰重建可信6.1 残差图是最快的质量检查每次跑完除了看事件表我会先拉一张残差图横轴是年份纵轴是去除背景和峰值之后的残留噪声。理想状态是残差均匀分布在一个窄带里没有明显趋势。如果某个时段残差持续偏高说明那里有未拆干净的信号事件表可能漏了事件如果残差在低浓度区持续偏低说明阈值可能太宽把噪声也选成了事件。残差图还能反映年龄模型的问题。如果在某个深度段残差出现突变而炭数据本身没有异常那大概率是年龄模型在该段出现拐点导致通量计算方式变化。这时候回去检查年龄模型而不是急着调参数。6.2 随机化检验可以增加可信度更严格的验证是用随机化方法。把每年的事件随机重排后重新计算火旋回长度与真实序列做对比看观测值落在随机分布的位置。这个过程不需要额外数据把年代-浓度序列复制成多条随机序列跑同样的流程即可成本低但说服力强。具体操作上我会在 MATLAB 里生成 1000 组随机序列保持总体炭通量分布不变只打乱时间顺序然后逐一提交给 CharAnalysis 的事件检测流程得到 1000 个“随机火旋回长度”。如果真实序列的火旋回短于 95% 的随机结果可以认为局地火灾具有显著的时间聚集性。这一步能让审稿人对你的重建更有信心。6.3 我个人的一个小习惯我现在每次跑之前都会先写一个 parameter_set.txt 存下来文件名带上岩芯编号和日期。这个习惯救过我两次一次是审稿人要参数一次是半年后自己忘了当时怎么调出来的结果。参数文件写好后先在诊断模式里跑一遍确认 SNI 过了再启动正式分析。宁可多花 10 分钟记录也不要在结果被质疑的时候翻车。希望帮到你。本文还有配套的精品资源点击获取