吴欣睿,那晓东,臧淑英
哈尔滨师范大学,寒区地理环境监测与空间信息服务黑龙江省重点实验室, 哈尔滨 150025
土壤湿度是监测土地生产力的重要指标之一,它不仅在水文循环过程中充当了大气和陆地之间能量交换的纽带[1- 2],而且土壤湿度同农作物及植被的生长发育状况密切相关,对生态景观结构和空间格局分布也具有决定性作用[3]。松嫩平原地形平坦、面积广阔,土壤肥沃,是中国主要的商品粮生产基地之一。然而随着全球气候的变化和人为干扰的加剧,松嫩平原河湖面积缩减、极端水文事件频发、土壤盐碱化程度日益加剧、农业干旱综合态势较为严峻,严重威胁国家粮食安全[4- 7]。亟需对松嫩平原的土壤湿度及其变化趋势开展长期连续的监测和评估,分析大尺度土壤湿度的变化对区域粮食产量的影响。不同于前人研究中常采用的线性趋势分析方法、相关分析、M_K突变检验法等趋势分析方法,本文引入适用于非线性波动数据特征分析的集合经验模态分解方法对土壤湿度变化趋势进行探索。该方法在信号分析领域被广泛发展和使用,近年来已有学者将其引入到对气象和生态因子的趋势分析中:罗那那[8]和吴燕峰[9]基于集合经验模态分解分析了降水多尺度变化特征;刘可[10]利用该方法研究了不同陆地生态系统的NDVI波动特征;胡悦[11]通过此方法得出了宁夏中部干旱带的变化趋势。以上研究均表明集合经验模态分解具有良好非线性分析能力。
本文在前期研究学者的基础上,以松嫩平原为研究区,做了以下3个方面工作:(1)利用温度植被干旱指数模型获取2000—2015年时空连续的土壤湿度数据,以补充松嫩平原土壤湿度等数据资料,为以后的科研工作提供数据参考;(2)引入集合经验模态分解对时间序列连续的土壤湿度数据进行分析,更加精准的探索土壤湿度变化趋势,以期为松嫩平原的合理灌溉、农作物种植工作做出有效的指导;(3)结合土壤湿度指标和农作物产量,进一步分析土壤湿度状况对农作物产量的影响,旨在为政府及农业部门防灾减灾、确保生态农业安全提供有力的理论支持。
图1 温度植被干旱指数原理示意图 Fig.1 Schematic diagram of temperature vegetation dryness indexTs:地表温度,Surface temperature;NDVI:归一化植被指数,Normalized difference vegetation index;TVDI:温度植被干旱指数,Temperature vegetation dryness index
TVDI能够表征土壤湿度的原理,是由于植被蒸腾和土壤水分蒸发能够使地表温度降低。很多学者以此为切入点对三者空间关系进行了研究。Price,Calson等研究发现,在植被覆盖由疏到密和土壤湿度由干旱到湿润连续变化的区域,将从遥感影像中获得的归一化植被指数(Normalized differential vegetation index,NDVI)和地表温度(Ts)数据分别作为横坐标和纵坐标建立坐标散点图,得到二者空间上呈三角形关系[12- 13](图1)。
Sandholt等人在观察Ts-NDVI空间特征中发现该空间特征中含有较多等直线,于是将其简化并提出了温度植被干旱指数(Temperature Vegetation Dryness Index,TVDI)用以监测旱情及监测土壤湿度[14]。表达式如下:
TVDI=(Ts-Tsmin)/(Tsmax-Tsmin)
(1)
式中,Ts为任意像元地表温度;Tsmin为最小地面温度,对应的是湿边;Tsmax为最大地面温度,对应的是干边。Nemani、Moran等在计算Ts-NDVI的空间特征关系中表明:Tsmin和Tsmax随着不同植被覆盖条件而变化[15- 16],Sandholt将Tsmin和Tsmax进行线性拟合得到方程[15]:
Tsmax=a1+b1×NDVI
(2)
Tsmin=a2+b2×NDVI
(3)
式中,a1、b1是干边的拟合系数,a2、b2是湿边的拟合系数,并将TDVI方程进行整合最终可得到TVDI计算方程。
集合经验模态分解[17- 18](Ensemble Empirical Mode Decomposition,EEMD)是一种在经验模态分解的基础上发展而来的数据减噪分析方法,此方法能够剔除数据中不稳定噪声信息并保留数据原本所具有的内在变化特征,因此在信号特征分析领域中被广泛应用。集合经验模态分解作为经验模态分解方法的优化算法,具有自适性,其关键在于白噪声的添加解决了尺度混合问题。白噪声的添加同原始信号之间的关系如下[19]:
(4)
式中,σ为标准偏差,k为添加白噪声值,N为样本数。
图2 研究区及气象站点分布Fig.2 Study area and meteorological stations
松嫩平原位于黑龙江省西南部和吉林省西北部(42°50′—49°18′N,119°45′—129°36′E),平原以大兴安岭东麓丘陵和台地为界,北部和东部以小兴安岭及长白山地外缘山麓台地为邻,南抵达松辽分水岭(图2)。研究区地处半干旱半湿润交界地带,属于典型的温带大陆性气候,自东向西大陆性逐渐增强。年平均降水量为400—500 mm,集中在夏季,并且由东南向西北递减,年平均蒸发量为1250—1650 mm,蒸发量大而降水少使得气候干旱。平原面积广阔,约23.75×104 km2,其中耕地面积占5.59×104 km2。
松花江、嫩江及其支流常年冲积是该平原形成的主要原因,加之平原地貌主要由平原周围漫岗低丘和岗丘间洼地组成,为多种土壤的形成提供了良好条件,不同土壤类型的理化性质及其生态环境状况详见表1。该区域主要盛产大豆、小麦、玉米、甜菜、亚麻、马铃薯等,是黑龙江省和国家重要的商品粮基地,其粮食商品率所占比例高达30%以上,因此,对松嫩平原进行土壤湿度变化趋势的监测对于推动农业发展显得极为重要。
研究所使用遥感数据来源于美国国家宇航局(NASA)数据中心,由于逐日晴空资料难以获取,本次工作选用2000—2015年MODIS的3级陆地标准产品:8天合成的1 km空间分辨率的包含LST波段的MODIS11A2和16天合成的1 km分辨率包含NDVI波段的MODIS13A2,其行列号为可以覆盖松嫩平原的h26v04、h27v04。并利用MRT、ENVI对原始数据进行拼接、转换格式、投影、掩膜裁剪等预处理操作。由于LST数据为8天合成,为了使LST数据和NDVI数据有一致的时间分辨率,对LST数据相邻的两个时相利用最大取值法进行合并得到时相为16天的LST数据。然后通过式(5)对16天的LST地表温度进行校正:
Hd=Ts+H×a
(5)
式中,Hd为校正后地表温度,Ts为校正前地表温度,H为对应像元高程值,a为地表温度随高程增加而改变的变化幅度,取常数0.6℃/100 m。
表1 松嫩平原土壤类型的生态环境及理化性质简表
将相同时期对应的NDVI数据和LST数据相对应,并提取某一NDVI值所对应的LST的最大值、最小值(图3),将NDVI和LST最大最小值进行线性拟合,得到干边和湿边对应的线性拟合方程。进而利用方程(1)计算每个像元所对应的TVDI值(图3)。因为已有研究表明[20],在植被覆盖度较低时,NDVI不能够很好的反映出植被生长状况,故本研究只取NDVI>0.2的区域拟合TVDI特征方程。
图3 2015年第113天NDVI-LST线性拟合及TVDI反演结果Fig.3 NDVI-LST linear fitting and TVDI inversion results in the 113th day of 2015TVDI:温度植被干旱指数,Temperature vegetation dryness index
根据历年的NDVI与LST构成的特征空间来看,NDVI与LST最大值构成的干边呈现出显著的负相关关系(P<0.05),能非常好的拟合出干边方程,而LST的最小值随NDVI的变化较小,相关关系相对较弱,但依然能够较好的拟合线性方程。
除此以外,本研究从中国气象数据共享网上获取了2000—2015年松嫩平原23个气象站点观测的0—10 cm土壤相对湿度数据(中国农作物生长发育和农田土壤湿度旬值数据集,中国气象局陆面数据通化系统CLDAS-V2.0产品数据集),由于观测数据集未经质量控制,经过剔除筛选,最终保留20个站点的数据,并采用SPSS 22.0的最大期望算法进行数据补全。
以16年中受干旱影响较为严重的2000年为例,利用2000年4月上旬至10月下旬的0—10 cm土壤相对湿度数据与相应位置的TVDI值进行相关性分析,选取的观测站点数据时间尽量与基于遥感数据提取的土壤干湿状况的时相吻合。从图4可以看出,TVDI值与土壤相对湿度之间表现出显著的相关性,TVDI值随着土壤相对湿度的增大显示出减小的趋势。并对TVDI与土壤湿度二者的拟合结果进行F检验,通过了置信度95%的检验,这表明将TVDI作为反映土壤湿润程度的指标具有一定的合理性。
图4 相对土壤湿度和温度植被干旱指数线性拟合Fig.4 Linear fit to relative humidity and temperature vegetation dryness index
从上图显示可看出个别数据点较为离散,这是因为即便在选取相对土壤湿度数据和建立TVDI模型所使用的Ts和NDVI数据在时间上是尽可能吻合的,但是由于遥感数据是合成数据,其与地表观测数据还是会存在一定的偏差,且由于各种不可抗因素,地表观测数据和遥感影像数据不能保证空间位置的绝对匹配,其可能存在的实际空间位置的偏差对遥感反演结果也有一定的影响。尽管如此,TVDI的分布仍然能够很好的反映土壤水分在空间上的分布状况。遂用温度植被干旱指数模型逐一计算2000—2015年的TVDI值,建立时间序列上的平均干旱指数TVDImean,即每年各个时相的平均TVDI值来表示该年份中土壤水分的总体情况,及时间序列上的最大干旱指数TVDImax,即每年中各个时相的最大TVDI值来表示该年份中的极端土壤缺水干旱情况。
本研究引入集合经验模态分解方法,对2000—2015年松嫩平原土壤湿度数据进行处理。利用该方法分解原始数据所得到的残余分量能够反映土壤湿度数据时间序列上的内在变化趋势。通过Matlab执行集合经验模态分解算法分析松嫩平原2000—2015年的土壤湿度数据时,输入辅助白噪声同原始信号的标准差比率为0.02,执行样本算法的总次数为100。并采用一元线性回归及皮尔森相关系数研究分析土壤水分的总体情况及土壤水分极端缺失情况两种状态对粮食产量的影响。
从历年监测结果来看,松嫩平原的土壤湿度情况存在明显的空间分异规律:总体上自西南向东北土壤湿度水平逐渐升高。为了更方便的描述松嫩平原土壤湿度的变化情况,根据研究区的干旱分布规律,将研究区按平行于右对角线方向划分为3个区域,分别命名为西南区、中部地区、东北区。在3个子研究区中,干旱发生频率和强度相对最高的是西南区,东北区虽然有过短时期的干旱情况,总体而言相较于其他两区较为湿润。以2004年为例,全年土壤旱情较为严重且干旱频率较高的地区为西南部地区,松嫩平原的东北部地区土壤长时期处于湿润状态,且松花江流域沿岸也未受到土壤干旱化的影响(图5)。由监测结果来看,自4月份起始,土壤干旱由中部研究区的哈尔滨一带向西南部地区、西部地区东北部地区扩展,由于研究区受到冬夏季风交替控制,春天土壤干旱,需灌溉出苗,5月份的人工补水对土壤干旱起到了缓解作用,但是旱情依然较为严峻。6月迎来了该地区的雨季,直至9月,降水频繁,土壤湿润程度明显增高,土壤缺水情况得到了很大改善,且夏季植被生长旺盛也对土壤水分的散失起到了阻碍作用。该地区2000年至2015年作物生长季的TVDI多年均值为0.499,最大为0.5491,表明生长季期间松嫩平原土壤水分含量一般保持在正常状态,由于松嫩平原主要为农业区,面积广阔地势平坦,其土壤水分的变化主要受降水量变化的影响。
图5 松嫩平原2004年TVDI反演结果Fig.5 TVDI inversion results in Songnen Plain in 2004
进一步根据吴黎[21]对TVDI干旱指标的划分标准,对研究区进行土壤湿度分级(表2),并制作不同等级土壤缺水频率空间分布图(图6),分析松嫩平原全局土壤湿度水平。
总体来看,土壤缺水的时间频率在0—32.6%之间,其中,中度缺水、重度缺水和特别重度缺水16年间发生的时间频率区间分别为0—30.8%,0—29.9%,0—22.3%,且土壤缺水程度越高,发生的时间频率越低。不同等级的土壤缺水状况在空间上分异格局明显,呈现出由西南向东北递减的趋势,松嫩平原西南部和中部地区常出现土壤中度缺水,重度缺水和特别重度缺水发生频繁的地区则多集中在平原西南部的白城、松原地区。松嫩平原的风沙土主要分布于其西南的白城、松原、四平[22],由于风沙土发育微弱、且持水能力很差,极易造成土壤缺水现象,且白城、松原地处内陆,具有很强的大陆性气候,其年均降水量低于平原其他地区,而蒸散量较大[23],加剧了该区域的土壤缺水程度;另外,此处主要以莫莫格、向海为代表的湿地生态系统为主,近年来,湿地面积锐减[24]使得湿地周围土壤盐碱化现象严重,土壤保水保肥能力大大降低,也是致使松嫩平原西南部土壤缺水时间频率较高的主要原因。
表2 土壤湿度干旱等级划分
图6 土壤缺水等级发生时间频率分布图Fig.6 Time frequency distribution of soil water shortage different levels
2000—2015年TVDImean和TVDImax的折线图呈波动起伏变化(图7),研究时间范围内出现3次土壤水分缺失严重的年份,分别是2000年、2004年和2007年,这与纪仰慧等人的研究结果具有很好的一致性[25]。对时间序列上的TVDImean和TVDImax进行集合经验模态分解,最终得到土壤湿度的残余分量(图7),残余分量反映了土壤湿度时间序列的内在发展趋势。2000—2015年期间,TVDImean呈现出下降趋势,且下降速率由慢转快,表明松嫩平原整体土壤湿度呈现出变湿趋势,且变湿的速率近年来呈增加趋势;但是TVDImax变化趋势呈现逐年增高,表明土壤极端缺水导致的干旱事件呈现出逐渐升高的趋势。由于松嫩平原长期受到人为活动的影响[26],平原的原生植被多被次生植被和单一的农田代替,生态环境具有一定的脆弱性,如过度放牧使得湿地草被退化,湿地面积锐减,导致土壤盐渍化加重,土壤持水能力减弱。人类垦荒、大量的工程建设等行为对松嫩平原生态景观结构的改变,干扰了水文循环过程,破坏了土壤结构,造成了植被的退化和单一化,弱化了对土壤水分的拦截能力,从而导致极端干旱事件近年来呈现不断增加趋势。
图7 2000—2015年松嫩平原年际TVDI及其残余分量变化特征Fig.7 Variability of TVDImean and TVDImax and residual components in the Songnen Plain, 2000—2015TVDImean:平均温度植被干旱指数,Mean temperature vegetation dryness index;TVDImax:最大温度植被干旱指数,Maximum temperature vegetation dryness index
该土壤湿度的表征指标是基于植被和地表温度建立的,其与植被的长势情况密切相关,为此,本研究对历年来松嫩平原粮食产量和春夏秋季TVDI值、TVDImean和TVDImax的关系进行研究分析。松嫩平原的历年粮食产量通过查阅中国统计年鉴网发布的区域国民经济和社会发展统计公报获取,为了保证数据的准确性,研究保留了2002—2015年的粮食产量数据,中间个别缺失数据采用前两年和后两年取平均算法补全获取。研究表明,粮食产量和TVDImean呈现显著相关关系(P<0.05),即年均土壤湿度水平对粮食产量的影响显著(图8),而粮食产量和TVDImax的线性拟合未通过P<0.05显著性检验,表明粮食产量受到极端土壤缺水事件的影响较小,即只有在作物受到长期缺水的威胁,才会导致农作物减产,该结果的得出与各界的农业干旱研究结论是一致的[27](图8);就粮食产量与春夏秋季的TVDI关系分析,其与夏季土壤湿度相关关系显著,表明夏季干旱极易导致松嫩平原粮食减产(图8),而春季、秋季土壤缺水也会在一定程度上致使粮食减产,且不同时期土壤缺水对粮食产量的影响程度由深及浅依次为夏季>春季>秋季(图8)。其原因和不同物候期作物的需水量密不可分,松嫩平原作物的关键物候期分别为5月出苗期、7月抽穗期、9月成熟期[28- 29]。抽穗期间,作物植株开始分化,其茎、叶等开始迅速发育,随着叶面积的增大,作物植株代谢旺盛,消耗水量较多。若该时期缺水,则易导致作物植株分化受阻,植株矮小,粮食产量低,因此夏季土壤水分对粮食产量影响最为显著;而在出苗期和成熟期,作物对土壤水分的需求量明显减小。虽然我们的研究发现不同季相的土壤湿度对粮食产量有不同程度影响,同时也应注意到光照、温度、土壤肥力对粮食产量的重要影响,要想对作物长势和产量进行科学合理的预测还需要对其他影响因素进一步深入研究。
图8 TVDI与粮食产量的关系Fig.8 Relationships between TVDI and crop productionTVDIspr:春季温度植被干旱指数,Spring temperature vegetation dryness index;TVDIsum:夏季温度植被干旱指数,Summer temperature vegetation dryness index;TVDIaut:秋季温度植被干旱指数,Autumn temperature vegetation dryness index
本文基于MODIS影像,利用松嫩平原区域可适用的温度植被干旱指数模型对土壤湿度进行反演,引入集合经验模态分解方法,分析了松嫩平原2000—2015年土壤湿度分布特征和变化趋势,并探讨了不同土壤湿度指标对农业粮食产量的影响。研究得到以下几点结论:
(1)松嫩平原土壤湿度在空间上呈现出自西南向东北逐渐变湿润的趋势,且16年来发生土壤缺水较为严重且时间频率较高的区域出现在平原西南部,达到重度缺水和特别重度缺水的时间比例分别高达22.9%、22.3%,中度缺水的发生频率较高的区域集中在平原中部,达到30.8%。
(2)近年来松嫩平原在2000、2004和2007年土壤缺水较为严重。从年尺度来看,近16年来,该平原平均土壤湿度有逐年变湿润的趋势,但是极端土壤干旱事件发生的频率却有升高趋势。
(3)短暂的土壤极端缺水现象对农作物产量并不明显,长期的土壤湿度水平对农作物的产量影响较为显著;季节性土壤湿度变化也会影响粮食产量,其中夏季土壤缺水极易致使粮食减产,次之是春季和秋季。
本研究仅从植被指数和地温指数计算温度植被干旱指数用以代替表示土壤含水量情况,来揭示松嫩平原2000—2015年土壤湿度的时空分布情况及变化趋势,存在要素单一的问题,在日后研究工作中,应进一步综合气象要素、农业活动干预要素等多种因素进行土壤湿度的时空变化以及其对粮食产量的影响研究。本研究虽然存在不足之处,但仍然能够反映松嫩平原长时期土壤湿度的时空变化特征及趋势,以及粮食产量对不同季相土壤湿度的响应,为农业部门的合理灌溉及科学管理提供了科学有力的理论依据。