基于核熵成分分析的綜合干旱指數的構建與應用—以黑河流域中上游為例

郭盛明，粟曉玲,2，吳海江，姜田亮，梁 箏，馮 凱

(1.西北農林科技大學水利與建筑工程學院，陜西 楊凌 712100；2. 西北農林科技大學旱區農業水土工程教育部重點實驗室，陜西 楊凌 712100)

旱災發生頻率高、影響范圍廣、災害損失大，被認為是世界上最嚴重的自然災害之一[1-2]。隨著全球氣候變化和社會經濟的高速發展，人類對自然環境的干預越來越大，降低了氣候系統的穩定性，加速了降水、蒸發等水循環過程，使得洪旱災害的發生風險進一步增大[3-4]。

干旱通常被分為氣象干旱、農業干旱、水文干旱、社會經濟干旱4種類型。常用影響干旱的氣象、水文要素構建干旱指數監測干旱。氣象干旱是指降水量與蒸發量之間長時間不平衡而出現的水分短缺現象[5]，常用帕爾默干旱指數(Palmer drought severity index，PDSI)、標準化降水指數(Standardized precipitation index，SPI)、標準化降水蒸散發指數(Standardized precipitation evapotranspiration index，SPEI)等表征[6]。農業干旱主要指土壤水分在一段時間內持續虧缺造成農作物不能正常生長發育的現象[7]，常用K指數、降水距平百分率、標準化土壤濕度指數(Standardized soil moisture index，SSMI)等表征[8]。水文干旱是指河川徑流滿足不了供水需求，使得可利用水量低于正常水平的現象[9]，常用徑流距平指數、地表供水指數、標準徑流指數(Standardized runoff index，SRI)等表征[10]。

干旱的發展及傳遞過程復雜，使用單變量干旱指數無法描述干旱的全部特征，因此，需要探究一種融合多變量的方法來構建綜合干旱指數[11]。常用的構建方法主要有水量平衡法、權重法和概率統計法[12]。其中，權重法包括熵權法、主成分分析法(PCA)等。王璽圳等[13]采用客觀賦權法聯合SPEI、NDVI、SSI和SGI構建綜合干旱指數CDI；Wang 等[14]基于PCA考慮了徑流變化提出了一種新的標準化降水蒸散徑流指數SPERI。概率統計法常用Copula函數，Hao等[15]基于Copula 函數聯合SSI和SPI提出了多變量標準化干旱指數MSDI；張迎等[16]采用Copula 函數聯合SPI和SRI構建綜合干旱指數MSDIp。熵權法、PCA等方法將若干變量線性組合在一起，可有效提取特征間的線性關系，但不能很好地反映變量間的非線性關系。Copula函數能連接多個不同分布的變量，但計算方法繁雜，尤其是連接3個及以上變量時。核熵成分分析(Kernel entropy component analysis，KECA)是對PCA的深層次改進，其在特征空間進行熵成分分析以實現數據變換，具有良好的非線性處理能力且計算相對簡單[17]。目前，KECA已被廣泛應用于目標分類、故障檢測等領域[18-20]，但在干旱領域鮮有應用。本文基于SPEI、SSMI和SRI，采用KECA構建綜合干旱指數SMDI，并探討其時空演變特征和適用性，以期為干旱監測提供科學指導。

1 研究區概況和數據來源

1.1 研究區概況

黑河是我國第二大內陸河，位于河西走廊中部，東以山丹縣境內的大黃山為源，西以嘉峪關境內的黑山為界，地理坐標為97°50′～102°00′E，37°50′～42°40′N[21]。黑河干流發源于祁連山系的走廊南山和冷龍嶺，全長821km，由鶯落峽和正義峽兩斷面將其劃分為上、中、下游。上游為主要產流區，植被較好，為牧業區；中游為徑流利用區，也是重要農業區，主要依靠黑河供水，干旱缺水；下游植被稀少，為徑流消耗區[22]。

1.2 數據來源

氣象數據為黑河流域中上游12個氣象站1966—2016年逐月降水、風速等氣象資料，來源于中國氣象科學數據共享服務網(http://data.cma.cn/)。水文數據為黑河干流鶯落峽和正義峽兩個水文站1966—2016年的逐月徑流資料，來源于國家自然科學基金委員會“黑河數據計劃管理中心”。黑河流域中上游1966—2016年逐月土壤濕度數據來源于全球陸地數據同化系統(GLDAS)數據集，空間分辨率為0.25°×0.25°。通過ArcGIS將土壤濕度數據重采樣為0.125°×0.125°柵格數據集，與氣象、水文聯合進行黑河流域中上游柵格計算。流域氣象站及水文站分布見圖1。

2 研究方法

2.1 單變量干旱指數的研究方法

用SPEI、SSMI和SRI3種指數分別表征氣象干旱、農業干旱和水文干旱。SPEI的計算方法詳見文獻[23]，其中潛在蒸散發量(ET0)采用Penman-Monteith公式計算[24]。

土壤濕度的邊緣累積概率P(xi)使用非參數的Gringorten[25]經驗位置劃分公式來計算，公式如下：

P(xi)=(i-0.44)/(n+0.12)

(1)

式中，i為土壤濕度序列升序時的次序，n為樣本序列長度。

將經驗累積概率P(xi)進行逆標準化即得到非參數化的SSMI，公式如下：

SSMI=φ-1(P)

(2)

式中,φ為標準正態分布函數，P為經驗累積概率。

SRI的計算方法與SPI類似，但與SPI假定降水服從Gamma分布不同，SRI假定徑流服從對數正態分布。在求得徑流的累積概率之后將其逆標準化即得到SRI。

2.2 核熵成分分析(KECA)基本原理

KECA通過構造數據的Renyi熵，以高維特征空間能最大程度保持原始空間Renyi熵的坐標軸為投影方向[26]。

給定N維樣本x，其概率密度函數為p(x)，則p(x)的Renyi 熵可表示為：

(3)

由于對數函數具有單調性，所以只需關注真數部分，真數部分用V(p)表示：

(4)

引入Parzen窗函數估計V(p)值，表達式如下：

(5)

式中，K(xi,xi)為Parzen窗函數(核函數)，xi為中心，為寬度參數。

可以獲得V(p)的估計值：

(6)

式中，K是一個N階的核矩陣，K(i,j)=K(xj,xi)；I是一個N×1向量(向量中的每個元素均為1)。

樣本集的Renyi熵可通過核函數進行估計。將核矩陣K進行特征值分解：K=EDET，其中D=diag(λ1,λ2,…，λN)是特征值所構成的對角線矩陣，E=[α1,α2,…,αN]是由相應的特征向量組成的特征矩陣，λ和α為Renyi熵的估計值按降序排列后的特征值和特征向量。將D和E代入式(6)可得到Renyi熵值V(p)，表達式如下：

(7)

(8)

目前常用的核函數包括徑向基(Radial basis function, RBF) 核函數、多項式核函數、感知器核函數等。其中RBF核函數應用廣泛且具有良好的適用性[27]，本研究選用RBF核函數，表達式為：

(9)

式中,k(xi,xj)為核矩陣元素；θ為寬度參數；‖· ‖為2范數。

θ的選取原則是使用較少的主成分來達到顯著的降維效果，即使得ξ值較大。通過比較不同的參數進行擇優，這里θ取值范圍為0～5。基于KECA將3種單變量干旱指數聯合后，即可得到綜合干旱指數MDI。

2.3 標準化綜合干旱指數SMDI計算方法

為與SPEI、SSMI及SRI比較，對MDI標準化處理。選取正態、對數正態、Weibull、Gamma及Logistic 5種具有較強適用性的理論函數擬合MDI序列，選取極大似然法估計分布函數參數，采用K-S檢驗進行擬合優度檢驗，并用AIC準則進行最優擬合函數選擇，之后將所得累積概率逆標準化，得到標準化綜合干旱指數SMDI。若5種分布都不滿足精度要求，則選擇非參數法計算SMDI。

本研究的干旱等級劃分參考SPI的統一標準(表1)。

表1 干旱指數等級劃分Table 1 Drought severity grades and their corresponding D-scale

2.3.1 K-S檢驗法 K-S檢驗是一種擬合優度檢驗法，用于判別樣本是否服從假定的理論分布。比較樣本的累積分布曲線與假設的理論頻率曲線并計算兩者的差值，若差值的絕對值在規定范圍內，表明該樣本服從假設分布。

以H0、H1分別表示樣本序列服從理論分布和不服從理論分布，構造統計量D：

D=max|F1(x)-F2(x)|

(10)

式中，F1(x)表示樣本累積分布概率，F2(x)表示理論分布概率。

選定顯著性水平α=0.05，當D>D(n,α)時，接受H1，否則接受H0。

2.3.2AIC準則AIC準則是一種常用的最優擬合函數選擇方法，能找出適配數據最優且自由參數個數較少的模型。AIC值越小，擬合效果就越好[29]。

AIC=2k-lnL

(11)

式中,k為模型自由參數個數，L為擬合的極大似然函數。

2.4 Mann-Kendall(M-K)趨勢檢驗法

M-K趨勢檢驗是一種分析時間序列趨勢的非參數檢驗方法，具有檢測范圍廣、定量化程度高等優點[30]。公式如下：

(12)

式中,S為構造的統計量；n為樣本序列長度；xj和xi是第j年和第i年對應數值。兩者之差大于0時，sgn值為1，小于0時sgn值為-1，等于0時sgn值為0。

當n大于8時服從正態分布，均值E(S)為0，方差Var(S)通過下式求得。

Var(S)=n(n-1)(2n+5)/18

(13)

構造標準化的檢驗統計量ZMK:

(14)

檢驗統計量ZMK無量綱單位，當ZMK的值為正(負)時，表示樣本序列隨時間的增加有增大(減少)的趨勢。當|ZMK|大于1.96時，認為在α=0.05的顯著性水平下變化趨勢明顯，反之則不明顯。

2.5 小波周期分析

小波方差圖可以確定一個序列中各種時間尺度相對擾動強度的大小，進而反映樣本序列的主要周期[31]。在經濟學、水文學等諸多需要分析時間序列周期性的領域被廣泛應用[32]。計算公式如下：

(15)

式中，Wf(a)為小波方差，a為尺度因子，b為平移因子，Wf(a,b)為小波系數。

3 結果與分析

當核函數的寬度參數θ為2時，黑河流域中上游各柵格第一主成分的累積貢獻率ξ>99%，所以取參數θ為2。本研究選取季時間尺度(3個月，即3—5月為春季，6—8月為夏季，9—11月為秋季，12月—次年2月為冬季)來分析黑河流域中上游干旱指數四季的趨勢變化、年代變化、周期變化和干旱發生頻率變化，并與歷史旱情事件比較以分析SMDI的適用性。

3.1 不同季節干旱指數序列分析

用黑河流域中、上游各柵格的不同干旱指數均值來表征中、上游的干旱狀況，并按季節分別比較各干旱指數間的差異(圖2～圖3)。

由圖2、圖3可知，同時考慮氣象、農業、水文3種要素的綜合干旱指數SMDI與分別只反映氣象干旱、農業干旱、水文干旱的單變量干旱指數SPEI、SSMI及SRI的年際變化不同，且在不同季節各干旱指數值亦有較大差異。單變量干旱指數值在同一時段有時相差很大，而綜合3種變量的標準化綜合干旱指數SMDI則在3者之間波動變化。在有些年份SMDI和3種單變量干旱指數值趨勢相反，這是因為3種單變量干旱指數都是按照同一分布計算得到。而SMDI則是使用5種參數分布和一種非參數法進行優選得到，各個柵格的分布類型不盡相同，這樣可在黑河流域中上游空間分布上達到最優。但由于不同柵格SMDI有正有負，在取柵格平均值時會抹去一些信息，從而表現為SMDI值偏大或偏小。2005年以后，SMDI值整體要比其他3種單變量干旱指數值偏小，說明近年來黑河流域中上游的干旱程度有所加劇。

3.2 干旱趨勢演變分析

分析黑河流域中上游4種不同干旱指數的四季變化趨勢(圖4)，根據干旱加重或減輕及變化是否顯著將其分為4類，分別用干旱顯著加重(T1)、干旱不顯著加重(T2)、干旱不顯著減輕(T3)、干旱顯著減輕(T4)來表示。

由圖4可知，氣象干旱除上游部分區域在夏、秋兩季干旱不顯著減輕，流域整體呈現干旱加重趨勢，且春、冬兩季尤為顯著；農業干旱除中游邊界地區有部分區域干旱加重，整體呈現干旱減輕趨勢；水文干旱上游表現為干旱顯著減輕趨勢，中游在夏、秋兩季有干旱不顯著減輕趨勢，而春、冬兩季呈現干旱不顯著加重趨勢。綜合了氣象、農業、水文的SMDI在趨勢演變上與氣象干旱相近，但變化趨勢不顯著，除秋季中游部分區域有干旱顯著加重趨勢，其他區域在四季變化上整體呈現干旱不顯著加重趨勢。統計流域干旱不顯著加重趨勢的格點占比，春季85.8%，夏季62.9%，秋季80.2%，冬季80.5%，全年77.6%，所以流域有77.6%的區域呈現干旱不顯著加重趨勢。

考慮SPEI是由降水與潛在蒸散發的差值序列擬合分布得到，分析黑河流域中上游12個氣象站四季降水和潛在蒸散發趨勢變化見表2。

表2 黑河流域中上游四季降水量(P)、潛在蒸散發量(ET0)趨勢變化Table 2 The tendency change of precipitation (P) and potential evapotran spiration (ET0)in the upper and middle Heihe River Basin

由表2可以看出，12個氣象站在51年中，除個別站點在某些季節降水有不顯著減小趨勢外，整體呈現增大趨勢，這與陳潔等[33]發現的近年來西北地區降水增多相符，但增多趨勢不顯著。而流域所有站點的潛在蒸散發量在各季節都呈增大趨勢，且多數站點都呈顯著增大趨勢。降水量的增加使得徑流量增多，表現為水文干旱逐漸減弱。但是降水量的增大不能完全抵消潛在蒸散發量的增大，因此氣象干旱表現為在流域大部分地區仍然加重。全球氣溫升高導致蒸散發量增大，同時也加速了冰雪融水，使得徑流量增多，表現為水文干旱逐漸減輕。由于近年來灌溉技術的不斷發展且抽取地下水進行農業灌溉，農業干旱有減輕趨勢。根據趨勢演變結果可見，單變量干旱指數在反映干旱變化趨勢時表征結果不相同。而SMDI則兼顧了氣象、農業、水文三者的特征，是對干旱狀況的全面反映，其趨勢變化結果表明黑河流域中上游整體仍呈干旱加重趨勢。

3.3 干旱年代變化分析

統計不同年代干旱的發生頻率，可進一步比較不同年代旱情的嚴重程度。本文研究序列為1966—2016年，由于1966—1969年只有4 a的數據，將其合并于20世紀70年代。將研究序列分為4段：1966—1979年、1980—1989年、1990—1999年、2000—2016年，代表20世紀70年代、20世紀80年代、20世紀90年代、21世紀初，統計得到黑河流域中上游SMDI在這4個不同時期的干旱發生頻率(圖5)。

由圖5可知黑河流域中上游在1980s旱情最輕，干旱發生頻率在季節尺度上均低于30%，其次是1970s，僅黑河流域上游在春季干旱的發生頻率較高。而1990s以及21世紀初流域旱情相對嚴重，同時這兩個時期在不同季節干旱的發生頻率又有所不同，1990s表現為夏、秋兩季干旱的發生頻率高，這與該時期黑河流域中上游城市化進程及工業發展迅速使得耗水量急劇增加進而導致嚴重水文干旱相符。而21世紀初則是在春、冬兩季干旱的發生頻率高。這是因為自2000年起黑河流域實施分水政策[34]等一系列積極舉措，使得夏、秋兩季干旱發生頻率降低。但春、冬兩季本身降水少，所以在21世紀初表現為春、冬兩季干旱發生頻率較高。

3.4 干旱空間特征分析

將干旱分為輕中旱和重特旱兩類(表1)，對SMDI統計得到黑河流域中上游不同季節干旱的發生頻率(圖6)。

由圖6可知，四季中，重特旱的發生頻率均比輕中旱低，表明黑河流域中上游極端干旱事件的發生頻率較低。其中春、夏兩季輕中旱的發生頻率要高于秋、冬兩季，夏季輕中旱的發生頻率最高。而重特旱的發生頻率則表現為夏、秋兩季高于春、冬兩季，秋季重特旱的發生頻率最高。整體來看，黑河流域中上游在夏季干旱的發生頻率最高，中游干旱的發生頻率整體上要高于上游。中游為流域重要的農業區，需要重點關注中游地區防旱減災工作。

3.5 干旱周期變化分析

采用小波函數對黑河流域中上游SMDI序列做周期分析，得到SMDI的小波方差圖(圖7、8)。

小波方差圖峰值顯示對應主要時間尺度，根據圖7、8將峰值從高到低排，可以得到黑河流域中上游SMDI四季主周期值(表3)。

由表3可知，黑河流域上游和中游的主周期在四季較為接近，但是在排序上略有區別。黑河流域中上游周期集中表現為43 a的長周期，15～23 a的中周期以及3～8 a的短周期。這與鐘鋒等[35]研究的黑河流域干旱存在19～23 a和14～15 a的中長振蕩周期及7～10 a和4～5 a的短振蕩周期結論基本一致。同時與太陽黑子活動周期約22 a、厄爾尼諾事件周期2～6 a相對應[36]，這表明黑河流域干旱的周期變化可能是太陽黑子活動及厄爾尼諾事件的綜合驅動結果。其中僅上游春季第一主周期為43 a的長周期，其余階段第一主周期皆為中周期和短周期，這與黑河流域中上游干旱發生頻繁相對應。根據周期情況可預測未來干旱大致出現時間。

表3 黑河流域中上游SMDI四季主周期值Table 3 The main cycle values of SMDI in four seasons in the upper and middle Heihe River Basin

3.6 典型干旱事件驗證

黑河流域中上游在行政劃分上屬于青海和甘肅兩省。根據《中國氣象災害大典》[37-38]：青海卷和甘肅卷，青海、甘肅兩省水資源公報以及中國水旱災害公報記錄：1969年甘肅全省大部分地區旱災嚴重，春夏連旱，重旱災涉及16縣，同年青海省海南、海西州發生嚴重干旱；1997年甘肅全省均發生不同程度的干旱，是近60年來干旱范圍最大、旱情最重的一年，同年青海省的祁連、野牛溝地區旱災嚴重；2009年冬季西北大部區域發生干旱。

為進一步展現SMDI的干旱監測能力，選取1969年春、1997年秋和2009年冬與同時期的其他3種單變量干旱指數干旱監測情況進行比較(圖9)。

由圖9可知，1969年春，SPEI與SRI的指數值表現為無旱，與旱情記錄的實際情況不符，而SMDI的監測結果與SSMI相近，都可以監測到較大區域的干旱。1997年秋4種干旱指數值在流域大部分區域均表現為發生干旱，但SPEI的監測結果在部分區域表現為無旱。而SMDI和SSMI及SRI識別區域基本一致，但在紅框所示區域不同，此區域SMDI值比3種單變量干旱指數小，即SMDI所表征的干旱等級高于3種單變量干旱指數所表征的干旱等級。歷史旱情記錄1997年秋是近60年來旱情最嚴重的一年，可以看出SMDI與本年旱情更為相符。2009年冬，西北地區大部區域發生干旱，但SRI和SSMI在黑河流域中上游表現為無旱，而SPEI識別出的發生干旱的流域面積沒有SMDI識別出來的大。與以熵權法聯合SPEI及SSI構建的氣象農業綜合干旱指數AMDI[39]相比，AMDI在旱情驗證上不能很好反映氣象和農業之間的非線性關系，表現為在春旱驗證上與歷史旱情事件沒有較好對應。而基于KECA構建的SMDI則能夠反映不同類型干旱之間的非線性關系，旱情驗證結果優于AMDI，體現了KECA在反映不同變量間非線性關系上具有優越性。說明SMDI比單變量干旱指數及熵權法構造的AMDI在黑河流域中上游的干旱監測適用性更好。

4 結 論

本文以黑河流域中上游為研究區，使用新的降維方法KECA聯合了氣象、農業、水文3種變量來構建綜合干旱指數SMDI，并對黑河流域中上游干旱的時空特征進行分析，同時與典型歷史旱情事件作驗證比較，主要得到以下結論：

1)黑河流域中上游除部分區域表現為干旱不顯著減輕外，整體仍表現為干旱加重趨勢，且有些區域加重趨勢顯著；在20世紀90年代夏、秋兩季和進入21世紀以來春、冬兩季旱情最為嚴重；干旱存在43 a的長周期，15～23 a的中周期以及3～8 a的短周期。

2)四季中，重特旱的發生頻率均比輕中旱的小，且干旱的發生頻率在上游低于中游，綜合來看，夏季干旱的發生頻率最高。

3)歷史典型旱情驗證結果表明，基于KECA構建的SMDI在干旱的捕獲區域以及極端干旱事件的反映上優于其他3種單變量干旱指數及熵權法構建的綜合干旱指數，能夠反映多種干旱變量之間的非線性關系，在黑河流域中上游干旱監測中有著更好的適用性。

在氣候變化與人類活動影響更劇烈的大環境下，氣象、農業、水文3種干旱之間的關系更為復雜，傳統的單變量干旱指數難以全面表征干旱的綜合特征，綜合干旱指數的發展與應用，能夠彌補傳統單變量干旱指數的不足，為決策者提供有益參考。本文使用KECA構建的SMDI，可表征氣象、農業和水文干旱之間的非線性關系，適用于黑河流域中上游干旱的監測。該方法是否存在區域性差異，在其他區域是否適用，需作進一步驗證。