童晟飛，王正憶，陳 星

（浙江大學(xué) 理學(xué)部 物理系，浙江 杭州310027）

1 引 言

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

2 基本原理

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

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

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

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

其中，

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

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

3 薄膜無吸收情況的處理

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

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

所以本程序使用可變步長來進(jìn)行n和d的循環(huán)．根據(jù)對大量實際計算的總結(jié)，發(fā)現(xiàn)n的改變對ψ和Δ的影響很大，而d的改變產(chǎn)生的影響相對較小，因此對d使用可變步長．

3．1 第1步循環(huán)

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

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

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

圖2 第1步循環(huán)

3．2 第2步循環(huán)

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

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

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

圖3 第2步循環(huán)

3．3 厚度周期的計算

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

3．4 最終結(jié)果

第2次循環(huán)得到的2組解分別加上整數(shù)個周期即為方程最終的解．

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

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

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

表1 薄膜樣品的計算結(jié)果比較數(shù)據(jù)

3．5 計算結(jié)果的檢驗

用上述算法開發(fā)的軟件對在單晶硅襯底上經(jīng)過熱氧化形成的SiO2／Si薄膜樣品1和2進(jìn)行分析［3－4，6］，對不同的數(shù)據(jù)處理方法進(jìn)行比較．選定相關(guān)參量襯底如下：樣品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èi)肷浣堑姆椒ǎ畼悠凡捎蒙鲜鰞刹降h(huán)逼近計算，可得材料折射率n和厚度d的值．樣品4的最后結(jié)果界面如圖4所示，整個計算過程大約持續(xù)了2 s．表1為本軟件與國外儀器軟件及國內(nèi)其他軟件［4，6－8］的比較．本文算法計算的折射率和厚度已分別精確到第3位和第2位的穩(wěn)定解，表明只要實驗測得的（ψ，Δ）足夠精確，就能夠得到足夠精確的（n，d）．

圖4 樣品4的最后結(jié)果界面

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

4 薄膜有吸收情況的處理

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

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

4．1 粒子群算法

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

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

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

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

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

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

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

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

3）計算評價參量

其中1／g也稱為創(chuàng)造系數(shù)，可以用β表示．4）刷新位置，

U（0，1）屬于（0，1）上區(qū)間的隨機(jī)數(shù)．

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

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

將評價函數(shù)設(shè)置為系統(tǒng)的“能量”，即為模擬退火算法，該方法具有極強(qiáng)的全局搜尋能力，可以有效避免陷入局部最優(yōu)的問題［8］．

在SA中，從一個狀態(tài)過渡到另一個狀態(tài)的概率為

其中－Δf在本問題中為適應(yīng)度之差．當(dāng)溫度從初始溫度逐漸降低到零點時，可以認(rèn)為系統(tǒng)已經(jīng)到達(dá)最優(yōu)狀態(tài)，求解完成．

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

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

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

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

圖5 SA－QPSO算法流程圖

4．4 計算結(jié)果的檢驗

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

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

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

表2 TiO2模型3組數(shù)據(jù)計算結(jié)果比較

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

表3 SiO2／Si模型2組數(shù)據(jù)計算結(jié)果比較

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

5 結(jié) 論

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

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

［2］ 劉細(xì)成．透射光譜法測量薄膜參數(shù)的研究［D］．成都：四川大學(xué)，2003．

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

［4］ 陳星，童晟飛，王正憶，等．橢圓偏振儀測量薄膜折射率及周期厚度解的分析［J］．實驗技術(shù)與管理，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］ 王洪濤．橢圓偏振法測量薄膜參量的數(shù)據(jù)處理［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］ 廖清君，王植恒，王磊，等．模擬退火法在吸收薄膜的橢偏反演算法中的應(yīng)用［J］．光學(xué)學(xué)報，2002，22（6）：683－687．

［10］ 劉靜．粒子群優(yōu)化算法研究及其在優(yōu)化理論中的應(yīng)用［D］．無錫：江南大學(xué)，2007．

［11］ 王黨社，張建科，徐均琪．薄膜光學(xué)常數(shù)的粒子群算法［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］ 周建華，游佰強(qiáng)，洪志哲．反射式橢偏測試系統(tǒng)的全局優(yōu)化算法［J］．激光與紅外，2005，35（3）：214－216．