基于小波和MK檢驗(yàn)的董志塬年降水量分析

趙雙慶，范 文，于寧宇

(長(zhǎng)安大學(xué) 地質(zhì)工程與測(cè)繪學(xué)院，陜西 西安 710054)

董志塬是中國(guó)黃土高原區(qū)面積最大、保存最完整的黃土塬，但干旱和水土流失限制了該地區(qū)社會(huì)經(jīng)濟(jì)的可持續(xù)發(fā)展。董志塬地區(qū)以雨養(yǎng)農(nóng)業(yè)，農(nóng)業(yè)發(fā)展與降雨密切相關(guān)，降雨是當(dāng)?shù)剞r(nóng)業(yè)生產(chǎn)的主要限制因素。引起水土流失最主要的外部動(dòng)力因素是降水，由其引發(fā)的水土流失問題是黃土高原地區(qū)最嚴(yán)峻的環(huán)境問題，大面積的水土流失導(dǎo)致了土地資源的喪失和退化[1-4]。因此，研究董志塬區(qū)的降水時(shí)間分布及變化趨勢(shì)，對(duì)于定量估算董志塬區(qū)區(qū)域土壤侵蝕量、預(yù)測(cè)未來降水變化趨勢(shì)、合理安排農(nóng)業(yè)生產(chǎn)以及科學(xué)制定水土保持措施具有重要的指導(dǎo)意義。降水在時(shí)間和空間上的分布規(guī)律及其對(duì)地表生態(tài)系統(tǒng)的影響是水文生態(tài)學(xué)研究的關(guān)鍵問題之一[5]。一直以來，降水規(guī)律都是國(guó)內(nèi)外專家學(xué)者廣泛關(guān)注的研究課題。Silva[6]使用MK方法對(duì)巴西東北部的氣候變量的時(shí)間序列進(jìn)行分析, 指出該地區(qū)的相對(duì)濕度和降水量呈下降趨勢(shì)。Andreo等人[7]通過頻譜分析、相關(guān)分析和連續(xù)小波分析研究了二十世紀(jì)以來伊比利亞半島南部的氣候和水文變化，指出該地區(qū)年平均降水存在5年的恒定周期，并預(yù)測(cè)整個(gè)21世紀(jì)該地區(qū)降水量呈下降趨勢(shì)。左洪超等人[8]對(duì)中國(guó)近50年來降水量進(jìn)行分析，表明我國(guó)年平均降水量沒有明顯的變化趨勢(shì)但年降水量的變化存在著明顯的區(qū)域特征。

不同地區(qū)降水變化的時(shí)間特征不同，對(duì)不同流域的降水時(shí)間特征進(jìn)行評(píng)價(jià)仍然有很重要的價(jià)值，有很多專家學(xué)者應(yīng)用MK檢驗(yàn)和小波分析方法，研究了黃土高原地區(qū)降水量的趨勢(shì)、周期、突變性等降水時(shí)間特征, 對(duì)黃土高原降水變化規(guī)律有了更深層次的認(rèn)識(shí)[9-15]。目前，對(duì)黃土地區(qū)的降水規(guī)律研究大多是大范圍整體性研究，對(duì)董志塬地區(qū)的降水規(guī)律研究較少，所以本文對(duì)該地區(qū)的年降水規(guī)律進(jìn)行深入分析，將對(duì)董志塬的水土保持、生態(tài)建設(shè)以及災(zāi)害預(yù)防起到一定的指導(dǎo)意義。

1 研究區(qū)概況

董志塬位于甘肅省慶陽市中南部，北起驛馬鎮(zhèn)北部的桐川鄉(xiāng)沿蔡家廟溝至慶城縣一線，南止涇河谷底，東西介于馬蓮河與蒲河之間，地理坐標(biāo)為N35.25—36.06,E107.46—107.96[16-18]。董志塬地區(qū)屬于大陸性半干旱氣候，降水分布不均勻，時(shí)有伏旱、暴雨、冰雹災(zāi)害；秋季降溫迅速，潮濕多雨；冬季寒冷漫長(zhǎng)，降雨稀少，年降雨量500～600 mm。涇河、馬蓮河、蒲河、黑河4條河流的373條一級(jí)支溝從四面八方伸向塬面腹地，塬區(qū)周邊受沖溝切蝕呈樹枝狀，形成塬高谷深，溝壑縱橫的隴東黃土塬侵蝕地貌。慶城縣附近溝壑區(qū)侵蝕地貌見圖1。董志塬區(qū)塬面和溝谷地層劃分為：第四系全新世Q4(沖積土層)、晚更新世Q3(馬蘭黃土)、中更新世Q2(離石黃土)、早更新世Q1(午城黃土)、上、中新世N2(藍(lán)田組紅粘土)，工程地質(zhì)圖見圖2。

圖1 董志塬上的溝壑地形Fig.1 The gulliesofthe Dongzhi Tableland

圖2 董志塬北部研究區(qū)工程地質(zhì)圖Fig.2 The engineering geological map of the research area in northern of the Dongzhi Tableland

2 數(shù)據(jù)來源與研究方法

2.1 數(shù)據(jù)來源

本文降水資料來源于中國(guó)氣象局董志塬區(qū)內(nèi)的西峰站點(diǎn)及南部附近的長(zhǎng)武站點(diǎn)1959—2016年的逐日監(jiān)測(cè)數(shù)據(jù)。氣象監(jiān)測(cè)站點(diǎn)位置，見圖3。為了得到董志塬的降雨變化的總體趨勢(shì)，將各站日降雨量累計(jì)，得到1959—2016年各年的年降雨數(shù)據(jù)。

圖3 董志塬數(shù)字高程模型及氣象站的位置Fig.3 The numerical elevation model and location of weather station in the Dongzhi Tableland

為了保證數(shù)據(jù)質(zhì)量，本研究?jī)H使用中國(guó)氣候參考網(wǎng)(CCRN)的數(shù)據(jù)。由中國(guó)氣象局管理的CCRN由最先進(jìn)的站點(diǎn)組成，這些站點(diǎn)配備了高質(zhì)量的儀器，并遵循嚴(yán)格的觀測(cè)協(xié)議，特別是用于檢測(cè)氣候變化信號(hào)[19]。

2.2 Mann-Kendall檢驗(yàn)法

本文采用MK檢驗(yàn)方法來分析董志塬降水站點(diǎn)的年降水量,得到董志塬地區(qū)年降水量變化趨勢(shì)以及降水量突變年份，具體分析方法簡(jiǎn)述如下。

2.2.1 MK趨勢(shì)檢驗(yàn)

(1)構(gòu)建董志塬降水站點(diǎn)的年降水量時(shí)間序列(x1，x2，…，xn)，假設(shè)H0表示該時(shí)間序列中的數(shù)據(jù)是獨(dú)立的，隨機(jī)變量同分布的樣本，沒有線性變化趨勢(shì)。

(2)計(jì)算檢驗(yàn)的統(tǒng)計(jì)變量S:

(1)

(3)對(duì)于有10個(gè)或者更多數(shù)據(jù)點(diǎn)的時(shí)間序列，使用正態(tài)近似。統(tǒng)計(jì)值S的方差計(jì)算如下：

(2)

其中年降水量相同的樣本的個(gè)數(shù)記為i，i值相同的樣本中降水量不同的樣本的個(gè)數(shù)記為ei。

(4)構(gòu)造的標(biāo)準(zhǔn)正態(tài)分布變量Z計(jì)算如下式：

(3)

(5)在一定的α置信水平上,如果Z≥|Z1-α/2|則拒絕原假設(shè),即在α置信水平上, 年降水量時(shí)間序列存在明顯的趨勢(shì)。

(6)其變化趨勢(shì)的大小用β表示,計(jì)算如下：

(4)

若β>0，表示呈上升趨勢(shì)；若β<0，表示呈下降趨勢(shì)。

2.2.2 MK突變檢驗(yàn)

根據(jù)董志塬降水站點(diǎn)的年降水量時(shí)間序列(x1，x2，…，xn)，然后我們構(gòu)造一秩序列：

(5)

其中

(6)

定義統(tǒng)計(jì)量：

(7)

其中，sk的均值和方差分別為E(sk)，var(sk)，時(shí)間序列獨(dú)立同分布的情況下，它們的表達(dá)式如下：

(8)

把時(shí)間序列x逆序排列xn，xn-1，…，x1，再通過(7)式求出UFk，同時(shí)令UBk=-UFk(k=n，n-1，…，1)，UB1=0。得到的統(tǒng)計(jì)序列UFk和UBk，UFk的正負(fù)表明xk的變化趨勢(shì)，正值表明序列呈上升趨勢(shì)，負(fù)值表明呈下降趨勢(shì)；當(dāng)UFk和UBk超過顯著性臨界直線時(shí)，說明趨勢(shì)顯著。UFk和UBk曲線的交點(diǎn)如果在臨界線之間，則交點(diǎn)對(duì)應(yīng)的時(shí)刻為突變發(fā)生的時(shí)間[20]。

2.3 Morlet小波分析

在分析時(shí)間序列時(shí)，需要獲得平滑連續(xù)的小波振幅以及時(shí)間序列的振幅和相位的信息，因而選擇非正交的復(fù)值小波。Morlet小波是非正交性的復(fù)值小波，因此國(guó)內(nèi)外專家學(xué)者大多采用Morlet小波對(duì)氣象時(shí)間序列進(jìn)行周期分析[21]。曹曉敏等[9]應(yīng)用小波分析，得出黃土高原中部降雨量具有2～3 a和7 a左右的年際振蕩，其中以2～3 a的年際振蕩尤為明顯。郝立生等[21]通過Morlet小波分析發(fā)現(xiàn)，太陽活動(dòng)存在14.5 a的顯著周期。李淼等[23]采用Morlet小波函數(shù)，對(duì)該地區(qū)近300 a 來降水的年際變化時(shí)間序列進(jìn)行了小波分析，得出該地區(qū)年降水量時(shí)間序列存在21 a、35 a 和85 a 左右的主周期，并指出85 a 周期是第一主周期。本研究采用氣象學(xué)分析時(shí)間序列常用的Morlet小波函數(shù)，對(duì)董志塬降水站點(diǎn)的降水時(shí)序進(jìn)行分析，并用紅噪聲或白噪聲標(biāo)準(zhǔn)譜對(duì)小波功率譜的顯著性進(jìn)行檢驗(yàn)，從而得到董志塬年降水量周期性變化規(guī)律。

Morlet小波函數(shù)及離散小波變換公式為：

φ(t)=π-1/4e-t2/2eiabt

(9)

(10)

式中，ω0是角頻率，Torrence等人認(rèn)為取6比較合適[21-22],ωf(a,b)為小波變換后的小波系數(shù)，f(iδt)為進(jìn)行分析的降水時(shí)間序列，δt為時(shí)間序列的時(shí)間間隔；ωf(a,b)為卷積后的小波系數(shù)，a和b分別為尺度因子(表示小波的周期長(zhǎng)度)和平移因子(表示時(shí)間上的平移)，“*”表示復(fù)共軛。本文對(duì)年降水量進(jìn)行分析，降水時(shí)序的時(shí)間間隔為1年，δt取1。

小波功率譜Ea,b定義為：

(11)

小波功率譜的顯著性檢驗(yàn)一般應(yīng)用紅噪聲或白噪聲標(biāo)準(zhǔn)譜來實(shí)現(xiàn)。求得降水時(shí)間序列滯后1的自相關(guān)系數(shù)r(1)，如果大于0.1則用紅噪聲譜進(jìn)行檢驗(yàn)，若小于0.1，則令r(1)=0，用白噪聲譜進(jìn)行檢測(cè)。Torrence[21]指出小波功率譜服從χ2分布。在一定的顯著水平上，計(jì)算得到紅噪聲或白噪聲的顯著性理論譜P。如果理論譜小于小波功率譜，則說明這種尺度對(duì)應(yīng)的周期具有顯著性。

理論功率譜：

(12)

Pa為紅噪聲或白噪聲譜：

(13)

若Ea,b>P，說明小波功率譜是顯著的。

3 結(jié)果與分析

3.1 年降水量變化趨勢(shì)分析

通過前述方法得到的西峰、長(zhǎng)武站點(diǎn)MK趨勢(shì)檢驗(yàn)統(tǒng)計(jì)參數(shù)見表1、表2；西峰、長(zhǎng)武站點(diǎn)年降水量隨年份的變化趨勢(shì)線見圖4、圖5。

表1 西峰站點(diǎn)MK趨勢(shì)檢驗(yàn)統(tǒng)計(jì)參數(shù)

表2 長(zhǎng)武站點(diǎn)MK趨勢(shì)檢驗(yàn)統(tǒng)計(jì)參數(shù)

圖4 西峰站點(diǎn)年降水量時(shí)序及變化趨勢(shì)線Fig.4 The timing and the trend line of annual precipitation at Xifeng station

圖5 長(zhǎng)武站點(diǎn)年降水量時(shí)序及變化趨勢(shì)線Fig.5 The timing and the trend line of annual precipitation at Changwu station

應(yīng)用MK檢驗(yàn)法檢驗(yàn)分析水文時(shí)間序列時(shí)，在0.05的顯著性水平下，Z≥|Z1-α/2|拒絕原假設(shè)具有統(tǒng)計(jì)學(xué)意義。顯著性水平α取0.05時(shí)，拒絕原假設(shè)說明存在顯著的變化趨勢(shì)；顯著性水平α取0.01時(shí)，拒絕原假設(shè)說明存在非常顯著的變化趨勢(shì)[24]。

西峰站點(diǎn)，MK分析中位數(shù)β=0.553>0，年降水量變化成上升趨勢(shì)，上升幅度為5.53 mm/10 a，西峰站點(diǎn)年降水的Mann-Kendall趨勢(shì)檢驗(yàn)值Z=0.550<1.96，未通過置信度95 %的顯著性趨勢(shì)檢驗(yàn)，西峰站點(diǎn)年降水量呈上升趨勢(shì)，但未表現(xiàn)出顯著性特點(diǎn)。長(zhǎng)武站點(diǎn)，MK分析中位數(shù)β=1.42>0，年降水量的變化成上升趨勢(shì)，上升幅度為14.2 mm/10 a，長(zhǎng)武站點(diǎn)年降水的Mann-Kendall趨勢(shì)檢驗(yàn)值Z=1.542<1.96，未通過置信度95 %的顯著性趨勢(shì)檢驗(yàn)，長(zhǎng)武站點(diǎn)年降水量呈上升趨勢(shì)，趨勢(shì)不顯著。

綜上分析，董志塬地區(qū)1959—2016年年降水量呈不顯著上升趨勢(shì)，年降水量以5.54～14.2 mm/10 a的速率增加。

3.2 MK突變檢驗(yàn)結(jié)果與分析

通過前述MK突變檢驗(yàn)方法，繪出的UFk和UBk曲線圖如下：

圖6 西峰站點(diǎn)年降水量MK統(tǒng)計(jì)曲線圖Fig.6 The MK statistical curve of annual precipitation at Xifeng station

圖7 長(zhǎng)武站點(diǎn)年降水量MK統(tǒng)計(jì)曲線Fig.7 The MK statistical curve of annual precipitation at Changwu station

從西峰站點(diǎn)年降雨量MK檢驗(yàn)的UFk趨勢(shì)線中可以清晰看出，20世紀(jì)年降水量變化趨勢(shì)：20世紀(jì)60年代初期下降，然后上升，到70年代初呈下降趨勢(shì)，70年代后至90年代初一直呈現(xiàn)上升趨勢(shì)，90年代末至00年代初呈下降趨勢(shì)，之后又呈現(xiàn)上升趨勢(shì)。從UFk和UBk的交點(diǎn)可以看出年降水量突變現(xiàn)象，降水量突變年份分別為1962年、1970年、1985年、2003年、2009年。從長(zhǎng)武站點(diǎn)年降水量MK檢驗(yàn)的UFk趨勢(shì)線中可以清晰看出，20世紀(jì)60年代后，長(zhǎng)武年降水量一直處于上升趨勢(shì)。從UFk和UBk的交點(diǎn)可以看出年降水量突變現(xiàn)象發(fā)生在1963年、1971年、1979年、1995年。

綜上得出，董志塬降水量突變起始年份分別為1962年、1970年、1985年、1995年、2003年、2009年。最小與最大降水的年份為1995年、2003年，這正好與年降水量突變起始年份1995年與2003年相對(duì)應(yīng)。

3.3 小波分析結(jié)果

應(yīng)用Morlet小波對(duì)董志塬氣象站點(diǎn)年降水時(shí)序進(jìn)行小波分析，其中年降水時(shí)序小波功率譜中的黑色實(shí)線為影響錐(COI)，COI考慮了小波函數(shù)與信號(hào)兩端數(shù)據(jù)卷積所產(chǎn)生的計(jì)算結(jié)果影響范圍，其內(nèi)部為可信分析結(jié)果。黑色虛線內(nèi)則為顯著性90%的功率譜。小波系數(shù)實(shí)部時(shí)頻分布圖可以得到降水時(shí)序的時(shí)間、頻率變化特征。在時(shí)間尺度一定的情況下，降水量在該尺度下的變化特征通過小波系數(shù)實(shí)部隨時(shí)間的變化來反映。正值對(duì)應(yīng)降水量偏多的情況，負(fù)值對(duì)應(yīng)降水量偏少的情況。年降水量時(shí)間序列的主要時(shí)間尺度(主周期)，可以通過降水時(shí)序小波功率譜來確定。

3.3.1 西峰年降水量的小波分析結(jié)果

西峰的情況如圖8及圖9所示，可見其存在三種時(shí)間尺度變化周期，分別為1～2年、7～8年、9～14年時(shí)間尺度。從圖9可以看出，影響錐內(nèi)的可信分析結(jié)果，時(shí)間尺度7～8年在能量上占了主要部分，是影響未來該地區(qū)降水量的主要時(shí)間尺度，該時(shí)間尺度在20世紀(jì)80年代趨勢(shì)比較明顯，震蕩中心顯示其枯水與豐水期的交替周期約為8年。而1～2年的周期主要存在于20世紀(jì)60年代早期一段時(shí)間，9～14年的周期主要存在于20世紀(jì)80年代末至90年代末。

圖8 西峰站點(diǎn)年降水時(shí)序小波系數(shù)實(shí)部時(shí)頻分布圖Fig.8 Real part time-frequency distribution of wavelet coefficients of annual precipitation at Xifeng station

圖9 西峰站點(diǎn)年降水時(shí)序小波功率譜Fig.9 Wavelet power spectrum of annual precipitation time-series at Xifeng station

3.3.2 長(zhǎng)武年降水量的小波分析結(jié)果

對(duì)于長(zhǎng)武地區(qū)，根據(jù)圖10、圖11可見其存在三種時(shí)間尺度變化周期，分別為1～2年、4～5年、7～9年時(shí)間尺度。時(shí)間尺度7～9年在能量上占了主要部分，是影響未來該地區(qū)降水量的主要時(shí)間尺度，該時(shí)間尺度在20世紀(jì)70年代初到80年代末趨勢(shì)顯著，震蕩中心顯示其枯水與豐水期的交替周期約為8年。1～2年左右的周期，主要存在于20世紀(jì)60年代初和70年代末。而4～5年左右的周期主要存在于20世紀(jì)80年代末和21世紀(jì)00年代末。

圖10 長(zhǎng)武站點(diǎn)年降水時(shí)序小波系數(shù)實(shí)部時(shí)頻分布圖Fig.10 Real part time-frequency distribution of wavelet coefficients of annual precipitation at Changwu station

圖11 長(zhǎng)武站點(diǎn)年降水時(shí)序小波功率譜Fig.11 Wavelet power spectrum of annual precipitation time-series at Changwu station

綜合分析得到，董志塬地區(qū)年降水量受2年、4年和8年的尺度波動(dòng)所影響，時(shí)間尺度8年在能量上占了主要部分，是影響未來該地區(qū)降水量的主要時(shí)間尺度，主要出現(xiàn)在20世紀(jì)70年代到80年代末；董志塬地區(qū)年降水量時(shí)間序列的主周期為8年，以8年為周期振蕩，震蕩中心在1980年。董志塬地區(qū)降水的豐枯交替明顯，董志塬地區(qū)整個(gè)時(shí)間序列上的年降水量呈現(xiàn)出上升下降的循環(huán)交替特征。

4 結(jié)論

1)董志塬年降水量呈現(xiàn)上升趨勢(shì)，但趨勢(shì)不顯著。

2)董志塬年降水量存在突變性，突變起始年份分別為1962年、1970年、1985年、1995年、2003年、2009年，突變年份間隔為10年左右。

3)董志塬年降水量主要受2年、4年和8年的尺度周期波動(dòng)影響，第一主周期為8年。

4)董志塬年降水在整個(gè)研究的時(shí)域上豐水期與枯水期交替變化明顯，以8年為周期，年降水量呈現(xiàn)上升下降周期性變化規(guī)律。