趙春蘭,李 屹,何 婷,武 剛，王 兵

(1.西南石油大學理學院,四川 成都 610500；2.西南石油大學人工智能研究院，四川 成都 610500;3.中國石油天然氣股份有限公司大港油田分公司,天津 300280,4.西南石油大學計算機科學學院，四川 成都 610500)

1 引言

岷江是長江上游的重要支流之一，也是成都平原最重要的水資源，其水質與人民的生活緊密相關。由于工業(yè)廢水和生活污水的排放，以及超量化肥和農藥的使用，岷江的污染也日益嚴重。因此，對未來一段時間水質進行較為準確的動態(tài)預測是該流域水資源管理的當務之急。

國內外關于水質預測有許多模型，可以分為機理模型和非機理模型2類。機理模型是依據(jù)物質能量守恒和質量守恒原理，通過流體力學中的連續(xù)方程、運動方程和能量方程推導得出[1]。1925年美國工程師Streeter和Phelps提出的氧平衡S-P模型開創(chuàng)了水質模擬模型的先河，隨后又出現(xiàn)了其他一系列水質模擬模型，如QUAL、WASP、BASINS和MIKE模型等[2]。機理模型考慮了生物化學等因素，預測精度較高，但其模型往往較為復雜，且所需采集的數(shù)據(jù)非常多。而非機理模型是指利用統(tǒng)計方法或人工智能方法建立的水質預測數(shù)學模型，只需對水質指標數(shù)據(jù)進行研究分析。在過去幾十年中，時間序列預測法、回歸分析預測法、灰色模型法、人工神經網絡預測法和支持向量機等非機理模型已成功地應用于水質預測。如Sarkar等[3]利用人工神經網絡模型對馬圖拉下游的溶解氧濃度進行了預測，取得了較好的預測效果；Kisi等[4]考慮水質數(shù)據(jù)的非線性特性，分別利用最小二乘支持向量機和多元自適應回歸樣條2種模型對印度Yamuna河的水質進行了預測，認為非線性預測模型有較好的預測性能；Avila等[5]使用回歸分析、分類樹和馬爾科夫鏈等統(tǒng)計模型對游樂場水質中的大腸桿菌進行了預測研究，發(fā)現(xiàn)貝葉斯模型在處理缺失數(shù)據(jù)和異常值方面很靈活，并且允許實時的連續(xù)更新，是一種有參考價值的預測方法。

由于單獨的預測模型大都存在一定的短板，如回歸分析對數(shù)據(jù)量需求較大，泛化能力不足；BP(Back Propagation)神經網絡存在收斂速度慢、容易陷入局部最優(yōu)的缺點。目前國內外學者為了提高水質預測精度，通常將多種方法聯(lián)合使用，以此來克服單獨使用一種方法的缺點[6]。如Faruk[7]采用ARIMA(AutoRegressive Integrated Moving Average)模型與人工神經網絡模型相結合的方法，提取水質數(shù)據(jù)的線性和非線性趨勢，對BuyukMenderes河的水溫、硼和溶解氧進行了預測，取得了比單一模型更高的預測精度；Liao等[8]將神經網絡與決策樹模型相結合，建立了一種改進的決策樹學習IDTL(Improved Decision Tree Learning)模型，對中國最大的淡水湖之一巢湖的水質數(shù)據(jù)進行了預測，在犧牲了較小分類精度的情況下，能夠為分類提供一些明確的啟發(fā)式方法，提高了水質監(jiān)測水平和預測能力；Liu等[9]提出了一種將實值遺傳算法RGA(Real-value Genetic Algorithm)與支持向量回歸SVR(Support Vector Regression)結合的混合模型RGA-SVR，對宜興市水產養(yǎng)殖場水質數(shù)據(jù)進行了預測?；旌夏Ｐ屠脤嵵颠z傳算法優(yōu)化了SVR的參數(shù)，故RGA-SVR模型比BP神經網絡和傳統(tǒng)SVR模型更適合于水質預測。以上研究表明，混合方法的預測精度通常高于單獨方法的，故將多種方法聯(lián)合使用進行水質預測，是一種可行的提高預測精度的辦法。

考慮到水質時間序列數(shù)據(jù)一般都具有季節(jié)性、非線性和模糊性等多種復雜特性，又因為傳統(tǒng)的時序預測對數(shù)據(jù)的完整性和準確性要求較高，本文提出了一種基于動態(tài)隸屬度的模糊時間序列模型。首先，將水質時間序列處理成平穩(wěn)序列；其次，根據(jù)水質時間序列的特點以及模糊方法在處理不確定性問題和模糊問題的有效性，將數(shù)據(jù)模糊化，利用模糊C均值FCM(Fuzzy C-Means)聚類算法得到聚類中心和動態(tài)子隸屬度序列；隨后對每個子隸屬度序列建立ARIMA乘積季節(jié)模型，在提取季節(jié)效應的同時預測下一個時刻的隸屬度；最后去模糊化，得到水質預測結果。為了驗證模型的預測效果，本文收集了岷江某斷面2011年1月至2017年8月逐月監(jiān)測的水質污染指標的數(shù)據(jù)，對水質進行了未來3個月的預測分析，并與ARIMA乘積季節(jié)模型和經典模糊時間序列FTS(Fuzzy Time Series)模型進行了對比分析。

2 水質預測模型

2.1 ARIMA乘積季節(jié)模型

由于水質數(shù)據(jù)表現(xiàn)出明顯的季節(jié)性特征，故本文采用ARIMA乘積季節(jié)模型作為對照。簡單的季節(jié)模型認為序列的季節(jié)效應、長期趨勢和隨機效應之間是很容易分開的；但乘積季節(jié)模型認為序列的長期趨勢、季節(jié)效應和隨機效應之間相互影響，且它們之間的關系比較復雜[10]。

ARIMA乘積季節(jié)模型的建模思路為：對序列進行d階差分和D階以周期S為步長的差分運算；當序列仍然存在短期相關性時，使用低階ARIMA(p,d,q)提取短期相關性，p和q分別表示非季節(jié)自回歸階數(shù)和移動平均階數(shù)；當序列還存在季節(jié)效應時，則使用以周期步長為單位的ARIMA(P,D,Q)模型提取相關關系，P和Q分別表示季節(jié)自回歸階數(shù)和移動平均階數(shù)。

此時擬合的模型為2個模型的乘積，如式(1)～式(5)所示：

(1)

θ(B)=1-θ1BS-θ2BS-…-θqBS

(2)

(3)

ΘS(B)=1-Θ1BS-Θ2BS-…-ΘQBS

(4)

ΦS(B)=1-Φ1BS-Φ2BS-…-ΦPBS

(5)

ARIMA模型常用于時間序列預測，經典的ARIMA模型預測流程可分為3步：首先對原始數(shù)據(jù)進行ADF(Augmented Dickey-Fuller)平穩(wěn)性和白噪聲檢驗，若非平穩(wěn)，可進行差分處理；其次根據(jù)AIC(Akaike Information Criterion)最小準則進行模型的識別與定階；最后對ARIMA(p,d,q)模型進行參數(shù)估計與檢驗,模型的建立與預測流程如圖1所示。參數(shù)估計常用的方法有最大似然估計、無約束最小二乘法和條件最小二乘法[11]。本文采用條件最小二乘估計法對模型的參數(shù)進行估計。

Figure 1 The modeling &predition of ARIMA model圖1 ARIMA模型建模和預測流程

2.2 經典模糊時間序列預測模型

1994年，Song等[12,13]運用模糊集合、隸屬函數(shù)等概念，對時間序列數(shù)據(jù)建立了經典的模糊時間序列預測模型(FTS)，隨后該模型及其改進模型被廣泛地應用于各個領域?？紤]水質數(shù)據(jù)存在模糊性的特點，本文將水質數(shù)據(jù)模糊化后再進行預測分析。

定義1[12,14]設U為論域，且U={u1,u2,…,un} 是給定U上的一個次序分割集，定義A為U上的語義變量集(模糊集)，如式(6)所示：

(6)

或表達為向量形式，如式(7)所示：

A={fA(u1),fA(u2),…,fA(un)}

(7)

其中,fA(·)是模糊集A上的模糊隸屬度函數(shù)，且fA(·):U→[0,1];fA(ui)是ui在模糊集A上的模糊隸屬度值，且有fA(ui)∈[0,1],1≤i≤n。

定義2[12,14]設Z(t)(t=0,1,2,…)為實數(shù)域R的一個子集，Z(t)上定義了一組模糊集fi(t)，假設F(t)={f1(t),f2(t),…}，稱F(t)為Z(t)上的模糊時間序列。

這里的F(t)為語言變量，fi(t)為F(t)中可能的語言值，即定義1中的fA(ui)。因為對于不同的時刻t，F(xiàn)(t)是不同的，故F(t)為t的函數(shù)。

定義3[12,14]設R(t,t-1)為從F(t-1)到F(t)的模糊關系，且滿足F(t)=F(t-1)°R(t,t-1)，則稱F(t)是由F(t-1)通過模糊關系R(t,t-1)推導得到的，其中“°”代表合成運算，關系R(t,t-1)稱為定義在F(·)上的一階模糊關系。

定義4[12,14]設F(t)是一個模糊時間序列。令F(t-1)=Ai，F(xiàn)(t)=Aj，這2個連續(xù)的觀測值F(t-1)和F(t)之間可以用一個模糊邏輯關系來表示，記為Ai→Aj，其中Ai和Aj分別被稱為“左件”和“右件”，Ai對應第i個模糊概念。

定義5[12,14]定義關系矩陣R的運算方法如式(8)所示：

R=R(t,t-1)=∪i,jRi,j

(8)

經典的模糊時間序列建模步驟一般為：

(1)根據(jù)序列確定論域U，對U進行劃分。為方便計算，基于序列的最小值和最大值，對最小值向下減一個合適的正數(shù)，最大值向上加一個合適的正數(shù)，將得到的2個數(shù)作為區(qū)間的2個端點，該區(qū)間也即論域，隨后對論域進行等間距的模糊劃分。

(2)定義模糊集和數(shù)據(jù)模糊化。論域U被劃分成了k個子區(qū)間，相對應地就有k個模糊概念Ai。根據(jù)相應的隸屬函數(shù)，求出數(shù)據(jù)對每個模糊集的隸屬度向量，并將數(shù)據(jù)模糊化為最大隸屬度所對應的模糊集。

(3)建立模糊邏輯關系并確定模糊關系矩陣?；跀?shù)據(jù)的模糊概念及其出現(xiàn)的先后順序，建立模糊邏輯關系。根據(jù)定義5可求得模糊關系矩陣R。

(4)去模糊化后進行預測。選擇合適的去模糊化方法，將輸出去模糊化后得到的預測值。

本文參照國家的水質分類標準，將序列等分為5個模糊區(qū)間，按照I和II 類的水質為優(yōu)、III 類的水質為良、IV 類的水質為輕度污染、V 類的水質為中度污染、劣 V 類的水質為重度污染進行劃分，隨后根據(jù)上述步驟進行預測分析。

3 基于動態(tài)隸屬度的模糊時間序列模型

3.1 基本思想

時間序列ARIMA模型基礎理論完善，方法簡單易行，但模型的建立依賴于數(shù)據(jù)，要求數(shù)據(jù)完整且較為精準，而且在預測中主要考慮了時間因素，忽略了引起序列本身變化的一些不確定因素。但是，某些水質數(shù)據(jù)存在固有的不確定性問題，例如測量不準確、觀測集不完整或難以獲得測量值等[13]。模糊時間序列預測模型作為常用的不確定信息處理方法，為水質預測開辟了新路徑。但是，經典的模糊時間序列模型對實現(xiàn)步驟的要求較高，其預測步驟所應用的方法不同也會對預測結果造成很大的影響。根據(jù)國內外學者將不同模型結合以提高模型性能和預測精度的經驗，以及水質數(shù)據(jù)本身的模糊性和季節(jié)性特性，本文結合時間序列分析和模糊理論，提出了一種基于動態(tài)隸屬度的模糊時間序列預測模型。該模型以模糊時間序列預測模型為基礎，將ARIMA 模型加入其中，對隸屬度進行動態(tài)的季節(jié)性預測，簡化運算步驟的同時還提高了預測精度。

3.2 模糊C均值聚類算法

(9)

其中，N為樣本總數(shù);k為聚類數(shù)目，在本文中即k個模糊概念;m是控制模糊重疊程度的模糊劃分矩陣指數(shù)，且m>1;‖xt-ch‖表示t時刻的樣本值xt到聚類中心ch的歐幾里得距離;uth表示樣本xt對第h個聚類的隸屬度。引入拉格朗日參數(shù)λ，利用拉格朗日方法對以上問題求解，令拉格朗日函數(shù)如式(10)所示：

(10)

令?F/?uth=0,?F/?λ=0,通過一系列計算，可以得到聚類中心,如式(11)所示：

(11)

隸屬度函數(shù)如式(12)所示：

(12)

由此可以得到FCM算法的步驟為：

(1)隨機初始化聚類成員值uth；

(2)利用式(11)計算聚類中心，按最小距離的原則進行聚類；

(3)根據(jù)式(12)更新隸屬度uth的值，通常情況下，m的值取為2；

(4)計算目標函數(shù)Jm的值；

(5)重復步驟2～步驟4，直到Jm的值小于指定的最小閾值，或者達到最大迭代次數(shù)之后停止迭代，得到最后的聚類結果;

(6)確定聚類個數(shù)k后，利用上述步驟，便可以得到k個聚類中心。

3.3 模型預測步驟

模型的預測步驟主要包括：原始數(shù)據(jù)平穩(wěn)化、構建論域并劃分子論域、模糊化數(shù)據(jù)、子隸屬度序列時間序列預測、去模糊化得到預測值。圖2展示了其預測流程,具體執(zhí)行步驟如下所示：

(1)原始數(shù)據(jù)平穩(wěn)化。

判斷序列是否為平穩(wěn)序列，若非平穩(wěn)，則需要通過差分化處理得到平穩(wěn)的序列用于后續(xù)分析，以確保后續(xù)分析的有效性和準確性。

(2)構建論域并劃分子論域。

找到建模數(shù)據(jù)中的最大值DMax和最小值DMin，為方便劃分論域，選擇2個合適的正實數(shù)d1和d2，則得到論域U={DMin-d1,DMax+d2}。

利用式(10)確定聚類個數(shù)：

(13)

其中，N為樣本總數(shù);xt和xt-1分別為t時刻和t-1時刻的樣本值。

將聚類中心按從小到大排序，可得到聚類中心c1,c2,c3,…,ck,設聚類中心ch和ch+1之間的中點為bh，那么bh就為子論域的邊界值。假設子隸屬度序列為Y={Y1,Y2,…,Yk}，則它們的論域分別為：u1:[DMin-d1,b1],u2:[b1,b2],u3:[b2,b3],…,uk:[bk-1,DMax+d2]。

(3)模糊化數(shù)據(jù)。

FCM聚類算法在得到聚類中心的同時還得到一個隸屬度矩陣(uht)k×N(1≤h≤k,1≤t≤N)，如式(14)所示：

(14)

在隸屬度矩陣中，每一行代表每個子論域所對應的子隸屬度序列。假設u12是第2列中數(shù)值最大的值，就代表原始序列的第2個值隸屬于第一個論域的程度最大，矩陣中其他數(shù)值的意義依次類推，即：

Y1={u11,u21,…,uN1},

Y2={u12,u22,…,uN2},

?

Yk={u1k,u2k,…,uNk}

(15)

其中,uth(t=1,2,…,N;h=1,2,…,k)是樣本xt對區(qū)間uh的隸屬度，且uth∈[0,1]。

(4)子隸屬度序列的時間序列預測。

平穩(wěn)序列X={x1,x2,…,xN}通過FCM聚類算法中的k個隸屬度函數(shù)模糊化后，得到k個子隸屬度序列Y={Y1,Y2,…,Yk}。分別對每個子隸屬度序列進行時間序列預測，并利用ARIMA乘積季節(jié)模型提取季節(jié)趨勢,若乘積季節(jié)模型無法擬合，則選擇傳統(tǒng)的ARIMA(p,d,q)模型進行隸屬度的預測。

(5)去模糊化得到預測值。

FCM算法可得到k個聚類中心c1,c2,c3,…,ck，則利用如式(16)所示的隸屬度加權平均法得到預測值G(t)：

(16)

其中，fth表示t時刻對應的數(shù)值隸屬于第h個模糊集的程度，即隸屬度。

Figure 2 Flow chart of model prediction圖2 模型預測流程圖

4 實證研究

4.1 預測指標選取

本文對岷江某斷面2011年1月至2017年11月共83個月逐月監(jiān)測的6項指標的水質數(shù)據(jù)，按照國家規(guī)定的《地表水環(huán)境質量標準》(GB3878-2002)中相應的評價標準進行評價，利用最常用的單因子法對其逐月評價完畢之后發(fā)現(xiàn)，其中溶解氧與總磷的污染較為嚴重，于是確定再對這2個水質指標進行預測分析，以達到污染防治的目的。原始數(shù)據(jù)如表1所示，其中，前80個樣本點作為模型訓練集，后3個樣本點作為模型預測集。

Table 1 Water quality indicator data表1 水質指標數(shù)據(jù)

Figure 3 Time series fitting diagrams of dissolved oxygen sub-membership series圖3 溶解氧數(shù)據(jù)子隸屬度序列時間序列擬合圖

4.2 模型預測性能評價指標

為了更科學準確地評價預測效果，并比較不同模型的優(yōu)劣，本文選用式(17)～式(19)所示的3個指標來評價模型的預測性能。

(1)均方根誤差(RMSE)：

(17)

(2)平均絕對誤差(MAE):

(18)

(3)平均絕對百分誤差(MAPE):

(19)

4.3 實例分析

本節(jié)以溶解氧數(shù)據(jù)為例展示模型的預測過程，總磷的預測過程同理。將溶解氧數(shù)據(jù)進行一階差分得到平穩(wěn)數(shù)據(jù)后，計算聚類個數(shù)，并利用FCM聚類算法中的隸屬度函數(shù)，得到4個子隸屬度序列如表2所示。隨后對這4個子隸屬度序列Y1,Y2,Y3和Y4依次進行傳統(tǒng)的時間序列分析。

Table 2 Dynamic sub-membership sequence表2 動態(tài)子隸屬度序列

對4個子隸屬度序列按照2.1節(jié)中的流程進行預測，分別擬合如下模型：

Y1:ARIMA(3,1,0)

Y2:ARIMA(3,1,0)

Y3:ARIMA(0,1,1)×(0,1,1)12

Y4:ARIMA(1,1,0)×(1,1,0)12

可得到4個子隸屬度序列的擬合效果,如圖3所示。

通過對4個子隸屬度序列的預測，得到每個子隸屬度序列最后3個月的預測值以及去模糊化之后的值,如表3所示。表3中每一行的值表示同一個預測點隸屬于不同子論域的程度，每一列的值表示不同預測點隸屬于同一個子論域的程度。利用式(16)對每一行進行去模糊化操作，得到的值如表3最后一列所示，該值表示每一個點去模糊化后的真實預測值。將預測值還原到差分前的值，可得到原始數(shù)據(jù)后3個月的預測值。

Table 3 Predicted values of membership表3 隸屬度預測值

4.4 結果評測與分析

為比較本文所提的新模糊時間序列模型與其他模型對水質指標的預測效果，表4給出了不同模型對不同水質指標的預測結果。對預測結果進行可視化分析，得到模型預測的走勢圖,分別如圖4和圖5所示。

Table 4 Model prediction results表4 模型預測結果

Figure 4 Prediction trend of dissolved oxygen圖4 溶解氧預測走勢圖

Figure 5 Prediction trend of total phosphorus圖5 總磷預測走勢圖

從圖4和圖5可以看出，對2個水質指標未來3個月的預測中，基于動態(tài)隸屬度的新模糊時間序列模型的預測值與真實值的偏差最小，說明該模型擬合最好。

為更科學貼切地反映預測效果，利用前文所述的3個評價指標，分別對上述3個模型進行評價,如表5所示。3個指標的評價結果表明：本文構建的基于動態(tài)隸屬度的模糊時間序列預測模型的預測性能最好，預測誤差指標均優(yōu)于其他對比模型的。

Table 5 Comparison of prediction error indicators of various models表5 各模型預測誤差指標的比較

將評價指標的結果進行可視化分析如圖6所示。從圖6可以看出，基于動態(tài)隸屬度的新模糊時間序列模型的預測精度最高，表明其預測效果最好，其次由于ARIMA乘積季節(jié)模型能提取較為復雜的季節(jié)效應，其預測效果次之，經典的模糊時間序列預測模型效果最差。就3種評價指標而言，本文提出的模型預測的溶解氧的RMSE相對于FTS和ARIMA乘積季節(jié)模型分別提高64.1%和30.1%，MAE分別提高68.0%和26.5%，MAPE相對于FTS和ARIMA乘積季節(jié)模型分別提高68.4%和27.5%。而預測的總磷的RMSE、MAE和MAPE都遠小于其他對比模型的。故本文根據(jù)水文時間序列的特點，將數(shù)據(jù)模糊化后再利用ARIMA乘積季節(jié)模型進行動態(tài)隸屬度的預測，并盡量建立ARIMA乘積季節(jié)模型提取隸屬度序列的季節(jié)趨勢，成功地提高了預測精度。

Figure 6 Comparison of evaluation results of dissolved oxygen prediction accuracy圖6 溶解氧預測精度評價結果柱狀圖對比

5 結束語

針對復雜多變的水質時間序列，本文基于時間序列分析方法和模糊數(shù)學基礎理論，提出了一種基于動態(tài)隸屬度的模糊時間序列預測模型。首先采用模糊C均值聚類將平穩(wěn)的序列進行論域劃分，既避免了經典的模糊時間序列中等分導致的數(shù)據(jù)分布不均勻問題，也考慮到了水質數(shù)據(jù)屬于各個污染類別時存在的亦此亦彼性；其次，將數(shù)據(jù)模糊化后得到的動態(tài)子隸屬度序列進行傳統(tǒng)的時間序列預測，充分結合了水質數(shù)據(jù)的模糊性特點和傳統(tǒng)時間序列預測法在處理動態(tài)時間序列上的優(yōu)點，克服了傳統(tǒng)時間序列預測方法要求數(shù)據(jù)較為完整和精確的弊端,與經典的模糊時間序列模型相比，無需建立模糊邏輯關系矩陣，大大簡化了計算過程。將新模型對污染較為嚴重的總磷和溶解氧的預測結果與FTS、ARIMA乘積季節(jié)模型的預測結果分別進行了比較，預測結果表明，新模糊時間序列模型的預測精度相對于其他2種模型的有較大的提高，是一種有效的短期預測模型，可為水質預警與保護提供有價值的參考。另外，是否有更加適合水質數(shù)據(jù)的模糊化方法?其他人工智能方法(如神經網絡和支持向量機)能否與該方法進行結合，從而進一步提高預測模型的精度?今后將在這方面進行更深入的探索分析。