碳化硅外延层厚度的确定（优秀范文二：频域反演法）

本题为 2025 年全国大学生数学建模竞赛 B 题。本站范文逐题手写发布，本文为范文二，视角为频域（波数域）反演法，与范文一（时域非线性最小二乘拟合法）形成互补。全文约 2700 字、含 8 张 SVG 配图与可运行 Python 附录，正文、配图、附录、真源四路数字严格一致。

---

摘要

针对碳化硅（4H‑SiC）外延层厚度测定问题，本文提出一种基于波数域（频域）相位周期图的反演方法。在空气 / 外延层 / 衬底三层薄膜的干涉模型下，反射率 $R(\lambda)$ 对波长的振荡在波数 $\kappa=1/\lambda$ 域表现为近均匀的周期信号，其周期满足 $\Delta\kappa=1/(2n d)$。据此，本文先将实测反射率变换到波数域并去除慢变包络，得到纯振荡分量 $\mathrm{osc}(\kappa)$；随后在候选厚度网格上构建"相位周期图"——对每个候选厚度 $d$，以 $\arg=4\pi n(\kappa)d\kappa$ 为载波做线性最小二乘拟合，取振荡功率 $P(d)=a^2+b^2$ 的峰值作为厚度估计，并经抛物线插值与局部精扫细化。对一条含高斯噪声（标准差 $\sigma=0.006$）的仿真反射率曲线，本方法给出 $d=3.2001\ \mu\mathrm{m}$，与真值 $3.20\ \mu\mathrm{m}$ 的相对误差仅 $3\times10^{-5}$；蒙特卡洛 220 次重复得到均值 $3.2006\ \mu\mathrm{m}$、标准差 $0.0044\ \mu\mathrm{m}$。独立的"极值间距"交叉验证给出中位数估计 $3.2129\ \mu\mathrm{m}$，与频域结果及范文一时域拟合法（$3.1971\ \mu\mathrm{m}$）三者互相印证。本文同时量化了不确定度来源：测量噪声贡献约 $4.4\ \mathrm{nm}$，而忽略材料色散将引入约 $99.9\ \mathrm{nm}$ 的系统偏差——本方法通过逐点代入色散 $n(\kappa)$ 将其消除。

关键词：薄膜干涉；波数域；相位周期图；外延层厚度；蒙特卡洛不确定度

---

一、问题重述

已知空气（折射率 $N_0=1.0$）、碳化硅外延层（折射率随波长色散）、衬底 4H‑SiC（$N_2=2.70$）构成的三层薄膜结构，在波长 $400\sim1000\ \mathrm{nm}$ 范围内测得反射率 $R(\lambda)$（含测量噪声，标准差 $\sigma=0.006$）。试由该反射率曲线反演外延层厚度 $d$、并给出不确定度。外延层折射率服从 Cauchy 色散 $n_1(\lambda)=A+B/\lambda^2$，文献先验 $A=2.57,\ B=0.015$；真值厚度 $d=3.20\ \mu\mathrm{m}$。

范文一采用时域非线性最小二乘直接拟合 $R(\lambda)$，本文改从频域切入：利用干涉振荡在波数域的均匀周期性，将厚度反演转化为一个"寻峰"问题，从而在原理上规避了时域拟合对初值的依赖，并天然兼容色散。

---

二、模型假设

1. 薄膜为平行于表面的均匀平面层，入射光垂直入射，忽略吸收与散射；
2. 衬底折射率 $N_2=2.70$ 已知且无色散，空气折射率 $N_0=1.0$；
3. 外延层折射率服从 Cauchy 色散 $n_1(\lambda)=A+B/\lambda^2$，其中 $A, B$ 取文献先验，反演中视为已知；
4. 测量噪声为加性零均值高斯白噪声，标准差 $\sigma=0.006$，与波长无关；
5. 波长采样为 $400\sim1000\ \mathrm{nm}$ 等间隔 201 点（步长 $3\ \mathrm{nm}$）。

---

三、符号说明

符号	含义	符号	含义
$\lambda$	波长（nm）	$\kappa$	波数 $\kappa=1/\lambda$（$\mu\mathrm{m}^{-1}$）
$d$	外延层厚度（$\mu\mathrm{m}$）	$n_1(\lambda)$	外延层折射率（Cauchy）
$R(\lambda)$	反射率	$\mathrm{osc}(\kappa)$	去包络后的振荡分量
$\delta$	单程相位厚度	$P(d)$	相位周期图功率
$\sigma$	测量噪声标准差	$\Delta\kappa$	振荡在波数域的周期

---

四、模型建立

4.1 为何转向波数域

三层薄膜的双程干涉给出反射率（见范文一正演模型）：
$$R(\lambda)=|r|^2,\qquad \delta=\frac{2\pi n_1(\lambda)d}{\lambda}.$$
其中干涉项以 $\cos(2\delta)=\cos(4\pi n_1 d/\lambda)$ 振荡。在波长域，相邻反射率极值满足 $\Delta\lambda\approx \lambda^2/(2n d)$——即振荡间距随 $\lambda$ 增大而显著拉长，是非均匀的（图 1 上）。这给"数周期"式的波长域分析带来麻烦：长波端周期稀疏、短波端密集，离散化误差随波长变化。

若改用波数 $\kappa=1/\lambda$（单位 $\mu\mathrm{m}^{-1}$），则振荡相位变为
$$\arg = \frac{4\pi n_1(\kappa)d}{\lambda}=\frac{4\pi n_1(\kappa)d}{\kappa^{-1}}=4\pi n_1(\kappa)d\,\kappa,$$
在 $\kappa$ 域表现为近均匀的周期信号（图 1 下），周期为
$$\Delta\kappa=\frac{1}{2n\,d}\quad\Longrightarrow\quad d=\frac{1}{2n\,\Delta\kappa}.$$
这一步坐标变换，把"随波长变化的周期"变成了"近似恒定的频率"，从而可以用成熟的频域方法稳健提取厚度。

4.2 去包络与振荡分量提取

反射率曲线同时包含"慢变包络"（反射率整体水平随波长的缓变）与"快变干涉振荡"。为得到纯净的振荡分量，本文对 $\lambda$ 做三阶多项式拟合
$$R(\lambda)\approx c_0+c_1\lambda+c_2\lambda^2+c_3\lambda^3,$$
以其作为包络估计并相减：
$$\mathrm{osc}(\lambda)=R(\lambda)-\mathrm{trend}(\lambda).$$
随后按 $\kappa=1/\lambda$ 重新排序，得到波数域振荡序列 $\mathrm{osc}(\kappa)$。图 2 显示，$\mathrm{osc}(\kappa)$ 是一条振幅基本恒定、周期近均匀的准正弦波，正是频域分析的理想输入。

4.3 相位周期图与厚度反演

对某个候选厚度 $d$，构造载波 $\arg_i=4\pi n_i d\kappa_i$（其中 $n_i=n_1(\lambda_i)$ 为逐点色散值），用线性最小二乘拟合
$$\mathrm{osc}_i\approx a\cos(\arg_i)+b\sin(\arg_i),$$
得到系数 $(a,b)$。定义振荡功率
$$P(d)=a^2+b^2,$$
它衡量"以厚度 $d$ 为载波频率时，数据能被正弦成分解释的程度"。当 $d$ 接近真值时，拟合残差最小、$P(d)$ 最大。因此厚度估计即周期图功率曲线的峰值位置（图 3）。

为兼顾效率与精度，采用"粗扫 + 细化"两步：

1. 粗扫：在 $d\in[0.5,8.0]\ \mu\mathrm{m}$ 以步长 $0.02$ 计算 $P(d)$，定位峰值 $d^*$；
2. 细化：在 $d^*$ 邻域以步长 $0.001$ 作局部精扫，并用相邻三点做抛物线插值，最终得到 $d=3.2001\ \mu\mathrm{m}$（图 4）。峰位与真值 $3.20\ \mu\mathrm{m}$ 几乎重合，验证了频域反演的准确性。

4.4 极值间距交叉验证

作为与周期图互相独立的检验，本文另用"极值间距法"：在 $\mathrm{osc}(\kappa)$ 中检测局部极大值，对每对相邻极大值测其波数间距 $\Delta\kappa$，由 $d=1/(2\bar n\,\Delta\kappa)$ 反推一个厚度估计。为抑制噪声诱发的伪极大值，先做窗口 5 轻平滑，并以最小间隔（4 个采样点）贪心筛除过近的次级峰。共得到 21 对独立估计，其中位数 $3.2129\ \mu\mathrm{m}$、均值 $3.1234\ \mu\mathrm{m}$，散布区间 $1.54\sim4.86\ \mu\mathrm{m}$（图 5）。中位数与时域、频域结果高度一致，说明三种方法在系统层面互相印证；而间距法的较大散布也提示它作为单点法对噪声更敏感——这正是周期图"整体拟合"更稳健的原因。

---

五、模型求解与结果分析

5.1 频域反演流程

综合上述步骤，频域反演流程为：① 读取反射率 $R(\lambda)$；② 三阶多项式去包络；③ 变换至波数域得 $\mathrm{osc}(\kappa)$；④ 相位周期图粗扫定位峰值；⑤ 抛物线 + 局部精扫细化得最终 $d$；⑥ 极值间距法独立交叉验证。该流程无需时域拟合那样的多参数初值猜测，对初值完全不敏感。

5.2 反演结果与多方法对比

将频域法与时域法（范文一的非线性最小二乘）、包络法（波长域极值间距初值）对同一条含噪曲线反演，结果如下（图 6）：

方法	厚度估计（$\mu\mathrm{m}$）	说明
包络法	3.8268	波长域初值，受噪声上漂明显
时域拟合（范文一）	3.1971	非线性最小二乘，需初值
频域相位（本文）	3.2001	波数域寻峰，免初值、兼容色散
真值	3.20	—

频域法估计 $3.2001\ \mu\mathrm{m}$ 与真值仅差 $0.0001\ \mu\mathrm{m}$；包络法因波长域噪声导致极值定位上漂（偏差约 $0.63\ \mu\mathrm{m}$），凸显了频域变换在稳健性上的优势。

5.3 残差与稳健性

频域法本质是对 $\mathrm{osc}(\kappa)$ 的整体正弦拟合，残差主要来自测量噪声而非模型失配；由于它利用了全部 201 个采样点的相位信息，对个别异常点的鲁棒性优于仅依赖少数极值点的间距法。后文蒙特卡洛实验进一步量化了这一稳健性。

---

六、灵敏度分析

6.1 蒙特卡洛不确定度

固定真值 $d=3.20\ \mu\mathrm{m}$ 与噪声水平 $\sigma=0.006$，重复 220 次"加噪反射率 → 频域反演"实验。得到的厚度估计分布（图 7）均值为 $3.2006\ \mu\mathrm{m}$、标准差 $0.0044\ \mu\mathrm{m}$，95% 置信区间 $[3.1914,\ 3.2092]\ \mu\mathrm{m}$，偏差（均值减真值）仅 $0.0006\ \mu\mathrm{m}$。分布对称、无系统偏移，说明频域反演具有良好的一致性与精度。

6.2 噪声与色散的影响

不确定度预算分解（图 8）显示：① 测量噪声（$\sigma=0.006$）是随机不确定度的主要来源，贡献约 $4.4\ \mathrm{nm}$；② 网格离散经细化后仅余约 $0.5\ \mathrm{nm}$，可忽略；③ 若忽略色散、用常折射率近似，将引入约 $99.9\ \mathrm{nm}$ 的系统偏差——本文在周期图载波中逐点代入 $n(\kappa)$，将该系统误差消除至可忽略水平。可见，材料色散是比采样网格更值得重视的误差源，必须在模型中显式处理。

---

七、模型评价

优点：① 原理清晰，将厚度反演转化为频域寻峰，物理直观；② 对初值完全不敏感，仅需单次扫描即可收敛，工程实现简单；③ 天然兼容材料色散（逐点 $n(\kappa)$），规避了常折射率近似的百纳米级系统偏差；④ 与极值间距法互为独立交叉验证，结果互相印证；⑤ 对测量噪声稳健，蒙特卡洛标准差仅 $4.4\ \mathrm{nm}$。

局限与改进：① 周期图分辨率受波长范围与采样点数制约，更厚样品（更多周期）精度更高、更薄样品（周期数少）分辨率下降；② 当前仅反演厚度，若 $A,B$ 亦未知，可扩展为 $(d,A,B)$ 二维/三维周期图或频域–时域混合反演；③ 极值间距法对噪声较敏感，可改用 Hilbert 包络或加窗谱估计进一步提升交叉验证的精度。

---

八、结论

本文将薄膜干涉反射率分析从波长域拓展到波数域，建立了基于相位周期图的碳化硅外延层厚度频域反演方法。对含噪仿真数据的反演结果为 $d=3.2001\ \mu\mathrm{m}$，相对真值误差约 $3\times10^{-5}$；蒙特卡洛 220 次给出不确定度 $0.0044\ \mu\mathrm{m}$。该方法无需初值猜测、兼容色散、稳健性强，与范文一的时域拟合法、波长域包络法三者结论一致，为外延层厚度的光学测定提供了一条互补且可靠的途径。

---

九、参考文献

1. 2025 高教社杯全国大学生数学建模竞赛 B 题《碳化硅外延层厚度的确定》赛题.
2. Born M., Wolf E. Principles of Optics. Cambridge University Press, 1999.（薄膜干涉与特征矩阵）
3. 徐叙瑢等. 薄膜光学与镀膜技术. 科学出版社, 2018.（包络法与色散模型）
4. Press W. H., et al. Numerical Recipes: The Art of Scientific Computing. Cambridge, 2007.（周期图与最小二乘）
5. 全国大学生数学建模竞赛组委会. 数学建模方法与分析. 高等教育出版社.

---

附录：核心 Python 实现

以下代码复用本题真源模块 gen_cumcm2025b_2.py（前向模型与实测数据复用 gen_cumcm2025b.py），可直接运行并复现正文全部关键数字。

import os, sys, math
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_2 as G2          # 频域反演真源
import gen_cumcm2025b_2 as G2mod

# 1) 读取同一条含噪仿真反射率（与范文一、正文同源）
Rmeas, clean = G2.G.build_meas()

# 2) 变换到波数域并去包络
kappa, osc, n_arr, trend = G2.detrend(Rmeas)

# 3) 相位周期图反演厚度（粗扫 + 抛物线 + 局部精扫）
d_fd, P_fd = G2.fd_invert(osc, n_arr, kappa, coarse=0.02)
print("频域反演 d = %.4f um  (振荡功率 P = %.4f)" % (d_fd, P_fd))

# 4) 极值间距法独立交叉验证
spacing = G2.extreme_spacing_d(osc, kappa, n_arr)
sp_med = sorted(spacing)[len(spacing) // 2]
print("极值间距法: 对数=%d  median=%.4f um  mean=%.4f um"
      % (len(spacing), sp_med, sum(spacing) / len(spacing)))

# 5) 与范文一时域拟合法、波长域包络法对比
x_td, f_td, it = G2.G.fit(Rmeas, G2.LAM, seed=2026)
env = G2.G.envelope_d(Rmeas)
print("三方法对比: 包络法=%.4f  时域拟合=%.4f  频域相位=%.4f"
      % (env, x_td[0], d_fd))

# 6) 蒙特卡洛不确定度（220 次加噪重复反演）
est, mean, std, lo, hi, bias = G2.fd_monte_carlo()
print("MC(220): mean=%.4f std=%.5f CI=[%.4f, %.4f] bias=%.5f"
      % (mean, std, lo, hi, bias))