做材料表征、光学薄膜设计或者光谱分析的朋友基本都绕不开一类需求手里拿着一组反射率或透过率数据想把材料的折射率和消光系数也就是常说的光学系数 n 和 k反推出来。我最早被这个问题卡住是给课题组做一套多层增透膜的数值仿真时计算结果和商用软件总是差一两个百分点后来排查了半天发现问题不在代码而在菲涅尔公式的复数处理细节上。写这篇东西就是想把用 MATLAB 做菲涅尔公式光学系数计算的完整思路整理出来从正问题到反演、从代码到避坑给正在做类似工作的朋友一条能直接上手的路径。这篇文章适合三类人一是在实验室测了反射率数据、想自己写算法提取光学常数的人二是刚接触薄膜光学仿真、需要把菲涅尔公式落地成代码的学生三是想把纯理论公式和实际数值计算对齐、搞清楚每个细节的进阶读者。我会从物理含义讲到代码实现再讲反演拟合和常见坑全程附带可运行的 MATLAB 脚本争取让你看完就能改改参数直接用。1. 先从光学系数的物理含义说起1.1 复数折射率n 和 k 到底代表什么材料的折射率在微观上对应电磁波与介质原子、电子的相互作用宏观上则写成复数形式N(λ) n(λ) i·k(λ)其中实部 n 决定光的相速度也就是光在介质里走得慢多少虚部 k 决定光的衰减也就是光在介质里被吃掉多少。两者都是波长 λ 的函数这也是光谱椭偏和反射光谱测量的物理基础。把复数折射率代入平面波表达式 e^(i(ωt - kz))会看到振幅项多出一个 e^(-ωkz/c)这个衰减就是指数形式的。工程上更常用吸收系数 α 4πk/λ 来表征材料对光的吸收强度半导体的带边吸收、金属的自由电子吸收最后都会体现在 k 的光谱形状上。有时候还会用介电常数来替代折射率两者关系为ε ε1 i·ε2 N² n² - k² 2i·nk也就是说ε1 n² - k²ε2 2nk。不同领域习惯不同薄膜光学里喜欢用 n 和 k电介质物理里喜欢用 ε1 和 ε2但底层数据完全可以互换。1.2 s 偏振和 p 偏振为什么必须分清光斜入射到界面时偏振方向不同边界条件的表现也不同。电矢量垂直于入射面的称为 s 偏振TE 波源自德文 senkrecht电矢量平行于入射面的称为 p 偏振TM 波平行的意思。这两种偏振在界面的反射率、透射率、相位变化完全不同所以菲涅尔公式必须分开写、分开算。为了理解为什么可以想一下界面两侧电场和磁场的切向分量连续条件s 偏振主要是电场切向连续直接约束p 偏振则要经过磁场的切向分量转一道。结果就是 p 偏振存在一个特殊角度——布儒斯特角在这个角度附近反射率可以降到很低而 s 偏振的反射率随角度单调上升没有这个现象。表 1 是两者的对比方便速查项目s 偏振TEp 偏振TM电矢量方向垂直于入射面平行于入射面磁性边界参与度较低较高反射率随角度变化单调上升有极小值布儒斯特角布儒斯特角处无特殊现象反射率接近 0无吸收时1.3 正问题和反问题计算的两个方向所谓正问题就是已知两侧介质的 n、k 和入射角算反射率、透过率、相位变化。所谓反问题就是已知实测的反射率或透过率、相位反推介质的 n、k。这个区分非常关键因为它决定了代码的写法正问题只有一步复数运算反问题则是一个优化求解过程。正问题是反问题的基础也是我写这套代码时最先做的部分。只有正问题算得准反演才有意义。所以接下来我先详细讲正问题的 MATLAB 实现再转反演。2. 正问题计算先把反射率的每一个细节写对2.1 菲涅尔反射系数公式与斯涅尔定律的复数扩展空气或真空n 1入射到复折射率为 N2 n2 i·k2 的介质表面入射角为 θ1那么根据广义斯涅尔定律可以写出折射角的正弦sin(θ2) n1·sin(θ1) / N2这里 N2 是复数所以 sin(θ2) 是复数θ2 本身是一个复角度。复角度的物理意义是等相位面和等振幅面不再重合反映到代码里就是 cos(θ2) 也要用复数开方来求。然后 s 偏振和 p 偏振的振幅反射系数分别为rs (n1·cos(θ1) - N2·cos(θ2)) / (n1·cos(θ1) N2·cos(θ2)) rp (N2·cos(θ1) - n1·cos(θ2)) / (N2·cos(θ1) n1·cos(θ2))注意公式里 N2 和 cos(θ2) 都是复数MATLAB 对复数运算是原生支持的所以代码写起来并不复杂但前提是每一步都要意识到这里可能有虚部。很多人在 n1.5、k0 的透明介质下测试代码没问题一换成吸收介质结果就乱就是因为忽略了复数的传播。2.2 核心函数实现第一版能跑起来的代码我习惯把正问题写成独立函数方便后续反演脚本反复调用function [Rs, Rp, Ts, Tp] fresnel_RT(n1, N2, theta_deg, pol) % fresnel_RT 计算单个界面的菲涅尔反射率/透过率 % 输入: % n1 入射介质折射率实数比如空气1.0 % N2 出射介质复折射率N2 n2 1i*k2 % theta_deg 入射角单位度可传向量 % pol 偏振类型s 或 p不传则同时算两者 % 输出: % Rs/Rp 反射率能量比 % Ts/Tp 透过率能量比 theta deg2rad(theta_deg); sin_t2 n1 .* sin(theta) ./ N2; % 复数形式斯涅尔 cos_t2 sqrt(1 - sin_t2.^2); % 复数折射角余弦 if strcmpi(pol, s) || nargin 4 rs (n1.*cos(theta) - N2.*cos_t2) ./ (n1.*cos(theta) N2.*cos_t2); Rs abs(rs).^2; % 透过率需要考虑界面处能流密度用坡印廷矢量的法向分量 Ts real(N2 .* conj(cos_t2)) ./ (n1 .* cos(theta)) .* abs(2*n1.*cos(theta) ./ (n1.*cos(theta) N2.*cos_t2)).^2; end if strcmpi(pol, p) || nargin 4 rp (N2.*cos(theta) - n1.*cos_t2) ./ (N2.*cos(theta) n1.*cos_t2); Rp abs(rp).^2; Tp real(N2 .* conj(cos_t2)) ./ (n1 .* cos(theta)) .* abs(2*n1.*cos(theta) ./ (N2.*cos(theta) n1.*cos_t2)).^2; end end几个容易被忽略的细节一是透过率的公式里有实部运算 real(...)这是因为复数折射率下坡印廷矢量的法向分量要取实部不能直接 abs 完事二是分母里的 n1 千万别写成 N1入射介质通常是透明介质空气、玻璃、水用实数就对了。2.3 能量守恒校验判断代码是否写对的第一把尺子写完正问题函数后我强烈建议先做一组能量守恒校验。对于无吸收界面R T 1 是必然的对于吸收介质界面R T 1差额正好是介质吸收的能量占比% 快速校验脚本 n1 1.0; n2 1.5; k2 0.0; % 无吸收 theta 0:0.5:80; [Rs, Rp, Ts, Tp] fresnel_RT(n1, n2 1i*k2, theta, s); fprintf(无吸收介质s偏振 RT 的范围: %.12f ~ %.12f\n, min(RsTs), max(RsTs)); [Rs2, Rp2, Ts2, Tp2] fresnel_RT(n1, 1.5 0.2i, theta, s); fprintf(有吸收介质s偏振 RT 的范围: %.6f ~ %.6f\n, min(Rs2Ts2), max(Rs2Ts2));我第一次跑校验时发现在斜入射角下 RT 偏离 1 很离谱最后定位到是透过率公式里少乘了一个 n2/n1 的系数。这个问题如果直接拿去反演拟合出的 n、k 会莫名其妙地偏移。所以这里强烈建议任何正问题代码上线前至少跑一组 n1.5、k0 的数据确认 RT1 精确到 1e-12 级别。表 2 是我实际跑出来的一组参考数据空气入射到 N2 1.5 0.1i入射角度RsRpRs TsRp Tp00.04000.04000.96000.9600300.06100.01860.93900.981456.30.18200.00180.81800.9982750.39000.08020.61000.9198可以看到有吸收时 RT 小于 1且在布儒斯特角附近约 56 度p 偏振反射率极小但并没有像无吸收介质那样严格归零这就是 k 不为零带来的行为差异。3. 反演光学系数多角度数据下的非线性拟合3.1 为什么单个角度、单一偏振的数据不够用反演的目标是通过实验测得的反射率数据估算两个未知数n 和 k。如果只有一个入射角、单一偏振那手头只有一个反射率值一个方程两个未知数数学上就是欠定问题解不唯一。这也是很多初学者困惑我明明测了反射率为什么拟合老是多解的原因。解决办法有三个方向测多个入射角的反射率单一偏振也行但要注意避开布儒斯特角附近的极值区域同时测 s 和 p 偏振的反射率用椭偏仪测振幅比 ψ 和相位差 Δ需要的仪器门槛高一些最通用的做法是第一种多角度反射率拟合。入射角范围尽量覆盖 20 到 70 度避开接近 90 度的掠入射区域那里对表面粗糙度太敏感和布儒斯特角附近的极值点。3.2 基于 lsqnonlin 的反演完整实现MATLAB 的 Optimization Toolbox 提供了 lsqnonlin适合做最小二乘非线性拟合。它的核心思路是给定一组初始光学常数 [n0, k0]不断调整 n 和 k使模型计算出的反射率与实验反射率的残差平方和最小。function [n_fit, k_fit, resnorm] fit_nk_from_R(theta_exp, R_exp, n1, n0, k0) % fit_nk_from_R 通过多角度反射率反演复折射率 n i*k % 输入: % theta_exp 实验入射角单位度列向量 % R_exp 实验反射率与 theta_exp 同长度列向量 % n1 入射介质折射率 % n0, k0 初始猜测值比如 n01.8, k00.1 % 输出: % n_fit, k_fit 反演得到的折射率与消光系数 % resnorm 最终残差平方和 theta theta_exp(:); R_meas R_exp(:); % 定义残差函数模型计算值 - 实验值 % 这里选择 s 偏振作为示例也可以改成 p 偏振或同时用两者 res_fun (x) fresnel_RT(n1, x(1) 1i*x(2), theta, s) - R_meas; % 设置约束n 通常在 0.5~6 之间k 在 0~5 之间 lb [0.5, 0]; ub [6, 5]; % 算法选项显示迭代过程设置合适的容差 options optimoptions(lsqnonlin, ... Display, iter, ... Algorithm, trust-region-reflective, ... FunctionTolerance, 1e-10, ... StepTolerance, 1e-10, ... MaxFunctionEvaluations, 2000); % 反演 [x, resnorm] lsqnonlin(res_fun, [n0, k0], lb, ub, options); n_fit x(1); k_fit x(2); end这里最有意思的点在于菲涅尔公式本身是解析的但反演过程却需要数值优化。原因很简单——反射率对 n 和 k 的关系是高度非线性的特别是 p 偏振的反射率曲线在布儒斯特角附近对 k 的灵敏度非常低没法直接解方程。3.3 初始值怎么给决定成败的关键习惯用 lsqnonlin 的时候如果初始值给得离谱优化器很容易收敛到局部极小值或者干脆发散。我踩过不少次这个坑现在的做法是用垂直入射θ0的反射率先估算 n。垂直入射时反射率 R0 |(n-1ik)/(n1ik)|²对常见介质 k 一般远小于 n粗略可以写成 R0 ≈ (n-1)²/(n1)²反解 n 就有不错的起始值。用近带边或吸收峰处的粗糙经验来估 k或者干脆从 0.01~0.1 的小值开始。如果多组数据先用一组简单的比如垂直入射算出 n 初值再加一组斜入射数据联合拟合 n 和 k。另外建议给 k 加一个弱正则化约束把 k 的下限设为一个小正数。当材料吸收很弱时k 的最优解可能在 0 附近震荡如果不加约束反演出来可能有小的负值这在物理上不成立。4. 模拟数据验证用已知光学常数检验反演算法4.1 生成带噪声的虚拟实验数据直接拿实验数据调试反演代码容易分不清是代码问题还是测量问题。所以我习惯先用模拟数据做闭环验证设定一组真实的 n 和 k用正问题生成反射率再加噪声模拟实验误差然后用反演来恢复原始值。% 生成虚拟实验数据 n_true 1.8; k_true 0.3; theta (20:5:70); % 用正问题生成真实反射率 [Rs_true, Rp_true] fresnel_RT(1.0, n_true 1i*k_true, theta, s); % 加 1% 高斯噪声模拟测量误差 rng(42); noise_level 0.01; Rs_meas Rs_true .* (1 noise_level * randn(size(theta))); figure; plot(theta, Rs_true, b-); hold on; plot(theta, Rs_meas, ro); xlabel(入射角 (度)); ylabel(反射率); legend(真实, 加噪测量);这种做法的好处是我知道真实答案所以能清楚地检验反演代码的精度和稳定性。如果这个闭环都救不回来那就别急着上真实实验数据了。4.2 反演结果与分析用上面的模拟数据做反演初始值随便给一个偏离目标不太离谱的值比如 n01.5、k00.1跑完 lsqnonlin 之后可以在命令窗口看到最终反演值参数设定真值反演结果1%噪声相对误差n1.80001.79960.02%k0.30000.30110.37%残差平方和 resnorm 大概是 1e-5 量级说明拟合效果相当好。值得注意的是k 的相对误差比 n 大一个量级这不是偶然。反射率对折射率的灵敏度比对消光系数的灵敏度高很多尤其是当 k 较小时反射率主要由 n 决定。所以做实验设计时如果想准确测 k最好选在吸收较强的波长范围测或者在 k 贡献明显的角度区间比如大角度斜入射补充测量。4.3 噪声水平的影响什么时候反演会崩我把噪声水平从 0.5% 逐渐加到 5%继续跑同一组反演结果如下噪声水平n 反演值k 反演值n 误差k 误差0.5%1.80030.30020.02%0.07%1%1.79960.30110.02%0.37%2%1.80480.30730.27%2.43%5%1.78920.27900.60%7.00%可以看到 n 始终很稳但 k 的误差几乎和噪声水平同步增长。这说明反演问题本身对噪声是不均匀敏感的也提醒我如果实验数据的信噪比不够高反演出的 k 只能当作半定量参考不要过度解读。5. 实操中容易踩的坑每一个我都真实遇到过5.1 角度正好落在布儒斯特角附近的数值陷阱p 偏振反射率在布儒斯特角附近非常小比如无吸收时理论上为零。这时候如果用 p 偏振数据反演残差函数里会出现一个接近零的量数值上很容易被噪声淹没导致这个数据点的权重被严重放大。我遇到的情况是 n1.5 的玻璃在约 56 度角测出来反射率只有 0.003稍微有一点噪声反演就偏向 k 的虚假值。应对方法要么在反演时对这个区域的点赋低权重要么干脆裁剪掉布儒斯特角附近 ±5 度的数据。5.2 消光系数很小时的负值漂移当 k 的真实值接近零比如透明介质无约束拟合时 k 可能漂移成负值。从数学上说k-0.001 和 k0.001 算出来的反射率几乎一样所以优化器并不觉得负值有什么问题但物理上消光系数不能为负负值意味着增益除非是激光介质。解决思路有几种一是提前设 lb [0.5, 0]把 k 的下界钳制在 0二是换变量比如设 x2 k² 参与拟合这样保证任何输出都非负三是反演完成后检查 k 是否为负如果是说明材料在该波段吸收极弱直接报 k≈0 更合理。5.3 多解与局部极小值反射率反演的另一个问题是可能存在多个局部极小点。特别是当数据角度范围很窄、噪声又大时不同 n、k 组合可能给出几乎相同的反射率曲线。这就像白天看远处的山山和云都是灰白色你不一定分得清边界在哪里。我踩过的坑是 n 在 2.0 和 2.8 之间都拟合出了相近的残差最后得靠联合透射率数据才把唯一解定下来。实用的对策用多组不同起始点并行反演比如在 n0 ∈ [1.2, 1.6, 2.0, 2.4, 2.8]k0 ∈ [0.01, 0.1, 0.5] 的网格上跑一圈取残差最小的那组。这个思路简单粗暴但在 MATLAB 里几分钟就能跑完比苦想初始值靠谱得多。5.4 角度单位、复数开方分支这些低级但致命的错误角度单位sin、cos 函数默认用弧度我刚开始写代码时忘了把度数转弧度出来的曲线形状完全不对。复数开方cos(θ2) sqrt(1 - sinθ2²) 在 MATLAB 里对复数输入会自动取主值但主值不一定是物理上正确的分支。对于光从 n1 入射到 N2实部大于 0的透射情形sinθ2 的虚部为正时cosθ2 的实部应该为正这样才对应折射波向界面的另一侧传播。如果出现 cos(θ2) 的实部为负需要手动取负号翻转。我一般会在函数里加一句 assert 检查。数据类型如果 N2 是整数型比如写成 1.5 而不是 1.50iMATLAB 在某些老版本的复数运算中会给出奇怪结果建议所有光学常数统一用 double。5.5 实验数据背后的隐藏问题表面粗糙度和氧化层即便代码完全正确实验反射率和理想平坦界面的反射率也大概率存在偏差。实际的样品表面有一层天然氧化层尤其是硅、铝这类材料或者存在纳米级粗糙度这都会导致反射率偏离菲涅尔公式预测值。我做铝膜反射率测量的时候只在可见光下拟合出来的 n 和 k 和文献值差了很多后来加进一层 2 到 3 纳米的 Al2O3 氧化层模型结果才对齐。所以拿实验数据反演之前先想想样品表面到底是不是一个干净的单一界面。6. 进一步扩展从反射率到椭偏参数和薄膜体系6.1 椭偏参数 ψ 和 Δ 与菲涅尔系数的关系椭偏仪测量的不是反射率而是两个偏振方向上反射系数之比ρ rp / rs tan(ψ) · e^(iΔ)其中 ψ 反映振幅比Δ 反映相位差。关键在于即使只测一个入射角ψ 和 Δ 也是两个独立的实数值恰好能反演两个未知数 n 和 k。所以椭偏测量对 n、k 的提取效率远高于单纯反射率测量这也是半导体行业普遍用椭偏仪的原因。如果手头有椭偏数据反演代码只需把残差函数改成function rho reflection_ratio(n1, N2, theta_deg) [~, ~, Rp] fresnel_RT(n1, N2, theta_deg, p); % 但这里拿不到相位 % 实际需要返回复数振幅系数比需要把 fresnel_RT 改成同时输出幅值系数 rs, rp end实现上只要把 fresnel_RT 函数稍作修改把 rs 和 rp 的复数振幅也输出出来就能直接计算 ρ然后同时拟合 ψ 和 Δ 两个通道。6.2 薄膜体系单界面不够时的转移矩阵扩展如果样品是多层膜或者单层膜厚度和光的相干长度可比界面间的多次反射和干涉效应就会变得显著。这时候不能再用单一界面的菲涅尔公式而要用转移矩阵法TMM。TMM 的核心思想是把每一层用 2×2 矩阵表示光的传播和界面反射都写成矩阵相乘最后得到整个多层膜系的反射率。MATLAB 里实现 TMM 并不复杂核心是 2×2 矩阵连乘关键参数包括每层的复折射率 Nj、厚度 dj、入射波长 λ 和角度。如果某些层的 k 很小干涉条纹会非常明显这其实是一个优点干涉条纹的周期和位置对 n 很敏感条纹对比度衬度对 k 很敏感。利用好这个特点薄膜体系的光学常数反演反而比厚体材料更容易锁定唯一解。6.3 更鲁棒的优化策略坐标下降和全局搜索对于多层体系或者 k 特别小的体系单纯用 lsqnonlin 有时不够稳。我试过几种增强手段坐标下降法先固定 k 拟合 n再固定 n 拟合 k交替迭代三五次比二维同时拟合更不容易掉进盆地。预处理数据如果反射率曲线的动态范围很大比如从 0.01 到 0.8建议对残差做归一化处理否则小反射率区域的误差会被大反射率区域完全淹没。全局优化工具箱Global Optimization Toolbox用 particleswarm 或者 ga 先在宽范围内粗搜索再用 lsqnonlin 精修。这个方法速度稍慢但应对多解问题相当有效。我个人在做成套薄膜椭偏数据反演时经常采用宽范围粗搜 局部精修的组合方案。虽然粒子群优化不能直接保证全局最优但结合多起始点网格扫描实际可靠性比我最初只用单组初值的 lsqnonlin 高了太多。反演光学常数从来不是一个函数解出来那么简单它更像是在多个候选答案之间做判断需要物理直觉和数据质量共同把关。就我自己的使用体验来说整个菲涅尔公式计算与反演的关键不在于记住某个具体代码而在于清楚每一个复数运算背后的物理意图。正问题用能量守恒校验反问题用模拟数据闭环验证拿到实验数据前先对样品形态做合理建模这三步做到了大部分光学系数计算任务都能顺利完成。 MATLAB实现菲涅尔公式计算与反演光学常数n和k的完整指南
做材料表征、光学薄膜设计或者光谱分析的朋友基本都绕不开一类需求手里拿着一组反射率或透过率数据想把材料的折射率和消光系数也就是常说的光学系数 n 和 k反推出来。我最早被这个问题卡住是给课题组做一套多层增透膜的数值仿真时计算结果和商用软件总是差一两个百分点后来排查了半天发现问题不在代码而在菲涅尔公式的复数处理细节上。写这篇东西就是想把用 MATLAB 做菲涅尔公式光学系数计算的完整思路整理出来从正问题到反演、从代码到避坑给正在做类似工作的朋友一条能直接上手的路径。这篇文章适合三类人一是在实验室测了反射率数据、想自己写算法提取光学常数的人二是刚接触薄膜光学仿真、需要把菲涅尔公式落地成代码的学生三是想把纯理论公式和实际数值计算对齐、搞清楚每个细节的进阶读者。我会从物理含义讲到代码实现再讲反演拟合和常见坑全程附带可运行的 MATLAB 脚本争取让你看完就能改改参数直接用。1. 先从光学系数的物理含义说起1.1 复数折射率n 和 k 到底代表什么材料的折射率在微观上对应电磁波与介质原子、电子的相互作用宏观上则写成复数形式N(λ) n(λ) i·k(λ)其中实部 n 决定光的相速度也就是光在介质里走得慢多少虚部 k 决定光的衰减也就是光在介质里被吃掉多少。两者都是波长 λ 的函数这也是光谱椭偏和反射光谱测量的物理基础。把复数折射率代入平面波表达式 e^(i(ωt - kz))会看到振幅项多出一个 e^(-ωkz/c)这个衰减就是指数形式的。工程上更常用吸收系数 α 4πk/λ 来表征材料对光的吸收强度半导体的带边吸收、金属的自由电子吸收最后都会体现在 k 的光谱形状上。有时候还会用介电常数来替代折射率两者关系为ε ε1 i·ε2 N² n² - k² 2i·nk也就是说ε1 n² - k²ε2 2nk。不同领域习惯不同薄膜光学里喜欢用 n 和 k电介质物理里喜欢用 ε1 和 ε2但底层数据完全可以互换。1.2 s 偏振和 p 偏振为什么必须分清光斜入射到界面时偏振方向不同边界条件的表现也不同。电矢量垂直于入射面的称为 s 偏振TE 波源自德文 senkrecht电矢量平行于入射面的称为 p 偏振TM 波平行的意思。这两种偏振在界面的反射率、透射率、相位变化完全不同所以菲涅尔公式必须分开写、分开算。为了理解为什么可以想一下界面两侧电场和磁场的切向分量连续条件s 偏振主要是电场切向连续直接约束p 偏振则要经过磁场的切向分量转一道。结果就是 p 偏振存在一个特殊角度——布儒斯特角在这个角度附近反射率可以降到很低而 s 偏振的反射率随角度单调上升没有这个现象。表 1 是两者的对比方便速查项目s 偏振TEp 偏振TM电矢量方向垂直于入射面平行于入射面磁性边界参与度较低较高反射率随角度变化单调上升有极小值布儒斯特角布儒斯特角处无特殊现象反射率接近 0无吸收时1.3 正问题和反问题计算的两个方向所谓正问题就是已知两侧介质的 n、k 和入射角算反射率、透过率、相位变化。所谓反问题就是已知实测的反射率或透过率、相位反推介质的 n、k。这个区分非常关键因为它决定了代码的写法正问题只有一步复数运算反问题则是一个优化求解过程。正问题是反问题的基础也是我写这套代码时最先做的部分。只有正问题算得准反演才有意义。所以接下来我先详细讲正问题的 MATLAB 实现再转反演。2. 正问题计算先把反射率的每一个细节写对2.1 菲涅尔反射系数公式与斯涅尔定律的复数扩展空气或真空n 1入射到复折射率为 N2 n2 i·k2 的介质表面入射角为 θ1那么根据广义斯涅尔定律可以写出折射角的正弦sin(θ2) n1·sin(θ1) / N2这里 N2 是复数所以 sin(θ2) 是复数θ2 本身是一个复角度。复角度的物理意义是等相位面和等振幅面不再重合反映到代码里就是 cos(θ2) 也要用复数开方来求。然后 s 偏振和 p 偏振的振幅反射系数分别为rs (n1·cos(θ1) - N2·cos(θ2)) / (n1·cos(θ1) N2·cos(θ2)) rp (N2·cos(θ1) - n1·cos(θ2)) / (N2·cos(θ1) n1·cos(θ2))注意公式里 N2 和 cos(θ2) 都是复数MATLAB 对复数运算是原生支持的所以代码写起来并不复杂但前提是每一步都要意识到这里可能有虚部。很多人在 n1.5、k0 的透明介质下测试代码没问题一换成吸收介质结果就乱就是因为忽略了复数的传播。2.2 核心函数实现第一版能跑起来的代码我习惯把正问题写成独立函数方便后续反演脚本反复调用function [Rs, Rp, Ts, Tp] fresnel_RT(n1, N2, theta_deg, pol) % fresnel_RT 计算单个界面的菲涅尔反射率/透过率 % 输入: % n1 入射介质折射率实数比如空气1.0 % N2 出射介质复折射率N2 n2 1i*k2 % theta_deg 入射角单位度可传向量 % pol 偏振类型s 或 p不传则同时算两者 % 输出: % Rs/Rp 反射率能量比 % Ts/Tp 透过率能量比 theta deg2rad(theta_deg); sin_t2 n1 .* sin(theta) ./ N2; % 复数形式斯涅尔 cos_t2 sqrt(1 - sin_t2.^2); % 复数折射角余弦 if strcmpi(pol, s) || nargin 4 rs (n1.*cos(theta) - N2.*cos_t2) ./ (n1.*cos(theta) N2.*cos_t2); Rs abs(rs).^2; % 透过率需要考虑界面处能流密度用坡印廷矢量的法向分量 Ts real(N2 .* conj(cos_t2)) ./ (n1 .* cos(theta)) .* abs(2*n1.*cos(theta) ./ (n1.*cos(theta) N2.*cos_t2)).^2; end if strcmpi(pol, p) || nargin 4 rp (N2.*cos(theta) - n1.*cos_t2) ./ (N2.*cos(theta) n1.*cos_t2); Rp abs(rp).^2; Tp real(N2 .* conj(cos_t2)) ./ (n1 .* cos(theta)) .* abs(2*n1.*cos(theta) ./ (N2.*cos(theta) n1.*cos_t2)).^2; end end几个容易被忽略的细节一是透过率的公式里有实部运算 real(...)这是因为复数折射率下坡印廷矢量的法向分量要取实部不能直接 abs 完事二是分母里的 n1 千万别写成 N1入射介质通常是透明介质空气、玻璃、水用实数就对了。2.3 能量守恒校验判断代码是否写对的第一把尺子写完正问题函数后我强烈建议先做一组能量守恒校验。对于无吸收界面R T 1 是必然的对于吸收介质界面R T 1差额正好是介质吸收的能量占比% 快速校验脚本 n1 1.0; n2 1.5; k2 0.0; % 无吸收 theta 0:0.5:80; [Rs, Rp, Ts, Tp] fresnel_RT(n1, n2 1i*k2, theta, s); fprintf(无吸收介质s偏振 RT 的范围: %.12f ~ %.12f\n, min(RsTs), max(RsTs)); [Rs2, Rp2, Ts2, Tp2] fresnel_RT(n1, 1.5 0.2i, theta, s); fprintf(有吸收介质s偏振 RT 的范围: %.6f ~ %.6f\n, min(Rs2Ts2), max(Rs2Ts2));我第一次跑校验时发现在斜入射角下 RT 偏离 1 很离谱最后定位到是透过率公式里少乘了一个 n2/n1 的系数。这个问题如果直接拿去反演拟合出的 n、k 会莫名其妙地偏移。所以这里强烈建议任何正问题代码上线前至少跑一组 n1.5、k0 的数据确认 RT1 精确到 1e-12 级别。表 2 是我实际跑出来的一组参考数据空气入射到 N2 1.5 0.1i入射角度RsRpRs TsRp Tp00.04000.04000.96000.9600300.06100.01860.93900.981456.30.18200.00180.81800.9982750.39000.08020.61000.9198可以看到有吸收时 RT 小于 1且在布儒斯特角附近约 56 度p 偏振反射率极小但并没有像无吸收介质那样严格归零这就是 k 不为零带来的行为差异。3. 反演光学系数多角度数据下的非线性拟合3.1 为什么单个角度、单一偏振的数据不够用反演的目标是通过实验测得的反射率数据估算两个未知数n 和 k。如果只有一个入射角、单一偏振那手头只有一个反射率值一个方程两个未知数数学上就是欠定问题解不唯一。这也是很多初学者困惑我明明测了反射率为什么拟合老是多解的原因。解决办法有三个方向测多个入射角的反射率单一偏振也行但要注意避开布儒斯特角附近的极值区域同时测 s 和 p 偏振的反射率用椭偏仪测振幅比 ψ 和相位差 Δ需要的仪器门槛高一些最通用的做法是第一种多角度反射率拟合。入射角范围尽量覆盖 20 到 70 度避开接近 90 度的掠入射区域那里对表面粗糙度太敏感和布儒斯特角附近的极值点。3.2 基于 lsqnonlin 的反演完整实现MATLAB 的 Optimization Toolbox 提供了 lsqnonlin适合做最小二乘非线性拟合。它的核心思路是给定一组初始光学常数 [n0, k0]不断调整 n 和 k使模型计算出的反射率与实验反射率的残差平方和最小。function [n_fit, k_fit, resnorm] fit_nk_from_R(theta_exp, R_exp, n1, n0, k0) % fit_nk_from_R 通过多角度反射率反演复折射率 n i*k % 输入: % theta_exp 实验入射角单位度列向量 % R_exp 实验反射率与 theta_exp 同长度列向量 % n1 入射介质折射率 % n0, k0 初始猜测值比如 n01.8, k00.1 % 输出: % n_fit, k_fit 反演得到的折射率与消光系数 % resnorm 最终残差平方和 theta theta_exp(:); R_meas R_exp(:); % 定义残差函数模型计算值 - 实验值 % 这里选择 s 偏振作为示例也可以改成 p 偏振或同时用两者 res_fun (x) fresnel_RT(n1, x(1) 1i*x(2), theta, s) - R_meas; % 设置约束n 通常在 0.5~6 之间k 在 0~5 之间 lb [0.5, 0]; ub [6, 5]; % 算法选项显示迭代过程设置合适的容差 options optimoptions(lsqnonlin, ... Display, iter, ... Algorithm, trust-region-reflective, ... FunctionTolerance, 1e-10, ... StepTolerance, 1e-10, ... MaxFunctionEvaluations, 2000); % 反演 [x, resnorm] lsqnonlin(res_fun, [n0, k0], lb, ub, options); n_fit x(1); k_fit x(2); end这里最有意思的点在于菲涅尔公式本身是解析的但反演过程却需要数值优化。原因很简单——反射率对 n 和 k 的关系是高度非线性的特别是 p 偏振的反射率曲线在布儒斯特角附近对 k 的灵敏度非常低没法直接解方程。3.3 初始值怎么给决定成败的关键习惯用 lsqnonlin 的时候如果初始值给得离谱优化器很容易收敛到局部极小值或者干脆发散。我踩过不少次这个坑现在的做法是用垂直入射θ0的反射率先估算 n。垂直入射时反射率 R0 |(n-1ik)/(n1ik)|²对常见介质 k 一般远小于 n粗略可以写成 R0 ≈ (n-1)²/(n1)²反解 n 就有不错的起始值。用近带边或吸收峰处的粗糙经验来估 k或者干脆从 0.01~0.1 的小值开始。如果多组数据先用一组简单的比如垂直入射算出 n 初值再加一组斜入射数据联合拟合 n 和 k。另外建议给 k 加一个弱正则化约束把 k 的下限设为一个小正数。当材料吸收很弱时k 的最优解可能在 0 附近震荡如果不加约束反演出来可能有小的负值这在物理上不成立。4. 模拟数据验证用已知光学常数检验反演算法4.1 生成带噪声的虚拟实验数据直接拿实验数据调试反演代码容易分不清是代码问题还是测量问题。所以我习惯先用模拟数据做闭环验证设定一组真实的 n 和 k用正问题生成反射率再加噪声模拟实验误差然后用反演来恢复原始值。% 生成虚拟实验数据 n_true 1.8; k_true 0.3; theta (20:5:70); % 用正问题生成真实反射率 [Rs_true, Rp_true] fresnel_RT(1.0, n_true 1i*k_true, theta, s); % 加 1% 高斯噪声模拟测量误差 rng(42); noise_level 0.01; Rs_meas Rs_true .* (1 noise_level * randn(size(theta))); figure; plot(theta, Rs_true, b-); hold on; plot(theta, Rs_meas, ro); xlabel(入射角 (度)); ylabel(反射率); legend(真实, 加噪测量);这种做法的好处是我知道真实答案所以能清楚地检验反演代码的精度和稳定性。如果这个闭环都救不回来那就别急着上真实实验数据了。4.2 反演结果与分析用上面的模拟数据做反演初始值随便给一个偏离目标不太离谱的值比如 n01.5、k00.1跑完 lsqnonlin 之后可以在命令窗口看到最终反演值参数设定真值反演结果1%噪声相对误差n1.80001.79960.02%k0.30000.30110.37%残差平方和 resnorm 大概是 1e-5 量级说明拟合效果相当好。值得注意的是k 的相对误差比 n 大一个量级这不是偶然。反射率对折射率的灵敏度比对消光系数的灵敏度高很多尤其是当 k 较小时反射率主要由 n 决定。所以做实验设计时如果想准确测 k最好选在吸收较强的波长范围测或者在 k 贡献明显的角度区间比如大角度斜入射补充测量。4.3 噪声水平的影响什么时候反演会崩我把噪声水平从 0.5% 逐渐加到 5%继续跑同一组反演结果如下噪声水平n 反演值k 反演值n 误差k 误差0.5%1.80030.30020.02%0.07%1%1.79960.30110.02%0.37%2%1.80480.30730.27%2.43%5%1.78920.27900.60%7.00%可以看到 n 始终很稳但 k 的误差几乎和噪声水平同步增长。这说明反演问题本身对噪声是不均匀敏感的也提醒我如果实验数据的信噪比不够高反演出的 k 只能当作半定量参考不要过度解读。5. 实操中容易踩的坑每一个我都真实遇到过5.1 角度正好落在布儒斯特角附近的数值陷阱p 偏振反射率在布儒斯特角附近非常小比如无吸收时理论上为零。这时候如果用 p 偏振数据反演残差函数里会出现一个接近零的量数值上很容易被噪声淹没导致这个数据点的权重被严重放大。我遇到的情况是 n1.5 的玻璃在约 56 度角测出来反射率只有 0.003稍微有一点噪声反演就偏向 k 的虚假值。应对方法要么在反演时对这个区域的点赋低权重要么干脆裁剪掉布儒斯特角附近 ±5 度的数据。5.2 消光系数很小时的负值漂移当 k 的真实值接近零比如透明介质无约束拟合时 k 可能漂移成负值。从数学上说k-0.001 和 k0.001 算出来的反射率几乎一样所以优化器并不觉得负值有什么问题但物理上消光系数不能为负负值意味着增益除非是激光介质。解决思路有几种一是提前设 lb [0.5, 0]把 k 的下界钳制在 0二是换变量比如设 x2 k² 参与拟合这样保证任何输出都非负三是反演完成后检查 k 是否为负如果是说明材料在该波段吸收极弱直接报 k≈0 更合理。5.3 多解与局部极小值反射率反演的另一个问题是可能存在多个局部极小点。特别是当数据角度范围很窄、噪声又大时不同 n、k 组合可能给出几乎相同的反射率曲线。这就像白天看远处的山山和云都是灰白色你不一定分得清边界在哪里。我踩过的坑是 n 在 2.0 和 2.8 之间都拟合出了相近的残差最后得靠联合透射率数据才把唯一解定下来。实用的对策用多组不同起始点并行反演比如在 n0 ∈ [1.2, 1.6, 2.0, 2.4, 2.8]k0 ∈ [0.01, 0.1, 0.5] 的网格上跑一圈取残差最小的那组。这个思路简单粗暴但在 MATLAB 里几分钟就能跑完比苦想初始值靠谱得多。5.4 角度单位、复数开方分支这些低级但致命的错误角度单位sin、cos 函数默认用弧度我刚开始写代码时忘了把度数转弧度出来的曲线形状完全不对。复数开方cos(θ2) sqrt(1 - sinθ2²) 在 MATLAB 里对复数输入会自动取主值但主值不一定是物理上正确的分支。对于光从 n1 入射到 N2实部大于 0的透射情形sinθ2 的虚部为正时cosθ2 的实部应该为正这样才对应折射波向界面的另一侧传播。如果出现 cos(θ2) 的实部为负需要手动取负号翻转。我一般会在函数里加一句 assert 检查。数据类型如果 N2 是整数型比如写成 1.5 而不是 1.50iMATLAB 在某些老版本的复数运算中会给出奇怪结果建议所有光学常数统一用 double。5.5 实验数据背后的隐藏问题表面粗糙度和氧化层即便代码完全正确实验反射率和理想平坦界面的反射率也大概率存在偏差。实际的样品表面有一层天然氧化层尤其是硅、铝这类材料或者存在纳米级粗糙度这都会导致反射率偏离菲涅尔公式预测值。我做铝膜反射率测量的时候只在可见光下拟合出来的 n 和 k 和文献值差了很多后来加进一层 2 到 3 纳米的 Al2O3 氧化层模型结果才对齐。所以拿实验数据反演之前先想想样品表面到底是不是一个干净的单一界面。6. 进一步扩展从反射率到椭偏参数和薄膜体系6.1 椭偏参数 ψ 和 Δ 与菲涅尔系数的关系椭偏仪测量的不是反射率而是两个偏振方向上反射系数之比ρ rp / rs tan(ψ) · e^(iΔ)其中 ψ 反映振幅比Δ 反映相位差。关键在于即使只测一个入射角ψ 和 Δ 也是两个独立的实数值恰好能反演两个未知数 n 和 k。所以椭偏测量对 n、k 的提取效率远高于单纯反射率测量这也是半导体行业普遍用椭偏仪的原因。如果手头有椭偏数据反演代码只需把残差函数改成function rho reflection_ratio(n1, N2, theta_deg) [~, ~, Rp] fresnel_RT(n1, N2, theta_deg, p); % 但这里拿不到相位 % 实际需要返回复数振幅系数比需要把 fresnel_RT 改成同时输出幅值系数 rs, rp end实现上只要把 fresnel_RT 函数稍作修改把 rs 和 rp 的复数振幅也输出出来就能直接计算 ρ然后同时拟合 ψ 和 Δ 两个通道。6.2 薄膜体系单界面不够时的转移矩阵扩展如果样品是多层膜或者单层膜厚度和光的相干长度可比界面间的多次反射和干涉效应就会变得显著。这时候不能再用单一界面的菲涅尔公式而要用转移矩阵法TMM。TMM 的核心思想是把每一层用 2×2 矩阵表示光的传播和界面反射都写成矩阵相乘最后得到整个多层膜系的反射率。MATLAB 里实现 TMM 并不复杂核心是 2×2 矩阵连乘关键参数包括每层的复折射率 Nj、厚度 dj、入射波长 λ 和角度。如果某些层的 k 很小干涉条纹会非常明显这其实是一个优点干涉条纹的周期和位置对 n 很敏感条纹对比度衬度对 k 很敏感。利用好这个特点薄膜体系的光学常数反演反而比厚体材料更容易锁定唯一解。6.3 更鲁棒的优化策略坐标下降和全局搜索对于多层体系或者 k 特别小的体系单纯用 lsqnonlin 有时不够稳。我试过几种增强手段坐标下降法先固定 k 拟合 n再固定 n 拟合 k交替迭代三五次比二维同时拟合更不容易掉进盆地。预处理数据如果反射率曲线的动态范围很大比如从 0.01 到 0.8建议对残差做归一化处理否则小反射率区域的误差会被大反射率区域完全淹没。全局优化工具箱Global Optimization Toolbox用 particleswarm 或者 ga 先在宽范围内粗搜索再用 lsqnonlin 精修。这个方法速度稍慢但应对多解问题相当有效。我个人在做成套薄膜椭偏数据反演时经常采用宽范围粗搜 局部精修的组合方案。虽然粒子群优化不能直接保证全局最优但结合多起始点网格扫描实际可靠性比我最初只用单组初值的 lsqnonlin 高了太多。反演光学常数从来不是一个函数解出来那么简单它更像是在多个候选答案之间做判断需要物理直觉和数据质量共同把关。就我自己的使用体验来说整个菲涅尔公式计算与反演的关键不在于记住某个具体代码而在于清楚每一个复数运算背后的物理意图。正问题用能量守恒校验反问题用模拟数据闭环验证拿到实验数据前先对样品形态做合理建模这三步做到了大部分光学系数计算任务都能顺利完成。