基于季节性时间序列模型的西藏那曲牧草生长预测研究

汪书乐 雒伟群 赵益民 嘎桑久美

1.西藏民族大学信息工程学院，咸阳 712082；2.西藏民族大学医学院，咸阳 712082

草原是西藏生态环境大整体中关键的一个部分，草原也是牧民生活的重要依靠，目前归一化植被指数（NDVI）预测已成为一个可行的办法来评估草原环境保护的方法。《西藏自治区国家生态文明高地建设条例》要求我们保护生态环境，保护好我们的绿水青山。充分了解草场的情况，采用合理的放牧可以帮助我们更好地保护生态环境，NDVI 的预测让我们对草原的植被情况有着更好地了解，从而来更好的判断草场的生长情况，来预估草原的保护是否有效。

NDVI 是利用遥感技术对地表覆盖进行定量评估的指标之一，它是通过计算红外波段和可见光波段的比值来反映地表植被覆盖状况的指标。NDVI 在气候变化、环境监测、农业生产等领域中具有重要的应用价值。因此，如何精确预测NDVI的变化趋向，对于相关领域的研究和实践具有重要的意义。植被生长情况是一个关键的环境指标，它与气候变化、生态环境和农业生产等密切相关。对植被生长情况进行准确的预测和监测，对于环境保护、生态建设和农业生产都具有关键的意义。归一化植被指数（NDVI）是衡量植被覆盖程度的重要指标，它通过反映地表反射率中红光和近红外光的比值来反映植被的生长状况。NDVI数据能够经过遥感技术获取，对植被的生长状态进行监测和预测。

过去的研究中，在NDVI的研究方面，国内外学者采用了多种方法来预测NDVI的变化趋势。例如人工神经网络在NDVI 时间序列预测中的应用［1］，利用深度学习模型对归一化植被指数(NDVI)实现精准预测［2］，通过空间分析、统计分析和机器学习等方法，探讨了草地NDVI 与降水量、气温和日照时数的相关关系［3］，以青海湖流域2000—2014 年归一化植被指数作为输入数据(8 天，500 米)，利用时间序列模型生成预测NDVI。采用不同时间序列模型:趋势移动平均法、二次指数平滑法、差分自回归滑动平均模型(ARIMA)预测2015 年NDVI 的空间分布并进行精度比较［4］，这些方法并不能很好地解决季节性变化的问题。因此，本文采用季节性自回归移动平均模型(SARIMA)来预测西藏那曲市NDVI 的变化趋势[5-15]。本文的主要贡献在于采用了一种基于季节性时间序列模型的NDVI预测方法，该方法能够很好地解决西藏那曲NDVI 季节性变化的问题。本文将该方法应用于西藏那曲市的NDVI 预测中，并通过实验验证了该方法的有效性和准确性。

1 相关技术

研究区域：那曲市属于高海拔草原生态系统，这种类型的生态系统的特点是气候寒冷干燥，由于海拔高而缺氧。白天和黑夜温度差距大，经常刮大风。每年的5 月至9 月有一个暖和的时期，被认为是草原牧草的重要生长期。在此期间，降雨量最丰富，高山草甸和高山草甸草原类型茁壮成长。该生态系统每年的绿色植物生长期约为100 天，并且集中在这个温暖的季节。

遥感数据：MOD13Q1 和MYD13Q1 是两种种植被指数产品，源自Terra 和Aqua 中分辨率成像光谱辐射计（MODIS）传感器获取的数据。它提供了全球尺度植被活动的估计值，空间分辨率为250m。MOD13Q1和MYD13Q1 的NDVI 数据产品每16 天生成一次，两者的数据产品时间之间还会相差8 天，我们采用的数据是陆地植被产品中的NDVI，它提供了植被冠层绿度的一致时空比较，这是叶面积、叶绿素和冠层结构的复合特性。植被指数产品广泛用于各种应用，如监测植被动态、干旱评估和作物产量估算。

植被指数NDVI：NDVI 用于估计给定区域的绿色植被情况。它基于健康的植被吸收更多红光并反射更多近红外光的原理。NDVI 公式将反射的近红外光和红光之间的差值除以两者的总和：NDVI=(NIR-R)/(NIR+R)。其中NIR是植被反射的近红外光量，R是植被反射的红光量。NDVI 公式在分析植被覆盖方面具有重要意义，因为它可以区分茂密植被区域、稀疏植被区域和裸露地面区域。NDVI 的值范畴为-1～+1，其中接近-1 的值表示裸露的土壤或水，接近0 的值表示植被稀疏或没有绿色植被的区域，接近+1 的值表示茂密的绿色植被。因此，通过使用NDVI，研究人员可以识别具有不同植被覆盖水平的区域，并监测植被随时间的变化，这对于研究生态系统动态、土地利用和气候变化非常有用。综上所述，NDVI 能够评估牧草的生长情况，NDVI 公式是通过遥感影像数据分析植被覆盖的有力工具。它的意义在于它能够区分不同层次的植被覆盖，可用于研究生态系统动态、土地利用和气候变化。

时间序列模型：时间序列模型是功能强大的工具，通过可识别表现出随时间变化的数据中的模式和趋势，以及预测未来值。时间序列模型被宽泛作用于经济、气象、环境等领域的数据分析和预测中。在时间序列模型中，ARIMA 模型对于非季节性时间序列数据建模特别有用，而SARIMA 模型能够更好地处理表现出随着时间和季节性变化而变化数据。ARIMA模型通过对时间序列数据进行差分来消除非平稳性，进而通过对差分后的序列进行模型选择、参数估计和模型诊断等步骤来预测未来的数据变化趋势。ARIMA 模型，需要选择自回归项的阶数（AR 部分的参数p）、差分的次数（d 参数）和移动平均项的阶数（MA 部分的参数q）。这些参数的选择可以通过观察自相关性和偏自相关性函数来进行。参数估计能够运用最大似然估计或最小二乘法来实现。完成参数估计后，还需要进行模型诊断，以确定模型是否满足假设和是否需要进一步调整。SARIMA 模型的显示为SARIMA(p,d,q)(P,D,Q)s，其中，p：非季节性自回归项的阶数（即AR 部分的参数）；d：进行非季节性差分的次数；q：非季节性移动平均项的阶数（即MA 部分的参数）；P：季节性自回归项的阶数；D：进行季节性差分的次数；Q：季节性移动平均项的阶数；s：季节性的周期长度，要将非平稳时间序列（例如表现出趋势和季节性模式的NDVI）转换为平稳时间序列，必须使用微分方法。平稳时间序列的特点是具有恒定的平均值且没有趋势变化，而非平稳时间序列则表现出向上或向下的趋势。确保时间序列的稳定性至关重要，应用适当的变换对于实现这一目标至关重要。季节性时间序列通常表现出循环模式，需要通过微分来消除周期性变化。为此，使用了可以处理季节性时间序列数据的SARIMA 模型。在季节性时间序列模型中，时间序列数据可以被看作由一个长期趋势、一个季节性周期和一个随机波动组成的线性组合。长期趋势是时间序列数据随着时间的变化而整体呈现出的变化趋势，季节性周期是时间序列数据在特定时间范围内的重复变化趋势。

精度检验方法：

（1）MAPE:平均绝对百分比误差（Mean Absolute Percentage Error）

平均绝对百分比误差（MAPE）是用于评估预测模型的准确性或测量实际值和预测值之间的误差的指标。它表示为百分比，并计算为实际值和预测值之间的绝对百分比差异的平均值。

（2）RMSE：均方根误差（Root Mean Square Error）

均方根误差（RMSE）是用于评估回归模型或预测模型性能的广泛使用的指标。它以与原始数据相同的单位测量预测值与实际值的平均偏差。

真实值：y={yi,y2,…yn}，n是模型中的参数数量。

2 时间序列模型的实现

本部分主要介绍实验设计和结果分析。

季节性时间序列模型的构建主要分为以下几个步骤：数据收集、数据预处理、模型选择、模型训练与评估、结果分析等步骤，如图1所示。

图1 时间序列模型构建步骤流程图Figure 1 Flow chart of time series model construction steps

数据收集：从NASA官网收集2012年到2020年遥感数据产品(如MOD13Q1 和MYD13Q1)中NDVI（归一化差异植被指数）。数据包含2012 至2020 共9 个年份，每年12×4=48 项数据，数据集由9×48=432 项数据组成。

数据预处理：对收集到的NDVI 时间序列数据进行预处理，包括缺失值处理，异常值，数据平滑等。应用ArcGIS、ENVI 对其中的MOD13Q1 和MYD13Q1 指数产品处理，处理完成后将16天最大值为输入数据。

缺失值处理：ArcGIS 最大值合成法（MVC）方法用于用相邻单元的最大值替换栅格数据集中的缺失值。使用ArcGIS中的“像元统计数据”工具来应用此方法。选择两张相邻8 天的MOD13Q1 和MYD13Q1 数据作为输入栅格数据，配置叠加统计为MAXIMUM，使用最大合成法补全缺失值。得到周期为16 天的NDVI数据集。

异常值处理：利用ENVI 中band math 的工具进行筛选出大于-1 小于1，输入公式：(b1 lt 0)×0+(b1 ge 0 and b1 le 10000)×b1×0.0001，式中b1 为NDVI 波段，lt意思为小于，ge 意思为大于等于，le 意思为小于等于，整体上式子为小于0 的都乘0，大于等于0、小于等于10000 的都保留本身并乘以0.0001，因此，该公式通过将NDVI 值固定在［-1,1］范围内，并在有效范围内重新缩放值，有效地处理异常值。

数据平滑：利用Savitzky-Golay（S-G）滤波算法于NDVI（归一化植被指数）数据进行平滑处理，简单介绍S-G 滤波算法的示例表达式：smoothed_nevi=savitzky_golay(b1,window_size,polynomial_order) 其 中smoothed_nevi是平滑后的NDVI数据。b1是原始NDVI数据。window_size 是S-G 滤波的窗口大小，它用于确定在每个时间点应用的数据点数量。polynomial_order是多项式阶数，它决定了S-G滤波所使用的多项式的次数。图2 为展示批量经过MVC 合成为30 天最大值处理，异常值处理后，2012年那曲市各月NDVI数据。

图2 2012年那曲NDVI 时间序列图（250m，30天）Figure 2 Naqu NDVI time series chart in 2012(250m,30days)

模型选择：初步选择为差分自回归模型，主要是看加了季节性因素后，两个模型是否为合适的ARIMA模型和SARIMA 模型。利用赤池信息量准则(AIC)或贝叶斯方法(BIC)对所选的ARIMA 模型和SARIMA 模型进行参数估计，较小的AIC 和BIC 能够说明模型具有较好的拟合效果，经过比较AIC和BIC值比较之后，我们选择了SARIMA 模型，因为它在这两个信息准则下都具有较小的数值。这意味着SARIMA 模型更能有效地描述和预测我们的数据，从而成为最终的模型选择。AIC的计算公式：AIC=-2×ln(L)+2×k，BIC的计算公式：BIC=-2×ln(L)+k×ln(n)，其中，L是似然函数的最大值。k是模型中的参数数量。n是数据点的数量。

模型训练与评估：将预先处理后的NDVI 时间序列数据归类为训练集和测试集，采用训练集对选定的模型进行训练，并应用测试集进行模型的分析。使用评估指标如平均绝对百分比误差（MAPE），均方根误差（RMSE）来分析和评估模型的预测性能。预测NDVI的变化趋势，将模型预测结果与实际观测值进行比较，制作预测曲线和实际观测值曲线，查看它们之间的拟合水平图3、图4所示。再经过计算RMSE、MAPE来评测所选季节性时间序列模型的预测精度，表1 所示。较低的RMSE和MAPE值说明模型具备较高的预测准确性。

表1 ARIMA模型和SARIMA模型预测精度评价表Table 1 ARIMA model and SARIMA model prediction accuracy evaluation table

图3 基于ARIMA模型的预测图Figure 3 Forecast chart based on ARIMA model

图4 基于SARIMA模型的预测图Figure 4 Prediction chart based on SARIMA model

结果分析：SARIMA 模型是适合拥有季节性波动的数据来分析其时间序列的建模办法。首先，确定了数据集是由16天为一个周期的数据组成，一年便是时间周期为s=24，并提取了非季节性成分来调整每个季节周期内部值的变化。建立模型流程如下：第一步，取2012—2020 年NDVI 数据为数据集，对数据集进行差分，平稳时间序列；第二步，利用ADF测试确定序列是否平稳，并通过PACF方法确定了SARIMA模型的p和q值（p=2，q=1）；第三步，利用赤池信息量准则（AIC）或贝叶斯方法(BIC)标准来估算所选差分自回归平均模型（ARIMA）和季节性自回归移动平均模型（SARIMA）的参数，并在选择最佳模型时权衡了模型的复杂度和拟合数据的好坏；第四步，进行残差测试，利用RMSE、MAPE等指标来评估模型的准确性，ARIMA 模型中MAPE=0.045，SARIMA 模型中MAPE=0.043，ARIMA 模型中RMSE=0.0115，SARIMA 模型中RMSE=0.0108，结果表明SARIMA模型的残差误差比ARIMA模型更小，RMSE和MAPE的值较低说明SARIMA模型具有更高的预测准确性。

3 总结与展望

综上所述，实验设计和结果分析应包括数据收集、数据预处理、模型选择、模型训练与评估、结果分析等步骤，以全面评估基于季节性时间序列模型的NDVI 预测研究的有效性和可靠性。本文探索了两种常见的时间序列模型来探索NDVI 的预测精度，选择ARIMA 模型和SARIMA 模型这两种模型进行理论分析，通过NDVI 的训练集进行对比实验，实验证明SARIMA 模型具备较高的拟合精度，也因为加入了季节性参数，SARIMA 模型在具有季节性数据预测方面，精度更高。这为草原牧草生长和环境监测等方面提供了重要的参考和决策依据。

基于季节性时间序列模型的NDVI预测研究仍有很大的发展空间。一方面，可以进一步优化和改进季节性时间序列模型，如采用多元季节性时间序列模型、基于深度学习的季节性时间序列模型，基于神经网络相关的季节性时间序列模型研究等以提高不同下预测精度和预测效率。另一方面，可以将季节性时间序列模型与其他遥感数据进行融合，如降雨数据、温度数据等，用来取得更为精确和充分的结果。总之，基于季节性时间序列模型的NDVI 预测研究有着广阔的应用前景和研究价值，有望将来实现更多的功能和应用。