一種基于ITA改進(jìn)的水文氣象序列趨勢分析法

吳國棟，劉廷璽，薛河儒

(1.內(nèi)蒙古農(nóng)業(yè)大學(xué)計算機(jī)與信息工程學(xué)院，內(nèi)蒙古 呼和浩特 010021；2.內(nèi)蒙古農(nóng)業(yè)大學(xué)理學(xué)院，內(nèi)蒙古 呼和浩特 010018； 3.內(nèi)蒙古農(nóng)業(yè)大學(xué)水利與土木建筑工程學(xué)院，內(nèi)蒙古 呼和浩特 010018)

由于氣候的變化，在水資源規(guī)劃、管理等研究中，對氣溫、徑流、降水、蒸散發(fā)等數(shù)據(jù)序列在不同時間尺度(如年度、月度、季度等)的趨勢識別、檢測和評價成為近年來重要的研究內(nèi)容[1-2]。目前，國內(nèi)外學(xué)者就全球不同地區(qū)的水文氣象要素開展了大量的趨勢分析研究，基于不同的理論提出和發(fā)展了多種趨勢分析方法[3-6]?？傮w而言，針對水文氣象要素的時間序列在時域上的趨勢分析方法包括兩大類：參數(shù)型趨勢檢測法和非參數(shù)型趨勢檢測法。參數(shù)型趨勢檢測方法主要有線性回歸、滑動平均、累積距平等方法，而Sen斜率估計、Mann-Kendall秩次檢驗等屬于非參數(shù)型趨勢檢測方法。

近30年來，相關(guān)研究人員采用Mann-Kendall秩次檢驗法(MK法)開展了很多關(guān)于環(huán)境大氣、水文、氣象和農(nóng)業(yè)等方面的研究。MK法是一種應(yīng)用廣泛的非參數(shù)型趨勢檢測法[7-10]，不僅可以探測時間序列所蘊含趨勢的增減性，還可以從統(tǒng)計學(xué)角度分析趨勢的顯著性水平，具有簡單易操作、結(jié)果可信度高的特點，而且不受限于數(shù)據(jù)的概率分布。2012年，Sen[11]提出一種改進(jìn)的非參數(shù)型趨勢檢測法，即ITA(innovative trend analysis)法。ITA法是從圖形上直觀地觀察數(shù)據(jù)序列的趨勢。由于ITA的直觀性和對數(shù)據(jù)獨立性及分布規(guī)律無要求，所以受到廣泛關(guān)注和應(yīng)用[9，11-15]。ITA法不僅可以從整體上觀察數(shù)據(jù)序列，還可以分開研究數(shù)據(jù)在低值區(qū)、中值區(qū)和高值區(qū)等不同數(shù)據(jù)區(qū)域的變化趨勢[12]。Sen[16]在2017年補充了ITA法的趨勢指標(biāo)及其顯著性檢驗方法，使得ITA法的功能更加完善。

雖然ITA法是一種易用、直觀的趨勢檢測方法，但由于其為參數(shù)型檢測方法，這使得其結(jié)果會受到數(shù)據(jù)序列分布特征的影響。本文基于ITA法設(shè)計一種新的趨勢統(tǒng)計量，采用非參數(shù)統(tǒng)計方法——自舉法[17-18](Bootstrap)進(jìn)行顯著性檢驗。采用數(shù)值模擬方式分析多種人工數(shù)據(jù)序列，將改進(jìn)的趨勢檢測法分別與經(jīng)典的MK法和ITA法進(jìn)行檢測結(jié)果的一致性比較，驗證其可行性。最后應(yīng)用改進(jìn)的趨勢檢測法對黑河上游年徑流總量、日本福岡每年風(fēng)暴天數(shù)、北京最大日降水量和瓊海年均氣溫等4種時間序列數(shù)據(jù)進(jìn)行趨勢檢測和顯著性分析。

1 研究方法

圖1 ITA方法中的5種趨勢類型Fig.1 Five kinds of trend according to ITA

ITA法是一種較為直觀的趨勢檢測方法，通過它可以觀察數(shù)據(jù)序列單調(diào)或非單調(diào)增加或減少趨勢，共有5種情況：單調(diào)增加、非單調(diào)增加、單調(diào)減少、非單調(diào)減少和無趨勢。ITA法的具體過程為：首先將數(shù)據(jù)序列平均分成相等的兩部分，并分別按照升序排列，記為{y1i}和{y2i}(i=1，2，…，n/2)。然后在直角坐標(biāo)系內(nèi)以{y1i}為橫坐標(biāo)，{y2i}為縱坐標(biāo)畫出散點圖并與無趨勢直線(1∶1)進(jìn)行比較(圖1)。如果散點分布在45°線的上方，說明存在單調(diào)增加的趨勢；如果散點分布在45°線的下方，說明存在單調(diào)減少的趨勢；如果散點分布于45°線的附近，說明無趨勢。而散點越遠(yuǎn)離45°線，說明趨勢越明顯。

2017年，Sen[13]又提出基于ITA法的趨勢顯著性檢驗方法，其統(tǒng)計量為

(1)

雖然ITA法可以在圖形上直觀地判斷趨勢的增減性和單調(diào)性，甚至可以觀測在不同大小區(qū)間內(nèi)數(shù)據(jù)序列的趨勢，但是該方法在趨勢的顯著性檢驗中回到了參數(shù)性檢驗過程，趨勢顯著性結(jié)果必然會受到數(shù)據(jù)序列概率分布和數(shù)字特征的影響。而水文氣象要素的變化往往存在較大的隨機(jī)性，所以相應(yīng)的時間序列服從何種概率分布難以確定，尤其在序列長度較小時。另外，水文氣象要素時序數(shù)據(jù)記錄的不可重復(fù)性，也增加了確定其概率分布的難度。Sen[13]提出ITA方法時，并沒有給出非常嚴(yán)格的數(shù)學(xué)證明來說明按照標(biāo)準(zhǔn)正態(tài)分布對ZITA進(jìn)行檢驗的合理性，因此考慮將ITA中的趨勢顯著性檢驗過程由參數(shù)型檢驗改為效果等同或更好的非參數(shù)型檢驗過程。

自舉法是非參數(shù)統(tǒng)計應(yīng)用中對統(tǒng)計量進(jìn)行區(qū)間估計的重要方法之一，是現(xiàn)代統(tǒng)計學(xué)中依托計算機(jī)高速發(fā)展并被廣泛應(yīng)用的一種重抽樣技術(shù)。其基本原理是在一個樣本中進(jìn)行多次有放回的均勻采樣，生成一系列能夠代表研究總體的多個樣本，然后再根據(jù)抽出的樣本計算統(tǒng)計量。按此過程可以生成關(guān)于待檢驗統(tǒng)計量的一個經(jīng)驗分布，便可估計其置信區(qū)間來進(jìn)行統(tǒng)計推斷。自舉法的優(yōu)點是對樣本的概率分布沒有任何要求，且無需增加更多的樣本信息就可以對統(tǒng)計量進(jìn)行準(zhǔn)確的統(tǒng)計推斷，這與水文氣象要素序列的隨機(jī)性和數(shù)據(jù)的不可重復(fù)性完全契合。

另一方面，根據(jù)ITA法的原理，散點分布越遠(yuǎn)離無趨勢直線，趨勢越明顯。為了更好地說明散點偏離無趨勢直線的程度，采用相對距離作為衡量趨勢及顯著性的新指標(biāo)。

綜上，在ITA法的基礎(chǔ)上，引入一個新指標(biāo)P檢測趨勢的增減性，并用自舉法確定其顯著性：

(2)

具體過程描述如下：

a.將時間序列數(shù)據(jù){xi}(i=1，2，…，n)平分為兩段，按升序排列后記為{y1i}和{y2i}，并計算P。

b.在{xi}(i=1，2，…，n)中進(jìn)行有放回的重采樣，即從序列{xi}中依次等概率地隨機(jī)抽取n個數(shù)據(jù)(其中有的數(shù)據(jù)可能會被多次抽到)構(gòu)成一個新時間序列樣本，詳細(xì)計算過程可參考文獻(xiàn)[18]。重復(fù)采樣M次后生成M個長度為n的新序列，對每個新序列重復(fù)步驟a，并將得到的新指標(biāo)按照升序排列：P1，P2，…，PM。

c.考慮置信水平α，計算

L=Mα/2H=M(1-α/2)

(3)

則置信區(qū)間為[PL，PH]。如，M=1 000，α=0.05，則L=25，H=975。從升序排列P1，P2，…，P1 000中找到第25和第975位置的兩個數(shù)P25和P975，則原始序列{xi}的新指標(biāo)P在置信水平為0.05上的置信區(qū)間為[P25，P975]。

d.如果{xi}(i=1，2，3，…,n)的P落入置信區(qū)間，則認(rèn)為在α置信水平上趨勢不具有顯著性，反之，認(rèn)為有顯著性，且P>0代表增趨勢，P<0代表減趨勢。

自舉法是在未知序列分布規(guī)律的情況下進(jìn)行顯著性檢驗的有效辦法，Mailhot等[19]指出，M=1 000可以保證結(jié)果的可靠性，本文選擇M=5 000進(jìn)行計算。

2 數(shù)值模擬驗證

為了驗證改進(jìn)的趨勢檢測法的可行性，人工生成一階自回歸隨機(jī)過程數(shù)據(jù)序列進(jìn)行蒙特卡羅模擬。按式(4)生成長度為n的數(shù)據(jù)序列：

(4)

式中：μ為平均值；ρ為一階序列相關(guān)系數(shù)；εi為服從標(biāo)準(zhǔn)正態(tài)分布的隨機(jī)擾動。然后根據(jù){Xi+d·i}(i=1，2，…，n)嵌入斜率為d的線性趨勢。

利用改進(jìn)的趨勢檢測法與MK法(具體過程可參考文獻(xiàn)[7]或[9])和ITA法分別對人造序列進(jìn)行趨勢和顯著性檢驗，并記錄結(jié)果。在相同的參數(shù)設(shè)定下，重復(fù)多次后統(tǒng)計改進(jìn)的趨勢檢測法與MK法和ITA法檢測結(jié)果的一致率。如，當(dāng)μ=20，σ=5，按序列的相關(guān)系數(shù)ρ=0、±0.1、±0.5、±0.9以及植入的線性趨勢斜率d=0.01、±0.05、±0.09組合生成共35種參數(shù)數(shù)據(jù)序列，每種參數(shù)序列各1 000條。對該1 000條數(shù)據(jù)序列分別進(jìn)行MK法、ITA法和改進(jìn)的趨勢檢測法的趨勢檢驗，統(tǒng)計改進(jìn)的趨勢檢測法分別與MK法和ITA法在趨勢類型(增或減)、10%置信水平的趨勢顯著性檢驗結(jié)果的一致率，以及綜合兩方面檢驗結(jié)果的一致率，結(jié)果見表1。

檢驗結(jié)果表明：

a.改進(jìn)的趨勢檢測法與MK法和ITA法一致率不受序列均值的影響。

b.序列中蘊含的趨勢越強(d的絕對值越大)，3種方法的檢測結(jié)果一致程度越大。

c.在趨勢類型一致率方面，改進(jìn)的趨勢檢測法與MK法的結(jié)果最大誤差小于5%，與ITA法的結(jié)果最大誤差小于1%。

d.在顯著性一致率方面，只有當(dāng)趨勢較弱(d絕對值很小)且前后兩個子序列存在強正相關(guān)時，出現(xiàn)一致率下降的現(xiàn)象。如，當(dāng)d=±0.01、ρ=0.9時，一致率下降到90%以下。

在數(shù)值模擬過程中發(fā)現(xiàn)，ITA法的檢驗指標(biāo)會受到數(shù)據(jù)序列的長度和標(biāo)準(zhǔn)差的影響，而改進(jìn)的趨勢檢測法受其影響較小。總體而言，數(shù)據(jù)序列越長，改進(jìn)的趨勢檢測法和ITA法的趨勢顯著性檢驗結(jié)果與MK法的結(jié)果越一致，且標(biāo)準(zhǔn)差越大時越需要更長的數(shù)據(jù)序列才能保證較高的一致率。相比而言，用改進(jìn)的趨勢檢測法檢驗同樣標(biāo)準(zhǔn)差的數(shù)據(jù)序列，達(dá)到指定的一致率所需的數(shù)據(jù)序列長度較小。換而言之，如果數(shù)據(jù)序列長度較小時(如小于100)，改進(jìn)的趨勢檢測法的檢驗效果要優(yōu)于ITA法。

表1 利用改進(jìn)的趨勢檢測法與MK法和ITA法模擬結(jié)果比較

如當(dāng)生成1 000條ρ=0、μ=50、σ=30、n=100的序列，并加入0.05的趨勢斜率，分別用MK法和ITA法進(jìn)行趨勢檢驗。結(jié)果表明兩種方法的趨勢類型檢驗結(jié)果一致率為84.3%，而顯著性檢驗結(jié)果的一致率只有25.9%。而針對同樣的數(shù)據(jù)序列，改進(jìn)的趨勢檢測法與MK法的趨勢類型檢驗結(jié)果一致率可達(dá)到90.2%。圖2反映了模擬過程中的1個數(shù)據(jù)序列，經(jīng)MK法和改進(jìn)的趨勢檢測法檢驗無顯著趨勢，但是ITA法的檢驗結(jié)果是在0.01的置信水平上顯著負(fù)趨勢。究其原因，可能是ITA法中的s與時間序列的平均值有關(guān)，當(dāng)樣本容量不大時均值易受數(shù)據(jù)序列中異常值的影響，造成均值不能很充分地代表樣本的中心趨勢。另外數(shù)據(jù)序列的標(biāo)準(zhǔn)差增大后，不僅加劇了異常值對均值的影響，而且也會影響ITA法顯著性檢驗中的置信區(qū)間。事實上，根據(jù)該數(shù)據(jù)序列的正態(tài)分布Q-Q圖可以確定其為非正態(tài)分布(圖2(c))。由此可見，ITA法會受到時間序列分布特征和數(shù)據(jù)統(tǒng)計特征的影響而產(chǎn)生錯誤結(jié)果。但是改進(jìn)的趨勢檢測法的趨勢顯著性檢驗過程完全是非參數(shù)化過程，顯著性不易受到來源于數(shù)據(jù)序列的分布特征和統(tǒng)計特征的影響，所以檢驗結(jié)果與經(jīng)典的MK法一致程度高。

圖2 一列ITA法與MK法結(jié)果差異大的數(shù)據(jù)Fig.2 A series with big differences between ITA and MK results

3 實例應(yīng)用

將改進(jìn)的趨勢檢測法應(yīng)用于不同地區(qū)、不同長度、不同水文氣象要素的數(shù)據(jù)序列?？紤]數(shù)據(jù)的可獲得性、序列長度、數(shù)據(jù)的完整性和水文氣象要素的種類，對黑河上游(1891—2008年)的年徑流量、日本福岡站(1948—2019年)年風(fēng)暴天數(shù)、北京氣象站(1952—2019年)年最大日降水量和瓊海氣象站(1954—2019年)年平均氣溫共4種時間序列進(jìn)行分析。黑河上游的年徑流量數(shù)據(jù)來源于國家冰川凍土沙漠科學(xué)數(shù)據(jù)中心 (http：//www.ncdc.ac.cn)；日本福岡站的數(shù)據(jù)來源于Tutiempo Network(https：//en.tutiempo.net)，個別缺失數(shù)據(jù)經(jīng)3次樣條插值做了補充；北京氣象站和瓊海氣象站的數(shù)據(jù)來源于國家氣象科學(xué)數(shù)據(jù)中心(http：//data.cma.cn)。4個時間序列的基本信息和統(tǒng)計特征信息見表2、表3。

表2 時間序列基本信息

表3 時間序列統(tǒng)計特征信息

為了直觀反映各數(shù)據(jù)序列的平均趨勢，將每個數(shù)據(jù)序列及其線性回歸直線合并成圖3，同時將檢測結(jié)果匯總，見表4。

圖3 4個時間序列的數(shù)據(jù)分布及回歸斜率Fig.3 Data distribution and regression slope of four time series

表4 4個時間序列的趨勢檢測結(jié)果

黑河上游年徑流量在0.05的顯著性水平上具有顯著的增加趨勢，解陽陽等[20-21]應(yīng)用MK法也得出了相同的結(jié)論，并指出黑河上游徑流量增加的主要原因是降水量和平均氣溫的增加。Zhang等[22]研究了近50年北京市暴雨趨勢及其空間差異，采用MK法檢驗得到的結(jié)論表明，北京市年最大日降雨量只在0.1的顯著性水平上呈現(xiàn)顯著下降的趨勢，該結(jié)論與改進(jìn)的趨勢檢測法的檢驗結(jié)果完全一致。孫瑞等[23]同樣采用MK法研究了1959—2013年海南島氣象要素的時空分布特征，結(jié)果顯示海南島平均氣溫的年際增長趨勢通過了0.01顯著性水平檢驗，這與本文檢驗的瓊海氣象站年均氣溫記錄在0.01顯著性水平上具有顯著增加趨勢的結(jié)論一致。日本福岡站的年風(fēng)暴天數(shù)時間序列也表現(xiàn)出與瓊海氣象站年均氣溫相同顯著性水平的增加趨勢(表4)。

4 結(jié) 論

a.改進(jìn)的趨勢檢測法是一種非參數(shù)型趨勢檢測方法，檢測結(jié)果與經(jīng)典的MK法保持高一致率，表明了改進(jìn)的趨勢檢測法的可行性。

b.改進(jìn)的趨勢檢測法不受限于數(shù)據(jù)序列的分布特征，在數(shù)據(jù)序列方差較大或數(shù)據(jù)系列長度較小時仍能保證與MK法良好的一致性，優(yōu)于ITA法。

c.在數(shù)據(jù)序列蘊含的趨勢斜率很小且正相關(guān)系數(shù)很大時，檢測結(jié)果的一致性會受到一定的影響，原因需進(jìn)一步研究。

d.經(jīng)實際應(yīng)用，黑河上游的年徑流量呈高顯著性水平(5%)的增加趨勢，日本福岡每年發(fā)生風(fēng)暴的天數(shù)呈現(xiàn)高顯著性水平的增加趨勢，北京的最大日降水量是中等顯著性水平(10%)的下降趨勢，瓊海的年平均氣溫表現(xiàn)為高顯著性水平的上升趨勢。