基于地表水循环遥感观测的黑河流域水平衡分析

闫柏琨，李文鹏，甘甫平，郑跃军，祁晓凡，白 娟，郭 艺，吴艳红，王龙凤，马燕妮

（1.中国自然资源航空物探遥感中心,北京 100083；2.自然资源部航空地球物理与遥感地质重点实验室,北京 100083；3.中国地质环境监测院,北京 100081；4.中国科学院空天信息创新研究院,北京 100094）

西北干旱内流区甘肃黑河流域南起祁连山脉、北至荒漠戈壁，孕育了上游寒区高山草甸、中游农田绿洲、下游荒漠绿洲的生态系统格局，是河西走廊主要经济、文化、生态走廊，是国家重点粮油、生态、经济保护区。但中下游干旱少雨、中游农业开发水资源消耗量较大，水资源矛盾突出。由于中游水土资源开发利用程度提高，造成了下游湖泊萎缩消亡，西居延海1961年干涸，东居延海1992年干涸。为缓解下游水资源紧张状况，从2000年开始通过水量统一调度实施黑河分水方案，规定了不同保证率条件下进入下游的水量，从2003年开始东居延海实现连年不干涸，水域面积常年维持在36 km2以上，下游生态系统健康状况明显好转[1-5]。

2000年以来，黑河流域自然与社会水文条件经历了很大变化。由黑河干流上游径流数据（莺落峡水文站）可知，2000年后，年径流量总体呈显著增加趋势[6-8]，2006—2018年期间年平均径流量较往年增加了27.4%[9]。由于上游径流量显著增加[10-12]，向下游分水量增加，这是下游缺水得以缓解的重要原因之一。2001年实施水资源治理工程以来，流域下游狼心山水文站年径流量总体呈增加趋势，气候变化是主因，治理工程的实施也发挥了重要的积极作用[11]。尽管如此，尚未达到国务院97 分水方案的要求，而且未来转入枯水年或平水期，上游径流量可能减少，向下游的下泄水量势必下降，下游天然生态健康存在较大的逆转风险[4]。黑河流域未来是否仍面临“中下游用水矛盾凸显的可能”，需加强流域水循环研究与水平衡分析，促进中下游水资源合理配置与协调发展。

地表水循环遥感观测技术快速发展，正逐步成为流域水平衡分析、水资源合理配置的重要手段。黑河流域具有西北内流盆地特有的水文条件，可作为水循环遥感观测技术研发应用的重要基地。多手段集成的天—空—地综合观测试验[13]、陆面蒸散发与水面蒸发估算方法研究[14-17]、山区水文模拟[18]、农田水文模拟[19-20]等水文生态过程监测、模拟、分析方法的研究，促进了流域水循环监测技术的进步。目前尽管各种地表水循环遥感观测技术快速发展，但仍以单项要素反演方法研究为主。本次研究以流域水资源问题为导向，综合应用多种遥感观测技术，开展地表水循环分析与水平衡分析。

本文以2000—2019年黑河流域水文显著变化期为研究时段，综合应用TRMM 与GPM 卫星数据观测的降水量、遥感估算的蒸散发量数据，结合气象站点、水文站点等观测数据进行了研究，以期对流域地表水资源时空变化特征进行分析，为地表水资源协调利用提供依据，促进地表水循环遥感观测技术的发展，明确技术瓶颈与发展方向。

1 研究区概况

黑河流域地处青海、甘肃、内蒙古三省（区）交界，是我国第二大内陆河，发源于南部青海祁连山中段，北至中蒙边界（图1）。

图1 黑河流域范围及位置Fig.1 Location of the Heihe River Basin

根据地表水力联系可分为东、中、西3 个相对独立的子水系。东部子水系包括黑河干流、梨园河及20多条支流。中部子水系为酒泉马营河至丰乐河诸小河流水系，为浅山短流，归宿于肃南县明花区至高台盐池盆地。西部子水系为酒泉洪水坝河至讨赖河水系，多为浅山短流，只有洪水坝河与讨赖河可贯穿酒泉盆地，讨赖河经嘉峪关市后改称北大河并经鸳鸯池水库进入北部金塔盆地。

黑河干流发源于青海省祁连县，从祁连山发源地到东居延海，全长约928 km。莺落峡以上为上游，海拔高，气候严寒湿润，年降水量为250～500 mm，年蒸发量为700～800 mm，是全流域的产流区。莺落峡—正义峡之间为中游，绿洲、荒漠、戈壁、沙漠断续分布，地势平坦，年降水量为110～370 mm，年蒸发量为1 200～2 200 mm，是河西走廊重要的灌溉农业区。正义峡以下为下游，地势开阔平坦，分布有东居延海等湖盆洼地和广阔的沙漠戈壁，属荒漠干旱区和极端干旱亚区，年降水量为40～54 mm，年蒸发量为2 200～2 400 mm[21]。

2 数据与方法

2.1 数据源

数据源主要包括气象站点数据、卫星数据（GPM降水卫星、MODIS 陆地观测卫星）、陆面模式产品、土地覆盖产品、径流数据（表1）。

表1 主要数据源及说明Table 1 Descriptions of the main data

2.2 研究方法

以水循环要素遥感观测为基础，进行了不同分区水量平衡关系与地表水循环分析、黑河干流上游径流变化原因分析、黑河中游主要土地覆盖类型蒸散发水量消耗等水平衡分析。研究方法与流程主要包括：

（1）综合降水卫星数据及气象与遥感数据，估算月尺度降水量、陆面蒸散发量、水面蒸散发量、潜在蒸散发量等水循环观测参量；

（2）根据河流水系、地形高程、地下水分区及水文站点位置进行水循环分区的划分，并明确各分区的水量平衡关系；

（3）计算不同水循环分区的多年平均水循环量（流入量、流出量、降水量、蒸散发量）；

（4）基于Budyko 水热模型计算黑河干流上游降水、潜在蒸散发对径流变化的贡献率，并分析原因；

（5）计算黑河中游主要土地覆盖类型的蒸散发耗水量；

（6）针对黑河流域的现状，重点探讨了黑河干流2000年以来径流量增加的原因及可持续性、中游农业用水与下游湖泊蓄水矛盾、蒸散发量估算的不确定性及原因。

3 结果

3.1 水循环遥感观测

3.1.1 降水量

根据卫星降水（TRMM/GPM）与站点降水对比研究，发现在月尺度上卫星降水与站点降水存在很好的线性关系，可达0.96[23]。由于流域内国家基本气象降水站点稀疏（有4 个站点），且降水量空间变异大，内插误差大，采用线性校正法以全国基本气象站月降水量数据为基准逐月校正生成了2000—2019年全国卫星月降水量数据并裁剪得到黑河流域月降水量，发挥了站点数据单点精度高与卫星数据区域分布精度高的优势。对比同期TRMM、GPM 与站点监测降水量的相关性，发现GPM 数据融合了雷达降水探测精度高与红外高轨卫星数据时间分辨率高的优势，相关性更高。2000年1-5月由于无GPM 数据，采用TRMM数据，2000年6月—2019年12月采用GPM 数据。与流域内4 个站点月降水量数据对比，二者一致性高，均方根误差为8.39 mm（图2）。全流域2000—2017年卫星降水数据（经站点降水数据校正后）与0.5°×0.5°网格站点插值降水数据多年平均降水量分别为88.86，89.14 mm，二者各年相对偏差最小值、最大值、平均值分别为1.00%、20.23%、8.80%。

图2 卫星月降水量与站点实测月降水量对比Fig.2 Comparison between the monthly measured precipitation and monthly satellite precipitation

3.1.2 陆面实际蒸散发量

基于互补相关模型估算了陆面实际蒸散发，该方法认为在地表水分供应不充足时潜在蒸散发越大，则实际蒸散发越小，当地表水分供应逐渐充足时，二者趋于相同。基于新发展的非线性互补相关法，估算了2000—2019年月陆面蒸散发，经与13 个涡度相关站点实测数据对比，全国范围内均方根误差为4.9～16.2 mm[24]。考虑到陆面蒸散发过程影响因素多，估算精度低于降水量的估算精度，假设在全流域尺度上多年平均降水量等于多年平均蒸散发量，以降水量为基准对蒸散发数据进行核校，系数为1.22，二者对比见图3。

图3 流域年降水量、陆面年蒸散发量变化Fig.3 Changes of annual precipitation and land surface evapotranspiration in the basin

3.1.3 湖泊水面蒸发量与流域潜在蒸散发量

水面蒸发量估测方法有实测法、模型估算法。如果蒸发皿直径较小或架设于陆面之上，因蒸发皿与周围环境的蒸发气象条件差异显著，测量的蒸发量与湖泊实际蒸发量差别较大，可高达40%，本文采用架设于开阔水面的大型蒸发皿（如E601）实测的蒸发数据，可准确测量湖泊实际蒸发量[16]。本文采用FAO 参考作物蒸散发模型模拟开阔湖面蒸发[25]，经与湖面E601蒸发皿2014年与2015年4—9月实测数据对比，二者一致性高，均方根误差为39.33 mm，相对于湖面1 178.9～1 183.7 mm/a 的蒸发量，相对误差约为3.3%。

根据每年的湖面面积（基于LandSat、Sentinel-2 中分辨率系列卫星提取）与湖面蒸发量模拟值，计算了多年湖面蒸发量，多年年均蒸发量为0.65×108m3，其中2000—2002年因湖面干涸、2012年因水面分布数据缺失，无水面蒸发量数据（图4）。

图4 东居延海年蒸发量Fig.4 Annual evaporation of the Eastern Juyan Lake

3.2 不同分区水量平衡关系

为了分析流域内水资源平衡，需对地表水循环进行分区，以便于统计分析各分区内水资源的补给、消耗、排泄量。分区主要考虑因素为：（1）区内产汇流及水资源消耗条件相对趋同；（2）现有水文站点分布位置。

产流条件分析的主要依据是地形、高程，分区边界需与自然分水岭一致。

汇流条件分析的主要依据是水系、流向、连通关系，分区需包含主干河流及其各汇入支流。

水资源消耗条件分析的主要依据为：（1）蒸发的地理与气象条件，避免分区大面积横跨山区与平原区，因为山区与平原区气象条件（辐射、气温、气压、风速、湿度）差异较大，潜在蒸散发量差别明显，水循环条件不同；（2）因为该流域内农业灌溉用水、农田蒸发、作物蒸腾量较大，分区尽可能涵盖完整的大型农业灌溉区。

此外，水文站点分布是分区的必要约束因素，分区边界尽可能与重要水文站重合。

按照以上分区原则，黑河流域内划分出15 个分区：祁连山区5 个，山前平原区3 个，下游荒漠戈壁区7 个（图5）。各分区水资源补给、产流、排泄、消耗、水储量及变化等，各分区水资源流向关系，水平衡公式详见表2。

图5 黑河流域地表水循环分区方案Fig.5 Scheme of surface water cycle zoning in the Heihe River basin

流域水资源平衡分析中，水储量变化是需要考虑的变量。对于山区，无大型的储水盆地，地下水主要通过河川基流的方式排泄，且黑河上游祁连山区无大型农灌区抽取山间盆地地下水，可认为水储量在多年尺度上不变。流域中下游盆地地下水储变量是根据2000—2019年地下水水位测量数据求算的平均年度变化量。土壤中含水量多年尺度也可以认为保持不变[26]。

3.3 不同分区水循环分析

因不同分区水文地理条件不同，遥感反演的蒸散发量偏差不同，需基于各分区水平衡公式（表2）对各分区蒸散发量进行校准。

表2 黑河流域地表水循环分区说明Table 2 Descriptions of the zones of surface water cycle

（1）将分区1 水平衡计算的蒸散发与遥感蒸散发量之比作为所有山区（分区1～4、13）遥感蒸散发校正系数，将山区校正前后差值部分计入所有平原区（分区5～12、14～15）蒸散发量，得到校正后的平原区蒸散发量。

（2）将分区12 水平衡计算的蒸散发（假设分区12 的多年平均降水量与蒸散发量相同）与遥感蒸散发量之比作为分区11、12 遥感蒸散发校正系数，将分区5 水平衡计算的蒸散发与遥感蒸散发量之比作为分区5、6、14 遥感蒸散发校正系数（依据是分区11 与12，分区5、6 与14 的地形高程相似，蒸散发估算的偏差量及偏差趋势相似），根据校正系数对各区蒸散发量进行校正。

（3）将上一步校正前后差值部分按照剩余分区面积比，计入其余各分区。最终各分区的多年平均水循环量见表3。

由于黑河干流上游永久冰川分布面积较小、冰川消融在径流中的占比较小[27]，且缺乏冰川消融量准确数据，所以分析中忽略该量。流域7 号水循环区向额济纳盆地（8 号分区）供给水资源约14.25×108m3（7 号分区的流出量），与2015—2019年哨马营水文断面实测年均下泄水量（据黑河流域管理局网站公开数据计算，为10.6×108m3）相近，说明了各区循环量的合理性。

根据各分区多年平均地表水循环量可知（表3），祁连山区是径流主要产流区，向中游下泄约45.11×108m3/a（1～4、13 号分区流出量之和），其中黑河干流为19.00×108m3/a（约占43%），中游是上游来水的主要消耗区，上游来水中约66%（29.92×108m3/a，上游来水减去5 号分区与14 号分区的流出量）用于中游消耗，约34%（15.19×108m3/a，5 号分区与14 号分区流出量）用于补充下游；中游（5～6、14 号分区）消耗水资源约58.68×108m3/a，其中临泽、张掖、民乐一带是中游水资源主要消耗区，年均消耗的上游来水和当地降水量约为43.97×108m3/a（5 号分区蒸散发量，约占中游消耗量的75%），高台、酒泉、嘉峪关一带次之，消耗量约为14.71×108m3/a（6 号与14 号分区蒸散发量之和，约占中游消耗量的25%）。

表3 各地表水循环分区多年平均水循环量（2000—2019年）Table 3 Average surface water cycle pattern in each water division from 2000 to 2019

3.4 黑河干流上游径流变化原因

应用Budyko 水热（降水与蒸发）平衡模型[11，28]对黑河干流上游降水、蒸发、径流关系进行建模，并对径流变化进行归因分析。

式中：Q——流域年径流深/mm；

P——年降水量/mm；

ET0——潜在年蒸散发量/mm；

n——流域下垫面特征的概化参数。

经对比，相对于GPM 卫星降水数据，0.5°×0.5° 网格降水数据输入Budyko 模型模拟的年径流量与实测径流量更接近，因此选用0.5°×0.5° 网格降水数据进行黑河干流上游径流变化的原因分析。黑河流域上游人类活动干扰少，可认为下垫面特征概化参数n保持不变。经计算，n取值0.7 时，2000—2019年实测与模拟径流最为接近，均方根误差为1.96×108m3，相对于2000—2019年年均径流量，相对误差为10.3%（图6），表明Budyko 模型可以很好地解释黑河干流上游径流量的变化。

图6 黑河干流上游年实测与Budyko 模拟径流量Fig.6 Annual measured and Budyko simulated runoff in the upper reaches of the main stream of the Heihe River

对式（1）进行求导，得出年降水、年潜在蒸散发对年径流变化的贡献分别为[11]：

式中：dQP——降水导致的年径流变化/mm；

dQET0——潜在蒸散发导致的年径流变化/mm；

P——年降水量/mm；

B——模拟的实际年蒸散发/mm；

QB——径流深/mm；

ET0——年潜在蒸散发量/mm。

除降水与潜在蒸发外，下垫面、降水强度变化等其它因素均会导致径流量改变。实际径流变化量减去式（2）（3）中dQP、dQET0之和，即为其它因素导致的径流变化量。2000—2019年黑河干流上游，降水量增加是径流量增加的主因，导致年均径流量增加0.35×108m3，约占96%；其它因素是次因，导致年均径流量增加0.014 7×108m3，约占4%，潜在蒸散发对径流量的影响更小，导致年均径流量减少仅为0.001 2×108m3（图7）。

图7 黑河干流上游径流变化原因分析图Fig.7 Attribution of runoff variation in the upper reaches of the main stream of the Heihe River

3.5 黑河中游主要土地覆盖类型蒸散发水量消耗

统计流域中游主要土地覆盖类型区实际蒸散发量，见表4。干流中游的民乐—张掖盆地（5 号分区）农田蒸散发量最大，达372～442 mm；盐池与酒泉盆地次之，为236～330 mm。张掖、盐池、酒泉盆地（分别为5、6、14 号分区）种植结构相似，均以夏玉米为主[29-30]；民乐县以小麦/燕麦为主。民乐—张掖盆地单位面积农田、草地、裸地实际蒸散发比例约为3.2∶2.6∶1，盐池盆地单位面积农田、裸地实际蒸散发比例约为1.9∶1.2∶1，酒泉盆地单位面积农田、裸地实际蒸散发比例约为1.5∶1.2∶1。综合表3、表4 可知，流域中游（张掖、盐池、酒泉盆地）水资源蒸散发消耗中，农田、草地、裸地年均水资源消耗量分别约20.3×108，21.5×108，16.4×108m3，各占消耗总量的35%、37%、28%。

表4 黑河流域中游主要土地覆盖类型蒸散发量Table 4 Evapotranspiration of main land covers in the middle reaches of the Heihe River basin

4 讨论

4.1 黑河干流上游径流量增加的可持续性

前人[6，9]对黑河干流上游径流变化特征与原因已有一定研究，尽管研究时段不完全相同，认识基本一致：（1）2000年以来，径流显著增加，增加幅度达27.4%[9]；（2）径流增加与降水增加、冰川积雪融化、潜在蒸散发减小等自然因素有关，与人类活动及下垫面变化无关。

本文研究进一步表明，降水是径流增加主因，降水对径流增加的贡献率为96%，其它因素只占4%。其它因素包括下垫面变化、人类取水、冻土消融、降水强度等。根据Budyko 模型，黑河干流上游降水量每减少10%，将引起径流量减少14%。如未来年降水量回落至292 mm（2000—2005年平均），莺落峡年径流量将减少至15.8×108m3，相比2015—2019年年均径流量减少5.73×108m3。考虑到降水量的波动变化，黑河干流上游径流量难以维持在较高水平，中下游水资源利用规划应以枯水年为基准。

4.2 黑河中游农业用水与下游湖泊蓄水矛盾分析

根据表3、4，民乐—张掖盆地是黑河中游水资源消耗的主要区域，约占中游75%，其中农田消耗占全盆地的39%，农田水资源消耗以引自黑河干流的灌溉用水为主。因此该盆地农田蒸散发耗水量是分析中游农业用水与下游湖泊蓄水矛盾的关键。根据张掖甘州区农灌区涡度相关仪观测与遥感蒸散发估算结果[31]，农灌区年蒸散发量约为600 mm，除本地年均降水量135 mm 外，另需约465 mm 的年均灌溉用水量。如因黑河干流上游径流量减少至2000—2005年平均水平，且不减少向下游的下泄水量，中游农业灌溉年均亏缺量约为5.73×108m3，每年约1 232 km2的农田灌溉水源无法得到保证。

黑河流域中游农业用水与下游生态保护用水矛盾仍存在较大的显现风险。根据水资源利用多目标规划模型模拟结果[32]可知，只有通过灌区节水措施降低灌溉定额才能有效解决黑河干流中游与下游的用水矛盾。

4.3 蒸散发估算的不确定性与原因探讨

实际蒸散发既是水平衡计算、又是不同类型下垫面耗水量统计的关键数据。黑河流域中游农田蒸散发是流域耗水量统计的基础。实际蒸散发量时空分布合理、不同类型下垫面蒸散发数值比例准确，且与多年降水量相平衡。张掖盆地实际蒸散发量分布（图8）显示，蒸散发空间分布总体合理，裸地、草地、农田蒸散发量由小变大，裸地位于该区北部，紧邻荒漠区，蒸散发量与降水量基本相当，草地分布于农田外围与祁连山北麓之间，水资源供应除降水外，还有山区来水补给，水分供应相对充足，蒸散发量明显增加，人工灌溉农田蒸散发量最大，其中张掖甘州区尤甚。

图8 民乐—张掖盆地2019年实际蒸散发量Fig.8 Actual evapotranspiration of the Minle—Zhangye basin(water resources division 8) in 2019

实际蒸散发受气象、土地覆盖类型、植被长势、土壤水分供应能力等多方面因素的影响，且植被蒸腾、土壤蒸发物理过程复杂。蒸散发量在时间、空间进行全面精确评价超过本文范围。尽管据此得出了相对合理的不同类型下垫面蒸散发耗水的比例数据，但以降水、径流作为约束对蒸散发量进行核校方可完成水平衡分析。这表明实际蒸散发的估算仍是水循环观测中的主要难点。

张掖盆地2019年实际蒸散发量分布（图8）存在2 处估算偏高区，偏高区下垫面类型为分布于孤立山丘之上的非人工林地与草地。全年蒸散发量高于甘州玉米育种灌溉区，原因有待研究。根据互补相关法的基本原理，基于Priestley-Taylor 计算的湿环境（近地表空气湿度饱和）蒸散发量（ETwes）对估算结果影响较大，该数值越大计算的实际蒸散发越大。对照ETwes分布图可以发现在上述孤立山丘处ETwes偏高。由此可见该系数空间变异很大，易受地形影响。因此，Priestley-Taylor 湿环境蒸散发量计算公式中，乘性经验系数在全流域采用统一数值会引入较大误差。

通过对比研究区已公开的代表性的蒸散发方法反演的蒸散发量（表5），有助于进一步理解蒸散发估算的不确定性与原因。区域蒸散发反演方法可分为基于物理模型与基于统计模型2 种。基于物理模型的方法是对蒸散发物理过程进行参数化并估算蒸散发量。基于统计模型的方法需要基于一定数量的实测数据建立蒸散发与诸多影响参量间的统计关系。考虑到基于统计模型的方法机理难辨析，结果难推广，本次仅对基于物理模型的方法与结果进行对比分析。

表5 收集的用于对比分析的流域实际蒸散发量Table 5 Evapotranspiration data collected for comparisons

基于物理模型的方法可分为Penman-Menteith、互补相关、能量平衡3 类。Penman-Menteith（P-M）是一种直接计算蒸散发量的方法，Shuttleworth-Wallace 双源模型是基于P-M 模型估算植被冠层蒸散与冠层下/间土壤蒸发的方法[31,33-34]，PML（Penman-Monteith-Leuning）是以P-M 为基础融合了植被GPP（总初级生产力）的模型[35-36]。互补相关通过“实际蒸散发与潜在蒸散发间存在近似对称/互补关系”现象，根据潜在蒸散发量推算实际蒸散发量[24]。能量平衡法基于“地表感/显热通量、潜热通量、热通量三者之和等于地表净辐射通量”这一能量平衡原理，通过估算地表感/显热通量、热通量间接推算潜热通量与实际蒸散发量[26,37-38]。

多种蒸散发数据对比可发现2 个基本特点（图9）：（1）在全流域、张掖农田灌区（5 号分区）各数据间以及与水平衡估算的蒸散发均存在不可忽略的差异；（2）均存在山区高估，平原区低估的现象，其表现是全流域高于或略低于水平衡估算量，而平原农灌区远小于水平衡估算量。

图9 不同方法估算蒸散发量对比Fig.9 Comparisons of annual evapotranspiration data

由各模型计算结果可知，以P-M 模型为基础的Shuttleworth-Wallace 双源模型反演数据与水平衡估算量接近。

5 结论

（1）祁连山区是径流主要产流区，向中游下泄量约为45.11×108m3/a，其中消耗于中游的水量约为29.92×108m3/a，占66%；补充下游的水量约为15.19×108m3/a，占34%。张掖盆地是黑河中游水资源消耗的主要区域，消耗的上游来水和当地降水量达43.97×108m3/a，约占中游消耗量的75%。中游农田蒸散发量约20.3×108m3/a，占上游来水的45%。降水量增加是黑河干流上游径流增加的主因，对径流量增加的贡献率为96%，导致径流增加0.35×108m3/a。潜在蒸散发对径流增加几乎没有贡献。根据目前黑河干流上游径流量变化与中游水资源消耗现状，如果未来由于水文周期变化导致上游径流减少，中下游用水矛盾凸显的风险较大。

（2）地表水循环遥感观测可作为流域水平衡分析的技术之一，分析流域地表水资源的空间分布状况、揭示水资源变化趋势与原因，支撑水资源合理配置。陆面实际蒸散发量是水平衡分析中主要的不确定性因素，准确估测不同类型下垫面实际蒸散发量是提升分析可靠性的关键。

（3）基于互补相关法估算地表实际蒸散发时，发现用于计算湿环境蒸散发量的Priestley-Taylor 公式中乘性经验系数受地形影响空间变异很大，区域上采用统一数值会对结果造成不可忽视的影响。这是互补相关法在流域水平衡分析应用中需解决的首要问题。