陡河水庫(kù)年徑流系列變異點(diǎn)分析

和吉，燕琴，雷冠軍

（1.華北水利水電大學(xué)水利學(xué)院，河南鄭州450000；2.中國(guó)水利水電科學(xué)研究院水資源研究所，北京100038）

陡河水庫(kù)年徑流系列變異點(diǎn)分析

和吉1，燕琴1，雷冠軍2

（1.華北水利水電大學(xué)水利學(xué)院，河南鄭州450000；2.中國(guó)水利水電科學(xué)研究院水資源研究所，北京100038）

文中利用R/S分析法和滑動(dòng)平均系數(shù)法2種方法對(duì)陡河水庫(kù)流域1953—1999年47年的降雨、徑流序列進(jìn)行了變異點(diǎn)分析，可知陡河水庫(kù)年徑流的變異點(diǎn)為1975年和1985年，與陡河流域70年代大興水利工程、80年代降雨偏枯及地下水開采過(guò)量等人類活動(dòng)因素對(duì)其的影響相符，結(jié)果可靠，具有一定的實(shí)用價(jià)值。

R/S分析法；滑動(dòng)平均相關(guān)系數(shù)法；年徑流；陡河水庫(kù)

變異是指水文觀測(cè)序列在某個(gè)時(shí)刻突然變化，發(fā)生變化的時(shí)刻稱為變異點(diǎn)。傳統(tǒng)的水文序列變異點(diǎn)的研究方法有水量平衡法、對(duì)比分析法、流域水文模擬法、熵值法、相關(guān)系數(shù)法等[1]。這些方法大都以單純的數(shù)理統(tǒng)計(jì)變量來(lái)對(duì)系統(tǒng)進(jìn)行分析檢驗(yàn)，但序列值的基準(zhǔn)點(diǎn)的選擇受人為因素影響，使得變異點(diǎn)很難一致[2]。而R/S分析法與滑動(dòng)相關(guān)平均系數(shù)法的共同點(diǎn)是不直接對(duì)水文序列的均值、方差等系統(tǒng)狀態(tài)變量分析檢驗(yàn)，且R/S分析法比較強(qiáng)健，假定條件少，能將隨機(jī)序列與非隨機(jī)序列區(qū)分開；滑動(dòng)相關(guān)系數(shù)法研究?jī)蓚€(gè)相關(guān)序列之間相關(guān)的性變異，方法簡(jiǎn)單[3]。

1 R/S分析法

R/S分析法，也稱重標(biāo)極差分析法，是20世紀(jì)英國(guó)水文學(xué)家Hurst提出的，經(jīng)過(guò)Mandelbrot（1972，1975）、Mandelbrot、Wallis（1969）和Lo（1991）等人的研究逐步完善。R/S分析法認(rèn)為人類水事活動(dòng)能使流域水文過(guò)程序列的分形特征值、赫斯特指數(shù)、分形維數(shù)發(fā)生明顯變化，可通過(guò)赫斯特指數(shù)判斷水文過(guò)程的變異點(diǎn)[3]。

1）R/S分析法處理序列值

一時(shí)間序列：ξ1，ξ2，ξ3，…，ξn，對(duì)應(yīng)的時(shí)間為t1，t2，t3，…，tn。將此連續(xù)時(shí)間序列抽出一個(gè)時(shí)間尺度的序列ξ1，ξ2，ξ3，…，ξn，序列的時(shí)間跨度為τ=tN-t1，在時(shí)間τ內(nèi)，該序列的平均值為：

式中：N為時(shí)間序列的長(zhǎng)度。

在tj時(shí)刻，序列值ξi相對(duì)于其平均值＜ξ＞Ni的累積偏差為：

同一個(gè)N值所對(duì)應(yīng)的最大X（t）值和最小X（t）值之差為極差，并記為：

Hurst利用的標(biāo)準(zhǔn)偏差為∶

2）Hurst指數(shù)

以ln（R/S）為縱坐標(biāo)，lnτ為橫坐標(biāo)，作線性回歸圖，圖形是一簇斜率為H、過(guò)定點(diǎn)（lnc，0）的直線。通過(guò)ln（R（τ）/S（τ））=Hlnc+Hlnτ計(jì)算H。

分別以ξi（i=10，11，12，…，n-10）為分界點(diǎn)將數(shù)據(jù)分為兩部分，對(duì)這兩部分進(jìn)行R/S分析，相應(yīng)的H記為H1和H2，求出H1和H2之差的絕對(duì)值△H=|H1-H2|，△H最大的一點(diǎn)為變異最大的點(diǎn)，即可相應(yīng)找出變異年份[3]。

2 滑動(dòng)平均相關(guān)系數(shù)法

相關(guān)分析法反映2個(gè)系列間統(tǒng)計(jì)相依程度，以一個(gè)步長(zhǎng)在2個(gè)相關(guān)水文序列上滑動(dòng)得到相關(guān)序列，新序列若有變異點(diǎn)，則有一個(gè)或者兩個(gè)能發(fā)生變異，若沒有變異點(diǎn)，不能直接確定2個(gè)序列是否存在變異?；瑒?dòng)相關(guān)系數(shù)法是一種綜合診斷方法，求多組相關(guān)聯(lián)的序列的滑動(dòng)相關(guān)系數(shù)序列，診斷出新序列的變異點(diǎn)，變異點(diǎn)反映每組序列相關(guān)程度的變化情況，從而得到2個(gè)序列的可能變異點(diǎn)，綜合分析可得徑流序列的變異點(diǎn)[4]。

先求出2個(gè)序列（P1，P2，P3，…，Pm）與（W1，W2，W3，…，Wm）的相關(guān)系數(shù)r1，然后再求（P2，P3，P4，…，Pm+1）與（W2，W3，W4，…，Wm+1）的相關(guān)系數(shù)r2，…，rn，得到一個(gè)相關(guān)系數(shù)序列，其中m≥10。

對(duì)m取多個(gè)值（取m=a），可以更有效地尋找變異點(diǎn)，并最終得到a個(gè)相關(guān)系數(shù)序列，并以時(shí)間為橫坐標(biāo)，滑動(dòng)相關(guān)系數(shù)為縱坐標(biāo)繪制曲線，得到a條滑動(dòng)相關(guān)系數(shù)曲線。將a個(gè)相關(guān)系數(shù)序列求平均值得到另一組滑動(dòng)相關(guān)系數(shù)序列，繪制滑動(dòng)相關(guān)系數(shù)系列平均值與時(shí)間的曲線。分析觀察這a條曲線以及滑動(dòng)相關(guān)系列平均值與時(shí)間的曲線，便可得到變異點(diǎn)。相關(guān)系數(shù)采用皮爾遜相關(guān)系數(shù)，其求解方程式為：

相關(guān)系數(shù)r的取值在-1.0～+1.0之間。r＞0時(shí)，兩變量呈正相關(guān)，越接近于1.0，正相關(guān)越顯著；r＜0時(shí)，兩變量呈負(fù)相關(guān)，越接近于-1.0，負(fù)相關(guān)越顯著；r=0時(shí)，則兩變量相互獨(dú)立[4]。

3 實(shí)例

3.1 概況

陡河位于燕山南部，介于沙河、還鄉(xiāng)河之間，河道地勢(shì)落差大，形式陡峭，獨(dú)流入海，全長(zhǎng)120 km，流域面積1 340 km2。陡河水庫(kù)位于唐山市東北15 km處，控制流域面積530 km2，總庫(kù)容為5.152×108m3，是以防洪為主并為城市生活和工農(nóng)業(yè)供水的大型水庫(kù)，也是引灤入唐工程的終端水利樞紐工程。流域內(nèi)地形比較復(fù)雜，坡陡流急，陡河流域處于半干旱溫帶季風(fēng)氣候區(qū)內(nèi)，季節(jié)變化明顯，受季風(fēng)氣候影響，多年平均降水量671 mm，降水年際變化大，年內(nèi)分配不均勻，年降水量的80%集中在6至9月份[6]。

3.2 基本資料

本文采用了陡河水庫(kù)流域以上1953—1999年共47年的降雨和徑流資料，其中降雨序列是通過(guò)陡河流域17處雨量站點(diǎn)實(shí)測(cè)資料，結(jié)合使用泰森多邊形法計(jì)算得到的，徑流序列是根據(jù)陡河水庫(kù)管理中心提供的《陡河水庫(kù)以上流域年徑流分析》摘取的。徑流量用水量平衡法還原，還原公式為：

式中：W還為還原年徑流量；W出為陡河水庫(kù)實(shí)測(cè)出庫(kù)流量；W蓄為陡河水庫(kù)年蓄變量；W滲為陡河水庫(kù)年滲漏水量；W用為工農(nóng)業(yè)及生活用水量；W蒸為陡河水庫(kù)蒸發(fā)損失水量；W引為跨流域引水量[7]。見表1。

表1 陡河水庫(kù)1953—1999年還原總水量

3.3 R/S法分析降雨和徑流序列變異點(diǎn)

使用R/S法分別對(duì)陡河水庫(kù)流域以上1953—1999年共47年的資料序列進(jìn)行分析，得到的△H序列與時(shí)間變化的折線圖，見圖1，2。

圖1 陡河水庫(kù)降雨序列分段R/S分析圖

圖2 陡河水庫(kù)徑流序列分段R/S分析圖

圖3 陡河水庫(kù)降雨—徑流序6種步長(zhǎng)滑動(dòng)相關(guān)系數(shù)序列值

圖4 陡河水庫(kù)降雨—徑流序6種步長(zhǎng)滑動(dòng)相關(guān)系數(shù)序列平均值

1）從R/S法分析降雨的變異點(diǎn)來(lái)看，降雨的最大△H在1975年，次大△H分別在1967年和1980年；1963—1975年的△H的數(shù)值普遍偏大，且變化不大；1975—1990年，△H最大值在1980年，△H值普遍偏小，降雨系列比較明顯的變異點(diǎn)在1975年。

2）從R/S法分析徑流的變異點(diǎn)來(lái)看，徑流量的最大△H在1975年，次大△H分別在1967年和1966年；1963—1975年的△H的數(shù)值普遍偏大，△H的值都在0.2以上；1977—1990年△H普遍偏小，除1980年和1985年之外其余年份的△H都在0.2以下，其中1985年是1978—1990年的△H最大點(diǎn)，綜合分析徑流量的變異點(diǎn)在1975年以及1985年。

3.4 滑動(dòng)相關(guān)系數(shù)法分析降雨和徑流序列變異點(diǎn)

使用滑動(dòng)相關(guān)系數(shù)法對(duì)陡河水庫(kù)流域以上1953—1999年共47年的降雨—徑流進(jìn)行分析，以n=12，14，16，18，20，22為分段，得到6個(gè)降雨—徑流皮爾遜相關(guān)系數(shù)序列，以及這個(gè)6個(gè)相關(guān)系數(shù)的平均值，并以相關(guān)系數(shù)為縱坐標(biāo)曲線，以時(shí)間為橫坐標(biāo)做相關(guān)系數(shù)曲線圖以及相關(guān)系數(shù)平均值曲線圖，見圖3，4。

觀察陡河水庫(kù)1953—1999年的降雨—徑流6種步長(zhǎng)（n=12，14，16，18，20，22）滑動(dòng)相關(guān)系數(shù)序列值與時(shí)間的曲線圖，同時(shí)觀察降雨—徑流各步長(zhǎng)滑動(dòng)相關(guān)系數(shù)序列的平均值與時(shí)間的曲線圖，可以看出降雨—徑流的皮爾遜相關(guān)系數(shù)按時(shí)間分布被1975年和1985年分為3個(gè)取值范圍，1960—1975年，皮爾遜相關(guān)系數(shù)的值在0.925上下浮動(dòng)；1975—1985年，皮爾遜相關(guān)系數(shù)的值在0.875上下浮動(dòng)；1985—1999年，皮爾遜相關(guān)系數(shù)的值在0.82上下浮動(dòng)。綜合分析，陡河水庫(kù)流域以上降雨—徑流系列的變異點(diǎn)發(fā)生在1975年以及1985年。

3.5 結(jié)果分析

利用R/S法與滑動(dòng)平均系數(shù)法分析陡河水庫(kù)流域降雨系列、徑流系列、降雨—徑流相關(guān)系列，分析得出陡河水庫(kù)降雨系列的變異點(diǎn)是1975年；徑流系列的變異點(diǎn)是在1975年和1985年；降雨—徑流皮爾遜相關(guān)系數(shù)的變異點(diǎn)在1975年和1985年。

結(jié)合實(shí)際情況分析，為滿足工農(nóng)業(yè)生產(chǎn)不斷發(fā)展的需要，70年代初開始，水利工程建設(shè)掀起高潮，人們?cè)诙负恿饔蚪⒘嗽S多中、小型水庫(kù)，灌渠，水閘，揚(yáng)水站，提引水工程等，變異點(diǎn)發(fā)生在1975與實(shí)際情況吻合。

根據(jù)陡河水庫(kù)管理中心提供的資料，1980—1984年由于降水偏少，同時(shí)地下水開采過(guò)量，導(dǎo)致了1984年較1977年地下水位下降了10 m左右，后幾十年由于降水偏豐，地下水位有所回升。所以徑流系列，以及降雨—徑流相關(guān)系數(shù)的變異點(diǎn)發(fā)生在1985年也與實(shí)際情況吻合。

［1］成靜清.非一致性年徑流序列頻率分析計(jì)算［D］.陜西：西北農(nóng)林科技大學(xué)，2010：1—117.

［2］沈宏.天然徑流還原計(jì)算方法初步探討［J］.水利規(guī)劃與設(shè)計(jì).2003（3）：15—18，47.

［3］周蓓，劉俊民，等.RS法在徑流還原和預(yù)測(cè)中的應(yīng)用［J］.人民長(zhǎng)江，2008，39（15）：42—45.

［4］陳廣圣，王義民，等.基于滑動(dòng)相關(guān)系數(shù)法的流域徑流變異診斷［J］.自然災(zāi)害學(xué)報(bào)，2016，25（1）：11—17.

［5］劉志濤，程兵峰.近10年陡河水庫(kù)來(lái)水變化趨勢(shì)及其影響因素分析［J］.海河水利，2011（4）：36—38.

［6］劉素英.陡河流域雨、水情分析［J］.河北水利科技，1999，20（1）：12—16.

［7］魏茹生.徑流還原計(jì)算技術(shù)方法及其應(yīng)用研究［D］.陜西：西安理工大學(xué)，2008：1—103.

TV121

A

1002－0624（2017）04－0018－03

2016-11-09