import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.signal import find_peaks, savgol_filter, welch
from scipy.fft import fft, fftfreq, fftshift
from scipy.optimize import minimize_scalar, curve_fit
from scipy.interpolate import interp1d
import warnings
warnings.filterwarnings('ignore')
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
class EnhancedInterferenceAnalyzer:
"""增强的干涉测厚分析器 - 支持动态折射率和优化的傅里叶分析"""
def __init__(self):
print("增强干涉测厚分析器")
print("=" * 60)
print("新增功能:")
print("1. 动态折射率计算:n = f(λ, 载流子浓度)")
print("2. 优化傅里叶分析:窗函数、零填充、功率谱密度")
print("3. 多波段分析:分波段计算提高精度")
print("4. 不确定性评估:计算测量不确定度")
print("=" * 60)
# SiC基本参数
self.n0_sic = 2.6 # 基础折射率
self.wavelength_ref = 10.0 # 参考波长 (μm)
# 色散参数 (Cauchy方程: n = A + B/λ² + C/λ⁴)
self.cauchy_A = 2.6
self.cauchy_B = 0.05 # μm²
self.cauchy_C = 0.001 # μm⁴
def calculate_dynamic_refractive_index(self, wavelength_um, carrier_concentration=1e16):
"""
计算动态折射率
参数:
wavelength_um: 波长 (μm)
carrier_concentration: 载流子浓度 (cm^-3)
返回:
n: 折射率
"""
# 基础色散 (Cauchy方程)
n_base = (self.cauchy_A +
self.cauchy_B / (wavelength_um ** 2) +
self.cauchy_C / (wavelength_um ** 4))
# 载流子对折射率的影响 (Drude模型)
plasma_frequency = np.sqrt(carrier_concentration * 1.6e-19 ** 2 / (8.85e-12 * 0.5 * 9.11e-31))
wavelength_m = wavelength_um * 1e-6
frequency = 3e8 / wavelength_m
# 简化的载流子修正
carrier_correction = -plasma_frequency ** 2 / (2 * frequency ** 2) * 0.1
n_corrected = n_base + carrier_correction
return np.clip(n_corrected, 1.5, 4.0) # 限制在合理范围内
def enhanced_fourier_analysis(self, wavenumbers, reflectances):
"""
增强的傅里叶分析
改进:
1. 使用汉宁窗减少频谱泄露
2. 零填充提高频率分辨率
3. 功率谱密度分析
4. 多峰检测
"""
print("\n增强傅里叶分析:")
# 确保数据等间距
k_min, k_max = wavenumbers.min(), wavenumbers.max()
n_points = len(wavenumbers)
k_uniform = np.linspace(k_min, k_max, n_points)
# 插值到等间距网格
interp_func = interp1d(wavenumbers, reflectances, kind='cubic',
bounds_error=False, fill_value='extrapolate')
r_uniform = interp_func(k_uniform)
# 预处理:去趋势和去均值
r_detrended = self.detrend_data(k_uniform, r_uniform)
r_centered = r_detrended - np.mean(r_detrended)
# 应用汉宁窗
window = np.hanning(len(r_centered))
r_windowed = r_centered * window
# 零填充 (提高频率分辨率)
padding_factor = 4
n_padded = n_points * padding_factor
r_padded = np.zeros(n_padded)
r_padded[:n_points] = r_windowed
# FFT
fft_result = fft(r_padded)
freqs = fftfreq(n_padded, d=(k_max - k_min) / n_points)
# 功率谱密度
power_spectrum = np.abs(fft_result) ** 2
# 归一化
power_spectrum = power_spectrum / np.max(power_spectrum)
# 只考虑正频率
positive_mask = freqs > 0
positive_freqs = freqs[positive_mask]
positive_power = power_spectrum[positive_mask]
# 平滑功率谱以减少噪声
if len(positive_power) > 20:
positive_power_smooth = savgol_filter(positive_power, 11, 3)
else:
positive_power_smooth = positive_power
# 多峰检测
peak_indices, peak_props = find_peaks(
positive_power_smooth,
prominence=0.1, # 相对于最大值的10%
height=0.05, # 至少5%的高度
distance=len(positive_power_smooth) // 20 # 最小间距
)
dominant_frequencies = []
dominant_powers = []
if len(peak_indices) > 0:
# 按功率排序
sorted_indices = peak_indices[np.argsort(positive_power_smooth[peak_indices])[::-1]]
for i, peak_idx in enumerate(sorted_indices[:3]): # 取前3个主要峰
freq = positive_freqs[peak_idx]
power = positive_power_smooth[peak_idx]
period = 1 / freq if freq > 0 else 0
dominant_frequencies.append(freq)
dominant_powers.append(power)
print(f"主要频率 {i + 1}: {freq:.6f} cm, 功率: {power:.3f}")
print(f"对应周期: {period:.2f} cm^-1")
return dominant_frequencies, dominant_powers, positive_freqs, positive_power_smooth
def detrend_data(self, x, y):
"""去除数据趋势"""
# 使用二次多项式拟合趋势
coeffs = np.polyfit(x, y, 2)
trend = np.polyval(coeffs, x)
return y - trend
def multi_band_analysis(self, wavenumbers, reflectances, theta_deg):
"""
多波段分析
将光谱分为多个波段分别分析,提高精度
"""
print(f"\n多波段分析:")
# 转换为波长
wavelengths = 10000 / wavenumbers # cm^-1 to μm
# 定义波段
bands = [
(2.5, 5.0, "短波红外"),
(5.0, 10.0, "中波红外"),
(10.0, 25.0, "长波红外")
]
band_results = []
for wl_min, wl_max, band_name in bands:
# 选择波段数据
band_mask = (wavelengths >= wl_min) & (wavelengths <= wl_max)
if np.sum(band_mask) < 20: # 数据点太少
continue
band_wavenumbers = wavenumbers[band_mask]
band_reflectances = reflectances[band_mask]
band_wavelengths = wavelengths[band_mask]
print(f"\n{band_name} ({wl_min}-{wl_max} μm):")
print(f"数据点数: {len(band_wavenumbers)}")
# 计算该波段的平均折射率
avg_wavelength = np.mean(band_wavelengths)
avg_refractive_index = self.calculate_dynamic_refractive_index(avg_wavelength)
print(f"平均波长: {avg_wavelength:.2f} μm")
print(f"平均折射率: {avg_refractive_index:.3f}")
# 傅里叶分析
dominant_freqs, powers, freqs, power_spectrum = self.enhanced_fourier_analysis(
band_wavenumbers, band_reflectances
)
# 计算厚度
if dominant_freqs:
for i, freq in enumerate(dominant_freqs[:2]): # 取前两个主要频率
period = 1 / freq if freq > 0 else 0
if period > 0:
thickness = self.calculate_thickness_from_period_dynamic(
period, theta_deg, avg_refractive_index
)
band_results.append({
'band': band_name,
'wavelength_range': (wl_min, wl_max),
'avg_wavelength': avg_wavelength,
'refractive_index': avg_refractive_index,
'period': period,
'thickness': thickness,
'power': powers[i] if i < len(powers) else 0
})
print(f"厚度 (频率{i + 1}): {thickness:.3f} μm")
return band_results
def calculate_thickness_from_period_dynamic(self, period_cm_inv, theta_deg, refractive_index):
"""使用动态折射率计算厚度"""
theta_rad = np.radians(theta_deg)
# 计算层内传播角度
sin_theta1 = np.sin(theta_rad) / refractive_index
if sin_theta1 >= 1:
cos_theta1 = 0.1 # 避免全反射情况
else:
cos_theta1 = np.sqrt(1 - sin_theta1 ** 2)
# 厚度计算
thickness_cm = 1 / (2 * refractive_index * cos_theta1 * period_cm_inv)
thickness_um = thickness_cm * 1e4 # cm → μm
return thickness_um
def uncertainty_analysis(self, results_list):
"""不确定性分析"""
print(f"\n不确定性分析:")
if not results_list:
return None
thicknesses = [r['thickness'] for r in results_list]
powers = [r['power'] for r in results_list]
# 加权平均 (用功率作为权重)
total_power = sum(powers)
if total_power > 0:
weights = [p / total_power for p in powers]
weighted_mean = sum(t * w for t, w in zip(thicknesses, weights))
else:
weighted_mean = np.mean(thicknesses)
# 标准差
std_dev = np.std(thicknesses)
# 相对不确定度
relative_uncertainty = std_dev / weighted_mean * 100 if weighted_mean > 0 else 0
print(f"测量结果数: {len(thicknesses)}")
print(f"厚度范围: {min(thicknesses):.3f} - {max(thicknesses):.3f} μm")
print(f"加权平均: {weighted_mean:.3f} μm")
print(f"标准差: {std_dev:.3f} μm")
print(f"相对不确定度: {relative_uncertainty:.2f}%")
return {
'weighted_mean': weighted_mean,
'std_dev': std_dev,
'relative_uncertainty': relative_uncertainty,
'count': len(thicknesses)
}
def load_data(self, filename):
"""加载数据"""
try:
df = pd.read_excel(filename)
wavenumbers = df.iloc[:, 0].values # cm^-1
reflectances = df.iloc[:, 1].values # %
# 基本清理
valid_mask = np.isfinite(wavenumbers) & np.isfinite(reflectances)
wavenumbers = wavenumbers[valid_mask]
reflectances = reflectances[valid_mask]
# 排序
sort_idx = np.argsort(wavenumbers)
wavenumbers = wavenumbers[sort_idx]
reflectances = reflectances[sort_idx]
print(f"\n数据加载: {filename}")
print(f"数据点数: {len(wavenumbers)}")
print(f"波数范围: {wavenumbers.min():.1f} - {wavenumbers.max():.1f} cm^-1")
print(f"波长范围: {10000 / wavenumbers.max():.2f} - {10000 / wavenumbers.min():.2f} μm")
return wavenumbers, reflectances
except Exception as e:
print(f"数据加载失败: {e}")
return None, None
def comprehensive_analysis(self, filename, theta_deg):
"""综合分析"""
print(f"\n" + "=" * 80)
print(f"增强综合分析: {filename} (入射角: {theta_deg}°)")
print("=" * 80)
# 加载数据
wavenumbers, reflectances = self.load_data(filename)
if wavenumbers is None:
return None
# 多波段分析
band_results = self.multi_band_analysis(wavenumbers, reflectances, theta_deg)
# 不确定性分析
uncertainty_result = self.uncertainty_analysis(band_results)
# 整体傅里叶分析
print(f"\n整体傅里叶分析:")
dominant_freqs, powers, freqs, power_spectrum = self.enhanced_fourier_analysis(
wavenumbers, reflectances
)
# 使用平均折射率计算整体厚度
avg_wavelength = np.mean(10000 / wavenumbers)
avg_n = self.calculate_dynamic_refractive_index(avg_wavelength)
overall_thicknesses = []
if dominant_freqs:
for freq in dominant_freqs[:2]:
period = 1 / freq if freq > 0 else 0
if period > 0:
thickness = self.calculate_thickness_from_period_dynamic(
period, theta_deg, avg_n
)
overall_thicknesses.append(thickness)
# 结果汇总
results = {
'filename': filename,
'theta_deg': theta_deg,
'band_results': band_results,
'overall_thicknesses': overall_thicknesses,
'uncertainty': uncertainty_result,
'wavenumbers': wavenumbers,
'reflectances': reflectances,
'avg_refractive_index': avg_n
}
print(f"\n" + "=" * 60)
print("结果汇总:")
print("=" * 60)
if uncertainty_result:
print(f"多波段加权平均: {uncertainty_result['weighted_mean']:.3f} ± {uncertainty_result['std_dev']:.3f} μm")
print(f"相对不确定度: {uncertainty_result['relative_uncertainty']:.2f}%")
if overall_thicknesses:
print(f"整体分析结果: {np.mean(overall_thicknesses):.3f} μm")
print(f"平均折射率: {avg_n:.3f}")
# 绘图
self.plot_enhanced_results(results)
return results
def plot_enhanced_results(self, results):
"""绘制增强结果"""
fig = plt.figure(figsize=(16, 12))
wavenumbers = results['wavenumbers']
reflectances = results['reflectances']
wavelengths = 10000 / wavenumbers
# 1. 原始光谱 (双坐标轴)
ax1 = plt.subplot(3, 3, 1)
ax1.plot(wavenumbers, reflectances, 'b-', linewidth=1)
ax1.set_xlabel('波数 (cm^-1)')
ax1.set_ylabel('反射率 (%)', color='b')
ax1.set_title(f'{results["filename"]} - 原始光谱')
ax1.grid(True, alpha=0.3)
# 波长坐标轴
ax1_twin = ax1.twiny()
ax1_twin.plot(wavelengths, reflectances, 'r-', alpha=0.6)
ax1_twin.set_xlabel('波长 (μm)', color='r')
ax1_twin.tick_params(axis='x', labelcolor='r')
# 2. 动态折射率
ax2 = plt.subplot(3, 3, 2)
wl_range = np.linspace(wavelengths.min(), wavelengths.max(), 100)
n_range = [self.calculate_dynamic_refractive_index(wl) for wl in wl_range]
ax2.plot(wl_range, n_range, 'g-', linewidth=2)
ax2.set_xlabel('波长 (μm)')
ax2.set_ylabel('折射率')
ax2.set_title('动态折射率')
ax2.grid(True, alpha=0.3)
# 3. 多波段厚度结果
ax3 = plt.subplot(3, 3, 3)
if results['band_results']:
bands = [r['band'] for r in results['band_results']]
thicknesses = [r['thickness'] for r in results['band_results']]
colors = ['red', 'green', 'blue', 'orange', 'purple'][:len(bands)]
bars = ax3.bar(range(len(bands)), thicknesses, color=colors, alpha=0.7)
ax3.set_xticks(range(len(bands)))
ax3.set_xticklabels(bands, rotation=45)
ax3.set_ylabel('厚度 (μm)')
ax3.set_title('多波段厚度结果')
# 添加数值标签
for bar, thickness in zip(bars, thicknesses):
height = bar.get_height()
ax3.text(bar.get_x() + bar.get_width() / 2., height + 0.01,
f'{thickness:.3f}', ha='center', va='bottom', fontsize=8)
# 4. 傅里叶功率谱
ax4 = plt.subplot(3, 3, 4)
freqs, power_spectrum = self.enhanced_fourier_analysis(wavenumbers, reflectances)[2:4]
ax4.semilogy(1 / freqs[freqs > 0], power_spectrum[freqs > 0])
ax4.set_xlabel('周期 (cm^-1)')
ax4.set_ylabel('功率谱密度')
ax4.set_title('傅里叶功率谱')
ax4.grid(True, alpha=0.3)
# 5. 不确定性分析
ax5 = plt.subplot(3, 3, 5)
if results['uncertainty']:
unc = results['uncertainty']
mean_val = unc['weighted_mean']
std_val = unc['std_dev']
# 误差条图
ax5.errorbar([0], [mean_val], yerr=[std_val],
fmt='ro', capsize=10, capthick=2, markersize=8)
ax5.set_xlim(-0.5, 0.5)
ax5.set_ylabel('厚度 (μm)')
ax5.set_title(f'测量不确定度\n{mean_val:.3f} ± {std_val:.3f} μm')
ax5.grid(True, alpha=0.3)
ax5.set_xticks([])
# 6. 波段功率对比
ax6 = plt.subplot(3, 3, 6)
if results['band_results']:
wavelength_centers = [r['avg_wavelength'] for r in results['band_results']]
powers = [r['power'] for r in results['band_results']]
ax6.scatter(wavelength_centers, powers, s=100, alpha=0.7, c=colors[:len(wavelength_centers)])
ax6.set_xlabel('波长 (μm)')
ax6.set_ylabel('功率')
ax6.set_title('各波段功率分布')
ax6.grid(True, alpha=0.3)
# 7-9. 分波段光谱
band_plots = [ax for ax in [plt.subplot(3, 3, 7), plt.subplot(3, 3, 8), plt.subplot(3, 3, 9)]]
bands = [(2.5, 5.0, "短波"), (5.0, 10.0, "中波"), (10.0, 25.0, "长波")]
for i, ((wl_min, wl_max, band_name), ax) in enumerate(zip(bands, band_plots)):
band_mask = (wavelengths >= wl_min) & (wavelengths <= wl_max)
if np.sum(band_mask) > 0:
ax.plot(wavelengths[band_mask], reflectances[band_mask], 'b-', linewidth=1)
ax.set_xlabel('波长 (μm)')
ax.set_ylabel('反射率 (%)')
ax.set_title(f'{band_name}红外段')
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
def analyze_both_files(self):
"""分析两个文件"""
print("增强分析附件1和附件2")
print("=" * 60)
# 分析附件1
result1 = self.comprehensive_analysis("C:\\Users\\Virgo\\Desktop\\数学建模大赛\\B题资料\\cleaned_碳化硅_10度入射.xlsx", 10)
# 分析附件2
result2 = self.comprehensive_analysis("C:\\Users\\Virgo\\Desktop\\数学建模大赛\\B题资料\\cleaned_碳化硅_15度入射.xlsx", 15)
# 综合分析
if result1 and result2:
print(f"\n" + "=" * 80)
print("最终综合结果")
print("=" * 80)
unc1 = result1.get('uncertainty')
unc2 = result2.get('uncertainty')
if unc1 and unc2:
t1 = unc1['weighted_mean']
std1 = unc1['std_dev']
t2 = unc2['weighted_mean']
std2 = unc2['std_dev']
# 加权平均 (权重与不确定度成反比)
w1 = 1 / (std1 ** 2 + 1e-6)
w2 = 1 / (std2 ** 2 + 1e-6)
final_thickness = (w1 * t1 + w2 * t2) / (w1 + w2)
final_uncertainty = 1 / np.sqrt(w1 + w2)
print(f"附件1厚度: {t1:.3f} ± {std1:.3f} μm")
print(f"附件2厚度: {t2:.3f} ± {std2:.3f} μm")
print(f"加权平均厚度: {final_thickness:.3f} ± {final_uncertainty:.3f} μm")
thickness_diff = abs(t1 - t2)
print(f"厚度差异: {thickness_diff:.3f} μm")
print(f"相对差异: {thickness_diff / final_thickness * 100:.2f}%")
if thickness_diff / final_thickness < 0.05:
print("两个角度结果一致性优秀")
elif thickness_diff / final_thickness < 0.1:
print("两个角度结果一致性良好")
else:
print("两个角度结果存在差异,需进一步分析")
return result1, result2
def main():
"""主程序"""
print("问题2:碳化硅外延层厚度确定算法 (增强版)")
print("=" * 80)
analyzer = EnhancedInterferenceAnalyzer()
result1, result2 = analyzer.analyze_both_files()
print("\n增强算法完成!")
if __name__ == "__main__":
main()(对代码进行全面的优化和检查,得出一个全新的代码,)
最新发布