林博安

(廣東省水利電力勘測(cè)設(shè)計(jì)研究院有限公司，廣州 510635)

1 研究背景

水位流量關(guān)系作為水文資料整編的關(guān)鍵、河道水面線計(jì)算重要復(fù)核依據(jù)，其合理與否直接影響了水文資料整編的可靠性[1]。20世紀(jì)以來，受大規(guī)模人類活動(dòng)及上游來水來沙條件改變的影響，梅江干流沿程河道斷面變化劇烈[2-4]，以往水位流量關(guān)系曲線已不適用于現(xiàn)今河道斷面，且近些年梅江流量較小，無法根據(jù)實(shí)測(cè)資料確定水位流量關(guān)系曲線，需進(jìn)行斷面水位流量關(guān)系的延長(zhǎng)。

對(duì)斷面水位流量關(guān)系曲線的延長(zhǎng)，國(guó)內(nèi)外研究者已經(jīng)開展了大量工作，目前采用的方法大致可劃分為依靠河道水力學(xué)規(guī)律的公式法，不考慮河道水力學(xué)規(guī)律僅依靠水位流量關(guān)系曲線數(shù)學(xué)模型的方法，兩種方法應(yīng)用均較為廣泛。以往多數(shù)學(xué)者通過研究河道水力學(xué)規(guī)律，運(yùn)用水力學(xué)公式推算水位流量關(guān)系曲線，較為常用的是通過曼寧公式法與斯蒂文斯法[5-7]定線推流；門玉麗等[8]采用運(yùn)動(dòng)波方程及擴(kuò)散波方程來求解水位流量關(guān)系曲線，推出了水位流量關(guān)系雙向轉(zhuǎn)換的公式；袁帥等[9]為了解決欠缺實(shí)測(cè)水文控制斷面資料而無法進(jìn)行水位流量關(guān)系曲線經(jīng)驗(yàn)擬合的問題，提出了一種以圣維南方程為理論基礎(chǔ)的推導(dǎo)方法。

隨著數(shù)學(xué)的發(fā)展，越來越多優(yōu)化算法，使得水位流量關(guān)系曲線數(shù)學(xué)模型的參數(shù)優(yōu)化求解變得可能，越來越多學(xué)者通過水位流量關(guān)系曲線數(shù)學(xué)模型法直接擬合水位流量關(guān)系。康玲等[10]使用小波分析、舒棟才等[11]使用免疫進(jìn)化算法、李中志等[12]使用改進(jìn)BP神經(jīng)網(wǎng)絡(luò)算法擬合水位流量關(guān)系，均得到了良好擬合效果。

目前以上兩種方法均可應(yīng)用于河道斷面水位流量關(guān)系的擬合，但較少學(xué)者將兩種方法進(jìn)行性能比較?；诖?，本文選取梅江流域水口(二)站進(jìn)行分析，為其他水文站水位流量關(guān)系曲線延長(zhǎng)時(shí)方法選取提供參考。

2 方法介紹

2.1 水力學(xué)公式法

水力學(xué)公式法可以減少水位流量關(guān)系曲線的高水延長(zhǎng)中水位、流速順勢(shì)延長(zhǎng)的人為任意因素，本質(zhì)上是通過計(jì)算得出的流速值來輔助曲線延長(zhǎng)。

2.1.1曼寧公式法(Manning Formula，MF)

1889年美國(guó)水利學(xué)者曼寧提出曼寧公式，曼寧公式一般表達(dá)式為：

(1)

式中：

v——河道流速；

n——河道糙率；

R——水力半徑；

S——水面比降。

在對(duì)水位流量關(guān)系延長(zhǎng)時(shí)，通過實(shí)測(cè)資料計(jì)算得到流速，進(jìn)行水位面積關(guān)系曲線的延長(zhǎng)，從而實(shí)現(xiàn)水位流量關(guān)系的延長(zhǎng)。

因水力半徑R可通過斷面資料求得，故計(jì)算斷面流速關(guān)鍵在于確定河道水面比降和糙率。如水面比降S與糙率n均無法確定時(shí)，將n-1S1/2看成一個(gè)未知數(shù)，依據(jù)實(shí)測(cè)資料通過下式計(jì)算出n-1S1/2，點(diǎn)繪Z～n-1S1/2曲線，因高水部分n-1S1/2近似于一個(gè)常數(shù)，故可沿平行于縱軸方向延長(zhǎng)。

(2)

2.1.2斯蒂文斯法(Stevens，STE)

斯蒂文斯法利用謝才公式：

(3)

式中：

C——謝才系數(shù)(C=n-1R1/6)；

A——河道斷面過水面積；其余符號(hào)同曼寧公式。

2.2 水位流量數(shù)學(xué)模型法

天然河道里穩(wěn)定的水位流量關(guān)系可以抽象為一種非線性關(guān)系，概化后的數(shù)學(xué)模型形式可表示:

Q=a(H+b)c

(4)

式中:

Q——河道流量;

H——斷面水位;

a、b、c——待擬定的模型參數(shù)。

為了正確得到待擬定的模型參數(shù)，避免目標(biāo)函數(shù)陷入局部最優(yōu)解問題，本文引入模擬退火算法、遺傳算法及粒子群算法，綜合計(jì)算待擬定的模型參數(shù)，各算法簡(jiǎn)介如下。

2.2.1模擬退火算法(Simulated Annealing，SA)

模擬退火算法是由N.Metropolis等人于1953年提出，基于固體物質(zhì)的退火過程與一般組合優(yōu)化問題之間的相似性[13]。主要原理主要是通過高溫固體溫度的不斷下降，結(jié)合概率的突然跳躍特性，最終跳出問題的局部最優(yōu)解，求得全局最優(yōu)解。

2.2.2遺傳算法(Genetic Algorithm，GA)

遺傳算法是由Holland教授于1975年提出，起源于達(dá)爾文生物進(jìn)化論的自然選擇和遺傳學(xué)機(jī)理的生物進(jìn)化過程[14]。主要原理是通過解集的種群出發(fā)，不斷進(jìn)化，通過適者生存法則不斷產(chǎn)生出更優(yōu)解，每一代種群之間進(jìn)行組合交叉與變異，逐步適應(yīng)環(huán)境，進(jìn)而得到問題的最優(yōu)解。

2.2.3粒子群算法(Particle Swarm Optimization，PSO)

粒子群算法是在1995年由Eberhart和Kennedy博士提出，源于對(duì)鳥群捕食行為的研究[15]。主要原理是利用群體中的個(gè)體對(duì)信息的共享，使整個(gè)群體的運(yùn)動(dòng)在問題求解空間中產(chǎn)生從無序到有序的演化過程，從而獲得問題的最優(yōu)解。

a 水口(二)水位流量關(guān)系散點(diǎn)

3 案例分析

3.1 研究區(qū)概況

梅江是韓江的主干流，發(fā)源于廣東省河源市紫金縣七星崠，上游稱琴江，至五華縣水寨鎮(zhèn)與五華河匯合后稱梅江，至梅州市大埔縣與汀江匯合后稱韓江。梅江河長(zhǎng)為307 km，流域集雨面積為13 929 km2，自上而下建有近江、合江、龍上、三龍、西陽、丙村、單竹窩、蓬辣灘等8座梯級(jí)電站，水口(二)水文站位于合江及龍上電站之間。研究區(qū)域及水口(二)站位置見圖1。

圖1 研究區(qū)域及水口(二)位置示意

3.2 數(shù)據(jù)描述

本研究收集了梅江流域水口(二)站多年的實(shí)測(cè)洪水水文要素資料，資料可靠性強(qiáng)；由于本文研究結(jié)果考慮了人類活動(dòng)對(duì)水文要素的影響，因此，資料一致性可不用考慮；選取了水口(二)站多年資料，涵蓋了測(cè)站實(shí)測(cè)最大流量及最新資料，代表性強(qiáng)。

本次研究依據(jù)洪峰大小、年代遠(yuǎn)近等因素，分別選取了水口(二)站1986年、2013年以及2019年的河道大斷面及洪水要素?cái)?shù)據(jù)(見表1)，根據(jù)《韓江流域規(guī)劃》中水口(二)站設(shè)計(jì)流量成果(見表2)可知，1986年洪水約為30年一遇，2013年洪水接近30年一遇，2019年洪水不到5年一遇。為了能夠?qū)ρ娱L(zhǎng)的水位流量關(guān)系曲線進(jìn)行驗(yàn)證，因此，本次研究將各個(gè)年份流量較大值的序列進(jìn)行提取，作為未知數(shù)，對(duì)延長(zhǎng)的流量進(jìn)行實(shí)測(cè)數(shù)據(jù)對(duì)比，分析水位流量曲線延長(zhǎng)的效果，不同年份提取序列信息見表1，水口(二)站實(shí)測(cè)數(shù)據(jù)示意見圖2。

表1 水口(二)站不同年份數(shù)據(jù)主要信息 m3/s，m

表2 水口(二)站設(shè)計(jì)洪峰流量 m3/s

3.3 水力學(xué)公式法

表3 水口(二)站水位流量序列水力學(xué)公式法延長(zhǎng)結(jié)果

3.3.1曼寧公式法(Manning Formula，MF)

從計(jì)算的Z～n-1S1/2曲線(見圖3a)可知，在低水時(shí)n-1S1/2呈現(xiàn)波動(dòng)上升的趨勢(shì)，在高水時(shí)數(shù)值趨向平穩(wěn)。通過提取的高水序列與延長(zhǎng)后的序列的線性回歸統(tǒng)計(jì)參數(shù)可知，1986年、2013年模擬的高水序列與提取的高水序列基本一致，且呈現(xiàn)高水部分逐漸貼合的正?，F(xiàn)象，因此，可以認(rèn)為計(jì)算參數(shù)可用于水位流量關(guān)系曲線的延長(zhǎng)；2019年模擬的高水序列與提取的高水序列決定系數(shù)小于0，說明模擬的高水序列與提取的高水序列有很大不同。

a Z～n-1S1/2曲線

以上說明，若年內(nèi)觀測(cè)流量足夠大，使得計(jì)算得到的n-1S1/2位于平穩(wěn)區(qū)間，其值作為高水延長(zhǎng)的依據(jù)是可行的，并且能夠較為有效的反映現(xiàn)實(shí)高水部分水位流量關(guān)系；而若年內(nèi)觀測(cè)流量小，計(jì)算得到的n-1S1/2依舊位于低水不穩(wěn)定段，使用曼寧公式進(jìn)行水位流量關(guān)系的延長(zhǎng)，將會(huì)帶來較大誤差。

3.3.2斯蒂文斯法(Stevens，STE)

3.4 水位流量數(shù)學(xué)模型法

各方法擬合的水位流量關(guān)系曲線如圖4所示，結(jié)果見表4所示。

表4 水口(二)站水位流量數(shù)學(xué)模型法結(jié)果

圖4 水口(二)站水位流量數(shù)學(xué)模型法擬合效果示意

從模擬結(jié)果可知：

1)各方法擬合的數(shù)學(xué)公式結(jié)果說明，雖然各算法均能夠跳出局部最優(yōu)解去尋找全局最優(yōu)解，但針對(duì)3參數(shù)的優(yōu)化，每一個(gè)算法尋求的最優(yōu)解均不相同，得到的目標(biāo)函數(shù)值也均不同，而針對(duì)不同年份最優(yōu)算法也不同，因此，在運(yùn)用數(shù)學(xué)模型法進(jìn)行水位流量關(guān)系延長(zhǎng)時(shí)，需要綜合不同的算法進(jìn)行計(jì)算，定義合理的可能解取值域，對(duì)比分析最優(yōu)算法及最優(yōu)曲線。

2)1986年、2013年水位流量關(guān)系的延長(zhǎng)結(jié)果表明，通過水位流量數(shù)學(xué)模型法可以較好的延長(zhǎng)現(xiàn)有的水位流量關(guān)系，關(guān)聯(lián)度、均方根誤差及決定系數(shù)均表明模擬的高水序列與實(shí)測(cè)的高水序列基本一致，且相較于水力學(xué)公式，各年的最優(yōu)算法的統(tǒng)計(jì)參數(shù)均較優(yōu)，說明數(shù)學(xué)模型法在水位流量關(guān)系線性關(guān)系明顯時(shí)，計(jì)算結(jié)果優(yōu)于水力學(xué)公式法。

3)2019年結(jié)果表明，數(shù)學(xué)模型法用于擬合年內(nèi)流量較小的水位流量關(guān)系時(shí)，類似于斯蒂文斯法，能夠向上擬合一定區(qū)間內(nèi)的水位流量關(guān)系，但延長(zhǎng)較大的流量，就趨勢(shì)而言，均會(huì)陷入錯(cuò)誤的水位流量關(guān)系。

4)就模擬效果而言，不同的算法適用于不同年份的水位流量關(guān)系曲線，1986年最優(yōu)算法是粒子群算法(PSO)，2013年最優(yōu)算法是模擬退火算法(SA)，2019年最優(yōu)算法是粒子群算法。換言之，各算法在不同年份均不可避免的陷入局部最優(yōu)解的情況，以本文采用的計(jì)算代碼及計(jì)算案例而言，最優(yōu)算法應(yīng)為粒子群算法，在1986年、2019年目標(biāo)函數(shù)值最小，且雖然2013年目標(biāo)函數(shù)值不是最小，但模擬的高水序列為最接近實(shí)測(cè)的高水序列。

殘差絕對(duì)值和；Cov是數(shù)據(jù)關(guān)聯(lián)度；RMSE是均方根誤差；R2是決定系數(shù)。

3.5 綜合分析

通過擬合的水位流量關(guān)系曲線，求對(duì)應(yīng)設(shè)計(jì)頻率流量下的測(cè)站水位(見表5所示)可知：有實(shí)測(cè)高流量年份(1986年、2013年)的各方法求解的結(jié)果差別不大，100年一遇情況下，各方法求解結(jié)果最大差值僅為0.27 m；缺失實(shí)測(cè)高流量年份的(2019年)，各方法求解結(jié)果差別較大，100年一遇情況下各方法求解結(jié)果最大差值為2.21 m。

表5 水口(二)站水位流量延長(zhǎng)結(jié)果 m

將水力學(xué)公式法(曼寧公式法MF、斯蒂文斯法STE)擬合的曲線，水位流量數(shù)學(xué)模型法各優(yōu)化算法中擬合最好的曲線進(jìn)行比較(見圖5)。結(jié)果可知，針對(duì)擁有高水流量的年份(1986年、2013年)水力學(xué)公式法、數(shù)學(xué)模型法擬合效果均較好，從對(duì)提取的實(shí)測(cè)高水流量擬合結(jié)果而言，數(shù)學(xué)模型法對(duì)高水延長(zhǎng)效果更好；在水位流量關(guān)系曲線延長(zhǎng)方面，水力學(xué)公式法計(jì)算得到的對(duì)應(yīng)水位均高于數(shù)學(xué)模型法計(jì)算結(jié)果。當(dāng)年內(nèi)實(shí)測(cè)流量較小，如果僅通過小流量數(shù)據(jù)進(jìn)行延長(zhǎng)序列，將會(huì)陷入擬合較好但錯(cuò)誤的曲線中，此時(shí)應(yīng)通過曼寧公式法進(jìn)行高水序列的延長(zhǎng)，并結(jié)合往年水位流量曲線以及河道斷面變化進(jìn)行水位流量的矯正。

圖5 水口(二)站水位流量擬合效果示意

水力學(xué)公式法中的斯蒂文斯法及數(shù)學(xué)模型法涉及目標(biāo)函數(shù)的最值優(yōu)化問題，其中斯蒂文斯法為單參數(shù)的函數(shù)優(yōu)化求解，擬合過程較為簡(jiǎn)單；數(shù)學(xué)模型法涉及3參數(shù)的函數(shù)優(yōu)化求解，為避免陷入局部最優(yōu)解，可通過多種優(yōu)化算法進(jìn)行目標(biāo)函數(shù)最值求解。而對(duì)于曲線方程而言，目標(biāo)函數(shù)的確定則會(huì)很大程度的影響參數(shù)的優(yōu)化求解，本文采用高水水位、流量與擬合結(jié)果的殘差絕對(duì)值和最小作為目標(biāo)函數(shù)，通過提取的高水資料對(duì)比模擬的高水資料結(jié)果證明該目標(biāo)函數(shù)是可行的。但水位數(shù)學(xué)模型法涉及到的參數(shù)較多，雖然引入了優(yōu)化算法極大的減小了函數(shù)陷入局部最優(yōu)值的概率，但部分優(yōu)化算法進(jìn)行目標(biāo)函數(shù)的優(yōu)化求解時(shí)，依舊會(huì)陷入局部最優(yōu)解無法跳出，應(yīng)綜合多種優(yōu)化算法進(jìn)行計(jì)算，并且結(jié)合測(cè)站多年的參數(shù)，規(guī)定合理的取值范圍。

4 結(jié)語

本文選取了水力學(xué)公式法(曼寧公式法、斯蒂文斯法)及水位流量數(shù)學(xué)模型法，基于梅江水口(二)站3 a的實(shí)測(cè)水位流量數(shù)據(jù)，提取部分高水序列，對(duì)比分析了不同方法在測(cè)站水位流量關(guān)系曲線延長(zhǎng)方面的性能差異，進(jìn)一步探討針對(duì)不同來水條件下各方法的適用性。實(shí)例分析結(jié)果表明：

1)針對(duì)有實(shí)測(cè)高流量年份(1986年、2013年)的水位流量關(guān)系延長(zhǎng)，傳統(tǒng)水力學(xué)公式及水位流量數(shù)學(xué)模型法擬合效果均較好。傳統(tǒng)水力學(xué)公式方面，曼寧公式、斯蒂文斯法延長(zhǎng)的流量與實(shí)測(cè)流量的R2均大于0.98，相較于曼寧公式，斯蒂文斯法結(jié)果更優(yōu)；水位流量數(shù)學(xué)模型法方面，引入的SA、GA以及PSO優(yōu)化算法均能夠有效的優(yōu)化求解目標(biāo)函數(shù)，模擬的流量與實(shí)測(cè)流量R2最優(yōu)算法均大于0.99，其中PSO算法表現(xiàn)更為優(yōu)異。

2)針對(duì)有實(shí)測(cè)高流量年份的水位流量關(guān)系延長(zhǎng)，水位流量數(shù)學(xué)模型法相較于傳統(tǒng)水力學(xué)法而言，需要的數(shù)據(jù)量小且模擬效果更優(yōu)，且由于優(yōu)化算法的不斷發(fā)展，能夠較為有效的避免目標(biāo)函數(shù)陷入局部最優(yōu)值，針對(duì)3參數(shù)的目標(biāo)函數(shù)優(yōu)化求解變得容易，但研究結(jié)果證明：對(duì)3參數(shù)的優(yōu)化求解，依然需要引入多個(gè)算法進(jìn)行綜合對(duì)比求解。

3)針對(duì)沒有發(fā)生實(shí)測(cè)高流量年份的水位流量關(guān)系延長(zhǎng)，通過傳統(tǒng)水力學(xué)的斯蒂文斯法及水位流量數(shù)學(xué)模型法均容易陷入錯(cuò)誤水位流量關(guān)系，為確保延長(zhǎng)水位流量曲線，不出現(xiàn)明顯因水庫蓄水等人為因素影響低水水位流量關(guān)系而導(dǎo)致其余方法得到錯(cuò)誤水位流量曲線，此時(shí)，可通過曼寧公式法進(jìn)行水位流量曲線的延長(zhǎng)。