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 m
3
/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.
80 trang |
Chia sẻ: lvcdongnoi | Lượt xem: 2615 | Lượt tải: 0
Bạn đang xem trước 20 trang tài liệu 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ệ, để xem tài liệu hoàn chỉnh bạn click vào nút DOWNLOAD ở trên
độ 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.
Các file đính kèm theo tài liệu này:
- Á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ệ.pdf