基于ERHIs法的雅魯藏布江流域水文情勢變化研究

楊 靖，藺子琪，趙慶緒，謝 靜，廖衛(wèi)紅

(1.華能西藏水電安全工程技術研究中心，四川 成都 610093；2.四川大學 水利水電學院，四川 成都 610065；3.華中科技大學 土木與水利工程學院，湖北 武漢 430072；4.長江勘測規(guī)劃設計研究有限責任公司，湖北 武漢 430010；5.水利部水網(wǎng)工程與調(diào)度重點實驗室，湖北 武漢 430010；6.中國水利水電科學研究院 流域水循環(huán)模擬與調(diào)控國家重點實驗室，北京 100038)

1 研究背景

雅魯藏布江是世界上海拔最高、西藏最大的一條河流，也是我國西南地區(qū)非常重要的跨界河流之一。雅魯藏布江流域徑流主要由降水、冰雪融水和地下水補給組成[1]。由于流域面積廣闊，流域內(nèi)的降水量、蒸發(fā)量以及地形、冰雪覆蓋面積、植被等條件差異很大，流域內(nèi)的徑流補給類型復雜多樣[2-3]。精確描述徑流特征及其變化的基礎是獲得大量高精度流量數(shù)據(jù)，特別是在水文站點缺乏的地區(qū)。分布式水文模型為流域徑流數(shù)據(jù)的延展和加密提供了重要的支撐[4-5]，基于此，可量化評估人類活動、下墊面和氣候變化等變化條件對水文特性影響程度，從而可為區(qū)域水資源合理開發(fā)利用提供重要的依據(jù)[6]。

針對氣候變化和水資源開發(fā)利用而導致的水文改變度的分析一直為生態(tài)水文研究領域的熱點問題[7-8]，Richter等[9]在1996年提出了一套能夠評估水生態(tài)水文改變過程的指標體系(Indicators of Hydrologic Alteration，IHA)，目前被廣泛應用于生態(tài)水文領域[10-14]。但Olden等[15]學者分析結(jié)果表明，眾多水文變量之間存在十分突出的冗余性問題，而IHA指標與大多數(shù)水文變量之間存在較高的信息重疊現(xiàn)象。相對于早期研究中只考慮單一水文變量，盡管IHA體系較完善，但是水文變量仍然較多，IHA變量之間也存在著自相關性和嚴重的冗余性問題。IHA法目前被認為可以相對最全面、最系統(tǒng)的來評估水文情勢變化特征。此外，雅魯藏布江地區(qū)的水文情勢研究較多，但具體對于環(huán)境流量組分的變化分析較少[16]，本文采用改進后的IHA法[17]降低生態(tài)水文變量中的冗余性，保留原始數(shù)據(jù)中的關鍵信息，分析雅魯藏布江下游、尼洋河出口斷面受上游水庫影響的水文改變度，聚焦于對流量組分分級以及缺資料地區(qū)的水文情勢評估。以期為分析雅魯藏布江地區(qū)水資源管理和生態(tài)治理提供科學參考。

2 研究概況

2.1 研究區(qū)概況雅魯藏布江(28°N—31°N、82°E—97°E)發(fā)源于西藏自治區(qū)西南部喜馬拉雅山中段北麓的杰馬央宗冰川，是世界上海拔最高、西藏最大的一條河流。雅魯藏布江從源頭～入海口全長約3162 km，其中，中國境內(nèi)長約2057 km。入?？谔幙偭饔蛎娣e175萬km2，在中國境內(nèi)流域面積約33.20萬km2，雅魯藏布江徑流總體以降水補給為主，流域內(nèi)多年平均降水量300～500 mm，且年際變化較小，年內(nèi)分配不均。流域多年平均凈流量約1655億m3，年徑流深與年降水量的分布趨勢具有較好的一致性[18]。測站分布位置如圖1所示。

圖1 研究區(qū)地理位置和觀測站點示意圖Fig.1 Schematic diagram of geographical location and observation stations in the study area

2.2 研究數(shù)據(jù)奴下水文站位于雅魯藏布江流域干流，控制流域面積占其國內(nèi)總流域面積的57.6%。尼洋河出口斷面的更張測站為雅魯藏布江的上游控制站，是青藏高原重要的水文站點之一[19]。本文根據(jù)雅魯藏布江流域水文站點分布及各站基本資料條件，選用奴下、更張站作為徑流計算的依據(jù)站，進行分析評價采用的是1965—2018年的徑流日值數(shù)據(jù)。然而，由于資料獲取難度問題，本文分析的幾個斷面的數(shù)據(jù)均存在不同程度的資料缺失，如奴下站點缺1958—1959年、1969年的實測資料。因此，本文采用分布式水文模型EasyDHM[20-22]進行主要站點逐日流量數(shù)據(jù)的插補還原。

2.3 研究方法論文基于Richter等[9，23-24]提出的IHA指標進行分析，通過重新構(gòu)建IHA指標體系將降水時期劃分為枯水期(1—5月和11—12月)和豐水期(6月—10月)兩個時段，將IHA指標劃分為不同時段內(nèi)IHA指標，進一步優(yōu)化IHA體系計算生態(tài)最相關水文指標 ERHIs，從而達到減少水文指標的同時盡可能多地保留原數(shù)據(jù)關鍵信息[17]。然后基于ERHIs進一步優(yōu)選參數(shù)，采用6種突變檢驗診斷指標變化規(guī)律，最后應用水文改變度描述雅魯藏布江徑流變化情況。

2.3.1 水文改變指標(IHA) Richter提出的IHA體系，該方法定量分析了水文變異程度，目前該指標體系被普遍用來代表水文情勢變化特征，該指標體系包含33個指標，目前被認為可以相對最全面、最系統(tǒng)的來評估水文情勢變化特征。該體系共分為5組，詳細指標劃分可參照文獻[23]。由于研究時段內(nèi)雅魯藏布江未出現(xiàn)斷流現(xiàn)象，故本研究不考慮“斷流天數(shù)”這一變量[25]。

2.3.2 突變診斷 采用貝葉斯推斷[26]、曼-肯德爾(Mann-Kendall)法[27-28]、勒帕熱(Le-Page)法[28]、BG(Bernaola-Galvan)分割算法[29]以及有序聚類法[30]對改進后的指標突變點進行識別，用以排除單一突變診斷方法的局限性以及數(shù)據(jù)漂移問題。

(1)貝葉斯推斷

貝葉斯推斷[26]是利用貝葉斯定理來修正特定假設的可能性，該方法可以充分利用先驗信息和樣本信息，得到統(tǒng)計結(jié)論。計算公式為：

p(X，α)=p(X|α)p(α)=p(α|X)p(X)

(1)

(2)

p(α|Y)∝p(α)p(Y|α)

(3)

(2)M-K法

M-K法[27-28]作為常見的非參數(shù)統(tǒng)計檢驗方法，無需樣本遵從一定的分布，適用于類型變量和順序變量，計算簡便，可以明確突變開始的時間并指出突變區(qū)域。計算公式如下：

(4)

(5)

(6)

其中，ai是第i時刻數(shù)值個數(shù)的累計數(shù)，E(ai)，D(ai)分別是ai的均值和方差，統(tǒng)計量UFi為標準正態(tài)分布。按時間序列x逆序xn，xn-1，…，x1，再重復上述過程，同時使UBk=-UFk，其中k=n，n-1，…，1。

2）經(jīng)營健康穩(wěn)健，盈利質(zhì)量和運營效率提升。公司過去三年每年經(jīng)營活動現(xiàn)金流量凈額均高于凈利潤，業(yè)績質(zhì)量較高。公司運營效率這兩年也有很大改善，存貨周轉(zhuǎn)率不斷優(yōu)化，2016年和2017年存貨周轉(zhuǎn)率分別為1.69次和2.07次，今年上半年已經(jīng)達到1.07次，較17年繼續(xù)提升7%，運營效率的提升將會顯著減少存貨減值準備計提，對業(yè)績提升有很大幫助。

(3)勒帕熱法

勒帕熱(Le Page)法[31]是一種具有兩個樣本不受分布影響的非參數(shù)判別法，計算公式如下：

(7)

(8)

U為統(tǒng)計量，式(8)是勒帕熱統(tǒng)計量。注意到當樣本量足夠大時，WU漸近具有自由度為2的χ2分布。其余參數(shù)含義同上。

(4)BG分割算法

BG分割算法[29]主要針對非平穩(wěn)時間序列進行突變檢測，其核心思想是通過顯著性檢驗找到突變點的左右兩部分序列的差異最大的位置。計算公式如下：

(9)

(10)

(5)有序聚類法

有序聚類分析法[32]進行突變點的識別，該方法核心思想是在分類時不打亂次序的情況下，通過計算突變點前后點位的離差平方和尋求最優(yōu)分割點。計算公式為：

2.3.3 水文改變度 基于IHA指標體系，Richter等[24]在1997年提出了用于判斷IHA指標改變程度的變化范圍法(Range of Variability Approach，RVA)，并定義了RVA閾值概念，即在IHA的基礎上，將各指標發(fā)生幾率的75%和25%作為RVA的上下閾值，用來作為評價各指標受影響程度的標準[33-35]。為了量化ERHIs受干擾后的變化程度，采用水文改變度來評估[36]，其定義如下：

式中：Di為第i個ERHIs指標的水文變化程度；Nio為第i個ERHIs指標變化在之后(之前)仍在RVA閾值范圍內(nèi)的實際觀測年數(shù)；Nie為第i個ERHIs突變(前)后預計落在RVA閾值內(nèi)的年數(shù)，可以表示為r×NT，其中r為突變前(后)第i個ERHIs落在RVA閾值內(nèi)的比例(本文取50%)，NT是受突變(前)和后影響的時間序列中記錄的總年數(shù)。D0為流域整體的水文改變度，規(guī)定水文變化程度在0～<33%之間為無變化或低程度變化，33%～<67%為中度變化，67%～≤100%之間是高度變化。

3 模型驗證與指標確定

3.1 測站驗證依據(jù)流域和數(shù)據(jù)資料情況，論文對奴下、更張斷面的1999到2012年水文日模擬精度進行了參數(shù)率定，將優(yōu)化得到的參數(shù)對流域各主要斷面1999—2018年的逐日水文過程進行了模擬，各站Nash系數(shù)評價見表1，說明論文采用的EasyDHM模型用于日徑流資料的延展和加密是可行的。依據(jù)流域和數(shù)據(jù)資料情況，本研究對奴下、更張斷面挑選洪水場次進行了參數(shù)率定和驗證，并對精度進行了評價，針對流域特點，主要針對洪峰誤差精度對方案合格程度進行了評定，其中水庫入庫數(shù)據(jù)主要通過出庫和蓄量反推，本身就存在誤差，因此對其洪峰許可誤差做了適當放大。選取1998—2018年間的20組徑流過程進行率定，徑流過程擬合效果較好，洪峰流量誤差和峰現(xiàn)時間誤差基本符合要求。對1998—2018年間的20組徑流過程樣本進行評定，洪水場次合格率為90%。EasyDHM模型采用SCE-UA全局參數(shù)優(yōu)化算法進行參數(shù)率定[37-40]。

表1 EasyDHM模型模擬徑流Nash系數(shù)表Table 1 Nash coefficient of simulated runoff by EasyDHM model

3.2 IHA的冗余性分析IHA體系包含33個指標，由于研究時段內(nèi)雅魯藏布江未出現(xiàn)斷流現(xiàn)象，故本研究不考慮“斷流天數(shù)”這一變量[25]。論文對所選擇的32個IHA指標進行冗余性分析，主要針對指標間的相關系數(shù)，認為相關系數(shù)小于0.4為冗余度較小。從圖2可以看出，部分水文變量之前展現(xiàn)了較強的相關性，10月流量，11月流量，12月流量，年最小 1 日、3 日、7 日、30 日流量之間的相關性較強，高流量年內(nèi)發(fā)生次數(shù)，高流量年內(nèi)平均歷時，上升率，下降率之間相關性較強。從圖3可以看出，年最小 3 日、7 日、30 日流量，年最大 1 日、3 日、7 日、30 日、90日流量之間也具有較強的相關性，其相關系數(shù)甚至超過了0.9。結(jié)果表明，雅魯藏布江的站點斷面的IHA指標變量相互之間顯然存有高度自關性，即指標冗余性問題非常明顯。

圖2 奴下斷面IHA指標相關關系熱力圖Fig.2 Correlation heat map of IHA index under Nuxia section

圖3 更張斷面IHA指標相關關系熱力圖Fig.3 Correlation heat map of IHA index under Gengzhang section

3.3 IHA指標體系重建鑒于IHA變量之間存在的冗余性問題比較突出，本文通過重新構(gòu)建IHA指標體系將降水時期劃分為枯水期(1—5月和11—12月)和豐水期(6月—10月)兩個時段，將IHA指標劃分為不同時段內(nèi)IHA指標，進一步優(yōu)化IHA體系計算ERHIs，在緩解IHA變量冗余性時，減少水文指標，在保證全面描述流域水文情勢變化時盡可能多地保留原數(shù)據(jù)關鍵信息，通過涵蓋豐枯時段以及指標所在組別信息來確定。

流域的徑流量以及徑流特征受到季節(jié)變化的影響，由于不同季節(jié)降水量的不同，可將降水時期分為豐水期和枯水期。由于一年內(nèi)枯水期歷史長短與流域自然地理及氣象條件相關，根據(jù)豐水期和枯水期的特點將雅魯藏布江奴下斷面6月至10月劃分為豐水期，將1—5月、11—12月劃分為枯水期(表2)。

表2 重建后IHA指標Table 2 IHA indicator after reconstruction

為使篩選出的ERHIs指標具有較高的可靠性，本文通過選取樣本年并挑選各樣本年所對應的相似年，計算相似年各IHA指標的貢獻率，選取貢獻率排在前六位且重復出現(xiàn)次數(shù)最多的IHA指標為ERHIs指標。

相似年通常指在表征同一流域有兩年或多個年份具有相同的水文指標或徑流特征或其他流域特征時，都可以稱這兩年或多個年份之間在某一特征前提下存在相似性，即這兩年或多個年互為相似年[41-42]。相似年在各個領域的應用十分廣泛，在天文、氣象、水利等各行業(yè)的預報都有應用。本文從奴下、更張斷面57年數(shù)據(jù)中分別隨機選取十年作為相應斷面的樣本年，將能夠反應徑流特征的水文統(tǒng)計值作為選取樣本年所對應相似年的基準，運用二維歐氏距離篩選方法結(jié)合統(tǒng)計值對選取好的樣本年進行相似年份的選取[41-42]。

為對應不同時期IHA指標體系，將1962—2018年流量數(shù)據(jù)按枯水期(1—5月和11—12月)和豐水期(6月—10月)分開，分別計算57年數(shù)據(jù)的枯水期統(tǒng)計值和豐水期統(tǒng)計值。統(tǒng)計值包括：均值、方差、極差、中位數(shù)、25分位數(shù)、75分位數(shù)、變異系數(shù)、偏度和峰度。

3.4 ERHIs指標確定通過計算各樣本年所對應十個相似年枯水期和豐水期IHA指標的變異系數(shù)，將十個相似年中各年枯水期和豐水期變異系數(shù)前六個指標挑選出來并記錄各指標在十年中出現(xiàn)的重疊次數(shù)，挑選出枯水期和豐水期出現(xiàn)次數(shù)排在前六的指標，將這十二個指標再次按照出現(xiàn)次數(shù)排序，查看前六個指標在一級指標中的覆蓋度是否合理。通過計算各樣本年所對應十個相似年枯水期和豐水期IHA指標的變異系數(shù)，將十個相似年中各年枯水期和豐水期變異系數(shù)前六個指標挑選出來并記錄各指標在十年中出現(xiàn)的重疊次數(shù)，挑選出枯水期和豐水期出現(xiàn)次數(shù)排在前六的指標，將這十二個指標再次按照出現(xiàn)次數(shù)排序，查看前六個指標在一級指標中的覆蓋度是否合理。經(jīng)合理性確認，將前六個指標確定為ERHIs。在合理性判定的基礎上計算各指標的相關性，經(jīng)計算后發(fā)現(xiàn)只有兩個指標存在相關性，其他指標之間的相關系數(shù)都小于0.4為不相關，將這些指標確定為篩選指標。具體步驟見圖4。

圖4 ERHIs指標確定步驟Fig.4 ERHIs indicator determination steps

指標的變異系數(shù)反應該指標在所有IHA指標中的鑒別能力。變異系數(shù)越大，表明該指標的信息分辨能力就越強，從而保證了篩選出的指標對流域影響最大[43]。初步篩選出的ERHIs指標見表3、表4。

表3 更張斷面ERHIs指標Table 3 ERHIs index of Gengzhang section

表4 奴下斷面ERHIs指標Table 4 ERHIs index of Nuxia section

3.4.1 ERHIs指標的合理性分析 為驗證2個斷面選出的ERHIs指標相互獨立且具有較高的代表性，現(xiàn)對選出的ERHIs指標進行合理性分析，本文中主要通過ERHIs指標的覆蓋性分析和冗余性分析來進行合理性分析。在減少IHA變量冗余性的同時，以保證ERHIs指標可以均勻覆蓋不同時期以及不同組別，有效保留原始信息。

(1)ERHIs覆蓋性分析

覆蓋性分析通過對各ERHIs指標在組別中的覆蓋范圍以及各指標在不同時期的覆蓋范圍進行指標的合理性分析，研究斷面指標覆蓋范圍分析結(jié)果見表5、表6。

表5 更張斷面覆蓋性分析Table 5 Coverage analysis of Gengzhang

表6 奴下斷面覆蓋性分析Table 6 Coverage analysis of Nuxia

由表可以得出，更張斷面六個ERHIs指標涉及到了第一組、第三組至第五組且達到豐枯水期全覆蓋，奴下斷面的六個ERHIs指標涉及到了第一組至第五組指標且達到豐枯水期全覆蓋。可以看出雅魯藏布江下游的2個斷面挑選出的ERHIs指標在組別中的覆蓋率以及各指標在不同時期的覆蓋率都很高，得出挑選出的ERHIs指標具有一定的合理性。

(2)ERHIs冗余性分析

為驗證5個斷面選出的ERHIs指標相互獨立且具有較高的代表性，現(xiàn)對選出的ERHIs指標進行相關性分析，確定ERHIs指標各指標間的相關性，圖5展示了各指標間的相關系數(shù)。結(jié)果表明：大部分指標之間存在較小的相關性，即信息并不存在冗余性問題，則可最終確定為ERHIs指標。

圖5 更張斷面和奴下斷面ERHIs相關關系熱力圖Fig.5 Heatmap of the correlation relationship of ERHIs in Gengzhang and Nuxia section

4 計算結(jié)果分析

4.1 初步診斷初步診斷是對已經(jīng)挑選出的ERHIs指標進行指標變異強度的確定，將強變異指標確定為最終用于識別變異年份的指標。變異系統(tǒng)的初步診斷采用Hurst系數(shù)法[44]分別對雅魯藏布江的更張水文站、奴下水文站的ERHIs指標進行變異檢驗。經(jīng)檢驗得出尼洋河站的10月指標Hurst系數(shù)為0.90屬于強變異，奴下站的4月指標Hurst系數(shù)為0.86屬于強變異。計算結(jié)果見表7、表8。

表7 更張站ERHIs指標變異強度Table 7 Intensity of variation of ERHIs indicators in Gengzhang section

表8 奴下站ERHIs指標變異強度Table 8 Intensity of variation of ERHIs indicators in Nuxia section

4.2 詳細診斷根據(jù)初步診斷結(jié)果針對各斷面挑選出的強變異指標進行詳細診斷。本文利用水文變異診斷系統(tǒng)的跳躍診斷法和氣候突變檢測方法包括Mann-Kendall方法、有序聚類法、貝葉斯推斷方法、Le-Page方法和BG(Bernaola-Galvan)分割法等方法進行突變年份的確定。更張斷面有4種方法診斷10月指標在2000年發(fā)生顯著突變，有2種方法診斷其在2014年發(fā)生顯著突變。奴下斷面有3種方法診斷4月指標在2013年發(fā)生顯著突變。具體診斷結(jié)果見表9。

表9 雅魯藏布江各斷面強變異指標的突變診斷結(jié)果Table 9 Mutation diagnosis results of strong variation indicators in each section of the Yarlung Zangbo River

結(jié)果表明，不同突變方法對突變點的識別有所差異，由于LePage方法需要人為確定子序列長度，可能會由于突變點的漂移造成誤差，只用一種突變方法難以全面的評判指標突變的結(jié)果。綜上所述，運用多種突變檢驗方法分析更張、奴下站57年徑流數(shù)據(jù)得出以下結(jié)論：①對于更張斷面：10月指標突變時段為1990—2014年。②對于奴下斷面：4月指標突變時段綜合序列為1970—2016年。造成奴下斷面指標綜合突變時間較長的原因可能為上游支流拉薩河直孔及旁多水庫對下游的調(diào)蓄作用，2010年沖久水庫投入運行后，改變了下游河段的天然徑流過程。因而，奴下斷面通過3種突變方法中體現(xiàn)出這種突變狀態(tài)，其中Hurst系數(shù)指標顯示年最小流量出現(xiàn)時間、年最小1日流量指標均有中等程度的改變。

4.3 強變異指標變化程度本文通過判斷變異前后時間序列的長短來確定穩(wěn)定時段，通過以上對奴下和更張斷面水文情勢突變點的研究，認為時間序列較長一段的水文情勢更為穩(wěn)定[24]。將更張斷面1965—1989年作為穩(wěn)定時段，將奴下斷面1965—1969年作為穩(wěn)定時段。以Richter等[26]定義的閾值概

念即將穩(wěn)定時段強變異指標發(fā)生幾率的75%作為RVA的上閾值，發(fā)生幾率的25%作為RVA的下閾值作為標準。

論文分別計算各斷面強變異指標包括更張斷面10月指標、奴下斷面4月指標的水文改變度(表10)。結(jié)果顯示，兩個斷面的強變異指標都發(fā)生了高度改變。綜上得出，更張斷面和奴下斷面整體水文改變度為高度改變。

表10 雅魯藏布江2個斷面強變異指標改變度計算結(jié)果Table 10 Calculation results of the change degree of strong variation index of two sections of the Brahmaputra River

5 結(jié)論

(1)通過EasyDHM插補雅江流域站點數(shù)據(jù)是可行的，其中，參數(shù)優(yōu)化后的Nash系數(shù)最大為0.97，最小為0.67，均在合理范圍內(nèi)，在缺資料地區(qū)運用較好。此外，ERHIs方法有效的通過減少變量簡化水文響應推求過程，以獲得水文情勢變化情況，提高水資源調(diào)配效率，對于水資源管理和生態(tài)系統(tǒng)保護都具有重要的參考價值和借鑒意義。

(2)ERHIs不僅有效緩解了IHA的冗余性問題，還有利于抓住最關鍵的生態(tài)水文變量。從IHA指標相關性熱圖可以看出，IHA參數(shù)之間的相關性較高，說明指標間的冗余問題十分突出。ERHIs指標突變結(jié)果分析表明，運用歐式距離結(jié)合相似年篩選IHA指標與多種突變結(jié)果的橫向?qū)Ρ纫约癊RHIs自身的縱向?qū)Ρ?，都表明所選取的生態(tài)最相關指標是合理的。

(3)氣候變化和人類活動對環(huán)境組分流量和水文狀況的構(gòu)成都有影響，而大規(guī)模水利建設是其中的主導因素。針對更張斷面和奴下斷面的水文改變度分析可以得出，在受到多布水電站建設的影響，尼洋河下游水文改變度差異顯著，這和張冬艷、劉海平等人的研究結(jié)果相一致[45-46]。

(4)上游水電站建設運行對水文情勢有不同程度的影響，水文情勢變化總體屬高度改變，非汛期所受影響程度比汛期更大。干旱季節(jié)的徑流變化較大，除受到上游水電站影響外，氣溫的增溫效應導致冰川融雪的補給量增加也在影響雅魯藏布江的年內(nèi)徑流分配，但這種變化還需相關的資料和分析去驗證。