富春江入庫控制斷面水位計(jì)算方法研究及誤差分析

陳學(xué)秋，瞿思敏,2，李代華，石 朋,2，王軼凡，單 帥，勾建峰

(1.河海大學(xué)水文水資源學(xué)院，南京 210098；2.河海大學(xué)水文水資源與水利工程國家重點(diǎn)實(shí)驗(yàn)室，南京 210098；3.云南省水文水資源局文山分局，云南 文山 663099)

富春江上游蘭溪入庫站水位受富春江水庫水位頂托、開機(jī)開閘回水波、動(dòng)庫容和上游橡膠壩等人類活動(dòng)影響，水位變化劇烈，傳統(tǒng)的水位流量關(guān)系線無法應(yīng)用。為應(yīng)用水位比降法準(zhǔn)確估計(jì)蘭溪站入庫流量，為洪水調(diào)度提供科學(xué)依據(jù)，首先需獲取蘭溪站及其下游女埠站水位資料。在調(diào)查過程中發(fā)現(xiàn)由于實(shí)際條件局限，無法實(shí)時(shí)獲得女埠站的水位數(shù)據(jù)，因此本文采用水動(dòng)力學(xué)方法和經(jīng)驗(yàn)相關(guān)法分別對(duì)女埠站水位進(jìn)行模擬計(jì)算。

以圣維南方程組為理論基礎(chǔ)的水動(dòng)力學(xué)方法是河道洪水演算的常用方法之一[1]。圣維南方程組是描述河道一維非恒定流水力要素隨時(shí)空變化的一階雙曲線型擬線性偏微分方程組[2]，隨著計(jì)算機(jī)技術(shù)的快速發(fā)展，高精度的數(shù)值求解方法被廣泛應(yīng)用于河道洪水模擬，并發(fā)展到多維河道水流[3]運(yùn)動(dòng)及河網(wǎng)計(jì)算[4]，以此為基礎(chǔ)建立的水動(dòng)力學(xué)模型具有一定物理基礎(chǔ)，對(duì)水文資料要求較低，計(jì)算精度高，適用于平原地區(qū)河道洪水演算。此外，水動(dòng)力學(xué)模型受到初始條件[5]、邊界條件以及河道糙率[6]等影響，對(duì)水文地形資料的選取有一定要求。經(jīng)驗(yàn)相關(guān)法是水位計(jì)算研究中出現(xiàn)最早、應(yīng)用最廣泛的一種數(shù)理統(tǒng)計(jì)方法，其基于大量實(shí)測(cè)水文資料的統(tǒng)計(jì)分析建立水文要素間的相關(guān)關(guān)系，根據(jù)已知站點(diǎn)水位過程推求未知站點(diǎn)水位過程。

1 研究方法與模型原理

1.1 水動(dòng)力學(xué)模型

圣維南方程組是描述河道一維非恒定水流運(yùn)動(dòng)的基本方程，它刻畫了各斷面水力要素隨時(shí)間、空間位置的變化規(guī)律。圣維南方程組依據(jù)質(zhì)量守恒和動(dòng)量守恒2大基本定律[7]，由連續(xù)方程和動(dòng)量方程組成：

(2)

式中：x為距離，m；t為時(shí)間，s；B為過水?dāng)嗝娴讓?，m；A為過水?dāng)嗝婷娣e，m2；Q為流量，m3/s；Z為水位，m；qL為旁側(cè)入流，m2/s，入流為正，出流為負(fù)；vx為旁側(cè)入流沿水流方向的速度，m/s，若旁側(cè)入流方向垂直于水流方向，則vx=0；α為動(dòng)量修正系數(shù)；g為重力加速度，m/s2；C為謝才系數(shù)，m1/2/s；R為水力半徑，m。

圣維南方程組是一組一階雙曲線型擬線性偏微分方程，本文擬采用簡(jiǎn)化四點(diǎn)線性隱格式對(duì)圣維南方程組進(jìn)行求解[8]，差分格式如下：

(3)

式中：f為變量，代表水位、流量等；θ為加權(quán)系數(shù)，0≤θ≤1。

1.2 經(jīng)驗(yàn)相關(guān)法

經(jīng)驗(yàn)相關(guān)法是水位計(jì)算研究中出現(xiàn)時(shí)間最早、應(yīng)用最廣泛的方法。該方法對(duì)大量實(shí)測(cè)資料進(jìn)行統(tǒng)計(jì)分析，找出具有相關(guān)關(guān)系的水文要素作為自變量和因變量，利用數(shù)理統(tǒng)計(jì)方法對(duì)其進(jìn)行擬合，得到各變量之間的經(jīng)驗(yàn)公式或相關(guān)關(guān)系圖[9]。經(jīng)驗(yàn)相關(guān)法依據(jù)變量多少分為單因素和多因素相關(guān)關(guān)系。依據(jù)自變量和因變量的關(guān)系分為線性和非線性相關(guān)關(guān)系。因此要根據(jù)研究河段影響水位變化的因素不同，因地制宜地選擇合適的水文要素?cái)M定經(jīng)驗(yàn)公式，進(jìn)行水位計(jì)算。

多因素線性相關(guān)基本模型一般表達(dá)為：

Y=b0+b1X1+b2X2+…+bnXn+ε

(4)

式中：X1,X2,…,Xn為自變量；Y為是因變量；b0,b1,…,bn為系數(shù)；ε為隨機(jī)誤差，均值為0。

2 模型應(yīng)用

2.1 流域概況

富春江水庫位于浙江省錢塘江中游富春江七里垅峽谷的出口處。錢塘江位于東經(jīng)118°～121°，北緯28°～31°，流域面積49 930 km2，是浙江省第1大河。富春江流域處于亞熱帶地區(qū)，氣候溫潤(rùn)多雨，多年平均氣溫為17.1 ℃。流域多年平均降水量為1 659.2 mm，最大年降水量為2 271.9 mm。流域降水量時(shí)空分配不均，年內(nèi)分配以6月份最多，12月份最少，空間分布西多東少。富春江流域山地占61.5%，丘陵占23.7%，平原占14.8%，地形自西向東傾斜，河流的坡降以上游較大，衢江以上約為0.5%，衢縣至梅城為0.25%～0.5%，梅城至壩址約為0.15%，壩址以下坡降平緩，已受錢塘江強(qiáng)潮影響。

近年來富春江流域下墊面因子發(fā)生了重大改變，流域內(nèi)部人類活動(dòng)影響增加，特別是流域內(nèi)小水電站較前幾年大量增加，使富春江流域洪水預(yù)報(bào)的難度大大增大。富春江上游有蘭溪入庫站，控制流域面積的85%，是富春江入庫估計(jì)有效的參照站。本文研究對(duì)象為蘭溪及其下游三河、梅城、江南站之間河段。

2.2 水動(dòng)力學(xué)模型應(yīng)用

2.2.1 模型構(gòu)建

蘭溪斷面為富春江水庫入庫控制斷面，取蘭江蘭溪至三河段為研究河段，河段長(zhǎng)21 km，其間女埠站距上游蘭溪站6 km，河段示意見圖1。

圖1 研究河段站點(diǎn)分布Fig.1 Distribution of hydrometric stations

河道首、末斷面采用85國家高程基準(zhǔn)面作為絕對(duì)基面，蘭溪站凍結(jié)基面與絕對(duì)基面高差為20.135 m，三河站凍結(jié)基面與絕對(duì)基面高差為20.096 m，2站相對(duì)高差為0.039 m。蘭溪站至三河站之間通過線性內(nèi)插得到20個(gè)斷面，各斷面之間相差1 km，計(jì)算斷面示意見圖2。內(nèi)插斷面及三河站無實(shí)測(cè)斷面資料，考慮到計(jì)算河段較短，且兩岸多為鋼筋混凝土防洪堤壩，斷面形狀變化不大，因此所有斷面均使用蘭溪站平均實(shí)測(cè)大斷面數(shù)據(jù)進(jìn)行計(jì)算。

圖2 研究河段計(jì)算Fig.2 Calculation structure of Lanxi-Sanhe river section

基于水位數(shù)據(jù)的一致性原則，本文選取2013年5場(chǎng)洪水過程，分別是20130323、20130423、20130601、20130626及20131007次洪水，其中前4場(chǎng)洪水用于參數(shù)率定，第5場(chǎng)洪水用于驗(yàn)證。

初始條件包括各斷面計(jì)算開始時(shí)刻水位或流量[5]。本文已知上下邊界初始時(shí)刻水位，對(duì)其進(jìn)行線性內(nèi)插可得各斷面初始水位。各斷面水位均取85國家高程基面。

根據(jù)已有水文資料，本文采用蘭溪站水位過程作為上邊界條件，三河站水位過程作為下邊界條件，中間不考慮旁側(cè)入流，計(jì)算女埠站水位過程。

2.2.2 參數(shù)率定及模型計(jì)算結(jié)果

對(duì)圣維南方程組采用4點(diǎn)隱式格式進(jìn)行差分，需要率定的參數(shù)包括河道糙率和權(quán)重系數(shù)。權(quán)重系數(shù)θ在[0.5,1]內(nèi)無條件穩(wěn)定，率定結(jié)果為0.75。考慮到河道糙率的空間變異性，模型對(duì)所有計(jì)算河段糙率逐一率定，率定結(jié)果見表1。

表1 各河段糙率率定結(jié)果Tab.1 Calibration results of reach roughness

參數(shù)率定及驗(yàn)證過程水位計(jì)算結(jié)果見表2，計(jì)算過程見圖3、圖4。5場(chǎng)洪水洪峰誤差均小于實(shí)測(cè)變幅的20%，峰現(xiàn)時(shí)間均不超過±1 h，確定性系數(shù)均大于0.8，說明此次模擬結(jié)果滿足許可誤差，精度較高。率定期洪水過程線出現(xiàn)計(jì)算洪峰水位大于實(shí)測(cè)洪峰水位，計(jì)算峰現(xiàn)時(shí)間晚于實(shí)測(cè)峰現(xiàn)時(shí)間，退水段洪水過程線擬合不佳的現(xiàn)象?？偨Y(jié)為以下2點(diǎn)：①漲水段與退水段模擬不一致，總體看來漲水段模擬結(jié)果比較符合實(shí)測(cè)值，退水段模擬結(jié)果偏高；②計(jì)算過程線相比實(shí)測(cè)過程出現(xiàn)坦化，不能反映真實(shí)水位的動(dòng)蕩變化，驗(yàn)證期洪水過程線表現(xiàn)尤為明顯。

圖3 率定期洪水模擬水位過程Fig.3 Simulated water level process of floods in model calibration

圖4 驗(yàn)證期洪水模擬水位過程(20131007次洪水)Fig.4 Simulated water level process of floods in model validation(flood 20131007)

類型洪號(hào)實(shí)測(cè)洪峰水位/m計(jì)算洪峰水位/m峰值絕對(duì)誤差/m峰值相對(duì)誤差/%峰現(xiàn)時(shí)間誤差/h確定性系數(shù)率定2013032323.8623.880.020.0800.942013042323.9924.010.020.08-10.922013060125.1125.160.050.1900.972013062625.3925.400.010.0400.98驗(yàn)證2013100723.3823.400.020.0900.88

注：峰現(xiàn)時(shí)間誤差一欄中，“-”代表洪峰滯后，“+”代表洪峰提前，以下表格同。

綜合分析水動(dòng)力學(xué)模型計(jì)算誤差來源主要是模型結(jié)構(gòu)的不足及斷面資料缺乏：①水動(dòng)力學(xué)模型用差分方程逼近微分方程時(shí)引入了誤差，誤差項(xiàng)使計(jì)算相位發(fā)生變化，導(dǎo)致洪峰滯后。②模型對(duì)河道底坡進(jìn)行簡(jiǎn)化，可能與河底真實(shí)底坡相差較大，如果內(nèi)插所得斷面底坡高于真實(shí)值，則會(huì)導(dǎo)致水深起算點(diǎn)比實(shí)際偏高，計(jì)算過程坡降偏小，重力作用被削減，流速減小，計(jì)算水位偏高。③糙率是河道洪水坦化的內(nèi)在因素[7]，漲退水過程中的水力要素往往不一致，本文模型概化中漲水和退水段采取同一套糙率值進(jìn)行計(jì)算，不能反映其時(shí)間變異性，導(dǎo)致退水過程計(jì)算誤差較大。

2.3 經(jīng)驗(yàn)?zāi)Ｐ蛻?yīng)用

2.3.1 斷面資料

本文通過分析江南水位站、梅城水位站、三河水位站和蘭溪水位站的水位資料，建立相關(guān)關(guān)系，推求女埠站水位資料。各水位站水文資料統(tǒng)計(jì)見表3，各水位站的地理位置見圖1。

2.3.2 經(jīng)驗(yàn)?zāi)Ｐ陀?jì)算

首先通過分析各水位站水位資料，找出女埠站與其他水位站之間的函數(shù)關(guān)系。通過20130323次、20130423次、20130601次及20130626次洪水實(shí)測(cè)資料率定得到女埠站計(jì)算水位表達(dá)式(回歸統(tǒng)計(jì)數(shù)據(jù)見表4，回歸參數(shù)分析見表5)：

表3 各水位站水文資料統(tǒng)計(jì)Tab.3 Hydrological data of water stations

Z(女)=0.697 908Z(蘭)+0.033 926Z(三)+

0.277 292Z(梅)+0.028 699Z(江)-1.007 58

(5)

表4 回歸數(shù)據(jù)統(tǒng)計(jì)Tab.4 Regression data statistical table

表5 回歸參數(shù)分析Tab.5 Regression parameters statistical table

注：表中Z(蘭)、Z(三)、Z(梅)、Z(江)分別表示蘭溪站、三河站、梅城站和江南站水位。

將20131007次洪水各水位站實(shí)測(cè)資料代入經(jīng)驗(yàn)公式計(jì)算得到女埠站計(jì)算水位過程，與實(shí)測(cè)水位過程進(jìn)行比較分析(見表6及圖5)。本文總結(jié)的經(jīng)驗(yàn)公式在進(jìn)行水位過程計(jì)算時(shí)模擬程度較好，確定性系數(shù)較高，洪峰誤差及峰現(xiàn)時(shí)間誤差皆在許可誤差范圍內(nèi)。由洪水過程線可知，漲水過程計(jì)算結(jié)果略低于實(shí)測(cè)過程，退水過程則相反，說明計(jì)算過程較實(shí)測(cè)過程較為滯后。

表6 經(jīng)驗(yàn)相關(guān)法計(jì)算結(jié)果Tab.6 Calculation results of the empirical correlation method

圖5 20131007次洪水計(jì)算過程Fig.5 Simulated water level process of flood 20131007

基于數(shù)據(jù)分析的經(jīng)驗(yàn)方法對(duì)資料的數(shù)量要求較大，本文只選取4場(chǎng)洪水率定經(jīng)驗(yàn)公式，1場(chǎng)洪水進(jìn)行驗(yàn)證，雖然計(jì)算結(jié)果總體較好，但是所得經(jīng)驗(yàn)公式是否具有普適性，能否滿足不同年份、不同場(chǎng)次洪水過程對(duì)水位計(jì)算的要求，還需要進(jìn)一步驗(yàn)證。

3 結(jié) 語

本文選擇水動(dòng)力學(xué)方法及經(jīng)驗(yàn)相關(guān)法計(jì)算富春江入庫控制斷面水位，2種方法計(jì)算結(jié)果精度均滿足要求。在洪峰水位計(jì)算上水動(dòng)力學(xué)方法優(yōu)于經(jīng)驗(yàn)相關(guān)法，而在洪水過程的模擬上則是經(jīng)驗(yàn)相關(guān)法的確定性系數(shù)更高?？傮w而言，本文建立的水動(dòng)力學(xué)方法具有一定物理基礎(chǔ)，對(duì)水文資料要求不高，能夠較好地對(duì)富春江入庫控制斷面水位進(jìn)行模擬計(jì)算，反映河道水流運(yùn)動(dòng)的內(nèi)在規(guī)律。如果模型概化能夠考慮糙率的時(shí)間變異性，則計(jì)算結(jié)果會(huì)更加準(zhǔn)確，因此值得深入研究。

□

[1] 芮孝芳.水文學(xué)原理[M]. 1版.北京:中國水利水電出版社,2004:6-7.

[2] M 喬芬,A 雷麥斯,芮孝芳. 圣維南方程組的解析解和數(shù)值解[J]. 河海大學(xué)科技情報(bào),1987,(2):155-163.

[3] 王船海,李光熾.流域洪水模擬[J].水利學(xué)報(bào),1996,(3):44-50.

[4] 李光熾.流域洪水演進(jìn)模型及其參數(shù)反問題研究[D]. 南京：河海大學(xué),2001.

[5] 張小琴,包為民,王 濤.水動(dòng)力學(xué)模型初始條件問題淺議[J].水電能源科學(xué),2008,26(2):62-65.

[6] 張小琴,包為民.糙率對(duì)感潮河段水位預(yù)報(bào)的影響[J].水文,2009,29(3):19-23.

[7] 趙振興,何建京.水力學(xué)[M].2版.北京:清華大學(xué)出版社,2010:108-109.

[8] 胡慶云,王船海.圣維南方程組4點(diǎn)線性隱格式的穩(wěn)定性分析[J].河海大學(xué)學(xué)報(bào)(自然科學(xué)版),2011,39(4):397-401.

[9] 包為民.水文預(yù)報(bào)[M].4版.北京：中國水利水電出版社，2009:132-134.