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

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

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

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

1 數(shù)據(jù)來(lái)源

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

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

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

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

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

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

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

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

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

表1 各時(shí)段徑流時(shí)間序列特征值

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

2.2 突變分析

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

2.2.1 有序聚類(lèi)分析法

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

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

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

總離差平方和為

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

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

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

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

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

其中

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

其中

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

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

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

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

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

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

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

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

其中

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

T服從t（n1+n2-2)分布,與顯著性水平 α＝0.05時(shí)在 t分布表中查得的臨界值 t1-α／2對(duì)比,若 ｜T ｜＞t1-α／2,拒絕原假設(shè),則兩子序列存在顯著性差異；反之,兩子序列不存在顯著性差異［12］。本文子序列長(zhǎng)度n1和n2取8,選取α＝0.01顯著性水平檢驗(yàn)變異程度。運(yùn)用滑動(dòng)t檢驗(yàn)法對(duì)花園口水文站的年徑流量進(jìn)行突變分析,檢驗(yàn)樣本均值存在的差異是否顯著,結(jié)果見(jiàn)圖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é)果及成因分析

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

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

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

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

2.3 趨勢(shì)分析

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

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

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

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

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

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

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

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

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

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

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

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

其中

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

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

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

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

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

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

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

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

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

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

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

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

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

2.3.4 趨勢(shì)綜合分析

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

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

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

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

3 結(jié) 論

（1)綜合有序聚類(lèi)分析法、Mann-Kendall突變檢驗(yàn)法、滑動(dòng)t檢驗(yàn)法分析結(jié)果,花園口徑流時(shí)間序列在1933—1934年、1989—1990年發(fā)生了突變,突變?cè)蚯罢咧饕獮樽匀灰蛩兀鉁亍⒔邓?大尺度的變化,后者主要為人為因素（用水增加、下墊面改變)流域尺度的變化。

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

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