黃 琨，李 昂，劉融融，宗長榮，王 琪

（江蘇省水文水資源勘測局鹽城分局，江蘇 鹽城 224051）

1 研究區(qū)域

淮河流域[1]地處我國東部，介于長江、黃河兩大流域之間，面積約為27 萬km2?；春恿饔虻乩砦恢媒橛谀媳敝g，氣候條件相對復(fù)雜，流域地勢低平，平原遼闊，洼地大面積分布，受澇風(fēng)險大，加之歷史上黃河長期奪淮，淮河逐漸失去了獨(dú)立的入海通道，造成水系出現(xiàn)紊亂，環(huán)境不斷惡化，水旱災(zāi)害逐年加重，這些均導(dǎo)致淮河流域成為一個水旱災(zāi)害頻生的區(qū)域。因此，對淮河流域降水演變規(guī)律的研究是十分必要的。

2 研究數(shù)據(jù)

本文采用的降水量資料為中國氣象科學(xué)數(shù)據(jù)共享服務(wù)網(wǎng)提供的《中國地面氣候資料月值數(shù)據(jù)集》[2]，該數(shù)據(jù)集包括了平均本站氣壓平均氣溫、平均最高氣溫、平均最低氣溫、平均水汽壓、平均相對濕度等眾多氣象要素，本文選取降水量為研究數(shù)據(jù)進(jìn)行分析。全系列降水資料均經(jīng)過嚴(yán)格的質(zhì)量控制和檢查，錯誤記錄也進(jìn)行了修正，統(tǒng)計(jì)的結(jié)果通過了極值時間一致性檢驗(yàn)，人工抽查，較為可靠。

本文選取淮河流域內(nèi)46 個基本、基準(zhǔn)地面氣象觀測站及自動站1951 年～2016 年的月值數(shù)據(jù)集進(jìn)行整理，去除部分?jǐn)?shù)據(jù)使得每個站的數(shù)據(jù)均從某年的1 月份開始，于當(dāng)年的12 月份結(jié)束，并對缺測的數(shù)據(jù)進(jìn)行線性插值，以保證降水序列的完整性。按照各站點(diǎn)降水序列起始年份進(jìn)行降序排列，并根據(jù)站網(wǎng)的變化以及站點(diǎn)資料的長度將1951 年～2016 年動態(tài)地劃分為15 個區(qū)段（節(jié)點(diǎn)），采用動態(tài)泰森多邊形法計(jì)算得出每個站點(diǎn)每個不同年份的權(quán)重，最終得出淮河流域1951 年～2016 年的逐年降水量。

圖1 淮河流域1951 年（左）及1995 年（右）各站點(diǎn)泰森多邊形

3 淮河流域降水演變規(guī)律分析

3.1 趨勢分析

為更加準(zhǔn)確地對淮河流域降水序列進(jìn)行趨勢分析，本文選取三種方法分別進(jìn)行計(jì)算分析。

（1）Spearman 秩次檢驗(yàn)法

對給定的時間周期Y1……YN及其相應(yīng)值進(jìn)行X1……XN升序排列，統(tǒng)計(jì)檢驗(yàn)所采用的秩相關(guān)系數(shù)按下式計(jì)算[3]：

式中：Xi為周期1-N 按照流量值升序排列序號；Yi為按照時間序列排序的序號；di為Xi與Yi之差值；Rs是否不同于0，選用統(tǒng)計(jì)量t 檢驗(yàn)；H0為無趨勢。

如果T 小于0，那么序列呈趨勢上升；如果T 大于0，那么序列呈趨勢下降。

根據(jù)淮河流域1951 年～2016 年年降水量數(shù)據(jù)，通過Spearman 秩次相關(guān)檢驗(yàn)法，計(jì)算結(jié)果見表1。

表1 Spearman 秩次相關(guān)檢驗(yàn)結(jié)果

（2）Kendall 秩次檢驗(yàn)法

Kendall 秩次相關(guān)法[4]按以下步驟進(jìn)行：

1）H0：假定序列 X1，X2，…，XN無變化趨勢；

2） 計(jì)算 P 值，P 為檢驗(yàn) X1，X2，…，XN序列中所有對偶值（Xi，Yj）（i＜j）當(dāng)中Xi＜Yj）的數(shù)量；

3）構(gòu)造統(tǒng)計(jì)量u：

當(dāng) n→∞，u～N（0，1），即 u 符合標(biāo)準(zhǔn)化正態(tài)分布。

根據(jù)淮河流域1951 年～2016 年年降水量數(shù)據(jù)，通過Kendall 秩次相關(guān)檢驗(yàn)法，計(jì)算結(jié)果見表2。

表2 Kendall 秩次相關(guān)檢驗(yàn)結(jié)果

（3）線性傾向估計(jì)

給定某樣本時間序列為x1，x2，…，xn，采用直線方程線性地?cái)M合它的變化趨勢，即：

式中：a0為回歸常數(shù)；a1為回歸系數(shù)，也可稱為傾向值；用最小二乘法對a0和a1進(jìn)行估計(jì)。

a1的正負(fù)表征 x 的變化傾向，若a1＞0，則 x 趨勢上升，反之則趨勢下降。變量和時間的線性相關(guān)密切程度采用相關(guān)系數(shù)r進(jìn)行表示。

式中：n 為樣本度；xi表示第 i 年的降水量表示降水量時間序列的均值ˉ=（n+1）/2。

當(dāng) r＞0 時，降水量序列有線性增的趨勢；當(dāng) r＜0 時，降水量序列有線性減的趨勢。當(dāng)r=0 時，回歸系數(shù)也為0，表明x 的變化與時間無相關(guān)；越接近0，x 與t 之間的線性相關(guān)就越小，反之相關(guān)性就越密切。

根據(jù)淮河流域1951 年～2016 年年降水量數(shù)據(jù)，通過線性傾向估計(jì)法，計(jì)算結(jié)果見圖2。

圖2 淮河流域逐年降水量線性傾向估計(jì)趨勢過程線

淮河流域線性傾向估計(jì)計(jì)算結(jié)果總結(jié)見表3。

表3 線性傾向估計(jì)計(jì)算結(jié)果

3.2 突變分析

（1）累計(jì)距平法

距平在本文中的含義為降水量序列與多年均值之差，累積距平則是將序列中某一年之前的降水距平值進(jìn)行累加，并繪制成累積過程線，再通過過程線來判斷降水量的變化趨勢[5]。計(jì)算公式如下：

淮河流域1951 年～2016 年年降水量序列累計(jì)距平過程見圖3，累積距平曲線呈上升表示距平值增加，豐水年發(fā)生或者繼續(xù)；呈下降趨勢表示距平值減小，枯水年發(fā)生或者繼續(xù)；拐點(diǎn)表示豐枯年的轉(zhuǎn)折或者突變。

圖3 淮河流域逐年降水量累計(jì)距平曲線

表4 累計(jì)距平法計(jì)算結(jié)果

（2）M-K 檢驗(yàn)法

M-K 檢驗(yàn)法[6]是目前在時間序列趨勢分析的非參數(shù)檢驗(yàn)中應(yīng)用較為廣泛的方法，最早由Mann 于1945 年提出，后經(jīng)過Kendall 和Sneyers 等人的完善，該方法能夠大致判斷變化趨勢開始發(fā)生的時間節(jié)點(diǎn)，Gooddens 等人又將其應(yīng)用到反序列當(dāng)中，進(jìn)一步發(fā)展成為可以檢驗(yàn)時間序列突變的方法。

M-K 檢驗(yàn)法的實(shí)質(zhì)是通過數(shù)據(jù)序列的秩判斷兩個變量間的相關(guān)程度。對于樣本容量為n 的時間序列x，構(gòu)造一個秩序列為：

其中

從上式可以看出，sk表示第i 時刻數(shù)值大于第j 時刻數(shù)值個數(shù)的累加。

假定時間序列隨機(jī)獨(dú)立，定義統(tǒng)計(jì)量為：

其中UF1=0，E（sk）為sk的均值，Var（sk）為sk的方差。在x1，x2，…xn，相互獨(dú)立，并且具有相同連續(xù)分布式的條件下，兩者可利用下列公式計(jì)算得出：

UF1服從標(biāo)準(zhǔn)正態(tài)分布，是按時間序列 x 順序，x1，x2，…xn，計(jì)算得出的統(tǒng)計(jì)量序列，給定顯著性水平α，查看正態(tài)分布表，如果表示序列有顯著變化趨勢。

將序列 x 按照時間逆序排列為 xn，xn-1，…，x1，再重復(fù)以上過程，并使 UBk=－UFk，k=n，n-1，…，1，UB=0。

這種方法不僅計(jì)算簡便，而且可以明確判定突變發(fā)生的起始時間，并且能夠指出突變區(qū)域。

采用M-K 檢驗(yàn)法計(jì)算并繪制出淮河流域1951 年～2016 年降水量序列的M-K 統(tǒng)計(jì)曲線見圖4。

圖4 淮河流域逐年降水量曼- 肯德爾統(tǒng)計(jì)曲線

從淮河流域1951 年～2016 年逐年降水量M-K 統(tǒng)計(jì)曲線可以看出，UF 和UB 兩條曲線的交點(diǎn)都在臨界線之間，但是UF 曲線始終沒有超過信度線，說明淮河流域年降水序列無顯著突變。

（3）滑動 t 檢驗(yàn)法

滑動t 檢驗(yàn)[7]是通過兩組樣本的均值之差異的顯著性來檢驗(yàn)序列的突變，人為地設(shè)定一個基準(zhǔn)時間點(diǎn)，如果該時間基準(zhǔn)點(diǎn)前后樣本的均值差異超過了一定的顯著性水平，那么就認(rèn)為序列發(fā)生了突變。計(jì)算公式如下：

式中：n 是樣本容量，x1和x2表示時間基準(zhǔn)點(diǎn)前后的兩段子樣本，n1和n2分別是兩段子樣本的容量表示兩段樣本的均值為方差。對于統(tǒng)計(jì)量t，首先給定顯著性水平α，如果t 超過臨界值，就認(rèn)為在設(shè)置的基準(zhǔn)時間點(diǎn)處發(fā)生了突變，反之則無突變。

采用滑動t 檢驗(yàn)法計(jì)算并繪制出淮河流域1951 年～2016 年逐年降水量滑動t 檢驗(yàn)曲線，見圖5。

圖5 淮河流域逐年降水量滑動t 檢驗(yàn)曲線

計(jì)算結(jié)果見表5。

表5 線性傾向估計(jì)計(jì)算結(jié)果

3.3 周期分析

為了更好地研究時間序列問題，Morlet 在20 世紀(jì)80 年代初提出了具有時- 頻多分辨功能的小波分析方法，該方法可以清晰地揭示隱藏于時間序列當(dāng)中的多種變化周期，也能夠充分反映序列在不同時間尺度上的變化趨勢，同時能對序列未來發(fā)展趨勢進(jìn)行定性估計(jì)。

小波分析理論目前廣泛地應(yīng)用于信號處理、數(shù)值分析以及大氣科學(xué)非線性科學(xué)領(lǐng)域。在時間序列研究中，小波分析主要用來對時間序列進(jìn)行消噪與濾波，計(jì)算信息量系數(shù)與分形維數(shù)，突變點(diǎn)的監(jiān)測和周期成分的識別以及多時間尺度的分析等[8]。

（1）小波函數(shù)

小波分析的實(shí)質(zhì)是利用一簇小波函數(shù)系來表示或逼近某個信號或者函數(shù)。因此，小波函數(shù)在小波分析中尤為關(guān)鍵，它具有震蕩性、能夠迅速衰減到零的特性，即小波函數(shù)ψ（t）∈L2（R）并且滿足：

式中：ψ（t）為基小波函數(shù)，可以通過尺度的伸縮以及時間軸上的平移形成一簇函數(shù)系：

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

（2）小波變換

如果ψa,b（t）是（2）式中給定的子小波，那么對于能量有限信號f（t）∈L2（R），其連續(xù)小波變換為：

由以上可知小波分析的基本原理，即通過尺度因子a 的增加或減小來得到信號的高低頻信息，進(jìn)而分析信號的概況或者細(xì)節(jié)，實(shí)現(xiàn)對信號在不同時間尺度和空間局部特征的分析。

在實(shí)際研究當(dāng)中，由小波變換方程得到小波系數(shù)是最主要的步驟，然后再通過這些系數(shù)對時間序列的時頻變化特征進(jìn)行分析。

（3）小波方差

對小波系數(shù)的平方值在b 域上進(jìn)行積分，即可得到小波方差：

小波方差隨伸縮尺度a 變化而變化的過程即為小波方差圖。由式（20）可知，它可以反映信號波動能量隨著尺度因子a的分布情況。因此，可利用小波方差圖來確定信號中不同尺度擾動的相對強(qiáng)度以及存在的主要時間尺度，即主周期。

根據(jù)淮河流域1951 年～2016 年的降水?dāng)?shù)據(jù)，利用Morlet小波分析方法對其進(jìn)行周期性分析，結(jié)果見圖6。

圖5 小波方差圖、小波系數(shù)實(shí)部等值圖

對淮河流域1951 年～2016 年年降水量作Morlet 小波變化得到小波系數(shù)實(shí)部等值圖，由圖5 可以清晰看出該地區(qū)年降水量數(shù)據(jù)在不同時間尺度上的周期震蕩。由年降水量的Morlet小波系數(shù)時頻分布圖和小波方差圖可以看出年降水量存在如下的周期演變特征：

淮河流域年降水量變化可能存在4 a、14 a、27 a、55 a 的周期，其中55 a 為第一峰值，對應(yīng)實(shí)部圖中“多- 少- 多”的變化過程，震蕩最強(qiáng)，但是完整性不足，需要更長的降水序列進(jìn)行驗(yàn)證，只能初步判定55 a 可能是淮河流域年降水量變化的一個周期；27 a 的變化周期主要從20 世紀(jì)80 年代開始發(fā)生，之后一直存在；而14 a 的變化周期基本上一直存在；4 a 的變化周期持續(xù)發(fā)生之后與21 世紀(jì)初消失。

4 結(jié)論

（1）淮河流域1951 年～2016 年降水量序列整體呈下降趨勢，但是顯著性不強(qiáng)。突變分析得出的結(jié)果也不顯著，分割時間為20 世紀(jì)70 年代中期，可見淮河流域整體降水較為穩(wěn)定。

（2）周期分析結(jié)果得出淮河流域存在55 a 的較大第一主周期，對應(yīng)實(shí)部圖中“多- 少- 多”的變化過程，震蕩最強(qiáng)，但是完整性不足，需要更長的降水序列進(jìn)行驗(yàn)證。