劉潭秋，王巧玲

?

基于ARIMA模型的湘江流域DO和 NH4+–N含量貝葉斯預(yù)測(cè)

劉潭秋1，王巧玲2

(1.長(zhǎng)沙理工大學(xué)經(jīng)濟(jì)與管理學(xué)院，湖南長(zhǎng)沙 410114；2.長(zhǎng)沙環(huán)境保護(hù)職業(yè)技術(shù)學(xué)院，湖南長(zhǎng)沙 410004)

為實(shí)時(shí)把控湘江流域水質(zhì)的變化趨勢(shì)，采用污染比較嚴(yán)重的湘江流域長(zhǎng)沙段和益陽段水質(zhì)指標(biāo)溶解氧(DO)和氨氮(NH4+–N)含量的監(jiān)測(cè)數(shù)據(jù)，用貝葉斯方法推斷經(jīng)典的ARIMA時(shí)間序列模型，并用馬爾可夫鏈蒙特卡羅(MCMC)模擬方法對(duì)DO和NH4+–N含量進(jìn)行貝葉斯預(yù)測(cè)。結(jié)果表明，該模型的貝葉斯預(yù)測(cè)能實(shí)現(xiàn)對(duì)湘江流域長(zhǎng)沙段和益陽段水質(zhì)指標(biāo)DO和NH4+–N含量的精確點(diǎn)預(yù)測(cè)、區(qū)間預(yù)測(cè)和概率預(yù)測(cè)。

湘江流域；溶解氧(DO)；氨氮(NH4+–N)含量；貝葉斯預(yù)測(cè)；馬爾可夫鏈蒙特卡羅(MCMC)模擬方法

湘江是湖南省重要的飲用水源地，其流域內(nèi)集中了全省60%的人口，也承載了60%以上的污染，其中以重金屬污染最為嚴(yán)重，非點(diǎn)源污染，特別是因現(xiàn)代農(nóng)業(yè)生產(chǎn)(包括種植業(yè)和畜禽養(yǎng)殖業(yè))引起的耗氧性有機(jī)物污染和氮、磷等植物營(yíng)養(yǎng)物質(zhì)濃度增高所導(dǎo)致的“富營(yíng)養(yǎng)化”現(xiàn)象亦不容樂觀[1–3]。地方政府除了限制和規(guī)范沿岸的工業(yè)、生活污水排放以及農(nóng)業(yè)耕種、養(yǎng)殖業(yè)的發(fā)展外，還需實(shí)時(shí)把控流域水質(zhì)的變化趨勢(shì)，對(duì)流域相關(guān)水質(zhì)指標(biāo)進(jìn)行實(shí)時(shí)、精確的預(yù)測(cè)，以便采取措施防患于未然，保證沿岸人們的生產(chǎn)、生活用水安全。

溶解氧(DO)是反映水體質(zhì)量的重要指標(biāo)之一，能夠反映出水體受到耗氧性有機(jī)物污染的程度。氨氮(NH4+–N)含量是反映水體富營(yíng)養(yǎng)化的重要指標(biāo)之一。傳統(tǒng)水質(zhì)預(yù)測(cè)模型的數(shù)學(xué)表達(dá)式一般都較復(fù)雜，計(jì)算成本高且耗時(shí)，很難做到實(shí)時(shí)預(yù)測(cè)?；跀?shù)據(jù)或數(shù)據(jù)驅(qū)動(dòng)類的模型越來越受到重視，其中，時(shí)間序列模型因建模步驟簡(jiǎn)單且預(yù)測(cè)精度較高而備受關(guān)注[4–5]。該類模型通常使用Box and Jenkins提出的經(jīng)典推斷方法[6]，其模型參數(shù)被設(shè)定為常數(shù)且預(yù)測(cè)結(jié)果為點(diǎn)預(yù)測(cè)。水環(huán)境是一個(gè)充滿不確定性的復(fù)雜系統(tǒng)，僅僅用模型中所包含的隨機(jī)殘差項(xiàng)來處理這種不確定性是不充分的。水質(zhì)管理中概率預(yù)測(cè)比點(diǎn)預(yù)測(cè)更有意義。目前已有貝葉斯方法應(yīng)用于水質(zhì)預(yù)測(cè)[7–8]的報(bào)道，但這些研究很少與時(shí)間序列模型相結(jié)合，因此，筆者采用貝葉斯方法來推斷時(shí)間序列模型，并對(duì)湘江流域DO和NH4+–N含量進(jìn)行貝葉斯預(yù)測(cè)。

1　材料與方法

1.1　研究區(qū)概況

湘江為長(zhǎng)江中游南岸的重要支流，干流全長(zhǎng)856 km，流域面積9.46萬km2，沿途匯入大小支流1 300多條。湘江水域集飲用水、灌溉水、漁業(yè)用水、工業(yè)用水提供和航運(yùn)、納污等多種功能于一體，為區(qū)域內(nèi)居民生活及工農(nóng)業(yè)生產(chǎn)提供重要保障。在20世紀(jì)70年代，湘江的整體水質(zhì)仍為Ⅱ類或Ⅲ類，有“綠色湘江”之稱。近幾十年來，伴隨著沿岸經(jīng)濟(jì)的迅猛發(fā)展，工業(yè)廢水、城市生活污水和農(nóng)業(yè)面源污水排放急劇增加，湘江水環(huán)境污染形勢(shì)嚴(yán)峻。目前，湘江流域設(shè)有長(zhǎng)沙新港斷面、益陽萬家嘴斷面、岳陽城陵磯斷面、常德坡頭斷面、常德沙河口斷面水質(zhì)自動(dòng)監(jiān)測(cè)站。所監(jiān)測(cè)的水質(zhì)指標(biāo)中，反映污染情況的主要指標(biāo)有DO和NH4+–N含量這2個(gè)指標(biāo)。

1.2　樣本數(shù)據(jù)

據(jù)國(guó)家環(huán)保部官網(wǎng)公布，2016年第11周至2017年第11周，長(zhǎng)沙新港斷面和益陽萬家嘴斷面DO均有多次高于Ⅲ類水質(zhì)標(biāo)準(zhǔn)限值(《GB3838—2002》限值為 5.0 mg/L)；長(zhǎng)沙新港斷面NH4+–N含量出現(xiàn)過1次高于Ⅲ類水質(zhì)標(biāo)準(zhǔn)限值(《GB3838—2002》限值為1.0 mg/L)，其余3個(gè)斷面的這2個(gè)指標(biāo)均未出現(xiàn)高于Ⅲ類水質(zhì)標(biāo)準(zhǔn)限值的情況?？梢?，相較于岳陽和常德，長(zhǎng)沙和益陽段受到的污染相對(duì)嚴(yán)重，所以，本研究中選取污染比較嚴(yán)重的湘江流域長(zhǎng)沙斷面和益陽斷面2016年第11周至2017年第11周DO和NH4+–N含量實(shí)測(cè)值進(jìn)行分析。

由圖1可見，監(jiān)測(cè)期間長(zhǎng)沙斷面的DO含量大部分情況下均低于益陽斷面的，這表明湘江長(zhǎng)沙斷面的有機(jī)物污染比湘江益陽斷面嚴(yán)重。由圖2可見，監(jiān)測(cè)期間長(zhǎng)沙斷面的NH4+–N含量大部分情況下均小于益陽斷面的，只有少數(shù)幾次大于益陽斷面的，出現(xiàn)了1次高于1.0 mg/L的情況?？梢?，湘江長(zhǎng)沙斷面的氨氮污染程度比湘江益陽斷面的輕?；谝陨戏治?，本研究中對(duì)湘江流域長(zhǎng)沙斷面和益陽斷面的DO和NH4+–N含量進(jìn)行預(yù)測(cè)試驗(yàn)，并以2016年第11周至第40周為取樣時(shí)間(因?yàn)樵诖似陂g2個(gè)斷面均頻繁出現(xiàn)了DO超標(biāo)的水質(zhì)輕度污染情況)。樣本指標(biāo)DO和NH4+–N含量各取30個(gè)監(jiān)測(cè)值。為了檢驗(yàn)?zāi)Ｐ偷念A(yù)測(cè)能力，將整個(gè)樣本(30個(gè)監(jiān)測(cè)值)劃分為2個(gè)部分：前期20個(gè)監(jiān)測(cè)值用于模型的訓(xùn)練，后10個(gè)監(jiān)測(cè)值用于檢驗(yàn)?zāi)Ｐ偷念A(yù)測(cè)能力。在數(shù)據(jù)被應(yīng)用于模型訓(xùn)練之前先進(jìn)行數(shù)據(jù)的平穩(wěn)性檢測(cè)和處理。

圖1　監(jiān)測(cè)期間長(zhǎng)沙斷面和益陽斷面的DO實(shí)測(cè)值

圖2監(jiān)測(cè)期間長(zhǎng)沙斷面和益陽斷面的NH4+–N含量實(shí)測(cè)值

Fig.2Variations of NH4+–N at Changsha section and Yiyang section during the monitoring periods

1.3　預(yù)測(cè)模型

自回歸整合移動(dòng)平均(Auto–Regressive Integrated Moving–Average，ARIMA)模型是最經(jīng)典的時(shí)間序列模型，能簡(jiǎn)單、有效地識(shí)別數(shù)據(jù)中的復(fù)雜變化模式，被廣泛應(yīng)用于各領(lǐng)域的一元時(shí)間序列預(yù)測(cè)，其一般表達(dá)式為

1.4　貝葉斯推斷與預(yù)測(cè)

采用馬爾可夫鏈蒙特卡羅(Markov Chain Monte Carlo，MCMC)模擬法，各參數(shù)逐一從其滿條件分布取樣。在執(zhí)行MCMC模擬過程中監(jiān)控和檢驗(yàn)馬爾可夫鏈的收斂，以確保樣本是從平穩(wěn)分布抽取的。設(shè)經(jīng)歷燃燒階段后抽樣迭代的次數(shù)為，繼續(xù)抽樣迭代至結(jié)束所完成的迭代總次數(shù)為。為了確保模型參數(shù)貝葉斯估計(jì)的準(zhǔn)確性，必須拋棄燃燒階段迭代獲得的預(yù)測(cè)模擬值，且–必須是一個(gè)足夠大的正整數(shù)值，因?yàn)槔碚撋现挥挟?dāng)?shù)螖?shù)充分大時(shí)，這個(gè)馬爾可夫鏈的收斂才可能發(fā)生[10]。完成迭代抽樣后，模型參數(shù)均值估計(jì)的計(jì)算公式如下：

2　結(jié)果與分析

2.1　預(yù)測(cè)模型的確定

采用Box and Jenkins經(jīng)典推斷法，使用樣本數(shù)據(jù)可以確定2個(gè)斷面的DO模型結(jié)構(gòu)均為ARIMA(1，0，1)模型，或簡(jiǎn)寫為ARMA(1，1)模型

而2個(gè)斷面的NH4+–N含量模型結(jié)構(gòu)均為不含常數(shù)項(xiàng)的ARMA(1，1)模型

用MCMC方法，分別對(duì)長(zhǎng)沙斷面和益陽斷面的DO和NH4+–N含量ARMA (1，1)模型進(jìn)行參數(shù)估計(jì)，其中，燃燒期長(zhǎng)度設(shè)定為1 000次，抽樣迭代次數(shù)設(shè)定為6 000次。待自相關(guān)圖和追蹤圖顯示所有參數(shù)的馬爾可夫鏈都已很好地收斂，表明模型已被訓(xùn)練好，可用來進(jìn)行實(shí)際預(yù)測(cè)(鑒于論文版面限制，模型參數(shù)的貝葉斯估計(jì)結(jié)果不在這里給出，但備索)。

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

根據(jù)等式(11)和(12)所確定的ARIMA模型結(jié)構(gòu)及模型參數(shù)的貝葉斯估計(jì)結(jié)果，按照1.4節(jié)的貝葉斯推斷步驟，采用MCMC方法模擬5 000次，對(duì)2016年第31周至第40周長(zhǎng)沙新港斷面和益陽萬家嘴斷面DO和NH4+–N含量進(jìn)行貝葉斯預(yù)測(cè)，獲得各監(jiān)測(cè)點(diǎn)在各預(yù)測(cè)期的5 000個(gè)水質(zhì)指標(biāo)預(yù)測(cè)值，構(gòu)成相應(yīng)預(yù)測(cè)模擬序列。

2.2.1點(diǎn)預(yù)測(cè)結(jié)果

1) 由表1、表2可見，各點(diǎn)預(yù)測(cè)結(jié)果與實(shí)測(cè)值之間的差異均較小。長(zhǎng)沙斷面DO實(shí)測(cè)值與預(yù)測(cè)值均值的最大差值為0.10 mg/L，發(fā)生在2016年第34周，最小差值為0.00 mg/L，發(fā)生在2016年第38周；益陽斷面最大差值為0.22 mg/L，發(fā)生在2016年第34周，最小差值為0.02 mg/L，發(fā)生在2016年第31、32周。長(zhǎng)沙斷面NH4+–N含量實(shí)測(cè)值與預(yù)測(cè)值均值之間的最大差值為0.01 mg/L，發(fā)生在2016年第33、35、36、37周，最小為0.00 mg/L，發(fā)生在2016年第32、34周；益陽斷面最大值為0.01 mg/L，發(fā)生在2016年第33、37周，最小值為0.00 mg/L，發(fā)生在其余預(yù)測(cè)時(shí)期。

表1　DO的貝葉斯預(yù)測(cè)結(jié)果

表1(續(xù))

實(shí)測(cè)值帶“*”號(hào)表示該值小于國(guó)家標(biāo)準(zhǔn)《GB3838–2002》規(guī)定的III類水質(zhì)標(biāo)準(zhǔn)限值，是IV類水質(zhì)，屬輕度污染。

表2　NH4+–N含量的貝葉斯預(yù)測(cè)結(jié)果

2) 長(zhǎng)沙斷面和益陽斷面DO的誤差百分比絕對(duì)值均值(MAPE)分別為1.07%、1.93%， NH4+–N含量的誤差百分比絕對(duì)值均值(MAPE)分別為4.12%、0.39%，可見，2 個(gè)指標(biāo)的MAPE值都相當(dāng)小，其中的最大值也只有4.12%，表明點(diǎn)預(yù)測(cè)結(jié)果精確。

3) 根據(jù)DO預(yù)測(cè)值均值，長(zhǎng)沙斷面有8次、益陽斷面有2次高于III類水質(zhì)標(biāo)準(zhǔn)限值，表現(xiàn)為輕度污染(IV級(jí)水質(zhì))。該結(jié)果與實(shí)測(cè)值相符，表明可根據(jù)點(diǎn)預(yù)測(cè)結(jié)果對(duì)實(shí)際污染的發(fā)生進(jìn)行預(yù)測(cè)。

2.2.2區(qū)間預(yù)測(cè)結(jié)果

由表1和表2中預(yù)測(cè)模擬序列第2.5百分位數(shù)和第97.5百分位數(shù)可知，預(yù)測(cè)期間2個(gè)斷面的DO和NH4+–N含量均處于這2個(gè)百分位數(shù)之間，或者說模型在預(yù)測(cè)2個(gè)斷面的DO和NH4+–N含量的95%可信區(qū)間內(nèi)都包含了真實(shí)值，并且MC誤差相當(dāng)小，該值最大也只有0.01 mg/L，表明模型的區(qū)間預(yù)測(cè)結(jié)果精確。

2.2.3概率預(yù)測(cè)結(jié)果

與預(yù)測(cè)概率69.92%和99.82對(duì)應(yīng)的時(shí)段，長(zhǎng)沙斷面和益陽斷面均實(shí)際發(fā)生了DO超標(biāo)所導(dǎo)致的水污染事件；NH4+–N含量IV類水質(zhì)發(fā)生的預(yù)測(cè)概率除益陽斷面2016年第37周為1%以外，其余時(shí)間均接近于0，也就是說，據(jù)概率預(yù)測(cè)結(jié)果，2個(gè)斷面因NH4+–N含量超標(biāo)而導(dǎo)致水污染的事件幾乎不會(huì)發(fā)生?？梢姡怕暑A(yù)測(cè)結(jié)果與實(shí)際情況完全相符。

3　結(jié)論與討論

水質(zhì)建模研究中，各種不確定性所帶來的影響日益受到重視。本研究中將貝葉斯方法應(yīng)用于ARIMA模型推斷，將模型中所有參數(shù)與變量都看成是隨機(jī)變量，很好地順應(yīng)了這一研究形勢(shì)。本研究中對(duì)水質(zhì)指標(biāo)DO和NH4+–N含量進(jìn)行點(diǎn)預(yù)測(cè)、區(qū)間預(yù)測(cè)的結(jié)果和概率預(yù)測(cè)的結(jié)果均與實(shí)際測(cè)量數(shù)據(jù)相符，表明該方法可供水質(zhì)管理部門對(duì)污染事件進(jìn)行預(yù)報(bào)預(yù)警，并為應(yīng)急處置方案的制定提供參考依據(jù)。

[1] 鄧琳，李宏卓，廖達(dá)偉，等．南淝河沉積物有機(jī)物污染及其生態(tài)風(fēng)險(xiǎn)評(píng)價(jià)[J]．環(huán)境科學(xué)與技術(shù)，2011，34(8)：179–184．DOI:10.3969/j.issn.1003–6504.2011.08.040．

[2] 楊鑫芳，李欣欣，李晉超，等．農(nóng)業(yè)土地利用中磷流失研究進(jìn)展[J]．山西農(nóng)業(yè)大學(xué)學(xué)報(bào)(自然科學(xué)版)，2013，33(3)：255–261．DOI:10.3969/j.issn.1671–8151.2013.03. 015．

[3] 趙英，崔福義，郭亮，等．基于BP神經(jīng)網(wǎng)絡(luò)的天津于橋水庫(kù)COD Mn預(yù)測(cè)研究[J]．南京理工大學(xué)學(xué)報(bào)(自然科學(xué)版)，2008，32(3)：376–380．DOI:10.3969/j.issn. 1005–9830.2008.03.026．

[4] WEST D，DELLANA S．An empirical analysis of neural network memory structures for basin water quality forecasting[J]．International Journal of Forecasting，2011，27(3)：777–803．DOI:10.1016/j.ijforecast.2010.09. 003．

[5] EMAMGHOLIZADEH S，KASHI H，MAROFPOOR I，et al．Prediction of water quality parameters of Karoon River (Iran) by artificial intelligence–based models[J]. International Journal of Environmental Science and Technology，2013，11(3)：645–656．DOI:10.1007/ s13762–013–0378–x．

[6] BOX G E P，JENKINS G M．Time Series Analysis：Forecasting and Control[M]．San Francisco：Holden–Day Press，1976．

[7] YANG J，REICHERT P，ABBASPOUR K C，et al. Hydrological modelling of the Chaohe Basin in China： Statistical model formulation and Bayesian inference[J]. Journal of Hydrology，2007，340(3)： 167–182．DOI:10. 1016/j.jhydrol.2007.04.006．

[8] WELLEN Christopher，ARHONDITSIS George B，LONG Tanya，et al．Quantifying the uncertainty of nonpoint source attribution in distributed water quality models：a Bayesian assessment of SWAT’s sediment export predictions[J]．Journal of Hydrology，2014，519：3353–3368．DOI:10.1016/j.jhydrol.2014.10.007．

[9] BATEN W D．Book review：theory of probability by Harold Jeffreys[J]．Natl Math Mag，1940(3)：159．

[10] BARNETT G，KOHN R，SHEATHER S．Robust bayesian estimation of autoregressive–moving–average models[J]．Journal of Time Series Analysis，1997，18(1)： 11–28．DOI:10.1111/1467–9892.00036．

責(zé)任編輯：王賽群

英文編輯：王庫(kù)

Prediction the contents of DO and NH4+–N in Xiangjiang river basin using Bayesian approach based on the ARIMA model

LIU Tanqiu1, WANG Qiaoling2

(1.School of Economics and Management, Changsha University of Science & Technology, Changsha 410114, China; 2.Changsha Environmental Protection College, Changsha 410004, China)

To master the variation of water quality in case of water security event and to take measures in advance against that in Xiangjiang river basin, the monitoring data of DO and NH4+–N in Changsha section and Yiyang section, which are two serious pollution river sections in the basin, were adopted for predicting their contents through ARIMA model which infers a classical time series model using Bayesian approach, the model parameters and prediction results were simulated by employing Markov Chain Monte Carlo (MCMC) method. The results showed that Bayesian approach in the model could accurately predict contents of DO and NH4+–N at section level, interval level, and probability level in the two selected sections.

Xiangjiang river basin; DO; NH4+–N; Bayesian prediction; Markov Chain Monte Carlo (MCMC) method

S273.2；X52

A

1007-1032(2017)05-0575-06

2017–06–01

2017–09–12

全國(guó)統(tǒng)計(jì)科學(xué)研究計(jì)劃項(xiàng)目(2013LY027)

劉潭秋(1971—)，女，四川成都人，博士，主要從事時(shí)間序列計(jì)量經(jīng)濟(jì)學(xué)模型理論與應(yīng)用研究，ltq870530569@yahoo.com

投稿網(wǎng)址：http://xb.hunau.edu.cn