时间序列自回归模型在土壤水分预测中的应用研究

白冬妹，郭满才，郭忠升，陈亚楠

(1．西北农林科技大学 理学院，陕西 杨凌712100;2．西北农林科技大学 水土保持研究所，陕西 杨凌712100;3．中国科学院/水利部 水土保持研究所，陕西 杨凌712100)

(责任编辑 徐素霞)

黄土高原大部分地区水资源缺乏，土壤水分是限制植被建设的重要因素，一直受到人们的高度关注与重视。土壤水分含量受到很多因素的影响，如气候条件、土壤、前期土壤水分存储、土壤水分补给和消耗等［1］，具有很大的随机性。在黄土丘陵半干旱区，多年生柠条林地土壤旱化严重，为了实现土壤水资源可持续利用，需要依据土壤水分状况，调控柠条生长与土壤水关系［2］，因此监测和预报土壤水分是防旱抗旱、调控柠条生长与土壤水关系的基础。

时间序列分析主要是采用参数模型对所观测到的有序随机数据进行分析和处理的一种方法［3］，这种方法不仅可以分析因果关系，而且能根据时间序列过去的变化规律推断以后的趋势［4］。目前时间序列分析广泛应用于工农业生产、自然科学、社会科学、医学、金融及各种工程领域［4］。在土壤水分动态模拟方面，刘洪斌等［5］采用时间序列自回归建模方法对紫色土丘陵旱坡地土壤水分动态进行了模拟和预测，取得了较好的成效;杨绍辉等［3］应用ARIMA 时间序列模型对北京市大兴区芦城的土壤墒情建立模型并预测，认为ARIMA时间序列模型能很好地预测土壤墒情的变化趋势，有实际的应用价值;张和喜等［6］在干旱地区土壤墒情研究中用到了时间序列自回归模型，认为该模型可以为干旱研究及治理提供依据。本文在研究黄土丘陵半干旱地区柠条林地土壤含水量动态变化方面运用时间序列自回归模型，试图建立合理的模型并预测土壤含水量，为黄土丘陵区植物与土壤水分关系调控和植被建设提供技术支撑。

1 时间序列分析基本原理

1．1 自回归模型

自回归模型就是用自身做回归变量，q 阶自回归过程{Xt}满足以下方程［7］:

式中:μ 为时间序列的均值;φ1，φ2，…，φq为系数;时间序列{Xt}的当期值是自身最近q 阶滞后项和et的线性组合，其中et包括了序列在t 期无法用过去值来解释的所有新信息。

1．2 相关分析

在时间序列分析中常用自相关系数反映不同时期观测值之间的相互关系，t 时刻的观测值与滞后k 时刻的观测值之间的相关程度称为时间延迟k 的自相关系数，记为rk，即

一般情况下计算n/5 个自相关系数［8］，即k =1，2，…，［n/5］。

1．3 定阶及模型确定

模型定阶的准则主要有自相关函数、偏自相关函数、AIC 和BIC［7］。本文选用研究最多的赤池信息量准则(Akaike’s Information Criterion，AIC)确定自回归模型的阶数，计算公式为

式中:L 为似然函数;m 为参数的数量，如果模型包含截距项或常数项，则m=q+1，否则m=q。

根据土壤含水量时间序列的自相关分析，确定模型的可能最高阶，用最大似然估计法估计各阶模型的参数，取AIC 最小值对应的阶数作为模型的最佳阶数，最后确定AR 模型。

1．4 模型诊断

如果模型被正确识别，并且参数估计充分接近真值，那么拟合的残差应该是独立同分布的具有零均值和相同标准差的白噪声序列。拟合模型的适合性检验就是检验残差是否是白噪声序列，文中采用卡方检验法验证白噪声序列。

设{εt}为选定AR 模型估计出的残差序列，由式(2)得到残差序列样本自相关函数为

Ljung-Box 统计量［7］计算式为

2 研究区概况及土壤水分测定

2．1 试验区概况

试验区位于黄土丘陵半干旱区的上黄生态试验站(宁夏固原)，该站地理位置为北纬35°59'—36°02'、东经106°26'—106°30'。试验地坡度为0 ～15°，海拔约为1 650 m。该地区年内降水量分配不均，主要集中在6—9月，占年降水量的70%以上，年际降水量变化在634．7(1984年)～259．9 mm(1991年)之间，无霜期152 d，土壤为黄绵土，植被类型为灌丛草原。试验区主要植物有长芒草(Stipa bungeana Trin．)、阿尔泰狗娃花(Heteropappus altaicus)、茭蒿(Artemisia giraldii Pamp．)、百里香(Thymus mongolicus Ronn．)等。

2．2 土壤水分测定

试验小区为5 个立地条件基本相同的标准径流小区，每个小区长20 m、宽5 m，柠条林为2002年播种，各小区播种量分别为2．0、1．5、1．0、0．5、0 kg(对照)。在每个径流小区中部安置两个相距2 m 的PVD 套管，管长8 m，用CNC503A(DR)型智能中子水分仪测定土壤水分。测定前对中子仪进行标定，标定方程为y =55．76x+1．89，式中y 为容积含水量，x 为中子仪读数［10］。土壤水分测定深度为5—780 cm，每20 cm 测定并记录一次数据。用环刀取样，测定5 和20 cm 深处土壤水分，对表层土壤水分资料进行验证和校正。除表层5 cm 测点数值代表0—10 cm 平均值外，其他测点土壤水分值代表测点高度±10 cm 范围内土壤水分［10］。受当地气候影响，柠条4月中旬开始生长，10月份之后停止生长，直到第二年4月中旬继续生长，因此土壤含水量测定时间为:每年4月中旬到10月底每隔15 d 测定一次，11月到次年3月每月测定一次。本文选取2011年4月至2013年2月的土壤水分资料作为研究对象，对柠条林地土壤含水量动态变化进行分析及预测。

3 数据处理及结果分析

3．1 数据预处理

运用MATLAB 及R 软件处理数据。

3．1．1 数据标准化

通过计算各层土壤含水量的变异系数，发现各土层土壤变异系数存在以下规律:5、20、40 cm 土层土壤水分的变异系数范围为20% ～30%，属于比较剧烈的变异;60—120 cm 土层土壤水分的变异系数在10% ～20%之间，属于中等程度的变异;而140 cm 以下土层土壤水分的变异系数均小于10%，属于弱变异。这说明表层的土壤水分值变化较快，不稳定，而深层土壤水分值变化缓慢，比较稳定。考虑到各变异程度的土层土壤含水量变化，本研究选取20、80、160 cm 深处土层的土壤水分数据为例进行分析，建立时间序列模型。

为了满足数据分析的需要，同时提高运算精度，首先应对原始数据进行标准化处理。标准化公式为

式中:Yt为标准化后的时间序列;Xt为原始时间序列;X 是时间序列{Xt}的均值;σ 是时间序列{Xt}的标准差。以20 cm 深处土层土壤含水量序列为例，标准化后如图1。

3．1．2 平稳性

时间序列模型的分析都是建立在序列平稳的基础上的，一个平稳的随机过程有以下要求:均值不随时间变化;自相关系数只与时间间隔有关，而与所处的时间无关。实际上，自然界大多数时间序列都是不平稳的［11］，由图1 可以看出，20 cm 深处土层土壤含水量序列是非平稳的。本研究对各土层含水量序列进行了Dickey-Fuller 单位根检验，结果见表1。20、80 和160 cm 深处土层的含水量序列经Dickey-Fuller 检验后p值依次为0．504 9、0．399 5、0．551 1，均大于0．1，序列存在单位根，所以原始序列为不平稳的。通常经过差分使得序列变得平稳，差分公式为Wt=Yt-Yt-1。经过一阶差分之后，检验得到的p 值分别为0．035 8、0．014 3、0．051 0，均小于0．1，据此可以认为经过一阶差分之后的时间序列为平稳序列。

图1 标准化后的20 cm 深土层土壤含水量序列

表1 各土层含水量时间序列的平稳性检验

以20 cm 深处土层土壤含水量为例，差分后序列如图2，可以看出差分后的序列没有明显的时间趋势;有些时间点上序列波动较剧烈，这是由于20 cm 土层土壤含水量受降雨量影响较大，Dickey-Fuller 检验结果认为差分后的20 cm 土层土壤含水量序列为平稳序列，符合时间序列的建模要求。

图2 差分后的20 cm 土层土壤含水量序列

3．2 建立模型

3．2．1 模型的定阶及参数估计

应用2011年4月到2012年11月的39 组土壤含水量序列数据进行自相关分析，分析39/5≈8 个自相关系数，结果见表2。rk为滞后阶数k 时的自相关系数，当k 值取4 的时候，表2 中20、80、160 cm 三个土层的rk值都比较接近于0，故可以将模型参数估计时的最高阶数定为4。

表2 土壤含水量时间序列自相关系数

确定阶数后，运用最大似然估计法估计模型的参数，并用式(4)计算各阶次的AIC 值，结果见表3。选取AIC 最小值对应的阶次为模型的滞后阶数，可以看出20 cm 土层AIC 最小值为113．69，说明20 cm 土层含水量序列的模型为AR(1);80 cm 土层AIC 最小值为112．25，确定模型为AR(1);160 cm 土层AIC 最小值为77．35，确定模型为AR(2)。各层土壤含水量拟合模型见表4，Wt为差分后的土壤含水量时间序列。

表3 各土层含水量时间序列模型的参数估计结果

表4 各层土壤含水量时间序列模型的拟合结果

3．2．2 模型诊断

由式(1)可知AR 模型的残差为

结合式(6)对20、80、160 cm 三个土层土壤含水量序列所对应的模型残差序列进行白噪声检验，评价模型的拟合优度，模型诊断结果见表5。

表5 各层土壤含水量序列模型的诊断结果

从表5 可以看出，这3 个模型都通过了检验，可以认为构建时间序列自回归模型预测黄土丘陵区半干旱地区柠条林地土壤含水量是合理的。

3．3 模型验证

一个模型的拟合程度较好不仅在于模型诊断结果，还在于模型的实际验证。本研究选取2011年4月至2012年11月的39 组数据进行模型拟合，经模型诊断认为拟合较好，为了进一步检验模型的实用性，我们运用建立的模型预测2012年12月至2013年2月的6组土壤含水量数据。用Yt-Yt-1替换表4 中的Wt，得到预测模型，将式(6)变为Xt=σYt+X 得到原始时间序列{Xt}，20、80、160 cm 土层土壤含水量实测值与预测值及相对误差见表6。

由表6 可以看出，未来3 个月的预测数据与实测数据相对误差较小，20 cm 土层土壤含水量预测值与实测值的相对误差变化范围为-5．76% ～5．92%，80 cm 土层实测值与预测值的相对误差变化范围为0．05% ～5．28%，160 cm 土层预测值与实测值之间的相对误差范围为-5．71% ～2．66%，均小于10%，所以认为本研究建立的自回归时间序列模型能够有效地预测土壤含水量，对不同变异程度的土层土壤含水量预测效果也较好，可为黄土丘陵区土壤水分与植物关系调控和植被建设提供技术支撑。

表6 各层土壤含水量序列模型的验证结果 %

4 结 论

文中将3 个土层土壤含水量序列进行标准化及差分处理后，用Dickey-Fuller 单位根检验得出差分后的土壤含水量时间序列为平稳序列，符合时间序列的建模要求。对于时间序列自回归模型的阶数问题，本研究选用AIC 准则确定了3 个土层的阶数并建立模型，由卡方检验结果认为模型拟合得较好。最后用2012年12月至2013年2月的数据验证模型，3 个土层预测值与实测值的相对误差均小于10%，说明该模型对各变异系数的土层土壤含水量预测精度高，更进一步说明了时间序列自回归模型能够较好地预测黄土丘陵半干旱地区柠条林地土壤含水量，为该地区植物与土壤水关系的调控及植被建设提供技术支撑。

［1］尚松浩．土壤水分模拟与墒情预报模型研究进展［J］．沈阳农业大学学报，2004，35(5-6):455-458．

［2］郭忠升，李耀林．植物生长与土壤水关系调控起始期［J］．生态学报，2009，29(10):5721-5729．

［3］杨绍辉，王一鸣，郭正琴，等．ARIMA 模型预测土壤墒情研究［J］．干旱地区农业研究，2006，24(2):114-118．

［4］祖元刚，赵则海，于景华，等．非线性生态模型［M］．北京:科学出版社，2004．

［5］刘洪斌，武伟，魏朝富，等．AR 模型在土壤水分动态模拟中的应用［J］．山地学报，2003，22(1):121-125．

［6］张和喜，杨静，方小宇，等．时间序列分析在土壤墒情预测中的应用研究［J］．水土保持研究，2008，15(4):82-84．

［7］Cryer Jonathan D，Chan Kung-Sik．时间序列分析及应用:R语言［M］．第2 版．潘红宇，译．北京:机械工业出版社，2011:1．

［8］秦华光，李家才，穆丹，等．时间序列自回归模型预测茶园小绿叶蝉种群动态的探讨［J］．安徽农业大学学报，2008，35(4):564-570．

［9］于宁莉，易冬云，涂先勤．时间序列中自相关与偏自相关函数分析［J］．数学理论与应用，2007，27(1):54-57．

［10］郭忠升，邵明安．人工柠条林地土壤水分补给和消耗动态变化规律［J］．水土保持学报，2007，21(2):119-123．

［11］李庆雷，马楠，付遵涛．时间序列非平稳检测方法的对比分析［J］．北京大学学报:自然科学版，2013，49(2):252-260．