松花江干流依蘭-佳木斯段近40年徑流量時間序列分析

王春雷，邢貞相，付 強(qiáng)，閆丹丹，劉美鑫

(1.黑龍江省黑河水文局，黑龍江 黑河 164300；2.東北農(nóng)業(yè)大學(xué)水利與建筑學(xué)院，哈爾濱 150030)

徑流作為地貌形成的外營力之一，參與在地殼中的化學(xué)過程，影響土壤發(fā)育，植物生長以及湖泊、沼澤形成等。徑流量作為地區(qū)工農(nóng)業(yè)供水的重要來源，制約著地區(qū)社會經(jīng)濟(jì)的發(fā)展規(guī)模，在國民經(jīng)濟(jì)中具有十分重要的意義。

中國東北地區(qū)的松花江、烏蘇里江、黑龍江匯流而成的三江平原腹地，是國家重要的商品糧基地的重要組成部分，2014年實現(xiàn)糧食總產(chǎn)量實現(xiàn)“十一連增”，對保障我國糧食安全做出了巨大貢獻(xiàn)。近些年來，由于盲目擴(kuò)展水田面積，加之農(nóng)民節(jié)水意識淡薄，使得地下水位普遍下降，進(jìn)而誘發(fā)了一系列水資源短缺問題。

對流域年徑流量進(jìn)行時間序列分析[1-3]，對研究旱澇變化規(guī)律，水資源量評估及農(nóng)業(yè)水資源的可持利用都具有重要的意義。松花江干流依蘭-佳木斯段地處三江平原西部，流經(jīng)依蘭、曙光農(nóng)場、湯原農(nóng)場多個國有農(nóng)場，徑流量是水資源重要的自然補(bǔ)給來源和灌溉水源。為此，陸志華[4,5]曾對松花江干流中游段的年內(nèi)分配規(guī)律進(jìn)行了研究，結(jié)果發(fā)現(xiàn)該河段徑流年內(nèi)分配極不均勻。為進(jìn)一步獲取松花江干流年徑流總量的多年變化趨勢性和周期性規(guī)律，本文擬采用趨勢分析法、Mann-Kendall突變檢驗法和小波分析技術(shù)研究松花江依蘭-佳木斯段的年徑流量的變化特征。

1 研究區(qū)概況

松花江流域介于北緯41°42′～51°38′、東經(jīng)119°52′～132°31′，流域面積為55.68萬km2，約占黑龍江總流域面積的30.2%。流域西部以額爾古訥河、大興安嶺為界，海拔高度介于700～1 700 m之間；北部以黑龍江、小興安嶺為界，海拔高度介于1 000～2 000 m之間；東南部以張廣才嶺、完達(dá)山脈與烏蘇里江、圖們江等為界，海拔高度介于200～2 700 m之間；西南部以遼河、松花江的松遼分水嶺為界，海拔高度介于140～250 m之間。于同江附近匯入黑龍江的松花江，與黑龍江、烏蘇里江下游的廣袤土地共同構(gòu)成著名的三江平原。

本文主要以松花江干流依蘭站、佳木斯站為研究對象，探究其徑流量的時間變化特征，以期為今后工農(nóng)業(yè)生產(chǎn)規(guī)劃的制定提供依據(jù)(見圖1)。

圖1 研究區(qū)域概況Fig.1 The geographic location of the Sanjiang Plain

2 數(shù)據(jù)與方法

2.1 基本氣象數(shù)據(jù)

徑流量數(shù)據(jù)摘錄自《黑龍江省水文年鑒》。根據(jù)收集的已有年鑒整理結(jié)果，分析發(fā)現(xiàn)1972-2010年期間的依蘭、佳木斯兩站的徑流資料均包含豐水年13 a、平水年12 a、枯水年14 a，具有較好的代表性，因此選取該時段內(nèi)的徑流量觀測數(shù)據(jù)作為本次研究的基礎(chǔ)資料，通過進(jìn)一步分析和計算得到依蘭站和佳木斯站的年徑流量。

2.2 徑流量時間序列分析

本文從趨勢性、突變性和周期性三方面分別對徑流量時間序列分析。其中，趨勢性分析采用趨勢分析法，該方法能夠較好地獲取長序列總體變化趨勢特征，但無法分析序列內(nèi)部的具體突變特征，因此，本文在進(jìn)行徑流序列突變性分析時選用世界氣象組織(WMO)推薦的Mann-Kendall突變分析法[6-9]，該方法可獲取長序列內(nèi)部具體突變點(diǎn)和突變特征；以上兩種方法均無法確切地分析長序列周期性變化特征，故在分析徑流序列周期性時采用Morlet小波分析法[10-13]對其進(jìn)行綜合判斷，該方法能夠準(zhǔn)確地確定序列各個主周期，獲取序列的周期性變化特征。

2.2.1趨勢分析法

將樣本量為n的某氣候變量用xi表示，其所對應(yīng)的時間用ti表示，建立用以描述xi與ti關(guān)系的一元線性回歸方程：

xi=a+bi(i=1,2,…,n)

(1)

式中：a為回歸常數(shù)；b為回歸系數(shù)。a、b均用最小二乘法對其進(jìn)行估計。

對觀測數(shù)據(jù)及其所對應(yīng)的時間ti、回歸系數(shù)b、常數(shù)a進(jìn)行最小二乘估計，分別為：

(2)

(3)

(4)

其中b為趨勢分析系數(shù)。對于線性回歸的計算結(jié)果，主要分析回歸系數(shù)b的正負(fù)表示氣候變量的趨勢傾向，即b>0時x隨時間t的增加呈上升趨勢，b<0時x隨時間t的增加呈下降趨勢。b值的大小反映了上升或下降的傾向程度。

2.2.2突變分析法

設(shè)有一時間序列如下：x1,x2,…,xn構(gòu)造秩序列ri，ri表示當(dāng)xi>xj(1≤j≤i)時的樣本累積數(shù)。定義：

(6)

sk均值E(sk)與方差Var(sk)的定義如下：

(8)

假定時間序列隨機(jī)獨(dú)立，對統(tǒng)計量UFk作如下定義：

(9)

其中UF1=0。UFk為標(biāo)準(zhǔn)正態(tài)分布，k給定一顯著水平α，查正態(tài)分布表得到臨界值Uα，當(dāng)|UFk|>Uα，表明序列存在一個明顯的趨勢變化，UFk表示為UF。把此法應(yīng)用到反序列中，重復(fù)以上計算過程，將計算值乘以-1得到UBk，UBk用UB表示，令UB1=0。

通過統(tǒng)計序列UFk和UBk可對序列x的變化趨勢進(jìn)一步進(jìn)行分析，同時明確突變時間，指出突變區(qū)域。若UFk>0，則序列呈上升趨勢；若UFk<0，則序列呈下降趨勢；當(dāng)它們超過臨界直線時，則序列變化趨勢顯著。

2.2.3小波分析法

小波分析是將一簇頻率各異的振蕩函數(shù)作為窗口函數(shù)φ(t)對信號f(t)進(jìn)行掃描和平移。時間序列f(kt)(k=1,2,…,N；Δt為取樣的時間間隔)的小波變換：

(10)

式中：a為尺度因子；b為平移因子；Wf(a,b)稱為小波系數(shù)。

基本小波函數(shù)涵蓋墨西哥帽小波(Mexican hat)、Morlet小波和Wave小波等。本文采用Morlet小波，小波系數(shù)：

(11)

式中：c為常數(shù)；i為虛數(shù)。

它是周期函數(shù)經(jīng)由Gaussian函數(shù)平滑得到，其伸縮尺度a與Fourier分析中的周期T有如下關(guān)系：

(12)

因此當(dāng)取c=6.2時，周期T可以近似用a來替代。時間域上將與a相關(guān)的全部小波變換系數(shù)進(jìn)行平方然后積分，即為小波方差：

(13)

式中：Var(a)為小波方差；Wf(a,b)為小波系數(shù)。

小波方差反映的是波動能量隨尺度a的分布，可用以確定各尺度在一個時間序列中擾動的相對強(qiáng)度，與峰值處相對應(yīng)的尺度即為該序列的主周期。

3 結(jié)果分析

3.1 趨勢性分析

年際變化上，近40年來松花江依蘭站和佳木斯站年徑流量均呈下降趨勢(見圖2)，年際傾向率分別為-2.532 3、-5.142 4 億m3/a。依蘭站2010年徑流量為603.111億m3，同1972年徑流量相比，增加了7.72%；佳木斯站2010年徑流量為621.495億m3，同1972年徑流量相比，減少了10.40%。為消除周期變化對數(shù)據(jù)進(jìn)行五年滑動平均，結(jié)果表明近40年來徑流量無明顯的總體變化趨勢，呈波動狀態(tài)。

圖2 1972-2010年徑流量總體變化特征Fig.2 The overall variation characteristics of annual runoff from 1972 to 2010

3.2 突變性分析

為進(jìn)一步分析松花江依蘭-佳木斯段徑流量變化趨勢，現(xiàn)采用M-K方法深入分析徑流量的突變特征。M-K突變檢驗分析研究結(jié)果見圖3(圖3中虛線表示α=95%的顯著性水平的臨界值)。

根據(jù)圖3分析可知：依蘭站UF和UB兩條曲線可知，依蘭站徑流量在1973年之前為上升趨勢，1974-1983年為下降趨勢，1984-2004年為上升趨勢，2004年之后為下降趨勢；突變發(fā)生在1975、1976和2004年。佳木斯站UF和UB兩條曲線可知，佳木斯站徑流量在1973年之前為上升趨勢，1974-1984年為下降趨勢，1985-2002年為上升趨勢，2002年之后為下降趨勢；突變發(fā)生在1975、1976和2000年。

綜合兩站徑流的突變特征，兩站年徑流量突變點(diǎn)均發(fā)生在20世紀(jì)70年中期、80年代中期和21世紀(jì)初，這與20世紀(jì)70年大規(guī)模農(nóng)業(yè)開墾加劇流域蒸發(fā)而致使徑流量減少有關(guān)；而自20世紀(jì)80年代中期“家庭聯(lián)產(chǎn)承包責(zé)任制”實施以來，農(nóng)業(yè)種植積極性得到極大提高，中國農(nóng)業(yè)產(chǎn)業(yè)結(jié)構(gòu)調(diào)整初見成效，農(nóng)業(yè)生產(chǎn)力極大發(fā)展，工農(nóng)業(yè)水資源消耗量劇增，致使年徑流量恢復(fù)能力減弱；21世紀(jì)初，隨著黑龍江省糧食產(chǎn)能的連續(xù)提升，水資源用量尤其是水稻種植面積進(jìn)一步擴(kuò)大，加之水利工程建設(shè)日益增多，致使研究河段的年徑流復(fù)又下降，為此，黑龍江省正在大力推進(jìn)“節(jié)水增糧食”行動、積極推廣“水稻節(jié)水控制灌溉技術(shù)”和加快“兩大平原”現(xiàn)代農(nóng)業(yè)綜合配套改革實驗方案的實施，以提高水資源的利用效率。

圖3 年徑流量M-K突變分析結(jié)果Fig.3 The M-K mutation analysis on annual runoff

3.3 周期性分析

圖4為小波變換系數(shù)實部時頻分布圖，從圖4可以看出，近40 a來松花江依蘭站和佳木斯站包含了不同尺度的周期變化。圖5為松花江依蘭站、佳木斯站Morlet小波方差分布，可對圖4中所識別出的周期進(jìn)一步判斷。綜合分析圖4和圖5可知：18、8和5 a左右尺度波動相對較為明顯，存在徑流量偏多偏少的循環(huán)變化，且2010年后上述尺度的小波正在形成且為正值，因此預(yù)測在這3種尺度上未來徑流量均呈現(xiàn)偏多趨勢。

圖4 小波變換系數(shù)實部的時頻分布Fig.4 The time-frequency distribution of real parts of wavelet transform coefficient

圖5 小波變換方差Fig.5 The variance distribution of real parts of Morlet wavelet

本文得出的佳木斯站年徑流序列趨勢分析結(jié)果和突變分析結(jié)果與文獻(xiàn)[4]的研究成果具有較好的一致性。此外，與文獻(xiàn)[4]相比，本文延長了徑流序列長度，同時增加了依蘭站進(jìn)行對比分析，使本文結(jié)果可靠性更強(qiáng)。

4 結(jié) 語

本文運(yùn)用松花江干流依蘭站和佳木斯站1972-2010年逐年徑流量資料，采用趨勢分析法、Mann-Kendall突變分析法和Morlet小波分析法，對松花江依蘭-佳木斯段年徑流量進(jìn)行時間序列分析，得到出以下幾點(diǎn)結(jié)論。

(1)年際變化上，近40 a來松花江依蘭站和佳木斯站年徑流量表現(xiàn)出總體微弱下降趨勢，且伴隨明顯的波動特性。這與全球氣候變暖使蒸發(fā)加劇和水文過程的周期性有密切關(guān)系。

(2)時間突變上，兩站年徑流量均存在多個 突變點(diǎn)，且突變點(diǎn)具有較好的同步性，這與該地區(qū)的人類活動、社會生產(chǎn)力和上下游的水文聯(lián)系有極大關(guān)系。因此，欲控制年徑流的突變，應(yīng)從規(guī)范人類活動、推廣高效節(jié)水技術(shù)等方面入手，同時要注意考慮上下游的影響。

(3)周期變化上，松花江依蘭站和佳木斯站包含了不同尺度的周期變化，18、8和5 a左右尺度波動相對較為明顯，存在徑流量偏多偏少的循環(huán)變化，且2010年后上述尺度的小波正在形成且為正值，因此預(yù)測在這3種尺度上未來徑流量均呈現(xiàn)偏多趨勢。

□

[1] 王 英,曹明奎,陶 波,等. 全球氣候變化背景下中國徑流量空間格局的變化特征[J]. 地理研究,2006,25(6):1 031-1 040.

[2] 袁 喆,嚴(yán)登華,楊志勇,等. 1961-2010年中國400 mm和800 mm等雨量線時空變化[J]. 水科學(xué)進(jìn)展,2014,25(4):494-502.

[3] 徐利崗,周宏飛,梁 川,等. 中國北方荒漠區(qū)降水多時間尺度變異性研究[J]. 水利學(xué)報,2009,40(8):1 002-1 011.

[4] 陸志華，夏自強(qiáng)，于嵐嵐. 松花江佳木斯站徑流變化規(guī)律及演變趨勢分析[J].水電能源科學(xué)，2011,29(4)：14-17.

[5] 陸志華，夏自強(qiáng)，于嵐嵐，等.松花江干流中游段年內(nèi)分配變化規(guī)律[J].河海大學(xué)學(xué)報(自然科學(xué)版)，2012,40(1)：63-69.

[6] 胡 剛,宋 慧. 基于Mann-Kendall的濟(jì)南市氣溫變化趨勢及突變分析[J]. 濟(jì)南大學(xué)學(xué)報(自然科學(xué)版),2012,26(1):96-101.

[7] Z X Xu, K Takeuchi, H Ishidaira. Monotonic trend and step changes in Japanese precipitation[J]. Journal of Hydrology, 2003,279(3):144-150.

[8] Henry B Mann. Nonparametric tests against trend[J]. Econometrica, 1945, 13(3):245-259.

[9] M G Kendall. Rank correlation methods[M]. London: Charles Griffin, 1948: 202.

[10] 邵曉梅,許月卿,嚴(yán)昌榮. 黃河流域降水序列變化的小波分析[J]. 北京大學(xué)學(xué)報(自然科學(xué)版),2006,42(4):503-509.

[11] 祝青林,張留柱,于貴瑞,等. 近30年黃河流域降水量的時空演變特征[J]. 自然資源學(xué)報,2005,20(4):477-482.

[12] 劉 勤,嚴(yán)昌榮,張燕卿,等. 近50年黃河流域氣溫和徑流量變化特征分析[J]. 中國農(nóng)業(yè)氣象,2012,33(4):475-480.

[13] 李占杰,魚京善. 黃河流域降水要素的周期特征分析[J]. 北京師范大學(xué)學(xué)報(自然科學(xué)版)，2010，46(3):401-404.

[14] 邢貞相，閆丹丹，劉美鑫，等. 三江平原近60年降水量時空變異特征分析[J]. 農(nóng)業(yè)機(jī)械學(xué)報,2015,46(11):337-343.