基于年徑流序列的五種趨勢檢測方法性能對比

姜 瑤，徐宗學，王 靜

（1.北京師范大學 水科學研究院，北京 100875；2.城市水循環(huán)與海綿城市技術(shù)北京市重點實驗室，北京 100875；3.西藏自治區(qū)水文水資源調(diào)查局，西藏 拉薩 850000）

1 研究背景

在地球氣候系統(tǒng)顯著變化及人類活動影響下，全球水循環(huán)過程正發(fā)生著不同程度的改變[1-3]，認識和掌握變化環(huán)境下的流域水循環(huán)演變規(guī)律具有重要意義。趨勢變化是水循環(huán)過程的重要變化特性之一，利用趨勢識別技術(shù)理解和認識水文要素在不同時間尺度上的變化規(guī)律是變化環(huán)境下水循環(huán)演變研究的重要內(nèi)容[4-6]。

趨勢識別的主要內(nèi)容是在所關(guān)心的時間尺度之上確定水文時間序列趨勢的具體形狀以及評估趨勢的顯著性[7]。目前，針對全球不同地區(qū)、不同流域的水文氣象要素，國內(nèi)外研究者已開展了大量趨勢識別研究，應(yīng)用或發(fā)展了多種趨勢識別方法[4，8-12]。常用的趨勢識別方法有線性回歸方法、累積距平曲線、滑動t檢驗等參數(shù)統(tǒng)計法，以及Mann-Kendall（MK）秩次檢驗、Spearman 秩次檢驗、Sen’s斜率估計等非參數(shù)統(tǒng)計方法。這些方法多是時間域內(nèi)的趨勢分析[7]，其中線性回歸方法和MK方法在水文氣象序列的趨勢識別中應(yīng)用最為廣泛。線性回歸方法根據(jù)水文觀測數(shù)據(jù)與其對應(yīng)時間之間的回歸系數(shù)判斷水文時間序列的線性趨勢及變化率，然而由于自然系統(tǒng)的非線性、非平穩(wěn)性，水文時間序列在長時間尺度上往往表現(xiàn)出更為復雜的非單調(diào)變化趨勢[13]，因此線性趨勢常常不符合真實的變化過程。MK方法對數(shù)據(jù)的分布形態(tài)沒有特定要求，且不受部分數(shù)據(jù)缺測的影響[14-15]，特別適用于水文氣象數(shù)據(jù)的趨勢檢測，但在實際應(yīng)用中MK分析結(jié)果會受到序列自相關(guān)性、序列長度、周期波動以及趨勢幅度等因素影響，同時也無法給出趨勢變化的具體形狀[16-17]。

為更準確地識別水文時間序列的多時間尺度變化特性，一些時頻分析方法開始廣泛應(yīng)用于水文氣象要素的趨勢分析中，其中離散小波變換（Discrete wavelet transform，DWT）和經(jīng)驗?zāi)B(tài)分解（Em?piric mode decomposition，EMD）方法最為常用。在DWT方法應(yīng)用方面，Pathak 等[18]基于DWT法分析了美國中西部地區(qū)季節(jié)溫度、降水和徑流的變化趨勢；Palizdan 等[19]利用DWT法對Langat 河流域降雨的長期趨勢進行了分析；Sang 等[20]提出了基于DWT法的非單調(diào)趨勢識別及顯著性檢驗方法，并應(yīng)用該方法估計了氣溫與潛在蒸散發(fā)的變化趨勢。DWT法能更有效地識別水文序列的非單調(diào)變化趨勢，但其也受到小波函數(shù)選擇和分解水平選擇等因素的制約[20]。EMD法依據(jù)時間序列本身的特性對序列進行分解，具有直接、后驗和自適應(yīng)等突出特點[21-22]，在水文要素的周期及趨勢識別方面都有著廣泛的應(yīng)用[16，23-26]。

由于趨勢識別方法所基于的原理不盡相同，且不同方法自身均具有局限性，在具體應(yīng)用時往往存在分析能力和分析結(jié)果上的差異[5-6]。為了獲得準確的趨勢識別結(jié)果，一些研究綜合應(yīng)用多種方法進行對比分析，如Tabari 等[10]基于線性回歸、MK檢驗、Pettitt檢驗和Sen’s斜率估計等方法分析研究了伊朗地區(qū)的氣溫、降水及潛在蒸發(fā)量等要素的演變特征；Nalley 等[27]聯(lián)合利用MK方法和DWT方法估計了月尺度、季節(jié)尺度以及年尺度平均徑流量與降雨數(shù)據(jù)的變化趨勢；Kahya 等[28]采用多種非參數(shù)估計方法（MK檢驗、Sen’s斜率、Spearman檢驗和Seasonal Kendall檢驗）對土耳其26個流域的月徑流變化趨勢進行了分析。但目前關(guān)于趨勢的定義不甚統(tǒng)一，對趨勢的認識和識別仍具有局限性，實際研究中很難判斷哪種方法得到的檢驗結(jié)果更為準確和可靠，水文時間序列的趨勢識別仍然存在諸多不確定性[16，29]。

為此，本文選取5種原理不同的趨勢識別方法（線性回歸方法、累積距平方法、MK方法、EMD方法和DWT方法），基于已知成分的人工生成序列和多個流域的實測年徑流序列，綜合對比5種趨勢識別方法的性能差異，探討不同方法對年徑流序列趨勢識別的適用性，以期提高對常用水文要素趨勢識別方法的認識，為水文要素的趨勢分析及其方法選擇提供參考。

2 方法介紹

2.1 線性回歸法線性回歸法是構(gòu)建水文變量x與其時間t之間的一元線性回歸方程：

式中：x（t）為樣本量為n的某一水文變量；t為x（t）所對應(yīng)的時間；a為截距；b為斜率。

a和b可通過最小二乘法求出，其中，回歸系數(shù)b 即可反映水文變量的變化速率，其正負表示水文變量的增減變化趨勢。

回歸系數(shù)b的顯著性可通過t檢驗進行判斷，構(gòu)造T 統(tǒng)計量：

2.2 累積距平法累積距平法是一種由曲線直觀判斷變化趨勢的方法。對于樣本長度為n的水文變量x（t），其某一時刻t的累積距平值cm表示為：該方法的核心是當數(shù)據(jù)持續(xù)大于平均值時，累積距平值增加，曲線呈上升趨勢，反之則呈下降趨勢[6]。根據(jù)曲線的起伏變化，可以判斷時間序列的長期演變趨勢及發(fā)生突變的大致時間。

2.3 MK法Mann-Kendall（MK）檢驗[14-15]是基于原始時間序列秩次的一種非參數(shù)統(tǒng)計方法，其無需數(shù)據(jù)服從一定的分布，同時不受少數(shù)異常值的影響，對于非正態(tài)分布的水文氣象數(shù)據(jù)分析具有突出的適用性，在水文要素的趨勢識別中得到廣泛應(yīng)用。對于一個長度為n的水文時間序列x（t），原假設(shè)H0認為原數(shù)據(jù)序列是n個獨立同分布的隨機變量，沒有趨勢存在。MK檢驗的統(tǒng)計量按下式計算：

其中，

統(tǒng)計量S的正（負）值表示序列的上升（下降）趨勢。當數(shù)據(jù)長度n＞8時，S近似服從標準正態(tài)分布[14-15]，其期望E（S）和Var（S）可表示為：

式中：m為數(shù)據(jù)序列內(nèi)含有的相等數(shù)據(jù)的組數(shù)；kl為第l組內(nèi)相等數(shù)據(jù)的個數(shù)。

標準化的統(tǒng)計檢驗量Z 按下式計算：

2.4 EMD法經(jīng)驗?zāi)B(tài)分解（Empirical model decomposition，EMD）方法的本質(zhì)是設(shè)法識別并分離原始序列固有的多尺度振蕩特征，從而形成一組本征模態(tài)函數(shù)（Intrinsic mode function，IMF）分量和一個殘余項，各IMF表征了原序列不同時間尺度的振蕩變化，殘余項一定程度上表現(xiàn)原序列的總趨勢[21]。EMD法的分解是基于序列時間尺度的局部特性，具有自適應(yīng)性，同時物理意義明確，特別適用于非線性非平穩(wěn)序列的分析處理[23]。EMD 通過自適應(yīng)均值篩選過程實現(xiàn)，詳細描述見文獻[21，23]。經(jīng)過分解后的原序列可以表示為：

式中：N為從原序列x（t）分解得到的IMF總個數(shù)；Cj為第j個IMF分量；rn為殘余項，是一個單調(diào)函數(shù)或無法繼續(xù)分離IMF分量的唯一極值函數(shù)。

各IMF分量的顯著性檢驗可利用Wu 等[23]所建立的基于均勻分布白噪聲統(tǒng)計特征的顯著性檢驗方法。

2.5 DWT法令L（2R）是定義在實軸上、可測的平方可積函數(shù)空間，對于時間序列 f (t)∈L2(R)，其離散小波變換表示為：

其中，

式中：Wf(j,k)為離散小波系數(shù)；ψ*(t)為母小波ψ(t)的復共軛函數(shù)；a、b均為常數(shù)；j為分解水平；k為時間位移因子。

在實際分析中，常使用二進制離散小波變換，即設(shè)a=2和b=1[30-31]，將ψ(j,k)

t表示為：

通常，分解水平j(luò)的理論最大值M可根據(jù)時間序列f（t）的長度L求得[32]：

當選定合理的小波分解水平之后，可應(yīng)用離散小波分解和重構(gòu)將水文序列分解為不同分解水平上的子序列：

各分解水平下子序列的顯著性檢驗參考Sang 等[34]提出的基于噪聲能量分布的顯著性檢驗方法進行。

3 案例分析

3.1 數(shù)據(jù)描述綜合使用人工生成序列和實測徑流序列對比5種趨勢檢測方法的性能差異。3組人工生成序列分別考慮無趨勢、上升趨勢和下降趨勢，并在此基礎(chǔ)上疊加周期項和隨機波動后生成（圖1）。各序列長度均為200，其基本組成見表1。實測序列采用3個不同水文站點的年徑流量序列進行對比分析，分別為雅魯藏布江干流奴下水文站年徑流量實測序列（RN序列）、黑河流域干流鶯落峽水文站年徑流量實測序列（RY序列）和黃河干流上游石嘴山水文站年徑流量實測序列（RH序列），序列長度均為60年（1956—2015年）（表1）。從直觀上看，3組實測數(shù)據(jù)（RN序列、RY序列和RH序列）分別呈不明顯變化、明顯下降和明顯上升趨勢（表1）。為便于對比分析，人工生成序列和實測年徑流序列均采用平均值及標準差進行標準化處理。

圖1 3組人工生成序列實際變化及其真實趨勢

表1 3組人工生成序列與3組實測徑流序列的主要信息

表2 人工生成序列的線性回歸法趨勢檢測結(jié)果

圖2 3組人工生成序列的MK分析結(jié)果

3.2 人工生成序列趨勢識別結(jié)果表2為人工生成序列的線性回歸法趨勢檢測結(jié)果，線性回歸結(jié)果表明，3組人工序列（S1、S2和S3）的平均變化率分別為-0.0046、0.0077和-0.0099，總體上分別呈減小、增加和減小的顯著變化趨勢。

MK法得到3組序列的趨勢檢測結(jié)果見圖2，統(tǒng)計量UFi=200分別為-3.632（下降）、6.439（上升）和-8.327（下降），平均變化率分別為-0.0035、0.0083和-0.0104（圖2），與線性回歸法結(jié)果一致。然而，兩種方法下S1序列的檢驗結(jié)果與其真實趨勢不一致，而S2和S3序列的趨勢檢測結(jié)果雖然正確，但僅是序列的總體變化趨勢，未反映其真實的非線性變化過程。

圖3為3組人工生成序列的累積距平曲線，其呈現(xiàn)了3組人工序列的趨勢變化過程及其拐點。圖3與圖2中MK法所得結(jié)果相似，但由于受序列周期變化的影響，其描述與序列的真實變化均有差異，尤其S1序列，其累積距平值受異常點和周期變化的影響而與實際變化差異較大。

圖3 3組人工生成序列的累積距平曲線

圖4 3組人工生成序列的EMD分析結(jié)果及其顯著性

圖5 3組人工生成序列的DWT分析結(jié)果及其顯著性

應(yīng)用EMD法分別對3組人工序列進行分解，得到其不同時間尺度下的變化特征曲線（即IMF1-5分量和殘余項），其中殘余項的變化代表序列在最大時間尺度下的趨勢。圖4描述了3組序列分解后殘余項與IMF5分量的變化特征及其顯著性檢驗結(jié)果。EMD法識別出S1序列總體表現(xiàn)為減小但趨勢和周期變化均不顯著（圖4（a）（d）），但由于序列隨機波動及個別異常點的影響，所得S1序列變化仍與其真實情況不完全一致。對于S2和S3序列，EMD法所得趨勢項與其真實趨勢基本一致，但識別出的序列周期變化仍存在不平滑的問題，與序列的實際變化存在一定偏差（圖4（b）（c）（d））。

圖5給出了基于DWT方法得到的不同分解水平下的子序列變化及其顯著性檢驗結(jié)果，可分別代表原序列不同尺度下的波動特征。對比EMD法，DWT法檢測出S1序列的趨勢不顯著但周期變化顯著，與S1序列的實際變化基本一致（圖5（a）（b））。同時，DWT方法識別出了S2序列和S3序列的線性變化趨勢及其真實的非線性變化趨勢，且趨勢均顯著，得到的序列變化曲線也與其實際變化更為接近（圖5（b）（c））。

總體上，累積距平法只能描述序列的變化過程及拐點，無法判斷序列的總體變化趨勢，且結(jié)果受周期變化和異常點的影響較大，很難得到準確的變化拐點。線性回歸和MK法能給出序列的總體變化趨勢、平均變化率及其顯著性，對于變化趨勢顯著的序列具有較準確的檢驗結(jié)果，但對周期變化明顯而趨勢不顯著的序列，其識別結(jié)果很容易受周期變化的影響而出現(xiàn)偏差。EMD和DWT方法能對原始序列的不同成分進行有效分解，從而較準確地識別出原序列的周期與趨勢變化，因此更適用于具有復雜變化特征的水文時間序列趨勢分析，其中DWT法的趨勢識別效果較EMD法更優(yōu)。

3.3 實測序列分析分別應(yīng)用上述5種方法對3個流域的實測年徑流序列（RN、RY和RH）進行分析，進一步對比5種趨勢識別方法在徑流變化趨勢分析中的性能差異。3組徑流序列分別具有不同的趨勢和周期變化特征（表1），表3統(tǒng)計了基于5種方法得到的3組實測徑流序列的趨勢檢測結(jié)果。線性回歸法和MK法的結(jié)果均表明，RN、RY和RH序列在60年尺度上的趨勢分別為上升但不顯著、顯著上升和顯著下降（圖6和圖7）。兩種方法計算得到的平均變化速率也相差不大，線性回歸結(jié)果為0.0015/a、0.0268/a和-0.0234/a，MK法為0.0015/a、0.0281/a和-0.0191/a（表3）。同時，由圖7可知，RN序列在1964年和2000年可能存在突變點，這與實測RN序列在1964年和1998—2000年左右的峰值點相對應(yīng)，而RY和RN序列分別在1975年和1978之后變?yōu)轱@著增加趨勢（圖7）。相較于線性回歸法，MK方法定量化地顯示出更多序列變化的細節(jié)。累積距平曲線顯示，RN序列總體有2個明顯的趨勢變化拐點，分別在1964年和1998年，而對于RY和RH序列，其累積距平曲線變化的拐點分別在2000年和1985年左右，該結(jié)果進一步驗證了MK方法所得突變點（圖8）。

表3 3組實測年徑流序列趨勢檢測結(jié)果

圖6 3組實測徑流序列的線性回歸法趨勢分析結(jié)果

圖7 3組實測徑流序列的MK分析結(jié)果

圖8 3組實測徑流序列的累積距平曲線

EMD和DWT方法的趨勢識別結(jié)果與前述3種方法（線性回歸、MK法和累積距平法）所得結(jié)果差異較大，尤其對于RN序列。EMD法識別結(jié)果顯示，在60年尺度上，RN和RH序列呈顯著下降趨勢，RY序列為顯著增加趨勢（表3），同時RN序列在第5分量上（IMF5）表現(xiàn)為顯著非單調(diào)變化（圖9）。與EMD法相似，DWT法將原序列分解為6個水平下的子序列，其中各實測年徑流序列在分解水平6下的子序列變化趨勢與EMD結(jié)果總體一致（圖9和圖10）。然而顯著性分析表明，RN序列的單調(diào)下降（D6）不顯著，在第5分解水平上的子序列為確定成分，表現(xiàn)出非單調(diào)變化趨勢；RY序列的單調(diào)增加趨勢（D6）顯著，為確定的趨勢成分；而RH序列的單調(diào)減小趨勢（D6）和弱非單調(diào)變化趨勢（D5）均顯著（圖10和表3）。

圖9 3組實測徑流序列的EMD分析結(jié)果

圖10 3組實測徑流序列的DWT分析結(jié)果

綜上所述，累積距平法直觀地描述了年徑流量的變化特征，識別出的變化拐點與MK法結(jié)果相似。對于RY和RH序列，由于其趨勢變化非常明顯，因此5種方法中除累積距平法外，均識別出了相似的變化趨勢且趨勢顯著（表3）。然而，對于RN序列，DWT法和EMD法結(jié)果均表明該序列存在顯著的非單調(diào)變化趨勢，受其影響，基于線性回歸法和MK法得到的序列整體表現(xiàn)為增加趨勢但不顯著的結(jié)果，這表明僅線性回歸或MK方法很難準確地識別出RN序列的真實變化趨勢。

3.4 討論基于上述分析，5種趨勢識別方法的特點及適用性如表4所示。從特點上看，累積距平法、線性回歸法和MK法均為時域內(nèi)分析方法，能描述出序列整體隨時間的變化特征，但也易受序列內(nèi)在復雜變化成分（如不同尺度周期變化）的影響。其中，累積距平法僅能反映序列某一區(qū)間的變化特征（上升或者下降），線性回歸法可描述出時間序列總體的線性趨勢，MK法可定量給出序列的總體變化趨勢、拐點及顯著性，但無法確定趨勢的具體形狀。EMD法和DWT法均為時頻域內(nèi)分析方法，可對序列不同時間尺度的變化特征進行有效分解，同時從時域和頻域揭示時間序列的局部特性，從而較準確地識別出序列的真實趨勢成分及其變化的具體形狀（線性和非線性）。

表4 5種趨勢識別方法特點及適用性

從適用性上看，對于變化趨勢明顯的年徑流序列（如RH序列和RY序列），5種方法中除累積距平法外，均得到了較一致的趨勢識別結(jié)果，適用于演變特征相對簡單且趨勢明顯的年徑流序列的趨勢檢測。其中，線性回歸法和累積距平法簡單直觀，可用于年徑流序列變化趨勢的初步分析，MK法與線性回歸法的趨勢識別結(jié)果差異不大，但MK法受序列缺測值、異常值影響較小，且能定量地給出序列趨勢變化的顯著性及突變點，性能優(yōu)于線性回歸法。但MK法與線性回歸法均易受序列周期或非線性特征的影響，不適用于具有明顯周期波動成分的水文時間序列的趨勢識別（如月尺度徑流序列），而EMD法和DWT法均能對原始序列的不同成分進行有效分解，從而較準確地識別出原序列的周期與趨勢變化及其形狀（線性和非線性），其性能明顯優(yōu)于其它3種方法，可用于具有多時間尺度變化、非線性、非平穩(wěn)特性的水文時間序列分析，其中DWT法表現(xiàn)更優(yōu)。例如在案例分析中，由于人工生成序列S1的周期變化影響，其線性回歸和MK法識別結(jié)果均與真實趨勢不符，而EMD和DWT法不僅識別出了S2和S3序列真實的非線性變化趨勢，也較好地識別出了S1序列的趨勢和周期，其中DWT法的結(jié)果與真實變化相符更好，而EMD法仍會受個別異常值的影響而與真實變化產(chǎn)生一定偏差。因此，對于演變機理及趨勢變化復雜（不明顯或非單調(diào)）的年徑流序列（如RN序列），僅用線性回歸或MK方法很難得到準確的趨勢判斷，有必要在序列分解基礎(chǔ)上進一步分析，此時EMD法和DWT法相對更適用。

5種趨勢檢驗方法是水文要素趨勢分析的常用方法，在國內(nèi)外不同流域的降水、徑流、蒸散發(fā)等的變化規(guī)律分析中均有應(yīng)用。研究基于具有不同變化趨勢的人工序列和實測序列，評價了各趨勢檢驗方法對不同趨勢特征的識別能力，所得結(jié)論也可以從以往相關(guān)研究中得到驗證[6，16-17]，可為不同流域水文要素的趨勢分析及其方法選擇提供一定的參考。但同時，由于本文僅以3組年徑流序列為實例，從應(yīng)用角度對比分析了不同方法在水文時間序列趨勢分析中的性能差異，對不同方法的性能評價存在局限性，有待進一步深入研究。

4 結(jié)論

本文選取5種基本原理不同的趨勢識別方法（線性回歸、累積距平、MK、EMD和DWT方法），基于3組已知成分的人工生成序列對比了5種方法在趨勢檢測方面的性能差異，并以3組具有不同趨勢變化特征的實測年徑流序列為例，進一步探討了不同方法在水文時間序列趨勢識別中的適用性。

實例分析結(jié)果表明，累積距平法、線性回歸法和MK法能在時域內(nèi)描述出序列整體的變化特征，但也易受序列內(nèi)在復雜變化成分（如不同尺度周期變化）的影響，適用于趨勢單調(diào)且無明顯周期波動的水文時間序列，對于演變特征相對簡單且趨勢明顯的年徑流序列具有較準確的趨勢檢測結(jié)果。其中，累積距平法僅能表明序列某段變化趨勢，且結(jié)果對序列的異常值較敏感；MK法與線性回歸法的趨勢識別結(jié)果差異不大，但MK法受序列異常值影響較小，且能直觀、定量地給出序列趨勢變化的顯著性及突變點，性能優(yōu)于線性回歸法。EMD法和DWT法可基于時頻分析準確地識別并分離出序列自身含有的不同尺度波動成分，從而避免序列周期、非線性等變化特征對趨勢識別的影響，對于具有多時間尺度變化特征及非平穩(wěn)特性的水文要素序列的趨勢檢測性能更優(yōu)，其中DWT法所得結(jié)果與真實變化更相符。因此，對于演變機理及變化特征復雜的年徑流序列，僅用線性回歸或MK方法很難得到準確的趨勢判斷，有必要在序列成分分解基礎(chǔ)上進行進一步分析，以避免序列趨勢程度、周期、非線性等特征的影響，而其中DWT方法應(yīng)為首選方法，EMD方法則是另一種較為可行和值得推薦的方法。