从一块岩石到一个年龄

LA-ICP-MS 锆石 U-Pb 定年：物理原理、每一步算法、以及结果该怎么读

面向第一次接触这套流程的人 · 对应 DRUID 2.8.0 与 ADEPT 1.3.0 · 所有公式与数值均取自源码

一、时间是怎么被"读"出来的

1.1 一块自己会走快的表

想象一个沙漏，但它的沙子不是匀速漏的：剩下的沙子越多，漏得越快。每一粒沙子的漏出概率是固定的，所以一开始沙子多、漏得猛，越到后面越慢，永远漏不完，只是越来越接近空。

这正是放射性衰变的行为。一个²³⁸U 原子在任何时刻衰变的概率是固定的，不受它已经等了多久、周围有什么影响——这叫恒定衰变律，是这套定年方法的全部物理基础。写成数学，母体原子数随时间的变化是：

其中 λ 叫衰变常数，单位是"每年"的倒数。它不是理论算出来的，是实验室反复测出来的物理常数。²³⁸U 的 λ 是 1.55125×10⁻¹⁰ /年（Jaffey 等，1971）。这个数字小得出奇——意思是：一个 ²³⁸U 原子平均要等 64 亿年才衰变一次。

把一个 λ 换算成更直观的"半衰期"（母体剩一半需要的时间，T½ = ln2/λ）：

1.2 为什么必须用"比值"

上面第二个公式看起来可以直接用：量出样品里积累了多少 ²⁰⁶Pb，除以母体 ²³⁸U，再反解 t 就行了。但有一个问题——你永远无法知道"初始有多少 ²³⁸U"。四十四亿年前这个晶体开始生长时有多少铀，已经无从考证了。

所以我们换一个角度。²³⁸U 少了多少，就等于 ²⁰⁶Pb 多了多少（假设没有铅跑出去）。设当前铀的原子数是 N：

N₀ 被消掉了，只剩下当前测得到的两个量之比。这就是为什么 U-Pb 定年测的是比值而不是含量。

这个技巧还有一个额外的好处：质谱仪对"某个元素测得准不准"依赖很多不稳定因素（灵敏度漂移、基体效应）。但如果你测的是同一个元素的两个同位素之比（比如 ²⁰⁷Pb/²⁰⁶Pb），那些不稳定的因素大多会同时乘在分子和分母上，相除就抵消了。同位素比值比绝对含量稳健得多——这是贯穿整套流程的一个主题，第 3 章还会反复遇到。

1.3 为什么锆石是最好用的那块表

理论上任何含铀矿物都能当时钟。实践中锆石（ZrSiO₄）几乎是唯一的选择，原因有三个，每一个都对应一个具体的麻烦：

• 它含铀，而且含量适中。通常几百到几千 ppm。太少则信号弱、测不准；太多则辐射损伤会让晶体变成"蜕晶质"（metamict），铅就容易跑掉。
• 它几乎不含铅，但很容易容纳铀。这是关键。晶体结构（配位）排斥半径较大的 Pb²⁺，却接受 U⁴⁺（U⁴⁺ 与 Zr⁴⁺ 半径接近，可以替代）。所以晶体形成时"钟的读数是零"这个前提基本成立。如果矿物一开始就带着一堆铅（叫"普通铅"），这个时钟就不是从零开始走的——第 3.4 节几乎全部在对付这件事。
• 它非常结实。熔点高、化学性质惰性，能扛过搬运、风化、甚至一定程度的变质。所以岩浆里结晶出的锆石，很可能原封不动地保存到今天。

最后一条还有另一个侧面，它既是优点也是麻烦：因为结实，一个锆石晶体可以记录多次事件。如果岩浆在几百万年甚至几亿年后又重新活动，新的锆石会在旧晶体的外面继续生长一层。于是你得到一颗有"核"和"边"的晶体，它们的年龄不同。

1.4 一个晶体里的两个年龄

传统做法是"躲"：在显微镜下挑出均一的区域，把激光对准那里。但一个 20 μm 的坑对一个只有几微米宽的过渡带来说太粗了，你没法用眼睛保证躲开。更常见的结局是——你以为打的是边，其实打进去了核，两条年龄混在一起，最后得到一个"协和度很差"的数字（第 4.4 节讲这个指标），然后整批数据被丢掉。

另一种做法是"用"。既然剥蚀过程中连续测了三十秒、每秒十次，为什么不看它随深度的变化？如果剖面在一段时间里稳定在 240 Ma，然后干净地跳到 400 Ma，两个年龄就都在数据里，只是需要把它们分别提取出来。这就是 DRUID（提取逐深度窗口的年龄剖面）和 ADEPT（在剖面上自动找出这些"稳定的段"）的分工。

这套思路成立的前提是：信号真的能分辨出深度的变化。下一章讲这个前提在物理上有多牢靠。

二、数据是怎么产生的：一台质谱仪里的七十五秒

2.1 从激光到计数

要读时钟，得先把晶体的一小块"搬进"仪器。这一步由激光完成，之后的每一步都是在把这团物质转换成电信号：

1. 激光剥蚀。一束聚焦的紫外激光打在锆石表面，瞬间把几十微米范围内的物质加热到汽化。被剥蚀下来的物质不是气体，而是极细的气溶胶——大量亚微米到微米级的颗粒。
2. 载气搬运。氦气（有时是氦氩混合）从样品池里把气溶胶吹出来，沿着管路送进质谱仪。这段管路有几十厘米到一米长，气溶胶在里面的传输效率本身就是一个会漂移的量。
3. 等离子体电离。气溶胶进入氩等离子体炬——温度大约 8000 K，比太阳表面还热。颗粒在这里被完全分解、原子化，然后被电离成带电离子。这一步的效率对不同的元素是不同的，叫"元素分馏"。
4. 质谱分离与计数。离子被电场加速、通过质量分析器（按质量/电荷比分离），最后打到检测器上。检测器每收到一个离子就输出一个电脉冲，电子学部分数这些脉冲。

输出就是每个质量数、每个时间切片上收到了多少个脉冲。仪器通常把它换算成 cps（counts per second）——每秒计数。

2.2 一次测点的时间序列

2.3 cps 不是"计数"——这件事很要紧

看起来把一个时间段的 cps 加起来就得到"总计数"，噪声就是 √N。但这个想法有两个漏洞，两个都会在后面的不确定度里出现：

漏洞一：cps 已经除过积分时间，且做了死时间校正

检测器在高计数率下会"数漏"——两个离子在极短时间内先后到达，电子学会把它们算成一个。这叫死时间效应，计数率越高漏得越多，是典型的系统性偏差。校正公式（非延长型）是：

本实验室实测值是 14.7645 ns（记在 druid/deadtime_ns.txt），但这个数字默认不启用——因为一旦默认打开，历史上所有未校正的结果就不可比了。要用得显式传参。

漏洞二：真实噪声不是泊松的

泊松统计假设"离子到达是完全独立的随机事件"。但 ICP-MS 上还有一种占主导的噪声源：等离子体闪烁噪声（flicker noise，也叫 1/f 噪声）。它来自等离子体本身的温度与密度起伏，其大小与信号强度成正比，而不是与 √信号成正比。

后果很实际：把积分时间加长，泊松噪声按 √N 下降，但闪烁噪声不怎么下降。所以"多点几个点就更准"在某个点之后就没有收益了。DRUID 的做法是干脆不用解析公式，改用刀切法（jackknife）从数据本身估计离散度——它对两种噪声一视同仁（reduction/ratios.py，第 3.8 节细讲）。

2.4 八个通道分别在测什么

一次测点通常同时采集下面这些质量数。它们的角色很不一样，值得记住：

2.5 从"一串数字"到"一个年龄"，中间隔着什么

到这里我们有了：一次测点、七个通道、几百个时间点、每个点一个 cps 值。目标是：这一点的年龄，以及它的不确定度。

距离还很远。接下来要处理的，是下面这些每一件都会单独污染结果的东西：

下一章按这个顺序逐个讲。每一步都会给出为什么（物理原因）和怎么做（实际公式），并且在必要时说明"这一步做错了会导致什么方向的偏差"。

三、从 cps 到比值：八个必须做的校正

下面八步按实际执行顺序排列。每一步都回答同一个问题：有没有什么东西正在污染这个比值？有些步骤是必须的，有些是可选的（比如死时间），但即使在"可选"的地方，也说明了为什么默认不启用——因为把一个校正关掉比错误地打开它更安全。

3.1 死时间校正（默认关闭）

检测器在高计数率下会漏计数，前面 2.3 节给了公式。这里补充三件事：

• 为什么是系统性的而不是随机的：漏计数只发生在高计数率上，所以它总是把信号往低的方向拉，不会两边对称。这是偏差，不是噪声。
• 为什么对 Pb/U 影响不对称：²³⁸U 的计数率通常远高于 ²⁰⁶Pb，所以 U 被压低的比例更大。结果是 Pb/U 比值被系统性抬高，年龄偏老。
• 为什么默认关闭：死时间是每台仪器、每个检测器自己的参数，必须用标准物质单独标定。本实验室实测 14.7645 ns，记在 druid/deadtime_ns.txt。但它不是默认值——一旦默认开启，历史上所有未校正的结果就不再可比了。要用就显式传参：

python -m druid.cli.reduce_batch --dir 某批次 --deadtime-ns 14.7645

3.2 找到剥蚀区间：为什么要两级阈值

目标很简单：找出激光真正打开的那一段。做法是拿 ²³⁸U 的信号（最强、形状最干净）去比对阈值。麻烦在于阈值的选取：

• 只用一个"本底 + k×噪声"的阈值不够。实测 iCAP RQ 在激光关闭后 10 秒，拖尾还有 100–200 cps——远高于本底噪声（本底可能只有几十 cps）。用一个统计阈值去卡它，会把窗口一直拉到记录末尾，把大量零信号混进来。
• 所以再加一个"峰高的百分比"阈值。本流程用峰值的 3%：thr2 = 本底 + 0.03 × (峰值 − 本底)。这个阈值不问"这是不是噪声"，只问"这是不是峰的一部分"，能干净地砍掉拖尾。
• 两个阈值取最大值，保证任何情况下阈值都不会落进本底噪声里（万一某个点信号极弱，百分比阈值会退化成接近本底的水平，那时靠统计阈值兜住）。

两个细节值得注意。第一，"峰值"用的是 98 分位数而不是最大值——最大值可能只是一次随机尖刺（比如一个包裹体被打爆），用它定阈值会让门槛虚高。第二，本底的 σ 用 MAD 估计而不是标准差（σ ≈ 1.4826 × MAD）。原因是本底段里可能混有零星尖峰，标准差会被它们放大，中位数绝对偏差不会。

信号在阈值附近的抖动会产生若干段短暂的"超阈值"区间，所以最后取最长的一段连续区间作为剥蚀段，而不是简单地取第一个超阈值点。

3.3 气体空白：为什么要取"紧挨着剥蚀之前"的那一段

激光关着时测到的信号，来自载气里的杂质（主要是氩气瓶和管路里的汞）、样品池的残留、以及仪器本身的电子本底。这部分必须在所有通道上逐通道扣除。

取哪一段？直觉上取文件开头最方便，但这里取的是剥蚀开始前 15 秒（并且留 0.5 秒的安全间隔，躲开激光实际起跳与电信号之间的延迟）。原因是：气体空白会缓慢漂移——氩气瓶里的汞含量、管路温度、进样系统的状态都在变。时间上离剥蚀最近的那一段，才代表剥蚀那几十秒里的真实本底水平。

扣除之后，净信号允许是负数。这不是错误，是统计学上正常的：信号恰好等于本底时，扣完就在零附近随机波动。这里绝对不能 clip 到零——那会把一半的随机涨落抹掉，引入一个系统性正偏差。ablation.py 的注释专门强调了这一点。

3.4 汞干扰与"这个铅是真的吗"——整条链路最关键的一步

这是整份文档里技术含量最高、也最容易做错的地方。它值得单独一节。

²⁰⁴Pb 是唯一能直接指示"普通铅"的同位素，而普通铅必须扣掉（下一节讲为什么）。但 ²⁰⁴ 这个质量数上坐着两位乘客：

Hg 的 ²⁰²/²⁰⁴ 丰度比 4.35 是自然界常数，而 ²⁰² 这个质量数上没有别的干扰，可以直接读。所以公式本身没问题。问题出在数量级上：

于是有了这套流程里我认为最重要的一个决定：先问这个 ²⁰⁴Pb 是不是真的存在，再决定要不要扣。

如果不管显著不显著，一律扣除，会发生什么？

后果是方向明确的，而且很严重。假设某个测点其实没有普通铅，但你按噪声估计出的一个正的 i204 去扣：

两项相除得到 ²⁰⁷Pb/²⁰⁶Pb。²⁰⁷Pb 被多扣了，所以这个比值偏小，而 ²⁰⁷Pb/²⁰⁶Pb 年龄随之大幅偏年轻。你会得到一批"看起来很精确、但系统性错了"的年龄，而且因为 ²⁰⁷Pb/²⁰⁶Pb 在协和度检验里占重要位置（4.4 节），这种偏差还能伪装成"数据不一致"，让你去怪样品。

3.5 普通铅校正：一个鸡生蛋问题的迭代解法

为什么要扣普通铅：晶体生长时可能裹进了一点点周围环境的铅（来自地壳或地幔的铅，不是铀衰变出来的）。这部分铅的存在会让"钟的读数是零"这个前提不成立，测出的年龄偏老。

要扣掉它，需要知道普通铅的同位素组成——即 ²⁰⁶Pb/²⁰⁴Pb、²⁰⁷Pb/²⁰⁴Pb、²⁰⁸Pb/²⁰⁴Pb 三个比值。而普通铅的组成取决于它形成时的年龄：地球的铅同位素组成随时间演化（²³⁸U 和 ²³⁵U 一直在衰变，所以古老的地球铅里 ²⁰⁶Pb/²⁰⁴Pb 更低）。这个演化模型叫 Stacey-Kramers 两阶段模型，可以查表得到任意年龄对应的铅组成。

于是问题来了：要算年龄，需要知道铅组成；而铅组成需要知道年龄。

这是标准的不动点迭代。初值取 500 Ma（随便一个中间年龄即可），实测 3–4 轮必然收敛，所以代码里迭代上限直接写死为 4 次（reduction/ratios.py）。如果某轮迭代出现了非物理的结果（比值非正），就放弃校正、走"不扣"的路线——又一次"宁可少扣不可错扣"。

3.6 滑动窗口：为什么窗口要重叠

到这里我们得到的是"整个剥蚀段一个年龄"。但要提取深度剖面，需要把它切成一串更细的窗口。

窗口参数不是随便定的。太小则单个窗口内只剩几个 cycle（4 s 才 13 个），各通道的计数统计都变差，先撑不住的是计数本来就少的 ²⁰⁴Pb；太大则把核-边的过渡带糊在一起，看不出结构。本批实测的合理值是 win = 4 s / step = 1 s——每个窗口里 12~13 个采样周期（cycle 周期实测 0.31871 s）。

3.7 剥蚀坑越打越深，比值会自己漂移

这是深度剖面方法必须解决的问题，否则整件事都不成立。

解决办法是用标样。标样的年龄是已知的（"推荐值"），所以在同一台仪器、同一次运行里，标样测出来的比值偏离推荐值多少，就反映了当下的分馏程度。关键在于要在同一个深度上比——因为分馏是深度的函数，不是常数。

3.8 不确定度：为什么不用 √N，为什么要两个来源相加

一个年龄数字如果不带不确定度，是没有意义的——因为"453.1 Ma"和"453.1 ± 8 Ma"是两种完全不同的结论。DRUID 给出的逐点不确定度由两部分正交相加：

为什么必须给出相关系数 ρ

最后还有一个容易被忽略的输出：²⁰⁶Pb/²³⁸U 与 ²⁰⁷Pb/²⁰⁶Pb 之间的相关系数 ρ（典型值 0.7–0.9）。它们之所以强相关，是因为两者共用同一个 ²⁰⁶Pb 测量值——²⁰⁶Pb 测高了，第一个比值偏大、第二个比值偏小。

这个 ρ 在算年龄的时候用不上（两个比值各自换算成年龄），但在画协和度图时必须用：数据点的误差椭圆倾角由 ρ 决定。忽略它，椭圆会被画成"正的"，视觉上完全失真，可能让人以为某个点偏离协和曲线，其实它只是沿着协和线方向有一大段误差。

刀切法（jackknife）在这里怎么用

做法很直接：窗口里有 n 次采样，就依次去掉第 j 次、用剩下 n−1 次重算一遍比值，得到 n 个"伪值"。这 n 个伪值的离散程度就是不确定度的估计：

为什么不用解析公式？因为刀切法不需要假设噪声的分布形式。前面说过真实噪声里混着泊松成分和闪烁成分，比例还随信号强度变化，解析式很难写对。刀切法直接从数据本身估计离散度，两类噪声都会被算进去。代价是计算量（每个窗口要重算 n 次），但这是值得的。

四、从比值到年龄：一个公式，和三套时钟

到这里，每一个窗口都有了两组数：比值和它的不确定度。剩下的数学反而简单了——因为这部分的物理只有一条定律，没有别的干扰项。

4.1 有解析解的两个：一行公式

把第 1.2 节那个式子反解，就是年龄：

两个实现细节值得提。第一，代码里用 log1p(r) 而不是 log(1+r)——当 r 很小时（年轻样品），这两个在数学上等价但在浮点数里不等价，log1p 专门为这个情况设计，能保住有效数字。

第二，r 可能是负数。扣完气体空白和普通铅之后，净信号是允许在零附近随机波动的（第 3.3 节）。如果 r ≤ −1，对数发散。代码里的处理是把 r 钳到 −0.999999 而不是返回 NaN——因为返回一个接近零的负年龄，比返回 NaN 要好得多：NaN 会在后面的加权平均里"传染"，把一个窗口的 NaN 变成一个坪的 NaN，而一个负年龄只会让那个点在平均里被稀释。这是"选择在哪一层承担错误"的典型例子。

4.2 没有解析解的那个：为什么用二分法

²⁰⁷Pb/²⁰⁶Pb 就麻烦了。它的正演公式是两个指数项之比：

要反解 t，需要解一个超越方程——写不出闭合形式。工程上有两条路：用牛顿法数值求根，或者用泰勒级数近似公式。这套代码选的是二分法，理由写得很实在：因为放射性累积是单调的，函数在区间上单调递增，二分法保证收敛；而牛顿法在某些区间会"跳飞"到无意义的解上去。在科研代码里，可预测的慢比不可预测的快更值钱。

实现上，先确定一个包含真解的区间 [1 Ma, 4567 Ma]（下界取 1 Ma 是因为方程在 t→0 附近退化，上界取地球年龄），然后每次取中点、判断正演比值偏大还是偏小、收缩区间。R 侧迭代 60 次，Python 侧 80 次——都远超精度需要（区间长度除以 2⁶⁰ 早已小于 1 年），多迭代没有收益。

4.3 三套时钟，各有各的擅长区间

为什么要有三套？因为它们对年龄的灵敏度分布完全不同。

4.4 协和度：三套时钟互相检验

同一个样品的三套年龄理论上应该相等——因为它们是同一个时钟的三个不同读数。实际上它们会有差异，而差异本身就是信息：

为什么会有不协和？最主要的原因是铅丢失：如果晶体后来受过热或变形，一部分放射成因的铅可能扩散丢失。因为不同同位素在不同矿物中的扩散行为不同，铅丢失对三套年龄的影响程度也不一样，于是三个数字就分开了。另一个原因是多期混合——这正是本方法要解决的（第 2 章），而它也会表现为不协和。第三个原因是普通铅扣得不干净（第 3.4 节），未扣净的普通铅会拉低 ²⁰⁷Pb/²⁰⁶Pb。

通常的判断标准是：协和度落在 90–110% 之内算协和，可以考虑使用；超出这个范围的点需要解释。

还有一张图在 U-Pb 文献里很常见，这里简单说明它的含义：协和曲线图（concordia diagram）。横轴是 ²⁰⁶Pb/²³⁸U、纵轴是 ²⁰⁷Pb/²³⁵U，把所有可能年龄的理论值连成一条曲线。协和的样品点落在这条曲线上；铅丢失的样品点因为两个比值按不同比例降低，会落在曲线下方。这也是为什么前面强调那个相关系数 ρ：数据点的误差椭圆倾角由它决定，画协和曲线图时必须用对。

五、从剖面到坪：在一条抖动的线上找台阶

到这里我们有了：每个测点一条"年龄随深度变化"的剖面，每个点带不确定度。剩下的问题是把这条线读成一到两个年龄。这件事看起来简单——人眼一眼就能看出台阶在哪——但要让程序每次都做出同样的判断，需要把它拆成明确的规则。

5.1 为什么不能整段平均

回到图 2 里的例子：一个坑从上到下穿过了 240 Ma 的边和 400 Ma 的核。如果对整段做一次积分平均，得到 330 Ma——这个数字不对应任何一次地质事件，它是两个人的平均身高。更糟的是它会带着一个看似合理的协和度，让人觉得"这是个有效数据"。用一个不存在的事件去讲一个故事，是这类数据最危险的用法。

所以要找的是剖面里那些稳定的水平段——在钻探行话里叫"坪"（plateau）。一个坪意味着一件事：在这段深度范围内，测到的年龄没有系统性变化，这通常对应晶体生长时的一个稳定阶段。

5.2 分段：PELT 在做什么

"找出水平段"可以写成一道最优化题。假设把序列切成 k 段，每段内部的"不平整程度"用段内残差平方和衡量（就是各点到本段均值的距离平方和）：

第一项让每段尽可能平，第二项惩罚"切太多段"。第二项不能少：没有它，最优解永远是"每个点自成一段"，因为那样每段的残差平方和恰好是零。这个加惩罚的想法和 BIC 准则同源——你想让模型更贴合数据，但要为复杂度付出代价。

求解这个问题的算法叫 PELT（Pruned Exact Linear Time，Killick、Fearnhead & Eckley，2012），ADEPT 里的实现来自 R/pelt.R。它的关键是"剪枝"：如果某个候选变点已经不可能成为最优解的一部分，就直接剔除，不必再考虑它。这让平均计算量从 O(n²) 降到接近 O(n)。

代价函数用的是前缀和技巧：要算"第 s 到第 t 个点的残差平方和"，不需要遍历这 t−s+1 个点，而是用两个预先算好的累加数组直接相减：

惩罚项取 (4/100) × log(n)。那个 4 不是随手写的：它来自 changepoint 包对"单变量均值模型"给出的 AIC 惩罚值（2 个参数 × 2）。ADEPT 把这个常数硬编码进来，就不必依赖那个第三方包的特定版本——因为包的版本一变，这个内部常数可能跟着变，而结果会跟着变。

5.3 分段之前要先平滑——以及一个危险的例外

直接对原始剖面跑 PELT 会出问题：噪声本身就有很多局部起伏，算法会把它们当成变点，于是分出大量没有意义的短段。所以标准流程是先用 LOESS（局部加权回归）把剖面磨平一点。

LOESS 的做法是：对每个点，用它周围一段范围内的点做加权线性回归，用拟合值代替原值。span = 0.15 意味着每个点用邻近 15% 的数据来拟合。这样保留了大尺度的趋势（真正的域边界），抹掉了小尺度的抖动（噪声）。

5.4 标准化：为了一个技术性的理由

平滑之后，剖面值会被线性缩放到 [0, 1] 区间（min-max 标准化）。这一步不是为了好看，而是因为 PELT 的惩罚项是一个固定数值：如果数据本身的尺度是几十秒的年龄差，惩罚太小不起作用；如果尺度是 0.01，惩罚又过大。把尺度固定下来，惩罚项才有可比的意义。

顺带说明代码里的一个细节：标准化用的是 loess 值的最小/最大值作为参照，而原始年龄列也用同一组参照值来缩放。这样两条曲线（平滑的和未平滑的）落在同一个尺度上，可以直接比较——这正是 smooth = "none" 时它们能完全重合的原因。

5.5 四步筛选：从"切出来的段"到"能用的坪"

PELT 会切出很多段，其中大部分是噪声或过渡带。接下来要过滤。过滤是有顺序的，每一步针对一类特定的"不该要的段"：

1. 年龄范围：坪年龄必须落在 [0, 4540] Ma 之内。这一步看起来多余（年龄怎么可能是负的？），但它会抓到前面步骤的真实错误——比如普通铅扣错、比值列搞反。超出这个范围的段一律丢弃。
2. 段内方差：把每段的值归一化到 [0,1] 之后算方差，阈值 0.1192。方差大的段是过渡带——它"跨越"两个域，本来就不该被当成坪。这个阈值是从原始发布的流程沿用的经验值，不是从数据推导出来的，属于可调参数。
3. 最短时长：一段的积分时长必须 ≥ 5 秒。太短的段虽然平，但点数太少、不确定度大，不足以支撑一个年龄。第 1 段豁免这个规则——它带着剥蚀起始的信号，常常被 lower_ablation_time 截断，如果一视同仁就永远会被丢掉。
4. 方向性：见下一节。这是四个步骤里唯一一个考虑了"年龄之间的相互关系"的。

5.6 方向性：从"全有或全无"到"保留最大的自洽子集"

锆石的深度剖面通常有一个方向性：如果坑从边打向核，年龄应该单调变化（一直增或一直减）。所以那些"和趋势相反"的段可能不可靠。

但这里有个真实的困难：噪声会让局部的顺序偶尔颠倒。三个坪的年龄是 453.0、453.2、453.1 Ma——这显然是一条上升趋势，只是中间那个点略微回落了一点。怎么处理？

这个改动是这套工具一个重要的演进：从"全有或全无"变成"保留最大的自洽子集"。算法上，求最长单调子序列是标准的两层循环动态规划（O(n²)，但 n 是段数、通常只有几十，代价可以忽略）。

为什么一定要保留方向性判断？因为一个剖面里出现的"台阶"有可能不是地质事件，而是仪器扰动。要求坪之间保持一致的年龄变化方向，是一种利用地质常识来排除随机噪声的手段。filter_direction 参数决定要"递增"（Forward）还是"递减"（Reverse），direction_tolerance（默认 2%）控制允许多大的反褶。

5.7 坪年龄怎么算：两种情形

最后一步是把坪内那些点合成一个年龄。这里有两条路径，取决于输入里有没有不确定度：

加权平均的直觉是：误差小、更可信的点应该有更大的话语权。权重取方差的倒数，这样一个 σ 小一半的点，权重会大四倍。理论上这也给了我们一个可以验证的性质：如果所有点的 σ 都一样，权重就均匀，加权平均必须退化成普通的算术平均。实测确实如此——同一批数据用 σ = 0.5 和 σ = 5（都是常数）算出来的坪年龄完全相同（21.80575，差在 1e-10 量级）。这种"退化测试"是验证加权实现有没有写错的最直接办法。

六、MSWD：判断"这个坪自洽吗"

6.1 一个必须回答的问题

假设一个坪里有 30 个点，它们的加权平均是 453.06 ± 1.83 Ma。这个数字的可信度取决于一件没有出现在数字里的事：这 30 个点彼此之间一致吗？

有两种完全不同的情形会给出同样的均值和同样的标准误：

• 情形 A：30 个点都落在 452–454 之间，每个点的误差是 ±2 Ma。散布和误差相称——这是一个真实的单一年龄。
• 情形 B：15 个点在 445 附近，15 个点在 462 附近，每个点的误差是 ±2 Ma。均值同样是 453，标准误同样很小——但这根本不是同一个年龄，是两个年龄混在一起。

情形 B 的均值看起来很精确，但它毫无意义。只看均值和标准误，你无法区分这两个情形。这就是 MSWD 存在的理由。

6.2 公式与直觉

读这个公式的方式：分子是"实际散布"，分母是"如果误差是诚实的，应该有多少散布"。

更严格地看：如果每个点的误差 σᵢ 都是真实可信的、且所有点来自同一个年龄，那么标准化残差 (xᵢ − x̄)/σᵢ 应该服从标准正态分布，它们的平方和服从自由度 n−1 的卡方分布。MSWD 就是这个卡方值除以自由度，期望值恰好是 1。

所以 MSWD 不需要记住公式，只需要记住一句话：它是"实际散布 ÷ 预期散布"，理想值 1。

6.3 三种情形分别意味着什么

6.4 概率读法：比 1 大多少才算大

MSWD 等于 1.2 算不算问题？取决于有几个点。所以除了 MSWD 本身，还应该给出概率：

常规的警示线是 p < 0.05。这个概率的含义要说清楚：它不是"这个年龄有 95% 的概率是对的"，而是"如果这些点真的来自同一个年龄且误差可信，那么出现这么大离散的概率只有 5%"。它是对一致性假设的检验，不是对年龄本身的检验。

pf(mswd, n-1, Inf, lower.tail=FALSE)  ==  chi2_sf(mswd * (n-1), n-1)
                                            ↑ 注意要乘自由度

6.5 一个我们必须承认的局限：重叠窗口

MSWD 的推导假设各个点互相独立。而我们的剖面点不独立——第 3.6 节说过，窗口宽 4 秒、步长 1 秒，相邻窗口共享 3 秒的数据。

这会让 MSWD 系统性偏小，我们可以估算偏差有多大。假设相关系数 ρ ≈ 0.75（共享 3/4 数据），则有效点数约为：

于是 DRUID 内部用的一条判据 MSWD < 1 + 2√(2/(k−1))（用来决定一段是"真实年龄域"还是"过渡带"）会偏宽松：

方向是漏杀而非错杀——该被剥掉的混合窗口没剥干净，但不会把好的域误判成坏的。这个偏差在我们全部报告里都存在，属于已知且已量化的局限，不是错误。诚实地说：ADEPT 报出的 MSWD 也应该按这个方向理解——它偏乐观。

6.6 这两个工具在这件事上的分工

最后值得一提的是：MSWD 在这两个工具里的角色并不相同。

算法是一致的（都是反比方差加权平均 + 卡方上尾概率），所以两边的数字可以直接并列对照。但它们回答的问题不同：DRUID 问"这个深度范围内是不是一个年龄"，ADEPT 问"选出来的这个坪自不自洽"。两个都通过，才说明这个年龄值得写进文章。

6.7 为什么"整段平均"不是一个年龄

到这里有一个自然而然的疑问：既然每个窗口都能算一个年龄，那为什么还要挑坪、分域？ 把所有窗口平均一下不就行了？这个问题值得正面回答——它正是"逐深度"这套做法存在的理由， 也是后面那张"不分域年龄"表存在的理由。

先澄清一个容易混的地方："整段"其实有两种算法，它们不是一回事。 一种是把整段的全部计数合起来、算一个比值再代入年龄方程（下面叫"整段积分"，结果表的 年龄206_238 就是它）；另一种是每个窗口各算一个年龄，再按 1/σ² 加权平均。 两者在单一年龄的测点上很接近，但在真有结构的测点上可以差 20% 以上——因为它们对"哪一段计数更重要" 的取舍不同：积分让计数率高的窗口占更大权重，加权平均让年龄不确定度小的窗口占更大权重。 不要把它们当成同一个数。

再回到那个关键问题。把整段当一个域做加权平均，得到的确实还是一个数， 但这个数有没有意义，取决于一件事：整段的 MSWD 是不是大约 1。 如果剖面本来就有两个年龄（外圈老内圈新，或者反过来），加权平均出来的值落在两者之间， 它不对应任何一个真实时刻——它是一个混合值。而 MSWD 会把这个事实报出来。

那为什么还要专门出一张"不分域"的表？因为对比需要同一把尺子。 分域之后，每个测点得到的域数、域的位置都不一样（示例批次里 33 个点在 2~3 个域之间，另 15 个是均一的）， 测点之间没法直接比；而"不做任何分域"是所有测点都走同一条路径，于是它就是那把尺子—— 整段值、分域主域、整段积分三个数并排，两级之间的差额精确地告诉你"这一步动了多少"。 把这三步的差异摊开，比只报告最后一个数更容易发现自己在哪一步被骗了。

七、结果表里每个数是什么意思

前面六章讲了数字是怎么算出来的。这一章反过来：拿到一个结果文件，从哪里开始看。

7.1 DRUID 的输出：六张表

结果表里的关键列

7.2 把剖面画出来：一个真实测点

7.3 ADEPT 的输出

ADEPT 返回两个表（以及可导出的 Excel）：

• $full——每一个 PELT 切出的段一行，包含全部中间量。调试和复现时看这张：你能看到哪一段在哪一步被过滤掉、以及为什么。
• $summary——只包含最终确认的坪。写文章时看这张。

7.4 三个判据：这批数据能不能用

每次拿到一批结果，先看这三个数字。三项都过，数据可用；任何一项不过，先查原因而不是先改结论。

这三个判据的共同点是：它们都是"内部一致性"检查，不能证明年龄是对的。它们能发现的是"这批数据内部有问题"。真正的正确性要靠外部证据——与已发表年龄的对照、不同方法的交叉验证、地质关系的合理性。软件能保证的只是"没有算错"，不能保证"结论是对的"。这个界限在第八章还会再讲一次。

八、我们怎么知道这些数字是对的

这一章是这份文档存在的主要理由。前面七章解释了每一步"应该"怎么做，但一个流程最难的部分不是设计它，而是确认它真的在按设计运行。特别是当它开始产出看起来合理的数字之后。

8.1 锚点：不许悄悄改数

这套工具最核心的防线是把已经核对过的结果硬编码进测试。具体说：

# tests/testthat/test-adept.R —— 每次运行都会被检查
expect_equal(ages, c(21.61723849, 22.56383385, 23.60464105))

这三个数字来自项目的早期版本，是用真实数据一次一次核对出来的。它们的意义不是"正确"（没有任何单一数字能证明正确），而是"没有变"。任何重构、优化、参数调整，只要动了任何一个数值计算环节，这三个数字就会变，测试立刻失败。

这条防线在这次开发中真的起作用了。我们做了三次较大的重构（拆分 700 行的主函数、把 20 万次循环改成向量化、加入不确定度加权），每次都靠这个锚点确认"默认路径一个数没动"。有一次改动确实让锚点变了——那正是它抓到的一个真实 bug（见 8.4 节第 6 条）。

8.2 两套独立实现互相对照

U-Pb 的计算在 R 和 Python 里各写了一遍（ADEPT 用 R，DRUID 用 Python），这不是浪费——它们可以互相做交叉验证。

有价值的对照发生在统计函数上，因为这类函数最容易"看起来对"却差一点点：

还有一个独立的第三方对照：PELT 分段的实现与 R 的 changepoint 包做过比对（在 n = 60/161/411、1–6 个真实变点的多组序列上，变点集合完全一致），之后才把那个依赖去掉。去掉的理由是：第三方包的内部常数可能随版本变化，而结果不该跟着变。

8.3 端到端：用真实数据跑完整条链

单元测试只能保证"每个零件是对的"，不能保证"装起来是对的"。所以还有一个端到端检查，跑的是仓库里自带的匿名示例批次：

python tests/check_example_batch.py

它会完整跑一遍 85 个测点（约 1 分钟），然后逐项核对：各表的行数、标样 QC 的偏差、样品年龄的分布、深度结构的分布，以及不分域整段结果的判定。这些期望值都硬编码在测试文件里，和 8.1 节的锚点同样的道理。

这个检查在 CI 上也会跑——所以每次代码改动，都会有一台干净的机器用真实数据验证一遍数值没变。

8.4 已知的局限：一份诚实清单

一份只讲优点的文档没有价值。下面是这套流程已知的、我们没打算隐瞒的问题：

第 6 条背后有一个这次开发中发生的真实故事，值得记下来。在给剖面的加权平均做重构之后，锚点年龄变了，而且变化的方式很微妙：剖面里完全平坦的段会得到"未定义"的方差，于是被静默丢弃。原因是一段数学上正确的链条：段内归一化时，零极差的段返回 NaN；NaN 经过一个 pmax 变成 NA；NA 让方差检查返回缺失；于是整段被当作"不合格"扔掉。

而"完全平坦的段"恰恰是最好的坪——它只会在一种输入上出现：做过分馏校正、又关掉了平滑（即 smooth="none"）的数据。这正是为 MSWD 准备的那种输入。换句话说，一个数学上正确的中间步骤，把最重要的一类数据静默地丢掉了，而且丢得没有任何提示。这个 bug 是被锚点测试抓到的；如果没有锚点，它会带着"结果略少几个坪"的症状一直存在下去。

8.5 我们验证不了什么

最后要划清一条线。以上所有验证，保证的都是"计算过程没有错"，不是"这个年龄是对的"。

• 如果输入的 cps 数据本身有问题（比如激光打到了别的矿物、样品表面有污染），流程不会知道。它只会给出一个自洽的、协和度良好的错误年龄。
• 如果某个校正的前提在某个样品上不成立（比如这个锆石的普通铅组成不服从 Stacey-Kramers 模型），流程也不会知道。
• 如果地质解释错了（把捕获锆石当成岩浆锆石），软件更不会知道。

这也是为什么第 7.4 节三个判据的措辞是"内部一致性检查"。软件的责任是把每一步做对、把不确定度诚实地说出来、并且在可能的时候提醒你哪里可疑。判断这个年龄意味着什么，永远是人的工作。

九、如果你要参与

这套工具是在解决一个具体问题时长出来的，所以它的结构反映的是问题的形状，而不是某种"优秀架构"。这意味着它有真实的缺口，而且缺口的位置是可预测的。下面按"需要的技术门槛从低到高"排列。

9.1 最有价值的贡献：报告一个说不通的结果

这件事不需要写代码，但对项目的价值可能超过任何代码贡献。因为开发者在自己的数据上永远看不到别人数据里的毛病。

一个有用的报告长这样：

1. 哪一步出了问题。是 DRUID 的还原阶段（"标样 QC 看起来不对"），还是 ADEPT 的分段阶段（"坪切得位置不对"）？
2. 期望是什么、实际是什么。"这个测点的剖面在 15–25 秒之间明显是平的，但 ADEPT 没有在这里切出坪" 比 "结果不对" 有用一百倍。
3. 具体数值。哪个测点、哪一秒、期望多少、实际多少。如果是数值差异，请说明差多少量级（1% 和 10 倍是完全不同的问题）。
4. 能不能最小复现。如果只用一个测点就能重现，问题就解决了一半。

9.2 加数据：不同仪器、不同批次、已知年龄的样品

这套流程的很多参数是从本实验室的单台仪器（iCAP RQ）和有限批次上标定出来的：搜索窗口宽度、剥蚀定位的阈值倍数、方差阈值、外部重现性的上下限。它们在别的仪器上大概率需要调整。

所以下面这些数据特别有价值：

• 别的仪器/实验室的数据——最能暴露"参数是仪器相关的"这类问题。
• 已知年龄的样品——比如已经用 ID-TIMS 测过的标样。这是唯一能验证绝对准确度的方式（8.4 节第 5 条）。
• 刻意构造的异常案例——比如普通铅很高的样品、或者年龄跨度很大的批次。边界情况是最容易出错的地方。

加数据的时候有一件事必须做：如果数据要进仓库，先匿名化。做法在 AGENTS.md 里写着——样品名替换成顺序编号，但标样名必须保留（它们是国际标准物质，而且角色判定靠名字匹配，改了名流程就找不到主标了）。

9.3 加测试：一个我们自己踩过的坑

写测试的时候，最容易犯的错误是凭印象写出一个"应该是"的值。这件事发生在我们身上，值得详细讲。

在给卡方上尾概率写测试时，我"记得"χ² = 200、自由度 10 的概率大约是 3.0559×10⁻³⁸，就把它写进了断言。测试立刻失败，实测值是 1.6139×10⁻³⁷——差了五倍。

这时候有两种反应：去改实现（因为"测试说不对"），或者去核对参考值。我们去 R 里算了一遍：R 给出的正是 1.6139305336977309×10⁻³⁷。也就是说实现是对的，是我编的期望值错了。

同样的原则也适用于加新的端到端检查：期望值要么来自同一批次的历史输出，要么来自手工验证过的计算。测试的价值完全取决于它的期望值有多可信。

9.4 加代码：先看这两份文件

如果要在代码层面动手，有两个文件必须先读，因为它们记录的是读代码看不出来的东西：

改完之后必须做的三件事：静态检查通过、自检通过、以及跑一次端到端对照确认数值没变。第三条不能省——它才是真正在保护这套东西的那一条。

9.5 关于"开源能带来什么"

这份文档最后想说一件不完全是技术的事。

这套工具的价值不取决于它有多聪明，而取决于有多少人真的用它、并且把遇到的问题说出来。原因很直接：软件的缺口是由使用方式决定的，而使用方式是开发者想象不出来的。

举两个真实的例子，都来自这次开发过程：

• 有人把结果输出重定向到文件，于是发现中文进度信息在 Windows 上会让程序崩溃——这个问题在交互式控制台下永远不会出现。
• 有人问"为什么样品年龄中位数和我手算的不一样"，于是发现那个统计量把 43 个标样测点也算进去了——这个 bug 活了很久，因为打印出来的 n=48 让人以为它只统计了样品。

这两个都不是"代码写得不好"，而是"没人想到会这样用"。用户反馈在这里不是帮开发者找 bug，而是帮开发者看到自己没看到的用例。

所以如果你在用这套东西，遇到任何说不通的地方——哪怕只是"我觉得这个数字怪怪的"——那都值得说出来。这份文档里所有的坦承（包括 8.4 节那份局限清单）都是为了让这种交流更容易发生：承认自己不知道，才有人愿意告诉你答案。