基于希爾伯特-黃變換與小波分析的降雨序列多時(shí)間尺度研究

賀軍奇,郭鑫佳,陳云飛,劉秀花,高萬德,龍婷

(1.長(zhǎng)安大學(xué)水利與環(huán)境學(xué)院,西安 710054;2.旱區(qū)地下水與生態(tài)效應(yīng)教育部重點(diǎn)實(shí)驗(yàn)室,西安 710054;3.水利部旱區(qū)生態(tài)水文與水安全重點(diǎn)實(shí)驗(yàn)室,西安 710054)

降雨是一個(gè)復(fù)雜的自然過程,受到諸多因素的影響,其序列往往是非線性、非平穩(wěn)的,并伴隨著一定的隨機(jī)性和周期性[1]。特別是近年來全球氣候變暖和人類活動(dòng)的疊加影響,不僅加快了區(qū)域尺度上的水文循環(huán)過程,也加劇了降雨的區(qū)域差異性和不確定性[2]。針對(duì)降雨序列表現(xiàn)出的多時(shí)間尺度、多頻率的變化特性,水文周期分析方法經(jīng)歷了從時(shí)間域到時(shí)頻域的發(fā)展歷程[3]。其中,Fourier變換因其時(shí)頻分析特性自1965年來便得到廣泛應(yīng)用,但受Heisenberg不確定性原理的限制[4],Fourier變換不能刻畫任意時(shí)刻的頻率成分,會(huì)在某種尺度上掩蓋序列特征[5]。隨后以其為基礎(chǔ)發(fā)展的小波分析被更多地用于研究降雨序列的多尺度變化規(guī)律和演變趨勢(shì),雖然小波分析使用可調(diào)時(shí)頻窗解決了時(shí)間和頻率分辨率的矛盾,具有良好的時(shí)頻局部化特征和多分辨率能力,但其本質(zhì)上仍是一種窗口可調(diào)的Fourier變換,依然存在Fourier變換的局限[6]。

希爾伯特-黃變換(Hilbert-Huang Transform,HHT)是Huang等在1998年提出的一種時(shí)頻分析方法,主要由經(jīng)驗(yàn)?zāi)B(tài)分解(Empirical Mode Decomposition,EMD)和希爾伯特譜分析(Hilbert Transform,HT)兩部分組成:通過EMD 對(duì)非線性、非平穩(wěn)信號(hào)逐級(jí)進(jìn)行線性化和平穩(wěn)化處理,在保留數(shù)據(jù)特性的同時(shí)把不同尺度的波動(dòng)和變異分離出來;隨后通過HT 求解瞬時(shí)頻率進(jìn)一步強(qiáng)化數(shù)據(jù)的局部特征,從而精確地給出頻率與時(shí)間的對(duì)應(yīng)關(guān)系,使HHT 不再受不確定性原理的限制[7]。因此該方法在目前的水文序列多尺度研究中得到了較多應(yīng)用。如周生輝等[8]通過EMD 分析了1959—2019年海流兔河流域周邊3個(gè)氣象站的降雨周期變化特征;Luo 等[9]利用HHT 識(shí)別徑流與產(chǎn)沙量的多時(shí)間尺度特征,并通過描述兩個(gè)變量在不同尺度上的關(guān)聯(lián)性闡明其內(nèi)在的產(chǎn)沙機(jī)制;邵駿等[10]利用HHT 對(duì)雅魯藏布江干支流上的10個(gè)水文站1956—2000年的天然徑流量進(jìn)行分析,探討了雅魯藏布江流域年徑流變化的近似周期及其演變趨勢(shì)。此外,為提高周期特征的識(shí)別精度,近年來一些研究也嘗試將改進(jìn)后的HHT 引入水文分析中。如范琳琳等[11]將基于T 檢驗(yàn)改進(jìn)的HHT 用于長(zhǎng)江宜昌水文站50 a日徑流量的分析中,識(shí)別出徑流在不同時(shí)間尺度上的年際變化周期;余世鵬等[12]為了消除HHT 可能存在的模態(tài)混疊現(xiàn)象,利用基于集合經(jīng)驗(yàn)?zāi)B(tài)分解(Ensemble Empirical Mode Decomposition,EEMD)改進(jìn)的H HT 剖析了濱海圍墾灘涂區(qū)域降水的多尺度周期特征,并指出改進(jìn)的HHT 能獲取更準(zhǔn)確的周期尺度特征。

然而,目前大多數(shù)研究都側(cè)重于HHT 的實(shí)例應(yīng)用,缺少與傳統(tǒng)周期識(shí)別方法的對(duì)比,同時(shí)也缺乏分析區(qū)域尺度上降雨周期演變特征的可靠方法。因此,本文基于多元經(jīng)驗(yàn)?zāi)B(tài)分解(Multivariate Empirical Mode Decomposition,MEMD)獲取多變量共同尺度的優(yōu)勢(shì),提出將其應(yīng)用于區(qū)域尺度降雨周期演變特征的提取。在此基礎(chǔ)上,采用2組已知成分的人工序列對(duì)比分析基于MEMD 改進(jìn)的HHT(后稱HHT)和傳統(tǒng)小波分析在周期識(shí)別準(zhǔn)確性上的差異;隨后通過毛烏素沙地11個(gè)氣象站點(diǎn)的實(shí)測(cè)降雨序列,進(jìn)一步探討HHT 在推求區(qū)域降雨序列周期演變特征上的適用性。本研究不僅能為區(qū)域水文周期分析方法的選擇提供更多參考,同時(shí)也通過研究毛烏素沙地降雨周期的演變特征,為旱澇災(zāi)害應(yīng)對(duì)和區(qū)域水資源管理提供一定的科學(xué)依據(jù)。

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

1.1 HHT的基本理論及特性

EMD 作為原始HHT 的基本理論,可以將復(fù)雜的原始數(shù)據(jù)分解為若干固有經(jīng)驗(yàn)?zāi)B(tài)函數(shù)(Intrinsic Mode Function,IMF)和殘差項(xiàng)(Residual)的有限組合。在本研究中,采用MEMD 取代原有的EMD 以改進(jìn)HHT。MEMD 是標(biāo)準(zhǔn)EMD 的多元擴(kuò)展,通過引入映射概念,將含有多個(gè)變量的信號(hào)組在不同空間中沿不同方向投影,生成多維包絡(luò)線[13]。然后計(jì)算多元變量的局部平均值,確保多個(gè)變量的IMF 分量在振蕩頻率和數(shù)量上的同步性,從而實(shí)現(xiàn)不同頻率尺度下多變量的共同模式分解[14]。

假設(shè)代表降雨的n維矢量數(shù)據(jù)S(t)={s1(t),s2(t),…,sn(t)}是數(shù)據(jù)t的函數(shù),Xθk={xk1,xk2,…,xkn}是角度矢量θk={θk1,θk2,…,θkn-1}(k=1,2,…,N)定義的沿著不同方向的矢量數(shù)據(jù),其中N是所有方向的個(gè)數(shù)。

非線性的k個(gè)空間數(shù)據(jù)的MEMD具體步驟如下[15]:

(1)獲取合適的n維方向矢量數(shù)據(jù)集X;

(2)計(jì)算時(shí)間序列S(t)沿著給定方向Xθk的投影pθk(t);

(3)識(shí)別每一個(gè)方向向量上pθk(t)極值出現(xiàn)的位置tθki;

(4)通過多元樣條插值函數(shù)插值〔tθki,S(tθki)〕,獲取多元包絡(luò)曲線eθk(t);

(5)計(jì)算包絡(luò)線均值M(t),公式如下:

(6)通過D(t)=S(t)-M(t)得到IMF。假如D(t)滿足多元IMF 的判斷標(biāo)準(zhǔn),則D(t)作為IMF保存,將S(t)=S(t)-D(t)作為步驟(2)的輸入變量。如果不滿足,則將S(t)=D(t)作為步驟(2)的輸入變量,繼續(xù)進(jìn)行迭代運(yùn)算(圖1)[16]。

圖1 MEMD流程示意圖Fig.1 Schematic diagram of MEMD process

HT 作為HHT 的核心內(nèi)容,其目的是得到能充分反映S(t)時(shí)—頻—能量分布的瞬時(shí)頻率,使分解得到的每個(gè)IMF 都有著明確的物理意義(即時(shí)間尺度周期)[17]。具體計(jì)算過程如下:對(duì)IMFi(t)進(jìn)行HT 轉(zhuǎn)換(公式2);提取瞬時(shí)頻率和振幅(公式3);最終得到能反映信號(hào)能量在時(shí)間和頻率上分布規(guī)律的Hilbert譜(公式4):

式中:P.V.表示Cauchy主值積分,ωi(t)和ai(t)分別為瞬時(shí)頻率和振幅。

1.2 小波分析的基本理論及特性

小波是一種傳統(tǒng)的時(shí)頻多分辨率分析方法,在進(jìn)行時(shí)間序列分析時(shí),通常選用可以較好表達(dá)相位信息的復(fù)值小波。其中,Morlet小波在時(shí)頻空間的局部性較好,因此常被用于降雨徑流時(shí)間序列分析[18]。本研究使用的Morlet小波解析形式為[19]:

式中:i為虛部單位;w0為無量綱頻率,t為時(shí)間。

對(duì)于時(shí)間序列f(t)∈L2(R),其連續(xù)小波變化為:

在時(shí)間域上不同尺度a的所有小波系數(shù)Wf(a,b)的平方值對(duì)位移因子b進(jìn)行積分,得到小波方差的計(jì)算公式:

為減小序列邊界端點(diǎn)效應(yīng)的影響,利用MATLAB小波工具箱中的信號(hào)延伸(Signal Extension)功能,采用對(duì)稱、周期延伸法對(duì)距平后的降雨數(shù)據(jù)兩端進(jìn)行延展。然后,計(jì)算其復(fù)小波系數(shù),刪除延伸數(shù)據(jù)的系數(shù)后,計(jì)算Mrolet復(fù)小波系數(shù)的實(shí)部、模、模方和方差。

1.3 配對(duì)樣本t檢驗(yàn)

配對(duì)樣本檢驗(yàn)是針對(duì)配對(duì)數(shù)據(jù)的t檢驗(yàn),可以用來檢驗(yàn)兩組數(shù)據(jù)是否有差別。設(shè)X1,X2表示兩個(gè)隨機(jī)的成對(duì)樣本,兩個(gè)變量的n對(duì)觀測(cè)值形成配對(duì)樣本,它們可以表示為(X11,X21),(X12,X22),…,(X1n,X2n),計(jì)算每一對(duì)樣本觀察值之差,記為dj,進(jìn)一步得到其均值和標(biāo)準(zhǔn)差Sd。假設(shè)H0:μ1-μ2=μd=0,H1:μd≠0。此時(shí)檢驗(yàn)統(tǒng)計(jì)量為[20]:

在原假設(shè)成立的條件下,它服從自由度為〔n-1〕的t分布。對(duì)應(yīng)的伴隨概率小于設(shè)定的顯著性水平,則拒絕兩組樣本均值相等的假設(shè)。將人工序列的預(yù)設(shè)周期分別與不同方法的計(jì)算結(jié)果進(jìn)行配對(duì)樣本t檢驗(yàn),以分析不同方法在周期識(shí)別方面的準(zhǔn)確性。本次研究顯著性水平設(shè)置為0.05。

1.4 數(shù)據(jù)資料

使用人工序列和實(shí)測(cè)降雨序列綜合對(duì)比HHT和小波分析在多周期尺度分析的準(zhǔn)確性和適用性。參考姜瑤等的研究案例[21],創(chuàng)建2組已知成分的人工序列數(shù)據(jù)測(cè)試對(duì)比H HT 和小波分析在周期識(shí)別方面的準(zhǔn)確性,并探討造成不同方法性能差異的主要原因。第1組數(shù)據(jù)考慮不同的趨勢(shì)變化,并在此基礎(chǔ)上疊加相同的周期和隨機(jī)波動(dòng);第2組數(shù)據(jù)無趨勢(shì)變化,在此基礎(chǔ)上疊加不同的周期和隨機(jī)波動(dòng)。所有序列的長(zhǎng)度均為120 a,詳細(xì)內(nèi)容見表1和圖2。

表1 人工序列數(shù)據(jù)的主要信息Table 1 Main information of artificial sequence data

圖2 人工序列數(shù)據(jù)示意圖Fig.2 Schematic diagram of manual sequence data

實(shí)測(cè)序列采用1982—2020年毛烏素沙地11個(gè)氣象站點(diǎn)的逐日降雨數(shù)據(jù)計(jì)算各站點(diǎn)的逐年降雨量(包含冬季降雪),數(shù)據(jù)源自中國(guó)氣象數(shù)據(jù)網(wǎng)(http:∥data.cma.cn)。獲取數(shù)據(jù)后,統(tǒng)一對(duì)所有站點(diǎn)的降雨數(shù)據(jù)進(jìn)行質(zhì)量檢查與校正,主要包括:檢查降雨數(shù)據(jù)是否有缺失,對(duì)缺失數(shù)據(jù)進(jìn)行插補(bǔ);刪除多余的重復(fù)數(shù)據(jù);對(duì)異常值進(jìn)行修正等。

如圖3所示,毛烏素沙地39 a年均降雨量在170～440 mm:沙地西部降雨較少,年平均降雨量為170～200 mm,包含陶樂站(S6),銀川站(S7)以及吳忠站(S9);中部次之,年平均降雨量為270～350 mm,包含定邊站(S3),鹽池站(S8)以及鄂托克旗站(S10);東部降雨較多,年平均降雨量為360～440 mm,包含榆林站(S1),神木站(S2),靖邊站(S4),橫山站(S5)以及伊金霍洛旗站點(diǎn)(S11)。

圖3 1982-2020年年均實(shí)測(cè)降雨空間分布圖Fig.3 Spatial distribution map of annual average measured rainfall from 1982 to 2020

2 結(jié)果與分析

2.1 人工序列周期識(shí)別結(jié)果

圖4和圖5分別為HHT 對(duì)2組人工序列的分解結(jié)果。括號(hào)內(nèi)的數(shù)字表示各個(gè)IMF和殘差分量占原始序列的方差相對(duì)貢獻(xiàn)率,該值表征了IMF 分量對(duì)原始序列的影響程度,方差貢獻(xiàn)率越大,IMF影響越大,代表性越強(qiáng)。由于數(shù)據(jù)驅(qū)動(dòng)的自適應(yīng)性,第1組人工序列(M1,M2及M3)經(jīng)MEMD 分解后得到了7個(gè)IMF和1個(gè)殘差分量;第2組人工序列(M1,M4及M5)經(jīng)分解后得到了6個(gè)IMF和1個(gè)殘差分量。從IMF1到殘差分量,數(shù)據(jù)振蕩幅度逐漸減小,波長(zhǎng)逐漸增大,體現(xiàn)了從高頻—低頻的轉(zhuǎn)變,反映了各分量所代表的時(shí)間尺度的周期變化特征。以M1為例,經(jīng)HT 變換后得到各模態(tài)(IMF1-7)對(duì)應(yīng)的時(shí)間表征尺度分別為2.9 a,4.7 a,12.1 a,24.1 a,30.4 a,56.4 a以及119.1 a(表2)。通過對(duì)比各序列的所有IMF分量,可以發(fā)現(xiàn)IMF3的方差貢獻(xiàn)率占據(jù)了原始序列的80%以上,這表明兩組人工序列的周期以IMF3分量的周期為主,進(jìn)而得到M1,M2和M3的周期分別為12.1 a,12.2 a和12.2 a;M1,M4和M5的周期分別為12.1 a,14.9 a和10.1 a,與表1的預(yù)設(shè)周期幾乎一致。此外,圖中的殘差分量均能很好地展示圖2原始數(shù)據(jù)的趨勢(shì)變化,這表明MEMD 能很好地提取2組數(shù)據(jù)的局部特征和趨勢(shì)變化?？偟膩碚f,H HT 不僅能準(zhǔn)確識(shí)別5個(gè)變量各個(gè)尺度下的周期特征,還具備很好的趨勢(shì)檢測(cè)性能。

表2 基于人工序列使用HHT得到的各固有模態(tài)函數(shù)周期表Table 2 The periodic table of each intrinsic mode function obtained by HHT based on the artificial generated sequence a

圖4 基于第1組人工序列M1,M2,M3使用HHT分解得到的固有模態(tài)函數(shù)(IMF1-7)和殘差項(xiàng)Fig.4 The intrinsic mode function(IMF1-7)and residual term decomposed by HHT based on the first set of artificially generated sequences M1,M2 and M3

圖5 基于第2組人工序列M1,M4,M5使用HHT分解得到的固有模態(tài)函數(shù)(IMF1-6)和殘差項(xiàng)Fig.5 The intrinsic mode function(IMF1-6)and residual term decomposed by HHT based on the first set of artificially generated sequences M1,M4 and M5

圖6A 和6B分別為人工序列采用不同延拓方式(對(duì)稱和周期延拓)獲取到的小波方差以及對(duì)應(yīng)的小波實(shí)部。小波方差反映了時(shí)間序列的波動(dòng)能量隨時(shí)間尺度的分布,根據(jù)峰值的個(gè)數(shù)可以確定時(shí)間序列中存在的主周期數(shù)量;小波實(shí)部則能反映出不同特征時(shí)間尺度信號(hào)在不同時(shí)間上的分布和相位等信息。結(jié)果表明:(1)當(dāng)采用對(duì)稱延拓的小波分析時(shí),M1,M2和M3的小波方差最大值均于19處取得,識(shí)別出的周期均為12.6 a;M4和M5的小波方差最大值分別于23,15處取得,識(shí)別出的周期為16 a和10.4 a,對(duì)周期的識(shí)別總體上相對(duì)偏大(表1)。(2)相較于前者,采用周期延拓的小波分析識(shí)別出的M1-M5周期與人工序列的預(yù)設(shè)周期基本一致,分別為12.2 a,12.2 a,12.2 a,14.9 a以及10.1 a。

圖6 基于人工序列使用小波分析方法得到的小波方差(A)和小波實(shí)部(B)圖Fig.6 Wavelet variance(A)and wavelet real part(B)obtained by using wavelet analysis method based on artificially generated sequence

為進(jìn)一步檢驗(yàn)各方法的周期結(jié)果與人工序列預(yù)設(shè)周期的精度差異,進(jìn)行了配對(duì)樣本t檢驗(yàn)。統(tǒng)計(jì)結(jié)果表明:HHT 的顯著性為Sig.=0.142＞0.05,采用周期延拓的小波分析的顯著性為Sig.=0.109＞0.05,兩者的周期識(shí)別結(jié)果與預(yù)設(shè)周期一致,不存在顯著差異,且HHT 的周期準(zhǔn)確性要略高于周期延拓的小波分析;相比之下,采用對(duì)稱延拓的小波分析的顯著性為Sig.=0.003＜0.05,識(shí)別出的周期結(jié)果偏大,與預(yù)設(shè)周期存在顯著差異(表3)。總的來說,H HT 與小波分析(周期延拓和對(duì)稱延拓)相比,周期識(shí)別精度最高。

表3 配對(duì)樣本t檢驗(yàn)結(jié)果Table 3 Paired sample t test results

2.2 實(shí)測(cè)序列分析

HHT 將毛烏素沙地11個(gè)降雨站點(diǎn)的時(shí)間序列分解為5個(gè)IMF和1個(gè)殘差分量。IMF的方差貢獻(xiàn)率計(jì)算結(jié)果表明,11個(gè)站點(diǎn)的IMF1和IMF2方差貢獻(xiàn)率最大,兩者貢獻(xiàn)了原始序列41.36%～74.93%的方差(圖7)。其中,毛烏素東南部3 個(gè)站點(diǎn)S1(圖7A)、S4(圖7D)和S5(圖7E)IMF2的振幅均在2000年之后表現(xiàn)出了明顯增加,這表明這3個(gè)站點(diǎn)的降雨波動(dòng)經(jīng)歷了平緩—?jiǎng)×业倪^程,且當(dāng)前正處在降雨量變化較為強(qiáng)烈的階段,年際變化不穩(wěn)定。而其余站點(diǎn)主周期對(duì)應(yīng)的IMF 振幅總體較平緩,無異常波動(dòng)。此外,圖7L還展示了11個(gè)站點(diǎn)1982—2020年的降雨趨勢(shì)變化特征:西部站點(diǎn)(S6,S7和S9)的殘差項(xiàng)低于200,變化趨勢(shì)不明顯;中部站點(diǎn)(S3,S8和S10)的殘差位于300左右,有輕微上升趨勢(shì);而東部站點(diǎn)(S1,S2,S4,S5 及S11)的殘差項(xiàng)位于300～500 之間,上升趨勢(shì)明顯?？傮w上來看,殘差項(xiàng)表現(xiàn)為:西部＜中部＜東部,與前文所述的降雨量空間分布格局一致。通過對(duì)實(shí)測(cè)降雨的分析發(fā)現(xiàn),H HT 方法不僅可以了解降雨量隨時(shí)間的波動(dòng)情況,還能通過殘差項(xiàng)準(zhǔn)確把握序列的變化趨勢(shì)。

圖7 基于實(shí)測(cè)降雨序列使用HHT得到的主要固有模態(tài)函數(shù)(IMF1-2)和殘差項(xiàng)Fig.7 Main intrinsic mode function(IMF1-2)and residual term decomposed by HHT based on measured rainfall series

表4展示了HHT和小波分析計(jì)算出的各站點(diǎn)主周期。(1)HHT結(jié)果表明:11個(gè)站點(diǎn)的第一主周期范圍為2.7～5.5 a,第二主周期為2.8～5.6 a。(2)小波分析(對(duì)稱延拓)的結(jié)果表明:11個(gè)站點(diǎn)中的4個(gè)站點(diǎn)識(shí)別出了3個(gè)主周期(S2,S6,S9和S10),3個(gè)站點(diǎn)識(shí)別出了2個(gè)主周期(S1,S3和S11),剩余4個(gè)站點(diǎn)僅識(shí)別出了1個(gè)主周期(S4,S5,S7和S8),且站點(diǎn)第一主周期范圍為4.6～13.0 a。(3)而小波分析(周期延拓)的結(jié)果表明:除S10識(shí)別出2 個(gè)主周期外,其余站點(diǎn)均識(shí)別出了3個(gè)主周期,且第一主周期在5.6～21.0 a之間,第二主周期在5.6～19.5 a之間。圖8A從3種方法識(shí)別出的周期來看,采用不同延拓方式的小波分析得出的周期結(jié)果在空間分布上無明顯規(guī)律,無論是采用對(duì)稱延拓還是周期延拓,絕大多數(shù)站點(diǎn)的小波分析主周期依然偏大,且在空間上差異大,尤其在相鄰站點(diǎn)出現(xiàn)了異常,如S8與S10的對(duì)稱延拓結(jié)果,以及S5,S2與S1的周期延拓結(jié)果。而H HT可以識(shí)別出全部站點(diǎn)相同數(shù)量的主周期,這主要是因?yàn)镸EMD克服了多變量數(shù)據(jù)中的模式對(duì)齊問題,能獲取多個(gè)站點(diǎn)的共同尺度,更適合區(qū)域尺度降雨周期的求解。

表4 基于實(shí)測(cè)降雨序列使用HHT和小波分析得到的主周期表Table 4 The main periodic table obtained by using HHT and wavelet analysis based on the measured rainfall series

圖8 基于實(shí)測(cè)降雨序列的主周期空間分布Fig.8 Spatial distribution of rainfall main period based on measured rainfall series

此外,采用H HT 還發(fā)現(xiàn)毛烏素沙地的降雨主周期在空間上具有明顯的區(qū)域差異,即腹地區(qū)域的降雨主周期大于沙地邊緣地區(qū)。由圖8B 可知,毛烏素腹地站點(diǎn)(包含S1,S4,S5,S10)的第一主周期均為5.5 a,第二主周期在2.8～3.6 a之間;而沙地邊緣或臨近黃河站點(diǎn)的降雨序列第一主周期為2.7～3.3 a,第二主周期在4.9～5.6 a之間。

3 討論

本研究利用人工序列和實(shí)測(cè)降雨序列,對(duì)H HT和小波分析在求解區(qū)域尺度降雨周期演變特征的適用性上進(jìn)行了測(cè)試、對(duì)比和分析。人工序列的結(jié)果可為周期識(shí)別的準(zhǔn)確性提供依據(jù),而實(shí)測(cè)降雨數(shù)據(jù)為研究區(qū)域多尺度的周期演變提供應(yīng)用實(shí)例,兩者結(jié)果相輔相成。

在測(cè)試對(duì)比中,HHT 能準(zhǔn)確識(shí)別人工序列的周期,同時(shí)也能很好地反映數(shù)據(jù)的趨勢(shì)變化;而相較于對(duì)稱延拓,采用周期延拓的小波分析能更準(zhǔn)確地識(shí)別人工序列周期,并能有效減小端點(diǎn)效應(yīng)造成的周期識(shí)別誤差。造成上述現(xiàn)象的原因主要有以下兩個(gè)方面:一是小波分析相較于其他分解方法(經(jīng)驗(yàn)?zāi)B(tài)分解EMD,變分模態(tài)分解VMD 等),受Heisenberg測(cè)不準(zhǔn)原理影響,只能分析出近似周期,缺乏自適應(yīng)性[22];二是在本研究中,人工序列的預(yù)設(shè)周期是以正弦函數(shù)為基底,有明顯規(guī)律波動(dòng),這使得序列端點(diǎn)處按照周期延拓更加合理,由此得出的結(jié)果也幾乎和預(yù)設(shè)周期一致。證實(shí)了不同延拓方式的選取在一定程度上會(huì)導(dǎo)致數(shù)據(jù)周期結(jié)果失準(zhǔn),這也從側(cè)面論證了若要使用小波分析得出準(zhǔn)確的周期結(jié)果,許多參數(shù)選取(如小波基、延拓方式等)都需要慎重考慮。除了上文闡述的主要原因外,本研究還發(fā)現(xiàn)當(dāng)小波方差含有多個(gè)主峰或主峰不明顯時(shí),很難準(zhǔn)確識(shí)別主周期下的周期尺度,這可能也是小波分析周期偏大的主要原因[23]。總的來說,測(cè)試分析的結(jié)果表明相較于傳統(tǒng)的小波分析,HHT 因其自適應(yīng)性分解能力使之在水文時(shí)間序列的周期識(shí)別和趨勢(shì)檢測(cè)方面具有較為廣闊的應(yīng)用前景。

在實(shí)例應(yīng)用中,HHT 和小波分析(對(duì)稱延拓)識(shí)別出的周期結(jié)果與前人的研究基本一致,而小波分析(周期延拓)則偏差較大[24-25]。與人工序列不同,造成小波分析(周期延拓)識(shí)別效果不佳的原因可能是因?yàn)閷?shí)測(cè)降雨序列往往沒有明顯的波動(dòng)規(guī)律,常呈現(xiàn)多時(shí)間尺度、多頻率的變化特性,致使周期延拓會(huì)以整個(gè)時(shí)間序列為單位進(jìn)行延展,進(jìn)而忽略序列內(nèi)部蘊(yùn)藏的短時(shí)間尺度周期。相比之下,H HT 和小波分析(對(duì)稱延拓)識(shí)別出的降雨主周期也存在一定差異,主要體現(xiàn)在兩方面:一是H HT 識(shí)別出毛烏素沙地所有站點(diǎn)均有2個(gè)降雨主周期;而小波分析(對(duì)稱延拓)得出的各站點(diǎn)降雨主周期數(shù)量不等,在1～3之間。二是從計(jì)算結(jié)果的數(shù)值來看,HHT 計(jì)算得到的毛烏素沙地降雨主周期平均為3.9 a,且各站點(diǎn)的主周期較為接近;而采用對(duì)稱延拓的小波分析得到的降雨主周期平均為7.1 a,但不同站點(diǎn)之間差別較大,且個(gè)別站點(diǎn)還高達(dá)10.3 a(S6)。這些差異主要依賴于MEMD方法在尋找多變量共同尺度上的優(yōu)勢(shì),這不僅有助于實(shí)現(xiàn)區(qū)域尺度上降雨周期的求解,同時(shí)也能避免區(qū)域尺度上因個(gè)別站點(diǎn)的極端變化導(dǎo)致的周期失準(zhǔn)問題[26-27]。此外,本研究還發(fā)現(xiàn)毛烏素東南緣站點(diǎn)(S1,S4和S5)的降雨量波動(dòng)在2000年后變得更為劇烈,這可能是這些區(qū)域在2000年后大降雨事件頻發(fā)所導(dǎo)致的[28]。在空間分布上,毛烏素腹地站點(diǎn)的降雨主周期大于邊緣站點(diǎn),這可能與這些站點(diǎn)的地理位置有關(guān),這些區(qū)域相對(duì)遠(yuǎn)離黃河、植被蓋度較低且土地利用類型較為單一,局地水汽循環(huán)更新緩慢,從而導(dǎo)致區(qū)域上存在明顯的周期差異。并且賈路等的研究[29]也表明了大氣環(huán)流和地形因素對(duì)水汽輸送過程的影響也是改變區(qū)域降雨周期的重要因素?？偟膩碚f,毛烏素沙地的降雨周期演變差異是復(fù)雜下墊面條件下不同緯度大氣環(huán)流連貫性調(diào)整的結(jié)果[30]。

雖然本研究展示了HHT 在周期識(shí)別和趨勢(shì)檢測(cè)方面的性能,可為水文時(shí)間序列的多尺度周期特征分析以及非線性趨勢(shì)檢驗(yàn)提供一定的參考。但同時(shí),本研究?jī)H以11個(gè)站點(diǎn)39 a的降雨序列為實(shí)例,從應(yīng)用角度對(duì)比分析了HHT 方法和采用不同延拓方式的小波分析在區(qū)域降雨周期演變特征上的差異,對(duì)方法性能的評(píng)價(jià)可能還存在一些局限性。此外,針對(duì)降雨過程在區(qū)域尺度上呈現(xiàn)出的復(fù)雜性和不確定性,可以考慮在之后的研究中將“分解—預(yù)測(cè)—重構(gòu)”的耦合思路用于建立混合預(yù)測(cè)模型,利用分解技術(shù)將數(shù)據(jù)分解為多個(gè)相對(duì)平穩(wěn)的時(shí)間序列,不僅能有效降低降雨的復(fù)雜度和不確定性,而且還能使模型更好地捕捉水文序列的變化特征,之后可以再結(jié)合機(jī)器學(xué)習(xí)或深度學(xué)習(xí)技術(shù),進(jìn)一步提高水文預(yù)測(cè)的精度。目前,HHT 方法在水文領(lǐng)域的應(yīng)用仍處于初步探索階段,尚需在后續(xù)研究工作中進(jìn)一步深化。

4 結(jié)論

本研究基于2組已知成分的人工序列對(duì)HHT方法和小波分析在周期準(zhǔn)確性上進(jìn)行了測(cè)試對(duì)比,并基于毛烏素沙地1982—2020年11個(gè)站點(diǎn)實(shí)測(cè)降雨序列進(jìn)行多尺度周期分析,現(xiàn)得出以下主要結(jié)論:

(1)對(duì)于波動(dòng)規(guī)律較為明顯的人工序列,HHT的周期識(shí)別精度最高。經(jīng)統(tǒng)計(jì)檢驗(yàn)證明,HHT 與小波分析(周期延拓)的周期識(shí)別結(jié)果與預(yù)設(shè)周期一致,不存在顯著差異;而小波分析(對(duì)稱延拓)結(jié)果大于預(yù)設(shè)周期。

(2)對(duì)于實(shí)測(cè)降雨序列,H HT 和小波分析(對(duì)稱延拓)可以較為準(zhǔn)確地識(shí)別降雨周期特征,而小波分析(周期延拓)則有較大偏差;小波延拓方式的選取會(huì)影響實(shí)測(cè)降雨序列的周期識(shí)別結(jié)果。

(3)HHT對(duì)區(qū)域降雨序列多尺度分析有很好的適用性,該方法計(jì)算出的降雨主周期存在明顯的區(qū)域差異:沙地腹地站點(diǎn)的降雨主周期(5.5 a)遠(yuǎn)大于沙地邊緣站點(diǎn)的主周期(2.7～3.3 a);且沙地降雨量在時(shí)間上呈明顯的上升趨勢(shì),在空間上由東到西逐漸減少。

(4)本研究認(rèn)為H HT 相較于小波分析,因其自適應(yīng)性特征,無需考慮眾多參數(shù)的選取,在保持精度的同時(shí)可以提取更多水文信息,在降雨序列多尺度分析上更具優(yōu)勢(shì),因此具有重要的實(shí)際應(yīng)用價(jià)值。