基于貝葉斯方法的流量數據的不確定性分析*

顧西輝，張 強

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

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

但是以往的研究往往局限于對水位流量關系曲線的擬合方法、擬合線型的研究，并且在流量不確定性方面局限于沖淤、變動回水、洪水漲落、水生植物、結冰等因素對擬合精度的影響，對水位流量關系的變異性引起的不確定性研究較少，而水位流量關系變異是多種因素影響下的結果，同時也是流域水資源管理等諸多實踐應用的基礎。基于此，本次采用貝葉斯方法來估計水位流量冪律關系式的參數，定量界定水位流量關系變異性引起的平均低流量、平均中流量和平均高流量的不確定性[11]。東江流域是珠江流域主要支流之一，主要擔負珠三角地區大城市，如廣州、東莞等，以及香港80%的用水，對于保障區域水資源安全與社會穩定具有極其重要的作用和意義。研究東江流域流量不確定性對于東江流域水資源管理具有重要理論與現實意義。

1 研究區域和數據

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

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

2 研究方法

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

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

(1)

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

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

(2)

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

f(θ,σ2|D)=

(3)

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

2.2 流量不確定性計算公式

(4)

3 分析結果

3.1 貝葉斯方法的擬合優度檢驗

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

表1 貝葉斯方法擬合的參數值和擬合優度指標

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

3.2 水位流量關系的時間變異性

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

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

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

站點參數均值離勢系數最大差值河源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

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

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

3.3 水位流量關系的變異性引起的不確定性

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

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

圖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

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

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

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

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

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

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

4 結 論

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

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

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

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

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

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

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

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

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

[6] 李中志. 基于改進 BP 神經網絡的水位流量關系擬合[J].中國農村水利水電,2008(10):30-35.

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

[8] 施秋紅, 王聯國. 基于簡化的人工魚群算法的水位流量關系擬合[J].甘肅農業大學學報,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] 曹深西，陳子燊.廣東沿海的極值風速概率分布研究[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.