不同水文序列突變檢測方法在漳河觀臺站徑流分析中的對比研究

鮑振鑫，張建云，王國慶，劉翠善，嚴(yán)小林，劉 晶，劉 悅

(1.南京水利科學(xué)研究院 水文水資源與水利工程科學(xué)國家重點實驗室，南京 210029； 2.水利部應(yīng)對氣候變化研究中心，南京 210029；3.河海大學(xué)，南京 210098)

0 引 言

隨機性是水文序列的重要特性之一。受多種因素的影響和制約，長時間尺度水文序列的演變規(guī)律表現(xiàn)出顯著的隨機性，很難用確定的數(shù)學(xué)物理方程來描述[1]。統(tǒng)計方法是研究水文序列隨機性的重要工具之一。利用統(tǒng)計方法分析水文序列的均值、方差等統(tǒng)計特征，有一個基本前提，即統(tǒng)計樣本服從獨立同分布[2]。但是近些年來由于環(huán)境變化的影響，水文序列的一致性受到破壞，對傳統(tǒng)統(tǒng)計方法研究隨機水文特征帶來了挑戰(zhàn)[3]。因此，識別水文序列的突變特征，尋找可用于統(tǒng)計分析的水文序列段，對于科學(xué)認(rèn)識水文要素的演變規(guī)律，分析環(huán)境變化對水文過程的影響，具有重要的科學(xué)意義和實用價值。

突變檢測方法包括參數(shù)檢驗法和非參數(shù)檢驗法等，被廣泛應(yīng)用于相關(guān)流域水文序列的突變點分析，取得了大量的研究成果[4-6]。例如，田小靖等[7]利用有序聚類分析法、Mann-Kendall檢驗法、Pettitt法等9種突變檢測方法，研究了黃河中游頭道拐站和龍門站的輸沙變化。袁滿[8]等同時考慮類內(nèi)部的離差較小和類間的離差較大原則，提出了改進的有序聚類分析法，識別了年平均流量序列的突變點，結(jié)果更合理。梁欣陽等[9]利用秩和檢驗法和TFPW-MK檢驗法分析了黃河上游蘭州站的水文序列突變特征，結(jié)果表明1985年前后水文指標(biāo)發(fā)生了顯著變化。黃麗娜[10]利用空間差異突變診斷模型研究了淮河上游降水和徑流的突變特征，結(jié)果表明降水和徑流突變發(fā)生在20世紀(jì)80年代中期。張敬平[11]等利用漳澤水庫的降水和徑流序列，對比分析了有序聚類法和啟發(fā)式分割算法在水文序列突變檢測中的應(yīng)用，結(jié)果表明啟發(fā)式分割法能更好地適用于非線性、非平穩(wěn)水文序列的變異診斷。

以往主要側(cè)重于利用突變檢測方法分析水文序列的突變點，關(guān)于多種突變檢測方法有何異同，特別是不同方法的適用范圍，及其檢測結(jié)果的物理意義有何異同等方面的研究較少。本文主要對比分析不同水文序列突變檢測方法，研究其適用范圍及突變點的物理意義，為水文序列演變規(guī)律分析提供技術(shù)支撐。

1 研究區(qū)與資料

以海河流域漳河水系觀臺水文站以上流域為研究對象，流域概況見圖1。流域地處中國北方，位于北緯35.9°～37.6°，東經(jīng)112.4°～114°，地跨山西、河北和河南三省，流域面積17 800 km2。流域以山區(qū)地貌為主，處于半濕潤區(qū)，屬于溫帶大陸性季風(fēng)氣候，四季分明，多年平均氣溫13.2 ℃，降水560 mm，徑流深53 mm，降水-徑流系數(shù)小于0.1。流域內(nèi)支流眾多，水系呈扇形分布，有北部的清漳河與南部的濁漳河兩大水系。流域內(nèi)有漳澤等3座大型水庫，12座中型水庫和93座小型水庫，總庫容約14 億m3，并建有4個大型跨流域引水工程[12]。利用觀臺水文站觀測的1951-2015年徑流序列，分析漳河流域水文特征的突變性。年徑流數(shù)據(jù)來源于水利部刊印的中國水文年鑒，該資料經(jīng)過水文資料整編等多重手續(xù)，結(jié)果可靠。

圖1 觀臺水文站以上漳河流域概況Fig.1 The basic information of the Zhanghe River Basin above the Guantai hydrologic station

2 研究方法

水文序列的突變檢測方法種類繁多，包括經(jīng)驗方法、統(tǒng)計學(xué)方法和水文模擬法等，其中統(tǒng)計學(xué)方法種類最多，應(yīng)用較廣泛[13]。本文選用常用的經(jīng)驗方法中的降水-徑流雙累積曲線，統(tǒng)計學(xué)方法中的有序聚類突變檢測方法和Mann-Kendall突變檢測方法，以及水文模擬法中的基于水文模擬的突變檢測方法等多種方法對比分析水文序列突變檢測結(jié)果的異同。

2.1 降水-徑流雙累積曲線

降水-徑流雙累積曲線是一種廣泛用于分析流域降水-徑流特性是否發(fā)生變化的經(jīng)驗方法[14]，其橫坐標(biāo)為累積年降水量(PSt)，縱坐標(biāo)為累積年徑流量(RSt)，計算公式如下：

(1)

式中：P和R分別為年降水量和年徑流量；n為序列長度。

通過查看降水-徑流雙累積曲線是否發(fā)生轉(zhuǎn)折來判斷序列的突變點。一般為了排除偶然因素的干擾，只有在雙累積曲線的坡度變化顯著，且轉(zhuǎn)折后有連續(xù)5年以上的觀測資料時，才判定降水-徑流序列發(fā)生了突變。該方法主要依賴于專家經(jīng)驗判斷，未通過統(tǒng)計變量檢測。

2.2 有序聚類突變檢測方法

有序聚類突變檢測方法最早由丁晶等人[15]提出，廣泛應(yīng)用于水文序列的突變分析。該方法的基本思想是，有一個樣本長度為n的水文序列，首先任意假定某一點(τ)為突變點，將水文序列分成兩段，則這兩段水文序列分別代表兩類水文樣本，統(tǒng)計兩段樣本的離差平方和(Snτ)，其計算公式為：

(2)

根據(jù)同類水文樣本間的方差較小，類之間的方差較大這一原則，找到Snτ的最小值，則相應(yīng)的為水文樣本的突變點。

2.3 Mann-Kendall突變檢測方法

Mann-Kendall檢測方法是一種非參數(shù)檢驗方法，最早由H B Mann和M G Kendall提出[16]。首先根據(jù)樣本序列計算xi>xj(i>j)的個數(shù)(Sk)，其計算公式為：

(3)

(4)

式中：xi是樣本長度為n的水文序列。

再統(tǒng)計Si的期望值[E(Si)]和方差[Var(Si)]：

(5)

則Mann-Kendall統(tǒng)計量UFi可由下式估計：

(6)

對于某一個顯著性水平α，查得正態(tài)分布臨界值Uα/2。如果|UFi|

2.4 基于水文模擬的突變檢測方法

考慮到傳統(tǒng)基于統(tǒng)計方法的突變檢測技術(shù)只針對單一的水文序列如何改變，即將人類活動和氣候變化對水文過程的影響混合在一起，為了更科學(xué)地分析人類活動改變降水-徑流特性，進而引起水文過程改變的機理，Wang等[17]提出了基于水文模擬的突變檢測技術(shù)。首先根據(jù)研究流域的實際情況，選擇未受人類活動影響的天然時期，利用天然時期的水文數(shù)據(jù)率定水文模型參數(shù)，滿足天然時期的觀測值和模擬值沒有系統(tǒng)偏差，即模擬值圍繞觀測值波動。天然時期率定的模型參數(shù)反映流域在天然時期未受人類活動干擾的下墊面特性。利用天然時期率定的模型參數(shù)，模擬整個時期的徑流過程。如果模擬徑流與實測徑流出現(xiàn)系統(tǒng)偏差，則將系統(tǒng)偏差的起始點認(rèn)為是人類活動對流域水文過程干擾的起始時刻，即流域降水-徑流特性發(fā)生改變的突變點。定義模擬徑流和實測徑流的偏差(Ki)與累積偏差(Km)如下：

(7)

式中：QSi和QOi分別為模擬徑流和實測徑流。

在天然時期，累積偏差(Km)在0值附近波動，將其持續(xù)大于0或者小于0的起始點判定為流域降水-徑流特性發(fā)生改變的突變點。

目前水文模型有很多種，可以根據(jù)流域特征和徑流模擬的需要來選擇合適的水文模型。本文選用VIC(Variable Infiltration Capacity)模型來模擬觀臺水文站的徑流過程。VIC模型是一個基于正交網(wǎng)格的分布式水文模型，它能夠同時模擬地表間的能量平衡和水量平衡，廣泛應(yīng)用于全球很多流域的徑流模擬[18, 19]。VIC模型將土層劃分為3層，總徑流由上兩層土壤產(chǎn)生的地表徑流和第3層土壤產(chǎn)生的基流兩部分組成。借用了新安江模型中的蓄水容量曲線來描述土壤含水量的空間分布不均勻性，計算地表徑流過程[20]。基流部分采用Arno模型中的原理來計算，即當(dāng)土壤含水量在某一閾值以下時，基流是線性消退的；而高于此閾值時，基流過程是非線性的[21]。在模型計算時，首先分別在每個網(wǎng)格內(nèi)獨立進行降水-蒸發(fā)-產(chǎn)流過程的計算，然后再統(tǒng)一匯流到流域出口斷面形成流量過程。VIC模型中共有6個參數(shù)需要利用觀測的流量資料來進行率定。包括第2層和第3層的土層厚度，描述土壤蓄水容量曲線的形狀參數(shù)，以及控制基流的3個參數(shù)。模型參數(shù)的率定主要考慮2個目標(biāo)函數(shù)：Nash-Sutcliffe效率系數(shù)和相對誤差。其中，Nash-Sutcliffe效率系數(shù)反映了模擬的流量過程和觀測的流量過程之間的吻合程度，其值越接近1，表示模擬效果越好；而相對誤差是一個水量平衡指標(biāo)，它反映了模擬總徑流量和觀測總徑流量之間的相對誤差，其值越接近0，則表示模擬效果越好。

3 結(jié) 果

3.1 觀臺水文站年徑流序列變化特征

1951-2015年觀臺水文站觀測的年徑流序列見圖2。從圖2可以看出，1951年以來，觀臺水文站觀測的年徑流序列呈現(xiàn)出持續(xù)性的下降趨勢，其中最大值發(fā)生在1963年，為259 mm，最小值發(fā)生在1999年，僅為1.67 mm，最大值是最小值的155倍。20世紀(jì)50和60年代處于豐水期，年徑流在100 mm左右波動；從70年代初開始，年徑流持續(xù)下降，到70年代末80年代初，年徑流降低到10 mm以下；隨后該流域進入連續(xù)枯水期，年徑流在16 mm左右波動。從觀臺水文站實測年徑流序列的演變特征來看，徑流序列的突變點發(fā)生在20世紀(jì)70年代。

圖2 觀臺水文站年徑流序列Fig.2 The observed annual streamflow at the Guantai hydrologic station

3.2 觀臺水文站年徑流序列突變檢測分析

觀臺水文站以上流域年降水量和徑流量雙累積曲線見圖3。從圖3可以看出，降水-徑流關(guān)系的突變點發(fā)生在1978年。在同樣的降水條件下，1978年以后的年徑流量明顯小于1977年之前的徑流量。例如1951-1977年，年降水徑流系數(shù)為0.158，但是1978-2015年的降水徑流系數(shù)僅為0.034。除了1978年這一明顯的轉(zhuǎn)折點以外，在1971年降水量和徑流量雙累積曲線開始有所轉(zhuǎn)折，但是不是很顯著。這表明降水-徑流關(guān)系從1971年開始有所變化，但是變化程度不大；從1978年開始降水-徑流關(guān)系發(fā)生了突變。

圖3 觀臺水文站年降水量和徑流量雙累積曲線Fig.3 The double-accumulation curve of precipitation and runoff at the Guantai hydrologic station

觀臺水文站實測年徑流序列的有序聚類突變檢測結(jié)果見圖4。從圖4可以看出，年徑流序列離差平方和的最小值位于1977年，則可以將整個序列分成1951-1977和1978-2015兩段。這一結(jié)果和降水-徑流雙累積曲線的分析結(jié)果相一致。離差平方和的次小值是1973年，其值與1977年的值相差很小，而且從20世紀(jì)60年代后期到70年代末，年徑流序列離差平方和的值都比較小，在谷值附近波動。這表明年徑流序列的突變不一定是在某一個年份突然發(fā)生的，而是在70年代緩慢進行的。這一結(jié)論和年徑流的演變特征相吻合，即在20世紀(jì)70年代，徑流持續(xù)下降，而不是在某一年突然躍變。

圖4 觀臺水文站年徑流序列有序聚類突變檢測Fig.4 The breakpoint of annual streamflow at the Guantai hydrologic station detected by the orderly clustering breakpoint detection method

觀臺水文站實測年徑流序列的Mann-Kendall突變檢測結(jié)果見圖5。從圖5可以看出，Mann-Kendall統(tǒng)計量為-6，表明年徑流呈現(xiàn)出顯著的下降趨勢，達到了1%的顯著性水平。UF和UB曲線的交點位于1973和1974年之間，表明年徑流序列的突變點發(fā)生在1974年，則可以將整個序列分成1951-1973和1974-2015兩段。這一結(jié)果略早于前兩種方法檢測的突變點，但是仍處于年徑流持續(xù)下降的20世紀(jì)70年代之中。

圖5 -觀臺水文站年徑流序列Mann-Kendall突變檢測Fig.5 The breakpoint of annual streamflow at the Guantai hydrologic station detected by the Mann-Kendall’s test methodology

觀臺水文站實測與VIC模型模擬的1955-1970年月徑流過程見圖6。從圖6可見VIC模型具有較好的徑流模型效果，Nash-Sutcliffe效率系數(shù)超過了0.7，相對誤差小于5%，可用于觀臺站的徑流模擬。觀臺水文站實測年徑流序列基于水文模擬的突變檢測結(jié)果見圖7。從圖7可以看出，模擬年徑流和實測年徑流的累積差值在20世紀(jì)50年代和60年代圍繞0值附近波動，在1971年以后持續(xù)偏大。這表明從70年代初開始，由于人類活動的影響，實測徑流小于模擬徑流，即降水-徑流關(guān)系開始發(fā)生改變。這與實測徑流下降的起始時刻相吻合。

圖6 觀臺水文站實測與模擬流量過程Fig.6 The observed and simulated streamflow at the Guantai hydrologic station

圖7 觀臺水文站年徑流序列基于水文模擬的突變檢測Fig.7 The breakpoint of annual streamflow at the Guantai hydrologi c station detected by the hydrological simulation based method

3.3 多種突變檢測方法對比分析

根據(jù)上文分析結(jié)果，利用降水-徑流雙累積曲線、有序聚類突變檢測方法、Mann-Kendall突變檢測方法、基于水文模擬的突變檢測方法等多種方法，分析得到觀臺水文站年徑流序列的突變點分別是1978、1978、1974和1971年?？梢娀谒哪M的突變檢測方法分析得到的突變點最早，Mann-Kendall突變檢測方法結(jié)果其次，而降水-徑流雙累積曲線和有序聚類突變檢測方法的結(jié)果最晚。

氣候要素的分析結(jié)果表明，漳河流域降水在20世紀(jì)50和60年代較大，70-90年代持續(xù)下降，但是在21世紀(jì)有所回升。從60年代初開始，流域內(nèi)修建了大量的水庫和引調(diào)水工程。1978年中國開始實行改革開放，經(jīng)濟快速發(fā)展，對水資源的需求迅速增加，極大地改變了流域的降水-徑流特性。結(jié)合徑流的演變特征和流域內(nèi)的氣候變化與人類活動情勢，4種突變檢測方法檢測的徑流突變點物理意義有所差別。降水-徑流雙累積曲線檢測的突變點反映的是流域降水-徑流關(guān)系發(fā)生改變的時刻，在本案例中從1971年開始流域降水徑流關(guān)系開始發(fā)生變化但是程度不大，從1978年開始發(fā)生了顯著變化。有序聚類突變檢測方法檢測的突變點反映的是徑流序列的聚類特征，在本案例中以1978年為界分成前后兩類，但是從20世紀(jì)60年代后期到70年代末的離差平方和都較小，表明以其中任意一年為突變點都具有一定的意義。Mann-Kendall突變檢測方法檢測的突變點反映的是徑流序列的突變特征，在本案例中其分析結(jié)果和有序聚類突變檢測方法的次小值相一致?；谒哪M的突變檢測方法檢測的突變點反映的是人類活動對流域降水徑流關(guān)系干擾的起始點，即本案例中的1971年，因此其結(jié)果早于其他方法的檢驗結(jié)果。綜合分析四種突變檢測方法計算結(jié)果的物理意義和流域?qū)嶋H氣候、水文、人類活動的狀況相符合。

4 主要結(jié)論

科學(xué)診斷水文序列的突變特性，對于認(rèn)識水文循環(huán)的演變規(guī)律，研究環(huán)境變化對水文過程的影響具有重要的科學(xué)意義和實用價值?；谡暮佑^臺水文站觀測的1951-2015年徑流序列，利用降水-徑流雙累積曲線、有序聚類突變檢測方法、Mann-Kendall突變檢測方法、基于水文模擬的突變檢測方法等多種方法分析了其突變特征。對比分析了各種突變檢測方法的異同及適用范圍。主要結(jié)論如下。

(1)1951年以來，觀臺水文站實測年徑流序列呈現(xiàn)出持續(xù)性的下降趨勢，20世紀(jì)50和60年代處于豐水期，從70年代初開始年徑流持續(xù)下降，80年代進入連續(xù)枯水期。

(2)降水-徑流雙累積曲線、有序聚類突變檢測方法、Mann-Kendall突變檢測方法和基于水文模擬的突變檢測方法等多種方法，檢測得到觀臺水文站年徑流序列的突變點分別是1978、1978、1974和1971年。

(3)相比而言，降水-徑流雙累積曲線檢測的突變點反映的是流域降水-徑流關(guān)系發(fā)生明顯改變的時刻；有序聚類突變檢測方法檢測的突變點反映的是徑流序列的聚類特征；Mann-Kendall突變檢測方法檢測的突變點反映的是徑流序列的綜合突變特征；基于水文模擬的突變檢測方法檢測的突變點反映的是人類活動對流域降水徑流關(guān)系干擾的起始點。

在分析流域水文序列的突變特性時，應(yīng)根據(jù)研究的需要合理選用相應(yīng)的突變檢測方法，科學(xué)識別水文序列的變異情勢。

□