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

鮑振鑫，張建云，王國慶，劉翠善，嚴小林，劉 晶，劉 悅

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

0 引 言

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

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

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

1 研究區與資料

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

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

2 研究方法

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

2.1 降水-徑流雙累積曲線

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

(1)

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

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

2.2 有序聚類突變檢測方法

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

(2)

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

2.3 Mann-Kendall突變檢測方法

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

(3)

(4)

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

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

(5)

則Mann-Kendall統計量UFi可由下式估計：

(6)

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

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

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

(7)

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

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

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

3 結 果

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

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

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

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

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

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

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

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

觀臺水文站實測年徑流序列的Mann-Kendall突變檢測結果見圖5。從圖5可以看出，Mann-Kendall統計量為-6，表明年徑流呈現出顯著的下降趨勢，達到了1%的顯著性水平。UF和UB曲線的交點位于1973和1974年之間，表明年徑流序列的突變點發生在1974年，則可以將整個序列分成1951-1973和1974-2015兩段。這一結果略早于前兩種方法檢測的突變點，但是仍處于年徑流持續下降的20世紀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效率系數超過了0.7，相對誤差小于5%，可用于觀臺站的徑流模擬。觀臺水文站實測年徑流序列基于水文模擬的突變檢測結果見圖7。從圖7可以看出，模擬年徑流和實測年徑流的累積差值在20世紀50年代和60年代圍繞0值附近波動，在1971年以后持續偏大。這表明從70年代初開始，由于人類活動的影響，實測徑流小于模擬徑流，即降水-徑流關系開始發生改變。這與實測徑流下降的起始時刻相吻合。

圖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 多種突變檢測方法對比分析

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

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

4 主要結論

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

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

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

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

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

□