基于MODIS数据的黄河源区土壤 干湿状况时空格局变化*

刘 馨，宋小宁†，冷 佩，夏 龙

(1 中国科学院大学资源与环境学院， 北京 100049； 2 中国农业科学院农业资源与农业区划研究所， 北京 100081)

干旱是严重的自然灾害之一，持续的干旱灾害会直接影响农业生产和经济发展，长期干旱甚至会造成植被减少、河水断流、沙漠化现象[1]。有证据表明，自20世纪70年代以来，干旱等极端气候现象频次不断增加[2]。因此，加强干旱监测，对于农业管理、旱灾防治、水资源管理及保护都有着极强的现实意义[3]。

黄河源是黄河重要的水源补给区，近30年来由于人类活动及气候变化的影响，黄河源区气温上升，下游频繁断流，极大地影响着中下游地区的生态环境和生产生活。因此对黄河源区进行区域尺度上的干湿状况时空分析对于水资源管理和生态环境保护具有重要意义。

干旱作为一种复杂的现象难以直接观测，因此通常采用干旱指标进行评估。基于观测手段的不同，干旱监测评估主要可分为两类：基于站点的观测指标和基于遥感的观测指标[2]。相比传统基于站点的观测，遥感技术能够实时动态地获取地表参数的时空分布信息，从而表征地表干湿状况，具有覆盖范围广，多时相等优点，是目前干旱监测的重要手段之一[4]。

干旱过程与土壤干湿状况以及作物水分亏缺有着紧密联系。20世纪90年代开始，基于可见光-近红外波段反演的植被指数和基于热红外反演的地表温度形成的三角形或梯形的特征空间逐渐发展用于旱情遥感监测。由于相同光照条件下，植被温度在空间上没有变化，三角形内的温度变化只反映土壤表面的干燥程度。随着植被覆盖增加，温度变化范围逐渐缩小，因此温度/植被特征空间内的散点随着植被增加而倾斜形成三角形或梯形。三角形边界反映物理边界限制，即裸露土壤或完全的植被覆盖，以及完全干燥或完全湿润的土壤状况[5]。三角形概念最初由Price[6]提出，后由Carlson等[7-8]详细阐述。Moran等[9]提出用三角形空间内作物水分胁迫指数等值线评估水分亏缺。随后，学者们陆续利用温度/植被特征空间研究土壤水分反演[10-12]。Sandholt等[12]在水分亏缺指数WDI (water deficit index) 基础上发展了温度植被干旱指数TVDI (temperature vegetation dryness index)。该指数假设土壤水分是温度变化的主要来源，指数的计算无需地面的辅助数据，通过光学和热红外遥感数据的反演就可直接计算。该指数已被广泛用于土壤水分反演和干旱状况的监测[12-16]，其结果表明，TVDI表征的干旱程度及范围与实地监测气象数据和土壤水分数据吻合[17-18]。因此，本研究基于LST-NDVI特征空间通过计算TVDI来研究黄河源区土壤相对干湿状况，进而分析干旱时空变化，为有效对黄河源地区进行水资源管理和生态环境保护提供科学依据。

1 研究区与数据源

1.1 研究区

黄河源区(95°50′E～103°30′E，32°30′N～36°10′N)为河源至唐乃亥之间的汇水区域(图1)，流经青海、甘肃、和四川3省[19]。该区域集水面积为1.22×105km2，占黄河流域面积的15.3%，多年平均径流量为2.051×1011m3，占全流域的34.5%，是黄河流域重要的水源地和产流区，在黄河流域有着不可或缺的重要作用[18]。源区平均海拔4 000 m以上，为低山谷和湖盆组合地貌。植被类型以高寒草甸、高寒沼泽草甸和高寒草原为主，土壤类型包括高山草甸土、沼泽化草甸土、高山荒漠土、高山草原土等[20]。

图1 黄河源地区示意图Fig.1 Location of the source area of Yellow River(SAYR)

黄河源区气候属典型高原大陆性高寒气候，空气稀薄、辐射强烈且日照时间长。受西南季风影响比较明显，气温和降水从东南向西北呈递减趋势。源区年均气温在5 ℃左右，源区多年平均降水量为426 mm，年内变化较大且分布极度不均匀，主要集中在6—9月，约占全年降水量的75%～90%甚至以上[20]。

1.2 数据源及预处理

本研究使用遥感数据为2007—2016年7—9月(时相第193～257天)黄河源区的地表温度和NDVI的MODIS影像数据，分别为16天合成的1 km分辨率的MOD13A2的NDVI产品数据和8天合成的1 km分辨率的MOD11A2的LST产品数据。数据处理过程中使用MRT对上述数据进行几何校正、投影变换的批量处理，使用IDL将数据批量转换成实际地表温度和NDVI的数值，同时对NDVI数据水体进行掩膜处理，以消除水体对特征空间可能造成的影响。考虑到植被指数的值对植被生长状况的敏感性问题，只选取NDVI处于0.1～0.9之间的数据进行干湿边反演计算。此外，利用Arcgis裁剪出产品数据中所需的研究区域部分，剔除异常值。同时，为保持时相一致性，将8天合成的MOD11A2的LST产品数据合成为每16天的LST数据。

为了对构建的TVDI作为旱情指标进行可信度检验，本研究利用黄河源玛曲地区(玛曲站点连续原位土壤水分与土壤温度的测量——CEOP-AEGIS项目)站点提供的表层土壤水分观测数据与TVDI进行相关性分析。图1(b)为这20个站点的分布位置(具体信息见表1)，站点监测频率为每15 min一次。光学-热红外遥感所在波段穿透性弱，因此研究选择表层土壤水分进行TVDI相关性分析。众多研究选取的表层土壤水分深度不尽相同，大多数研究选择0～10 cm深度范围，而有的研究选取0～20 cm深度范围内的土壤水分数据，均与TVDI有较好相关性[18，21-22]。本研究参照黄河源区的土壤水分相关研究，同时考虑到实测数据的可获取性，选择5 cm深度的土壤水分数据。由于卫星过境的时间在11时左右，为保持实测土壤水分与MODIS数据时相一致，选取对应的每16天11时的土壤水分测量值计算平均土壤水分。

2 研究方法

2.1 温度植被干旱指数

LST和NDVI 特征空间可以表征地表植被和干湿状况，从而实现旱情监测，在国内外有了广泛的研究[8]。基于热红外遥感的地表温度和基于可见光-近红外波段的植被指数改善了单一使用表征地表干湿状况的不足，可以有效地判断作物长势和干旱程度，同时解决了植物在水分亏缺时表征现象的滞后性。2002年Sandholt等[12]提出TVDI(温度-植被干旱指数)估测地表土壤水分状况，由地表温度和NDVI计算而得(计算示意图如图2)，其形成的三角形或梯形的空间特征表征地表的干湿状况。计算公式如下

(1)

表1 实测站点信息Table 1 In-situ site information

式中：Ts为像元温度；Tsmin为相同NDVI值的最小地表温度，对应特征空间的湿边；a+bNDVI对应特征空间的干边，为相同NDVI值的最大地表温度，a和b为干边拟合线性关系的系数。干边对应的TVDI为1，表示极度干旱；湿边对应的TVDI为0，表示极度湿润。TVDI范围在0～1.0之间，地表干旱程度与TVDI数值呈正比[12,23]。

图2 TVDI计算示意图[13]Fig.2 TVDI triangle algorithm[13]

2.2 LST/NDVI特征空间构建

本研究使用批量重投影、拼接、裁剪等方法处理LST和NDVI产品数据，并利用处理好的产品数据构建出LST/NDVI特征空间。计算TVDI时，根据像元的NDVI值确定对应的温度变化，从而确定该像元对应的干边和湿边，再根据像元的地表温度值以及确定的干湿边计算出TVDI。干湿边的确定参照Garcia等[23]的方法，干边对应公式(1)中(a+bNDVI)，湿边与X轴平行，对应NDVI所在像元最低温度的均值。如前所述，研究区植被覆盖和土壤水分范围涵括所有极端情况且像元足够多时，LST/NDVI特征空间呈三角形，但实际情况中到达模拟干边的温度时，植被为减少水分流失会主动关闭气孔，因此土壤水分湿度不会因为高温而迅速减少。在NDVI高植区域的干边像元数值也不可能为0，特征空间内的散点多呈梯形分布特征。

图3为基于MODIS的地表温度和NDVI产品数据拟合的TVDI梯形特征空间，选取2007—2016年间的部分拟合结果为例。由图3可知，干湿边的拟合效果较好，R2均为0.9以上，最高的甚至达到0.99。算法假定当特征空间内植被覆盖度及土壤水分含量变化范围较大时，特征空间呈梯形。空间内干湿边为直线。LST/NDVI梯形特征空间对应的4个顶点在图3中分别为：1)水分饱和的裸土；2)干燥裸土；3)水分充足的全覆盖植被；4)水分亏缺的全覆盖植被。

2.3 干旱等级划分标准

本研究获取TVDI以后，根据相关文献提出的基于TVDI的干旱等级划分标准[16, 19, 24-26]，将干旱程度做了等级划分，便于后续对黄河源区域的干湿状况进行时空分析。具体划分指标如表2。

图3 LST/NDVI特征空间Fig.3 LST/NDVI space

等级类型TVDI1极湿润0～0.22湿润0.2～0.43正常0.4～0.64干旱(轻度)0.6～0.85重旱0.8～1.0

3 结果与分析

3.1 TVDI作为旱情指标的可信度评价

由TVDI的物理意义可知，TVDI与土壤水分具有负相关性；众多研究也表明TVDI与土壤水分存在负相关的线性关系[12-15, 27-28]。本研究利用TVDI作为干旱指标对黄河源区进行10年的干湿状况时空分析，因此在分析前利用站点所测的表层土壤水分(5 cm)与TVDI进行相关性分析,有利于检验TVDI作为干旱指标表征黄河源区域尺度土壤干湿状况变化的有效性。

由于云、水等因素影响导致部分LST和NDVI产品数据缺失，加上极端天气、设备条件等因素影响站点的实测数据，研究选取时间尺度上横向和纵向两组数据：2010年D193～D257的5个时相的土壤水分数据，以及2008—2010年D209的3年土壤水分数据与对应的TVDI进行相关性分析，增强了验证的可信度。由图4可知，TVDI与土壤水分呈显著的负相关线性关系，其中图4 (a)为TVDI与2008—2010年期间D209时相的土壤水分的线性关系，图4 (b)为TVDI与2010年D193～D257的5个时相的土壤水分的线性关系。图4 (a)中线性拟合的R为0.71，图4 (b)线性拟合的R为0.70，说明TVDI与土壤水分具有较好的相关关系。图4 (a)、4(b)两组线性关系略有差异，可能是由于时相选取不同，不同时相具有的不同天气状况(如云层、降雨等)会对遥感反演结果造成一定的影响；且气候条件和设备条件也会影响实地站点的数据精度，这些因素都使得TVDI作为干旱指标可能会与实际地表干湿状况有一定的偏差。但是总体来看，图4两组线性关系良好，说明TVDI很够有效表征土壤表层干湿状况。

3.2 黄河源区土壤干湿状况时空格局分析

3.2.1 土壤湿度时间格局分析

图5(e)为10年来TVDI自D193至D257五个时相的变化规律图。逐时相来看，除2008年D257外，其余时相的趋势大致一致。10年间的TVDI均在0.4～0.6，属干旱分级中的正常状况。各时相旱情大体上波动趋势一致，说明在旱情严重的年份，不同时相旱情都有一定程度的加重，旱情缓解的年份，不同时相的旱情也有一定程度的缓解。各月旱情在10年期间的波动变化各有区别，也是由于不同时相当地的天气变化，如降雨、云层的影响，造成对旱情变化的短期影响。

图4 TVDI与土壤水分线性关系Fig.4 Linear relationship between TVDI and soil moisture

图5 2007—2016年TVDI时间格局分析Fig.5 TVDI analysis of temporal patterns from 2007 to 2016

此外，根据中国科学院资源环境科学数据中心发布的2010年土地利用分类图，本研究选取高、中、低覆盖度草地区域，逐时相分析TVDI 在不同植被状况下对土壤干湿状况的指示。为控制降雨量、地形、地表覆盖类型等变量因素，选取东部地区相离较近的实测站点所在区域。东部地区地势平坦且植被覆盖度较广，选取的3个实测站点(CST_01, NST_01, NST_07)地形均为河谷区域，地表覆盖类型为草地。TVDI变化规律如图5(a)～5(d)所示。图5(a)、(b)、(c)分别为中、低、高覆盖度草地的TVDI趋势图，图5 (d)为10年内D241时相上TVDI在3种不同覆盖度草地的变化趋势图。东北部地区整体来看处于干旱-重旱。2007年—2016年间，在中低草地覆盖区，7月上旬旱情最为严重，中覆盖区有5年时间TVDI数值都在0.8及以上，为重旱程度，10年期间TVDI均在0.6以上，为干旱程度；低覆盖区有7年时间TVDI达0.8以上，为重旱程度，10年期间TVDI均在0.6以上，为干旱程度。其次是9月上旬，中覆盖区有4年时间TVDI数值在0.8以上，为重旱程度，10年期间TVDI均在0.6以上，为干旱程度；低覆盖区有6年时间达0.8及以上，为重旱程度；10年期间TVDI均在0.6以上，为干旱程度。中低植被覆盖区10年来旱情均较为严重，TVDI数值大致均在0.6～1间浮动，其中低植被覆盖区所有时相的TVDI在10年内数值均在0.6及以上，旱情持续。各时相旱情在中、低覆盖度区域TVDI年际趋势相似，各月旱情变化基本具有一致性，而浮动差异可能是由于短期天气变化造成的一定影响。图5 (d)可以看出，NST_01与CST_01站点所在中、低植被覆盖区域与NST_07所在高植被覆盖区域10年间TVDI趋势有所差异。中、低植被覆盖区域在10年间TVDI变化大体相似，而高植被覆盖区域TVDI波动较大，在2008、2012、2016年与中、低植被覆盖区域的TVDI变化呈现出相反的趋势。高植被覆盖区域7月上旬有5年旱情严重，TVDI达0.8以上，其余时相旱情较轻，均仅有2～3年TVDI达0.8以上。尽管严重干旱的情况较中低覆盖区域缓解，但是干旱依然在10年间持续发生。D193和D257在10年间TVDI均在0.6～1之间小幅波动，D209仅有1年TVDI略低于0.6；D225和D241分别有3年和2年的TVDI数值低于0.6，但也均大于0.5。干旱与植被覆盖度互为因果，植被具有水分涵养的作用，植被覆盖度低的地方易水土流失，引发旱情，而干旱持续发生时亦会造成植被覆盖减少的情况[1]。

3.2.2 土壤湿度空间格局分析

本研究采用空间统计的分析方法，计算各等级旱情(参照表2)区域占黄河源区总面积的比例，从而分析旱情的时空分布特征。图6(a)为2007—2016年间D193～D257时相的旱情占比频谱图，横坐标为2007—2016年的D193～D257顺序各时相，纵坐标为各旱情面积占比；图6(b) 为2007—2016年间D193～D257时相的旱情占比直方图，横坐标为2007—2016年的D193～D257顺序各时相，纵坐标表示各旱情面积占比及其总和。图6(a)、6(b)结合分析可得，10年来7—9月正常干湿状况的区域所占面积最多，最高占比达48%；干旱区域和湿润区域面积交替占比第二，在10%～40%之间浮动；而极端情况即极湿润和重旱区域占比均较少，在0～17%之间浮动。各旱情面积占比波动规律较为明显，而干旱与湿润等级的区域面积波动有着明显的相反趋势。

本文中以2016年7月上旬—9月上旬(D193～D257时相)的TVDI图像为例进行说明。旱情的空间变化特征显著：每年7月上旬西部地区旱情较轻，东北部地区和东南部受旱区域多为严重干旱。7月下旬，干旱区域逐渐向西部迁移，东南部旱情缓解，东北部仍然遭受严重干旱。8月上旬东北部旱情缓解，东南、西北两部旱情有加重趋势。8月下旬开始，西北部旱情缓解，且逐渐湿润，旱情逐渐向东南部转移。9月上旬西北部旱情继续缓解，基本处于湿润-正常的地表状态，而东北、东南区域保持较严重的干旱状态。中部区域在7—9月一直处于湿润-正常的稳定状态。分布趋势基本吻合中西部水分整体水分较充足，干旱比较严重的区域主要在东北部和东南部的分布规律。

4 结论与讨论

本研究以黄河源为研究区，利用2007—2016年D193～D257的MODIS的1 km分辨率的地表温度(LST)和植被指数(NDVI)的产品数据反演10年温度-植被干旱指数(TVDI)。分析结果表明：

1)本研究提取NDVI对应的最大、最小地表温度构建LST/NDVI特征空间，10年间D193～D257时相的LST/NDVI散点图均构成梯形特征空间，符合前人的特征空间研究成果。

图6 2016年7—9月(D193～D257)黄河源区旱情变化图Fig.6 Variation of drought situation from July to September (D193-D257) in 2016

2)通过研究区站点监测的实地土壤水分与反演的TVDI建立线性关系进行验证，结果表明实地土壤水分与TVDI具有良好的线性负相关关系。2008—2010年实测土壤水分与TVDI线性关系的相关系数达到0.7，说明TVDI能作为干旱指标有效地指示黄河源区的干湿状况。

3)10年间D193～D257的TVDI均值属干旱分级中的正常状况；各时相TVDI大体上波动趋势一致，说明在旱情严重的年份，不同时相上旱情都有一定程度的加重，旱情缓解的年份，不同时相的旱情也有一定程度的缓解。东北部地区整体来看处于干旱-重旱，年年有旱情发生，7月上旬(D193)的干旱情况均最为严重。中、低植被覆盖区域的干旱趋势较为相似，10年来旱情均较为严重，其中低植被覆盖区所有时相在十年间旱情持续。高植被覆盖区域严重干旱的情况较中低覆盖区域缓解，但是干旱依然在10年间持续发生。

4)旱情的空间变化特征显著，7月上旬至9月上旬期间(D193～D257)，西部尤其是中部区域旱情较轻，严重旱情多集中于东北部和东南部区域，分布趋势基本吻合中西部整体上土壤水分较充足，干旱比较严重的区域主要在东北部和东南部的分布规律。

根据中国科学院资源环境科学数据中心(http:∥www.resdc.cn/data.aspx?DATAID=125)提供的中国生态地理分区图可知，黄河源区并非处于干旱区，而TVDI的旱情分级也只是提供黄河源区的相对干湿状况信息，并非绝对的旱情。本研究以TVDI为指标，研究黄河源区土壤相对干湿状况，为有效对黄河源地区进行水资源管理和生态环境保护提供科学依据。TVDI结合基于热红外遥感的地表温度和基于可见光-近红外波段的植被指数两种指示因子，改善了单一使用表征地表干湿状况的不足，同时解决了植物在水分亏缺时表征现象的滞后性。由于本研究使用NDVI、LST以及监测站点土壤水分观测值的时相不同，最终将所有数据统一对应NDVI时相，3种数据在时间对应上存在偏差，从而让反演结果产生误差。实地站点仅有2008—2010年3年的实地土壤水分数据，且由于极端天气，站点设备等原因，许多实测值不具有对比验证意义，可能也对遥感反演的数据验证存在一定的影响。后续将对TVDI反演结果进行进一步的研究，如参照Hoffmann等[29]提出的方法，考虑改进干湿边的拟合，为提高遥感影像反映地区干湿状况进行更深一步的工作。