吳國(guó)棟，劉廷璽，薛河儒

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

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

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

1 研究方法

圖1 ITA方法中的5種趨勢(shì)類(lèi)型Fig.1 Five kinds of trend according to ITA

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

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

(1)

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

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

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

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

(2)

具體過(guò)程描述如下：

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

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

c.考慮置信水平α，計(jì)算

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位置的兩個(gè)數(shù)P25和P975，則原始序列{xi}的新指標(biāo)P在置信水平為0.05上的置信區(qū)間為[P25，P975]。

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

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

2 數(shù)值模擬驗(yàn)證

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

(4)

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

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

檢驗(yàn)結(jié)果表明：

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

b.序列中蘊(yùn)含的趨勢(shì)越強(qiáng)(d的絕對(duì)值越大)，3種方法的檢測(cè)結(jié)果一致程度越大。

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

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

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

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

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

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

3 實(shí)例應(yīng)用

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

表2 時(shí)間序列基本信息

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

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

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

表4 4個(gè)時(shí)間序列的趨勢(shì)檢測(cè)結(jié)果

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

4 結(jié) 論

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

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

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

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