基于最小二乘法拟合碳化硅外延层厚度

摘要

随着国家大力发展半导体为代表的新兴产业，碳化硅外延层厚度作为第三代半导体器件性能的关键参数之一，对器件性能有着重要影响。因此，外延层厚度的测定有其必要性。本文围绕红外干涉法在外延层厚度测量中的应用展开研究，系统建立了单次反射与多光束干涉情形下的厚度反演模型，提出了基于最小二乘法与相邻极值点差法的交叉验证算法，利用了附件1 – 4完成模型的验证与精度分析。
针对问题1，基于双光束干涉效应，考虑外延层 - 衬底界面单次反射情形，构建了反射率关于掺杂载流子浓度、红外光谱的波长的Drude-Lorentz模型。依据菲涅耳公式得出理论反射率的表达式，通过非线性最小二乘拟合确定最优载流子浓度，将其回代并利用相邻极值点差法进行交叉验证，从而确保厚度测量的可靠性。

针对问题2，基于问题1建立的数学模型，首先针对数据进行预处理，有效克服了红外光谱中噪声与拼接误差带来的干扰，选取1818cm^-1以后的高波数段。将该模型应用于附件1与附件2中，反演出外延层厚度为8.5383um。通过残差重采样评估，结果表明相对误差约为0.066%，具有良好的重复性与可靠性。

针对问题3，进一步考虑多光束干涉效应，推导得出形成显著多光束干涉的必要条件：界面反射系数需足够大且光学厚度满足相长干涉条件。理论分析表明，多光束干涉会引入高频振荡分量，导致厚度反演结果系统性偏大。基于此，构建了包含多次反射项的扩展模型。应用于附件3与附件4，判定其均存在显著多光束干涉，经修正后反演厚度为3.19um，最后将复折射率中的消光系数考虑后代入模型对问题2结果计算得到8.5392um。

关键词：碳化硅；Drude-Lorentz模型；介电函数；最小二乘；相邻极值点差法

一、问题重述

1.1问题背景

随着国家大力发展半导体为代表的新兴产业，碳化硅作为第三代新型半导体材料，因其优越的综合性能备受瞩目。碳化硅外延层厚度是外延材料的关键参数之一，对器件性能有着重要影响。因此，建立一套完善的测试标准是极其必要的。

1.2问题重述

红外干涉法是一种常见的无损测量外延层厚度的方法，其基本原理是，由于外延层与衬底存在掺杂载流子浓度的差异，导致其折射率不同，外延层的折射率通常不是常数，它与掺杂载流子浓度、红外光谱的波长等参数有关。红外光入射到外延层后，一部分从外延层表面反射出来，另一部分从衬底表面反射回来，这两束光在一定条件下会产生干涉条纹。通过分析红外光谱的波长、外延层的折射率和红外光的入射角等参数确定外延层的厚度。

问题1：假设外延层和衬底界面只有一次反射、透射所产生的干涉条纹时：
（1）建立确定外延层厚度的数学模型。

问题2：依据第1问的模型：
（1）设计测量外延层厚度的算法。
（2）根据题目中所提供的碳化硅晶圆片的光谱实测数据，计算外延层厚度，分析结果的可靠性。

问题3：假设光波可以在外延层界面和衬底界面产生多次反射、透射产生多光束干涉时：
（1）探讨产生多光束干涉的必要条件，以及多光束干涉对外延层厚度计算精度可能产生的影响。
（2）依据产生多光束干涉的必要条件，分析题目中提供的硅晶圆片的测试结果是否产生多光束干涉，设计测量外延层厚度的数学模型和算法，并计算相应结果。如若多光束干涉的出现影响到了外延层厚度计算精度，尽可能消除其影响，计算消除影响后的结果。

二、问题分析

2.1问题一的分析

问题1主要要求我们给出测量碳化硅的外延层厚度的模型，探究外延层厚度与红外光谱的波长、外延层的折射率和红外光的入射角之间的关系。根据题目条件，已知红外光谱的波长和红外光的入射角，问题关键在于求出外延层折射率。由于外延层折射率与掺杂载流子浓度、红外光谱的波长等参数有关，红外波段下，碳化硅的光学性质受自由载流子与晶格振动的共同影响，Drude-Lorentz 模型能较好地描述这种复杂的介电响应，取复折射率的实部即为碳化硅外延层折射率，依据菲涅耳公式得出理论反射率的表达式，通过最小二乘法拟合得到最优载流子浓度和厚度的有效方法。最后，依据实际反射率的极值点差法计算厚度，并与拟合厚度对比，利用相邻极值点差法进行交叉验证，从而确保厚度测量的可靠性。

2.2问题二的分析

问题2主要是根据问题1建立的数学模型，计算外延层厚度，并给出可靠性分析。通过分析原始数据的特性，将数据进行归一化、Savitizky-Golay滤波去噪、PCHIP插值重采样的预处理，设定外延层厚度和载流子浓度初值，通过最小二乘法拟合得到最优载流子浓度，分析理论反射率与实验反射率之间的对误差，将拟合得到的最优载流子浓度作为已知参数，利用相邻极值点差法得到平均厚度，进一步使用残差重采样得到外延层厚度。验证三种方法误差是否在95%置信区间以内，保障了模型的可靠性。

2.3 问题三的分析

问题3主要是推导产生多光束干涉的必要条件，探讨其出现是否会影响到外延层厚度计算精度。根据查阅相关文献，多光束干涉应该先满足产生干涉现象的条件：反射率足够大，任意两相邻反射光满足相位差，满足相干长度条件，满足弱吸收条件。再根据问题2设计的算法，验证其是否对外延层厚度计算精度的影响。

三、模型假设

碳化硅外延层表面洁净，具有良好的光学表面。
忽略测试环境对测试结果可能的影响。
实验测得的数据真实可靠。
碳化硅衬底与外延层载流子掺杂浓度差异足够大。
碳化硅衬底厚度远大于外延层厚度。

四、符号说明

符号	说明
	空气介质折射率
$$n_1$	外延层的折射率
$$n_2$	衬底的折射率
	外延层几何厚度
	光程差
	相位差
	S偏振反射系数
	P偏振反射系数
	理论反射率

（剩余符号在后续有做说明）

五、模型准备

为分析光在“空气-外延层-衬底”多层结构中的传播与干涉行为，本文基于波动光学理论，建立相应的光学模型。本节主要涵盖角度计算、光程差与相位差的计算以及偏振状态的求解，具体如下：

5.1角度的确定

红外光在空气、外延层、衬底三层介质中传播时，其传播方向遵循斯涅耳定律，该定律描述了光从一种介质入射至另一种介质时，入射角与折射角的关联。

入射光由空气进入外延层，入射角为 theta_0，折射角为 theta_1，则有：

由此可得：

入射光由外延层（折射率为 n_1）进入衬底（折射率为 n_2），入射角与进入外延层的折射角相等为 theta_1，折射角为 theta_2，则有：

由此可得：

5.2光程差、相位差的计算

在红外干涉法测量外延层厚度的过程中，“空气→外延层”界面反射形成的反射光1以及“外延层→衬底”界面反射形成的反射光2。由于反射光2较反射光1额外穿过厚度为 d 的外延层，二者之间就产生了光程差，进而导致相位差，最终决定干涉增强或减弱。
1.光程差：
两束参与干涉的反射光之间的光程差由两束光的传播路径差异与半波损失带来的额外光程差一起决定。
首先需分析半波损失的影响：当光从光疏介质入射至光密介质发生反射时,反射光会发生相位突变，带来 \lambda/2 的额外光程差。由于两束光都是由光疏介质进入光密介质，所以都存在半波损失并产生 \lambda/2 的额外光程差，不过因其大小相等且方向一致，二者相互抵消。因此，不考虑半波损失带来的影响。
由下图可得，总光程差为反射光2在外延层内的往返路径与反射光1在空气中的路径差：

2.相位差：
两束光通过折射率不同的介质时，因光程差会产生相位变，即：

当光从一种介质入射到另一个介质时，会形成一个入射平面，其是由入射光线的传播方向与界面法向量所确定的平面。
根据光的电场方向，存在s偏振和p偏振两种偏振，由于光矢量振动方向的本质差异，分别垂直和平行于入射平面。它们在两种介质分界面发生反射时，反射光与入射光的复振幅比值由菲涅耳公式推导得出[1]，界面反射系数正是用于描述这两种偏振光在反射过程中振幅变化规律的关键物理量，然后对于空气与外延层上界面依据菲涅耳公式有：

将其统一到一个式子里，即：

对于分量，有

而对于分量，有

其中为有效折射率。为方便计算，后文将简写为。

各界面的菲涅尔反射系数和透射系数分别记为、、、和、、、。下标表示光从一个介质传播到另一个介质的方向，例如：

表示光从折射率为的介质入射到折射率为的介质时的反射系数，
其他符号可依此类推。

有了菲涅尔反射比和透射比，那我们就可以求出各反射光波的振幅，假设入射光的振幅为1，因此我们可以求出第1条到第\ j\ 条反射光的振幅分别为

六、模型的建立与求解

6.1 问题一模型的建立

由于折射率同时受掺杂载流子浓度和红外光波长的影响，载流子浓度不能作为定量引入，本节基于 Drude-Lorentz 理论建立了折射率与载流子浓度、红外光波长的模型。考虑到红外光的偏振特性，分别针对S偏振和P偏振光，用菲涅耳公式计算各界面的反射系数，进而得到理论反射率。随后通过最小二乘法，将理论反射率与实验反射率进行拟合，以此确定最优的载流子浓度和外延层厚度。再将最优载流子浓度重新代入上述模型，拟合得到折射率，依据波长与反射率曲线中干涉条纹的相邻极值点，利用极值点差法计算外延层的厚度。最后，将通过统计得到的实际外延层厚度与拟合得到的外延层厚度进行对比，验证计算结果的准确性，上述数学模型的流程图如图3：

6.1.1 Drude - Lorentz 模型

由题目可知，碳化硅外延层的折射率并非常数，其数值同时依赖于掺杂载流子浓度与红外光波长。若采用固定折射率值，忽略其影响将导致厚度计算结果出现偏差。因此，本节引入能够同时描述载流子浓度\ \rho\ 、波长\ \lambda\ 以及折射率三者关联的Drude - Lorentz 模型。
该模型是描述介质光学响应的经典物理模型，反应了折射率与波长、掺杂载流子浓度之间的关系，能够拟合和预测材料在宽光谱范围的光学性质，主要受两种机制影响：
1.自由载流子效应Drude项：掺杂引入的自由电子或空穴在外加光场作用下发生振荡，对介电响应产生贡献。该效应随波长增加而愈发显著，表现为折射率随波长增大而减小。
2.晶格振动与束缚电子效应Lorentz项：材料本身的晶格振动和电子在能带间的跃迁也会影响介电响应，通常在特定共振频率附近表现出强烈的色散和吸收。
综上，经典介电函数式[2][3]为：

其中，

，角频率与波长的关系；
，等离子体频率平方表达式；

各符号含义如下：

：真空中的光速；
：高频极限介电常数；
：等离子体频率；
：载流子浓度；
：元电荷；
：真空介电常数；
：载流子的有效质量；
：Drude 阻尼系数；
：第个 Lorentz 振子的共振频率；
：第个 Lorentz 振子的阻尼系数。

上述常数和系数在已有文献中给定 [4]。

通过麦克斯韦方程组与介质本构关系推导 [3]，复折射率可由介电函数得到：

其中：

实部为碳化硅外延层的折射率；
虚部为消光系数。

6.1.2 光的反射率模型

针对光在空气 - 外延层 - 衬底三层介质中的反射，需结合界面反射系数与双光束干涉原理进行分析，最终建立反射率与外延层厚度、波长、载流子浓度等参数的关系。
1.单界面反射系数
对于空气 - 外延层界面，反射系数r_{01}表达式为：
空气–外延层界面（0→1）的反射系数为：

外延层–衬底界面（1→2）的反射系数为：

2. 双光束干涉

考虑到为双光束干涉，总反射系数只需考虑第一束与第二束光的振幅比例与相位差。分别对 S 偏振和 P 偏振进行分析，可得总振幅反射系数、和总反射率、的表达式如下：

S 偏振：

由于，代入后得：

P 偏振：

同理可得：

3.总反射率的计算
实验光通常是非偏振光，它等概率包含S偏振和P偏振[5]。因此总的理论反射率为：

6.1.3 最小二乘优化模型

基于理论反射率与附件给定的反射率，我们构建最小二乘优化问题进行数值优化，通过拟合可得到最优载流子浓度以及外延层厚度，进而确定复折射率随波长的分布，目标函数如下：

其中，表示在第个波长点处测得的实验反射率数据。

6.1.4 相邻极值点差法交叉验证模型

为了进一步验证拟合得到的外延层厚度 d 的可靠性，利用最优载流子浓度计算折射率，再采用相邻极值点差法计算厚度，并与拟合结果对比。该方法基于干涉条纹的周期性特征，通过相邻峰值对应的波数差推导厚度，是红外干涉法中经典的简化测量手段，具体步骤如下：
1.最优载流子浓度求折射率
得到最优载流子浓度后，将其代入(1)式得到：

$$
\varepsilon(\lambda) = \varepsilon_\infty

\frac{ \dfrac{\rho e^2}{\varepsilon_0 m^\ast} }{ \left( \dfrac{2\pi c}{\lambda} \right)^2 + i \gamma_D \left( \dfrac{2\pi c}{\lambda} \right) }

\sum_j \frac{ S_j \omega_{0j}^2 }{ \omega_{0j}^2 - \left( \dfrac{2\pi c}{\lambda} \right)^2 - i \gamma_j \left( \dfrac{2\pi c}{\lambda} \right) }
$$

通过复折射率与介电函数的关联(2)式得到：

$$
\widetilde{n}_1^(\omega) = \sqrt{\varepsilon(\omega)} = n_1^(\omega) + i k_1(\omega)
$$

提取复折射率的实部，此为碳化硅外延层实际折射率。

2. 干涉极值条件

双光束干涉通过波峰与波谷的简单叠加，会产生相长干涉和相消干涉。设第级峰值对应的波长为，对应的外延层等效厚度为。

相长干涉：
当总光程差为波长的整数倍时，两束光同相叠加，反射率达到极大值：

相消干涉：
当总光程差为波长的半整数倍时，两束光反相叠加，反射率达到极小值：

4. 实际厚度推导

对于两个相邻的第和级峰值对应的波长和，有：

由于与相邻，波长差异较小，假设在相邻极值点之间折射率变化很小，因此我们可将它们取近似值即相邻极值点对应数据的平均值，具体如下：

$$
n_1^(\lambda_m) \approx n_1^(\lambda_{m+1}) = \bar{n}1^* = \frac{n_1^(\lambda_m) + n_1^(\lambda{m+1})}{2}
$$

联立公式 (3) 与 (4) 可得：

将式 (5) 代入式 (3)，可得第对相邻极值对应的等效厚度：

其中：

，取绝对值是为了避免因波长随级数增大而减小导致负差值；
为相邻波长处折射率的平均值；
为对应角度余弦的平均值（见前文近似）。

5. 统计得到最终的实际厚度 d

对于每一对相邻极值，代入公式 (6) 可计算出对应的等效厚度：

其中为测量到的总极值点数（即波长峰值或谷值个数），因此共有个相邻间隔。

最终，取这些厚度的算术平均值作为外延层的实际厚度：

6. 相对误差交叉验证
为验证计算解决的可靠性，将统计得到最终的实际厚度 d^\prime 和拟合得到的厚度 d 进行一个对比。由于相对误差更能反映测量结果的可信程度，本节引入了相对误差这一概念进行交叉验证。
通过计算拟合值与实际值之差与实际值的比值两者之间的偏差，其定义为：

由于测量噪声及模型近似，计算得到的厚度存在波动。通过相对误差评估结果的合理性。若相对误差在10%以内，则认为厚度计算是可靠的。
综上，上述数学模型的流程图为：

6.2 问题二模型的建立与求解

根据问题一的数学模型，可以通过编写程序得到结果，其算法流程图为

6.2.1数据分布分析

将数据导入并绘图，根据图5可知，附件1和附件2中均存在一定的噪声；存在反射率为0的数据；附件2中存在反射率大于100%的数据。部分数据显然是不符合物理意义的。因此，有必要进行数据预处理。

我们选取行业内常用的1818cm^-1以后的高波数区间[7]。该区段干涉条纹分辨率清晰，信噪比较高，存在稳定的周期，能够在减少低波数段不确定性影响的同时，提高碳化硅外延层厚度和载流子浓度拟合结果的准确性。

6.2.2数据预处理

STEP1：数据归一化
数据导入后，去除反射率为 0 的数据点，并将反射率限制在(0,\ 100]之间，为了后续计算，将反射率统一归一化到(0,\ 1]区间，以保证与理论反射率计算的一致性。

STEP2：Savitizky-Golay滤波去噪
为了抑制高频噪声并突出干涉条纹的主要特征，我们对归一化后的曲线进行平滑处理，具体采用自适应的Savitizky-Golay滤波方法，该方法通过局部多项式回归实现噪声滤除，从而避免削弱条纹峰谷结构。

STEP3：PCHIP插值重采样
由于我们去除了一些异常值，因此我们进一步查看当前数据的采样步长，发现附件一的采样步长约为0.482，附件2的采样步长为0.500，可以发现，两者的采样步长非常接近但不是一致的。为了便于对比与统一拟合，我们采用PCHIP插值方法对平滑后的两组数据以0.5cm^-1为步长重新采样。

6.2.3算法流程

Drude–Lorentz模型
首先我们根据Drude–Lorentz模型，求得复折射率，接着使用菲涅耳公式分别计算 s 偏振和 p 偏振下的界面反射系数，通过双干涉模型得到理论反射率曲线，对于实验中常见的非偏振光情形，反射率取s、p 偏振两种偏振结果的平均。
最小二乘拟合模型
然后使用最小二乘拟合，将理论反射率和实验中的反射率差值的平方作为优化目标，优化会得到拟合载流子浓度和拟合厚度，我们将得到的拟合载流子浓度回带入Drude–Lorentz模型求得折射率。
相邻极值点差法模型
紧接着通过相邻极值点差法计算厚度，方法核心是：
（1）提取反射率极值点：条纹的极值对应不同干涉阶次，每两个相邻极值对应的波数间隔与厚度相关。
（2）计算波数间隔：对极值波数求差，计算得到相邻波数的间隔、计算平均间隔和标准差。
（3）折射率修正与厚度计算：将拟合得到的载流子浓度代入折射率计算公式得到修正的折射率，通过厚度计算公式得到每个相邻极值的厚度，最后取相邻厚度的平均值作为实际厚度。

由于测量噪声及模型近似，计算得到的厚度存在波动。通过相对误差评估结果的合理性。若相对误差在10%以内，则认为厚度计算是可靠的。
在相邻极值点差法中，对于区间范围的选择，依据选取干涉条纹明显的、信噪比高的、极值间隔规律的区域。综合上述依据，选择1818-4000{cm}^{-1}波数区间，这个区间既保证条纹清晰、极值点充分，又能减少低波数段强吸收带来的不确定性，能够得到稳定且可信的厚度值。

进一步分析
为了进一步评估结果的可靠性，我们采用残差重采样的方法对拟合结果进行不确定度评估，在最小二乘拟合时，实验测得的反射率-模型计算得到的反射率就会得到一个残差，对残差做有放回的抽样，叠加到模型预测上得到伪观测，然后对伪观测重新拟合，从而获得一组拟合参数分布。
该方法的原理在于：残差可视为实验数据中未被模型解释的噪声部分，其统计特征反映了测量与模型误差。通过对残差进行重采样并叠加到模型预测值上，可在保持模型结构不变的前提下构造新的样本，进而模拟不同可能的实验观测情形。多次迭代后，拟合参数的分布即可反映模型在不确定性下的稳健性。

6.2.4模型求解

我们使用python求解计算得到载流子浓度、外延层厚度、折射率

表1 问题2的求解

载流子浓度 ()	外延层厚度 ()	折射率
	8.54	2.184

拟合得到的归一化相对误差约为6.60*10^-3，即理论反射率与实验反射率之间的平均相对偏差仅约为 0.660%。这一极小的误差表明，所建立的Drude-Lorentz单层膜模型能够准确描述外延层的光学响应，拟合结果在整个波段范围内与实验测量值高度一致，验证了模型在反映载流子浓度、外延层厚度及折射率等参数方面的可靠性和有效性。
进一步，将实验反射率与理论反射率进行对比，如下图7：

观察到实验数据与理论数据绝对反射率偏差最大不超过3%，并且两者相位是几乎重合的，说明模型拟合得很好。
（1）相邻极值点差法得到的结果：
通过相邻极值点差法计算厚度得到平均厚度为：8.348um+-2.580um。
（2）残差重采样得到的结果：
通过残差重采样方法对厚度进行统计分析。
分析残差重采样样本中载流子浓度和厚度的分布特性，如下图所示

由于残差重采样样本的分布偏离正态，采用百分位法计算 95% 置信区间，以更准确地反映参数的不确定性。
因此，得到外延层厚度为：8.68um，95%置信区间为：[6.14-11.57].
三个结果高度一致，均在可接受的误差范围内，进一步验证了外延层厚度的可靠性。

6.3 问题三模型的建立与求解

6.3.1 多光束干涉必要条件推导

1. 光学基本公式

首先给出菲涅耳反射系数和透射系数。设空气、外延层、衬底的折射率分别为、、，外延层几何厚度为。各界面的菲涅耳反射系数和透射系数分别记为：

反射系数：、、、
透射系数：、、、

下标表示光从一个介质传播到另一个介质的方向，例如表示从折射率为的介质入射到折射率为的介质时的反射系数，其余类推。
因此有：

有了菲涅尔反射比和透射比，那我们就可以求出各反射光波的振幅，假设入射光的振幅为1，因此我们可以求出第1条到第j条反射光的振幅分别为

任意两相邻反射光之间的相位差为：

由于各反射光相干，总反射振幅为所有反射波的叠加：

该级数为等比级数，公比，故可求和：

利用能量守恒关系（无吸收时）：

并结合、，可进一步简化得：

此式给出振幅为 1 的入射光经单层膜系反射后所得反射光的合振幅，单层膜的反射率则为：

2.必要条件
（1）反射率足够大
一个无穷级数

只有在收敛的情况下，总反射振幅 r 能够得到一个具体的值，即需要满足

由于，上式等价于：

（2）相干长度条件
相干长度L_c是描述光源时间相干性的物理量，定义为光波在传播方向上能保持相位相关性的最大距离差，它与光源的频谱带宽∆v存在以下傅里叶变换关系：

要观察到清晰的阶干涉条纹，第 1 束光与第束光相遇时必须仍能发生干涉。这就要求它们之间的光程差：

必须小于相干长度。

对于高阶干涉（即），这就要求相干长度必须足够长。不过，保证至少两束光能发生干涉的基本条件是：

（3）弱吸收条件
光在介质中传播的强度衰减遵循比尔-朗伯定律：

要使多光束干涉显著，必须保证足够多的高阶光束仍有足够的强度参与叠加。如果吸收太强，指数衰减因子会使高阶项迅速衰减到零,从而退化为双光束干涉。因此，弱吸收更利于多光束显著。

6.3.2外延层厚度测量的精度

在碳化硅外延层的光学测厚中，由于衬底与外延缓层界面反射率差异显著，入射光会在高反射率界面间发生多次反射，形成复杂的多光束干涉效应。不仅导致干涉条纹峰的相位发生偏移，还会引起条纹的混合叠加以及强度重分布，使得传统的双光束干涉模型无法准确反演外延层厚度。若忽视这种多光束干涉的影响，厚度计算结果会出现系统性误差，影响碳化硅外延层厚度测量的精度和可靠性。因此，需要考虑消光系数对厚度的影响。

6.3.3 模型的求解

将数据导入后，同样进行了归一化，将反射率归一化到(0,1] 区间，保证与理论反射率计算的一致性；采用Savitizky-Golay滤波去噪，抑制高频噪声、突出干涉条纹；进行PCHIP插值重采样，在保持数据单调性和主要趋势的前提下，避免插值过程中产生的伪振荡。
与问题2相比，对于Drude-Lorentz模型中的折射率，我们进一步考虑其虚部，采用复折射率的Drude-Lorentz 模型进行描述，将得到的外延层复折射率代入菲涅耳公式中得到理论反射率，理论反射率与实验反射率做最小二乘得到拟合载流子浓度和外延层厚度。

表2 硅晶外延层厚度

载流子浓度 ()	外延层厚度 ()	折射率实部	消光系数
	3.19	2.628	0.036

将复折射率中的消光系数考虑后代入模型，消除影响后，重新对问题2结果进行计算得到：

表3 碳化硅精确外延层厚度结果

载流子浓度 ()	外延层厚度 ()	折射率
	8.5392	2.184

七、模型的评价与改进

7.1 模型的评价

本文基于干涉原理和Drude–Lorentz模型构建了外延层厚度与反射率的数学模型，利用最小二乘进行参数反演。结果表明，该模型能够较好地反映实验数据中的主要干涉条纹特征，拟合得到的厚度利用多种方法进行交叉验证，均在可接受的误差范围内，说明模型具有较好的物理合理性和适用性。
从结果上来看，模型在高波数段拟合准确，能够提取厚度，残差较小。各类模型参数均有实际的物理背景，保证了模型等可解释性，并且该模型不仅适用于碳化硅的外延层，还能够推广到硅晶的外延层等其他不同材料的外延层，具有可扩展性和普适性。

7.2 模型的改进

尽管本文在实验数据拟合和厚度提取方面取得了较为理想的结果，但仅利用了高波数段计算外延层厚度，由于材料晶格振动和探测器响应的影响，模型拟合精度不足，未能充分考虑该区域的厚度信息。最小二乘拟合可能会导致模型陷入局部最优，可以利用遗传算法、粒子群优化等全局优化方法，结合局部搜索策略来提高稳健型，并且还可以引入贝叶斯优化或正则化约束来缓解参数之间的耦合，提高结果的稳健型和可靠性。

参考文献

[1] 母国光, 战元龄. 光学（第2版）[M]. 北京: 高等教育出版社, 2009.
[2] KAZAN M. Oxygen behavior in aluminum nitride[J]. Journal of Applied Physics, 2005, 98: 103529. DOI:10.1063/1.2130713.
[3] KAZAN M, OTTAVIANI L, MOUSSAED E, et al. Effect of introducing gettering sites and subsequent Au diffusion on the thermal conductivity and the free carrier concentration in n-type 4H-SiC[J]. Journal of Applied Physics, 2008, 103: 053707. DOI:10.1063/1.2883475.
[4] 马格林.一种新的SiC外延材料质量评估方法[D].西安电子科技大学,2011.
[5] Handbook of Optics. Vol. 1: Geometrical and Physical Optics, Polarized Light, Components and InstrumentsM. 3rd ed. New York: McGraw-Hill Education, 2010.
[6] Kato H S, Sato Y, Yamamoto M, et al. Surface-Specific Vibrational Spectroscopy of the Graphene Electrode/Aqueous Electrolyte Interface[J]. Journal of Physical Chemistry C, 2020, 124(21): 11447–11455. DOI:10.1021/acs.jpcc.0c01941
[7] 中关村天合宽禁带半导体技术创新联盟.4H-碳化硅同质外延层厚度的红外反射测量方法:T/IAWBS007—2018[S].北京:中国标准出版社,2018.

附录

附录A支撑材料文件列表

文件名	文件内容
`AI工具使用详情`	关于 AI 工具的使用说明
`problem1.ipynb`	问题 1 和问题 2 的求解代码与结果
`problem3.ipynb`	问题 3 的求解代码与结果
`Proview.ipynb`	4 个附件的预处理结果
`附件1.csv`	处理后的附件 1 数据
`附件2.csv`	处理后的附件 2 数据
`附件3.csv`	处理后的附件 3 数据
`附件4.csv`	处理后的附件 4 数据

附录B程序代码

import numpy as np
import pandas as pd
from scipy.optimize import minimize
from numpy.lib.scimath import arcsin

e = 1.602e-19
eps0 = 8.854e-12
m = 0.3 * 9.11e-31
c = 3e8
eps_inf = 6.5
gamma_D = 1e13

omega_0 = 1.5e14
gamma_L = 1e12
S_L = 1.0

# 测试 3个振子分别测试
lorentz_params = [
    (1.5e14, 1e12, 1.0),
    #(2.5e14, 2e12, 0.5),
    #(3.5e14, 1.5e12, 0.8)
]

def pre(df):
    lambda_m = 1e-2 / df["波数 (cm-1)"].values
    R = df["反射率 (%)"].values
    mask = R <= 1.0
    lambda_m = lambda_m[mask]
    R = R[mask]
    idx = np.argsort(lambda_m)
    return lambda_m[idx], R[idx]

df1 = pd.read_csv("附件1.csv")
df2 = pd.read_csv("附件2.csv")

# 测试 不同波数范围
#df1 = df1[df1["波数 (cm-1)"] <= 1400]
df1 = df1[df1["波数 (cm-1)"] > 1818]
#df2 = df2[df2["波数 (cm-1)"] <= 1400]
df2 = df2[df2["波数 (cm-1)"] > 1818]

lambda1, R1 = pre(df1)
lambda2, R2 = pre(df2)

theta1 = np.radians(10)
theta2 = np.radians(15)
n0 = 1.0

def drude_lorentz(lambda_m, N):
    omega = 2 * np.pi * c / lambda_m

    omega_p2 = N * e**2 / (eps0 * m)
    eps = eps_inf - omega_p2 / (omega**2 + 1j*gamma_D*omega)

    for (omega_0, gamma_L, S_L) in lorentz_params:
        eps += S_L * omega_0**2 / (omega_0**2 - omega**2 - 1j*gamma_L*omega)
    n_c = np.sqrt(eps)
    n_real = np.real(n_c)
    return n_real

def Fresnel_S_P(N, d, theta, lambda_array, n_sub=2.6, n0=1.0):
    n_c = drude_lorentz(lambda_array, N)
    theta_t = arcsin(n0 * np.sin(theta) / n_c)
    theta_s = arcsin(n_c * np.sin(theta_t) / n_sub)

    r01_s = (n0 * np.cos(theta) - n_c * np.cos(theta_t)) / (n0 * np.cos(theta) + n_c * np.cos(theta_t))
    t01_s = 2 * n0 * np.cos(theta) / (n0 * np.cos(theta) + n_c * np.cos(theta_t))
    t10_s = 2 * n_c * np.cos(theta_t) / (n0 * np.cos(theta) + n_c * np.cos(theta_t))
    r12_s = (n_c * np.cos(theta_t) - n_sub * np.cos(theta_s)) / (n_c * np.cos(theta_t) + n_sub * np.cos(theta_s))

    r01_p = (n_c * np.cos(theta) - n0 * np.cos(theta_t)) / (n_c * np.cos(theta) + n0 * np.cos(theta_t))
    t01_p = 2 * n0 * np.cos(theta) / (n_c * np.cos(theta) + n0 * np.cos(theta_t))
    t10_p = 2 * n_c * np.cos(theta_t) / (n_c * np.cos(theta) + n0 * np.cos(theta_t))
    r12_p = (n_sub * np.cos(theta_t) - n_c * np.cos(theta_s)) / (n_sub * np.cos(theta_t) + n_c * np.cos(theta_s))

    delta = 2 * np.pi * n_c * d * np.cos(theta_t) / lambda_array

    rtot_s = r01_s + t01_s * t10_s * r12_s * np.exp(2j * delta)
    rtot_p = r01_p + t01_p * t10_p * r12_p * np.exp(2j * delta)

    Rs = np.abs(rtot_s) ** 2
    Rp = np.abs(rtot_p) ** 2
    return Rs, Rp

def obj(p_log):
    N = 10**p_log[0]
    d = 10**p_log[1]
    Rs1, Rp1 = Fresnel_S_P(N, d, theta1, lambda1)
    Rs2, Rp2 = Fresnel_S_P(N, d, theta2, lambda2)
    Rth1 = 0.5*(Rs1 + Rp1)
    Rth2 = 0.5*(Rs2 + Rp2)
    
    err1 = np.mean(((R1 - Rth1) / (R1 + 1e-12))**2)
    err2 = np.mean(((R2 - Rth2) / (R2 + 1e-12))**2)
    return err1 + err2
    
res = minimize(obj, x0=[22, -5], bounds=[(2,26), (-7, -2)])
N_fit = 10**res.x[0]
d_fit = 10**res.x[1]

print(f"拟合载流子浓度 N: {N_fit:.3e} m^-3")
print(f"拟合外延层厚度 d: {d_fit*1e6} 微米")
print(f"折射率 n ≈ {drude_lorentz(10e-6, N_fit)}")

df1 = pd.read_csv("附件1.csv")
df2 = pd.read_csv("附件2.csv")

lambda1, R1 = pre(df1)
lambda2, R2 = pre(df2)
def wavenumber_to_lambda(lambda_m):
    return 1e-2 / lambda_m

def lambda_m_from_wavenumber_cm(nu_cm):
    return 1e-2 / nu_cm

# 做一半发现多余了 shit山 就这么放在这里了
def smooth(nu_cm, R, window=None, poly=None, baseline_deg=None):
    return R, R, np.zeros_like(R)

def findPeaks(nu_cm, R_corr, smooth_window=11, poly=3, height=None, distance=None, prominence=None):
    if nu_cm[0] > nu_cm[-1]:
        nu_cm = nu_cm[::-1]; R_corr = R_corr[::-1]
    R_for_peak = savgol_filter(R_corr, window_length=smooth_window, polyorder=poly) if smooth_window>3 else R_corr
    peaks, props = find_peaks(R_for_peak, height=height, distance=distance, prominence=prominence)
    return peaks, props, R_for_peak

def computeSpacing(nu_cm, peaks_idx):
    nu_peaks = nu_cm[peaks_idx]
    dnu = np.diff(nu_peaks)
    return nu_peaks, dnu, dnu.mean() if len(dnu)>0 else (nu_peaks, dnu, None)

def thickness(delta_nu_cm, n_avg, theta_rad):
    delta_nu_m = delta_nu_cm * 100.0
    d_m = 1.0 / (2.0 * n_avg * np.cos(theta_rad) * delta_nu_m)
    return d_m

def temp01(lambda_m, R, N_value, theta_rad,
                            nu_min_cm=None, nu_max_cm=None,
                            smooth_win_peak=15, peak_prominence=None, peak_distance=None):
    """
    返回： d_est_m, d_std_est_m (基于条纹间隔 std), info dict
    """
    nu_cm = wavenumber_to_lambda(lambda_m)
    if nu_cm[0] > nu_cm[-1]:
        nu_cm = nu_cm[::-1]; R = R[::-1]; lambda_m = lambda_m[::-1]
    if nu_min_cm is None: nu_min_cm = nu_cm.min()
    if nu_max_cm is None: nu_max_cm = nu_cm.max()
    mask = (nu_cm >= nu_min_cm) & (nu_cm <= nu_max_cm)
    nu_sel = nu_cm[mask]; R_sel = R[mask]; lam_sel = lambda_m[mask]
    if len(nu_sel) < 10:
        return None, None, {"error":"selected band too small"}
    R_corr, R_smoothed, baseline = smooth(nu_sel, R_sel, window=15, poly=3, baseline_deg=2)
    peaks_idx, props, R_for_peak = findPeaks(nu_sel, R_corr, smooth_window=smooth_win_peak,
                                                     height=None, distance=peak_distance, prominence=peak_prominence)
    if len(peaks_idx) < 3:
        peaks_idx, props, R_for_peak = findPeaks(nu_sel, R_corr, smooth_window=7,
                                                         height=None, distance=peak_distance, prominence=peak_prominence)
    if len(peaks_idx) < 2:
        return None, None, {"error":"not enough peaks found", "n_peaks":len(peaks_idx)}
    nu_peaks, dnu_array, dnu_mean = computeSpacing(nu_sel, peaks_idx)
    dnu_std = dnu_array.std() if len(dnu_array)>0 else 0.0
    n_c_s = drude_lorentz(lambda_m, N_value)
    lam_peaks = lambda_m_from_wavenumber_cm(nu_peaks)
    n_real_interp = np.real(np.interp(lam_peaks, lambda_m, n_c_s))
    n_avg = np.mean(n_real_interp)
    d_est = thickness(dnu_mean, n_avg, theta_rad)
    rel_err = np.sqrt( (dnu_std/dnu_mean if dnu_mean!=0 else 0.0)**2 + (np.std(n_real_interp)/n_avg if n_avg!=0 else 0.0)**2 )
    d_std = d_est * rel_err
    info = {"nu_peaks":nu_peaks, "dnu_array":dnu_array, "dnu_mean":dnu_mean, "dnu_std":dnu_std,
            "n_at_peaks":n_real_interp, "n_avg":n_avg, "peaks_idx":peaks_idx, "R_for_peak":R_for_peak,
            "nu_sel":nu_sel, "R_sel":R_sel, "lam_sel":lam_sel}
    return d_est, d_std, info

nu_min = 1000
nu_max = 4000
d_est_m, d_std_m, info = temp01(lambda1, R1, N_fit, theta1, nu_min_cm=nu_min, nu_max_cm=nu_max,
                                                 smooth_win_peak=21, peak_prominence=0.002, peak_distance=4)
if d_est_m is not None:
    print("条纹法估计 d = {:.3f} 微米 ± {:.3f} 微米".format(d_est_m*1e6, 3*d_std_m*1e6))
else:
    print("错误：", info)

def computer_R(N, d, theta, lambda_array, n_sub=2.6, n0=1.0):
    n_c = drude_lorentz(lambda_array, N)

    theta_t = arcsin(n0 * np.sin(theta) / n_c)

    theta_s = arcsin(n_c * np.sin(theta_t) / n_sub)

    r01 = (n0 * np.cos(theta) - n_c * np.cos(theta_t)) / (n0 * np.cos(theta) + n_c * np.cos(theta_t))
    r12 = (n_c * np.cos(theta_t) - n_sub * np.cos(theta_s)) / (n_c * np.cos(theta_t) + n_sub * np.cos(theta_s))

    delta = 2 * np.pi * n_c * d * np.cos(theta_t) / lambda_array

    r_total = (r01 + r12 * np.exp(2j * delta)) / (1 + r01 * r12 * np.exp(2j * delta))

    return np.abs(r_total) ** 2

def bootstrap(res, nboot=200):
    ests = []
    resids1 = R1 - computer_R(N_fit, d_fit, theta1, lambda1)
    resids2 = R2 - computer_R(N_fit, d_fit, theta2, lambda2)
    
    for i in range(nboot):
        r1 = computer_R(N_fit, d_fit, theta1, lambda1) + np.random.choice(resids1, size=len(resids1), replace=True)
        r2 = computer_R(N_fit, d_fit, theta2, lambda2) + np.random.choice(resids2, size=len(resids2), replace=True)

        def obj(p):
            N_try, d_try = 10**p[0], 10**p[1]
            err1 = np.mean((r1 - computer_R(N_try, d_try, theta1, lambda1))**2)
            err2 = np.mean((r2 - computer_R(N_try, d_try, theta2, lambda2))**2)
            return err1 + err2
        
        try:
            x0 = res.x + np.random.normal(0, 0.05, size=2)
            rr = minimize(obj, x0=x0, bounds=[(20,26),(-6,-2)])
            if rr.success:
                N_try, d_try = 10**rr.x[0], 10**rr.x[1]
                if 1e2 < N_try < 1e24 and 0.1e-6 < d_try < 50e-6:
                    ests.append((N_try, d_try))
        except:
            continue
    
    return np.array(ests)
ests = bootstrap(res, nboot=500)

if len(ests) > 0:
    N_mean, d_mean = ests.mean(axis=0)

    N_lower, N_upper = np.percentile(ests[:,0], [2.5, 97.5])
    d_lower, d_upper = np.percentile(ests[:,1], [2.5, 97.5])

    print(f"N_fit: {N_mean:.3e} (95% CI: {N_lower:.3e} - {N_upper:.3e})")
    print(f"d_fit: {d_mean*1e6:.2f} 微米 (95% CI: {d_lower*1e6:.2f} - {d_upper*1e6:.2f} 微米)")
else:
    print("无")

problem3.ipynb
import numpy as np
import pandas as pd
from scipy.optimize import minimize
from scipy.signal import savgol_filter
from numpy.lib.scimath import arcsin

e = 1.602e-19
eps0 = 8.854e-12
m = 0.3 * 9.11e-31
c = 3e8
eps_inf = 11.7
gamma_D = 1e13

omega_0 = 1.5e14
gamma_L = 1e12
S_L = 1.0

# 3个振子测试
lorentz_params = [
    (1.5e14, 1e12, 1.0),
    #(2.5e14, 2e12, 0.5),
    #(3.5e14, 1.5e12, 0.8)
]

def pre(df):
    lambda_m = 1e-2 / df["波数 (cm-1)"].values
    R = df["反射率 (%)"].values
    mask = R <= 1.0
    lambda_m = lambda_m[mask]
    R = R[mask]
    idx = np.argsort(lambda_m)
    return lambda_m[idx], R[idx]

df3 = pd.read_csv("附件3.csv")
df4 = pd.read_csv("附件4.csv")

# 测试 不同波数段
#df1 = df1[df1["波数 (cm-1)"] <= 1400]
df3 = df3[df3["波数 (cm-1)"] > 500]
#df2 = df2[df2["波数 (cm-1)"] <= 1400]
df4 = df4[df4["波数 (cm-1)"] > 500]

lambda1, R1 = pre(df3)
lambda2, R2 = pre(df4)


theta1 = np.radians(10)
theta2 = np.radians(15)
n0 = 1.0

def drude_lorentz(lambda_m, N):
    omega = 2 * np.pi * c / lambda_m

    omega_p2 = N * e**2 / (eps0 * m)
    eps = eps_inf - omega_p2 / (omega**2 + 1j*gamma_D*omega)

    for (omega_0, gamma_L, S_L) in lorentz_params:
        eps += S_L * omega_0**2 / (omega_0**2 - omega**2 - 1j*gamma_L*omega)
    n_complex = np.sqrt(eps)
    return n_complex

def Fresnel_S_P(N, d, theta, lambda_array, n_sub=3.42, n0=1.0):
    n_complex = drude_lorentz(lambda_array, N)
    theta_t = arcsin(n0 * np.sin(theta) / n_complex)
    theta_s = arcsin(n_complex * np.sin(theta_t) / n_sub)

    r01_s = (n0 * np.cos(theta) - n_complex * np.cos(theta_t)) / (n0 * np.cos(theta) + n_complex * np.cos(theta_t))
    r12_s = (n_complex * np.cos(theta_t) - n_sub * np.cos(theta_s)) / (n_complex * np.cos(theta_t) + n_sub * np.cos(theta_s))

    r01_p = (n_complex * np.cos(theta) - n0 * np.cos(theta_t)) / (n_complex * np.cos(theta) + n0 * np.cos(theta_t))
    r12_p = (n_sub * np.cos(theta_t) - n_complex * np.cos(theta_s)) / (n_sub * np.cos(theta_t) + n_complex * np.cos(theta_s))

    delta = 2 * np.pi * n_complex * d * np.cos(theta_t) / lambda_array

    rtot_s = (r01_s + r12_s * np.exp(2j * delta)) / (1 + r01_s * r12_s * np.exp(2j * delta))
    rtot_p = (r01_p + r12_p * np.exp(2j * delta)) / (1 + r01_p * r12_p * np.exp(2j * delta))

    Rs = np.abs(rtot_s) ** 2
    Rp = np.abs(rtot_p) ** 2
    return Rs, Rp

def Fresnel_S_P_1(N, d, theta, lambda_array, n_sub=3.42, n0=1.0):
    n_complex = drude_lorentz(lambda_array, N)
    theta_t = arcsin(n0 * np.sin(theta) / n_complex)
    theta_s = arcsin(n_complex * np.sin(theta_t) / n_sub)

    r01_s = (n0*np.cos(theta) - n_complex*np.cos(theta_t)) / (n0*np.cos(theta) + n_complex*np.cos(theta_t))
    r12_s = (n_complex*np.cos(theta_t) - n_sub*np.cos(theta_s)) / (n_complex*np.cos(theta_t) + n_sub*np.cos(theta_s))
    r01_p = (n_complex*np.cos(theta) - n0*np.cos(theta_t)) / (n_complex*np.cos(theta) + n0*np.cos(theta_t))
    r12_p = (n_sub*np.cos(theta_t) - n_complex*np.cos(theta_s)) / (n_sub*np.cos(theta_t) + n_complex*np.cos(theta_s))

    delta = 2*np.pi * n_complex * d * np.cos(theta_t) / lambda_array

    rtot_s = (r01_s + r12_s * np.exp(2j*delta)) / (1 + r01_s * r12_s * np.exp(2j*delta))
    rtot_p = (r01_p + r12_p * np.exp(2j*delta)) / (1 + r01_p * r12_p * np.exp(2j*delta))

    Rs = np.abs(rtot_s)**2
    Rp = np.abs(rtot_p)**2
    return Rs, Rp

def obj(params_log):
    N = 10**params_log[0]
    d = 10**params_log[1]
    Rs1, Rp1 = Fresnel_S_P(N, d, theta1, lambda1)
    Rs2, Rp2 = Fresnel_S_P(N, d, theta2, lambda2)
    Rth1 = 0.5*(Rs1 + Rp1)
    Rth2 = 0.5*(Rs2 + Rp2)

    err1 = np.mean(((R1 - Rth1) / (R1 + 1e-12))**2)
    err2 = np.mean(((R2 - Rth2) / (R2 + 1e-12))**2)
    return err1 + err2
    
res = minimize(obj, x0=[22, -5], bounds=[(1,26), (-9, -1)])
N_fit = 10**res.x[0]
d_fit = 10**res.x[1]

print(f"拟合载流子浓度 N = {N_fit:.3e} m^-3")
print(f"拟合外延层厚度 d = {d_fit*1e6:.2f} μm")
print(f"折射率 n ≈ {drude_lorentz(10e-6, N_fit)}")

import numpy as np
import pandas as pd
from scipy.optimize import minimize
from scipy.signal import savgol_filter
from numpy.lib.scimath import arcsin

e = 1.602e-19
eps0 = 8.854e-12
m = 0.3 * 9.11e-31
c = 3e8
eps_inf = 6.5
gamma_D = 1e13

omega_0 = 1.5e14
gamma_L = 1e12
S_L = 1.0

lorentz_params = [
    (1.5e14, 1e12, 1.0),   # 振子1
    #(2.5e14, 2e12, 0.5),
    #(3.5e14, 1.5e12, 0.8)
]

def pre(df):
    lambda_m = 1e-2 / df["波数 (cm-1)"].values
    R = df["反射率 (%)"].values
    mask = R <= 1.0
    lambda_m = lambda_m[mask]
    R = R[mask]
    idx = np.argsort(lambda_m)
    return lambda_m[idx], R[idx]


df3 = pd.read_csv("附件1.csv")
df4 = pd.read_csv("附件2.csv")

#df1 = df1[df1["波数 (cm-1)"] <= 1400]
df3 = df3[df3["波数 (cm-1)"] > 1818]
#df2 = df2[df2["波数 (cm-1)"] <= 1400]
df4 = df4[df4["波数 (cm-1)"] > 1818]

lambda1, R1 = pre(df3)
lambda2, R2 = pre(df4)


theta1 = np.radians(10)
theta2 = np.radians(15)
n0 = 1.0

def drude_lorentz(lambda_m, N):
    omega = 2 * np.pi * c / lambda_m
    omega_p2 = N * e**2 / (eps0 * m)
    eps = eps_inf - omega_p2 / (omega**2 + 1j*gamma_D*omega)
    for (omega_0, gamma_L, S_L) in lorentz_params:
        eps += S_L * omega_0**2 / (omega_0**2 - omega**2 - 1j*gamma_L*omega)
    n_complex = np.sqrt(eps)
    return n_complex

def Fresnel_S_P(N, d, theta, lambda_array, n_sub=2.6, n0=1.0):
    n_complex = drude_lorentz(lambda_array, N)
    theta_t = arcsin(n0 * np.sin(theta) / n_complex)
    theta_s = arcsin(n_complex * np.sin(theta_t) / n_sub)

    r01_s = (n0 * np.cos(theta) - n_complex * np.cos(theta_t)) / \
            (n0 * np.cos(theta) + n_complex * np.cos(theta_t))
    r12_s = (n_complex * np.cos(theta_t) - n_sub * np.cos(theta_s)) / \
            (n_complex * np.cos(theta_t) + n_sub * np.cos(theta_s))

    r01_p = (n_complex * np.cos(theta) - n0 * np.cos(theta_t)) / \
            (n_complex * np.cos(theta) + n0 * np.cos(theta_t))
    r12_p = (n_sub * np.cos(theta_t) - n_complex * np.cos(theta_s)) / \
            (n_sub * np.cos(theta_t) + n_complex * np.cos(theta_s))

    delta = 2 * np.pi * n_complex * d * np.cos(theta_t) / lambda_array

    rtot_s = (r01_s + r12_s * np.exp(2j * delta)) / (1 + r01_s * r12_s * np.exp(2j * delta))
    rtot_p = (r01_p + r12_p * np.exp(2j * delta)) / (1 + r01_p * r12_p * np.exp(2j * delta))

    Rs = np.abs(rtot_s) ** 2
    Rp = np.abs(rtot_p) ** 2
    return Rs, Rp

def Fresnel_S_P_1(N, d, theta, lambda_array, n_sub=2.6, n0=1.0):
    n_complex = drude_lorentz(lambda_array, N)
    theta_t = arcsin(n0 * np.sin(theta) / n_complex)
    theta_s = arcsin(n_complex * np.sin(theta_t) / n_sub)

    r01_s = (n0*np.cos(theta) - n_complex*np.cos(theta_t)) / (n0*np.cos(theta) + n_complex*np.cos(theta_t))
    r12_s = (n_complex*np.cos(theta_t) - n_sub*np.cos(theta_s)) / (n_complex*np.cos(theta_t) + n_sub*np.cos(theta_s))
    r01_p = (n_complex*np.cos(theta) - n0*np.cos(theta_t)) / (n_complex*np.cos(theta) + n0*np.cos(theta_t))
    r12_p = (n_sub*np.cos(theta_t) - n_complex*np.cos(theta_s)) / (n_sub*np.cos(theta_t) + n_complex*np.cos(theta_s))

    delta = 2*np.pi * n_complex * d * np.cos(theta_t) / lambda_array

    rtot_s = (r01_s + r12_s * np.exp(2j*delta)) / (1 + r01_s * r12_s * np.exp(2j*delta))
    rtot_p = (r01_p + r12_p * np.exp(2j*delta)) / (1 + r01_p * r12_p * np.exp(2j*delta))

    Rs = np.abs(rtot_s)**2
    Rp = np.abs(rtot_p)**2
    return Rs, Rp

def obj(params_log):
    N = 10**params_log[0]
    d = 10**params_log[1]
    Rs1, Rp1 = Fresnel_S_P(N, d, theta1, lambda1)
    Rs2, Rp2 = Fresnel_S_P(N, d, theta2, lambda2)
    Rth1 = 0.5*(Rs1 + Rp1)
    Rth2 = 0.5*(Rs2 + Rp2)
    err1 = np.mean(((R1 - Rth1) / (R1 + 1e-12))**2)
    err2 = np.mean(((R2 - Rth2) / (R2 + 1e-12))**2)
    return err1 + err2

res = minimize(obj, x0=[22, -5], bounds=[(1,26), (-9, -1)])
N_fit = 10**res.x[0]
d_fit = 10**res.x[1]

print(f"拟合载流子浓度 N ≈ {N_fit:.3e} m^-3")
print(f"拟合外延层厚度 d ≈ {d_fit*1e6} μm")
print(f"折射率 n ≈ {drude_lorentz(10e-6, N_fit)}")

智浩的Blog

基于最小二乘法拟合碳化硅外延层厚度

基于最小二乘法拟合碳化硅外延层厚度

摘要

一、问题重述

1.1问题背景

1.2问题重述

二、问题分析

2.1问题一的分析

2.2问题二的分析

2.3 问题三的分析

三、模型假设

四、符号说明

五、模型准备

5.1角度的确定

5.2光程差、相位差的计算

六、模型的建立与求解

6.1 问题一模型的建立

6.1.1 Drude - Lorentz 模型

6.1.2 光的反射率模型

6.1.3 最小二乘优化模型

6.1.4 相邻极值点差法交叉验证模型

6.2 问题二模型的建立与求解

6.2.1数据分布分析

6.2.2数据预处理

6.2.3算法流程

6.2.4模型求解

6.3 问题三模型的建立与求解

6.3.1 多光束干涉必要条件推导

6.3.2外延层厚度测量的精度

6.3.3 模型的求解

七、模型的评价与改进

7.1 模型的评价

7.2 模型的改进

参考文献

附录