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

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

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

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

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

1 研究方法與模型原理

1.1 水動力學模型

圣維南方程組是描述河道一維非恒定水流運動的基本方程，它刻畫了各斷面水力要素隨時間、空間位置的變化規律。圣維南方程組依據質量守恒和動量守恒2大基本定律[7]，由連續方程和動量方程組成：

(2)

式中：x為距離，m；t為時間，s；B為過水斷面底寬，m；A為過水斷面面積，m2；Q為流量，m3/s；Z為水位，m；qL為旁側入流，m2/s，入流為正，出流為負；vx為旁側入流沿水流方向的速度，m/s，若旁側入流方向垂直于水流方向，則vx=0；α為動量修正系數；g為重力加速度，m/s2；C為謝才系數，m1/2/s；R為水力半徑，m。

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

(3)

式中：f為變量，代表水位、流量等；θ為加權系數，0≤θ≤1。

1.2 經驗相關法

經驗相關法是水位計算研究中出現時間最早、應用最廣泛的方法。該方法對大量實測資料進行統計分析，找出具有相關關系的水文要素作為自變量和因變量，利用數理統計方法對其進行擬合，得到各變量之間的經驗公式或相關關系圖[9]。經驗相關法依據變量多少分為單因素和多因素相關關系。依據自變量和因變量的關系分為線性和非線性相關關系。因此要根據研究河段影響水位變化的因素不同，因地制宜地選擇合適的水文要素擬定經驗公式，進行水位計算。

多因素線性相關基本模型一般表達為：

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

(4)

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

2 模型應用

2.1 流域概況

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

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

2.2 水動力學模型應用

2.2.1 模型構建

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

圖1 研究河段站點分布Fig.1 Distribution of hydrometric stations

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

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

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

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

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

2.2.2 參數率定及模型計算結果

對圣維南方程組采用4點隱式格式進行差分，需要率定的參數包括河道糙率和權重系數。權重系數θ在[0.5,1]內無條件穩定，率定結果為0.75。考慮到河道糙率的空間變異性，模型對所有計算河段糙率逐一率定，率定結果見表1。

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

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

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

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

類型洪號實測洪峰水位/m計算洪峰水位/m峰值絕對誤差/m峰值相對誤差/%峰現時間誤差/h確定性系數率定2013032323.8623.880.020.0800.942013042323.9924.010.020.08-10.922013060125.1125.160.050.1900.972013062625.3925.400.010.0400.98驗證2013100723.3823.400.020.0900.88

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

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

2.3 經驗模型應用

2.3.1 斷面資料

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

2.3.2 經驗模型計算

首先通過分析各水位站水位資料，找出女埠站與其他水位站之間的函數關系。通過20130323次、20130423次、20130601次及20130626次洪水實測資料率定得到女埠站計算水位表達式(回歸統計數據見表4，回歸參數分析見表5)：

表3 各水位站水文資料統計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 回歸數據統計Tab.4 Regression data statistical table

表5 回歸參數分析Tab.5 Regression parameters statistical table

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

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

表6 經驗相關法計算結果Tab.6 Calculation results of the empirical correlation method

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

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

3 結 語

本文選擇水動力學方法及經驗相關法計算富春江入庫控制斷面水位，2種方法計算結果精度均滿足要求。在洪峰水位計算上水動力學方法優于經驗相關法，而在洪水過程的模擬上則是經驗相關法的確定性系數更高。總體而言，本文建立的水動力學方法具有一定物理基礎，對水文資料要求不高，能夠較好地對富春江入庫控制斷面水位進行模擬計算，反映河道水流運動的內在規律。如果模型概化能夠考慮糙率的時間變異性，則計算結果會更加準確，因此值得深入研究。

□

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

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

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

[4] 李光熾.流域洪水演進模型及其參數反問題研究[D]. 南京：河海大學,2001.

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

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

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

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

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