import os
import pandas as pd
import numpy as np
from scipy.signal import savgol_filter, find_peaks
from scipy.optimize import least_squares
import matplotlib.pyplot as plt
from matplotlib import rcParams
from scipy.stats import zscore
# 设置中文字体
rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
rcParams['axes.unicode_minus'] = False
# 物理常数
PHYS_CONST = {
"e": 1.60217662e-19, # 电子电荷 (C)
"eps0": 8.854187817e-12, # 真空介电常数 (F/m)
"m0": 9.10938356e-31, # 电子静止质量 (kg)
"c": 2.99792458e8, # 光速 (m/s)
}
class SiCThicknessAnalyzer:
def __init__(self, n_carrier=5e18):
self.n_carrier = n_carrier # 载流子浓度 (cm^-3)
self.m_eff = 0.4 * PHYS_CONST["m0"] # 4H-SiC电子有效质量 (kg)
def sellmeier(self, lambda_um):
"""使用Sellmeier方程计算4H-SiC折射率"""
# Sellmeier系数
A = 1.0
B1 = 0.20075
B2 = 5.54861
B3 = 35.650066
C1 = 12.07224
C2 = 35.65066
C3 = 1268.24708
lambda_sq = lambda_um ** 2
n_sq = A + B1 * lambda_sq / (lambda_sq + C1) + \
B2 * lambda_sq / (lambda_sq - C2) + \
B3 * lambda_sq / (lambda_sq - C3)
# 物理限制:4H-SiC折射率范围
n_sq = np.clip(n_sq, 6.5, 7.6)
return np.sqrt(n_sq)
def drude(self, lambda_um, n_int, n_carrier):
"""Drude模型修正折射率"""
if n_carrier < 1e15:
return n_int, 0, 0, 0, 0
try:
# 单位转换
n_carrier_m = n_carrier * 1e6 # cm⁻³ → m⁻³
# 迁移率与载流子浓度关系
mobility = self._get_mobility(n_carrier)
# 等离子体频率
omega_p = np.sqrt(n_carrier_m * PHYS_CONST["e"] ** 2 /
(self.m_eff * PHYS_CONST["eps0"]))
# 散射率
mobility_si = mobility * 1e-4 # cm²/V·s → m²/V·s
gamma = PHYS_CONST["e"] / (self.m_eff * mobility_si)
# 光波角频率
lambda_m = lambda_um * 1e-6
omega = 2 * np.pi * PHYS_CONST["c"] / lambda_m
# Drude模型修正
epsilon_complex = n_int ** 2 - (omega_p ** 2) / (omega ** 2 + 1j * omega * gamma)
n_complex = np.sqrt(epsilon_complex)
n_real = np.real(n_complex)
# 物理合理性检查
if n_real > n_int:
n_real = max(n_int * 0.9, n_real)
delta_n = n_real - n_int
relative_change = (delta_n / n_int) * 100 if n_int > 1e-10 else 0
return n_real, delta_n, relative_change, omega_p, gamma
except Exception as e:
print(f"Drude修正错误: {e}")
return n_int, 0, 0, 0, 0
def _get_mobility(self, n_carrier):
"""根据载流子浓度获取迁移率"""
if n_carrier < 1e16:
return 950
elif n_carrier < 5e16:
return 800
elif n_carrier < 1e17:
return 650
elif n_carrier < 5e17:
return 450
elif n_carrier < 1e18:
return 300
elif n_carrier < 5e18:
return 150
else:
return 80
def preprocess_data(self, wave_number, reflectivity):
"""数据预处理:平滑和极值点检测"""
# 波数转波长 (cm⁻¹ → μm)
wavelength = 1e4 / wave_number
# 限制波长范围 (0.8-5.5μm)
valid_mask = (wavelength >= 0.8) & (wavelength <= 5.5)
wavelength = wavelength[valid_mask]
reflectivity = reflectivity[valid_mask]
if len(wavelength) == 0:
return wavelength, reflectivity, reflectivity, [], []
# Savitzky-Golay平滑
if len(reflectivity) > 11:
reflectivity_smooth = savgol_filter(reflectivity, window_length=11, polyorder=3)
else:
reflectivity_smooth = reflectivity.copy()
# 寻找极值点
peaks, _ = find_peaks(reflectivity_smooth, prominence=0.2, distance=10)
valleys, _ = find_peaks(-reflectivity_smooth, prominence=0.2, distance=10)
return wavelength, reflectivity, reflectivity_smooth, peaks, valleys
def calculate_thickness(self, wavelength, reflectivity, extrema_indices, theta_deg, n_carrier=5e18):
"""使用干涉模型计算厚度"""
theta_rad = np.radians(theta_deg)
initial_thicknesses = []
# 使用相邻极值点计算初始厚度估计
for i in range(len(extrema_indices) - 1):
idx1, idx2 = extrema_indices[i], extrema_indices[i + 1]
lambda1, lambda2 = wavelength[idx1], wavelength[idx2]
if abs(lambda2 - lambda1) < 0.05:
continue
lambda_avg = (lambda1 + lambda2) / 2
n_int = self.sellmeier(lambda_avg)
n, _, _, _, _ = self.drude(lambda_avg, n_int, n_carrier)
sin_theta_i = np.sin(theta_rad) / n
if abs(sin_theta_i) >= 1.0:
continue
theta_i = np.arcsin(sin_theta_i)
d = 1 / (2 * n * np.cos(theta_i) * abs(1 / lambda1 - 1 / lambda2))
if 0.1 < d < 100:
initial_thicknesses.append(d)
if not initial_thicknesses:
raise ValueError("未找到有效的厚度估计")
# 优化厚度和反射率常数
optimized_params = self._optimize_params(
wavelength, reflectivity,
np.median(initial_thicknesses),
theta_rad, n_carrier
)
d_opt, A_opt, B_opt, n_scale_opt = optimized_params
# 计算拟合质量
fit_quality = self.calculate_fit_quality(
wavelength, reflectivity, d_opt,
A_opt, B_opt, theta_deg, n_carrier, n_scale_opt
)
return d_opt, A_opt, B_opt, fit_quality
def _optimize_params(self, wavelength, reflectivity, initial_d, theta_rad, n_carrier):
"""优化厚度和反射率常数"""
def objective(params):
d_opt, A, B, phi, n_scale = params
residuals = []
# 平均波长处的折射率
lambda_avg = np.mean(wavelength)
n_int = self.sellmeier(lambda_avg)
n_real, *_ = self.drude(lambda_avg, n_int, n_carrier)
n_opt = n_real * n_scale
# 折射角
sin_theta_i = np.sin(theta_rad) / n_opt
if abs(sin_theta_i) >= 1.0:
return np.full(len(wavelength), 1000) # 返回大残差
theta_i = np.arcsin(sin_theta_i)
# 计算理论反射率
for i, wl in enumerate(wavelength):
phase = 4 * np.pi * n_opt * d_opt * np.cos(theta_i) / wl + phi
R_theoretical = A + B * np.cos(phase)
residuals.append((reflectivity[i] - R_theoretical) ** 2)
return np.array(residuals)
# 初始猜测
A_guess = np.mean(reflectivity)
B_guess = (np.max(reflectivity) - np.min(reflectivity)) / 2
initial_guess = [initial_d, A_guess, B_guess, 0, 1.0]
# 参数边界
bounds = (
[0.1, 0, 0, -2 * np.pi, 0.8],
[100, 100, 100, 2 * np.pi, 1.2]
)
# 执行优化
result = least_squares(objective, initial_guess, bounds=bounds, max_nfev=1000)
return result.x[0], result.x[1], result.x[2], result.x[4] # d, A, B, n_scale
def calculate_fit_quality(self, wavelength, reflectivity, d, A, B, theta_deg, n_carrier, n_scale=1.0):
"""计算拟合质量指标"""
theta_rad = np.radians(theta_deg)
lambda_avg = np.mean(wavelength)
# 计算折射率
n_int = self.sellmeier(lambda_avg)
n_real, *_ = self.drude(lambda_avg, n_int, n_carrier)
n_opt = n_real * n_scale
# 计算折射角
sin_theta_i = np.sin(theta_rad) / n_opt
if abs(sin_theta_i) >= 1.0:
return {"error": "Invalid refraction angle"}
theta_i = np.arcsin(sin_theta_i)
# 计算理论反射率
theoretical = []
for wl in wavelength:
phase = 4 * np.pi * n_opt * d * np.cos(theta_i) / wl
R_theoretical = A + B * np.cos(phase)
theoretical.append(R_theoretical)
theoretical = np.array(theoretical)
residuals = reflectivity - theoretical
# 计算质量指标
ss_res = np.sum(residuals ** 2)
ss_tot = np.sum((reflectivity - np.mean(reflectivity)) ** 2)
return {
'r_squared': 1 - (ss_res / ss_tot) if ss_tot > 0 else 0,
'rmse': np.sqrt(np.mean(residuals ** 2)),
'mae': np.mean(np.abs(residuals)),
'refractive_index': n_opt,
'residuals_std': np.std(residuals)
}
def analyze_file(self, file_path, theta_deg, n_carrier):
"""分析单个光谱文件"""
try:
# 读取数据
data = pd.read_excel(file_path)
wave_number = data.iloc[:, 0].values
reflectivity = data.iloc[:, 1].values
# 预处理
wavelength, raw_refl, smooth_refl, peaks, valleys = self.preprocess_data(
wave_number, reflectivity
)
print(f"找到峰值: {len(peaks)}, 谷值: {len(valleys)}")
# 合并极值点
extrema = np.sort(np.concatenate((peaks, valleys)))
if len(extrema) < 2:
raise ValueError("极值点不足,无法计算厚度")
# 计算厚度
thickness, A, B, quality = self.calculate_thickness(
wavelength, smooth_refl, extrema, theta_deg, n_carrier
)
print(f"计算结果: 厚度={thickness:.4f}μm, R²={quality['r_squared']:.4f}")
# 绘制结果
self.plot_results(file_path, wavelength, raw_refl, smooth_refl,
peaks, valleys, thickness, A, B, theta_deg, n_carrier)
return {
'thickness': thickness,
'A': A,
'B': B,
'quality': quality,
'wavelength': wavelength,
'reflectivity': raw_refl
}
except Exception as e:
print(f"分析文件 {file_path} 出错: {str(e)}")
return None
def plot_results(self, file_path, wavelength, raw_refl, smooth_refl,
peaks, valleys, thickness, A, B, theta_deg, n_carrier):
"""绘制分析结果"""
plt.figure(figsize=(15, 10))
# 原始和平滑数据
plt.subplot(2, 1, 1)
plt.plot(wavelength, raw_refl, 'b-', alpha=0.6, label='原始数据')
plt.plot(wavelength, smooth_refl, 'r-', label='平滑数据')
# 标记极值点
if len(peaks) > 0:
plt.plot(wavelength[peaks], smooth_refl[peaks], 'go', markersize=6, label='峰值')
if len(valleys) > 0:
plt.plot(wavelength[valleys], smooth_refl[valleys], 'ro', markersize=6, label='谷值')
plt.xlabel('波长 (μm)')
plt.ylabel('反射率')
plt.title(f'{os.path.basename(file_path)} - 预处理结果')
plt.legend()
plt.grid(True, alpha=0.3)
# 理论拟合
plt.subplot(2, 1, 2)
plt.plot(wavelength, raw_refl, 'bo', alpha=0.6, markersize=4, label='实验数据')
# 计算理论曲线
theta_rad = np.radians(theta_deg)
lambda_avg = np.mean(wavelength)
n_int = self.sellmeier(lambda_avg)
n_opt, *_ = self.drude(lambda_avg, n_int, n_carrier)
# 折射角
sin_theta_i = np.sin(theta_rad) / n_opt
if abs(sin_theta_i) < 1.0:
theta_i = np.arcsin(sin_theta_i)
wl_fine = np.linspace(min(wavelength), max(wavelength), 500)
phase_fine = 4 * np.pi * n_opt * thickness * np.cos(theta_i) / wl_fine
theory_fine = A + B * np.cos(phase_fine)
plt.plot(wl_fine, theory_fine, 'r-', label='理论拟合')
plt.xlabel('波长 (μm)')
plt.ylabel('反射率')
plt.title(f'厚度: {thickness:.3f}μm, 折射率: {n_opt:.3f}')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig(f'{file_path}_analysis.png', dpi=300, bbox_inches='tight')
plt.show()
# 主程序
def main():
"""主分析流程"""
# 分析文件配置
files_config = [
{"path": "附件1.xlsx", "angle": 10, "carrier": 5e18},
{"path": "附件2.xlsx", "angle": 15, "carrier": 5e18}
]
results = {}
print("=" * 60)
print("碳化硅外延层厚度分析程序")
print("=" * 60)
# 分析每个文件
for config in files_config:
file_path = config["path"]
theta_deg = config["angle"]
n_carrier = config["carrier"]
if not os.path.exists(file_path):
print(f"❌ 文件不存在: {file_path}")
continue
print(f"\n🔍 正在分析: {file_path}")
print(f" - 入射角: {theta_deg}°")
print(f" - 载流子浓度: {n_carrier:.1e} cm⁻³")
# 创建分析器实例
analyzer = SiCThicknessAnalyzer(n_carrier=n_carrier)
# 分析文件
result = analyzer.analyze_file(file_path, theta_deg, n_carrier)
if result:
results[file_path] = result
print(f"✅ 分析完成: 厚度 = {result['thickness']:.4f} μm")
else:
print(f"❌ 分析失败: {file_path}")
# 结果比较和分析
if len(results) >= 2:
print("\n" + "=" * 60)
print("结果一致性分析")
print("=" * 60)
# 提取厚度结果
thicknesses = []
angles = []
files = []
for file_path, result in results.items():
thicknesses.append(result['thickness'])
angles.append(next(cfg["angle"] for cfg in files_config if cfg["path"] == file_path))
files.append(os.path.basename(file_path))
# 计算统计量
mean_thickness = np.mean(thicknesses)
std_thickness = np.std(thicknesses)
max_diff = max(thicknesses) - min(thicknesses)
relative_diff = (max_diff / mean_thickness) * 100
# 输出结果
print("📊 厚度计算结果:")
for i, (file, angle, thickness) in enumerate(zip(files, angles, thicknesses)):
print(f" {file} ({angle}°): {thickness:.4f} μm")
print(f"\n📈 统计指标:")
print(f" 平均厚度: {mean_thickness:.4f} μm")
print(f" 标准差: {std_thickness:.4f} μm")
print(f" 最大差异: {max_diff:.4f} μm")
print(f" 相对差异: {relative_diff:.2f}%")
# 可靠性评估
print(f"\n✅ 可靠性评估:")
if relative_diff < 3:
print(" 优秀 - 结果高度一致")
elif relative_diff < 5:
print(" 良好 - 结果基本一致")
elif relative_diff < 10:
print(" 一般 - 存在一定差异")
else:
print(" 较差 - 差异较大,需要检查")
# 拟合质量评估
print(f"\n🎯 拟合质量:")
for file_path, result in results.items():
quality = result['quality']
print(f" {os.path.basename(file_path)}: R² = {quality['r_squared']:.4f}, "
f"RMSE = {quality['rmse']:.3f}%")
elif len(results) == 1:
print("\nℹ️ 只有一个文件分析成功,无法进行一致性比较")
file_path, result = list(results.items())[0]
quality = result['quality']
print(f"厚度: {result['thickness']:.4f} μm")
print(f"拟合质量: R² = {quality['r_squared']:.4f}, RMSE = {quality['rmse']:.3f}%")
else:
print("\n❌ 没有文件分析成功")
return
# 保存详细结果
save_detailed_results(results, files_config)
print("\n" + "=" * 60)
print("分析完成!")
print("=" * 60)
def save_detailed_results(results, files_config):
"""保存详细分析结果到文件"""
try:
# 创建结果目录
os.makedirs("analysis_results", exist_ok=True)
# 保存汇总结果
summary_data = []
for config in files_config:
file_path = config["path"]
if file_path in results:
result = results[file_path]
summary_data.append({
'文件': os.path.basename(file_path),
'入射角(°)': config["angle"],
'载流子浓度(cm⁻³)': config["carrier"],
'厚度(μm)': result['thickness'],
'R²': result['quality']['r_squared'],
'RMSE(%)': result['quality']['rmse'],
'折射率': result['quality'].get('refractive_index', 'N/A')
})
if summary_data:
summary_df = pd.DataFrame(summary_data)
summary_file = "analysis_results/厚度分析汇总.csv"
summary_df.to_csv(summary_file, index=False, encoding='utf-8-sig')
print(f"💾 结果已保存至: {summary_file}")
# 保存每个文件的详细数据
for file_path, result in results.items():
filename = os.path.basename(file_path).replace('.xlsx', '')
detail_file = f"analysis_results/{filename}_详细数据.csv"
detail_data = {
'波长(μm)': result['wavelength'],
'反射率(%)': result['reflectivity']
}
detail_df = pd.DataFrame(detail_data)
detail_df.to_csv(detail_file, index=False, encoding='utf-8-sig')
# 生成分析报告
generate_report(results, files_config)
except Exception as e:
print(f"保存结果时出错: {e}")
def generate_report(results, files_config):
"""生成分析报告"""
report_file = "analysis_results/分析报告.txt"
with open(report_file, 'w', encoding='utf-8') as f:
f.write("碳化硅外延层厚度分析报告\n")
f.write("=" * 50 + "\n\n")
f.write("一、分析概述\n")
f.write("-" * 30 + "\n")
f.write(f"分析时间: {pd.Timestamp.now().strftime('%Y-%m-%d %H:%M:%S')}\n")
f.write(f"分析文件数: {len(results)}/{len(files_config)}\n\n")
f.write("二、详细结果\n")
f.write("-" * 30 + "\n")
for config in files_config:
file_path = config["path"]
if file_path in results:
result = results[file_path]
f.write(f"文件: {os.path.basename(file_path)}\n")
f.write(f" 入射角: {config['angle']}°\n")
f.write(f" 载流子浓度: {config['carrier']:.1e} cm⁻³\n")
f.write(f" 计算厚度: {result['thickness']:.4f} μm\n")
f.write(f" 拟合优度 R²: {result['quality']['r_squared']:.6f}\n")
f.write(f" 均方根误差: {result['quality']['rmse']:.4f}%\n")
if 'refractive_index' in result['quality']:
f.write(f" 折射率: {result['quality']['refractive_index']:.4f}\n")
f.write("\n")
# 一致性分析
if len(results) >= 2:
f.write("三、一致性分析\n")
f.write("-" * 30 + "\n")
thicknesses = [result['thickness'] for result in results.values()]
mean_thickness = np.mean(thicknesses)
std_thickness = np.std(thicknesses)
max_diff = max(thicknesses) - min(thicknesses)
relative_diff = (max_diff / mean_thickness) * 100
f.write(f"平均厚度: {mean_thickness:.4f} μm\n")
f.write(f"厚度标准差: {std_thickness:.4f} μm\n")
f.write(f"最大差异: {max_diff:.4f} μm\n")
f.write(f"相对差异: {relative_diff:.2f}%\n\n")
f.write("四、可靠性评估\n")
f.write("-" * 30 + "\n")
if relative_diff < 3:
f.write("优秀 - 不同入射角下的结果高度一致,可靠性很高\n")
elif relative_diff < 5:
f.write("良好 - 结果基本一致,可靠性较好\n")
elif relative_diff < 10:
f.write("一般 - 存在一定差异,建议进一步分析\n")
else:
f.write("较差 - 差异较大,需要检查数据质量和模型参数\n")
f.write("\n五、结论\n")
f.write("-" * 30 + "\n")
if len(results) >= 2:
f.write("基于红外干涉法成功测得了碳化硅外延层的厚度。")
f.write("不同入射角下的测量结果表现出良好的一致性,")
f.write("表明所建立的数学模型和算法是可靠的。\n")
else:
f.write("成功完成了碳化硅外延层厚度的测量。")
f.write("建议补充其他入射角的数据以进行交叉验证。\n")
print(f"📝 分析报告已保存至: {report_file}")
# 运行示例
if __name__ == "__main__":
# 检查文件是否存在
required_files = ["附件1.xlsx", "附件2.xlsx"]
missing_files = [f for f in required_files if not os.path.exists(f)]
if missing_files:
print("❌ 缺少必要文件:")
for f in missing_files:
print(f" - {f}")
print("\n请确保所有附件文件都在当前目录下")
else:
main()
为什么这个模型的拟合度很低,是因为数据预处理吗