我正在编辑【python】代码,遇到了 【Unresolved attribute reference 'analyze_file' for class 'SiCThicknessAnalyzer'】
,请帮我检查并改正错误点。我的原始代码如下:
【import os
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import rcParams
import core
# 设置中文字体
rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
rcParams['axes.unicode_minus'] = False
def plot_results(self, file_path, wavelength, raw_refl, smooth_refl,peaks, valleys, thickness, A, B, theta_deg, n_carrier, quality):
"""绘制分析结果(包含质量指标)"""
plt.figure(figsize=(15, 12)) # 增加高度以容纳质量指标
# 原始和平滑数据
plt.subplot(3, 1, 1) # 改为3行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(3, 1, 2)
plt.plot(wavelength, raw_refl, 'bo', alpha=0.6, markersize=4, label='实验数据')
# 计算理论曲线(使用实际波长点计算)
theta_rad = np.radians(theta_deg)
n_int = self.sellmeier(wavelength) # 改为向量化计算
n_opt, *_ = self.drude(wavelength, n_int, n_carrier) # 改为向量化计算
# 折射角计算(向量化)
sin_theta_i = np.sin(theta_rad) / n_opt
valid_idx = np.abs(sin_theta_i) < 1.0
theta_i = np.arcsin(np.where(valid_idx, sin_theta_i, 0))
# 理论曲线计算(使用新的反射率公式)
R1 = ((1 - n_opt) / (1 + n_opt)) ** 2
phase = 4 * np.pi * n_opt * thickness * np.cos(theta_i) / wavelength
numerator = R1 + self.R2 - 2 * np.sqrt(R1 * self.R2) * np.cos(phase)
denominator = 1 + R1 * self.R2 - 2 * np.sqrt(R1 * self.R2) * np.cos(phase)
theory_curve = numerator / denominator
theory_curve *= 100 # 转换为百分比
# 只绘制有效的理论曲线部分
plt.plot(wavelength[valid_idx], theory_curve[valid_idx], 'r-', label='理论拟合')
plt.xlabel('波长 (μm)')
plt.ylabel('反射率')
plt.title(f'厚度: {thickness:.3f}μm, 折射率: {np.mean(n_opt):.3f}, R2: {self.R2:.3f}')
plt.legend()
plt.grid(True, alpha=0.3)
# 新增质量指标展示区域
plt.subplot(3, 1, 3)
# 隐藏坐标轴
plt.axis('off')
# 展示质量指标表格
quality_text = (
f"拟合质量评估:\n\n"
f"• 相关系数 R² = {quality['r_squared']:.4f}\n"
f"• 均方根误差 RMSE = {quality['rmse']:.4f}%\n"
f"• 折射率 = {quality['refractive_index']:.4f}\n"
f"• 入射角 = {theta_deg}°\n"
f"• 载流子浓度 = {n_carrier:.1e} cm⁻³\n"
f"• R2 (常数近似) = {self.R2:.4f}"
)
plt.text(0.5, 0.5, quality_text,
fontsize=12,
ha='center', va='center',
bbox=dict(facecolor='white', alpha=0.8))
plt.tight_layout()
# 创建规范化的保存路径
save_dir = os.path.dirname(file_path)
save_name = os.path.splitext(os.path.basename(file_path))[0] + '_analysis.png'
save_path = os.path.join(save_dir, "analysis_results", save_name)
os.makedirs(os.path.dirname(save_path), exist_ok=True)
plt.savefig(save_path, dpi=300, bbox_inches='tight')
plt.close() # 关闭图形以避免内存泄漏
def analyze_file(self, file_path, theta_deg, n_carrier, plot=True):
"""分析单个光谱文件"""
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, 折射率={quality['refractive_index']:.4f}, R²={quality['r_squared']:.4f}")
# 绘制结果
if plot:
self.plot_results(file_path, wavelength, raw_refl, smooth_refl,
peaks, valleys, thickness, A, B, theta_deg, n_carrier, quality)
return {
'thickness': thickness,
'A': A,
'B': B,
'quality': quality,
'wavelength': wavelength,
'reflectivity': raw_refl
}
except Exception as e:
print(f"分析文件 {file_path} 出错: {str(e)}")
import traceback
traceback.print_exc()
return None
# 主程序
def main():
"""主分析流程"""
files_config = [
{"path": "附件1.xlsx", "angle": 10, "carrier": 1e18, "R2": 0.20388},
{"path": "附件2.xlsx", "angle": 15, "carrier": 1e18, "R2": 0.20388}
]
results = {}
print("=" * 60)
print("碳化硅外延层厚度分析程序(基于精确干涉模型)")
print("=" * 60)
for config in files_config:
file_path = config["path"]
theta_deg = config["angle"]
n_carrier = config["carrier"]
R2 = config["R2"]
if not os.path.exists(file_path):
print(f"文件不存在: {file_path}")
continue
# 创建分析器实例,传入R2值
analyzer = core.SiCThicknessAnalyzer(n_carrier=n_carrier, R2=R2)
print(f"分析文件: {file_path}")
print(f" - 入射角: {theta_deg}°")
print(f" - 载流子浓度: {n_carrier:.1e} cm⁻³")
print(f" - R2 (常数近似): {R2:.3f}")
# 分析文件
result = analyzer.analyze_file(file_path, theta_deg, n_carrier)
if result:
results[file_path] = result
# 保留4位小数,符合题目精度要求
print(
f"分析完成: 厚度 = {result['thickness']:.4f} μm, 折射率 = {result['quality']['refractive_index']:.4f}")
else:
print(f"分析失败: {file_path}")
# 结果比较和分析
if len(results) >= 2:
print("\n" + "=" * 60)
print("结果一致性分析")
print("=" * 60)
# 提取厚度结果
thicknesses = []
refractive_indices = []
angles = []
files = []
for file_path, result in results.items():
thicknesses.append(result['thickness'])
refractive_indices.append(result['quality']['refractive_index'])
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 if mean_thickness > 0 else 0
mean_refractive = np.mean(refractive_indices)
std_refractive = np.std(refractive_indices)
# 输出结果
print(" 厚度计算结果:")
for i, (file, angle, thickness, refractive) in enumerate(zip(files, angles, thicknesses, refractive_indices)):
print(f" {file} ({angle}°): {thickness:.4f} μm, 折射率: {refractive:.4f}")
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" 平均折射率: {mean_refractive:.4f}")
print(f" 折射率标准差: {std_refractive:.4f}")
# 可靠性评估
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"折射率: {quality['refractive_index']:.4f}")
print(f"拟合质量: R² = {quality['r_squared']:.4f}, RMSE = {quality['rmse']:.3f}%")
else:
print("\n 没有文件分析成功")
return
# 灵敏度分析
print("\n" + "=" * 60)
print("灵敏度分析")
print("=" * 60)
sensitivity_results = {}
# 分析载流子浓度对厚度的影响
if "附件1.xlsx" in results:
print("\n分析附件1的厚度对载流子浓度的敏感性...")
n_carrier_range = np.logspace(17, 19, 15) # 从1e17到1e19,取15个点(对数尺度)
analyzer = core.SiCThicknessAnalyzer(n_carrier=1e18, R2=0.20388)
carrier_sensitivity = analyzer.sensitivity_analysis(
"附件1.xlsx",
theta_deg=10,
n_carrier_range=n_carrier_range,
param='carrier'
)
sensitivity_results['carrier'] = carrier_sensitivity
# 绘制灵敏度曲线
plt.figure(figsize=(12, 8))
# 厚度变化曲线
plt.subplot(2, 1, 1)
plt.semilogx(carrier_sensitivity['param_values'], carrier_sensitivity['thickness'], 'bo-')
plt.axvline(1e18, color='r', linestyle='--', label=f'标称值 (n=1e18 cm⁻³)') # 修改标签,避免使用上标
plt.xlabel('载流子浓度 n_carrier (cm⁻³)')
plt.ylabel('计算出的厚度 d (μm)')
plt.title('厚度对载流子浓度的灵敏度分析')
plt.grid(True, which="both", ls="--", alpha=0.5)
plt.legend()
# 折射率变化曲线
plt.subplot(2, 1, 2)
plt.semilogx(carrier_sensitivity['param_values'], carrier_sensitivity['refractive_index'], 'go-')
plt.axvline(1e18, color='r', linestyle='--', label=f'标称值 (n=1e18 cm⁻³)') # 修改标签,避免使用上标
plt.xlabel('载流子浓度 n_carrier (cm⁻³)')
plt.ylabel('折射率 n')
plt.title('折射率对载流子浓度的灵敏度分析')
plt.grid(True, which="both", ls="--", alpha=0.5)
plt.legend()
plt.tight_layout()
plt.savefig("analysis_results/载流子浓度灵敏度分析.png", dpi=300)
plt.show()
# 分析入射角度对厚度的影响
if "附件2.xlsx" in results:
print("\n分析附件2的厚度对入射角度的敏感性...")
angle_range = np.linspace(5, 25, 11) # 从5°到25°,取11个点
analyzer = core.SiCThicknessAnalyzer(n_carrier=1e18, R2=0.20388)
angle_sensitivity = analyzer.sensitivity_analysis(
"附件2.xlsx",
theta_deg=15,
n_carrier_range=angle_range,
param='angle'
)
sensitivity_results['angle'] = angle_sensitivity
# 绘制灵敏度曲线
plt.figure(figsize=(12, 8))
# 厚度变化曲线
plt.subplot(2, 1, 1)
plt.plot(angle_sensitivity['param_values'], angle_sensitivity['thickness'], 'bo-')
plt.axvline(15, color='r', linestyle='--', label=f'标称值 (θ=15°)')
plt.xlabel('入射角度 θ (°)')
plt.ylabel('计算出的厚度 d (μm)')
plt.title('厚度对入射角度的灵敏度分析')
plt.grid(True, ls="--", alpha=0.5)
plt.legend()
# 拟合优度变化曲线
plt.subplot(2, 1, 2)
plt.plot(angle_sensitivity['param_values'], angle_sensitivity['r_squared'], 'mo-')
plt.axvline(15, color='r', linestyle='--', label=f'标称值 (θ=15°)')
plt.xlabel('入射角度 θ (°)')
plt.ylabel('拟合优度 R²')
plt.title('拟合质量对入射角度的灵敏度分析')
plt.grid(True, ls="--", alpha=0.5)
plt.legend()
plt.tight_layout()
plt.savefig("analysis_results/入射角度灵敏度分析.png", dpi=300)
plt.show()
# 保存灵敏度分析结果
print("\n保存灵敏度分析结果...")
for param, result in sensitivity_results.items():
df = pd.DataFrame({
result['param_name']: result['param_values'],
'厚度(μm)': result['thickness'],
'折射率': result['refractive_index'],
'R²': result['r_squared']
})
df.to_csv(f"analysis_results/{param}_灵敏度分析.csv", index=False)
# 保存详细结果
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.get("angle", "N/A"),
'载流子浓度(cm⁻³)': config.get("carrier", "N/A"),
'R2': config.get("R2", "N/A"),
'厚度(μm)': result.get('thickness', "N/A"),
'折射率': result.get('quality', {}).get('refractive_index', "N/A"),
'R²': result.get('quality', {}).get('r_squared', "N/A"),
'RMSE(%)': result.get('quality', {}).get('rmse', "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.splitext(os.path.basename(file_path))[0] # 安全移除扩展名
detail_file = f"analysis_results/{filename}_详细数据.csv"
detail_data = {
'波长(μm)': result.get('wavelength', []),
'反射率(%)': result.get('reflectivity', [])
}
detail_df = pd.DataFrame(detail_data)
detail_df.to_csv(detail_file, index=False, encoding='utf-8-sig')
print(f"文件 {filename} 详细数据已保存")
# 生成分析报告
generate_report(results, files_config)
except Exception as e:
print(f"保存过程中发生错误: {str(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" R2 (常数近似): {config['R2']:.4f}\n")
f.write(f" 计算厚度: {result['thickness']:.4f} μm\n")
f.write(f" 折射率: {result['quality']['refractive_index']:.4f}\n")
f.write(f" 拟合优度 R²: {result['quality']['r_squared']:.6f}\n")
f.write(f" 均方根误差: {result['quality']['rmse']:.4f}%\n")
f.write("\n")
# 一致性分析
if len(results) >= 2:
f.write("三、一致性分析\n")
f.write("-" * 30 + "\n")
thicknesses = [result['thickness'] for result in results.values()]
refractive_indices = [result['quality']['refractive_index'] 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 if mean_thickness > 0 else 0
mean_refractive = np.mean(refractive_indices)
std_refractive = np.std(refractive_indices)
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")
f.write(f"平均折射率: {mean_refractive:.4f}\n")
f.write(f"折射率标准差: {std_refractive:.4f}\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()import numpy as np
from scipy.signal import savgol_filter, find_peaks
from scipy.optimize import least_squares
# 定义可能会用到的物理常数
PHYS_CONST = {
"e": 1.60217662e-19, # 电子电荷 (C) - 基本电荷单位
"eps0": 8.854187817e-12, # 真空介电常数 (F/m) - 真空中的电容量度
"m0": 9.10938356e-31, # 电子静止质量 (kg) - 电子的质量
"c": 2.99792458e8, # 光速 (m/s) - 真空中的光速
}
# 碳化硅外延层厚度计算模型
# 基于红外干涉法和物理光学模型计算4H-SiC外延层的厚度
class SiCThicknessAnalyzer:
# 初始化函数,设定默认参数
# n_carrier: 载流子浓度,单位cm⁻³,默认5×10¹⁸
# R2: SiC-SiC界面反射率,默认0.20338(根据文献资料)
def __init__(self, n_carrier=5e18, R2=0.20338):
self.n_carrier = n_carrier # 载流子浓度 (cm^-3)影响Drude模型修正
self.m_eff = 0.4 * PHYS_CONST["m0"] # 4H-SiC电子有效质量 (kg)
self.R2 = R2 # 常数近似法中的R2值(SiC-SiC界面反射率)- 衬底反射率
# Sellmeier方程计算4H-SiC折射率
# lambda_um: 波长,单位微米(μm)
# 返回: 对应波长的折射率
def sellmeier(self, lambda_um):
# Sellmeier系数(均为特定条件下的常数)
# 这些系数通过实验数据拟合得到,适用于4H-SiC材料
A = 1.0
B1 = 0.20075
B2 = 5.54861
B3 = 35.650066
C1 = 12.07224
C2 = 35.65066
C3 = 1268.24708
# Sellmeier方程计算折射率平方
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折射率范围(2.3~2.8)
# 确保计算结果在物理合理范围内,避免数值误差导致的异常值
n_sq = np.clip(n_sq, 2.3 ** 2, 2.8 ** 2)
return np.sqrt(n_sq) # 返回折射率
# Drude模型修正折射率
# 考虑载流子浓度对折射率的影响(自由载流子吸收效应)
# lambda_um: 波长数组,单位微米
# n_int: 本征折射率数组(Sellmeier方程计算结果)
# n_carrier: 载流子浓度数组
# 返回: 修正后的折射率及相关参数
def drude(self, lambda_um, n_int, n_carrier):
# 将输入转换为数组,确保处理标量和数组的一致性
lambda_um = np.asarray(lambda_um)
n_int = np.asarray(n_int)
n_carrier = np.asarray(n_carrier)
# 初始化输出数组
n_real = n_int.copy() # 修正后的折射率实部
delta_n = np.zeros_like(n_int) # 折射率变化量
relative_change = np.zeros_like(n_int) # 相对变化百分比
omega_p = np.zeros_like(n_int) # 等离子体频率
gamma = np.zeros_like(n_int) # 散射率
# 找出需要Drude修正的点
# 只有载流子浓度足够高时才需要修正(n_carrier >= 1e15)
mask = n_carrier >= 1e15
if np.any(mask):
# 单位转换 cm⁻³ → m⁻³
# 1 cm⁻³ = 10⁶ m⁻³(因为1 m³ = 10⁶ cm³)
n_carrier_m = n_carrier[mask] * 1e6
# 迁移率与载流子浓度关系
# 迁移率随载流子浓度增加而减小
mobility = np.zeros_like(n_carrier_m)
for i, nc in enumerate(n_carrier_m):
mobility[i] = self._get_mobility(nc / 1e6) # 转换回cm⁻³单位输入
# 等离子体频率计算
omega_p_mask = 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_mask = PHYS_CONST["e"] / (self.m_eff * mobility_si)
# 光波角频率计算
lambda_m = lambda_um[mask] * 1e-6 # μm → m
omega = 2 * np.pi * PHYS_CONST["c"] / lambda_m
# Drude模型修正介电函数
epsilon_complex = n_int[mask] ** 2 - (omega_p_mask ** 2) / (omega ** 2 + 1j * omega * gamma_mask)
# 计算复折射率
n_complex = np.sqrt(epsilon_complex)
n_real_mask = np.real(n_complex) # 取实部(折射率)
# 物理合理性检查:如果修正后的折射率大于本征折射率,则进行限制
# 选择将修正后的折射率限制在本征折射率的90%以上
# Drude修正通常使折射率减小,但数值计算可能产生异常
n_real_mask = np.maximum(n_int[mask] * 0.9, n_real_mask)
# 将修正后的值赋给输出数组
n_real[mask] = n_real_mask
delta_n[mask] = n_real_mask - n_int[mask]
relative_change[mask] = (delta_n[mask] / n_int[mask]) * 100
omega_p[mask] = omega_p_mask
gamma[mask] = gamma_mask
return n_real, delta_n, relative_change, omega_p, gamma
# 根据载流子浓度获取迁移率
# 迁移率随载流子浓度增加而减小(散射增强)
# n_carrier: 载流子浓度,单位cm⁻³
# 返回: 迁移率,单位cm²/V·s
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 # 高浓度时迁移率低
# 数据预处理(筛选有效干涉条纹)
# wave_number: 波数数据,单位cm⁻¹
# reflectivity: 反射率数据,单位%
# 返回: 处理后的波长、反射率数据及极值点位置
def preprocess_data(self, wave_number, reflectivity):
# 波数转波长
wavelength = 1e4 / wave_number
# 限制波长为3~8μm(碳化硅透明窗口)
# 此范围外数据可能受吸收影响,不适合干涉分析
valid_mask = (wavelength >= 3.0) & (wavelength <= 8.0)
wavelength = wavelength[valid_mask]
reflectivity = reflectivity[valid_mask]
if len(wavelength) == 0:
raise ValueError("无有效波长数据(需3~8μm范围)")
# 平滑处理 - 使用Savitzky-Golay滤波器去除噪声
# 保持信号特征的同时减少随机噪声
n_points = len(reflectivity)
window_length = max(11, int(n_points * 0.05)) # 动态窗口大小
if window_length % 2 == 0:
window_length += 1 # 确保窗口长度为奇数
reflectivity_smooth = savgol_filter(reflectivity, window_length=window_length, polyorder=3)
# 极值点提取 - 寻找干涉条纹的峰值和谷值
# 这些极值点对应干涉相长和相消的条件
peaks, _ = find_peaks(
reflectivity_smooth,
prominence=0.2, # 最小 prominence 确保是显著峰值
distance=10, # 最小 peak 间距
width=3 # 最小 peak 宽度
)
valleys, _ = find_peaks(
-reflectivity_smooth, # 找反射率的谷值(即-reflectivity的峰值)
prominence=0.5, # 谷值需要更大的 prominence
distance=20, # 谷值间距通常较大
width=3
)
return wavelength, reflectivity, reflectivity_smooth, peaks, valleys
# 使用干涉模型计算厚度
# wavelength: 波长数组
# reflectivity: 反射率数组
# extrema_indices: 极值点索引数组
# theta_deg: 入射角度,单位度
# n_carrier: 载流子浓度
# 返回: 优化后的厚度和拟合质量信息
def calculate_thickness(self, wavelength, reflectivity, extrema_indices, theta_deg, n_carrier=5e18):
theta_rad = np.radians(theta_deg) # 角度转弧度
initial_thicknesses = [] # 存储初步估算的厚度
refractive_indices = [] # 存储对应的折射率
# 双光束干涉的相位差修正
phase_correction = True
# 遍历所有极值点对,计算初步厚度估计
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.1:
continue
# 计算平均波长处的折射率(Sellmeier+Drude)
lambda_avg = (lambda1 + lambda2) / 2
n_int = self.sellmeier(lambda_avg)
n, _, _, _, _ = self.drude(lambda_avg, n_int, n_carrier)
refractive_indices.append(n) # 记录折射率
# 斯涅尔定律计算折射角
sin_theta_i = np.sin(theta_rad) / n
if abs(sin_theta_i) >= 1.0: # 全反射检查
continue
theta_i = np.arcsin(sin_theta_i)
# 相位差修正 - 计算厚度
# 基于相邻极值点波长差计算厚度
if phase_correction:
# 考虑相位突变的修正公式
d = (lambda1 * lambda2) / (2 * n * np.cos(theta_i) * abs(lambda2 - lambda1))
else:
# 标准公式
d = 1 / (2 * n * np.cos(theta_i) * abs(1 / lambda1 - 1 / lambda2))
# 碳化硅外延层工艺厚度范围(2~20μm)
# 过滤不合理的结果
if 2 < d < 20:
initial_thicknesses.append(d)
if not initial_thicknesses:
raise ValueError("未找到有效的厚度估计")
# 计算平均折射率
avg_refractive_index = np.mean(refractive_indices) if refractive_indices else 2.65
# 优化厚度和反射率常数
# 使用最小二乘法进一步优化初步估计结果
optimized_params = self._optimize_params(
wavelength, reflectivity,
np.median(initial_thicknesses), # 使用中位数作为初始值
theta_rad, n_carrier
)
d_opt, n_scale_opt = optimized_params
# 计算拟合质量
fit_quality = self.calculate_fit_quality(
wavelength, reflectivity, d_opt,
theta_deg, n_carrier, n_scale_opt
)
# 添加折射率到质量评估中
fit_quality['refractive_index'] = avg_refractive_index * n_scale_opt
return d_opt, fit_quality
# 优化厚度和反射率常数
# 使用最小二乘法优化参数,使理论反射率最接近实验值
def _optimize_params(self, wavelength, reflectivity, initial_d, theta_rad, n_carrier):
def objective(params):
# 参数: [厚度d, 相位φ, 折射率缩放因子n_scale]
d_opt, phi, n_scale = params
residuals = [] # 残差数组
# 每个波长单独计算折射率(避免平均波长导致的色散误差)
for i, wl in enumerate(wavelength):
n_int = self.sellmeier(wl) # 本征折射率
n_real, *_ = self.drude(wl, n_int, n_carrier) # Drude修正
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)
# 计算R1(空气-外延层界面反射率)
# 菲涅尔反射公式
R1 = ((1 - n_opt) / (1 + n_opt)) ** 2
# 使用精确的反射率公式(考虑多次反射干涉)
phase = 4 * np.pi * n_opt * d_opt * np.cos(theta_i) / wl + phi + np.pi
numerator = R1 + self.R2 - 2 * np.sqrt(R1 * self.R2) * np.cos(phase)
denominator = 1 + R1 * self.R2 - 2 * np.sqrt(R1 * self.R2) * np.cos(phase)
R_theoretical = numerator / denominator
# 转换为百分比(实验数据通常是百分比形式)
R_theoretical *= 100
# 反射率物理限制:0~100%
R_theoretical = np.clip(R_theoretical, 0, 100)
# 计算残差平方(最小二乘法的目标函数)
residuals.append((reflectivity[i] - R_theoretical) ** 2)
return np.array(residuals)
# 初始值:基于碳化硅特性设定
A_guess = np.clip(np.mean(reflectivity), 15, 35) # 反射率基线15~35%(符合空气-SiC反射特性)
B_guess = np.clip((np.max(reflectivity) - np.min(reflectivity)) / 2, 5, 20) # 振幅5~20%
initial_guess = [initial_d, A_guess, B_guess, 0, 1.0] # [厚度, A, B, 相位, 折射率缩放]
# 参数边界:严格贴合碳化硅工程范围
bounds = (
[2, 10, 2, -2 * np.pi, 0.95], # 下限:厚度≥2μm,A≥10%,B≥2%,n_scale≥0.95
[20, 40, 25, 2 * np.pi, 1.05] # 上限:厚度≤20μm,A≤40%,B≤25%,n_scale≤1.05
)
# 优化算法:使用最小二乘法
result = least_squares(
objective,
initial_guess,
bounds=bounds,
max_nfev=3000, # 最大函数评估次数
ftol=1e-8, # 函数容差
xtol=1e-8 # 参数容差
)
return result.x[0], result.x[1], result.x[2], result.x[4] # 返回优化后的参数
# 计算拟合质量指标
# 评估理论模型与实验数据的吻合程度
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:
# 计算R1(空气-外延层界面反射率)
R1 = ((1 - n_opt) / (1 + n_opt)) ** 2
# 使用精确的反射率公式
phase = 4 * np.pi * n_opt * d * np.cos(theta_i) / wl
numerator = R1 + self.R2 - 2 * np.sqrt(R1 * self.R2) * np.cos(phase)
denominator = 1 + R1 * self.R2 - 2 * np.sqrt(R1 * self.R2) * np.cos(phase)
R_theoretical = numerator / denominator
# 转换为百分比
R_theoretical *= 100
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, # 决定系数,越接近1越好
'rmse': np.sqrt(np.mean(residuals ** 2)), # 均方根误差,越小越好
'mae': np.mean(np.abs(residuals)), # 平均绝对误差,越小越好
'refractive_index': n_opt, # 优化后的折射率
'residuals_std': np.std(residuals) # 残差标准差
}
】