水文變異條件下潦河生態流量計算研究

王 強， 夏 瑞， 鄒 磊， 陳 焰， 張 遠， 劉成建,4， 陳明昊,5

1.中國環境科學研究院水生態保護修復研究室, 北京 100012 2.武漢大學, 水資源與水電工程科學國家重點實驗室, 湖北 武漢 430072 3.中國科學院地理科學與資源研究所, 陸地水循環及地表過程重點實驗室, 北京 100101 4.西北大學城市與環境學院, 陜西省地表系統與環境承載力重點實驗室, 陜西 西安 710127 5.中國科學院重慶綠色智能技術研究院三峽生態環境研究所, 重慶 400714

氣候變化和高強度人類活動加速了流域水文循環過程，進而改變了降水和地表水資源的時空分布，導致徑流序列的一致性遭到了破壞[1-3]. 地表水文循環過程的顯著變化，改變了流域產匯流過程的原有屬性，導致徑流序列出現變異，給流域生態系統帶來嚴重影響. 此外，隨著社會經濟的發展和人類日益增長的用水需求，高強度的人類活動(如跨流域調水、水利工程建設、城市化建設等)對流域水文過程的影響將進一步加劇，給變化環境下水資源的適應性管理帶來了嚴峻挑戰[4-6]. 面對水資源量時空分布不均、水質污染嚴重、水生態環境日益惡化等問題，科學、合理估算河流生態流量將為保障河流生態系統健康良性發展、促進流域生態建設和水資源合理配置提供重要依據[7-9].

自20世紀七八十年代以來，生態流量的相關研究開始蓬勃發展. 據不完全統計，全球生態流量的計算方法超過200種，大致可分為4類，包括利用歷史流量資料的水文學方法(如Tennant法[10]、Texas法、NGPRP法等)、根據河道斷面水力特性的水力學法(如濕周法[11]、R2CROSS法等)、以指示物種所需的水力條件為依據的棲息地法(如WUW法、IFIM法[12]等)以及綜合考慮研究區生態環境整體的綜合法(如BBM法[13]等). 雖然各種方法都有各自的生態學理論基礎，但均有一定的適用條件和不足，目前針對河道生態流量的計算還并沒有統一的標準，各種方法仍在不斷發展中，不同方法的適用性也存在較大的爭議.

由于水文資料易于獲取、不需要開展現場測定數據，能夠利用長序列徑流資料分析統計特性，且具有較強的通用性，因此，水文學方法在全世界范圍內得到了廣泛的推廣和應用. 如朱才榮等[14]采用9種水文學方法估算了漢江中下游基本生態流量，提出了表征生態需水動態特征的生徑比概念，對不同水文學方法生態流量進行生徑比分析；李紫妍等[15]采用變化范圍法(RVA)、多年日流量資料排頻法等分析了漢江子午河生態流量過程，研究該河道維持多種生態環境目標所需的生態流量；龍凡等[16]通過建立年月豐枯遭遇的Copula聯合分布，求得不同典型年條件下各月豐平枯的條件概率，利用改進的FDC法計算了豐平枯典型年的生態流量過程；王俊釵等[17]以生徑比為標準，運用流速法、月最小徑流法、逐月頻率設定法選取動態生態流量過程；張濤等[18]采用多項式各時段的流量頻率曲線，建立流量與頻率的函數關系，以某枯水期頻率計算該時段河道的最小生態流量. 但上述研究均未充分考慮水文序列的一致性條件是否發生變化. 隨后，國內外部分學者開展水文變異條件下的生態需水研究，如肖才榮等[19]根據秩和檢驗法分析東江水文變異，并以變異點前各月平均流量序列進行分析推求生態流量；李劍鋒等[20]采用滑動秩和檢驗法分析黃河干流水文變異，并對河道內生態流量進行研究；劉劍宇等[21]分析了鄱陽湖流域的水文變異及其成因，并估算了鄱陽湖流域五河的生態流量；孟鈺等[22]分析了淮河干流的水文變異，以典型魚類長吻鮠為保護對象，采用改進的FLOWS法對水文變異前后長吻鮠的生態環境進行對比分析.

由于受到氣候變化和人類活動的影響，流域自然水文氣象序列往往容易發生突變，從而導致水文序列的時空變異，顯著影響當地生態系統完整性與河流健康狀態. 然而，目前針對水文變異條件下的生態流量研究大多是針對變異點前水文序列進行分析，但往往變異點前水文序列數據有限，限制了以長序列資料為數據基礎的水文學方法應用. 鄱陽湖是我國第一大淡水湖，與洞庭湖并稱“長江雙腎”，對長江中游水環境和水生態保護具有十分重要的作用. 在氣候變化和強人類活動作用下，鄱陽湖流域的管理和保護愈發重要，是當前長江保護修復的關鍵環節之一. 基于此，該研究以江西省鄱陽湖西北部的潦河為例，通過實測水文數據分析水文序列的變異性，提出了基于水文變異條件下河道內生態流量的計算方法，旨在為探索潦河生態需水標準的制定提供數據參考和技術支撐.

1 研究區域與數據來源

潦河流域位于江西省鄱陽湖區西北部，為修水的最大支流，又稱為奉新江. 潦河流域以萬家埠水文站位出口點，地理位置為114°87′E～115°38′E、28°38′N～29°04′N，流域面積 3 548 km2. 潦河流域地勢西高東低，中游為高山叢林地貌，植被良好，下游為丘陵崗地，植被稀少，水土流失嚴重. 流域內降水量較為充沛，屬于半濕潤地區，當地產生暴雨的天氣系統主要是臺風、高空槽、低壓冷鋒等，特大暴雨往往由多種天氣系統配合形成，多年平均年降水量為 1 600 mm，但受地形影響，降雨分布極不均勻，上游山區最大年降雨量記錄超過 2 000 mm. 潦河流域屬于暴雨多發區，為江西省四大暴雨中心之一，加上河床坡降較大，洪水匯流速度快，因其特殊的地理位置、氣候因素和洪水特點，洪災頻繁發生.

該研究收集了潦河流域及其周邊11個氣象站1963—2014年的逐日氣象觀測數據，站點分布較為均勻，能夠基本反映流域氣象水文變量的時空變化，部分缺測數據通過與相鄰氣象站點相對應的氣象變量序列建立回歸關系進行插補. 徑流數據來自于萬家埠水文站1953—2014年的水文要素摘錄表，包括日均流量和月均流量，用于對模型進行參數率定和模型驗證. 水文模型采用的數字高程模型(DEM)、土地覆蓋和土壤類型數據均來自中科院資源環境科學數據中心(http:www.resdc.cn)，流域位置及站點分布如圖1所示.

圖1 潦河流域位置及站點分布Fig.1 The location of Liaohe River Basin and stations

2 研究方法

2.1 水文變異診斷方法

徑流變異點檢測方法眾多，常見的有基于貝葉斯理論的里海哈林法、Mann-Whitney非參數檢驗法、有序聚類法、RS分析法、F檢驗法、T檢驗法、滑動秩和檢驗法和Mann-Kendall非參數統計檢驗法、Pettitt突變點檢驗法、SNHT(Standard Normal Homogeneity Test)方法、Buishand檢驗法和Buishand秩和突變檢驗法等. 單一檢測方法在徑流突變點檢驗中存在局限性，因此該研究采用Mann-Kendall非參數統計檢驗法、Pettitt突變點檢驗法、SNHT方法、Buishand檢驗法和Buishand秩和突變檢驗法對突變點進行綜合評估. 其中，Mann-Kendall非參數統計檢驗法所需樣本無需遵從預先假定的分布，且不受少數異常值的干擾，計算簡便實用，廣泛應用于水文氣象時間序列的趨勢和突變檢驗. 該方法具體計算步驟：對于具有n個樣本的水文氣象時間序列X，其觀察值序列為{xi,i=1,2,3,…,n}，假設該時間序列相互獨立且有相同連續分布，則按照時間序列順序和逆序構建統計量UF和UB.

(1)

(2)

其中，

(3)

式中，xi為第i個水文氣象變量，n為序列長度，Sk為秩序列，E(Sk)為方差，VarSk為期望，UF和UB為統計量. 通過進一步分析統計量序列UFk和UBk可以獲取時間序列的趨勢變化情況，可以明確突變時間.

2.2 分布式時變增益模型

由于氣候變化和高強度人類活動的影響，目前水文站實測徑流過程對于表征天然徑流情況存在較大局限，因此需要對測站以上受人類活動影響的水量進行還原，即為天然徑流還原. 目前，采用天然徑流還原方法有多種，其中應用最為廣泛的是分項調查法、蒸發差值法、降雨徑流關系法和水文模型法等[23-24]. 其中，水文模型是對自然界中復雜降雨徑流水文循環過程的概化和近似描述，以流域產匯流理論為基礎，可以綜合考慮氣象因素、地形、土地利用變化、土壤類型等的空間不均勻性，具有較為完善的物理基礎，且隨著模型軟件的不斷成熟，其獲取方式更為便捷，應用更加廣泛.

該研究采用分布式時變增益模型(Distributed Time-Variant Gain Model, DTVGM)[25-27]在徑流突變檢驗的基礎上進行徑流還原. DTVGM模型是將XIA等[28-30]提出的水文非線性時變增益模型(TVGM)通過DEMGIS平臺進行拓展，進而推廣到流域水文變量時空模擬的分布式水文模型. 該模型通過建立不同下墊面狀況(土地利用、土壤類型等)與水文非線性產流參數之間的關系，一方面具有分布式模型概念性模擬的特點，另一方面具有水文非線性系統分析適應能力強的優點，該模型已經在國內外諸多流域水文模擬方面進行了應用和驗證. 針對分布式模型水文模擬，首先根據潦河流域數字高程模型(Digital Elevation Model，DEM)對流域進行子流域劃分，進而選取兩種評價指標——Nash-Sutchliffe效率系數(NSE)[31]和相關系數(r)作為評估模型性能的評價指標對模型進行參數率定和驗證. 其中，NSE反映了觀測值和對應模擬值之間的擬合程度，可用式(4)計算：

(4)

(5)

2.3 生態需水計算方法

該研究主要采用逐月頻率計算法估算生態流量，逐月頻率計算法即對多年各月的流量序列進行排頻和適線分析，一般取90%的保證率作為河道內生態流量. 該方法一方面可以考慮河流生態系統對水量的要求，另一方面還能夠根據不同月份徑流的變化特征，考慮不同月份生態環境需水的不同要求. Green等[32]研究指出，任何一種概率分布函數只可能對某種分布或子樣本容量的檢驗效果比較好，并不存在一種占絕對優勢的概率分布函數. 因此，該研究采用兩參數正態分布(ND)、伽瑪分布(GD)，三參數廣義極值分布(GEV)、皮爾遜Ⅲ型分布(P-Ⅲ)、對數正態分布(LN3)函數對經模型還原后的各月徑流序列進行分析，各分布詳細介紹如表1所示. 為了尋求最適宜于描述月徑流過程的概率分布線型，該研究采用Kolmogorov-Smirnov(K-S)、Anderson Darling(A-D)和概率點據相關系數(PPCC)3種檢驗法[33]進行擬合優度檢驗分析. 其中，PPCC檢驗法是假設待檢驗樣本服從某種分布，排序后的觀測值與假設的理論分布值的概率點據呈近似線型相關，相關系數r即可用來比較不同假設分布函數線型的優劣.

表1 該研究采用的概率分布函數

3 結果與討論

3.1 突變檢驗與成因分析

采用Mann-Kendall非參數統計檢驗法在顯著性水平α=0.05下對潦河萬家埠站1953—2014年年徑流時間序列進行趨勢突變分析，結果(見圖2)表明，萬家埠站徑流在1972年發生了突變. 突變前，萬家埠站徑流在1953—1972年間呈顯著下降趨勢，1972年后徑流量迅速增加. 突變點前后多年平均徑流深分別為910和1 020 mm，相比突變點前，突變點后的流域年徑流量增加了12%. 采用Pettitt突變點檢驗法、SNHT方法、Buishand檢驗法和Buishand秩和突變檢驗法對萬家埠站水文序列進行突變點診斷，均檢測出1972年為突變點，進一步增加了突變點存在的可信度，該結果與肖麗英等[34]得到的變異檢驗結果一致.

圖2 潦河流域年徑流的變化趨勢及突變點分析Fig.2 Trends analysis and break-point detection of annual runoff in Liaohe River Basin

潦河位于江西鄱陽湖流域，該區域自1950s以來，氣候呈現出一定程度的波動性變化特征. 潦河流域降水量呈增加趨勢，降水量在20世紀60年代初期和90年代初期發生突變，而蒸發量呈現下降趨勢，夏季尤為明顯[35]. 鄱陽湖流域20世紀60年代氣溫處于下降趨勢，70年代初期氣溫有所回升，尤其在90年代后呈現顯著增溫時期. 潦河降水量的突變增加和蒸發量的突變減少以及氣溫的增加，導致徑流在20世紀70年代初期發生突變，徑流呈現增加趨勢[36-37].

3.2 水文模擬與徑流還原

根據潦河萬家埠站水文序列突變檢驗結果，選擇突變點前的1953—1971年為基準期，用于率定和驗證DTVGM模型參數[38]. 考慮到部分雨量站數據缺測，選擇1963—1968年(6年)作為模型率定期，1969—1971年(3年)作為模型驗證期. 采用SCE-UA算法[39]，以NSE及r作為目標函數，對所構建水文模型進行率定. 表2列出了模型最優參數，表3列出了模型在率定期和驗證期的模擬精度評估，圖3給出了模型的徑流模擬結果.

表2 DTVGM模型參數意義及最優參數結果

表3 潦河流域模型率定期和驗證期結果

圖4 潦河流域1963—2014年月徑流模擬過程Fig.4 Simulated monthly hydrograph during 1963-2014 in Liaohe River Basin

圖3 潦河流域率定期和驗證期月徑流模擬過程Fig.3 Simulation of monthly runoff during calibration and validation period in Liaohe River Basin

從表3可以看出：針對日徑流過程，模型率定期和驗證期NSE分別為0.78和0.85，r分別為0.88和0.92；針對月徑流過程，模型率定期和驗證期NSE分別為0.85和0.91，r分別為0.92和0.96. 可見，率定期和驗證期模擬徑流和實測徑流擬合程度較高，從而驗證了DTVGM模型. 能夠模擬潦河流域水資源的天然分布狀況[40]，為進一步開展生態流量研究奠定了堅實的基礎. 固定DTVGM模型參數，輸入1960—2014年流域降水量，將1960—1962年設置為模型預熱期，即可獲取基于水文模型的1963—2014年潦河流域的月徑流過程(見圖4).

3.3 生態流量計算結果

該研究采用GD、ND、GEV、P-Ⅲ及LN3五種分布對潦河流域的月徑流過程進行適線分析，采用線性矩法獲取各分布參數[41]. 由于不同分布針對同一樣本系列(如7月徑流序列)的擬合存在顯著差異，因此需分析不同分布的擬合優度，對比不同分布的擬合效果，進而選擇最優分布進行生態流量估算. 以多年7月徑流序列為例，將5種分布及樣本序列繪制在一起(見圖5).

圖5 潦河流域萬家埠站7月徑流序列5種概率分布函數擬合結果Fig.5 The fitting diagrams of 5 probability distribution functions of runoff in July at Wanjiabu Station, Liaohe River Basin

由于生態流量主要取90%保證率條件下的月徑流過程，因此該研究重點聚焦于頻率為0.5～1.0之間的適線效果. 從圖5(b)可以明顯看出，5種分布中，GEV、P-Ⅲ與LN3分布在低徑流中擬合效果較好，GD分布雖常用于月降水的擬合[42]，但擬合效果明顯劣于GEV和P-Ⅲ分布. 5種分布中，ND分布擬合效果最差，這也從側面證明了水文序列為偏態分布，與已有研究結果[43]相一致. 為了定量分析不同分布對月徑流過程的擬合效果，該研究采用K-S檢驗、A-D檢驗及PPCC檢驗法對5種分布的適線性進行檢驗. 以潦河萬家埠站7月徑流序列為例，選用3種檢驗方法對5種分布進行擬合優度檢驗，結果如表4所示. 綜合3種擬合優度檢驗方法的分析結果表明，GEV分布為7月徑流序列的最優概率分布函數. 同時，繪制5種分布下理論值與觀測值之間的相關關系(見圖6)，結果表明，P-Ⅲ、GEV與LN3分布對樣本擬合效果最好，GD分布與ND分布擬合效果較差，這與表4擬合優度分析結果相一致.

表4 潦河流域萬家埠站7月徑流序列概率 分布函數擬合優度檢驗結果

選用5種分布對潦河流域逐月徑流序列進行分析，采用3種檢驗方法對頻率分布函數進行擬合優度綜合檢驗. 檢驗分析表明，針對不同月份、不同檢驗指標其最優分布函數存在差異，如對萬家埠站7月徑流序列分布函數進行檢驗，其中K-S、A-D檢驗下GEV分布為最優，而PPCC檢驗法下LN3分布為最優(見表4). 因此，需綜合分析多種方法對概率分析函數的擬合優度進行檢驗更為合理，且每個月份的水文序列相互獨立，也需要分別進行概率分析. 采用相同方法對經還原的潦河萬家埠站逐月徑流序列進行擬合優度檢驗分析，即可獲取相應水文序列的最優分布函數，結果如表5所示. 由表5可見，不同月份徑流的最優分布不盡相同，其中GEV分布在6個月中表現為最優分布，GD、P-Ⅲ和LN3均在2個月中為最優分布. 選取逐月最優分布頻率曲線，以90%保證率對應的流量作為生態基流，獲取如表5所示的生態流量結果.

將最優分布的逐月頻率計算法與Tennant法、最枯月平均流量法和7Q10法等常用方法進行對比，各方法所求生態流量及水量如表5所示. 其中Tennant法是較為傳統且經典的水文學方法，該方法將一年分為兩個階段，根據鄱陽湖流域主要魚類的繁殖規律，選定10月—翌年3月為一般用水期，4—9月為魚類產卵期. 根據Tennant法對河流流量的描述，選定“非常好”對應的生態流量作為參照標準，即選用多年平均流量的50%作為魚類產卵期的生態流量，一般用水期采用多年平均流量的30%作為河道內生態流量[44]. 最枯月平均流量法以河流最小月平均實測徑流量的多年平均值作為河流的基本生態環境需水量. 7Q10法根據各年連續7日最低徑流流量繪制理論頻率關系分布曲線，但由于該標準對國內河流要求比較高，因此采用近10年最枯月平均流量作為生態流量. 對比表5中各種方法得到的生態流量可知，該研究采用最優分布的逐月頻率計算法推求的年生態流量最大，為1 548.5×106m3；處于“非常好”水平的Tennant法計算的年生態流量次之，為1 282.7×106m3；7Q10法和最枯月平均流量法獲取的年生態流量最低，分別為706.8×106和620.6×106m3. 從生態需水的季節性特征分析而言，最枯月平均流量法和7Q10法所給出的每月生態流量均相同，無法反映出不同月份生態流量之間的差異，其中最枯月平均流量法獲取生態流量的Tennant法等級位于“差”的水平，而7Q10法獲取生態流量的Tennant法等級位于“中等”水平. 基于“非常好”水平的Tennant法推求的生態需水過程，將一年分為汛期和枯水期，雖然考慮了生態需水在汛期和枯水期的差異，但仍未考慮水文變量的季節性變化特征，因此對于水資源量年內分布極度不均的流域仍需進一步分析[45]. 該研究采用最優分布的逐月頻率計算法進行推求生態流量，根據逐月90%保證率條件下徑流作為月生態流量，所獲取的生態流量能夠考慮生態需水的動態變化，該方法所求得年生態流量與“非常好”等級下Tennant法計算所得年生態流量相近，且遠大于其他兩種水文學方法得到的年生態流量，表明該生態流量能夠滿足水生生物的正常需求，計算結果具有確定性和合理性[18,20-21].

圖6 潦河流域萬家埠站7月徑流序列5種頻率曲線理論值與觀測值散點圖Fig.6 The scatter plot of observed and theoretical values of 5 probability distribution functions for runoff in July at Wanjiabu Station, Liaohe River Basin

表5 不同方法下潦河流域生態流量計算結果比較

4 結論

a) 采用多種水文序列突變檢驗方法對潦河萬家埠站進行檢驗，發現潦河萬家埠站徑流在1972年發生變異，主要是由于降水量的突變增加和蒸發量的突變減少，導致徑流在20世紀70年代初期發生突變，徑流量呈現增加趨勢.

b) 基于突變點前基準期對DTVGM模型進行參數率定和驗證，率定期和驗證期模擬徑流和實測徑流擬合程度較高，從而驗證了DTVGM模型能夠較好地模擬潦河流域水資源的天然分布狀況.

c) 基于多種擬合優度指標對分布函數進行分析，選用最優分布函數開展生態流量計算，生態流量能夠滿足水生生物的正常需求. 與Tennant法、最枯月平均流量法、7Q10法等方法比較，基于最優分布函數的生態流量結果更具確定性與合理性，可為變化環境下水資源的適應性管理提供重要依據.