Kalman Filter For Sensor Signal Processing

Nghiên Cứu Bộ Lọc Kalman Áp Dụng Cho Bài Toán Cảm Biến Nguyễn Quốc Đính [email protected] Tp. Hồ Chí Minh, Tháng 8 năm 2008

Tóm tắt nội dung Trong một số ứng dụng thực tế, giá trị cảm biến đưa về chứa rất nhiều nhiễu, có thể bao gồm nhiễu của môi trường, của nguồn, của quá trình xử lý . . . , vấn đề đặt ra là làm sao loại bỏ được các nhiễu này để có được giá trị gần đúng với giá trị thật của biến cần đo nhất. Bộ lọc Kalman được giới thiệu ở đây như là một bộ lọc thích nghi thông thấp để ước lượng giá trị trung bình của giá trị, đồng thời sự thỏa hiệp giữa tốc độ hội tụ và sự giao động của giá trị ước lượng cũng được xem xét.

Mục lục 1

Giới Thiệu Chung Về Bộ Lọc Kalman Cho Hệ Tuyến Tính Rời Rạc

Mô hình của đối tượng tuyến tính rời rạc được biểu diễn thông qua hệ phương trình trạng thái sau: xk = Axk−1 + Buk−1 + wk−1

(1)

zk = Hxk + vk

(2)

Trong đó: - x: biến trạng thái. - u: biến đầu vào. - z: trạng thái đầu ra (đo lường được). - w, v: nhiễu quá trình và nhiễu đo lường tương ứng. Giả sứ các nhiễu quá trình w(t) và nhiễu đo lường v(t) là nhiễu trắng, có phương sai tương ứng là Q và R. p(w) ∼ N (0, Q)

(3)

p(v) ∼ N (0, R)

(4)

Với hệ thống như vậy, bộ lọc Kalman tuyến tính cho phép ta ước lượng một cách tốt nhất giá trị biến trạng thái xk , kí hiệu là x ˆk , sao cho hiệp phương sai của chúng Pk đạt giá trị nhỏ nhất. Pk = E[ek ek T ]

(5)

ek ≡ xk − x ˆk

(6)

1

Thuật toán Kalman cho hệ rời rạc như sau:

Như vậy ta có thể nhận thấy rằng thuật toán này gồm hai bước: bước ước lượng dự đoán (time update) và bước làm chính xác dự đoán này (measurement update) dựa trên thông tin đầu vào đo được zk . ˆk−1 Thuật toán này Chú ý rằng x ˆ− k là giá trị dự đoán được cập nhật từ giá trị ước lượng x được trình bày trên sơ đôi khối như sau:

Hình 1: Mô hình bộ lọc Kalman cho hệ tuyến tính

2

Bài Toán Cảm Biến

Cho bài toán như sau: có một thông số được đưa về từ cảm biến, tuy nhiên thông số này chịu tác động của nhiễu (nguồn có thể là từ nhiễu của quá trình, nhiễu đo lường, nhiễu do nguồn điện cung cấp, do ADC . . . ). Ta sẽ thiết kế bộ lọc Kalman để ước lượng được gía trị tốt nhất giá trị cần đo này. Mô hình của bài toán này: xk = Axk−1 + Buk−1 + wk = xk−1 + wk

(7)

Với giá trị thu được z ∈ < zk = Hxk + vk = xk + vk Với: - giá trị cần đo là x(t) 2

(8)

- tín hiệu thu được từ cảm biến là z(t) - nhiễu của quá trình là w(t). - nhiễu đo lường là v(t). Các phương trình cho bộ lọc Kalman: • Phương trình cho quá trình "Time Update": x ˆ− ˆk−1 k =x Pk− = Pk−1 + Q

(9)

• Phương trình cho quá trình "Measurement Update": Pk− Kk = − Pk + R x ˆk = x ˆ− ˆ− k + Kk (zk − x k) Pk = (1 − Kk )Pk−

(10)

Một số giả thiết đưa ra: Tín hiệu cần đo là x(t) ∼ (0, 50). Với bài toán này thì nhiễu hệ thống là bằng không, tuy nhiên để hệ thống đáp ứng được với sự thay đổi nhạy hơn, ta có thể cho nó có một phương sai nhỏ. Giả sử w(t) ∼ (0, 1e − 3) ↔ Q = 1e − 3. Với nhiễu của đo lường, trong bài toán mô phỏng dưới đây khi lấy giá trị từ cảm biến ta đã cho nó một nhiễu khoảng 10% giá trị thật. Việc chọn giá trị phương sai R(của giá trị nhiễu v(t)) ảnh hưởng đến tốc độ ước lượng của hệ thống. Từ phương trình tính hệ số Kalman Kk , ta thấy Kk tỉ lệ nghịch với R, giá trị R lớn thì tốc độ ước lượng chậm hơn, giá trị ước lượng có vẻ như ít tin tưởng hơn vào giá trị đo được, và đương nhiên kết quả ước lượng sẽ phẳng. Còn với giá trị R nhỏ thì ngược lại, tốc độ ước lượng nhanh hơn, bộ lọc tin tưởng hơn vào giá trị đo được, đương nhiên kết quả ước lượng ít phẳng hơn. Như vậy việc chọn thông số R khá quan trọng, có thể chọn R cứng bằng một giá trị cố định, hoặc là giá trị mềm có thể thay đổi, ví dụ như thay đổi theo tốc độ thay đổi thông số cần đo chẳng hạn. Chọn thông số ước lượng ban đầu: x(0) = 0, P0− = 5. Ta sẽ nhận thấy giá trị ban đầu x(0) là không quan trọng, vì nó sẽ tự thích nghi từng bước để đạt đến giá trị cần xác định như trong các mô phỏng bên dưới. Còn P0− thì chỉ cần khác không là được (nếu không thì giá trị ước lượng x(t) sẽ không tự thích nghi được!), và ta sẽ nhận thấy giá trị này nhỏ dần, tương ứng với giá trị ước lượng xấp xỉ với trị kì vọng. Chương trình mô phỏng được viết bằng Matlab được trình bày ở phần phụ lục A. Phân tích các kết quả thu đươc từ mô phỏng với những giá trị khác nhau của R Kết quả mong muốn là đường liền màu đen, giá trị từ cảm biến là các điểm màu xanh, Kết quả ước lượng là nét liền màu đỏ. Với R = 1:

3

Hình 2: Kết quả với R = 1 Nhận thấy tuy có tồn tại nhiễu khá lớn nhưng kết quả thu được từ bộ lọc gần bằng giá trị mong muồn không có nhiễu. Giá trị ban đầu x(0) tuy chọn khá xa so với giá trị thực, nhưng quá trình là thích nghi cho thấy nó tự biến đổi để đạt đến giá trị tốt nhất. Hiệp phương sai P giảm dần theo thời gian cho thấy kết quả ước lượng gần giống với giá trị thực của thông số cần đo. Sau khoảng 50 bước thì bắt đầu hội tụ,bằng 0.031, so với giá trị gán ban đầu là 1.

Hình 3: R = 1, đồ thị của Covariance Pk Pk giảm dần từ giá trị ban đầu là 1 về 0.031. Với R = 10: Bây giờ tăng giá trị R lên 10 lần, K = 10, thấy rằng bộ lọc chậm đáp ứng hơn. Tuy nhiên đường ước lượng phẳng hơn. Điều này cho thấy lựa chọn R lớn với những thông số có tốc độ thay đổi chậm. Giá trị P hội tụ bằng 0.095.

4

Hình 4: R = 10 Với R = 0.1: Bây giờ giảm R xuống 100 lần, K = 0.1, từ kết quả ta thấy kết quả ước lượng không được phẳng cho lắm, nhưng tốc độ ’học’ khá nhanh. R nhỏ thích hợp cho thông số có tốc độ thay đổi nhanh. Giá trị P hội tụ bằng 0.0095.

Hình 5: R = 0.01 Giá trị đầu vào thay đổi: Bây giờ ta không cho đầu vào là một hằng số cố định nữa, mà có sự thay đổi. Ví dụ dưới đây cho nó biến đổi theo hàm bậc nhất. Ta vẫn giữ nguyên mô hình như trên (mô hình ở trên mô tả cho hệ thống mà đầu vào là hằng số). Ta thấy với các giá trị R nhỏ thì kết quả ước lượng vẫn bám theo được giá trị mong muốn. Ví dụ với R = 0.1:

5

Hình 6: Với giá trị nhỏ R. Với R = 10: thì kết quả không tốt như vậy:

Hình 7: Với giá trị lớn R. Sở dĩ ở đây tôi không chọn lại một mô hình cho bài toán vì trong bài toán mà chúng ta cần giải quyết mức nhiên liệu thay đổi không theo một quy luật tuyến tính cố định nào cả. Nhận xét: • Bộ lọc Kalman cho thấy nó có khả năng giải quyết tốt với nhiễu cho mô hình của bài toán đã biết trước được mô hình. • Việc chọn thông số R là quan trọng, đòi hỏi chúng ta phải dung hòa giữa tốc đô hội tụ và độ phẳng của kết quả ước lượng.

6

Tài liệu [1] Leslie Lamport, An Introduction to the Kalman filter. Department of Computer Science, University of North Carolina at Chapel Hill, 2006 [2] http://en.wikipedia.org/wiki/Kalman_filter

7

Phụ Lục A Chương trình Matlab mô phỏng bộ lọc Kalman function kalman( number) %input for i=1:number input(i) = 10; end %process noise and measurement noise w = 1e-5; %process noise, may it equal to Zero v = 1e0; %mesurement noise %equivelent covariance Q = 1e-3; R = 1; %initial value xpre = 0; Ppre =1;

%chosing this value is not importance. %must different from Zero

%to plot the result time = []; % timing pos = []; % measurement result, from sensor posest = []; % measurement estimate, after Kalman filter Pest = []; %loop for i=1:number %making input and output x = input(i) + w*randn; z = x + v*randn; %measurement update (correct) K = Ppre* inv(Ppre + R); xpre = xpre + K*(z - xpre); P = (1 - K)*Ppre; %update for ploting time = [time i]; pos = [pos z]; posest = [posest xpre]; Pest = [Pest P]; %time update (predict for next step) xpre = xpre; Ppre = P + Q; end clf plot(time, pos, ’*’) hold on 8

plot(time,posest,’LineWidth’,3,’MarkerEdgeColor’,’r’, ’MarkerFaceColor’,’g’,’MarkerSize’,20,’Color’,’r’); plot(time, input,’LineWidth’,2,’MarkerEdgeColor’,’r’, ’MarkerFaceColor’,’g’,’MarkerSize’,20,’Color’,’k’); grid on

9