湘江流域金屬污染的貝葉斯預(yù)測(cè)實(shí)證研究

劉潭秋　孫湘海

摘要：水環(huán)境是一個(gè)充滿不確定性的復(fù)雜巨系統(tǒng)，傳統(tǒng)水質(zhì)模型很難體現(xiàn)重金屬污染物在河流中遷移的隨機(jī)性。本文選擇ARIMA模型作為重金屬預(yù)測(cè)模型，運(yùn)用貝葉斯相關(guān)理論進(jìn)行分析、參數(shù)估計(jì)和預(yù)測(cè)，從而不僅獲得點(diǎn)預(yù)測(cè)結(jié)果，而且獲得區(qū)間估計(jì)和概率預(yù)測(cè)結(jié)果。實(shí)例分析證實(shí)，基于貝葉斯方法的ARI-MA模型能夠獲得很好的點(diǎn)預(yù)測(cè)和區(qū)間表現(xiàn)。

關(guān)鍵詞：時(shí)間序列模型；貝葉斯理論；河流重金屬污染；預(yù)測(cè)

中圖分類號(hào)：TP391.9 文獻(xiàn)標(biāo)識(shí)碼：A

1引言

水環(huán)境是一個(gè)充滿不確定性的復(fù)雜巨系統(tǒng)，根據(jù)水流流速、污染物的彌散系數(shù)、污染物的降解系數(shù)與污染物濃度的關(guān)系構(gòu)造的傳統(tǒng)水質(zhì)模型很難體現(xiàn)重金屬污染物在河流中遷移的隨機(jī)性和不確定性。

由Box and Jenkins（1970）提出的自回歸整合移動(dòng)平均（autoregressivintegratedmoving-aver-age，ARIMA）模型作為最典型的時(shí)間序列預(yù)測(cè)技術(shù)，已廣泛應(yīng)用于各個(gè)領(lǐng)域的時(shí)序預(yù)測(cè)研究，水質(zhì)管理領(lǐng)域亦不例外。此外，貝葉斯統(tǒng)計(jì)學(xué)作為處理不確定性的一種恰當(dāng)方法，自上個(gè)世紀(jì)九十年代起已廣泛被應(yīng)用于建模。目前，貝葉斯方法廣泛被應(yīng)用于簡(jiǎn)單集總式概念性水文模型的參數(shù)和不確定性估計(jì)，但很少有研究應(yīng)用這些方法去估計(jì)其他水文模型不確定性，原因是輸入數(shù)據(jù)和計(jì)算時(shí)間有限。但是水文研究中，水文建模中的不確定性已經(jīng)成為一個(gè)不可避免的主體。雖然以前只有參數(shù)不確定性被關(guān)注，但是越來越多的人意識(shí)到在建模方法中其它來源不確定性的重要性。

將ARIMA模型與貝葉斯理論相結(jié)合，這樣既能夠充分利用ARIMA模型點(diǎn)預(yù)測(cè)精度較高的優(yōu)點(diǎn)，又能運(yùn)用貝葉斯分析和推斷的優(yōu)勢(shì)，實(shí)現(xiàn)重金屬污染的概率預(yù)測(cè)。概率預(yù)測(cè)不同于點(diǎn)預(yù)測(cè)，前者不僅能夠給出具體的數(shù)值，而且能夠給出出現(xiàn)這個(gè)數(shù)值結(jié)果的可能性（概率），以及在給定的可能性（概率）下，預(yù)測(cè)結(jié)果是什么樣的數(shù)值范圍。事實(shí)上，自從1969年美國(guó)國(guó)家氣象局開始制作并發(fā)布概率降水預(yù)測(cè)至今，概率水文預(yù)測(cè)的概念在國(guó)外得到了廣泛的關(guān)注，水文學(xué)家對(duì)此進(jìn)行了大量的研究，提出了許多方法與模型，極大地豐富了水文預(yù)測(cè)不確定性研究的理論與實(shí)踐，其中貝葉斯預(yù)測(cè)是一種重要的概率預(yù)測(cè)方法。

用貝葉斯方法來分析時(shí)間序列ARIMA模型，國(guó)外己有不少學(xué)者做過相關(guān)研究，但是將其應(yīng)用于水質(zhì)預(yù)測(cè)，甚或是重金屬污染濃度預(yù)測(cè)，鮮有人涉足。因此，我們將進(jìn)行嘗試性研究，這無疑在理論上和實(shí)踐上都是一次有益的探索。

2ARIMA模型的貝葉斯推斷

2.1ARIMA模型結(jié)構(gòu)

一個(gè)ARIMA（p，o，q）模型的一般表達(dá)式為：

其中{X_t}是被研究時(shí)間序列，u是常數(shù)項(xiàng)，φ_i和θ_i是待估參數(shù)，ε_i是隨機(jī)誤差項(xiàng)，p和q分別是被研究時(shí)間序列和誤差項(xiàng)時(shí)間序列的滯后階數(shù)且均為正整數(shù)。而ARIMA（p，o，q）中“o”表示{X_t}是平穩(wěn)序列，否則需通過差分處理將其變?yōu)槠椒€(wěn)序列。

假定被研究變量有T個(gè)觀察值，記為X=（X₁，…，X_T）。模型對(duì)應(yīng)的似然函數(shù)f（X|ψ）的計(jì)算表達(dá)式為：

在對(duì)模型參數(shù)完全未知情況下，可以根據(jù)Jeffrey提出的非信息先驗(yàn)分布設(shè)定方法，將模型中參數(shù)的先驗(yàn)分布都被設(shè)置為均勻分布，其中模型方差σ²的非信息先驗(yàn)分布就是f（σ²）∝1/σ²。通常假設(shè)模型中各個(gè)參數(shù)變量是相互獨(dú)立的，因此這個(gè)聯(lián)合先驗(yàn)密度π（ψ）是各個(gè)參數(shù)先驗(yàn)密度的乘積，即π（ψ）=π（φ，θ，u，σ²）=π（φ）·π（θ）·π（u）·π（σ²）∝1/∝，其中-∞<φ，θ，u<-∞，0<σ²<∞。

根據(jù)貝葉斯定理，待估參數(shù)的后驗(yàn)密度為：

π（ψ|X）∝f（X|ψ）·π（ψ） （3）

將每個(gè)參數(shù)的均勻分布式代入上式，就形成模型參數(shù)的聯(lián)合后驗(yàn)分布。

在具體的推斷實(shí)施過程中，這里采了用Gelf andan dSmith（1990）提出的Gibbs取樣方法，逐一從與各個(gè)參數(shù)相聯(lián)系的滿條件分布取樣。具體步驟為：

6）這樣，就獲得全部參數(shù)的第1次完整抽樣，重復(fù)上述1）-5）步驟，獲得全部參數(shù)的第2次完整抽樣，然后再多次重復(fù)上述1）-5）步驟，直至獲得全部參數(shù)的第N次完整抽樣。這樣，就形成這個(gè)ARIMA模型參數(shù)的樣本點(diǎn)序列{φ₁^（j），…，φ_p^（j），θ₁^（j），…，θ_q^（j），μ^（j），σ^{2，（j）}}_j^N=1。

在實(shí)施上述抽樣步驟中，同時(shí)監(jiān)控和檢驗(yàn)馬爾可夫鏈的收斂，以確保樣本是從平穩(wěn)分布抽取。目前在實(shí)際操作中對(duì)于馬爾可夫鏈?zhǔn)諗颗c否的常用監(jiān)測(cè)方法是樣本圖形分析法、追蹤圖法、自相關(guān)圖法，以及一些正式收斂測(cè)試法，例如Geweke測(cè)試法，Gelman-Rubin測(cè)試法等。

2.2ARIMA模型的貝葉斯估計(jì)

在完成上述迭代抽樣之后，直接采用如下公式獲得模型參數(shù)的均值估計(jì)值：

其中，N表示迭代總次數(shù)，M表示馬爾可夫鏈?zhǔn)諗壳暗臉颖军c(diǎn)。為了避免初始值對(duì)最終估計(jì)結(jié)果的影響，這前M個(gè)樣本點(diǎn)必須被拋棄，而這些被拋棄的抽樣迭代部分被稱作燃燒階段。此外，為了保證模型參數(shù)貝葉斯估計(jì)的準(zhǔn)確性，N-M必須是一個(gè)足夠大的正整數(shù)值，以確保該馬爾可夫鏈?zhǔn)諗俊?

ARIMA模型的貝葉斯估計(jì)，不僅能夠通過4）式獲得模型參數(shù)的均值估計(jì)值，而且能夠獲得參數(shù)其他統(tǒng)計(jì)特征值，例如百分位數(shù)、標(biāo)準(zhǔn)差等。

2.3ARIMA模型的貝葉斯預(yù)測(cè)

在貝葉斯方法中，預(yù)測(cè)是基于所研究變量未來值向量X^F的一個(gè)概率分布的構(gòu)建，并且以過去（被觀察到的）值X向量為條件，并考慮這個(gè)預(yù)測(cè)模型參數(shù)ψ的后驗(yàn)分布。后驗(yàn)預(yù)測(cè)密度表達(dá)式是：

f（X^F|X）=∫f（X^F|X，ψ）π（ψ|X）dψ （5）

其中f（X^F|X，ψ）是條件預(yù)測(cè)密度。這個(gè)后驗(yàn)預(yù)測(cè)密度的MCMC解為：

其中，{ψ^（i）}^Gi=1，是從參數(shù)的滿條件分布抽樣獲得，G是馬爾可夫鏈的迭代次數(shù)。

這個(gè)后驗(yàn)預(yù)測(cè)密度f(wàn)（X^F|X）能夠被解釋為來自條件預(yù)測(cè)分布f（X^F|X，ψ）的一次平均，權(quán)重值是這些參數(shù)的后驗(yàn)概率。貝葉斯預(yù)測(cè)結(jié)果就是預(yù)測(cè)的可信區(qū)間，其反映了所研究現(xiàn)象的不確定性。

3湘江流域重金屬污染的貝葉斯預(yù)測(cè)實(shí)證研究

3.1樣本數(shù)據(jù)與ARIMA模型的確定

以湘江中下游某個(gè)重金屬監(jiān)測(cè)點(diǎn)獲得的一組鎘污染濃度數(shù)據(jù)為樣本進(jìn)行實(shí)證分析。這組鎘濃度數(shù)據(jù)是2007-2010年期問每個(gè)月固定某個(gè)時(shí)點(diǎn)鎘在水中濃度的均值，而不是鎘污染濃度隨時(shí)間變化的均值，共計(jì)40個(gè)樣本數(shù)據(jù)。這些觀測(cè)值是等時(shí)間間隔連續(xù)被記錄的，時(shí)間間隔單位是月。選擇前35個(gè)樣本數(shù)據(jù)對(duì)神經(jīng)網(wǎng)絡(luò)模型進(jìn)行訓(xùn)練，后5個(gè)樣本數(shù)據(jù)則用于與預(yù)測(cè)結(jié)果進(jìn)行比較，從而驗(yàn)證模型的預(yù)測(cè)能力。本次研究所使用的數(shù)據(jù)與筆者前期開展的使用ARIMA模型來預(yù)測(cè)重金屬污染濃度研究所使用的數(shù)據(jù)相同，這樣選擇樣板數(shù)據(jù)也是便于將本次研究與之前研究工作進(jìn)行對(duì)比。

根據(jù)前期研究可知，ARIMA模型的具體表達(dá)式為一個(gè)ARIMA（1，0，1）模型：

X_t=u+φ₁X_t-1+θ₁ε_t-1+ε_t （7）

其中，X_t表示第t時(shí)刻的鎘污染濃度，u是常數(shù)項(xiàng)，ε_t是模型殘差項(xiàng)且服從均值為0且方差為常數(shù)的學(xué)生t分布，φ₁和θ₁是模型需要被估計(jì)的參數(shù)，這個(gè)模型被記為服從學(xué)生t分布的ARIMA（1，0，1）模型，或稱為服從學(xué)生t分布的ARMA（1，1）模型，建模過程詳見文獻(xiàn)。

3.2模型的貝葉斯分析

基于上述模型的殘差項(xiàng)分布設(shè)定，模型的似然函數(shù)為：

鑒于本次實(shí)證研究對(duì)象應(yīng)用本文所討論方法的相關(guān)研究很少，這里采用Jeffrey提出的非信息先驗(yàn)分布設(shè)定方法，在使用傳統(tǒng)ARIMA建模方法獲得參數(shù)估計(jì)結(jié)果基礎(chǔ)上，可將模型的非信息先驗(yàn)分布設(shè)定如下：

π（φ）～Uniform（0，1）

π（θ）～Uniform（-1，0） （9）

π（σ²）～Uniform（0，160）

根據(jù)貝葉斯定理，給定數(shù)據(jù)，模型參數(shù)的聯(lián)合后驗(yàn)分布正比于參數(shù)聯(lián)合先驗(yàn)分布密度與似然函數(shù)的乘積。按照上述對(duì)這個(gè)ARIMA（1，0，1）模型參數(shù)先驗(yàn)分布的設(shè)定，以及顯然這些先驗(yàn)分布密度相互獨(dú)立，因此模型參數(shù)的聯(lián)合后驗(yàn)密度就直接正比于模型的似然函數(shù)與每個(gè)參數(shù)先驗(yàn)密度的乘積，如等式（3）所示。

相應(yīng)地，參數(shù)φ的滿條件分布為：

顯然，這個(gè)滿條件分布不服從某個(gè)標(biāo)準(zhǔn)的概率統(tǒng)計(jì)分布，因此在實(shí)施σ²取樣時(shí)可采用Metropolis-Hastings算法（參見文獻(xiàn)）。

3.3貝葉斯估計(jì)結(jié)果

在對(duì)所構(gòu)建的鎘污染濃度預(yù)測(cè)模型實(shí)施完貝葉斯分析后，采用OpenBUGS軟件對(duì)模型實(shí)施貝葉斯估計(jì)。在估計(jì)過程中，燃燒期長(zhǎng)度選擇為5000次（即，M=5000），觀察自相關(guān)圖和追蹤圖可以發(fā)現(xiàn)，所有參數(shù)的馬爾可夫鏈都很好地收斂了。為了保證估計(jì)值的準(zhǔn)確性，在燃燒期之后又迭代取樣50000次（即，N=50000）所獲得的樣本點(diǎn)用于（7）式以計(jì)算每個(gè)參數(shù)的馬爾可夫積分，即參數(shù)均值估計(jì)值。具體結(jié)果見表10

3.4貝葉斯預(yù)測(cè)結(jié)果

利用已經(jīng)完成參數(shù)估計(jì)的ARIMA（1，0，1）模型來進(jìn)行預(yù)測(cè)。首先需要計(jì)算預(yù)測(cè)值的后驗(yàn)預(yù)測(cè)密度，根據(jù)等式（5），模型的后驗(yàn)預(yù)測(cè)密度為

顯然直接對(duì)上式積分是很困難的，同樣需要采用MCMC方法來獲得模型的預(yù)測(cè)結(jié)果，即貝葉斯概率預(yù)測(cè)結(jié)果，這樣就很容易預(yù)測(cè)未來湘江流域鎘污染濃度變化范圍，以及發(fā)生可能性，甚至能夠預(yù)測(cè)極端事件發(fā)生的概率。其中貝葉斯均值預(yù)測(cè)值是通過等式（4）計(jì)算得到，其中為了保證馬爾可夫鏈的收斂，G取值了一個(gè)相當(dāng)大的值，即50000。具體預(yù)測(cè)結(jié)果如下一系列表。

由表2可知：貝葉斯概率預(yù)測(cè)結(jié)果表明，在百分位5到百分位95之間所對(duì)應(yīng)的90%可信區(qū)間內(nèi)包含了真實(shí)值；貝葉斯概率預(yù)測(cè)的均值預(yù)測(cè)結(jié)果與ARIMA模型采用傳統(tǒng)建模方法所獲得的點(diǎn)預(yù)測(cè)結(jié)果非常接近真實(shí)值，但第二個(gè)預(yù)測(cè)點(diǎn)除外，這也能夠從圖1直觀地感受到。采用常用衡量預(yù)測(cè)精度的統(tǒng)計(jì)量——誤差百分比絕對(duì)值均值（MAPE）來進(jìn)行對(duì)比，我們發(fā)現(xiàn)貝葉斯均值預(yù)測(cè)和ARIMA模型的傳統(tǒng)點(diǎn)預(yù)測(cè)結(jié)果所對(duì)應(yīng)的MAPE值分別為0.0542和0.0813，顯然在單個(gè)值預(yù)測(cè)上ARIMA傳統(tǒng)預(yù)測(cè)方法似乎更優(yōu)。然而，貝葉斯預(yù)測(cè)的最大優(yōu)勢(shì)在于其不斷改進(jìn)的能力，即將實(shí)踐經(jīng)驗(yàn)不斷融入先驗(yàn)信息從而不斷改進(jìn)貝葉斯預(yù)測(cè)結(jié)果，這是傳統(tǒng)ARIMA模型預(yù)測(cè)所缺乏的。鑒于本研究在應(yīng)用領(lǐng)域的開創(chuàng)性，因此缺乏相關(guān)經(jīng)驗(yàn)積累，所以本次貝葉斯預(yù)測(cè)中所使用的先驗(yàn)分布采用非信息先驗(yàn)分布，這是導(dǎo)致這一結(jié)果產(chǎn)生的原因。因此，為了進(jìn)一步改進(jìn)模型的貝葉斯預(yù)測(cè)能力，需要在以后的實(shí)踐工作中不斷總結(jié)預(yù)測(cè)結(jié)果，不斷改進(jìn)模型參數(shù)信息先驗(yàn)分布的設(shè)定，使其不斷逼近真實(shí)統(tǒng)計(jì)分布。

4結(jié)論

將貝葉斯理論和方法與重金屬時(shí)間序列預(yù)測(cè)ARIMA模型相耦合，并對(duì)該模型所對(duì)應(yīng)的貝葉斯分析、參數(shù)估計(jì)和預(yù)測(cè)理論進(jìn)行了探討，并在此基礎(chǔ)上實(shí)施湘江重金屬污染濃度貝葉斯預(yù)測(cè)。實(shí)例預(yù)測(cè)結(jié)果表明，獲得了較好的均值預(yù)測(cè)結(jié)果和區(qū)間預(yù)測(cè)，而對(duì)于現(xiàn)有預(yù)測(cè)結(jié)果的改進(jìn)有賴于模型參數(shù)先驗(yàn)分布選擇，這需要在預(yù)測(cè)實(shí)踐基礎(chǔ)上不斷改進(jìn)。