牧馬河流域年徑流變化特征分析

夏煥清，霍秀秀，曹苗苗

（陜西省漢中水文水資源勘測局，陜西 漢中 732000）

1 研究區(qū)概況與數(shù)據(jù)來源

1.1 研究區(qū)概況

牧馬河為漢江干流一級支流，源于城固、南鄭、西鄉(xiāng)3縣交界處米倉山北麓的白熊山老鷹崖，流域位于大巴山北部，山地高在2500 m～3000 m之間，大致東向流經(jīng)城固縣、西鄉(xiāng)縣，在西鄉(xiāng)城南三花石鄉(xiāng)回龍灣匯入漢江。西鄉(xiāng)水文站是牧馬河重要控制站，集水面積1224 km2，多年平均降雨量為877.8 mm，7月～9月降水量占全年的50%以上，多年平均徑流量為9.04億m3。

1.2 數(shù)據(jù)來源及研究方法

采用的徑流量數(shù)據(jù)為西鄉(xiāng)水文站(1974年～2018年)的年徑流量數(shù)據(jù)，見表1。

表1 西鄉(xiāng)水文站(1974年～2018年)的年徑流量數(shù)據(jù) 單位：億m3

研究方法為：采用線性回歸、五年滑動(dòng)平均對牧馬河流域徑流變化趨勢進(jìn)行分析，M-K檢驗(yàn)、小波分析分別進(jìn)行年徑流量突變檢驗(yàn)與周期分析。數(shù)據(jù)處理工具：小波分析和M-K檢驗(yàn)用Matlab，顯著性檢驗(yàn)用皮爾遜相關(guān)性分析。

1.3 研究方法

1.3.1 線性回歸

通過建立牧馬河流域西鄉(xiāng)站年徑流與時(shí)間序列之間的線性回歸方程，通過回歸方程的斜率以及變化趨勢的顯著程度從而判斷因變量趨勢變化的方法，表達(dá)式：

式中：Xt為t時(shí)刻年徑流值；t為時(shí)間；b為方程截距；a為回歸方程斜率。

用10 a表示研究對象時(shí)間序列趨勢變化程度，當(dāng)a＞0時(shí)，表明年徑流具有上升趨勢；當(dāng)a＜0時(shí)，具有下降趨勢。此外，通過皮爾遜相關(guān)分析法來判斷年徑流變化趨勢是否顯著，當(dāng)p＜0.05時(shí)，年徑流在時(shí)間序列上具有顯著性上升或下降趨勢；當(dāng)p＞0.05時(shí)，年徑流不具有顯著變化趨勢，本文通過SPSS17.0對年徑流變化趨勢進(jìn)行顯著性檢驗(yàn)。

1.3.2 滑動(dòng)平均

滑動(dòng)平均是可以直觀反映年徑流時(shí)間序列變化趨勢的一種方法，通過計(jì)算西鄉(xiāng)站年徑流系列（1974年～2018年）44年徑流項(xiàng)連續(xù)5項(xiàng)序列的平均值。從第1項(xiàng)至5項(xiàng)的算術(shù)平均為滑動(dòng)平均第1項(xiàng)，最后1項(xiàng)是原來44年徑流序列的第40至第44項(xiàng)的算術(shù)平均。若西鄉(xiāng)站年徑流序列具有趨勢成分，則5年滑動(dòng)平均能夠?qū)⑽鬣l(xiāng)站年徑流的趨勢直觀的顯示出來。

1.3.3 Mann-Kendall突變檢驗(yàn)

對于西鄉(xiāng)站年徑流樣本時(shí)間序列 X1、X2、X3、X4…Xn，mi表示第i個(gè)樣本Xi大于Xj（1≤j≤i）的累計(jì)數(shù)，則構(gòu)造統(tǒng)計(jì)量：

在序列相互獨(dú)立的假設(shè)條件下，的均值和方差為：

將公式標(biāo)準(zhǔn)化得到以下計(jì)算式：

UFk計(jì)算結(jié)果組成UF曲線，同樣以此方法引用到反序列中，計(jì)算得到 UBk，使得 UBk=-UFk。

給定顯著性水平0.05，統(tǒng)計(jì)量臨界值±1.96。通過繪制UFk、曲線圖UBk，如果兩條曲線在置信度0.05水平內(nèi)相交，則定義為突變點(diǎn)。如果UFk或UBk的值大于0，則表明序列具有上升趨勢，反之，序列呈下降趨勢。當(dāng)曲線超出置信度區(qū)間時(shí)，表明序列具有顯著的上升或下降趨勢，該過程通過Matlab工具調(diào)用編輯的西鄉(xiāng)站年徑流M-K突變檢驗(yàn)程序，來實(shí)現(xiàn)相關(guān)過程計(jì)算及作圖。

1.3.4 小波分析

小波分析方法是一種在空間固定的基礎(chǔ)上反映時(shí)頻局部發(fā)生變化的分析方法。基本原理用一簇小波函數(shù)系來表示或逼近某一信號或函數(shù)。小波分析最關(guān)鍵的部分是小波函數(shù)，小波函數(shù)φ（t）∈L2（R），且有：

φ（t）為基小波函數(shù)，它可通過尺度的伸縮和時(shí)間軸上的平移構(gòu)成一簇函數(shù)系：

式中：φa,b（t）為子小波；a為尺度因子，反映小波的周期長度；b為平移因子，反應(yīng)時(shí)間上的平移。

對西鄉(xiāng)站年徑流的研究選用的是Morlet小波函數(shù)，其表達(dá)式為：

式中：ω0為常數(shù)（ω0≥5）；i為虛數(shù)。

對于任意函數(shù)f（X），其小波變換可以定義為：

牧馬河流域西鄉(xiāng)站年徑流周期分析計(jì)算與出圖過程，通過Matlab工具調(diào)用編輯的西鄉(xiāng)站年徑流周期分析程序來實(shí)現(xiàn)小波分析。

3 結(jié)果與分析

3.1 趨勢分析

圖1中圓點(diǎn)折現(xiàn)圖為年徑流量圖，三角折線圖為五年滑動(dòng)平均圖，方程式為年徑流量圖趨勢線方程式。由圖可知，年徑流趨勢線k為-0.1505，即年徑流量趨勢為減?。焕闷栠d相關(guān)性分析，在顯著性水平為0.01以下，皮爾遜相關(guān)性系數(shù)為-0.459，兩個(gè)變量為負(fù)相關(guān)且顯著相關(guān)；五年平均滑動(dòng)圖趨勢線k為-0.1346，五年平均徑流量趨勢為減小。

圖1 年徑流及五年滑動(dòng)平均圖

3.2 突變分析

圖2 UF小于0表明呈下降趨勢，當(dāng)UF和UB兩條曲線出現(xiàn)交點(diǎn)，且交點(diǎn)在臨界直線之間，交點(diǎn)對應(yīng)的時(shí)間即為突變開始時(shí)刻。從圖2可知，西鄉(xiāng)站年徑流量突變時(shí)間為1990年。

圖2 M-K突變檢驗(yàn)圖

3.3 周期風(fēng)險(xiǎn)

在不同的時(shí)間尺度下，年徑流量的變化有周期性和時(shí)間域上豐、枯變化的特性。

圖3 牧馬河流域年徑流小波系數(shù)模方圖

圖3 為小波模方圖，由圖3來分析年徑流在小波變化域中波動(dòng)能量大小的變化特性，進(jìn)而反映哪些能量聚集中心主導(dǎo)年徑流在時(shí)間域上的波動(dòng)變化的情形，分析可知，在5 a、10 a、23 a左右的時(shí)間尺度上能量聚集中心閉合較好并且分布密集，與圖5比較沒有30 a以上的，原因?yàn)闀r(shí)間序列尺度不夠長。

圖4 牧馬河流域年徑流小波變換實(shí)部圖

圖4 為小波系數(shù)實(shí)部的時(shí)頻分布圖，表示年徑流在小波變化域中以基準(zhǔn)面上下起伏的波動(dòng)情況，從而反映年徑流在時(shí)間尺度上豐、枯變化的特性，依圖可知西鄉(xiāng)站年徑流豐枯變換情況，豐水年在1980年、1989年、2011年左右，枯水年在1986年、2001年左右。

圖5 牧馬河流域年徑流小波變化方差圖

由圖5可知，小波方差在8 a、22 a、38 a左右的時(shí)間尺度上存在峰值，22 a左右的周期震蕩最強(qiáng)，是流域的第一主周期，第二、三周期依次約為8 a、38 a，由于西鄉(xiāng)水文站建站時(shí)間較短，年徑流數(shù)據(jù)序列較短，因此無法準(zhǔn)確預(yù)測更長時(shí)間尺度的周期。

總體來看，西鄉(xiāng)水文站年徑流周期為22 a～23 a，三個(gè)豐水年之間的間隔為 9 a、22 a、31 a 與圖 5 分析的周期 8 a、22 a、38 a左右，契合度較高。但在長周期上有較大差別，分析主要原因是隨著時(shí)間序列的增加，小時(shí)間尺度上的突變點(diǎn)逐漸退化為非突變點(diǎn)，隨著時(shí)間尺度的降低，突變點(diǎn)增多，大時(shí)間尺度上的一些平常態(tài)，在低層次上轉(zhuǎn)化為突變態(tài)，而本文因時(shí)間序列較短，在長周期上判斷不夠準(zhǔn)確，進(jìn)一步說明年徑流周期隨時(shí)間尺度的不同而變化。

4 結(jié)論

采用采用線性回歸、滑動(dòng)平均法、MK檢驗(yàn)、皮爾遜相關(guān)分析和小波分析方法分析牧馬河西鄉(xiāng)水文站天然年徑流時(shí)間序列進(jìn)行多時(shí)間尺度變化特性的分析。結(jié)果顯示:西鄉(xiāng)水文站的年徑流呈減小趨勢，年徑流突變時(shí)間為1990年左右，呈周期性變化，主周期為22 a～23 a左右，它主導(dǎo)著牧馬河年徑流長期變化的特性，同時(shí)也存在著8 a～10 a、38 a左右的較強(qiáng)周期變化。在1980年、1989年、2011年左右為豐水年，在1986年、2001年左右為枯水年。