Thử nghiệm dự báo lượng mưa ngày bằng phương pháp dùng mạng thần kinh nhân tạo hiệu chỉnh sản phẩm mô hình số

Tạp chí Khoa học đhqghn, KHTN & CN, T.xxII, Số 1PT., 2006 Thử nghiệm dự báo l−ợng m−a ngày bằng ph−ơng pháp dùng mạng thần kinh nhân tạo hiệu chỉnh sản phẩm mô hình số Hồ Thị Minh Hà, Nguyễn H−ớng Điền Khoa Khí t−ợng-Thuỷ văn và Hải d−ơng học Tr−ờng Đại học Khoa học Tự nhiên, ĐHQG Hà Nội 1. Mở đầu Trên thế giới có rất nhiều nhà khí t−ợng đã quan tâm tới mạng thần kinh nhân tạo (MTKNT) vì nó có khả năng học và xử lý tính phi tuyến của các quá trình biến đổi phức tạp mà các ph−ơng pháp dự báo khác không thực hiện đ−ợc. MTKNT đã đ−ợc ứng dụng trong rất nhiều lĩnh vực khoa học nh− điều khiển tự động, nhận dạng dấu vân tay, ... và cũng đ−ợc các nhà khí t−ợng sử dụng trong dự báo thời tiết, khí hậu. Có thể kể ra một số công trình nghiên cứu đã sử dụng MTKNT nh− dự báo giáng thuỷ hạn dài ở California, mô hình hoá ôzôn, dự báo m−a đá nguy hiểm, dự báo giáng thủy hạn ngắn,... Tuy nhiên, trong n−ớc mới chỉ có rất ít các nghiên cứu về lĩnh vực này và ch−a có công trình nào đ−ợc áp dụng vào thực tiễn. ở đây chúng tôi sử dụng MTKNT để thử nghiệm hiệu chỉnh sản phẩm của mô hình số, cụ thể là l−ợng m−a của mô hình phân giải cao HRM, theo số liệu quan trắc trong 3 tháng mùa m−a (6, 7, 8) năm 2004 nhằm nâng cao chất l−ợng của sản phẩm dự báo số. Tr−ớc hết ta sẽ xem xét cấu trúc và hoạt động của một mạng thần kinh đơn giản cũng nh− ph−ơng pháp học giảm gradient của MTKNT, sau đó sẽ tiến hành hiệu chỉnh sản phẩm của mô hình số và đánh giá kết quả. 2. Khái niệm về MTKNT và ph−ơng pháp học giảm gradient 2.1. Cấu trúc và hoạt động của MTKNT Tổng có trọng số của các giá trị đầu vào Giá trị đầu vào 1 Giá trị đầu vào n Hàm truyền (a) (b) Hình 1. Sơ đồ mạng thần kinh sinh học (a) và cấu trúc của MTKNT một nút ẩn (b). - Mạng thần kinh sinh học bao gồm nhánh thần kinh (dendrites), tế bào thần kinh (cell body), trục thần kinh (axon) và các xung thần kinh (electrical spike) t−ơng ứng với 3 phần chính bao gồm đầu vào, lớp ẩn và kết xuất của MTKNT. 1 Hồ Thị Minh Hà, Nguyễn H−ớng Điền 2 - Đầu vào của MTKNT bao gồm N giá trị xi và các trọng số t−ơng ứng. Kết quả lấy tổng có trọng số của các x ija i sẽ đ−ợc nhập vào lớp ẩn : 1 n j o i i u a x a = = + ij∑ (1) - Lớp ẩn bao gồm hàm truyền, có thể là hàm sigma, hàm tang-hypebol,…Ng−ời ta th−ờng sử dụng hàm sigma có dạng: ( ) ( ) 1 1 exp x x σ = (2) + − Thay giá trị của u vào hàm truyền sẽ đ−ợc đầu ra của lớp ẩn, ký hiệu là yj j (j=1,.., H). - Các giá trị đầu ra của lớp ẩn là đầu vào của kết xuất, thực hiện t−ơng tự nh− đối với lớp ẩn đ−ợc giá trị kết xuất, ký hiệu là zk (k = 1,.., K). - Kết xuất đích mà ta muốn mạng học đ−ợc là các giá trị tk (k = 1,.., m). xi ∑ = += n i ijioj axau 1 ( )jj uy σ= yj ∑ = += h j jkjok bybz 1 jkbija Hình 2. Quá trình lan truyền tiến - Một quá trình đi từ đầu vào, qua lớp ẩn và đến kết xuất là một quá trình lan truyền tiến (feed-forward) của mạng (Hình 2). Quá trình này cũng t−ơng tự nh− quá trình nhận, xử lý thông tin và truyền thông tin đến não bộ của một tế bào thần kinh sinh học. 2.2. Ph−ơng pháp học giảm gradient (giảm dốc nhất) Bản chất của ph−ơng pháp này là cực tiểu hoá hàm lỗi (sai số trung bình bình ph−ơng) giữa kết xuất zk và kết xuất đích tk . Hàm lỗi E là hàm của trọng số: ( ) ( )∑ = −= m k kk ztE 1 2 2 1ωr (3) với ωr là vectơ các trọng số, ứng với mỗi vectơ ωr ta có một giá trị sai số và các giá trị này lập thành một mặt lỗi trong không gian. Để dễ hình dung ta giả sử có 2 trọng số w1 và w2 và ta có thể biểu diễn hàm E là một mặt lỗi theo 2 biến này. Hình 3. Mặt lỗi w1 w2 E Thử nghiệm dự báo l−ợng m−a ngày bằng ph−ơng pháp... 3 Mục tiêu của mạng là tìm đ−ợc điểm trũng nhất trên mặt lỗi, nơi đó zk và tk gần nhau nhất. Ban đầu, khi chọn một cặp trọng số bất kỳ ta đ−ợc một điểm trên mặt lỗi và mạng phải đi theo một đ−ờng nào đó để tìm đ−ợc điểm trũng nhất. H−ớng đi của mạng qua từng b−ớc là h−ớng làm giảm sai số của b−ớc tr−ớc nó. Biểu diễn toán học của quá trình cực tiểu hoá sai số đ−ợc trình bày d−ới dây. Thiết kế một MTKNT gồm đầu vào bao gồm n mẫu, mỗi mẫu gồm các giá trị input xi (i=1..N), lớp ẩn gồm các nút yj (j=1..H) và kết xuất là các giá trị đích tk (k=1..K). Ký hiệu lại hàm truyền sigma là g(x): ( ) ( )xxg −+= exp1 1 (4) x1 x2 … xN a0 a11 a22 a1H aN2 y1 y2 … yH b0 b11 b21 b1K bH1 b2k … u1 u2 … uH v1 … vK z1 … … zK t1 … … tK g(uj) g(vk) Đ ầ u v à o Lớp ẩ n đ ầ u r a đ íc h Hình 4: Các thành phần chủ yếu của một mẫu của MTKNT. Đối với mỗi mẫu, tr−ớc hết ta tính tổng trọng số các giá trị đầu vào: (5) ∑ = += N i iijj xaau 1 0 . Hàm truyền sigma đ−ợc áp dụng cho lớp ẩn: ( )jjj uugy − 1 +== exp1)( ) (6) Lấy tổng trọng số các nút ẩn yj : (7) ∑ = += H j jjkkk ybbv 1 0 . áp dụng hàm truyền sigma cho lớp xuất: ( kkk vvgz −+== exp1 1)( (8) B−ớc (8) đ−ợc thực hiện nếu nh− giá trị đích là biến nhị phân. Nếu là các giá trị thực thì zk đ−ợc gán trực tiếp bằng vk t−ơng ứng. Sai số tổng cộng của mỗi mẫu là: Hồ Thị Minh Hà, Nguyễn H−ớng Điền 4 ∑ = −= K k kk tzE 1 2)( 2 1 (9) Để sai số này đạt cực tiểu thì đạo hàm của E theo các trọng số phải tiến dần đến 0. Ta tính đạo hàm riêng của E theo từng trọng số. Tr−ớc hết là theo trọng số của lớp xuất: jk k k k kjk b v v z z E b E ∂ ∂ ∂ ∂ ∂ ∂=∂ ∂ .. , (10) trong đó kk k tz z E −=∂ ∂ (11) )1.( kk k k zz v z −=∂ ∂ (12) và j jk k y b v =∂ ∂ (13) Khi đó: )).(1.(. kkkkj jk tzzzy b E −−=∂ ∂⇒ (14) Đặt: )).(1.( kkkkk tzzzp −−= (15) Vậy kj jk py b E .=∂ ∂⇒ (16) Tiếp theo là đạo hàm theo trọng số của lớp ẩn: ij i i j K k jk k k k kij i i j jij a u u y b v v z z E a u u y y E a E ∂ ∂ ∂ ∂ ⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ ∂ ∂ ∂ ∂ ∂ ∂=∂ ∂ ∂ ∂ ∂ ∂=∂ ∂ ∑ = ...... 1 (17) Sử dụng ph−ơng trình (16) ta đ−ợc: ij i i j K k jkk ij a u u y bp a E ∂ ∂ ∂ ∂⎟⎠ ⎞⎜⎝ ⎛=∂ ∂ ∑ = ... 1 , (18) trong đó )1.( jj i j yy u y −=∂ ∂ (19) và i ij i x a u =∂ ∂ (20) Vậy ijj K k jkk ij xyybp a E ).1.(. 1 −⎟⎠ ⎞⎜⎝ ⎛=∂ ∂⇒ ∑ = (21) Thử nghiệm dự báo l−ợng m−a ngày bằng ph−ơng pháp... 5 Nếu g(x) là hàm hybepol hay tang hypebol thì đạo hàm của E theo các trọng số sẽ khác. Các trọng số aij và bjk đ−ợc cập nhật sao cho gradient tổng cộng của sai số theo mọi trọng số trên tất cả các mẫu giảm đi. Gọi d là tổng gradient sai số theo mọi trọng số (ký hiệu chung là ω ) của b−ớc tr−ớc nó: i n i Ed ∑ = ⎟⎠ ⎞⎜⎝ ⎛ ω∂ ∂= 1 (22) với /∂ ∂ωE là đạo hàm riêng tổng cộng của E theo trọng số của mẫu thứ i, n là số l−ợng mẫu. H−ớng giảm dốc nhất là h−ớng ng−ợc lại của h−ớng d. Nh− vậy, mỗi trọng số sau một b−ớc sẽ đ−ợc cập nhật mới: ωωω ∆+=+ kk 1 (23) với: d.εω −=∆ (24) trong đó ε là hệ số học, là hệ số thực nghiệm do ng−ời sử dụng tự chọn, ở đây chọn hệ số học 05.0=ε . Theo dõi sự biến thiên của sai số E theo thời gian đến khi E giảm ít hoặc hầu nh− không giảm nữa thì dừng lại và chọn số thế hệ luyện phù hợp. Quá trình luyện đ−ợc thực hiện nhiều lần để tránh rơi vào các cực tiểu địa ph−ơng trên mặt lỗi hay còn gọi là bị rơi vào “bẫy”. Nh− vậy, trên lý thuyết MTKNT có thể tìm đ−ợc sai số cực tiểu, có thể bằng 0 và kết xuất sẽ hoàn toàn trùng khớp với kết xuất đích. Nh−ng có một số vấn đề làm cho MTKNT không tìm đến đ−ợc điểm trũng nhất của mặt lỗi. Việc học của MTKNT hoàn toàn độc lập với con ng−ời, con ng−ời không tác động đến bất cứ một hoạt động nào của nó nh−ng việc chọn số nút, số lớp ẩn, số thế hệ quá nhiều hoặc quá ít lại có thể làm cho MTKNT không tìm đ−ợc sai số cực tiểu. Ngoài ra mạng còn th−ờng gặp phải một vấn đề nữa là vấn đề quá khớp. Khi kiểm nghiệm cho thấy đối với tập số liệu luyện, sai số rất nhỏ nh−ng với tập số liệu độc lập dùng để đánh giá thì sai số lại tăng lên. Điều này phần lớn là do mạng đã luyện với quá nhiều thế hệ. Những vấn đề này chỉ có thể hạn chế đ−ợc nhờ việc luyện nhiều lần và lấy tổ hợp các kết quả nhận đ−ợc hoặc tiến hành kiểm nghiệm chéo (không trình bày chi tiết ở đây). Sau đây sẽ sử dụng MTKNT để thử nghiệm hiệu chỉnh cho l−ợng m−a dự báo của mô hình phân giải cao HRM. 3. Thử nghiệm hiệu chỉnh l−ợng m−a dự báo của mô hình HRM 3.1. Số liệu sử dụng - Số liệu bao gồm các giá trị l−ợng m−a ngày dự báo tr−ớc 24h bởi mô hình HRM với độ phân giải ngang 14km và 31 mực thẳng đứng, đ−ợc nội suy về 314 trạm trên khu vực Việt Nam và l−ợng m−a quan trắc tại các trạm đó trong thời gian từ tháng 6 đến tháng 8 năm 2004. - Hai tập số liệu này đ−ợc xử lý loại bỏ những ngày không đủ số liệu quan trắc hoặc dự báo. Sau khi luyện mạng xong sẽ phục hồi các vị trí thiếu giá trị bằng –9999 và đánh giá chất l−ợng dự báo theo các chỉ số, chủ yếu là l−ợng m−a trung bình ngày, sai số trung bình toàn ph−ơng (RMSE), sai số hệ thống (BIAS) và độ lệch tuyệt đối trung Hồ Thị Minh Hà, Nguyễn H−ớng Điền 6 bình (MAE) cho từng trạm, từng vùng (Bắc Bộ, Trung Bộ, Nam Bộ) và toàn bộ Việt Nam. 3.2. Chạy mô hình Phần xử lý số liệu đầu vào, số liệu đầu ra đ−ợc thực hiện bằng lập trình Fortran trên hệ điều hành Linux. Phần luyện mạng, kiểm nghiệm và hiển thị đ−ợc thực hiện nhờ phần mềm Matlab. a/ Phần luyện mạng: Số liệu đầu vào (input) của mạng là chuỗi giá trị l−ợng m−a ngày dự báo đã đ−ợc nội suy về 314 trạm trên Việt Nam. Đích cần học là các giá trị m−a quan trắc trên các trạm t−ơng ứng. Mạng đ−ợc thiết kế bao gồm 2 lớp ẩn và 1 lớp xuất. Lớp ẩn thứ nhất có 3 nút ẩn, lớp thứ 2 có 5 nút ẩn, số thế hệ học là 300 thế hệ. Số thế hệ này đ−ợc chọn sau một vài lần luyện mạng. Các bộ trọng số tại 2 lớp ẩn và 1 lớp xuất đ−ợc l−u lại để kiểm nghiệm. b/ Phần tái tạo, kiểm nghiệm và đánh giá: - Tái tạo: Nhân bộ trọng số với chính số liệu input đ−ợc đ−a vào mạng học, kết quả so sánh với đích học xem mức độ trùng khớp đ−ợc bao nhiêu. - Kiểm nghiệm: Nhân bộ trọng số ở trên với chuỗi dự báo của tập số liệu kiểm nghiệm, kết quả so sánh với chuỗi quan trắc của tập số liệu kiểm nghiệm. - Đánh giá: Hiển thị giá trị m−a trung bình theo thời gian cho 314 trạm đối với cặp số liệu input-tái tạo và input-kiểm nghiệm. Tính các chỉ số RMSE, BIAS, MAE đối với quá trình luyện và quá trình kiểm nghiệm. N 2 i i i 1 1 RMSE (fcst obs ) N = = −∑ ; N i i 1 1 MAE fcst obs N = = −∑ i ; N i i i 1 1 BIAS (fcst obs ) N = = −∑ trong đó N là số ngày của chuỗi, fcsti là dự báo thứ i, obsi là quan trắc thứ i t−ơng ứng. 4. Kết quả thử nghiệm Sau đây là một số kết quả dự báo l−ợng m−a ngày bằng ph−ơng pháp dùng MTKNT hiệu chỉnh sản phẩm của mô hình HRM. Số liệu luyện và dự báo bao gồm cả 314 trạm trên Việt Nam sau đó tách ra làm 3 chuỗi gồm 159 trạm trên Bắc Bộ (19.5-24 độ Bắc, 100-110 độ Đông), 95 trạm trên Trung Bộ (12-19.5 độ Bắc, 100-110 độ Đông) và 57 trạm trên Nam Bộ (7-12 độ Bắc, 100-110 độ Đông), loại ra 3 trạm nằm ngoài đảo và trên n−ớc Lào. Trên mỗi hình vẽ biểu diễn l−ợng m−a ngày trung bình theo thời gian của các trạm trên từng khu vực, đ−ờng chấm chấm là số liệu dự báo 24h của HRM, đ−ờng liền Thử nghiệm dự báo l−ợng m−a ngày bằng ph−ơng pháp... 7 nét có dấu hoa thị là số liệu đ−ợc hiệu chỉnh bởi MTKNT và đ−ờng liền nét là số liệu quan trắc t−ơng ứng. Hình 5. L−ợng m−a ngày (mm/ngày) trung bình tháng 6 và tháng 7 của Bắc bộ, tập số liệu input và tái tạo Hình 6. L−ợng m−a ngày (mm/ngày) trung bình tháng 8 của Bắc bộ, tập số liệu input và kiểm nghiệm Hình 7. L−ợng m−a ngày (mm/ngày) trung bình tháng 6 và tháng 7 của Trung bộ, tập số liệu input và tái tạo Hồ Thị Minh Hà, Nguyễn H−ớng Điền 8 Hình 8. L−ợng m−a ngày (mm/ngày) trung bình tháng 8 của Trung bộ, tập số liệu input và kiểm nghiệm Hình 9. L−ợng m−a ngày (mm/ngày) trung bình tháng 6 và tháng 7 của Nam bộ, tập số liệu input và tái tạo Hình 10: L−ợng m−a ngày (mm/ngày) trung bình tháng 8 của Nam bộ, tập số liệu input và kiểm nghiệm Thử nghiệm dự báo l−ợng m−a ngày bằng ph−ơng pháp... 9 Bảng 1: Đánh giá sai số tổng quát của quá trình luyện mạng và quá trình kiểm nghiệm BIAS (mm/ngày) MAE (mm/ngày) RMSE (mm/ngày) Trung bình dự báo (mm/ngày) Trung bình quan trắc (mm/ngày) HRM 6.585 7.547 9.835 13.654 Việt Nam MTKNT 0.002 2.710 3.563 7.071 7.069 HRM 9.692 10.265 12.516 17.872 Bắc Bộ MTKNT 0.276 2.577 3.609 8.457 8.180 HRM 3.337 4.721 6.065 9.014 Trung Bộ MTKNT -0.288 2.840 3.456 5.390 5.678 HRM 3.445 4.843 5.759 9.881 Luyện (06- 07/2004) Nam Bộ MTKNT -0.341 2.912 3.666 6.090 6.432 HRM 6.179 7.663 9.884 13.898 Việt Nam MTKNT -0.499 3.007 4.058 7.220 7.719 HRM 10.306 10.704 12.677 20.265 Bắc Bộ MTKNT -0.673 2.823 3.939 9.287 9.959 HRM 0.980 4.153 5.384 5.942 Trung Bộ MTKNT -0.518 3.454 4.462 4.444 4.961 HRM 3.469 5.249 6.355 9.730 Kiểm nghiệm (08/2004) Nam Bộ MTKNT -0.077 2.824 3.744 6.184 6.261 Từ hình 5 đến hình 10 và bảng 1 ta thấy MTKNT đã hiệu chỉnh sản phẩm của HRM về gần với quan trắc hơn, sai số ở tập luyện nhỏ hơn sai số trên tập kiểm nghiệm. Nhìn chung, tính trung bình cho cả chuỗi thời gian (2 tháng đối với tập luyện và 1 tháng đối với tập kiểm nghiệm) thì các đ−ờng biểu diễn l−ợng m−a trung bình theo thời gian tại 314 trạm của MTKNT và quan trắc khá sát nhau, trong khi đó sản phẩm của HRM th−ờng cho l−ợng m−a lớn hơn nhiều so với quan trắc. Điều này thể hiện rõ hơn trong các chỉ số ở bảng 1. Tính cho tất cả các ngày và tất cả các trạm thì l−ợng m−a dự báo trung bình của HRM khoảng hơn 13mm/ngày trong khi l−ợng m−a trung bình của quan trắc và MTKNT chỉ khoảng 7-8mm/ngày. Do đó, các sai số RMSE, BIAS và MAE của MTKNT nhỏ hơn hẳn so với HRM. Tuy nhiên, những ngày m−a lớn trong tập số liệu luyện rất ít so với số ngày không m−a hoặc m−a rất nhỏ, nên MTKNT học đ−ợc ít thông tin về m−a lớn. Vì vậy, MTKNT không hiệu chỉnh tốt đối với những ngày có m−a lớn. Nh− thấy trên hình 5 đến hình 10, tại một vài trạm th−ờng có m−a lớn trong 3 tháng mùa hè năm 2004, đ−ờng mô phỏng của MTKNT không gần với đ−ờng quan trắc và chỉ số BIAS có giá trị âm, trong khi đó dự báo của HRM lại nắm bắt tốt. 5. Kết luận - MTKNT đã hiệu chỉnh sản phẩm của HRM về gần với quan trắc hơn, sai số ở tập luyện nhỏ hơn sai số trên tập kiểm nghiệm. - MTKNT đã thiết lập có thể hiệu chỉnh tốt đối với những số liệu m−a vừa và m−a nhỏ trong khoảng 5-15mm, ch−a hiệu chỉnh tốt đối với những ngày m−a lớn và những trạm th−ờng có m−a lớn. Hồ Thị Minh Hà, Nguyễn H−ớng Điền 10 Để cải tiến kết quả hiệu chỉnh cần thay đổi một số tham số của mạng nh− hệ số học, số lớp ẩn, số nút ẩn của mỗi lớp ẩn và luyện mạng nhiều lần để chọn đ−ợc một MTKNT phù hợp nhất. Phức tạp hơn, ta có thể thiết lập nên MTKNT với mục tiêu học không chỉ là cực tiểu hoá sai số trung bình bình ph−ơng mà còn tăng thêm trọng số đối với những quan trắc đích có giá trị lớn để đ−a thêm thông tin về những ngày m−a lớn cho mạng học. Sau khi chọn đ−ợc mạng phù hợp nhất ta có thể l−u trữ lại mạng để dùng cho các mô phỏng tiếp sau mà không cần luyện lại mạng. Bài báo đ−ợc hoàn thành với sự hỗ trợ kinh phí của đề tài NCCB 733104. Tài liệu tham khảo 1. Caren Marzban and Arthur Witt, A Bayesian Neural Network for Severe-Hail Size Prediction, Weather and Forecasting, Volume 16, 2000, p600-610. 2. David Silverman and John A. Dracup, Artificial Neural Networks and Long-Range Precipitation Prediction in California, Monthly weather review, JANUARY 2000, p.57-66. 3. Narasimhan et al., Ozone Modeling Using Neural Networks, Monthly weather review, MARCH 2000, p291-296. 4. Robert J. Kuligowski and Ana P. Barros, Experiments in Short-Term Pricipitation Forecasting Using Artifical Neural Networks, Monthly weather review, 2000, Volume 126, p.470-482. 5. Tom M. Mitchell, Machine Learning, Chapter 4, 1996, p81-127. VNU. JOURNAL OF SCIENCE, Nat., Sci., & Tech., T.xXII, n01AP., 2006 An Experiment in Daily Rainfall Forecasting Using Artificial Neural Network to Adjust Numerical Model Output Ho Thi Minh Ha, Nguyen Huong Dien Department of Hydro-Meteorology & Oceanography College of Science, VNU Artifical Neural Network (ANN) is one kind of Decision Support Systems. It has the similar structure and behaviours to that of biological neuron. It can learn and simulate the information as a brain through a learning course. The most common learning method is folowing the gradient descent rule. In this paper it’s used for HRM rainfall output adjustment. The input data is rainfall field predicted by HRM (with the horizontal resolution14km and 31 vertical levels) and the observation data, respectively in the period of time between 01th June 2004 and 31th August 2004 at 314 stations in Vietnam. Data is divided into 2 parts. One part is used for neural network learning and the other one for prediction and test. The average rainfall, root mean square error (RMSE), mean absolute error (MAE) and mean bias error (BIAS) are the indices used to verify the results. All the indices show that the rainfall adjusted by ANN is closer to observation than that of HRM output in general. In particle, however the rainfall of ANN is smaller than that of HRM in the heavy rainfall days. The best results are get in the North of Vietnam and they’re very not good in the South of Vietnam. The reason is maybe the lack of station in this area.