孫文宇，姚 成，劉志雨，何 健，黃鵬年

(1.河海大學(xué)水文水資源學(xué)院，江蘇 南京 210098； 2.水利部信息中心，北京 100053；3.江蘇省水文水資源勘測局，江蘇 南京 210029； 4.南京信息工程大學(xué)水文與水資源工程學(xué)院，江蘇 南京 210044)

中小河流是我國當(dāng)前洪水防控的重點(diǎn)薄弱環(huán)節(jié)[1-2]，建立流域水文水動(dòng)力模型及相應(yīng)的洪水預(yù)報(bào)與調(diào)度系統(tǒng)是流域防汛的重要非工程措施，有助于進(jìn)行必要的洪水模擬分析和評(píng)價(jià)，輔助防洪減災(zāi)決策方案的制定。秦淮河流域[3-6]有著封閉圩區(qū)的地形特征[7-8]和感潮河網(wǎng)[9-10]的水力條件，汛期排水不暢[11]，洪水多發(fā)，迫切需要建立充分考慮流域?qū)嶋H情況的水文水動(dòng)力模型，國內(nèi)外已有很多流域進(jìn)行了相應(yīng)的研究[12-18]。水文水動(dòng)力模型中涉及的相關(guān)參數(shù)通常由實(shí)測降雨徑流資料反求或根據(jù)流域特征進(jìn)行推算，只能反映流域降雨徑流過程的平均情況，用于實(shí)際預(yù)報(bào)時(shí)難以避免誤差；用于計(jì)算斷面徑流過程的實(shí)測資料精度有限，流域內(nèi)人為抽排水等因素干擾難以量化，存在系統(tǒng)性預(yù)報(bào)誤差。相較大型流域，各種誤差因素對(duì)中小河流徑流過程的影響更為顯著。實(shí)時(shí)校正[19]是對(duì)預(yù)報(bào)誤差的一種有效補(bǔ)救手段，可以減少不確定性帶來的誤差，提高預(yù)報(bào)精度，其中反饋法、K最近鄰算法(K-nearest neighbor algorithm，KNN法)運(yùn)用廣泛，如劉開磊等[20]將實(shí)時(shí)校正方法用于淮河流域水動(dòng)力學(xué)模型中，韓通等[21]將KNN法運(yùn)用于前期雨量指數(shù)模型來預(yù)報(bào)山區(qū)小流域洪水，徐杰等[22]等將多種實(shí)時(shí)校正方法用于半濕潤流域洪水預(yù)報(bào)，均取得較好的結(jié)果。

本文構(gòu)建適用于秦淮河流域的水文水動(dòng)力模型，并采用實(shí)時(shí)校正方法提高徑流模擬精度，可用于秦淮河流域洪水模擬預(yù)報(bào)，也可為其他中小河流流域防洪減災(zāi)工作提供參考。

1 模型構(gòu)建

秦淮河流域范圍、水系、主要站點(diǎn)及計(jì)算斷面情況見圖1。

圖1 秦淮河流域示意圖Fig.1 Map of Qinhuai River Basin

1.1 流域概況及資料收集

秦淮河流域位于長江下游南京段南岸，四周為丘陵山區(qū)，腹地為低洼圩區(qū)。水系總體流向自東南往西北，發(fā)源于溧水河、句容河，干流自三汊河匯入長江。秦淮河下游汛期洪水水位流量受長江潮位影響較大，長江南京下關(guān)潮位站實(shí)測潮位在夏季較高，是時(shí)秦淮河排洪受江潮頂托，排水不暢，易出現(xiàn)長期高水位現(xiàn)象。秦淮河流域內(nèi)有多處水庫、水閘等水利工程設(shè)施，下游干流屬感潮河段，降雨徑流及洪水演進(jìn)模擬相對(duì)復(fù)雜。

水文水動(dòng)力模型模擬所需的主要資料有流域下墊面資料和水文氣象資料兩大類。

采用地理信息數(shù)據(jù)云網(wǎng)站(www.gscloud.cn)的高分辨率數(shù)字高程資料，處理分析后可以得到流域范圍、高程、坡度、河道、下墊面類型等信息。結(jié)合水系概況、山丘區(qū)和圩區(qū)分布、水利工程分布，將秦淮河流域劃分為62個(gè)流域水文模型計(jì)算單元，再以水面、旱地、水稻田和城鎮(zhèn)道路4種下墊面類型，分別建立產(chǎn)匯流模型。

選用江蘇省水文水資源勘測局提供的2009—2017年水文氣象資料，主要包括降雨、蒸發(fā)、水位等實(shí)測數(shù)據(jù)。收集并整理了秦淮河流域范圍內(nèi)及周邊共23個(gè)雨量站日雨量數(shù)據(jù)，部分缺測數(shù)據(jù)采用臨近站替代；蒸發(fā)數(shù)據(jù)采用東山站、句容站實(shí)測資料，缺失部分根據(jù)流域蒸發(fā)特性(表1)估算；各水文站水位數(shù)據(jù)用于確定計(jì)算起始時(shí)刻的初始水位。

表1 秦淮河流域蒸發(fā)特性

1.2 河道一維水流模擬

用一維圣維南方程組描述河道漸變非恒定水流運(yùn)動(dòng)[23-24]：

(1)

式中：x、t分別為距離和時(shí)間坐標(biāo)；Q、A、B、Z分別為河道斷面的流量、過水面積、河寬和水位；qL、vx分別為旁側(cè)入流流量及其沿水流方向的流速分量；α為動(dòng)量校正系數(shù)；g為重力加速度；K為流量模數(shù)。

感潮河段水流運(yùn)動(dòng)同時(shí)受上游來水影響和下游潮位頂托，需采用動(dòng)力波方程描述，必須直接采用數(shù)值方法如雙追趕法求解圣維南方程組。

時(shí)段平均水力要素以時(shí)段初的水力要素替代，以首、末斷面水位為基本未知量，令：

(2)

式中：α、β、ζ、θ、η、γ為追趕系數(shù)；下標(biāo)i、j代表任意中間斷面，f、l代表首、末斷面。導(dǎo)出逆推及順推公式，由遞推公式可得首、末斷面流量與節(jié)點(diǎn)水位線性組合的方程組：

(3)

求得首、末節(jié)點(diǎn)水位后，由式(2)取同一中間斷面消去Qi，可求得任意節(jié)點(diǎn)水位：

(4)

得到Zi后，回代即可求得Qi。

1.3 邊界條件處理

1.3.1 流量上邊界計(jì)算

通過水文學(xué)方法，對(duì)降雨徑流過程進(jìn)行模擬，各支流起始斷面的計(jì)算流量作為流量上邊界。降雨徑流過程模擬包括產(chǎn)流計(jì)算和匯流計(jì)算兩部分，產(chǎn)流依據(jù)4種土地利用類型分別計(jì)算[25-27]。

a.水面，如湖泊水庫等大水體，凈雨深即為徑流深：

Rw=P-K1Em

(5)

式中：Rw為水面產(chǎn)生的徑流深；P為降水量；Em為蒸發(fā)皿觀測值；K1為大水體水面蒸發(fā)與蒸發(fā)皿觀測值的比率。

b.旱地的產(chǎn)流量用新安江模型中蓄滿產(chǎn)流[25]的方法計(jì)算，包氣帶土壤含水量達(dá)到田間持水量前不產(chǎn)流，達(dá)到田間持水量后產(chǎn)生的徑流深等于凈雨深。

c.水稻田產(chǎn)流考慮水稻不同生長階段的田間水深，每日進(jìn)行一次判斷，田間水深低于階段適宜水深下限則灌水至適宜水深，高于階段最大耐淹水深則排水至適宜水深上限，排水量為水稻田產(chǎn)流量。水稻生長期以外的時(shí)間水稻田作為旱地計(jì)算。

d.城鎮(zhèn)道路區(qū)域中不透水部分采用以下產(chǎn)流公式計(jì)算：

R1=KcP-Dc

(6)

式中：R1為不透水面產(chǎn)流量；Kc為徑流轉(zhuǎn)化系數(shù)，一般取0.8～0.95；Dc為不透水面填洼損失，一般取2 mm。

e.采用概化水庫法計(jì)算平原圩區(qū)匯流。平原圩區(qū)的河溝、塘壩有一定蓄水能力，概化成以圩區(qū)面積為底的平底水庫后，通過閘門、泵站調(diào)節(jié)蓄水量：

We=Wb+Rs-Wo

(7)

式中：Wb、We分別為概化庫區(qū)時(shí)段初、末蓄水量；Rs為概化庫區(qū)產(chǎn)流量；Wo為出庫水量；Wpl、Wys分別為排澇水量和引水量；Wmax、Wmin分別為概化庫區(qū)興利蓄水量的上、下限。

f.采用匯流曲線法計(jì)算平原非圩區(qū)匯流，匯流曲線可由實(shí)測資料率定或根據(jù)經(jīng)驗(yàn)曲線確定。

g.采用無因次單位線法計(jì)算山丘區(qū)子流域匯流，計(jì)算公式為

(8)

式中：Qt為t時(shí)刻出流；qt-k+1為t-k+1時(shí)刻調(diào)蓄入流；rk為k時(shí)刻無因次時(shí)段單位線的值，無因次時(shí)段單位線可通過綜合瞬時(shí)單位線法，基于流域下墊面特征和實(shí)測資料分析得出。

h.山丘區(qū)子流域內(nèi)的池塘和小水庫的調(diào)節(jié)作用采用一個(gè)虛擬水庫模擬，水量平衡方程為

(9)

式中：Vb、Ve分別為時(shí)段初、末水庫蓄水量，Qib、Qie分別為時(shí)段初、末水庫入流量；Qob、Qoe分別為時(shí)段初、末水庫出流量；QP為供水量。

1.3.2 水位下邊界處理

秦淮河于河定橋分流后，老秦淮河在下關(guān)潮位站處匯入長江，秦淮新河在下三山潮位站舊址處匯入長江。下三山潮位站在下關(guān)潮位站上游約14.5 km處，特征潮出現(xiàn)時(shí)間較下關(guān)潮位站滯后 30～45 min，僅有1951—1961年逐日特征潮位資料。

根據(jù)現(xiàn)有潮位資料，擬合下三山潮位站與下關(guān)潮位站特征潮位相關(guān)方程：

Zgs=1.055Zgg-0.212

(10)

Zds=1.054Zdg-0.070

(11)

式中：Zgs、Zds分別為下三山潮位站高、低潮位；Zgg、Zdg分別為下關(guān)潮位站高、低潮位。

1.4 節(jié)點(diǎn)堰閘流量計(jì)算

建立節(jié)點(diǎn)水位方程，考慮水量平衡，有：

(12)

式中：Qul為第l條河道流向節(jié)點(diǎn)u的流量；Au、Zu分別為節(jié)點(diǎn)u的蓄水面積和水位；L為匯入河道數(shù)。根據(jù)計(jì)算時(shí)段建立方程，用收斂迭代法求解即可。

堰閘過流有3種情形，關(guān)閘時(shí)，閘上、閘下作為兩條單一河道，閘上、下游河道斷面流量均為0；自由出流時(shí)，過閘流量不受閘下水位影響，閘上、閘下仍可視為兩條單一河道，不同在于河道斷面的流量為過閘流量：

(13)

(14)

式中：Qz為過閘流量；Zd為堰頂高程；Z1、Z2分別為堰上、下游水位；B為閘門開啟寬度；m為自由出流系數(shù)，取0.320～0.385；φ為淹沒出流系數(shù)，取1.0～1.08。

考慮連續(xù)性，堰上、下游流量均等于過閘流量。

2 實(shí)時(shí)校正方法

根據(jù)實(shí)測的水位或者流量資料，通過反饋機(jī)制修正預(yù)報(bào)值的過程稱為實(shí)時(shí)校正。中小河流的洪水模擬預(yù)報(bào)易受人為因素?cái)_動(dòng)，采用水文模型中常用的KNN法和反饋法對(duì)模擬結(jié)果進(jìn)行實(shí)時(shí)校正，以水文流量的連續(xù)性減少人為因素帶來的誤差，提高模擬精度。

2.1 KNN法

利用實(shí)測期內(nèi)各時(shí)刻預(yù)報(bào)誤差序列生成訓(xùn)練樣本庫，再基于誤差回歸分析的思想，假定t+e時(shí)刻，即當(dāng)前預(yù)報(bào)計(jì)算時(shí)刻的誤差值與實(shí)測期各時(shí)刻預(yù)報(bào)誤差值存在相關(guān)關(guān)系：

pt+e=f(pt，pt-1，…，pt-s+1)

(15)

式中：pt+e、pt，pt-1，…，pt-s+1為不同時(shí)刻預(yù)報(bào)誤差值；e為預(yù)見期；s為特征向量維數(shù)，即與計(jì)算時(shí)刻預(yù)報(bào)誤差值具有相關(guān)關(guān)系的預(yù)報(bào)誤差值個(gè)數(shù)。

KNN法可以看作一種歷史匹配方法，將當(dāng)前預(yù)報(bào)時(shí)刻的特征向量(pt，pt-1，…，pt-s+1)與訓(xùn)練樣本庫中各樣本的特征向量進(jìn)行比對(duì)，按歐氏距離從小到大來匹配近鄰樣本序列，用反距離權(quán)重法求得當(dāng)前預(yù)報(bào)時(shí)刻誤差估計(jì)值，校正值即為原預(yù)報(bào)結(jié)果加上預(yù)報(bào)誤差估計(jì)值。

2.2 反饋法

反饋法充分利用了最近時(shí)刻實(shí)測和預(yù)報(bào)值的誤差計(jì)算結(jié)果，先計(jì)算實(shí)測序列相鄰差值和預(yù)報(bào)序列相鄰差值，根據(jù)預(yù)報(bào)序列相鄰流量差值判斷漲水段及退水段，再進(jìn)行反饋校正計(jì)算。

a.當(dāng)流量序列處于漲水段時(shí)，實(shí)時(shí)校正流量為

(16)

式中：N為實(shí)測流量序列長度；n為序列號(hào)；Qfk、Qobs、Qcal分別為校正流量、實(shí)測流量和計(jì)算流量；F為反饋校正計(jì)算因子。

b.當(dāng)流量序列處于退水段時(shí)，直接進(jìn)行校正計(jì)算：

(17)

3 模擬及校正結(jié)果

選取秦淮河流域2009—2017年每年汛期典型洪水過程，采用構(gòu)建的秦淮河流域水文水動(dòng)力模型進(jìn)行模擬，計(jì)算時(shí)段取1 h。各計(jì)算斷面中，前垾村(秦)水文站斷面位于秦淮河干流上部，句容河與溧水河交匯處，具有多年汛期實(shí)測資料，可作為代表斷面，用于水文水動(dòng)力模型模擬及實(shí)時(shí)校正結(jié)果分析；武定門閘下斷面位于分流后老秦淮河下游，有實(shí)測日流量資料，可用于評(píng)價(jià)水文水動(dòng)力模型模擬效果。

選取洪峰相對(duì)誤差(δRPE)和納什效率系數(shù)(CNSE)用于分析模擬結(jié)果精度，洪峰相對(duì)誤差絕對(duì)值越小、納什效率系數(shù)越接近1則模擬效果越好。

3.1 代表斷面模型模擬結(jié)果

表2為前垾村(秦)和武定門閘下斷面實(shí)測洪峰流量(qobs,max)、計(jì)算洪峰流量(qcal,max)和洪峰相對(duì)誤差模擬結(jié)果。

表2 代表斷面徑流模擬結(jié)果

由表2可知，實(shí)時(shí)校正前，秦淮河流域水文水動(dòng)力模型的模擬結(jié)果較好，發(fā)生較大洪水時(shí)，模擬洪峰流量的精度較高。洪峰許可誤差取20%，9場洪水中，前垾村(秦)斷面有7場洪峰相對(duì)誤差絕對(duì)值小于20%，斷面洪峰預(yù)報(bào)合格率達(dá)78%，武定門閘下斷面合格率為67%。由于武定門閘下斷面位于干流分流后的河道，徑流量絕對(duì)值較小，易放大徑流模擬誤差，所以斷面模擬效果不如前垾村(秦)斷面。

3.2 預(yù)見期6 h模擬效果對(duì)比

以6 h預(yù)見期[22]為例，分析前垾村(秦)斷面實(shí)時(shí)校正前后的徑流模擬效果。圖2為部分年份前垾村(秦)斷面實(shí)測流量、模擬流量及兩種實(shí)時(shí)校正方法校正后的流量對(duì)比，圖3為實(shí)時(shí)校正前后各場次洪水模擬的洪峰相對(duì)誤差和納什效率系數(shù)分布。

圖2 前垾村(秦)斷面部分場次洪水模擬結(jié)果Fig.2 Partial flood simulation results of the Qianhan Village (Qin)

圖3 校正前后洪峰相對(duì)誤差、納什效率系數(shù)分布Fig.3 Distribution of flood peak error and Nash-Sutcliffe efficiency coefficient before and after real-time correction

實(shí)時(shí)校正前，在洪峰模擬較好的基礎(chǔ)上，9場洪水中有7場納什效率系數(shù)大于0.7，但實(shí)測洪峰陡漲陡落、峰形尖瘦的場次，模擬誤差相對(duì)偏大，主要原因可能在于原始資料，雨量數(shù)據(jù)為日雨量數(shù)據(jù)，插值后得到時(shí)段雨量數(shù)據(jù)，插值結(jié)果難以準(zhǔn)確匹配實(shí)際的降雨過程，坦化了短歷時(shí)高強(qiáng)度的暴雨過程，易放大漲水段流量及峰值流量的模擬誤差。人為因素的干擾也會(huì)對(duì)模擬精度帶來影響，模型對(duì)閘壩調(diào)度、農(nóng)田灌排、城區(qū)取用水過程的模擬可能存在偏差，易造成系統(tǒng)性誤差，且在洪水徑流較小的場次，模擬誤差影響更大。

經(jīng)由6 h預(yù)見期實(shí)時(shí)校正后，各場次洪水的模擬效果得到明顯提高，KNN法和反饋法校正結(jié)果中，斷面洪峰預(yù)報(bào)合格率達(dá)0.89，且9場洪水校正結(jié)果中KNN法納什效率系數(shù)均大于0.7，反饋法納什效率系數(shù)均大于0.8，反饋法表現(xiàn)略優(yōu)于KNN法，但KNN法能夠輸入更多參考樣本、調(diào)整預(yù)熱期、根據(jù)需求調(diào)節(jié)樣本向量維數(shù)或近鄰樣本個(gè)數(shù)，實(shí)際應(yīng)用更靈活。

3.3 不同預(yù)見期校正結(jié)果對(duì)比

依次選取2 h、4 h、6 h、8 h、12 h作為預(yù)見期，不同預(yù)見期內(nèi)兩種實(shí)時(shí)校正方法的洪峰相對(duì)誤差見圖4，納什效率系數(shù)見圖5。

圖4 不同預(yù)見期洪峰相對(duì)誤差分布Fig.4 Distribution of relative error of flood peak under different forecast periods

圖5 不同預(yù)見期納什效率系數(shù)分布Fig.5 Distribution of Nash-Sutcliffe efficiency coefficient under different forecast periods

隨著預(yù)見期減小，兩種實(shí)時(shí)校正方法洪峰相對(duì)誤差分布范圍逐漸減小，模擬精度明顯提高。在預(yù)見期4 h或更小時(shí)，反饋法所有場次洪水洪峰相對(duì)誤差絕對(duì)值均小于20%，KNN法所有場次洪水洪峰相對(duì)誤差絕對(duì)值均小于30%。但根據(jù)計(jì)算原理，預(yù)見期較大時(shí)，洪峰陡漲的場次模擬誤差會(huì)進(jìn)一步放大，實(shí)際校正結(jié)果也出現(xiàn)了洪峰相對(duì)誤差明顯放大的異常值?？紤]洪峰相對(duì)誤差分布情況，KNN法校正后，僅有一場洪水誤差較大，其余誤差分布較為集中，而反饋法誤差分布稍分散。分析不同預(yù)見期下納什效率系數(shù)的分布，兩種方法實(shí)時(shí)校正后，納什效率系數(shù)均明顯提高，預(yù)見期越小校正結(jié)果越佳，預(yù)見期小于8 h時(shí)所有場次洪水納什效率系數(shù)均能達(dá)到0.7以上。

以校正后洪峰相對(duì)誤差大小進(jìn)行比較，反饋法校正效果更好；以校正后納什效率系數(shù)結(jié)果進(jìn)行比較，兩種校正方法表現(xiàn)接近。KNN法根據(jù)實(shí)際需求，可以靈活調(diào)節(jié)計(jì)算參數(shù)以提高校正效果；反饋法充分依賴實(shí)測值及預(yù)報(bào)值，校正結(jié)果不易出現(xiàn)突變的異常值，但在預(yù)見期較大時(shí)，校正結(jié)果近似等于最新時(shí)刻實(shí)測值與預(yù)見期前后預(yù)報(bào)值差值之和，此時(shí)穩(wěn)定性下降，存在理論性不足的缺點(diǎn)。

4 結(jié) 論

a.秦淮河流域水文水動(dòng)力模型考慮了感潮河段、復(fù)雜平原河網(wǎng)的徑流模擬需要，模擬效果良好，雨洪較大的場次模擬效果更佳。

b.KNN法和反饋法兩種實(shí)時(shí)校正方法都能有效提高秦淮河流域水文水動(dòng)力模型的模擬精度，且隨著預(yù)見期減小，模擬精度越來越高。結(jié)合實(shí)時(shí)校正方法的秦淮河流域水文水動(dòng)力模型可應(yīng)用于流域的洪水預(yù)報(bào)、輔助分析決策等工作。