基于突變檢驗的黃河上游生態(tài)水文變異分析

梁欣陽，盧玉東，孫東永，張洪波，楊亞慧(1.長安大學(xué)環(huán)境科學(xué)與工程學(xué)院，西安 710054;2.長安大學(xué) 旱區(qū)地下水文與生態(tài)效應(yīng)教育部重點實驗室，西安 710054)

0 引 言

近些年，在全球氣候變化與人類活動共同作用下，河川徑流的演變規(guī)律受到各國學(xué)者的關(guān)注[1]。研究表明，變化環(huán)境下世界上很多河流的徑流都表現(xiàn)出不同程度的下降趨勢[2]，這一現(xiàn)象在我國也不例外，黃河上游表現(xiàn)尤為明顯。近年來黃河上游徑流量總體呈下降趨勢，到1990年后平均水量只有270～290 億m3，平均較常年偏少15.6%[3]。在黃河上游徑流量的下降趨勢過程中，針對水文變異情況作進一步分析對于具有操作性的河流生態(tài)水文調(diào)度具有深遠的意義。如樊輝等[4]采用TFPW-MK法檢驗?zāi)陱搅髋c年輸沙長期變化趨勢，并利用貝葉斯變點分析方法診斷時序數(shù)據(jù)突變點，分析了黃河干、支流徑流量年際變化特征；李劍鋒等[5]在考慮水文變異的基礎(chǔ)上，系統(tǒng)提出河道內(nèi)生態(tài)需水的計算方法，并對黃河流域生態(tài)需水問題進行分析；于茜等[6]應(yīng)用變化范圍法(RVA)對黃河上游蘭州站的生態(tài)水文特征變異進行分析，給出了不同水庫運行模式對黃河蘭州站生態(tài)水文特征變異的影響程度情況。然而以系統(tǒng)性的思路對黃河上游蘭州斷面徑流序列進行變異診斷、分級和指標識別的研究還較少。本文基于黃河流域水文生態(tài)變異問題的重要性及其特殊意義，對研究方法及思路作了進一步的考慮：采用TFPW-MK突變檢驗法和秩和檢驗法2種方法分析黃河上游流域代表性斷面徑流變異，避免單一變異分析方法的缺陷，增加變異分析結(jié)果的可信度；進一步采用IHA水文指標改變法分析變異前后水文生態(tài)情勢。

1 研究區(qū)概況

黃河上游區(qū)從源頭至內(nèi)蒙古托克托縣河口鎮(zhèn)，河道長3 472 km，流域面積約38.6 萬km2，匯入的較大支流有43條，徑流量占全河的60%，其水文生態(tài)的變化對區(qū)域的可持續(xù)發(fā)展具有重要的影響。根據(jù)已有成果[7-9]中水文分區(qū)與生態(tài)分區(qū)的特點，考慮不同位置的調(diào)蓄水庫的影響控制范圍以及干流控制斷面的分布，將黃河上游劃分為黃河源頭-蘭州河段、蘭州-河口鎮(zhèn)段、蘭州-河口鎮(zhèn)段。本文選取劉家峽水電站下游的蘭州斷面，其流量的變化可以反映龍羊峽、劉家峽聯(lián)合調(diào)度對流量特征的影響，將對下游河段的供水灌溉、水沙關(guān)系以及生態(tài)保護等產(chǎn)生關(guān)鍵性的影響。

2 方 法

2.1 去趨勢預(yù)置白Mann-Kendall檢驗法

基于秩的Mann-Kendall檢驗法[10,11]是一種非參數(shù)統(tǒng)計檢驗方法，除了可以進行趨勢檢驗，也可用于變異點的檢驗分析。由于其不需要遵從一定的分布，不受少數(shù)異常值的干擾，且結(jié)構(gòu)簡單，計算方便，使其也成為水文序列的變異點分析中常用方法之一，但其檢驗受到趨勢項的影響。章誕武等[12]通過對比分析普通Mann-Kendall檢驗法和去趨勢預(yù)置白Mann-Kendall檢驗法得出：序列中存在的正(負) 自相關(guān)性將增大( 減少) 其方差，從而在某一顯著性水平下，增大( 降低) 拒絕原假設(shè)的概率；同時，序列趨勢項會增大其自相關(guān)系數(shù)的計算值，對于有明顯趨勢的序列，應(yīng)先剔除趨勢項后估算其自相關(guān)系數(shù)。

去趨勢預(yù)置白Mann-Kendall檢驗 (TFPW-MK, Mann-Kendall test with trend-free pre-whitening)[13,14]： 考慮數(shù)據(jù)序列明顯趨勢對自相關(guān)系數(shù)估計的影響，通過去趨勢預(yù)置白方法，剔出趨勢對自相關(guān)系數(shù)估計的影響，更準確地對數(shù)據(jù)序列進行檢驗。步驟如下。

(1) 假定序列由線性趨勢和AR(1)組成，采用Theil(1950年)和Sen(1968年)提出的方法計算新樣本數(shù)據(jù)的線性趨勢β：

β=median( ) ?j

(1)

(2)去除趨勢項，形成不含趨勢項的序列Yt：

Yt=Xt-Tt=Xt-βt

(2)

(3)計算序列Yt的一階自相關(guān)系數(shù)r1，剔除序列中的自相關(guān)項：

Y′t=Yt-r1Yt-1

(3)

(4)補還趨勢項Tt，得到不含自相關(guān)影響的新序列：

Y″t=Y″t+Tt

(4)

該序列將不再受自相關(guān)性的影響。

對經(jīng)過去趨勢預(yù)白化后形成的新序列使用Mann-Kendall檢驗法進一步檢驗序列變異，定義統(tǒng)計量為：

(5)

(6)

(7)

E(sk)=k(k+1)/4

(8)

(9)

將時間序列x按降序排列，再按上式計算，同時使：

UFk=-UB′k,k′=n+1-k,k=1,2,…,n

(10)

通過分析統(tǒng)計序列UFk和UBk，不僅可以進一步分析序列x的趨勢變化，還可以明確突變的時間，指出突變的區(qū)域。若UFk值大于0，則表明序列呈上升趨勢，小于0則呈下降趨勢。2統(tǒng)計序列構(gòu)成的曲線分別標記為UF和UB，如果2條曲線超過臨界直線，則表明上升或下降趨勢顯著；如果2條曲線出現(xiàn)交點，且交點在臨界直線之間，那么交點對應(yīng)的時刻就是突變開始的時刻。

此方法與別的檢驗方法不同的是，在檢驗時，其2條曲線的交點通常為與2個年份之間，即突變點不指向某一個確定的點。我們需要采取一定的方法從2個點中選取一點作為突變點。

2.2 滑動秩和檢驗(Mann-Whitney U檢驗)

秩和檢驗法，也稱Mann-WhitneyU檢驗法[15,16]，是非參數(shù)統(tǒng)計方法，它不依賴于總體分布類型，不以推斷總體參數(shù)為目的，而旨在檢驗2種或2種以上觀察變量的分布有無顯著差異。由于秩和檢驗法對變異點診斷有適應(yīng)性強且效率高(均值和Cv、Cs)的特點，本文采用該法來驗證TFPW-MK檢驗的結(jié)果是否合理?；瑒又群蜋z驗法計算步驟如下[15-17]。

(1)考察點前后2序列總體的分布函數(shù)Fpre(x)和Fpost(x)，從總體Fpre(x)和Fpost(x)中分別抽取容量為npre和npost的2個樣本，要求檢驗原假設(shè)Fpre(x)=Fpost(x)。

(2)將2組數(shù)據(jù)分別由小到大排列，再將2組數(shù)據(jù)由小到大統(tǒng)一編秩。如有原始數(shù)據(jù)相同時，可按相同數(shù)據(jù)全在同組內(nèi)不用求平均秩次，不同組間有相同數(shù)據(jù)需求平均秩次的原則處理。

(3)求秩和并確定檢驗統(tǒng)計量。當2樣本容量不同時，以容量小者容量為n1，容量大者容量為n2，小容量樣本的秩和為統(tǒng)計檢驗量T1，大容量樣本的秩和為統(tǒng)計檢驗量T2。若2樣本容量相同，則n1=n2，任取一組的秩和檢驗統(tǒng)計量為T1，另一組為T2。

(4)計算和作出推斷結(jié)論。統(tǒng)計量U值計算公式如下：

(11)

或

(12)

式中：U為統(tǒng)計量；n1為較小容量樣本容量；n2為較大容量樣本容量；T1為小容量樣本的秩和；T2為大容量樣本的秩和。

如果n2≤8，則分別用式(11)和式(12)計算，取較小值作為U值。然后查閱Mann-WhitneyU檢驗中U的相伴概率表得相伴概率p。

如果n2>8，無法查閱相伴概率標，可把正態(tài)近似用于U的抽樣分布來檢驗：

(13)

式中：z為統(tǒng)計檢驗量；其他參數(shù)意義同前。

2.3 水文變異的影響域識別方法

水文生態(tài)變異指標的定義：在有關(guān)生態(tài)水文聯(lián)系的研究中，人們通常用指標體系來表征和評估一個流域或區(qū)域的生態(tài)水文系統(tǒng)健康狀態(tài)。指標體系一般由一系列具有生態(tài)意義的水文要素或特征所組成，并作為生態(tài)水文系統(tǒng)對水需求的外在表現(xiàn)形式。

水文變異法(RVA, Range of Variability Approach)[18]是對擾動河流進行水文情勢變異評估的一種方法，其基本原理在于以變異點位置為依據(jù)，進行水文階段劃分，結(jié)合水文變異指標(IHA, Indicators of Hydrologic Alteration)，識別變異前后水文序列的變異領(lǐng)域(變異點前稱為參考序列，變異點后稱為變異序列)，對水文變異進行細節(jié)解析，找出發(fā)生強烈變異的水文指標。其以與生態(tài)相關(guān)的流量特征的統(tǒng)計分析為基礎(chǔ)，從量、時間、頻率、延時和變化率5個方面對32個河流水文特征(指標)進行提取(見表1)。指標的選取應(yīng)遵循以下原則。

(1)立足流域河情和生態(tài)水文特點，充分考慮河流的水流特性，選取能夠反映生態(tài)水文特征的指標。

(2)考慮已有的開發(fā)既定事實。指標的選取不僅僅以滿足河流生態(tài)需求為標準，同時還要兼顧生防洪、供電的需要。

(3)重視指標的代表性。不同斷面的水文特征不同，選取的指標要具有獨立性和代表性，同時要避免不必要的重復(fù)。

(4)提高指標的可操控性，為后續(xù)的生態(tài)調(diào)控奠定基礎(chǔ)。

(5)考慮指標的獲取難度，保證在一定時期能夠獲取全部指標用以評價或表征。

將參考序列的提取結(jié)果作為擾動閾值， 通過對比評價序列的河流水文特征落入擾動閾值的數(shù)量， 來反映河流流量受擾動的影響程度，一般以32個IHA 指標以及綜合影響程度的等級來表示，其中各個IHA指標代表的生態(tài)環(huán)境影響不同。水文指標是通過其參數(shù)的改變來反映影響生態(tài)變化的程度，水文指標的高程度改變無疑會導(dǎo)致相應(yīng)的生態(tài)體系發(fā)生相對應(yīng)的高強度擾動。表1為各個指標具體含義及其對河流生態(tài)系統(tǒng)的影響。

表1 河流的主要生態(tài)水文改變指標(IHA)及其對河流生態(tài)系統(tǒng)的影響Tab.1 Summary of IHA parameters in rivers and their ecosystem influences

2.4 方法耦合

本文采用TFPW-MK檢驗法為主要變異檢測方法，結(jié)合滑動秩和檢驗方法驗證檢測結(jié)果，得到變異年份后作為RVA法的間隔年份對IHA變異指標進行分析，見圖1。

圖1 方法耦合Fig.1 Method coupling

3 實例分析

選取蘭州斷面實測日徑流資料，資料來源于黃委水文資料匯編，序列長度為1950-01-01-2004-12-31，共55 a。本文將對黃河上游的生態(tài)水文特征變化包括徑流年際變化、汛期變化、年內(nèi)分配、極端流量變化、水文極值發(fā)生時間、高流量和低流量、水流條件變化頻率與速率、大小洪水進行分析。圖2為蘭州斷面年徑流變化趨勢。

3.1 變異點的識別分析

鑒于變異點檢驗的特殊性和復(fù)雜性，采用 TFPW-MK變異點檢驗對蘭州斷面實測年徑流序列進行變異點研究。檢驗結(jié)果如圖2所示，可以看到，所得趨勢線UFk數(shù)值基本都為負值，表示徑流序列整體呈現(xiàn)連續(xù)的下降趨勢，且在1985年徑流序列有顯著下降趨勢。根據(jù)文獻分析[19]，從20世紀80年代后期到90年代末，黃河上游氣溫逐步升高，降水量也呈現(xiàn)明顯減少趨勢， 降水量減少自然會導(dǎo)致徑流減少，且徑流變化比例大于降水量變化比例。因此，黃河上游徑流趨勢和突變特征在一定程度上是對流域內(nèi)降水變化的響應(yīng)，特別是年徑流序列20世紀80年代中期與流域降水突變時間相吻合。

圖2 蘭州站年天然徑流量變化趨勢線Fig.2 Annual natural runoff change trend line of Lanzhou station

依據(jù)TFPW-MK突變檢驗結(jié)果，UFk和UBk曲線交點位于1985年，即初步認為蘭州斷面徑流變異年份為1985年。

圖3 TFPW-MK變異檢驗Fig.3 Mann-Kendall test with trend-free pre-whitening

3.2 滑動秩和檢驗

依照時間順序，對水文序列逐點進行秩和檢驗，即逐點作為分割點，以此把水文序列分成2個樣本。設(shè)定2個假設(shè)，即：①H0，2個樣本的分布相同，即無顯著差異；②H1，2個樣本的分布不同，即具有顯著差異。

取置信水平α=0.05 ，采用雙側(cè)分位數(shù)檢驗，查得α=0.05 時，Z1-α/2=1.96 ，若|z|≤Z1-α/2，則接受H0，2個樣本無顯著差異；若|z|>Z1-α/2，則接受H1，2個樣本有顯著差異，即發(fā)生變異，并以|z|最大處作為水文變異點。表2為蘭州斷面滑動秩和檢驗成果。

表2 蘭州斷面首個變異點滑動秩和檢驗成果Tab.2 The moving rank sum test results of the first change point of Lanzhou

步驟1：對蘭州整個序列逐年進行秩和檢驗，計算出首個變異點。序列分段為1950-1985、1986-2004 時，統(tǒng)計檢驗量|z|最大，且|z|>1.96，檢驗結(jié)果為拒絕H0，2樣本分布有顯著差異。

步驟2：接著找出所有可能變異點。對1950-1985年進行滑動秩和檢驗，序列分段為1950-1963和1964-1985時，統(tǒng)計檢驗量最大，檢驗結(jié)果為拒絕H0，2樣本分布有顯著差異；對1986-2004 年進行滑動秩和檢驗，該時段內(nèi)不存在變異點?，F(xiàn)在，蘭州序列分段為1950-1963，1964-1985，1986-2004。

步驟3：每相鄰兩段序列整合成一個序列，進行變異點分析，對1964-2004年進行滑動秩和檢驗，序列分段為1964-1985、1986-2004時，統(tǒng)計檢驗量最大，檢驗結(jié)果為拒絕H0，2樣本分布有顯著差異；對1950-1985年進行滑動秩和檢驗，序列分段為1950-1963、1964-1985時，統(tǒng)計檢驗量最大，檢驗結(jié)果為拒絕H0，2樣本分布有顯著差異。每相鄰兩段序列合并成為一個序列進行檢驗，其變異點均沒有發(fā)生變化。需要再次用步驟2對1950-1963年進行滑動秩和檢驗，該時段內(nèi)不存在變異點；對1964-1985年進行滑動秩和檢驗，該時段內(nèi)不存在變異點；對1986-2004年進行滑動秩和檢驗，該時段內(nèi)不存在變異點，各序列均沒有變異點。已經(jīng)找出所有可能的變異點，最后得出唐乃亥的水文變異分段為：1950-1963、1964-1985、1986-2004。

根據(jù)圖3所示TFPW-MK檢驗結(jié)果與秩和檢驗結(jié)果對照分析，在1963年，UFk和UBk曲線沒有交點，且趨勢線UFk變化較為平緩，則認為在1963年未發(fā)生變異。而在1985年，TFPW-MK檢驗結(jié)果與秩和檢驗結(jié)果吻合，均顯示發(fā)生水文變異，則確定1985年為黃河上游蘭州斷面徑流序列的突變年份。

3.3 IHA指標分析

通過以上分析，確定了水文序列的變異時間、變異強度，但是各水文指標所發(fā)生的變異程度仍舊未知，故本文通過RVA法評價了變異點前后序列的變異范圍(指標)以及強度，結(jié)果如圖4所示。

圖4 指標改變度分析Fig.4 Alteration analyzing

通過分析可知，蘭州斷面的生態(tài)水文特征變化較大，整體改變度為0.61，屬高度改變。32個指標中發(fā)生高度改變的近1/2，并且大多集中在7、8、9、10月平均流量、年最大流量平均值、每年高流量次數(shù)和延時以及每年流量逆轉(zhuǎn)次數(shù)。高度改變的水文指標改變度如表3所示。

表3 高度改變的水文指標Tab.3 Highly changed alterations

以水文指標的改變?yōu)橐罁?jù)，來分析其對生態(tài)系統(tǒng)的影響。

(1)7、8、9、10月屬于汛期，植物生長茂盛，蒸發(fā)強度大，因此對水的需求量大，1985年后7、8、9、10月平均流量高度改變，尤其是9月份從1 583 m3/s減少到863 m3/s，改變度達到-0.1，這對下游依賴洪水流量及其創(chuàng)造的棲息地環(huán)境而生存的生物物種而言，影響是破壞性的。

(2)最小化指標整體趨勢趨于平緩，最大化指標則受其削減洪峰的影響，在1985年后出現(xiàn)了銳減，說明大洪水的洪峰一般只持續(xù)1～5 d， 最大流量平均值的高度改變會在一定程度上影響河流渠道地貌和自然棲息地的構(gòu)建。

(3)高脈沖流量通過與漫灘和高地的連通,可大量地輸送營養(yǎng)物質(zhì)并塑造漫灘多樣化形態(tài)，維系河道并育食河岸生物，從而提高河流生物的生物量和多樣性。1985年后高脈沖流量的次數(shù)和延時均有明顯下降趨勢，將會導(dǎo)致生物多樣性的下降。

(4)每年流量逆轉(zhuǎn)次數(shù)明顯增多，而水流條件的變化頻率對干旱地區(qū)河流生物物種影響顯著，陡漲陡落的流量變化會導(dǎo)致水生生物被沖刷或擱淺，洪水的暴落導(dǎo)致生物幼苗種群不能建立。由以上分析可以看出1985年后蘭州斷面流量的主要變化可以概括為汛期流量下降、洪峰削減、高脈沖流量減少和水流條件變化頻率加快，這些改變與蘭州斷面受龍羊峽、劉家峽水庫不同組合調(diào)節(jié)聯(lián)系密切。1985年以后，龍羊峽水庫建成，龍劉聯(lián)合調(diào)度使得防洪能力和調(diào)蓄能力進一步增強,蓄豐補枯，汛期水流被削平，枯水流量被填補，峰谷缺失，以致水流情勢改變嚴重。7-10月份流量過程趨于平緩，最大化流量平均值發(fā)生較大改變。同時水庫兼有供水灌溉、發(fā)電、防洪任務(wù),相機供水的次數(shù)增多，水流流量起落變化頻繁。

4 結(jié) 語

通過對黃河上游蘭州斷面的實測徑流水文序列變異點的識別和分析可得出以下結(jié)論。

(1)本文所采用的方法經(jīng)過耦合從多方面對水文序列進行了檢驗，較全面地反映了時間序列的變異特性，因而檢驗分辨率較高。

(2)通過變異診斷分析，發(fā)現(xiàn)蘭州斷面水文生態(tài)指標在龍、劉兩庫運行后發(fā)生較高程度改變，這無疑非常不利于區(qū)域水文生態(tài)安全，建議在黃河生態(tài)水文指標體系的構(gòu)建和水庫生態(tài)調(diào)度中給予充分考慮。

(3)本研究以黃河流域蘭州斷面為例，研究了黃河上游生態(tài)水文變異情況，取得了一定的成果，但在實踐應(yīng)用層面仍存在問題，比如如何使用生態(tài)水文變異分析結(jié)果指導(dǎo)區(qū)域的水資源管理實踐，在后續(xù)研究中需重點考量。作者認為生態(tài)水文變異分析應(yīng)該成為河流水資源管理和水利工程規(guī)劃中的重要限制性因素。因此，建立生態(tài)水文聯(lián)系變異指標與水事活動間的關(guān)系，并依據(jù)其關(guān)系設(shè)置限建規(guī)則，將成為生態(tài)水文變異分析應(yīng)用于流域水資源管理實踐的重要手段。

□

[1] IPCC.The physical science basic[C]∥Contribution of working group I to the fourth assessment report of the intergovernmental panel on climate change(IPCC). Paris, 2007.

[2] Kundzewicz Z W,Graczyk D. Trend detection in river flow series:1. Annual maximum flow [J]. Hydrological Sciences Journal, 2005,50(5):797-810.

[3] 王金花， 康玲玲，余 輝，等. 氣候變化對黃河上游天然徑流量影響分析[J]. 干旱區(qū)地理, 2005,8(3):88-291.

[4] 樊 輝，楊曉陽. 黃河干、支流徑流量與輸沙量年際變化特征[J]. 泥沙研究, 2011,(4):11-15.

[5] 李劍鋒，張 強，陳曉宏，等. 考慮水文變異的黃河干流河道內(nèi)生態(tài)需水研究[J]. 地理學(xué)報, 2011,66(1):99-110.

[6] 于 茜，黃 強，張洪波. 不同水庫運行模式對黃河蘭州站生態(tài)水文特征變異的影響[J]. 干旱區(qū)地理, 2010,33(5):747-755.

[7] 水利部黃河水利委員會.黃河流域及西北內(nèi)陸河水功能區(qū)劃[R]. 黃河水利委員會,2001.

[8] 黃河流域氣候編寫組.黃河流域氣候[R].黃河水利委員會水文局,1987.

[9] 沈珍瑤，楊志峰.黃河流域水資源可再生性評價指標體系與評價方法[J].自然資源學(xué)報,2002,17(2):188-197.

[10] Kendall M G.Rank correlation measures[M]. London:Charles Griffin,1976:110-120.

[11] Mann H B.Non-parametric tests against trend[J].Econometrica,1945,13:245-259.

[12] 章誕武，叢振濤，倪廣恒. 基于中國氣象資料的趨勢檢驗方法對比分析[J]. 水科學(xué)進展, 2013,24(4):490-496.

[13] Sheng Yue,Wang Chunyuan.The mann-kendall test modified by effective sample size to detect trend in serially correlated hydrological series[J].Water Resources Management,2004,18:208-218.

[14] Sheng Yue,Paul P,Bob P,et al.The influence of autocorrelation on the ability to detecttrend in hydrological series[J].Hydrological Processes,2002,16(1):1 807-1 829.

[15] Seigel S. Nonparametric statistics for the behavioral sciences[Z]. New York: McGraw-Hill Book Company, 1956.

[16] Mann H B, Whitney D R. On a test of whether one of two random variables is stochastically larger than the other[J]. The Annals of Mathematical Statistics, 1947,18(1):50-60.

[17] 拜存有. 渭河流域關(guān)中段徑流過程變異點診斷研究[D]. 陜西楊凌: 西北農(nóng)林科技大學(xué),2008.

[18] Brian D Richter,Jeffrey V Baumgartner,Jennifer Powell,et al.A method for assessing hydrologic alteration within ecosystems[J].Conservation Biology,1996,10(4):1 163-1 174.

[19] 趙芳芳，徐宗學(xué). 黃河蘭州以上氣候要素長期變化趨勢和突變特征分析[J].氣象學(xué)報，2006，64(2)：246-255.

[20] 馮 波，閆佰忠，章光新. 松花江流域水文干旱聯(lián)合概率分布特征研究[J]. 節(jié)水灌溉，2014,(5):38-42.