Luận văn Áp dụng phương pháp ước lượng bất định khả năng (GLUE) cho dự báo lũ trên lưu vực sông Vệ

độ dẫn thủy lực tƣơng đƣơng với lƣợng ẩm trung bình ở thời gian t(mm/h), t bƣớc thời gian (h), Ki,s độ dẫn thủy lực bão hoà trong ô thứ i (mm/h), S,i độ rỗng (m3/m3), r,i lƣợng ẩm còn thừa trong ô thứ i(m3/m3), A chỉ số tính không liên kết giữa các lỗ rỗng tính toán từ phƣơng trình A=(2+3B)/B với B là chỉ số phân bố kích thƣớc độ rỗng của ô. Dòng sát mặt là thành phần quan trọng của cân bằng nƣớc trong đất. Nó là lƣợng nƣớc thấm xuống lớp đất mặt và di chuyển theo phƣơng ngang đến khi gia nhập vào kênh. Các yếu tố ảnh hƣởng đến dòng sát mặt là: các thuộc tính vật lý và độ dày của lớp đất: cấu trúc đất thô làm dòng chảy theo phƣơng thẳng đứng còn cấu trúc đất mịn hay lớp đất chống lại dòng chảy theo phƣơng thẳng đứng và dòng sát mặt có thể diễn ra nhanh chóng; thảm phủ: trực tiếp liên quan đến khả năng ngấm và ảnh hƣởng của vật chất vô cơ; địa hình: gradien độ dốc là yếu tố quan trọng quyết định vận tốc và lƣợng dòng chảy sát mặt; địa chất và khí hậu của khu vực nghiên cứu. Trong WetSpa cải tiến, dòng chảy sát mặt đƣợc giả định xảy ra sau quá trình ngấm và ngừng khi lƣợng ẩm của đất thấp hơn khả năng trữ. Lƣợng dòng chảy sát mặt đƣợc tính toán từ định luật Darcy và sóng động học có nghĩa là gradien thủy lực bằng với độ dốc ở mỗi ô: iiiiSi W/t)t(KSDc)t(RI 1 (2.16) Trong đó : RIi(t) (mm) là lƣợng dòng sát mặt chảy ra từ ô thứ i ở mỗi bƣớc thời gian (h), Di là độ sâu rễ ở ô thứ i (m), Si là độ dốc ở ô thứ i (m/m), Wi là chiều rộng của ô thứ i (m), CS là hệ số tỉ lệ phụ thuộc vào thảm phủ. 39  Lượng trữ nước ngầm và dòng chảy cơ sở Lƣợng trữ nƣớc ngầm là lƣợng nƣớc trong đới bão hoà. Thông thƣờng lƣu lƣợng dòng ngầm hình thành một dòng chảy cơ sở ở cửa ra của lƣu vực. m sgs tSGctQG 1000/)()( (2.17) Trong đó : QGs(t) dòng ngầm trung bình ở cửa ra của lƣu vực con (m 3 /s), SGs(t) là lƣợng trữ nƣớc ngầm ở lƣu vực con ở thời điểm t (mm), m là số mũ: m=1 cho hồ chứa tuyến tính và m = 2 cho hồ chứa phi tuyến, cg là hệ số triết giảm dòng ngầm phụ thuộc vào diện tích, hình dạng, thể tích của lỗ hổng và khả năng truyền của lƣu vực con. Cho mỗi tiểu lƣu vực, cân bằng lƣợng ngầm dƣới dạng: A ttQG tEG A AtRG tSGtSG ss s N i ii ss s 1000 )( )( )( )1()( 1 (2.18) Trong đó: SGs(t) và SGs(t-1) là lƣợng nƣớc ngầm của lƣu vực con ở bƣớc thời gian t và t-1 (mm), Ns là số lƣợng ô trong lƣu vực con, Ai là diện tích ô (m 2 ), As là diện tích tiểu lƣu vực (m2), EGs(t) là lƣợng bốc thoát hơi trung bình từ lƣợng nƣớc ngầm của lƣu vực con (mm), QGi(t) là lƣu lƣợng dòng ngầm (m 3 /s). Một phƣơng trình tuyến tính sử dụng trong mô hình liên quan giữa bốc thoát hơi dƣới tầng sâu với PET và lƣợng nƣớc ngầm: (t)SEtEDtEIEPcctEG iiivdi )()()( (2.19) Trong đó : EGi(t) là lƣợng bốc thoát hơi trung bình từ lƣợng nƣớc ngầm (mm), EP là PET (m), Cd là một biến, tính toán từ SGi(t)/SGs,o với SGs,o là lƣợng nƣớc ngầm của lƣu vực con ở thời điểm t (mm) và SGs,o là khả năng trữ nƣớc ngầm của lƣu vực con (mm).  Dòng chảy tràn và diễn toán trong kênh Trong WetSpa cải tiến, diễn toán dòng chảy tràn và dòng chảy trong kênh dùng phƣơng pháp sóng khuếch tán tuyến tính. Phƣơng pháp này phù hợp mô 40 phỏng dòng chảy ở mức độ nhất định và một trong những thuận lợi là nó có thể giải quyết theo phƣơng pháp giải tích, tránh tính toán bằng phƣơng pháp số và xác định các điều kiên biên một cách chính xác. Giả định ô lƣới là một đoạn sông với dòng chảy không ổn định một chiều và bỏ qua các số hạng chuyển động trong phƣơng trình động lƣợng St. Venant, quá trình chảy trong ô có thể đƣợc mô phỏng bởi phƣơng trình sóng khuếch tán (Miller và Cunge, 1975 ): 0 2 2 x Q d x Q c t Q ii (2.20) Q là lƣu lƣợng ở thời điểm t (s) và vị trí x (m), ci là tốc độ sóng khuếch tán ở ô thứ i(m/s), di là hệ số khuếch tán của sóng ở ô thứ i(m 2 /s). Xem xét một hệ thống đƣợc giới hạn bởi biên trên và biên dƣới, giải phƣơng trình (3.20) ở cửa ra của ô lƣới, khi vận tốc dòng chảy và hệ số khuếch tán không đổi có thể giải bằng quá trình phân bố mật độ thời gian lần chuyển tiếp đầu tiên của một chuyển động Brown: 2 3 ( ) ( ) exp 42 i i i i ii l c t l u t d td t (2.21) Với ui(t) là hàm phản ứng xung của ô lƣới (1/s) và li là kích cỡ ô lƣới (m). Hai tham số ci và di đƣợc tính toán theo công thức Manning (Henderson, 1966) : ii vc 3 5 (2.22) i ii i S Rv d 2 (2.23) Trong đó : Ri là bán kính thủy lực ở ô thứ i (m), Si là độ dốc ở ô thứ i (m/m), Vi là vận tốc dòng chảy (m/s). Bán kính thủy lực đƣợc tính theo công thức: 41 pb pi AaR )( (2.24) Trong đó Ai là diện tích lƣu vực ở thƣợng lƣu (km 2 ), ap là hằng số, bp là số mũ theo tỉ lệ hình dạng, tất cả phụ thuộc vào tần suất lũ. Vận tốc dòng chảy tính toán theo phƣơng trình Manning: 2 1 1 i 3 2 i i i SR n v (2.25) với ni là hệ số nhám Manning phụ thuộc vào loại thảm phủ và đặc điểm của kênh. Dƣới giả định hệ thống diễn toán tuyến tính, hàm phản ứng xung ở cuối dòng chảy là kết quả từ một xung đơn vị đầu vào đến một ô riêng lẻ, có thể đƣợc tính toán mà không cần sự can thiệp đến các ô khác. Dọc theo hƣớng phản ứng của dòng chảy gồm xung di chuyển qua nhiều ô, mỗi ô có một hàm phản ứng cấp đơn vị khác nhau. Trong quá trình diễn toán này, đầu ra của bất kì ô lƣới nào lại trở thành đầu vào của ô lƣới kế tiếp và phân bổ đầu vào gốc đƣợc thay đổi liên tục bởi quá trình động lực trong các ô, đó chính là các hàm phản ứng xung. Phản ứng theo hƣớng dòng chảy là tổng của các phản ứng xung liên tiếp nhau. )()( 1 tutU N j ji (2.26) Trong đó: Ui(t) là hàm phản ứng theo hƣớng dòng chảy (s -1 ) Chỉ số i liên quan đến ô nơi bắt đầu tính là đầu vào, j là số chuỗi số và N là tổng số ô dọc theo hƣớng dòng chảy. Mô hình phƣơng trình sóng khuếch tán thoả mãn phƣơng trình (2.26) cho các ô lƣới, điều đó có nghĩa rằng nó cho phép khả năng phân tích theo chiều dọc. Bởi vìa các hàm phản ứng xung là bất biến theo thời gian, do đó phƣơng trình (2.26) cũng bất biến theo thời gian và do đó có một mối quan hệ tuyến tính giữa phản ứng theo hƣớng dòng chảy và lƣợng đầu vào.Giả định rằng hàm phản ứng theo hƣớng dòng chảy Ui(t) cũng là một phân bố thời gian di chuyển đầu tiên, De smedt [8] đƣa ra một cách giải gần đúng phƣơng trình (2.26) bằng 42 phƣơng pháp số liên quan đến lƣu lƣợng ở cuối của một hƣớng dòng chảy đến nơi bắt đầu của một hƣớng dòng chảy tiếp theo : ii i ii tt tt tt tU /2 )( exp /2 1 )( 2 2 332 (2.27) Trong đó : ti là thời gian chảy trung bình từ ô thứ i đến dòng chảy (s), 2 i là phƣơng sai của thời gian chảy (s2). Tham số ti , 2 i là các tham số phân bố theo không gian, và có thể tính đƣợc bằng tích phân chập dọc theo địa hình quyết định dòng chảy nhƣ là một hàm của tốc độ sóng và hệ số phân tán : và j N j j i l c t 1 1 (2.28) j N j j j i l c d 1 3 2 2 (2.29) Hàm phản ứng ở cuối dòng chảy đến một đầu vào tuỳ ý ở ô bắt đầu, có thể đƣợc tính toán bằng cách cộng dồn tổng lƣợng dòng chảy đầu vào bởi hàm phản ứng xung đơn vị. Từ cách nhìn nhận trên, điều này tƣơng đƣơng với việc phân tích đầu vào thành các hàm phản ứng xác định và cộng tất cả các hàm phản ứng thành một phản ứng đơn. Do đó, quá trình lƣu lƣợng từ một lƣợng đầu vào bất kì đƣợc tính theo công thức : )()()( 0 tUVtQ t t iii (2.30) Trong đó : Qi(t) la lƣu lƣợng ở nơi kết thúc của một hƣớng dòng chảy sinh bởi một lƣợng đầu vào bất kì trong ô thứ i; Ui(t- ) tƣơng đƣơng với đƣờng đơn vị tức thời (IUH) sử dụng trong quy ƣớc của ngành thủy văn và là thời gian trễ (s) ; Vi( ) là tổng lƣợng dòng chảy đầu vào ở ô thứ i và ở thời gian gồm có dòng chảy 43 mặt và dòng chảy sát mặt, cộng thêm dòng chảy ngầm nếu ô thứ i nằm ở cửa ra của tiểu lƣu vực. Xem xét quá trình phân tích thực tế trong một hệ thống diễn toán tuyến tính, hàm phản ứng của dòng chảy trên lƣu vực có thể đƣợc quyết định là tổng các hàm phản ứng của các nhân tố đƣợc đóng góp từ tất cả các ô lƣới. Do đó, hàm phản ứng ở cấp độ lƣu vực có thể đƣợc tính toán là : wN i i tQtQ 1 )()( (2.31) Nw là số các ô trên toàn bộ lƣu vực, Q(t) là lƣu lƣợng ở cửa ra của lƣu vực. 2.2 Phƣơng pháp ƣớc lƣợng bất định khả năng - GLUE Theo Yang, Reichert, Abbaspour, Xia, and Yang [17] các mô hình cấu trúc lƣu vực ngày càng đƣợc sử dụng nhiều nhằm hỗ trợ trong vấn đề quản lý nguồn nƣớc vì các mô hình này đã đƣợc hiệu chỉnh một cách cẩn thận và đã đƣợc phân tích tính bất định. Tính bất định xuất hiện bởi sai số trong quá trình chọn mẫu, điều kiện ban đầu, số liệu đầu vào, đầu ra, và độ chính xác của các thông số (Blasone, Vrught, Madsen, Rosbjerg, Robinson, và Zyvolosky, 2008). Beven và Benley (2001) đã đƣa chỉ ra rằng phƣơng pháp GLUE có thể xử lý hoàn toàn những sai số yếu tố này. 2.2.1 Cơ sở lý thuyết phương pháp GLUE Phƣơng pháp ƣớc lƣợng bất định khả năng – Generalized Likehood Uncertainty Estimate (GLUE), là phƣơng pháp đƣợc đề xuất bởi Beven [12] nhằm tính toán đến tính bất định của mô hình, đây là một biện pháp để khắc phục các sai số yếu tố và tính toán khả năng. Phƣơng pháp GLUE trở thành 1 công cụ để phân tích độ nhạy và tính toán ƣớc lƣợng bất định khả năng bằng cách sử dụng các kết quả mô phỏng Monter Carlo. Phƣơng pháp GLUE đã đƣợc nghiên cứu và sử dụng làm lý thuyết chung về ƣớc lƣợng bất định khả năng ứng dụng trong thủy văn nhƣ Beven (1993), Romanowicz et al. (1994, 1996), Freer et al. (1996), Franks et al. (1997) và Zak et al. (1997) đã trình bày. Blasone et al (2008), và Binley [17] đã cho thấy GLUE đƣợc ứng dụng rộng rãi để đánh giá độ bất định đối với các mô hình 44 nhƣ mô hình mƣa rào dòng chảy, mô hình xói mòn đất, mô hình nƣớc dƣới đất, mô hình lũ lụt và mô hình thủy văn phân phối. Beven và Freer (2001) đã kết luận rằng, phƣơng pháp Glue có tính đến các yếu tố bất định: bất định đầu vào, bất định về cấu trúc, bất định về bộ thông số,.. Xuất phát từ ý tƣởng cho rằng: luôn tồn tại sai số nào đó ở cấu trúc của mô hình và trong quá trình xác định bộ thông số cũng tồn tại sai số. Vì vậy, không thể khẳng định rằng các giá trị tham số của mô hình xác định đƣợc (dù tốt nhất) là tối ƣu. Do đó phƣơng pháp GLUE không tập trung xác định một bộ thông số tối ƣu, mà xác định một tập hợp bộ thông số khả năng. Phƣơng pháp này dựa trên mô phỏng Monter Carlo bằng cách sử dụng những bộ thông số đƣợc lựa chọn ngẫu nhiên. Mỗi một mẫu lựa chọn đƣợc so sánh với khoảng giá trị giới hạn đƣợc mặc định bởi mô hình. Những mẫu không nằm trong khoảng giới hạn đó bị loại bỏ xem nhƣ không có. Những mẫu đƣợc giữ lại kết hợp với chỉ tiêu khả năng dựa trên cơ sở đánh giá mức độ phù hợp của số liệu đối với mô hình. Điêu kiện duy nhất của chỉ tiêu là những giá trị không phù hợp xem nhƣ bằng 0 và giá trị tăng dần với mức độ phù hợp, cụ thể nhƣ chỉ tiêu dựa theo phƣơng pháp thống kê xác xuất truyền thống, ngoài ra cũng có thể đƣợc lựa chọn bởi một định nghĩa khác. Ngoài ra phƣơng pháp ƣớc lƣợng bất định khả năng còn có khả năng cập nhật các giá trị khi có bộ số liệu mới, và tính toán xác định giá trị các bộ dữ liệu mới. Nhƣ vậy phƣơng pháp GLUE đƣợc xác định bởi các bƣớc sau đây: - Lựa chọn hay định nghĩa một chỉ tiêu để đánh giá độ phù hợp - Xác định khoảng giá trị và hàm phân bố của các tham số - Thiết lập quy trình sử dụng chỉ tiêu đánh giá phù hợp để tính toán khoảng bất định - Thiết lập quy trình cập nhật độ phù hợp khi có thêm số liệu 45 - Đánh giá giá trị của chuỗi số liệu bổ sung đối với thay đổi giá trị khoảng bất định. 2.2.2 Xác định chỉ tiêu đánh giá độ phù hợp Đây là bƣớc đầu tiên để đƣa ra một chỉ tiêu đánh giá mức độ phù hợp cho bộ tham số. Chỉ tiêu đánh giá cho biết mức độ phù hợp của mô phỏng (với mỗi bộ tham số) so với thực tế. Chỉ tiêu đánh giá phải tuân thủ một số đặc điểm nhất định. Giá trị của chỉ tiêu nên bắt đầu từ giá trị 0 đối với tất cả các mô phỏng cho kết quả hòan toàn không phù hợp với thực tế và đơn điệu tăng khi mức độ phù hợp giữa kết quả mô phỏng và thực tế tăng. Đặc tính này có thể thỏa mãn bởi nhiều công thức, do đó ngƣời sử dụng mô hình có thể lựa chọn nhiều chỉ tiêu đánh giá phù hợp. Các nghiên cứu từ trƣớc đã sử dụng các chỉ tiêu phù hợp khác nhau, và chúng đều bao gồm hai thành phần: công thức xác định chỉ tiêu và giá trị ngƣỡng loại bỏ. - Trong phƣơng pháp GLUE thƣờng sử dụng chỉ tiêu Nash, đƣợc xác định bởi công thức sau: 2 1 2 ,1 )( )( 1 avej M j jji M j i QoQo QoQs NS (2.32) Trong đó: i = 1, 2, 3....,N là số lần mô phỏng NSi là chỉ số phù hợp của lần mô phỏng thứ i j = 1, 2, 3, ...., M là bƣớc của của mô phỏng Qsi,j là lƣu lƣợng tính toán của lần mô phỏng thứ i tại thời điểm của bƣớc thời gian j Qoj là lƣu lƣợng quan trắc tại bƣớc thời gian j Qoave là lƣu lƣợng trung bình quan trắc đƣợc Chỉ tiêu thứ hai là chỉ tiêu hiệu quả mô hình (ME) đã đƣợc sử dụng nhiều trong GLUE . Công thức xác định hiệu quả mô hình nhƣ sau: )exp( 2 0 2 i i WL (2.33) 46 Trong đó: i = 1, 2, ... N là số lần mô phỏng Li là khả năng mô phỏng của lần thứ i σi là phƣơng sai của số dƣ của lần mô phỏng thứ i σo là phƣơng sai của các giá trị quan trắc W là trọng số có thể điều chỉnh đƣợc Trong nghiên cứu này W có thể tăng từ giá trị 1, 5,...100. Theo Blassone (2008) thì với W = 5 là hợp lý đối với tính toán bất định. Theo Beven và Binley (1992) sử dụng chỉ tiêu đánh giá sai số phƣơng sai EV (error variance), đƣợc tính nhƣ sau: V iiL )( 2 (2.34) Trong đó: i = 1, 2, ... N là số lần thực hiện mô phỏng Li là độ phù hợp của mô phỏng thứ i σ2i là phƣơng sai của số dƣ của lần mô phỏng thứ i V là trọng số Trong phạm vi nghiên cứu này V tăng từ 1, 5, 10. Với giá trị V = 5 đƣợc chỉ ra là phù hợp với đánh giá độ bất định - Giá trị ngƣỡng loại bỏ: là một giá trị dùng để phân biệt các mô phỏng đƣợc chấp nhận và không phù hợp. Các mô phỏng không phù hợp có trị số các chỉ tiêu đánh giá bằng 0 và các mô phỏng này bị loại bỏ trong quá trình ƣớc lƣợng khoảng bất định. Trong thực tế giá trị ngƣỡng loại bỏ này thƣờng là một giá trị xác định của một chỉ tiêu đánh giá (ví dụ chỉ tiêu NS > 0.8). Với chỉ tiêu đánh giá thứ nhất, Andersen, Refsgaard, và Jensen (2001) đã sắp xếp mức độ mô phỏng từ kém, trung bình và tốt với giá trị giới hạn là 0,7: NS nhỏ hơn 0,7 đƣợc đánh giá là mô phỏng kém, từ 0,7 trở lên đƣợc đánh giá là trung bình đến tốt (NS càng cao càng tốt) [4]. 47 Với chỉ tiêu đánh giá thứ hai, giá trị ngƣỡng đƣợc đƣa ra dựa trên thử nghiệm của Lamb, Beven và Myrabo (1998) là 10% các mô phỏng cho giá trị tốt nhất. Theo Beven và Binley thì với chỉ tiêu đánh giá thứ ba không có giá trị giới hạn của chỉ tiêu, nghĩa là các mô phỏng đều đƣợc đƣa vào trong ƣớc lƣợng khoảng bất định. 2.2.3 Xác định khoảng giá trị và hàm phân bố của các thông số Trong phƣơng pháp GLUE, việc xác định khoảng giá trị của các tham số là cần thiết. Độ rộng của dải giá trị phải phù hợp. Nếu rộng quá sẽ dẫn đến những mô phỏng không cần thiết, ngƣợc lại sẽ bỏ qua nhiều giá trị của tham số. Trong bƣớc này cần phải chú ý: - Lựa chọn tham số: cần xem xét các tham số nào có ảnh hƣởng thực sự đối với điều kiện thực tế áp dụng. Có thể dùng phân tích độ nhạy để lựa chọn các tham số này. - Xác định khoảng giá trị và hàm phân bố các tham số: thông thƣờng các khoảng giá trị của các tham số đƣợc xác định từ các nghiên cứu trƣớc đây và đối với hàm phân bố thƣờng đƣợc lấy là hàm phân bố đều khi ta không biết nhiều về giá trị của chúng. - Phƣơng pháp chọn mẫu: có hai phƣơng pháp chính là chọn mẫu ngẫu nhiên (Monte Carlo) và chọn mẫu theo phƣơng pháp siêu lập phƣơng Latin (Latin Hypercube Sampling - LHS). 2.2.4 Thiết lập quy trình sử dụng chỉ tiêu đánh giá độ phù hợp để tính toán khoảng bất định Sau khi xác định chỉ tiêu đánh giá và khoảng giá trị ban đầu của tham số, sử dụng phƣơng pháp phân tích Monte Carlo tính toán với nhiều bộ thông số. Trong thực tế nếu thời gian mô phỏng dài (đối với các mô hình phân phối), ngƣời ta thƣờng sử dụng phƣơng pháp chọn LHS để tăng hiệu quả của quá trình tính toán. Thực chất LHS là phƣơng pháp Monte Carlo cải tiến. Trong luận văn này, phƣơng pháp LHS đƣợc sử dụng và số mẫu mô phỏng là 200. Với mỗi bộ thông số đƣợc tạo ra bởi phƣơng pháp LHS, mô hình WetSpa sẽ tính toán đƣợc lƣu lƣợng dòng chảy ra. Từ đó giá trị của chỉ tiêu đánh giá đã chọn đƣợc tính toán. Có thể sử dụng một trong ba chỉ tiêu đánh giá nêu trên. Trong nội dung thực hiện luân văn cả ba chỉ tiêu đánh giá đƣợc thử nghiệm. Do quy trình tính 48 toán chỉ áp dụng với một chỉ tiêu nên sau đây là mô tả các bƣớc tiếp theo khi sử dụng một chỉ tiêu duy nhất. Khi các giá trị của chỉ tiêu đánh giá đƣợc xác định, các mô phỏng đƣợc chấp nhận (NSi > 0.7) đƣợc giữ lại để tính toán khoảng bất định bằng việc sử dụng giá trị ngƣỡng loại bỏ. Giá trị của các chỉ tiêu của các mô phỏng đƣợc chấp nhận sau đó đƣợc biến đổi sao cho tổng của chúng bằng một đơn vị theo công thức sau: RLi = Li/ (L1 + L2 +… + Ln) (2.35) Trong đó: RLi là giá trị biến đổi chỉ số độ phù hợp của mô phỏng thứ i, Li giá trị của chỉ số phù hợp của mô phỏng thứ i, L1 and L2 tƣơng ứng là giá trị của các chỉ số phù hợp của các mô phỏng đƣợc chấp nhận thứ nhất và thứ hai LN là giá trị của chỉ số của mô phỏng cuối cùng đƣợc đánh giá là phù hợp khi dùng giá trị ngƣỡng loại bỏ. Ở mỗi bƣớc thời gian, giá trị lƣu lƣợng ứng với 5% và 95% của hàm phân bố lũy tích các chỉ số phù hợp đƣợc sử dụng làm khoảng bất định của giá trị dự báo. Lƣu lƣợng Q5% và Q95% đƣợc xác định bởi công thức (5). )(%% nnbnna nnbnna nnbn nnbn QQ CLCL CLCL QQ (2.36) Trong đó: Qn% là lƣu lƣợng tƣơng ứng với n% của hàm phân bố lũy tích các chỉ số phù hợp; CLn% là giá trị của chỉ số phù hợp tƣơng ứng với n% của hàm phân bố lũy tích các chỉ số phù hợp; CLnna, CLnnb tƣơng ứng là giá trị của các chỉ số phù hợp ngay trên và dƣới giá trị CLn%; Qnna và Qnnb tƣơng ứng là giá trị lƣu lƣợng ngay trên và dƣới giá trị Qn%. 49 2.2.5 Thiết lập quy trình cập nhật độ phù hợp khi có thêm số liệu Trong quy trình tính toán bất định theo phƣơng pháp GLUE, có thể cập nhật giá trị chỉ tiêu đánh giá độ phù hợp khi có dữ liệu mới. Sau đó các giá trị này có thể đƣợc cập nhật, sử dụng phƣơng trình Bayes: L( y) = L( |y)L( ) (2.37) Trong đó: L( ) là phân phối các chỉ tiêu phù hợp của tập các bộ tham số trƣớc khi cập nhật L( |y) là phân phối các chỉ tiêu phù hợp khi có các số liệu mới (trƣớc khi áp dụng ngƣỡng) L( y) là phân phối các chỉ tiêu phù hợp của tập các bộ tham số sau khi cập nhật (trƣớc khi áp dụng ngƣỡng loại bỏ) Phân phối các chỉ tiêu phù hợp đƣợc cập nhật sau đó đƣợc sử dụng để cập nhật khoảng bất định. Chú ý rằng khi sử dung các chỉ tiêu ME, EV, giá trị ngƣỡng loại bỏ có thể giữ nguyên. Còn khi sử dụng chỉ tiêu NS thì ngƣỡng giá trị lọa bỏ cần đƣợc tính lại theo công thức sau: NS > 0,7 n (2.38) với n là số tập số liệu đƣợc bổ sung. 2.2.6 Chế độ mô phỏng và chế độ dự báo Các quy trình tính toán và cập nhật chỉ tiêu phù hợp có thể sử dụng trong cả chế độ mô phỏng và chế độ dự báo. Đối với chế độ mô phỏng thì số liệu mƣa và dòng chảy đã có sẵn, và từ những dữ liệu này các bƣớc trong phƣơng pháp GLUE thực hiện tính toán ƣớc lƣợng và cập nhật bất định. Do số liệu về lƣu lƣợng thực đo đã có nên đƣờng quá trình thực đo cùng các khoảng giá trị bất định có thể cùng vẽ lên trên một biểu đồ để xem xét tính phù hợp. Đối với chế độ dự báo thì số liệu dòng chảy chƣa có. Số liệu mƣa đƣa vào quy trình là mƣa thiết kế hay mƣa dự báo, từ đó mô hình WetSpa sẽ sử dụng bộ thông số ban đầu để tính toán dòng chảy và sử dụng các chỉ số phù hợp đƣợc xác định trƣớc cộng với lũ mô phỏng để tính toán 50 khoảng bất định. Và nhƣ vậy, quy trình chỉ có thể áp dụng đƣợc ở chế độ dự báo sau khi đã áp dụng nó ở chế độ mô phỏng. 51 CHƢƠNG 3: ÁP DỤNG PHƢƠNG PHÁP ƢỚC LƢỢNG BẤT ĐỊNH CHO DỰ BÁO LŨ TRÊN SÔNG VỆ 3.1 Thu thập và xử lý dữ liệu 3.1.1 Số liệu không gian Trong mô hình WetSpa cải tiến sử dụng 3 bản đồ số là: DEM, bản đồ đất và sử dụng đất. Ngoài ra để so sánh và tính toán các đặc trƣng lƣu vực còn cần sử dụng đến bản đồ mạng lƣới sông suối và bản đồ mạng lƣới trạm thủy văn trên lƣu vực. Bản đồ DEM với kích thƣớc 50x50 m dùng để tính toán các tham số liên quan đến địa hình. Các loại thảm phủ trên lƣu vực sông Vệ đƣợc chuyển đổi sao cho phù hợp với các thuộc tính trong mô hình và đƣa về cùng kích cỡ ô lƣới giống DEM. Từ bản đồ này ta có đƣợc các tham số về hệ số dòng chảy tiềm năng và khả năng trữ của các vùng trũng. Các loại đất trên lƣu vực sông cũng đƣợc thay đổi cho phù hợp với các loại đất trong mô hình và đƣa về kích cỡ ô lƣới 50x50m. Từ bản đồ này, các tham số về khả năng về độ rỗng, độ dẫn thủy lực, độ ẩm dƣ... đƣợc đƣa vào mô hình. Hình 3.1: Bản đồ số độ cao DEM và bản đồ sử mạng lƣới trạm lƣu vực sông Vệ tính đến trạm An Chỉ 52 Hình 3.2: Bản đồ sử dụng đất và bản đồ đất lƣu vực sông Vệ tính đến trạm An Chỉ 3.1.2 Số liệu khí tượng thủy văn Số liệu mƣa tại 4 trạm An Chỉ, Sơn Giang, Giá Vực, Ba Tơ đƣợc sử dụng để tính toán dòng chảy trên lƣu vực. Trong đó, có hai trạm Ba Tơ và An Chỉ là trạm nằm trong lƣu vực, hai trạm đo còn lại Sơn Giang và Giá Vực nằm ngoài lƣu vực, đƣợc sử dụng để vẽ đa giác Thiessen và nội suy số liệu mƣa trên toàn bộ lƣu vực. Số liệu mƣa quan trắc 6 giờ một lần đƣợc quy về mƣa hàng giờ. 3.2 Tính toán trong Arcview Sử dụng phần mềm Arcview tính toán nội suy các bản đồ thô thủy văn, thủy lực tại từng ô lƣới trên lƣu vực: - Tính toán hƣớng dòng chảy và tích tụ dòng chảy - Tính toán mạng lƣới sông suối, với kích thƣớc ô lƣới ban đầu cấp sông suối - Tính toán độ dốc với giá trị độ dốc nhỏ nhất là 0.01% - Tính toán bán kính thủy lực cho từng ô lƣới ứng với tần suất l - Phân chia lƣu vực thành 13 lƣu vực con với giá trị ô lƣới là 40 - Tính toán độ dẫn thủy lực - Tính toán độ rỗng của đất - Tính toán khả năng trữ - Tính toán lƣợng ẩm dƣ - Tính toán chỉ số phân bố kích cỡ độ rỗng của đất - Tính toán giai đoạn héo của thực vật 53 - Tính toán độ sâu của đới rễ cây - Tính toán khả năng ngƣng tụ - Tính toán hệ số Manning sử dụng bảng tra của mô hình - Tính toán hệ số dòng chảy tiềm năng với giả thiết bề mặt không thấm là 70% - Tính toán khả năng trữ - Tính toán vận tốc dòng chảy và thời gian chảy trong ô lƣới và lƣu vực con - Phân chia đa giác Thiessen từ các trạm mƣa đã có Các bản đồ thông số này đƣợc xuất ra dƣới định dạng ASCII phù hợp với đầu vào của ngôn ngữ lập trình Fortran. 3.3 Các thông số toàn cục của mô hình A. Bahremand và F. De Smedt (2007) đã chỉ ra 11 thông số chính có thể hiệu chỉnh đƣợc trong mô hình WetSpa, các thông số còn lại đƣợc tính toán tự động trong Arcview và không thể thay đổi đƣợc. 11 thông số này cùng cũng đồng thời là 11 trong số 12 thông số toàn cục có trong mô hình, bao gồm: hệ số dòng sát mặt, hệ số dòng ngầm, độ ẩm đất, thông số hiệu chỉnh bốc thoát hơi nƣớc khả năng, lƣợng trữ nƣớc ngầm ban đầu, lƣợng trữ nƣớc ngầm cực đại, nhiệt độ tuyết tan, hệ số nhiệt độ, hệ số mƣa, hệ số dòng chảy mặt, cƣờng độ mƣa tƣơng ứng với hệ số dòng chảy mặt bằng 1, và một thông số nữa không thể đƣa vào phân tích độ nhạy là thời gian. Những tham số này mang tính chất vật lý quan trọng trong kiểm soát quá trình sản sinh dòng chảy và lƣu lƣợng ở cửa ra lƣu vực, nhƣng lại rất khó xác định chính xác trên từng ô lƣới. Do đó, việc hiệu chỉnh những tham số toàn cục này dựa vào số liệu dòng chảy thực đo là cần thiết đối với mô hình phân phối này. * Thời gian tính toán dt(h) Trong mô hình, ta có thể sử dụng chuỗi số liệu tính toán theo bƣớc thời gian giờ hay ngày. Giá trị của dt phụ thuộc vào chuỗi số liệu đo đạc đƣợc, dt bằng 1 đối với chuỗi quan trắc hàng giờ, dt bằng 24 đối với chuỗi quan trắc hàng ngày. Thông số này không thể hiệu chỉnh đƣợc. * Thông số dòng sát mặt Ci 54 Dòng sát mặt là thành phần quan trọng của cân bằng nƣớc trong đất. Nó là lƣợng nƣớc thấm xuống lớp đất mặt và di chuyển theo phƣơng ngang đến khi gia nhập vào dòng chính. Đây là thành phần dòng chảy chủ yếu ở vùng nhiệt đới ẩm, đặc biệt là ở các vùng dốc và có độ che phủ tốt. Trong mô hình WetSpa cải tiến, dòng sát mặt đƣợc giả định xảy ra sau quá trình ngấm khi độ ẩm đất vƣợt quá khả năng chứa nƣớc và gradient thủy lực đủ lớn để làm nƣớc chảy, và ngừng khi lƣợng ẩm của đất thấp hơn khả năng trữ. Hệ số này góp phần quyết định dòng chảy ra của lƣu vực. Theo A. Bahremand và F. De Smedt (2007) Ci thƣờng có giá trị lớn hơn 1. * Hệ số triết giảm dòng ngầm Cg Cg phản ánh chế độ triết giảm nƣớc ngầm cho toàn bộ lƣu vực. Cg là hệ số triết giảm dòng ngầm phụ thuộc vào diện tích, hình dạng, thể tích của lỗ hổng và khả năng truyền nƣớc của lƣu vực con. Để đơn giản mô hình, một giá trị chung cho hệ số triết giảm dòng ngầm đƣợc xác định tại cửa ra lƣu vực trong 1 file đầu vào. Hệ số này đặc trƣng cho mỗi lƣu vực con dựa vào diện tích thoát nƣớc và độ dốc trung bình, trong đó giá trị cao hơn đƣợc ấn định cho lƣu vực con có diện tích thoát nƣớc lớn và độ dốc lớn, và giá trị thấp hơn cho lƣu vực con có diện tích thoát nƣớc nhỏ và độ dốc thoải. * Tiêu chuẩn độ ẩm đất Kss Kss liên quan tới sức chứa tối đa và độ ẩm đất ban đầu. Độ ẩm đất là thành phần then chốt trong mô hình điều khiển các quá trình thủy văn sản sinh dòng chảy mặt, bốc thoát hơi nƣớc, thấm và dòng sát mặt. * Thông số hiệu chỉnh bốc thoát hơi nƣớc khả năng Kep Tốc độ bốc thoát hơi nƣớc khả năng liên quan đến bề mặt nƣớc hay cỏ bao phủ trên diện rộng. Để tính đến những ảnh hƣởng này, cần sử dụng một hệ số hiệu chỉnh, hệ số này thƣờng gần bằng 1, thƣờng đƣợc đƣa vào mô hình dƣới dạng một hằng số giả định cho toàn lƣu vực trong suốt chuỗi thời gian và có thể đƣợc hiệu chỉnh bởi mô hình thông qua mô phỏng cân bằng nƣớc hạn dài. * Lƣợng trữ nƣớc ngầm ban đầu G0 55 G0 chính là độ sâu tầng nƣớc ngầm (mm). Trong WetSpa cải tiến, cân bằng nƣớc ngầm đƣợc duy trì trên quy mô lƣu vực con và coi lƣợng trữ nƣớc ngầm hiệu quả là một phần của lƣợng trữ trong tầng ngậm nƣớc góp phần vào dòng chảy mặt. Giá trị này có thể đƣợc điều chỉnh trong suốt quá trình hiệu chỉnh bằng cách so sánh giữa dòng chảy tính toán và thực đo cho giai đoạn ban đầu, là một trong những điều kiện ban đầu quan trọng trong mô hình. * Lƣợng trữ nƣớc ngầm cực đại Gmax Lƣợng trữ nƣớc ngầm cực đại là độ sâu cực đại của tầng nƣớc ngầm trên lƣu vực. Mặc dù là thông số rất quan trọng nhƣng Gmax thƣờng không thể đo đạc đƣợc trực tiếp mà phải giả định mang một giá trị nào đó, giá trị này có thể đƣợc thay đổi trong quá trình hiệu chỉnh mô hình. * T0 T0 là giá trị nhiệt độ cơ sở (°C) mà tại đó giáng thủy thay đổi trạng thái từ mƣa thành tuyết, giá trị điển hình thƣờng bằng 0 hoặc rất gần với 0. * Ksnow Ksnow là hệ số tan chảy tuyết theo nhiệt độ hàng ngày (mm/°C/ngày). Phạm vi hệ số nhiệt độ ngày điển hình là 1.8 - 3.7mm/0C/ngày đối với mƣa tự do. Thông thƣờng, hệ số nhiệt độ ngày thay đổi cả theo thời gian và không gian, nhƣng trong mô hình WetSpa nó đƣợc coi là một thông số toàn cục và mang một giá trị cố địnhcho toàn bộ lƣu vực. * K_rain Hệ số nhiệt độ - mƣa ngày (mm/mm/°C/day) quyết định tốc độ tuyết tan gây ra bởi sự ngƣng tụ của không khí ẩm trên bề mặt tuyết và sự truyền nhiệt bình lƣu cho tuyết thông qua giáng thủy, và đƣợc sử dụng để tính toán tuyết tan thêm vào lƣợng mƣa rơi. Giá trị hệ số mƣa ngày thƣờng rất nhỏ, điển hình khoảng 0.01 *Thông số dòng chảy mặt Krun Krun là hệ số dòng chảy mặt đối với cƣờng độ mƣa rất nhỏ gần bằng 0. Đây là một hệ số quan trọng, phản ánh dòng chảy mặt cơ sở trên lƣu vực. 56 *Cƣờng độ mƣa tƣơng ứng với thành phần dòng chảy mặt bằng 1 Pmax Thông số này là cƣờng độ mƣa giới hạn (hay chính là cƣờng độ mƣa vƣợt ngƣỡng) với đơn vị là mm/h hay mm/ngày phụ thuộc vào bƣớc thời gian của quá trình mô phỏng, với hệ số dòng chảy mặt bằng 1. Trên thực tế giá trị 2 thông số này thay đổi theo không gian, phụ thuộc vào đặc tính từng ô lƣới, nhƣ loại đất, sử dụng đất, độ dốc..., để đơn giản trong mô hình giả thiết Krun và Pmax là hằng số. Việc hiệu chỉnh 2 thông số này có thể đƣợc thực hiện bằng cách so sánh lƣợng dòng chảy mặt và đỉnh lũ tính toán và thực đo. 3.4 Xây dựng quy trình dự báo lũ có tính đến độ bất định của bộ thông số Trong nội dung này trình bày các bƣớc trong quy trình đƣợc thực hiện bằng cách lập trình trong ngôn ngữ Matlab và Fortran. Trình tự nhƣ sơ đồ sau: Hình 3.3: Quy trình tính toán bất định 57 3.4.1 Lựa chọn thông số Các tham số trong mô hình WetSpa đƣợc chia thành hai phần: các thông số trong quá trình thiết lập thời gian ArcView và các thông số toàn cầu (đã trình bày ở mục trên). Không phân tích tính toán ƣớc lƣợng bất định với các thông số trong ArcView do các thiết lập trong ArcView thực hiện một cách thủ công, do đó với khối lƣợng mô phỏng lớn thì việc phân tích tính bất định là không khả thi. Vì vậy, xem xét 12 tham số toàn cầu xác định đƣợc 7 tham số đƣợc đƣa vào tính toán bất định. Bƣớc thời gian là tham số đầu tiên đó không đƣa vào tính toán do bƣớc thời gian mặc định trong dự báo lũ là 1 giờ. Ba thông số, T0, K_snow và K_rain, chỉ đƣợc sử dụng khi có tuyết tan. Đối với lƣu vực sông Vệ, không có hiện tƣợng tuyết tan nên các tham số này không ảnh hƣởng đến lƣu lƣợng lũ do đó không xét đến. Tham số thứ năm không đƣợc vào tính toán là K_ep - một yếu tố liên quan đến bốc thoát hơi. Trên thực tế, trong thời gian lũ, lƣợng bốc hơi xem nhƣ bằng 0. Tham số K_ep không có thể ảnh hƣởng kết quả đầu ra mô hình vì thế không xét đến đối với phân tích bất. Vì vậy, chỉ có 7 tham số đƣợc đƣa vào tính toán bất định trong nghiên cứu này. 3.4.2 Khoảng bất định của các thông số Trong quá trình xác định các tham số của mô hình cho thấy có một số bộ tham số khác nhau cho kết quả giống nhau. Từ đó cho thấy trong một giới hạn nào đó của bộ tham số có thể cho kết quả phù hợp. Vì vậy để phân tích tính bất định cần phải xác định khoảng bất bất định của các tham số. Xác định khoảng bất bất định của các tham số đƣợc dựa trên bản chất vật lý đặc trƣng của tham số trong mô hình. Nó phụ thuộc hoàn toàn vào chủ quan của ngƣời sử dụng mô hình. Trong trƣờng hợp kinh nghiệm bị hạn chế thì khoảng bất định đƣợc lựa chọn một cách thống nhất từ các nghiên cứu đã có là phù hợp (Beven và Binley, 1992) Để xác định khoảng bất định của các thông số dựa trên nghiên cứu về phân tích độ nhạy của các tham số. Việc phân tích độ nhạy các thông số trong mô hình WetSpa cải tiến đã từng đƣợc thực hiện bởi Bahremand và De Smedt (2008) bằng phƣơng pháp PEST. Hai tác giả này đề xuất để nghiên cứu sâu hơn, nên sử dụng 58 một phƣơng pháp tổng thể để phân tích độ nhạy. Kết quả phân tích 8 thông số còn lại của Bahremand và De Smedt (2008) [1] bằng phƣơng pháp PEST tính toán đƣợc độ nhạy của chúng nhƣ trong bảng và khẳng định mối quan hệ nghịch biến giữa Kss và G0. Bảng 3.1: Phân tích độ nhạy của các thông số theo phƣơng pháp PEST Tuy nhiên phƣơng pháp cục bộ, nhƣ PEST, chỉ sử dụng các dữ liệu cục bộ không phản ánh đƣợc toàn bộ không gian thông số, do đó kết quả chỉ đúng với độ nhạy và độ bất định cục bộ [1]. Trong luận văn này sử dụng kết quả của Tom Doldersum [9] sử dụng phƣơng pháp Morris để phân tích độ nhạy tổng thể cho các thông số trong mô hình WetSpa cải tiến, trong đó có xét đến sự khác biệt khi áp dụng đối với các lƣu vực ở Việt Nam và Châu Âu. Kết quả nhƣ trong bảng sau: Bảng 3.2: Khoảng bất định của các thông số Thông số Ki Kg Kss G0 Gmax Krun Pmax Giá trị nhỏ nhất 2 0.002 0 0 50 0 0 Giá trị lớn nhất 11 0.06 1.5 50 150 10 500 3.4.3 Phương pháp lấy mẫu Nhƣ trong nội dung cơ sở lý thuyết đã trình bày, có hai phƣơng pháp lấy mẫu: lấy mẫu ngẫu nhiên và lấy mẫu siêu lập phƣơng LHS. Lấy mẫu ngẫu nhiên tạo ra một số lƣợng lớn các bộ thông số của mô hình một cách ngẫu nhiên theo xác suất (Saltelli, Chan and Scott, 2000) và các mẫu tạo ra sẽ đƣợc tính toán. Phƣơng pháp này phù hợp với những mô hình bán phân phối, thời gian tính toán ngắn. Tuy nhiên, mô hình phân phối đầy đủ-có thời gian tính toán lớn hơn, do đó phƣơng pháp ngẫu nhiên lấy mẫu không hiệu quả. Vì vậy cần phải có phƣơng pháp lấy mẫu phù hợp với mô hình phân phối đầy đủ (theo lập luận của Uhlenbrook và Sieber (2005)). Uhlenbrook và Sieber (2005) đã sử dụng phƣơng pháp LHS, phù hợp với mô hình phân phối đầy đủ và tính toán trong phƣơng pháp GLUE. 59 LHS là một phƣơng pháp lấy mẫu phân tầng mà hiệu quả ƣớc lƣợng số liệu thống kê khoảng giá trị. Các phân phối xác suất của mỗi tham số đƣợc chia thành N dải giá với một xác suất bằng nhau xảy ra (1 / N). Mỗi một giá trị mô phỏng đƣợc lấy ngẫu nhiên 1 lần trong một dải giá trị. Trình tự lựa chọn các dãy là ngẫu nhiên và mô hình đƣợc thực hiện với một tổ hợp ngẫu nhiên của các giá trị tham số từ khoảng xác định mỗi dải giá trị giới hạn của các tham số (Yu, Yang, and Chen, 2001) Trong nội dung thực hiện của luận văn, N đƣợc đặt tới 5 và số lƣợng mẫu lựa chọn là 200 bộ thông số trên dải phân bố của các giá trị. Hình 3.4: File 200 bộ thông số đƣợc chọn theo phƣơng pháp LHS 3.4.4 Tính toán với mô hình WetSpa Với các mẫu đã lựa chọn ở trên tính toán với mô hình WetSpa với mỗi số liệu đầu vào (số liệu khí tƣợng) tính toán lƣu lƣợng dòng chảy (output) tƣơng ứng 60 với 200 bộ thông số sau khi thực hiện thuật toán tính lặp 200 lần bằng ngôn ngữ Fortran trong mã nguồn của mô hình WetSpa. Thời gian xử lý, tính toán đối với 200 mẫu khoảng 2 giờ. Để tính toán với phƣơng pháp Glue cần phải thực hiện mô phỏng trƣớc khi tiến hành dự báo. Sử dụng trực tiếp trận lũ 11/99 và 10/2003 để thực hiện mô phỏng. Từ kết quả thu đƣợc là 200 đƣờng quá trình lũ có thể thực hiện tính toán ƣớc lƣợng bất định theo các trình tự đã nêu trong chƣơng 2. Hình 3.5: Kết quả tính toán lƣu lƣợng lũ tháng 10/2003 tƣơng ứng với 200 bộ thông số 3.4.5 Lựa chọn chỉ tiêu Trong luận văn này, có 3 chỉ tiêu đƣợc đƣa ra để đánh giá độ phù hợp của mẫu (hay bộ thông số) đã lựa chọn ở trên. Trong bƣớc mô phỏng tiến hành đánh giá với cả ba chỉ tiêu: NS, ME và EV. Từ đó phân tích lựa chọn một chỉ tiêu dùng cho bƣớc dự báo. Sử dụng các chƣơng trình đƣợc viết bằng ngôn ngữ Matlab thực hiện tính toán đánh giá bằng các chỉ tiêu. 61 3.4.6 Tính toán khả năng Dựa vào kết quả tính toán theo các chỉ tiêu đã thực hiện ở trên bƣớc này tính toán thƣớc đo khả năng của mẫu mô phỏng. Các mẫu mô phỏng đƣợc giữ lai là mẫu có kết quả tính toán theo các chỉ tiêu vƣợt qua ngƣỡng loại bỏ đối với mỗi loại chỉ tiêu. Trong bƣớc tính toán này có tính quyết định đối với tính toán bất định. Do số lƣợng mẫu tính toán là 200, nên trong phần trình bày sẽ đƣa một ví dụ thu nhỏ về kết quả tính toán khả năng nhƣ sau: Bảng 3.4: Kết quả tính toán khả năng Hàng/ cột Các bộ thông số Thực đo 1 2 3 .. .. .. .. Tinh toán chỉ tiêu 1 0.0519 7 0.7306 7 0 0.1456 0.212 4 0.8457 7 0.5804 8 Tính toán khả năng 1 0 0.7306 7 0 0 0 0.8457 7 0 Giá trị mẫu (bộ thông số) 0 5.043 2.077 3.917 4.755 5.79 3.252 3.627 0 0.042 0.004 0.063 0.069 0.06 0.005 0.013 0 1.207 0.06 1.011 0.433 0.213 0.55 0.962 0 25.004 7.204 42.343 12.037 14.02 2 29.798 35.562 0 83.418 54.12 95.933 89.324 70.46 5 117.95 9 90.131 0 0 0 0 0 0.001 0 0 0 153.70 8 26.497 341.99 1 266.53 7 495.4 3 21.456 375.83 3.4.7 Tính toán bất định (UE) Bao gồm hai chế độ: chế độ mô phỏng và dự báo. a) Mô phỏng Đối với chế độ mô phỏng, số liệu cần để tính toán là lƣu lƣợng thực đo, lƣu lƣợng tính toán từ mô hình WetSpa và tính toán khả năng từ chỉ tiêu đƣợc lựa chọn để đánh giá. Kết quả của tính toán bất định là những chuỗi lƣu lƣợng có các giá trị đánh giá theo các chỉ tiêu là đạt (chỉ tiêu so sánh giữa tính toán và thực đo, giá trị căn cứ theo ngƣỡng loại bỏ). Từ đó tính toán đƣợc dải giá trị lƣu lƣợng 62 phù hợp mà giá trị thực đo nằm trong khoảng giới hạn đó. Luận văn thực hiện mô phỏng với 2 trận lũ: tháng 11/1999 và tháng 10/2003 Hình 3.6: Kết quả tính toán bất định chế độ mô phỏng với trận lũ tháng 10/2003(chỉ tiêu NS) 63 Hình 3.7: Kết quả tính toán bất định chế độ mô phỏng với trận lũ tháng 10/2003 (chỉ tiêu EV) Hình 3.8: Kết quả tính toán bất định chế độ mô phỏng với trận lũ tháng 10/2003 (chỉ tiêu ME) 64 Hình 3.9: Kết quả tính toán bất định chế độ mô phỏng với trận lũ tháng 11/1999 (chỉ tiêu NS) Hình 3.10: Kết quả tính toán bất định chế độ mô phỏng với trận lũ tháng 11/1999 (chỉ tiêu EV) 65 Hình 3.11: Kết quả tính toán bất định chế độ mô phỏng với trận lũ tháng 11/1999 (chỉ tiêu ME) b) Dự báo Đối với chế độ dự báo, số liệu lƣu lƣợng thực đo là chƣa biết, lƣợng mƣa có thể sử dụng là lƣợng mƣa dự báo. Do đó, sử dụng mô hình WetSpa để tính toán các chuỗi lƣu lƣợng tƣơng ứng với bộ thông số từ số liệu mƣa. Trận lũ mô phỏng này sẽ đƣợc tính toán khoảng tất định dựa trên chỉ số phù hợp cũ thu đƣợc từ các mô phỏng trƣớc. Các quá trình tính toán nhƣ đã trình bày nhƣng ở bƣớc tính toán bất định lựa chọn chế độ mô phỏng. Trong nội dung luận văn đã tính toán dự báo đối với trận lũ tháng 12/99, số liệu mƣa là số liệu thực đo và lƣu lƣợng đƣợc xem nhƣ chƣa biết. Lựa chọn chỉ tiêu NS và tính toán theo trình tự thu đƣợc kết quả nhƣ sau: 66 Hình 3.12: Kết quả dự báo trận lũ tháng 11/1999 (Đƣờng liền nét mảnh là đƣờng thực đo) 67 KẾT LUẬN VÀ KIẾN NGHỊ Công tác dự báo lũ đóng một vai trò quan trọng trong việc ổn định phát định phát triển kinh tế xã hội. Dự báo lũ phục vụ cho mục đích phòng chống giảm nhẹ các thiên tai do nƣớc gây ra và phục vụ cho công tác vận hành các công trình khai thác nguồn nƣớc. Vì vậy, việc đƣa ra một kết quả dự báo chính xác có ý nghĩa đối với cộng đồng và mọi ngành kinh tế. Nghiên cứu phƣơng pháp ƣớc lƣợng bất định là một hƣớng đi mới trong nghiên cứu công nghệ dự báo lũ. Trên thế giới hiện nay đã có một số nghiên cứu về ƣớc lƣợng bất định nhƣng nghiên cứu mới chỉ mang tính tham khảo chƣa đƣợc đƣa vào ứng dụng. Do đó, vấn đề phân tích tính toán độ bất định và khoảng dự báo trong dự báo lũ là vấn đề khá mới mẻ ở Việt Nam. Luận văn này đã nghiên cứu cơ sở lý thuyết của phƣơng pháp ƣớc lƣợng bất định GLUE. Từ đó, kết hợp với mô hình WetSpa xây dựng chƣơng trình tính toán khoảng dự báo bằng phần mềm Matlab. Kết quả là sơ đồ quá trình tính toán khoảng dự báo đƣợc thiết lập, bao gồm các thủ tục con dễ hiểu, dễ hiệu chỉnh đầu vào và kết quả áp dụng thử nghiệm cho dự báo lũ trên sông Vệ Tính đến trạm An Chỉ. Trong quá trình thực hiện đề tài, học viên rút ra một số vấn đề sau: 1) Về mô hình WetsSpa Mô hình tính toán với chuỗi số liệu đầu vào liên tục. Do đó, trong giai đoạn chuẩn bị dữ liệu phải thực hiện kiểm tra tính liên tục và độ tin cậy của số liệu. Các giá trị âm trong chuỗi số liệu đại diện cho trƣờng hợp các dữ liệu thực đo bị thiếu phải đƣợc thay thế bằng các giá trị nội suy. Cách phân chia các loại sử dụng đất không rõ ràng gây khó khăn cho ngƣời sử dụng. Các giá trị đƣợc gán cho mỗi ô lƣới biểu hiện giá trị trung bình trên diện tích mỗi ô. Sự biến thiên trên mỗi ô lƣới càng lớn, sai số sẽ càng tăng. Do đó, kích thƣớc ô lƣới nên đƣợc xác định rõ ràng. Kích cỡ ô lƣới nhỏ có thể biểu hiện tốt hơn sự thay đổi các đặc điểm vật lý trên lƣu vực, nhƣng dẫn đến việc giả định thời gian và tốn bộ nhớ hơn trong suốt thời gian mô phỏng, đặc biệt cho những lƣu vực lớn. - Bƣớc thời gian trong mô hình là ngày hoặc giờ sẽ không khả thi khi dự báo lũ cho một lƣu vực rất nhỏ, nơi lƣợng nƣớc thừa có thể chảy ra ngoài ngay ở bƣớc thời gian đầu tiên. Phần diện tích không thấm ở khu vực đô thị đƣợc đƣa vào mô hình một cách chủ quan, phụ thuộc vào kích cỡ ô lƣới. 68 - Trong một ô lƣới kích cỡ 50x50 m thì 30% diện tích không thấm đƣợc gán vào khu vực dân cƣ, 70% cho khu vực công nghiệp và thƣơng mại và 100% cho bãi đỗ xe, đƣờng chính… Điều này không phản ánh thực tế và mang đến nhiều sai số cho kết quả mô hình. - Mô hình sử dụng nhiều hệ số kinh nghiệm đƣợc mặc định qua nội suy và hiệu chỉnh từ các nghiên cứu trƣớc đây và sử dụng cho toàn bộ lƣu vực. Do phạm vi dao động quá lớn, nhiều tham số nhƣ độ dẫn thuỷ lực, hệ số nhám…có thể thay đổi lớn khi ứng dụng mô hình đến những địa điểm khác với môi trƣờng hoàn toàn khác. Do đó, việc hiệu chỉnh mô hình là cần thiết và điều này mang đến những khó khăn cho quá trình tham số hóa của mô hình ở lƣu vực không có trạm đo. WetSpa sử dụng nhiều loại ngôn ngữ lập trình phức tạp nhƣ ArcView Avenue, Fortran và Visual Basic, gây khó khăn cho ngƣời dùng khi muốn thay đổi mô hình cho phù hợp với nhu cầu sử dụng. 2) Về phƣơng pháp GLUE: nghiên cứu phƣơng pháp ƣớc lƣợng bất định GLUE nghiên cứu mới trong công nghệ dự báo lũ. Trong đó đã tính toán, lƣợng hóa các yếu tố bất định và thể hiện bằng kết quả đƣa ra là một dải giá trị thay vì một giá trị duy nhất theo truyền thống. Dải giá trị đƣợc đƣa ra từ những tính toán với nhiều bộ thông số của mô hình WetSpa có khả năng cho kết quả mô phỏng, dự báo tƣơng đối tốt khác với những phƣơng pháp truyền thống (chỉ tính toán với 1 bộ thông số). Trong luận văn, phƣơng pháp GLUE đã đƣợc khai triển theo các quy trình đã trình bày ở chƣơng 2 để áp dụng cho dự báo lũ sử dụng mô hình WetSpa. Ngôn ngữ lập trình cấp cao Matlab có hệ thống tính toán khoa học kỹ thuật, thƣ viện hàm phong phú có khả năng mô phỏng, vẽ đồ thị, biểu đồ và phân tích dữ liệu … Vì vậy, sử dụng ngôn ngữ Matlab để thực hiện các bƣớc trong phƣơng pháp GLUE, bao gồm các bƣớc chính nhƣ sau: chọn mẫu LHS; các chỉ tiêu NS, ME, EV: từ kết quả mô phỏng của mô hình sẽ thu đƣợc giá trị các chỉ tiêu phù hợp từ các chỉ tiêu đánh giá; các mô phỏng đƣợc chấp nhận; tính toán khoảng bất định. Trong đó, lựa chọn mẫu LHS và xác định các chỉ tiêu đƣợc đánh giá là quan trọng nhất. Đối với lựa chọn mẫu LHS, vấn đề xác định các thông số và khoảng phân bố của các thông số ảnh hƣởng trực tiếp đến kết quả dự báo. Trong nội dung luận văn này đã tham khảo một số nghiên cứu và lựa chọn đƣợc 7 thông số chính và khoảng phân bố của các thông số đã đƣợc kiểm nghiệm trong quá trình thực hiện. Khoảng 69 phân bố càng nhỏ thì mức độ thể hiện của các mẫu đƣợc chi tiết hơn, dải dự báo càng hẹp nên độ chính xác càng cao. Về việc xác định các chỉ tiêu đánh giá, hiện nay trên thế giới đã có một vài nghiên cứu mới về việc đƣa ra những chỉ tiêu phù hợp hơn trong phƣơng pháp GLUE. Mức độ phù hợp đƣa ra dựa trên quan điểm của ngƣời sử dụng nên còn phụ thuộc nhiều vào yếu tố chủ quan. Trong nội dung luận văn, học viên lựa chọn 3 chỉ tiêu NS, EV, ME đƣợc xem là phổ biến và đƣợc thử nghiệm trong nhiều nghiên cứu. Áp dụng ba chỉ tiêu đã nêu cho thấy chỉ số NS đƣợc sử dụng cho kết quả tốt nhất. Tuy nhiên, việc xây dựng các chỉ tiêu đánh giá độ phù hợp luôn là quan trọng bậc nhất trong quy trình nhƣng chƣa thể hiện đƣợc mức độ phù hợp trong dự báo và cần đƣợc nghiên cứu thêm trong các giai đoạn tiếp theo. 3) Về kết quả tính toán Kết hợp phƣơng pháp GLUE và mô hình WetSpa xây dựng đƣợc sơ đồ tính toán độ bất định nhằm đƣa ra một dải kết dự báo thay vì một giá trị duy nhất. Kết quả mô phỏng tƣơng đối phù hợp với đƣờng quá trình lũ thực đo của hai trận lũ tháng 11/1999 và tháng 10/2003 cho thấy sử dụng chỉ tiêu NS để tính toán bất định cho kết quả phù hợp hơn nên sử dụng chỉ tiêu NS để thực hiện bƣớc dự báo. Khi sử dụng quy trình ở chế độ dự báo (trận lũ tháng 12/1999), kết quả cho thấy có sự chênh lệch giữa giá trị thực đo và các biên trên (đối với đỉnh lũ thứ 2) và biên dƣới (đối với đỉnh lũ thứ 3) của khoảng bất định dự báo. Tuy nhiên mức độ chênh lệch là không lớn. Sự sai khác ở đỉnh lũ thứ 3 có thể đƣợc giải thích một phần là do với cấp lƣu lƣợng lớn (Q>3500 m3/s), tại An Chỉ có hiện tƣợng nƣớc tràn bờ và trao đổi với lƣu vực sông Trà Khúc. Do đó lƣu lƣợng đo đƣợc không phải là tổng lƣu lƣợng thu đƣợc từ thƣợng lƣu chảy về. Nếu chỉ sử dụng đƣờng trung bình là đƣờng dự báo thì sai số sẽ lớn hơn. Từ sự chênh lệch giữa khoảng dự báo và giá trị thực đo có thể cần tiếp tục nghiên cứu các vấn đề sau: i) áp Áp dụng quy trình cập nhật độ phù hợp sử dụng nhiều trận lũ quá khứ; ii) Áp dụng quy trình hiện tại với các lƣu vực khác mà ở đó không có sự trao đổi nƣớc với lƣu vực xung quanh khi có lũ lớn. 70 71 TÀI LIỆU THAM KHẢO Tiếng Việt [1]. Phạm Thị Phƣơng Chi (2009), Sử dụng phương pháp Morris đánh giá độ nhạy các thông số trong mô hình WetSpa, Luận văn thạc sỹ khoa học ngành Thủy văn, Trƣờng Đại học Khoa học Tự nhiên, Hà Nội. [2]. Nguyễn Lan Châu và các cộng sự, Các bài toán trong ứng dụng mô hình thủy văn Marine để mô phỏng và dự báo lũ sông Đà, Tạp chí Khí tƣợng Thủy văn (2005) 1. [3]. Nguyễn Anh Đức (2005), Hiệu chỉnh, áp dụng công thức SCS và mô hình sóng động học phương pháp phần tử hữu hạn mô phỏng quá trình lũ lưu vực sông Vệ - trạm An Chỉ, Khóa luận tốt nghiệp ngành Thủy văn, Trƣờng Đại học Khoa học Tự nhiên, Hà Nội. [4]. Nguyễn Tiền Giang, Daniel van Putten, Phạm Thị Thu Hiền, Công nghệ dự báo lũ khi xét đến tính bất định của mô hình thủy văn: Cơ sở lý thuyết , Tạp chí Khoa học ĐHQGHN, Khoa học Tự nhiên và Công nghệ 25, Số 3S (2009) 403‐411. [5]. Nguyễn Tiền Giang, Nguyễn Thị Thủy, Khai thác mô hình WetSpa phục vụ dự báo lũ các lưu vực sông quốc tế: tính bất định số liệu, tham số, cấu trúc mô hình và đề xuất các giải pháp, Tạp chí Khoa học ĐHQGHN, Khoa học Tự nhiên và Công nghệ 25 (1S), (2009) 35. [6]. Hồ Thị Minh Hà (2008), Nghiên cứu khả năng mô phỏng mùa các yếu tố khí tượng trên lãnh thổ Việt Nam bằng phương pháp thủy động và thống kê , Luận án Tiến sỹ ngành Khí tƣợng, Trƣờng Đại học Khoa học Tự nhiên, Hà Nội. [7]. Nguyễn Thanh Sơn (2008), Nghiên cứu mô phỏng quá trình mưa - dòng chảy phục vụ sử dụng hợp lý tài nguyên nước và đất một số lưu vực sông thượng nguồn miền Trung, Luận án Tiến sỹ ngành Địa lý, Trƣờng Đại học Khoa học Tự nhiên, Hà Nội. 72 [8]. Nguyễn Thị Thủy (2008), Ứng dụng mô hình WetSpa cải tiến dự báo lũ cho lưu vực sông Cả tính đến trạm Dừa, Khóa luận tốt nghiệp ngành Thủy văn, Trƣờng Đại học Khoa học Tự nhiên, Hà Nội. [9]. Viện Quy hoạch Thủy lợi (2003), Quy hoạch sử dụng tổng hợp nguồn nước lưu vực sông Trà Khúc - Tỉnh Quảng Ngãi, Hà Nội. [10]. Bofu Yu (2004), Báo cáo Thủy văn và Hình thái địa hình bồi tích các sông Trà Bồng, Trà Khúc và sông Vệ tại Quảng Ngãi, Việt Nam Tiếng Anh [11]. Aronica G., Bates P. D., Horritt M. S. (2002), Assessing the uncertainty in distributed model predictions using observed binary pattern information within GLUE, Hydrological processes, 16, 2001-2016. [12]. K. J. Beven, Uncertainty in Predictions of Floods and Hydraulic Transport, Publs. Inst. GeoPhys. Pol. Acad. Sc., E-7 (2007) 401. [13]. K.J. Beven, A.M. Binley, The future of distributed models: model calibration and uncertainty prediction, Hydrological Process 6 (1992) 279. [14]. Beven K., Binley A. (1992), The future of distributed models [15]. Beven Keith (2001), How far can we go in distributed hydrological modelling, Hydrology and Earth System Sciences, 5, 1-12. [16]. Campolongo F., Saltelli A. (1997), Sensitivity analysis of an environmental model: an application of different analysis methods, Reliability Engineering & System Safety, 57, 49-69. [17]. Daniel Van Puten (2009), Estimating and updating uncertainty with the GLUE methodology, Bachelor thesis Twente University, Enschede, The Netherlands. [18]. FAO (2006), World reference base for soil resources 2006, Italia. [19]. FAO (2006), FAO Soil Unit, Italia. 73 [20]. Nguyen, T. G., De Kok J. (2006), Systematic testing of an integrated systems model for coastal zone management using sensitivity and uncertainty analyses, Environmental Modelling & Software 22 (2007), 1572-1587. [21]. V.B. Liu, F.De Smedlt, Document and user manual WetSpa extension, Belgium 2004. [22]. Y. B. Liu et al., A diffusive transport approach for flow routing in GIS- based flood modeling, Journal of Hydrology 283 (2003) 91. [23]. Granger Morgan, Max Herion, Mitchell Small (1990), Uncertainty, Cambridge University Press, The United Stated of America. [24]. Iman R.L., Helton J.C. (1988), An investigation of uncertainty and sensitivity analysis techniques for computer models, Risk Analysis 8 (1), 71-90. [25]. NSW Department of Commerce Manly Hydraulics Laboratory (2006), Review and Assesment of Hydrologic/Hydralic Flood Models. [28]. Nurmohamed, R., Naipal, S., De Smedt, F. (2006), Hydrologic modeling of the Upper Suriname River basin using WetSpa and ArcView GIS, Journal of spatial Hydrology, 6, 1-17. [26]. Roberta-Serena Blasone, Jasper A. Vrugt, Henrik Madsen, Dan Rosbjerg, Bruce A. Robinson, George A. Zyvoloski (2008), Generalized likelihood uncertainty estimation (GLUE) using adaptive Markov Chain Monte Carlo sampling, Water Resources, 31, 630-648. [27]. Ryan Fedak (1999), Effect of Spatial Scale on Hydrologic Modeling in a Headwater Catchment, Master Thesis. [28]. Tom Doldersum (2009), Global sensitivity analysis of the WetSpa model, Bachelor thesis, Twente University, Enschede, The Netherlands. 74 [29]. S. Uhlenbrook et al., Prediction uncertainty of conceptual rainfall- runoff models caused by problems in identifying model parameters and structures, Hydrological Sciences Bulletin 44 (5), (1999) 779. [30]. T. Wagener, H.V. Gupta, Model identification for hydrological forecasting under uncertainty, Stochastic Environmental Research and Risk Assessment 19 (2005) 378.