基于協整理論的淮河流域上游洪水預報實時校正方法

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

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

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

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

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

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

1 研究區概況與數據來源

淮河流域地處我國東部，介于長江與黃河流域之間，流域面積18.7萬km2。本文以淮河魯臺子以上流域為研究區域，面積為8.86萬km2，氣溫由北向南升高，蒸發量南小北大，年平均蒸散發為900～1 500 mm。淮河流域為典型半干旱半濕潤流域，多年平均降水量為911 mm，多年平均徑流深為231 mm，多年平均徑流系數為0.25，其徑流年內分布不均，主要集中在汛期。流域水系如圖1所示。研究數據來源于水文年鑒(2003—2014年)逐日逐時段實測降水、流量和蒸散發數據(逐時段蒸散發數據由逐日數據平均得到)。

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

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

2 分布式垂向混合產流模型構建

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

3 實時校正方法

3.1 ECM模型

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

yt=b0+b1xt+εt

(1)

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

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

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

(2)

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

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

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

(3)

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

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

Δ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時刻實測流量數據一階差分；ΔQsim,t-1為t-1時刻模擬流量數據一階差分。

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

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

(5)

式中b2為系數。

3.2 AR模型

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

(6)

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

3.3 ARECM模型

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

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

(7)

式中α0、α1為系數。

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

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

(8)

式中β0、β1、β2為系數。該模型顯示，t時刻誤差值et可由et-1和et-2估計得到。由于誤差序列可能是非平穩的，因此不能直接應用普通最小二乘法(ordinary least square，OLS)計算參數取值。根據式(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)/λ

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

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

(11)

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

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

(12)

3.4 評價指標

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

修正效果評價系數為

(13)

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

4 實時校正模型應用效果分析

4.1 平穩性檢驗

采用廣泛應用于時間序列平穩性檢驗的Augment Dickey-Fuller (ADF)檢驗法對10場洪水誤差系列進行平穩性檢驗。基于ADF檢驗的回歸方程為

(14)

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

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

表2 ADF檢驗結果Table 2 Results of ADF test

4.2 協整檢驗

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

表3 協整檢驗結果的P值Table 3 P values of cointegration test results

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

4.3 模型效果分析

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

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

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

表4 修正效果評價系數統計結果Table 4 Statistics of correction effect evaluation coefficient

表5 確定性系數統計結果Table 5 Statistics of deterministic coefficient

表6 洪峰相對誤差統計結果 單位：%Table 6 Statistics of relative error of flood peaks unit: %

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

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

表7 峰現時差統計結果 單位：hTable 7 Statistics of error of time to flood peaks unit: h

表8 徑流深相對誤差統計結果 單位：%Table 8 Statistics of relative error of runoff depth unit: %

(a) 20030620次洪水

5 結 論

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

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

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