基于貝葉斯方法的流量數(shù)據(jù)的不確定性分析*

顧西輝，張 強

(中山大學(xué) 水資源與環(huán)境系∥華南地區(qū)水循環(huán)與水安全廣東省普通高校重點實驗室，廣東 廣州 510275)

水文資料整編規(guī)范(2000)指出[1]：水位流量關(guān)系受到?jīng)_淤、變動回水、洪水漲落、水生植物、結(jié)冰等因素的影響，會產(chǎn)生不確定性。國內(nèi)外對于這種不確定性的分析進行了大量研究。程偉等[2]對三峽水庫蓄水后下游近壩段水位流量關(guān)系中的研究指出：河道侵蝕基點條件是影響水位流量關(guān)系的主要因素; 戴凌全等[3]用多次冪多項式函數(shù)作為水位流量關(guān)系的擬合公式，用最小二乘法進行擬合；夏軍強等[4]提出一維水動力模型對黃河下游平攤流量進行分析; 董曉華等[5]基于最小二乘法對繩套型水位流量關(guān)系曲線進行擬合。此外采用改進的BP神經(jīng)網(wǎng)絡(luò)、混合禁忌搜索算法、簡化的人工魚群算法等算法進行水位流量關(guān)系曲線的擬合[6-8]。文獻[9]對意大利Po河流域下游河段進行上述各個來源不確定性定量分析，得出流量總的不確定性在6.2%～42.8%之間；José-Luis Guerrero等[10]用廣義似然估計和蒙特卡洛方法分析水位流量關(guān)系的時間變異性。上述研究對于水位流量關(guān)系的時空變異和流量數(shù)據(jù)的不確定性大小及原因有著重要的科學(xué)意義。

但是以往的研究往往局限于對水位流量關(guān)系曲線的擬合方法、擬合線型的研究，并且在流量不確定性方面局限于沖淤、變動回水、洪水漲落、水生植物、結(jié)冰等因素對擬合精度的影響，對水位流量關(guān)系的變異性引起的不確定性研究較少，而水位流量關(guān)系變異是多種因素影響下的結(jié)果，同時也是流域水資源管理等諸多實踐應(yīng)用的基礎(chǔ)?；诖耍敬尾捎秘惾~斯方法來估計水位流量冪律關(guān)系式的參數(shù)，定量界定水位流量關(guān)系變異性引起的平均低流量、平均中流量和平均高流量的不確定性[11]。東江流域是珠江流域主要支流之一，主要擔(dān)負珠三角地區(qū)大城市，如廣州、東莞等，以及香港80%的用水，對于保障區(qū)域水資源安全與社會穩(wěn)定具有極其重要的作用和意義。研究東江流域流量不確定性對于東江流域水資源管理具有重要理論與現(xiàn)實意義。

1 研究區(qū)域和數(shù)據(jù)

圖1 河源、嶺下和博羅3個測站的地理分布圖Fig.1 The geographic distribution of three stationsHeyuan，Lingxia and Boluo

東江屬于珠江三大水系之一，位于中國南方濕潤地區(qū)，長度523 km，流域總面積35 636 km2，平均坡降為0.35‰。本文以東江流域河源、嶺下和博羅3個測站為研究對象(圖1)。采用河源站1981-2009年水位流量數(shù)據(jù)，嶺下站1956-2009年水位流量數(shù)據(jù)，博羅站1956-2005年水位流量數(shù)據(jù)進行研究。

2 研究方法

2.1 水位流量曲線擬合的貝葉斯公式

本文采用有著位置參數(shù)的冪律關(guān)系式，進行單一曲線法推流[1]:

(1)

公式中，Q是流量，h是水位。(A1,A2,…,An)是尺度參數(shù)，(b1,b2,…,bn)是形狀參數(shù)，h0是位置參數(shù)，代表河道橫斷面最低點的高程，(h1,h2,…,hn-1)是不同片段水位流量曲線的分割點，例如：河道的深槽、邊灘、低灘、高灘等由于水力性質(zhì)、徑流條件等因素的差異，公式(1)中的各種參數(shù)是不同的。由于各個片段在擬合時是相互獨立的，所以本文只考慮一條單一的水位流量曲線進行研究。文獻[12]以假設(shè)河道橫斷面寬度-深度呈指數(shù)關(guān)系作為基礎(chǔ)，結(jié)合曼寧公式推導(dǎo)上述滿足冪律關(guān)系的水位流量關(guān)系式，給出了公式(1)詳細的水力學(xué)方程式推導(dǎo)過程。單一片段的對數(shù)回歸模型為：

qi=a+blg(hi-h0)+εi

(2)

公式(2)中qi=lg(Qi)，a=lg(A)，為尺度參數(shù)A的自然對數(shù)值。εi為模型殘差，假設(shè)εi～N(0,σ2)。構(gòu)造向量θ=(a,b,h0)，根據(jù)貝葉斯概率公式，參數(shù)(θ,σ2)的后驗概率密度函數(shù)為：

f(θ,σ2|D)=

(3)

公式(3)中，D代表實測的水位-流量值，L是似然函數(shù)，π是先驗分布。

2.2 流量不確定性計算公式

(4)

3 分析結(jié)果

3.1 貝葉斯方法的擬合優(yōu)度檢驗

參數(shù)a和b分別是河道的尺度參數(shù)和形狀參數(shù)，共同為河道水力性質(zhì)和橫斷面幾何形狀所決定，所以給予一個二維聯(lián)合正態(tài)分布作為先驗分布。位置參數(shù)h0給予一個無信息先驗分布，噪聲(公式(2)對數(shù)回歸模型殘差)參數(shù)σ2給予一個逆伽馬分布。作者在另一篇文章詳細介紹了參數(shù)集(a,b,h0,σ2)的先驗分布、后驗分布以及推斷方法。

表1 貝葉斯方法擬合的參數(shù)值和擬合優(yōu)度指標

圖2 嶺下站2005年水位流量關(guān)系的貝葉斯擬合效果及置信區(qū)間Fig.2 The fitting effect and confidence interval of stage-discharge relationship fitted by Bayesian method in Lingxia station

3.2 水位流量關(guān)系的時間變異性

公式(2)中的尺度參數(shù)a、形狀參數(shù)b和位置參數(shù)h0反映出了河道水力特性和橫斷面幾何形狀的變化情況。用貝葉斯方法分別分析1981-2009年河源站, 1956-2009年嶺下站和1956-2005年博羅站逐年的水位-流量數(shù)據(jù)，研究尺度參數(shù)a、形狀參數(shù)b和位置參數(shù)h0的變化情況(圖3，表2)。

從圖3中可以看出，以1980年為分界點，嶺下站的尺度參數(shù)a呈下降趨勢，形狀參數(shù)b呈上升趨勢，博羅站正好相反。河源站在1980年之后，變化情況規(guī)律性不明顯，總體來說，尺度參數(shù)a和形狀參數(shù)b均呈下降趨勢。對于位置參數(shù)h0，3個測站均呈下降趨勢，這種趨勢在1980年之后表現(xiàn)的更加顯著。并且博羅站在3個測站中的下降趨勢最顯著。這些變化可能受到1980年代東江流域開始興起的大量采砂的影響。Luo等[14]指出：從1986-2003年珠江流域共采掘沙子超過8.7×108m3，其中導(dǎo)致東江流域河道平均下切1.77～6.48 m。河道采沙的另一后果是導(dǎo)致河道形狀向窄深方向發(fā)展，即河道寬深比下降[15]。參數(shù)h0持續(xù)下降反應(yīng)了東江流域河床下切的現(xiàn)象，參數(shù)a、b反映了河道形狀—寬深比的變化。

圖3 河源站、嶺下站和博羅站貝葉斯估計的系數(shù)歷年變化及趨勢圖Fig.3 The coefficients changing with time estimated by Bayesian method and their trends of Heyuan，Lingxia and Boluo stations

站點參數(shù)均值離勢系數(shù)最大差值河源a4.9520.1132.042b1.6220.191.467h029.6110.0354.44嶺下a4.0720.1772.903b2.1060.1331.177h011.4130.0592.615博羅a4.8950.1342.745b1.7730.1481.107h03.7130.3124.638

從河源、嶺下和博羅站各個參數(shù)變化情況來看，洪水對尺度參數(shù)a和形狀參數(shù)b的趨勢也起到了顯著的影響。比如1959年的洪水，嶺下站尺度參數(shù)a由下降趨勢變成上升趨勢，同樣的情況也出現(xiàn)在1966年洪水。1994、2005年2場洪水則把尺度參數(shù)a由上升趨勢變成下降趨勢。所以洪水對河道水力特性以及幾何形狀的影響比較復(fù)雜，變化沒有明顯統(tǒng)一的趨勢。

從表2中可以看得出來，地理位置相對靠近河口的嶺下站和博羅站比相對遠離河口的河源站，其離勢系數(shù)相對較大，這表明東江流域河道水力特性和幾何形狀的變化越靠近河口，變化越劇烈，其水位流量關(guān)系的不確定性也應(yīng)該相應(yīng)的增加。這種性質(zhì)也可以從參數(shù)的最大差值中進一步的反映出來。

3.3 水位流量關(guān)系的變異性引起的不確定性

東江流域河源、嶺下、博羅3個測站每一年的水位數(shù)據(jù)作為一個時間序列。使用這個時間序列的水位值，用貝葉斯推斷的水位流量關(guān)系曲線來估計流量值，同時計算相關(guān)的95%置信區(qū)間的上界和下界。公式(4)用近似的求和公式來代替，求取平均的低流量、中流量和高流量的不確定性(圖4)。圖4中可以看出3個測站1956-1974年，低流量、中流量和高流量的不確定性總體趨勢均在減小。珠江流域20世紀60-70年代大規(guī)模的聯(lián)圍筑閘、整治河道使河道形狀趨于穩(wěn)定，同時河道內(nèi)流量趨于穩(wěn)定，這一時期水位流量關(guān)系受到干擾較小，相對穩(wěn)定。90年代末期，低流量、中流量和高流量不確定性總體趨勢趨于穩(wěn)定(2009年除外，2009年珠江流域遇到大旱，流量的急劇減小是不確定性升高的主要原因)；河源站、嶺下站和博羅站，低流量、中流量和高流量的不確定性趨于吻合，個別年份有較大差異。主要因為90年代末期東江流域主要水利工程已經(jīng)修建完畢、河道大量采沙得到有效管理和治理，河道形狀趨于穩(wěn)定，并且東江流域3大水庫(楓樹壩、新豐江、白盆珠)對東江流域徑流量的有效調(diào)節(jié)，綜合影響下，東江流域低流量、中流量和高流量流量大小差距顯著變小，水位流量關(guān)系趨于穩(wěn)定。

70年代中期，嶺下站低流量不確定性有一個明顯增加的跳躍性變化。這一時期楓樹壩和新豐江兩座水庫的建成和運行使得嶺下站河道枯水期流量增加，低流量水位流量關(guān)系產(chǎn)生明顯異變。70年代中期到90年代初，嶺下站低、中、高流量普遍高于其他時期(圖5(b))。這一時期河道采沙大量增加，對河道形狀變化影響劇烈，水位流量關(guān)系非常不穩(wěn)定。

圖4 河源站、嶺下站和博羅站低流量、中流量和高流量不確定性Fig.4 The uncertainty in low flow、medium flow and high flow of Heyuan、Lingxia and Boluo stations

圖5 河源、嶺下和博羅3站歷年流量不確定性Fig.5 The flow uncertainty over the years of Heyuan,Lingxia and Boluo stations

單個測站來說，河源站、嶺下站和博羅站，表現(xiàn)的共同規(guī)律是：高流量不確定性最大，中流量不確定性次之，低流量不確定性最小(個別年份除外)。以博羅站為例(圖5c)，絕大部分年份，高流量要比中流量、低流量的不確定性要大，最高達到108%。由于博羅站位于東江流域的干流，上有秋香江等支流的匯流，因此水量充沛，洪峰流量曾經(jīng)達到12 800 m3/s(1959年)。高流量條件下，洪水對河道的沖刷、河道斷面幾何形狀以及河道水力特性的改變更加劇烈，相比低、中流量，高流量還涉及到河漫灘的形狀、水力性質(zhì)以及粗糙度的變化?；谏鲜鲈?，高流量不確定性相對較大。但是，圖5(c))顯示，70年代初期之后，高流量和低、中流量不確定性的差異顯著在減小。1974年新豐江水庫和楓樹壩水庫的建成和使用，對博羅站的消峰作用很明顯，有效的調(diào)節(jié)了博羅站高、低流量之間的差距。兩大水庫對穩(wěn)定博羅站流量的不確定性起到了很大作用。

河源站、嶺下站和博羅站每一年的水位-流量數(shù)據(jù)作為一個時間序列，共有133個時間序列。文獻[11]提出了一個質(zhì)量等級系統(tǒng)：不確定性在0～9%，10%～19%，20%～39%，40%～79%和>80%被分別評為優(yōu)、良、中、差和極差。統(tǒng)計這133個時間序列不確定性在各個等級中的比例(圖6)?？梢钥闯觯瑬|江流域低流量不確定性集中在0～20%，占到63%，表現(xiàn)最好。這一點與文獻[11]在挪威測站的研究是不同的：他認為挪威測站低流量不確定性表現(xiàn)最差。而東江流域高流量不確定性相對較大，集中在20%～40%，屬于中等水平。中流量不確定性居于低、高流量之間，不確定性低于40%的比例達到78%?？傮w來說，東江流域流量不確定性在中等偏上水平。

從均值和離勢系數(shù)來看(表3)，三個測站共同表現(xiàn)的規(guī)律是：流量越大，其不確定性的均值越大，離勢系數(shù)越小；同一類別的流量，越靠近河口，不確定性的均值相對較大，離勢系數(shù)相對較小。

圖6 133個時間序列流量不確定性各等級統(tǒng)計比例Fig.6 The statistical proportions of each grade of 133 time series flow uncertainty

表3 河源站、嶺下站和博羅站不確定性統(tǒng)計值Table 3 The statistical values of uncertainty in Heyuan，Lingxia and Boluo stations

這說明：同一個測站，低流量的不確定性相對較小，變化幅度相對較大；同一類別流量，越靠近河口，流量的不確定性相對較大，變化幅度相對較小。

4 結(jié) 論

本文通過用貝葉斯方法推斷水位流量冪律關(guān)系式的參數(shù)，繼而求出水位流量關(guān)系的變異性引起的不確定性，得出以下結(jié)論：

2)東江流域水位流量關(guān)系時間變異性顯著，河源站和博羅站河道尺度參數(shù)呈下降趨勢，形狀參數(shù)呈上升趨勢，嶺下站則相反。3個測站的位置參數(shù)均呈下降趨勢。3個測站參數(shù)的變化趨勢以1980年為分界點，前后變化顯著。

3)東江流域20世紀60-70年代大規(guī)模的聯(lián)圍筑閘、整治河道對減小3個測站的流量不確定性有顯著的作用；20世紀70年代中期到90年代初期，大規(guī)模的河道采沙對嶺下站流量不確定性影響最顯著。90年代之后，河源站和嶺下站流量不確定性在穩(wěn)定減少，博羅站流量不確定性呈增加趨勢；3個測站低、中、高流量不確定性之間的差距減小。

4)東江流域高流量不確定性相對較大，集中在20%～40%，屬于中等水平。中流量不確定性居于低、高流量之間，不確定性低于40%的比例達到78%?？傮w來說，東江流域流量不確定性在中等偏上水平。

[1] 水利部長江水利委員會水文局.水文資料整編規(guī)范[S],2000.

[2] 程偉,陳立,許文盛,等. 三峽水庫蓄水后下游近壩段水位流量關(guān)系[J].武漢大學(xué)學(xué)報：工學(xué)版,2011,44(4):434-444.

[3] 戴凌全,戴會超,蔣定國,等. 基于最小二乘法的河流水位流量關(guān)系曲線推算[J].人民黃河,2010,32(9):37-39.

[4] 夏軍強,吳保生,李文文. 黃河下游平灘流量不同確定方法的比較[J].泥沙研究,2009(3):20-29.

[5] 董曉華, 薄會娟, 鄧霞. 基于最小二乘法的繩套型水位流量關(guān)系最優(yōu)定線研究[J].中國農(nóng)村水利水電,2010(1):51-55.

[6] 李中志. 基于改進 BP 神經(jīng)網(wǎng)絡(luò)的水位流量關(guān)系擬合[J].中國農(nóng)村水利水電,2008(10):30-35.

[7] 汪嘉楊,李祚泳,倪長健,等. 基于混合禁忌搜索算法的水位流量關(guān)系擬合[J].系統(tǒng)工程,2006,24(6):107-110.

[8] 施秋紅, 王聯(lián)國. 基于簡化的人工魚群算法的水位流量關(guān)系擬合[J].甘肅農(nóng)業(yè)大學(xué)學(xué)報,2010,45(2):147-151.

[9] BALDASSARRE G D, MONTANARI A. Uncertainty in river discharge observations: a quantitative analysis[J].Hydrology and Earth System Sciences Discussions, 2009(6):39-61.

[10] GUERRERO J L, WESTERBERG I K, HALLDIN S， et al. Temporal variability in stage-discharge relationships[J].Journal of Hydrology, 2012, 446-447: 90-102.

[11] PETERSEN-?VERLEIR A, SOOT A, REITAN T. Bayesian rating curve inference as a streamflow data quality assessment tool[J].Water Resour Manage, 2009, 23(9):1835-1842.

[12] REITAN T, PETERSEN-?VERLEIR A. Dynamic rating curve assessment in unstable rivers using Ornstein-Uhlenbeck processes[J].Water Resources Research, 2011, 47(2):1-14.

[13] 曹深西，陳子燊.廣東沿海的極值風(fēng)速概率分布研究[J].海洋通報，2013,32(1):12-18.

[14] LUO X L, ZENG E Y, JI R Y, et al. Effects of in-channel sand excavation on the hydrology of the Pearl River Delta, China[J].Journal of Hydrology, 2007, 343(3/4):230-239.

[15] 季榮耀，陸永軍，左利欽. 東江下游博羅河段人類活動影響下的河床演變[J].泥沙研究, 2010(5):48-54.