基于NDVI分区的全国草地生态系统覆盖区土壤水分遥感反演

2022-05-04郭兴国

水土保持研究 2022年3期

杨娜, 郭兴国

(1.济南市勘察测绘研究院, 济南 250100; 2.三峡大学水利与环境学院, 湖北宜昌 443002)

草地生态系统作为陆地生态系统中的敏感成员，在发展畜牧业、维持生物多样性、维持生态系统平衡等方面发挥着重要作用，但其发育与演替受气候、人类活动等影响较为严重[1-2]。据气象报告显示过去几十年间气候经历了明显变化，这势必会对我国(尤其在干旱半干旱地区)草地生长产生显著影响，进而影响全球的物质循环和能量流动，对生态系统的稳定性产生重要影响[3-4]。土壤水分是表示土壤干旱程度的重要物理量，是地表和大气—水热过程的重要纽带，是监测土地退化的一个重要指标[5-8]，其含量虽仅占全球水循环分量的0.005%左右，却影响着地—气界面间的物质循环和能量交换[9]，特别对干旱事件和高温热浪的形成和发展及其严重程度有重大驱动作用[10]。传统的土壤水分测量方法(如烘干称重法)是利用探针或重量测量法测定不同深度的土壤水分[11-13]，虽可得到较为精确的土壤水分数据，但数据代表性差、获取过程劳民伤财，获得的基于站点的数据已不能准确刻画大尺度范围内土壤水分的时空变化[13-16]。20世纪70年代中期，遥感技术因具有时效快、范围广、长期动态监测及高时空分辨率等优点，被广泛应用在土壤水分的反演工作中[17-21]。利用遥感监测土壤水分的方法有可见光—近红外法(如反射率法、植被指数法)、热红外法(如表观热惯量法)、微波法(被动微波遥感法和主动微波遥感法)和多源遥感法等[21-30]。其中，可见光—近红外法虽具有高空间和时间分辨率、大范围长时间动态监测的优势，但因其穿透能力较弱，且易受到天气、植被覆盖类型等的影响，使反演精度较低[5]。热红外法仅适用于裸土和植被覆盖度低的区域，在植被覆盖度高的区域因植被覆盖了土壤的信息，严重影响土壤水分反演精度[10,31-32]。土壤水分遥感反演方法很多，但对不同下垫面(可谓裸土、部分植被覆盖和全植被覆盖区域)需使用不同的反演模型[33]。虽目前已有很多学者对土壤水分遥感反演做了大量工作，但因研究数据、研究区、植被类型和起始年份不同而存在差异[30-33]。而土壤水分研究作为草地生态系统动态监测的一项长期性、经常性的工作，被列为生态保护与修复工程的重要任务之一，引入土壤水分探究在生态保护与修复工程措施下草地覆盖动态及其对气候变化的响应情况，对掌握和评价生态保护与修复等工程的所取得的效果有重大意义。

鉴于此，本文以全国草地生态系统为研究区，鉴于表观热惯量和温度植被干旱指数模型优缺点，基于2012年、2016年已有的实测数据与反演结果在精度验证的基础上确定NDVI阈值，建立土壤水分遥感反演模型。并使用该模型反演2016年全国草地生态系统0—10 cm土壤水分，分析土壤水分时空分布格局，以期为我国生态环境改善、畜牧业管理措施及未来生态环境建设提供科学依据。

1 材料与方法

1.1 研究区概况

草地类型数据源于中国生态系统评估与生态安全数据库(http:∥www.ecosystem.csdb.cn/ecosys/)提供的草地生态系统数据。由中国科学院生态环境研究中心依据《中国生态系统》的分类方法，采用6级分类单位将草地生态系统分为了草原生态系统(分布于我国温带地区，如东北，华北地区，青藏高原，黄土高原等)、草丛生态系统(主要分布在在热带和亚热带，向北可以分布到华北地区。)、草甸生态系统(普遍分布于我国北方各省区及青藏高原)、高寒生态系统(广泛分布在青藏高原以及外缘山地，北起祁连山，南至横断山、天山北坡、玛纳斯河上游的大中大坂)四大类。中国草地生态系统总面积达3.9亿hm2，占世界草地总面积的13%，占全国国土面积的41%(图1)。

1.2 数据来源

DEM数据：DEM数据源于地理空间数据云(http:∥www.gscloud.cn/sources/accessdata/310?pid=302)的SRTM产品数据，分辨率为90 m。

土壤水分实测数据：2012年土壤水分实测数据由地球大数据科学工程数据共享服务系统(http:∥data.casearth.cn/)提供，多为人工观测，观测频率为5～15 d一次，观测深度有0—10，0—20，0—40，0—80，0—120 cm。2016年土壤水分实测数据由寒区旱区科学数据中心(http:∥westdc.westgis.ac.cn/data/)提供的大沙龙站和峨堡站及国家青藏高原科学数据中心(http:∥www.tped-atabase.cn/portal/MetaDataInfo.jsp?MetaDataId=249 456)提供包含寒冷半干旱的那曲站、寒冷潮湿的玛曲站、寒冷干旱的阿里站和狮泉站，每隔15 min观测一次距地面5 cm，10 cm，30 cm，50 cm，80 cm处的土壤水分。

土壤湿度产品数据：由中国气象数据网(http:∥data.cma.cn/ data/ cdcdetail/ dataCode)提供的“CLDASV2.0”产品。土壤相对湿度产品数据有0—10 cm，0—20 cm，0—50 cm观测深度数据。利用中国区域业务的质量控制后的土壤水分自动站观测资料对其评估，结果表明：CLDAS-V2.0土壤相对湿度产品与地面实际观测相关系数0.8左右，均方根误差小于0.1，偏差为0.1。

遥感数据：遥感数据包括MOD09GA(时间分辨率为1 d，空间分辨率为500 m的地表反射率数据，提供1～7波段的反射率、质量评价等数据)，MOD11A1和MYD11A1(时间分辨率为 1 d、空间分辨率为 1 km的陆地表面温度，包含昼夜地表面温度和对应的质量控制数据)，所有数据均由Google Earth Engine平台处理并下载。投影为Albers正轴等面积双标准纬线圆锥投影，并利用质量控制文件对MOD09GA和MOD11A1数据进行质量控制，最大限度减少云等其他噪声影响。使用MOD09GA数据提供的波段分别计算得到NDVI和EVI。NDVI，EVI计算公式见(1)，(2)：

NDVI=(NIR-R)/(NIR+R)

(1)

EVI=2.5(NIR-R)/(NIR+6R-7.5B+1)

(2)

式中:NIR为MOD09GA数据对应的近红外波段;R为红光波段;B为蓝光波段。

2 研究方法

2.1 遥感数据处理方法

2.1.1 地表温度的地形订正不同海拔高度下，地表温度受气温和大气湍流的影响。当研究区有明显地形起伏时，随高程增加，气温降低，地温也降低，不经高程订正，高程高区地表温度低，会造成对该区域土壤湿度的高估。对地表温度数据(LST)使用DEM数据对其进行订正，订正见公式(1)[7,16,26]。

Td=Ts+aH

(3)

式中:Td为经DEM订正后的地表温度值;Ts为原始的地表温度值;H数字高程(DEM)值;a表示地表温度随海拔高度增加而降低的程度,结合相关参考文献,本文的a取值为-0.6/100℃。

2.1.2 遥感数据缺失值填充因NASA提供的 MODIS LST标准产品采用的是劈窗算法和昼夜地表温度反演算法，在晴空条件下经度可达1 K，但热红外波段易受云和其他干扰因素影响，使获得的数据存在大量数据缺失，需对数据做填充。本文是采用上午星(Terra，从北向南于地方时10∶30左右通过赤道)获得的MOD11A1数据和下午星(Aqua，从南向北于地方时13∶30左右通过赤道)获得的MYD11 A1数据，在质量控制文件的控制下，分别对两种数据的对应的白天数据和夜晚数据进行最大值合成，实现数据填充。

2.2 土壤水分遥感反演方法

2.2.1 温度植被干旱指数(TVDI) Price[34]、Sandholt[35]等研究发现植被覆盖度与土壤水分之间的变化范围较大时，NDVI与LST间构成的散点图是一个三角形。并在大量研究下提出了温度植被干旱指数(TVDI)的概念，指出TVDI是利用Ts-NDVI特征空间提取的水分胁迫指标来估算陆面表层土壤水分的一种方法，Moran等[36]对Ts-NDVI特征空间的研究中发现地表最低温度与植被覆盖度和植被覆盖类型存在一定的关系。TVDI的计算公式如下所示：

(4)

式中:Ts为地表温度;Tsmax为最高地表温度;Tsmin为最低地表温度。干湿边方程如下所示。

Tsmax=a1+b1×NDVI

(5)

Tsmin=a2+b2×NDVI

(6)

式中:a1,b1,a2,b2分别为干湿边拟合方程的系数;NDVI为归一化植被指数。

TVDI的数据获取与处理过程较为简单，但NDVI对土壤背景变化较为敏感，当植被覆盖度过低，NDVI难以指示区域内植物生物量，植被覆盖度过高，NDVI增加延缓而呈饱和状态，对植被检测灵敏度下降[37]，导致反演的土壤水分与实际值偏差较大。因此，本文在NDVI高于某个值之后，采用基于增强型植被指数(EVI)的TVDI反演该区域土壤水分。

2.2.2 表观热惯量模型(ATI) 土壤热惯量是土壤的一种热特性，是引起土壤表层温度变化快慢的内在因素，与土壤水分有密切关系[37-39]。Price[34]于1985年提出了表观热惯量(ATI)的概念，考虑到在一定条件下入射的太阳辐射可视为常数，公式如下:

ATI=1-A/ΔT

(7)

式中:ATI为表观热惯量;A为全波段地表反照率;ΔT为最大地表温差。本研究采用马里大学Liang的宽波段反射率算法[38],公式为:

A=0.16a1+0.291a2+0.243a3+0.116a4+

0.112a5+0.08a7-0.0015

(8)

式中:ai(i=1,2,3,4,5,7)为MODIS数据产品对应的波段地物反射率。采用ATI模型的方法反演土壤相对湿度简单方便,但只适合于裸土和低植被覆盖区域,没有考虑到地表蒸发的影响[39]。

2.3 土壤水分关系模型的建立

本文采用的土壤目前最为广泛的线性模型，建立每天采样点的实测数据与对应的遥感反演的TVDI值或ATI值的关系模型，进而反演研究区每天的土壤水分值，然后以每天的土壤湿度值为基础合成2016年全国草地的各月内每旬的土壤水分数据[39,40-42]。

W=d·X+c

(9)

式中:W为各观测深度的土壤水分值;d为回归模型系数；c为回归模型常数项;X为表观热惯量值或温度干旱植被指数值。

3 结果与分析

3.1 确定NDVI阈值

为选取出适合本研究的NDVI阈值，参考前人阈值设置[37,43-45]，选取343组土壤水分实测数据。选取NDVI小于0.15的采样点22个、0.15～0.18的20个、0.18～0.2的40个、0.2～0.215的采样点35个、0.215～0.25的40个、0.7～0.75的51个、0.75～0.78的45个、0.78～0.8的60个、0.8～0.82的20个、0.82～0.85的10个。在不同的阈值设置下，使用不同模型进行反演，计算站点土壤水分含水量与不同深度实测土壤水分的相对误差。结果发现NDVI≤0.2时，不同深度的土壤水分采用ATI反演的土壤水分与实测值之间的相对误差小于基于NDVI的温度干旱植被模型的相对误差，且随着土层深度的加深，相对误差呈现增加趋势，即ATI模型适合于地表的土壤水分反演，随着深度的加深，ATI和TVDI的相对误差均呈现增加趋势，因此在NDVI≤0.2时采用ATI模型反演精度较高；NDVI≥0.78时，采用基于EVI的TVDI与实测数据的相对误差较小，且在土层深度小于40 cm时，其相对误差并没有随着深度的加深而加大，但当深度大于40 cm之后，基于EVI的TVDI反演的土壤水分与实测土壤水分之间的误差迅速增加。NDVI∈[0.2,0.78]时，采用基于NDVI的TVDI模型反演土壤水分精度较高(图2)。根据以上确定的阈值，对NDVI<0.2的82个点、0.20.78的90个点的土壤水分实测数据与对应的ATI、基于NDVI计算的TVDI和基于EVI的TVDI建立一元线性回归模型，公式见表1。

表1 土壤水分模型

3.2 土壤水分反演模型

通过每日的中国草地覆盖区域的NDVI值，借助于GEE(Google Earth Engine)平台通过设定不同的NDVI阈值，分别获得NDVI<0.2，0.2≤NDVI≤0.78，NDVI>0.78的遥感影像，并使用NDVI<0.2的影像对MOD09GA的B1-B7波段进行裁剪，以便于使用ATI模型反演得到土壤水分。然后对遥感反演得到的土壤水分采用平均值法将日数据合成月数据，并采用最小二乘法点对点的回归拟合各月ATI、基于NDVI计算的TVDI和基于EVI计算的TVDI和实测数据的0—10 cm的土壤水分，由于存在短暂的降雨等情况，土壤水分的变化可能不连续，在拟合过程中剔除极端值，计算得到不同NDVI值得情况下的土壤水分反演模型。可以看出，中国草地覆盖区域土壤水分于ATI，TVDI呈现除不同程度的相关性，并对各月回归拟合得到的土壤水分反演模型进行统计检验，检验结果表明各月拟合结果均通过了0.01的假设检验，具备统计学意义，模型反演效果较好(表2)。

表2 土壤水分与 ATI、基于NDVI的TVDI、基于EVI的TVDI回归拟合模型

3.3 土壤水分反演值、CLDASV2.0产品与实测值散点图

从散点图(图3)可以看出遥感反演土壤水分与不同来源地面实测土壤水分比较，本文研究结果与其具有很好的一致性，但CLDASV2.0土壤水分产品数据与地面实测数据之间偏差较大：(1) 地面实测的土壤水分集中在5%～35%，遥感反演的土壤水分也集中在5%～35%，与地面实测土壤水分的R2=0.7853(p<0.01)， RMSE=3.9841。但在25%～30%间时，大部分遥感反演值比地面观测值高，在土壤水分大于30%时，地面观测土壤水分比遥感反演值低，但大部分站点与地面实测值误差小于5%，这一误差仍在可接受范围之内；(2) CLDASV2.0土壤水分产品值集中在5%～25%，与地面实测土壤水分的R2=0.5553(p<0.01)， RMSE=7.321 8%。在5%～10%间时，大部分CLDASV2.0土壤水分产品值低，土壤水分介于15%～20%间时，CLDASV2.0土壤水分产品值远低于地面观测值。CLDASV2.0土壤水分产品值与地面实测值之间的差距较大。

图3 土壤水分遥感反演值、CLDASV2.0土壤水分产品值与站点实测数据比较

3.4 遥感反演与CLDASV2.0产品的空间尺度分析

分析2016年草地遥感反演的土壤水分含水量年均值分布图(图4A)：全年草地土壤水分集中在0%～40%。总体来看，内蒙古土壤水分高于新疆、青海、西藏和四川。土壤水分最低值分布在昆仑山脉的西藏与青海交界处、天山山脉北坡，为0～10%，新疆与西藏交界处土壤水分分布在30%～40%。同是昆仑山脉，但地域差异性尤其明显。中国南方地区土壤水分集中在30%～40%，比北方草地普遍高，这与南北降水、气温差异有关，中国南部的土壤水分常年比北方土壤水分高。与“CLDASV2.0”产品的0—10 cm的土壤水分(图4B)相比，“CLDASV2.0”产品的土壤水分在北方地区集中在10%～20%，只有西藏少部分地区、新疆塔里木盆地分布在0～10%。总体来看，遥感反演的土壤水分较“CLDASV2.0”产品数据在数值上存在较小差异，但在空间分布上保持较高的一致性。

3.5 土壤水分空间分布

通过平均法对每日土壤水分空间数据合成月土壤水分空间分布数据(图5)，发现：全国草地土壤水分存在明显的时空差异性。从全局来看，在2016年1月全国草地土壤水分由西北向东南增加，最低值(0～10%)出现在青海、西藏及新疆和内蒙古地区，最高值集中在南方各省的大部分区域。从局部分布来看，在喜马拉雅山脉、昆仑山脉、冈底斯山脉、横断山脉、天山山脉、祁连山山脉及阴山山脉地区的土壤水分较同一省份内的其他区域的土壤水分高，高出约10%。究其原因可能是该地区海拔较高，在9月、10月份，气温较低，土壤蒸发量较少，土壤保湿度比其他地区较低，1月份仍存有往年的土壤水分和积雪，因此土壤水分较高。

图4 全国草地0-10 cm土壤水分年均值分布

2月的土壤水分较1月平均高出约20%，空间变化与1月份较为相似，但很多地区的土壤水分缓慢增加，增加不明显；3月全国草地土壤水分迅速降低，4月变化区域基本是在3月的基础上缓慢增长。5月到6月内蒙古区域的土壤水分显著提高，从20%～30%提高到70%～90%，其他地区的土壤水分也从0～40%提高到30%～70%。据中国气象数据网提供的降水量在线产品数据显示，2016年5月降水量比3月、4月份的降水量增加很多；6月到7月的土壤水分下降了约10%；在7月到9月中国北方草地土壤水分普遍提高，据调查数据显示，2016年8月北方地区降雨量比往年偏多近7%，夏、秋季降水分别偏多了10%和6%，夏季暴雨过程频繁，其中，黑龙江、四川、甘肃、青海、宁夏等地区属于异常丰年，因此导致土壤水分在8月、9月大大提高。9月到11月，土壤水分显著降低，11月到12月土壤水分集中在0～50%，与前几月相比较低。

3.6 时间分布特征

通过比较土壤水分遥感反演值与CLDASV2.0土壤水分产品值的月平均值变化特征发现二者虽存在一定差异，但月际变化基本保持一致的变化趋势：(图6)全国草地土壤水分年内呈现由“升—降”的变化周期与产品数据、平均降水量变化趋于一致，特别在夏季，降水量对土壤水分有直接影响。7月、8月土壤水分比其他月高，1—6月与9—12月土壤水分偏低。其中“升—降”的第一个周期为1—7月为上升趋势，1—6月变化较小，基本在20%之内，7月、8月突增到25%；第二周期为8—12月呈下降趋势，下降趋势比1—7月的上升趋势显著。从季节变化特征来看，土壤水分最大值出现在8月，最小值出在3月。按季节看，全国整体土壤水分的最高值出现在夏季(7月、8月份)，其次是秋(9月、10月、11月份)和冬季(12月、1月、2月份)，最低值出现在是春季(3月份)。

4 讨论与结论

4.1 讨论

针对NDVI对不同植被覆盖状况下敏感性不同的问题，本文对比分析了不同NDVI阈值设置下采用不同的土壤水分反演模型对土壤湿度状况的监测效果，得到了合理的结论，NDVI≤0.2时，采用ATI模型监测土壤水分效果更好，0.20.78时，采用基于EVI的TVDI监测土壤水分精度较高。这一研究结果与前人发现的NDVI对土壤背景的变化较为敏感，当NDVI<0时，可认为地表湿度为100%，但是当植被覆盖度<20%时，NDVI对区域内的植被很难有指示意义，植被覆盖度>80%时，NDVI呈现饱和状态，NDVI值得增加出现缓慢得延迟趋势，对植被检测得灵敏度下降的研究结论基本保持一致。因此，使用单一的TVDI指数获取土壤水分干湿信息时，会由于NDVI对植被检测的灵敏度的原因而对植被覆盖较高区域效果较好，但在植被稀疏、荒漠或者裸土地区、植被覆盖度很高区域效果可能不佳。张顺谦[46]、卢远[47]、杨曦[48]等提出用EVI取代NDVI，降低NDVI易于饱和对TVDI的影响，结果发现 EVI-LST构建的TVDI特征空间更具有优势，但在植被覆盖度较低的裸土、荒漠地区的模拟精度提高并不是很好。戴睿等[49]研究发现TVDI-NDVI的相关系数总体高于TVDI-EVI。杨茹等[50]基于LANDSAT数据计算不同植被指数建立TVDI反演土壤水分，发现TVDI-NDVI的反演精度高于TVDI-EVI反演精度，因此，在直接使用EVI替代NDVI建立TVDI模型反演土壤水分的精度上仍存在争议。而ATI被很多研究证明在植被覆盖度较低的地区反演的土壤水分精度较高[37]。

因此，本研究综合考虑了ATI和TVDI模型优势，反演2012年、2016年草地土壤水分，并以实测土壤水分和反演土壤水分的相对误差为判断指标，确定NDVI阈值，建立2016年全国草地0—10 cm土壤水分反演模型，并基于实测数据对遥感反演土壤水分、CLDASV2.0产品数据的精度对比并分析2016年全国草地生态系统覆盖区域土壤水分时空变化特征、发现：基于本文提出的模型反演的土壤水分与地面实测土壤水分的R2=0.7853(p<0.01)，RMSE=3.9841，且集中分布范围较为一致；CLDASV2.0土壤水分产品与地面实测土壤水分的R2=0.5553(p<0.01)，RMSE=7.321 8%，远低于地面观测值。并从遥感反演土壤水分产品与 “CLDASV2.0”产品数据对比发现，土壤水分值上存在较小差异，但在空间分布上保持较高的一致性，两者的时间变化趋势基本保持一致，相关系数较高。就草地类型来看，草原生态系统、草甸生态系统和高寒稀疏植被和冻原生态系统分布区域，“CLDASV2.0”产品值明显低于遥感反演值。崔园园等[51]通过研究CLDAS在东北地区的适用性发现CLDAS在对极值的描述上有所欠缺；刘欢欢[52]、于月明[53]等研究发现CLDAS的土壤水分监测结果具有较高的使用价值，但在不同地区的偏差不尽相同，系统性的偏干/偏湿站点分布区域明显，呈东部偏干、西部偏湿的态势；甄熙等[54]研究发现CLDAS能够很好地反映土壤相对湿度的空间分布状况，但在数值上存在一定偏差；黄晓龙等通过实测站点对CLDAS对比发现CLDAS在某一些地区与实测值存在较大差异，误差来源主要与辐照度为受数据采集时间、局地环境和云影响有较大关系。这些研究结论能够解释本研究遥感反演产品与CLDAS在数值上存在的差异(如草原生态系统、草甸生态系统和高寒稀疏植被和冻原生态系统分布区域，究其原因可能是该区域地形复杂，云量影响等引起)，而两者在时间变化趋势和空间分布特征较为一致。

图5 全国草地生态系统覆盖区0-10 cm土壤水分空间分布

图6 全国草地土壤相对湿度月均值与各月降水量均值分布

本研究在不同植被覆盖状况下选取不同的模型反演土壤水分，虽在一定程度上提高土壤湿度状况反演的准确性，但是土壤湿度受多种因素的影响，下一步可以考虑使用数据同化、机器学习建立综合多种影响因素的土壤湿度状况监测模型，进一步提高监测的准确性。在本研究过程中主要使用MODIS产品的数据，虽能够监测大尺度的土壤湿度状况，也取得了较好的结果，但MODIS影响容易受云、雨雪天气、地形等影响，且分辨率较低，后期可以考虑基于大数据处理平台GEE上实现利用高分辨率的雷达影像建立土壤水分模型，已达到更为细致的研究。