賀冉冉，陳元芳，黃 琴，吳 晶

(1. 河海大學(xué)水文水資源學(xué)院，江蘇 南京 210098； 2.蚌埠學(xué)院材料與化學(xué)工程學(xué)院，安徽 蚌埠 233030)

基于Generalized LASSO模型的長(zhǎng)江寸灘站年均徑流突變研究

賀冉冉1,2，陳元芳1，黃 琴1，吳 晶1

(1. 河海大學(xué)水文水資源學(xué)院，江蘇 南京 210098； 2.蚌埠學(xué)院材料與化學(xué)工程學(xué)院，安徽 蚌埠 233030)

為了區(qū)分水文時(shí)間序列的趨勢(shì)和跳躍變異，將基于L1范數(shù)正則化技術(shù)的generalized LASSO模型應(yīng)用于水文序列變異識(shí)別。經(jīng)過(guò)識(shí)別，發(fā)現(xiàn)長(zhǎng)江寸灘水文站年平均流量序列在1969年發(fā)生了向下的均值躍變。此外，趨勢(shì)分析表明無(wú)論是躍變前后的子序列還是剔除跳躍成分的整個(gè)序列，均未檢出顯著的趨勢(shì)，這說(shuō)明對(duì)寸灘水文站年平均流量序列的跳躍變異假設(shè)是合理的?；趃eneralized LASSO模型的結(jié)果與其他突變檢測(cè)方法結(jié)果進(jìn)行比較，所得結(jié)論是一致的。

水文時(shí)間序列；變異識(shí)別；跳躍點(diǎn)；趨勢(shì)分析；generalized LASSO；長(zhǎng)江上游寸灘水文站

近幾十年來(lái)，由于氣候變化和人類(lèi)活動(dòng)引起下墊面性質(zhì)的改變，導(dǎo)致了徑流序列發(fā)生變異的可能性增大，這也使得對(duì)水文序列的變異分析成為水文時(shí)間序列分析的熱點(diǎn)問(wèn)題之一[1-2]。水文序列變異，是指水文變量的分布形式和分布參數(shù)隨著時(shí)間發(fā)生了變化[2]。在對(duì)水文序列進(jìn)行頻率分析時(shí)，水文序列的非平穩(wěn)性會(huì)導(dǎo)致計(jì)算結(jié)果的不確定性增加[3-4]。對(duì)水文序列的變異進(jìn)行診斷可以使人們更深入地了解水文序列的變化規(guī)律，并能夠合理地對(duì)未來(lái)變化趨勢(shì)進(jìn)行評(píng)估。

在多種類(lèi)型的水文序列變異中，水文變量的均值水平發(fā)生變異是最基本的也是最重要的變異類(lèi)型。水文序列的均值變異有2種可能的表現(xiàn)形式，即趨勢(shì)或者跳躍。從函數(shù)逼近的角度來(lái)看，所謂的趨勢(shì)是指可以用一個(gè)對(duì)時(shí)間的線性函數(shù)來(lái)描述均值水平隨著時(shí)間的變化規(guī)律；而跳躍則是指可以用一個(gè)分段常數(shù)函數(shù)來(lái)描述均值的變化規(guī)律。區(qū)分這2種變異類(lèi)型有著重要的意義，因?yàn)橼厔?shì)相對(duì)于跳躍，更加肯定未來(lái)變異的持續(xù)性[5]。將歷史水文序列診斷為趨勢(shì)或者跳躍變異將會(huì)導(dǎo)致不同的外推結(jié)論。確定水文序列的變異是趨勢(shì)還是跳躍得到一些研究者的重視[6]。

目前有許多水文序列變異識(shí)別方法。針對(duì)趨勢(shì)識(shí)別，最常用的有Mann-kendall趨勢(shì)檢驗(yàn)；而針對(duì)跳躍點(diǎn)進(jìn)行識(shí)別或檢驗(yàn)的方法包括滑動(dòng)t檢驗(yàn)、有序聚類(lèi)法、Pettitt檢驗(yàn)等[7-9]。由于水文序列的復(fù)雜非線性特征和普遍不長(zhǎng)的觀測(cè)年份(一般為幾十年)，變異識(shí)別在應(yīng)用中常將多種識(shí)別檢驗(yàn)方法聯(lián)用，例如通過(guò)綜合的變異診斷系統(tǒng)來(lái)進(jìn)行變異的識(shí)別檢驗(yàn)[10]。多種方法聯(lián)用可以降低結(jié)果的不確定性，這也說(shuō)明開(kāi)發(fā)新的有效的變異點(diǎn)識(shí)別方法是有益的。

LASSO(least absolute shrinkage and selection operator)是近年來(lái)應(yīng)用較為廣泛的變量選擇方法[11]，其原理是對(duì)回歸系數(shù)的L1范數(shù)進(jìn)行約束，可以將一部分系數(shù)收縮為0，進(jìn)而獲得稀疏的模型結(jié)構(gòu)。Tibshirani等[12]提出的generalized LASSO是對(duì)LASSO的拓展，通過(guò)將不同類(lèi)型的懲罰項(xiàng)表達(dá)成統(tǒng)一形式，不僅可以獲得稀疏的模型結(jié)構(gòu)，而且可以獲得特定的模型結(jié)構(gòu)特征。盡管generalized LASSO模型并未在水文時(shí)間序列變異檢驗(yàn)中有過(guò)應(yīng)用，但是與其相關(guān)的fused LASSO模型已經(jīng)在水文時(shí)間序列突變點(diǎn)檢驗(yàn)中有所應(yīng)用[13]。筆者將基于generalized LASSO的識(shí)別方法引入水文序列變異的識(shí)別，其思路是將變異識(shí)別看作是對(duì)水文時(shí)間序列的函數(shù)逼近過(guò)程。將generalized LASSO模型應(yīng)用于函數(shù)逼近時(shí)，通過(guò)對(duì)逼近序列構(gòu)造特定的約束，可以獲得一個(gè)分段常數(shù)逼近序列或者分段一次函數(shù)序列，通過(guò)調(diào)整正則化參數(shù)可以識(shí)別一系列可能的變異點(diǎn)?？梢曰诮徊骝?yàn)證技術(shù)根據(jù)最小的驗(yàn)證集誤差來(lái)選擇最優(yōu)的正則化參數(shù)值，進(jìn)而獲得最優(yōu)的分?jǐn)啾平蛄泻推鋵?duì)應(yīng)的變異點(diǎn)。此外，交叉驗(yàn)證誤差還可以用來(lái)選擇最優(yōu)的假設(shè)，即是用分段常數(shù)序列還是線性函數(shù)序列來(lái)逼近水文序列，進(jìn)而可以確定變異的類(lèi)型是均值跳躍還是趨勢(shì)變異。

本文首先敘述generalized LASSO模型的原理，并介紹其求解過(guò)程和交叉驗(yàn)證的方法，然后以長(zhǎng)江上游寸灘水文站的年平均流量序列為例來(lái)探討基于generalized LASSO的變異識(shí)別方法的應(yīng)用。

1 Generalized LASSO原理

1.1 Generalized LASSO

LASSO模型的實(shí)質(zhì)是在線性回歸最小二乘法的基礎(chǔ)上增加對(duì)系數(shù)絕對(duì)值之和(即系數(shù)向量的L1范數(shù))的約束，用來(lái)對(duì)模型系數(shù)進(jìn)行收縮，進(jìn)而獲得稀疏的模型。LASSO的原理可以參考相關(guān)文獻(xiàn)[11]，本文不再贅述。在LASSO被提出之后，出現(xiàn)了一些類(lèi)似的同樣是基于L1范數(shù)約束的擴(kuò)展模型[14]。在這些模型中，通過(guò)修改模型L1范數(shù)的形式可以獲得一系列特定的模型性質(zhì)。Tibshirani[12]在2011年提出的generalized LASSO模型將多個(gè)不同類(lèi)型的L1范數(shù)正則化問(wèn)題封裝成為同樣的形式。Generalized LASSO模型應(yīng)用于對(duì)水文序列的逼近時(shí)，其形式為

(1)

式中：y——長(zhǎng)度為n年的水文序列，n×1矩陣；β——對(duì)序列y的逼近序列，n×1矩陣；D——某種形式的懲罰矩陣；λ——正則化參數(shù)。在式(1)中，通過(guò)選擇某種類(lèi)型的懲罰矩陣D，可以使Dβ的稀疏性與β的某種性質(zhì)聯(lián)系起來(lái)。

如果假設(shè)對(duì)水文序列逼近函數(shù)的形式為分段常數(shù)函數(shù)(即假設(shè)水文序列的變異類(lèi)型為跳躍型)，則式(1)中的矩陣D可以被設(shè)置為(n-1)×n矩陣D0：

(2)

(3)

式中：βj——β的第j個(gè)元素。也就是說(shuō)，通過(guò)對(duì)β相鄰元素之差的絕對(duì)值之和進(jìn)行懲罰，可以獲得對(duì)原序列y的階梯函數(shù)逼近。在式(3)中，對(duì)逼近序列β跳躍點(diǎn)個(gè)數(shù)的控制通過(guò)設(shè)定λ來(lái)實(shí)現(xiàn)。當(dāng)λ為0時(shí)逼近序列β與原序列相等，而當(dāng)λ無(wú)窮大時(shí)逼近序列β為常數(shù)。在λ從0變?yōu)闊o(wú)窮大的過(guò)程中，β相鄰的元素逐漸融合(即大小變得相等，β序列跳躍點(diǎn)逐漸消失)。這個(gè)模型被稱(chēng)為fused LASSO[14]。

如果假設(shè)逼近函數(shù)的形式為分段線性函數(shù)，則可以設(shè)置矩陣D為(n-2)×n矩陣D1：

(4)

(5)

因此，對(duì)D1β的懲罰就成了對(duì)β序列的離散二階導(dǎo)數(shù)的絕對(duì)值之和進(jìn)行懲罰，而所獲得的逼近序列β就是分段線性函數(shù)。在一些文獻(xiàn)中，這個(gè)過(guò)程被稱(chēng)為趨勢(shì)濾波[15]。

1.2 求解方法

對(duì)(1)式的求解，選取不同的λ值會(huì)有不同的解β，可以通過(guò)所謂求解路徑的方法來(lái)求解所有λ值對(duì)應(yīng)的解。以fused LASSO問(wèn)題為例[12]，求解思路是首先將原問(wèn)題轉(zhuǎn)化為對(duì)偶問(wèn)題進(jìn)行求解：

(6)

式中：u——拉格朗日乘子。

可以證明對(duì)于任意的λ，原問(wèn)題和對(duì)偶問(wèn)題的解存在如下的關(guān)系：

(7)

式中：βλ——正則化參數(shù)取值為λ時(shí)對(duì)應(yīng)的解；uλ——正則化參數(shù)取值為λ時(shí)對(duì)應(yīng)的拉格朗日乘子。

所謂的求解路徑，就是依次求解λ從變化到0的所有值對(duì)應(yīng)的解。由于對(duì)偶問(wèn)題的解向量uλ在λ的變化過(guò)程中是分段線性變化的，也就是λ在相鄰拐點(diǎn)的區(qū)間內(nèi)部變化時(shí)，uλ路徑的方向是不變的，而變化只發(fā)生在若干個(gè)拐點(diǎn)之處。由于式(7)的關(guān)系，可以知道βλ也是分段線性變化的。這樣只需要給出λ在若干關(guān)鍵節(jié)點(diǎn)處對(duì)應(yīng)的解βλ，即完成了整個(gè)求解路徑的求解。矩陣D取值為D0和D1時(shí)求解方法的證明過(guò)程以及計(jì)算程序可以參考相關(guān)文獻(xiàn)[12]。

1.3 交叉驗(yàn)證方法

交叉驗(yàn)證的目的是為了評(píng)估假設(shè)的合理性和選擇最優(yōu)的λ值，其程序簡(jiǎn)述如下[12]。交叉驗(yàn)證首先需要對(duì)樣本進(jìn)行分類(lèi)，其方法是循環(huán)將時(shí)間序列元素依次編入指定的類(lèi)別。具體來(lái)說(shuō)，假設(shè)序列的長(zhǎng)度為8，記為y1,…,y8，進(jìn)行3折交叉驗(yàn)證，將樣本類(lèi)別記為c1、c2、c3，則將樣本依次循環(huán)編入指定類(lèi)別：y1(×),y2(c1),y3(c2),y4(c3),y5(c1),y6(c2),y7(c3),y8(×)，其中兩端括號(hào)中的×表示兩端的元素(即y1和y8)不歸入任何類(lèi)別，始終作為交叉驗(yàn)證時(shí)建模的訓(xùn)練樣本。在進(jìn)行交叉驗(yàn)證時(shí)，如需將c2類(lèi)中的樣本作為測(cè)試集，只需將其余2類(lèi)樣本(也包括兩端的樣本)作為訓(xùn)練集，對(duì)指定的λ計(jì)算逼近序列。任意一個(gè)c2類(lèi)樣本點(diǎn)的預(yù)測(cè)值都可以通過(guò)相鄰的點(diǎn)(因?yàn)榭隙ú辉赾2類(lèi)中)的逼近值取平均求出。這樣，通過(guò)計(jì)算不同λ產(chǎn)生的逼近序列在測(cè)試集上的誤差，可以作為選擇最優(yōu)λ值的參考。筆者在交叉驗(yàn)證時(shí)選擇進(jìn)行10折交叉驗(yàn)證。

圖1 寸灘站年平均流量序列及多年均值水平Fig.1 Annual mean streamflow series at Cuntan station and multi-year mean value

1.4 計(jì)算的軟件實(shí)現(xiàn)

本文的所有計(jì)算和繪圖都是在R語(yǔ)言的環(huán)境下進(jìn)行的，其中對(duì)generalized LASSO模型的求解和交叉驗(yàn)證使用了R語(yǔ)言擴(kuò)展包genlasso[16]。

2 實(shí) 例 分 析

2.1 基于generalized LASSO的寸灘站年平均徑流量序列的變異識(shí)別

以長(zhǎng)江上游重要控制站寸灘站1893—2012年的年平均徑流量序列(圖1)為例，對(duì)其進(jìn)行變異類(lèi)型的確定和變異點(diǎn)的識(shí)別和檢驗(yàn)。從圖1可見(jiàn)，20世紀(jì)60年代以來(lái)，年平均流量總體上處于歷史平均水平之下。

在求解generalized LASSO模型時(shí)，不同的λ值對(duì)應(yīng)著不同的逼近序列。較大的λ對(duì)應(yīng)的逼近序列較為粗糙，而較小的λ對(duì)應(yīng)的逼近序列則較為細(xì)致。圖2分別給出了逼近多項(xiàng)式階數(shù)為0(即分段常數(shù)函數(shù))和階數(shù)為1(即分段線性函數(shù))時(shí)求解路徑中最大的3個(gè)λ所對(duì)應(yīng)的逼近序列，其對(duì)應(yīng)的變異點(diǎn)數(shù)量分別為0個(gè)、1個(gè)和2個(gè)。從圖2可見(jiàn)，隨著正則化參數(shù)的減小，變異點(diǎn)的個(gè)數(shù)在逐漸增加。這些變異點(diǎn)可以作為具有優(yōu)先次序的潛在變異點(diǎn)。

圖2 求解路徑中最大的3個(gè)λ值對(duì)應(yīng)的逼近序列及對(duì)應(yīng)的變異點(diǎn)位置Fig.2 Approximation series and change points corresponding to three maximum values of λ in solution path

如果限制變異點(diǎn)的個(gè)數(shù)為1個(gè)，對(duì)應(yīng)于階數(shù)為0的情況，在1969年出現(xiàn)了一個(gè)向下的跳躍點(diǎn)(圖2(c))；對(duì)應(yīng)于階數(shù)為1的情況，則在1955年出現(xiàn)了趨勢(shì)轉(zhuǎn)折(圖2(d))，在其之后徑流的下降速率明顯增大。當(dāng)減小正則化參數(shù)λ的值，出現(xiàn)2個(gè)變異點(diǎn)時(shí)，在跳躍變異的假設(shè)下會(huì)檢出1927年出現(xiàn)了微弱的向下跳躍(圖2(e))；而在趨勢(shì)轉(zhuǎn)折的假設(shè)下(圖2(f))，在1934年會(huì)出現(xiàn)趨勢(shì)轉(zhuǎn)折，從1934—1955年有短期的上升趨勢(shì)。

考慮到過(guò)多的變異點(diǎn)難以通過(guò)統(tǒng)計(jì)檢驗(yàn)驗(yàn)證，結(jié)合圖2的信息，主要考慮3種備選的變異假設(shè)，分別稱(chēng)為假設(shè)A(全局趨勢(shì)，見(jiàn)圖2(b))、假設(shè)B(1955年為轉(zhuǎn)折點(diǎn)的趨勢(shì)變異，見(jiàn)圖2(d))和假設(shè)C(1969年為跳躍點(diǎn)的均值變異，見(jiàn)圖2(c))。

2.2 通過(guò)交叉驗(yàn)證來(lái)判斷變異類(lèi)型

圖3 不同變異點(diǎn)數(shù)量對(duì)應(yīng)的交叉驗(yàn)證誤差Fig.3 Cross-validation errors corresponding to different numbers of change points

如1.3節(jié)所述，為了確定寸灘年平均流量序列的變異類(lèi)型，需要進(jìn)行交叉驗(yàn)證。筆者同時(shí)給出了多項(xiàng)式階數(shù)為0和1時(shí)，不同變異點(diǎn)個(gè)數(shù)對(duì)應(yīng)的交叉驗(yàn)證誤差(均方誤差)，見(jiàn)圖3。在圖3中，標(biāo)注了2.1節(jié)中所述假設(shè)A、假設(shè)B和假設(shè)C對(duì)應(yīng)的點(diǎn)。考慮到較大的變異點(diǎn)個(gè)數(shù)沒(méi)有實(shí)際意義，在圖3中只顯示變異點(diǎn)個(gè)數(shù)小于等于10的情況。在計(jì)算交叉驗(yàn)證誤差之前，首先將流量序列標(biāo)準(zhǔn)化處理為無(wú)量綱數(shù)，這樣圖3中的交叉驗(yàn)證誤差也是無(wú)量綱的。

從圖3中可以看出，在變異點(diǎn)個(gè)數(shù)為0的情況下，階數(shù)為1相對(duì)于階數(shù)為0來(lái)說(shuō)有著明顯小的交叉驗(yàn)證誤差，說(shuō)明相對(duì)于無(wú)均值變異的假設(shè)來(lái)說(shuō)，全局的線性趨勢(shì)假設(shè)是個(gè)更優(yōu)的假設(shè)。但是從變異點(diǎn)數(shù)量從1到10來(lái)看，0階逼近相對(duì)于1階逼近一直有著較小的誤差。其中在變異點(diǎn)個(gè)數(shù)為1的0階逼近(即假設(shè)C)取得了最小的誤差。這說(shuō)明了存在一個(gè)均值跳躍點(diǎn)的分段常數(shù)函數(shù)逼近是最優(yōu)逼近，也就是說(shuō)，交叉驗(yàn)證的結(jié)果支持一個(gè)跳躍點(diǎn)變異的假設(shè)。

2.3 對(duì)整個(gè)序列及子序列的趨勢(shì)檢驗(yàn)

表1 對(duì)整個(gè)序列及子序列的趨勢(shì)檢驗(yàn)

為了進(jìn)一步考察上述3種假設(shè)的合理性，以上面3種變異類(lèi)型為基礎(chǔ)，分別對(duì)子序列和整個(gè)序列的趨勢(shì)進(jìn)行Mann-Kendall檢驗(yàn)和線性回歸趨勢(shì)檢驗(yàn)，結(jié)果見(jiàn)表1。表1中Z值是Mann-Kendall檢驗(yàn)的標(biāo)準(zhǔn)化統(tǒng)計(jì)量。從表1可以看出，檢驗(yàn)整個(gè)時(shí)間范圍(1893—2012年)時(shí)，2種趨勢(shì)檢驗(yàn)方法都表明有顯著下降的趨勢(shì)(p<0.01)，這樣假設(shè)A就無(wú)法排除。如果檢驗(yàn)假設(shè)B的趨勢(shì)轉(zhuǎn)折變異點(diǎn)前后的子序列，則都沒(méi)有顯著的趨勢(shì)，說(shuō)明假設(shè)B不合理，即序列沒(méi)有發(fā)生趨勢(shì)轉(zhuǎn)折。對(duì)于假設(shè)C，分別檢驗(yàn)跳躍點(diǎn)前后的子序列和剔除躍變成分后的整個(gè)序列，都沒(méi)有發(fā)現(xiàn)顯著的趨勢(shì)存在，這與假設(shè)C不矛盾。這樣，通過(guò)一系列趨勢(shì)檢驗(yàn)，假設(shè)B就被排除，剩下假設(shè)A和假設(shè)C這2種競(jìng)爭(zhēng)性的假設(shè)。

從上面的分析可以看出，趨勢(shì)檢驗(yàn)無(wú)法在2種競(jìng)爭(zhēng)性的假設(shè)(即全局線性下降的趨勢(shì)和一個(gè)跳躍點(diǎn)的變異)之間進(jìn)行選擇，這是因?yàn)槠錄](méi)有能力區(qū)分跳躍和趨勢(shì)這2種變異類(lèi)型[6]，而上述的基于generalized LASSO的交叉驗(yàn)證技術(shù)可以在一定程度上提供另外的證據(jù)。

2.4 與其他跳躍點(diǎn)識(shí)別方法的對(duì)比

為了與其他的跳躍點(diǎn)識(shí)別檢驗(yàn)方法相對(duì)照，筆者進(jìn)行了4種突變點(diǎn)識(shí)別檢驗(yàn)(即有序聚類(lèi)法、滑動(dòng)t檢驗(yàn)法、Pettitt檢驗(yàn)和累積距平法)。在進(jìn)行滑動(dòng)t檢驗(yàn)時(shí)，依次選擇1903—2002年作為分割點(diǎn)，將其前后序列作為2個(gè)樣本進(jìn)行t檢驗(yàn)。結(jié)果表明所有的方法都給出了一致的結(jié)論，即1969年寸灘年平均流量發(fā)生了向下均值躍變。其中滑動(dòng)t檢驗(yàn)和Pettitt檢驗(yàn)表明檢出的變異在0.01水平上是顯著的。這說(shuō)明基于generalized LASSO模型識(shí)別出的跳躍點(diǎn)(假設(shè)C)是合理的。

3 討論與結(jié)論

關(guān)于長(zhǎng)江上游寸灘站的徑流變異特征，已經(jīng)有多個(gè)文獻(xiàn)的研究所涉及。王順久[17]通過(guò)線性趨勢(shì)分析和Mann-Kendall檢驗(yàn)得出寸灘的年徑流存在顯著下降的趨勢(shì)。夏軍等[18]也通過(guò)Mann-Kendall檢驗(yàn)發(fā)現(xiàn)寸灘站年徑流下降的趨勢(shì)，并且通過(guò)累積距平法發(fā)現(xiàn)寸灘站年徑流在1968年之前增加，而在1993年之后減小。根據(jù)李林等[19]的研究，長(zhǎng)江上游河流可以認(rèn)為是雨水補(bǔ)給型河流，降水是影響徑流最主要的因素。馮亞文等[20]指出，長(zhǎng)江上游的年徑流量呈現(xiàn)減少的趨勢(shì)，并且造成這種現(xiàn)象的主要原因仍是降水。夏軍等[18]也指出氣候變化是導(dǎo)致寸灘站徑流發(fā)生了變異的主要因素。根據(jù)相關(guān)的研究及筆者的分析，寸灘站年徑流量確實(shí)在近幾十年是減小的，而難點(diǎn)在于如何評(píng)估其變異類(lèi)型。如前所述，現(xiàn)今主流的趨勢(shì)分析方法并不能區(qū)分全局趨勢(shì)和跳躍變異，所以即使趨勢(shì)檢驗(yàn)的結(jié)果給出了顯著變異的結(jié)論，寸灘年平均流量序列是否存在全局性趨勢(shì)仍然是值得商榷的。本文通過(guò)多種方法的對(duì)比研究，認(rèn)為寸灘站年平均流量序列呈現(xiàn)向下的跳躍變異?？紤]到徑流序列的隨機(jī)性，這個(gè)結(jié)果還存在著一定的不確定性，而且其物理機(jī)制還需要深入研究。

需要指出的是，趨勢(shì)分析和突變分析經(jīng)常被同時(shí)使用。盡管多種方法聯(lián)用可以增加結(jié)果的全面性，但是趨勢(shì)和跳躍之間存在的邏輯關(guān)系卻是一個(gè)需要厘清的問(wèn)題，而這個(gè)問(wèn)題在以往的研究中并未得到重視。如前所述，趨勢(shì)是指均值在整個(gè)時(shí)間段內(nèi)均值緩慢的變化，而跳躍突變是指均值在跳躍點(diǎn)前后從一個(gè)均值狀態(tài)變?yōu)榱硪粋€(gè)狀態(tài)。盡管由于水文序列的長(zhǎng)度較短、噪聲較大導(dǎo)致區(qū)分這2種變異類(lèi)型有難度，但是在邏輯上，這2種結(jié)論是不相容的。以本文對(duì)寸灘流量的分析為例，當(dāng)檢測(cè)出跳躍點(diǎn)時(shí)，即表明全局趨勢(shì)檢驗(yàn)的結(jié)果已經(jīng)受到跳躍點(diǎn)的干擾，所以需要對(duì)跳躍點(diǎn)前后的子序列或者是去除跳躍成分的整個(gè)序列重新進(jìn)行趨勢(shì)檢驗(yàn)。這種先進(jìn)行跳躍點(diǎn)檢驗(yàn)再進(jìn)行趨勢(shì)檢驗(yàn)的方法被一些文獻(xiàn)所強(qiáng)調(diào)[21]。此外，還需要強(qiáng)調(diào)的是，對(duì)跳躍點(diǎn)前后子序列的趨勢(shì)檢驗(yàn)也是對(duì)跳躍點(diǎn)檢驗(yàn)程序前提的驗(yàn)證，例如滑動(dòng)t檢驗(yàn)假設(shè)分割點(diǎn)前后的序列值來(lái)自于2個(gè)隨機(jī)樣本，也就是說(shuō)分割點(diǎn)前后序列的均值應(yīng)該是平穩(wěn)的。

從方法上說(shuō)，基于generalized LASSO模型的水文變異識(shí)別方法的有效性還需要進(jìn)一步探討。需要強(qiáng)調(diào)的是，由于水文序列的長(zhǎng)度較短，并且存在較大的噪聲，根據(jù)最小化交叉驗(yàn)證誤差獲得的逼近序列可能存在過(guò)擬合現(xiàn)象，所以其識(shí)別的變異點(diǎn)應(yīng)該作為“可能的變異點(diǎn)”，最終還需要通過(guò)統(tǒng)計(jì)檢驗(yàn)來(lái)對(duì)其進(jìn)行確認(rèn)。總之，generalized LASSO 模型可作為水文時(shí)間序列可能變異點(diǎn)的識(shí)別程序，其識(shí)別的潛在變異點(diǎn)可以作為后續(xù)統(tǒng)計(jì)檢驗(yàn)的參考。

[ 1 ] MILLY P, JULIO B, MALIN F, et al. Stationarity is dead: whither water management?[J]. Ground Water News & Views, 2007, 4(1): 6-8.

[ 2 ] 謝平, 許斌, 劉媛, 等. 水資源序列變異的時(shí)間尺度主因分析方法[J]. 水力發(fā)電學(xué)報(bào), 2013, 32(6): 24-30. (XIE Ping, XU Bin, LIU Yuan, et al. Key factor analysis of alteration in the hydrological series on time scale[J]. Journal of Hydroelectric Engineering, 2013, 32(6): 24-30.(in Chinese))

[ 3 ] 胡義明, 梁忠民, 趙衛(wèi)民, 等. 基于跳躍性診斷的非一致性水文頻率分析[J]. 人民黃河, 2014, 36(6): 51-53. (HU Yiming, LIANG Zhongmin, ZHAO Weimin, et al. Study on frequency analusis method of non-stationary observation based on jump analysis[J]. Yellow River, 2014, 36(6): 51-53.(in Chinese))

[ 4 ] 謝平, 陳廣才, 夏軍. 變化環(huán)境下非一致性年徑流序列的水文頻率計(jì)算原理[J]. 武漢大學(xué)學(xué)報(bào)(工學(xué)版), 2005, 38(6): 6-9. (XIE Ping, CHEN Guangcai, XIA Jun. Hydrological frequency calculation principle of inconsistent annual runoff series under Changing environments[J]. Engineering Journal of Wuhan University, 2005, 38(6): 6-9.(in Chinese))

[ 5 ] MCCABE G J, WOLOCK D M. A step increase in streamflow in the conterminous United States[J]. Geophysical Research Letters, 2002, 29(24): 38-1-38-4.

[ 6 ]ROUG C, GE Yan, CAI Ximing. Detecting gradual and abrupt changes in hydrological records[J]. Advances in Water Resources, 2013, 53(53): 33-44.

[ 7 ] 周園園, 師長(zhǎng)興, 范小黎, 等. 國(guó)內(nèi)水文序列變異點(diǎn)分析方法及在各流域應(yīng)用研究進(jìn)展[J]. 地理科學(xué)進(jìn)展, 2011, 30(11): 1361-1369. (ZHOU Yuanyuan, SHI Changxing, FAN Xiaoli, et al. Advances in the research methods of abrupt changes of hydrologic sequences and their applications in drainage basins in China[J]. Progress in Geography, 2011, 30(11): 1361-1369.(in Chinese))

[ 8 ] 孫東永,黃強(qiáng), 王義民，等. 滑動(dòng)近似熵在徑流序列突變性分析中的應(yīng)用[J]. 水力發(fā)電學(xué)報(bào), 2014, 33(4): 1-6. (SUN Dongyong, HUANG Qiang, WANG Yimin, et al. Application of moving approximate entropy to mutation analysis of runoff time series[J]. Journal of Hydroelectric Engineering, 2014, 33(4): 1-6.(in Chinese))

[ 9 ] 李彬彬, 謝平, 李析男, 等. 基于Hurst系數(shù)與Bartels檢驗(yàn)的水文變異聯(lián)合分析方法[J]. 應(yīng)用基礎(chǔ)與工程科學(xué)學(xué)報(bào), 2014, 22(3): 481-491. (LI Binbin, XIE Ping, LI Xi’nan, et al. Joint analysis method for hydrological variation based on Hurst coefficient and Bartels test[J]. Journal of Basic Science and Engineering, 2014, 22(3): 481-491.(in Chinese))

[10] 謝平, 陳廣才, 雷紅富, 等. 水文變異診斷系統(tǒng)[J]. 水力發(fā)電學(xué)報(bào), 2010, 29(1): 85-91. (XIE Ping, CHEN Guangcai, LEI Hongfu, et al. Hydrological alteration diagnosis system[J]. Journal of Hydroelectric Engineering, 2010, 29(1): 85-91.(in Chinese))

[11] TIBSHIRANI R. Regression shrinkage and selection via the lasso: a retrospective[J]. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 2011, 73(3): 273-282.

[12] TIBSHIRANI R J, TAYLOR J. The solution path of the generalized lasso[J]. Annals of Statistics, 2011, 39(3): 1335-1371.

[13] JEON J J, SUNG J H, CHUNG E S. Abrupt change point detection of annual maximum precipitation using fused lasso[J]. Journal of Hydrology, 2016, 538：831-841.

[14] TIBSHIRANI R, SAUNDERS M, ROSSET S, et al. Sparsity and smoothness via the fused lasso[J]. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 2005, 67(1): 91-108.

[15] KIM S J, KOH K, BOYD S, et al. L1 Trend Filtering[J]. SIAM Review, 2009, 51(2): 339-360.

[16] ARNOLD T B, TIBSHIRANI R J. Genlasso: path algorithm for generalized lasso problems[EB/OL]. [2014-09-15]. https://CRAN.R-project.org/package=genlasso.

[17] 王順久. 長(zhǎng)江上游川江段氣溫、降水及徑流變化趨勢(shì)分析[J]. 資源科學(xué), 2009, 31(7): 1142-1149. (WANG Shunjiu. Changing pattern of the temperature, precipitation and runoff in Chuanjiang section of the Yangtze River[J]. Resources Science, 2009, 31(7): 1142-1149.(in Chinese))

[18] 夏軍, 王渺林. 長(zhǎng)江上游流域徑流變化與分布式水文模擬[J]. 資源科學(xué), 2008, 30(7): 962-967. (XIA Jun, WANG Miaolin. Runoff changes and distributed hydrologic simulation in the upper reaches of Yangtze River[J]. Resource Science, 2008, 30(7): 962-967.(in Chinese))

[19] 李林, 王振宇, 秦寧生, 等. 長(zhǎng)江上游徑流量變化及其與影響因子關(guān)系分析[J]. 自然資源學(xué)報(bào), 2004, 19(6): 694-700. (LI Lin, WANG Zhenyu, QIN Ningsheng, et al. Analysis of the relationship between runoff amount and its impacting factor in the upper Yangtze River[J]. Journal of Natural Resources, 2004, 19(6): 694-700.(in Chinese))

[20] 馮亞文, 任國(guó)玉, 劉志雨, 等. 長(zhǎng)江上游降水變化及其對(duì)徑流的影響[J]. 資源科學(xué), 2013, 35(6): 1268-1276. (FENG Yawen, REN Guoyu, LIU Zhiyu, et al. rainfall and runoff trends in the upper Yangtze River[J]. Resources Science, 2013, 35(6): 1268-1276.(in Chinese))

[21] LI Dingfang, XIE Huantian, XIONG Lihua. Temporal change analysis based on data characteristics and nonparametric test[J]. Water Resources Management, 2014, 28(1): 227-240.

Abrupt change of annual mean streamflow at Cuntan Station on Yangtze River based on generalized LASSO

HE Ranran1,2, CHEN Yuanfang1, HUANG Qin1, WU Jing1

(1.CollegeofHydrologyandWaterResources,HohaiUniversity,Nanjing210098,China; 2.SchoolofMaterialandChemistryEngineering,BengbuUniversity,Bengbu233030,China)

In order to distinguish trends and abrupt changes in hydrological time series, the generalized LASSO model based on L1-norm regularization was used to detect changes in hydrological series. It is found that the annual mean streamflow series at the Cuntan Station shows a downward step change in 1969. Furthermore, trend analysis shows that both the sub-series before/after 1969 and the whole series after removing step change components show no significant trend, indicating that the assumption of the step change in the annual mean streamflow series at the Cuntan Station is reasonable. The conclusions obtained from the generalized LASSO method are consistent with those obtained from other change detection methods.

hydrological time series; change detection; jump point; trend analysis; generalized LASSO; Cuntan Station in upper reaches of Yangtze River

10.3876/j.issn.1000-1980.2017.04.012

2016-07-23

國(guó)家自然科學(xué)基金(51479061)

賀冉冉(1983—)，男，安徽蚌埠人，講師，博士研究生，主要從事水文不確定性分析。E-mail：heranran2006@163.com

P333

A

1000-1980(2017)04-0358-07