蒸散发（evapotranspiration，ET）主要包括地表土壤水分的蒸发和植被的蒸发蒸腾^[1]，作为全球水循环的重要组成部分，它是连接水分运移和能量平衡的重要纽带，也是造成土壤水分损失、导致干旱的一个主要因素。

遥感方法可以实现区域尺度蒸散量的估算^{[2 − 3]}，遥感蒸散发数据产品的出现，也为学者的研究带来了便利^{[4 − 5]}。但是，由于遥感影像的重访周期较长，且通过遥感手段估算的地表蒸散量一般都是卫星过境时刻的瞬时值，难以获得连续时间序列的蒸散量，因此通过遥感手段估算蒸散量的方法还存在一定的局限性。

蒸发强度不仅受到太阳辐射、温湿度、风速和植被等多种因素的影响^{[6 − 8]}，土壤水分作为蒸散发的直接供给来源，在蒸散发的过程中也起着极其重要的控制作用^[9]，土壤水分的胁迫导致实际蒸散发强度低于潜在蒸散发强度。因此，结合土壤水分数据进行蒸散量的遥感估算，在改进蒸散量的估算精度方面具有重要意义^[10]。

联合国粮农组织（Food and Agriculture Organization of the United Nations，FAO）在P-M公式的基础上提出了单作物系数法，将参考作物蒸散量乘以作物系数、土壤水分胁迫系数推求实际蒸散量^{[11 − 12]}，结合遥感数据可以实现该方法在区域尺度上的扩展。余昭君等^[13]基于单作物系数法，结合气象条件对葡萄中期和后期的作物系数进行了校正，计算了葡萄园的蒸散量，结果与波文比能量平衡法测定的结果具有很好的一致性。李毅等^[12]建立单作物系数与砂石覆盖量的关系，对砂石覆盖下的冬小麦的蒸散量进行了计算，获得了良好的效果。但是目前对单作物系数法的研究多集中在作物系数方面，对土壤水分胁迫系数的研究仍有很大的发展空间。

Richards方程是土壤水分运移的基本方程，可以模拟不同条件下连续时间尺度的土壤水分运动过程，为研究土壤水分对蒸散发的影响机理提供了理论依据。张堂堂等^[14]基于此方程结合遥感数据对黄土高原典型塬区的蒸散量进行了估算，为区域蒸散量的估算提供了一种新的思路。涡度相关法计算蒸散量被认为是各种方法中较精密而可靠的方法，王韦娜等^[15]将蒸渗仪法和涡度相关法测定的蒸散量进行对比，发现两者具有较高的相关性。在蒸散发的相关研究中，常将此作为一种标准参考方法。

针对目前蒸散发研究在土壤水分胁迫和连续性方面的不足，本文旨在以涡度相关法计算的蒸散量作为实际蒸散量，结合FAO推荐的单作物系数法，将土壤含水量信息引入P-M公式，实现土壤水分胁迫下的日蒸散量估算。同时基于Richards方程，模拟蒸发条件下连续时间尺度的一维土壤水分运动，实现连续日蒸散量的估算，并结合遥感数据进行区域尺度的扩展，为区域蒸散量的估算以及揭示水分的传输与蒸散发过程提供参考。

1.   研究区概况

黑河流域位于河西走廊中部，发源于祁连山中段，流域气候干燥，降水稀少且集中，多大风，日照充足，昼夜温差大。年均降水量为113 mm，年均气温为8 °C，年均蒸散量为407 mm。本次研究区为黑河流域处于甘州区境内的部分区域（图1）^{[16 − 17]}，研究区面积为3150.57 km²，南北高中间低，中部多为耕地且河流分布较为密集，南北两侧分布多为林地和草地。研究区内包含涡度相关仪站点大满超级站。该站点位于甘肃省张掖市大满灌区内，生长季下垫面种植作物为玉米，为排除降雨和生长季农业灌溉的影响，研究时间段为2020年4月27日—5月3日。

图  1  研究区范围及土地利用类型

Figure  1.  Study area and land use type

下载: 全尺寸图片 幻灯片

2.   材料与方法

2.1   数据

本研究数据来自国家青藏高原科学数据中心（https://data.tpdc.ac.cn/home），包括大满超级站气象要素梯度观测系统的气象数据^{[18 − 19]}、3个站点涡度相关仪观测的通量数据^{[20 − 22]}、高分辨率国家土壤信息网格基本属性数据集^{[23 − 24]}、遥感土壤水分数据集^[25]和地表日蒸散量数据集^[26]。

其中气象数据包括空气温湿度、风速和降雨等，采样时间间隔为10 min；通量数据包括感热通量、潜热通量、净辐射通量和土壤热通量，采样时间间隔为30 min；土壤数据包括地下埋深5，15，30，60 cm的土壤容重和土壤颗分数据以及地下埋深2，4，10，20，40 cm的土壤体积含水量数据以及1 km分辨率逐日全天候地表土壤水分遥感数据集。

2.2   研究方法

2.2.1   蒸散量计算

本研究采用FAO改进的P-M公式对潜在蒸散量（potential evapotranspiration，PET）进行估算：

式中：${I}_{\mathrm{P}\mathrm{E}\mathrm{T}}$——潜在蒸散量/（mm·d⁻¹）；

$\Delta$——蒸汽压曲线的斜率/（kPa·°C⁻¹）；

${R}_{\mathrm{n}}$——作物表面的净辐射/（MJ·m⁻²·d⁻¹）；

G——土壤热通量/（MJ·m⁻²·d⁻¹）；

T——平均温度/°C；

$\gamma$——干湿计常数/（kPa·°C⁻¹）；

${u}_{2}$——2 m高度处的风速/（m·s⁻¹）；

${e}_{\mathrm{s}}$——饱和水汽压/kPa；

${e}_{\mathrm{a}}$——实际水汽压/kPa。

涡度相关技术能通过测量各种属性的湍流脉动值直接测量通量。通用的方法多是建立在经验关系的基础上，或从其它气象参量推论而来，该方法不受平流条件限制，具有较高的准确性和可靠性。本研究以涡度相关法计算的蒸散量作为实际蒸散量：

式中：${I}_{\mathrm{E}\mathrm{T}\mathrm{V}}$——涡度相关法计算的蒸散量/（mm·d⁻¹）；

${I}_{\mathrm{L}\mathrm{E}i}$——潜热通量观测值的校正值/（W·m⁻²）；

$i$——潜热通量的观测次数；

$\lambda$——水的汽化潜热/（J·kg⁻¹），为2.5×10⁶ J/kg。

2.2.2   土壤水分胁迫下的实际蒸散量估算

在实际应用中，由于土壤水分和作物生长状态等的影响，实际蒸散量低于潜在蒸散量。FAO 推荐单作物系数法计算土壤水分胁迫条件下的实际蒸散量：

式中：$I_{\mathrm{E}\mathrm{TA}}$——土壤水分胁迫下的实际蒸散量/（mm·d⁻¹）；

${K}_{\mathrm{c}}$——作物系数；

${K}_{\mathrm{\omega }}$——土壤水分胁迫系数。

将实测土壤含水量数据作为胁迫因子引入，实现土壤水分胁迫下的实际日蒸散量的估算。由于所选时间段为连续日尺度的短时间段，作物生长状态不具有较大的变化，所以认为作物系数${K}_{\mathrm{c}}$为常数，因此将作物系数法中的${K}_{\mathrm{c}}$和${K}_{\mathrm{\omega }}$统一考虑为土壤含水量的函数$f\left(\omega \right)$：

式中：$\omega$——土壤体积含水量/（cm³·cm⁻³）；

$f\left(\omega \right)$ ——土壤体积含水量函数；

d——常数/（mm·d⁻¹）。

针对土壤体积含水量函数，本研究采取线性、指数、对数、Sigmoid函数4种形式进行拟合：

2.3   蒸发条件下的Richards方程

土壤水分运移的基本方程大多采用Richards方程。Richards方程是根据Darcy定律和质量守恒方程推导而来^{[11, 27]}，是一个二阶非线性偏微分方程，很难求得解析解，常用数值解法模拟土壤水分运动过程。当蒸发强度已知时，蒸发条件下的一维土壤水分运动的数学模型为：

式中：$K\left(\omega \right)$——土壤导水率/（m·s⁻¹）；

$D\left(\omega \right)$——土壤水分扩散率/（m²·s⁻¹）；

$E\left(t\right)$ ——水分蒸发强度/（mm·d⁻¹）；

${\omega }_{\textit{z}}$——不同深度的初始土壤体积含水率/（cm³· cm⁻³）；

$t$——模拟的时间步长/s;

z——土壤的深度/cm;

L——蒸发下边界的深度/cm。

式（9）可以建立矩形差分网格，将计算区域离散化，利用追赶法求解，可得到土壤剖面上诸点的$\omega$值，按时间层推进，即可模拟土壤水分在剖面上的动态变化。

3.   结果与分析

3.1   潜在/实际蒸散量

本文将涡度相关仪所观测到的潜热通量进行筛选、插值以及能量强制闭合，以此计算大满站的实际日蒸散量。同时用FAO推荐的P-M公式计算对应日期的潜在日蒸散量。对比4月27日—5月3日的潜在日蒸散量和实际日蒸散量（表1）可以看出，实际蒸散量普遍小于潜在蒸散量，潜在日蒸散量约为实际日蒸散量的1.7倍，但两者总体上具有相似的变化规律，相关性达到0.918。

表  1  潜在/实际日蒸散量计算结果

Table  1.  Calculationresults of potential/actual daily evapotranspiration

日期	潜在日蒸散量/mm	实际日蒸散量/mm
4月27日	5.597	4.211
4月28日	5.020	3.040
4月29日	5.175	3.336
4月30日	5.226	3.454
5月1日	5.385	3.462
5月2日	4.522	2.260
5月3日	3.055	1.775

下载: 导出CSV 

| 显示表格

3.2   土壤水分胁迫下的实际蒸散量估算

由于蒸发蒸腾作用首先从土壤表层提取水分，深层土壤水分通过毛细管的过程较慢^[28]，且对蒸散发的影响通过影响表层土壤水分实现。本次仅对地下埋深10 cm以浅的土壤体积含水量进行研究。选取不同埋深处（2，4，10 cm）的土壤含水量的实测值，以及0～4 cm、2～4 cm、4～10 cm、2～10 cm、0～10 cm共5个深度区间土壤含水量的平均值进行分析。为了验证模型的准确性和可行性，选取6月10—13日的数据对模型进行精度验证。根据式（4）进行土壤水分胁迫下的实际日蒸散量估算，并与基于涡度相关仪法计算的日蒸散量进行对比，以均方根误差（root mean square error，RMSE）作为模型的评价指标，结果见表2。

表  2  不同土壤水分胁迫函数构建的蒸散量估算模型精度对比

Table  2.  Accuracy comparison of evapotranspiration estimation models constructed by different soil water stress functions

埋深/cm	均方根误差/（mm·d−1）
	线性函数	指数函数	对数函数	Sigmoid函数
2	0.496	0.433	0.437	0.393
4	0.413	0.329	0.332	0.232
10	0.416	0.331	0.337	0.223
0～4	0.428	0.347	0.350	0.133
2～4	0.424	0.310	0.311	0.157
4～10	0.414	0.330	0.335	0.226
2～10	0.409	0.324	0.327	0.217
0～10	0.407	0.314	0.316	0.195

下载: 导出CSV 

| 显示表格

从不同函数模型的评价指标可以看出，在不同的函数模型中，线性函数模型的精度最低，指数和对数函数模型的精度相当。Sigmoid函数在高值和低值时变化较为缓慢，可以灵活的调节土壤含水量在低值和高值情况下的函数模型，能更准确地反映土壤水分对土壤蒸散发的胁迫作用，因而其相较于其他模型具有更好的拟合效果。经过比较，选择基于0～4 cm深度平均土壤体积含水量（${\omega }_{0-4}$）的Sigmoid函数作为最终的土壤水分胁迫函数，构建蒸散量的估算模型：

3.3   基于Richards方程的连续日蒸散量估算

Richards方程是非饱和渗流理论的基本方程，本研究基于式（9）对蒸发条件下一维垂向的土壤水分运动进行数值模拟。上边界为蒸发边界，蒸发量由式（10）计算；下边界为恒定含水率边界，通过对实测土壤含水量波动情况进行分析，发现地下埋深40 cm处的土壤含水量波动强度为0.17%～0.47%，波动程度较小，因此选择40 cm处为模拟过程的下边界。

具体计算方法为：用初始k时刻的土壤含水量${\omega }^{k}$计算初始蒸发强度，并用k时刻的${\omega }^{k}$值计算K($\omega$)和D($\omega$)作为k+1时刻的预报值，然后用迭代法进行误差校正，求得k+1时刻的${\omega }^{k+1}$值，重复上述步骤即可求得任一时段的土壤体积含水量和蒸散量。

本次研究以4月27日0时的土壤含水量为初始值，用上述方法模拟4月28日—5月3日的土壤水分运动，得到每日0时地下不同节点处（2，4，10，20，40 cm）土壤含水量的模拟值（图2）。从图中可以看出，随着蒸散发的进行，地下埋深40 cm以上的土壤含水量不断减少，且每日剖面土壤含水量的分布具有相同的分布规律。对比模拟与实测效果可以看出，表层的土壤水分运动较为复杂，实测值和模拟值相差较大，深度越深，模拟值和实测值越为接近。此次试验将地表埋深40 cm以内的土壤作为均质进行模拟，但是实际上不同埋深的土壤性质存在一定差异，这也是误差产生的原因。

图  2  4月28日—5月3日0 时土壤含水量模拟值与实测值对比

Figure  2.  Comparison between simulated and measured values of soil water content from April 28 to May 3

下载: 全尺寸图片 幻灯片

结合每日土壤含水量的模拟结果和潜在日蒸散量，根据式（10）计算每日的实际蒸散量，实现顾及土壤水分胁迫的连续日蒸散量的估算。并对每天模拟值与实测值的误差进行分析（表3），结果表明，模拟的日蒸散量与实测值较为接近，平均绝对误差为0.260 mm/d，均方根误差为0.288 mm/d，预测结果较好。

表  3  日蒸散量模拟值与实测值的误差分析

Table  3.  Error analysis of simulated and measured evapotranspiration

日期	实测值/mm	模拟值/mm	绝对误差/mm
4月28日	3.040	3.191	0.151
4月29日	3.336	3.037	−0.299
4月30日	3.454	3.684	0.231
5月1日	3.462	3.917	0.455
5月2日	2.260	2.185	−0.074
5月3日	1.775	1.428	−0.346

下载: 导出CSV 

| 显示表格

3.4   区域尺度日蒸散量的遥感反演

基于点尺度的研究，结合遥感数据实现区域尺度上的扩展。研究过程中所需水力学参数（$\omega_{\mathrm{s}}、\omega_{\mathrm{r}}、K_{\mathrm{s}}$和$\alpha、n$）由土壤的容重、砂粒、粉粒和黏粒含量数据通过Rosetta模型预测得出。土壤含水量数据采用赵天杰等^[25]提供的“中国1 km分辨率逐日全天候地表土壤水分数据集（2003—2022）”。该产品是通过在36 km处对基于AMSR-E和AMSR-2的SSM进行降尺度生成的。

本文构建模型用到的土壤含水量为地表0～4 cm土壤含水量的平均值，而遥感反演的是土壤表层的土壤含水量，因此在使用之前需要对土壤水分数据集的适用性进行验证。将大满站点的实测土壤含水量与土壤水分数据集的结果进行对比分析（图3）。结果显示，土壤水分数据集的结果略小于实测土壤含水量，但是两者具有显著的相关性，相关系数达到0.922；线性拟合的拟合优度也达到0.849，且两者的均方根误差为0.0445 cm³/cm³，在该数据集的验证误差范围之内。因此，该数据集可用于本次研究区域尺度的扩展。

图  3  实测土壤含水量和土壤水分数据集拟合关系

Figure  3.  Fitting relationship between measured soil water content and soil water data set

下载: 全尺寸图片 幻灯片

以4月27日的土壤含水量数据作为初始数据，结合Richards方程和式（10）对蒸发条件下区域尺度的土壤水分运动进行模拟，并实现区域尺度蒸散发的遥感反演，得到甘州区内黑河流域4月28日—5月3日的日蒸散量的反演结果（图4）。

图  4  4月28日—5月3日研究区日蒸散量估算

Figure  4.  Estimation of daily evapotranspiration in the study area from April 28 to May 3

下载: 全尺寸图片 幻灯片

结果显示，4月28日—5月3日研究区的日蒸散量具有相同的空间分布特征，即蒸散量的高值都出现在研究区的中部，西南和东北地区的蒸散量相对较小。研究区地势南北高中间低，河流多向中部聚集，导致中部地区土壤含水量较高，蒸发强度较大。在时间序列上，研究时间段内4月30日的日蒸散量最大，5月3日的日蒸散量最小，通过与点尺度日蒸散量的实测结果进行对比分析，发现区域尺度上连续日蒸散量的遥感反演结果与站点的实测结果在时间序列上表现出一致的变化规律，也验证了此方法的可行性。

4.   讨论

4.1   区域蒸散发遥感反演精度评价

结合研究区内大满站、花寨子站和张掖湿地站的观测数据对区域尺度日蒸散量遥感反演结果的精度进行验证，结果见表4。

表  4  三个站点日蒸散量遥感反演结果的精度评价

Table  4.  Accuracy evaluation of remote sensing inversion results of daily evapotranspiration at three stations

站点	实测平均蒸 散量/（mm·d−1）	遥感反演平均蒸 散量/（mm·d−1）	平均绝对误 差/（mm·d−1）	均方根误差 /（mm·d−1）
大满站	2.787	2.907	0.283	0.287
花寨子站	1.732	1.599	0.334	0.308
张掖湿地站	4.881	5.070	0.307	0.343

下载: 导出CSV 

| 显示表格

花寨子站的下垫面为荒漠，对应的日蒸散量为3个站点中最小；张掖湿地站的下垫面为湿地，对应的日蒸散量最大。站点的土壤含水量与日蒸散量表现出正相关的关系。

选取姚云军等^[26]提供的“祁连山地区基于MODIS的逐日地表蒸散发数据（2020）”对本文区域尺度日蒸散量遥感反演结果的空间分布进行验证。对研究区外接一个矩形，并将其等间隔的分为55行×55列，取格网点为验证的数据采样点。通过对比发现，两者的相关系数达到0.940（图5），日蒸散量具有相似的空间分布，高值集中在研究区中部的水域和耕地区域（由于城市用地多为不透水面，蒸散发过程复杂，本文不对其进行分析）。不同土地利用类型的平均日蒸散量大小为：水域>耕地>林地>草地>未利用土地（表5）。

图  5  遥感反演日蒸散量和基于MODIS的逐日地表蒸散发数据对比

Figure  5.  Comparison of daily evapotranspiration retrieved from remote sensing and daily surface evapotranspiration data based on MODIS

下载: 全尺寸图片 幻灯片

表  5  不同土地利用类型下遥感反演日蒸散量与基于MODIS的逐日地表蒸散发数据的精度评价

Table  5.  Accuracy evaluation of daily evapotranspiration retrieved by remote sensing and daily surface evapotranspiration data based on MODIS under different land use types

土地利用 类型	基于MODIS的平均 蒸散量/（mm·d−1）	遥感反演平均蒸 散量/（mm·d−1）	平均绝对误 差/（mm·d−1）	均方根误 差/（mm·d−1）
耕地	3.152	3.233	0.157	0.173
林地	2.281	2.426	0.127	0.143
草地	1.971	1.980	0.103	0.153
水域	5.771	5.571	0.114	0.165
未利用土地	1.571	1.620	0.133	0.225

下载: 导出CSV 

| 显示表格

4.2   模型适用性分析

本文基于大满站农田下垫面进行了研究，但是在花寨子站，下垫面为荒漠，植被覆盖较少，土壤含水量较低，模型的均方根误差由大满站的0.287 mm/d下降至花寨子站的0.308 mm/d（表4）。在张掖湿地站，下垫面为芦苇湿地，土壤含水量较高，植被覆盖更为密集，模型的均方根误差为0.343 mm/d。当下垫面发生变化时，土壤的性质发生了变化，尽管模型的估算精度出现下降，但整体仍维持在较为理想的水平，因为Richards方程提供了土壤水分运动的理论基础，可以模拟土壤中不同条件下的水分运动，为单作物系数法提供了更为准确的土壤水分状态，且模型中的Sigmoid函数在土壤水含量较低和较高时变化较为缓慢，可以更好地表达不同土壤水分条件对蒸散发的胁迫作用。

通过与MODIS地表日蒸散发数据集的对比，发现基于Richards方程的单作物系数法可以适应不同类型的土壤和植被，在不同植被覆盖类型下以及水域都有较好的估算结果（表5），均方根误差维持在0.172 mm/d左右。但是模型对土壤和气象数据的要求较高，且Richards方程对土壤参数和初始条件较为敏感，使得在一些地区或数据受限的情况下限制了此方法的应用。此外，当进行更长时间序列的蒸散量的估算时，还需要考虑植被的生长状态对结果的影响。

5.   结论

（1）基于P-M公式计算的潜在日蒸散量与基于涡度相关法计算的实际日蒸散量具有相似的变化规律，两者的相关性为0.918。由于受到土壤水分等因素的限制，实际蒸散量普遍小于潜在蒸散量，潜在日蒸散量约为实际日蒸散量的1.7倍。

（2）通过构建实际蒸散量与土壤含水量和潜在蒸散量之间的线性、指数、对数、Sigmoid函数关系，发现引入土壤含水量信息后的P-M公式能更准确地估算日尺度的蒸散量，且基于0～4 cm深度平均土壤体积含水量的Sigmoid函数作为土壤水分胁迫函数所构建的蒸散量估算模型精度最好，估算的日蒸散量的均方根误差为0.133 mm/d。

（3）用Richards方程模拟蒸发条件下一维垂向土壤水分的运移过程，研究时间段内，土壤水分模拟值的平均均方根误差为0.0299 cm³/cm³。用每日土壤水分的模拟值结合每日的潜在日蒸散量，计算土壤水分胁迫下每日的实际日蒸散量，所估算的日蒸散量与实际日蒸散量较为接近，均方根误差为0.288 mm/d。

（4）基于点尺度上的研究结果，结合遥感土壤水分数据和土壤容重数据实现区域尺度土壤水分胁迫下的日蒸散量的估算。受地形南北高、中间低的影响，研究区日蒸散量的高值集中在中部的水域和耕地区域，西南和东北地区的蒸散量相对较小，不同土地利用类型下的平均日蒸散量水域>耕地>林地>草地>未利用土地，且日蒸散量在时间尺度上的变化规律与站点的实测结果保持一致。

由于下垫面的异质性以及不同深度土壤性质的差异，今后仍然需要结合下垫面的类型和不同深度土壤的属性，同时考虑植物的生长状态，对复杂下垫面下区域尺度日蒸散量的遥感估算进行深入探讨。