基于一種時間序列模型的河流重金屬污染濃度預測研究

摘要：水環境是一個充滿不確定性的復雜巨系統，傳統水質模型很難體現重金屬污染物在河流中遷移的隨機性，因此經典的時間序列模型——ARIMA模型被應用于河流重金屬污染濃度的預測。實例分析證實，通過采用將獲得的最新數據不斷地添加到用于模型設定的樣本中，并再此基礎上獲得最近向前一個時期預測值的動態預測方法，ARIMA模型能夠獲得很好的預測表現，尤其是在充分考慮模型殘差統計分布特征的情況下，采用具有學生t分布的模型預測更精確。



關鍵詞：時間序列模型；河流重金屬污染；預測

中圖分類號：TP391.9 文獻標識碼：A



Forecast Study on Forecasting Pollutant Concentration of Heavymetal Contaminants in Streams



LIU Tanqiu1，SHEN Xinping2，WANG Hanhua1

（1.School of Management and Economics，Changsha University of Science Technology， Changsha410114， China； 

2. Dongting Lake Water Resources Administration Bureau of Hunan Province， Changsha410007， China）

Abstract：Traditional stream waterquality models are hardly able to describe stochastic behavior of heavymetal contaminants in water， due to stream environment influenced by various uncertainties. Therefore， a classic time series model， namely autoregressive integrated moving average （ARIMA） model， is used to predict pollutant concentration of heavymetal contaminants in streams. An empirical analysis evaluates the forecasting performance of two ARIMA models with different statistical distribution errors using a dynamic forecast approach. The results indicate that the two ARIMA models both perform very well， especially the one with student t distribution.



Key words： the time series model， heavy—metal contaminants in streams， forecasting



1引言

隨著近年來經濟建設的快速發展以及各種生產活動的擴大，各類水環境中重金屬污染日趨嚴重，河流作為飲用水的重要來源，其重金屬污染問題尤為引人矚目。由于大多數重金屬污染物在水環境中不易被微生物所降解，卻能被水生生物累積富集，并且重金屬在水體中的吸附、解吸、絡合、沉降、絮凝等化學過程及物理過程非常復雜，且其含量和濃度隨著水體的物理情況和化學條件的變化而變化，因此建立一個完整的河流重金屬水質模型非常困難。

傳統水質模型的建立通常是根據水流流速、污染物的彌散系數、污染物的降解系數與污染物濃度的關系，構造一個河段的水質模型，然后通過對該河段的長期監控、實驗獲得長期數據來確定模型相關參數。然而，水環境是一個充滿不確定性的復雜巨系統，傳統水質模型很難體現重金屬污染物在河流中遷移的隨機性。因此，包含隨機項的一種時間序列模型——自回歸整合移動平均（AutoRegressive Integrated MovingAverage， ARIMA）模型獲得了研究者們的關注，并被用于研究水質變量的變化規律[1—4]。由Box and Jenkins（1970）開發的ARIMA模型作為最典型的時間序列預測技術，已廣泛應用于各個領域的時序預測研究，水質管理領域亦不例外。水質變量時序觀測值往往表現出很強的序列相關性能夠被ARIMA模型很好地描述[5]。雖然因水文過程復雜且變量之間存在非線性關系和時變動態性的特征，近年來一些研究者更偏好使用神經網絡模型進行水質預測，但是也有研究者通過比較研究發現，在對污水過程變量進行向前多步預測時，神經網絡模型并不總是有更好的預測精度[4]。此外，相對于神經網絡模型，ARIMA模型有明確的數學函數關系表達式，能夠獲得水質變量與時間相關變化的信息。因此，這里我們嘗試將ARIMA模型應用于河流重金屬污染濃度的預測，事實上迄今為止很少有研究者將該模型應用于河流重金屬污染管理研究，雖然河流重金屬污染濃度無疑是一種重要的水質變量。利用河流中重金屬污染監測獲得的數據確定相應ARIMA模型，實現對污染水體的模擬分析和實時預測，行政主管部門制定事故應變決策提供科學依據。2ARIMA建模方法和技術

對河流中的某種重金屬污染物濃度變化過程{yt}建立一個ARIMA（p， d， q）模型，被表示為：

yt=φ1yt—1+…+φpyt—p+θ1εt—1+

…+θqεt—q+εt （1）

其中t表示時間序號，φi（i=1， …，p）和θi（i=1， …，q）表示模型參數，{εt}是一個均值為0且方差為σ2ε白噪音過程。使用后移算子B，上式可以寫為：

φ（B）yt=θ（B）εt（2）

其中φ（B）和θ（B）分別是階p和q的多項式：

φ（B）=1—φ1B—…—φpBp （3）

和

θ（B）=1+θ1B+…+θqBq（4）

計算技術與自動化2012年9月

第31卷第3期劉潭秋等：基于一種時間序列模型的河流重金屬污染濃度預測研究

這個模型要求yt序列是平穩的，否則需進行d階（≥0）差分處理。Box and Jenkins為ARIMA模型建立了一套完整的建模步驟[6]，包括模型的識別、估計和診斷方法，這也是該模型如此之流行的一個重要原因。

識別步驟是為了確定模型的整、自回歸項AR階的p值和移動平均項MA階的q值。首先，確定所研究的重金屬污染物濃度時序數據是否平穩，以及在不平穩情況下需經過幾階差分實現平穩，從而確定模型的“整”，即d值。這是通過單位根測試實現的，通常采用Dickey and Fuller（1979）建立的一種統計測試，即DF測試。1981年他們進一步提出ADF測試，對如下等式

△yt=μ+τt+φyt—1+γ1Δyt—1+

…+γkΔyt—k+εt（5）

進行回歸分析，其中△y表示對被研究變量的一階差分處理，參數φ的t統計量的顯著性可查Dickey and Fuller所制定臨界值表獲得[7]。如果t統計量小于表中的臨界點，那么這意味著這個原始數據是平穩的，否則是非平穩的。包含的滯后步長k是通過信息標準來選擇的。對于沒有一種明確可識別的趨勢的數據，趨勢項τt可以省略。為了測試二階差分，用Δy代替y以及Δ2y代替Δy重復這個測試過程。二階測試的回歸等式中通常不包括τt項。對于時間序列中的季節趨勢，Box and Jenkins（1976）建議采用一種季節差分。如果一個變量的季節差分Δ4yt=yt—yt—4是平穩過程而yt不是，那么稱這個變量是季節整。雖然Box and Jenkins（1976）所作的季節差分決策是基于對相關圖的肉眼觀察基礎上，但是Hylleberg et al.（1990）建立了一種統計測試，即HEGY測試，它的構造原理與DF測試相似。然后，觀察原始數據序列的自相關函數（ACF）圖和偏自相關函數（PACF）圖初步確定自回歸項AR滯后步長p值和移動平均項MA滯后步長q值，方法是分別觀察圖中的截斷（cutoff）點，即ACF（PACF）統計量表現為統計顯著的最大滯后步長對應的點。由于ACF和PACF函數都以幾何速率接近0，因此僅僅從對ACF圖和PACF圖的視覺觀察很難確定p值和q值。研究者們往往通過一些信息標準來比較一組被備選模型，其中最有名的是Akaike信息準則（Akaike Information Criterion， AIC）：

AIC=—2TL+2T×（p+q） （6）

等式中L表示最大似然估計值所對應的對數似然函數，T是樣本容量，（p+q）表示模型參數個數，如果ARIMA模型中包含有常數項，那么參數個數在原有基礎上加1。AIC公式中的第二項為懲罰函數，是通過所使用的參數個數來懲罰備選模型，其中最小AIC值對應的模型是被選定的模型。雖然這種方法會高估p和q值，但是這不會造成有偏預測，并且會在建模的最后階段——診斷測試階段得到糾正。

估計階段是對等式（1）中的參數φ=（φ1，…，φp）''、θ=（θ1，…，θq）''和σ2ε進行估計，常采用最大似然（QML）方法進行估計，其中假設εt序列服從均值為0和方差為σ2ε的正態分布，因此以yt和εt的初始值，即y0≡（y0，y—1，…，y—p+1）''和ε0≡（ε0，ε—1，…，ε—q+1）′，為條件構建ARMA（p， q）過程的條件對數似然函數為：

L（φ，θ，σ2ε|y0，ε0）=—T2log （2π）—

T2log （σ2ε）—∑Tt=1ε2t2σ2ε（7）

其中對等式（1）進行轉換有

εt=yt—φ1yt—1—…—φpyt—p—

θ1εt—1—…—θqεq—1 （8）

其中令s=0，—1，…，—p+1時ys=0，以及令s=0，—1，…，—q+1時εs=0。利用上式進行迭代可最終獲得參數QML估計值。研究證實，即使誤差εt的真實分布不是正態分布，這樣獲得的QML估計值也是一致的。

診斷階段是后驗識別模型參數個數，避免p和q值被高估或低估。如果p和q被高估。首先，檢驗估計獲得的參數p（q）是否統計顯著。如果p不顯著，那么應假設φp=0，那么自回歸階最大應為p—1而不是p。然后，重復采用估計步驟的方法對ARMA（p—1， q）模型參數進行估計。重復這兩個步驟直到AR項參數估計值統計顯著。然后，將同樣的方法應用于移動平均階。

如果p和q被低估。在這種情況下，模型不能抓住數據中所有當期依賴性。這種錯誤設定能夠從對模型殘差分析中檢驗出來，因為此時殘差不再滿足白噪音的條件。理論上，殘差能夠被一個無窮的自回歸表達式定義，即

φ（B）θ（B）yt=εt （6）

但實際上，殘差是通過將被估計系數帶入模型并計算出

t=（B）θ^（B）yt （7）

對應于樣本之前的觀察值被簡單地處理為等于0。然后，計算殘差的自相關函數

h（ε）=1T∑Tt=h+1tt—h∑Tt=12t （8）

其中h為整數。當殘差序列表現像一個白噪音，那么這個自相關函數會近似等于0。研究者還采用一種更正式的方法進行判斷，即Q統計量：

QH=T∑Hh=12h（ε） （9）

當p和q被恰當選擇且觀察值總數T很大時，這個統計量服從一個自由度為Hpq的卡方分布。如果該統計量小于臨界值，接受p和q值是正確的原假設；否則，拒絕原假設。

3實例分析

為了檢驗ARIMA模型對河流中重金屬污染濃度預測的有效性，實證分析某河流一個監測點實測獲得的重金屬污染濃度數據。限于篇幅，筆者僅以2007～2010年期間40個重金屬鎘濃度均值的歷史數據作為樣本進行研究，這里均值是指每個月固定某個時點鎘在水中濃度的均值，而不是鎘污染濃度隨時間變化的均值。這些觀測值是等時間間隔連續被記錄的，時間間隔單位是月，如圖1所示。

圖1樣本時序數據圖

按照上述ARIMA模型的建模程序，首先需檢驗數據的平穩性。應用Matlab R2012a軟件相關軟件包實施ADF單位根檢驗，在測試等式分別包含常數項，常數項+時間趨勢項和無常數項和時間趨勢項情況下，獲得t統計量（分別為—4.0968、—4.0193和—2.3567）在5%統計顯著水平下均小于對應的臨界值（分別等于—2.9402、—3.5303和—1.9478），這意味著在5%統計顯著水平ADF測試拒絕單位根原假設，即該濃度序列是平穩的，不存在單位根現象。因此，可直接用該數據序列建模而無需進行差分處理。

接下來，根據圖2所示的相關圖（correlogram）初步確定模型中AR項滯后階p和MA項滯后階q值。圖中滯后一階的相關系數和偏自相關系數值明顯較大，因此可以初步確定p和q值分別為1。更精確地確定這兩個值的方法是通過選擇最大滯后階數，這里取樣本容量的十分之一，即4，然后分別令p和q在0～4范圍內取整數值，組合成32個不同的ARIMA模型，并應用樣本數據估計每個模型和計算AIC值，其中ARIMA（1， 0， 1）模型有最小的AIC值（注：篇幅所限無詳細結果列出，備索）。這意味著對于這個樣本數據，最終確定的預測模型是ARIMA（1， 0， 1）模型。使用樣本數據對模型估計，結果見表1。圖2樣本時序數據的相關圖

表1模型參數估計結果

參數

估計值

標準誤

t統計量

p值

φ1

0.9766

0.0540

18.0774

0.0000

θ1

—0.5057

0.2241

—2.2567

0.0120

殘差測試

檢驗統計量

臨界值

p值

Q檢驗

17.3113

31.4104

0.6327

JB檢驗

26.0896

4.7481

0.0018 注：Q檢驗和JB檢驗的臨界值所對應的統計顯著性水平均為5%。

顯然，模型的AR項參數估計值φ1和MA項的參數估計值θ1都是統計顯著不為0的，這直接可由每個參數估計值對應的p值遠小于5%觀察到。我們對模型殘差進行了兩項測試，一項是LjungBox的Q檢驗，用于測試殘差序列中是否存在序列相關從而診斷模型是否穩健，另一項則是用于檢驗殘差序列是否服從正態分布的JarqueBera檢驗，簡稱JB檢驗。Q檢驗統計量（=17.3113）小于其臨界值（=31.4104），意味著在5%統計顯著性水平下接受這個ARIMA（1，0，1）模型的殘差序列不存在序列相關的原假設，這也由其對應的p值（=0.6327）>5%得到驗證。而對于JB檢驗，其檢驗統計量（=26.0896）大于對應的臨界值（=4.7418），這表明這個ARIMA（1，0，1）模型的殘差序列在5%的統計顯著水平下拒絕服從正態分布的原假設，這同樣在其p值（=0.0018） <5%得到驗證。圖3所示的模型殘差直方圖顯示該殘差序列表現出非常明顯的尖峰厚尾特征，這是一種比較典型的學生t分布的分布特征。事實上，為了計算方便研究者們大都假設ARIMA模型殘差項服從正態分布，并認為即使在小樣本情況下這一假設并不合理，但是隨著樣本容量的增大，根據大數定理，模型殘差項會漸近服從正態分布。為了了解這一觀點是否會影響模型的預測能力，這里將采用兩個ARIMA（1， 0， 1）模型作為鎘污染濃度預測模型，其中一個模型的殘差序列被假設為服從正態分布，另一個則被假設服從學生t分布。并且，根據ARIMA（1，0，1）模型參數個數，這里設定該學生t分布的自由度為3。

圖3正態分布ARIMA（1，0，1）模型殘差直方圖

為了驗證模型的預測能力，采用動態預測方法：將整個樣本數據序列劃分為兩個部分，其中前30個數據作為最初的模型識別和估計數據，而后10個數據作為預測驗證數據。因此，首先利用前30個觀察值數據來建立ARIMA模型，并在這個模型基礎上預測向前一個時期的數據。然后，將這個預測時期的真實數據加入到這個包含30個觀察值的數據序列中來，并用這31個觀測值對模型進行識別和估計，然后再在此基礎上預測向前一個時期的鎘污染濃度值。重復這個步驟，直到預測完最后一個時期的鎘污染濃度值。

圖4具有不同分布的ARIMA（1，0，1）模型預測

觀察由這兩個具有不同分布類型的ARIMA（1，0，1）模型的預測圖4可發現，具有學生t分布的ARIMA（1，0，1）模型的預測效果要比正態分布ARIMA（1，0，1）模型更好，這還能夠從常用于衡量模型的總體預測表現的統計指標值獲得驗證。三個常用指標分別是預測誤差平方均值（MSE），預測誤差絕對值均值（MAE）和預測誤差百分比絕對值均值（MAPE），當使用它們來比較不同模型的預測能力時，統計值越小則表明對應模型的預測水平越好。表2顯示t分布模型的三個預測衡量指標值都比正態分布模型的對應值要小幾乎一半，這表明總體而言前者確實有更好的預測能力。

表2模型預測結果

模型

MSE

MAE

MAPE

正態分布ARIMA（1，0，1）模型

3.3789

2.4156

0.0807

t分布ARIMA（1，0，1）模型

2.4291

1.6015

0.0542

4結論

通過經典的ARIMA時間序列模型能夠在不考慮重金屬含量和濃度隨著水體的物理情況和化學條件的變化而發生的復雜變化機理情況下，僅根據歷史數據來確定河流重金屬污染變化趨勢，并在此基礎上進行未來一個時期的預測。通過不斷將最新獲得的重金屬濃度數據添加到用于設定模型的樣本中去，以保證模型不斷更新從而能夠描述最新的重金屬污染變化趨勢，以便獲得準確預測結果。實例分析證實，ARIMA模型能夠很好地預測河流重金屬濃度變化，并且在充分考慮模型殘差的分布特征后，選擇學生t分布而不是常用的正態分布后，ARIMA模型預測水平能夠進一步得到提高。

參考文獻

[1]Nunnari， G.， Nucifora， M.， Randieri， C. The application of neural techniques to the modelling of time—series of atmospheric pollution data [J]. Ecological Modelling 111（1998）： 187—205.

[2]Chaloulakou， A.， Assimacopoulos， D.， Lekkas， T. Forecasting daily maximum ozone concentrations in the Athens basin [J]. Environmental Monitoring and Assessment， 56（1999）： 97—112.

[3]Lehmann， A. and Rode， M.： Long—term behaviour and crosscorrelations analysis of water quality parameters of the Elbe river at Magdeburg， Germany [J]. Water Research， 35（2001）： 2153–2160.

[4]Dellana， Scott A.， David West. Predictive modeling for wastewater applications： Linear and nonlinear approaches [J]. Environment Modelling Software， 24（2009）： 96—106.

[5]Yurekli， K.， Kurunc， A.， Cevik， O. Simulation of drought periods using stochastic models [J]. Turkish Journal of Engineering and Environmental Sciences 28 （3）， 2004： 181—190.

[6]Box， G.E.P. Jenkins，G.M. Time series analysis forecasting and control [C]. San Francisco： HoldenDay， 1970.

[7]Dickey， D. A. and Fuller， W. A. Distribution of the estimators for autoregressivetime series with a unit root [J]. Journal of the American Statistical Association， 74（ 1979）： 427—31.