1. 引言
地表实际蒸散发是水文循环和地表能量平衡的重要环节，也是植被生存必须的生理过程和地表物质循环的重要载体 [1] 。目前，已发展了从传统方法，模拟方法到遥感方法等很多种方法用于估算一个区域的蒸散发量。
传统方法包括水量平衡法、器测法、经验公式法、模式计算法等。传统方法的劣势在于：蒸散发的时空分布特征与区域气候条件、土壤水分状况、植被覆盖等因素密切相关，传统的方法都是以点的观测资料为基础来开展蒸散发的研究 [2] ，因此，很难直接运用于蒸散发研究所需要的区域尺度。遥感方法以其能够获知大范围的地表特征信息的优势为较为准确估算地表蒸散发提供了可能。基于遥感信息的蒸散估算方法在国内外得到迅速发展，主要有统计经验法、特征空间法和能量平衡法，代表模型有基于Penman-Monteith公式的VITT模型、基于Priestly-Taylor公式的三角形模型以及基于能量平衡的SEBAL模型、SEBS模型和TSEB模型等。但遥感方法也存在一些问题，例如，地表参数反演困难，模型有局限性及误差，在区域尺度上检验困难 [3] 。
近年来，国内从遥感蒸散发模型的发展，到地表蒸散发观测试验，以及典型生态系统类型的蒸散发模拟都取得了很多的成果 [4] [5] [6] [7] 。然而，由于遥感蒸散发模拟的复杂性，区域蒸散发的模拟主要应用在一些特定的区域或流域，针对更大空间尺度的研究不多。目前，SEBS模型、Priestley-Taylor模型、NOAH陆面模式和BEPS模型已应用于整个中国区域的陆面蒸散发模拟研究。
吴桂平、位贺杰和范建忠等人已率先运用MOD16产品分别运用于鄱阳湖流域、渭河流域及陕西省的地表实际蒸散发的时空分布特征研究，并且都取得了良好的、具有实际意义的分析结果 [8] [9] 。这些研究的过程及结果都对本课题的研究提供了有力的理论基础和参考依据 [10] [11] 。
2011年美国NASA研究团队在MODIS遥感数据蒸散反演算法上取得了重要研究成果，并通过NASA地球观测系统发布了全球MODIS陆地蒸散产品数据MOD16 [12] ，该产品不仅提供了地表蒸散发量的参数，还具有高时间分辨率以及免费获取的特点，因此利用MOD16产品来反映嫩江流域地表蒸散量的时空分布特征及其与土地利用类型的关系具有一定的优势。
本文将主要使用MOD16A2和MOD16A3产品，作为本课题研究的基础。NASA EarthData提供的MOD16产品有MOD16A2、MOD16A3两种，均提供包括蒸散发(ET)、潜热通量(LE)、潜在蒸散发(PET)、潜在潜热通量(PLE)等要素，MOD16A2产品的时间分辨率为8天，空间分辨率为500 m；MOD16A3产品的时间分辨率为1年，空间分别率为500 m。本文首先探讨运用该产品研究嫩江流域地表蒸散发量的可行性。然后，在结合嫩江流域土地利用数据，分析嫩江流域在不同时间尺度下地表蒸散发的数量特征和空间格局特。最后，再更进一步对不同土地利用类型下蒸散发量进行统计分析，从而可以为嫩江流域水文分析计算、水资源管理及高效利用提供参考。
2. 研究区概况及数据来源
2.1. 研究区概况
嫩江流域位于中国东北部，西侧以大兴安岭分水岭为界，北侧以伊勒呼里山为界，南侧以霍林河南部的分水岭为界，东侧大部分以嫩江为界，坐标为119˚12′~127˚54′，北纬44˚02′~51˚42′，土地总面积为28.30万平方米，地跨黑龙江、内蒙古、吉林3省区。该流域属于寒温带半湿润大陆性气候，春季多风少雨干旱，夏季短促高温多雨，秋季降温急剧且常有霜冻，冬季漫长严寒干燥。多年平均气温2℃~4℃，冬季嫩江冰封期达150天左右，冰厚约为1米。研究区范围：以大赉站控制范围作为本论文的研究区范围。
嫩江流域多年平均径流深为76.5 mm，多年平均降雨447 mm，天然年径流总量227.3亿m3。流域的径流主要来自于降水，春季会有一定量由于冰雪融化带来的融水。径流量的年内和年际变化均比较大，每年的1、2、3、11、12月份是枯水、枯冰期，降水量会明显减少，河川径流主要靠地下水补给；4、5、6、10月份是一年的平水期；7、8、9月份是丰水期，常发生大、暴雨，且易形成洪水。
嫩江流域北部、西部和南部三面地势较高，东部地势较低，形成著名的松嫩平原。土地利用类型多样，有林地、草地、耕地、水域、湿地等类型。林地、草地主要分布在北部、西部高山地区，东部和南部零星分布；耕地和水域主要分布在东部平原区，湿地主要有位于流域东南部的扎龙、莫莫格湿地群，位于南部的向海湿地，如图1。嫩江是我国重要的农牧业，工业，能源生产基地。近些年来，气候变化与人类对水土资源不合理开发利用，引起下游湿地入流量减少和蒸散量增加等。
2.2. 数据来源
本研究使用的数据具体来源见表1。MOD16-ET产品的算法是Mu等 [12] 在Penman-Monteith公式的基础上改进的，算法将地表分为裸土和冠层两种组分，又将冠层分为湿润和干燥两种类型，分别根据不同下垫面的生物物理过程计算各自的空气动力学阻抗和冠层阻抗，然后结合MODIS数据产品估算 [12] 。模型主要输入的变量及其来源：叶面积指数(LAI)和植被光合有效辐射分量(FPAR)从MODIS遥感数据中的MOD15产品获取；陆地覆盖类型(LAND COVER)来自MODIS遥感数据MOD12Q1产品；相关气象数据来自GMAO气象网站。根据MODIS数据轨道号的排列规律及嫩江流域所在的地理位置，选择的卫星轨道号为h25v03、h25v04、h26v03和h26v04。下载到的数据从2000年至2014年共15年的数据。根据嫩江流域的土地利用方式和分析精度要求，将研究区域的土地利用类型归并为6个一级地类，包括林地、农田、草地、城镇、裸地和水域。
Figure 1. Schematic diagram of land use types in the Nenjiang River Basin
图1. 嫩江流域土地利用类型示意图
Table 1. Data used in the study
表1. 本研究所用的数据资料
3. 研究方法
3.1. 数据预处理
对于MOD16数据的处理方法的具体如下：
(1) 原始的MOD16产品采用的是HDF数据格式，投影类型为正弦曲线投影(sinusoidal)。因此，需要用到NASA提供的MRT软件进行格式转换，将MOD16-ET产品的HDF文件转化为WGS-1984经纬度坐标系统下的GeoTiff格式文件，在转换格式的同时进行重采样等操作。
(2) 根据嫩江流域的地理位置，将不同卫星轨道的数据(h25v03, h25v04, h26v03, h26v04)用ArcGIS软件进行轨道镶嵌，制成一个栅格数据集，然后利用ArcGIS软件对经过投影转换、轨道镶嵌之后的MOD16数据基于嫩江流域的边界进行掩膜截取，从而可以得到研究区域的每一种时间尺度下地表蒸散发量的空间分布图。通过统计可以计算得到在不同时间尺度下，不同空间尺度下蒸散发量的平均值、最大值、最小值等数据。
(3) 为了能够利用MOD16产品进一步分析研究嫩江流域地表蒸散发量的变化特点和时空分布特点，利用ArcGIS软件和配准后的嫩江流域土地利用类型图，掩膜截取五种不同土地利用类型的区域，再通过上述步骤计算出年尺度下各种土地利用类型的蒸散发量的平均值，最后分析不同土地利用类型蒸散发量的变化特征。
3.2. 数据精度检验
MOD16产品的精度是本文研究的基础，为了验证MOD16产品数据可以用于嫩江流域的地表蒸散发量的研究，首先需要对流域MOD16产品地表蒸散发量的有效值进行检验，检验的依据是流域水量平衡方程。
(1)
公式(1)中，P为降雨量(mm)、E为蒸散量(mm)、R为径流量(mm)。在年尺度上可以用水量平衡方程式对流域内的MOD16地表产品的适用性进行评估和验证。
3.3. 技术路线图
本文技术路线图见图2，共分为3个步骤，1) 数据处理，包括针对源数据的格式转换以及淹没截取等；2) 精度检验，利用水量平衡原理对遥感产品的结果进行评价；3) 分析嫩江流域蒸散发时空分布特征。
Figure 2. Technology roadmap of this article
图2. 本文技术路线图
4. 结果分析
4.1. MOD16产品精度检验
利用嫩江流域大赉水文站2005年~2012年的降雨径流数据，估算出年均降雨量与年均径流深，从表2结果来看，总体上呈现出了较好的一致性。
Table 2. MOD16 product accuracy inspection form in the Nenjiang River Basin from 2005 to 2012
表2. 2005年~2012年嫩江流域MOD16产品精度检验表
图3展示的是嫩江流域多年蒸散发量的平均值空间分布，嫩江流域多年平均蒸散发量呈现由北向南逐渐降低的趋势，这也与图1显示的嫩江流域土地利用类型分布规律相吻合，所以MOD16地表蒸散产品在嫩江流域研究区范围内的精度总体上是符合要求的，所以我们可以继续利用MOD16产品数据对流域蒸散发量的时空分布规律进行后续的分析研究。
Figure 3. Spatial distribution of multi-year ET average values in the Nenjiang River Basin
图3. 嫩江流域多年ET平均值空间分布示意图
4.2. 蒸散发量年平均值变化特征
对MOD16/Annual ET产品数据经过上文介绍过的方法处理，得到嫩江流域地表2000年至2014年逐年蒸散发量的数值。下图显示了嫩江流域地表蒸散量的逐年变化过程。
Figure 4. Change process of surface evapotranspiration in the Nenjiang River Basin from 2000 to 2014
图4. 嫩江流域2000~2014年地表蒸散量变化过程
由图4可知，年均蒸散发量的波动范围在359.5~432.6 mm/a之间，而且各年间的平均蒸散发量平均值呈现出一定的波动。可以计算得出嫩江流域多年平均蒸散发量为396 mm/a (如图中虚线所示)。超出多年蒸散发量平均值的年份出现在2002年、2004年、2005年、2008年、2012年、2013年、2014年，其中，2004年接近多年平均蒸散发量，其余年份均低于多年平均蒸散发量。其中最大值出现在2014年，蒸散发量为432.6 mm/a，最小值出现在2009年，蒸散发量为359.5 mm/a。
4.3. 蒸散发量多年月平均值时空分布特征
根据MOD16/Monthly产品数据统计分析2000~2014年嫩江流域月平均蒸散发量的结果，如图5所示，总体来看，嫩江流域地表蒸散发量从12月到2月的变化较为平缓，在3月份会出现一个小的峰值，最大值出现在7月份或者8月份。蒸散发量主要集中在6月、7月、8月，这三个月的蒸散发量能占到全年蒸发量的50%以上。2000~2014年各年6、7、8月的蒸散发量占到一年中蒸散发量的比例，如表3所示。
Table 3. Proportion of evapotranspiration in June, July and August from 2000 to 2014
表3. 2000年~2014年6、7、8月蒸散发量所占比例
Figure 5. Distribution of surface evapotranspiration in different months in the Nenjiang River Basin from 2000 to 2014
图5. 嫩江流域2000年~2014年不同月份地表蒸散发量分布
Figure 6. Annual distribution of evapotranspiration in the Nenjiang River Basin from 2000 to 2014
图6. 嫩江流域2000~2014年蒸散发量年内分布
15年中，6、7、8三个月的总蒸散发量所占的比例变化范围在54%到61%之间，表明嫩江流域的蒸散发量受到流域的温度和降雨的影响很大。6月之后，嫩江流域的温度显著提升，大气循环加快，降雨也达到一年中的最大值，可以推断降雨是嫩江流域蒸散发最重要的来源。2000~2014年，各月平均蒸散发平均值的年内分布如图6所示。
分析出现这些现象的原因是：12月、1月、2月嫩江流域气温低，日照时间短，降水主要以降雪为主，地表水主要以固态的形式储存，不利于蒸发。所以呈现出一种平缓的趋势。3月份，气温开始回升，大量以固态储存的水开始大量转化为液态形成地表径流，液态水相对于固态水更利于蒸发，同时植被开始发育，加大了植物散发，因此3月份蒸散发量会出现一个小峰值。到了6月，气温继续上升，降雨量也大大增加，供水量大，日照充足，再加上风速大，因此地表蒸散发将大大增加，到了7、8月达到峰值。9月之后，气温开始缓慢降低，降雨量也开始下降，各种因素又向着不利于蒸发的方向发展，因此蒸散发量降到一年中的最低值。
4.4. 不同土地利用类型的蒸散发量分布特征
根据MOD16/Annual产品数据,结合嫩江流域土地利用类型，可以统计分析在2000年~2014年间，5种土地利用类型的年均蒸散发量并绘制年际变化曲线，如图7所示。可以看出5种土地利用类型的年际波动变化在趋势大体上是一致的。但是在数值上有很大的差异。
Figure 7. Average annual evapotranspiration process curve for different land use types
图7. 不同土地利用类型的年均蒸散发量过程曲线
Figure 8. Average annual evapotranspiration for different land use types
图8. 不同土地利用类型的年均蒸散发量
在一个区域环境中，由于不同土地利用类型本身的生理生态特性，会影响平均蒸散发量。为了能够分析嫩江流域上不同土地利用类型的蒸散发特点，运用ArcGIS软件的掩膜截取和空间统计功能，按照不同的土地类型(农田、林地、裸地、城镇、草地)进行裁剪，然后统计计算嫩江流域内不同土地利用类型的年均蒸散量，如图8所示，由图中可以看出不同土地利用类型跟蒸散量密切相关，蒸散量的大小按照“林地 > 草地 > 农田 > 城镇 > 裸地”的顺序递减。这也印证了植被覆盖好的地方，蒸散发量高。其中林地的年均蒸散发量最高，达到480.6 mm/a。城镇和裸地的年均蒸散发量最小，城镇的年均值是341.1 mm/a，裸地的年均值是340.8 mm/a。在鄱阳湖流域的蒸散发研究中，各土地利用类型的年均蒸散发量的顺序为“林地 > 草地 > 农田 > 裸地 > 城镇”。
从图8中可以看出，在11、12、1、2月份各种土地类型的蒸散发量区别不大，在某些月份林地的蒸散发量相比其他土地类型的蒸散发量要小，这与林地年均蒸散发量最大正好相反。原因可能是冬天的时候嫩江流域林地都被冰雪所覆盖，也很少受到人类活动影响，十分不利于蒸散发。不过到了6、7、8月份，林地的蒸散发量要远远大于其他土地利用类型的蒸散发量。
结合图3，可以发现每种土地利用类型的年内分布在4至9月份的差别很大，分别计算每种土地利用类型的变异系数，变异系数可以反映一组数据的离散程度。由表4，可以看出林地的变异系数最大，裸地的变异系数最小。这表明林地的月均值离散程度最大，这是因为林地的蒸散发量在一年中对温度和降雨量的影响最敏感。
Figure 9. Intra-year distribution of evapotranspiration by different land use types from 2008 to 2014
图9. 2008年~2014年不同土地利用类型的蒸散发量的年内分布
从图9中可以看出，在11、12、1、2月份各种土地类型的蒸散发量区别不大，在某些月份林地的蒸散发量相比其他土地类型的蒸散发量要小，这与林地年均蒸散发量最大正好相反。原因可能是冬天的时候嫩江流域林地都被冰雪所覆盖，也很少受到人类活动影响，十分不利于蒸散发。不过到了6、7、8月份，林地的蒸散发量要远远大于其他土地利用类型的蒸散发量。
Table 4. Coefficient of variation of annual changes in different land use types
表4. 不同土地利用类型年内变化的变异系数
结合图3，可以发现每种土地利用类型的年内分布在4至9月份的差别很大，分别计算每种土地利用类型的变异系数，变异系数可以反映一组数据的离散程度。由表4可知，林地的变异系数最大，裸地的变异系数最小。这表明林地的月均值离散程度最大，这是因为林地的蒸散发量在一年中对温度和降雨量的影响最敏感。
5. 结论
(1) 通过嫩江流域水文站的降雨、径流资料统计分析，利用水量平衡方程，实测降水量与实测径流量估算的流域的蒸散发量与MOD16产品测算出的流域蒸散发量有较好的一致性，所以MOD16产品的精度和有效性能够满足用于研究嫩江流域地表蒸散发量时空变化分析的需求。
(2) 通过统计分析年尺度的MOD16产品，从时间序列上来看，在2000年~2014年中，最大值出现在2014年，最小值出现在2009年，多年平均蒸散发量为396 mm/a。从空间分布上来看，流域地表蒸散量呈现出由北向南逐渐减小的趋势。
(3) 嫩江流域蒸散发量年内分布由于受到气温、降雨量、风速、太阳辐射等影响，各月间会出现较大的波动。大致呈现双峰型的分布趋势。在三月会出现一个小峰值，在七月或八月出现一年中的极大值。而一年中超过一半的蒸发量出现在6、7、8月。
(4) 通过2000年~2014年的年尺度的MOD16产品结合嫩江流域的土地利用类型，可以分析得到林地的年均蒸散量最大，达到480.6 mm/a。之后依次是草地、农田、城镇、裸地。
(5) 通过2008年~2014年的月尺度的MOD16产品结合嫩江流域的土地利用类型，可以分析得到各种土地利用类型的蒸散发量的年内分布于整个流域的蒸散发量的年内分布大致吻合，但是林地的月均蒸散发量变异系数最大，裸地的最小。这种现象是各个土地利用类型受到自身特性、气象条件、农作物特性和人类活动等综合作用的结果。
基金项目
中国长江三峡集团有限公司科研项目资助(0799251, 0704189)。