電腦模擬與實驗驗證

在本章中，第一節利用時間序列分析法建立半導體混貨製程模型，第二節說明利 用ANOVA 模型分析截距並使用擴張卡曼濾波器作預測之流程，第三節分別利用擴張 卡曼濾波器、D-EWMA 與 time-varying D-EWMA 控制方法比對效能。

5-1 電腦模擬

5-1.1 製程干擾模型

為了驗證ANOVA 模型的效能，因此本文以(3.1)式作為基本模型。混貨製程模擬 [42]假設為兩個機台(1,2)，每個機台包含三個項目(1,2,3)，如(5.1)式所示：

( ) ( ) ( ) ( ) 1, 2 , ( ) 1, 2,3

i i i i i i j ki i i i

y k −b u k = +t p +

η

k i= j k = (5.1) 其中每一批次執行序的選擇由各項出現的機率決定，共四個執行序。各個項目出現的 機率如圖4。

圖 4 產品與機台選擇機率圖 P1：0.3

P3：0.4

T1：0.5

T2：0.5

P1T1：0.15 P2T1：0.15

P2T2：0.15 P3T2：0.20 P2：0.3

P1T2：0.15 P3T1：0.2

其中斜率 b 與製程的輸入 u 為常數，設定為：

1 2 1 2

[ , ] [1,1.3] , [ , ] [1,1]b b = u u = (5.2) 機台截距t 與產品截距_{i ( )}

j ki

p 設定如下圖 5 所示：

圖 5 產品與機台各項截距圖

並且假設機台製程干擾

η

_i( )k_i 滿足於ARIMA(0,1,1)。(3.2)式之 ARIMA(p,d,q)隨機程式 可改寫成下式：

0 1 1

(1−

φ

_iB )(1−B)

η

_i( ) (1k_i = −

θ

_iB ) ( )

ε

_i k_i (5.3)

其中

ε

_i( )k_i ∈N(0,0.2 )² 且[ , ] [0.5,0.2]

θ θ

_{1 2} = 。完成所有設定後，產生500 個批次的輸出 值作為歷史資料，並預測後100 個批次的輸出值。

5-1.2 ANOVA 模型

依上文所述，可得到混貨模擬製程之量測值。在得到量測值之後，依據第二章 所述之ANOVA 模型分析方法，可將(5.1)式表示如下：

( ) ( ) ( ) ( ) ( )

ij ki y ki i b u ki i i i pj ki i ki

α

= − = + +

μ τ

+

ξ

(5.4)

P1：10

P2：17

T1：5

T2：7

P1T1：11 P2T1：15

P1T2：13 P2T2：17 P1：6

P3T1：22

P3T2：24

製程原始截距如表1。ANOVA 分析出之

μ

、

τ

_i、 _{( )}

j ki

p 如表2、表 3 所示。

表 1 製程原始截距

分析截距 產品1 產品2 產品3

機台1 11 15 22

機台2 13 17 24

表 2 各因子之效果

A 因子效果(機台)：τ

τ

₁

τ

₂

τ

₃

-4.91129 -1.02495 5.93624

B 因子效果(產品)：p p₁ p ₂

-1.11029 1.11029

處理總平均值：μ 17.53468

表 3 ANOVA 分析截距

產品1 產品2 產品3

分析截距

分析 原始 誤差 分析 原始 誤差 分析 原始 誤差 機台1 11.51 11 4.6% 15.39 15 2.6% 22.36 22 1.6%

機台2 13.73 13 5.6% 17.62 17 3.6% 24.58 24 2.4%

在表 2 中，得知量測值經由 ANOVA 分析後，可得到各因子之效果。然後將之 整理成表3，其中，分析項為 ANOVA 分析後之截距，原始項為原始設定之截距，誤 差 項 為 分 析 截 距 與 原 始 截 距 的 平 均 絕 對 誤 差 率(mean absolute percentage error,

MAPE)。為了判定 ANOVA 分析後所得知截距是否為可信任的，故以平均絕對誤差率 作為判定。平均絕對誤差率表示如下：

MAPE=

1 n

i i

i i

D D D n

=

∑ −

(5.5)

其中

D

_i為截距實際值，

D

_i為截距分析值，n 為樣本數。

對於平均絕對誤差率值之準確率，Lewis [53]並提出表4的評估準則。

表 4 MAPE 評估標準 平均絕對誤差率(%) 準則說明

<10 高精確預測

10-20 優良預測

20-50 合理預測

>50 不準確預測

透過表 4 之評估準則，表 3 中各項截距之平均絕對誤差率皆小於 10%，因 此其 ANOVA 分析之截距乃是可以信任的。而在本文中，使用 ANOVA 分析出 之截距乃是一個初步的估計值，因此對於 ANOVA 的前提假設(常態性假設、

同質性假設、獨立性假設)，在此並不多加探討。

分析出模型之截距後，即可得(3.25)式，由(3.25)式可知其製程模型的製程干擾

i( )ki

ξ

。其製程干擾之數列如圖6、圖 7 所示。為了分析並建立製程干擾之數學模型，

依據第三章之建立時間序列模型方法，先將製程干擾歷史資料作一次差分後，由ACF 與 PACF 圖形以及計算各種不同階次組合之 AIC 值來確定製程干擾歷史資料之最佳 模型。數據處理過程如圖8 所示。

圖 6 機台 1 製程干擾數列圖(紅線為 ( )

η

_i k_i ，藍線為

ξ

_i( )k_i )

圖 7 機台 2 製程干擾數列圖(紅線為 ( )

η

_i k_i ，藍線為

ξ

_i( )k_i )

圖 8 擬合時間數列流程圖

圖 8 中令 ( )

ξ

_i k_i ，

ξ

_i(k_i− ， (1)

ξ

_i k_i− …為製程干擾數列，將其作一次差分，使數2)

列成為平穩數列，差分後的數列可表示為 W k ，_i( )_i W k_i( _i− ， (1) W k_{i i}− 2)

其中

W k_i( )_i =

ξ

_i( )k_i −

ξ

_i(k_i− = −1) (1 B) ( )

ξ

_i k_i

將差分數列W k 減去_i( )_i W k 之平均值_i( )_i

μ

_W 後，再利用時間序列擬合，即可得到製程干 擾時間數列之模型。

由上述所得之歷史資料擬合後的階次、係數與 AIC 值整理在表 6 表 7。表 5 中機

台1 製程干擾資料批次數為 236 個，表 6 中機台 2 製程干擾資料批次數為 264 個。

製程干擾歷史資料

ξ

_i( )k_i

i( )i

W k 一次差分

減去平均值

μ

_W

i( )i

W k

利用時間數列 擬合數列W k_i( )_i

表 5 機台 1 製程干擾數列擬合階次、係數與 AIC 值 機台1

φ

_{1 θ AIC}

(0,1,1) 0.48756 -65.47 (1,1,1) -0.06887 0.43281 -65.26 資料批次數：236

表 6 機台 2 製程干擾數列擬合階次、係數與 AIC 值

機台2

φ

₁

θ

AIC

(0,1,1) 0.21832 -48.44 (1,1,1) -0.42274 0.23993 -50.44 資料批次數：264

依據2-4 節所知，擬合後之數列已選擇 AIC 越接近 0 越為合適。由表 5、表 6 中 可知，機台1 擬合後之模式以選擇 ARIMA(1,1,1)較為合適，機台 2 擬合後之模式以 選擇ARIMA(0,1,1)較為合適。

5-1.3 預測

在半導體方面，已有一些學者[26]將 EWMA 控制器用來預測濺鍍製程，進而控

制沈積時間，使輸出厚度在目標範圍內。在處理製程漂移現象上D-EWMA 已被證明 比EWMA 控制器之效能來的好[54]，因此將 D-EWMA 作為比較性能的控制器之一，

並將D-EWMA 與擴張卡曼濾波器兩者做一比較，如表 7 所示。其中若製程模型為非 線性時變方程，則擴張卡曼濾波器可以適用；在權重的部份，卡曼增益功能類似於 D-EWMA 之權重，但擴張卡曼濾波器能夠隨製程狀態調整大小，D-EWMA 則否。另 一個最大特點為擴張卡曼濾波器是以狀態空間方程式衍生出之狀態估計方法，所以將 製程干擾與量測干擾都加入考慮，而D-EWMA 只考慮量測干擾，與實際製程情況有

所出入。除了D-EWMA 外，為了增加其可信度，另外再加上 time-varying D-EWMA 控制器作比較。time-varying D-EWMA 與 D-EWMA 的差異性在於 time-varying D-EWMA 在權重的選擇上多了一個折扣因子，藉由使用折扣因子以增強 D-EWMA 對參數的估計，使預測值較為快速的得到目標值，並且更快的達到控制的效果。

表 7 D-EWMA 與 Kalman filter 之比較

估計器類型 D-EWMA extended Kalman filter

適用範圍 線性方程式 線性或非線性方程式

參數類型 系統參數為固定值 系統參數可變動

gain 權值(weighting)為固定 隨製程而改變

誤差項 只有一個量測誤差項 考慮系統干擾以及量測誤差

參數估計 無 有

上節已擬合製程干擾數列

ξ

_i( )k_i 之時間序列，在此以機台1 為例。機台 1 之擬合 製程干擾數列

ξ

_i( )k_i 時間序列模型為ARIMA(1,1,1)，因此可以將(2.12)式改寫為下式：

1 1 1

1 1

(1−

φ

_iB )(1−B) ( ) (1

ξ

_i k_i = −

θ

_iB ) (

ε

_i k_i−1) (5.6) 假設製程干擾

ξ

_i( )k_i =X( ) { ( ),k_i = X k_{1 i} X k₂( )}_i ，並將(3.14)式與(3.15)式改寫成下二式：

X k₁( ) (1_i = +

φ ξ

₁) (_i k_i− +1) X (₂ k_i− +1)

ε

_i( )k_i (5.7)

2( )_{i i}1 _i( _i 1) _i1 _i( )_i

X k = −

φ ξ

k − −

θ ε

k (5.8)

因此可以將(5.6)式轉換成狀態空間形式：

¹

1 1

1 1 1

( ) ( 1) ( )

0

i

i i i i

i i

k

φ

k k

φ θ ε

⎡ + ⎤ ⎡ ⎤

= ⎢ ⎣ − ⎥ ⎦ − + ⎢ ⎣ − ⎥ ⎦

X X

(5.9)

ξ

ˆ ( ) [1 0] ( )

_i

k

_i

= X k

_i

+ v ( ) k

_i (5.10) 式中X( )k 為第_i i 個機台的第

k

批次製程干擾之狀態，ξ

ˆ ( )

_i

k

_i 為第 i 個機台的第

k

批次處

理過後之製程干擾，

ε

_i( )k_i 、v( )k_i 分別為系統與量測干擾，兩者皆屬於white-noise。

再將ARIMA 模型係數當作被估計狀態放入狀態方程式中，並依據擴張卡曼濾波器需 求，將狀態X( )k 擴充為_i { ( ),X k_{1 i} X k₂( ), ( ), ( )}_i

φ

_i1 k_i

θ

_i1 k_i ^T，再以4-2 節線性化方法可得到 線性化之狀態空間方程：

( )k_i = _i(k_i−1) +

ε

_i( ) k_i

X FX Γ (5.11)

ˆ ( )

_i

k

_i

X( ) k

_i

( ) k

_i

ξ

= H + v

(5.12)

其中

1 i1 1

2 1 1

1 1

( ) 1+ 1 ˆ (0) 0 ˆ

( ) - 0 (0) 0 X( ) , F=

( ) 0 0 1 0

( ) 0 0 0 1

i

i i

i

i i

i i

X k X

X k X

k k

k

φ φ φ

θ

⎡ ⎤

⎡ ⎤

⎢ ⎥

⎢ ⎥

⎢ − ⎥

⎢ ⎥

=⎢ ⎥ ⎢ ⎥

⎢ ⎥

⎢ ⎥

⎢ ⎥

⎣ ⎦ ⎣ ⎦

¹ 1

Γ , H=[1,0,0,0]

0 0

θ

i

⎡ ⎤

⎢− ⎥

⎢ ⎥

=⎢ ⎥

⎢ ⎥

⎣ ⎦

利用(5.12)式與擴張卡曼濾波器方法，即可計算下一批次的製程干擾，即預測形式可 寫為：

( ) ˆ

ˆ ( ) ˆ ( ) ( ) ( )

ij ki y ki i b u ki i i ti pj ki i ki

α

= − = + +

μ

+

ξ

(5.13)

為了測試擴張卡曼濾波器預測之效能，以兩個機台的製程干擾歷史資料為例，使 用D-EWMA 與 time-varying D-EWMA 兩個方法作為比對。其中 D-EWMA 控制器之 預測形式可寫為：

( ) ( ) ( ) _{( )} ( )

ij ki y ki i b u ki i i ti pj ki i ki

α

= − = + +

μ

+

ξ

(5.14)

更新製程干擾與干擾漂移項為：

1 1

[ ] [ ] (1 )( [ 1 1] [ 1 1])

i k ki i wi i ki wi i ki ki p ki i ki

ξ

=

ξ

+ −

ξ

− − + − − (5.15)

2 2

[ ] ( [ ] [ 1 1]) (1 ) [ 1 1])

i i i i i i i i i i i i i

p k k =w

ξ

k −

ξ

k − k − + −w p k − k − (5.16)

[ 1 ] [ ] [ ]

i ki ki i k ki i p k ki i i

ξ

+ =

ξ

+ (5.17)

其中為ξ 估計干擾，_i p 為干擾漂移量，_i w 為權重。 _i

其利用估計 k 時期製程之干擾與漂移量，將之當作為下一個時期(k+1)的製程干 擾與干擾漂移量。

而在time-varying D-EWMA 方面，其方法雖與 D-EWMA 方法相同，但在估計權 重時，使用一個簡單卻有效的時變性權重調整法，可以增強D-EWMA 對參數的估計，

其中time-varying D-EWMA 控制器之預測形式可寫為：

( ) ( ) ( ) ( ) ( )

ij ki y ki i b u ki i i ti pj ki i ki

α

= − = + +

μ

+

ξ

(5.18)

其調整模式為[55]：

1( ) max{ ,1 2}

i i i i

w k = w w (5.19)

2( ) min{ ,1 2} ^k

i i i i i

w k = w w +

γ

(5.20)

其中

γ

_i(0≤ < 表示為折扣因子。加入折扣因子，可以使初始暫態得到補償，快速

γ

_i 1) 的收斂至穩態。

time-varying D-EWMA 之設計是利用

γ

_i來加速干擾漂移量的收斂，亦即，當切換 執行序時，利用加大權重的方式來使D-EWMA 能夠快速地收斂至穩態，以降低在切 換執行序時所造成之暫態誤差；故在此time-varying 調整模式中，隨著批次量的增加，

2( )

i i

w k 會趨近於穩態之最佳權重。

在決定權重方面，選用均方誤差(mean square error，MSE)最小作為指標，其計算 方式為

MSE=

2 1

( ( ) ( ))

n

ij i ij i

i

k k

n

α α

=

∑

−

(5.21) 為了求得合適之權重，在得到製程歷史資料後，依下列設定求得w 與_i1 w ： _i2

1 0.1,0.2, ,0.9 wi =

2 0.1,0.2, ,0.9 wi =

並依(5.21)式，求得與歷史資料之 MSE 最小時之權重。如圖 9、圖 10 所示。

圖 9 各權重之 MSE 圖(機台一)

圖 10 各權重之 MSE 圖(機台二)

由圖9、圖 10 可以得知，機台一在w_i1=0.7,w_i2 =0.2時之 MSE 最小，機台二在

1 0.8, 2 0.1

i i

w = w = 時之MSE 最小。

然後以此求得之權重，加入

γ

_i做計算，再與歷史資料做比對，並以最小之 MSE 作為判定標準，藉此選定合適的

γ

_i。在圖 11、圖 12 中可以得知，機台一在

γ

_i = 時0 之MSE 最小，機台二在

γ

_i =0.83時之 MSE 最小。 此時所得之

γ

_i，為預測時折扣因 子之初始值，當開始預測後，

γ

_i^k隨著批次數的增加，會漸漸趨近於0，w_i2( )k 亦會趨_i 近於穩態之最佳權重。

圖 11 各折扣因子之 MSE 圖(機台一)

圖 12 各折扣因子之 MSE 圖(機台二)

擴張卡曼濾波器所用之參數取自表 5、表 6 中之參數，D-EWMA 之權重設定如 表8，time-varying D-EWMA 之權重設定如表 9。

表 8 D-EWMA 權重表

D-EWMA w i1 w _i2

機台1 0.7 0.2

機台2 0.8 0.1

表 9 time-varying D-EWMA 權重表

time-varying D-EWMA w i1 w _i2

γ

_i

機台1 0.7 0.2 0

機台2 0.8 0. 1 0.83

在比較預測效能方面，同樣選用MSE 作為比較方法。三種預測方法之各比較值 列於表10。機台 1 預測結果如圖 13，機台 2 預測結果如圖 14。機台 1 製程干擾預測 如圖15，機台 2 製程干擾預測如圖 16。

表 10 三種預測方法比較值

MSE (機台 1) MSE (機台 2) ARIMA(0,1,1)-EKF 0.047469 0.037871

D-EWMA 0.054652 0.040008 time-varying D-EWMA 0.054652 0.041150

圖 13 三種預測方法與原始資料之預測結果比較圖(機台一)

圖 14 三種預測方法與原始資料之預測結果比較圖(機台二)

圖 15 三種預測方法與原始資料之製程干擾預測比較圖(機台一)

圖 16 三種預測方法與原始資料之製程干擾預測比較圖(機台二)

由表 10 之各比較值可以得知，擴張卡曼濾波器之效能確實比 D-EWMA 和 time-varying D-EWMA 之效能還要好。

在上文例子中，機台製程干擾滿足於 ARIMA(0,1,1)。為了進一步驗證擴張卡曼 濾波器之效能，故加入機台製程干擾

η

_i( )k_i 滿足ARIMA(1,1,1)之例子作為比較。其中 加入[ , ] [0.3,0.4]

φ φ

_{1 2} = ，其他設定與上文相同。計算過程同上文所述。

ANOVA 分析截距列於表 11，機台製程干擾數列擬合階次、係數與 AIC 值列於 表 12、表 13，D-EWMA 與 time-varying D-EWMA 之權重列於表 14 和表 15。

ARIMA(1,1,1)三種預測方法之比較值列於表 16。機台 1 預測結果如圖 17，機台 2 預 測結果如圖18。機台 1 製程干擾預測如圖 19，機台 2 製程干擾預測如圖 20。

表 11 ARIMA(1,1,1)：ANOVA 分析截距

產品1 產品2 產品3

分析截距

分析 原始 誤差 分析 原始 誤差 分析 原始 誤差 機台1 11.92 11 8.3% 15.80 15 5.3% 22.76 22 3.4%

機台2 14.20 13 9.2% 18.08 17 6.4% 25.04 24 4.3%

表 12 ARIMA(1,1,1)：機台 1 製程干擾數列擬合階次、係數與 AIC 值 機台1

φ

_{1 θ AIC}

(0,1,1) 0.48421 -67.51 (1,1,1) -0.05225 0.44592 -64.46 資料批次數：236

表 13 ARIMA(1,1,1)：機台 2 製程干擾數列擬合階次、係數與 AIC 值 機台2

φ

1 θ AIC

(0,1,1) 0.21014 -50.22 (1,1,1) 0.17955 0.387 -47.33 資料批次數：264

表 14 ARIMA(1,1,1)：D-EWMA 權重表 D-EWMA w i1 w _i2

機台1 0.6 0.2

機台2 0.8 0.1

表 15 ARIMA(1,1,1)：time-varying D-EWMA 權重表 time-varying D-EWMA w i1 w _i2

γ

_i

機台1 0.6 0.2 0

機台2 0.8 0.43 0.33

表 16 ARIMA(1,1,1)：三種預測方法比較值 MSE (機台 1) MSE (機台 2) ARIMA(1,1,1)-EKF 0.049369 0.040039

D-EWMA 0.053389 0.040204 time-varying D-EWMA 0.053389 0.040364

圖 17 ARIMA(1,1,1)：三種預測方法與原始資料之預測結果比較圖(機台一)

圖 18 ARIMA(1,1,1)：三種預測方法與原始資料之預測結果比較圖(機台二)

圖 19 ARIMA(1,1,1)：三種預測方法與原始資料之製程干擾預測比較圖(機台一)

圖 20 ARIMA(1,1,1)：三種預測方法與原始資料之製程干擾預測比較圖(機台二)

5-2 實驗驗證

上節所述之內容，是以電腦模擬所得之結果，為了驗證本研究是否能在實際製程 中使用，因此以CMP 製程數據，作為驗證本研究之資料來源。隨著半導體製造技術 對高積集度與微縮化的不斷追求，在目前最尖端的半導體科技中，CMP 已經成為眾 所矚目的關鍵製程。不論在原始矽晶圓的拋光、IC 絕緣層及金屬層的形成、記憶體 的電容製作或新一代的鑲嵌法(damascene)及 SOI (silicon on insultor)的相關應用上，

CMP 不但是目前最直接也最有效的平坦化技術，也是突破下一世代半導體元件製造 瓶頸不可或缺的利器。

本文所分析的資料為某半導體廠CMP 製程的製程資料，機台為應材(applied materials)所製造，機型為 Reflexion，研磨對象為 300mm 晶圓，研磨層材質為 PE-TEOS oxide。機台有 4 個晶圓承載具(carrier)，同時有 3 個晶圓進行研磨，另一個承載具負 責填充或移出晶圓。每一待研磨晶圓必須通過3 個研磨平台進行研磨。研磨墊固定在 研磨平台上，和晶圓承載具以相反方向旋轉。

分析資料的第一步為測得晶圓的研磨率，晶圓研磨率的演算法如下式

PR = ( - ) /A A′ s (5.22)

其中PR 為研磨率，A 為研磨前量測晶圓的厚度，即前測值， A′為研磨後量測晶圓的 厚度，即後測值，s 為研磨時間。

本案例之製程為4 產品 1 機台，以研磨墊 1 的數據，共 69 個批次的研磨率作為 歷史資料，預測下一個研磨墊的研磨率(研磨墊 2)，然後重設初始值，再預測研磨墊 3 與研磨墊 4。研磨墊 2、3、4 共有 143 個批次的研磨率。研磨墊歷史資料如圖 21 所示。

在文檔中 中 華 大 學 (頁 33-60)