中国东北黑土区地表土壤含水量时空分异特征研究

李雪冬,红 英,费 龙,刘昀昊,王 浩,薛 泽

(1.长春师范大学地理科学学院,吉林 长春 130032;2.内蒙古师范大学地理科学学院,内蒙古 呼和浩特 010000)

0 前言

土壤含水量(Soil Water Content,SWC)强烈控制着显热和潜热通量之间的能量再分配,是地表与大气界面的重要状态参量,直接影响土壤的性质和植物生长状况,间接控制着植物分布以及生态系统小气候的变化[1]。土壤水分是气候、水文、生态、农业等领域衡量土壤干旱程度的重要指标,更是陆面过程模型、水文模型与生态过程模型的一个必要变量[2-3]。土壤含水量的时空分布与变化研究有助于进一步明确植被与土壤的关系,便于理解陆气间热量平衡与温度变化,同时也可为土地覆被空间格局、作物估产、抗旱救灾以及水土资源利用等提供决策支持。

SWC时空格局特征与旱情监测通常涉及较大区域,具有持续性和偶然性的特点,同时受外界环境作用导致空间上具有较强的变异性。传统的土壤水分测量方法[4],如烘干法、土壤时域反射仪(TDR)测量法等,虽然能够非常精确地获得SWC信息,但在短时间内获取大范围的土壤水分信息时,会受到人力与财力方面的限制。遥感技术具有时效性强、监测范围广且成本较低等优点,凭借较高的时空分辨率,可以实现长时间大区域的动态监测。众多学者基于土壤水分与植被生长状况和地表温度的关系理论,寻求基于遥感的土壤水分监测和旱情判别方法,温度植被干旱指数(Temperature Vegetation Dryness Index,TVDI)是最常用方法之一[5-6],基于MODIS数据的TVDI干旱监测方法在国际上得到了广泛的应用和认可,也是我国气象部门干旱监测业务化应用的主要方法之一。PRICE、CARLSON、GOETZ和SANDHOLT等认为NDVI和Ts特征空间呈三角形关系,高中灵等[7]在此基础上对植被指数进行对数变换,建立了ATVDI方法;HOLZMAN等[8]建立了LST与其他植被指数的关系,发现植被指数对土壤水分反演结果影响较大;刘立文等[9]引入多种植被指数,并对比了各种指数之间的关系,证实了上述结论,同时还发现利用DEM对Ts进行地形校正,反演TVDI具有更高的精度。

合成孔径雷达(Synthetic Aperture Radar,SAR)凭借穿透性强的特点,使其获取数据不受云层和光线的影响,并且对土壤表面水分高度敏感,在SWC反演中被广泛使用。丰富的模型为微波遥感反演土壤水分信息提供了物理基础和多样化的工具,但在植被覆盖条件下,微波信号的组成则变得十分复杂[10],反演植被覆盖下的SWC重点在于如何有效地分离出植被对微波散射信号的相对贡献[11]。借助植被的生物、物理特性挖掘植被散射信号的特征,并借助光学遥感数据定量反演植被参数信息是近年来研究热点[12]。BAGHDADI等[13]利用植被指数来减少植被层的散射贡献,降低土壤水分反演的不确定性;BAO等[14]使用归一化水体指数校正水云模型,减少了植被对估算SWC的影响,在植被特征参数表达后向散射信号的能力评价、模型参数的识别以及整体求解方案等方面的研究取得了一定的成果。微波数据虽然在反演SWC时表现出较强的能力,但数据获取时间和资源上还不及可见光数据。

我国东北黑土区被称为世界三大黑土区之一,主要分布在松嫩平原以及周围低台地区域,具有土层厚、有机质含量高、肥力高等特点,素有“土中之王”的美称,对于我国粮食安全有着举足轻重的意义。本研究从土壤水分含量对植被生长过程健康程度与农作物产量影响的相对重要性出发,协同Sentinel-1的合成孔径雷达数据与Sentinel-2多光谱数据,结合修正的水云模型,实现东北黑土区的SWC定量反演,分析SWC的空间分异特征;为了弥补Sentinel数据量大、获取时间序列数据困难的不足,借助时间序列MOD11A2和MOD13A3产品计算TVDI,构建黑土区SWC时间序列数据集,分析不同土地利用类型的SWC变化趋势,进而分析土壤水分的时空分布特征与气象因子的关系,挖掘土壤水分变化的主导因子。研究结果可为黑土地保护、抗旱监测、作物估产以及研究土地覆被格局特征等提供数据支撑。

1 数据源与方法

1.1 研究区

由于对“黑土”土壤类型的看法不同而导致东北黑土区的界限划定存在不同方案。本研究选取的“黑土区”为刘宝元等以土壤发生学分类的黑土、黑钙土、栗钙土和灰色森林土作为黑土区标志、采用“中心引力集聚法”反复迭代产生的东北典型黑土区[15],该区域面积为33.32万km2,共涉及138个县级行政区(表1),分为松嫩黑土亚区、蒙东黑土亚区以及三江黑土亚区。松嫩黑土亚区主要分布在松嫩平原及其四周台地低丘区,蒙东黑土亚区主要分布在内蒙古东部地区,三江黑土亚区为向东延伸出的部分(图1)。本研究选取的黑土区土壤结构良好,土层疏松软绵,呈中性或微酸性,理化性能好,素有“土中之王”的美称,对我国粮食安全起着关键作用,是重要的商品粮基地。该黑土区气候冬季寒冷,1月平均气温-6～-32 ℃,春秋两季干旱多风,年降雨量一般在400～1 200 mm,夏季雨量充沛,主要集中于6月至9月(全年降水量的70%以上)。地貌分布规律性强且类型多样,自南向北跨越了暖温带、中温带、寒温带3个不同的自然地带,从东到西横穿湿润区、亚湿润区、亚干旱区等不同的气候区,自然植被由东南向西北依次为落叶阔叶林、落叶阔叶及针叶混交林、针叶林和草原,植被生长茂盛,根系发达,在地上与地下积累了大量的有机质。

图1 研究区位置图

表1 东北典型黑土区面积统计

1.2 数据源

“哨兵”系列卫星是欧洲哥白尼空间部分(GSC)的专用卫星系列,Sentinel-1搭载了C波段SAR传感器,涉及波模式、干涉宽幅模式、条带模式和超宽幅模式等4种数据获取模式,时间分辨率为12 d,极化方式为VV和VH。利用欧空局提供的SNAP软件对Sentinel-1A的SAR数据进行斑点滤波、辐射定标和地理编码,提取相应的微波后向散射系数。Sentinel-2多光谱数据的时间分辨率为5 d,包含13个波段,本研究采用空间分辨率为20 m的植被红边波段和短波红外波段计算NDWI用于反演植被含水量。土地利用类型数据为Sentinel-2卫星图像构建的10 m分辨率数据。本研究还用到了MODIS产品数据,其中包括地表温度产品(MODIS11A2,LST)和植被指数产品(MODIS13A3,NDVI),空间分辨率1 km,数据来源于美国国家航空航天局(https://modis.gsfc.nasa.gov),研究区内产品数据的可信度均在95%以上。DEM数据来源于地理空间数据云(http://www.gscloud.cn),图像预处理包括拼接、裁剪以及坏道处理等,空间分辨率重采样至1 km,将高程与坡度利用等间隔进行分割,便于后期分析。气象数据(月平均气温和降水数据)来源于中国气象局气象数据中心(http://data.cma.cn),气温数据对海拔变化较为敏感,插值考虑到了海拔因素的影响,采用ANUSPLIN气象插值获得,这样在很大程度上消除山地区域海拔等因素对温度的影响,因此空间插值的准确率得到了极大提高。降水采用克里金插值,插值空间分辨率与MODIS产品数据一致。

1.2.2 地面实测数据

地面采样时间为2020年8月1～5日(图1),由于Sentinel-1A卫星8月2日从研究区过境,尽量保障了地面实测与卫星过境时间的同步性,样点主要分布在吉林省和黑龙江省境内。在研究区内共设置29个样方,每个样方在2 m×2 m的正方形范围内,选择对角线与中心5个位置进行取样,使用GPS记录坐标,每次取样选择土壤表面0～5 cm左右的土壤,充分混合后取杂质较少的土壤放入铝盒中称重,之后在105 ℃的烘干箱中进行充分烘干,待土壤恒重后进行称重记录,测得土壤含水量[16]。

1.3 研究方法

1.3.1 土壤含水量反演模型

1.3.1.1 TVDI计算方法

研究表明植被指数与地表温度构建的特征空间呈现三角形或者梯形分布。SANDHOLT等利用简化的LST-NDVI特征空间提出了TVDI[17],其中TVDI的表达式为:

在地质遗迹调查中发现的冰臼大多属于侧洞冰臼，就是冰臼的开口朝向旁侧，与山顶部位开口向上的冰臼相区别。朗乡花岗岩石林地质公园是侧洞冰臼分布最多的地区。

(1)

式中,T代表TVDI,L代表的是地表温度(℃);Lmax代表的是给定植被指数下的最大地表温度(℃);Lmin为给定植被指数下的地表温度的最小值(℃)。Lmax和Lmin的表达式分别为:

Lmax=a+b×ND,

(2)

Lmin=c+d×ND.

(3)

式(1)为干、湿边的拟合方程,T数值在0～1之间,其中T值越大,则代表土壤含水量越低,旱情越严重,在0.4之内为地表湿润状态,无干旱现象,当T在0.4～0.6之间为轻旱,而T大于0.6为较干旱现象,ND表示植被指数NDVI。

1.3.1.2 修正水云模型

1978年ATTEMA和ULABY[18]在估算农作物覆盖下土壤水分时提出了水云模型(Water-Cloud model),该模型假设植被覆盖层为均值的散射体,不考虑植被与地表的多次散射,因此将获取的散射信号看作是农作物的直接反射以及农作物双次衰减后的地表散射之和,因此微波后向散射可以表示为公式(4):

σC=σV+γ2σS,

(4)

式中,σC为植被覆盖地表获取的总的微波后向散射系数,σV与σS分别为植被层的后向散射系数和土壤层直接后向散射系数,γ2为微波穿透植被层的双层衰减因子。它们的表达式分别为:

σV=Amvcosθ(1-γ2),

(5)

γ2=exp(-2Bmvsecθ),

(6)

式中,mv(kg/m2)为植被含水量,θ为入射角,γ2为植被观测透过率。A和B为模型的经验参数,水云模型参数采用BINDLISHI等[19]在不同土地覆被方式中的参数结果。将式(5)、式(6)代入式(4),并变形可获得土壤层直接后向散射系数为:

(7)

植被含水量是水云模型的重要输入参数,在植被覆盖区域的土壤水分反演中具有重要的作用,研究表明植被含水量与光谱指数之间具有良好的拟合关系,本研究采用经验模型法对不同的植被含水量进行反演,其中模型参数借助JACKSON等[20]的研究结果。

1.3.1.3 模型评价

本研究反演的土壤水分数据用野外实测的数据进行校验,反演结果精度用决定系数(R2)和均方根误差(Re)进行评价,决定系数(R2)用于表征实测值和模拟值之间的拟合优度,值越接近于1,代表拟合性越好,Re反映了模拟值与实测值间的偏离程度,值越小模型的反演精度越高。

(8)

式中,Qi为实测值,Pi为模型反演值,n为验证样本数量。

1.3.2 分析方法

1.3.2.1 相关分析

本研究中的相关系数采用的是Pearson相关系数,主要用于分析TVDI与温度和降水间的相关关系,相关系数可用来表示两个变量之间的紧密程度,公式如下:

(9)

式中,x表示TVDI,y表示温度或降水,n为样本数量。rxy为所求变量相关系数,绝对值越大代表相关程度越大。采用t检验判断rxy是否显著。

1.3.2.2 趋势分析

本文应用一元线性回归法,分析黑土区2000—2020年SWC变化趋势,其变化斜率采用最小二乘法计算,公式为:

(10)

式中,n为研究时段年数;Ti为某像元第i年的TVDI均值;θs为该像元TVDI年际变化斜率,θs大于零时TVDI为增大趋势,则旱情加重,反之表示旱情减轻。

2 结果分析

2.1 模型适用性评价

本研究分别用到了基于Sentinel微波数据的修正水云模型和MODIS产品数据的TVDI两种方式反演研究区土壤水分含量,并且对比了Sentinel微波数据VH/VV不同极化方式下的后向散射系数,反演结果用实测的29个采样点对其进行验证,结果如图2所示。图2(a)、图2(b)分别为VH和VV极化方式下,总后向散射系数和考虑植被后的后向散射系数对比;图3为MODIS的TVDI反演结果精度。通过图2(a)、图2(b)可以发现,实测的数据与借助Sentinel反演的数据均在拟合值附近,说明Sentinel在反演土壤水分时具有一定的潜力,分别分析图2(a)、图2(b)可以发现,在考虑到植被对后向散射系数的影响后,修正的水云模型在VV/VH不同极化方式下的地表后向散射反演结果精度均较高,肯定了植被与地表的交叉散射在总后向散射中的贡献,在考虑植被影响后的水云模型均较原始的总后向散射系数值反演的结果有明显提升。对比图2(a)与图2(b)可以发现,无论是否考虑到植被影响的条件,VV极化条件下土壤后向散射系数的反演精度均比VH极化条件下得到的土壤后向散射系数的结果要好,其中VV在考虑植被影响的前提下,反演精度最高,R2和Re分别为0.75和0.015。基于Sentinel反演的土壤水分数据虽然可以得到较为理想的效果,但该数据的应用受数据获取时间的限制,因此本研究采用MODIS计算的TVDI,进行时间序列土壤水分数据反演。由图3分析可知,基于MODIS的TVDI反演出的土壤水分数据与实测数据也具有较强的相关关系,其中R2和Re分别为0.71和0.053。

图2 基于Sentinel不同极化方式反演土壤水分(SWC)精度验证

图3 TVDI反演土壤水分(SWC)精度验证

2.2 黑土区地表水分的空间变化特征

为了分析研究区内土壤水分的空间分异情况,本研究基于MODIS产品2000—2020年生长季TVDI的均值以及Sentinel的VV极化方式下基于修正水云模型反演的2020年土壤含水量分别生成了黑土区土壤水分的空间分布图(图4)。由图4可知基于Sentinel的反演结果与MODIS的TVDI反演结果具有相似的空间分布特征,黑土区的土壤水分情况呈现较强的空间异质性。从空间特征而言,内蒙古黑土亚区的土壤水分含量相对较低,干旱现象较为明显,在松嫩平原区腹地干旱现象较为突出,低值区域主要分布在该区域的耕地区域,在松嫩平原的周边区域,土壤湿度相对较高。整个黑土区的土壤水分含量相对较低,呈现干旱现象。

图4 黑土区土壤水分空间分布特征

2.3 黑土区土壤含水量的时间变化特征

为了讨论黑土区多年来土壤水分的变化情况以及时间序列特征,本研究分别获取黑土区2000—2020年各土地利用类型的土壤水分变化趋势,如图5所示,(a)表示的是整个黑土区,(b)～(f)分别代表的土地利用类型依次为水田、旱地、有林地、灌木林地、草地。由图5(a)可知黑土区的TVDI值随着年份的变化呈现下降趋势并且具有较强的显著性特征,在各个土地利用类型中TVDI值均随着年份均呈现下降趋势,其中旱地的TVDI值最高为0.84,灌木林地和草地的最大TVDI也均大于0.8,而水田的TVDI值虽然最高值小于0.75,但是在变化趋势上水田的下降趋势却最为明显,其次为草地和灌木林地,有林地相对变化较弱。不同植被类型区TVDI的平均值由大到小排序为旱地、灌木林地、草地、有林地、水田,这不仅反应了不同土地利用类型TVDI的差异,也强调了空间尺度的重要性。

图5 不同土地利用类型土壤水分变化情况

由于年均的TVDI值并不能全面反映出研究区内的土壤水分的时空演变特征,因此本研究借助时间序列趋势分析方法,获得了研究区内2000—2020年的逐象元的土壤水分变化趋势的空间分布,得到的趋势分布和P值(显著性检验)如图6所示。由图可知TVDI值呈现升高趋势的区域面积较小,而呈现下降趋势的区域则相对较大,研究区内的部分区域趋向干旱,其中变化最为明显的为松嫩平原的南部以及内蒙古亚黑土区的北部,而三江平原亚区则呈现波动下降的趋势。

图6 2000—2020年黑土区土壤水分变化趋势图

2.4 影响因素

为了探讨气候因素对土壤水分的影响,本研究将土壤水分数据与温度和降雨数据进行了相关分析(图7),结果表明,研究区内的TVDI 与降水量呈负相关,有50.2%的地区TVDI与降水量的负相关关系通过显著性检验,尤其是在内蒙古亚区,土壤水分与降雨的相关性显著特征较为明显。研究期间整体黑土区的土壤水分情况与气温的相关性特征并不显著,其中有80.2%的区域没有通过显著性检验,说明随着气温的升高,TVDI 的变化并不明显。就整个区域而言,在松嫩平原尤其是南部区域,土地利用类型相对复杂,并且耕地相对较多,人为干扰较为严重,从而导致干旱情况不受单一因素控制,是多种因素共同作用的结果,降水量只在一定程度上影响这些地区的干旱情况。

图7 TVDI变化影响因素分析

3 结论

本文从黑土区退化严重以及该区域对我国粮食安全的重要性出发,基于Sentinel-1微波数据,结合Sentinel-2数据,借助水云模型实现了研究区内的土壤水分反演,基于实测数据验证了Sentinel-1不同极化方式的反演精度,探讨了Sentinel在土壤水分反演过程中的潜力。因受Sentinel数据获取时间限制,本研究进一步借助MOD11A2 LST、MOD13A2 NDVI数据集,计算了TVDI干旱监测指数,并分析了黑土区的时空分布格局,讨论了不同土地利用类型长时间序列的土壤水分变化情况。研究的主要结论如下:

第一,Sentinel-1数据和MODIS数据集在监测黑土区土壤水分过程中均具有一定的潜力。在Sentinel-1微波反演土壤水分过程中,植被对后向散射系数影响较大,在去除植被影响后,Sentinel的VV/VH极化方式反演精度均有所提升,且VV极化方式表现出更大的潜力。对比Sentinel-1数据和MODIS数据可以发现,Sentinel反演的土壤水分值要低于MODIS的TVDI反演结果,主要原因可能是MODIS数据并不能有效去除植被的影响。

第二,黑土区内植被生长季的TVDI平均值为0.6,具有较强的空间异质性,其中干旱现象主要分布在松嫩平原腹地的旱地区域。分析2000—2020年的黑土区土壤水分含量的时间变化趋势可以发现,研究区内的干旱情况较为显著,但干旱有缓解的趋势,其中旱地、草地区域旱情较为明显。

第三,黑土区内的土壤水分情况受到气温和降雨的共同作用,研究区内TVDI与气温呈现正相关与降雨呈现负相关,分析研究区内的逐象元的控制因子可知,降雨相对于气温对研究区内的干旱影响更为直接,其中受降雨影响的区域较大,而受温度影响的区域较小。整体而言在松嫩平原中部受降雨的影响较为严重,而在内蒙古黑土亚区受温度的影响较为严重。

黑土SWC的时空变化特征研究对于黑土地保护、资源合理利用、农业的可持续发展、农业灾害的防治等具重要意义。借助微波的土壤水分反演可以提供较高的反演精度,但反演的过程复杂,同时受到地表粗糙度、植被覆盖度以及土壤质地等因素影响严重。微波数据的获取受时间的限制难以实现长时间序列的推演,因此如何实现快速、便捷的时间序列土壤含水量反演依旧是未来研究的热点。相对于微波数据基于MODIS的TVDI反演方法更加简单,但是基于MODIS反演的数据空间分辨率较低。本研究反演结果的验证主要是借助实测数据,就研究区的面积而言,验证点的数量略显不足,并且验证点主要分布在黑龙江和吉林省境内,因此在验证的过程中存在一定的局限性。另外,在黑土区的土壤含水量反演结果与气象因素的相关分析中,气象因子的插值过程考虑到了地形的影响,但并未单独分析地形因素对土壤含水量的作用。