单层膜反射谱与膜厚反演：从 FFT 初值到模型拟合

薄膜的反射谱里有干涉条纹，条纹的疏密编码了膜的光学厚度。从谱反推厚度叫膜厚反演（film thickness retrieval）。基本思路是：写出正向模型，调整厚度让模型谱贴近测量谱。本文用一个完全可复现的合成数据例子，把整个流程以及它的陷阱走一遍。

1. 正向模型#

考虑垂直入射、一层均匀无吸收的薄膜（厚度 $d$、折射率 $n_1$）在衬底（折射率 $n_2$）上，入射介质是空气（$n_0 = 1$）。界面的菲涅尔反射系数

$$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 d/\lambda$，总反射系数与反射率

$$r = \frac{r_{01} + r_{12}\,e^{\,2i\delta}}{1 + r_{01}\,r_{12}\,e^{\,2i\delta}},\qquad R(\lambda) = |r|^2 .$$

这里取 $n = n + ik$ 的约定（时间因子 $e^{-i\omega t}$），衬底可以有吸收（$k > 0$）。

import numpy as np

def reflectance(lam, d, n_film, n_sub):
    """Normal-incidence reflectance of one lossless film on a substrate in air.
    lam, d in the same length unit; n_film(lam), n_sub(lam) are callables."""
    n0, n1, n2 = 1.0, n_film(lam), n_sub(lam)
    r01 = (n0 - n1) / (n0 + n1)
    r12 = (n1 - n2) / (n1 + n2)
    e = np.exp(2j * 2 * np.pi * n1 * d / lam)          # exp(2 i delta)
    r = (r01 + r12 * e) / (1 + r01 * r12 * e)
    return np.abs(r) ** 2

下文的例子用的是玩具材料（只为演示，不要拿去当真实光学常数）：

• 膜：柯西色散 $n_1(\lambda) = 1.450 + 0.0035/\lambda^2$（$\lambda$ 单位 μm），类似氧化硅；
• 衬底：$n_2 = 3.6 + 0.08/\lambda^2 + 0.02\,i$，类似硅；
• 光谱：400–800 nm，512 个点；
• 真实厚度 $d_{\mathrm{true}} = 1.234\,\mu\mathrm{m}$，加 0.2% 的高斯噪声（绝对反射率）。

2. 目标函数与局部极小#

反演就是最小化残差：

$$\hat d = \arg\min_d \sum_i \big[\, R_{\mathrm{model}}(\lambda_i; d) - R_{\mathrm{meas}}(\lambda_i) \,\big]^2 .$$

麻烦在于 $R$ 随 $d$ 振荡：$d$ 变一点点，条纹就整体平移，$d$ 变 $\sim\bar\lambda/2n$（约 0.2 μm）时，条纹又"差不多对上了"，所以目标函数有很多局部极小。

图 b 里，真值处的 rms 残差是 0.2%（噪声水平），而邻近的局部极小（距真值约 ±0.19 μm）高达约 10%。实测从不同初值出发的收敛范围：初值偏离真值 ±0.09 μm 以内都能收敛到正确结果，更远就掉进邻近的坑。所以拟合前必须有一个靠谱的初值。

3. 用 FFT 给初值（厚膜）#

反射谱里的条纹分量 $\cos(4\pi n_1 d\,\sigma)$ 在 $\sigma = 1/\lambda$ 域是单一频率，正是 光谱域深度谱 里处理的对象：把 $R(\lambda)$ 重采样到均匀 $\sigma$、扣均值、加窗、做 FFT，峰在 $z_{\mathrm{app}} = n_g d$，所以

$$d_0 = \frac{z_{\mathrm{peak}}}{n_g},\qquad n_g = n - \lambda\frac{\mathrm{d}n}{\mathrm{d}\lambda} .$$

对柯西色散 $n = A + B/\lambda^2$，$n_g = A + 3B/\lambda^2$，本例 $n_g\approx 1.485$（550 nm）。FFT 给出 $d_0 = 1.2434\,\mu\mathrm{m}$，距真值只差 0.009 μm，远在 ±0.09 μm 的收敛范围内。用这个初值做最小二乘（scipy.optimize.least_squares），得到 $\hat d = 1.2340\,\mu\mathrm{m}$，残差 0.199%，等于噪声水平。

import numpy as np
from scipy.optimize import least_squares

def fft_guess(lam, R, n_g, z_min=0.25):
    sigma = 1.0 / lam
    s = np.linspace(sigma.min(), sigma.max(), len(R))
    Ru = np.interp(s, sigma[::-1], R[::-1])
    x = (Ru - Ru.mean()) * np.hanning(len(R))
    P = np.abs(np.fft.rfft(x, n=16 * len(R)))
    z = np.fft.rfftfreq(16 * len(R), d=s[1] - s[0]) / 2      # z_app = OPD / 2
    m = z > z_min                                            # skip the DC leakage
    return z[m][np.argmax(P[m])] / n_g

def fit_thickness(lam, R, d0, n_film, n_sub):
    f = lambda p: reflectance(lam, p[0], n_film, n_sub) - R
    return least_squares(f, [d0], x_scale=[0.01])

最后的精度可以从雅可比矩阵 $J$ 估计：$\mathrm{cov} = s^2\,(J^{\mathsf T}J)^{-1}$，其中 $s^2$ 是残差方差。

4. 薄膜：没有条纹，FFT 就失效了#

膜很薄（光学厚度小于一个条纹周期）时，整个光谱里只有半个周期不到的起伏，FFT 没有峰可找。例子里取 $d = 0.12\,\mu\mathrm{m}$：$z_{\mathrm{app}} = 0.18\,\mu\mathrm{m}$，低于这套光谱的分辨单元（400–800 nm 下 $\delta z\approx 0.4\,\mu\mathrm{m}$）。FFT 峰搜索只会停在搜索下限附近（0.25 μm），那不是信号。

此时只能靠模型拟合。图 c 是残差对厚度的扫描：只有一个明确的极小，落在 0.1198 μm，精修后得到 0.1200 μm。这里能成功，是因为模型假设（折射率、色散）完全正确。薄膜里厚度和折射率强烈相关，条件稍有偏差，结果就会漂（下一节）。

5. 精度的真相：噪声很小，系统误差很大#

图 d 给出蒙特卡洛结果（每个噪声水平 150 次）：

0.03 nm 的"精度"显然不可信，这只是统计精度：噪声被 512 个点平均掉了，而且模型与数据完全一致。真实测量的上限由系统误差决定：

• 折射率偏了，结果跟着偏。 数据用真实折射率生成，拟合时假设膜折射率高了 0.01（0.7%），得到 $\hat d = 1.2257\,\mu\mathrm{m}$，比真值小 8.3 nm（−0.67%），几乎正好等于 $n$ 的相对误差。因为条纹主要由乘积 $n\,d$ 决定。
• 残差几乎不报警。 上面这个错误的拟合，rms 残差从 0.199% 只升到 0.284%，肉眼看拟合曲线仍然完美贴合。
• 把 $n$ 也放开，$d$ 和 $n$ 几乎完全相关。 给 $n$ 加一个整体偏移 $\delta n$ 作为自由量，两者的相关系数是 $-0.997$，厚度的统计误差放大约 11.6 倍（0.2% 噪声下 0.032 nm → 0.37 nm）。

6. 实用建议#

1. 厚膜（出现多个条纹）：先 FFT 找初值，再最小二乘精修。初值需要在目标函数的"吸引盆"内（本例 ±0.09 μm，约 7% 的厚度），所以 $n_g$ 估到几个百分点就够。
2. 薄膜（条纹少于一个周期）：只能在一个宽范围内做网格扫描找全局极小，再精修；同时要承认厚度与折射率相关，只在 $n$ 可靠时才报厚度。
3. 看残差的结构，不只看它的大小。残差应该像白噪声；若有与波长相关的起伏，说明模型有问题（色散、衬底常数、入射角、粗糙度、膜不均匀）。
4. 衬底是硅这类强色散材料时，必须用实测的 $n(\lambda)$、$k(\lambda)$ 表，不能用常数。
5. 汇报结果时，把"统计误差"和"系统误差"分开写：前者来自噪声和拟合，后者来自折射率、波长标定和模型假设，通常后者大得多。