单源能量平衡模型反演灌溉牧草场实际蒸散量及其时空特征分析

Single-Source Energy Balance Model Inversion of Actual Evapotranspiration from Irrigated Pasture and Analysis of Its Spatio-Temporal Characteristics

摘要: 水资源短缺是制约干旱地区经济发展的重要因素，所以探究干旱地区生长季人工灌溉牧草实际蒸散量时空变化规律非常必要。本研究使用2013~2021年Landsat-8OIL影像基于pySEBAL模型结合气象站数据反演锡林浩特市沃原奶牛场种植基地人工灌溉牧草实际蒸散量。结果表明：由于人工灌溉牧草水分充足，该区域实际蒸散量较高，介于0~10 mm/d、50~200 mm/m之间。研究区牧草生长季(5~9月)实际蒸散量变化规律整体呈“单峰型”，其中峰值日为2014年7月25日(5.46 mm/d)，峰值月为7月(155.71 mm/m)，种植区界限随牧草蒸散量升高而清晰。此外，实际蒸散量变化规律与降水频率和降水量也有很好的相关性。模型反演所得实际蒸散量与FAO Penman-Monteith算法结果比较(R2 = 0.7504、RMSE = 1.2575 mm/d、MRE = 0.9366 mm/d)，模型整体精度较高，能够较为准确地反映该区牧草实际蒸散情况。反演结果基本符合实际情况，可以为人工灌溉牧草用水管理提供方案，为政府合理政策制定提供依据。

Abstract: Water shortage is an important factor restricting economic development in arid areas, so it is very necessary to explore the spatiotemporal changes in actual evapotranspiration of artificially irrigated forage during the growing season in arid areas. This study uses Landsat-8 OIL images from 2013 to 2021 based on the pySEBAL model and weather station data to invert the actual evapotranspiration of artificially irrigated pasture at the Woyuan Dairy Farm planting base in Xilinhot City. The results show that due to the sufficient water content of artificially irrigated pasture, the actual evapotranspiration in this area is relatively high, ranging from 0 to 10 mm/d and 50 to 200 mm/m. The change pattern of actual evapotranspiration during the pasture growing season (May-September) in the study area showed a “single peak” overall, with the peak day being July 25, 2014 (5.46 mm/d) and the peak month being July (155.71 mm/d). The boundaries of planting areas become clearer as the evapotranspiration of pasture increases. In addition, the change pattern of actual evapotranspiration is also well correlated with precipitation frequency and precipitation amount. Comparing the actual evapotranspiration obtained by the model inversion with the results of the FAO Penman-Monteith algorithm, the overall accuracy of the model is higher and can more accurately reflect the actual evapotranspiration of forage in the area (R2 = 0.7504, RMSE = 1.2575 mm/d, MRE = 0.9366 mm/d). The inversion results are basically consistent with the actual situation, and can provide solutions for artificial irrigation pasture water management and provide a basis for the government to formulate reasonable policies.

1. 引言

气候变化和人口持续增加给地球上有限的淡水资源带来巨大压力 [1] [2] 。水资源短缺已成为世界上大多数干旱地区发展的限制因素 [3] 。为了应对干旱缺水，在牧草生长过程中进行人工灌溉的情况逐渐增多。人工灌溉可在降雨较少的情况为牧草补充水分、改善草地生态环境和促进牧草生长。国内外科研和生产实践表明，人工灌溉草地产草量相比天然草地提高20~40倍 [4] 。锡林郭勒盟位于内蒙古自治区中部，是我国重要的畜产品基地和青储基地 [5] 。该区域内可利用水资源有限，近年的干旱导致牧草减产、草地退化、沙化和区域生态环境不断恶化；另一方面由于人类活动导致地下水枯竭和用水需求增加，这就要求当地需要通过更先进的牧草灌溉管理方案来提高水的利用效率，优化农业用水管理。

蒸散发(Evapotranspiration)是连接水、能源和碳循环的主要水文通量，在水文、气象和农业水管理中发挥着重要作用 [6] ，而实际蒸散量(ET_a)通常被认为是作物需水量的代表 [7] 。传统的监测技术如波文比法、涡流协方差法、称重蒸渗仪法和土壤水分平衡法等，在小范围内可以相对准确测量ET_a，但不便用于大面积时空特征分析 [8] [9] 。基于遥感的蒸散模型具有反演高时空分辨率ET_a能力，能准确有效地反映ET_a时空连续性和动态性。通常使用基于遥感的单源能量平衡模型反演ET_a，包括陆地表面能量平衡模型(SEBAL) [10] 、高分辨率内化校正蒸散发模型(METRIC) [11] 、和地表能量平衡模型(SEBS)等 [12] 。长期以来，SEBAL模型因为不需要事先了解作物类型、土壤和管理措施等现场条件，一直被认为是最适用于估算ET_a的遥感模型 [2] [3] 。SEBAL模型采用迭代算法，在研究区内选定冷热极端像素计算显热通量和潜热通量，并使用遥感影像信息计算每个像素ET_a。pySEBAL由IHE-Delft研究所开发(https://ocw.un-ihe.org/)，在SEBAL模型基础上嵌入python自动化像素选择程序，可以有效减少人工选择极端像素时间、消除主观选择极端像素所引起的误差。长期以来，SEBAL模型因为不需要事先了解作物类型、土壤和管理措施等现场条件，一直被认为是最适用于估算ET_a的遥感模型。Xue等人利用pySEBAL模型成功在加利福尼亚州农场反演番茄、马铃薯、杏树苗的实际蒸散量 [13] ；苏婷婷等人基于SEBAL模型以较高精度估算了土默特平原农田蒸散量 [14] 。但就目前来说，SEBAL模型对于估算人工灌溉牧草ET_a的适用性还需进一步验证。

基于遥感的蒸散发模型研究已经十分成熟，在世界各地反演各种作物ET_a均有成功案例，但对人工灌溉牧草ET_a反演的研究较少。因此，本文使用遥感影像基于pySEBAL模型结合气象站数据反演人工灌溉牧草的ET_a；检验此模型反演人工灌溉牧草ET_a的适用性；分析牧草生长季ET_a时空变化规律。通过本研究对人工灌溉草场的用水管理提出建议，为政府以及相关部门政策制定提供理论依据。

2. 材料与方法

2.1. 研究区概况

研究区位于内蒙古自治区锡林郭勒盟锡林浩特市北郊沃原奶牛场种植基地，是全国首家现代化畜牧业科技示范园区、国家综合开发草原建设及万亩节水灌溉高产饲草料基地，每年为自治区提供亿斤饲草料。研究区概况图如图1所示，该区域内主要种植作物为牧草，种植面积约15.5 km²，研究区概况如图1所示 [15] 。属于温带干旱、半干旱大陆性气候区，年平均气温0℃~3℃，无霜期90~120 d [16] 。年平均降水量295 mm，降雨多集中在7、8、9三个月内 [17] 。年平均相对湿度在60%以下，年蒸发量在1500~2700 mm之间。

Figure 1. Overview of the study area

图1. 研究区概况图

2.2. pySEBAL模型原理

SEBAL (Surface Energy Balance Algorithm for Land)模型是以遥感为基础建立起来的蒸散量反演模型，1998年由Bastiaanssen最先提出 [18] 。pySEBAL模型是Hessels等人在Python环境中开发的SEBAL算法版本，基于地表能量平衡原理通过估算净辐射通量、土壤热通量和感热通量来计算潜热通量，推导出瞬时ET，计算公式如下 [19] ：

$\lambda ET=R_n-G-H$ (1)

式中，λET为潜热通量(W/m²)；λ为水的汽化潜热；R_n为地表净辐射通量(W/m²)；G为土壤热通量(W/m²)；H为显热通量(W/m²)。净辐射通量R_n (W/m²)为地表能量、动量和水运输和交换的主要能量来源，计算公式如下：

$R_n=\left(1-\alpha \right)R_{S\downarrow }+{\epsilon }_0{R}_{L\downarrow }-{\epsilon }_0\sigma T_s^4$ (2)

式中， $R_{S\downarrow }$ 到达地表太阳短波辐射(W/m²)，即入射短波辐射； $R_{L\downarrow }$ 到达地表太阳长波辐射，即入射长波辐射(W/m²)。α为表面反照率(−)。ε₀通过NDVI和叶面积指数(LAI)的半经验关系得到的地表比辐射率，可从红和近红外波段检索。σ为史蒂芬玻尔兹曼常数，为5.67*10⁻⁸ (W/m²∙K⁴)，T_S为地表温度(土壤和植被综合辐射温度) (K)。pySEBAL利用Bastiaanssen开发的经验方程，将土壤热通量G计算为R_n的一部分 [18] ：

$G=T_{s,datum}\left(0.0038+0.007\alpha \right)\left(1-0.98{\text{NDVI}}^4\right)\times R_n$ (3)

式中，T_s,datum为基于研究区数字高程地图(DEM)坡度和角度修正后的地表温度(T_S)。显热通量也叫做感热通量，是指由于温度变化而引起大气与下垫面之间发生湍流形式的热交换。其计算过程比较复杂，计算公式如下：

$H=\frac{\rho \times C_p\times dT}{r_{ah}}$ (4)

式中：ρ为空气密度(kg/m³)；C_p为空气定压比热容(J/kg∙K)；dT为距离地面高度Z₁和Z₂处温度差，通常取0.1 m和2 m；r_ah为热量传输空气动力学阻抗(s/m)。公式5中系数“a”和“b”为针对冷像素和热像素迭代确定，因此它们对于每张影像都是特定的 [20] [21] 。

$dT=a+b\times T_{s,datum}$ (5)

$a=\frac{d{T}_{hot}-d{T}_{cold}}{T_{s,datum,hot}-d{T}_{s,datum,cold}}$ (6)

$b=\frac{d{T}_{hot}-a}{T_{s,datum,hot}}$ (7)

在pySEBAL中，通过识别三个像素来确定极端像素，即(1)植物像素、(2)水体像素和(3)热点像素。植物像素被自动识别为场景中具有最大NDVI的像素。如Jaafar和Ahmad所述，使用大气顶部反射率带对水体像素进行分类，水体像素为非冻结温度和负NDVI组合，冷像元T_s在植物像素和水体像素最小值中选择，热像元被识别为NDVI值在0.03和0.2之间的像元 [22] 。以上各通量是基于影像计算得到，均为瞬时通量。因此，需要通过转换得出日通量。基于卫星过境时间的瞬时R_n、H和G，可以计算瞬时蒸发分数(EF_i) (公式8)，并通过使用平流因子Ω将其转换为每日蒸发分数(EF₂₄)(公式9)，该因子用于减少下午ET_a增加引起的误差 [23] ：

$E{F}_i=\frac{R_n-H-G}{R_n-G}$ (8)

$E{F}_{24}=\Omega \times E{F}_i$ (9)

Ω通过以下公式计算：

$\Omega =1+0.985\times E{F}_i\times \left\{\exp \left[0.08\times \left(e_s-e_a\right)\right]-1\right\}$ (10)

式中，e_s为冠层参考高度上方气温饱和蒸汽压，e_a为冠盖高度上方实际蒸汽压。计算每日ET_a公式如下：

$E{T}_{24}=\frac{8.64\times {10}^7\times \Omega \times E{F}_i\times \left(R_{n24}-G_{24}\right)}{\lambda \times {\rho }_w}$ (11)

式中，ET₂₄卫星过境当日实际蒸散量(mm/d)，λ为潜热通量(J/kg)，ρ_w为水密度(kg/m³)。G₂₄为日平均土壤热通量(W/m²)，假设土壤和植被表面为0。R_n24为当天的平均净辐射(W/m²)，可计算为：

$R_{n24}=\left[\left(1-\alpha \right)\times R_a-110\right]\times {\tau }_{sw}$ (12)

式中，R_a为每日地外太阳辐射(W/m²)。τ_sw为大气透射率，受空气中湿度、灰尘和其他污染物影响。关于pySEBAL算法更多细节，详见Hessels等人pySEBAL3.3.7手册 [24] 。1998年，联合国粮食及农业组织推荐使用FAO Penman-Monteith方法计算作物蒸散发，其基本原理为：

$E{T}_0=\frac{0.408\Delta \left(R_n-G\right)+\gamma \frac{900}{T+273}{u}_2\left(e_s-e_a\right)}{\Delta +\gamma \left(1+0.34{\mu }_2\right)}$ (13)

$E{T}_a=K_c\times E{T}_0$ (14)

式中：ET₀为潜在蒸散发(mm)；K_c为作物系数；R_n为作物表层净辐射(W/m²)；G为土壤热通量(W/m²)；γ为干湿表常数(kPa/℃)；T为日平均气温(℃)；u₂为2米高度处风速(m/s)；Δ为饱和水汽压曲线斜率(kPa/℃)；e_s为饱和水汽压(kPa)；e_a为实际水汽压(kPa)。

2.3. 数据与处理

(1) 遥感数据

Landsat-8 OLI遥感影像来源于美国地质调查局官网(http://glovis.usgs.gov/)，空间分辨率30 m，重返周期为16天。如图2所示，选取2013~2021年生长季(5-9月)无云或少云影像共25幅，用于估算研究区牧草ET_a。

太阳辐射数据来自美国家航空航天局提供的大气再分析数据MERRA-2 (The second Modern-Era Retrospective analysis for Research and Applications version 2)，其空间分辨率为0.625˚ × 0.5˚，时间分辨率为1小时。

DEM数据来源于地理空间数据云(http:/www.gscloud.cn/)，空间分辨率为30 m。

将Landsat-8 OLI影像裁剪到研究区大小，以缩短pySEBAL模型的计算时间。所有的图像预处理都是使用QGIS和ArcGIS，pySEBAL模型自动执行大气校正处理。

(2) 气象数据

气象数据来自美国国家海洋和大气管理局(NOAA)国家环境信息中心网站(NCEI) (https://www.ncei.noaa.gov/)国际交换数据集。选取数据集中内蒙古锡林浩特市气象站2013~2021年生长季(5~9月)的温度、相对湿度、日照时数、风速、降水量等气象数据。

本研究使用的K_c来自于侯琼等人计算所得的不受水分胁迫内蒙古典型草原作物系数，如表1所示 [25] 。

Figure 2. 25 Landsat-8 images from 2013 to 2021 used for pySEBAL irrigated pasture ET_a estimation

图2. 2013~2021年25幅Landsat-8影像用于pySEBAL受灌溉牧草ET_a估算

Table 1. Crop coefficient table

表1. 作物系数表

3. 结果与分析

3.1. pySEBAL模型ET_a反演分析

由于可用影像数量较少，仅2014年生长季各月有连续可用影像。故以2014年为例，分析研究区生长季日、月ET_a时空变化特征。

图3是2014年研究区牧草生长季日实际蒸散量时空分布和频率分布图，ET_a范围为0~10 mm/d，图3(c)、图3(d)两图中蒸散发较高的圆形区域为牧草种植区。图3(a)中牧草处于生长期早期，研究区蒸散量均值为3.51 mm/d，大部分区域在2~5 mm/d之间。图3(b)中6月蒸散量明显增加，均值为4.1 mm/d。图3(c)中种植区与非种植区蒸散量差距较大，种植边界清晰，种植区ET_a在6~7 mm/d之间，靠近牧场边界的非种植区域小于4 mm/d，均值6.17 mm/d为生长季内最高水平。图3(d)中牧草ET_a为4~6 mm/d，均值为4.38 mm/d，进入生长季后期，ET_a逐渐降低，非种植区降低显著。9月牧草刈割完成，图3(e)中ET_a均值为2.19 mm/d，空间分布均匀，为2014年生长季各月最低水平。

图4为2013~2021年研究区生长季(5~9月)日实际蒸散量时空分布图和频率分布直方图。5月研究区蒸散发水平较低，空间分布均匀，多集中在2~4 mm/d，均值2.98 mm/d。6月牧草快速生长，ET_a均值升高到4.31 mm/d。7月种植区与非种植区蒸散发差异显著，牧草ET_a在5~7 mm/d之间，均值达到5.46 mm/d，而靠近研究区边界的非种植区ET_a小于3 mm/d。8月牧场ET_a降低至3~5 mm/d，均值为4.12 mm/d。9月由于牧草刈割完成，种植区与非种植区ET_a均处于较低水平1~3 mm/d，且空间分布均匀。图3、图4中ET_a空间分布、像元频率、发展趋势基本一致，遵循“5~6月增长，7月达到峰值，8~9月下降”的规律，但图4表现出长时间序列上生长季各月ET_a空间分布更均匀、时间变化更平稳。

Figure 3. Spatiotemporal distribution and frequency distribution histogram of actual daily evapotranspiration in the 2014 growing season

图3. 2014年生长季日实际蒸散量时空分布和频率分布直方图

Figure 4. Spatiotemporal distribution and frequency distribution histogram of actual daily evapotranspiration during the growing season from 2013 to 2021

图4. 2013~2021年生长季日实际蒸散量时空分布和频率分布直方图

本研究对2014年日ET_a在时间尺度上进行聚合处理，得到2014年生长季月ET_a时空分布以及频率分布如图5所示，ET_a在50~200 mm/m范围内。2014年5月研究区绝大部分区域ET_a大于100 mm/m，均值为128.76 mm/m。6月研究区蒸散发量增加，大部分区域ET_a高于120 mm/m，均值升高至138.67 mm/m。进入7月草场ET_a达到峰值，在牧草ET_a在160~180 mm/m之间，均值为155.71 mm/m。8月牧场ET_a开始降低，非种植区域下降明显，多数区域ET_a在140~160 mm/m之间，均值降低至133.6 mm/m，略低于6月。9月为牧草刈割时段，牧草蒸散量降至120~140 mm/m之间，均值117.25 mm/m为生长季月ET_a最低值。

月实际蒸散量空间分布和时间动态与日值基本一致，5月份为生长季早期，种植区与非种植区植被蒸散发水平趋同；6月牧草进入快速生长阶段，圆形种植区在空间分布图中逐渐显现，蒸散量迅速升高；7月植被呼吸作用最强，牧草实际蒸散量提升至最大值；8~9月进入生长季后期，非种植区蒸散量快速下降，牧草种植区蒸散量逐渐降低。

Figure 5. Spatiotemporal distribution and frequency distribution histogram of actual monthly evapotranspiration in the growing season of 2014

图5. 2014年生长季月实际蒸散量时空分布和频率分布直方图

Figure 6. Bar chart of daily precipitation (blue) and line chart of daily actual evapotranspiration (orange) during the 2014 growing season

图6. 2014年生长季(蓝色)日降水量柱状图和(橙色)日实际蒸散量折线图

本文通过锡林浩特气象站2014年5~9月降雨数据来探究区降雨与ET_a的关系，如图6所示。日降水量最高(27.6 mm)发生在2014年6月16日(第170天)，降水量和降水频率在研究区生长季中期(6~7月)较多，在早期和后期(5月和9月)较少，这与牧草ET_a变化规律基本一致。明显可以观察到降水事件发生时，牧草ET_a下降；而未发生降水的时段，牧草ET_a又逐渐升高，这与空气湿度对ET_a响应的自然规律基本一致。

3.2. pySEBAL模型精度分析

本文以锡浩特气象站气象数据通过第二章第二节中Penman-Monteith算法公式(13) (14)，检验pySEBAL模型反演人工灌溉牧草实际蒸散发的精度。选择研究区模拟结果均值与Penman-Monteith算法结果进行比较，结果如图7所示，R² = 0.7504，RMSE = 1.2575 mm/d，MRE = 0.9366 mm/d。与Penman-Monteith算法结果相比pySEBAL模型略微高估了牧草生长季ET_a，但模型总体能够比较准确反演人工灌溉牧草ET_a变化趋势。

Figure 7. Accuracy verification results

图7. 精度验证结果

4. 讨论与结论

通过分析研究区灌溉牧草生长季实际蒸散量时空变化可以得出：研究区实际蒸散量在生长季初期(5~6月)蒸散量逐渐升高，中期(7月)到达顶峰，后期(8~9月)逐渐降低。研究区牧草种植边界在5~7月逐渐清晰，又从7~9月缓慢消失，这表明受灌溉和其他因素的影响，研究区种植区牧草与非种植区杂草的实际蒸散量有很大差异，从侧面反映出灌溉对牧草生长发育的影响巨大，对提高产量有显著效果。

pySBEAL模型能够反映牧草实际蒸散量在整个生长季时空变化规律，与降水事件发生也具有很高相关性。该模型反演的研究区实际蒸散量均值与Penman-Monteith算法计算结果相比R²达到了0.75，说明两种算法的结果有很高的相关性；RMSE = 1.26 mm/d、MRE = 0.94 mm/d体现出两组结果的之间的差距较小，表明pySEBAL模型能以较高精度反演人工灌溉牧草实际蒸散量。

SEBAL模型应用于估算作物ET_a方面成功的案例有很多，使用场景也十分广泛。Jingyuan Xue等人利用pySEBAL模型在加利福尼亚州农村分别反演番茄、马铃薯、杏树苗的实际蒸散量 [13] ；苏婷婷等人基于SEBAL模型估算土默特平原农田蒸散量 [14] ，都具有较高精度。李宝富等人也成功应用SEBAL模型估算塔里木河干流域蒸散发 [26] 。

本研究仍存在一些不足之处，由于2013~2021年期间研究区内仅有25景影像符合无云或少云的要求，影像数量较少可能会对研究结果造成偏差。2014年各月只有一景影像可用，较少的作物系数插值可能会出现较高的不确定性。另外，研究区非种植区也占有很大面积，在计算区域平均值时会导致ET_a偏低，使用实测蒸散量数据在点位上验证可能会得到更高的相关性。

本研究基于Landsat-8遥感影像，结合锡林浩特市气象站气象数据，利用pySEBAL模型对锡林浩特市沃原奶牛场种植基地人工灌溉牧草从实际蒸散量时空变化、和模型精度验证两方面分析pySEBAL模型估算锡林郭勒牧场生长季蒸散量适用性。得出以下结论：

(1) 研究区牧草蒸散量随生长季时间变化明显，ET_a在0~10 mm/d和50~200 mm/m之间。5~6月研究区牧草实际蒸散量逐渐增加，7月达到峰值(5.46 mm/d、155.71 mm/m)，8~9月逐渐减少。

(2) 将Penman-Monteith算法计算结果作为验证值对pySEBAL模型反演实际蒸散量结果进行精度验证分析，R² = 0.75、RMSE = 1.26 mm/d、MRE = 0.94mm/d，表明在研究区使用pySEBAL模型反演实际蒸散量可信度较高，能够描述研究区实际蒸散量时空变化情况。

目前对灌溉牧草蒸散量的研究较少，但随着农业用水管理制度的完善，以及越来越多牧场通过灌溉来增加牧草产量，牧草灌溉管理在未来会进一步受到重视。本研究基于遥感模型反演实际蒸散量相较于地面测量具有更好的时效性，以及更节省人力物力。这为今后的牧草灌溉管理提供了新的思路，为政府制定合理的政策提供有力依据。