碳化硅外延层厚度的确定（2025 国赛 B）优秀范文（三）

摘要

题目要求由 400–1000 nm 波段测得的反射率谱 $R(\lambda)$ 反演 4H-SiC 外延层厚度 $d$，并量化结果的不确定度。本题范文一（时域最小二乘点估计）与范文二（波数域频域点估计）已分别给出单一"最优"厚度值及其误差棒；本文换一个互补的视角——贝叶斯反演，不再只追求一个点估计，而是对参数 $(d,A)$ 建立完整的后验分布，刻画 $d$–$A$ 强耦合带来的不确定性结构，并给出可信区间与先验敏感性。我们复用范文一的三层薄膜双程干涉正演模型，构造高斯噪声似然与弱信息正态先验，用自适应 Metropolis-Hastings（Haario 协方差自适应）采样 80000 次得到后验。最终后验均值 $d=3.1968\ \mu\text{m}$、$A=2.5708$，95% 可信区间 $d\in[3.1880,3.2066]\ \mu\text{m}$（半宽 $0.00932\ \mu\text{m}$）；$d$–$A$ 后验相关系数高达 $0.9977$，从统计上再次印证联合反演的病态本质。改变 $A$ 的先验宽度（$\sigma_A$ 从 $0.50$ 到 $0.02$），后验均值 $d$ 仅在 $3.1970\sim3.1972$ 间波动，说明结果由数据（似然）主导、对主观先验稳健；剖面似然与边缘后验峰值一致，表明频率派与贝叶斯结论互通；后验预测带完整覆盖实测曲线。本文与范文一、二共同构成"点估计—频域交叉验证—贝叶斯不确定度"三位一体的厚度测定方案。

一、问题重述

给定一片在衬底上外延生长的 4H-SiC 外延片，已测得反射率随波长变化的曲线 $R(\lambda)$（波段 400–1000 nm，201 个采样点，步长 3 nm）。衬底折射率 $n_2=2.70$、空气折射率 $n_0=1$ 视为已知；外延层折射率随波长按 Cauchy 式 $n_1(\lambda)=A+B/\lambda^2$ 变化（$B$ 固定为文献先验 $B_0=0.015$）。需要由该光谱反演外延层厚度 $d$，并给出结果的精度与不确定度。问题本质是由观测光谱反演薄膜光学参数的不适定反问题：干涉条纹的周期唯一确定 $d$，但折射率色散未知时 $d$ 与 $n_1$ 间存在强耦合，单一波长拟合会陷入病态；贝叶斯方法的长处恰恰在于把这种耦合显式表达为后验的协方差结构，而不是强行忽略它。

二、模型假设

1. 单层均匀膜：外延层厚度均匀、各向同性，折射率仅随波长变化。
2. 界面理想、无吸收：两界面平整无散射，所考波段内外延层与衬底均无吸收，反射仅由折射率失配引起。
3. 色散形式：Cauchy 经验式 $n_1(\lambda)=A+B/\lambda^2$（$\lambda$ 取微米）。
4. 衬底/空气折射率已知：$n_0=1.0$、$n_2=2.70$ 作为固定先验。
5. 测量噪声模型：实测反射率由无噪模型叠加加性高斯白噪声得到，标准差 $\sigma=0.006$。
6. 贝叶斯设定：似然取高斯噪声模型；先验取弱信息正态先验 $d\sim\mathcal N(3.5,1.0^2)$、$A\sim\mathcal N(2.57,0.05^2)$，仅施加覆盖物理合理范围的温和约束，避免引入强主观偏向。

三、符号说明

符号	含义	取值/单位
$\lambda$	波长	400–1000 nm
$n_0,n_1,n_2$	空气/外延层/衬底折射率	1.0 / 未知 / 2.70
$d$	外延层厚度（反演目标）	$\mu\text{m}$
$A,B$	Cauchy 色散常数	$A\approx2.57,\ B=0.015$
$\delta$	单程相位厚度	$\delta=2\pi n_1 d/\lambda$
$R$	总反射率	$R=|r|^2\in[0,1]$
$\sigma$	测量噪声标准差	0.006
$p(d,A\mid R)$	参数后验分布	贝叶斯推断目标

四、模型建立

4.1 三层薄膜正演模型

体系视为空气 ($n_0$)/外延层 ($n_1$)/衬底 ($n_2$) 三层结构。垂直入射时上下界面的菲涅尔振幅反射系数分别为

$$r_{01}=\frac{n_0-n_1}{n_0+n_1},\qquad r_{12}=\frac{n_1-n_2}{n_1+n_2}.$$

外延层引入单程相位厚度 $\delta=2\pi n_1(\lambda)d/\lambda$，由特征矩阵得总反射振幅

$$r=\frac{r_{01}+r_{12}e^{-2i\delta}}{1+r_{01}r_{12}e^{-2i\delta}},$$

总反射率 $R(\lambda)=|r|^2$。代入真值 $d=3.20\ \mu\text{m}$、$A=2.57$、$B=0.015$，模型反射率在 $[0.1806,0.2111]$ 间呈现规则干涉条纹，与实测形态一致，正演模型成立。

4.2 病态本质：$d$–$A$ 强耦合

由于色散常数 $A$ 同时进入 $n_1$ 而 $n_1$ 又进入相位 $\delta$，$d$ 与 $A$ 并非独立可辨：略微改变 $A$ 并相应调整 $d$ 可使反射谱几乎不变。这正是联合反演 $(d,A)$ 病态的根源。在贝叶斯语言里，这一耦合表现为后验在 $(d,A)$ 平面上沿一条倾斜"脊"分布——沿脊方向参数不确定，而垂直于脊的方向被数据牢牢锁定。理解这一结构，是选择采样策略与解释结果的关键。

4.3 贝叶斯框架：似然 + 先验 + 后验

给定实测谱 $\hat R_i$，高斯噪声似然为

$$p(\hat R\mid d,A)\propto\exp\!\Bigl[-\frac{1}{2\sigma^2}\sum_i\bigl(R_{\text{model}}(\lambda_i;d,A,B_0)-\hat R_i\bigr)^2\Bigr].$$

联合弱信息正态先验 $p(d,A)$，由贝叶斯定理得后验

$$p(d,A\mid \hat R)\propto p(\hat R\mid d,A)\,p(d,A).$$

与逐点优化的方法不同，后验同时编码了"最可能的值"和"可能的散布"，可直接读出可信区间，而无需额外做蒙特卡洛扰动实验。

4.4 自适应 Metropolis-Hastings 采样

后验无解析形式，采用 Haario 协方差自适应的随机游走 Metropolis-Hastings：以时域最小二乘解 $(d_{\text{LS}},A_{\text{LS}})$ 为初值（避免沿脊滑入相邻周期别名 basin），前 100 步用较小固定协方差保证混合，之后用历史样本的经验协方差 $\hat\Sigma$ 构造 proposals，步长取 $(2.38^2/\dim)\hat\Sigma$，使采样自动沿脊方向拉伸、克服强耦合导致的混合困难。共迭代 80000 次，舍去前 12000 次燃烧期，按间隔 68 薄化，保留 1000 个有效后验样本。

五、模型求解与结果分析

5.1 收敛诊断

抽样轨迹（图 1）显示 $d$ 与 $A$ 两条链在燃烧期后迅速平稳，无趋势性漂移或周期跳变，表明采样已充分遍历后验。参数平面上的后验散点（图 2）呈明显的倾斜脊状结构，$d$ 与 $A$ 沿脊强相关；最大后验点 MAP $(3.1972,2.5705)$ 与时域最小二乘点 LS $(3.1971,2.5706)$ 几乎重合，说明三种独立方法得到的点估计彼此一致，脊的方向正是病态耦合的几何体现。

5.2 后验推断

厚度 $d$ 的边缘后验（图 3）以真值 $3.20$ 为对称中心，峰值与后验均值 $3.1968$、MAP $3.1972$ 高度吻合，95% 可信区间 $[3.1880,3.2066]\ \mu\text{m}$ 完整覆盖真值。色散常数 $A$ 的边缘后验（图 4）被数据压缩到极窄区间，峰值 $2.5708$ 与真值 $2.5700$ 一致——相比先验 $\sigma_A=0.05$，后验方差被大幅收缩，说明数据对 $A$ 的信息远多于先验，所谓"先验被数据强烈修正"正体现于此。

5.3 频率派与贝叶斯的一致性

固定 $d$ 并对 $A$ 做一维最小二乘优化，得到剖面似然（profile likelihood）；它与贝叶斯边缘后验（图 5）在 $d\approx3.197$ 处几乎重合于同一峰值，真值竖线落在峰内。这从方法论上说明：当似然足够强、先验足够弱时，贝叶斯后验与点估计的置信结构殊途同归，二者并非对立，而是同一信息的两种表述。

5.4 三方法反演结果对比

将本文贝叶斯后验均值与范文一时域最小二乘、范文二频域相位的结果并列（图 7）：时域 LS $d=3.1971\pm0.0066$、频域 $d=3.2001\pm0.0087$、贝叶斯 $d=3.1968\pm0.0093$（± 为 95% 半宽），三者误差棒高度重叠且全部覆盖真值 $3.20\ \mu\text{m}$。三套独立数学框架给出彼此印证的解，是该厚度测定方案可靠性的最强证据。

六、灵敏度分析

6.1 先验敏感性

贝叶斯结果是否受主观先验左右，是评审最关心的质疑。我们取三组 $A$ 的先验宽度（$\sigma_A=0.50,0.05,0.02$）分别重跑 MCMC，后验均值 $d$ 依次为 $3.1970$、$3.1972$、$3.1970\ \mu\text{m}$（图 6），相对变化不足 $0.01\%$。这说明 $d$ 的推断几乎完全由似然（数据）决定，弱信息先验仅起定义域约束作用，结论对先验选择稳健。

6.2 后验预测能力

用 1000 个后验样本在每条波长上预测反射率，取 $\pm1\sigma$ 形成后验预测带（图 8）。该带完整覆盖实测曲线，且后验均值预测（绿虚线）与实测（红）几乎贴合，表明后验不仅给出参数不确定度，也正确传播为可观测量的预测不确定度——模型无系统遗漏。若预测带明显偏窄而实测逸出，则提示噪声模型或正演假设需修正。

七、模型评价

优点：(1) 提供完整的后验分布而非单点值，天然给出可信区间，不确定度量化更直接；(2) 自适应 MCMC 沿脊采样，显式刻画 $d$–$A$ 强耦合结构，物理含义清晰；(3) 先验敏感性分析证明结果由数据主导，对主观设定稳健；(4) 剖面似然与边缘后验一致，桥接频率派与贝叶斯两套话语；(5) 后验预测带提供对实测的覆盖检验，模型可信度高；(6) 全部计算仅用标准库，可完整复现。

局限：(1) MCMC 采样较点估计耗时，需合理设置燃烧期与薄化以保证有效样本量；(2) 仍依赖衬底折射率 $n_2$ 与色散常数 $B$ 的先验，若 $B$ 未知需与 $(d,A)$ 联合推断，耦合进一步加剧；(3) 假设无吸收、界面理想，对吸收性或高粗糙样品需引入损耗项；(4) 当条纹稀少（外延层过薄）时后验可能变宽乃至多峰，需补充光谱范围或先验信息。

八、结论

本文以贝叶斯反演视角解决 4H-SiC 外延层厚度测定问题，在复用三层干涉正演的基础上，用自适应 Metropolis-Hastings 对 $(d,A)$ 采样得到后验。后验均值 $d=3.1968\ \mu\text{m}$、95% 可信区间 $[3.1880,3.2066]\ \mu\text{m}$，与真值 $3.20\ \mu\text{m}$ 一致；$d$–$A$ 后验相关 $0.9977$ 从统计上揭示了联合反演的病态根源。先验敏感性（后验均值 $d$ 在 $3.1970\sim3.1972$ 间）与三方法互证（时域/频域/贝叶斯全部覆盖真值）共同表明结论稳健、数据主导。相比范文一、二的点估计，本篇补充了"不确定度从何而来、对假设多敏感"的贝叶斯回答，三者合起来构成一套可复现、可辩护、可工程落地的外延层厚度测定方案。

九、参考文献

[1] 全国大学生数学建模竞赛组委会. 2025 高教社杯全国大学生数学建模竞赛 B 题：碳化硅外延层厚度的确定[R]. 2025.

[2] Tompkins H G, Irene E A. Handbook of Ellipsometry[M]. William Andrew, 2005.

[3] Gelman A, Carlin J B, Stern H S, et al. Bayesian Data Analysis[M]. 3rd ed. CRC Press, 2013.

[4] Haario H, Saksman E, Tamminen J. An adaptive Metropolis algorithm[J]. Bernoulli, 2001, 7(2): 223–242.

[5] Born M, Wolf E. Principles of Optics[M]. 7th ed. Cambridge University Press, 1999.

附录：核心 Python 实现

import os, sys, math, random
# 定位并导入真源模块（与配图同源，四路一致）
try:
    _HERE = os.path.dirname(os.path.abspath(__file__))
except NameError:
    _HERE = os.getcwd()
sys.path.insert(0, os.path.abspath(os.path.join(_HERE, "..", "..", "..", "tools")))
import gen_cumcm2025b_3 as G3   # 贝叶斯反演真源

# ---- 复算关键数字（固定种子 → 与正文 / 配图逐字吻合）----
Rmeas, _ = G3.G.build_meas()
x_ls, _, _ = G3.G.fit(Rmeas, G3.G.LAM, seed=2026)          # 时域 LS 定位 basin
samples, acc = G3.mcmc(Rmeas,
                       prior={"d_mu": 3.5, "d_sd": 1.0, "A_mu": 2.57, "A_sd": 0.05},
                       n_iter=80000, burn=12000, thin=68, seed=20260822,
                       init=list(x_ls))
stat = G3.summarize(samples, Rmeas)

print("后验均值 d=%.4f μm   A=%.4f" % (stat["mean_d"], stat["mean_A"]))
print("后验中位数 d=%.4f   A=%.4f" % (stat["med_d"], stat["med_A"]))
print("MAP d=%.4f   A=%.4f" % (stat["map_d"], stat["map_A"]))
print("95%% 可信区间 d: [%.4f, %.4f]  半宽=%.5f μm" % (stat["lo_d"], stat["hi_d"], 0.5 * (stat["hi_d"] - stat["lo_d"])))
print("d-A 后验相关系数=%.4f   MCMC 接受率=%.3f   有效样本=%d" % (stat["corr"], acc, stat["n"]))

# 三方法点估计对比（复用范文一 / 二真源接口）
est_ls, mc_ls_mean, mc_ls_std, mc_ls_lo, mc_ls_hi, bias_ls, _ = G3.G.monte_carlo()
kappa, osc, n_arr, _ = G3.G2.detrend(Rmeas)
d_fd, _ = G3.G2.fd_invert(osc, n_arr, kappa, coarse=0.02)
print("三方法 d: 时域LS=%.4f  频域=%.4f  贝叶斯=%.4f  (真值 3.20)" % (x_ls[0], d_fd, stat["mean_d"]))

# 先验敏感性（3 组弱信息先验）
settings = [{"d_mu": 3.5, "d_sd": 1.0, "A_mu": 2.57, "A_sd": 0.50},
            {"d_mu": 3.5, "d_sd": 1.0, "A_mu": 2.57, "A_sd": 0.05},
            {"d_mu": 3.5, "d_sd": 1.0, "A_mu": 2.57, "A_sd": 0.02}]
sens = G3.prior_sensitivity(Rmeas, settings, init=list(x_ls))
print("先验敏感性(后验均值 d): 宽=%.4f  基准=%.4f  紧=%.4f" % tuple(sens))