Isoclock 保姆级解读：从零看懂一套 LA-ICP-MS U-Pb 定年流程

面向第一次接触 U-Pb 定年、或者第一次打开这份源码的人 —— 每一节都先讲人话，再指给你看代码在哪一行

图解版　·　56 张图　·　行号对应重构后的 3676 行 Isoclock2.0.py　·　图里所有数字由脚本实算，未手写

---

写在最前面

这份文档不是官方手册的重写，也不是 API 文档。它的来历是：有人把 Isoclock2.0.py 从头到尾读了一遍，把里面每一处"为什么"整理出来，写成一条从不懂到懂的路径。

所以先说清楚三件事：

1. 它不能代替原始手册。软件的操作说明在仓库里的两份 PDF 里，本文只做补充和解释。
2. 行号对应的是本仓库当前版本（重构后的 3676 行 Isoclock2.0.py）。上游原版行号不同，如果你用的是原版，按函数名搜索即可。
3. 本文对科学逻辑保持原样。代码里有些写法是故意保留的（哪怕看起来别扭），因为它们会影响已经发表的结果。凡是这种地方，我都会明确标注"原样保留"，并说明为什么不能顺手改。

阅读建议：

• 完全的新手：从第 1 章顺着读到第 5 章。
• 只想知道"普通铅到底是什么"：直接跳第 3 章。
• 只想知道"这套流程每一步在干嘛"：跳第 4 章。
• 想知道"我的结果怎么让审稿人接受"：跳第 5 章。
• 想知道"这软件哪儿不行"：跳第 6 章。

---

---

第 1 章　先弄懂在测什么：U-Pb 定年的三条腿

这一章只讲最基础的三件事：我们到底在测什么、一个比值凭什么能变成年龄、以及为什么一次测量能给四个年龄。这三件事是后面所有内容的底座。

1.1 三个计时器

自然界的铀和钍会慢慢变成铅，这个"变"的速度（衰变常数）是物理常数，不随温度、压力改变。这就是铀铅定年能用的根本原因。

母体	子体	衰变常数 λ（/年）	半衰期
²³⁸U	²⁰⁶Pb	1.55125 × 10⁻¹⁰	4.468 Ga
²³⁵U	²⁰⁷Pb	9.8485 × 10⁻¹⁰	0.704 Ga
²³²Th	²⁰⁸Pb	4.9475 × 10⁻¹¹	14.01 Ga

这三个数字在代码里到处都是，但写法有点"分裂"，第一次看很容易被绕晕，提前说清楚：

• 第 470 行写的是 0.000000000155125，单位是 每一年（1.55125 × 10⁻¹⁰）。
• 第 116-118 行（SK2model）写的是 0.155125，单位是 每百万年（1.55125 × 10⁻¹⁰ × 10⁶）。
• 两者是同一个常数，只差 10⁶。看到 0.155125 就知道年龄以 Ma 为单位，看到 0.000000000155125 就知道年龄以年为单位 —— 仅此而已，不是笔误。

第四个数字是 ²³⁸U/²³⁵U = 137.818（现代值）。它出现在几乎所有比值换算里：因为 ²³⁵U 和 ²³⁸U 是同一个元素的同位素，你测不出它们的绝对量，只能测它们的比值，而这个比值假设是常数。

1.2 年龄公式：为什么长成 ln(r+1)/λ 这个样子

假设矿物结晶那一刻，里面只有母体 U，一点子体 Pb 都没有（后面会看到，这个假设往往是错的，这就是整篇文章的主题）。过了时间 t 之后：

母体剩余：U = U₀ · e^(−λt)
子体生成：Pb* = U₀ · (1 − e^(−λt))

两式相除，U₀ 消掉，得到比值：

Pb*/U = e^(λt) − 1

反解出 t：

t = ln(Pb*/U + 1) / λ

这就是代码里那个反复出现的模式（第 470 行为例）：

log(abs(date_all[2][i]*f206_238+1))/(0.000000000155125)/1000000
#   └─ 校正后 206Pb/238U ──┘   └── λ₂₃₈（每年）──┘  └ 换成 Ma ┘

理解了三件事，这一行就通了：

• +1 是因为 e^(λt) − 1 反解时要把 −1 移过去；
• 除法用的是 λ 而不是别的，因为指数反函数是对数；
• 最后 /1000000 是把"年"换成"百万年"。

注意 abs()。这个绝对值函数体现了一种"不抛错"的工程态度：如果某个点的比值算出来是负数（信号太差、背景扣过头），数学上取对数会报错，程序直接崩。加了绝对值，程序能跑完，但结果是没有物理意义的垃圾值。这是这套代码的一个通病：宁可给你一个假数字，也不肯中断。后面第 6 章会回到这一点。

同一条公式换了母体就得到另外三个年龄：

206Pb/238U 年龄: log(r+1) / λ₂₃₈        # 第 470 行
207Pb/235U 年龄: log(r+1) / λ₂₃₅        # 第 472 行
208Pb/232Th 年龄: log(r+1) / λ₂₃₂       # 第 474 行

1.3 为什么单次测量能出四个年龄

一次激光打点，你同时得到 ²⁰⁶Pb、²⁰⁷Pb、²⁰⁸Pb、²³²Th、²³⁸U 五个信号，于是能算三对"母-子"关系和一堆比值。理论上，如果这个点从结晶到现在一直是封闭体系：

• ²⁰⁶Pb*/²³⁸U 给出的年龄
• ²⁰⁷Pb*/²³⁵U 给出的年龄
• ²⁰⁷Pb*/²⁰⁶Pb 给出的年龄

三者必须一致。 把它们画在一张图上（横轴 ²⁰⁷Pb/²³⁵U，纵轴 ²⁰⁶Pb/²³⁸U），一致的点全部落在一条曲线上 —— 这条曲线叫协和线（concordia）。

不一致的两种情况，含义完全不同：

偏离方向	常见原因	后果
沿一条指向原点的弦排开（落在协和线下方）	放射成因铅丢了（后期变质、流体、退晶质化）	叫"不协和"，需要做铅丢失模型
沿一条近水平的混合线向右漂（落在协和线右下方）	混进了普通铅	叫"普通铅污染"，需要做普通铅校正

这两件事长得像，但处理方式完全相反。第二章会看到，Isoclock 存在的全部理由，就是第二种情况。

1.4 ²⁰⁷Pb/²⁰⁶Pb 年龄为什么特殊

前面三个年龄里，后一个特别值得单独说，因为它在代码里是一个没有解析解的迭代函数。

Age76Pb()（第 122 行）要解决的问题是：已知测到的 ²⁰⁷Pb/²⁰⁶Pb 比值，反推年龄。

而理论关系是：

²⁰⁷Pb*/²⁰⁶Pb = (1/137.818) · (e^(λ₂₃₅·t) − 1) / (e^(λ₂₃₈·t) − 1)

（对应第 130 行的 Rap。）

这个式子里 t 同时出现在分子分母的指数上，没办法用代数解出来，只能数值求解。代码用的是"割线法"（secant method）：

while abs(Rap76 - Rap) > 0.00005 and N < 10:   # 第 133 行
    delta = Rap76 - Rap
    if Rap < Rap76:
        Tmin = Tm          # 猜小了，把下界提上来
    else:
        Tmax = Tm          # 猜大了，把上界压下去
    Rapi = (t=Tmin 时算出的比值)     # 第 138 行
    Raps = (t=Tmax 时算出的比值)     # 第 139 行
    Tm = Tmin + (Tmax - Tmin) * (Rap76 - Rapi) / (Raps - Rapi)
    Age76Pb = Tm
    N = N + 1

逻辑本身没问题。但这里有一个真实的缺陷，我必须在第 6 章单独讲：循环条件是 abs(Rap76 - Rap) > 0.00005，而 Rap 是在进入循环之前算好的、循环里从来没有更新过。所以这个条件一旦第一次为真，就永远为真 —— 循环实际固定跑 10 次就退出。

这不是致命的（10 次割线法对多数年龄已经够收敛），但它让"收敛判据"变成了一句空话。1.5–2.2 Ga 的老年龄会因此偏高最多约 35 Ma（实测在 2.1 Ga 处是 +34.77 Ma）。这一条已经写进仓库 README 的 Known issues，第 6 章有完整的实测数据。

---

第 2 章　分水岭：为什么 Isoclock 不是给锆石写的

这一章回答你问的最关键的那个问题。

2.1 先看锆石过的是什么日子

锆石（ZrSiO₄）在 U-Pb 定年里是"天选之子"，原因有两个：

第一，它排斥铅。 锆石的晶格是给 Zr⁴⁺ 和 Si⁴⁺ 准备的，Pb²⁺ 半径太大、电荷不对，几乎进不去。所以锆石结晶时抓到的"初始铅"极少。

第二，它富含铀。 锆石里的 U 常是几十到几千 ppm 的量级。

两个因素叠加，结果就是：锆石的 ²⁰⁶Pb/²⁰⁴Pb 常常大于 1000，也就是说，你测到的铅里，99.9% 以上都是铀衰变自己长出来的，外面混进来的那点普通铅可以忽略。

所以锆石的常规处理流程是：

1. 不做普通铅校正（或者只在极少数情况下做）；
2. 主要工作是处理铅丢失导致的"不协和"；
3. ²⁰⁷Pb/²⁰⁶Pb 比值在锆石里主要当成一个"体检指标"，用来看协和度。

如果你拿锆石数据跑 Isoclock，正确的选择是：计算模式选 "Standard not Corr."，样品普通铅校正选 "None"。界面里那个 Zircon(91500) 按钮指的是用 91500 当标样，不是"这块样品是锆石所以要走特殊流程"，这两件事别混。

2.2 换一批矿物，情况就反过来了

Isoclock 面向的是磷灰石、榍石、黑钨矿、锡石这一类矿物。它们的共同特征，恰恰和锆石相反：

	锆石	磷灰石	榍石	黑钨矿 / 锡石
典型 U 含量	几十–几千 ppm	< 50 ppm	几十–几百 ppm	低，变化大
晶格对 Pb 的排斥	极强	弱	中等	弱
²⁰⁶Pb/²⁰⁴Pb	常 > 1000	常 < 200	几十–几百	常常很低
普通铅占比	< 1%	可达 10–90%	常 10–50%	常很高
主要矛盾	铅丢失	普通铅	普通铅	普通铅 + 基体

2.3 用"误差传递"把这个困难量化

这是理解整件事最有用的一把尺子。

普通铅的 ²⁰⁶Pb/²⁰⁴Pb 大约是 18.7（第 3 章会给出来历）。假设你的样品实测 ²⁰⁶Pb/²⁰⁴Pb = R，那么普通铅在总 ²⁰⁶Pb 里的占比就是：

f_c = 18.7 / R

而你要做的那个减法，是拿"普通铅的量"去减。这份量本身的相对误差会原封不动地传到放射成因铅上，并且被放大：

放射成因 ²⁰⁶Pb 的相对误差 = σ_204 × f_c / (1 − f_c)
                                    └─ 204 通道的相对误差

样品 ²⁰⁶Pb/²⁰⁴Pb	普通铅占比 f_c	误差传递因子 f_c/(1−f_c)	若 204 通道有 5% 相对误差
1000（典型锆石）	1.9%	0.02	→ 年龄误差 0.10%
200	9.3%	0.10	→ 0.52%
100	18.7%	0.23	→ 1.15%
50	37.4%	0.60	→ 2.99%

两行对照着看，结论就出来了：

• 锆石（R ≈ 1000）：就算 ²⁰⁴Pb 通道有 5% 的误差，传到年龄上也只剩 0.1%。所以普通铅在锆石里是"扣也行、不扣也行"。
• 黑钨矿（R ≈ 50）：同样 5% 的 ²⁰⁴Pb 误差，传到年龄上是 3%。如果这块矿物是 200 Ma，那就是额外的 ±6 Ma —— 而这还只是 ²⁰⁴Pb 一个通道带来的。

这就是为什么低 U/Pb 矿物的数据，误差条天生就比锆石宽。不是测得更差，而是年龄信息被普通铅稀释了。

2.4 一句话总结这个差异

这句话解释了 Isoclock 为什么长成现在这样：

• 为什么它要提供四种普通铅校正路径（第 3 章）—— 因为换一种矿物，能用的通道就不一样；
• 为什么它有一个 "Standard Corr. / Standard not Corr." 开关 —— 因为是否需要对标样本身也做普通铅校正，取决于标样含不含普通铅（MAD-NEW 磷灰石含，91500 锆石不含）；
• 为什么输出表里专门留出 Tera-Wasserburg 图的两列数据（见第 5 章）—— 因为对这类矿物，最能让审稿人信服的证据不是"我扣了普通铅"，而是"我用等时线图证明我扣对了"。

2.5 还有一件比普通铅更根本的事

这一点必须说，因为它比技术细节重要：

普通铅校正是"数据整理"，不是"地质论证"。

如果你测的矿物在结晶之后经历过后期流体改造、重结晶、或者多期生长的混合，那么它的 U-Pb 体系可能已经被重置或者部分重置。这时候，无论你选哪条校正路径、参数调得多精细，你算出来的"年龄"都可能没有对应的地质事件 —— 它只是一个数学上的平均值。

所以拿到一个漂亮年龄时，先问自己三个问题：

1. 这个点和同一块样品其他点一致吗？
2. 这个矿物有独立的证据（结构观察、共生关系、其他定年方法）证明它是一期形成的吗？
3. 如果我把校正方法换一种，年龄会变多少？

第 5 章会给出具体的操作方法。

---

第 3 章　普通铅：把这件事彻底说透

这一章是全文最核心的部分。

3.1 普通铅到底是什么

矿物结晶的那一刻，它的晶格会从周围环境里"抓"一些铅进去。这些铅的来源五花八门：岩浆本身的、围岩带进来的、流体搬运来的。

关键在于：这些铅的量和这块矿物的年龄没有任何关系。 它们是"混进来"的，不是"长出来"的。

但是质谱仪不知道这件事。它只告诉你："这个样品里有 ²⁰⁶Pb、²⁰⁷Pb、²⁰⁸Pb、²⁰⁴Pb。"

所以你测到的量永远是两部分之和：

实测 Pb = 放射成因 Pb*（来自本矿物 U/Th 衰变，携带年龄信息）
        + 普通 Pb_c  （结晶时混进来的，不携带年龄信息）

定年的全部工作就是：**想办法估出 Pb_c，从实测值里减掉，只留下 Pb*。**

3.2 ²⁰⁴Pb：普通铅的指纹

铅有四个稳定同位素。其中三个会从 U 和 Th 变出来：

• ²⁰⁶Pb ← ²³⁸U
• ²⁰⁷Pb ← ²³⁵U
• ²⁰⁸Pb ← ²³²Th

只有 ²⁰⁴Pb 不是任何一条主要衰变链的产物。（²³²Th 链有一个极稀有的分支会产生极微量 ²⁰⁴Pb，可以忽略。）

这是一个极其有用的性质。它意味着：

有了这个指纹，再看代码里那些比值就很好懂了：

符号	含义	代码位置
Q8	普通铅的 ²⁰⁶Pb/²⁰⁴Pb	第 3177 行
R8	普通铅的 ²⁰⁷Pb/²⁰⁴Pb	第 3178 行
S8	普通铅的 ²⁰⁸Pb/²⁰⁴Pb	第 3179 行
Pbc	普通铅的 ²⁰⁷Pb/²⁰⁶Pb（= R8/Q8）	第 3180 行
P382	放射成因的 ²⁰⁷Pb/²⁰⁶Pb	第 3175 行

这几个变量名看着不友好，但记住 8 代表 208Pb、4 代表 204Pb、6 代表 206Pb、7 代表 207Pb 就能读懂了。Q8 就是 "²⁰⁶/²⁰⁴"，R8 是 "²⁰⁷/²⁰⁴"，以此类推。

3.3 可是 ²⁰⁴Pb 特别难测

如果 ²⁰⁴Pb 这么有用，为什么还需要别的校正方法？因为它有一个致命弱点：

丰度太低。 自然界的铅里，²⁰⁴Pb 只占 1.4%，而 ²⁰⁶Pb 占 24%、²⁰⁸Pb 占 52%。你要在低 U 矿物本来就微弱的信号里，去精确测量那 1.4%。

再加上一个专门的干扰：

²⁰⁴Hg 的干扰。 汞和铅的原子质量都是 204（²⁰⁴Hg vs ²⁰⁴Pb），单靠质量分辨率分不开。而汞在 Ar 等离子体和样品气溶胶里并不少见。所以测"²⁰⁴Pb"的通道，实际上同时接收了 ²⁰⁴Pb 和 ²⁰⁴Hg。

代码试图处理这个问题：loaddata() 在读文件时额外读了一个 202Hg 通道（第 280 行），就是留着用来反推 ²⁰⁴Hg 的量。

但是这段计算实际上没有生效。 在 _net_signal() 里（第 1569-1571 行）有一段注释：

# 204 通道：原实现直接用原始 204 减背景，未做 Hg 干扰扣除
# （源码里那段 Hg204cal_s 恒为 0、m204_s 从未被引用）
'pb204': np.average(sig['Hg204_s']) - np.average(sig['Hg204_b']),

也就是说，软件直接把 204 通道的原始计数当成普通铅的量，没有扣 Hg。如果仪器导出数据时已经做过 Hg 校正，这样是对的；如果没做，那 ²⁰⁴Pb 就是被高估的，所有基于 204 的校正都会偏。这条也在 README 的 Known issues 里。

你需要做的：确认你的质谱仪导出数据时是否已经做了 Hg 校正（Thermo 软件的 "Hg correction" 之类）。这是你自己的数据你自己的责任。

3.4 四条校正路径：各用什么当"比例尺"

既然 ²⁰⁴Pb 又弱又脏，就得有替代方案。于是有了四条路。

先看一个统一的思路。假设我们想知道"混合进来的普通铅贡献了多少 ²⁰⁶Pb"。如果我们知道普通铅的三个关键比值（Q8、R8、Pbc 等），那么任何一条能反映"普通铅占比"的观测都能反算出这个比例。

四条路的区别，就是"用什么当比例尺"：

路线一：²⁰⁷Pb 法（界面：207Pb）

比例尺：实测 ²⁰⁷Pb/²⁰⁶Pb 比放射成因值高出了多少。

原理：你要量的是 ²⁰⁶Pb 里有多少是普通铅带来的。虽然你测不出 ²⁰⁴Pb，但你知道普通铅的 ²⁰⁷/²⁰⁶ 是多少（Pbc）。而实测的 ²⁰⁷/²⁰⁶ 是普通铅和放射成因铅的加权混合，混合比例直接可以从"实测值偏离放射成因值多少"里解出来。

代码里的实现（_ratios_std207，第 1600 行）：

f206 = (n7 / n6 - P382) / (Pbc - P382)          # 第 1604 行：解出普通铅中 206 的占比
r['r76'] = (n7 - n6 * Pbc * f206) / (n6 * (1.0 - f206))
r['r68'] = (n6 / nU) * (1.0 - f206)

f206 就是"²⁰⁶Pb 里普通铅占的比例"（在代码里叫 f206，其实就是文献里的 f₂₀₆）。看出来了吗 —— 当 f206 = 0（没有普通铅）时，r68 就退化成 n6/nU，也就是最朴素的那个比值。整个校正就是"把混合物解回纯放射成因端元"。

优点：不需要测 ²⁰⁴Pb，信号最强的两个通道（206、207）就够了。 代价：你必须知道普通铅的 ²⁰⁷/²⁰⁶ 是多少（下一节讲这个数字从哪来），而且这个假设错了会直接错到年龄上。

路线二：²⁰⁸Pb 法（界面：208Pb）

比例尺：²⁰⁸Pb 与 Th 的关系。

这一路来自 Zack et al. (2011)，原本是为金红石设计的。它不用 ²⁰⁴Pb，而是利用 ²⁰⁸Pb 和普通铅的 ²⁰⁶/²⁰⁸、²⁰⁷/²⁰⁸ 比值。

代码在 _ratios_std208（第 1647 行），核心是这一行：

r['r76'] = (pb207_t - pb208_t * Pbc_78) / (pb206_t - pb208_t * Pbc_68)
r['r68'] = (pb206_t - pb208_t * Pbc_68) / net['nU']

其中 Pbc_78、Pbc_68 就是由 SK2model() 算出来的模式铅 ²⁰⁷/²⁰⁸ 和 ²⁰⁶/²⁰⁸（第 119-120 行）。逻辑是"用 ²⁰⁸Pb 当参照物，把属于普通铅的那部分 206、207 减掉"。

优缺点：适合 Th/U 比较稳定的矿物；如果矿物本身 Th 含量很低或者波动很大，²⁰⁸Pb 反而更不可靠。

路线三：²⁰⁴Pb 法（界面：204Pb）

比例尺：直接测的 ²⁰⁴Pb。

最直接、最少假设的一条路：既然 ²⁰⁴Pb 全部是普通铅，那测到多少 ²⁰⁴Pb 就知道有多少普通铅，再乘以 Q8、R8 就知道该扣掉多少 206、207。

代码在 _ratios_std204（第 1629 行）：

m74 = (avg(Pb207_s) - avg(Pb207_b)) / net['pb204']    # 第 1639 行：实测 207/204
m64 = (avg(Pb206_s) - avg(Pb206_b)) / net['pb204']    # 第 1640 行：实测 206/204
r['r76'] = ((Pb207_s - avg(Pb207_b) * ((m74 - R8) / m74)) / n6 * ((m64 - Q8) / m64))
r['r68'] = n6 * ((m64 - Q8) / m64) / nU

优点：物理上最直接，假设最少。 代价：²⁰⁴Pb 信号弱（见 3.3），误差大；而且如果 ²⁰⁴Pb 的净信号为零，这里会直接除零崩溃 —— 这是 README 里的 Known issue 第 2 条。代码没有做任何防护。

路线四：标定 ²⁰⁴Pb 法（界面：Cal 204Pb）

比例尺：用标样实测的 204 通道反过来标定。

这一路来自 Chew et al. (2014)，是专为"有普通铅的副矿物标样"设计的。它的巧妙之处在于：不追求 ²⁰⁴Pb 的绝对准确，只用它的相对变化。

在 _ratios_std_cal204（第 1613 行）里，先算一个 n383：

n383 = (r76_std * Pb206_s - avg(Pb206_b) * r76_std - Pb207_s + avg(Pb207_b)) \
       / (Q8 * r76_std - R8)

这一大串的意思是：用标样"应该有的"放射成因 ²⁰⁷/²⁰⁶（P382）去反推，标样每个时间点在 204 通道上"多出来"的那部分信号（含普通铅和 Hg）。然后把这个 n383 按 Q8、R8、S8 分配到 206、207、208 上扣掉。

这一路的好处：它让 204 通道的系统性偏差（比如 Hg 干扰、绝对灵敏度标定不准）在"标样 vs 样品"的比较中被抵消掉一大半。这也是为什么 README 里建议"标样含普通铅时选 Standard Corr. + Cal 204Pb"。

四条路的对照表

界面选项	代码方法名	比例尺	关键假设	什么时候用
207Pb	_METHOD_207	²⁰⁷/²⁰⁶ 的偏离	知道普通铅 ²⁰⁷/²⁰⁶	204 通道信号差，但 206/207 信号好
208Pb	_METHOD_208	²⁰⁸Pb 与 Th	知道普通铅 ²⁰⁶/²⁰⁸、²⁰⁷/²⁰⁸	Th 含量稳定的矿物
204Pb	_METHOD_204	直接测 ²⁰⁴Pb	204 通道干净、H 已校正	204 信号足够好
Cal 204Pb	_METHOD_CAL204	标样标定的 204	标样与样品基体相近	标样本身含普通铅时的首选

注意界面上的 Standard Corr. / Standard not Corr.（第 3375-3376 行）是管标样的，上面这四条路既管标样也管样品 —— 样品那一组由右侧的 None / 207Pb / 204Pb / 208Pb / Cal 204Pb 按钮（第 3409-3413 行）决定。

3.5 那个"普通铅比值"从哪来：Stacey–Kramers 两阶段模型

现在回到最关键的问题：普通铅的 ²⁰⁶/²⁰⁴、²⁰⁷/²⁰⁴、²⁰⁸/²⁰⁴ 是你不知道的，你怎么可能知道？

答案是一个模型。

地球的铅同位素比值不是恒定的 —— 它在几十亿年里一直在被 U 和 Th 衰变"喂养"而升高。Stacey & Kramers (1975) 提出了一个广被接受的两阶段模型来描述这件事：地球的 U/Pb 比值在 3.7 Ga 前后发生过一次跃变，之后按新的比值演化到今天。

代码里有这个模型的全部常数。看 Cal_age()（第 3170 行）：

def Cal_age(age):
    a = exp(0.000000000155125 * age * 1000000) - 1      # 206Pb*/238U
    b = exp(0.00000000098485  * age * 1000000) - 1      # 207Pb*/235U
    c = exp(0.000000000049475 * age * 1000000) - 1      # 208Pb*/232Th
    P382 = 1/137.818 * (b/a)                             # 放射成因 207Pb/206Pb
    Q8 = 18.700 - a * 9.735                              # 普通铅 206Pb/204Pb
    R8 = 15.628 - b * (9.735/137.818)                    # 普通铅 207Pb/204Pb
    S8 = 38.630 - c * 36.837                             # 普通铅 208Pb/204Pb
    Pbc = R8 / Q8                                        # 普通铅 207Pb/206Pb

逐个常数解释：

常数	含义	来历
18.700	今天地球铅的 ²⁰⁶Pb/²⁰⁴Pb	Stacey & Kramers (1975) 现代值
15.628	今天地球铅的 ²⁰⁷Pb/²⁰⁴Pb	同上
38.630	今天地球铅的 ²⁰⁸Pb/²⁰⁴Pb	同上
9.735	μ，即 ²³⁸U/²⁰⁴Pb（第二阶段）	同上，约 9.74
36.837	μ × κ，即 ²³²Th/²⁰⁴Pb（κ ≈ 3.78 为 ²³²Th/²³⁸U）	SK2model 里写作 36.84
11.152 / 12.998 / 31.23	3.7 Ga 时地球铅的三个比值	见 SK2model（第 116-118 行）
3.7	两阶段模型的分界时间（Ga）	同上

读法：Q8 = 18.700 − a × 9.735 的意思是

R8、S8 同理。然后 Pbc = R8/Q8 就给出了"T 时刻模式铅的 ²⁰⁷/²⁰⁶"。

SK2model() 是同一个模型的另一种写法，从 3.7 Ga 的初始值往后加：

Pbc_64 = 11.152 + 9.74 * (exp(0.155125*3.7) - exp(0.155125*age/1000))     # 第 116 行
Pbc_74 = 12.998 + (9.74/137.818) * (exp(0.98485*3.7) - exp(0.98485*age/1000))
Pbc_84 = 31.23 + 36.84 * (exp(0.049475*3.7) - exp(0.049475*age/1000))

两者数学上等价。实测对照（同一 T 下两个函数各自算出的 ²⁰⁶Pb/²⁰⁴Pb）：

T	Cal_age 的 Q8	SK2model 的 Pbc_64	差异
1 Ma	18.6985	18.7018	0.003
100 Ma	18.5478	18.5510	0.003
1000 Ma	17.0664	17.0689	0.002
3700 Ma	11.1526	11.1520	0.001

两条曲线全程重合，残余不超过 0.003（约 0.02%，符号在 3 Ga 附近翻转一次），来自 μ 与 κ 在两处取了不同的舍入位数。它们说的是同一件事。

（顺带一提：Cal_age 用 36.837、SK2model 用 36.84，Cal_age 用 9.735、SK2model 用 9.74 —— 这是从原始文献不同舍入位数抄来的，差异小于 0.01%，可以忽略，但它确实说明这些常数是手工录入的，将来换模型或换引用时必须逐处核对。）

3.6 为什么用"标样年龄"来算这些参数

这是初学者最容易卡住的一点。

看 skmethod()（第 3188 行，界面选 "S & K" 时触发）：

Standard_age = tk.simpledialog.askfloat(..., initialvalue='485')   # 问你标样年龄
Standard_name = tk.simpledialog.askstring(..., initialvalue='MAD')
R8, Q8, S8, P382, a, b, c, Pbc = Cal_age(Standard_age)              # 用标样年龄算参数

为什么是标样年龄，不是样品年龄？

因为标样里的普通铅，是标样结晶的那个时刻混进去的。MAD-NEW 磷灰石标样结晶于约 483 Ma，它当时抓到的地球铅就是 483 Ma 前的地球铅 —— 其 ²⁰⁶/²⁰⁴ 是 18.700 − (e^(λ·0.483Ga) − 1) × 9.735 ≈ 17.94，而不是今天的 18.700。

所以：标样的普通铅参数要按标样年龄算。

而样品的普通铅参数是另一套，由 PbCorrSamples()（第 2883 行）在界面上单独问你 —— 要么直接输入同位素比值，要么输入"样品的估计年龄"让程序用 S&K 模型算。

这就回答了"为什么那么设置"：普通铅的 ²⁰⁷/²⁰⁶ 不是常数，它随时间变化。你必须告诉软件"这块矿物的普通铅大概是哪个年代的"，或者"它长什么样"。

3.7 那些默认值是怎么来的（这一段能解答很多疑惑）

PbCorrSamples() 里每个分支都给了默认值。逐个验证一下，你会发现它们不是随便写的：

207Pb 法（第 2900、2905-2906 行）

initialvalue = '100'          # 样品估计年龄，默认 100 Ma
initialvalue = '0.048015'     # 放射成因 207Pb/206Pb
initialvalue = '0.842185'     # 初始（普通）207Pb/206Pb

把 age = 100 Ma 代进第 2902-2903 行：

• Radioactive_Pb = (1/137.818)(e^(λ₂₃₅·100Ma)−1)/(e^(λ₂₃₈·100Ma)−1)。算出来 = 0.048036
• common_Pb = (15.628 − (e^(λ₂₃₅·100Ma)−1)·μ/137.818) / (18.700 − (e^(λ₂₃₈·100Ma)−1)·μ)。算出来 = 0.842185

所以那两个"神仙数字"就是 100 Ma 时的 S&K 模式铅值。 手动输入 0.048015 和 0.842185，等价于告诉程序"我的样品普通铅是 100 Ma 的模式铅"。

（0.842185 与精确值逐位吻合；0.048015 与精确值 0.048036 差 0.04%，这是原作者手工取整留下的痕迹。这类"看起来随意、其实有出处"的默认值，照抄之前一定要先反算一遍它对应什么假设。）

204Pb 法（第 2926-2930 行）

initialvalue = '15.6207'   # 初始 207Pb/204Pb
initialvalue = '18.5478'   # 初始 206Pb/204Pb

同样代 age = 100 Ma：

• 18.700 − (e^(0.0155125) − 1) × 9.735 = 18.700 − 0.1522 = 18.548 ✓
• 15.628 − (e^(0.098485) − 1) × 9.735/137.818 = 15.628 − 0.00731 = 15.621 ✓

完全对上。

208Pb 法（第 2914-2919 行）

initialvalue = '0.48240'   # 初始 206Pb/208Pb
initialvalue = '0.40628'   # 初始 207Pb/208Pb

• 18.548 / 38.447 = 0.48242 ✓
• 15.621 / 38.447 = 0.40629 ✓

（其中 38.447 = 38.630 − (e^(0.0049475)−1) × 36.837，也是 100 Ma 的值。）

三个分支的默认值全部指向同一个年龄：100 Ma。 理解了这一点，你就知道这些数字该怎么改了 —— 它们应该对应你样品的真实估计年龄，而不是照抄默认值。

3.8 必须承认的事：模型铅只是"平均猜测"

上面整套东西有一个隐含前提：你的样品的普通铅，长得像那个年代的地球平均铅。

这个前提经常不成立：

• 如果样品来自地壳重熔的岩浆，普通铅可能明显富 ²⁰⁷Pb 和 ²⁰⁸Pb（地壳的 U/Pb、Th/Pb 更高）；
• 如果样品周围有方铅矿等硫化物的铅混入，普通铅可能极端富放射成因；
• 如果样品经历了多期流体，普通铅可能本身就是一个混合物。

模型铅给的是"统计意义上最可能的值"，不是"你这块样品的真实值"。

所以有一条原则要记住：

具体来说：

1. 如果你的 ²⁰⁴Pb 信号足够好，优先用 204 或 Cal 204Pb 路线，因为它是"实测"的；
2. 如果研究区有已知的普通铅同位素组成（比如同一矿区的方铅矿数据），直接输入那些比值，比用 S&K 模型强得多；
3. 最重要的一条：用 Tera-Wasserburg 图（第 5 章）让数据自己说出普通铅组成，而不是由你规定。这是唯一能让审稿人真正放心的做法。

---

第 4 章　从原始文件到年龄：每一步在干什么

这一章把整条流水线拆成 6 步，跟着数据走一遍：从仪器导出的原始文件，到最终输出表里的每一列。

4.1 全局视野

界面上的九个步骤（README 里的 (1)–(9)）对应的就是这条链。下面逐步展开。

4.2 第 0 步：仪器必须导出哪 8 个通道

loaddata()（第 273 行）返回固定的 8 个通道：

序号	通道	用途
0	Time	时间轴，用于把秒换算成数组下标
1	202Hg	监控 Hg 干扰（注意：这段计算目前没生效，见 3.3）
2	204Pb	普通铅指纹（实际是"204 通道"，含 Hg）
3	206Pb	主要年龄信息
4	207Pb	主要年龄信息
5	208Pb	Th-Pb 体系
6	232Th	208Pb/232Th
7	238U	分母

少任何一个都会报错。这是 README 里"读这个再开始"四条注意事项的第 4 条。

为什么必须齐全：因为不同校正路径要用不同的通道组合。哪怕你打算用 207Pb 法（不用 204），程序在设计上也要求 204 通道存在 —— 这是软件的强制项，不是科学上的必需。

4.3 第 1 步：读文件（loaddata）

三种仪器的导出格式不一样，代码用三个分支处理（第 275 行起的 mode = ele.get()）：

界面选项	mode	处理方式
Element	0	跳过 13 行，列名硬编码为 Time, 202Hg, 204Pb, ...
Agilent	1	先试跳 3 行，失败再试跳 2 行；列名按"纯数字"匹配（204Pb 匹配 204）
Thermo	2	跳过 7 行，列名按原名匹配

这一段经过一次重构，从"把整个文件读成字符串矩阵再逐列转换"改成"先读表头定位列，再用 np.loadtxt(usecols=...) 只读需要的那 8 列并直接按 float 解析"。对 3250 行的 Agilent 文件，从 39 ms 降到 9.5 ms（4.1 倍）。

这不是科学改动，是纯工程改动：输入输出仍然逐位相同（24 个数组、0 处不一致）。

新手最容易踩的坑：输入目录里不能有别的东西。 程序会把目录里所有文件都当数据文件处理，多一个无关的 CSV 就报错。

4.4 第 2 步：划分积分窗口

激光剥蚀是一个时间过程：先开激光打样品（信号段），然后关激光等背景（背景段），或者反过来。

界面上选中一个样品会画出信号曲线（showplt，第 1761 行），然后你框选两段：

• 背景段（background）：激光没打样品时测到的计数。它不是零，因为等离子体本身、载气里的杂质都会产生信号。
• 信号段（signal）：激光打样品时的计数。

为什么要扣背景：因为质谱仪测到的是"样品信号 + 环境底色"。你要的是样品那部分。

积分窗口存在 Time_setting.npy 里，dataprocess() 开头（第 1409-1417 行）有一个自检：

if numn[name][3] - numn[name][2] < 0 or numn[name][4] - numn[name][2] < 0:
    wrong_samples.append(name)          # 窗口越界 → 记入错误列表
elif numn[name][3] + 1 > numn[name][4] - 1:
    numn[name][3] = numn[name][4]       # 首尾太近 → 强行对齐
    numn[name][4] = numn[name][3]

这是一个"带病运行"的设计：窗口画错了，程序不中断，而是把一个明显不合理的窗口（起点=终点）硬塞进去算，最后给你一个垃圾结果。所以一定要看 wrong_samples 相关的日志信息。

4.5 第 3 步：背景扣除与比值（dataprocess + reduce_sample）

背景扣除

_net_signal()（第 1561 行）：

'n6': Pb206_s - np.average(Pb206_b),     # 净 206Pb 信号
'n7': Pb207_s - np.average(Pb207_b),
'n8': Pb208_s - np.average(Pb208_b),
'n2': Th232_s - np.average(Th232_b),
'nU': U238_s - np.average(U238_b),
'pb204': np.average(Hg204_s) - np.average(Hg204_b),   # 204 通道

注意一个细节：前四个是"逐点数组减一个标量"（信号段每个时间点减去背景段均值），而 pb204 是"先各自求平均再相减"（得到一个标量）。这个不一致是从原始代码继承下来的，后面会带来麻烦（第 4.7 节的 r84）。

六个比值

_ratios_sample()（第 1588 行，不校正的基准版本）：

比值	表达式	输出列
r76	n7 / n6	²⁰⁷Pb/²⁰⁶Pb
r68	n6 / nU	²⁰⁶Pb/²³⁸U
r82	n8 / n2	²⁰⁸Pb/²³²Th
r86	n8 / n6	²⁰⁸Pb/²⁰⁶Pb
r26	Th232_s / Pb206_s	²³²Th/²⁰⁶Pb（注意用的是原始信号，不是净信号）
r84	n8 / pb204	²⁰⁸Pb/²⁰⁴Pb（分母是标量）

r26 和 r84 这两处的"不走寻常路"是原样保留的。改造时必须照抄，否则结果会变。

2σ 离群过滤

激光剥蚀的信号常有几个时间点被包裹体、裂纹、或者等离子体抖动污染。代码用最朴素的方法剔除（_filter_2s，第 1525 行）：

return np.where(np.abs(x - np.average(x)) < 2.0 * np.std(x), x, nan)

把偏离均值超过 2 倍标准差的时间点替换成 nan（注意：不是删掉，是标记）。

这里有一个只有做过数值对照才会知道的细节：np.nanmean / np.nanstd 在 NumPy 里是把 nan 替换成 0 之后在整个原长数组上求和，而不是"压缩掉无效元素再求"。两种做法的浮点求和顺序不同，末位会漂。代码忠实复刻了 NumPy 的行为（第 1537 行、第 1533-1535 行的注释），所以不要"顺手优化"成 x[~isnan].mean() —— 那会让结果与原来不一致。

输出两个文件

• Mean_Cps.csv：背景段 7 通道均值 + 信号段 7 通道均值（共 16 个数）
• result_all.csv：23 列的比值与计数（列定义见 4.7 节）

4.6 第 4 步：分馏与基体效应校正（SSB）

为什么需要这一步

激光剥蚀出来的物质，轻的同位素比重同位素更容易被带走（分馏）。而且随着剥蚀坑越打越深，坑底的成分会变，分馏也会变（这叫 downhole fractionation）。

结果就是：你测到的 ²⁰⁶Pb/²³⁸U 不等于样品真实的 ²⁰⁶Pb/²³⁸U，它有一个系统性的偏移。

怎么办：Sample-Standard-Bracketing

思路很朴素：

代码里就是一行（factor，第 384 行）：

def factor(stand_v, i):
    return stand_v / date_all[i][standard].mean()

• stand_v：标样理论比值（由标样年龄算出）
• date_all[i][standard].mean()：标样实测比值的平均
• 相除得到分馏因子 f

然后（第 441-448 行）：

result_cal_207Pb_206Pb = date_all[0] * f207_206      # 注意：这是数组整体相乘
result_cal_206Pb_238U  = date_all[2] * f206_238
result_cal_207Pb_235U  = date_all[4] * f207_235
result_cal_208Pb_232Th = date_all[6] * f208_232

这就是核心：样品实测值 × f = 校正后真值。

Average 与 Liner 的区别

	Average（Age_Calculate_average）	Liner（Age_Calculate）
分馏因子来源	整批标样求一个平均值	相邻两个标样之间分段求值
代码	第 418-430 行	第 1024 行的 factor(SEP, ...)
能否校正时间漂移	不能	能（分段近似）
适用	仪器很稳定、批次短	批次长、仪器有漂移

Liner 的实现思路（第 961-992 行）：把标样的位置序列切成段，每一段标样"夹住"一组样品，这组样品就用这段标样的平均因子。

注意：Liner 是分段，不是插值，更不是回归。它把两个标样之间的所有样品当成同一段处理，中间不做时间趋势外推。

代码里其实写了 4 种回归模型（regression()，第 152 行：线性、对数、开方、幂函数，用 scipy.curve_fit 拟合），但这个函数从未被调用 —— 全局变量 regression_method 只在 main() 里被赋值（第 2998 行 regression_method = 1），界面上没有任何控件能改它。这是一个"写了但没接线"的功能，不是死代码，第 6 章会讨论怎么把它接上。

4.7 第 5 步：算年龄

四个年龄

（Age_Calculate_average，第 470-476 行）

输出列	公式	代码
²⁰⁶Pb/²³⁸U 年龄	ln(r68_corr + 1) / λ₂₃₈	第 470 行
²⁰⁷Pb/²³⁵U 年龄	ln(r75_corr + 1) / λ₂₃₅	第 472 行
²⁰⁸Pb/²³²Th 年龄	ln(r82_corr + 1) / λ₂₃₂	第 474 行
²⁰⁷Pb/²⁰⁶Pb 年龄	迭代求解（Age76Pb）	第 476 行

误差传播与那个 ρ

这一步值得仔细看，因为它决定了你报告的误差有没有意义。

第 455-461 行：

w = sqrt((2s_r68 / r68)² + excess_V²) * 100       # 206/238 的相对误差（%）
s = sqrt((2s_r75 / r75)² + excess_V²) * 100       # 207/235 的相对误差（%）
o = sqrt((2s_r76 / r76)² + 0.0001) * 100          # 207/206 的相对误差（%）
if abs((w² + s² − o²) / (2ws)) < w/s:
    rho = abs((w² + s² − o²) / (2ws))
else:
    rho = w/s

三件事：

第一，excess_V 是"额外误差"。 它是界面上的一个输入框（第 3533 行），默认值 3（第 3516 行 var5.set(3)），单位是百分比。它代表"测量本身之外的不确定性" —— 标样的长期重现性、基体差异、仪器日间波动等等。

这是整个软件里最需要你自己判断的一个数字。 默认 3% 是一个偏保守的通用值。如果你这批数据的标样重现性只有 1%，那么用 3% 会高估误差；如果标样重现性有 5%，用 3% 就是低估，审稿人一眼就能看出来。

第二，²⁰⁷Pb/²⁰⁶Pb 那一项用的是硬编码的 0.0001（等于 1% 的平方），不是 excess_V²。这是原始代码里的不一致，属于"原样保留"。实际影响：不管你把 excess_V 设成多少，207/206 的额外误差永远按 1% 算。

第三，rho 是 ²⁰⁶/²³⁸ 与 ²⁰⁷/²³⁵ 两个误差的相关系数。 为什么需要它？因为：

²⁰⁷Pb/²⁰⁶Pb ≈ (²⁰⁷Pb/²³⁵U) / (²⁰⁶Pb/²³⁸U) × (常数)

三个比值不是独立的，它们的误差满足 o² ≈ w² + s² − 2ρws。反解就得到代码里那个式子。这个 ρ 必须输出给画协和图的软件（IsoplotR 等），否则在图上画误差椭圆时会把两个不独立的误差当成独立的，椭圆画错，年龄也跟着错。

输出表里第 40 列就是 rho（第 879 行），这个设计是对的。

4.8 输出文件清单

文件	内容	什么时候产生
result_all.csv	每样品 23 列比值与计数	第 (8) 步 Fractionation Correction
Mean_Cps.csv	每样品 16 个通道均值	同上
cal_age_result_A*.xls	年龄表（Average）	第 (9) 步 Export (Average)
cal_age_result_L*.xls	年龄表（Liner）	第 (9) 步 Export (Liner)
信号图	每个样品选中时实时生成	第 (7) 步

年龄表（xls）的列很密，值得说明（Age_Calculate_average 第 806-888 行定义）：

列号	内容
0-1	文件名、样品名
2-15	六个比值 + 各自的 2s
17-18	²⁰⁶Pb/²³⁸U 年龄 + 2s
19-20	²⁰⁷Pb/²³⁵U 年龄 + 2s
21-22	²⁰⁸Pb/²³²Th 年龄 + 2s
23-24	²⁰⁷Pb/²⁰⁶Pb 年龄 + 2s
25-29	校正后年龄（按所选普通铅路径）
31-34	Tera-Wasserburg 图用：²³⁸U/²⁰⁶Pb + 2s、²⁰⁷Pb/²⁰⁶Pb + 2s
36-40	作图用：²⁰⁷Pb/²³⁵U、²⁰⁶Pb/²³⁸U、各自的 2s、ρ
42-45	²⁰⁸Pb/²³²Th、²⁰⁸Pb/²⁰⁶Pb + 2s
46-49	微量元素 U、Th、Pb（ppm）

31-34 这四列是这套软件最有价值的输出之一，第 5 章会讲怎么用。

关于微量元素换算

第 352-379 行有一个 coefficient() 函数，用 NIST 玻璃标样的已知含量把 cps 换算成 ppm：

标样	U (ppm)	Th (ppm)	Pb (ppm)
NIST610	461.5	457.2	426
NIST612	37.38	37.79	38.57
NIST614	0.832	0.748	2.32

用的是 系数 = 标样已知含量 / 标样实测 cps 均值，然后样品 cps × 系数。

注意：这里假设"微量元素的外标"和"U-Pb 的标样"是两回事 —— 微量元素走 NIST 玻璃，U-Pb 走矿物标样。如果两者混用，换算出来的 U、Th、Pb 含量只能当半定量看。

U 含量这一列（第 18 列）还有一个用处：校验数据合理性。如果你测的是磷灰石，U 含量算出来是 5000 ppm，那大概率哪里出问题了。

---

第 5 章　怎么让别人接受你的结果

这一章是实用的。审稿人和同行不会因为你用了某个软件就相信你，他们会问三个问题。

5.1 审稿人会问的三个问题

1. "你怎么知道普通铅扣对了？"
2. "你怎么知道仪器分馏校正对了？"
3. "你报的误差是内部误差还是真实可重复性？"

三个问题各有对策。

5.2 对策一：把你的标样表现亮出来

你的结果的可信度上限，就是你这批数据里标样的表现。

因为 SSB 的逻辑是"用标样校准样品"，如果标样自己在这批数据里就飘了 5%，那你不可能比 5% 更准。

好消息是：代码本来就已经算了这些量（第 423、428 行）：

std_68_2s = date_all[2][standard].std() / date_all[2][standard].mean()
std_75_2s = date_all[4][standard].std() / date_all[4][standard].mean()

但算完就丢掉了，没有输出到任何文件、也没有打印。这是软件的一个实际缺陷（第 6 章列为改进项）。

你可以怎么补救：

• 打开生成的 xls，找标样的行，看"²⁰⁶Pb/²³⁸U 年龄"那一列。如果标样名义年龄 483 Ma，而你这批数据里的标样年龄在 470–495 Ma 之间波动，那么你的外部误差就不该小于 ±2.5%。
• 论文里建议明确写一句："本批数据的标样 MAD-NEW 重复测量 ²⁰⁶Pb/²³⁸U 年龄为 483 ± 12 Ma（n=8，2SD），与推荐值一致。" —— 这一句话能省掉很多质疑。

5.3 对策二：用 Tera-Wasserburg 图自证（最重要的一招）

这是最能说服审稿人的一步，而软件已经把需要的数据准备好了。

什么是 Tera-Wasserburg 图

一种特殊的等时线图：

• 横轴：²³⁸U/²⁰⁶Pb（注意是倒过来的）
• 纵轴：²⁰⁷Pb/²⁰⁶Pb

把所有数据点画上去，如果它们来自同一个年龄、同一批普通铅，它们会大致落在一条直线上。这条直线：

• 下截距（y 轴上，即 ²³⁸U/²⁰⁶Pb → 0 的位置）= 普通铅的 ²⁰⁷Pb/²⁰⁶Pb
• 上截距与协和线的交点 = 年龄

为什么它能解决问题

因为它把"普通铅是多少"从假设变成了测量结果。

前面第 3 章说过：S&K 模型给的普通铅组成只是"平均猜测"。而在 TW 图上，你用的是"数据自己拟合出来的截距"，不需要假设普通铅是什么。

这就直接回答了审稿人的第一个问题：不是"我相信模型铅"，而是"我让 20 个数据点自己告诉我普通铅长什么样，而且它们线性度很好（MSWD 小），说明这个假设是自洽的"。

怎么操作

1. 从 cal_age_result_*.xls 里取 31-34 列（²³⁸U/²⁰⁶Pb + 2s、²⁰⁷Pb/²⁰⁶Pb + 2s）；
2. 也可以取 36-40 列（²⁰⁷Pb/²³⁵U、²⁰⁶Pb/²³⁸U、ρ）画普通协和图；
3. 把这些数据贴进 IsoplotR（在线版免费：isoplotr.london-geochron.com）或者 IsoplotR 的 R 包；
4. 选择 "Tera-Wasserburg" 模式，让它拟合下截距和上截距。

判据：

• 拟合出的下截距应该和你选的普通铅 ²⁰⁷/²⁰⁶ 差不离。如果差很多，说明你的假设错了 —— 这是好事，你提前发现了问题。
• MSWD 应该在 1 附近（20 个点的话，1.5 以内都算可以）。MSWD 太大说明有一期以上的普通铅，或者有铅丢失，或者数据里有离群点。

这一招还有一个额外好处：如果下截距显著不等于模型值，你可以直接把 TW 图拟合出的截距当作普通铅组成输入回 Isoclock 重算一遍。这就形成了一个闭环。

5.4 对策三：把外部误差算对

回顾 4.7 节：误差 = 内部误差 ⊕ 外部误差（excess_V）。

• 内部误差：来自这一次测量的计数统计。信号越强、积分时间越长，它越小。它是"这个点的精度"。
• 外部误差：来自两次测量的差异。它包含标样不确定性、基体差异、仪器漂移。它是"这个点的准确度"。

审稿人对只报内部误差的结果通常是不接受的，因为那会让误差条小得离谱（往往 ±0.5% 都不到），而实际上不同实验室测同一块标样能差 2–5%。

怎么定 excess_V：

情况	建议
有长期的标样重现性数据（同一实验室、同一仪器）	用那个 2SD 值
标样有推荐不确定度（如 91500 的 ²⁰⁶Pb/²³⁸U 推荐值有 ±1% 左右）	用它
两者都有	取大的，或者按误差传播合成
都没有	保守取 2–5%，并在文中说明

并且一定要在论文里写明白："年龄不确定度含外部误差 X%（依据为 …）"。含含糊糊地只报一个 2s，最容易被抓。

5.5 对策四：做敏感性检验（换一条路再跑一遍）

这是最省力、也最有力的稳健性检验。

具体做法：同一批数据，用不同的普通铅校正路径各跑一遍：

跑法	界面设置
A（主结果）	样品 Cal 204Pb，标样 Cal 204Pb(Std)，Standard Corr.
B	样品 207Pb，标样 207Pb(Std)，Standard Corr.
C	样品 208Pb，标样 208Pb(Std)，Standard Corr.
D（对照）	样品 None，Standard not Corr. —— 看不校正会偏多少

然后对比：

• 如果 A、B、C 的年龄在误差范围内一致 → 非常好，说明你的结果不依赖校正方法的选择，可以直接写进论文，附上一条"三种校正路径结果一致"的说明。
• 如果它们差很多（比如差 30 Ma，超过误差）→ 这是一个警示信号。说明你的普通铅组成假设不稳定，或者数据本身有问题。这时候不要挑一个最好看的报，要做的是回头查原因（是不是有多个普通铅端元？TW 图上是不是散开了？）。
• D 与 A 的差值，就是你这次校正到底改了多少。这个差值应该和你样品的 ²⁰⁶Pb/²⁰⁴Pb 反向相关 —— 如果不对，逻辑上就有问题。

这个检验不花什么时间（软件跑得很快），但它在审稿意见里非常有说服力。 建议直接做成附表。

5.6 一份"可以直接抄"的附录清单

论文方法部分或附录里，建议把下面这些写清楚。少写一条，就可能被追问一轮。

【数据采集】
- 仪器型号、激光系统、剥蚀策略（光斑、频率、能量密度）
- 采集的气体、载气流量
- 每次分析的前冲洗/后冲洗时间

【标样】
- 标样名称（如 MAD-NEW）
- 标样所用年龄及其来源文献（如 483 Ma, 某文献）
- 本批数据中标样的重复性（n、均值、2SD）

【数据处理】
- 软件名称与版本（Isoclock v2.0）
- 计算模式：Standard Corr. / Standard not Corr.
- 标样普通铅校正路径：Cal 204Pb(Std) / 207Pb(Std) / ...
- 样品普通铅校正路径：Cal 204Pb / 207Pb / ...
- 普通铅参数的来源：S&K 模型（并给出所用年龄）或实测输入值
- excess_V 的取值与依据
- 是否做了 Hg 校正（仪器自动 / 未做）

【结果校验】
- Tera-Wasserburg 图的下截距与 MSWD
- 不同校正路径的结果对比（附表）
- 标样与参考值的偏差

5.7 一句提醒

不要为了得到"好看的年龄"而去挑校正方法。 这不仅是学术道德问题，也是实际风险问题 —— 挑出来的结果在后续被别的方法交叉验证时很容易被发现。

换句话说：第 5.5 节的敏感性检验，应该在你决定最终报哪个数之前做，而不是之后。

---

第 6 章　局限与可以改进的方向

这一章我尽量说全，并且按"改了会不会影响已有结果"来分类。这是诚实的部分 —— 一个软件的真实能力边界，比它的功能列表更重要。

6.1 三条需要留意的数值行为（已写入 README 的 Notes for users）

缺陷 1：Age76Pb 的收敛判据失效

位置：第 122-145 行，特别是第 133 行。

Rap = (exp(0.00098485 * Tm) - 1) / (exp(0.000155125 * Tm) - 1) / 137.818   # 第 130 行，循环外

while abs(Rap76 - Rap) > 0.00005 and N < 10:    # 第 133 行
    ...
    Tm = Tmin + (Tmax - Tmin) * (Rap76 - Rapi) / (Raps - Rapi)
    Age76Pb = Tm
    N = N + 1                                    # Rap 从未更新

Rap 在进入循环前算了一次，循环体内没有任何地方重新计算它。所以 abs(Rap76 - Rap) 的值在整个循环里是恒定的 —— 条件要么一开始就是假（循环不执行），要么永远为真（循环跑满 10 次）。

后果：迭代次数被写死成 10。实测（把这段源码逐字复制出来，与二分法高精度解对比，输入为各年龄对应的理论 ²⁰⁷Pb/²⁰⁶Pb）：

目标年龄	Age76Pb 返回	精确解	偏差
100 Ma	100.00	100.00	−0.00
500 Ma	500.00	500.00	+0.00
1000 Ma	1000.03	1000.00	+0.03
1500 Ma	1502.26	1500.00	+2.26（+0.15%）
2100 Ma	2134.77	2100.00	+34.77（+1.66%）
2500 Ma	2499.31	2500.00	−0.69

这组数字与 README 里"约 +35 Ma（1.7%）@ 2.1 Ga"完全吻合。值得注意的是，误差不是单调的：1.5 Ga 以下几乎为零，2.1 Ga 附近最差，2.5 Ga 又回落到零 —— 这正是"迭代次数被写死"的典型症状（不同年龄需要不同的迭代次数才能收敛，10 次恰好对某些区间够、对另一些不够）。

影响范围：只影响 ²⁰⁷Pb/²⁰⁶Pb 年龄这一列（输出第 23-24 列）。²⁰⁶Pb/²³⁸U、²⁰⁷Pb/²³⁵U 年龄有解析解，不受影响。所以对年轻样品（< 1 Ga）来说，这条缺陷实际可以忽略；只在老样品里才需要警惕。

README 建议 1.5–2.2 Ga 的年龄谨慎使用。

修复难度：极低。把 Rap 的更新移进循环即可。但这会改变数值输出，所以需要作者确认。

更稳妥的替代：直接改用 SciPy 的 brentq 求根，几十年内不会再有收敛问题。

缺陷 2：²⁰⁴Pb 路径的除零

位置：_ratios_std204，第 1639-1640 行。

m74 = (...) / net['pb204']
m64 = (...) / net['pb204']

如果某个时间点的净 ²⁰⁴ 信号为 0（或整段 204 通道的背景恰好等于信号），这里直接 ZeroDivisionError，整个批次的处理中断。

后果：不是算错，是跑不完。所以危害等级反而比缺陷 1 低（至少它不会骗你）。

修复难度：低。加一个保护即可，但要想清楚"该怎么处理"（跳过该点？该点标记为无效？）—— 这需要作者决定语义。

缺陷 3：²⁰⁴Hg 校正恒为 0

位置：_net_signal，第 1569-1571 行。

# 204 通道：原实现直接用原始 204 减背景，未做 Hg 干扰扣除
# （源码里那段 Hg204cal_s 恒为 0、m204_s 从未被引用）
'pb204': np.average(sig['Hg204_s']) - np.average(sig['Hg204_b']),

原始代码里有一段试图计算 Hg 对 204 的贡献，但表达式是 np.average(x) - np.average(x) —— 恒等于 0，而且它的结果从未被引用。所以软件完全没有做 Hg 校正。

后果：如果仪器导出数据时已经做了 Hg 校正，无害；如果没有，那么：

• ²⁰⁴Pb 被高估 → 普通铅被高估 → 年龄偏年轻（因为扣掉了过多的放射成因铅）
• 208Pb/204Pb 列（输出第 14-15 列）整体偏低

危害等级：中高，因为它安静地给出错误结果，不报错、不警告。

你能做的：去确认仪器导出时到底有没有做 Hg 校正。这是数据层面的问题，软件层面改不了。

修复难度：中。需要 ²⁰⁴Hg/²⁰²Hg 的天然比值（约 0.23）和仪器的质量歧视校正因子，后者需要专门的 Hg 标定。

6.2 方法层面的局限

这些不是 bug，是这套方法本身的边界。

局限 4：标样与样品"同龄假设"的传递

Cal_age(Standard_age) 用标样年龄定出 Q8/R8/S8，供标样校正使用 —— 这是对的。但样品侧的普通铅参数（common_Pb、common_Pb207_204 等）来自用户输入的一个年龄或一组比值，是一个全局常数，对所有样品一视同仁。

问题：如果这批样品跨越了较宽的年龄区间，或者普通铅来源不同，用一个常数去校正所有人，会把系统误差引入到每个点上。

改进方向：允许按样品组（而不是整批）指定普通铅参数；或者做迭代 —— 先用初值算年龄，再用算出的年龄反推普通铅，反复两三次。（这是很多文献里的标准做法，代码里目前没有。）

局限 5：SSB 假设"标样与样品分馏行为相同"

外标法的前提是标样和样品的基体足够接近。MAD-NEW 磷灰石标样配磷灰石样品是合理的；但用 NIST 玻璃做标样去测黑钨矿，两者的剥蚀行为、气溶胶传输效率、电离效率都不同，分馏因子不能直接套用。

这是方法层面的固有局限，不是软件问题。你能做的是：选基体尽量接近的标样，并在论文里说明。

局限 6：Liner 只是分段，不插值、不回归

前面说过，Age_Calculate 把标样分成若干段，每段一个因子。如果两个标样之间隔了 10 个样品，这 10 个样品用的是同一个因子，中间的时间趋势被完全忽略。

regression()（第 152 行）实现了 4 种漂移模型，但没有接到界面上，从未执行。

改进方向：把 regression() 接线（工作量小），或者用更标准的时间插值（对每个样品，用它前后两个标样的因子做线性内插）。后者是当前文献里的普遍做法。

局限 7：没有等时线拟合、没有 MSWD

Isoclock 只做"逐点校正 + 逐点年龄"，不做统计拟合：

• 不给 TW 图的截距和误差
• 不给协和年龄（concordia age）和 MSWD
• 不给加权平均年龄

这意味着你必须把数据导出到 IsoplotR（或其他软件）做第二步。这不是缺陷，但使用者必须知道：Isoclock 的输出是"中间产品"，不是最终结果。

好消息是它输出的列（31-40）恰好是为这个目的设计的，导出很方便。

局限 8：²⁰⁸Pb/²³²Th 年龄通常不可靠

软件会输出这个年龄（第 21-22 列）。但要提醒：Th 和 U 的地球化学行为差异大，²³²Th 的信号强度和稳定性通常不如 ²³⁸U，而且初始 Th/Pb 假设更难满足。

建议：除非你研究的是独居石这类富 Th 矿物，否则把 ²⁰⁸Pb/²³²Th 年龄当作参考，不要作为主要结论。

局限 9：abs() 掩盖了失效

前面第 1.2 节提过：所有年龄计算都套了 abs()。这避免了负数取对数导致崩溃，但也意味着当数据病态时，你得到的不是报错，而是一个看似正常的数字。

判断方法：算完年龄后，回头检查输入的比值。如果某点校正后的 ²⁰⁶Pb/²³⁸U 是负数，或者年龄是 0（Age76Pb 在比值小于 0.0460455 时返回 0，第 128、144 行），那这个点就是无效的，必须剔除。

6.3 工程层面的局限

问题	影响	改进方向
输入/输出目录不能相同	必须开两个目录	加一个检查，或者允许同目录
输入目录里有无关文件就报错	容易踩坑	按扩展名过滤
每批数据必须重启程序	效率低	状态管理重构
输出表靠"列的位置"索引	改动列顺序就全错	引入列名 schema
状态文件用 pickle（Time_setting.npy 等）	安全与兼容性风险	改用 JSON
转换脚本曾用 eval() 解析文件名	任意代码执行风险	已修复（改为 int(float(...))）
异常静默失败	打包成 exe 后尤其严重	已改进（日志带 traceback、Tk 回调异常已接管）
缺自动化测试	改动容易引入回归	本 fork 已补（网页端 67 项测试 + 桌面端对照验证）

6.4 如果只能改三件事

按"性价比"排序（性价比 = 影响 ÷ 改动量）：

1. 修 Age76Pb 的收敛判据。 一行改动，解决一个影响 1.7% 的系统偏差。但会改变输出，需要作者确认。
2. 把标样重复性输出出来。 现在 std_68_2s、std_75_2s 算完就丢。多写两行 print 或加两列到 xls，用户立刻就能看到自己数据的可信度上限。
3. 加一个"可疑值标记"机制。 现在 abs() 和兜底逻辑会把病态点变成正常数字。加一列 flag，标出"比值负数""年龄为 0""204 信号为零"的点。

再往后的优先级：Hg 校正（需要仪器配合）、regression() 接线、按样品组指定普通铅、JSON 状态文件、列名 schema。

6.5 一句话总结局限

---

第 7 章　代码地图：新手该从哪读起

Isoclock2.0.py 有 3676 行，但真正需要理解的科学代码不到 300 行。剩下的是界面代码。

7.1 核心函数索引

函数	行号	职责	建议阅读顺序
SK2model	40	S&K 模式铅（3.7 Ga 起算写法）	6
Age76Pb	48	²⁰⁷Pb/²⁰⁶Pb 年龄（迭代）	7
regression	78	漂移回归模型（未接线）	—
loaddata	199	读三种仪器文件	1
Cal_age	3096	S&K 参数（现代值起算写法）→ Q8/R8/S8/P382/Pbc	5
Age_Calculate_average	244	平均分馏因子 + 算年龄 + 出 xls	8
Age_Calculate	852	分段分馏因子（Liner）	9
dataprocess	1331	主处理流程	2
_filter_2s	1451	2σ 离群标记	3
_mean_2sem	1456	均值与 2 倍标准误	4
_net_signal	1487	背景扣除	3
_ratios_std207	1526	207Pb 校正	10
_ratios_std_cal204	1539	Cal 204Pb 校正	11
_ratios_std204	1555	204Pb 校正	12
_ratios_std208	1573	208Pb 校正	13
reduce_sample	1607	还原成 23 列输出行	14
PbCorrSamples	2809	询问样品普通铅参数	15
instructure0/1/2	1781 / 2064 / 2348	三个界面页	—
main	2885	主界面与全局变量	—

7.2 三条推荐阅读路径

路径 A：想理解科学逻辑（推荐先走这条）

SK2model (114)  →  Cal_age (3170)  →  Age76Pb (122)
       ↓
  _net_signal (1561)  →  _ratios_std_cal204 (1613)
       ↓
  factor (384)  →  Age_Calculate_average 第 418-484 行

这条路上的每个函数都不到 40 行，加起来不到 200 行，但包含了整个方法的核心。

路径 B：想理解数据处理

loaddata (273)  →  dataprocess (1405)  →  _load_channels (1375)
       ↓
  _mean_cps_row (1393)  →  reduce_sample (1681)

重点是看数据从"文件里的数字"变成"比值"的每一步。

路径 C：想改界面

main (2959)  →  instructure0 (1855)  →  showplt (1761)

界面代码占了全文件 80% 以上，都是 Tk 的布局与回调，读的时候可以把里面所有的 Frame、Entry、Button 先跳过，只找 command= 后面接的函数名。

7.3 读这份代码的三个注意事项

1. 全局变量满天飞。 像 standard、Standard_age、Q8、excess_V、PbCorrS 这些都在函数之间通过 global 传递。读的时候先用搜索功能找到它在哪里被赋值。
2. 很多写法是"故意保留的不一致"。 看到两处数学等价但写法不同的代码，不要急着统一 —— 代码注释里通常写明了为什么（浮点路径不同会改变末位）。
3. 留意 abs() 和 except。 这是两个"消灭错误信息"的机制。看到它们，就要提高警惕。

---

第 8 章　术语速查表

术语	英文	一句话解释
协和	concordant	²⁰⁶Pb/²³⁸U 和 ²⁰⁷Pb/²³⁵U 给出的年龄一致
不协和	discordant	两者不一致，通常是放射成因铅丢失
协和线	concordia	所有协和点连成的曲线，是年龄的"标尺"
普通铅	common Pb (Pb_c)	矿物结晶时混进来的铅，不含年龄信息
放射成因铅	radiogenic Pb (Pb*)	由本矿物内 U/Th 衰变产生的铅，含年龄信息
初始铅	initial Pb	同普通铅，强调"结晶时刻"
模式铅	model Pb	用地球演化模型算出的理论普通铅组成
μ	mu	²³⁸U/²⁰⁴Pb 比值，S&K 模型取约 9.74
κ	kappa	²³²Th/²³⁸U 比值，S&K 模型取约 3.78
分馏	fractionation	仪器造成同位素比值偏离真值的系统偏差
剥蚀坑分馏	downhole fractionation	激光越打越深导致的分馏变化
SSB	Sample-Standard-Bracketing	用已知年龄标样校准样品的方法
外标	external standard	与样品分开测量的校准物质
Tera-Wasserburg	TW 图	横轴 ²³⁸U/²⁰⁶Pb、纵轴 ²⁰⁷Pb/²⁰⁶Pb 的等时线图
下截距	lower intercept	TW 图上给出普通铅 ²⁰⁷/²⁰⁶ 的截距
2s	2 sigma	2 倍标准差，约 95% 置信水平
内部误差	internal error	单次测量的计数统计误差
外部误差	external error	重复测量之间的重现性误差
MSWD	mean square weighted deviation	衡量数据点与拟合线一致程度的指标，1 附近最好
ρ (rho)	error correlation	²⁰⁶Pb/²³⁸U 与 ²⁰⁷Pb/²³⁵U 误差的相关系数
cps	counts per second	计数率，质谱仪原始输出单位

---

附录　参考文献

本文中提到的所有方法出处：

• Chew, D.M., Petrus, J.A., Kamber, B.S., 2014. U–Pb LA–ICPMS dating using accessory mineral standards with variable common Pb. Chemical Geology 363, 185–199.

→ 对应代码里的 Cal 204Pb 法（_ratios_std_cal204）。

• Stacey, J.S., Kramers, J.D., 1975. Approximation of terrestrial lead isotope evolution by a two-stage model. Earth and Planetary Science Letters 26, 207–221.

→ 对应代码里的 SK2model 与 Cal_age 中的全部模式铅常数。

• Zack, T., Stockli, D.F., Lander, R.H., McRae, C., Stockli, L.D., 2011. In situ U-Pb rutile dating by LA-ICP-MS: ²⁰⁸Pb correction and prospects for geological applications. Contributions to Mineralogy and Petrology 162, 515–530.

→ 对应代码里的 208Pb 法（_ratios_std208）。

• Tera, F., Wasserburg, G.J., 1972. U-Th-Pb systematics in three Apollo 14 basalts and the problem of initial Pb in lunar rocks. Earth and Planetary Science Letters 14, 281–304.

→ Tera-Wasserburg 图的原始出处。

• 软件说明书中引用的标样推荐年龄（AY-4 = 154 Ma、MAD-NEW ≈ 483 Ma、91500 = 1062 Ma、wc-1 = 254.4 Ma）请以各标样原始文献为准，代码中的数值只是默认建议值。

---

最后

这份文档的写法，是尽量把"为什么"讲清楚，而不是只讲"怎么点按钮"。

因为 U-Pb 定年里最容易出问题的从来不是操作 —— 是操作背后的那一步假设。软件把假设做成了默认值，默认值用起来毫无摩擦，于是大家就忘了它们存在。

记住三件事就够了：

1. 普通铅的参数不是常数，是你告诉软件的一个假设；
2. 模型铅只是平均猜测，能用实测就用实测，能用等时线就用等时线；
3. 报年龄的时候，顺便报一下你对普通铅做了哪些假设、以及换一种假设年龄会怎么变。

做到第三条，你的结果就比大多数已发表的结果更站得住。