時(shí)間序列自回歸模型在土壤水分預(yù)測(cè)中的應(yīng)用研究

白冬妹，郭滿(mǎn)才，郭忠升，陳亞楠

(1．西北農(nóng)林科技大學(xué) 理學(xué)院，陜西 楊凌712100;2．西北農(nóng)林科技大學(xué) 水土保持研究所，陜西 楊凌712100;3．中國(guó)科學(xué)院/水利部 水土保持研究所，陜西 楊凌712100)

(責(zé)任編輯 徐素霞)

黃土高原大部分地區(qū)水資源缺乏，土壤水分是限制植被建設(shè)的重要因素，一直受到人們的高度關(guān)注與重視。土壤水分含量受到很多因素的影響，如氣候條件、土壤、前期土壤水分存儲(chǔ)、土壤水分補(bǔ)給和消耗等［1］，具有很大的隨機(jī)性。在黃土丘陵半干旱區(qū)，多年生檸條林地土壤旱化嚴(yán)重，為了實(shí)現(xiàn)土壤水資源可持續(xù)利用，需要依據(jù)土壤水分狀況，調(diào)控檸條生長(zhǎng)與土壤水關(guān)系［2］，因此監(jiān)測(cè)和預(yù)報(bào)土壤水分是防旱抗旱、調(diào)控檸條生長(zhǎng)與土壤水關(guān)系的基礎(chǔ)。

時(shí)間序列分析主要是采用參數(shù)模型對(duì)所觀測(cè)到的有序隨機(jī)數(shù)據(jù)進(jìn)行分析和處理的一種方法［3］，這種方法不僅可以分析因果關(guān)系，而且能根據(jù)時(shí)間序列過(guò)去的變化規(guī)律推斷以后的趨勢(shì)［4］。目前時(shí)間序列分析廣泛應(yīng)用于工農(nóng)業(yè)生產(chǎn)、自然科學(xué)、社會(huì)科學(xué)、醫(yī)學(xué)、金融及各種工程領(lǐng)域［4］。在土壤水分動(dòng)態(tài)模擬方面，劉洪斌等［5］采用時(shí)間序列自回歸建模方法對(duì)紫色土丘陵旱坡地土壤水分動(dòng)態(tài)進(jìn)行了模擬和預(yù)測(cè)，取得了較好的成效;楊紹輝等［3］應(yīng)用ARIMA 時(shí)間序列模型對(duì)北京市大興區(qū)蘆城的土壤墑情建立模型并預(yù)測(cè)，認(rèn)為ARIMA時(shí)間序列模型能很好地預(yù)測(cè)土壤墑情的變化趨勢(shì)，有實(shí)際的應(yīng)用價(jià)值;張和喜等［6］在干旱地區(qū)土壤墑情研究中用到了時(shí)間序列自回歸模型，認(rèn)為該模型可以為干旱研究及治理提供依據(jù)。本文在研究黃土丘陵半干旱地區(qū)檸條林地土壤含水量動(dòng)態(tài)變化方面運(yùn)用時(shí)間序列自回歸模型，試圖建立合理的模型并預(yù)測(cè)土壤含水量，為黃土丘陵區(qū)植物與土壤水分關(guān)系調(diào)控和植被建設(shè)提供技術(shù)支撐。

1 時(shí)間序列分析基本原理

1．1 自回歸模型

自回歸模型就是用自身做回歸變量，q 階自回歸過(guò)程{Xt}滿(mǎn)足以下方程［7］:

式中:μ 為時(shí)間序列的均值;φ1，φ2，…，φq為系數(shù);時(shí)間序列{Xt}的當(dāng)期值是自身最近q 階滯后項(xiàng)和et的線(xiàn)性組合，其中et包括了序列在t 期無(wú)法用過(guò)去值來(lái)解釋的所有新信息。

1．2 相關(guān)分析

在時(shí)間序列分析中常用自相關(guān)系數(shù)反映不同時(shí)期觀測(cè)值之間的相互關(guān)系，t 時(shí)刻的觀測(cè)值與滯后k 時(shí)刻的觀測(cè)值之間的相關(guān)程度稱(chēng)為時(shí)間延遲k 的自相關(guān)系數(shù)，記為rk，即

一般情況下計(jì)算n/5 個(gè)自相關(guān)系數(shù)［8］，即k =1，2，…，［n/5］。

1．3 定階及模型確定

模型定階的準(zhǔn)則主要有自相關(guān)函數(shù)、偏自相關(guān)函數(shù)、AIC 和BIC［7］。本文選用研究最多的赤池信息量準(zhǔn)則(Akaike’s Information Criterion，AIC)確定自回歸模型的階數(shù)，計(jì)算公式為

式中:L 為似然函數(shù);m 為參數(shù)的數(shù)量，如果模型包含截距項(xiàng)或常數(shù)項(xiàng)，則m=q+1，否則m=q。

根據(jù)土壤含水量時(shí)間序列的自相關(guān)分析，確定模型的可能最高階，用最大似然估計(jì)法估計(jì)各階模型的參數(shù)，取AIC 最小值對(duì)應(yīng)的階數(shù)作為模型的最佳階數(shù)，最后確定AR 模型。

1．4 模型診斷

如果模型被正確識(shí)別，并且參數(shù)估計(jì)充分接近真值，那么擬合的殘差應(yīng)該是獨(dú)立同分布的具有零均值和相同標(biāo)準(zhǔn)差的白噪聲序列。擬合模型的適合性檢驗(yàn)就是檢驗(yàn)殘差是否是白噪聲序列，文中采用卡方檢驗(yàn)法驗(yàn)證白噪聲序列。

設(shè){εt}為選定AR 模型估計(jì)出的殘差序列，由式(2)得到殘差序列樣本自相關(guān)函數(shù)為

Ljung-Box 統(tǒng)計(jì)量［7］計(jì)算式為

2 研究區(qū)概況及土壤水分測(cè)定

2．1 試驗(yàn)區(qū)概況

試驗(yàn)區(qū)位于黃土丘陵半干旱區(qū)的上黃生態(tài)試驗(yàn)站(寧夏固原)，該站地理位置為北緯35°59'—36°02'、東經(jīng)106°26'—106°30'。試驗(yàn)地坡度為0 ～15°，海拔約為1 650 m。該地區(qū)年內(nèi)降水量分配不均，主要集中在6—9月，占年降水量的70%以上，年際降水量變化在634．7(1984年)～259．9 mm(1991年)之間，無(wú)霜期152 d，土壤為黃綿土，植被類(lèi)型為灌叢草原。試驗(yàn)區(qū)主要植物有長(zhǎng)芒草(Stipa bungeana Trin．)、阿爾泰狗娃花(Heteropappus altaicus)、茭蒿(Artemisia giraldii Pamp．)、百里香(Thymus mongolicus Ronn．)等。

2．2 土壤水分測(cè)定

試驗(yàn)小區(qū)為5 個(gè)立地條件基本相同的標(biāo)準(zhǔn)徑流小區(qū)，每個(gè)小區(qū)長(zhǎng)20 m、寬5 m，檸條林為2002年播種，各小區(qū)播種量分別為2．0、1．5、1．0、0．5、0 kg(對(duì)照)。在每個(gè)徑流小區(qū)中部安置兩個(gè)相距2 m 的PVD 套管，管長(zhǎng)8 m，用CNC503A(DR)型智能中子水分儀測(cè)定土壤水分。測(cè)定前對(duì)中子儀進(jìn)行標(biāo)定，標(biāo)定方程為y =55．76x+1．89，式中y 為容積含水量，x 為中子儀讀數(shù)［10］。土壤水分測(cè)定深度為5—780 cm，每20 cm 測(cè)定并記錄一次數(shù)據(jù)。用環(huán)刀取樣，測(cè)定5 和20 cm 深處土壤水分，對(duì)表層土壤水分資料進(jìn)行驗(yàn)證和校正。除表層5 cm 測(cè)點(diǎn)數(shù)值代表0—10 cm 平均值外，其他測(cè)點(diǎn)土壤水分值代表測(cè)點(diǎn)高度±10 cm 范圍內(nèi)土壤水分［10］。受當(dāng)?shù)貧夂蛴绊懀瑱帡l4月中旬開(kāi)始生長(zhǎng)，10月份之后停止生長(zhǎng)，直到第二年4月中旬繼續(xù)生長(zhǎng)，因此土壤含水量測(cè)定時(shí)間為:每年4月中旬到10月底每隔15 d 測(cè)定一次，11月到次年3月每月測(cè)定一次。本文選取2011年4月至2013年2月的土壤水分資料作為研究對(duì)象，對(duì)檸條林地土壤含水量動(dòng)態(tài)變化進(jìn)行分析及預(yù)測(cè)。

3 數(shù)據(jù)處理及結(jié)果分析

3．1 數(shù)據(jù)預(yù)處理

運(yùn)用MATLAB 及R 軟件處理數(shù)據(jù)。

3．1．1 數(shù)據(jù)標(biāo)準(zhǔn)化

通過(guò)計(jì)算各層土壤含水量的變異系數(shù)，發(fā)現(xiàn)各土層土壤變異系數(shù)存在以下規(guī)律:5、20、40 cm 土層土壤水分的變異系數(shù)范圍為20% ～30%，屬于比較劇烈的變異;60—120 cm 土層土壤水分的變異系數(shù)在10% ～20%之間，屬于中等程度的變異;而140 cm 以下土層土壤水分的變異系數(shù)均小于10%，屬于弱變異。這說(shuō)明表層的土壤水分值變化較快，不穩(wěn)定，而深層土壤水分值變化緩慢，比較穩(wěn)定?？紤]到各變異程度的土層土壤含水量變化，本研究選取20、80、160 cm 深處土層的土壤水分?jǐn)?shù)據(jù)為例進(jìn)行分析，建立時(shí)間序列模型。

為了滿(mǎn)足數(shù)據(jù)分析的需要，同時(shí)提高運(yùn)算精度，首先應(yīng)對(duì)原始數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化處理。標(biāo)準(zhǔn)化公式為

式中:Yt為標(biāo)準(zhǔn)化后的時(shí)間序列;Xt為原始時(shí)間序列;X 是時(shí)間序列{Xt}的均值;σ 是時(shí)間序列{Xt}的標(biāo)準(zhǔn)差。以20 cm 深處土層土壤含水量序列為例，標(biāo)準(zhǔn)化后如圖1。

3．1．2 平穩(wěn)性

時(shí)間序列模型的分析都是建立在序列平穩(wěn)的基礎(chǔ)上的，一個(gè)平穩(wěn)的隨機(jī)過(guò)程有以下要求:均值不隨時(shí)間變化;自相關(guān)系數(shù)只與時(shí)間間隔有關(guān)，而與所處的時(shí)間無(wú)關(guān)。實(shí)際上，自然界大多數(shù)時(shí)間序列都是不平穩(wěn)的［11］，由圖1 可以看出，20 cm 深處土層土壤含水量序列是非平穩(wěn)的。本研究對(duì)各土層含水量序列進(jìn)行了Dickey-Fuller 單位根檢驗(yàn)，結(jié)果見(jiàn)表1。20、80 和160 cm 深處土層的含水量序列經(jīng)Dickey-Fuller 檢驗(yàn)后p值依次為0．504 9、0．399 5、0．551 1，均大于0．1，序列存在單位根，所以原始序列為不平穩(wěn)的。通常經(jīng)過(guò)差分使得序列變得平穩(wěn)，差分公式為Wt=Yt-Yt-1。經(jīng)過(guò)一階差分之后，檢驗(yàn)得到的p 值分別為0．035 8、0．014 3、0．051 0，均小于0．1，據(jù)此可以認(rèn)為經(jīng)過(guò)一階差分之后的時(shí)間序列為平穩(wěn)序列。

圖1 標(biāo)準(zhǔn)化后的20 cm 深土層土壤含水量序列

表1 各土層含水量時(shí)間序列的平穩(wěn)性檢驗(yàn)

以20 cm 深處土層土壤含水量為例，差分后序列如圖2，可以看出差分后的序列沒(méi)有明顯的時(shí)間趨勢(shì);有些時(shí)間點(diǎn)上序列波動(dòng)較劇烈，這是由于20 cm 土層土壤含水量受降雨量影響較大，Dickey-Fuller 檢驗(yàn)結(jié)果認(rèn)為差分后的20 cm 土層土壤含水量序列為平穩(wěn)序列，符合時(shí)間序列的建模要求。

圖2 差分后的20 cm 土層土壤含水量序列

3．2 建立模型

3．2．1 模型的定階及參數(shù)估計(jì)

應(yīng)用2011年4月到2012年11月的39 組土壤含水量序列數(shù)據(jù)進(jìn)行自相關(guān)分析，分析39/5≈8 個(gè)自相關(guān)系數(shù)，結(jié)果見(jiàn)表2。rk為滯后階數(shù)k 時(shí)的自相關(guān)系數(shù)，當(dāng)k 值取4 的時(shí)候，表2 中20、80、160 cm 三個(gè)土層的rk值都比較接近于0，故可以將模型參數(shù)估計(jì)時(shí)的最高階數(shù)定為4。

表2 土壤含水量時(shí)間序列自相關(guān)系數(shù)

確定階數(shù)后，運(yùn)用最大似然估計(jì)法估計(jì)模型的參數(shù)，并用式(4)計(jì)算各階次的AIC 值，結(jié)果見(jiàn)表3。選取AIC 最小值對(duì)應(yīng)的階次為模型的滯后階數(shù)，可以看出20 cm 土層AIC 最小值為113．69，說(shuō)明20 cm 土層含水量序列的模型為AR(1);80 cm 土層AIC 最小值為112．25，確定模型為AR(1);160 cm 土層AIC 最小值為77．35，確定模型為AR(2)。各層土壤含水量擬合模型見(jiàn)表4，Wt為差分后的土壤含水量時(shí)間序列。

表3 各土層含水量時(shí)間序列模型的參數(shù)估計(jì)結(jié)果

表4 各層土壤含水量時(shí)間序列模型的擬合結(jié)果

3．2．2 模型診斷

由式(1)可知AR 模型的殘差為

結(jié)合式(6)對(duì)20、80、160 cm 三個(gè)土層土壤含水量序列所對(duì)應(yīng)的模型殘差序列進(jìn)行白噪聲檢驗(yàn)，評(píng)價(jià)模型的擬合優(yōu)度，模型診斷結(jié)果見(jiàn)表5。

表5 各層土壤含水量序列模型的診斷結(jié)果

從表5 可以看出，這3 個(gè)模型都通過(guò)了檢驗(yàn)，可以認(rèn)為構(gòu)建時(shí)間序列自回歸模型預(yù)測(cè)黃土丘陵區(qū)半干旱地區(qū)檸條林地土壤含水量是合理的。

3．3 模型驗(yàn)證

一個(gè)模型的擬合程度較好不僅在于模型診斷結(jié)果，還在于模型的實(shí)際驗(yàn)證。本研究選取2011年4月至2012年11月的39 組數(shù)據(jù)進(jìn)行模型擬合，經(jīng)模型診斷認(rèn)為擬合較好，為了進(jìn)一步檢驗(yàn)?zāi)Ｐ偷膶?shí)用性，我們運(yùn)用建立的模型預(yù)測(cè)2012年12月至2013年2月的6組土壤含水量數(shù)據(jù)。用Yt-Yt-1替換表4 中的Wt，得到預(yù)測(cè)模型，將式(6)變?yōu)閄t=σYt+X 得到原始時(shí)間序列{Xt}，20、80、160 cm 土層土壤含水量實(shí)測(cè)值與預(yù)測(cè)值及相對(duì)誤差見(jiàn)表6。

由表6 可以看出，未來(lái)3 個(gè)月的預(yù)測(cè)數(shù)據(jù)與實(shí)測(cè)數(shù)據(jù)相對(duì)誤差較小，20 cm 土層土壤含水量預(yù)測(cè)值與實(shí)測(cè)值的相對(duì)誤差變化范圍為-5．76% ～5．92%，80 cm 土層實(shí)測(cè)值與預(yù)測(cè)值的相對(duì)誤差變化范圍為0．05% ～5．28%，160 cm 土層預(yù)測(cè)值與實(shí)測(cè)值之間的相對(duì)誤差范圍為-5．71% ～2．66%，均小于10%，所以認(rèn)為本研究建立的自回歸時(shí)間序列模型能夠有效地預(yù)測(cè)土壤含水量，對(duì)不同變異程度的土層土壤含水量預(yù)測(cè)效果也較好，可為黃土丘陵區(qū)土壤水分與植物關(guān)系調(diào)控和植被建設(shè)提供技術(shù)支撐。

表6 各層土壤含水量序列模型的驗(yàn)證結(jié)果 %

4 結(jié) 論

文中將3 個(gè)土層土壤含水量序列進(jìn)行標(biāo)準(zhǔn)化及差分處理后，用Dickey-Fuller 單位根檢驗(yàn)得出差分后的土壤含水量時(shí)間序列為平穩(wěn)序列，符合時(shí)間序列的建模要求。對(duì)于時(shí)間序列自回歸模型的階數(shù)問(wèn)題，本研究選用AIC 準(zhǔn)則確定了3 個(gè)土層的階數(shù)并建立模型，由卡方檢驗(yàn)結(jié)果認(rèn)為模型擬合得較好。最后用2012年12月至2013年2月的數(shù)據(jù)驗(yàn)證模型，3 個(gè)土層預(yù)測(cè)值與實(shí)測(cè)值的相對(duì)誤差均小于10%，說(shuō)明該模型對(duì)各變異系數(shù)的土層土壤含水量預(yù)測(cè)精度高，更進(jìn)一步說(shuō)明了時(shí)間序列自回歸模型能夠較好地預(yù)測(cè)黃土丘陵半干旱地區(qū)檸條林地土壤含水量，為該地區(qū)植物與土壤水關(guān)系的調(diào)控及植被建設(shè)提供技術(shù)支撐。

［1］尚松浩．土壤水分模擬與墑情預(yù)報(bào)模型研究進(jìn)展［J］．沈陽(yáng)農(nóng)業(yè)大學(xué)學(xué)報(bào)，2004，35(5-6):455-458．

［2］郭忠升，李耀林．植物生長(zhǎng)與土壤水關(guān)系調(diào)控起始期［J］．生態(tài)學(xué)報(bào)，2009，29(10):5721-5729．

［3］楊紹輝，王一鳴，郭正琴，等．ARIMA 模型預(yù)測(cè)土壤墑情研究［J］．干旱地區(qū)農(nóng)業(yè)研究，2006，24(2):114-118．

［4］祖元?jiǎng)?，趙則海，于景華，等．非線(xiàn)性生態(tài)模型［M］．北京:科學(xué)出版社，2004．

［5］劉洪斌，武偉，魏朝富，等．AR 模型在土壤水分動(dòng)態(tài)模擬中的應(yīng)用［J］．山地學(xué)報(bào)，2003，22(1):121-125．

［6］張和喜，楊靜，方小宇，等．時(shí)間序列分析在土壤墑情預(yù)測(cè)中的應(yīng)用研究［J］．水土保持研究，2008，15(4):82-84．

［7］Cryer Jonathan D，Chan Kung-Sik．時(shí)間序列分析及應(yīng)用:R語(yǔ)言［M］．第2 版．潘紅宇，譯．北京:機(jī)械工業(yè)出版社，2011:1．

［8］秦華光，李家才，穆丹，等．時(shí)間序列自回歸模型預(yù)測(cè)茶園小綠葉蟬種群動(dòng)態(tài)的探討［J］．安徽農(nóng)業(yè)大學(xué)學(xué)報(bào)，2008，35(4):564-570．

［9］于寧莉，易冬云，涂先勤．時(shí)間序列中自相關(guān)與偏自相關(guān)函數(shù)分析［J］．?dāng)?shù)學(xué)理論與應(yīng)用，2007，27(1):54-57．

［10］郭忠升，邵明安．人工檸條林地土壤水分補(bǔ)給和消耗動(dòng)態(tài)變化規(guī)律［J］．水土保持學(xué)報(bào)，2007，21(2):119-123．

［11］李慶雷，馬楠，付遵濤．時(shí)間序列非平穩(wěn)檢測(cè)方法的對(duì)比分析［J］．北京大學(xué)學(xué)報(bào):自然科學(xué)版，2013，49(2):252-260．