秦超瑞，桑宇婷，趙雪花，祝雪萍

（太原理工大學 水利科學與工程學院，山西 太原 030024）

受自然因素及人類活動的影響，河流年徑流量會發(fā)生顯著變化［1-3］，徑流序列發(fā)生突變［4］，導致旱澇災害等一系列水文極值事件頻發(fā)，研究徑流序列的突變規(guī)律可以為流域治理與水資源合理開發(fā)利用提供科學依據(jù)。 目前，水文氣象方面常用的突變檢測方法主要有參數(shù)統(tǒng)計法和非參數(shù)秩檢測統(tǒng)計法兩大類。 其中參數(shù)統(tǒng)計法，如累積距平法［5］、滑動t 檢測法［6-7］及線性回歸法等，適用于滿足一定分布的數(shù)據(jù)序列，但在實際應用中無法對樣本序列的分布形態(tài)做出簡單假設，且統(tǒng)計值會受到原始數(shù)據(jù)序列長度、數(shù)據(jù)異常及缺失的影響；非參數(shù)秩檢測統(tǒng)計法，如Mann-Kendall 檢測法［8］和Pettitt 突變檢測法［9］等，無法去除樣本數(shù)據(jù)所固有的相關屬性，從而存在分析誤差。 多尺度直線擬合法可避免原始數(shù)據(jù)屬性及其固有相關屬性的影響，很好地克服上述方法的缺點，已被應用于電液伺服閥的故障診斷［10］，并取得了很好的效果。 除徑流突變規(guī)律外，突變前兆信號也很重要。 臨界慢化現(xiàn)象由于淥等［11］于1984年提出，具有揭示復雜動力系統(tǒng)是否趨于臨界災變的優(yōu)勢。 2009年，Scheffer 等［12］研究表明，當系統(tǒng)趨近于分岔點時，臨界慢化現(xiàn)象會導致動力學出現(xiàn)3 個可能的早期預警信號（擾動的恢復變慢、方差增大、自相關系數(shù)增大）。 由此為突變前兆信號研究提供了一種新途徑，并被廣泛用于多個領域，如晏銳等［13］將臨界慢化現(xiàn)象用于地震前兆信號檢測；顏鵬程等［14］和吳浩等［15］將臨界慢化現(xiàn)象應用于氣象突變前兆信號研究中，證明臨界慢化現(xiàn)象是氣候突變前的一個早期信號。 但臨界慢化現(xiàn)象在徑流序列中的應用仍較少，因此本文采用多尺度直線擬合法分析汾河上游水文站的徑流序列突變點并根據(jù)系統(tǒng)恢復速率及恢復力、方差及自相關系數(shù)的演變規(guī)律，探討前兆信號發(fā)生年份作為早期預警的可能性。

1 研究內(nèi)容與研究方法

1.1 研究區(qū)域概況及數(shù)據(jù)來源

汾河位于山西省境內(nèi)，是山西省第一大河流，流域面積為39471 km2。 汾河上游段從河源寧武縣東寨鎮(zhèn)管涔山到太原蘭村，流域面積為7705 km2，是太原市主要的水源涵養(yǎng)區(qū)和水土保持生態(tài)功能區(qū)。 上游流域位于黃土高原地區(qū)，屬于溫帶大陸性季風區(qū)，春冬干旱少雨，夏秋多雨，且大雨、暴雨多集中于7月、8月。 研究區(qū)地處山區(qū)，無大型城市群，城鎮(zhèn)化程度較低，工農(nóng)業(yè)處于落后水平，但礦產(chǎn)資源非常豐富。 20 世紀60年代開始，人類活動對徑流影響加劇，主要體現(xiàn)在水利工程建設及水資源開采。 流域內(nèi)有汾河水庫和汾河二庫兩座大型水庫，其中汾河水庫是山西省最大的水庫，位于汾河干流上游，地處婁煩縣杜交曲鎮(zhèn)，設計總庫容為7.21 億m3，是一座集防洪、灌溉、發(fā)電為一體的水利工程，于1956年攔洪，1961年竣工驗收后正式投入運行。 本文以汾河上游4 個水文站（上靜游站、汾河水庫站、寨上站及蘭村站）1956—2000年的年徑流量序列（資料來自山西省水文局）為研究資料。 4 個水文站年降水量最大值均發(fā)生于1967年，最小值均發(fā)生于1972年。 其中上靜游站位于汾河支流嵐河上，控制流域面積為1140 km2。 汾河水庫站、寨上站及蘭村站位于汾河干流上，控制流域面積分別為5268、6819、7705 km2。 4 個水文站在汾河上游流域內(nèi)的地理位置如圖1所示。

圖1 汾河上游流域水系及水文站分布

1.2 研究內(nèi)容

針對汾河上游流域4 個水文站的年徑流量序列進行突變點及突變前兆信號檢測。 具體分為如下3 個部分：采用多尺度直線擬合法檢測年徑流量序列突變點，根據(jù)系統(tǒng)的恢復速率及恢復力演變規(guī)律檢測突變前兆信號發(fā)生年份，根據(jù)序列的方差及自相關系數(shù)演變規(guī)律檢測突變前兆信號發(fā)生年份。

1.3 研究方法

1.3.1 多尺度直線擬合法

根據(jù)多尺度直線擬合法原理［10］，選取合適的尺度，用多條線段代替原序列，可利用多條擬合線段提取原序列的主要信息，以達到快速有效地反映原序列信息的目的。 初始尺度r1計算公式如下：

式中：n為徑流序列長度。

采用式（1）計算得到初始尺度r1后，可將原序列劃分成m段，每段用di（i取1 ～m）來表示，如最后一段長度小于r1，則將其并入前一段即可。 對于每段序列，均可用最小二乘法進行直線擬合。 根據(jù)最小二乘法原理可求得每段擬合線段的斜率ki，并依次求出相鄰兩段線段斜率差值的絕對值，找出最大值記為emax，說明原序列在該相鄰兩段線段間的變化最大，則突變點一定落在該區(qū)間內(nèi)，此時將r1變?yōu)閞2，r2=E（r1／2）。在r2尺度下對新區(qū)間采用上述方法逼近突變點，直至尺度rn=2 時，此時的區(qū)間長度為3，區(qū)間的中點即為原徑流序列的突變點。

1.3.2 基于系統(tǒng)的恢復速率及恢復力演變規(guī)律檢測突變前兆信號

徑流系統(tǒng)是一個復雜的動力學系統(tǒng)，徑流量在變化環(huán)境中演變時往往會發(fā)生突變，而徑流系統(tǒng)在發(fā)生突變前總是先演變到某個臨界閾值附近，一段時間后系統(tǒng)才發(fā)生突變。 因此，徑流系統(tǒng)在臨界閾值附近的一些特征可作為判斷徑流系統(tǒng)臨近突變的前兆信號。由臨界慢化現(xiàn)象［11］（動力系統(tǒng)由一種相態(tài)轉(zhuǎn)變?yōu)榱硪环N相態(tài)之前，系統(tǒng)在趨近于臨界點附近，表現(xiàn)出漲落幅度增大、擾動的恢復速率變慢以及恢復到舊相位能力變小的現(xiàn)象）可知，系統(tǒng)的恢復速率及恢復力可作為表征徑流系統(tǒng)突變前兆信號的特征量。 構(gòu)建徑流序列混沌模型，假設徑流量變化是從一個平衡態(tài)經(jīng)歷一次突變后向另一個平衡態(tài)轉(zhuǎn)移的過程，系統(tǒng)的恢復速率及恢復力是指徑流序列在受到突變干擾后重新恢復到平衡態(tài)的速率和能力。 系統(tǒng)的恢復速率可表示為

式中：x為徑流狀態(tài)變量；λ為徑流量的增長系數(shù)；λx為徑流量的增長速率；k為徑流量的減少系數(shù)；kx2為徑流量的減少速率；t為時間。

式（2）的物理意義是徑流狀態(tài)變量x隨著時間t增大（減?。┑乃俾?。 該模型中的初始值和參數(shù)λ、k只能取正值，為了便于討論，令λ／k=μ，由式（2）可知，當dx／dt=0 時，x=0 和x=μ為式（2）的兩個解，此時系統(tǒng)的狀態(tài)變量不變，即徑流量達到平衡狀態(tài)。 式（2）等號兩邊同時對t求導，可以得到系統(tǒng)的恢復力：

由式（3）可知，當x=0、x=μ／2 和x=μ時，系統(tǒng)的恢復力為零。

該模型描述的是徑流狀態(tài)變量由初始平衡態(tài)（x=x0）經(jīng)歷一次突變后趨向于另一個平衡態(tài)（x=μ）。 當徑流狀態(tài)變量開始發(fā)生演變時，系統(tǒng)的恢復速率和恢復力也發(fā)生演變，且變化幅度大于徑流狀態(tài)變量的變化幅度。 根據(jù)臨界慢化現(xiàn)象，在徑流狀態(tài)變量趨于突變點之前會出現(xiàn)系統(tǒng)的恢復速率及恢復力變小的現(xiàn)象，表現(xiàn)為存在某一時刻系統(tǒng)的恢復速率及恢復力同時為零，并發(fā)生相位轉(zhuǎn)變，該時刻能很好地反映徑流狀態(tài)變量偏離平衡態(tài)的演變過程，可作為徑流系統(tǒng)發(fā)生突變的前兆信號。

1.3.3 基于方差及自相關系數(shù)演變規(guī)律檢測突變前兆信號

臨界慢化現(xiàn)象除了導致系統(tǒng)的恢復速率及恢復力出現(xiàn)上述演變規(guī)律外，還將導致序列的方差和自相關系數(shù)增大［11］。 因此，徑流系統(tǒng)趨近閾值附近時序列的方差和自相關系數(shù)同時達到極小值并開始逐漸增大，表明徑流系統(tǒng)即將發(fā)生突變，也可作為徑流系統(tǒng)發(fā)生突變的前兆信號。

方差是概率論中用來度量隨機變量與其均值之間偏離程度的指標。 自相關系數(shù)度量的是同一變量在兩個不同時期的相關程度，即自身過去行為對現(xiàn)在的影響，自相關系數(shù)c的定義如下：

式中：j為滯后長度；xi為序列中第i個樣本值； ˉx為序列均值；s為序列均方差。

2 結(jié)果與分析

2.1 年徑流量序列突變點檢測

將多尺度直線擬合法應用于汾河上游4 個水文站年徑流量序列突變點檢測，結(jié)果如圖2 所示。

圖2 年徑流量序列突變點檢測結(jié)果

圖2（a）為上靜游站年徑流量序列的突變點檢測結(jié)果，經(jīng)計算可得分割點3、4 之間線段的斜率和分割點4、5 之間線段的斜率差值絕對值最大，由多尺度直線擬合法原理可知突變點一定落在分割點3 與分割點5 之間，因此將尋求突變點的范圍縮小至分割點3 與分割點5 之間。 采用同上的方法不斷縮小范圍和尺度，直至最后尺度rn=2 時，此時區(qū)間的中點即為原始序列的突變點。 圖中A點為最終檢測到突變點的位置，由此可知上靜游站年徑流量序列突變年份為1975年，年徑流量均值在1975年前后分別為0.57 億m3和0.44 億m3，突變后較突變前減少22.81%。 圖2（b）～（d）為汾河水庫站、寨上站和蘭村站年徑流量序列突變點檢測結(jié)果，經(jīng)計算可得3 個水文站初始尺度下擬合線段斜率差值絕對值最大值均為分割點2、3 之間的線段與分割點3、4 之間的線段之差，因此將尋求突變點的范圍縮小至分割點2 與分割點4 之間，采用同上的方法逐漸逼近突變點，最終A點為檢測到的突變點位置，由此可知汾河水庫站、寨上站和蘭村站年徑流量序列的突變年份均為1966年。 其中汾河水庫站年徑流量均值在1966年前后分別為4.00 億m3和3.21億m3，突變后較突變前減少19.75%；寨上站年徑流量均值在1966年前后分別為4.74 億m3和3.57 億m3，突變后較突變前減少24.68%；蘭村站年徑流量均值在1966年前后分別為4.72 億m3和3.25 億m3，突變后較突變前減少31.14%。 干流上3 個水文站徑流量的減少程度自上而下依次遞增。 通過降水徑流雙累積曲線斜率變化率計算［16］，可得降水和人類活動分別對4 個水文站徑流量序列突變的貢獻率，結(jié)果見表1。

表1 降水和人類活動對突變的貢獻率

由表1可知，汾河干流上的水文站年徑流量序列突變一致性較好，且人類活動是引起干流徑流量序列突變的主要因素。 而支流的年徑流量序列突變一致性較差，降水是引起上靜游站徑流量序列突變的主要因素。 降水是徑流產(chǎn)生的源泉，因此降水變化是影響徑流變化最敏感的因素。 4 個水文站年降水量均于1967年達到最大值，于1972年減到最小值。 由于短時期內(nèi)降水量發(fā)生較大幅度的波動且支流控制的流域面積遠小于干流，因此上靜游站1975年發(fā)生突變，而干流沿程不斷有支流匯入導致干流站控制的流域面積不斷增大以緩沖降水波動帶來的影響，汾河水庫站、寨上站和蘭村站受汾河水庫的調(diào)節(jié)作用使得這次降水量大幅度波動未對干流上的水文站徑流產(chǎn)生突變影響。

人類活動對支流的上靜游站影響較小，對干流徑流變化的影響主要體現(xiàn)在水利工程和水資源開采兩個方面。 水利工程主要為汾河水庫的修建，汾河水庫投入使用后防洪效益顯著，1965年全年來水量為2.01 億m3，而水庫給下游放水6.29 億m3，充分發(fā)揮了水庫的調(diào)節(jié)作用，保證了下游的用水需求［16］。 工業(yè)用水量和農(nóng)業(yè)用水量急劇增加，工業(yè)用水量在1949年約為180萬m3，從汾河水庫1960年投入運行后增至3051.11萬m3，1956年汾河灌區(qū)面積為8.32 萬hm2，在汾河水庫1960年投入運行后灌溉面積增至8.46 萬hm2［17］。因此，汾河水庫下游工農(nóng)業(yè)用水量急劇增加是干流水文站徑流量序列發(fā)生突變的原因之一。 汾河上游水土流失嚴重，汾河水庫的修建使上游來沙大量被攔蓄，下游輸沙量急劇減少（1965年水庫泥沙淤積為0.0125 億m3，1966年為0.1270 億m3）［18］，因此汾河水庫攔蓄泥沙也是造成干流水文站徑流量序列發(fā)生突變的原因之一。 隨著城鎮(zhèn)化程度提高，汾河上游流域內(nèi)地下水過量開采（1949年以前每日不超過1 萬m3，到1965年增至每日32.93 萬m3），導致地下水位下降，深層地下水漏斗從1965年起迅速擴大［19］，而取水口一般設在支流和干流匯合后的河段上，干流上3 個水文站均處于支流和干流匯合后的河段上，因此人類對地下水資源的過度開采也是導致干流徑流量序列發(fā)生突變的原因之一。 干流上3 個水文站徑流量減少程度自上而下依次遞增的原因有兩個：一是汾河水庫修建于汾河干流上游，自上而下對汾河水庫站、寨上站和蘭村站的調(diào)節(jié)作用遞減；二是沿著干流不斷有支流匯入，導致徑流量減少程度增大。

綜上，位于支流的上靜游站徑流量序列1975年發(fā)生突變的主要影響因素為降水，而位于干流的汾河水庫站、寨上站和蘭村站徑流量序列1966年發(fā)生突變的主要影響因素為人類活動。

2.2 基于系統(tǒng)的恢復速率及恢復力演變規(guī)律檢測突變前兆信號結(jié)果分析

為了更好地檢測徑流突變年份，將系統(tǒng)恢復速率及恢復力應用于汾河上游水文站徑流量序列突變前兆信號檢測中，如圖3 所示。 由圖3（a）和（b）可知上靜游站年徑流量序列的系統(tǒng)恢復速率在1973年為零并由正相位向負相位轉(zhuǎn)變（見圖3（a）中A點），恢復力在1973年為零并由正相位轉(zhuǎn)變?yōu)樨撓辔唬ㄒ妶D3（b）中B點），恢復速率和恢復力同時為零且發(fā)生相位轉(zhuǎn)變的點只有1973年，則上靜游站年徑流量序列的突變前兆信號發(fā)生于1973年。 由圖3（c）～（h）可知，汾河水庫站、寨上站和蘭村站年徑流量序列的系統(tǒng)恢復速率均在1958年附近為零并由負相位轉(zhuǎn)變?yōu)檎辔唬ㄒ妶D3（c）（e）（g）中A點），恢復力均在1958年附近為零并由正相位轉(zhuǎn)變?yōu)樨撓辔唬ㄒ妶D3（d）（f）（h）中B點），恢復速率和恢復力同時為零且發(fā)生相位轉(zhuǎn)變的點只有1958年，則干流水文站的突變前兆信號均發(fā)生于1958年。 綜上，根據(jù)徑流量序列的系統(tǒng)恢復速率及恢復力演變規(guī)律，得到上靜游站的突變前兆信號發(fā)生于1973年，汾河水庫站、寨上站和蘭村站的突變前兆信號均發(fā)生于1958年。

圖3 年徑流量序列恢復速率及恢復力演變過程

在適當?shù)臅r間范圍內(nèi)前兆信號發(fā)生年份較突變年份越早，系統(tǒng)的恢復速率及恢復力對引起突變的因子越敏感，可為預防突變發(fā)生提供充足的準備時間，因此預警效果越優(yōu)。 根據(jù)徑流狀態(tài)變量恢復速率及恢復力的演變規(guī)律，得到上靜游站年徑流量序列的突變前兆信號發(fā)生于1973年，比實際突變年份1975年提前2 a。 汾河水庫站、寨上站和蘭村站年徑流量序列的突變前兆信號均發(fā)生在1958年，比實際突變年份1966年提前8 a。 突變前兆信號均早于突變年份10 a 之內(nèi)，說明徑流狀態(tài)變量的恢復速率及恢復力的演變規(guī)律能為徑流量序列提供可靠的突變前兆信號。 由于前兆信號發(fā)生年份在支流上只比突變年份提前2 a，和突變年份較為接近，因此預警效果不如干流。 前兆信號發(fā)生年份在干流上比突變年份提前8 a，因此針對汾河上游流域，基于系統(tǒng)恢復速率及恢復力的演變規(guī)律更適合求取干流徑流量序列的突變前兆信號。 由于干流徑流量序列的突變受人類活動的影響較大，受降水影響較小，因此徑流狀態(tài)變量的恢復速率及恢復力對人類活動干擾更敏感。

2.3 基于方差及自相關系數(shù)演變規(guī)律檢測突變前兆信號結(jié)果分析

基于臨界慢化現(xiàn)象，將方差和自相關系數(shù)應用于汾河上游水文站年徑流量序列突變前兆信號檢測中，結(jié)果如圖4 所示。 由圖4（a）和（b）可知，上靜游站徑流量序列的方差和自相關系數(shù)同時達到極小值的點只有1968年（見圖4（a）和（b）中的A點和B點），之后開始逐漸增大，則上靜游站徑流量序列的突變前兆信號發(fā)生于1968年。 由圖4（c）～（h）可知，汾河水庫站、寨上站和蘭村站的方差和自相關系數(shù)同時達到極小值的點只有1960年（見圖4（c）～（h）中A點和B點），之后開始逐漸增大，則3 個水文站徑流量序列的突變前兆信號均發(fā)生于1960年。 綜上，根據(jù)徑流量序列的方差和自相關系數(shù)演變規(guī)律得到上靜游站徑流量序列的突變前兆信號發(fā)生于1968年，干流3 個水文站徑流量序列的突變前兆信號均發(fā)生于1960年。

圖4 年徑流量序列方差及自相關系數(shù)演變過程

根據(jù)臨界狀態(tài)時徑流量序列的方差及自相關系數(shù)的演變規(guī)律，得到上靜游站徑流量序列突變的前兆信號發(fā)生于1968年，比實際突變年份1975年提前7 a。汾河水庫站、寨上站和蘭村站3 個水文站徑流量序列突變的前兆信號均發(fā)生于1960年，比實際突變年份1966年提前6 a。 突變前兆信號均早于突變年份10 a之內(nèi)，說明徑流量序列的方差及自相關系數(shù)的演變規(guī)律能為徑流量序列提供可靠的突變前兆信號。 前兆信號發(fā)生年份在支流上較突變年份出現(xiàn)更早，預警效果優(yōu)于干流。 因此，基于方差及自相關系數(shù)的演變規(guī)律檢測突變前兆信號更適用于支流徑流量序列。 由于支流徑流量序列的突變受降水影響較大，受人類活動的影響較少，因此徑流量序列的方差和自相關系數(shù)對降水量變化更敏感。

3 討論與結(jié)論

3.1 討論

采用多尺度直線擬合法對汾河上游徑流量序列進行突變點檢測，計算降水和人類活動分別對突變的貢獻率并分析其原因。 根據(jù)臨界慢化現(xiàn)象檢測突變前兆信號發(fā)生年份，探討了前兆信號作為早期預警的可能性。 需要指出的是多尺度直線擬合法只能檢測出序列中突變程度最顯著的一個突變點，在今后的研究工作中，需改進檢測方法，使其可以檢測出序列中所有的突變點。 同時，對臨界慢化現(xiàn)象是否受突變類型及突變強度影響仍然需要進一步研究。

3.2 結(jié)論

（1）位于支流的上靜游站年徑流量序列的突變年份為1975年，年徑流量均值在1975年前后分別為0.57 萬m3和0. 44 萬m3， 突 變 后 較 突 變 前 減 少22.81%，降水對突變的貢獻率為79.67%，是突變的主要影響因素。 位于干流的汾河水庫站、寨上站和蘭村站年徑流量序列的突變年份均為1966年，其中汾河水庫站年徑流量均值在1966年前后分別為4.00 萬m3和3.21 萬m3，突變后較突變前減少19.75%，人類活動對突變的貢獻率為62.46%；寨上站年徑流量均值在1966年前后分別為4.74 萬m3和3.57 萬m3，突變后較突變前減少24.68%，人類活動對突變的貢獻率為64.97%；蘭村站年徑流量均值在1966年前后分別為4.72 萬m3和3.25 萬m3，突變后較突變前減少31.14%，人類活動對突變的貢獻率為68.94%，干流徑流量序列突變的主要影響因素為人類活動。

（2）根據(jù)臨界慢化現(xiàn)象，徑流量序列在發(fā)生突變前徑流狀態(tài)變量的恢復速率變慢和恢復力變小，得到上靜游站的突變前兆信號發(fā)生于1973年，汾河水庫站、寨上站和蘭村站的突變前兆信號均發(fā)生于1958年，且徑流狀態(tài)變量的恢復速率及恢復力對人類活動干擾更敏感，更適用于求取汾河上游干流的突變前兆年份。

（3）根據(jù)臨界慢化現(xiàn)象，徑流量序列在發(fā)生突變前將導致序列的方差和自相關系數(shù)增大，得到上靜游站年徑流量序列的突變前兆信號發(fā)生于1968年，汾河水庫站、寨上站和蘭村站的年徑流量序列的突變前兆信號均發(fā)生于1960年，且徑流量序列的方差和自相關系數(shù)對降水量變化更敏感，更適用于求取汾河上游支流的突變前兆年份。