土地退化是指在不利因素的影响下土地质量与生产力下降的过程^[1]，导致土地退化的因素有很多，大致可以分为自然因素、人为因素以及社会经济因素^[2]。然而由于人类的不合理开发，尤其是矿区开采等活动，土地退化问题已经严重威胁着人类的生存和农业的持续发展^[3]。作为鄂温克族自治区中的主要露天矿区之一，伊敏矿的建设以及开采导致矿区及其附近环境的沙漠化程度由轻度敏感转向高度敏感^[4]，导致土地退化的风险有所增加。因此针对该地区开展土地退化预测具有重要的意义。

由于土地退化是一个连续且复杂的过程，多数类型的土地退化具有不可逆性，因此选取合适的预测方法，对土地退化的监测和预警具有一定的指导意义。归一化植被指数(normalized difference vegetation index，NDVI)作为遥感数据中的一种能监测全球或区域植被状态与生态环境的重要指标^[4-8]，也是衡量土地退化状况的重要指标之一，草原植被类型较为单一，利用NDVI就可以很好地反映研究区土地退化情况。目前，已有大量学者利用不同预测模型对土地退化进行了深入研究。林楠等^[9]以吉林东部为研究区域，运用灰色系统G(1, 1)模型和马尔科夫矩阵方法对该区域的植被覆盖变化情况进行模拟，并比较了两种模型的预测精度；李云鹏等^[10]构建了植被覆盖度、土壤质地和裸沙百分比的三因子综合指标评价体系，确立了不同土地退化程度的评价标准，对内蒙古近20年风蚀土地退化情况进行了分析。

本文以伊敏矿为主要研究对象，采用灰色Verhulst模型、自适应滤波法、三次样条插值法，以2000~2018年缓冲区NDVI值为基础数据，对该研究区域的土地退化趋势进行了预测分析，探讨适合矿区的预测模型，并进行土地退化等级评价，从而为西部重点矿区土地退化预测和评估提供借鉴，为人类活动提供指导。

1 研究区概况及数据来源 1.1 研究区域概况

本文主要研究区域为鄂温克族自治旗内的伊敏露天矿。鄂温克族自治旗地处内蒙古自治区东北部、大兴安岭西侧、呼伦贝尔大草原东南部^[11]，其草原面积占全旗总面积的68.9%，如图 1中蓝色框选区域所示。以伊敏矿为研究中心，以10 km为半径建立缓冲区进行分析，得出最终的研究区域如图 1中黄色框选区域所示。伊敏露天煤矿区地处呼伦贝尔草原大兴安岭西坡，海拉尔盆地东部，伊敏河中下游地区，是呼伦贝尔主要的矿区之一，但长期的露天开采使原本脆弱的草原生态系统从“结果构型破坏”演变成“功能性破坏”^[12]。

1.2 数据来源

数据集采用的是分辨率为1 000 m的SPOT/VEGETATION NDVI卫星遥感数据，利用最大值合成法生成2000~2018年的年度植被指数数据集，以及2009~2018年的月度植被指数数据集。数据获取平台为中国科学院资源环境科学数据中心。以鄂温克族自治旗为边界，用ArcMap中的按掩模提取功能裁剪2000~2014年的年度NDVI数据，得到模型对比的原始数据。以伊敏矿区边界为中心，半径为10 km以内的区域为研究区域，使用ArcMap中的掩模提取功能裁剪2009~2018每年的1月、2月、5月、6月、7月、8月、11月、12月的月份NDVI数据，得到计算该时期潜在植被生长力的原始数据。以伊敏矿区边界为中心，半径为10 km以内的区域为研究区域，使用ArcMap中的掩模提取功能裁剪2000~2018年的年度NDVI数据，得到预测2019、2020年NDVI数据的原始数据。

2 方法 2.1 模型方法介绍 2.1.1 灰色Verhulst模型

灰色预测模型的建模过程主要为先将数据进行累加得到规律性较强的数列，然后对生成模型得到的数据进行累减得到原始数据的预测值^[13]，常用于人口预测。

灰色Verhulst模型表达式：

$ x^{(1)}(k+1)=\frac{a x^{(1)}(0)}{b x^{(1)}(0)+\left(a-b x^{(1)}(0)\right) \mathrm{e}^{a k}} $

(1)

还原预测值为：

$ \begin{array}{l} x^{(0)}(k+1)=x^{(1)}(k+1)-x^{(1)}(k), \\ k=0, 1, 2, \cdots \end{array} $

(2)

式中, x⁽⁰⁾、x⁽¹⁾分别为两数列; a、b分别为两系数；e是自然常数。

2.1.2 自适应滤波法

自适应滤波预测法能较好地反映过去的规律及将来的趋势，近年来在社会经济现象的定量分析和预测中得到了广泛应用。其基本原理如下：

$ \begin{array}{c} \hat{y}_{t+1}=w_{1} y_{t}+w_{2} y_{t-1}+\cdots+w_{N} y_{t-N+1} \\ =\sum\limits_{i=1}^{N} w_{i} y_{t-i+1} \end{array} $

(3)

式中，y₁、y₂、…、y_t是某种按照时间单位依次记录下来的观测值序列；$ {\hat y_{t + 1}}$为第t+1期的预测值；w_i为第t-i+1期的观测值权数；N为权数个数；y_t-i+1为第t-i+1期的观测值。

$ w_{i}^{\prime}=w_{i}+2 k \times e_{i+1} y_{t-i+1} $

(4)

式中，i=1, 2, …, N; t=N，N+1, …, M(M为序列数据的个数，N为权数的个数)；w_i和w′_i分别为调整前后的第i个权数；k为学习常数，又称为收敛因子或自适应常数，k的大小决定权数调整的速度；$ e_{i+1}=y_{t+1}-\hat{y}_{t+1}$为t+1期的预测误差；y_t-i+1为第t-i+1期的观测值。

2.1.3 三次样条插值法

三次样条插值是广泛采用的插值技术，相比二次样条插值，三次样条插值法具有更好的拟合效果。在筑堤沉降预测研究中，Kanayama等^[14]使用三次样条插值校正数据，显著提高了预测效果。

样条曲线S(x)是一个分段定义的公式。给定n+1个数据点，n个三次多项式分段为：

$ \begin{array}{c} S_{i}(x)=a_{i}+b_{i}\left(x-x_{i}\right)+C_{i}\left(x-x_{i}\right)^{2}+ \\ d_{i}\left(x-x_{i}\right)^{3}, i=0, 1, \cdots, n-1 \end{array} $

(5)

式中，a_i、b_i、c_i、d_i为样条曲线系数。

2.2 土地退化分析方法

1) 利用ArcGIS中的Raster Calculator模块中的最大值合成法对研究区每年1月、2月、11月、12月的数据进行处理，从这4幅图中选取每个像元的最大值重新生成一幅图，作为非生长季的植被最佳生长状态图，共得到2009~2018年十幅非生长季植被最佳生长状态图，利用最大值合成法对其进行处理，得到一幅植被非生长季状态图。具体运算过程如下:

$ \begin{array}{c} \mathrm{NDVI}_{\text {非生长季 }}=\min \left(\mathrm{NDVI}_{\max (2010)}, \right. \\ \left.\mathrm{NDVI}_{\max (2011)}, \cdots, \mathrm{NDVI}_{\max (2018)}\right) \end{array} $

(6)

2) 用最大值合成法对2009~2018年研究区每年5月、6月、7月、8月的数据进行处理，共生成十幅生长季的植被最佳生长状态图。再利用最大值合成法对生成的十幅图进行处理，生成一幅草原生长季状态图。具体运算过程如下：

$ \begin{array}{c} \mathrm{NDVI}_{\text {最佳 }}=\max \left(\mathrm{NDVI}_{\max (2010)}\right., \\ \left.\mathrm{NDVI}_{\max (2011)}, \cdots, \mathrm{NDVI}_{\max (2018)}\right)-\mathrm{NDVI}_{\text {非生长季 }} \end{array} $

(7)

因为每年5月、6月、7月、8月为当地植物生长旺季，而每年11月、12月、1月、2月这4个月份为当地植物生长淡季，所以用这几个月的数据进行处理分别得到植被生长季状态图和植被非生长季状态图。

3) 通过以下公式计算某一年的相对退化率(%)：

$ \begin{array}{c} \mathrm{NDVI}_{\text {潜在 }}=\mathrm{NDVI}_{\text {最佳 }}-\mathrm{NDVI}_{\text {非生长季 }} \end{array} $

(8)

$ \text { 相对退化率 }=100 \% \times \frac{\left(\mathrm{NDVI}_{\text {潜在 }}-\mathrm{NDVI}_{\text {某年 }}\right)}{\mathrm{NDVI}_{\text {潜在 }}} $

(9)

本文参照刘钟龄等^[15]划分土地退化的等级标准，通过研究区的NDVI_潜在值计算相对退化率，得出土地退化等级。按照不同植被覆盖情况分成4个退化等级(未退化区、轻度退化区、中度退化区、重度退化区)，以此作为土地退化状况的表征。

3 数据分析与结果 3.1 模型精度评定

为了提高3个模型的分析精度以及模型的普适性，将预测区域扩大到伊敏矿区所在的鄂温克族自治旗。鄂温克族自治旗地处内蒙古自治区东北部、大兴安岭西侧、呼伦贝尔大草原东南部，占呼伦贝尔市总面积的7.39%。其草原面积占全旗总面积的68.9%，选取该区域的NDVI数据作为基础数据，代入土地退化预测系统进行NDVI的预测，分别从相对误差与绝对误差两个方面进行精度评定。

1) 绝对误差检验。利用$e(k)=\left|x^{(0)}(k)-\hat{x}^{(0)}(k)\right| $计算出2010~2015年每年的预测值与真值的差值的绝对值，因为误差均值最能反映预测值与真值之间的偏离程度，所以将计算出的6年的绝对误差求和后取平均值作为该模型的绝对误差。

2) 相对误差检验。利用公式$\Delta_{k}=\left|\frac{e(k)}{x^{(0)}(k)}\right| \times $100%计算出2010~2015年每年的相对误差，后续处理与绝对误差检验相同，将计算得到的6年的相对误差求和后取均值作为该模型的相对误差。

将上述计算得到的模型预测精度进行对比，结果如图 2、图 3所示。

3种模型预测精度差异比较明显，模型对比主要集中在自适应滤波法与三次样条插值法两个模型之间。其中自适应滤波法的相对误差平均值与绝对误差平均值较小，预测精度较高，平均相对误差为92.8%；三次样条插值法次之，为92.3%；自适应滤波法的预测精度高出三次样条插值法0.5%。自适应滤波法预测的平均绝对误差为0.052 317，比三次样条插值法小0.002 3。2010~2014年3个模型的预测误差处于正常范围，而2015年灰色Verhulst模型预测的相对误差高达0.838 150，绝对误差为0.707 156，该模型连续循环预测导致误差呈指数增长，也可能是灰色Verhulst模型比较适于描述具有饱和状态的过程，即S形过程。综合分析，自适应滤波法为最优模型，因此选用自适应滤波法进行下一步土地退化预测的研究。

3.2 土地退化类型分析

根据土地退化等级划分标准，将土地退化状态进行分级，得到土地退化程度年际变化情况及退化类型分布情况(图 4)。2018~2020年，土地退化比例逐年提高，土地呈退化趋势。

从图 4、图 5可以看出，2018年退化区域占比不到总面积的12%，鄂温克族自治旗草原区域保持着较为完好的原生植被状态，中度或重度退化区域较少。通过模型预测，2019年轻度退化和中度退化区域都不断增加，轻度退化占比高达17.56%，相比上年增幅达7%。2020年退化区域大幅度增加，土地退化率达到42%。其中轻度退化占比40.29%，相比上年增加23%，中度退化区域略微增加。从伊敏矿及矿区周围土地退化程度年际变化情况(图 5)可以看出，研究区在2018~2020年内，重度退化区域面积几乎无变化，保持平稳趋势，集中分布在研究区的东北部；伊敏矿的采矿坑周围长时间处于中度退化状态，退化面积2019年、2020年连续增加，说明土地退化严重区域正在蔓延，急需采取有效整治措施；轻度退化区域面积总体呈上升趋势，上升趋势逐渐加快，说明土地退化程度依然持续加剧。

4 结束语 4.1 结果分析

本文利用Visual Studio 2013和ArcGIS Engine 10.2进行土地退化预测系统的开发，模块功能主体为3种预测模型，分别为灰色Verhulst模型、自适应滤波法以及三次样条插值法。研究利用预测系统预测出2010~2015年共6年的NDVI数据，将其与实际的NDVI数据进行比较，选用绝对误差和相对误差对模型进行精度评定。然后利用精度较高的自适应滤波法对2019年、2020年伊敏矿及其周边植被覆盖率进行预测与土地退化评估。

1) 伊敏矿区周边土地退化情况不容乐观，中轻度退化区域不断增加。重度退化区域保持平稳且集中分布在东北部，采矿坑周围长时间处于中度退化，预测2020年退化率高达42%。

2) 自适应滤波法预测的相对误差和绝对误差普遍小于灰色Verhulst模型和三次样条插值法，预测结果较准确。经过平均绝对误差的比较发现，自适应滤波法预测的平均绝对误差最小为0.052 317，比三次样条插值法小0.002 3，比灰色Verhulst模型小0.655。从平均相对误差可以得到，自适应滤波法预测精度最高，为92.8%，三次样条插值法次之，为92.3%，灰色Verhulst模型误差较大，预测精度不高。通过3种模型的对比可知，自适应滤波法能更好地反映和预测研究区的植被NDVI值的变化。

3) 标准的Verhulst模型曲线呈S型，它反映了生物生长由加速增长到减速增长到最终趋于稳定的过程。由图 2、图 3可知，这一模型2013年之前的NDVI预测较好，胜过其他模型，但由于2013年个别地区NDVI值突变导致预测值误差极大，因而从总体预测结果来看，精度较差。

4.2 讨论

1) 自适应滤波法利用每一年的真值与预测值之间的误差不断对权数进行调整，考虑到数据的变化趋势，其得出的预测值更符合真值的整体走向，针对突变的情况有较好的适应性，所以预测误差最小。整体来看，自适应滤波法的预测精度比较平稳，效果较好。自适应滤波法应用时首先要确定权数的个数，该模型的预测精度还与模板的长度有关。经过实际研究发现，在数据量一定的情况下，权数的个数要根据数据的特性来决定，如果数据中存在周期性变动，权数的个数应根据周期数进行调整，让其为周期数或周期数的整数倍。

2) 三次样条插值法利用三次多项式函数在两两数据点之间连接，相邻小区间的局部平滑曲线串联成整体曲线。三次样条插值法利用已知数据插值出已知数据区间内的未知量，而此次研究针对预测结果的特点对该算法进行调整，根据已知数据插值出已知数据区间外的未知量。由于三次样条插值法的基本特性是两两数据进行连接，所以插值出的数据主要是由最后两个数据决定的，数据利用率小，而且如果最后两个数据之间存在突变的情况会对最后的结果造成较大的影响。本文研究结果显示三次样条插值法的预测精度一般，但是总体预测精度较平稳。

3) 灰色Verhulst模型的预测精度波动较大，误差最小为0.032 87，最大达4.221 302，2011年、2012年、2013年误差稍有增大，但整体来说不算明显，而2013年、2014年、2015年误差突增。对2012年、2013年、2014年、2015年的部分地区数据进行分析发现，造成最终误差突增的主要原因是地区的NDVI值变化幅度较大，而灰色Verhulst模型对变化幅度较大的数据比较敏感，因而会产生预测值猛增的现象。对于此次研究区的数据来说，灰色Verhulst模型预测精度不稳定，预测效果不是很理想。