黃河近百年徑流演變特征分析

胡慧杰，崔 凱，曹 茜，李姝蕾，常肖杰，沈麗娜

（1.河南黃河水文勘測設(shè)計(jì)院，河南鄭州450004；2.黃河勘測規(guī)劃設(shè)計(jì)研究院有限公司，河南鄭州450003；3.河南省許昌水文水資源勘測局，河南許昌461000）

受氣候變化的影響和人類活動的干預(yù),流域水循環(huán)規(guī)律逐漸發(fā)生變化,徑流演變規(guī)律也隨之改變。河川徑流量是流域水資源科學(xué)開發(fā)利用、合理調(diào)度配置的基礎(chǔ)［1］,因此研究徑流演變特性對流域水資源的高效利用、合理規(guī)劃具有重要意義。

黃河是我國的第二大河,承擔(dān)著流域供水、灌溉、防洪、防凌、生態(tài)等多方面任務(wù)。黃河水情的不利形勢將對沿黃地區(qū)的工農(nóng)業(yè)生產(chǎn)、生態(tài)環(huán)境造成極大影響,研究黃河徑流的演變過程,分析其演化特性,有助于深刻認(rèn)識水資源變化規(guī)律,科學(xué)控制水資源開發(fā)利用過程。近幾年很多學(xué)者對上游唐乃亥站徑流有較深入研究［2］,而對黃河下游花園口水文站徑流研究不多。花園口水文站是國家級水文站,地處黃河中下游的節(jié)點(diǎn)部位［3］,發(fā)揮著為黃河下游防洪調(diào)度與水資源管理提供重要水文信息的作用?；▓@口水文站以上的匯流面積占黃河流域總面積的97%,來水量占黃河總水量的96.6%,因此花園口水文站來水量大小基本反映黃河水量的豐、枯狀況,分析花園口水文站近百年（1919—2016年)的徑流演變可以反映黃河流域的徑流演變特性,從而為黃河水資源配置調(diào)度及管理規(guī)劃提供重要的參考依據(jù)。

1 數(shù)據(jù)來源

花園口水文站1938年建站,1949年1月起有連續(xù)的觀測資料,建站之前的徑流數(shù)據(jù)根據(jù)陜縣站1919年以來的徑流量插補(bǔ)延長［4］,最終得到花園口水文站1919—2016年各水文年（7月至次年6月)的天然徑流量。本文采用經(jīng)過還原的天然徑流時間序列分析黃河徑流的演變規(guī)律。

2 數(shù)據(jù)分析

徑流時間序列是水文時間序列的一種,一般包含確定成分和隨機(jī)成分。為深刻反映黃河流域徑流變化規(guī)律,本文運(yùn)用水文時間序列分析方法,對徑流數(shù)據(jù)的突變點(diǎn)及突變前后的趨勢性進(jìn)行分析。

2.1 統(tǒng)計(jì)特征分析

水文學(xué)常用的統(tǒng)計(jì)特征值包括平均值、極值、極差、均方差、變差系數(shù)、偏態(tài)系數(shù)、年際極值比等［5］。其中,變差系數(shù)和極值比反映的是徑流的相對變化程度,這兩個特征值越大表明徑流的年際變化越劇烈,徑流資源的利用難度越大,對水資源管理的要求越高［6］,反之,則表示徑流的年際變化較為平緩,有利于水資源的開發(fā)利用。

花園口1919—2016年各水文年度的徑流量變化見圖1,其中：1964—1965年度的天然徑流量最大,為938.7億m3；2002—2003年度的天然徑流量最小,為257.4億m3。徑流量最大值為最小值的3倍多,極差為681.3億m3,說明黃河來水年際間相差較大。

圖1 花園口水文站年天然徑流量變化過程

1919—2016年多年平均徑流量為532.5億m3,標(biāo)準(zhǔn)差為135.4億m3,變差系數(shù)為0.25。第一個黃河流域規(guī)劃于1955年通過,此后開展了黃土高原水土流失治理,建設(shè)了一大批淤地壩；流域第一座大型年調(diào)節(jié)水庫龍羊峽水庫1986年建成并投入使用；1999年開始實(shí)施全河統(tǒng)一水量調(diào)度。按以上重要年份將徑流時間序列分為不同的時段,各時段的特征值見表1。其中,1956—2016年的標(biāo)準(zhǔn)差較大,其序列值較為分散,即1956年以來徑流量年際變化較為劇烈,徑流資源的利用難度較大。

1956—1986年徑流量平均值為606.1億m3,變差系數(shù)為0.23,相對1956—2016年的變化較為集中。1987—2016年徑流量均值為450.6億m3,變差系數(shù)為0.20,偏態(tài)系數(shù)、極差相對1956—1986年均較小?？梢?1987年以來黃河徑流相對集中,極差為373.0億m3,而1999年以來黃河來水較少并出現(xiàn)負(fù)偏現(xiàn)象。

表1 各時段徑流時間序列特征值

從特征值分布來看,1956—1986年徑流量平均值最大,1999—2016年徑流量平均值最小,表明1956—1986年黃河來水相對較多。1999—2016年徑流量標(biāo)準(zhǔn)差及變差系數(shù)最小,其次是1987—2016年,即近30 a的黃河來水變化相對較為平緩。1956—2016年徑流量偏態(tài)系數(shù)最大,正偏幅度最大；各時段中只有1999—2016年偏態(tài)系數(shù)為負(fù),出現(xiàn)負(fù)偏。1919—2016年與1956—2016年兩時段的極差最大,極差最小的時段為 1999—2016年,其次為1987—2016年,可見近30 a黃河流域的來水變化幅度相對較小。

2.2 突變分析

突變分析即對數(shù)據(jù)序列的變異進(jìn)行診斷,突變主要包括均值突變、方差突變、趨勢突變、概率突變、譜突變、模型參數(shù)突變等幾種類型［7］。水文時間序列的突變主要指均值突變或趨勢突變,即時間序列從一個均值狀態(tài)或者趨勢突然跳躍到另一個均值狀態(tài)或趨勢的突變［8］。突變分析常用的方法主要包括有序聚類分析法、Mann-Kendall突變檢驗(yàn)法、滑動t檢驗(yàn)法等。

2.2.1 有序聚類分析法

有序聚類分析法是通過統(tǒng)計(jì)分析提取序列突變點(diǎn)的有效方法,主要思想是利用離差平方和判斷突變點(diǎn),認(rèn)為同類之間的離差平方和最小,而類與類之間的離差平方和最大［9］。采用該方法可以分析水文時間序列中徑流量或降水量的突變點(diǎn)。

設(shè)有水文時間序列xt（x1,x2,…,xk-1),突變點(diǎn)為τ,則突變前后的離差平方和分別為

式中：ˉxτ、ˉxn-τ分別為τ前后兩部分的均值。

總離差平方和為

取S（τ) 的極小值,則當(dāng)S＝minSn（τ) （2≤τ≤n-1)時,分割點(diǎn)兩側(cè)的離差平方和最小,即τ為最優(yōu)二分割（突變點(diǎn))。

運(yùn)用有序聚類分析法對花園口水文站的年徑流量時間序列進(jìn)行突變點(diǎn)分析,計(jì)算與τ對應(yīng)的Sn（τ),點(diǎn)繪Sn（τ)變化曲線（見圖 2)。 由圖 2可知,Sn（τ)最小值出現(xiàn)在1990—1991年度,該年度為最可能的突變點(diǎn)。1933—1934年度Sn（τ)局部最小,故該年度也可能是突變點(diǎn)。因此,花園口水文站年徑流量可能在1990—1991年、1933—1934年發(fā)生了突變。

圖2 花園口水文站年徑流有序聚類S變化曲線

2.2.2 Mann-Kendall突變檢驗(yàn)法

Mann-Kendall突變檢驗(yàn)法廣泛應(yīng)用于非參數(shù)檢驗(yàn),是時間序列分析常用方法之一［10］。該方法由Mann和Kendall于1954年提出,近年來被廣泛應(yīng)用于分析徑流、氣溫、降水和水質(zhì)等要素時間序列的突變情況。Mann-Kendall突變檢驗(yàn)法不要求被分析樣本遵從一定分布,也不受其他異常值的干擾,適用于氣象、水文等非正態(tài)分布數(shù)據(jù),計(jì)算十分簡便。應(yīng)用Mann-Kendall突變檢驗(yàn)法進(jìn)行序列突變檢驗(yàn)時,對時間序列xt（t＝1,2,…,n),構(gòu)造一秩序列［11］：

其中

定義統(tǒng)計(jì)量：

其中

式中：E（Sk)和Var（Sk)分別為累計(jì)數(shù)Sk的均值和方差。

將時間序列xt（t＝1,2,…,n)按降序排列,再按式（8)計(jì)算,同時使

當(dāng)k＝1時,UF1為0,服從標(biāo)準(zhǔn)正態(tài)分布。

運(yùn)用Mann-Kendall突變檢驗(yàn)法檢驗(yàn)時間序列的突變性時,將UFk和UBk兩條統(tǒng)計(jì)量序列曲線和顯著性水平為0.05的正態(tài)分布值U0.05＝±1.96兩條臨界線繪制在同一張圖上。如果UFk和UBk兩條曲線在兩條臨界線U0.05＝±1.96之間相交,則說明該序列有突變,交點(diǎn)對應(yīng)的時刻便是突變開始的時間。按照上述方法,得出花園口水文站年徑流量在1922—1923年度、1993—1994年度可能發(fā)生了突變（見圖3)。

圖3 花園口水文站年徑流量Mann-Kendall突變檢驗(yàn)曲線

2.2.3 滑動t檢驗(yàn)法

滑動t檢驗(yàn)法是根據(jù)檢驗(yàn)樣本均值存在的差異是否顯著的方法來檢驗(yàn)突變的,其基本思想是把時間序列中兩段子序列的均值是否存在顯著差異,看成來自總體的均值有無顯著差異來進(jìn)行檢驗(yàn)。如果兩段子序列的均值之間的差異超過了一定的顯著性水平,就可認(rèn)為有突變發(fā)生。

設(shè)滑動點(diǎn)前后兩個序列總體的分布函數(shù)分別為F1（x)和F2（x),從中分別抽取容量為n1和n2的兩個樣本,檢驗(yàn)原假設(shè),若F1（x)＝ F2（x),則該抽樣分布函數(shù)為

其中

式中：ˉx1、ˉx2分別為前后兩個樣本的均值；、分別為兩個樣本的方差；Sw為兩個樣本方差的線性組合。

T服從t（n1+n2-2)分布,與顯著性水平 α＝0.05時在 t分布表中查得的臨界值 t1-α／2對比,若 ｜T ｜＞t1-α／2,拒絕原假設(shè),則兩子序列存在顯著性差異；反之,兩子序列不存在顯著性差異［12］。本文子序列長度n1和n2取8,選取α＝0.01顯著性水平檢驗(yàn)變異程度。運(yùn)用滑動t檢驗(yàn)法對花園口水文站的年徑流量進(jìn)行突變分析,檢驗(yàn)樣本均值存在的差異是否顯著,結(jié)果見圖4。由圖4可知,花園口水文站年徑流量在1929—1939年、1968年、1985年、1988—1991年、2002年發(fā)生變異,其中1930—1934年、1989—1991年為0.10顯著性水平下的變異點(diǎn)。

2.2.4 突變點(diǎn)檢驗(yàn)結(jié)果及成因分析

對以上3種檢驗(yàn)法得到的突變點(diǎn)進(jìn)行匯總,見表2。

表2 年徑流量突變點(diǎn)檢驗(yàn)結(jié)果

由表2可以看出,花園口水文站年徑流量在水文年1933—1934年度、1989—1990年度發(fā)生了突變。因此,將花園口的徑流時間序列分為1919—1933年、1934—1989年、1990—2016年3個階段。

徑流形成過程是多重因素相互作用的復(fù)雜自然現(xiàn)象。資料顯示,黃河流域1933年發(fā)生特大洪水,受災(zāi)面積超過1萬km2,而1972—1989年黃河連年斷流。20世紀(jì)30年代,流域用水規(guī)模不大,自然因素起主導(dǎo)作用,1933年突變的主要因素為全球大尺度的降水、氣溫變化；80—90年代用水量激增,流域尺度的小循環(huán)使得徑流時間序列發(fā)生突變,1989年的突變現(xiàn)象主要為人為因素導(dǎo)致的用水量增加和流域下墊面改變。

2.3 趨勢分析

趨勢分析即分析水文時間序列變化總體的傾向［13］,分析結(jié)果一般有增大、減小、不變3種。本文采用線性趨勢回歸檢驗(yàn)法、Mann-Kendall趨勢檢驗(yàn)法、Spearman秩次相關(guān)檢驗(yàn)法對花園口水文站徑流量變化趨勢進(jìn)行分析。

2.3.1 線性趨勢回歸檢驗(yàn)法

對于檢驗(yàn)自變量與因變量之間是否存在線性關(guān)系,序列xt（t＝1,2,…,n)的線性趨勢成分Yt與時序t的線性回歸方程為

其中a和b的估計(jì)式為

式中：Yt為因變量；a為常數(shù)；b為回歸系數(shù)；＾a和＾b分別為a和b的估計(jì)值；t為時序；ξt為殘差項(xiàng)；ˉx、ˉt分別為xt和t的均值。

構(gòu)造統(tǒng)計(jì)量：

式中：d為服從t分布的函數(shù)；為的方差為殘差項(xiàng)平方和與自由度的比值。

式（16)服從自由度為（n-2)的t分布［14］。 假設(shè)原序列無趨勢性,按照式（16)～式（19)計(jì)算d,選擇顯著水平α＝0.05,在t分布表中查出臨界值與計(jì)算的d值比較,當(dāng)時,拒絕原假設(shè),序列趨勢性顯著。按照上述方法對突變點(diǎn)劃分的3個階段進(jìn)行線性趨勢回歸檢驗(yàn),結(jié)果見表3。

表3 年徑流量線性趨勢回歸檢驗(yàn)結(jié)果

由表3可知,數(shù)據(jù)序列總體變化趨勢較為顯著,通過0.01的顯著性檢驗(yàn),其中：1919—1933年變化趨勢相對顯著,通過0.05的顯著性檢驗(yàn)；1934—1989年、1990—2016年變化趨勢平穩(wěn)。

2.3.2 Mann-Kendall趨勢檢驗(yàn)法

運(yùn)用Mann-Kendall趨勢檢驗(yàn)法進(jìn)行趨勢分析時,同樣不要求被分析樣本遵從一定分布,該方法是關(guān)于觀測值序列的秩次和時序的秩次相關(guān)檢驗(yàn),原假設(shè)H0為時間序列xt（t＝1,2,…,n)n個獨(dú)立的隨機(jī)變量同分布的樣本,Mann-Kendall檢驗(yàn)法中的統(tǒng)計(jì)變量S計(jì)算公式為

其中

式中：Ri和Rj分別為xj和xi的秩次；n為序列長度。

當(dāng)n＞8,假設(shè)實(shí)測數(shù)據(jù)服從獨(dú)立且同分布,統(tǒng)計(jì)量S服從正態(tài)分布,其均值和方差滿足：

標(biāo)準(zhǔn)正態(tài)分布統(tǒng)計(jì)量計(jì)算式為

在趨勢檢驗(yàn)中,選擇顯著性水平α＝0.10,當(dāng)｜Z｜＜時,接受原假設(shè),即趨勢不顯著；反之,則趨勢顯著。統(tǒng)計(jì)變量Z＞0時,序列存在上升趨勢；Z＜0時,序列存在下降趨勢。另選取顯著性水平α＝0.05及α＝0.01檢驗(yàn)趨勢顯著程度。根據(jù)上述方法分析花園口水文站年徑流量變化趨勢,結(jié)果見表4。

表4 年徑流量Mann-Kendall趨勢檢驗(yàn)結(jié)果

從表4可以看出,數(shù)據(jù)序列總體趨勢顯著,通過0.05的顯著性檢驗(yàn),統(tǒng)計(jì)值小于零,總體呈現(xiàn)下降趨勢,其中：1919—1933年變化趨勢相對顯著,通過0.10的顯著性檢驗(yàn),表現(xiàn)為下降趨勢；1934—1989年、1990—2016年變化趨勢平穩(wěn)。

2.3.3 Spearman秩次相關(guān)檢驗(yàn)法

Spearman秩次相關(guān)檢驗(yàn)法是通過分析時間序列xt與其時序t的相關(guān)性來檢驗(yàn)其趨勢性的方法［15］。xt用其秩次Rt（Rt為將xt從大到小排序所對應(yīng)的序號)代表,秩次相關(guān)系數(shù)計(jì)算式為

式中：dt為變量從大到小排序的序號與變量對應(yīng)的時序之間的差值,dt＝Rt-t。

秩次Rt與時序t越接近,則dt的值就越小,秩次相關(guān)系數(shù)r就越大,趨勢性就越顯著。相關(guān)系數(shù)是否異于0,用t檢驗(yàn)法,統(tǒng)計(jì)量計(jì)算式為

式（26)服從自由度為（n-2)的t分布。假設(shè)序列無趨勢,將計(jì)算出的T與在t分布表中查得的顯著性水平α＝0.05對應(yīng)的臨界值則拒絕原假設(shè),序列趨勢顯著；反之,序列趨勢不顯著,另選取顯著性水平α＝0.01檢驗(yàn)趨勢顯著程度。根據(jù)上述方法分析花園口水文站年徑流量變化趨勢,結(jié)果見表5。

表5 年徑流量Spearman秩次檢驗(yàn)結(jié)果

從表5可以看出,數(shù)據(jù)序列總體趨勢顯著,通過0.10的顯著性檢驗(yàn)。根據(jù)對突變點(diǎn)劃分的3個階段分析,趨勢表現(xiàn)平穩(wěn),均未通過0.10的顯著性檢驗(yàn)。

2.3.4 趨勢綜合分析

對以上趨勢檢驗(yàn)分析方法得到的結(jié)果進(jìn)行綜合,假如某種方法檢驗(yàn)趨勢顯著,則記為“1”,否則記為“-1”,將各種檢驗(yàn)方法的加和作為徑流時間序列的綜合趨勢（見表6)。

表6 年徑流量各方法檢驗(yàn)結(jié)果

從表6可以看出：總體趨勢較為顯著；1919—1933年減小趨勢顯著,兩種方法的結(jié)果均通過0.10的顯著性檢驗(yàn)；1934—1989年、1990—2016年的趨勢較為平穩(wěn)。綜合以上分析,花園口年徑流量的突變形式為跳躍突變,1934—1989年的平均年徑流量比1919—1933年增長122億m3,而1990—2016年的平均年徑流量顯著下降,比1934—1989年約減小163億m3（見圖5)。

圖5 花園口水文站年徑流量變異前后綜合趨勢檢驗(yàn)

3 結(jié) 論

（1)綜合有序聚類分析法、Mann-Kendall突變檢驗(yàn)法、滑動t檢驗(yàn)法分析結(jié)果,花園口徑流時間序列在1933—1934年、1989—1990年發(fā)生了突變,突變原因前者主要為自然因素（氣溫、降水)大尺度的變化,后者主要為人為因素（用水增加、下墊面改變)流域尺度的變化。

（2)綜合線性趨勢回歸檢驗(yàn)法、Mann-Kendall趨勢檢驗(yàn)法、Spearman秩次相關(guān)檢驗(yàn)法分析結(jié)果,花園口徑流時間序列在1919—1933年下降趨勢顯著（通過了0.10的顯著性檢驗(yàn))；1934—1989年、1990—2016年徑流量變化相對穩(wěn)定。徑流時間序列的突變類型為跳躍突變,按突變點(diǎn)劃分的3個階段年徑流量平均值相差較大。

黃河流域面積較大,河川徑流的演變比較復(fù)雜,在以后的研究中,應(yīng)結(jié)合上中游其他重要水文站,分析黃河不同河段的徑流演變特性。