基于遥感光谱的干旱区土地退化评价体系构建

2019-06-20孙强强孙丹峰孙敏轩刘浩田尤淑撑刘爱霞

农业工程学报 2019年9期

张平，孙强强，孙丹峰※，孙敏轩，刘浩田，尤淑撑，刘爱霞

（1. 中国农业大学土地科学与技术学院，北京 100193；2. 自然资源部国土卫星遥感应用中心，北京 100035）

0 引言

干旱区土地退化（荒漠化）对全球粮食安全、环境质量和人类福利有着重要影响[1-4]，是21世纪以来全球所面临的重要生态环境问题之一[1]。《联合国防治荒漠化公约》第十三次缔约方大会确立了“携手防治荒漠，共谋人类福祉”的主题[2]，旨在共议公约新战略框架，维护全球生态安全。作为受荒漠化影响较严重的国家之一，2014年中国荒漠化土地面积达261.16万km2[3]。因此，及时、准确地对干旱区土地退化状况及危害程度进行评价，可以指导决策者制定优化的荒漠化防治办法和修复措施，对早日实现土地退化零增长[4]的目标具有重要意义。

相比于实地样点观测，遥感技术因其覆盖范围大、时效性强、准确度高，成为揭示干旱区土地退化状态和动态的有效工具[5]。目前，遥感直接分类[6-8]、指标指数综合[9-14]、景观格局分析[15]等是土地退化遥感评价最常用的方法。其中，运用非监督分类、机器学习等监督分类直接划分荒漠化类型和等级能够提供较为精确的荒漠化信息，但基于光谱特征的传统分类评价所获取的知识通常既不可解释，又难以传递；以植被、土壤为监测对象，进行植被指数NDVI、土壤质量监测等的单因素评价或指标聚合，可以提供一定程度的土地退化表征，却忽略了驱动力-过程-症状的内在联系，导致指标体系庞大而无法重复，限制了遥感早期预警的优势；景观格局指数对土地退化的发展和逆转有良好的指示意义，但很难充分挖掘斑块尺度的变化过程信息。而土地退化是人与自然耦合因素协同作用下土地利用/覆被类型、数量、结构以及功能的改变而引起的生态服务价值降低，具有非线性、阈值、多尺度等复杂性和动态性[16-18]，退化症状诊断需要结合关键的生物物理变量（植被-土壤的互动变化特征）来揭示荒漠化状态-机理-效应[19]。

近年来，光谱混合分解广泛应用于干旱区土地覆被分类与评价。基于多季相宽波段遥感光谱混合分解获取植被-土壤组分参数，可以简单地、直接明了地建立地表过程物理模型，进行土地利用/覆被分类制图。如：Sun等[20]采用沙、盐、植被和暗色物质 4端元丰度值构建基于先验知识的决策树，有效模拟了民勤绿洲的景观要素分布。在此基础上，通过物理意义明确的荒漠化端元丰度季相变异信息，量化下垫面的质量功能特征，为土地荒漠化监测评价提供保障。如：张熙川等[21]利用不同时相 TM 影像，采用光谱混合分解模型，以沙地、裸地、农耕地、牧草地作为基本端元组分提取反映土地退化的定量特征。卢远等[22]选用植被、盐碱地、裸沙地、沼泽图、苇地作为光谱端元进行土地退化程度的划分。Sun[23]利用端元丰度值差值监测评价了研究区的沙化、盐化和植被退化/恢复的动态变化，并结合社会经济等因素提出退化综合征分析框架。姜宛贝等[24]尝试将土地覆被/利用分类后比较和端元丰度值差值相结合，综合评价民勤生态工程第二阶段初期（2010－2015年）的治理效果。综合以上研究成果，宽波段遥感标准端元空间是将多光谱传感器光谱空间转换为具有实际物理意义的端元空间，标准端元光谱具有传感器之间传递和在其他区域可比较的优势，为不同时空条件下土地退化评价提供统一可比的基准[25]。

因此，本文在以往研究基础上，进一步建立基于标准光谱端元空间的干旱区土地退化状态评价体系，从土地利用/覆被结构和空间特征，植被-土壤组分季相变异开展结构-功能-质量间的转化规则，综合评价退化类型和退化程度状态。并以民勤2015年GF-1/WFV时间序列遥感影像进行方法体系验证，为区域荒漠化发生机理深入探究、客观稳定的荒漠化状态预警提供有效借鉴，服务于针对性、差异化的干旱区土地退化修复政策调控。

1 研究区概况和数据预处理

1.1 研究区概况

民勤县（103°03′－104°02′E；38°05′－39°06′N）地处甘肃省河西走廊东北部，位于石羊河流域下游，西、北、东三面被腾格里和巴丹吉林两大沙漠所包围。所占总面积1.59×104 km2，属典型的温带大陆性气候，气候干燥。根据国家气象信息中心（http://data.cma.cn)）统计，2015年民勤降水量为122 mm左右，多集中在夏末秋初（7－9月份），平均温度10.2 ℃，与该地区2000年以来多年平均降水量（125.75 mm）和多年平均温度（9.44 ℃）相差不大，因此所选年份可以代表常规年份。大风天数较多，环境恶劣。中部有石羊河冲积形成狭长而平坦的绿洲，绿洲区面积占16.32%，是当地水资源最丰富的区域，绝大多数人口和主要的社会经济活动都集中在中部的民勤绿洲，年内总用水量3.58亿立方米，比上年增长0.12%。青土湖水域面积达到 22.36 km2，地下水位埋深升至3.14 m，形成旱区湿地106 km2。民勤县是全国乃至世界荒漠化最严重的区域之一，也是中国北方沙尘暴四大发源地之一，土地荒漠化是限制当地农业和经济发展的主要因素之一。图1为民勒研究区。

图1 民勤研究区Fig.1 Minqin study area

1.2 数据预处理

本研究以2015年GF-1 WFV （16m）作为基础数据源，通过中国资源卫星应用中心陆地观测卫星数据服务平台（http://www.cresda.com/CN/）下载无云或少量云（＜15%）、且选取云及阴影对研究区没有影响的多光谱数据。民勤研究区由2-4景高分影像所覆盖，因此本次研究共下载43景影像（表1）。

由于GF-1 WFV Level 1级产品为辐射校正影像产品，本研究利用 ENVI 5.1中辐射定标工具（Radiometric Calibration）进行辐射定标，该工具通过读取 GF-1 WFV 元数据文件中每个波段的 Gain，offset值将 DN值转化为辐射亮度值。其次，运用FLAASH大气校正模块进行大气校正，其所采用的算法为MODTRAN5辐射传输模型，输入参数包括传感器和图像目标信息、大气模型、气溶胶模型等，参数信息可从GF-1 WFV元数据文件中获取。最后，利用无缝镶嵌工具交互式地生成一条接边线，将有地理坐标同一时期的几景图像拼接，并根据民勤行政界线通过图像裁剪提取研究区域，最终得到质量较高的的16期GF-1反射率数据（表1）。

表1 GF-1 WFV数据源列表Table 1 GF-1 WFV images data list

此外，分别于2016年9月11至9月20日、2017年8月21至8月26在研究区采集共计146个样点数据（图1），采集了土壤样品并准确的记录了各个地类的覆被类型和综合景观照片。其中，2017年50个样点用于分类评价决策树训练样本，2016年95个样点用于精度评价检验。为避免地面验证数据与获取的遥感影像时间滞后影响，采样过程中以2015年遥感影像为基准期，若采样点土地利用/覆被类型发生变化，则用其邻近的土地利用/覆被未发生变化的地点替代。

2 研究方法

图 2是本文建立的反映土地退化生物物理过程的干旱区土地退化状态遥感评价方法流程。首先从地类结构、退化类型、退化程度三个层次构建评价体系。基于标准光谱端元的多季相分析和多时相分类结果，一方面运用土地利用/覆被分类宏观表示下垫面质量结构，同时将植被-土壤互动特征辅助荒漠化发生过程及程度等土地退化状态评价。最后，运用地面采样点建立以实地样地为基础的土地退化评价基准，进行土地退化状态评价的绝对定标和交叉验证。

图2 土地退化状态遥感评价方法流程Fig.2 Flow of land degradation assessment using remote sensing

2.1 标准光谱端元的多季相分析

干旱区地物类型复杂，混合像元普遍存在。光谱混合分解将混合像元分解成不同地物的组合，分解所得丰度值为像元构成地物在混合像元中所占的百分比[26-28]，具有一定的物理指示意义。前期 TM 系列研究表明，沙（SL）、盐（SA）、植被（GV）、暗色物质（DA）4类端元与干旱区地表覆被的物理组成成分一致，是干旱区物质的基本组成成分。沙、盐和暗色物质端元指示土壤及其生境状况，其中沙、盐隐含荒漠化发生过程中裸沙、盐分特征信息，暗色物质包括阴影以及深色的土壤和岩石，而植被组分可以用来描述当地植被的盖度，4类端元可以作为标准端元库，获得时空条件下一致可比的标准端元光谱曲线。由于GF-1 WFV宽波段影像仅有4个波段（n≤b-1），同时考虑到干旱区地表的复杂性和空间异质性，本文拟采用多端元光谱混合分解。多端元光谱混合分解通过将像元的光谱信息和空间信息相结合，逐像元调整端元的类型和数量动态地进行混合像元分解，比固定端元线性光谱混合模型更合理且分解精度更高。

在 Sun等[20]、刘娜等[29-30]多季相线性光谱混合分解方法的基础上，本研究运用时间序列 GF-1 WFV宽波段影像进行多端元光谱混合分解的方法流程如下。

1）选取与干旱区地表物质组成一致的标准端元（沙、盐、植被、暗色物质）构建标准端元库，任意 3类端元组合，构成4种标准端元候选集（如图3）。考虑到地物端元光谱时间序列变异性和可比性，选择每种纯净端元的代表性季相（2-3期）的平均光谱曲线作为标准端元空间的光谱曲线，建立包含 4类端元组合候选集光谱曲线的标准端元光谱库。

图3 标准端元候选集Fig.3 Candidate dataset of standard endmembers

2）运用所建立的标准端元光谱库（包含4种标准端元光谱组合），逐像元进行全约束线性光谱混合分解，并且基于 RMSE最优拟合原则逐像元确定最优分解模型，通过python编程实现自动化的最优模型分解。最终得到时间序列各端元丰度图和相应的残差均方根图。

3）选取春、夏、初冬典型三季相（4月，8月，11月）三季相合成R、G、B三色度图像，分析得到各端元丰度值的空间分布情况和年内端元丰度值的季相变异规律，为土地利用/覆被分类与土地退化状态评价提供基础。

2.2 标准光谱端元的多时相分类

土地利用/覆被系统结构-功能的变化势必会带来生态服务价值的改变。本文首先基于多时相标准光谱端元空间的分类要求和土地退化遥感的精细化监测需要，建立了三级的分类体系。其中一级反映气候水文地貌控制下的下垫面性质和主要生态服务类型。二级根据土地利用/覆被季相对比特征表征自然/人工植被群落。三级续分出主导植被/作物类型，并能反映植被-土壤生境的互动过程（表2）。其次，构建了4端元16期年内时间序列变化曲线，采用典型季相的端元组分进行三原色色度空间分析获取端元组分季相特征，和时间序列分析方法获取更详细的端元季相分类知识，提取典型地类端元光谱时间序列特征参数，例如：时间序列曲线最大值、最小值以及对于地物分类敏感的时间窗口（例如生长期）均值、极值、变化率、差值等，量化旱地系统“结构-状态-功能”用以揭示植被及其土壤生境的生态互动过程，建立分类知识库（表2）。在分类知识挖掘基础上，进行决策树的构建，从根部到树枝逐步从原始影像中分离出每一种土地利用/覆被类型。该分类决策树基于具有实际物理意义的端元丰度值建立，具有较好的稳定性和可移植性。最终共获得22种3级地类分类结果，分类结果如图4所示。结合2016年和2017年实地采样点位（图1）和高分辨率Google Earth影像，利用验证样本（约3800像元）对分类结果进行精度验证，总体精度达88.31%，为后期的土地退化状态评价奠定基础。

表2 多时相土地利用/覆被分类知识库Table 2 Knowledge database for multi-temporal land use/cover classification

图4 2015年民勤县标准光谱端元的多时相分类结果Fig.4 Multi-temporal land use/cover classification results based on standard EMs spectra of Minqin in 2015

2.3 土地退化状态评价及其定标验证

2.3.1 评价体系构建

干旱区土地退化状态评价分别从土地利用/覆被分类、退化类型和退化程度三个层面建立干旱区土地退化遥感状态评价体系。

1）第一级根据土地利用目的的不同，从土地利用/覆被的角度对干旱区土地退化进行分类与评价，不同土地覆被类型生态功能不同，表征其生态系统的稳定性和下垫面的承载力。

2）第二级以下垫面土壤属性发生过程为依据，进行土地退化类型的划分。一般来说，可以将下垫面退化发生过程分为未退化、沙化、盐化、沙-盐化四种过程。

3）第三级以土地退化的轻重程度为指标，对同一演替序列、不同演替水平的土地单元进行的阶段划分，一般分未退化、轻度退化、中度退化、重度退化四级。

2.3.2 土地退化类型状态评价

精细的土地利用/覆被分类是土地退化状态评价的基础和前提。以决策树构建的基于先验知识的土地利用/覆被分类体系以各端元（植被-土壤组分）为分割变量，可以系统展示各个地类生态属性特征（季相特征、植被土壤特征），为区别不同利用/覆被条件下的主要退化类型提供依据。因此，利用分类知识（表2）得到识别不同土地利用/覆被条件下的主导控制过程。对于存在不同退化类型的土地利用/覆被类型，本文基于土地退化状态三级评价体系，重点分析各土地利用类型下在典型季相（5,8,10为春、夏、冬的代表性季相）特征指标的状态及变异规律，运用该退化类型相对应的时期和植被组分的时间窗口选择阈值进行不同退化类型的划分（表3）。由于本研究所选用的状态评价变量为具有实际物理意义的端元丰度值，因此可以利用先验知识选取目标地类下主导退化类型的特征端元作为相应节点的分割变量，例如：盐化主导类型为盐端元，具体的阈值确定根据训练样本获得。最终运用决策树方法进行退化类型的划分，确定未退化、沙化、盐化和沙-盐化4种主导退化类型。

表3 土地退化类型评价体系Table 3 Evaluation system for types of land degradation

从表3中可以看出，不同土地利用/覆被条件下的土地退化类型状态各不相同。除林地外，受自然因素和人为活动耦合作用的地类，如耕地、园地、因长期不合理的利用，也会导致一定程度的退化过程，其主要退化状态类型为沙-盐化和盐化。人类干扰小的地类，受当地自然资源条件禀赋限制，如荒草地根据其所在生境不同，退化状态也不同；沙地、裸地主要表现为沙化，盐碱地主要表现为盐化等。

2.3.3 土地退化程度状态评价

不同的土地利用方式由于受人为干预程度不同，所表征的土地退化发生过程不同，这就决定了不同退化类型的程度评价指标阈值存在差异。本文根据民勤土地退化研究现状[31]，考虑到极重荒漠化和重度荒漠化差别不大，将沙化、盐化、沙-盐化按照轻度、中度、重度分级组合，在土地退化类型评价的基础上，建立了 16种不同退化发生过程的荒漠化等级评价指标及其评价规则（表4）。具体以表征植被生长期覆盖度和指示生境的沙/盐端元丰度的指标互为辅助，根据实地样点定标后确定各退化等级划分阈值，运用决策树方法实现退化程度的划分。

表4 土地退化程度评价体系Table 4 Evaluation system for degree of land degradation

2.3.4 土地退化评价结果验证

将2016年9月在民勤设置的95个地面样点作为主要结果验证数据。每个样点控制30 m×30 m的样方，将其划分为2×2个子样方，并取4个子样方的表面土壤（0～20 cm）混合作为样点土壤，测量其电导率（EC）和土壤酸碱度（pH值）指示土壤的盐化和碱化程度。同时获取了立地景观照片以供土地利用/覆被分类结果和土地退化评价综合描述与验证。民勤地类复杂多样，土地退化影响因子较多，本研究参考FAO/UNEP荒漠化状态评价指标体系并结合 Sun等（2016）研究成果[32]，建立民勤基于实地样点的荒漠化状态指标体系（表5），并与宽波段遥感标准光谱端元空间的评价进行交叉验证。

3 研究结果与分析

3.1 土地退化类型状态评价结果

基于土地退化类型状态评价体系和决策树分类法得到 2015年民勤不同土地利用/覆被类型控制下的退化状态类型分布图（图5）。由图可知，沙化过程（面积占比85.81%）、沙-盐化过程（面积占比 5.14%）是民勤地区主导的土地退化过程。沙化过程主要分布在腾格里沙漠和巴丹吉林沙漠及其与绿洲边缘交界处，农田集约化过程使得处于交界处的边缘脆弱区在非耕种区土壤裸露，同时林地向耕地转换、种植结构（春季作物向夏季作物转换）等土地利用行为也会引起沙化风险。沙-盐化位于西北部民勤洪积扇和东北部下游低洼区域，由于平均气温增高和蒸发量的增加加快了岩石的风化速度，降水量增加使得盐分随水流和风力搬运至山脚，在洪积扇处和下游低洼处大量积累，同时下游东北部由于农田灌溉和水流运输，盐分积累使之形成盐湖。在沙和盐分随季节的增减互动下，两个区域形成沙-盐化的状态格局。盐化过程在气候和地形的综合作用下，除西北部民勤洪积扇和东北部下游低洼区域盐分主导形成盐化集中区域外，部分由于土地利用不合理井水灌溉导致农田区域出现盐化现象（面积占比0.81%）。未退化地区（面积占比8.24%）大多位于民勤绿洲。

图5 民勤县土地退化类型状态评价分布图Fig.5 Land degradation types distribution mapping of Minqin county

对每一土地利用方式下的土地退化类型进行面积统计，揭示各不同土地利用方式下的荒漠化主导发生过程，为针对性地指导土地退化防控和差别化的政策管理提供参考。由图 6可知，耕地、园地、林地、水域等生态服务较积极的地类，以未退化占主导，但由于长期耕作利用，致使耕地部分出现沙-盐化，而水域和林地部分出现沙化风险。未利用地以沙化过程为主导，覆被类型大多以沙地、盐碱地、盐渍化沙地、戈壁等构成，生境条件恶劣，难以利用。灌草地虽然以沙化过程为主导，但由于生长灌木、小灌木、半灌木等荒漠植被，具有改善地区环境的能力。

图6 2015年民勤县各土地利用类型下退化类型面积比例Fig.6 Land degradation type area distribution map for every land use types of Minqin in 2015

3.2 土地退化程度状态评价结果

由土地退化程度分布图（图7）及相应面积统计可知，2015年民勤县轻度沙化（面积占比60.30%）、中度沙化（23.64%）。其中，轻度沙化、中度沙化主要分布在腾格里沙漠处，此外，未退化（面积占比8.24%）、轻度沙化轻度盐化（面积占比2.67%）主要分布在民勤绿洲及其周围区域，由于当地降水少（约100 mm）和地下水位较高（约15～18 m），绿洲边缘自然植被覆盖度较低，以荒草地、灌木、半灌木和小灌木等不稳定状态存在，加之部分地区不合理的耕作（抛荒、弃耕等），导致土壤存在一定裸露而被风蚀。

图7 土地退化程度分布图Fig.7 Distribution of land degradation degree

从绿洲周围区域轻度盐化，到洪积扇处和下游低洼处逐渐扩大的盐渍化的空间格局是民勤水资源空间再分配的结果。温度升高和降水量的增加使得山脉和平原交界处和下游低洼处进一步增加盐渍化的风险，而人类的土地利用活动（大面积夏季/秋季作物的种植）和管理措施（关井、地表水灌溉等）有利于土壤盐分随水流迁移至更深的地下深层和下游的低洼处（青土湖），导致农田及其边缘盐渍化风险大大降低。Sun等[31]的研究结果显示，97.35%的农田土壤为非盐化土壤，对于民勤作物生长期，占轻度盐化的农田仅为2.62%，与上述研究结果基本一致。同时相比于上游地区，下游地区更依赖盐碱度高的地下水灌溉，导致下游地区退化程度大于上游地区。

3.3 精度评价与验证

精度评价是土地退化状态评价的保障，本文利用实地采样点结合地面立地景观照片，运用混淆矩阵（表）对评价结果进行验证。结果显示，利用遥感技术开展基于土地利用/覆被分类和植被-生境参数的干旱区土地退化状态评价识别土地退化类型的能力为87.5%，识别土地退化程度的能力达到78.5%（表6），具有挖掘土地退化驱动力-过程-症状的潜力。

表6 精度评价混淆矩阵Table 6 Confusion matrix of accuracy assessment

4 讨论

本研究首次将光谱混合分解模型拓展到国产 GF-1 WFV传感器，并发展了适用于此卫星传感器多端元混合分解模型，克服了旱地系统景观局部变异性大的问题。通过对比姜宛贝等[30]对2015年OLI端元空间的分析，发现GF-1系列影像相比Landsat8在时间和空间分辨率上均有较大提升，能够得到年内高频次端元时间序列曲线[33]，这对后期获取更详细的端元季相分类知识具有重要作用。基于标准端元空间，本文充分挖掘斑块尺度上植被和土壤的时间序列互动变化过程，建立了表达年内退化过程时空差异的关键知识库，并运用决策树对不同土地利用/覆被结构进行可视化表达。相比于 Sun[23]、姜宛贝等[30]对民勤2008年和2015年的分类结果，本文在区分荒漠植被、耕地作物水平能够获取稳定可靠的结果，实现了更精细、准确的土地利用/覆被分类。相比于传统的依靠庞大指标体系土地退化状态评价，本文通过构建“退化类型-退化程度”三层次土地退化状态评价体系，运用土地退化诊断关键指标（植被端元、沙端元和盐端元）实现了不同土地退化类型和土地退化程度的快速、有效地识别和评价，所获得的评价结果具有较强可解释性和可对比性，可以指导土地退化修复的精细管理和高效决策。

本研究还存在一些不足之处，有待进一步提升完善。1）决策树阈值划分过分依赖于专家知识和训练样本，“自下而上”和“自上而下”自动化调控水平有待提升。2）本文仅基于一年数据进行土地退化状态的相对评价，该技术方法在更大范围的时空尺度拓展应用将成为今后研究的重点。