諸葛燕，方海泉

(1.浙江經(jīng)濟(jì)職業(yè)技術(shù)學(xué)院，浙江杭州 310018；2.中國航天系統(tǒng)科學(xué)與工程研究院，北京 100048)

我國實行最嚴(yán)格水資源管理制度，健全水資源監(jiān)控體系是實施這一制度的重要舉措。國務(wù)院《關(guān)于實行最嚴(yán)格水資源管理制度的意見》要求，加強(qiáng)取水、排水、入河湖排污口計量監(jiān)控設(shè)施建設(shè)。我國各地普遍開展了取用水監(jiān)測，積累了大量監(jiān)測數(shù)據(jù)[1-5]。為了有效發(fā)揮取用水監(jiān)測數(shù)據(jù)在水資源管理中的作用，需要對監(jiān)測數(shù)據(jù)進(jìn)行整編、清洗、加工，辨識異常數(shù)據(jù)，保留正常數(shù)據(jù)，從而對一個地區(qū)取用水監(jiān)測數(shù)據(jù)的有效性進(jìn)行評估。

1 取用水監(jiān)測數(shù)據(jù)異常值分析

取用水監(jiān)測數(shù)據(jù)的有效性評估要以數(shù)據(jù)真實可靠為基礎(chǔ)，為此，首先要對監(jiān)測數(shù)據(jù)進(jìn)行異常值分析，包括監(jiān)測數(shù)據(jù)總量、監(jiān)測數(shù)據(jù)時序變化兩個方面。

1.1 監(jiān)測數(shù)據(jù)總量異常

將1 個取用水戶同一天所有取水管的監(jiān)測數(shù)據(jù)求和，得到該戶的日取水量(以下將每個受監(jiān)測的取用水戶簡稱為“監(jiān)測點(diǎn)”)；通過計算一段時間的日取水量監(jiān)測數(shù)據(jù)，推算出1 年的取水量，即監(jiān)測年取水量。如果監(jiān)測年取水量與實際年取水量偏差太大，說明監(jiān)測數(shù)據(jù)可能存在異常，包括以下兩種情況。

1.1.1 監(jiān)測年取水量高于實際年取水量

用S1表示1 個監(jiān)測點(diǎn)的監(jiān)測年取水量，S2表示該監(jiān)測點(diǎn)實際年取水量。把監(jiān)測年取水量除以實際年取水量的比值稱為正監(jiān)測率，記為g.

計算出每一個監(jiān)測點(diǎn)的正監(jiān)測率得到數(shù)組｛gm｝，m 表示監(jiān)測點(diǎn)的個數(shù)。若g 越大，說明監(jiān)測年取水量高于實際年取水量，則該監(jiān)測點(diǎn)監(jiān)測數(shù)據(jù)異常的可能性就越大。

1.1.2 監(jiān)測年取水量低于實際年取水量

把正監(jiān)測率的倒數(shù)稱為倒監(jiān)測率，記為d.計算出每一個監(jiān)測點(diǎn)的倒監(jiān)測率得到數(shù)組｛dm｝。若d 越大，說明監(jiān)測年取水量低于實際年取水量，則該監(jiān)測點(diǎn)監(jiān)測數(shù)據(jù)異常的可能性就越大。d 異常大的一種原因可能是該監(jiān)測點(diǎn)只監(jiān)測了部分取水管，還有其它取水管沒有監(jiān)測到，即監(jiān)測不全面。

1.2 監(jiān)測數(shù)據(jù)時序變化異常

對區(qū)域內(nèi)每一個監(jiān)測點(diǎn)日取水量數(shù)據(jù)繪制時間序列圖，根據(jù)圖形變化趨勢可能存在三類異常值現(xiàn)象：某日數(shù)值異常大、某個數(shù)據(jù)出現(xiàn)頻率異常高、某日數(shù)據(jù)為負(fù)值，對每一類異?，F(xiàn)象抽象出數(shù)學(xué)算法，通過MATLAB 編程可以快速、準(zhǔn)確分析出每一類異常情況分別存在于哪些監(jiān)測點(diǎn)。

1.2.1 某日數(shù)值異常大

監(jiān)測點(diǎn)日取水量時間序列數(shù)據(jù)記為｛an｝，n 表示時間序列的天數(shù)，該數(shù)組中大于0 的數(shù)構(gòu)成的數(shù)組稱為正數(shù)組，記為｛ank｝，ank＞0，｛ank｝的中位數(shù)記為median ｛ank｝，｛an｝的最大值記為max ｛an｝。監(jiān)測點(diǎn)日取水量時間序列數(shù)據(jù)的最大值與中位數(shù)比值，即最大偏離率，記為。r 越大，則該監(jiān)測點(diǎn)時間序列數(shù)據(jù)的最大值為異常值的可能性越大。計算出每一個監(jiān)測點(diǎn)的r 得到1 個數(shù)組｛rm｝。

1.2.2 某個數(shù)據(jù)出現(xiàn)頻率異常高

把時間序列數(shù)據(jù)｛an｝中每一個不重復(fù)的數(shù)構(gòu)成的集合記為，集合中每個數(shù)在｛an｝中出現(xiàn)的頻率記為P(ant)，該監(jiān)測點(diǎn)中出現(xiàn)頻率最高的數(shù)據(jù)記為h，P(h)=max(P(ant))(簡記為P)，即1 個監(jiān)測點(diǎn)的最高頻率。P 越大，則該監(jiān)測點(diǎn)監(jiān)測數(shù)據(jù)中某個數(shù)據(jù)出現(xiàn)的頻率越高，該監(jiān)測點(diǎn)數(shù)據(jù)異常的可能性就越大。計算出每一個監(jiān)測點(diǎn)的P 得到一個數(shù)組｛Pm｝。

1.2.3 某日數(shù)據(jù)為負(fù)值

對于時間序列數(shù)據(jù)｛an｝，若an＜0，則該日取水量數(shù)據(jù)為異常值，｛an｝中小于0 的數(shù)據(jù)的個數(shù)記為w，即1 個監(jiān)測點(diǎn)的負(fù)值數(shù)，計算出每一個監(jiān)測點(diǎn)的得到1 個數(shù)組｛wm｝。

1.3 異常情況綜合分析

對每一個監(jiān)測點(diǎn)的取用水監(jiān)測數(shù)據(jù)，綜合考慮監(jiān)測數(shù)據(jù)總量和監(jiān)測數(shù)據(jù)時序變化兩類異常情況，據(jù)此建立異常情況綜合值的數(shù)學(xué)模型。

其中k01，k02，k03，k04，k05為權(quán)重系數(shù)分別為g，d，r，p，w 標(biāo)準(zhǔn)化之后的值。權(quán)重系數(shù)的確定采用層次分析法中的“兩兩比較法”[6]。

計算出每一個監(jiān)測點(diǎn)的異常情況綜合值x 得到數(shù)組｛xm｝。異常情況綜合值x 越大，該監(jiān)測點(diǎn)監(jiān)測數(shù)據(jù)異常的可能性也就越大。采用Grubbs 檢驗法檢驗｛xm｝中的異常大值[7]。再進(jìn)一步判斷異常數(shù)據(jù)是否可以校正，對于單個時間點(diǎn)數(shù)據(jù)異常大或者出現(xiàn)負(fù)值的數(shù)據(jù)可以采用時間序列預(yù)測方法校正異常值[8]，但對于長時間數(shù)據(jù)異常導(dǎo)致正監(jiān)測率或者倒監(jiān)測率異常大的數(shù)據(jù)是無法校正的。

2 取用水監(jiān)測數(shù)據(jù)有效性評估

2.1 有效性評估指標(biāo)體系

對一個地區(qū)的取用水監(jiān)測數(shù)據(jù)有效性進(jìn)行評估，不僅要考察監(jiān)測數(shù)據(jù)正常的監(jiān)測點(diǎn)的數(shù)量，而且要考察監(jiān)測數(shù)據(jù)正常的監(jiān)測點(diǎn)所監(jiān)測水量的大小。因此要建立包括監(jiān)測點(diǎn)比率、監(jiān)測覆蓋率、監(jiān)測水量比率的指標(biāo)體系。

監(jiān)測點(diǎn)比率y1是指某地區(qū)監(jiān)測點(diǎn)的監(jiān)測數(shù)據(jù)經(jīng)過異常值分析，排除有異常且無法校正的數(shù)據(jù)，剩下監(jiān)測數(shù)據(jù)正常的監(jiān)測點(diǎn)的個數(shù)占該地區(qū)有上報監(jiān)測數(shù)據(jù)的監(jiān)測點(diǎn)個數(shù)的比例。

監(jiān)測覆蓋率y2是指一個地區(qū)監(jiān)測數(shù)據(jù)正常的監(jiān)測點(diǎn)的年監(jiān)測水量占該地區(qū)年取水總量的比例。

監(jiān)測水量比率y3是指一個地區(qū)監(jiān)測數(shù)據(jù)正常的監(jiān)測點(diǎn)的年監(jiān)測水量占該地區(qū)有上報監(jiān)測數(shù)據(jù)的監(jiān)測點(diǎn)年取水量的比例。

本研究所考慮的監(jiān)測點(diǎn)不包括直流冷卻型火電廠監(jiān)測點(diǎn)。

2.2 有效性評估模型

取用水監(jiān)測數(shù)據(jù)有效性評估需要綜合考慮監(jiān)測點(diǎn)比率y1、監(jiān)測覆蓋率y2、監(jiān)測水量比率y33 個指標(biāo)，據(jù)此建立以下數(shù)學(xué)模型。

其中k1，k2，k3為屬于(0,1)的權(quán)重系數(shù)，且k1+k2+k3=1，權(quán)重系數(shù)的確定同樣采用層次分析法中的“兩兩比較法”。由于y1，y2，y3都屬于(0,1)，可知z∈（0,1）.z越大，則說明該地區(qū)的取用水監(jiān)測數(shù)據(jù)越有效。

取用水監(jiān)測數(shù)據(jù)經(jīng)過異常值分析可得到監(jiān)測數(shù)據(jù)正常的監(jiān)測點(diǎn)，為取用水監(jiān)測數(shù)據(jù)有效性分析打下基礎(chǔ)，監(jiān)測數(shù)據(jù)有效性分析建模流程圖如圖1 所示。

圖1 數(shù)據(jù)有效性分析流程圖Fig.1 Flow chart of data validity analysis

3 實例研究

3.1 研究對象

選取G 市的所有監(jiān)測點(diǎn)的取用水監(jiān)測數(shù)據(jù)作為研究對象。G 市有上報監(jiān)測數(shù)據(jù)的監(jiān)測點(diǎn)有27 個。不包括直流冷卻型火電廠。

每個監(jiān)測點(diǎn)的實際年取水量數(shù)據(jù)難以獲取，本研究采用2011 年水利普查時的年取水量(簡稱水利普查數(shù)據(jù))作為實際年取水量，對少部分在2011 年沒有被普查到的監(jiān)測點(diǎn)，采用該監(jiān)測點(diǎn)的取水許可證許可水量(簡稱許可證數(shù)據(jù))推算出年取水量。理論上，取水許可證的許可水量與實際用水量在統(tǒng)計學(xué)上有一定相關(guān)性，由此可以根據(jù)同時具備水利普查數(shù)據(jù)與許可證數(shù)據(jù)的監(jiān)測點(diǎn)，找到二者的關(guān)系，從而推算出沒有水利普查數(shù)據(jù)的監(jiān)測點(diǎn)在2011 年的取水量，以此作為該監(jiān)測點(diǎn)的年取水量。這里假定取水戶每年的取水量維持在相對穩(wěn)定的水平，也即假定取水戶在2011 年的取水量與2016 年的取水量變化不大。

3.2 監(jiān)測數(shù)據(jù)預(yù)處理

對于G 市所有的監(jiān)測數(shù)據(jù)，整理出每一個監(jiān)測點(diǎn)的日取水量時間序列數(shù)據(jù)?？疾鞎r間范圍從2016 年2 月1 日到2016 年7 月31 日，共計182 d，相當(dāng)于半年時間，以此半年的取水量的2 倍近似為1 年的監(jiān)測年取水量。

3.3 監(jiān)測數(shù)據(jù)總量異常值分析

3.3.1 監(jiān)測數(shù)據(jù)、水利普查數(shù)據(jù)、許可證數(shù)據(jù)整理

對于G 市有上報監(jiān)測數(shù)據(jù)的27 個監(jiān)測點(diǎn)，分別計算出每一個監(jiān)測點(diǎn)半年總用水量S0，進(jìn)而推算該監(jiān)測點(diǎn)監(jiān)測年取水量S1=2S0，監(jiān)測年取水量列于表1 的第2 列。

根據(jù)監(jiān)測點(diǎn)名稱找出對應(yīng)的水利普查數(shù)據(jù)和許可證數(shù)據(jù)。把查找到的水利普查數(shù)據(jù)和許可證數(shù)據(jù)分別列于表1 的第3 列和第4 列，其中有4 個監(jiān)測點(diǎn)沒有水利普查數(shù)據(jù)。

表1 監(jiān)測點(diǎn)的三類數(shù)據(jù)Tab.1 Three types of data of monitoring station

3.3.2 用許可證數(shù)據(jù)擬合出水利普查數(shù)據(jù)

把水利普查數(shù)據(jù)和許可證數(shù)據(jù)都存在的監(jiān)測點(diǎn)單獨(dú)列出來，由于數(shù)量級差別較大，為了提高擬合準(zhǔn)確度，把水利普查數(shù)據(jù)分為兩部分進(jìn)行擬合，以10 000×104m3為分界點(diǎn)。選出水利普查數(shù)據(jù)小于10 000×104m3的監(jiān)測點(diǎn)，計算出水利普查數(shù)據(jù)與許可證數(shù)據(jù)的相關(guān)系數(shù)為0.91，說明二者有較強(qiáng)的相關(guān)性。圖2 是水利普查數(shù)據(jù)小于10 000×104m3的監(jiān)測點(diǎn)的許可證數(shù)據(jù)與水利普查數(shù)據(jù)的關(guān)系圖，直線上的星號是表示根據(jù)許可證數(shù)據(jù)擬合出空缺的水利普查數(shù)據(jù)。

圖2 許可證數(shù)據(jù)與水利普查數(shù)據(jù)的關(guān)系Fig.2 The relationship between license data and census data

同理作出水利普查數(shù)據(jù)大于10 000×104m3的監(jiān)測點(diǎn)的許可證數(shù)據(jù)與水利普查數(shù)據(jù)的關(guān)系圖，也可以根據(jù)許可證數(shù)據(jù)擬合出空缺的水利普查數(shù)據(jù)。表1 第3 列中的加粗字體是表示由許可證數(shù)據(jù)擬合出的水利普查數(shù)據(jù)。

3.3.3 計算正監(jiān)測率和倒監(jiān)測率

根據(jù)表1 的監(jiān)測年取水量和水利普查數(shù)據(jù)，計算出G 市27 個監(jiān)測點(diǎn)的正監(jiān)測率g 和倒監(jiān)測率d，計算結(jié)果列于表2 的第2 列和第3 列。

表2 監(jiān)測數(shù)據(jù)異常情況的參數(shù)值和綜合值Tab.2 The parameter value and comprehensive value of abnormal monitoring data

3.4 監(jiān)測數(shù)據(jù)時序變化異常值分析

應(yīng)用MATLAB 編程計算出每個監(jiān)測點(diǎn)的最大偏離率r、最高頻率p、負(fù)值數(shù)w，結(jié)果發(fā)現(xiàn)w 都為0，也即該市的監(jiān)測數(shù)據(jù)沒有出現(xiàn)負(fù)值。最大偏離率、最高頻率的計算結(jié)果列于表2 的第4 列和第5 列.

3.5 計算監(jiān)測數(shù)據(jù)異常情況綜合值

采用“兩兩比較法”確定正監(jiān)測率g、倒監(jiān)測率d、最大偏離率r、最高頻率p 的權(quán)重系數(shù)，經(jīng)過咨詢水資源專家得到比較矩陣A1.

求出權(quán)重系數(shù)

經(jīng)過計算能夠通過一致性檢驗，求得的權(quán)重系數(shù)可以使用。

綜合計算公式如下：

計算出監(jiān)測點(diǎn)監(jiān)測數(shù)據(jù)異常情況綜合值x 列于表2 的第6 列，表2 的行排列順序就是按照x 值從大到小排序的。對x值采用Grubbs 檢驗得到前3 個監(jiān)測點(diǎn)監(jiān)測數(shù)據(jù)異常情況綜合值為異常大值，可初步判斷這3 個監(jiān)測點(diǎn)監(jiān)測數(shù)據(jù)屬于異常。經(jīng)過分析發(fā)現(xiàn)第一個監(jiān)測點(diǎn)“M1 自來水公司”監(jiān)測水量過小，實際用水量是監(jiān)測水量的12 倍多。第二個監(jiān)測點(diǎn)“M2 紡織印染有限公司”監(jiān)測水量過大，是實際用水量的22 倍多。第三個監(jiān)測點(diǎn)“M3 水利投資有限公司”監(jiān)測數(shù)據(jù)存在異常大值，如果是單個數(shù)據(jù)異常大可以校正，但是從該監(jiān)測點(diǎn)的監(jiān)測數(shù)據(jù)時間序列圖(圖3)可以看出不只是單個數(shù)據(jù)異常大，該監(jiān)測點(diǎn)前1 個月的數(shù)據(jù)較大，后5 個月的數(shù)據(jù)相對較小，導(dǎo)致總量低于擬合出的水利普查數(shù)據(jù)?？梢源_定這3個監(jiān)測點(diǎn)監(jiān)測數(shù)據(jù)屬于異常的可能性非常大，并且這3 個監(jiān)測點(diǎn)的監(jiān)測數(shù)據(jù)都無法校正。

圖3 監(jiān)測數(shù)據(jù)時序變化圖Fig.3 The time series change of monitoring data

3.6 計算監(jiān)測數(shù)據(jù)有效性評估值

篩選出數(shù)據(jù)異常的監(jiān)測點(diǎn)，從而確定出監(jiān)測數(shù)據(jù)正常的監(jiān)測點(diǎn)，計算出以下三個評價指標(biāo)值。

3.6.1 監(jiān)測點(diǎn)比率

對于G 市，排除以上3 個監(jiān)測數(shù)據(jù)異常的監(jiān)測點(diǎn)，還剩下24 監(jiān)測數(shù)據(jù)正常的監(jiān)測點(diǎn)。監(jiān)測點(diǎn)比率y1=24/27=0.888 9.

3.6.2 監(jiān)測覆蓋率

監(jiān)測覆蓋率y2=7.5043 億m3/43.45 億m3=0.172 7.

分子為排除3 個監(jiān)測數(shù)據(jù)異常的監(jiān)測點(diǎn)之后剩下的24 個監(jiān)測數(shù)據(jù)正常的監(jiān)測點(diǎn)監(jiān)測年水量之和，用表1第2 列后24 個監(jiān)測點(diǎn)的監(jiān)測年取水量。分母為G 市2015 年取用水量減去直流式火(核)電用水，該數(shù)據(jù)來源G 市所在省的2015 年水資源公報。

3.6.3 監(jiān)測水量比率

監(jiān)測水量比率y3=7.5043 億m3/22.2696 億m3=0.3369.

分子與監(jiān)測覆蓋率的分子一致，分母為G 市27 個有上報監(jiān)測數(shù)據(jù)的監(jiān)測點(diǎn)實際年取水量之和，用表1 第3 列的水利普查數(shù)據(jù)作為實際年取水量。

3.6.4 取用水監(jiān)測數(shù)據(jù)有效性評估

采用“兩兩比較法”確定y1，y2，y3權(quán)重系數(shù)，經(jīng)過咨詢水資源專家得到比較矩陣A2.

求出權(quán)重系數(shù)

經(jīng)過計算能夠通過一致性檢驗，求得的權(quán)重系數(shù)可以使用。

取用水監(jiān)測數(shù)據(jù)有效性評估值的計算

3.7 討論

已計算出G 市取用水監(jiān)測數(shù)據(jù)有效性評估值為0.67，距離理想值1 還有提升空間，若要提高G 市的取用水監(jiān)測數(shù)據(jù)的有效性，建議從三個方面加以改進(jìn)，一是對于監(jiān)測不全面的監(jiān)測點(diǎn)要確保監(jiān)測到所有的取水管；二是對于監(jiān)測數(shù)據(jù)大于或小于真實取水量數(shù)據(jù)的監(jiān)測點(diǎn)應(yīng)該檢查監(jiān)測設(shè)備以及通訊設(shè)施以確保監(jiān)測數(shù)據(jù)真實可靠；三是在一定程度上增加監(jiān)測點(diǎn)數(shù)量，并且盡量選擇取用水大戶進(jìn)行監(jiān)測。

4 結(jié)論

對一個地區(qū)的取用水在線監(jiān)測數(shù)據(jù)從數(shù)據(jù)的總量和時序變化兩個方面進(jìn)行異常值分析，能夠有效篩選出異常監(jiān)測數(shù)據(jù)。進(jìn)一步對取用水在線監(jiān)測數(shù)據(jù)進(jìn)行有效性評估，評估結(jié)果從一定程度上能夠反映出該地區(qū)監(jiān)測數(shù)據(jù)的質(zhì)量，從而能夠發(fā)現(xiàn)存在的問題并有針對性地改進(jìn)，為取用水在線監(jiān)測數(shù)據(jù)服務(wù)于水資源管理業(yè)務(wù)提供更加可靠的保障，也可以為將來智慧水利、智慧流域的監(jiān)測點(diǎn)建設(shè)提供借鑒[9-11]。