基于多源遥感数据和水量平衡原理的灌溉用水量分析

卢诗卉，赵红莉，蒋云钟，郝 震，张象明，陈根发

（中国水利水电科学研究院，北京 100038）

1 研究背景

随着经济社会发展、人口增长，水资源需求快速增加，水资源短缺问题日益严重［1-2］。农业用水占我国用水总量的60%以上［3］，灌溉用水占农业用水的90%以上。灌溉用水的精确估算对总用水量统计的准确性有着重要的影响，也对支撑水资源精细化管理有重要的意义。但农业用水分布广泛且分散，计量难度相对工业和生活用水较大，区域农业用水总量估算受数据条件限制，较多依靠人工经验，存在较大的不确定性，亟需结合新数据条件研究改善。

现行的灌溉用水估算方法大体分为两类，一是典型调查和定额推演法，二是水量平衡推算法。前者主要根据灌溉定额与实灌面积数据进行估算［4］。灌溉定额的确定首先要进行典型调查，然后进行定额推演，根据《灌溉用水定额编制导则》（GB/T 29404—2012）要求确定省级分区、典型县和水文年型，收集有关数据，整理分析资料，合理调整并确定省级分区主要作物灌溉用水定额。在灌溉定额的确定上，总体存在农业用水计量设施不完善、复杂灌区统计难度大等问题［5］。此外，灌溉用水定额是一个动态指标，但在农业灌溉用水推算的实际工作中，很难按照水文气象变化、田间水分状况及作物长势等进行动态调整，导致该方法不能准确计算农业灌溉用水量［6］。

牛文臣等［7］提出根据区域水量平衡原理计算农业用水量的方法，将耕地和非耕地区域的蒸腾蒸发消耗量作为农业用水量，在长时间尺度上回避区域内部地表水、土壤水和地下水等各部分水量相互转化的复杂关系，通过区域水量平衡关系推求农业用水量，其中耕地区域的蒸腾蒸发量对应为灌溉水量。徐建新等［8］沿用了其对农业用水量的定义，采用水量平衡原理分析计算多个典型灌区的实际农业用水量，并分析影响因素及其变化。这些早期的研究虽然从宏观尺度上给出了基于水量平衡的农业用水匡算方法，但存在定义不准确、没有解释天然蒸散发与灌溉蒸散发关系的问题。詹同涛等［9］基于水量平衡原理提出长序列农田灌溉定额测算方法，给出了相对准确的田间水量平衡关系描述，但计算需要用到农田区域精准的出入境水量和蓄变量数据，只有在试验农田区域才具备计算条件。可见，水量平衡方法虽然在理论上可以推求农业灌溉水量，但在实际应用中，还需要解决水量平衡方程构建与数据可获取性的匹配问题。

遥感技术是最为有效的对地观测技术和信息获取的手段之一［10］。随着各类高空间、时间、光谱分辨率民用卫星的出现，定量遥感技术进一步发展，遥感与地理信息系统、全球导航技术及物联网等技术不断融合，其在农业领域的应用广度和深度不断扩展［11］。遥感蒸散发、农作物长势、重力卫星陆水储量等遥感数据产品不断丰富，有望填补水量平衡方程中大范围耗水、蓄水量变化的监测空白，但由于各类产品存在概念不一致、通用性不强和精度不高等问题，在不断改进遥感产品精度的前提下，还需要研究遥感数据与水量平衡要素在概念上的对应关系，尝试结合多种遥感产品，通过相互校验和组合来支撑灌溉用水估算。

本文尝试通过两级水量平衡方程联合方法来控制多源数据的误差，找到相对合理的蓄变量、耗水量等的数据组合，推求农业灌溉水量。首先在完整区域尺度构建水量平衡方程，利用多源数据推算区域总耗水量，从中分离出与遥感蒸散发概念相对应的耗水量，以此为约束，遴选本区域适用的遥感蒸散发产品；其次在农田区域构建田间水分补给与耗散的水量平衡方程，利用逐日降水和遥感蒸散发数据分析农田降水耗水量与灌溉耗水量，进而估算出灌溉水量。最后以山东省济南市为例，通过收集济南市2012—2015年的降水量、出入境水量、跨流域调水量、工业与生活耗水量、遥感蒸散发量和储水变化量等数据，基于上述方法对济南市2012—2015年灌溉用水量的合理范围进行推算和分析。

2 研究区及数据

2.1 研究区概况本文以济南市为研究区（图1），该市位于山东省中西部，范围为东经116°11′—117°44′、北纬36°01′—37°32′，总面积8151 km2。南依泰山，北跨黄河，分别与西南部的聊城、北部的德州和滨州、东部的淄博、南部的泰安等交界。

图1 研究区概况

研究区地处鲁中南低山丘陵与鲁西北冲积平原的交接带上，地势南高北低。区域内有三大水系，即黄河、小清河以及徒骇河马颊河水系。黄河水系主要有玉符河、北大沙河、南大沙河、浪溪河、玉带河等河流；小清河水系主要有巨野河、绣江河、漯河等河流；徒骇河马颊河水系主要有徒骇河、德惠新河等河流。研究区地处华北中纬度地带，属暖温带半湿润大陆性季风气候区。其主要气候特征是季风明显，四季分明，雨量集中。年平均降水量648.0 mm，降水空间分布不均，总的分布趋势是由东南往西北递减。夏季（6月—8月）降水量为367～499 mm，集中了全年降水量的60%以上。

2.2 研究区数据收集为计算研究区农业用水量，本文收集了研究区及周边2012—2015年60个水文和气象雨量站逐日降水、跨流域调水量、非农业用水量与18个水文监测站点的行政区出入境水量等地面观测统计数据，以及中分辨率成像光谱仪（moderate-resolution imaging spectroradiometer ，MO⁃DIS）8日遥感蒸散发产品MOD16和16日遥感NDVI 产品MOD13、全球陆面同化数据（GLDAS）逐月遥感蒸散发产品GLDAS-NOAH、重力卫星（Gravity Recovery and Climate Experiment，GRACE）逐月陆地水储量数据产品、中国科学院遥感与数字地球研究所提供的2015年土地利用数据。如表1所示：

表1 数据来源

3 研究方法

3.1 基于两级水量平衡联合的灌溉水量估算方法

（1）区域水量平衡。水量平衡原理是指在一定的时空范围内，水分运动保持质量守恒，亦即输入水量与输出水量的差额等于系统内蓄水变化量。建立区域水量平衡方程，如式（1）：

式中：P为时段内区域降水量，可由地面观测数据计算得到；WI为计算时段内流入研究区域的水量，WO为计算时段内流出研究区域的水量，WD为计算时段内区域跨流域调水量，均可由地面径流观测或调水量计量得到；W耗为计算时段内区域内消耗的总水量；ΔW为计算时段内区域的蓄水变化量。

区域耗水量是指区域内自然消耗或人类在生产和生活过程中直接或间接消耗的水量，简单考虑来源主要包括太阳能、矿物能和生物能驱动的耗水量［12］。

式中：ET可以分为有灌溉农田的蒸散发ET农和自然下垫面无灌溉的蒸散发ET自；Wm为人类工业活动驱动的耗水；Wh为人类生活活动驱动的耗水。考虑工业和生活耗水具有较完备的监测计量体系，其耗水量统计相对准确，直接采用《济南市水资源公报》公布的工业与生活耗水量统计数据。

由上述分析可知，ΔW与ET为传统水量平衡方程构建时的数据获取难点。本文选用重力卫星陆水储量数据表示ΔW、遥感蒸散发数据表示ET，完成区域水量平衡方程构建。由于遥感数据存在不同程度的误差和区域适用性，选取以遥感为数据源的ΔW和ET时，须以区域水量平衡为约束，满足平衡方程的ΔW和ET数据组对，才能用于计算农田区域的耗水量ET农。

（2）农田水量平衡。对于农田区域，可以将耗水分解为来自降水的蒸散发和来自灌溉的蒸散发，记为有效降水P有效和有效灌溉水量W有效灌：

有效降水是指在旱作条件下消耗于农作物蒸散过程的降水量［13］，美国农业部土壤保持局所推荐的有效降水量方法是目前比较公认和得到推广的有效降水量计算方法之一［14-15］，有学者将此方法用于京津冀地区取得较好效果［16］，研究区与京津冀地区同属华北平原，气候条件相似，并且此方法以日为计算时段，能较好扣除大降水的产流量，具体公式如下：

田间有效灌溉水量可认为是到达田间的灌溉水量中用于蒸散发消耗的水量。一般来说，田间灌溉水量略大于田间有效灌溉水量，在高效节水灌溉的模式下，田间有效灌溉水量接近田间灌溉水量。

有效灌溉水量可认为是灌溉水量中用于蒸散发消耗的水量。一般来说，灌溉水量大于有效灌溉水量，灌溉水利用效率α越高，有效灌溉水量越接近灌溉水量。

3.2 区域降水量与出入境水量计算方法反距离权重法是气象要素插值最为常用的方法之一，该方法简单易行，插值效率高［17-19］。本研究的降水数据共选用常年雨量站60 个（见图1），平均站网密度248 km2/站。通过反距离权重法将点雨量转变为面雨量，计算区域降水总量。

区域流入、流出水量通过18个出入境控制站的流量监测数据统计得到。出入境控制站位置如图1所示。区域调入水量通过对调水工程调入水量统计得到。

3.3 陆地总蓄水变化量计算及分析水量平衡方程中的区域蓄水变化量是指区域地表、地下和土壤中蓄存水量的变化量，单纯依靠地面站点监测，难以直接获得；依靠水文模型连续模拟计算分析时段水蓄变量，一方面计算量庞大，另一方面又受模型精度影响，尤其是无资料地区水文模拟的不确定性干扰，难以实用。

有大量学者研究在大空间尺度上直接使用GRACE 卫星数据估算陆地总蓄水量，并通过水量平衡、实测数据和水文模型等方法进行验证，发现估算结果与验证数据具有较高的相关性［20-23］。本研究选用美国空间研究中心（CSR）提供的2007—2015年的CSR RL05 Mascon陆地水储量产品（本产品来自美国空间研究中心：http：//www.csr.utexas.edu），时间分辨率为月，空间分辨率为0.5°［24-25］。该数据以等效水柱高（地球水相关质量变化与水密度作商）表示陆地水储量变化，乘以单元格面积可转换为水量，在概念上与水量平衡方程中的区域水蓄变量一致。

由于目前重力卫星陆水储量产品的空间分辨率还较粗，用于较小空间范围时，需要做好数据的可用性分析。本文的研究区济南市总体面积不大，为检验重力卫星数据的可用性，一方面将陆水储量数据与济南市地下水观测数据进行对比分析，如图2所示，二者具有较好的一致性，可认为重力卫星在济南市范围内的数据可以反映蓄水量变化趋势；另一方面利用水量平衡方程来判别，如果存在一种蒸散发数据，二者能使水量平衡方程成立，则认为这组重力卫星和蒸散发数据在该区域可用。

图2 研究区2007—2015年地下水位及陆地总蓄水量变化

3.4 区域耗水量分解与灌溉水量估算区域总耗水量主要包括太阳能驱动的耗水量、工业能驱动耗水量和生活耗水量，可以由水量平衡方程推算得到，见式（1）；太阳能驱动的耗水量可以通过式（2）计算得到，同样能通过遥感数据反演得到，即遥感蒸散发。水量平衡得到的耗水量更符合实际情况，但是没有空间分布，可以利用遥感数据的空间分布，通过土地利用类型将总太阳能耗水量分解为农田区域和非农田区域耗水，得到区域农田耗水量。

现有遥感蒸散发产品比较多，由于其计算原理的不同，对不同区域和不同水文气象情况的适用性不一样［26-27］。本文收集了应用较多的MOD16和GLDAS两种蒸散发产品，根据3.1节所述原理，分析蒸散发产品与陆地水储量产品在水量平衡方程中的匹配性，依据水量平衡选取较为合理的蒸散发产品，与区域土地利用数据（图1）结合，计算农田蒸散发量。利用式（5）计算农田区域的有效降水（或称降水耗水），利用式（4）计算农田区域的田间有效灌溉水量（或称灌溉耗水），实际田间灌溉水量会大于田间有效灌溉水量。

4 结果分析

4.1 区域水量平衡分析与蒸散发产品选择研究区水量平衡各分量计算结果如表2所示。在研究选取的4年中，2012 与2015年接近多年平均降水条件（648 mm），2013年为丰水年，2014年为干旱年。

表2 研究区水量平衡计算结果表

利用重力卫星陆地水储量数据计算得到济南市蓄水量变化数据表明，2012年基本持平，2013年略减少，2014年明显减少，2015年明显增加，如图2所示。

根据水量平衡原理，表2中区域总耗水量（7）=（1）+（3）-（4）+（5）-（6）。采用2012—2015年济南市水资源公报统计发布的工业和生活耗水量作为工业能和生物能驱动的耗水量，从区域总耗水量中扣除该量，得到能与遥感蒸散发概念相对应的太阳能驱动的蒸散发量。

表3 研究区遥感蒸散发量比选

与两类遥感蒸散发产品GLDAS-NOAH、MOD16 进行对比可以发现：2012、2013 和2015年相对湿润的年份，MOD16 遥感蒸散发产品计算结果更接近水量平衡方程要求，相对误差小于13%；2014年是降水较少的干旱年份，GLDAS-NOAH 遥感蒸散发产品计算结果更接近水量平衡方程要求，相对误差为12%。这个对比结果与多篇文献对MOD16、GLDAS-NOAH 遥感蒸散发产品在我国适用性评价的结论基本一致，即正常来水年份，MOD16精度尚可［28-30］，干旱年份MOD16的蒸散发量明显偏小［26］；而GLDAS-NOAH的蒸散发产品在干旱年份表现较好［31］，正常年份普遍偏高［27］。本文选取2012、2013和2015年的MOD16产品和2014年的GLDAS-NOAH产品组合参与后续的灌溉用水量分析。

4.2 计算灌溉用水量按土地利用类型，提取出农田区域的蒸散发和降水量。利用插值后逐日网格降水量，计算农田区域的有效降水。从农田区域蒸散发总量中扣除有效降水量，得到农田区域的灌溉耗水量，参考公报的用水效率计算得到农田区域的灌溉水量，见表4。

由表4中的计算结果可见，2014年农田蒸散发量最大，但降水量和有效降水量均最小，农田区域必然存在大量灌溉水分补充，来支撑较高的蒸散发消耗。其余3年的农田蒸散发量较为接近，2015年降水量虽明显小于2012、2013年，但由于其降水的年内分配较为均匀，产生了较多的有效降水，推算2015年灌溉水量仍与2012、2013年相当。

表4 研究区农业用水计算结果

4.3 计算结果合理性分析将本文计算得到的灌溉耗水量与济南市水资源公报统计的灌溉耗水量进行对比，发现在2012年、2013年和2015年，来水条件接近多年平均或偏丰，本文计算的灌溉耗水量与公报中统计灌溉耗水量较为接近，本文的结果总体偏大。在较为干旱的2014年，本文的结果比公报统计的结果明显偏大。为分析2014年灌溉耗水量计算结果的合理性，从降水量、农作物长势、区域水蓄变量几方面进行年际差异的比较。

（1）降水量年际差异分析。通过站点插值计算得到研究区2012—2015年降水量空间分布图，如图3所示，区域降水量年际变异性较大，总体降水量和农田区域的降水量均是2013年>2012年>2015年>2014年，从图4中可以看出2012年和2013年降水集中于夏季，降水大部分产流消耗，2014和2015年降水较均匀，故2015年的有效降水大于2012年。

图3 研究区2012—2015年降水量空间分布

图4 研究区2012—2015年降水量时间过程分布

（2）作物长势年际差异分析。利用遥感归一化叶面积指数（NDVI）产品，分析研究区农作物长势的年际差异。图5为主要作物冬小麦每年第129天NDVI的空间分布，图6为NDVI与降水时间序列统计对比。由图5、图6 可见，整体上冬小麦的长势在时间和空间上均未表现出明显的年际差异，虽然2014年降水较少，但作物保持了正常的长势，这与4.2节中2014年农田区域出现较高的蒸散发也保持了一致。从支撑作物生长和蒸散发的水源上考虑，在2014年应该存在较多的灌溉补水。

图5 研究区（第129天）冬小麦NDVI空间分布

图6 研究区NDVI与降水时间序列统计结果

（3）区域水蓄变量年际差异分析。利用重力卫星的陆地水储量数据辅助分析区域蓄水量变化。由图2和表2第6列可见，2014年陆地水储量存在明显减少，且剧烈下降发生在春灌期间。结合降水量和作物长势的对比分析，可认为2014年消耗了较多的蓄水量为作物生长提供灌溉。

综上，干旱的2014年应存在较大的灌溉水量，本文分析得到的14.59亿m3灌溉耗水量，与降水量、蓄变量、作物长势等指标的年际相对性更为匹配，2014年灌溉水量远大于2012、2013、2015年应是更接近实际情况的结果。

5 结论与展望

本文以多源遥感、水文水资源地面观测和统计数据为基础，构建基于区域水量平衡和田间水量平衡联合的灌溉水量推算方法，梳理了遥感监测数据与水量平衡要素的对应关系，构建基于区域水量平衡方程对多源遥感数据进行比选、组合的方法，获取满足区域水量平衡约束的蒸散发空间分布，再通过田间水量平衡方程对农田总蒸散发进行分解，得到灌溉耗水或用水量的合理范围，使水量平衡方法用于区域灌溉水量分析的可操作性明显增强。

通过收集济南市2012—2015年相关数据，应用本文方法进行了灌溉水量分析。结果表明，济南市在2012—2015年降水有明显的年际差异，2014年干旱，区域蓄水量明显减少。通过水量平衡分析推算出2012—2015年济南市的灌溉水量分别为12.03亿、11.20亿、17.71亿、10.88亿m3。这4年济南市水资源公报统计的灌溉用水量分别为8.73亿、8.18亿、8.23亿、7.26亿m3，整体呈下降趋势，未体现出2014年干旱对农业灌溉的影响。总体上，本文推算的结果在2012、2013、2015年均表现出与公报数据的一致性，2014年二者存在差异，经分析认为，本文结果更为合理。

研究发现：引入多源遥感数据构建区域和田间两级水量平衡方程，推算区域灌溉水量，具有以下优点：（1）重力卫星陆地水储量变化监测和遥感蒸散发监测数据为区域水量平衡方程的构建提供了新的数据源，如果遥感数据及反演模型具备足够的分辨率和精度，有望解决水量平衡方程中蓄变量和农业耗水量获取的难题，提高定量推算灌溉水量的可操作性；（2）利用区域水量平衡方程约束遥感数据的误差，在一定程度上保证了利用田间水量平衡推算灌溉水量的合理性；（3）由于遥感数据本身具有的空间分布特征，本文方法在推算得灌溉水量总数的同时，还可输出灌溉用水的空间分布。

引入卫星遥感数据后，基于水量平衡的灌溉用水推算方法虽然得到一些改进，但也还存在一些有待完善的问题：（1）虽然有了区域蓄变量和农田耗水量的遥感监测数据，但灌溉水量的推算还需要以区域出入境水量、工业和生活用耗水等数据具有较高的统计精度为假设前提，依赖于水量平衡方程中其他要素准确获取，实际上某些区域的出入境水量，尤其是地下水出入境量还难以准确获取；（2）无论是基于重力卫星的陆地水储量数据，还是遥感蒸散发数据，都还存在着空间时空分辨率不高和精度不稳定等问题，有待于通过卫星遥感传感器和反演算法的改进来加以解决；（3）农田区域有效降水的计算多采用经验公式方法，不同区域的通用性差，也还有待进一步研究完善。