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

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

0 引 言

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

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

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

1 研究區(qū)與資料

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

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

2 研究方法

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

2.1 降水-徑流雙累積曲線

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

(1)

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

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

2.2 有序聚類突變檢測(cè)方法

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

(2)

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

2.3 Mann-Kendall突變檢測(cè)方法

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

(3)

(4)

式中：xi是樣本長(zhǎng)度為n的水文序列。

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

(5)

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

(6)

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

2.4 基于水文模擬的突變檢測(cè)方法

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

(7)

式中：QSi和QOi分別為模擬徑流和實(shí)測(cè)徑流。

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

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

3 結(jié) 果

3.1 觀臺(tái)水文站年徑流序列變化特征

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

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

3.2 觀臺(tái)水文站年徑流序列突變檢測(cè)分析

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

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

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

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

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

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

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

圖6 觀臺(tái)水文站實(shí)測(cè)與模擬流量過(guò)程Fig.6 The observed and simulated streamflow at the Guantai hydrologic station

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

3.3 多種突變檢測(cè)方法對(duì)比分析

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

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

4 主要結(jié)論

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

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

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

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

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

□