与乔木林相比, 灌木林在整个生态系统中所占的比重较小, 但由于其对水分、养分条件要求不高, 一直以来在干旱、半干旱区有着大量的分布^[1]。因此, 灌木对于干旱生态系统的稳定与平衡起着十分重要的作用^[2]。准确掌握沙地灌木分布及动态具有重要的生态学意义。传统植被分布调查主要以地面观测为主, 该方法具有一定的局限性, 例如在人力、物力上的花费都比较大, 样地的外推也面临着较大的不确定性^[3]。随着遥感对地观测技术的不断发展, 使得大区域、大尺度的植被监测成为了可能, 而植被类型分类与覆盖度估算一直是植被遥感监测领域的热点, 已在全球不同尺度、不同地区得到了广泛的应用^[4-8]。目前, 绝大部分植被类型遥感分类以像元或若干像元的组合(对象)为基本单元, 将其划归为某一指定植被类型^[9-10], 并未考虑像元中不同植被类型混合存在的情况, 而这一问题对于中低分辨率遥感数据来说相当普遍, 特别是在地表空间异质性较强的区域^[11-12]。意识到这一问题, 估算每一像元特定植被类型覆盖度的研究得到了开展, 其中最为成熟、应用最为广泛的是全球MODIS FVC树木覆盖度产品^[13], 并且Sexton等人^[14]已通过尺度转换方法生产了全球30m空间分辨率树木覆盖度产品。值得指出的是, FVC产品中树木覆盖度主要反映的是高于5m的树木的覆盖度, 并将灌木与草本等定义为其他植被覆盖度, 显然无法满足灌木信息提取的需求, 另外其精度在以稀疏森林为主的干旱地区也存在较大疑问^[15]。针对干旱地区的特殊情况, 木本植被覆盖度(包括乔木/灌木)的估算研究得到了开展^[16], 基于时间序列MODIS数据提取的季节参数估算了整个萨赫勒地区的木本覆盖度, Higginbottom等^[17]基于Landsat数据利用机器学习的方法估算了南非较小区域的木本覆盖度。与FVC产品不同, 这些研究将乔木与灌木信息混合在了一起, 也无法满足灌木覆盖度提取的需求。因此, 发展一个大尺度、中高分辨率、适用于干旱地区沙地、可与高大乔木与草本区分的灌木覆盖度遥感估算方法具有重要的意义。

植被覆盖度估算最为常用的就是植被指数法, 一般而言植被指数与植被覆盖度之间具有较高的相关性, 但单纯利用植被指数往往无法区分不同植被类型, 并且受到荒漠地区植被稀疏及地表异质性较强的影响, 植被指数对于以根茎生长为主的灌木探测能力不足^[11]。混合像元分解法也得到了广泛应用, 该方法的优势是不需要采集大量地面样本, 然而合理的端元选取对估算精度影响较大, 导致该方法的估算结果可能存在较大误差^[18-19]。与上述方法相比, 机器学习方法运算速度更快、稳定性更高且不容易产生过拟合, 通过构建样本数据集, 结合地形、气候、植被指数、光谱信息等众多变量, 往往能构建出更为稳定的估算模型^{[14, 16, 20]}。考虑到大区域、多时相中高分辨率光学遥感数据的可用性及机器学习算法的计算需求, Google Earth Engine(GEE)遥感分析云平台已成为开展相应研究的最佳工具之一。GEE提供了校正后的全球Landsat系列数据, 并提供了机器学习算法接口^[21]。研究人员已经利用GEE开展了国家尺度植被动态监测^[22]、洲际尺度耕地制图^[23-24]、全球地表温度估算^[25]等等。然而, 尚没有针对沙地灌木覆盖度估算的相关研究。

为此, 本研究提出了充分利用GEE遥感云平台, 通过Collect Earth工具构建地面灌木覆盖度样本库, 以Landsat-8影像及其他辅助数据为变量, 通过机器学习方法开展大尺度时间序列灌木覆盖度估算的技术体系, 并以中国北方毛乌素沙地为示范区开展了应用, 并对其应用效果进行了评价, 以期为大尺度沙地灌木覆盖度估算提供参考。

1 研究区概况

毛乌素沙地位于鄂尔多斯高原东南部和陕北黄土高原以北的洼地, 覆盖内蒙古自治区的鄂尔多斯南部、陕西省榆林市的北部风沙区和宁夏回族自治区盐池县东北部小片区域^[26], 如(图 1)所示。地理位置为37°27′30″—39°22′30″N, 107°20′—111°30′E, 地处中国东部季风的尾闾地区, 属典型大陆性半干旱气候, 年降水量在250—440mm之间, 面积约为4.22万km², 在中国北方第二大沙地。毛乌素沙地流动沙地、半固定沙地分布较为广泛, 荒漠化进程迅速, 因此, 一直以来都是我国荒漠化治理的重点地区^[27]。土壤在毛乌素沙地分布具有明显的地带性变化特征。在梁地、沙地地区, 土壤主要以栗钙土和棕钙土为主；在滩地中土壤主要呈草甸土—盐化草甸土—盐土过渡的形式；在河谷阶地中, 土壤由淡栗钙土过渡为盐化草甸土。整体来看, 毛乌素沙地土壤偏沙性, 其中栗钙土、棕钙土、灰钙土、风沙土所占面积较大, 养分含量低^[28]。灌木植被是毛乌素沙地最具代表性, 同时也是分布最为广泛的植物群落。西鄂尔多斯荒漠地区大量分布超旱生灌木, 如低矮灌木狭叶锦鸡儿(Caragana stenophylla)、驼绒黎(Ceratoides latens)等, 中部及中东部干草原地区植被物种丰富多样, 并保有大量的原生灌木、半灌木植被, 这其中就包括油蒿(Artemisia ordosica)、白沙蒿(Artemisia sphaerocephala krasch)等蒿类植被以及杨柴(Hedysarum mongolicum)、柠条(Caragana korshinskii)等豆科植被^[29]。随着三北防护林、京津冀风沙源等一系列生态恢复工程的推进, 据第五次全国荒漠化和沙化监测最新数据显示, 毛乌素沙地作为重点治理区域, 其植被覆盖情况, 尤其是作为优势植被的灌木覆盖情况发生了显著变化^[30]。

2 研究方法

整体研究流程如(图 2)所示。首先, 数据准备工作包括两个部分, 其一是对多源遥感数据进行选取并进行预处理, 其二是通过两种采样法获取样本并构建样本数据集。其次, 提取样方对应多源数据, 包括光谱参数、植被指数、物候因子以及地形因子, 并对因子进行相关性分析。最后, 构建CART和SVM两种估算模型并进行精度验证, 最终产出灌木覆盖度产品。

2.1 输入因子的获取与预处理

(1) 卫星遥感数据

本研究基于GEE平台, 以陆地数据Landsat-5 TM以及陆地数据Landsat-8 OLI地表反射率影像为数据源, 利用GEE平台强大的运算能力对影像进行质量筛选、云掩膜及平均值计算, 获取毛乌素沙地2000—2017年(2012、2016除外)每年5—9月共80期地表反射率影像。在此基础上, 通过NDVI最大合成的方法生成了各年度生长季地表反射率影像。

除反射率外, 本研究计算了归一化植被指数(NDVI)、土壤调节植被指数(SAVI)、修正的土壤调节植被指数(MSAVI)、差值植被指数(DVI)^{[11, 32]}4个植被指数作为输入参数, 以更好估算灌木覆盖度。

(1)

(2)

(3)

(4)

式中, NIR为近红外波段反射率, R为红光波段反射率, L为降低土壤背景影响而设定的一个土壤调节参数, 该参数的值随植被的密度变化而变化, Huete^[33]认为在大面积稀疏植被覆盖区域L取值为0.5时对消除土壤背景影响效果最佳。

(2) 地形数据获取与处理

地形数据调用GEE平台上的SRTM数字高程模型, 地形因子提取包括高程、坡度、坡向3个因子。

(3) 气候数据获取与处理

气候数据调用TerraClimate全球陆地表面月平均气候数据集^[31]。气候因子提取包括降水、风速、帕尔默干旱指数、蒸散量4个因子作为输入参数。

输入变量汇总如表 1所示。四类因子共包含17个变量。

2.2 样本库构建

样本数据集构建是本研究的关键步骤。在此, 我们运用了Collect Earth样本收集器, 这是一款由联合国粮农组织(FAO)开发的专业影像解译软件, 主要用于陆地数据的收集、管理及分析。该软件可以轻松地访问多个免费的高时空分辨率卫星影像存档, 利用这些数据的协同作用, 实现了一种全新的土地监测方法, 称为增强视觉解译法^[34]。该方法可提升空间采样精度, 特别是在地表异质性较强的干旱、半干旱地区, 可解决植被稀疏所导致的采样困难问题。本研究对毛乌素沙地进行样本库构建, 样方以10km为间隔均匀分布在研究区内, 共415个大小为100m×100m的样方, 每个样方中等距分布有5×5, 共25个采样点, 每个采样点对应4%的地块, 采用等距取样法对每个采样方分别进行覆盖度提取, 如(图 3)所示, 通过这种方式实现覆盖度的量化。根据实际情况对样本进一步筛选, 将影像时效性较差、影像被云遮盖、影像分辨率较低的样方进行去除, 最终通过Collect Earth样本收集器获取376个有效样方。

2017年8月, 课题组依据标准一致的野外采样方案, 对毛乌素沙地开展灌木覆盖度调查, 共采集大小为30m×30m的样方24个。上述共400个样方数据为毛乌素沙地灌木覆盖度遥感估算打下了重要基础。

2.3 灌木覆盖度估算方法

为提升数据的时间匹配度, 本研究运用GEE平台, 分别对Landsat-8及Landsat-5影像进行时间序列选取、云量筛选、生长季NDVI最大值合成以及镶嵌等处理, 对2010至2015年5年间(2012年除外)获取样方对应影像波段光谱值、植被指数、地形因子以及气候因子进行提取, 在此基础上分别构建分类回归树、支持向量机两种估算模型, 并对模型估算精度进行评价。

2.3.1 CART(Classification And Regression Tree)模型

CART算法是Breiman^[35]于1984年提出的一种基于决策树构建的算法, 实质上也是归纳学习算法中的一种。其主要原理是从已有的样本库中归纳出整体的某些特征规律, 生成树形分类结构从而用于除样本之外的事例分类。基本原理是通过训练数据集构建预测模型, 采用二分递归分割方法, 使得生成的非叶子节点都有两个分支。将多个分支预测模型进行拟合, 通过预测连续变量获得树形结构模型。该算法既可用于分类, 也可用于连续变量的预测, 可有效处理大区域、大样本量和髙维度的数据。

CART算法在树形结构生长过程中, 采用经济学中基尼(Gini)系数作为选择最优测试变量和分割阈值的准则, 其数学定义如下所示：

(5)

(6)

式中, P(j|h)是从训练样本集中随机抽取的一个样本, 当某一测试变量值为h时属于第j类的概率, n_jh为训练样本中该测试变量值为h时属于第j类的样本个数, n(h)为训练样本中该测试变量值为h的样本个数, j为类别个数。

在本次研究中, 目标变量是灌木覆盖度, 测试变量是用于构建模型的各种辅助因子, 但如果仅按照上述方法生成的树形结构模型往往会出现“过拟合”的现象, 这是由于完整的树形结构模型对于样本数据集特征的描述“过于精确”其原因是受到了噪声影响, 失去了一般代表性从而无法对新数据进行下一步精确的分类, 针对该问题本研究采用K折交叉验证的方法进行叶子结点修剪, 在验证过程中引入一个“可调错误率”的概念, 即对某一树枝的全部叶子节点增加一个惩罚因子, 若该树枝仍然保持较低错误率, 则予以保留, 否则给予剪除。最终得出一颗适用于本研究且兼顾高复杂度与低错误率的最优二叉树^[36]。

2.3.2 SVM(Support Vector Machine)模型

支持向量机(SVM)通过引入核函数低维空间中线性不可分的问题映射到高维空间中, 其基本思想是通过使用位于类域边缘的支持向量, 在特征空间中的类之间构造分离超平面, 在高维空间中最大化不同类别之间的间隔寻求最佳超平面。

1997年, Vnpiak等人^[37]基于回归函数构建了用于预测连续变量的支持向量回归(SVR)模型, 即给定训练样本集T={x(x₁, y₁), (x₂, y₂), ... (x_N, y_N)}, y₁∈R, 目的是学习到一个f(x)使得其与y尽可能的接近, 其数学定义如下所示：

(7)

式中，w和b为回归函数f(x)的斜率和偏移, d为维度输入空间。

选择合适的核函数对于SVR模型的关系构建非常重要, 前人的研究中较为常见的核函数有4种：1)线性核函数：K(x_i, x_j)=x_i×x_j；2)多项式核函数：K(x_i, x_j)= (x_i×x_j+1) ^d；3)Sigmoid核函数：K(x_i, x_j) =tanh[a (x_i×x_j)+b]；4)径向基核函数(Radial Basis Function, RBF)：K(x_i, x_j) =exp(－g‖x_i－x_j‖²)。在不同核函数的作用下, 数据在高维特征空间中的分布结构将不同, 从而直接决定了SVM与核方法的最终性能^[38-39]。

本研究选择应用最广、普适性最强的RBF作为构建模型的核函数, 以基于样本点提取的各种辅助因子作为f(x)函数中的输入因子, 灌木覆盖度作为输出结果, 借助台湾大学林智仁(Lin Chih-Jen)教授开发的Libsvm工具箱^[40]进行参数选择, 该工具箱具有参数调节相对较少以及提供很多默认参数的特点, 通过网格搜索法对RBF核函数中g参数和C(惩罚)参数进行选择, 主要步骤包括：(1)按照Libsvm软件包要求的格式准备数据集；(2)对数据集进行归一化处理；(3)选用RBF函数作为核函数并运用网格搜索算法；(4)采用交叉验证方法选择最佳参数C与g；(5)采用最优参数惩罚因子C和参数g构建支持向量机模型。

2.4 精度评价

本研究样本数据库共包含400个样方, 每块样方均包含了乔木、灌木、草地及沙地四种覆盖类型。所有样方采集时间与影像信息提取时间基本保持一致, 最大限度的保证数据的时间匹配度。随机选取25%的样方对模型精度进行验证, 通过均方根误差(RMSE)、确定系数R²和估算精度(Estimation Accuracy, EA)^[41]3个指标来评价不同模型估算灌木覆盖度的精度, 计算公式为：

(8)

(9)

(10)

式中，y_i为实测灌木覆盖度, 为估算的灌木覆盖度, 为实测灌木覆盖度的平均值, 为估算灌木覆盖度的平均值, N为样方数。

3 结果与讨论 3.1 毛乌素沙地灌木覆盖情况

本研究采用Collect Earth样本收集器结合野外样地调查的方式构建样本数据集, 共获取400个样方。其中, 0—20%覆盖度的样方共205个, 占总样方数的51.3%；20%—40%覆盖度的样方共86个, 占总样方数的21.5%；40%—60%覆盖度的样方共63个, 占总样方数的15.8%；60%—80%覆盖度的样方共37个, 占总样方数的9.2%；80%—100%覆盖度的样方共9个, 占总样方数的2.3%。

从不同覆盖度样方空间分布情况来看, 如(图 4)所示, 覆盖度较高的区域主要集中在毛乌素沙地的东北缘, 有较多覆盖度大于60%的样方分布, 中东部地区主要分布覆盖度40%—60%的样方, 中部及中西部地区主要分布覆盖度低于20%的样方, 零星分布有覆盖度在20%—40%的样方。从整体上来看, 毛乌素沙地东部地区样方平均覆盖度要高于西部地区。

3.2 估算变量选择

本研究涉及到了Landsat-5与Landsat-8两种数据, 尽管Landsat系列数据具有较好地延续性, 但比较GEE平台上大气校正后Landsat-5与Landsat-8数据一致性仍具有重要意义。因此, 我们选取一块未曾变化的典型区域(包括乔木、灌木、草地、沙地), 分别利用2011年和2013年最大NDVI合成影像进行相关性分析(共计11746个像元)。通过相关性分析得到Landsat-5和Landsat-8得到的各对应波段及植被指数决定系数分别为0.94, 0.92, 0.88和0.92, 如(图 5)所示, 均达到极显著水平(P＜0.001), 表明两种传感器的反射率存在显著的相关^[42]。其中近红外和短波红外波段的斜率趋近于1, 红光波段在4种波段中斜率相对较小但也达到了0.84。考虑到数据年份不一致, 上述结果说明两种传感器对应波段的反射率具有较高的一致性, 用于长时间序列样本光谱信息提取具有可行性。

不同类型变量与灌木覆盖度相关性分析结果如(图 6)所示。光谱参数各波段均在P=0.01水平上显著相关, 其中近红外波段相关性最强；各种植被指数均在P=0.01水平上显著相关, 其中NDVI、MSAVI、SAVI相关性一致, DVI相关性相对较低；物候因子中蒸散量以及降水在P=0.01水平上显著相关, 而风速与干旱指数相关性低；地形因子中坡度及坡向P=0.05水平上显著相关, 而高程相关性低。考虑到覆盖度与遥感影像之间关系复杂, 单一变量不能满足估算要求, 因此对四类参数因子进行综合选择, 分别选取各类因子中与覆盖度相关性较高的因子。针对植被指数, NDVI、SAVI、MSAVI 3个指数与覆盖度相关性一致, 并且三者之间自相关性强, 为避免波段信息重叠对模型产生影响, 因此只选其中一个(即NDVI)作为模型参数。同时, 相关研究表明DVI对低覆盖植被具有较强探测能力, 因此也将其选作模型参数^[3]。最终选取近红外波段、NDVI、DVI、Precipitation、AET以及Slope作为构建模型的自变量。

3.3 灌木覆盖度估算模型构建及精度验证 3.3.1 CART灌木覆盖度估算模型

利用CART对由上述的自变量和因变量组成的样本进行分析, (图 7)展示了CART树形结构的一部分, 为了防止过拟合, 对模型进行K-fold交叉验证(K=5), 树形结构为：首先以凹凸点NIR≤2309.00为切割点构造二叉树, 当NIR≤2309.00时, 通过Precipitation构造下一颗二叉树；当NIR＞2309.00时, 通过NDVI构造下一颗二叉树, 树形结构不断生长, 最终通过剪枝得到一颗最优二叉树。

(表 2)统计了CART模型每个输入自变量对于构建模型的贡献度, 可以看出贡献度最大的是NIR波段, 其次是NDVI和DVI, 贡献度分别达到64.60和59.36, AET和Precipitation贡献度均为13.51而Aspect的贡献度最低仅为10.97。这个排名在一定程度上体现了每个自变量在构建模型中的权重。其中NIR波段贡献度最高, 原因是灌木属于木本植物, 主要生长方式是通过根和茎因增粗生长从而形成大量的木质结构, 因此, 与叶绿素含量较高的草本植物相比, 灌木对于红光波段的吸收能力有限, 反射的近红外光也越少。由此看来, NIR波段是区分沙地草本植物和灌木植物的最佳因子。

3.3.2 SVM灌木覆盖度估算模型

将选出的自变量代入SVM模型中, 调用Libsvm工具箱使用交叉验证的方法确定支持向量机模型参数gamma和c。如(图 8)所示, 研究中将参数gamma和c的变化范围设置在-10—10之间, 本研究选取树脂复杂度较小的RBF核函数, 它能将样本非线性地映射到一个更高维的空间, 与线性核不同, RBF核函数能够处理属性间的非线性关系。调用网格搜索参数算法并设置K-fold交叉验证(K=5), 最终得出最优惩罚系数c为512, gamma为1。

3.3.3 模型精度验证

模型验证结果如(图 9)所示。从图中可知, 两种机器学习模型均可实现一定精度的灌木覆盖度估算, CART与SVM模型的R²分别达到了0.73和0.72, RMSE分别为13.66%和13.73%, 估算精度分别为61.8%和61.6%。尽管两种模型估算精度总体相差不大, 对特定区间灌木覆盖度估算效果却有所差异。当灌木覆盖度低于50%时, SVM模型估算精度较高, 而当覆盖度大于50%时CART模型估算精度较高。

3.4 时间序列灌木覆盖度变化分析

本研究采用CART及SVM算法对2000—2015年间毛乌素沙地灌木覆盖度进行计算(2012年除外), 部分年份估算结果如(图 10)所示。从总体分布来看, 可以发现两种模型的估算结果一致性较强。对2015年估算结果进行分析表明, 灌木覆盖度高于60%的区域主要分布在毛乌素沙地的南北两缘以及东部地区, 西部地区主要以中低覆盖度为主, 中部及西南地区以流动沙地为主。统计结果显示, CART模型及SVM模型估算毛乌素沙地平均灌木覆盖度分别为28.6%和32.9%, 其中, 覆盖度低于25%的区域分别占总面积的53.6%和36.6%；覆盖度在25%—50%之间的区域分别占总面积的31.4%和34.3%；覆盖度在50%—75%之间的区域分别占总面积的11.5%和25.7%；覆盖度大于75%的区域分别占总面积的3.5%和3.4%。从2000—2015年间时间动态来看, 毛乌素沙地灌木覆盖度呈整体增长趋势, 其中南北两缘及东部地区增长趋势明显, 中部及中西部地区长期处于较低覆盖水平。

为说明本研究所发展产品特点, 我们以毛乌素沙地为例比较了灌木覆盖度与树木覆盖度产品。美国马里兰大学的Hansen团队发布了首套全球30m分辨率树木覆盖度产品^[43], 但该产品计算了高于5m的树木的覆盖度, 该产品目前已成为区域木本植被丰富性的主要数据。然而, 在以灌木覆盖为主的中国北方沙地, 该数据对大范围存在的灌木分布无法反映(图 11), 考虑到该类地区在全球及中国的广泛分布, 开展灌木覆盖度的科学估算及动态分析具有重要的应用意义。本研究提出的灌木覆盖度估算技术体系, 具备全球推广生成相应产品的潜力, 可以对树木覆盖度产品形成有效补充, 可为干旱地区生态系统管理、科学研究提供重要支撑。

4 结论

本文提出了一种利用Collect Earth进行地面样本采集、通过Google Earth Engine云平台与机器学习方法估算中高分辨率灌木覆盖度的方法, 并以毛乌素沙地为例开展了应用研究。主要结论如下：

(1) 本研究采用Collect Earth样本收集器进行等距采样的方式构建样本数据集, 解决了样本数据集偏少或样本提取精度较低的问题, 为今后更大尺度的估算提供了新的手段。

(2) 机器学习算法对毛乌素沙地灌木覆盖度估算精度均达到61%以上, 其中, 分类回归树模型确定系数R²为0.73, 均方根误差为13.66%, 支持向量机模型确定系数R²为0.72, 均方根误差为13.73%。两种模型估算精度基本持平, 但对中低和中高覆盖度估算能力各有所长, 分类回归树模型对中高覆盖度区域估算精度高于支持向量机模型, 而在中低覆盖度区域支持向量机模型估算精度高于分类回归树模型。

(3) GEE云平台存储了公开可用的全球尺度近40年的遥感影像以及其他遥感产品, 通过Google数据中心众多超级计算机的快速运算将过去需要花费几天甚至几周时间的数据处理流程在几小时甚至分秒间完成。BP级的数据存储提供了更多的数据选择性, 大大提高了数据质量, 不同的数据源也拉长了数据的时间序列, 使得今后更大尺度上的干旱、半干旱区稀疏植被估算成为了可能。