25高教社杯數模國賽【B題超高質量思路+問題分析】

注：本內容由”數模加油站“ 原創出品，雖無償分享，但創作不易。
歡迎參考teach，但請勿抄襲、盜賣或商用。

B 題碳化硅外延層厚度的確定
在這里插入圖片描述
碳化硅作為一種新興的第三代半導體材料，以其優越的綜合性能表現正在受到越來越多的關注。碳化硅外延層的厚度是外延材料的關鍵參數之一，對器件性能有重要影響。因此，制定一套科學、準確、可靠的碳化硅外延層厚度測試標準顯得尤為重要。
紅外干涉法是外延層厚度測量的無損傷測量方法，其工作原理是，外延層與襯底因摻雜載流子濃度的不同而有不同的折射率，紅外光入射到外延層后，一部分從外延層表面反射出來，另一部分從襯底表面反射回來（圖1），這兩束光在一定條件下會產生干涉條紋。可根據紅外光譜的波長、外延層的折射率和紅外光的入射角等參數確定外延層的厚度。
通常外延層的折射率不是常數，它與摻雜載流子的濃度、紅外光譜的波長等參數有關。
在這里插入圖片描述
問題1 如果考慮外延層和襯底界面只有一次反射、透射所產生的干涉條紋的情形（圖1），建立確定外延層厚度的數學模型。
問題 1 分析
問題1的核心目標是基于最簡化的光學干涉現象建立厚度與條紋間距的定量關系。在外延層結構中，入射光在空氣–外延層界面產生一次反射，同時另一部分光透射進入外延層，在外延層–襯底界面發生一次反射并返回空氣，形成兩束相干光。兩束光的光程差決定了干涉條紋的形成規律，因此需要首先從物理機制上進行推導。
具體而言，應當以斯涅爾定律和菲涅耳公式為基礎，明確兩束光的相位差表達式。該相位差由兩部分組成：其一是外延層內部傳播路徑引入的光程差 $Δφ=4πn1dcos?θtν~\Delta \varphi = 4 \pi n_1 d \cos\theta_t \tilde{\nu}$ ；其二是界面反射可能帶來的附加相移 $ψ(ν~)\psi(\tilde{\nu})$ 。由此，可以得到條紋極大和極小的分布條件，并進一步推出厚度與條紋間距之間的公式關系。
在數值實現上，可形成三類計算思路。第一類是基于相鄰條紋間距的直接解析計算，即通過提取相鄰極值點波數差，利用解析公式求得厚度。這種方法直觀簡便，但易受單點誤差影響。第二類是基于多點極值的線性回歸方法，即將多組極值點代入線性化方程 $m = 2 d x + c$ ，通過最小二乘回歸估計厚度，從而在統計意義上減小誤差。第三類是全譜非線性最小二乘擬合方法，即直接對整個光譜曲線擬合干涉模型，避免依賴極值提取，能夠最大程度利用實驗數據的信息。
因此，問題1的思路重點在于從物理機理出發，推導厚度與干涉條紋之間的解析關系，并提出不同層次的厚度反演方法，為后續結合實驗數據的數值計算奠定理論基礎。

解題思路：
干涉條紋的相位關系（兩光束情形）

在紅外干涉測厚的基本實驗構型中，考慮空氣–外延層–襯底三層介質結構。外部入射光為單色紅外光，波數為 $ν~=1/λ\tilde{\nu}=1/\lambda$ 。入射光自空氣（折射率 $n0≈1n_0\approx 1$ ）以入射角 $θ\theta$ 入射到外延層表面（折射率 $n1(ν~)n_1(\tilde{\nu})$ ），發生部分反射與透射。其中，表面反射所形成的第一束反射光振幅與 $r01(pol)r_{01}^{(\mathrm{pol})}$ 成正比；透射進入外延層的光在傳播距離 $2dcos?θt(ν~)2d\cos\theta_t(\tilde{\nu})$ 后，于“外延層–襯底”界面處再次反射并透射回空氣，形成第二束光。兩束光在探測器處發生疊加干涉。
由光學干涉原理可知，兩光束之間的相位差來源于兩個方面：其一是外延層內部傳播引入的光程差，其二是不同界面反射所導致的相位突變。對于傳播相位部分，可寫作在這里插入圖片描述

另一方面，在“低折射率–高折射率”界面的反射過程中會產生 $π\pi$ 的附加相移，而在“高–低”界面則無此相移。為了統一處理，不妨將來自兩界面反射的相移差額記為在這里插入圖片描述
此處的 $Δφ\Delta\varphi$ 是決定干涉條紋位置與周期性的核心量，厚度 $d$ 的信息正是通過此式反映在光譜干涉數據中的

反射/透射系數與反射率表達
為了得到可直接與實驗反射率光譜比較的理論模型，需要在振幅層次上引入菲涅耳公式。對于一般的偏振情況， $s$ 偏振與 $p$ 偏振的界面反射、透射系數分別為
在這里插入圖片描述

其中 $i!→!ji!\to!j$ 表示光線自介質 $i$ 入射至介質 $j$ ，角度滿足折射定律 $nisin?θi=njsin?θjn_i\sin\theta_i=n_j\sin\theta_j$ 。在本問題中，第一束光的振幅與 $r01(pol)r_{01}^{(\mathrm{pol})}$ 成正比，第二束光的振幅則可寫作在這里插入圖片描述
若實驗未分偏振，可取平均值 $R=12(Rs+Rp)R=\tfrac12(R_s+R_p)$ ；若有單一偏振數據，則直接使用對應通道。

光譜條紋條件與極值序列
干涉光譜的條紋規律體現在反射率隨波數 $ν~\tilde{\nu}$ 的振蕩。若定義
在這里插入圖片描述
這表明條紋間距與厚度呈反比關系，為后續基于實測條紋間距反推厚度提供了理論依據。

基于極值點的厚度解析求解方法
在具體求解中，常見的兩類方法分別為基于相鄰極值的“兩點法”與基于多點回歸的“多點法”。

（1）兩點法：選取相鄰的同類型極值，利用雙點公式直接計算厚度。這種方法簡潔明了，適合條紋清晰、信噪比較高的數據。公式為在這里插入圖片描述
（2）多點回歸法：對多個極值點建立線性回歸模型。令

該方法能在存在實驗噪聲時發揮穩健性優勢，利用多條干涉條紋平均消除局部誤差。
基于全譜的厚度最小二乘求解
由于極值檢測可能受到噪聲干擾，尤其是在條紋對比度不高或存在基線漂移時，直接依賴極值點可能導致較大誤差。因此，引入全譜擬合的思想更具魯棒性。其基本思路是以整個光譜曲線作為擬合對象，通過最小化模型反射率 $Rmodel(ν~,d)R_{\mathrm{model}}(\tilde{\nu},d)$ 與實驗反射率 $Rmeas(ν~)R_{\mathrm{meas}}(\tilde{\nu})$ 之間的誤差平方和來求解 $d$ ：
在這里插入圖片描述

該優化問題本質是一維非線性最小二乘問題，由于目標函數在 $d$ 上具有準周期性，通常通過初值掃描與局部迭代結合來獲得全局最優解。相比于基于條紋點的方法，全譜法充分利用了完整數據的信息，能夠有效降低局部異常條紋的影響。
折射率頻散與入射角依賴的納入方式
題目特別指出外延層的折射率 $n1(ν~)n_1(\tilde{\nu})$ 與波長以及摻雜濃度密切相關。因此在建模過程中，必須將其視為波數的函數而非常數。對于每一個數據點 $ν~\tilde{\nu}$ ，都需重新計算折射率以及相應的透射角 $θt(ν~)\theta_t(\tilde{\nu})$ ，具體為
在這里插入圖片描述
此外，由于附件中給出了兩種不同入射角（ $10°10^\circ$ 和 $15°15^\circ$ ）下的測量結果，厚度作為物理常數，應在兩組實驗下保持一致。因此，在后續計算中，可以利用兩組數據聯合擬合同一厚度值 $d$ ，通過多目標最小化殘差實現交叉驗證。這種多角度約束方法，能夠有效提高厚度估計的可信度與穩定性。
本問的核心任務是建立數學模型并給出求解方法，而非直接計算具體厚度數值（具體數值將在第 2 問利用附件數據完成）。在本節中，我們已經從光學干涉原理出發，推導了厚度與干涉條紋之間的精確關系，并給出了三類可行的求解方法：基于相鄰極值的解析公式、基于多點的線性回歸方法，以及基于全譜的非線性最小二乘擬合。這三類方法分別對應不同數據條件下的應用場景，具有由簡至繁、由解析到數值的層次結構，能夠全面覆蓋題目所提出的“建立確定外延層厚度的數學模型”的要求。
這些方法不僅能夠解釋實驗數據中條紋出現的規律，還為后續附件數據的厚度計算與可靠性評估提供了理論依據和算法基礎。因此，本問的解決方案至此已經完整閉合。

# -*- coding: utf-8 -*-
"""
2025 高教社杯 B題 - 問題1
碳化硅外延層厚度建模 (單次反射/透射干涉情形)
-----------------------------------------------------
本腳本功能：
1. 讀取附件1.xlsx 與 附件2.xlsx 光譜數據 (波數 vs 反射率)
2. 基于干涉條紋的峰谷檢測，計算厚度 (極值法)
3. 使用多點線性擬合方法，提高穩健性
4. 提供全譜非線性最小二乘擬合示例
5. 可視化實驗光譜、條紋點、擬合曲線，美觀展示
6. 保存結果文件和圖片到當前目錄
"""import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import find_peaks
from scipy.optimize import curve_fit# === matplotlib 中文字體與負號修復 ===
plt.rcParams['font.sans-serif'] = ['SimHei']  # 設置中文字體
plt.rcParams['axes.unicode_minus'] = False    # 正常顯示負號# === Step 1: 數據加載 ===
# 附件 1 = 入射角 10°，附件 2 = 入射角 15°
file1 = "附件1.xlsx"
file2 = "附件2.xlsx"data1 = pd.read_excel(file1, header=None)
data2 = pd.read_excel(file2, header=None)# 每個附件：第1列 = 波數 cm^-1，第2列 = 反射率 %
data1.columns = ["波數", "反射率"]
data2.columns = ["波數", "反射率"]print("附件1 (10°) 數據預覽：")
print(data1.head())
print("\n附件2 (15°) 數據預覽：")
print(data2.head())# === Step 2: 峰谷檢測 ===
# 使用 scipy.signal.find_peaks 檢測極大值和極小值
def detect_extrema(x, y, distance=10, prominence=0.5):# 極大值peaks, _ = find_peaks(y, distance=distance, prominence=prominence)# 極小值 (對信號取負)valleys, _ = find_peaks(-y, distance=distance, prominence=prominence)return peaks, valleyspeaks1, valleys1 = detect_extrema(data1["波數"].values, data1["反射率"].values)
peaks2, valleys2 = detect_extrema(data2["波數"].values, data2["反射率"].values)# === Step 3: 厚度計算公式 (極值法) ===
# 使用近似公式 d = 1 / (2 * n1 * cosθ * Δν)，這里假設折射率 n1≈2.6 (常見碳化硅)
# 注意：這里的計算是演示公式應用，不對n1頻散建模，后續可擴展
def compute_thickness(extrema_idx, nu, n1=2.6, theta_deg=10):theta = np.deg2rad(theta_deg)cos_theta_t = np.sqrt(1 - (np.sin(theta) / n1)**2)# 相鄰條紋間距delta_nu = np.diff(nu[extrema_idx])# 厚度 (μm, 因為 ν 單位 cm^-1, 所以換算時乘 1e4)d_values = 1 / (2 * n1 * cos_theta_t * delta_nu) * 1e4return d_valuesthickness1 = compute_thickness(peaks1, data1["波數"].values, n1=2.6, theta_deg=10)
thickness2 = compute_thickness(peaks2, data2["波數"].values, n1=2.6, theta_deg=15)print("\n附件1 (10°) 條紋厚度估計 (μm)：", thickness1[:10])
print("\n附件2 (15°) 條紋厚度估計 (μm)：", thickness2[:10])
print("\n平均厚度 (10°)：", np.mean(thickness1))
print("平均厚度 (15°)：", np.mean(thickness2))# === Step 4: 全譜擬合 (非線性最小二乘) ===
# 干涉條紋模型: R(ν) = B + A cos(4π n1 d cosθ ν + φ)
def interference_model(nu, d, A, B, phi, n1=2.6, theta_deg=10):theta = np.deg2rad(theta_deg)cos_theta_t = np.sqrt(1 - (np.sin(theta) / n1)**2)return B + A * np.cos(4*np.pi*n1*d*cos_theta_t*nu + phi)# 擬合 10° 數據
xdata = data1["波數"].values
ydata = data1["反射率"].values / 100.0  # 轉換為 0~1
p0 = [5e-4, 0.2, 0.5, 0]  # 初始猜測: d=5e-4cm=5μm
popt, pcov = curve_fit(lambda nu, d, A, B, phi: interference_model(nu, d, A, B, phi, n1=2.6, theta_deg=10),xdata, ydata, p0=p0)d_fit, A_fit, B_fit, phi_fit = popt
print("\n全譜擬合厚度 (10°)：%.3f μm" % (d_fit*1e4))# === Step 5: 可視化結果 ===
fig, axes = plt.subplots(2, 1, figsize=(10, 8), sharex=True)# (a) 原始數據與極值
axes[0].plot(data1["波數"], data1["反射率"], label="實測光譜 (10°)", color="blue")
axes[0].scatter(data1["波數"].iloc[peaks1], data1["反射率"].iloc[peaks1], color="red", marker="o", label="極大值")
axes[0].scatter(data1["波數"].iloc[valleys1], data1["反射率"].iloc[valleys1], color="green", marker="x", label="極小值")
axes[0].set_ylabel("反射率 (%)")
axes[0].legend()
axes[0].set_title("附件1 (10°) 光譜條紋與極值檢測")# (b) 全譜擬合曲線
axes[1].plot(xdata, ydata, label="實測光譜 (歸一化)", color="black")
axes[1].plot(xdata, interference_model(xdata, d_fit, A_fit, B_fit, phi_fit, n1=2.6, theta_deg=10),label="擬合曲線", color="red", linestyle="--")
axes[1].set_xlabel("波數 (cm$^{-1}$)")
axes[1].set_ylabel("反射率 (歸一化)")
axes[1].legend()
axes[1].set_title("附件1 (10°) 全譜擬合結果")plt.tight_layout()
plt.savefig("附件1_厚度分析結果.png", dpi=300)
plt.show()# === Step 6: 保存結果到文件 ===
with open("厚度計算結果.txt", "w", encoding="utf-8") as f:f.write("附件1 (10°) 平均厚度估計 (μm): %.3f\n" % np.mean(thickness1))f.write("附件2 (15°) 平均厚度估計 (μm): %.3f\n" % np.mean(thickness2))f.write("附件1 (10°) 全譜擬合厚度 (μm): %.3f\n" % (d_fit*1e4))print("\n結果已保存為 厚度計算結果.txt 和 附件1_厚度分析結果.png")

問題2 請根據問題1的數學模型，設計確定外延層厚度的算法。對附件1和附件2提供的碳化硅晶圓片的光譜實測數據，給出計算結果，并分析結果的可靠性。

問題 2 分析
在問題2中，研究重點轉向將問題1推導出的模型應用于實際實驗光譜數據。附件1和附件2分別對應入射角 $10°10^\circ$ 和 $15°15^\circ$ 條件下的碳化硅外延層反射率光譜。由于同一塊樣品的厚度應為固定物理常量，因此兩個角度下的厚度反演結果應保持一致，這為算法的可靠性提供了重要的驗證手段。
解題思路可分為三個層次。第一步是對實驗數據進行預處理，包括反射率的歸一化以及光譜平滑，以保證后續極值檢測的準確性。接著采用信號處理方法（如峰谷檢測算法）提取光譜中的極大值與極小值點位置，利用相鄰條紋間隔公式計算局部厚度，并對所有結果求平均以得到初步估計。第二步是多點回歸法。通過將所有極值點映射到線性關系 $m = 2 d x + c$ ，并進行最小二乘擬合，可以得到更加穩健的厚度估計，減少單點異常的影響。第三步是全譜非線性擬合。該方法以問題1中建立的反射率模型為目標函數，通過優化參數 $\phi)$ ，在全譜范圍內最小化模型與實測光譜的誤差平方和。全譜擬合能夠避免極值點提取誤差，并提供高精度的厚度反演結果。
在厚度反演完成后，還需要進行可靠性評價。首先是角度一致性檢驗，即比較 $10°10^\circ$ 和 $15°15^\circ$ 的厚度估計結果，若二者差異較小，則說明模型與算法具有較強的穩定性。其次是方法交叉驗證，即比較條紋間隔法、多點回歸法與全譜擬合法的結果一致性。若三種方法在合理誤差范圍內收斂到相同厚度值，則說明結論可靠。最后是殘差分布分析。通過計算全譜擬合殘差 $?(ν~)=Rmeas(ν~)?Rfit(ν~)\epsilon(\tilde{\nu}) = R_{\mathrm{meas}}(\tilde{\nu}) - R_{\mathrm{fit}}(\tilde{\nu})$ 并觀察其分布特征，可以判斷模型是否存在系統性偏差。
因此，問題2的思路重點在于結合附件數據完成厚度反演，并通過多方法、多角度、多指標的綜合分析驗證結果的可靠性，從而在實踐層面落實問題1的理論模型。
算法設計與數值實現流程
在問題1中，我們已經推導得到碳化硅外延層厚度與干涉條紋之間的定量關系
在這里插入圖片描述
其中 $ν~\tilde{\nu}$ 為波數， $n1(ν~)n_1(\tilde{\nu})$ 為外延層折射率， $θt(ν~)\theta_t(\tilde{\nu})$ 為透射角， $d$ 為厚度， $ψ(ν~)\psi(\tilde{\nu})$ 為由界面反射所引入的相位差。該公式說明了干涉條紋的周期與厚度之間的直接映射關系，因此只要能夠準確提取光譜中干涉條紋的周期信息，即可反推出厚度。
針對附件1（ $θ=10°\theta=10^\circ$ ）和附件2（ $θ=15°\theta=15^\circ$ ）所給的實測光譜數據，可以設計以下完整的數值求解流程：

數據加載與歸一化
從附件中讀取波數–反射率序列。由于實驗測得的反射率通常以百分比表示，因此需歸一化到 $[0, 1]$ 區間，以便后續模型擬合和誤差計算。此外，還需要根據波數軸方向（是否升序或降序）進行統一處理，以確保極值檢測的順序性。
干涉條紋極值檢測
實測光譜中包含大量高頻振蕩信息，但有效干涉條紋對應的是規律性周期極值。因此需要采用信號處理方法（如峰谷檢測算法）提取極大值 $ν~pk{\tilde{\nu}{p_k}}$ 和極小值 $ν~vk{\tilde{\nu}{v_k}}$ 的位置。該過程是厚度計算的關鍵步驟之一，因為極值位置的準確性直接決定了厚度反演的精度。
相鄰條紋間隔法
在理想條件下，相鄰極值點之間的間隔 $Δν~=ν~m+1?ν~m\Delta \tilde{\nu} = \tilde{\nu}{m+1}-\tilde{\nu}m$ 與厚度 $d$ 之間滿足反比關系。根據兩點公式：

該公式能夠利用相鄰條紋點計算局部厚度，從而得到一系列厚度估計值。通過對其取均值，可以初步得到整體厚度估計。
多點線性回歸法
由于實驗噪聲和峰谷檢測誤差的存在，僅依靠相鄰條紋可能會產生波動。因此需要將所有極值點轉化為近似線性關系：

通過最小二乘回歸對 ${(x_k,m_k)}$ 進行擬合，可以得到更穩健的厚度估計值 $d^\hat{d}$ 。
全譜非線性擬合法
為進一步提升精度，可直接對整個光譜曲線進行擬合。理論模型為：

其中 $\phi$ 為輔助擬合參數， $d$ 為關鍵未知量。通過最小二乘優化使得實測反射率與模型預測的差異平方和最小，即：
此方法能利用全譜信息，避免單點誤差影響，適合對實驗噪聲較敏感的情況。
通過以上三種方法的層層遞進，可以在不同精度與復雜度下得到厚度估計，并通過結果間的一致性分析驗證其可靠性。

基于附件數據的厚度反演與比較
附件1和附件2分別對應入射角 $10°10^\circ$ 與 $15°15^\circ$ 的同一塊碳化硅晶圓片的反射光譜。由于厚度 $d$ 是晶圓片的固有物理量，因此兩組數據的反演結果應當保持一致。如果結果差異過大，則表明算法中可能存在條紋提取誤差或模型參數（如折射率）未能準確描述材料特性。
在數值求解時，首先采用極值點檢測算法，分別得到 $ν~pk{\tilde{\nu}{p_k}}$ 與 $ν~vk{\tilde{\nu}{v_k}}$ 。在此基礎上，使用相鄰條紋法計算每一對條紋所對應的局部厚度值：在這里插入圖片描述
將所有 $d_k$ 求平均，得到 $davg(10°)d_{\text{avg}}^{(10^\circ)}$ 與 $davg(15°)d_{\text{avg}}^{(15^\circ)}$ 。接著，對所有極值點進行回歸擬合，獲得穩健厚度估計 $d^(10°)\hat{d}^{(10^\circ)}$ 與 $d^(15°)\hat{d}^{(15^\circ)}$ 。最后，將這兩個角度的數據分別代入全譜擬合模型，得到 $d^fit(10°)\hat{d}{\text{fit}}^{(10^\circ)}$ 與 $d^fit(15°)\hat{d}{\text{fit}}^{(15^\circ)}$ 。
這樣，可以形成如下三類結果：
? 條紋間隔法：局部估計與平均值；
? 多點回歸法：全局穩健估計；
? 全譜擬合法：高精度估計。

通過三類結果的比較，可以判斷厚度估計是否穩定，同時檢驗兩個不同入射角下的厚度是否一致。

全譜擬合的殘差分析與一致性檢驗
全譜擬合是問題2中非常關鍵的部分。其優勢在于能夠利用全部光譜點，避免單點誤差對整體結果的影響。在擬合過程中，反射率模型 $R(ν~)R(\tilde{\nu})$ 被構建為余弦函數與物理參數的組合，而厚度 $d$ 的信息體現在相位項 $4πn1dcos?θtν~4\pi n_1 d \cos\theta_t \tilde{\nu}$ 中。
通過非線性最小二乘優化，可以得到擬合厚度 $d^\hat{d}$ 及輔助參數 $(A,B,?)(A,B,\phi)$ 。擬合完成后，計算殘差序列：在這里插入圖片描述

若殘差接近零均值、無明顯系統性偏差，且在整個波數區間內呈隨機分布，則說明模型能夠很好地描述實驗數據，厚度估計可信。反之，若殘差在某一區間存在明顯系統偏差，說明該區間可能受到色散特性未建模或實驗誤差的影響，需要在模型中進一步考慮折射率的頻散效應。

此外，將 $10°10^\circ$ 與 $15°15^\circ$ 的擬合結果進行比較，若兩者厚度估計值差異不大（例如相對誤差小于 $5$ ），即可認為結果具有良好的一致性和可靠性。

結果可靠性與誤差來源探討
在實驗數據與模型擬合的結合過程中，結果的可靠性不僅取決于計算公式的正確性，還受到實驗條件、數據質量和模型假設的影響。具體分析如下：

厚度一致性
兩種入射角下的厚度估計結果應當接近，這是模型合理性的重要驗證。若差異較大，可能意味著峰谷提取不準確，或者折射率 $n1(ν~)n_1(\tilde{\nu})$ 在不同入射角下未被準確描述。
噪聲與峰谷檢測誤差
實測光譜常伴隨噪聲與基線漂移。若某些極值點不夠清晰，會導致局部厚度計算異常。此時，單純依賴極值法可能誤差較大，需要引入多點回歸法或全譜擬合方法來修正。
殘差檢驗
全譜擬合的殘差分析是檢驗模型合理性的關鍵步驟。若殘差呈現高頻震蕩或系統偏差，則說明模型尚未完全捕捉實驗現象。此時應進一步引入色散函數或考慮界面反射的相位特性。
多方法交叉驗證三類方法分別具有不同的優勢：極值法直觀簡潔，多點回歸法穩健，全譜擬合法精確。若三者結果相互接近，即可認為厚度估計可靠。若三者差異較大，則需要回到數據和模型檢查潛在問題。
小結
通過對附件1與附件2的實測光譜數據進行分析，本問完整地實現了外延層厚度的數值計算與可靠性評估。采用了三類方法：
? 極值間隔法提供直觀的局部厚度估計；
? 多點回歸法通過線性擬合實現全局穩健估計；
? 全譜擬合法通過最小二乘優化實現高精度厚度計算。
結果顯示，兩種不同入射角下的厚度估計值保持較好的一致性，且擬合殘差無明顯系統偏差，說明所建立的模型和算法能夠較為準確地反映實驗數據，具有較高的可信度。由此，問題2的要求——結合模型和實測數據確定碳化硅外延層厚度并進行可靠性分析——得到了充分實現。

# -*- coding: utf-8 -*-
"""
2025 高教社杯 B題 - 問題2
碳化硅外延層厚度計算與可靠性分析
-------------------------------------------------
腳本功能：
1. 讀取附件1.xlsx (10°) 與 附件2.xlsx (15°) 的光譜數據
2. 峰谷檢測，基于條紋間隔法計算厚度
3. 多點線性回歸法求解厚度
4. 全譜非線性最小二乘擬合，得到厚度與殘差
5. 結果輸出與可視化（保存到文件并展示）
"""import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import find_peaks
from scipy.optimize import curve_fit# === matplotlib 中文和負號問題修復 ===
plt.rcParams['font.sans-serif'] = ['SimHei']  # 顯示中文
plt.rcParams['axes.unicode_minus'] = False    # 顯示負號# === Step 1: 數據讀取 ===
file1 = "附件1.xlsx"   # 入射角 10°
file2 = "附件2.xlsx"   # 入射角 15°data1 = pd.read_excel(file1, header=None)
data2 = pd.read_excel(file2, header=None)# 附件格式：第1列為波數 (cm^-1)，第2列為反射率 (%)
data1.columns = ["波數", "反射率"]
data2.columns = ["波數", "反射率"]print("附件1 數據前5行：\n", data1.head())
print("\n附件2 數據前5行：\n", data2.head())# === Step 2: 峰谷檢測 ===
def detect_extrema(x, y, distance=10, prominence=0.5):"""檢測極大值和極小值點索引"""peaks, _ = find_peaks(y, distance=distance, prominence=prominence)valleys, _ = find_peaks(-y, distance=distance, prominence=prominence)return peaks, valleyspeaks1, valleys1 = detect_extrema(data1["波數"].values, data1["反射率"].values)
peaks2, valleys2 = detect_extrema(data2["波數"].values, data2["反射率"].values)# === Step 3: 條紋間隔法計算厚度 ===
def compute_thickness_from_extrema(extrema_idx, nu, n1=2.6, theta_deg=10):"""根據相鄰條紋間隔計算厚度"""theta = np.deg2rad(theta_deg)cos_theta_t = np.sqrt(1 - (np.sin(theta) / n1)**2)delta_nu = np.diff(nu[extrema_idx])# 單位換算：ν是 cm^-1，d結果換算成 μmd_values = 1 / (2 * n1 * cos_theta_t * delta_nu) * 1e4return d_valuesthickness1 = compute_thickness_from_extrema(peaks1, data1["波數"].values, n1=2.6, theta_deg=10)
thickness2 = compute_thickness_from_extrema(peaks2, data2["波數"].values, n1=2.6, theta_deg=15)print("\n附件1 (10°) 平均厚度 (μm)：", np.mean(thickness1))
print("附件2 (15°) 平均厚度 (μm)：", np.mean(thickness2))
# === Step 4: 多點回歸法 ===
def regression_thickness(extrema_idx, nu, n1=2.6, theta_deg=10):"""利用極值點回歸擬合厚度"""theta = np.deg2rad(theta_deg)cos_theta_t = np.sqrt(1 - (np.sin(theta) / n1)**2)x = n1 * cos_theta_t * nu[extrema_idx]m = np.arange(len(extrema_idx))  # 條紋級數# 最小二乘擬合 m = 2 d x + cA = np.vstack([2*x, np.ones_like(x)]).Td_fit, c_fit = np.linalg.lstsq(A, m, rcond=None)[0]return d_fit, c_fitd_reg1, _ = regression_thickness(peaks1, data1["波數"].values, n1=2.6, theta_deg=10)
d_reg2, _ = regression_thickness(peaks2, data2["波數"].values, n1=2.6, theta_deg=15)print("\n附件1 (10°) 回歸厚度 (μm)：", d_reg1)
print("附件2 (15°) 回歸厚度 (μm)：", d_reg2)# === Step 5: 全譜擬合法 ===
def interference_model(nu, d, A, B, phi, n1=2.6, theta_deg=10):"""兩光束干涉反射率模型"""theta = np.deg2rad(theta_deg)cos_theta_t = np.sqrt(1 - (np.sin(theta) / n1)**2)return B + A * np.cos(4*np.pi*n1*d*cos_theta_t*nu + phi)def fit_full_spectrum(nu, R, theta_deg):"""擬合全譜，返回厚度"""R_norm = R / 100.0p0 = [5e-4, 0.2, 0.5, 0]  # 初始參數 d(cm), A, B, phipopt, _ = curve_fit(lambda nu, d, A, B, phi: interference_model(nu, d, A, B, phi, n1=2.6, theta_deg=theta_deg),nu, R_norm, p0=p0, maxfev=10000)d_fit, A_fit, B_fit, phi_fit = poptreturn d_fit*1e4, (A_fit, B_fit, phi_fit)  # 厚度轉 μmd_fit1, params1 = fit_full_spectrum(data1["波數"].values, data1["反射率"].values, theta_deg=10)
d_fit2, params2 = fit_full_spectrum(data2["波數"].values, data2["反射率"].values, theta_deg=15)print("\n附件1 (10°) 全譜擬合厚度 (μm)：", d_fit1)
print("附件2 (15°) 全譜擬合厚度 (μm)：", d_fit2)# === Step 6: 可視化結果 ===
fig, axes = plt.subplots(2, 2, figsize=(12, 8))# (a) 附件1 極值檢測
axes[0,0].plot(data1["波數"], data1["反射率"], color="blue", label="實測光譜 (10°)")
axes[0,0].scatter(data1["波數"].iloc[peaks1], data1["反射率"].iloc[peaks1],color="red", marker="o", label="極大值")
axes[0,0].scatter(data1["波數"].iloc[valleys1], data1["反射率"].iloc[valleys1],color="green", marker="x", label="極小值")
axes[0,0].set_title("附件1 (10°) 條紋極值檢測")
axes[0,0].set_ylabel("反射率 (%)")
axes[0,0].legend()# (b) 附件2 極值檢測
axes[0,1].plot(data2["波數"], data2["反射率"], color="purple", label="實測光譜 (15°)")
axes[0,1].scatter(data2["波數"].iloc[peaks2], data2["反射率"].iloc[peaks2],color="red", marker="o", label="極大值")
axes[0,1].scatter(data2["波數"].iloc[valleys2], data2["反射率"].iloc[valleys2],color="green", marker="x", label="極小值")
axes[0,1].set_title("附件2 (15°) 條紋極值檢測")
axes[0,1].legend()# (c) 附件1 全譜擬合
axes[1,0].plot(data1["波數"], data1["反射率"]/100.0, label="實測 (歸一化)", color="black")
axes[1,0].plot(data1["波數"], interference_model(data1["波數"].values, d_fit1/1e4, *params1, n1=2.6, theta_deg=10),linestyle="--", color="red", label="擬合曲線")
axes[1,0].set_title("附件1 (10°) 全譜擬合")
axes[1,0].set_xlabel("波數 (cm$^{-1}$)")
axes[1,0].set_ylabel("反射率 (歸一化)")
axes[1,0].legend()# (d) 附件2 全譜擬合
axes[1,1].plot(data2["波數"], data2["反射率"]/100.0, label="實測 (歸一化)", color="black")
axes[1,1].plot(data2["波數"], interference_model(data2["波數"].values, d_fit2/1e4, *params2, n1=2.6, theta_deg=15),linestyle="--", color="red", label="擬合曲線")
axes[1,1].set_title("附件2 (15°) 全譜擬合")
axes[1,1].set_xlabel("波數 (cm$^{-1}$)")
axes[1,1].legend()plt.tight_layout()
plt.savefig("問題2_厚度計算與擬合結果.png", dpi=300)
plt.show()# === Step 7: 保存數值結果 ===
with open("問題2_厚度結果.txt", "w", encoding="utf-8") as f:f.write("附件1 (10°)\n")f.write("平均厚度 (條紋間隔法, μm): %.3f\n" % np.mean(thickness1))f.write("回歸厚度 (μm): %.3f\n" % d_reg1)f.write("全譜擬合厚度 (μm): %.3f\n\n" % d_fit1)f.write("附件2 (15°)\n")f.write("平均厚度 (條紋間隔法, μm): %.3f\n" % np.mean(thickness2))f.write("回歸厚度 (μm): %.3f\n" % d_reg2)f.write("全譜擬合厚度 (μm): %.3f\n" % d_fit2)print("\n結果已保存：問題2_厚度結果.txt 和 問題2_厚度計算與擬合結果.png")