Kalman Filter 簡介

Kalman Filter 最早是 1960 年 Kalman [1]所提出，其利用狀態空間方程式 (State Space Equation)之特性所發展出的遞迴估計系統狀態方法。整體而言 卡曼濾波器具有以下之特點[23]：

1.卡曼濾波器之解是一個適合在電腦上計算的遞回方程式。

2.對於數據可以逐一的即時處理，並基於目前與前一刻計算出的狀 態估計值運用遞迴方式算出下一個狀態的估計值。

3.系統方程可以是時變的動態系統，輸入之訊號或數據之型態可以 是非平穩的。

4.可應用於各種線性(Linear)或非線性(Non-Linear)問題。

以下將介紹狀態空間方程式及 Extended Kalman Filter(EKF)之架構。

4-1 狀態空間方程式簡介

一個線性離散時間系統之狀態空間方程式可由狀態方程式與輸出 方程式之組合來表示，即[23]

(k + =1) (k +1, ) ( )k k + (k +1, ) ( )k k

X F X Γ w (4.1) (k + =1) (k +1) (k + +1) (k+1)

Y H X v (4.2) (4.1)式為狀態方程式，代表系統中狀態的變化，(4.2)式為輸出方程式，

代表由系統狀態產生變化後系統的輸出。其中 X 代表系統狀態向量

(n× ，Y 為系統輸出向量(1) m× ，F 為 k 到 k+1 時期之狀態轉移矩陣1) (n n× ，Γ 為干擾輸入矩陣() n n× ，H 為觀測矩陣() m n× ；w 與 v 分別) 為系統干擾向量 (n× 以及量測干擾向量(1) m× ，並且為 White-Noise，1) 且與系統狀態X 三者互相獨立，並且有以下性質：

[ ( ) ^T( )] ( , ) 0, , 0

E w k w j =Q

δ

k j Q≥ ∀k j≥ (4.3) [ ( ) ^T( )] ( , ), 0, , 0

E v k v j =R

δ

k j R≥ ∀k j≥ (4.4) [ ( ) ^T( 1)] 0, , 0

E w k v j+ = ∀k j≥ (4.5) [ (0) ^T( )] 0, [ (0) ^T( 1)] 0 0

E X w k = E X v k + = ∀ ≥k (4.6) 則與(4.1)式與(4.2)式相對應之卡曼濾波器形式為

ˆ(k + =1) ˆ(k +1 )k + (k +1)[ (k + −1) ˆ(k +1 )]k

X X K Y Y (4.7) ˆ(k +1 )k = (k +1, ) ( )k ˆ k

X F X (4.8) ˆ(k +1 )k = (k +1) (ˆ k +1 )k

Y H X (4.9)

K ( k + = 1) P ( k + 1 ) k H

^T

( k + 1)[ ( H k + 1) ( P k + 1 ) k H

^T

( k + + 1) R ( k + 1)]

⁻¹

(4.10)

( k + 1 ) k = ( k + 1, ) ( ) k k

^T

( k + 1, ) k + ( k + 1, ) ( ) k k

^T

( k + 1, ) k

P F P F Γ Q Γ

(4.11)

(k + = −1) [ (k+1) (k +1)] (k +1 )k

P I K H P (4.12) 初始值設定為

ˆ (0) = E [ (0)] = (0), (0) = var (0) ∀ ≥ k 0

X X X P X

(4.13) 式中K 為卡曼增益(Kalman Gain)，遞迴計算流程如下圖 4，利用已給定

之初始值預測下一個時點狀態變異矩陣P(k+1|k)，接著計算卡曼增益，

利用這些已知數據去預測下一個時點(k+1)之系統狀態 ˆ (X k +1 )k 及輸出 觀測值 ˆ (Y k+1 )k 。此時新的觀測值 Y(k+1)進入後，就可以更新目前在 k+1 時點的狀態估計向量

X ˆ ( k + 1)

以及狀態變異矩陣 (P k+1)。如此可計 算下一個狀態變異矩陣P(k+2|k+1)之預測值，依次類推。

圖 4 Kalman Filter 運算流程圖[23]

4-2 Extended Kalman Filter

Extended Kalman Filter 是將 Kalman Filter 估算法更進一步的延伸，

其可將系統未知的參數視為狀態並加以估計，或是處理非線性製程狀態 估計。對於目前濺鍍製程，因沈積率的歷史資料並不完整，故利用 EKF

初始值 P(0)

X

(0) k=0

計算

P ( k + 1 ) k

計算

K k ( + 1)

計算

X ^ˆ ( k + 1 ), ( k Y ^ˆ k + 1 ) k

計算

X ˆ ( k + 1)

計算

P ( k + 1)

k=1,2,…n Q(k)

R(k+1)

Y(k+1)

來估計沈積率的狀態，並將系統參數加入狀態估計。一個狀態空間隨機 離散非線性系統可表示為[23]：

(k + =1) f[ ( ), ( ), ] k w k k ∀ ≥k 0

X X (4.14)

(k + =1) h[ (k +1), (v k+1),k+1] ∀ ≥k 0

Y X (4.15)

X 為系統狀態向量 (n× ，1)

f

與

h

為非線性向量函數，Y 為系統觀測向量

(m× ；w、v 與其他之假設皆與上一節相同。所以(4.14)式與(4.15)式在1) 每一個時期作泰勒展開後，取其一次項之線性化狀態空間方程可寫為：

(k + =1) ( ) ( )k k + ( ) ( )k w k

X F X Γ (4.16)

(k + =1) (k+1) (k + +1) v k( +1)

Y H X (4.17)

其中

( ) ˆ( ), ( ) 0

[ ( ), ( ), ]

( ) ( ) _{X k X k w k}

f k w k k

k k _{= =}

= ∂

∂ F X

X (4.18)

( 1) ˆ( 1 ), ( 1) 0

[ ( 1), ( 1), 1]

( 1)

( 1)

_{X k X k k v k}

h k v k k

k k

_{+ = + + =}

∂ + + +

+ = ∂ +

H X

X

^(4.19)

( ) ˆ( ), ( ) 0

[ ( ), ( ), ]

( ) ( )

_{X k X k w k}

f k w k k

k w k

_{= =}

= ∂

∂

Γ X

(4.20)

則相對應之 Extend Kalman Filter 遞迴方程為：

ˆ(k + =1) f[ ( ), ]ˆ k k + (k +1)[ (k + −1) h[ (ˆ k +1 ),k k+1]]

X X K Y X (4.21)

( k + = 1) ( k + 1 ) k

^T

( k + 1)[ ( k + 1) ( k + 1 ) k

^T

( k + + 1) R k ( + 1)]

⁻1

K P H H P H

(4.22)

( k + 1 ) k = ( ) ( ) k k

^T

( ) k + ( ) ( ) k k

^T

( ) k

P F P F Γ Q Γ

(4.23)

+ = − + + +

ˆ (0) 0, (0) var (0) = =

X P X

(4.25)

其計算流程與 4-1 節相同

4-3 Non-fixed Interval Extended Kalman Filter

事實上，目前量產產品之濺鍍製程中要做到等間距量測是一件困難 的工程，廠商無不想盡辦法提昇工作機台的利用率，為了要使控制器更 能符合實際之情況，就必須把非等間距量測之因素考慮進去。本論文中 是將未量測到之沈積率批次資料，利用預測 d 個時期（d 為未量測之批 次數）之後的沈積率預測值與新量測而得的沈積率資料來更新模型的狀 態及參數，以此類推，求得最佳的沈積率預測值，進而換算出更精確的 沈積時間。以下將說明預測方法及流程。

以 ARIMA(1,1,0)模型為例子，考慮製程中，製程干擾未必屬於 White-Noise，故提高 MA 階次，消除製程干擾之相關性，所以模型轉 變為 ARIMA(1,1,1)，即

( ) (1 ) ( 1) ( 2) ( ) ( 1)

X k = +

φ

X k− −

φ

X k− +a k −

θ

a k− (4.26)

( ) ( ) ( )

Y k = X k + v k

(4.27) 其中X(k)為沈積率狀態，Y(k)為輸出沈積率，a(k)與v(k)分別為濺鍍製程 干擾與量測干擾，ψ與θ為製程模型參數，由(4.26)式與(4.27)式可轉換 為狀態空間方程形式，並將模型參數加入估計。

( ) k = ( k − + 1) a k ( )

X FX Γ

(4.28)

( ) ( ) ( )

Y k = HX k + v k

(4.29) 其中

1

2

( ) ( ) ( )

( ) ( ) X k k X k

k k φ θ

⎡ ⎤

⎢ ⎥

⎢ ⎥

= ⎢ ⎥

⎢ ⎥

⎣ ⎦

X

_，

1 1

1 1 ˆ (0) 0 0 ˆ (0) 0 0 0 1 0 0 0 0 1

X X

φ

φ

⎡ + ⎤

⎢ ⎥

− −

⎢ ⎥

= ⎢ ⎥

⎢ ⎥

⎢ ⎥

⎣ ⎦

F ，

1

0 0 θ

⎡ ⎤

⎢ − ⎥

⎢ ⎥

= ⎢ ⎥

⎢ ⎥

⎣ ⎦

Γ

，

[

^{1 0 0 0}

]

=

H

將上述狀態方程利用卡曼濾波方程預測靶材之沈積率，對於非等間距問 題，處理方法如下：假設兩個非等間距量測點之間(k, k+d)經過了 d 個 批次，利用在第 k 批次時預測第 k+d 批次之沈積率狀態 ˆ (X k +d k)與狀 態變異

P ( k + d k )

並計算在(k+d)批次預測之沈積率 ˆ (Y k +d k)，待新的 沈積率量測值得到後，更新沈積率狀態，計算流程如圖 5 所示，其餘計 算方式皆與 4-2 節相同。

圖 5 Non-Fixed Interval EKF 計算流程 初始值 P(0)

X

(0) k=0

計算

P ( k + 1 ) k

計算

K ( k + 1)

計算

X ˆ ( k + d )

計算

P ( k + d )

k=1,2,…n

Q(k)

R(k+1)

Y(k+d)

計算

X ^ˆ ( k + d k ), ( P k + d k ), ( Y ^ˆ k + d k )

d 為未量測之 lot 數

預測 d 個 lot 以後之 狀態

X 與狀態誤差 P

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