基于突變檢驗的黃河上游生態水文變異分析

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

0 引 言

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

1 研究區概況

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

2 方 法

2.1 去趨勢預置白Mann-Kendall檢驗法

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

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

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

β=median( ) ?j

(1)

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

Yt=Xt-Tt=Xt-βt

(2)

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

Y′t=Yt-r1Yt-1

(3)

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

Y″t=Y″t+Tt

(4)

該序列將不再受自相關性的影響。

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

(5)

(6)

(7)

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

(8)

(9)

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

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

(10)

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

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

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

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

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

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

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

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

(11)

或

(12)

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

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

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

(13)

式中：z為統計檢驗量；其他參數意義同前。

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

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

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

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

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

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

(4)提高指標的可操控性，為后續的生態調控奠定基礎。

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

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

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

2.4 方法耦合

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

圖1 方法耦合Fig.1 Method coupling

3 實例分析

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

3.1 變異點的識別分析

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

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

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

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

3.2 滑動秩和檢驗

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

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

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

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

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

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

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

3.3 IHA指標分析

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

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

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

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

以水文指標的改變為依據，來分析其對生態系統的影響。

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

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

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

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

4 結 語

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

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

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

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

□

[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]. 干旱區地理, 2005,8(3):88-291.

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

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

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

[7] 水利部黃河水利委員會.黃河流域及西北內陸河水功能區劃[R]. 黃河水利委員會,2001.

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

[9] 沈珍瑤，楊志峰.黃河流域水資源可再生性評價指標體系與評價方法[J].自然資源學報,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]. 水科學進展, 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] 拜存有. 渭河流域關中段徑流過程變異點診斷研究[D]. 陜西楊凌: 西北農林科技大學,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] 趙芳芳，徐宗學. 黃河蘭州以上氣候要素長期變化趨勢和突變特征分析[J].氣象學報，2006，64(2)：246-255.

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