Tổng quan hệ thống đồng hóa lọc Kalman tổ hợp và ứng dụng cho mô hình dự báo thời tiết WRF

một trạng thái ban đầu mới và sai số của trạng thái ban đầu này cho quá trình dự báo tiếp theo. Chúng ta sẽ đi vào từng quá trình một cách chi tiết hơn trong các phần tiếp theo. Bước dự báo Giả thiết khí quyển tại một thời điểm i nào đó được đặc trưng bởi một trạng thái a ix với một sai số a iε . Đầu tiên chúng ta sẽ dự báo cho trạng thái đến thời điểm i + 1 sẽ cho bởi: )(1 a i f i M xx  (1) trong đó M là mô hình dự báo. Do mô hình này là không hoàn hảo, dự báo bằng mô hình này sẽ có một sai số nào đó kể cả khi điều kiện ban đầu là chính xác. Gọi sai số nội tại này của mô hình là , khi đó một cách lý thuyết giá trị sai số này sẽ được xác định như sau:  )(1 titi M xx (2) trong đó t ii )1( x là trạng thái thực của khí quyển tại thời điểm i (i + 1). Chúng ta sẽ giả thiết rằng sai số nội tại này là không lệch và ma trận sai số hiệp biến của nó được cho bởi một ma trận Q, nghĩa là  T Q;0 (3) Song song với dự báo trạng thái, chúng ta sẽ dự báo cả sai số từ thời điểm thứ i đến thời điểm thứ i + 1 sử dụng mô hình tiếp tuyến L được định nghĩa dựa trên dạng biến phân của phương trình (1) như sau: iii M xxLx x x x  )()(1     (4) Với mô hình tiếp tuyến L này, sai số của trạng thái tại thời điểm thứ i + 1 sẽ được cho bởi i a ii εxLε )(1  (5) Trong thực tế, chúng ta không bao giờ biết được sai số tuyệt đối thực i và như thế không thể dự báo được sai số cho bước tiếp theo. Tuy nhiên, trong đa số các trường hợp, chúng ta lại có thể biết hoặc xấp xỉ được đặc trưng thống kê của sai số được đặc trưng bởi ma trận sai số hiệp biến P  . Thêm vào đó, ma trận này cũng sẽ được sử dụng để đồng hóa cho bước tiếp theo. Do đó, chúng ta sẽ viết lại (5) cho ma trận sai số hiệp biến thay vì cho sai số tuyệt đối i. Lưu ý theo định nghĩa rằng t i f i f i xxε  , t i a i a i xxε  chúng ta sẽ có mối quan hệ sau QLLP η)εLη)εL xxxxεεP T     Ta i a i a i Tt i f i t i f i T ii f i (( ))(( 11111 (6) Chú ý thêm rằng chúng ta đã giả thiết là sai số mô hình  và sai số trạng thái a iε là không có tương quan với nhau. Như vậy, cho trước giá trị sai số mô hình Q, mô hình M, và mô hình tiếp tuyến L, phương trình (2) và (6) cấu thành một quá trình dự báo cơ bản trong bước dự báo K.Q. Chánh / Tạp chí Khoa học ĐHQGHN, Khoa học Tự nhiên và Công nghệ 27, Số 1S (2011) 17-28 20 theo đó trạng thái a ix và sai số a iε tại thời điểm i sẽ được dự báo đến thời điểm i + 1. Bước phân tích Trong bước phân tích tiếp theo, giả sử tại thời điểm i + 1, chúng ta có một bộ số liệu quan trắc yo với sai số quan trắc là o. Nhiệm vụ của chúng ta trong bước này là phải kết hợp được trạng thái dự báo f i 1x và sai số f i 1P với quan trắc để tạo được một bộ số liệu đầu vào mới tốt hơn tại thời điểm i + 1. Lưu ý rằng mặc dù a ix là ước lượng tốt nhất của trạng thái khí quyển tại thời điểm i, giá trị dự báo f i 1x tại thời điểm i + 1 lại không phải là tốt nhất do sai số của mô hình và của a ix . Do đó chúng ta cần phải đồng hóa tại thời điểm i + 1 để trạng thái dự báo không bị lệch khỏi trạng thái thực tại các thời điểm này. Một cách hình thức, chúng ta sẽ ước lượng trạng thái khí quyển mới tốt hơn tại thời điểm i + 1 như sau: )]([ 111 f i of i a i H   xyKxx (7) trong đó H là một toán tử quan trắc nội suy từ trường mô hình sang các giá trị điểm lưới, và K là ma trận trọng số. Một cách trực quan, ma trận K càng lớn, ảnh hưởng của quan trắc lên trường phân tích càng nhiều. Do đó, ma trận K rất quan trọng và phải được dẫn ra một cách tối ưu nhất có thể. Để thuận tiện cho việc suy dẫn K, chúng ta định nghĩa một vài biến sai số sau: t i a i a i xxε  , t i f i f i xxε  , )( ti oo i H xyε  (8) Để tìm ma trận K, chúng ta trước hết phải tính ma trận sai số hiệp biến Pa cho trạng thái phân tích a i 1x và sau đó cực tiểu hóa ma trận này. Theo định nghĩa:     T Ta i )( )1 t 1i a 1i t 1i a 1i a 1i a 1i x-x)(x-x (εεP (9) Thay (7) vào (9) và xắp xếp lại, chúng ta sẽ thu được:      Ta i )εH-εKε)εH-εKεP f 1i o 1i f 1i f 1i o 1i f 1i ((((1 (10) trong đó ma trận H là tuyến tính hóa của toán tử quan trắc H. Đặt   Tf i )1 f 1i f 1i (εεP ,   T)o 1i o 1i (εεR , và giả thiết trạng thái nền không có tương quan với trạng thái phân tích, chúng ta sẽ thu được từ (10) phương trình sau: TTf i a i KRKKHIPKHIP   )()( 11 (11) Ma trận trọng số K sẽ cực tiểu hóa vết của ma trận sai số a i 1P khi và chỉ khi 0))(( 1     a itrace P K (12) trong đó trace() ký hiệu vết của ma trận. Ở đây, đạo hàm theo ma trận sẽ được hiểu là đạo hàm từng thành phần của ma trận. Lý do cho việc cực tiểu hóa vết của ma trận thay vì trực tiếp ma trận là do tổng các thành phần trên đường chéo của ma trận a i 1P sẽ chính là bình phương của tổng sai số căn quân phương trong trường hợp các biến là không tương quan chéo. Do vết của một ma trận là bảo toàn trong các phép biến đổi trực chuẩn, chúng ta luôn có thể chéo hóa ma trận sai số a i 1P để đưa về một cơ sở mà trong đó tổng sai số căn quân phương sẽ là tổng của các thành phần đường chéo. Lấy đạo hàm vết của ma trận a i 1P , chúng ta khi đó sẽ thu được từ (11) và (12) 1 11 )(    Tf i Tf i HHPRHPK (13) Với giá trị ma trận trọng số K cho bởi (13) ở trên, giá trị cực tiểu của ma trận sai số hiệp biến phân tích khi đó sẽ thu được bằng cách thay (13) vào (11). Biến đổi tường minh chúng ta sẽ thu được: .)( 11 f i a i   PKHIP (14) K.Q. Chánh / Tạp chí Khoa học ĐHQGHN, Khoa học Tự nhiên và Công nghệ 27, Số 1S (2011) 17-28 21 Như vậy, ở bước phân tích này chúng ta đã thu được một ước lượng ban đầu mới tốt hơn từ một trạng thái dự báo (hay dự báo nền) và quan trắc cho trước. Sau khi thu được trạng thái mới a i 1x và ma trận sai số mới a i 1P , quá trình dự báo lại được lặp lại cho bước đồng hóa kế tiếp theo. Một cách tóm tắt, lọc Kalman được cho bởi minh họa trong hình 1. Hình 1. Minh họa hai bước chính của bộ lọc Kalman. Mặc dù có ưu điểm vượt trội so với các phương pháp đồng hóa biến phân khác, lọc Kalman cho bởi hệ các phương trình (1), (6), (7), (13), (14) lại rất khó áp dụng trực tiếp trong các mô hình thời tiết có tính phi tuyến cao và bậc tự do rất lớn. Ba khó khăn chính của bộ lọc Kalman ở trên là 1) xây dựng mô hình tiếp tuyến L; 2) lưu trữ và thao tác các ma các trận sai số với số chiều có kích thước quá lớn; và 3) sai số nội tại của mô hình Q không được biết đầy đủ. Khó khăn thứ nhất có thể được giải quyết bằng cách sử dụng một biến thể của bộ lọc Kalman, gọi là Kalman tổ hợp mở rộng (EnKF) được đề xuất ban đầu bởi Evensen năm 1994. Khó khăn thứ hai được khắc phục bằng cách địa phương hóa các số liệu quan trắc xung quanh từng điểm nút lưới (localization) hoặc đồng hóa lần lượt từng giá trị quan trắc theo chuỗi (serial). Có nhiều các biến thể của EnKF dựa trên việc địa phương hóa hay xử lý chuỗi. Hai trong số đó sẽ được trình bày trong mục sau. Về sai số nội tại của mô hình, đây là một hướng phát triển còn mở của bộ lọc Kalman trong thời gian gần đây và có rất nhiều phương pháp xử lý. Ví dụ kỹ thuật tăng cấp cộng tính, kỹ thuật tăng cấp nhân, kỹ thuật hiệu chỉnh độ lệch hệ thống, kỹ thuật cộng nhiễu ngẫu nhiên. Trong nghiên cứu gần đây, tác giả có trình bày một kỹ thuật khác, tạm gọi là kỹ thuật nhiễu lực [25]. Chi tiết hơn về các thuật toán hiệu chỉnh sai số mô hình có thể được tham khảo trong [26]. Như đã đề cập trong phần giới thiệu, quá trình đồng hóa phải bao gồm hai bước chính là phân tích khách quan và ban đầu hóa. Trong bước phân tích của bộ lọc Kalman trình bày ở trên, chúng ta thấy rằng dường như bộ lọc Kalman chưa có quá trình ban đầu hóa một cách cụ thể. Tuy nhiên, các phân tích chi tiết cho thấy trong thực tế, bộ lọc Kalman đã tính đến quá trình ban đầu hóa một cách nội tại trong bước dự báo. Điều này là do trong bước dự báo này, ma trận sai số hiệp biến nền sẽ được tích phân theo thời gian. Do đó, các tương quan chéo giữa các biến động lực sẽ được hiệu chỉnh theo thời gian. Ở một giới hạn đủ dài, ma trận sai số hiệp biến nền thu được từ bộ lọc này sẽ có khả năng phản ánh được các tương quan chéo giữa các biến động lực và như vậy thông tin quan trắc thu được của bất kỳ một biến nào cũng sẽ được cập nhất cho tất cả các biến mô hình khác. Đây chính là ưu điểm của bộ lọc Kalman, đặc biệt trong vùng vĩ độ thấp tại đó không tồn tại một ràng buộc lý thuyết tường minh cho các mối quan hệ động lực giữa các biến giống như trong vùng ngoại nhiệt đới. (Trong vùng ngoại nhiệt đới, mối quan hệ động lực giữa các biến thường được sử dụng là mối quan hệ địa chuyển, cân bằng thủy tĩnh, hay cân bằng gradient.) Chúng ta sẽ xem xét phần này chi tiết hơn trong các mục sau. 3. Lọc Kalman tổ hợp a) Lọc Kalman tổ hợp Do khả năng phát triển mô hình tiếp tuyến và tích phân ma trận sai số hiệp biến theo thời gian với mô hình tiếp tuyến là không thực tế K.Q. Chánh / Tạp chí Khoa học ĐHQGHN, Khoa học Tự nhiên và Công nghệ 27, Số 1S (2011) 17-28 22 trong các mô hình dự báo thời tiết, lọc Kalman phải được cải tiến để có thể áp dụng được cho các bài toán nghiệp vụ. Một cách tiếp cận phổ biến nhất dựa trên tích phân ngẫu nhiên Monte- Carlo theo đó một tập các đầu vào được tạo ra xung quanh một giá trị trường phân tích cho trước. Lưu ý rằng tập đầu vào này không phải được lấy bất kỳ mà được tạo ra dựa theo phân bố xác xuất cũng như giá trị sai số của trường phân tích aP tại từng thời điểm1. Ví dụ nếu phân bố của trường phân tích có dạng phân bố chuẩn Gauss, khi đó tập đầu vào của trường phân tích sẽ phải tuân theo phân bố: )()()( 2 1 2/12/ 1 ||)2( 1 )( aaaTaa a ep an a xxPxx P P x     (15) trong đó a x là giá trị trường phân tích trung bình tổ hợp thu được từ bước phân tích của lọc Kalman như trình bày trong phần 2. Với một tập K các đầu vào { a kx }k=1..K sinh ra từ phân bố (15), chúng ta có thể thu được ma trận sai số hiệp biến dự báo cho bước thời gian tiếp theo như sau: Tff k K k ff k f K )()( 1 1 1 xxxxP      (16) trong đó )]([)( 1 i a ki f k tMt xx  . Các nghiên cứu trong [4] đã chỉ ra rằng một tổ hợp khoảng 25-50 các thành phần là đủ để lọc Kalman tổ hợp phát huy tác dụng. Ngoài việc giản lược quá trình phát triển mô hình tiếp tuyến, cách tiếp cận EnKF có một vài ưu điểm nổi trội bao gồm 1) rút bớt các tính toán với các ma trận có số chiều lớn; 2) không cần tuyến tính hóa mô hình cũng như mô hình liên hợp (adjoint); 3) cung cấp một tổ hợp các nhiễu ban đầu tối ưu hóa cho việc dự báo các ma trận sai số hiệp biến. Mặc dù có nhiều điểm thuận lợi song đến tại thời điểm này, chỉ có _______ 1 Để đơn giản, trong phần này chúng ta sẽ chỉ làm việc với một thời điểm cụ thể trong bước phân tích của bộ lọc Kalman và do đó chỉ số thời gian i sẽ được bỏ qua. trung tâm dự báo thời tiết của Canada đã đựa lọc Kalman tổ hợp vào dự báo nghiệp vụ do khối lượng tính toán lớn. EnKF hiện được coi là một hướng đi phát triển mạnh nhất trong thời gian tới cho các bài toán nghiên cứu có tính dự báo thấp như bão nhiệt đới và ngoại nhiệt đới, mưa lớn, hay dông. Trong các phần tiếp theo, chúng tôi sẽ trình bày hai biến thể được phát triển phổ biến nhất hiện này của bộ lọc EnKF là bộ lọc EnKF chuỗi (Serial EnKF), và bộc lọc EnKF địa phương biến đổi (LETKF) b) Lọc Kalman tổ hợp chuỗi Về mặt bản chất, lọc EnKF chuỗi (SEnKF) là quá trình trong đó số liệu quan trắc được đồng hóa lần lượt từng giá trị một. Quá trình này sẽ làm giảm kích thước cũng như khối lượng tính toán. Trong phần này, chúng tôi sẽ theo trình bày trong [17] (viết tắt là SZ03) do cách tiếp cận của SZ03 có tính ứng dụng cao và đã được sử dụng trong một loạt các bài toán đồng hóa quy mô vừa. Để thảo luận được rõ ràng, chúng ta viết lại phương trình cập nhật trạng thái phân tích như sau: ][ fofa xHyKxx  (17) .)( fa PKHIP  (18) trong đó, 1)(  TfTf HHPRHPK (19) Lưu ý rằng quá trình cập nhật trạng thái phân tích này sẽ chỉ được áp dụng cho trường trung bình tổ hợp. Để tạo ra được các thành phần tổ hợp phân tích tiếp theo, cần phải tạo thêm một bộ nhiễu phân tích và cộng vào trường phân tích trung bình tổ hợp a x như sẽ trình bày trong các phần tiếp theo. Để đồng hóa các quan trắc một cách hiệu quả, chúng ta viết lại tích hiệp biến Tf HP như sau: ),(),( ffTffTf HxxCovHxxCovHP  (20) K.Q. Chánh / Tạp chí Khoa học ĐHQGHN, Khoa học Tự nhiên và Công nghệ 27, Số 1S (2011) 17-28 23 trong đó chúng ta đã giả thiết rằng ma trận sai số hiệp biến dự báo fP được xấp xỉ bằng một tổ hợp K thành phần{ f kx }k=1..K như cho bởi (16). Một cách tường minh, (20) có thể được viết như sau:      K 1k xHHxxxHP Tjjk jj k Tf K ))(( 1 1 (21) trong đó    K k f k f K 1 1 xx (22) Về mặt ý nghĩa vật lý, tích ma trận (20) được hiểu như là một ma trận tương quan của trạng thái dự báo giữa các điểm nút lưới mô hình ( f x ) và các điểm mà tại đó quan trắc được cho ( f Hx ). Trong trường hợp chúng ta đồng hóa lần lượt từng giá trị quan trắc một, phương trình (17) và tích (20) sẽ được đơn giản hóa rất nhiều. Thật vậy, với chỉ một giá trị quan trắc, f Hx sẽ trở thành một vô hướng, và do đó tích Tf HP sẽ trở thành một vector, tạm gọi là vector c. Cũng như vậy, Tf HHP sẽ chỉ còn là một vô hướng, và do đó có thể viết lại số hạng nghịch đảo của )( Tf HHPR  như là nghịch đảo của một giá trị vô hướng, tạm gọi là d. Như vậy, (17) sẽ rút gọn thành dycxx foi f i a i /][ xH (23) với chỉ số chiều của mô hình i chạy trên toàn các nút lưới mô hình (i = 1… N) được viết tường minh. Phương trình cập nhật (23) phát biểu rằng với mỗi giá trị quan trắc được đưa vào y o, tương quan của giá trị quan trắc này với tất các điểm nút lưới khác (được biểu diễn bởi ci) sẽ lan truyền ảnh hưởng của quan trắc này lên tất cả các điểm nút lưới (lưu ý rằng chỉ số i chạy trên tất cả các điểm nút lưới). Mức độ ảnh hưởng của quan trắc tới các điểm nút lưới mô hình được cho bởi tham số d; sai số quan trắc càng nhỏ, mức độ ảnh hưởng của quan trắc sẽ càng lớn và ngược lại. Với một giá trị quan trắc điểm duy nhất, ma trận sai số hiệp biến phân tích theo phương trình (18) cũng được đơn giản hóa rất đáng kể: d T fa cc PP  (24) Phương trình cập nhật (24) phát biểu rằng sai số phân tích sẽ được giảm một lượng tỷ lệ với tương quan của biến được quan trắc với các điểm nút của mô hình. Ma trận sai số hiệp biến phân tích P a đóng vai trò quan trọng trong việc xây dựng bộ nhiễu cho SEnKF. Như đã đề cập ở trên, bước phân tích chỉ được áp dụng cho trạng thái trung bình tổ hợp. Để tạo ra một tổ hợp các đầu vào cho bước tích phân tiếp theo, chúng ta cần phải thực hiện thêm một bước cuối cùng là tạo ra bộ nhiễu phân tích. Bộ nhiễu phân tích này không thể chọn tùy ý mà phải được xây dựng dựa trên trường sai số hiệp biến phân tích P a cho bởi (24). Có rất nhiều lựa chọn khác nhau, xong một lựa chọn phổ biến nhất là sở dụng phương pháp căn quân phương ma trận. Gọi tổ hợp bộ nhiễu là }|...||{ 21 a K aaa xxxX  , chúng ta sẽ giả thiết rằng: Taaa K )( 1 1 XXP   (25) Một cách hình thức, nhiễu phân tích sẽ thu được như sau:   2/1)1( aa K PX  (26) Trong trường hợp số liệu quan trắc được đồng hóa lần lượt từng giá trị một theo chuỗi, chúng ta thu được kết quả sau: fa d XH c IX         (27) trong đó 12/1 ])/(1[  dR và }|...||{ 21 f K fff xxxX  . Một cách tường minh, chúng ta có: f k a k d xH c Ix         (28) và do đó, K.Q. Chánh / Tạp chí Khoa học ĐHQGHN, Khoa học Tự nhiên và Công nghệ 27, Số 1S (2011) 17-28 24 )( ffk aa k aa k d xxH c Ixxxx         (30) Gọi N là kích thước của trạng thái nền và K là số thành phần tổ hợp, các bước đồng hóa của lọc SEnKF như vậy sẽ bao gồm các bước như sau: - Bước 1: chọn một số liệu quan trắc duy nhất yo trong tập các số liệu quan trắc tại thời điểm i; - Bước 2: tính trung bình tổ hợp của trường dự báo fx theo công thức (22), và xây dựng ma trận nhiễu nền }|...||{ 21 ff K fffff xxxxxxX  ; - Bước 3: xây dựng ma trận H (1  N) và tính vector Tf HPc   RN theo công thức (21) dựa theo tổ hợp dự báo nền có được tại thời điểm i; - Bước 4: tính giá trị )( HcRd và 12/1 ])/(1[  dR . Lưu ý rằng vì chúng ta làm việc với từng giá trị quan trắc một, ma trận R sẽ chỉ là một giá trị vô hướng đặc trưng giá trị phương sai sai số của quan trắc yo; - Bước 5: Tính trạng thái phân tích trung bình tổ hợp xa dựa theo phương trình (23); - Bước 6: Tính nhiễu tổ hợp phân tích Xa theo (27) - Bước 7: Cộng nhiễu tổ hợp phân tích Xa với trung bình tổ hợp phân tích a x để thu được một tổ hợp phân tích Kk a k ..1}{ x ứng với quan trắc yo; - Bước 8: Gán trường nền tổ hợp bằng giá trị phân tích tổ hợp thu được ở bước 7 ở trên, nghĩa là Kk a k f k ..1}:{  xx . Lưu ý rằng bằng việc gán liên tục trường tổ hợp nền bằng giá trị trường phân tích mới, chúng ta sẽ tiết kiệm được bộ nhớ rất nhiều. - Bước 9: Quay lại bước 1 với số liệu quan trắc mới cho đến khi đồng hóa toàn bộ tập số liệu quan trắc có được tại thời điểm i. Có thể chứng minh được tường minh rằng trong các bài toán mô hình tuyến tính, đồng hóa chuỗi từng số liệu một sẽ cho kết quả tương tự như đồng hóa đồng thời tất cả các số liệu. Tuy nhiên chúng tôi sẽ không trình bày ở đây do hạn chế về độ dài của bài báo. Một chú ý là quá trình tính toán sẽ được giảm rất nhiểu nếu để ý rằng số liệu quan trắc tại một điểm cho trước có xu hướng chỉ ảnh hưởng các quan trắc ở lân cận nó. Do đó, thay vì cập nhật toàn bộ các điểm nút lưới với từng quan trắc, có thể chọn một lân cận xung quanh điểm quan trắc để cập nhật. Điều này không chỉ làm giảm khối lượng tính toán mà nó còn rất cần thiết để loại bỏ các tương quan chéo giả do số lượng thành phần tổ hợp là hữu hạn. Quá trình địa phương hóa số liệu quan trắc này sẽ được xem xét chi tiết hơn trong thuật toán đồng hóa tiếp theo. c)Lọc Kalman tổ hợp địa phương Có thể nhận thấy rằng trong phương pháp đồng hóa chuỗi ở trên, toàn bộ quá trình đồng hóa có tính kế thừa và không thể song song hóa để tăng tốc độ tính toán. Do đó, lọc SEnKF nói chung sẽ xử lý số liệu quan trắc không hiệu quả, đặc biệt trong trường hợp số liệu quan trắc là rất lớn. Mặc dù điều này có thể được khắc phục bằng cách xử lý lựa chọn số liệu quan trắc hoặc kết hợp các quan trắc quá gần nhau bằng phương pháp “siêu ép” (super-ob) số liệu, lọc SEnKF vẫn chưa được tối ưu hóa hoàn toàn cho các bài toán đồng hóa lọc tổ hợp. Trong phần này, chúng ta sẽ xem xét một biến thể khác của lọc EnKF cho phép xử lý số liệu một cách hiệu quả hơn rất nhiểu, được gọi là phương pháp lọc Kalman tổ hợp địa phương hóa biến đổi (Local ensemble transform Kalman filter, LETKF). Phương pháp này được đề xuất ban đầu bởi Hunt và ccs năm 2005 và đã được đưa vào thử K.Q. Chánh / Tạp chí Khoa học ĐHQGHN, Khoa học Tự nhiên và Công nghệ 27, Số 1S (2011) 17-28 25 nghiệm trong một vài nghiên cứu với các mô hình toàn cầu. Đây cũng chính là phương pháp sẽ được lựa chọn và đưa vào mô hình thời tiết WRF trong phần sau của chúng tôi. Một cách cơ bản, lọc LETKF là một phương pháp theo đó tại mỗi điểm nút lưới, chúng ta sẽ chọn một lân cận với kích thước cho trước. Với không gian con này, chúng ta sẽ chọn ra tất cả các quan trắc cho được bên trong không gian này và tạo ra một vector quan sát lân cận. Sau đó sử dụng ma trận nhiễu tổ hợp nền để biến đổi từ không gian con căng bởi số điểm nút lưới địa phương sang không gian con căng bởi số thành phần tổ hợp. Điều này sẽ làm giảm đáng kế khối lượng tính toán ma trận vì không gian tổ hợp thường nhỏ hơn không gian địa phương rất nhiều. Do đó, các phép toán ma trận sẽ có độ chính xác cao hơn. Để minh họa thuật toán một cách rõ ràng, nhắc lại rằng ma trận nhiễu tổ hợp nền f X (có số chiều N  K) được định nghĩa như sau: }|...||{ 21 ff K fffff xxxxxxX  (31) trong đó fx là trung bình tổ hợp. Gọi w là một vector biến đổi trong không gian tổ hợp được định nghĩa như sau: wXxx bb  (30) Khi đó, hàm giá trong không gian tổ hợp địa phương sẽ chuyển thành ][ }])([)({)1()( 1 wXx wXXXXIww bb bTbbTbT J kJ     (31) trong đó ][ wXx bbJ  là hàm giá trong không gian mô hình. Hàm giá sẽ được cực tiểu hóa nếu w là trực giao với không gian con rỗng của toán tử Xf. Lấy đạo hàm của )(wJ  theo w, chúng ta sẽ thu được giá trị aw làm cực tiểu hóa hàm giá (31) như sau: )]([)( 01 fTfaa H xyRYPw    (32) trong đó )](),....,(),([ 21 ff K fffff HHH xxxxxxY  (33) và 11 ])()1[(  fTfa K YRYIP  (34) Như vậy, trong không gian tổ hợp, ma trận trọng số thu được từ (32) sẽ có dạng: 1)(ˆ  RYPK Tfa  và do đó ma trận trọng số K trở thành 1)(ˆ  RYPXKXK Tfaff  (35) Với ma trận K thu được ở trên, giá trị trạng thái phân tích trung bình tổ hợp tại điểm nút chúng ta đang quan tâm sẽ được cho bởi ][ fofa xHyKxx  (36) Cũng giống như trong bộ lọc SEnKF, nhiệm vụ cuối cùng của chúng ta là xây dựng bộ tổ hợp các trạng thái phân tích. Để làm điều đó, chúng ta chú ý rằng: Tfff K )( 1 1 XXP   ; và Taaa K )( 1 1 XXP   Kết hợp (18) và (35) với chú ý trên, chúng ta thu được: TffTfaafTaa KK ))()()ˆ((ˆ 1 1 )( 1 1 11 XYRYPPXXX      (37) Sử dụng (34), chúng ta sẽ thu được: TfafTaa K )(ˆ)1()( XPXXX  (38) và do đó, 2/1]ˆ)1[( afa K PXX  (39) Quá trình đồng hóa theo bộ lọc LETKF như vậy có thể được tóm tắt như sau: - Bước 1: tại mỗi điểm nút lưới, chọn một vùng thể tích lân cận bao xung quanh điểm nút đó để xây dựng ma trận nhiễu nền địa phương )](),....,[( 1 ff K fff xxxxX  ; - Bước 2: Trong mỗi thể tích lân cận, tìm tất cả các quan trắc bên trong thể tích lân cận này và xây dựng ma trận quan trắc nhiễu nền )](),....,([ 1 ff K fff HH xxxxY  (nếu H là toán tử tuyến tính, khi đó ff HXY  ). K.Q. Chánh / Tạp chí Khoa học ĐHQGHN, Khoa học Tự nhiên và Công nghệ 27, Số 1S (2011) 17-28 26 Đồng thời xây dựng ma trận sai số quan trắc R ứng với các quan trắc bên trong thể tích; - Bước 3: Tính ma trận sai số hiệp biến biến đổi a P  theo (34) và sau đó ma trận trọng số K theo (35); - Bước 4: Cập nhật giá trị trung bình tổ hợp địa phương a x theo (36); - Bước 5: Tính ma trận nhiễu phân tích Xa theo (39) và cộng vào a x để thu được tổ hợp phân tích lân cận mới ( a k aa k Xxx  ); - Bước 6: Chọn điểm giữa của vector tổ hợp vector phân tích lân cận a kx và gán điểm này cho điểm nút lưới chọn ở bước 1 - Bước 7: Quay trở lại bước 1 và lặp cho đến hết tất cả các điểm nút lưới. Có thể nhận thấy dễ dàng trong các bước tính toán ở phía trên rằng các điểm nút lưới khác nhau được thực thi một cách hoàn toàn độc lập với nhau. Đây là một ưu điểm của lọc LETKF vì chúng ta có thể song song hóa bộ lọc này một cách hiệu quả bằng cách chia các phần công việc độc lập cho các lõi tính toán khác nhau. Điều này cho phép tăng tính hiệu quả tính toán lên rất nhiều so với SEnKF. 4. Hệ thống dự báo WRF-LETKF Trong nghiên cứu này, chúng tôi xây dựng thử nghiệm bộ lọc LETKF cho mô hình thời tiết phổ dụng nhất hiện nay là mô hình WRF. Việc thiết kế hệ thống này, tạm gọi là WRF-LETKF, được phát triển theo quy trình ứng dụng nghiệp vụ chuẩn với tất cả các quá trình vào ra, cập nhật số liệu, xử lý đồng hóa, kiểm tra chất lượng quan trắc, tạo điều kiện biên tổ hợp, và dự báo tổ hợp được tiến hành một cách tự động và đồng bộ hóa theo thời gian thực. Sơ đồ thiết kế hệ thống được biểu diễn minh hóa như trong Hình 2. Hình 2. Minh họa sơ đồ hệ thống dự báo tổ hợp WRF-LETKF. Số liệu quan trắc đầu tiên sẽ được xử lý kiểm định chất lượng thông qua bộ chương trình chuẩn WRFDA cho trong mô hình WRF. Quá trình kiểm định chất lượng này sẽ xác định các sai số cho các mực và các biến quan trắc tương ứng. Số liệu quan trắc sau khi được kiểm định sẽ được kết hợp với số liệu dự báo tổ hợp hạn rất n