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 是两者的对比，方便速查：

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 对复数运算是原生支持的，所以代码写起来并不复杂，但前提是每一步都要意识到"这里可能有虚部"。很多人在 n=1.5、k=0 的透明介质下测试代码没问题，一换成吸收介质结果就乱，就是因为忽略了复数的传播。

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偏振 R+T 的范围: %.12f ~ %.12f\n', min(Rs+Ts), max(Rs+Ts)); [Rs2, Rp2, Ts2, Tp2] = fresnel_RT(n1, 1.5 + 0.2i, theta, 's'); fprintf('有吸收介质，s偏振 R+T 的范围: %.6f ~ %.6f\n', min(Rs2+Ts2), max(Rs2+Ts2));

我第一次跑校验时发现在斜入射角下 R+T 偏离 1 很离谱，最后定位到是透过率公式里少乘了一个 n2/n1 的系数。这个问题如果直接拿去反演，拟合出的 n、k 会莫名其妙地偏移。所以这里强烈建议：任何正问题代码上线前，至少跑一组 n=1.5、k=0 的数据，确认 R+T=1 精确到 1e-12 级别。

表 2 是我实际跑出来的一组参考数据（空气入射到 N2 = 1.5 + 0.1i）：

可以看到有吸收时 R+T 小于 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 初始猜测值，比如 n0=1.8, k0=0.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 的时候，如果初始值给得离谱，优化器很容易收敛到局部极小值或者干脆发散。我踩过不少次这个坑，现在的做法是：

1. 用垂直入射（θ=0）的反射率先估算 n。垂直入射时反射率 R0 = |(n-1+ik)/(n+1+ik)|²，对常见介质 k 一般远小于 n，粗略可以写成 R0 ≈ (n-1)²/(n+1)²，反解 n 就有不错的起始值。
2. 用近带边或吸收峰处的粗糙经验来估 k，或者干脆从 0.01~0.1 的小值开始。
3. 如果多组数据，先用一组简单的（比如垂直入射）算出 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 反演结果与分析

用上面的模拟数据做反演，初始值随便给一个偏离目标不太离谱的值，比如 n0=1.5、k0=0.1，跑完 lsqnonlin 之后可以在命令窗口看到最终反演值：

残差平方和 resnorm 大概是 1e-5 量级，说明拟合效果相当好。值得注意的是，k 的相对误差比 n 大一个量级，这不是偶然。反射率对折射率的灵敏度比对消光系数的灵敏度高很多，尤其是当 k 较小时，反射率主要由 n 决定。所以做实验设计时，如果想准确测 k，最好选在吸收较强的波长范围测，或者在 k 贡献明显的角度区间（比如大角度斜入射）补充测量。

4.3 噪声水平的影响：什么时候反演会崩

我把噪声水平从 0.5% 逐渐加到 5%，继续跑同一组反演，结果如下：

可以看到 n 始终很稳，但 k 的误差几乎和噪声水平同步增长。这说明反演问题本身对噪声是不均匀敏感的，也提醒我：如果实验数据的信噪比不够高，反演出的 k 只能当作半定量参考，不要过度解读。

5. 实操中容易踩的坑：每一个我都真实遇到过

5.1 角度正好落在布儒斯特角附近的数值陷阱

p 偏振反射率在布儒斯特角附近非常小，比如无吸收时理论上为零。这时候如果用 p 偏振数据反演，残差函数里会出现一个接近零的量，数值上很容易被噪声淹没，导致这个数据点的权重被严重放大。我遇到的情况是 n=1.5 的玻璃，在约 56 度角测出来反射率只有 0.003，稍微有一点噪声，反演就偏向 k 的虚假值。

应对方法：要么在反演时对这个区域的点赋低权重，要么干脆裁剪掉布儒斯特角附近 ±5 度的数据。

5.2 消光系数很小时的负值漂移

当 k 的真实值接近零（比如透明介质），无约束拟合时 k 可能漂移成负值。从数学上说，k=-0.001 和 k=+0.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.5+0i，MATLAB 在某些老版本的复数运算中会给出奇怪结果，建议所有光学常数统一用 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 高了太多。反演光学常数从来不是"一个函数解出来"那么简单，它更像是在多个候选答案之间做判断，需要物理直觉和数据质量共同把关。

就我自己的使用体验来说，整个菲涅尔公式计算与反演的关键，不在于记住某个具体代码，而在于清楚每一个复数运算背后的物理意图。正问题用能量守恒校验，反问题用模拟数据闭环验证，拿到实验数据前先对样品形态做合理建模，这三步做到了，大部分光学系数计算任务都能顺利完成。