Đề tài Ước lượng sai số mô hình trong bộ lọc Kalman bằng phương pháp lực nhiễu động

Tạp chí Khoa học ĐHQGHN, Khoa học Tự nhiên và Công nghệ 26, Số 3S (2010) 310‐316 310 _______ Ước lượng sai số mô hình trong bộ lọc Kalman bằng phương pháp lực nhiễu động Kiều Quốc Chánh* Khoa Khí tượng Thủy văn và Hải dương học, Trường Đại học Khoa học Tự nhiên, ĐHQGHN, 334 Nguyễn Trãi, Hà Nội, Việt Nam Nhận ngày 11 tháng 8 năm 2010 Tóm tắt. Trong bài báo này, một phương pháp xác định sai số mô hình trong bộ lọc đồng hóa Kalman sẽ được trình bày. Kiểm định phương pháp này trên mô hình Lorenz 40 biến chỉ ra rằng phương pháp mới có nhiều ưu điểm so với phương pháp tăng cấp nhân đơn thuần. Mở rộng của phương pháp này cho các hệ với bậc tự do lớn như trong các mô hình dự báo thời tiết nghiệp vụ cũng sẽ được thảo luận. Từ khóa: Đồng hóa số liệu, lọc Kalman, mô hình Lorenz, mô hình dự báo số 1. Mở đầu∗ Đồng hoá số liệu về bản chất là một quá trình trong đó số liệu quan trắc và một trường phỏng đoán nền được kết hợp với nhau một cách thống kê để thu được điều kiện ban đầu tốí ưu cho mô hình số (trong bài này thuật ngữ ‘mô hình’ ngụ ý một phương trình biểu diễn dưới dạng sai phân dùng để giải một bài toán phương trình đạo hàm riêng một cách xấp xỉ với điều kiện biên và điều kiện ban đầu cho trước). Đặc trưng thống kê của bài toán đồng hóa số liệu chính là cốt lõi của tất các thuật toán đồng hoá hiện tại. Nếu mô hình và các dữ liệu quan sát là hoàn hảo, bài toán đồng hóa số liệu khi đó sẽ đơn thuần chỉ là một bài toán nội suy (hay ngoại suy) tối ưu nhiều chiều. Nếu quan trắc là tuyệt đối nhưng mô hình ẩn chứa các sai số nội tại, bài toán đồng hóa sẽ không còn là một bài toán nội suy tối ưu đơn thuần vì khi đó điều kiện ban đầu chính xác sẽ không còn luôn được trông đợi (thậm chí ngay cả khi phép nội suy là chính xác) do các dữ liệu quan sát có thể tiềm ẩn các thành phần không cân bằng mà mô hình không cho phép tích phân. Ví dụ các sóng trọng trường có thể được kích thích và lan truyền rất nhanh, dẫn đến sự phá huỷ tính ổn định của mô hình. Nếu cả mô hình và quan trắc là không hoàn hảo thì rõ ràng sự bất định này phải được tính đến trong mô hình một cách thích hợp. Trong nghiên cứu này, chúng tôi sẽ tập trung chủ yếu vào sai số nội tại của mô hình, gọi tắt là sai số mô hình. Vấn đề sai số của dữ liệu quan trắc thiên về bài toán kiểm định chất lương quan trắc nghiệp vụ và sẽ không được xem xét ở đây. ∗ ĐT: 84-4-38584943. E-mail: [email protected] Trong thực tế, ước lượng sai số mô hình là một vấn đề rất khó của bài toán đồng hoá số liệu do nguồn lớn nhất của sai số mô hình lại chính là các quá trình vật lí không được hiểu biết đầy đủ. K.Q. Chánh / Tạp chí Khoa học ĐHQGHN, Khoa học Tự nhiên và Công nghệ 26, Số 3S (2010) 310‐316 311 Ví dụ như các lực rối, lực ma sát hay tham số hoá không đầy đủ các quá trình vật l ý vi mô của mô hình. Thêm vào đó, các xấp xỉ số học của mô hình cũng có thể là một nguồn sai số đáng kể của do các thuật toán sai phân hữu hạn có thể chứa các nghiệm phi vật lí hoặc trở nên mất ổn định khi vi phạm các điều kiện tích phân. Các kĩ thuật xử lí sai số mô hình trong các thuật toán đồng hóa số liệu hiện đại bao gồm kĩ thuật tăng cấp nhân [1], kĩ thuật tăng cấp cộng tính [2], hay phương pháp hiệu chỉnh sai số hệ thống [3]. Một sự giới thiệu tổng quan đầy đủ về các kĩ thuật xử lí sai số mô hình có thể được tìm thấy trong nghiên cứu [4]. Trong nghiên cứu này, một phương pháp khác dựa trên giả thiết rằng nguồn của sai số mô hình chủ yếu là do sự biểu diễn không đầy đủ của các quá trình vật lí sẽ được trình bày, tạm gọi là phương pháp lực nhiễu động. Trong phần tiếp theo, sự thiết lập cơ sở lý thuyết của phương pháp lực nhiễu động sẽ được thảo luận. Phần 3 mô tả các ứng dụng của phương pháp này đối với mô hình Lorenz 40- biến. Sự mở rộng của phương pháp này cho một hệ với bậc tự do lớn hơn như là mô hình dự báo thời tiết số sẽ được xem xét trong phần 4, và một vài kết luận sẽ được đưa ra trong phần cuối cùng. 2. Cơ sở lí thuyết Xem xét một phương trình mô tả sự tiến triển của một trạng thái x có dạng tổng quát như sau: )()( tM dt d Fxx += (1) trong đó x(t) ∈ ℜn là một vector trạng thái n- chiều phụ thuộc vào thời gian có phân bố xác suất ban đầu đặc trưng bởi ma trận hiệp biến Pf, M là một mô hình phi tuyến mô tả sự tiến triển của trạng thái, và F(t) ∈ ℜn là vector lực1. Giả thiết một tập hợp số liệu quan sát ∈ ℜ _______ 1 Để đơn giản các ký hiệu, không gian Eulerian với metric đơn giản sẽ được ngầm hiểu sao cho x và các phép tính o iy p p- chiều được cho trước tại các thời điểm gián đoạn {ti∈I }với một phân bố xác suất đặc trưng bởi ma trận sai số hiệp biến R. Khi đó, lọc Kalman toàn phần cho phép đồng hoá tập số liệu quan sát này sẽ được cho bởi các phương trình dưới đây (xem [5]) o iy i T i a ii f i QLPLP += −−− 111 (2) 1 11 )( − −− += RHHPHPK TfiTfi ))(( fi o i f i a i xHyKxx −+= f i a i 1)( −−= PKHIP trong đó và là ma trận sai số hiệp biến nền (hay dự báo) tại thời điểm i-1 và i, L là mô hình tiếp tuyến của mô hình M, và là ma trận sai số hiệp biến phân tích tại thời điểm i-1 và i, K là ma trận trọng số, Q f i 1−P f iP a i 1−P a iP i là ma trận sai số mô hình, và H là toán tử biến đổi từ không gian mô hình sang không gian quan trắc. Lọc Kalman sẽ được áp dụng tại từng thời điểm i cho mỗi chu trình đồng hóa và sau đó được tích phân tiếp theo đến thời điểm thứ i+1 tại đó quá trình phân tích với bộ lọc Kalman lại được lặp lại. Như được thảo luận ở trong phần giới thiệu, hai nguồn sai số chính của mô hình đặc trưng bởi ma trận Q là các xấp xỉ số học của phương trình (1) và các lực cưỡng bức không được hiểu biết đầy đủ F. Mặc dù loại sai số đầu tiên liên quan đến thuật toán tích phân mô hình có thể được khắc phục bằng cách thiết kế các thuật toán tính toán hợp lí, loại sai số thứ hai liên quan đến tính chất vật lý rất khó kiểm soát và có đóng góp lớn nhất đến sai số mô hình tổng cộng, đặc biệt trong các hệ phức tạp như là hệ thống khí quyển-đại dương. Để bài toán được thiết lập một cách tường minh, giả thiết rằng các thuật toán sai phân hữu hạn của phương trình (1) là đủ chính xác sao cho sai số của mô hình do các xấp xỉ số học có thể được tạm bỏ qua và chúng ta do đó có thể tập trung hoàn vào vector có thể được thực hiện với topo tương ứng. Các ký tự in đậm ngụ ý các vector trong các không gian mô hình hay không gian quan trắc một cách tương ứng. K.Q. Chánh / Tạp chí Khoa học ĐHQGHN, Khoa học Tự nhiên và Công nghệ 26, Số 3S (2010) 310‐316 312 loại sai số mô hình vật lí. Giả thiết rằng lực F(t) là một biến ngẫu nhiên với một phân bố xác suất cho trước, nhiệm vụ của chúng ta bây giờ là tìm một biểu diễn cho ma trận sai số mô hình Q với giả thiết này. Từ phương trình mô hình (1), dạng biến phân của nó có dạng )()()( t dt d FxxJx δδδ += (3) trong đó gradient J(x) được định nghĩa bởi ∂M/∂x. Trong trường hợp tổng quát nghiệm chính xác của phương trình (3) là không khả tích, vì với δF(t) phụ thuộc tường minh vào thời gian, nghiệm chính xác sẽ liên quan đến việc thừa số hoá các ma trận không khả nghịch. Tuy nhiên, chú ý rằng mặc dù F(t) phụ thuộc vào thời gian, phân bố thống kê của nhiễu lực lại có thể được giả thiết là không phụ thuộc vào thời gian với một phân bố xác xuất có biên độ cho trước, nghĩa là δF không phụ thuộc vào thời gian. Trong trường hợp này, nghiệm thuần nhất của phương trình (3) sẽ có dạng 0321321 212111 ...))(()(()(( )(()(()((( 21 111 1 111 xxJxJxJ xJxJxJIx δ δ ++ ++= ∫∫∫ ∫∫∫ −−− −−− t t t t t t t t t t t t h i i i i i i i i i i tttdtdtdt ttdtdttdt (4) trong đó δxo là nhiễu động ban đầu do điều kiện ban đầu không chính xác tại t = 0. Với nghiệm thuần nhất (4), nghiệm cuối cùng của phương trình (3) sẽ được cho bởi FxJxx δδδ ))(()()( 1 ttt h −−= (5) Nghiệm này có thể được viết ngắn gọn hơn bằng việc đưa vào toán tử sắp xếp thời gian T, được định nghĩa như là [6] ∑ −−= − σ σσσσ σσθθ ))(())...1(()()...()}()...({ )()1()2()1(1 nttttttT nnn HHHH (6) trong đó tổng σ chạy trên tất cả các giao hoán của (1....n), (tổng cộng có n! các giao hoán) và hàm Heaviside được định nghĩa bởi ⎩⎨ ⎧ < ≥=− ji ji ji tt tt tt 0 1 )(θ Một cách thực chất, toán tử sắp xếp thời gian sẽ sắp xếp lại tất cả các ma trận sao cho các ma trận với thời gian trễ nhất sẽ đứng ở phía ngoải cùng bên trái. Đây là một kĩ thuật rất quen thuộc trong bài toán lí thuyết trường lượng tử [6]. Với toán tử sắp xếp thời gian T, nghiệm (5) có thể được viết lại một cách cô đọng như sau FJxxJx δδττδ 10}])(([exp{)( 1 −−= ∫ − t ti dTt (7) trong đó hàm mũ của ma trận được định nghĩa như là ∑∞ = = 0 ! )exp( n n n AA . (8) Nghiệm (7) có thể được kiểm tra một cách dễ dàng bằng cách thay nó trực tiếp vào phương trình (3). Với nghiệm (7), mô hình tiếp tuyến L được định nghĩa trong lọc Kalman sẽ có dạng }])(([exp{ 1 ∫ − = t ti dT ττ xJL (9) và sai số mô hình bây giờ sẽ được cho bởi }))({( 11 TE FJFJQ δδ −−≡ (10) Để tính toán sai số mô hình Q chú ý rằng nếu chúng ta có một tập mẫu n phép thử với cùng một điều kiện ban đầu sao cho δxo = 0, rõ ràng là khi đó từ phương trình (7) tất cả các sai số sẽ được tạo ra chỉ bởi lực nhiễu động δF, nghĩa là δx = δxh. Như vậy, chúng ta có thể thu được ma trận Q theo hai cách khác nhau. 1. Tính toán trực tiếp ma trận Q bằng cách thống kê các vector (J-1δF). Điều này được thực hiện bằng phương pháp lấy mẫu n phép thử δF để tạo ra một mẫu n các vector (J-1δF). Từ đó, ma trận sai số mô hình Q có thể thu được một cách dễ dàng từ phương trình (10). 2. Cách tiếp cận thứ hai là thực hiện n phép tích phân mô hình với các lực F được làm nhiễu một cách ngẫu nhiên. Các tích phân mô hình K.Q. Chánh / Tạp chí Khoa học ĐHQGHN, Khoa học Tự nhiên và Công nghệ 26, Số 3S (2010) 310‐316 313 tuy nhiên sẽ được tiến hành với cùng một điều kiện ban đầu sao cho tất cả các sai số mô hình trong các đầu ra có thể được gán cho các lực bị làm nhiễu. Đầu ra của các phép chạy này bây giờ có thể được tính toán thống kê để thu được ma trận Q. Phương pháp tiếp cận thứ hai sẽ được chọn trong nghiên cứu này bởi vì nó có thể được mở rộng một cách dễ dàng đối với các mô hình nghiệp vụ hoặc trong các trường hợp tổng quát hơn, ví dụ như các điều kiện biên bị làm nhiễu như được thảo luận trong phần 4. 3. Thiết kế thí nghiệm Để xem xét một cách đầy đủ nhất có thể hiệu quả của thuật toán lực nhiễu động, mô hình Lorenz 40 biến sẽ được sử dụng như là một mô hình mẫu trong nghiên cứu này sao cho lọc Kalman toàn phần có thể được sử dụng. Cùng với điểm nổi bật của việc sử dụng lọc Kalman toàn phần, mô hình này có thể được tích phân một cách rất chính xác bằng việc sử dụng thuật toán Runge-Kutta bậc 4. Điều này sẽ làm tối thiểu hoá sai số mô hình do các phương pháp tính toán số và do đó cho phép xem xét một cách đầy đủ phương pháp lực nhiễu động. Để so sánh phương pháp mới với các cách tiếp cận khác, kĩ thuật thừa số tăng cấp nhân sẽ được thực hiện song song với phương pháp lực nhiễu động được trình bày trong phần 2. Một sự so sánh đầy đủ hơn với các kỹ thuật xử lý khác bao gồm tăng cấp cộng tính hay khử sai số hệ thống sẽ được đề cập đến trong các nghiên cứu tiếp theo. 3.1. Mô hình Mô hình Lorenz 40-biến được cho bởi (xem [7]) FxMFxxxx dt dx iiiii i +≡+−−= −+− )()( 211 (11) Trong đó F được lấy chính xác F = 8.0 cho trạng thái thực. Điều kiện ban đầu của trạng thái thực sẽ được chọn một cách ngẫu nhiên. Một khi đã chọn, điều kiện ban đầu này tuy nhiên sẽ được giữ không đổi trong tất cả các thí nghiệm tiếp theo. Trạng thái thực ở trên sẽ được tích phân 1000 bước thời gian và in ra tại từng bước. Điều kiện biên tuần hoàn cho xi được áp dụng tại i = 0 sao cho x0 = xN. Mô hình này có sự tiến triển hỗn loạn sau một thời gian chuyển tiếp khoảng 50 bước tích phân. Bước thời gian δt = 0.01 sẽ được sử dụng trong tất cả các thí nghiệm. 3.2. Mô hình tiếp tuyến Mô hình tiếp tuyến cho mô hình Lorenz có thể thu được trực tiếp từ phương trình (9). Do khối lượng tính toán lớn và các sai số làm tròn của tích phân ma trận, chúng tôi chỉ giới hạn các tính toán tại các xấp xỉ bậc một và bậc hai của phương trình (9). Với xấp xỉ bậc một, chúng ta thu được dạng quen thuộc ∑ = +≈ M i itxt 1 ))((JIL δ (12) ở đó M là số bước tích phân mà tại đó quan trắc sẽ được đồng hóa. Ví dụ, M = 1 tương ứng với việc đồng hóa tại tất cả các bước tích phân, M = 2 sẽ thực hiện đồng hóa tại từng 2 bước tích phân một. Với xấp xỉ bậc hai, L được cho bởi ]))(([))(( 11 ∑∑ == ++≈ i j j M i i txttxt JIJIL δδ (13) Vì chúng ta làm việc tường minh trong không gian ℜ40, mô hình liên hợp đơn giản là chuyển vị của mô hình tiếp tuyến (13). 3.3. Lực nhiễu động Ngoại trừ trạng thái thực trong đó lực tác dụng được biết chính xác với F = 8.0, một tổ hợp gồm n thành phần các tích phân mô hình sẽ không có giá trị lực tác dụng chính xác mà được lấy từ một tập hợp gồm n phép lấy ngẫu nhiên có phân bố chuẩn với độ lệch chuẩn δF = 1.0 được cộng tại từng bước tích phân. Như K.Q. Chánh / Tạp chí Khoa học ĐHQGHN, Khoa học Tự nhiên và Công nghệ 26, Số 3S (2010) 310‐316 314 được thảo luận trong phần 2, tất cả các phép thử phải có cùng một điều kiện ban đầu trong quá trình tích phân tổ hợp. Ma trận sai số hiệp biến mô hình Q sẽ thu được bắng cách lấy mẫu n đầu ra của tích phân tại các thời điểm đồng hóa. Đối với số liệu quan trắc cần thiết cho việc đồng hóa, các nhiễu động với phân bố Gauss và độ lệch chuẩn bằng 1.0 sẽ được cộng vào thành phần trạng thái thực của mô hình. 3.4. Kết quả Để đánh giá độ chính xác của các phương pháp khác nhau, sai số căn quân phương (rms) của sai số giữa trạng thái phân tích và trạng thái thực tại các thời điểm đồng hóa sẽ được sử dụng với định nghĩa như sau: 2/140 1 2)( 40 1 ⎭⎬ ⎫ ⎩⎨ ⎧ −= ∑ =k t k a k xxRMS Hình 1 chỉ ra một sự so sánh của sự tiến triển theo thời gian của rms cho các trường hợp không có hiệu chỉnh sai số mô hình (NOC), hiệu chỉnh sai số bằng phương pháp lực nhiễu động (PF20) với 20 thành phần, và phương pháp tăng cấp nhân điển hình (INF). Chúng ta có thể nhận thấy dễ dàng rằng đối với tất cả các cửa sổ đồng hóa M từ 1 đến 8 bước tích phân, PF20 cho một kết quả tốt hơn và rất ổn định so với INF. Cả PF20 và INF đều cho kết quả tốt hơn so với trường hợp sai số mô hình không được tính đến trong mô hình như nhìn thấy trong trường hợp NOC. Do các mô hình tiếp tuyến có độ sai lệch tích luỹ tăng theo khoảng đồng hóa M, có thể nhận thấy từ Hình 1 là với giá trị M lớn thì rms cũng tăng nhanh. Với M > 15, lọc Kalman sẽ phân kỳ trong tất cả các phương pháp. Sự ổn định của phương pháp PF20 so với INF là có thể hiểu được nếu chúng ta chú ý là PF20 cho phép tính đến sai số nội tại của mô hình trong khi INF chỉ phụ thuộc vào tần số đồng hóa. Nói một cách khác, INF sẽ giả thiết là sai số mô hình tỷ lệ với sai số của ma trận hiệp biến phân tích. Hình 1. Sự tiến triển theo thời gian của rms giữa trạng thái phân tích và trạng thái thực cho các thí nghiệm NOC (đường liền nhạt), INF với thừa số nhân 0.03 (đường chấm), và PF20 (đường liền đậm) với M = 1, 2, 4, và 8. Để xem xét thêm độ nhậy của phương pháp lực nhiễu động, một loạt các thí nghiệm đã được tiến hành trong đó số lượng các thành phần tổ hợp tăng dần từ 10 đến 100 với M cố định bằng 14 như được chỉ ra trong Hình 2. Mặc dù các thí nghiệm với nhiều thành phần tổ K.Q. Chánh / Tạp chí Khoa học ĐHQGHN, Khoa học Tự nhiên và Công nghệ 26, Số 3S (2010) 310‐316 315 hợp cho rms nhỏ hơn như mong đợi (Hình 2), có thể nhận thấy rằng sự giảm của rms dường như bão hóa rất nhanh chỉ với 20 thành phần tổ hợp. Điều này chỉ ra rằng chỉ cần với một số ít các thành phân tổ hợp cũng có thể nắm bắt tốt cấu trúc và đặc trưng của trường sai số mô hình, một ưu điểm rất có ý nghĩa đối với các tính toán thực tế trong đó khối lượng tính toán lớn của mô hình không cho phép chúng ta có nhiều thành phần tổ hợp. 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 NO PF10 PF20 PF40 PF60 PF80 PF100 MF Hình 2. Sai số rms lấy trung bình trong khoảng 1000 bước tích phân cho phương pháp lực nhiễu động với số thành phần tổ hợp là 10, 30, 50, 100 (xám), hiệu chỉnh tăng cấp (xám nhạt), và không có hiệu chỉnh sai số mô hình (xám đậm). 4. Ứng dụng mở rộng Như đã được thảo luận trong phần 2, khuôn khổ lí thuyết trong phương pháp lực nhiễu động chỉ có ý nghĩa đối với các hệ như được cho bởi phương trình (1) với một số bậc tự do nhỏ. Đối với các hệ phức tạp hơn như là hệ trái đất-khí quyển, sẽ gần như không thể sử dụng phương pháp lọc Kalman toàn phần do số chiều của mô hình là quá lớn. Do đó, lọc Kalman tổ hợp phải được sử dụng [8]. Đối với các hệ như vậy, một sự mở rộng tự nhiên của phương pháp lực nhiễu động là lấy mẫu một cách trực tiếp các đầu ra của một tích phân tổ hợp mà có cùng điều kiện ban đầu như đã được xem xét trong phần 3. Vấn đề duy nhất phải chú ý là khi tính toán các sai số mô hình với phương pháp lực nhiễu động này là làm thế nào để có thể tạo ra một bộ nhiễu thích hợp. Đây là một câu hỏi mà phụ thuộc rất nhiều vào mô hình mà chúng ta có và bài toán chúng ta cần thiết phải giải quyết. Một cách cụ thể, giả thiết rằng chúng ta có một mô hình bão khu vực mà chúng ta muốn nghiên cứu tính dự báo của các bản tin dự báo đựờng đi của bão. Các nghiên cứu trước đã chỉ ra rằng đường đi của bão phụ thuộc rất nhiều vào các yếu tố môi trường không tính được trong mô hình số ví dụ các sơ đồ tham số hóa đối lưu hay lớp biên. Trong trường hợp này, một cách rõ ràng nhất để tạo ra trường lực nhiễu động là sử dụng ngay các sơ đồ tham số khác nhau để tạo ra bộ nhiễu. Các tính toán liên tục của trường sai số mô hình với các sơ đồ tham số hóa này không đòi hỏi các mô hình tiếp tuyến hay mô hình liên hợp và do đó sẽ có ý nghĩa thực tế hơn. Với sai số mô hình ước lượng được bằng cách này, bộ lọc Kalman tổ hợp có thể được kết hợp để tạo ra bộ nhiễu trên cùng một tổ hợp thay vì chạy 2 tổ hợp riêng rẽ cho điều kiện ban đầu và cho sai số mô hình. 5. Kết luận Trong bài nghiên cứu này, phương pháp xác định sai số mô hình bằng cách tạo ra bộ nhiễu ngẫu nhiên của lực tác dụng đã được khảo sát l ý thuyết một cách tường minh. Phương pháp nhiễu lực được dựa trên giả thiết rằng nguồn gốc lớn nhất của sai số mô hình là do các hiểu biết không đầy đủ của các quá trình vật lý trong mô hình, đặc biệt trong các hệ phức tạp như khí quyển đại dương. Phương pháp lực nhiễu động trên đã được kiểm nghiệm trên mô hình Lorenz và đã chỉ ra một vài tính chất nổi bật bao gồm 1) sự ổn định của thuật toán đối với một khoảng rộng của cửa số đồng hóa, 2) độ chính xác cao hơn phương pháp tăng cấp bội thuần tuý , và 3) K.Q. Chánh / Tạp chí Khoa học ĐHQGHN, Khoa học Tự nhiên và Công nghệ 26, Số 3S (2010) 310‐316 316 Tài liệu tham khảo độ chính xác được duy trì tốt ngay cả với một số ít các thành phần tổ hợp. Điều này rất được trông đợi cho các ứng dụng thực tế trong đó khối lương tính toán rất lớn của mô hình nghiệp vụ không cho phép chúng ta có nhiều thành phần tổ hợp. Mở rộng của phương pháp lực nhiễu động cho hệ thống với nhiều bậc tự do cũng đã được thảo luận. Nghiên cứu và ứng dụng chi tiết hơn của phương pháp lực nhiễu động sẽ được trình bày trong nghiên cứu tới. [1] J. L. Anderson, and S. L. Anderson, A Monte Carlo implementation of the non-linear filtering problem to produce ensemble assimilations and forecasts. Mon. Wea. Rev., 127 (1999) 2741. [2] H. L. Mitchell, and P. L. Houtekamer A