周紹陽,王 悅,龔朝海

(長江水利委員會水文局 長江下游水文水資源勘測局,江蘇 南京 210011)

0 引 言

在水文測驗領(lǐng)域,相比水位觀測,流量測驗更為復(fù)雜且費時費力。水位和流量在客觀上存在顯著關(guān)聯(lián)性,可通過實時監(jiān)測水位推算斷面相應(yīng)流量。然而,在天然河道中,受洪水漲落、斷面沖淤、變動回水等因素影響,水位流量關(guān)系通常不具有單值函數(shù)的特性,呈不穩(wěn)定形態(tài),無法直接通過水位推算相應(yīng)流量,只能采用連時序等方法整編流量;而連時序法需要布置大量測次,工作量大。因此,為了減少流量測次,降低工作強度,有必要建立水位與流量的單值化關(guān)系[1-2]。

水位和流量單值化關(guān)系研究的主要方法包括校正因素法、抵償河長法、落差法等[3]。落差指數(shù)法是以落差法為基礎(chǔ),通過優(yōu)選落差指數(shù)β,建立水位與流量同落差的β次方之比例關(guān)系。該方法適用于斷面基本穩(wěn)定、受變動回水或變動回水及洪水漲落綜合影響的測站。自20世紀70年代以來,該方法逐步發(fā)展完善,已被廣泛應(yīng)用于水文站的資料整編[4]。在實際應(yīng)用中,落差指數(shù)和落差系數(shù)等參數(shù)的確定是影響落差指數(shù)法精度的主要原因之一[5]。

傳統(tǒng)參數(shù)確定方法通常是利用相關(guān)資料,通過試錯法或優(yōu)選法確定部分參數(shù),再通過點繪關(guān)系圖或試錯法求解其他參數(shù)。該方法存在以下不足之處:① 點繪關(guān)系圖工作量大,且受人工經(jīng)驗和技能水平影響,存在非唯一解的問題;② 各參數(shù)之間可能存在密切關(guān)聯(lián),傳統(tǒng)方法得到的參數(shù)不一定是最優(yōu)結(jié)果。近年來,隨著計算機技術(shù)和機器學(xué)習(xí)的快速發(fā)展,落差指數(shù)法的參數(shù)優(yōu)選問題得到更廣泛研究。李正最[6]提出了一種直接求解落差指數(shù)法參數(shù)的方法,該方法簡單易行,但計算量較大。喻娓厚等[7]則利用遺傳算法的全局尋優(yōu)能力對落差指數(shù)法中的參數(shù)進行優(yōu)化,該方法適用范圍廣,但實現(xiàn)較為困難,需要對遺傳算法有一定的理解。鄧才等[8]引入粒子群算法直接求解落差指數(shù)法中的參數(shù),該方法計算速度快,但算法實現(xiàn)較為困難,參數(shù)設(shè)置不當(dāng)可能會導(dǎo)致出現(xiàn)局部最優(yōu)解的情況。

本文研究了落差指數(shù)法原理,發(fā)現(xiàn)落差指數(shù)法的參數(shù)優(yōu)選可視為一個非線性約束優(yōu)化問題。廣義簡約梯度算法是求解非線性約束問題的最有效方法之一。相比其他優(yōu)化算法,該方法具有簡單實用、算法收斂快、計算精度高等優(yōu)點[9-10]。因此,本研究采用廣義簡約梯度算法對落差指數(shù)法中的參數(shù)進行優(yōu)選,以確定水位與流量的單值化關(guān)系。

1 研究方法

1.1 落差指數(shù)法

1972年,研究發(fā)現(xiàn)落差指數(shù)在0.5上下變動可明顯提高定線精度[3]。1972～1975年,在深入研究了落差指數(shù)和相應(yīng)的落差水尺優(yōu)選及誤差分析的基礎(chǔ)上,提出了落差指數(shù)法[4]。落差指數(shù)法可由曼寧公式導(dǎo)出[11]:

(1)

式中:q為單值化流量,m3/s;Q為實測流量,m3/s;β為綜合落差指數(shù);ΔZ為綜合落差,m,由式 (2) 計算:

ΔZ=α1ΔZ1+α2ΔZ2+…+αnΔZn

(2)

式中:ΔZ1,ΔZ2,…,ΔZn為與各參證站的同時落差,m;α1,α2,…,αn為落差系數(shù),且α1+α2+…+αn=1。

單值化流量與水位具有明確的函數(shù)關(guān)系,在實際水文整編工作中,常采用的函數(shù)關(guān)系主要有指數(shù)方程、對數(shù)方程和多項式方程[12],其中,多項式方程表達式為

q=k0+k1Z+k2Z2+…+kmZm

(3)

式中:k0,k1,k2,…,km為多項式系數(shù);Z為測站水位,m。

由式 (1)～ (3) 可知,整個求解過程中需要確定β,α1,α2,…,αn,k0,k1,k2,…,km等參數(shù)。利用傳統(tǒng)方法很難將所有參數(shù)一起優(yōu)化,一般先初定β值,再通過試錯法定出Z-q關(guān)系曲線,直至定線精度符合單一曲線的定線要求。而利用廣義簡約梯度算法可快速實現(xiàn)參數(shù)的全部優(yōu)化。

1.2 廣義簡約梯度算法

簡約梯度算法(Reduced Gradient Method,簡稱RG算法)是非線性規(guī)劃的單純形法應(yīng)用于僅具有線性約束的非線性約束規(guī)劃問題。后來,RG算法被推廣到非線性約束問題,得到廣義簡約梯度算法(Generalized Reduced Gradient Method,簡稱GRG算法)[13]。GRG算法是一種針對含有非線性約束的優(yōu)化問題的求解方法[14-15],優(yōu)化問題可簡化為

minF(X)X∈En

s.t.H(X)=0

(4)

L≤X≤UL,U∈En

式中:X為設(shè)計變量;F(X)為目標(biāo)函數(shù);En表示n維空間;H(X)為約束函數(shù),H(X)=[h1(X),h2(X),…,hm(X)]T;L為設(shè)計變量X的下限向量,L=[l1,l2,…,ln]T;U為設(shè)計變量X的上限向量,U=[u1,u2,…,un]T。

求解時,首先將X分為兩部分,即X=[XB,XN]T。其中,XB為基向量,m維;XN為非基向量,(n-m)維。相應(yīng)地,L=[LB,LN]T,U=[UB,UN]T。由隱函數(shù)存在定理可知,存在連續(xù)映射:

XB=V(XN)

(5)

式中:V(XN)為連續(xù)映射函數(shù),將XN的取值通過某種關(guān)系映射到XB的取值。

目標(biāo)函數(shù)F(X)可轉(zhuǎn)化為

F(X)=F(XB,XN)=f(XN)

(6)

式中:f(XN)為關(guān)于非基向量XN的目標(biāo)函數(shù)。

原來n個變量的目標(biāo)函數(shù)F(X)可以轉(zhuǎn)化為n-m個變量的函數(shù)f(XN),則f(XN)在Xk關(guān)于XN的梯度即為簡約梯度,可得到F(X)關(guān)于XN的簡約梯度為

?f(XNk)=?NF(Xk)-?NH(Xk)

[?BH(Xk)]-1?BF(Xk)

(7)

式中:?f(XNk)為簡約梯度,表示非基向量XN的目標(biāo)函數(shù)關(guān)于XN在第k步迭代中的梯度;?NF(Xk)為原始目標(biāo)函數(shù)關(guān)于非基向量XN在當(dāng)前步驟k的梯度;?NH(Xk)為原始約束函數(shù)關(guān)于非基向量XN在當(dāng)前步驟k的梯度;?BH(Xk)為原始約束函數(shù)關(guān)于基向量XB在當(dāng)前步驟k的梯度;?BF(Xk)為原始目標(biāo)函數(shù)關(guān)于基向量XB在當(dāng)前步驟k的梯度?？珊営洖?/p>

?f(XNk)=[r1,r2,…,rn-m]T

(8)

式中:r1,r2,…,rn-m表示簡約梯度向量?f(XNk)的分量。

定義Sk=[s1,s2,…,sn-m]T,其中:

(9)

令

(10)

Yc+1=Yc-H(Yc,XNk+1)[?BH(Y0,XNk+1)]-1

c=1,2,…

(11)

最后,求得滿足H(Yc+1,XNk+1)=0時的Yc+1,即可得到最優(yōu)值。

2 應(yīng)用案例

2.1 研究區(qū)域概況

大通水文站位于安徽省池州市貴池區(qū)梅龍街道,東經(jīng)117°37′,北緯30°46′,是國家一類水文站,也是長江下游干流最后一個徑流控制站。測驗河段上下游10 km范圍內(nèi)基本順直。上游有皖河、秋浦河等水匯入,下游有九華河、青弋江等匯入,均對斷面水流影響較小。影響水流的主要因素為上游干流來水和上游219 km處鄱陽湖出水。由于大通水文站位于潮區(qū)界點,距?？?24 km,枯季感潮顯著。多年來,該站水位流量關(guān)系雖然復(fù)雜,但總體穩(wěn)定。落差指數(shù)法在大通水文站資料整編中具有較好的適用性,單值化后的水位流量關(guān)系采用多項式函數(shù)效果更優(yōu)[12]。安慶水位站作為大通水文站上游參證站,位于安徽省安慶市沿江東路,距大通水文站約80 km。測站區(qū)位見圖1。

圖1 大通水文站區(qū)位Fig.1 Location of Datong Hydrological Station

2.2 數(shù)學(xué)建模

采用落差指數(shù)法進行參數(shù)優(yōu)選,主要精度評價指標(biāo)為實測點對關(guān)系線的隨機不確定度,計算公式為

(12)

(13)

結(jié)合大通水文站特性,公式 (2),(3),(1) 和(4) 可對應(yīng)改寫為下式 (14),(15),(16) 和(17)。

ΔZ=Z安慶-Z大通

(14)

q=k0+k1Z大通+k2Z大通2+k3Z大通3

(15)

Qci=q(ΔZ)β

=(k0+k1Z大通+k2Z大通2+k3Z大通3)·

(Z安慶-Z大通)β

(16)

(17)

式中:ΔZ為安慶與大通的落差水位,m;Z安慶為安慶站水位,m;Z大通為大通站水位,m;q為大通站單值化流量,m3/s;k0,k1,k2,k3為多項式系數(shù)。

2.3 算法實現(xiàn)

算法編程可以通過FrontlineSolvers、Matlab等商用軟件實現(xiàn)[16]?？紤]到大通水文站每年的實測數(shù)據(jù)只有40組左右,采用Excel規(guī)劃求解工具中GRG算法模塊實現(xiàn)參數(shù)優(yōu)選更加便捷,算法實現(xiàn)流程見圖2。節(jié)選大通水文站2018～2022年整編數(shù)據(jù)并代入數(shù)學(xué)模型計算,參數(shù)優(yōu)選結(jié)果見表1。

表1 參數(shù)優(yōu)選結(jié)果

圖2 算法實現(xiàn)流程Fig.2 Flow chart of implementation of algorithm

3 結(jié)果分析

3.1 成果檢驗

將表1中計算出的相關(guān)參數(shù)代入公式 (15),即求得2008～2022年大通水文站水位流量關(guān)系單值化曲線,見圖3。由圖3可知,關(guān)系曲線擬合良好,符號、適線、偏離數(shù)值檢驗合理。統(tǒng)計實測點流量與推算流量的系統(tǒng)誤差及隨機不確定度,統(tǒng)計結(jié)果見表2。系統(tǒng)誤差范圍為-0.13%～-0.02%,隨機不確定度范圍為2.84%～5.24%。根據(jù)SL/T 247-2020《水文資料整編規(guī)范》中水位流量關(guān)系定線精度指標(biāo)要求,采用水力因素法的一類精度的水文站,應(yīng)滿足系統(tǒng)誤差不超過±2%,隨機不確定度不大于10%。由此可見,基于GRG算法的定線精度滿足規(guī)范要求,且在典型年(2020年長江流域特大洪水,2022年長江流域極端干旱)的情況下,成果可靠。

表2 GRG算法與南方片軟件定線精度對比

3.2 精度對比

大通水文站現(xiàn)有落差指數(shù)法是基于南方片軟件實現(xiàn)的,經(jīng)驗落差指數(shù)則是根據(jù)1998～2012年共15 a的實測資料試算求得,β=0.6,且在2013～2022年流量資料整編中驗證了可靠性。表2為2018～2022年南方片軟件定線與GRG算法定線的精度統(tǒng)計。由表2可知,兩種方法均可滿足大通水文站定線要求,系統(tǒng)誤差無明顯區(qū)別,但GRG算法定線的隨機不確定度和最大相對誤差(除2020年)均小于南方片軟件定線?？傮w而言,基于GRG算法的定線精度更高。

4 結(jié) 論

(1) GRG算法計算速度快且易于實現(xiàn),是目前求解約束非線性最優(yōu)化問題有效的方法之一。本文將該方法應(yīng)用于大通水文站2018～2022年落差指數(shù)法參數(shù)優(yōu)選,并與傳統(tǒng)方法進行了對比。結(jié)果表明:基于GRG算法的落差指數(shù)法參數(shù)優(yōu)選是可行且實用的,不僅實現(xiàn)了參數(shù)的全部優(yōu)化,提高了定線精度,還簡化了傳統(tǒng)方法的工作步驟,且不需要人工干預(yù),為落差指數(shù)法的參數(shù)求解提供了一種新思路。

(2) 使用GRG算法的前提是落差指數(shù)法適用于測站。當(dāng)測站水流或斷面條件發(fā)生變化時,應(yīng)注意分析方法的適用性。比如,大通水文站2022年個別數(shù)據(jù)不滿足定線精度要求,其主要原因為該站受極端枯水及感潮等因素影響,水流受下游水位頂托影響嚴重,導(dǎo)致測站不符合落差指數(shù)法的適用條件,南方片軟件和GRG算法的定線誤差均較大。因此,遇到反常數(shù)據(jù)時應(yīng)結(jié)合測站特性及水情加以分析。