光學薄膜橢圓偏振數據反演的優化算法

童晟飛，王正憶，陳 星

（浙江大學 理學部 物理系，浙江 杭州310027）

1 引 言

材料薄膜的厚度、光學常量是膜系設計、材料研究和薄膜制備中不可缺少的參量．膜厚和光學常量不能由橢圓偏振儀直接給出，后期需要經過一系列的擬合計算．橢圓偏振法由于無須測定光強的絕對值，因而具有很高的精度和靈敏度．對超薄膜（可達1 nm以下）同樣具有很高的靈敏度．如果橢圓偏振儀的起偏器和檢偏器回轉角有0．01°的精度時，則膜厚0．007 nm的變化就可能被檢測出來．由于橢偏參量確立的方程是超越方程，無法直接從測量數據通過計算得到薄膜參量的解析解［1］，需要進行反演計算，可通過優化算法做數值反演解．本軟件運用C＋＋語言通過變步長迭代循環的優化算法開發數據處理程序，能夠直接根據輸入的數據得到結果，無需進行查表或作圖等人工操作，并且適應性較廣（可以接受任意的激光波長和基底材料），速度較快（平均運行時間2 s），精度較高（折射率精確到0．001，厚度精確到0．01 nm）．程序可以接受1組測量數據或者激光入射角不同的2組測量數據，輸入1組數據時程序會輸出所有符合方程的解（包括跨周期的解），要得到唯一的解，需要改變橢偏儀的入射光角度后再進行1次實驗，2種情況下重疊的解即為最終結果．對于有吸收的薄膜，即折射率為復數，用了結合模擬退火法的量子行為粒子群混合算法進行計算，QPSO算法包含比較當前粒子的適應值與其歷史最優值以不斷刷新歷史最優值的過程迭代更新，系統就能夠獲得擺脫局部能量極小點的機會，并找到一個更好的、更接近于整體的極小點［2］．

2 基本原理

若有一偏振光在樣品表面反射，如圖1，可以將其分解成為在2個互相垂直方向上的分量波：振動面平行入射面的線偏振光稱p波，振動面垂直入射面的線偏振光稱s波．

圖1 光波在單層薄膜中的傳播

在界面1和界面2菲涅耳反射系數分別為［3－5］：

其中，n0為空氣折射率（n0＝1．0），n為膜折射率，ns為襯底折射率（可以是復數折射率），θ0為入射角，r1p和r1s是界面1上p光和s光的反射系數，r2p和r2s是界面2上p光和s光的反射系數．薄膜的反射率可由菲涅耳反射系數表示為

其中，

為使計算方便，定義ψ和Δ使之滿足下列關系［2］：

所以理論上只要確定ψ和Δ后，就可解出n和d．

3 薄膜無吸收情況的處理

當薄膜無吸收時，比較簡單的求解程思路是：將n從0．01到5．00，d從1 nm到500 nm的數值逐一代入方程，求出相應的ψ和Δ，然后與輸入值進行比對，當誤差精度小于指定的范圍時，就將這一組n和d輸出．程序使用的是固定步長，利用固定步長計算將會導致程序運行時間過長，計算量浪費．

作出（ψi，Δi）關于（n，d）變化的圖，找到圖上表現為尖銳的峰，說明（ψi，Δi）隨（n，d）的變化較快．可見，對于絕大多數n和d，對應的ψ和Δ都與輸入值相去甚遠，只有在接近方程組的解時，ψ和Δ才會接近輸入值［4］．

所以本程序使用可變步長來進行n和d的循環．根據對大量實際計算的總結，發現n的改變對ψ和Δ的影響很大，而d的改變產生的影響相對較小，因此對d使用可變步長．

3．1 第1步循環

具體方法是：n從0．01到5．00以步長0．01循環，d從1到501以5為步長進行循環，當｜δψ｜＜1且｜δΔ｜＜3時，以1為步長對d－4到d＋4范圍內的d進行逐一計算．計算完成后恢復5為步長．

建立存放10組解的數組，存放n和d，把n，d代入方程作正解．計算出的ψ和Δ的誤差精度ε＝｜δψ｜2＋｜δΔ｜2．每計算1組n和d，將ε與數組中10組解對應項分別進行比較．如果計算得到的解的偏差小于已存放的偏差最大的解，則替換偏差最大的這組解．這樣，在完成循環時，存放在數組中的解就是δψ與δΔ相對最小的解，即誤差精度ε為最小．

程序算法的流程圖如圖2所示．

圖2 第1步循環

3．2 第2步循環

在進行第1次循環之后，得到了10組近似的解．這些解的精度還是不夠的，需要進行更精確的計算．實際上，第1次循環得到的10組解并不是孤立的，這些解的折射率n值通常都集中在1個數值附近．通過編寫1個函數可以將這些解歸納起來，得到1組近似解，并使這組解中的厚度d為相近的n值對應的所有d值中最小的．

對這組解分別再進行1次循環，其中n的循環從n－0．05開始，到n＋0．05終止，步長為0．001；對應的d的循環從d－2到d＋2，步長為0．01．建立1個數組，只存放1組解，對解的篩選方法和前面相似，只是這次只保留偏差最小的解．這樣的循環可以使n精確到小數點后3位，d精確到小數點后2位，從而達到精確度的要求．

程序算法的程序圖如圖3所示．

圖3 第2步循環

3．3 厚度周期的計算

厚度d具有周期性，這是由（5）～（7）式中的e指數決定的．e－2iδ這一項的δ以π為周期，所以d的周期為

3．4 最終結果

第2次循環得到的2組解分別加上整數個周期即為方程最終的解．

從上面的分析可以看到，僅通過1組實驗數據得到的解是不唯一的，要得到唯一的解，需要改變橢偏儀的入射光角度后再進行1次實驗，2種情況下重疊的解即為最終結果．

要測量膜厚過1個周期的真實厚度，可采用改變入射角的方法．從表1可以看出：樣品3和5薄膜在2個不同的入射角測量下，2組測量厚度基本相等，薄膜的厚度相差0．580 nm，0．260 nm，相應的周期厚度相差4．750 nm，1．180 nm，薄膜的真實厚度應該在1個周期厚度以內．從表1看出4號樣品薄膜在2個不同的入射角測量下厚度相差7．19 nm，相應的周期厚度相差6．69 nm．因為膜的厚度是唯一的，所以一定疊加周期厚度．對于重合解的判斷，需要控制好精度．如果精度過高可能會找不到重合解，而精度過低則容易導致周期判斷錯誤．一般地，2次入射角相差越大則厚度周期相差得也越大，此時可以降低重合的精度；而入射角相差較小時厚度周期也相差得比較小，此時就必須使用較高的重合精度．

本程序設置了相應的參量．在本程序中，利用2組數據分別得到的折射率n相差不能超過入射角之差（使用角度作為單位時的數值）的0．001倍（最小不低于0．005），否則視為誤差過大．對于2組數據的每個厚度周期，判斷厚度d相差是否超過入射角之差的0．1倍（最小不低于0．5 nm），如果超過，則認為薄膜不是位于這一周期，繼續搜索其他的周期厚度，直到找到重合的為止．

表1 薄膜樣品的計算結果比較數據

3．5 計算結果的檢驗

用上述算法開發的軟件對在單晶硅襯底上經過熱氧化形成的SiO2／Si薄膜樣品1和2進行分析［3－4，6］，對不同的數據處理方法進行比較．選定相關參量襯底如下：樣品1，2，4單晶硅襯底（n＝3．85－0．02i）． 樣 品 1 和 2 入 射 光 波 波 長635 nm．樣品4入射光波波長為632．8 nm．樣品3和5為K9玻璃（n＝1．515）上的TiO2薄膜．入射光波波長為632．8 nm．要測量膜厚過1個周期的真實厚度，可采用改變入射角的方法．樣品采用上述兩步迭代循環逼近計算，可得材料折射率n和厚度d的值．樣品4的最后結果界面如圖4所示，整個計算過程大約持續了2 s．表1為本軟件與國外儀器軟件及國內其他軟件［4，6－8］的比較．本文算法計算的折射率和厚度已分別精確到第3位和第2位的穩定解，表明只要實驗測得的（ψ，Δ）足夠精確，就能夠得到足夠精確的（n，d）．

圖4 樣品4的最后結果界面

因為鍍完膜的樣品折射率和厚度在理論上是唯一的，所以選擇2個相差不大的入射角分別計算應該得到同樣的折射率和厚度結果，由此可以從多周期的解中確定出薄膜的真實折射率和厚度，以及樣品是否跨周期的判斷．

4 薄膜有吸收情況的處理

對光有吸收的薄膜，除了折射率n和厚度d之外，還存在消光系數k，一共有3個未知數要求解．如果通過一定的實驗方法能夠確定這3個未知數中的1個，則可以通過前面的變步長循環迭代的算法進行求解．

如果3個未知數都不能預先確定，對于前面的變步長循環法，雖然已經有所優化，但用于處理這種情況仍然不適應，為此需要使用更加優化的算法．常見的優化算法有遺傳算法、粒子群算法、模擬退火算法等．

4．1 粒子群算法

粒子群優化（particle swarm optimization）算法是廣泛應用的優化算法，基本思想是，將每組（n，k，d）視為三維空間中的1個點，然后一群粒子隨機地在該空間中運動，各粒子之間獨立運動的同時保持著一定的聯系．當1個粒子接近方程解所對應的點時，其他的粒子就會向這個粒子的方向靠近，從而使所有粒子達到收斂點．

4．2 量子行為粒子群算法和模擬退火算法

在傳統的粒子群算法中，粒子按照軌道運動，而有限的軌道無法覆蓋解空間中的每個點．因此傳統PSO算法不能保證收斂到全局最優點．同時，如果1個粒子找到了較好的點，其他的粒子也會朝這個方向運動，一旦粒子都到達這個較好的點附近，就可能無法發現其他更好的點，導致早熟收斂．

為了避免上述問題，人們提出了量子行為粒子群優化（quantum behaved particle swarm optimization，QPSO）算法．

在QPSO中，粒子不再有軌道的概念，而是按照量子力學的原理在解空間中運動．在這一模型中，通過適應度函數判斷粒子在某個位置的優越程度．每個粒子每次迭代時記錄它的歷史最優位置pbest，并與其他所有粒子的pbest進行比較得到種群的歷史最優位置gbest．設p點為pbest與gbest之間的任意一點，將p點設為一個delta勢阱，利用量子力學原理計算粒子出現在某點的概率，通過蒙特卡羅模擬的方式來測量粒子的位置，從而對粒子進行刷新．

刷新粒子的具體方法是［7］：

1）對粒子的每一維，在pid和pgd之間得到1個隨機點，

2）計算所有粒子的平均最好位置mbest，

3）計算評價參量

其中1／g也稱為創造系數，可以用β表示．4）刷新位置，

U（0，1）屬于（0，1）上區間的隨機數．

QPSO雖然能夠保證算法的全局收斂性，但是在收斂的情況下，由于所有的粒子都向最優解的方向飛去，導致粒子的多樣性損失，使得后期收斂速度明顯變慢，容易陷入局部最優．

為了解決這個問題，需要考慮與其他優化算法結合使用．模擬退火算法（simulated annealing，SA）是一種全局搜索能力極強的算法．該方法起源于金屬冷卻退火這一自然現象．從能量的角度分析，對于一個在溫度為T的處于熱力學平衡的體系，其在某一微觀狀態內能為E時的概率分布服從玻耳茲曼分布；系統按照概率f分布于所有不同的能量狀態中，即使在很低的溫度下，系統也可能處于較高的能量狀態，因此，相應地系統就能夠獲得擺脫局部能量極小點的機會，并找到一個更好的、更接近于整體的極小點［1，9］．

將評價函數設置為系統的“能量”，即為模擬退火算法，該方法具有極強的全局搜尋能力，可以有效避免陷入局部最優的問題［8］．

在SA中，從一個狀態過渡到另一個狀態的概率為

其中－Δf在本問題中為適應度之差．當溫度從初始溫度逐漸降低到零點時，可以認為系統已經到達最優狀態，求解完成．

4．3 結合模擬退火法的量子行為粒子群的混合算法

SA的全局搜索能力很強，但收斂速度偏低．相比之下，QPSO在收斂速度方面更有優勢．為了更好利用2種算法各自的優點，可以將SA算法混入QPSO算法中．

QPSO算法包含比較當前粒子的適應值與其歷史最優值以不斷刷新歷史最優值的過程．這個選擇過程會指導其他粒子搜索適應值相對好粒子的附近空間，這一過程中低適應值的粒子將被拋棄．然而，那些低適應值的粒子可能具備潛在的更好的進化趨勢，這樣就降低了整個種群的進化．為了解決這一問題，將模擬退火的概念引入這一過程，低適應值的粒子將以一定概率被接受．這樣具有較好進化趨勢的低適應值粒子在搜索空間會以一定的概率繼續飛行，可以有效地避免陷入局部最小值［7，10］．

由此形成了一種混合算法，流程圖如圖5所示．在該算法中共有3個參量需要調節：創造系數β，退火初始溫度T0，退火常量CR．這些參量決定了整個收斂的準確度以及收斂速度的快慢，需要根據求解的問題來不斷嘗試．在本問題的求解中，設定β＝1．2，T0＝50，CR＝0．98．

圖5 SA－QPSO算法流程圖

4．4 計算結果的檢驗

將一些模擬計算來驗證這種算法用于橢偏反演計算的可行性及有效性．在玻璃襯底上的單層TiO2模型及SiO2／Si模型薄面樣品進行分析，TiO2模型的參量值為：n＝2．625，k＝0．015，d＝417 nm，ns＝1．523，ks＝ 0，入射激光波長λ＝370 nm，從這些數據可以計算得到（ψ，Δ）在入射角分別為45°，60°，70°，80°時的準確值，具體數據見表2［1，7－9，11－14］．SiO2／Si 模 型 參 量 值 為：n＝1．46，k＝0，d＝108 nm，ns＝3．85，ks＝ －0．02．計算在入射角分別為45°，60°，70°時的準確值，具體數據見表3［15］．

可根據實驗數據（ψi，Δi）作出（ψi，Δi）關于（n，d）變化的圖，找到圖上表現為尖銳的峰，估算k，n，d值的范圍．求解時設定消光系數搜索范圍為［0，0．2］，折射率范圍［0，4］，厚度范圍［30，500］．

利用混合算法編寫的程序數據進行處理，每次選取2組數據進行計算，得到的結果如表2～3所示．可以看出，程序具有較高的求解精度，求解誤差通常能達到10－7的水平，表明對于精確的測量值，程序可以給出比較準確的計算．

表2 TiO2模型3組數據計算結果比較

注：n′，k′，d′為其他軟件求解值［7－8，12－13］．

表3 SiO2／Si模型2組數據計算結果比較

從表2～3可以看出：對于有吸收薄膜樣品用結合模擬退火法的量子行為粒子群混合算法，其反演得到薄膜光學常數計算結果精確．

5 結 論

本文提出變步長循環的優化算法對沒有吸收薄膜樣品計算速度快（2 s左右），準確性高（n精確到0．001，d精確到0．01 nm）．程序計算過程中沒有使用隨機數，因此同一組數據，每次計算的結果都相同．利用2組數據確定了多周期真實厚度．對于有吸收薄膜樣品本文提出了利用加入模擬退火的量子行為粒子群優化算法．本程序改善粒子群的全局、局部搜索能力和收斂速度，是具有很強全局搜索能力的算法，可以在保證計算速度的情況下減少早熟收斂的發生．利用量子力學原理計算粒子出現在某一點的概率，通過蒙特卡羅模擬的方式來測量粒子的位置，從而對粒子進行刷新，可很快跳出局部極小收斂達到全局最優點．

［1］ 廖清君．單波長消光橢偏儀的數據處理研究［D］．成都：四川大學，2002．

［2］ 劉細成．透射光譜法測量薄膜參數的研究［D］．成都：四川大學，2003．

［3］ 張麗娟．透射光譜法測量光柵參數的研究與全抗反射膜的設計［D］．成都：四川大學，2005．

［4］ 陳星，童晟飛，王正憶，等．橢圓偏振儀測量薄膜折射率及周期厚度解的分析［J］．實驗技術與管理，2011，28（6）：42－46．

［5］ Azzam R M A，Bashara N M．Ellipsometry and poliarized light［M］．Amsterdam：North－Holland Publishing Co．，1977：269．

［6］ 王洪濤．橢圓偏振法測量薄膜參量的數據處理［J］．物理實驗，2001，21（7）：8－12，17．

［7］ Law B M，Pak H K．Ellipsometric imaging of surface drops［J］．J．Opt．Soc．Am．A，1996，13（2）：379－384．

［8］ Dobrowolski J A，Ho F C，Waldorf A．Determination of optical constants of thin film coating materials base on inverse synthesis［J］．Applied Optics，1983，22（20）：3191－3200．

［9］ 廖清君，王植恒，王磊，等．模擬退火法在吸收薄膜的橢偏反演算法中的應用［J］．光學學報，2002，22（6）：683－687．

［10］ 劉靜．粒子群優化算法研究及其在優化理論中的應用［D］．無錫：江南大學，2007．

［11］ 王黨社，張建科，徐均琪．薄膜光學常數的粒子群算法［J］．計算物理，2008，25（2）：208－212．

［12］ Ross T，Cormier G．Particle swarm optimization for ellipsometric data inversion of samples having an arbitrary number of layers［J］．J．Opt．Soc．Am．A，2010，27（2）：319－326．

［13］ Cormier G，Boudreau R．Genetic algorithm for ellipsometric data inversion of absorbing layers［J］．J．Opt．Soc．Am．A，2000，17（1）：129－134．

［14］ Comfort J C，Urban F K，Barton D．An algorithm for analyzing ellipsometric data taken with multiple angles of incidence［J］．Thin Solid Films，1996，290－291：51－56．

［15］ 周建華，游佰強，洪志哲．反射式橢偏測試系統的全局優化算法［J］．激光與紅外，2005，35（3）：214－216．