基于協(xié)整理論的淮河流域上游洪水預(yù)報實時校正方法

張旭旻，瞿思敏，李 倩，石 朋，嵇海祥，宋蘭蘭，王麒棟

(1.河海大學水文水資源學院，江蘇 南京 210098； 2.蘇州市網(wǎng)慧水利設(shè)計咨詢有限公司，江蘇 蘇州 215128；3.水利部南京水利水文自動化研究所，江蘇 南京 210008； 4.舟山市生態(tài)環(huán)境局，浙江 舟山 316021)

流域水文模型研究采用時不變非線性系統(tǒng)，通過歷史水文數(shù)據(jù)對模型參數(shù)進行率定，用于未來洪水預(yù)報，但這種預(yù)報方案往往得不到令人滿意的結(jié)果[1]。其原因在于，流域水文系統(tǒng)是一個時變非線性系統(tǒng)，在人類活動和氣候變化的影響下，水文規(guī)律隨時間變化。此外，流域水文系統(tǒng)是一個復(fù)雜系統(tǒng)，用模型進行模擬時，常進行一系列的假設(shè)和簡化，這些假設(shè)和簡化在外延時會帶來較大誤差。包紅軍等[2-4]通過構(gòu)建模擬效果更準確的短時定量降水序列和洪水預(yù)報智能模型嘗試解決模型外延時產(chǎn)生的誤差問題。但在現(xiàn)有技術(shù)水平下，還無法從根源解決上述問題，因此，需要通過實時校正技術(shù)對預(yù)報結(jié)果進行適時、適量的校正，以提高洪水預(yù)報精度。

實時校正方法根據(jù)校正對象的不同大體上可以分為兩類：一是過程誤差校正方法(process bias correction，PBC)，二是終端誤差校正方法(terminal bias correction，TBC)[5]。PBC主要對模型輸入、參數(shù)、狀態(tài)變量等進行校正，包括遞推最小二乘校正法、卡爾曼濾波方法、動態(tài)系統(tǒng)響應(yīng)曲線方法等。已有研究表明，動態(tài)系統(tǒng)響應(yīng)曲線方法能夠?qū)δM流量結(jié)果進行有效的校正[6-7]，但在實際應(yīng)用中會出現(xiàn)“振蕩”現(xiàn)象，可通過引入平穩(wěn)矩陣減緩該缺陷。王莉莉等[8]構(gòu)建了淮河中游的Preissmann四點隱式差分格式的一維水動力模型，并利用卡爾曼濾波對構(gòu)建的水動力模型進行校正，結(jié)果表明，經(jīng)過校正的模型對水位和流量的預(yù)報精度更高。TBC方法不考慮預(yù)報過程誤差，直接對最終流量結(jié)果進行校正，從而達到實時更新原預(yù)報值的目的。TBC方法中自回歸模型(autoregressive model，AR)應(yīng)用最為廣泛，該校正模型已被應(yīng)用于多個流域的洪水預(yù)報結(jié)果實時校正中。已有研究表明，AR模型能夠有效地對預(yù)報洪水結(jié)果進行實時校正，但對洪峰流量誤差校正效果不佳[5,9-11]。其原因在于，AR模型存在一個基本假定，認為誤差序列存在序貫相關(guān)性[12]，當這種序貫相關(guān)性因為序列突變而變化時，會導(dǎo)致較差的校正結(jié)果，而洪峰前后誤差序列不可避免地存在突變，所以AR模型校正效果較差。同時AR模型需滿足時間序列為平穩(wěn)序列的前提條件，而實際上，大多數(shù)水文要素時間序列都是非平穩(wěn)的?；诜瞧椒€(wěn)時間序列運用AR模型進行分析會導(dǎo)致“偽回歸”現(xiàn)象的發(fā)生。

Engle等[13]1987年提出的協(xié)整理論為非平穩(wěn)時間序列分析提供了有力的理論依據(jù)，常被用來處理同階單整的非平穩(wěn)時間序列問題[14-15]?；趨f(xié)整理論的誤差修正模型(error correction model，ECM)通過引入一階差分項消除了變量可能存在的趨勢因素，避免了“偽回歸”問題，同時通過引入誤差修正項削弱了序貫相關(guān)性的影響。ECM模型作為計量經(jīng)濟學中解決誤差修正問題的經(jīng)典模型，在水文學中應(yīng)用較少。針對河川徑流為平穩(wěn)系列的假設(shè)所造成的“偽回歸”問題，張利亞等[14,16-17]分別將ECM模型應(yīng)用于第二松花江流域、黃河流域和渭河流域，構(gòu)建了年徑流量預(yù)測模型，取得了較好的應(yīng)用效果。

已有研究表明，分布式水文模型能夠較好地考慮流域降水與下墊面空間分布不均的特點，有效提高洪水預(yù)報精度[18-19]，同時垂向混合產(chǎn)流模型綜合考慮了蓄滿、超滲產(chǎn)流模式，因此，在半干旱半濕潤地區(qū)具有良好的應(yīng)用效果。本文采用分布式垂向混合產(chǎn)流模型模擬結(jié)果作為本次實時校正算法研究的數(shù)據(jù)來源。本文以淮河流域上游為研究區(qū)域，嘗試將協(xié)整理論與ECM模型應(yīng)用于自回歸建模中，解決AR模型無法分析非平穩(wěn)時間序列以及序貫相關(guān)性的問題，并在分布式垂向混合產(chǎn)流模型基礎(chǔ)上，建立一階至三階AR模型、ECM模型、基于誤差修正的自回歸模型(autoregressive error correction model, ARECM)對模擬結(jié)果進行校正，對比各模型校正效果，以期提高淮河流域上游洪水預(yù)報精度，為防洪調(diào)度提供更有力的科學依據(jù)。

1 研究區(qū)概況與數(shù)據(jù)來源

淮河流域地處我國東部，介于長江與黃河流域之間，流域面積18.7萬km2。本文以淮河魯臺子以上流域為研究區(qū)域，面積為8.86萬km2，氣溫由北向南升高，蒸發(fā)量南小北大，年平均蒸散發(fā)為900～1 500 mm?；春恿饔驗榈湫桶敫珊蛋霛駶櫫饔?，多年平均降水量為911 mm，多年平均徑流深為231 mm，多年平均徑流系數(shù)為0.25，其徑流年內(nèi)分布不均，主要集中在汛期。流域水系如圖1所示。研究數(shù)據(jù)來源于水文年鑒(2003—2014年)逐日逐時段實測降水、流量和蒸散發(fā)數(shù)據(jù)(逐時段蒸散發(fā)數(shù)據(jù)由逐日數(shù)據(jù)平均得到)。

圖1 魯臺子以上流域水系與站點Fig.1 Schematic diagram of water system and stations in river basin above Lutaizi station

表1 魯臺子以上流域洪水模擬結(jié)果Table 1 Results of flood simulation in river basin above Lutaizi station

2 分布式垂向混合產(chǎn)流模型構(gòu)建

垂向混合產(chǎn)流模型由包為民等[20]于1997年提出，模型將超滲產(chǎn)流與蓄滿產(chǎn)流進行垂向組合，降水到達地面后，首先通過下滲能力分布曲線劃分為地面徑流和下滲水流。在土壤缺水量大的區(qū)域，下滲水流補充土壤含水量，不產(chǎn)流；在缺水量小的區(qū)域，下滲水流補充土壤含水量達到飽和后，產(chǎn)生地面以下徑流。由于同時考慮了超滲和蓄滿兩種產(chǎn)流模式，該模型既適用于濕潤地區(qū)產(chǎn)匯流過程模擬，也適用于半干旱半濕潤地區(qū)產(chǎn)匯流過程模擬[20-21]。為了更好地考慮降水和下墊面分布不均對流域產(chǎn)匯流過程的影響，本文基于垂向混合產(chǎn)流模型和DEM資料，構(gòu)建了柵格型分布式垂向混合產(chǎn)流模型，根據(jù)2003—2014年淮河魯臺子站實測資料，選取10場洪水對模型進行參數(shù)率定與檢驗，模型模擬結(jié)果如表1所示。結(jié)果表明，模型平均確定性系數(shù)為 0.74，模擬結(jié)果合格率為70%，能夠較好地模擬淮河流域上游洪水過程。

3 實時校正方法

3.1 ECM模型

ECM模型是一種具有特定形式的計量經(jīng)濟學模型。首先需要對自變量x與因變量y進行協(xié)整分析，以發(fā)現(xiàn)變量之間的協(xié)整關(guān)系(即長期均衡關(guān)系)，假設(shè)協(xié)整關(guān)系的回歸方程形式為

yt=b0+b1xt+εt

(1)

式中：b0、b1為系數(shù)；εt為t時刻誤差修正項，即回歸方程的殘差序列，反映t時刻因變量在短期波動中偏離長期均衡關(guān)系的程度。

將誤差修正項看作一個解釋變量，連同其他反映短期波動的解釋變量一起，建立短期模型，即ECM模型。誤差修正模型的主要形式為

Δyt=b0+b1Δxt+γεt-1+μt

(2)

式中：Δyt、Δxt為t時刻因變量與自變量的一階差分；εt-1為t-1時刻誤差修正項；μt為純隨機數(shù)列；γ為系數(shù)，多為負值。

根據(jù)實際應(yīng)用的場景，可對ECM模型的形式進行變化，首先將觀測流量作為因變量，模擬流量的一階滯后項作為自變量進行協(xié)整分析，得到

Qobs,t=b0+b1Qsim,t-1+εt

(3)

式中：Qobs,t為t時刻實測流量；Qsim,t-1為t-1時刻模擬流量。

基于式(3)構(gòu)建的短期模型為

ΔQ′obs,t=b0+b1ΔQsim,t-1+γεt-1+μt

(4)

其中

εt-1=Qobs,t-1-b0-b1Qsim,t-2

式中：ΔQ′obs,t為估計的t時刻實測流量數(shù)據(jù)一階差分；ΔQsim,t-1為t-1時刻模擬流量數(shù)據(jù)一階差分。

假設(shè)具有t時刻的實測流量數(shù)據(jù)以及t+1、…、t+n時刻的預(yù)報流量數(shù)據(jù)，在利用ECM模型對預(yù)報數(shù)據(jù)進行滾動校正時，首先計算出ΔQsim,t與εt，根據(jù)式(4)計算得出ΔQ′obs,t+1，將ΔQ′obs,t+1與Qobs,t相加即可得到t+1時刻校正后流量結(jié)果Q′obs,t+1；類似地，對t+n+1時刻預(yù)報數(shù)據(jù)進行校正時，根據(jù)ΔQsim,t+n與εt+n計算得出ΔQ′obs,t+n+1，將ΔQ′obs,t+n+1與Q′obs,t+n相加即可得到t+n+1時刻校正后流量結(jié)果Q′obs,t+n+1。因此，假設(shè)分布式垂向混合產(chǎn)流模型預(yù)見期為n，則ECM模型預(yù)見期為n+1。式(4)為一階ECM模型，二階ECM模型如式(5)所示，三階模型形式相似。

ΔQ′obs,t=b0+b1ΔQsim,t-1+b2ΔQsim,t-2+γεt-1+μt

(5)

式中b2為系數(shù)。

3.2 AR模型

AR模型基于變量的自相關(guān)性，認為t時刻變量的值可以用其之前各時刻的值來預(yù)測。在實時校正過程中，根據(jù)實測流量與預(yù)報流量計算誤差序列，在假設(shè)誤差序列為平穩(wěn)序列的前提下，建立誤差A(yù)R模型，估計下一時刻誤差值，從而達到校正效果。AR模型為

(6)

式中:et為t時刻預(yù)報值與實測值之間誤差；ak為系數(shù)；k為自回歸階數(shù)；ε′t為t時刻隨機誤差，為零均值獨立同分布誤差序列。

3.3 ARECM模型

依據(jù)Granger表述定理，如果變量x與y是協(xié)整的，則變量間的短期非均衡關(guān)系總能由ECM模型表述[13]。為便于解釋推導(dǎo)，假設(shè)根據(jù)實測流量與預(yù)報流量構(gòu)建的誤差序列存在如下一階自回歸形式的長期均衡關(guān)系：

et=α0+α1et-1+εt

(7)

式中α0、α1為系數(shù)。

實際上，單場洪水誤差序列很少處在均衡點上，因此根據(jù)單場洪水建立的誤差A(yù)R模型只能反映誤差序列的短期或非均衡關(guān)系。假設(shè)針對單場洪水具有如下二階自回歸形式模擬方程：

et=β0+β1et-1+β2et-2+μt

(8)

式中β0、β1、β2為系數(shù)。該模型顯示，t時刻誤差值et可由et-1和et-2估計得到。由于誤差序列可能是非平穩(wěn)的，因此不能直接應(yīng)用普通最小二乘法(ordinary least square，OLS)計算參數(shù)取值。根據(jù)式(8)可得

Δet=β0+β1et-1-et-1+β2et-2+μt

(9)

對式(9)進一步變換后可得

Δet=β1Δet-1-λ(et-1-γ0-γ1et-2)+μt

(10)

其中

γ0=β0/λγ1=(β1+β2)/λ

式中λ為系數(shù)。如果et為一階單整序列，則Δet為平穩(wěn)序列，可以應(yīng)用OLS法計算模型參數(shù)取值。令γ0=α0，γ1=α1，由式(10)可得：

Δet=β1Δet-1-λεt-1+μt

(11)

式(11)通過引入一階差分項消除了變量可能存在的趨勢因素，避免了“偽回歸”與可能存在的多重共線性問題，解決了非平穩(wěn)誤差序列的自回歸建模問題，同時通過引入誤差修正項，體現(xiàn)了變量水平值信息。式(11)為本文采用的一階ARECM模型，式(12)為二階ARECM模型，三階模型形式相似。

Δet=β1Δet-1+β2Δet-2-λεt-1+μt

(12)

3.4 評價指標

選取修正效果評價系數(shù)[1]與確定性系數(shù)作為評價校正效果優(yōu)劣的指標。并根據(jù)GB/T 22482—2008《水文情報預(yù)報規(guī)范》，針對洪峰與徑流量分別選取洪峰相對誤差、峰現(xiàn)時差、徑流深相對誤差對洪水要素校正效果進行詳細的評價。

修正效果評價系數(shù)為

(13)

式中Qraw,t為校正后t時刻流量。

4 實時校正模型應(yīng)用效果分析

4.1 平穩(wěn)性檢驗

采用廣泛應(yīng)用于時間序列平穩(wěn)性檢驗的Augment Dickey-Fuller (ADF)檢驗法對10場洪水誤差系列進行平穩(wěn)性檢驗?；贏DF檢驗的回歸方程為

(14)

式中：α、β、δ為系數(shù)；p為總滯后階數(shù)；ζt為白噪聲序列；ξi為滯后階數(shù)為i的線性時間趨勢項。

對淮河流域上游10場洪水的誤差序列進行ADF檢驗，結(jié)果如表2所示。10場洪水中除20120907與20140827次洪水誤差序列外，其余洪水誤差序列均為平穩(wěn)序列。本研究中，為了體現(xiàn)ARECM模型在非平穩(wěn)時間序列建模中的優(yōu)勢，對20120907與 20140827次洪水依舊構(gòu)建AR模型，對垂向混合產(chǎn)流模型模擬結(jié)果進行校正。需要強調(diào)的是，針對20120907與20140827次洪水構(gòu)建的AR模型為“偽回歸”模型，并沒有實用價值。

表2 ADF檢驗結(jié)果Table 2 Results of ADF test

4.2 協(xié)整檢驗

協(xié)整關(guān)系能夠反映變量間長期存在的一種穩(wěn)定均衡關(guān)系。對多變量應(yīng)用ECM模型時需滿足變量間具有協(xié)整關(guān)系的前提條件，因此，對變量進行協(xié)整檢驗是很有必要的。協(xié)整檢驗[13]核心在于對回歸方程殘差進行單位根檢驗(即平穩(wěn)性檢驗)，若殘差序列為單整序列(即平穩(wěn)序列)，則變量間存在協(xié)整關(guān)系。本文對誤差修正模型中的變量進行協(xié)整檢驗，各變量間協(xié)整檢驗結(jié)果的P值如表3所示。當P>0.05時，認為在5%置信度區(qū)間內(nèi)接受原假設(shè)；當P>0.1時，認為在10%置信度區(qū)間內(nèi)接受原假設(shè)，而原假設(shè)為不具有協(xié)整關(guān)系。

表3 協(xié)整檢驗結(jié)果的P值Table 3 P values of cointegration test results

在建立ECM模型的10場洪水中，對模型模擬值與實測值進行協(xié)整檢驗，其中7場洪水P值小于0.1，在10%置信度區(qū)間內(nèi)拒絕原假設(shè)，認為序列間具有協(xié)整關(guān)系，通過了協(xié)整檢驗。未通過10%置信度協(xié)整檢驗的3場洪水中，有2場洪水P值為 0.14，可以近似認為通過了協(xié)整檢驗；而20140827次洪水，P值為0.47，未通過10%置信度協(xié)整檢驗，因此，不采用ECM模型對該場洪水進行誤差校正。在一階ARECM模型建立中，對et與et-1進行協(xié)整檢驗，10場洪水中有8場洪水通過了5%置信度協(xié)整檢驗，其中，未通過協(xié)整檢驗的2場洪水的P值分別為0.11與0.14，可以近似認為通過10%置信度協(xié)整檢驗，能夠用于模型建立。針對二階與三階ARECM模型，所有變量在各場次洪水中均通過了5%置信度的協(xié)整檢驗，可用于建立ARECM模型。

4.3 模型效果分析

選取流域2003—2014年10場洪水資料以及分布式垂向混合產(chǎn)流模型模擬結(jié)果，分別建立一階至三階AR、ECM和ARECM實時校正模型對模擬結(jié)果進行校正，并對校正效果進行對比分析。

從校正結(jié)果的修正效果評價系數(shù)(表4)看，ECM與ARECM模型相較于AR模型具有明顯優(yōu)勢，3種方法的校正效果與模型階數(shù)基本成正比，但針對某一場洪水，建議依照最小信息量準則確定校正模型的階數(shù)。AR、ECM、ARECM模型平均修正效果評價系數(shù)為-0.77、0.76、0.98。AR模型表現(xiàn)效果較差，原因在于20120907與20140827次洪水誤差序列不滿足平穩(wěn)序列的假設(shè)，校正效果較差，刪除以上兩場洪水后，AR模型的平均修正效果評價系數(shù)提高到0.20。針對20120907與20140827次洪水，ARECM模型修正效果評價系數(shù)明顯優(yōu)于AR模型，其中ARECM模型修正效果評價系數(shù)分別達到了1.00與0.92遠高于AR模型的-9.15與-0.65(一階至三階模型的平均值)，證明了ARECM模型在針對非平穩(wěn)誤差序列建模方面具有很強的能力，同時證明了可以通過引入?yún)f(xié)整理論解決AR模型無法針對非平穩(wěn)誤差序列建模的問題。根據(jù)確定性系數(shù)的統(tǒng)計結(jié)果(表5)，ARECM模型校正效果明顯優(yōu)于ECM模型與AR模型，且ECM模型校正效果優(yōu)于AR模型，證明了ECM與ARECM模型在實時校正方面的適用性。

由表6可見，以洪峰相對誤差為評價指標，AR模型校正效果最差，ECM模型與ARECM模型校正后的洪峰相對誤差相較于AR模型分別減小了11.22%與17.41%。其原因在于ARECM模型通過引入誤差修正項，保證了變量的水平值沒有被忽視，減弱了誤差序列的序貫相關(guān)性對建模的影響，使得模型在洪峰處不僅依靠誤差序列的自相關(guān)性進行校正，同時也考慮了校正后誤差序列與實測誤差序列上一時刻的偏離程度，避免了AR模型在誤差突變點處產(chǎn)生較大誤差的現(xiàn)象。峰現(xiàn)時差(表7)與洪峰相對誤差結(jié)果相似，ARECM模型與ECM模型校正效果遠優(yōu)于AR模型，其中，20120907與20140827次洪水不滿足平穩(wěn)序列假設(shè)，但仍構(gòu)建了AR模型，導(dǎo)致校正效果較差，因此沒有在表7中計算這兩場洪水二、三階AR模型校正后的峰現(xiàn)時差。由表8可見，以徑流深相對誤差為評價指標，ECM模型與AR模型對模擬結(jié)果的校正效果相近，不如ARECM模型。

表4 修正效果評價系數(shù)統(tǒng)計結(jié)果Table 4 Statistics of correction effect evaluation coefficient

表5 確定性系數(shù)統(tǒng)計結(jié)果Table 5 Statistics of deterministic coefficient

表6 洪峰相對誤差統(tǒng)計結(jié)果 單位：%Table 6 Statistics of relative error of flood peaks unit: %

綜上所述，AR模型校正效果最差，其主要原因在于AR模型以誤差序列的自相關(guān)性為基礎(chǔ)進行校正，無法對洪峰處誤差進行較好的校正；ECM模型依賴實測值與預(yù)報值，因此在洪峰處模型校正效果較好，但徑流深方面，校正效果與AR模型相似。ARECM模型模擬效果較好，通過誤差修正項的引入，避免了AR模型在洪峰處校正效果差的缺陷，提高了洪峰誤差校正效果，同時通過引入?yún)f(xié)整理論解決了AR模型無法對非平穩(wěn)序列建模的問題。

為了直觀地展示校正效果，選取2場典型洪水(20030620和20100811次洪水)，對比3種三階實時校正模型校正后的流量過程線(圖2)。從圖2可以看出，2場洪水的流量模擬與實測結(jié)果相差較大；AR模型以實測與模擬流量的誤差序列為基礎(chǔ)進行校正，對模擬流量線型影響較??；ECM模型直接考慮實測流量的影響，能夠顯著改變模擬流量線型，使校正結(jié)果更貼近實測流量；ARECM模型校正效果最好，能夠?qū)δM洪峰流量進行準確的校正，顯著提高了淮河流域上游洪水預(yù)報精度。

表7 峰現(xiàn)時差統(tǒng)計結(jié)果 單位：hTable 7 Statistics of error of time to flood peaks unit: h

表8 徑流深相對誤差統(tǒng)計結(jié)果 單位：%Table 8 Statistics of relative error of runoff depth unit: %

(a) 20030620次洪水

5 結(jié) 論

a.AR模型能夠有效地對垂向混合產(chǎn)流模型洪水模擬結(jié)果進行校正，但由于序貫相關(guān)性的影響，在洪峰處校正效果較差，而在徑流深方面與ECM模型校正效果相似。

b.ECM模型能夠較好地解決洪水實時校正問題，模型校正效果優(yōu)于AR模型但劣于ARECM模型，從洪峰相對誤差角度，校正效果相較于AR模型改善了11.22%。

c.ARECM模型通過引入誤差修正項，減弱了誤差序列序貫相關(guān)性對建模的影響，并通過引入一階差分項，避免了“偽回歸”與可能存在的多重共線性問題，避免了AR模型對洪峰校正效果較差的缺陷。模型校正后平均洪峰相對誤差為1.28%，平均徑流深相對誤差為1.49%，從洪峰相對誤差角度看，ARECM模型校正效果相較于AR模型改善了17.41%，校正效果較好，可用于校正實時洪水預(yù)報結(jié)果。