水文序列趨勢及形態變化的表征方法

張姝琪，張洪波,2，辛 琛，南政年，李哲浩

(1.長安大學環境科學與工程學院，陜西 西安 710054； 2.長安大學旱區地下水文與生態效應教育部重點實驗室，陜西 西安 710054； 3.陜西省江河水庫管理局，陜西 西安 710018)

近年來，受氣候、土地覆蓋等自然要素變化和人類活動擾動的影響，流域的水文過程發生變異，水文序列不再滿足一致性要求，導致變化環境下基于一致性假設的水文頻率計算方法備受質疑[1-3]。在水文序列非一致性的研究方面，數據序列的趨勢變化檢驗一直是一項重要的研究內容，其對水文分析、模擬、預測以及水文時空變化規律都有著關鍵性影響，同樣也對變化環境下流域水資源管理與決策具有深遠意義[4]。

現階段，識別水文序列趨勢變化的方法有很多，如回歸分析法[5]、Spearman法[6]、滑動平均法[7]等。其中最常用的是回歸分析法，常用的回歸類型有一元回歸、多元回歸、指數回歸、對數回歸、多項式回歸等[8]。回歸分析法也存在一定的不足，如當水文序列變化的無序性較大時，回歸系數偏小，趨勢特征不明顯，將會難以判斷。為此，有學者將Mann-Kendall非參數檢驗法(M-K法)引入氣象水文序列的趨勢分析中，用以判斷氣象水文序列發生的趨勢性變化[9-10]。M-K法是國際氣象組織推薦的一種判別氣象水文等自然變量趨勢變化的方法，由于不要求分析樣本遵從一定分布，且不受個別異常值和缺失值的干擾，故在氣象水文領域應用十分廣泛。此外，M-K法原理簡單，計算簡便[11]，能夠很好地描述氣象水文數據序列的趨勢特征[12]，并在以往的實踐中已取得很好的應用效果，如邸擇雷等[13]基于呼倫貝爾草原新巴爾虎右旗1958—2016年的氣溫和降水數據資料，應用M-K法得出了年平均氣溫顯著上升、降水不顯著下降的趨勢檢驗結論。曹靜等[14]以科爾沁沙地西部為研究區域，利用1981—2016年氣象資料，使用M-K法對蒸發皿蒸發量的變化進行了趨勢研究，結果表明，在過去36年里該地區的蒸發皿蒸發量呈顯著下降趨勢。李佳秀等[15]以新疆地區為研究對象，使用M-K法對氣溫和降水資料進行分析，結果表明近50年來，降水呈增長趨勢。張巖等[16]采用M-K法分析了1957—2012年三江源區徑流的變化趨勢，結果表明47.3%區域的徑流變化趨勢不顯著。

與許多氣象變量不同，水文變量由于受到人類活動的高強度擾動，不再呈現純粹的遞增或遞減趨勢變化，而是更趨多元化變化，通過文獻梳理，可大體將這種變化總結為加劇、抵消以及紊亂3種情形。所謂加劇，指持續性的高強度人類活動使原本較為顯著的趨勢變化更加嚴重，如河道階梯性取水引發的下游徑流量的持續性衰減；抵消，指人類活動使原本較為顯著的水文趨勢變化變緩或消失，如枯水期河川徑流日益趨于干涸，但是上游建庫后豐水期的水被調蓄至枯水期，導致枯水期水量增加，原有徑流序列中的枯水減少趨勢隨之不顯著甚至消失；紊亂，指人類活動使得原本顯著的水文趨勢變化趨于無序和復雜，原本單調的趨勢變化開始震蕩，表現出愈來愈強的非線性特征，導致很難確定水文序列的整體走向。

眾所周知，在高強度的人類活動影響下，水文要素發生了不容忽視的變化[17]，水文序列不再是平穩序列，趨勢變化的特征也日趨復雜，不再是簡單的線性關系。盡管前文所提到的M-K法具有不需要樣本服從某一分布且不受異常值干擾等優勢，但受其原理(秩次)制約，檢驗結果始終與秩的大小變化息息相關[18]。當兩個水文序列恰好具有相同的秩次關系時，其M-K檢驗結果必定一致，從而獲得兩個序列具有相同的趨勢特征的結論。事實上，兩個序列可能還存在一些有關趨勢變化的個性化信息或其他重要形態信息，因為M-K檢驗而被丟失或未能體現。可見，M-K法在識別序列趨勢變化的細節特征或者說形態變化上尚有不足。此外， M-K法的統計檢驗值僅表征序列的秩次關系，并不能用于實際水文序列的空間對比分析，導致M-K法在空間應用的拓展上一直舉步維艱。為此，本文基于M-K法提出一套水文序列趨勢及形態變化分析的指標體系與方法，通過提取水文數據序列的趨勢及形態特征，可實現水文序列變化的細節刻畫和空間對比分析。本文以渭河流域為研究區開展案例研究，以驗證指標體系方法的適用性。

1 研究區概況及數據來源

渭河流域西起鳥鼠山，東至潼關，北起白于山，南至秦嶺，流域范圍為北緯33°50′～37°18′、東經104°00′～110°20′，面積13.48萬km2[19](圖1)。渭河干流全長818 km，平均比降1.3‰，占黃河流域面積的18%，是黃河的第一大支流[20]。流域地處黃土高原南緣，秦嶺北麓，位于干旱地區與濕潤地區的過渡地帶[21]。流域的降雨時空分布極不均勻，年均降雨量介于200～900 mm之間，且年內、年際變化較大，導致流域內的水資源分布亦呈現時空上的差異性變化[22-23]。渭河流域是通往西北、西南的咽喉要地，也是我國重要的糧棉油產區和工業生產基地之一。特殊的地理位置使其在西部大開發戰略中具有重要地位，其流域水資源狀況對黃河水資源的豐枯變化和開發利用具有重要影響[24]。然而隨著大規模的工業、農業生產活動對水資源需求的不斷增大，渭河徑流持續縮減，趨勢變化強烈且日趨復雜，已引發下游河道淤積、水土流失嚴重、洪澇災害突出、水污染嚴重等一系列生態環境問題，導致流域生態環境系統遭到破壞，生態災害頻繁發生[25]。因此，了解流域內徑流序列的形態變化特征，尤其是趨勢變化，不僅可科學地認識變化環境下渭河流域的水文演變規律，也對指導區域未來的經濟布局以及水資源安全保障有著極其重要的意義。

圖1 渭河流域水系及水文站分布圖

為更好地描述渭河流域徑流趨勢變化的空間差異性，本文選取渭河干流上的北道、林家村、魏家堡、咸陽、華縣5個干流水文站和張家山、狀頭站2個支流水文站的實測年徑流資料作為基礎數據資料。資料長度統一為51年(1959—2009年)，資料來源于黃河水利委員會水文局。

2 研究方法

2.1 去趨勢預置白

眾所周知， M-K法的使用前提是檢驗數據之間相互獨立，然而很多水文數據并不滿足該要求，故本文采用去趨勢預置白(trend-free pre-whitening，TFPW)方法[26-27]對序列進行去一階自相關處理，從而使待檢序列滿足M-K法的使用要求。具體步驟[28]為：設一待檢序列為Xt(t=1, 2, …，n，n為序列長度)，則有

Yt=Xt-βt

(1)

Yt′=Yt-γYt-1

(2)

Yt″=Yt′+βt

(3)

(4)

式中：Yt為扣除趨勢后的殘余序列；Yt′為去除自相關項后的獨立白噪聲序列；Yt″是經過TFPW處理后得到的新序列；β為待檢序列的坡度，可通過Theil-Sen法[29]即式(4)計算；γ1為Yt的自相關系數，階數因序列差異而有所不同，常規計算中多選一階系數。需要注意的是，如果γ1足夠小，可認為Yt是一個獨立序列，即可不處理，直接使用M-K法進行趨勢檢驗，否則需要帶入式(2)(3)計算，并最終將Yt″代入M-K檢驗。

2.2 M-K趨勢檢驗

M-K法已在徑流、氣溫、降水等氣象水文要素的序列趨勢檢驗上[30-32]廣被認可。

先定義檢驗統計量S：

(5)

式中sign()為符號函數。當(Xi-Xj)<0、≥0時，sign(Xi-Xj)分別為-1、0或1；當檢驗統計量S≥0、<0時，標準化統計量Z可計算：

(6)

式中Z為正值表示上升趨勢，為負值表示下降趨勢。Z的絕對值大于等于1.96時可認為其通過了置信度為95%的顯著性檢驗，即趨勢顯著。

梳理M-K法的檢驗原理，不難發現M-K法主要是通過秩將非正態分布的原序列轉化為服從正態分布的秩序列，并由此采用統計檢驗方法判斷序列的趨勢性變化。從應用范圍上看，M-K法因為假設前提較少，故其應用更廣，對未經處理的數據所得的檢驗結果也更加可信。然而，M-K法也存在一定的不足。M-K法雖然可以給出一個定量化的趨勢統計量，即在整個時域上對數據序列進行趨勢判斷是行之有效的，但并不能區分具有相同Z值的數據序列趨勢及形態之間的差異性變化。水文工作者或流域管理人員依據同一流域不同站點徑流資料的M-K檢驗結果，并不能準確了解徑流序列的變化，對于序列中某一時間段的趨勢行為或者是哪些變化引起整個時域上的趨勢演變，也無法有效掌握。鑒于此，本文提出構建一套表征趨勢及形態變化的指標識別體系與方法來彌補M-K法在趨勢分析及形態化細節識別上的不足。

盡管水文序列形態變化形式眾多，且可能存在掩蔽問題，但如果能通過指標體系宏觀描述序列的這些趨勢和形態變化特征，則可極大程度地幫助讀者對序列變化有一個相對直觀的了解與認識，并可通過指標值的差異較為準確地判別不同序列之間的差異性。基于以上分析，筆者結合M-K法，提出了表征水文序列趨勢和形態變化的指標體系，主要包括：①M-K統計量Z值，反映待檢序列的總體趨勢；②擬合拋物線系數a，反映待檢序列趨勢形態變化的具體表現形式；③平穩序列(待檢序列去除趨勢成分)的集中度與集中期，反映徑流波動的能量分布及最大能量分布的位置；④方差序列擬合的拋物線系數a′，描述徑流波動在時域上的漸變性特征。具體指標計算方法如下：①使用TFPW方法對實測水文序列進行去自相關處理，得到滿足M-K法應用假設的待檢序列；②對待檢序列進行M-K非參數趨勢檢驗，得到序列在整個時域上的總體趨勢變化指標，即Z值；③對待檢序列進行二次擬合，得到擬合方程y=ax2+bx+c中的拋物線系數a，通過a可判斷趨勢形態變化的非線性表現形式；④對待檢序列作趨勢成分剔除處理，得到平穩序列，并對其集中度、集中期進行統計計算，揭示徑流序列波動的能量分布特征；⑤以平穩序列為基礎，通過滑窗計算，得到均方差序列，并對其進行拋物線擬合，得到拋物線系數a′，表征徑流序列波動的漸變性特征。

3 結果分析

3.1 指標計算

以咸陽站和狀頭站實測徑流序列為例，闡述各指標的計算過程。

a. M-K統計量Z值與拋物線系數a。將咸陽站及狀頭站的實測徑流序列代入TFPW前置處理程序進行自相關性剔除，得到新的徑流序列，即待檢序列，其偏自相關性統計檢驗結果見圖2。由圖2可知，兩站實測序列經自相關剔除處理后，原序列中的顯著自相關性較好地被控制在可接受范圍內。

圖2 咸陽站和狀頭站徑流序列TFPW處理前后的偏自相關性檢驗結果

然后，對待檢序列實施M-K趨勢檢驗和二次擬合(拋物線)計算。M-K檢驗結果顯示咸陽站及狀頭站待檢徑流序列的統計量Z分別為-4.548 4和-1.840 3，即咸陽站待檢序列通過了95%的顯著性檢驗，這表明咸陽站的年徑流序列在整個時域上表現為顯著的下降趨勢，而狀頭站待檢序列為不顯著下降趨勢。將二次擬合結果與徑流序列一并繪制于圖3中，可發現二次擬合曲線也表明序列呈現為下降趨勢。其中，咸陽站及狀頭站待檢序列的二次擬合拋物線系數a分別為0.0171和-0.0038，數值均較小，即為不顯著的凹減和凸減形態。就a值來看，咸陽站年徑流序列呈現先減后增的非線性趨勢，而狀頭站年徑流序列呈現先增后減的非線性趨勢。

b. 均方差序列拋物線系數a′。由M-K法的檢驗效能分析可知，M-K趨勢檢驗對序列的離散程度或分段離散特征無法有效表征，使其在形態變化上的指代意義有所欠缺。為了反映序列形態變化中的離散特征(方差漸變性和能量分布特征)，本文采取滑窗的方式對剔除趨勢成分的平穩序列逐段進行方差統計，得到一個連續的均方差序列，并對均方差序列進行二次擬合，獲得均方差序列的拋物線系數a′，并以此描述序列方差的漸變性特征，所得結果見圖4。由圖4可知，咸陽站及狀頭站徑流序列的a′分別為0.017和0.000 7，這表明兩站的方差序列呈不顯著的凹減型變化，且咸陽站的方差曲率大于狀頭站。

圖3 咸陽站和狀頭站年平均流量序列及其二次擬合曲線

圖4 咸陽站和狀頭站均方差序列及其二次擬合曲線

c. 集中度與集中期。方差序列的拋物線系數a′主要體現了年徑流序列在時域上的非線性漸變波動特征(方差漸變性)，如凸增型、凹減型等，但是并不能體現波動變化的能量分布。湯奇成等[33]在研究年徑流量在一年中的分布情況時，曾借鑒表示降水量年內分配的向量法，提出用集中期(concentration period)和集中度(concentration ratio)兩個概念來反映徑流年內分配的各種特性，效果令人滿意，Zhang等[34]也利用集中度集中期反映了黃河流域降雨在年內的分布情況。本文將借鑒集中期和集中度的概念來描述年徑流波動的能量分布情況。這里，方差集中度的大小主要反映一個序列的方差變化是均勻分布還是集中在某一年或某幾年。方差集中度越大表示該序列徑流波動的能量越集中；反之，序列的波動能量更趨分散或平緩。而方差集中期則指示了最大的波動所對應的年份，即徑流波動的能量極值點。

集中度、集中期的分析對象為剔除趨勢成分的平穩序列。具體計算步驟如下：①將整個研究時段認為是一個圓周，可分為N份，N表示數據序列的總個數；②每一個數據位置(年)相對應的角度為360°/N；③根據文獻[33]中所列公式，進行Matlab編程計算。由計算結果可知，兩站方差序列的集中度均小于33%，說明方差變化分布較為均勻，且狀頭站更顯著，這與圖4所展示的結果是一致的。在能量極值點上，咸陽站出現在1967年，而狀頭站出現在1982年，表明這兩個時間段徑流序列能量波動較大。

d. 綜合分析。以咸陽站年徑流序列的趨勢及形態變化特征為例，開展綜合分析。序列主要特征為：①整個序列呈現顯著的減小趨勢(M-K統計量Z為-4.548 4，且絕對值大于1.96)；②減小的形式為不太顯著的凹減型(實測序列拋物線系數a為0.0171)；③徑流波動在時域上呈現凹減型，即波動存在一個先減小后小幅增加的過程(方差序列拋物線系數a′為0.017)；④徑流波動能量較為分散，即波動相對平穩(集中度小于33%)，且波動的高值區主要集中在1967年附近，綜合兩種信息以及其下降趨勢，認為在1967年附近應該存在一個大的向下跳躍。

當獲取了這些形態變化特征后，可根據這些特征所攜帶的信息，大體勾勒所描述的徑流序列的趨勢及形態變化，并對徑流序列的表現形式給出定性的分析。以咸陽站為例，具體過程見圖5。

由特征①、②可勾勒出咸陽站徑流序列的整體趨勢變化(圖5(a))；由特征③可勾勒徑流序列的波動應呈現為先強后弱，再微弱漸強的過程(圖5(b))；由特征④可修訂圖5(b)，即除先強后弱，再小幅漸強的過程外，還應體現整體波動差異不是特別大，且在1967年附近存在徑流波動的高能區的特點(圖5(c))。最后將圖5(a)和圖5(c)相融合，可獲得整體的重構結果(圖5(d))。

同理，可重構狀頭站徑流序列的整體趨勢變化，如圖6所示。需要說明的是，由于只是示意圖，因此數據關系可能不完全閉合。將其與圖3中的咸陽站和狀頭站的實測年徑流序列進行對比，可發現重構的示意性序列能較好地體現原序列的趨勢與形態特征。可見，應用本文所提出的指標體系(Z，a，a′，集中度，集中期)，可較好地描繪原系列的趨勢與形態特征，并將其量化，讀者可借由各指標值宏觀了解掌握實測系列的趨勢與形態變化。

圖5 基于指標特征信息的咸陽站徑流序列趨勢與形態變化特征重構示意圖

圖6 基于指標特征信息的狀頭站徑流序列趨勢與形態變化特征重構示意圖

3.2 指標體系在流域空間分析中的應用

為驗證以上指標體系的空間適用性，本文進一步分析了渭河流域其他各站點的序列趨勢與形態變化，并進行了空間對比分析。

應用上述方法對其他站點進行指標計算與分析，所得結果見圖7～9。

由圖7可知，TFPW方法有效消除了各站點序列中自相關性的影響，可滿足M-K法應用假設。將待檢序列分別進行M-K檢驗，結果分別為-3.752 5、-4.759 6、-4.970 8、-5.133 2、-3.443 8，由此可知，5站徑流序列均顯示出顯著的下降趨勢。圖8給出了各站點序列的二次擬合結果。由圖8可知，M-K檢測結果與二次擬合結果基本一致，a系數分別為0.013 9、0.018 4、0.013 8、-0.004 2、-0.022 6，即華縣、魏家堡及林家村站呈現凹減型下降趨勢，而北道站、張家山為凸減型下降趨勢變化，除張家山站外，其他站點非線性變化均不顯著。

圖7 不同站點年徑流序列TFPW處理前后的偏自相關性檢驗結果

圖8 不同站點年徑流序列及其二次擬合曲線

圖9 各站均方差序列及其二次擬合曲線

圖9顯示了各站的方差序列及其二次擬合曲線，其拋物線系數a′分別為0.0161、0.018、-0.0016、0.009、0.0047，其方差的變化依次是凹減下降末尾漸平、快速凹減下降末尾漸升、凸減、凹減末尾漸升、緩慢凹減。

各站徑流序列的集中度和集中期統計結果見表1。從表1可知，除華縣站外(集中度大于70%)，各站的徑流波動均較為平穩(集中度小于20%)。其中，華縣站的集中度最大，張家山站最小，其結果也與方差變化的實際情況(圖9)較為相符。而集中期則主要分布在4個時段，即1970年、1963年、1973年及1960年。

表1 渭河流域各站集中度和集中期

最后，依據各站不同指標所指代的信息進行綜合分析，并與實際徑流序列進行對比，結果見圖10。圖10顯示了華縣、魏家堡、林家村、北道、張家山5個水文站實測水文序列及形態變化，綠色區域表示該序列的集中期區段，即徑流波動最大的位置。由圖10可知，5個指標較好地反映了實際水文序列中的形態變化，由此可以看出本文所提出的方法是行之有效的，指標體系可較好地體現原序列的趨勢及形態變化特征。

鑒于以上定量化指標，將渭河流域7個站點的指標值分別統計于表2中，可以發現渭河流域各站點徑流序列均表現為較為明顯的下降趨勢(狀頭站為不顯著下降)，但其趨勢形態略有不同，呈現差異性的非線性特征。另外，序列波動特征也不一致，在方差漸變性及能量分布上有一定區別。通過結果的對比，可發現本文提出的面向趨勢和形態變化的指標體系是切實可行的，能夠分析出序列趨勢及形態特征的差異，有利于全面認識渭河流域徑流序列的變化，實現同一流域不同站點或不同流域水文序列之間的比較與分析。

圖10 渭河流域5站年徑流序列及形態變化分析

表2 渭河流域上游至下游及支流7個水文站指標值

4 結 論

a. 基于秩的M-K法因檢驗機理限制，對序列的形態變化表征不夠充分，可能遺漏重要細節趨勢信息，從而影響對水文序列規律的認識。

b. 通過表征指標體系中的5個指標(Z、a、a′、集中度、集中期)能簡單有效地重構示意性水文序列，實現對序列趨勢及形態變化的細致刻畫與表征。

c. 借助指標體系的定量化優勢，可以實現同一流域不同水文站點或不同流域水文站點間的序列對比分析。

目前研究對象僅局限于渭河流域，對于該方法適用范圍的確定還需進一步補充驗證。同時，本文提出的集中期僅指示序列波動最大點的位置，對于其波動方向以及影響范圍尚無法準確界定，仍需進一步完善。