汪昊燃,王 容,黃鵬年,何 健,姚 成
(1.河海大學(xué)水文與水資源學(xué)院,江蘇 南京 210098; 2.水利部水文水資源監(jiān)測(cè)預(yù)報(bào)中心,北京 100053;3.南京信息工程大學(xué)水文與水資源工程學(xué)院,江蘇 南京 210044; 4.江蘇省水文水資源勘測(cè)局,江蘇 南京 210029)
中小河流流域的洪水預(yù)報(bào)一直是防汛工作的重點(diǎn)[1-3]。分布式水文模型如柵格新安江模型等,對(duì)山地丘陵地形時(shí)空模擬精細(xì)化程度較高,在匯流階段一般采用馬斯京根法或擴(kuò)散波計(jì)算匯流,計(jì)算結(jié)果受河道地形、回水頂托和閘壩建筑等影響,導(dǎo)致模擬精度降低[4-11]。水力學(xué)模型對(duì)感潮河段潮位計(jì)算具有良好的模擬精度,但上游山地丘陵地區(qū)多為天然河道,地形復(fù)雜,斷面形狀不規(guī)則,造成水力學(xué)相關(guān)參數(shù)計(jì)算困難[12],因此有必要將水文學(xué)模型與水力學(xué)模型相結(jié)合,以此來(lái)建立河流與流域之間的水力聯(lián)系。殷駿[13]在南四湖區(qū)對(duì)水文和一維水力耦合預(yù)報(bào)模型的研究明確了該模型的可行性和穩(wěn)定性;在水文學(xué)和水力學(xué)模型結(jié)合的基礎(chǔ)上,李致家等[14]對(duì)比了一維水力學(xué)模型和二維水力學(xué)模型的適用情況,發(fā)現(xiàn)一維水力學(xué)模型計(jì)算精度在低水位情況下受各水體之間連通性影響較大;孫映宏等[15]根據(jù)山區(qū)、平原區(qū)、分洪區(qū)以及河道等不同的產(chǎn)匯流條件將水文學(xué)與水力學(xué)模型耦合后應(yīng)用在山前小流域,取得了良好的模擬效果。一維水力學(xué)模型在低水位計(jì)算時(shí),由于各水體之間水流交換減少,導(dǎo)致模擬結(jié)果精度受到影響,需采用實(shí)時(shí)校正方法對(duì)誤差結(jié)果進(jìn)行修正[16]。實(shí)時(shí)校正方法一般分為過(guò)程誤差校正方法和終端誤差校正方法,在一維水力學(xué)模型的實(shí)時(shí)校正中,可采用卡爾曼濾波法對(duì)區(qū)間入流、狀態(tài)變量等過(guò)程進(jìn)行校正[17-19],也可采用K最近鄰(K-nearest neighbor,KNN)法等對(duì)計(jì)算結(jié)果進(jìn)行校正。KNN法[20]不需要建立自回歸方程,不需要計(jì)算相關(guān)參數(shù),在實(shí)時(shí)校正中有較好的應(yīng)用前景[21]。劉開(kāi)磊等[22]采用KNN法對(duì)一維水力學(xué)模型計(jì)算結(jié)果進(jìn)行校正,發(fā)現(xiàn)該方法能較大程度、穩(wěn)定地提高模擬精度;反饋模擬法對(duì)于較短的預(yù)見(jiàn)期有良好的校正效果,且與預(yù)報(bào)模型無(wú)關(guān),同樣是終端誤差實(shí)時(shí)校正的常用方法之一[23]。
秦淮河流域河網(wǎng)密集,上游以山區(qū)為主,坡度較大;下游河段出口與長(zhǎng)江相連,其水位受長(zhǎng)江潮汐作用影響顯著,屬于感潮河段且同時(shí)受水閘調(diào)度的影響。本文采用柵格新安江模型與一維河網(wǎng)水力學(xué)模型結(jié)合的模型模擬水位。在流域上游山地丘陵地區(qū)采用柵格新安江模型計(jì)算出流,流域下游平原感潮河段采用一維水力學(xué)模型計(jì)算水位,并對(duì)計(jì)算結(jié)果實(shí)時(shí)校正,以期提高洪水水位模擬的精度。
秦淮河流域位于長(zhǎng)江下游,地形四周高、中間低,坡度較大,總面積2684km2;上游屬山地丘陵區(qū),下游地區(qū)河道窄、坡度小的地形條件導(dǎo)致流域上游匯流快,下游泄洪緩慢[24]。
東山站位于南京市江寧區(qū)東山大橋上游200m處,是秦淮河流域下游出口的水位控制站。秦淮河在東山站下游300m處分流入老秦淮河及秦淮新河。秦淮新河全長(zhǎng)16.8km,秦淮新河閘位于秦淮新河河口處;武定門閘位于老秦淮河中段,與東山站相距12km,見(jiàn)圖1。武定門閘上水位站、秦淮新河閘上水位站主要為調(diào)度作用,東山站為主要負(fù)責(zé)報(bào)汛的水位控制站,因此選取東山站作為水位模擬計(jì)算的研究水位站。東山站水位受到上游來(lái)水、武定門閘、秦淮新河閘以及長(zhǎng)江潮汐的影響,汛期水位上漲迅速,因此水位呈現(xiàn)繩套或非單一性復(fù)雜關(guān)系。
柵格新安江模型產(chǎn)流部分采用蓄滿產(chǎn)流模型,計(jì)算單元柵格冠層截留量、河道降水量以及蒸散發(fā)量以推求地表徑流、壤中流和地下徑流,再將3種水源按各柵格間的演算次序演算至流域出口;模型匯流部分采用主網(wǎng)格馬斯京根法,將3種產(chǎn)流水源依次演算至站點(diǎn)出口[25]。整理2010—2018年?yáng)|山站逐時(shí)水位資料,2014—2018年秦淮河流域的逐日降雨、徑流資料,以秦淮河水系上游前垾村(秦)站為水文模型及水力學(xué)模型的分界點(diǎn),流域上游采用柵格新安江模型進(jìn)行模擬。DEM網(wǎng)格大小為1km×1km,流域內(nèi)共有9座水庫(kù),出流均按主網(wǎng)格馬斯京根法演算至前垾村(秦)站;前垾村上游包含11個(gè)雨量站,均在流域范圍內(nèi)。
秦淮河干流采用一維河網(wǎng)水力學(xué)模型進(jìn)行計(jì)算。模型的上邊界采用由柵格新安江模型計(jì)算的流量過(guò)程,以流域下游出口處的長(zhǎng)江潮位作為河網(wǎng)水力學(xué)模型的下邊界。
前垾村(秦)站下游包括秦淮新河閘上、武定門閘上、東山站3個(gè)水位站,云臺(tái)山河、牛首山河和外港河3條主要支流;秦淮河干流在東山站分為人工河道秦淮新河和老秦淮河兩條河流匯入長(zhǎng)江(圖2)。前垾村(秦)站至下游出口共設(shè)有86個(gè)斷面,其中秦淮新河糙率取0.022,老秦淮河糙率取0.024,其余支流糙率取0.021。秦淮新河閘與武定門閘為調(diào)度水閘,調(diào)度規(guī)則均為當(dāng)水位低于6.5m時(shí)引水,當(dāng)水位高于8.5m時(shí)泄洪。前垾村(秦)站以下流域地勢(shì)低洼,存在多個(gè)圩區(qū),將圩區(qū)分為水面、水稻田、旱荒地和城市道路進(jìn)行計(jì)算;產(chǎn)流計(jì)算將土壤分為兩層,采用蓄滿產(chǎn)流方法計(jì)算;圩區(qū)內(nèi)采用單位線法計(jì)算匯流,匯流至圩區(qū)河道后,以圩區(qū)排澇模數(shù)進(jìn)行量化,作為旁側(cè)入流進(jìn)入秦淮河河道內(nèi)。
圖2 河網(wǎng)概化圖
采用的河網(wǎng)水力學(xué)方程為圣維南方程組:
(1)
(2)
式中:A為過(guò)水?dāng)嗝婷娣e;Q為流量;u為流速;Z為斷面水位;g為重力加速度;S0為水流沿程損失;Sf為摩阻比降。將式(1)用四點(diǎn)線性Preissmann隱式格式進(jìn)行差分,可以得到任一河段差分方程:
(3)
式中Cj、Dj、Ej、Fj、Gj、φj均為系數(shù),由初值計(jì)算。
河網(wǎng)水力學(xué)模型將河網(wǎng)分為樹(shù)狀河網(wǎng)和環(huán)狀河網(wǎng)。前垾村(秦)站至東山站之間的河道及匯入河道的支流呈樹(shù)狀河網(wǎng)水系,采用樹(shù)狀河網(wǎng)追趕法進(jìn)行河道水位計(jì)算,即對(duì)于已知的流量邊界條件,上斷面L1和下斷面L2可假設(shè)如下追趕關(guān)系:
(4)
式中:j為斷面編號(hào)(j=L1,L1+1,L1+2,…,L2-1);K、T、P、V為追趕系數(shù)。由式(2)和式(3)聯(lián)立可得遞推關(guān)系,依次求得Kj+1、Tj+1、Pj+1、Vj+1,最后得到:
QL2=PL2-VL2ZL2
(5)
式(5)與下邊界條件QL2=f(ZL2)聯(lián)解可得ZL2,依次回代可求得Zj、Qj。
環(huán)狀河網(wǎng)包括內(nèi)河道和外河道,外河道計(jì)算與樹(shù)狀河網(wǎng)計(jì)算類似,均由河道一端邊界條件已知的外節(jié)點(diǎn)推求另一端水力要素未知的內(nèi)節(jié)點(diǎn)。東山站以下的老秦淮河、秦淮新河構(gòu)成環(huán)狀河網(wǎng)水系,采用環(huán)狀河網(wǎng)雙追趕系數(shù)法進(jìn)行水位計(jì)算。內(nèi)河道的計(jì)算沒(méi)有端點(diǎn)邊界條件可用,河段差分方程依然可表示為式(2),首、末斷面水位為基本未知量,因此需要利用雙追趕方程求解。
KNN法校對(duì)計(jì)算主要包括4個(gè)步驟:①計(jì)算模擬時(shí)段前70%的誤差數(shù)據(jù)存放為全部歷史誤差數(shù)據(jù),用于以后時(shí)段的實(shí)時(shí)校正;②根據(jù)誤差自回歸理論,考慮前面多個(gè)時(shí)段(NT),設(shè)定預(yù)見(jiàn)期為ext,則可以從NT+ext時(shí)刻開(kāi)始至m-1+NT+ext建立m個(gè)歷史預(yù)報(bào)誤差向量集;③識(shí)別最近鄰樣本,以向量間的歐式距離為評(píng)價(jià)相似相關(guān)性的指標(biāo),找出與當(dāng)前預(yù)報(bào)誤差相關(guān)向量最相似的X個(gè)預(yù)報(bào)誤差向量;④對(duì)該X個(gè)預(yù)報(bào)誤差向量對(duì)應(yīng)的預(yù)報(bào)誤差加權(quán),使用反距離權(quán)重法對(duì)當(dāng)前誤差進(jìn)行估計(jì)。
反饋法基于洪水預(yù)報(bào)誤差的相似性,對(duì)洪水預(yù)報(bào)誤差進(jìn)行實(shí)時(shí)校正,在洪水預(yù)報(bào)實(shí)時(shí)校正方面已取得較好的結(jié)果[26]。反饋法通過(guò)計(jì)算實(shí)測(cè)流量與預(yù)報(bào)流量之間的反饋因子(F),將流量過(guò)程分為漲洪段和落洪段,進(jìn)而對(duì)預(yù)報(bào)流量進(jìn)行修正。
先從第2個(gè)時(shí)刻起計(jì)算各時(shí)刻與前一時(shí)刻的流量差:
ΔQobs,i=Qobs,i-Qobs,i-1
(6)
ΔQcal,i=Qcal,i-Qcal,i-1
(7)
式中:Qobs,i為本時(shí)刻實(shí)測(cè)流量;Qcal,i為本時(shí)刻預(yù)報(bào)流量。
通過(guò)各時(shí)刻與前一時(shí)刻的預(yù)報(bào)流量差和實(shí)測(cè)流量差,計(jì)算F:
(8)
可以得到漲水段反饋模擬實(shí)時(shí)校正計(jì)算公式為
Qobs,i=Qobs,i-1+ΔQcal,iF
(9)
退水段反饋模擬實(shí)時(shí)校正計(jì)算公式為
(10)
選取秦淮河流域2014—2018年9場(chǎng)典型洪水,采用柵格新安江模型計(jì)算前垾村(秦)站上游至該站的產(chǎn)匯流過(guò)程,選定前8場(chǎng)進(jìn)行次洪模型參數(shù)率定,最后1場(chǎng)用來(lái)檢驗(yàn)。部分分布參數(shù)由輸入的土壤類型、植被覆蓋類型等數(shù)據(jù)直接計(jì)算,如自由水蓄水容量SM、張力水蓄水容量WM等。部分敏感參數(shù)率定結(jié)果為:蒸散發(fā)折算系數(shù)1.00;上層張力水容量占比0.17;下層張力水容量占比0.50;地下水消退系數(shù)0.993;壤中流消退系數(shù)0.90;河網(wǎng)水流消退系數(shù)0.83。洪水模擬結(jié)果如表1所示。
表1 前垾村(秦)站洪水模擬成果統(tǒng)計(jì)
根據(jù)GB/T 22482—2008《水文情報(bào)預(yù)報(bào)規(guī)范》:模擬計(jì)算結(jié)果中,前垾村(秦)站以上流域8場(chǎng)洪水,徑流深合格率與洪峰合格率均為87.5%;預(yù)報(bào)檢驗(yàn)的1場(chǎng)洪水徑流深和洪峰都合格。所選9場(chǎng)洪水的徑流深誤差與洪峰誤差都達(dá)到乙等精度要求。由柵格新安江模型計(jì)算的前垾村(秦)站洪水流量過(guò)程可用于一維水力學(xué)模型上邊界流量輸入過(guò)程。
選取次洪模擬中的5場(chǎng)洪水進(jìn)行水位模擬計(jì)算。將前垾村(秦)站的模擬流量過(guò)程作為一維水力學(xué)模型上邊界流量輸入條件,下游出口處長(zhǎng)江潮位作為一維水力學(xué)模型的下邊界,并采用KNN法和反饋法進(jìn)行實(shí)時(shí)校正。對(duì)模擬水位的評(píng)價(jià)標(biāo)準(zhǔn)主要考慮其與實(shí)測(cè)水位的絕對(duì)誤差,因此對(duì)東山站水位計(jì)算結(jié)果的評(píng)價(jià)采用確定性系數(shù)(CNSE)和洪峰水位誤差(δPE)2種指標(biāo)。對(duì)未校正的柵格新安江-河網(wǎng)水力學(xué)模型計(jì)算成果、反饋模擬法校正的計(jì)算成果和KNN法校正的計(jì)算成果分別計(jì)算兩種指標(biāo),計(jì)算結(jié)果如圖3所示。部分模擬水位過(guò)程如圖4所示。
圖3 校正方法指標(biāo)對(duì)比
圖4 模擬水位過(guò)程
由圖3中可以看出:未校正的水位模擬計(jì)算結(jié)果存在較大的不確定性,經(jīng)過(guò)反饋法和KNN法校正后的模擬水位過(guò)程確定性系數(shù)均有較大提升,且KNN法校正效果相比反饋法更好。對(duì)于洪峰水位絕對(duì)誤差,未校正的模擬洪峰水位計(jì)算結(jié)果波動(dòng)較大,經(jīng)過(guò)校正后精度得到提高,且KNN法校正后的模擬洪峰水位誤差更加穩(wěn)定,誤差值在0.04~0.55m之間,模擬精度較高。2015062508號(hào)洪水未校正的模擬水位確定性系數(shù)為0.388,該場(chǎng)洪水的上邊界前垾村(秦)站的入流模擬結(jié)果較好,但經(jīng)過(guò)河網(wǎng)水力學(xué)模型計(jì)算后,東山站的模擬水位較實(shí)測(cè)水位相差較大,造成確定性系數(shù)偏小。經(jīng)過(guò)調(diào)查,本場(chǎng)洪水在南京市區(qū)造成嚴(yán)重內(nèi)澇,因此在三汊河口處人工開(kāi)鑿臨時(shí)泄洪道,對(duì)下邊界條件產(chǎn)生影響,導(dǎo)致高水位時(shí)模擬結(jié)果與實(shí)測(cè)結(jié)果誤差較大。同時(shí)根據(jù)調(diào)度規(guī)則,下游兩處水位站會(huì)在水位高于8.5m時(shí)泄洪,因此模擬水位在8.5m左右時(shí)會(huì)出現(xiàn)模擬調(diào)度實(shí)際未調(diào)度或模擬未調(diào)度實(shí)際調(diào)度的情況,造成在低水位時(shí)模擬結(jié)果不佳。經(jīng)過(guò)實(shí)時(shí)校正后,2015062508號(hào)洪水洪峰/水位模擬結(jié)果和確定性系數(shù)均大幅度提高,可以滿足模擬精度要求。
2017060920號(hào)洪水未校正的模擬水位誤差只有0.13m,但確定性系數(shù)僅為0.322。同時(shí)根據(jù)水位過(guò)程線可知,由于模擬的水位過(guò)程漲水階段相較于實(shí)測(cè)水位過(guò)程滯后約12 h,造成確定性系數(shù)偏小。分析前垾村(秦)站入流情況可知,在洪水發(fā)生初始階段該站流量為負(fù)值,即出現(xiàn)逆流情況,而柵格新安江模型模擬流量過(guò)程不會(huì)出現(xiàn)逆流情況,因此上邊界條件出現(xiàn)較大誤差,造成模擬水位過(guò)程與實(shí)際水位過(guò)程在時(shí)間尺度出現(xiàn)偏移。經(jīng)過(guò)KNN法校正后的水位模擬過(guò)程在最高水位時(shí)出現(xiàn)了一些精度損失,而峰現(xiàn)時(shí)差以及確定性系數(shù)有了較大的提高。
在實(shí)時(shí)校正前,柵格新安江模型與河網(wǎng)水力學(xué)模型結(jié)合的預(yù)報(bào)方法對(duì)水位的模擬計(jì)算結(jié)果受平原圩區(qū)地形影響較大,且因模型的局限性導(dǎo)致在低水位時(shí)模擬結(jié)果不佳,洪峰水位模擬精度較高但峰現(xiàn)時(shí)間有較大誤差。而對(duì)模擬結(jié)果采用KNN法和反饋法進(jìn)行實(shí)時(shí)校正后,均能有效提高東山站水位模擬精度。對(duì)比2種實(shí)時(shí)校正方法,反饋法計(jì)算簡(jiǎn)單,KNN法可以根據(jù)實(shí)際預(yù)報(bào)需要靈活調(diào)整參數(shù),且校正后的水位模擬結(jié)果相較于反饋法精度更高。
秦淮河流域地勢(shì)四周高中間低,形成了許多圩區(qū),導(dǎo)致上游匯流快,下游泄洪慢,且受水閘泵站和長(zhǎng)江潮汐的影響,洪水預(yù)報(bào)工作具有不確定性和復(fù)雜性。本文構(gòu)建以柵格新安江模型計(jì)算的前垾村(秦)站流量過(guò)程為一維河網(wǎng)水力學(xué)模型上邊界、以下游出口處長(zhǎng)江潮位為模型下邊界的水位模擬計(jì)算模型,選取秦淮河流域2014—2017年的5場(chǎng)次洪水進(jìn)行水位模擬計(jì)算。結(jié)果表明,河網(wǎng)水力學(xué)模型對(duì)平原圩區(qū)河道匯流計(jì)算具有較高的精度,結(jié)合柵格新安江模型的預(yù)報(bào)徑流過(guò)程作為上邊界輸入條件后,實(shí)現(xiàn)了對(duì)受水工建筑物影響的感潮河段水位的模擬計(jì)算。在經(jīng)過(guò)實(shí)時(shí)校正之后,解決了低水位情況下水位預(yù)報(bào)精度不高的問(wèn)題,其中KNN法校正后確定性系數(shù)為0.718~0.975,達(dá)到了乙等預(yù)報(bào)精度要求,可為秦淮河流域防汛提供參考。由于本文所選洪水場(chǎng)次較少,還需要在收集更多資料后對(duì)該模型進(jìn)行進(jìn)一步計(jì)算驗(yàn)證和校正,逐步形成預(yù)報(bào)可靠性較高的秦淮河流域柵格新安江與河網(wǎng)水力學(xué)結(jié)合的預(yù)報(bào)模型。