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

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

(1.長安大學水利與環境學院,西安 710054;2.旱區地下水與生態效應教育部重點實驗室,西安 710054;3.水利部旱區生態水文與水安全重點實驗室,西安 710054)

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

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

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

1 方法與數據

1.1 HHT的基本理論及特性

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

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

非線性的k個空間數據的MEMD具體步驟如下[15]:

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

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

(3)識別每一個方向向量上pθk(t)極值出現的位置tθki;

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

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

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

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

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

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

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

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

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

對于時間序列f(t)∈L2(R),其連續小波變化為:

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

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

1.3 配對樣本t檢驗

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

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

1.4 數據資料

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

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

圖2 人工序列數據示意圖Fig.2 Schematic diagram of manual sequence data

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

如圖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)以及伊金霍洛旗站點(S11)。

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

2 結果與分析

2.1 人工序列周期識別結果

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

表2 基于人工序列使用HHT得到的各固有模態函數周期表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分解得到的固有模態函數(IMF1-7)和殘差項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分解得到的固有模態函數(IMF1-6)和殘差項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分別為人工序列采用不同延拓方式(對稱和周期延拓)獲取到的小波方差以及對應的小波實部。小波方差反映了時間序列的波動能量隨時間尺度的分布,根據峰值的個數可以確定時間序列中存在的主周期數量;小波實部則能反映出不同特征時間尺度信號在不同時間上的分布和相位等信息。結果表明:(1)當采用對稱延拓的小波分析時,M1,M2和M3的小波方差最大值均于19處取得,識別出的周期均為12.6 a;M4和M5的小波方差最大值分別于23,15處取得,識別出的周期為16 a和10.4 a,對周期的識別總體上相對偏大(表1)。(2)相較于前者,采用周期延拓的小波分析識別出的M1-M5周期與人工序列的預設周期基本一致,分別為12.2 a,12.2 a,12.2 a,14.9 a以及10.1 a。

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

為進一步檢驗各方法的周期結果與人工序列預設周期的精度差異,進行了配對樣本t檢驗。統計結果表明:HHT 的顯著性為Sig.=0.142＞0.05,采用周期延拓的小波分析的顯著性為Sig.=0.109＞0.05,兩者的周期識別結果與預設周期一致,不存在顯著差異,且HHT 的周期準確性要略高于周期延拓的小波分析;相比之下,采用對稱延拓的小波分析的顯著性為Sig.=0.003＜0.05,識別出的周期結果偏大,與預設周期存在顯著差異(表3)?？偟膩碚f,H HT 與小波分析(周期延拓和對稱延拓)相比,周期識別精度最高。

表3 配對樣本t檢驗結果Table 3 Paired sample t test results

2.2 實測序列分析

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

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

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

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

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

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

3 討論

本研究利用人工序列和實測降雨序列,對H HT和小波分析在求解區域尺度降雨周期演變特征的適用性上進行了測試、對比和分析。人工序列的結果可為周期識別的準確性提供依據,而實測降雨數據為研究區域多尺度的周期演變提供應用實例,兩者結果相輔相成。

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

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

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

4 結論

本研究基于2組已知成分的人工序列對HHT方法和小波分析在周期準確性上進行了測試對比,并基于毛烏素沙地1982—2020年11個站點實測降雨序列進行多尺度周期分析,現得出以下主要結論:

(1)對于波動規律較為明顯的人工序列,HHT的周期識別精度最高。經統計檢驗證明,HHT 與小波分析(周期延拓)的周期識別結果與預設周期一致,不存在顯著差異;而小波分析(對稱延拓)結果大于預設周期。

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

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

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