Luận án Nghiên cứu xây dựng mô hình giám sát sự bốc, thoát hơi nước của lớp phủ khu vực Tây Bắc Việt Nam từ dữ liệu ảnh vệ tinh

như độ cao của các trạm quan trắc khí tượng trong khu vực. 2.5.1.4. Tia tới sóng ngắn (RS↓) Bức xạ sóng ngắn đi tới là dòng năng lượng bức xạ mặt trời khuếch tán trực tiếp đi tới bề mặt đất (W/m2). Nó được tính toán trong điều kiện bầu trời trong, là hằng số trong thời gian sử dụng ảnh: 𝑅𝑆↓ = 𝐺𝑠𝑐 × 𝑐𝑜𝑠𝜃 × 𝑑𝑟 × 𝜏𝑠𝑤 (2.10) Thông thường 𝑅𝑆↓ nằm trong khoảng 200-1000 W/m 2 phụ thuộc vào vị trí và thời gian chụp ảnh. Trong đó: Gsc - Hằng số mặt trời (1367 W/m2); τsw - Hệ số truyền dẫn qua khí quyển tính theo công thức (2.9); dr - Là nghịch đảo của khoảng cách tương đối giữa mặt trời và trái đất; θ - Góc cao mặt trời thời điểm chụp ảnh (SUN_ELEVATION). Giá trị dr được tính toán theo công thức đề xuất bởi Duffie and Beckman (1980), đồng thời cũng đưa ra trong tài liệu FAO 56 (Crop Evapotranspiration Allen và cộng sự 1998): 𝑑𝑟 = 1 + 0.033 cos (𝐷𝑂𝑌 2𝜋 365 ) (2.11) Trong đó: DOY là các ngày liên tục trong năm; và góc 𝐷𝑂𝑌 2𝜋 365 có đơn vị là radian; giá trị dr nằm trong khoảng 0,97 tới 1,03 và là đơn vị thứ nguyên. 63 2.5.1.5. Tia phát xạ sóng dài (RL↑) Tia phát xạ sóng dài là năng lượng bức xạ nhiệt phát xạ từ bề mặt của trái đất tới khí quyển (W/m2) được tính theo các bước sau: a. Các chỉ số được sử dụng tính toán gồm: Chỉ số thực vật NDVI (Normalized Difference Vegetation Index), Chỉ số thực vật có hiệu chỉnh ảnh hưởng của đất SAVI (Soil Adjusted Vegetation Index), và LAI (Leaf Area Index) chỉ số diện tích lá, được sử dụng tính toán hệ số phát xạ trong (F05). Chỉ số thực vật NDVI là tỷ số giữa sự khác biệt giữa phản xạ của kênh cận hồng ngoại (ρ 5 ) và kênh đỏ (ρ 4 ) đối với ảnh Landsat 8. NDVI = (ρ 5 – ρ 4 ) / (ρ 5 + ρ 4 ) (2.12) Trong đó: ρ 5 là band phổ cận hồng ngoại (Near InfraRed), ρ 4 là band phổ thuộc bước sóng màu đỏ. Thông thường chỉ số NDVI nằm trong khoảng từ -1 đến +1, đối với lớp phủ thực vật bề mặt NDVI nằm trong khoảng 0 đến +1, đối với nước và mây chỉ số NDVI ≤ 0. Chỉ số thực vật có điều chỉnh ảnh hưởng của đất SAVI. Đây là một trong các chỉ số được xây dựng theo hướng biến đổi công thức tính NDVI. Huete 1988, đã đưa chỉ số SAVI với việc thêm thông số L vào công thức ban đầu của NDVI như dưới đây: SAVI = (1 + L) (ρ 5 – ρ 4 ) / (L + ρ 5 + ρ 4 ) (2.13) Trong đó: L là hằng số SAVI; nếu L = 0 thì SAVI trở thành công thức tính NDVI. Thông thường giá trị L = 0,1 – 0,5, giá trị L có thể được tính toán và phân tích từ nhiều ảnh với thảm thực vật không thay đổi, còn độ ẩm đất thì thay đổi. Theo Huete 1988, L = 1 khi mật độ thực vật không đáng kể; L = 0,5 khi thực vật có mật độ trung bình; L=0,25 khi thực vật có mật độ cao. Chỉ số diện tích lá (LAI) là tỷ số giữa tổng diện tích lá trên cây trồng với diện tích mặt bằng để trồng cây đó. Đây là một chỉ số về sinh khối và diện tích lá. Chỉ số tán lá LAI tính toán cho vùng phía nam Idaho theo công thức kinh nghiệm như sau: 64 𝐿𝐴𝐼 = 𝑙𝑛 ( 𝑆𝐴𝑉𝐼𝑚𝑎𝑥 − 𝑆𝐴𝑉𝐼 0.59 ) 0.91 (2.14) Trong đó: Chỉ số SAVI được tính toán từ ảnh vệ tinh cho khu vực thực nghiệm được theo công thức (2.13) lấy giá trị L=0,25. Giá trị lớn nhất cho LAI tương ứng với giá trị lớn nhất của SAVI, giá trị cho SAVI "bão hòa" với sự gia tăng LAI không thay đổi nhiều. Mối quan hệ giữa SAVI và LAI phụ thuộc vào vị trí và loại cây trồng. b. Hệ số phát xạ bề mặt (ε) là tỷ số của năng lượng nhiệt bức xạ bởi bề mặt với năng lượng nhiệt phát ra bởi một vật đen ở nhiệt độ tương tự. Có 2 loại phát xạ bề mặt được sử dụng trong mô hình SEBAL. Đầu tiên là đại diện cho phát xạ nhiệt bề mặt với sự phát xạ nhiệt kênh 10 dải phổ hẹp tương đối của ảnh Landsat 8 (10,3 – 11,3 μm), biểu diễn bởi εNB, Thứ hai là một đại diện cho phát xạ nhiệt bề mặt với sự phát xạ nhiệt trong quang phổ nhiệt dải rộng (6 - 14 μm), biểu diễn bởi εo, giá trị εNB được sử dụng để tính toán nhiệt độ bề mặt (Ts), và giá trị εo sử dụng để tính toán tổng phát xạ năng lượng sóng dài từ bề mặt (F06). Theo (Allen và cộng sự 2002), hệ số phát xạ bề mặt tính toán trong (F06) sử dụng công thức kinh nghiệm, khi NDVI > 0: εNB = 0,97 + 0,0033LAI; với LAI < 3 (2.15) ε0 = 0,95 + 0,01LAI; với LAI < 3 (2.16) và εNB = 0,98, ε0= 0,98 với LAI ≥ 3. Đối với vùng nước, tuyết ta sử dụng bộ lọc trong mô hình cài đặt εNB và ε0. Với vùng nước; NDVI < 0 và α < 0,47, εΝΒ = 0,99 và εo = 0,985 Với vùng tuyết; NDVI < 0 và α ≥ 0,47, εΝΒ = 0,99 và εo = 0,985 c. Nhiệt độ bề mặt đất Ts được tính toán trong F07 dựa vào công thức biến đổi của Plank như sau: 𝑇𝑠 = 𝐾2 𝑙𝑛 ( 𝜀ΝΒ𝐾1 𝑅𝑐 + 1) (2.17) Trong đó: Rc - Bức xạ nhiệt bề mặt đất đã hiệu chỉnh, theo Allen và cộng sự 2002, giá trị Rc được lấy bằng Lλ; εNB - Sự phát xạ nhiệt trong giải hẹp; K1 và K2 hằng 65 số hiệu chỉnh đối với kênh hồng ngoại nhiệt của ảnh Landsat, sử dụng kênh 10 với giá trị K1 = 774,8853 và K2 = 1321,0789 (Cục địa chất Hoa Kỳ (USGS) 2013) d. Phát xạ sóng dài RL↑ được tính trong F08, được tính toán theo công Stefan- Boltzmann như sau: RL↑ = ε0 × σ × Ts 4 (2.18) Trong đó: σ hằng số Stefan - Boltzmann (5,67 × 10-8 W/m2/K4); Ts - Nhiệt độ bề mặt (0K); ε0 hệ số phát xạ bề mặt sử dụng để tính toán tổng phát xạ năng lượng sóng dài từ bề mặt. Giá trị RL↑nằm trong khoảng 200 – 700 W/m 2 phụ thuộc vào vị trí và thời gian chụp ảnh. 2.5.1.6. Chọn điểm nóng và lạnh của ảnh Trong mô hình SEBAL đã sử dụng hai điểm ảnh hiệu chỉnh điều kiện biên cho cân bằng năng lượng. Đây là những điểm ảnh "nóng" và "lạnh" và được đặc trưng bởi các khu vực quan tâm. Điểm ảnh lạnh được chọn là những vùng ẩm có thực phủ và được tưới nước đầy đủ. Nhiệt độ bề mặt và nhiệt độ gần mặt đất được cho là tương tự tại điểm ảnh này. Điểm ảnh nóng được chọn là những vùng khô, không có cây nông nghiệp, ở đó được coi là không có sự bốc thoát hơi nước (Allen và cộng sự 2002). Cả hai điểm ảnh nên đặt trong một khu vực rộng lớn đồng nhất có chứa ít nhất 6 điểm ảnh trong một kênh. Việc lựa chọn các điểm ảnh đòi hỏi kinh nghiệm. Độ chính xác tính toán ET trong SEBAL phụ thuộc vào sự lựa chọn cẩn thận của hai điểm ảnh này. 2.5.1.7. Bức xạ tới sóng dài RL↓ Bức xạ đi tới sóng dài RL↓là dòng bức xạ nhiệt đi tới từ khí quyển (W/m 2), được tính toán theo công thức Stefan - Boltzmann như sau: RL↓ = εa × σ × Ta 4 (2.19) Trong đó: σ hằng số Stefan - Boltzmann (5,67 × 10-8 W/m2/K4); εa – hệ số phát xạ của khí quyển; Ta nhiệt độ gần mặt đất. Theo công thức kinh nghiệm của Bastiaanssen (1995) hệ số phát xạ của khí quyển εa được áp dụng bằng cách sử dụng dữ liệu từ các cánh đồng cỏ linh năng ở Idaho: 66 εa = 0,85 × (-ln τsw)0,09 (2.20) Trong đó; τsw được tính toán từ công thức (2.9). Thay công thức (2.20) vào công thức (2.19) và sử dụng Tcold cho điểm lạnh theo giá trị nhiệt độ gần mặt đất Ta ta được công thức sau: RL↓ = 0,85 × (-ln τsw) 0,09× σ × Tcold4 (2.21) Tính toán này được thực hiện bằng cách sử dụng một bảng tính hoặc máy tính. Giá trị RL↓ nằm trong khoảng 200 - 500 W/m2 phụ thuộc vào vị trí và thời gian chụp ảnh. 2.5.1.8. Tính giá trị năng lượng bức xạ ròng đi tới bề mặt đất Rni Bức xạ ròng đi tới bề mặt đất Rni được tính toán dựa theo các giá trị về suất phân sai bề mặt (α) được tính toán theo F01, F02, F03 và F04; Tia phát xạ sóng dài (RL↑) tính toán theo F05, F06, F07 và F08; Hệ số phát xạ bề mặt (ε0) tính theo F05, F06; Bức xạ sóng ngắn đi tới (RS↓) theo công thức (2.10) và bức xạ sóng dài đi tới (RL↓) theo công thức (2.19). Thay các giá trị vào công thức (2.1) tính được giá trị Rni theo F09. Giá trị Rni nằm trong khoảng 100 – 800 W/m2 phụ thuộc vào bề mặt lớp phủ. 2.5.2. Tính giá trị bức xạ ròng trung bình ngày Rnd từ Rni được tính từ ảnh vệ tinh Landsat 8 Sau khi tính được năng lượng bức xạ ròng tại thời điểm i (Rni) từ ảnh vệ tinh Landsat 8 từ mô hình SEBAL. Yêu cầu đặt ra là phải tính được năng lượng bức xạ ròng trung bình theo ngày từ năng lượng bức xạ ròng tại thời điểm i. Mặt khác, theo kết quả nghiên cứu của (Jackson và cộng sự 1983), dữ liệu năng lượng bức xạ ròng mặt trời tại các thời điểm trong ngày được thể hiện phù hợp nhất với đồ thị hàm sin. Do đó, đã đề xuất công thức tính năng lượng bức xạ mặt trời tại thời điểm i (Rni) như sau: 𝑅𝑛𝑖 = 𝑅𝑛,𝑚𝑎𝑥 𝑠𝑖𝑛(𝜋𝑡/𝑁) (2.22) Trong đó: Rn,max – Bức xạ năng lượng mặt trời lớn nhất (vào khoảng thời gian giữa trưa, thời điểm 12 giờ trưa); t - Khoảng thời gian bắt đầu từ lúc mặt trời mọc tới 67 thời điểm i; N – Khoảng thời gian từ lúc mặt trời mọc đến lúc mặt trời lặn tính theo đơn vị của t. Vào thời điểm giữa trưa, thời điểm 12 giờ trưa thì t = N/2. Khi đó Rni = Rn,max. Để tính năng lượng bức xạ ròng trung bình theo ngày ta sử dụng công thức: 𝑅𝑛𝑑 = ∫ 𝑅𝑛,𝑚𝑎𝑥 𝑠𝑖𝑛(𝜋𝑡/𝑁)𝑑𝑡 = (2𝑁/𝜋). 𝑅𝑛,𝑚𝑎𝑥 𝑁 0 (2.23) Khi đó tỷ số giữa tổng bức xạ hàng ngày và bức xạ tức thời tại thời điểm t là: Rnd/Rni = 2N/[𝜋𝑠𝑖𝑛(𝜋𝑡/𝑁)] (2.24) Hay Rnd = J. Rni (2.25) Trong đó J là hệ số, J = 2N/[𝜋𝑠𝑖𝑛(𝜋𝑡/𝑁)]; Thời gian ban ngày N được tính theo công thức sau: N = a + b.sin2[π.(D + 10)/365] (2.26) Trong đó: D – ngày trong năm (theo ngày tháng chụp ảnh vệ tinh) Theo Neiburger và cộng sự 1973 sử dụng công thức để tính các hệ số a và b dưới dạng một hàm của vĩ độ như sau: a = 12 – 5,69 10-2L – 2,02 10-4L2 + 8,25 10-6L3 – 3,15 10-7L4 (2.27) b = 0,123.L– 3,10 10-4L2 + 8,0 10-7L3 + 4,99 10-7L4 (2.28) L – Vĩ độ của vị trí tính (đơn vị tính độ) Thay các giá trị N, a, b tính theo các công thức (2.26), (2.27), (2.28) vào công thức (2.25) sẽ tính được giá trị năng lượng bức xạ ròng trung bình ngày Rnd 2.5.3. Tính bức xạ ròng trung bình ngày Rnd từ số liệu khí tượng đo trực tiếp tại các trạm quan trắc theo mô hình FAO 56 – Penman - Monteith (Allen và cộng sự 1998), đã đề xuất công thức tính năng lượng bức xạ ròng theo ngày dựa trên số liệu đo đạc trực tiếp tại các trạm quan trắc khí tượng thủy văn như sau: Rnd = Rns – Rnl (2.29) Trong đó: Rns – Bức xạ ròng sóng ngắn mặt trời (MJ/m2/ ngày); Rns – Bức xạ ròng sóng dài mặt trời (MJ/m2/ ngày). Rns = (1 – α).Rs (2.30) 68 Trong đó: α – Suất phân sai bề mặt (albedo); Rs – Bức xạ mặt trời tại các điểm quan sát trên mặt đất (MJ/m2/ngày). 𝑅𝑠 = (𝑎 + 𝑏. 𝑛 𝑁 ) . 𝑅𝑎 (2.31) Trong đó: Ra – Năng lượng bức xạ ròng ở đỉnh khí quyển (MJ/m2/ngày); n - Số giờ nắng thực của ngày (giờ); N – Độ dài của ngày (từ lúc mặt trời mọc đến khi mặt trời lặn); a, b - là các hệ số thể hiện phần bức xạ ngoài từ trái đất (a = 0,25, b = 0,5) 𝑅𝑎 = 24, 60 𝜋 𝐺𝑠𝑐 . 𝑑𝑟 . [𝜔𝑠 . sin(𝜑) . sin(𝛿) + cos(𝜑) . cos(𝛿) . sin(𝜔𝑠)] (2.32) Trong đó: Gsc – Hằng số mặt trời 0,082 (MJ/m2/phút); dr – Tham số hiệu chỉnh khoảng cách giữa mặt trời và mặt đất; ωs – góc giờ mặt trời (rad); δ - góc nghiêng mặt trời (rad); φ – Vĩ độ của điểm quan sát (rad); 𝑅𝑠𝑜 = [(𝑎 + 𝑏) + (2. 10 −5 . 𝑍)] (2.33) Trong đó: a,b - là các hệ số thể hiện phần bức xạ ngoài từ trái đất (a = 0,25, b = 0,5); Z – Là độ cao tại các điểm quan trắc. 𝑅𝑛𝑙 = 𝜎 [ 𝑇𝑚𝑎𝑥,𝐾 4 + 𝑇𝑚𝑖𝑛,𝐾 4 2 ] . (0,34 − 0,14. √𝑒𝑎). (1,35. 𝑅𝑠 𝑅𝑠𝑜 − 0,35) (2.34) Trong đó: σ – Hằng số Stefan – Boltzmann = 4,903.10-9; Tmax, K - Nhiệt độ cao nhất (0K); Tmin, K - Nhiệt độ thấp nhất (0K); ea - áp suất hơi thực tế [kPa]. Các tham số dữ liệu đầu vào để xác định Rnd theo tiêu chuẩn FAO 56 – Penman – Monteith được xác định từ dữ liệu khí tượng đo trực tiếp tại các trạm quan trắc khí tượng thủy văn. 2.5.4. Tính giá trị nhiệt ẩn của quá trình bốc thoát hơi nước (λ) từ dữ liệu ảnh vệ tinh Landsat 8 Giá trị nhiệt ẩn của quá trình bốc thoát hơi nước (λ) là năng lượng cần thiết để thay đổi một đơn vị khối lượng nước từ thể lỏng thành hơi trong một quá trình có áp suất không đổi và nhiệt độ không đổi. Giá trị nhiệt ẩn của quá trình bốc thoát hơi nước (λ) được tính cho từng pixel ảnh theo công thức sau: λ = 2.501 – 0.002361ˣTs (2.36) 69 Trong đó: λ - Giá trị nhiệt ẩn của quá trình bốc thoát hơi nước (MJ/kg); Ts – Nhiệt độ bề mặt (0C). Giá trị nhiệt độ bề mặt được tính trực tiếp từ ảnh vệ tinh Landsat 8 theo các pixel tương ứng. 2.5.5. Tính giá trị hằng số Psychrometric (γ) từ dữ liệu ảnh vệ tinh Landsat 8 và thông tin độ cao từ DEM Hằng số Psychrometric được tính toán thông qua tiêu chuẩn FAO 56 Penman - Monteith như sau: 𝛾 = 𝐶𝑝 𝑃 𝜀 𝜆 (2.37) Trong đó: Cp – Giá trị nhiệt dung riêng ở áp suất không đổi (Cp = 1,013 10-3 (MJ/kg/0C)); ε – Tỷ lệ khối lượng phân tử của hơi nước/không khí khô (ε = 0,622); λ - Giá trị nhiệt ẩn của quá trình bốc thoát hơi nước (MJ/kg) được tính theo công thức (2.38); P – Áp suất khí quyển (kPa); Giá trị áp suất khí quyển P phụ thuộc vào độ cao và được tính trực tiếp dựa vào mô hình số độ cao (DEM) như sau: 𝑃 = 101,3 ( 293 − 0.0065 𝑍 293 ) 5.26 (2.38) Trong đó: Z – Độ cao các điểm tính từ mực nước biển (m), giá trị độ cao Z được tính theo từng pixel ảnh từ DEM 2.5.6. Tính giá trị độ dốc của đường cong áp suất hơi bão hòa (Δ) từ dữ liệu ảnh vệ tinh Landsat 8 Giá trị độ dốc đường cong áp suất hơi nước bão hòa ∆ phụ thuộc vào tham số nhiệt độ không khí, tiêu chuẩn FAO 56 Penman - Monteith đã xác định độ dốc của đường cong áp suất hơi nước bão hòa bằng công thức sử dụng nhiệt độ không khí như sau: Δ = 4098 [0,6108. 𝐸𝑋𝑃 ( 17,27 𝑇𝑆 𝑇𝑎 + 237,3 )] (𝑇𝑆 + 237,3)2 (2.39) Trong đó: TS – nhiệt độ không khí gần bề mặt (0C). Giá trị Ts được tính trực tiếp từ ảnh vệ tinh Landsat 8 theo các pixel ảnh. 70 2.6. Xác định hệ số a, b của mô hình Priestley - Taylor phù hợp với địa hình khí hậu khu vực Tây Bắc Việt Nam Theo công thức tính lượng bốc thoát hơi nước (1.3) mô hình Priestley - Taylor các hệ số a, b của mô hình được đề xuất theo kết quả tính thực nghiệm cho từng khu vực khác nhau. Để xác định hệ số a, b của mô hình Priestley - Taylor phù hợp với điều kiện địa hình chia cắt mạnh, chênh cao lớn và bề mặt lớp phủ với nhiều trạng thái cây trồng khác của khu vực Tây Bắc Việt Nam. Nghiên cứu sinh đã sử dụng chuỗi dữ liệu khí tượng, thủy văn đo trực tiếp tại 11 trạm khí tượng thủy văn (08 trạm khí tượng thủy văn tại tỉnh Hòa Bình và 03 trạm khí tượng tại tỉnh Sơn La), tại mỗi trạm các dữ liệu khí tượng được thu thập phục vụ tính toán gồm Nhiệt độ, Độ ẩm, Số giờ nắng thực, Vận tốc gió, lượng bốc thoát hơi nước thực tế tại thời điểm các ngày 01/7/2015, ngày 15/6/2016, ngày 04/6/2017, ngày 25/7/2018, ngày 09/7/2019, ngày 07/7/2020 và ngày 18/8/2021 do Đài khí tượng thủy văn tỉnh Hòa Bình cung cấp. Tại mỗi thời điểm sử dụng dữ liệu khí tượng tại 8/11 trạm quan trắc khí tượng thủy văn gồm (Khí tượng Hòa Bình, Khí tượng Mai Châu, Khí tượng Kim Bôi, Khí tượng Chi Nê, Khí tượng Lạc Sơn, Thủy văn Hòa bình, Khí tượng Phù Yên, Khí tượng Mộc Châu) để tính toán xác định hệ số a, b của mô hình Priestley – Taylor và 3/11 trạm quan trắc khí tượng thủy văn gồm (Thủy văn Hưng Thi, Thủy văn Lâm Sơn, Khí tượng Bắc Yên) để kiểm tra kết quả tính. Cụ thể tại 8 trạm khí tượng thủy văn sử dụng để xác định hệ số a, b của mô hình Priestley-Taylor sẽ áp dụng tiêu chuẩn FAO 56 để xác định giá trị năng lượng bức xạ ròng trung bình ngày Rnd theo công thức (2.29), giá trị nhiệt ẩn quá trình bốc thoát hơi nước λ theo công thức (2.36), giá trị hằng số psychrometric γ theo công thức (2.37) và giá trị độ dốc đường cong áp suất hơi nước bão hòa Δ theo công thức (2.39). Sau khí tính được các giá trị (Rnd, λ, γ, Δ) từ dữ liệu khí tượng tại 8 trạm khí tượng thủy văn sử dụng để xác định hệ số a, b của mô hình Priestley-Taylor, kết hợp với lượng bốc thoát hơi nước thực tế đo trực tiếp tại 8 trạm khí tượng thủy văn tương ứng tại thời điểm các ngày 01/7/2015, ngày 15/6/2016, ngày 04/6/2017, ngày 25/7/2018, ngày 09/7/2019, ngày 07/7/2020 và ngày 18/8/2021. Thay các giá trị (Rnd, λ, γ, Δ) được tính từ dữ liệu khí tượng và giá trị bốc thoát hơi nước thực tế đo trực tiếp tại các trạm khí tượng thủy văn vào phương trình (1.3). Khi đó tương ứng với 71 các thời điểm ngày 01/7/2015, ngày 15/6/2016, ngày 04/6/2017, ngày 25/7/2018, ngày 09/7/2019, ngày 07/7/2020 và ngày 18/8/2021 tại mỗi vị trí đặt trạm quan trắc khí tượng thủy văn sẽ lập được 01 phương trình dạng (1.3) mỗi phương trình khi đó chỉ có 2 ẩn số là a và b. Như vậy, tương ứng mỗi thời điểm ngày 01/7/2015, ngày 15/6/2016, ngày 04/6/2017, ngày 25/7/2018, ngày 09/7/2019, ngày 07/7/2020 và ngày 18/8/2021 với 8 trạm quan trắc khí tượng thủy văn được sử dụng để xác định hệ số a, b lập được 8 phương trình có dạng: Aiaj + Bibj + Li = 0 (2.40) Chỉ số j tương ứng với mỗi thời điểm, chỉ số i tương ứng với vị trí các trạm khí tượng thuỷ văn, hệ số 𝐴𝑖 = ∆𝑖 ∆𝑖−𝛾𝑖 × 𝑅𝑛𝑑𝑖 λ𝑖 , hệ số Bi = 1 và Li hệ số tự do. Giải hệ phương trình theo nguyên lý số bình phương nhỏ nhất sẽ xác định được hệ số aj, bj tương ứng với từng thời điểm ngày 01/7/2015, ngày 15/6/2016, ngày 04/6/2017, ngày 25/7/2018, ngày 09/7/2019, ngày 07/7/2020 và ngày 18/8/20211. Giá trị a, b của mô hình Priestley – Taylor đề xuất được cho là phù hợp với điều kiện địa hình, khí hậu khu vực Tây Bắc Việt Nam là giá trị trung bình của 7 thời điểm tính tương ứng với ngày 01/7/2015, ngày 15/6/2016, ngày 04/6/2017, ngày 25/7/2018, ngày 09/7/2019, ngày 07/7/2020 và ngày 18/8/2021. 2.7. Đề xuất mô hình, quy trình ước tính, giám sát lượng bốc thoát hơi nước sử dụng kết hợp mô hình SEBAL với dữ liệu ảnh vệ tinh Landsat 8 và mô hình Priestley - Taylor 2.7.1. Đề xuất mô hình ước tính giám sát lượng bốc thoát hơi nước bề mặt lớp phủ khu vực Tây Bắc Việt Nam Mỗi mô hình ước tính giám sát lượng bốc thoát hơi nước bề mặt lớp phủ khác nhau có những ưu nhược điểm khác nhau tùy thuộc vào điều kiện địa hình, khí hậu, chênh cao địa hình và thực tế bề mặt lớp phủ của các khu vực đó, phụ thuộc vào trang thiết bị, trình độ khoa học công nghệ và yêu cầu các tham số đầu vào của các mô hình đó. Khu vực Tây Bắc Việt Nam với điều kiện địa hình khí hậu đặc trưng như địa hình chia cắt mạnh, chênh cao lớn, nhiều tiểu vùng khí hậu, hệ thống sông suối, diện tích mặt nước lớn và đặc biệt bề mặt lớp phủ với nhiều trạng thái cây trồng khác nhau, không thuần loài. Ngoài ra dữ liệu ảnh vệ tinh Landsat 8 với độ phân giải không gian, 72 thời gian phù hợp, có thể tính toán các tham số phục vụ ước tính, giám sát lượng bốc thoát hơi nước thay cho việc sử dụng dữ liệu khí tượng đo trực tiếp tại các trạm quan trắc khí tượng đảm bảo độ chính xác, tin cậy, hiệu quả và phù hợp với điều kiện về trang thiết bị máy móc, trình độ khoa học công nghệ tại Việt Nam. Có thể tự động hóa trong công tác ước tính giám sát lượng bốc thoát hơi nước bề mặt lớp phủ nhằm xác định nhu cầu nước của cây trồng và cảnh báo hạn hán, thiên tai, cháy rừng. Do vậy, đề tài luận án đã đề xuất sử dụng kết hợp mô hình SEBAL với dữ liệu ảnh vệ tinh Landsat 8 và mô hình Priestley-Taylor với các hệ số a, b được tính toán từ thực nghiệm để ước tính giám sát lượng bốc thoát hơi nước bề mặt lớp phủ khu vực Tây Bắc Việt Nam. 2.7.2. Quy trình ước tính, giám sát lượng bốc thoát hơi nước sử dụng kết hợp mô hình SEBAL với dữ liệu ảnh vệ tinh Landsat 8 và mô hình Priestley - Taylor Quy trình ước tính, giám sát lượng bốc thoát hơi nước sử dụng kết hợp mô hình SEBAL với dữ liệu ảnh vệ tinh Landsat 8 và mô hình Priestley - Taylor với các hệ số a, b xác định từ thực nghiệm được tiến hành qua các bước chính như sau: 73 Hình 2.8. Quy trình ước tính, giám sát lượng bốc thoát hơi nước thực tế ETa sử dụng kết hợp mô hình viễn thám SEBAL với dữ liệu ảnh vệ tinh Landsat 8 và mô hình Priestley - Taylor Bức xạ phổ của các kênh ảnh Lλ (Spectran radiance) B1 Hệ số phản xạ của các kênh ảnh 𝜌λ (Reflectivity) B2 Chỉ số thực vật (NDVI) Chỉ số (SAVI) Chỉ số diện tích lá (LAI) B3 Suất phân sai tại đỉnh khí quyển αtoa (Albedo-top of atmosphere) B4 Hệ số phát xạ bề mặt εNB và ε0 (Surface emissivities) B5 Nhiệt độ bề mặt TS (Surface temperature) B6 Ảnh vệ tinh Landsat 8 (DN) Suất phân sai bề mặt α (Surface albedo) B8 Tia tới sóng ngắn 𝑅𝑆↓ (Incoming shortwave) B7 Tia phát xạ sóng dài 𝑅𝐿↑ (Out going longwave) B9 Tia tới sóng dài 𝑅𝐿↓ (Incoming longwave) B10 Tính năng lượng bức xạ ròng thời điểm i (Rni) SEBAL B11 Tính năng lượng bức xạ ròng trung bình ngày (Rnd) B12 Tính lượng bốc thoát hơi nước thực tế ETa theo mô hình Priestley – Taylor (a,b xác định bằng thực nghiệm) B14 Biên tập bản đồ bốc thoát hơi nước thực tế ETa B15 Thông tin độ cao từ DEM Nhiệt ẩn bốc thoát hơi nước (λ); Hằng số Psychrometric (γ); Độ dốc đường cong áp suất hơi nước bão hòa (Δ) B13 74 Theo (Cục địa chất Hoa Kỳ (USGS) 2013) ảnh vệ tinh Landsat 8 được download miễn phí từ websize https://earthexplorer.usgs.gov/, hệ tọa độ WGS-84 múi chiếu 60, sau khi tăng cường chất lượng hình ảnh được cắt theo ranh giới hành chính tỉnh Hòa Bình. Quy trình ước tính, giám sát lượng bốc thoát hơi nước lớp phủ từ ảnh vệ tinh được tiến hành qua các bước như sau: B1: Tính chuyển giá trị pixel từ dạng số DN sang bức xạ phổ của các kênh ảnh 10 và 11 phục vụ tính nhiệt độ bề mặt; B2: Chuyển đổi giá trị điểm ảnh từ dạng số (DN) sang giá trị phản xạ bề mặt (Reflectance) phục vụ tính các chỉ số NDVI, SAVI, LAI và suất phân sai tại đỉnh khí quyển αtoa B3: Tính các chỉ số NDVI, SAVI, LAI phục vụ tính hệ số phát xạ bề mặt dải rộng (ε0) và dải hẹp (εNB). B4: Tính suất phân sai tại đỉnh khí quyển αtoa phục vụ tính suất phân sai bề mặt đất (α) và tính hệ số phát xạ bề mặt dải rộng (ε0) và dải hẹp (εNB). B5: Tính hệ số phát xạ bề mặt dải rộng (ε0) và dải hẹp (εNB). Giá trị εNB được sử dụng để tính toán nhiệt độ bề mặt (Ts), và ε0 sử dụng sau để tính toán tổng phát xạ năng lượng sóng dài từ bề mặt. B6: Tính nhiệt độ bề mặt (Ts) phục vụ tính các năng lượng phát xạ sóng dài 𝑅𝐿↑; năng lượng tới sóng dài 𝑅𝐿↓; Nhiệt ẩn bốc thoát hơi nước (λ) và Độ dốc đường cong áp suất hơi nước bão hòa (Δ) B7, B8, B9, B10: Tính các giá trị năng lượng tới sóng ngắn 𝑅𝑆↓; Suất phân sai bề mặt α; Tia phát xạ sóng dài 𝑅𝐿↑; Tia tới sóng dài 𝑅𝐿↓ phục vụ tính tính năng lượng bức xạ ròng thời điểm i (Rni) theo mô hình SEBAL. B11: Tính năng lượng bức xạ ròng trung bình ngày Rnd từ năng lượng bức xạ ròng thu được tại thời điểm i từ ảnh Landsat 8 để phục vụ tính lượng bốc thoát hơi nước theo mô hình Priestley – Taylor. B12: Tính năng lượng bức xạ ròng trung bình ngày từ dữ liệu ảnh vệ tinh Landsat 8 75 B13: Tính giá trị Nhiệt ẩn bốc thoát hơi nước (λ); Hằng số Psychrometric (γ); Độ dốc đường cong áp suất hơi nước bão hòa (Δ) từ nhiệt độ bề mặt và thông tin độ cao từ DEM phục vụ tính lượng bốc thoát hơi nước theo mô hình Priestley – Taylor. B14: Tính lượng bốc thoát hơi nước theo mô hình Priestley – Taylor từ dữ liệu tính được trong B12 và B13 (a, b xác định bằng thực nghiệm). B15: Biên tập bản đồ bốc thoát hơi nước thực tế ETa theo hệ tọa độ WGS-84 múi chiếu 60 tại các thời điểm thu nhận ảnh vệ tinh. Để có thể tính toán lượng bốc thoát hơi nước theo qui trình trên được nhanh chóng, hiệu quả tại bất kỳ thời điểm nào có ảnh vệ tinh Landsat 8 cần thiết phải xây dựng phần mềm tính lượng bốc thoát hơi nước. 2.8. Xây dựng chương trình ước tính, giám sát lượng bốc thoát hơi nước sử dụng kết hợp mô hình SEBAL và mô hình Priestley-Taylor trên nền Google Earth Engine 2.8.1. Khái quát về Google Earth Engine Google Earth Engine (GEE) là một nền tảng dựa trên kỹ thuật điện toán đám mây để phân tích không gian địa lý trên quy mô toàn cầu. GEE mang đến khả năng tính toán khổng lồ của Google để giải quyết nhiều vấn đề xã hội bao gồm phá rừng, hạn hán, thảm họa, bệnh tật, an ninh lương thực, quản lý nước, giám sát khí hậu và bảo vệ môi trường. GEE là nền tảng duy nhất trong lĩnh vực này, nó là nền tảng tích hợp (sử dụng siêu máy tính và tài nguyên điện toán đám mây quy mô lớn) được thiết kế cung cấp năng lực tính toán không chỉ dành riêng cho các nhà khoa học viễn thám, mà còn rộng lớn hơn trong các lĩnh vực khác nhau. Earth Engine bao gồm một danh mục dữ liệu sẵn sàng phân tích (petabyte) cùng với một dịch vụ tính toán song song hiệu năng cao, được truy cập và kiểm soát thông qua giao diện lập trình ứng dụng (API) có thể truy cập Internet và môi trường phát triển tương tác dựa trên web cho phép tạo mẫu nhanh và trực quan hóa kết quả. 2.8.2. Những lợi ích của Google Earth Engine trong việc xây dựng chương trình ước tính, giám sát lượng bốc thoát hơi nước Google Earth Engine cho phép người dùng tính toán trên hàng Petabytes (một triệu tỉ byte) dữ liệu một cách nhanh chóng mà không cần phải điều hướng sự phức 76 tạp của tính năng chạy song song cloud. Tăng cường khả năng tiếp cận toàn diện đ