水文時間序列高階自相關性識別及其對水文分析的影響研究

，，2，，

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

水文時間序列是指某種水文特征值隨時間而變的一系列觀測值，在社會水文學的視點下，它也可以被看作是由氣候條件、下墊面條件和人類活動等因素綜合作用結果的表征[1]。一般水文時間序列可分解為確定性和隨機性兩種成分，確定性成分主要表現(xiàn)為水文現(xiàn)象的趨勢變化和周期變化等，其主要受人類活動、氣候因素和下墊面因素的宏觀尺度驅(qū)動，使水文時間序列產(chǎn)生緩慢的漸變或劇烈的突變，亦或有規(guī)律的波動；而隨機性成分則表現(xiàn)為水文現(xiàn)象的相依性和純隨機變化，其主要受氣候、下墊面等各種隨機因素影響，并驅(qū)動水文時間序列的隨機性變化和不確定性變化[2]。

近年來，受全球氣候變化及高強度人類活動(如水庫蓄水、淤地壩攔水、引水灌溉及城鄉(xiāng)建設等)的影響，流域的水循環(huán)過程和水文物理形成機制都發(fā)生了顯著的改變[3-4]。很多地區(qū)的實測年徑流、降水、泥沙等水文時間序列都表現(xiàn)出相對較大的自相關性[5-6]，使得一些檢驗結果的準確性受到一定影響。目前，關于水文序列自相關性的研究工作主要集中在序列的一階自相關性上。如在國內(nèi)外相關領域廣泛應用的由Von Storch和NaVarra提出的預置白[7](Pre-whitening，PW)方法和Yue等提出的去趨勢預置白[8](Trend-free pre-whitening，TFPW)方法等，其主要作用都是消除或降低序列中的一階自相關性，以滿足數(shù)據(jù)獨立性的要求。然而，RazaVi和Vogel[9]的研究表明，在對樹木年輪序列進行預置白(Pre-whitening，PW)處理后會使時間序列中包含的重要信息丟失，同時會導致對氣候和水文極端條件的低估，并建議在使用預置白方法時需要考慮其在長期自相關性上的陷阱。但對于序列高階自相關性的剔除方法或影響，這些文獻均未有過多探討。目前，水文領域有關高階自相關問題的論述，并不多見。張洪波等[10]提出基于EMD的改進RBF方法來控制高頻分量的預測誤差，進而實現(xiàn)提高徑流預測精度的方法。其中作者根據(jù)各分量的自相關分析結果，篩選存在顯著自相關的時滯，完成選階后構建神經(jīng)網(wǎng)絡的輸入向量，盡管考慮了各顯著階數(shù)序列，但未說明其對徑流預測效果產(chǎn)生的影響。Sen[11]提出的一種新的用于統(tǒng)計檢驗的自相關處理方法，即過白化處理法(over-whitening，OW)，其通過在時間序列的標準化序列上添加一個白噪聲，可實現(xiàn)一階和高階自相關性的綜合去除，為高階自相關移除提供了手段。此外，Julian[12]對德國1978—2009年17個觀測站的降水序列中18O、2H穩(wěn)定同位素進行趨勢變化研究，結果表明在建模過程中考慮高階自相關性的影響時，有3個站的序列表現(xiàn)出顯著的趨勢變化，而廣泛應用的僅考慮一階自相關性的趨勢分析方法未能充分考慮到穩(wěn)定同位素序列中高階自相關性的存在。再者，張洪波等[13]在傳統(tǒng)過白化處理法的基礎上，提出了改進措施，采用分段過白化的處理手段實現(xiàn)了趨勢低損情況下的高階顯著自相關剔除。綜上所述，在對水文序列進行分析時，僅僅考慮序列的一階自相關性是不夠充分的。與此同時，盡管水文序列的高階自相關議題在一些文獻中出現(xiàn)過，但其并未引起廣泛的關注或普遍的共識。因此，非常有必要對水文序列中高階自相關性做以深入的解析并探討其對現(xiàn)有技術手段的影響。

為此，本文分別選取了陜西境內(nèi)漢江流域的2個水文站、渭河流域的6個水文站以及無定河流域的11個水文站的實測年徑流序列為研究對象，通過分段去趨勢方法，剔除序列中的趨勢成分，再對去趨勢序列進行自相關性分析，以確定序列高階自相關性在不同流域的空間分布。同時，文章進一步以BP神經(jīng)網(wǎng)絡為基礎，建立降水-徑流預測模型，旨在探討加入高階降雨序列對徑流預測效果的影響，進而提出水文時間序列分析中高階自相關性的處理方法。

1 研究區(qū)概況與數(shù)據(jù)來源

1.1 研究區(qū)概況

1.1.1漢江流域

漢江是陜南地區(qū)最大的河流，也是長江較大的一級支流，發(fā)源于陜西省寧強縣潘家山，由西向東流經(jīng)漢中、安康兩市進入湖北省境。漢江在陜西境內(nèi)河長654 km，集水面積4.86萬km2，多年平均徑流量247億m3。徑流量地區(qū)分布不均勻，總的趨勢是南岸多于北岸，即大巴山多，秦嶺少。陜南地區(qū)降水充沛，時空分布不均，年均降水量1 050 mm，秦嶺南坡年均降水量在600～800 mm之間，巴山北坡米倉山區(qū)年降水量最大，可達1 400 mm以上。流域內(nèi)河川徑流的補給主要以降雨形成的地表徑流為主，地下水補給為次[16]。

漢江沿岸的地貌特征以丘陵、盆地為主，東西長約400 km，南北寬約3～60 km，漢江流貫中部，形成以漢中盆地和石泉-安康盆地為主的多個盆地區(qū)[17]。

1.1.2渭河流域

渭河是黃河的一級支流，起源于甘肅省渭源縣，流經(jīng)甘肅、寧夏、陜西，于潼關匯入黃河，流域總面積約13.43萬km2，全長818 km，陜西境內(nèi)河流總長約502 km，流域面積6.71萬km2。渭河流域年平均氣溫12 ℃～14 ℃，年平均降雨量為400～800 mm，降雨主要集中在6—9月，降雨量占全年的60%左右[13-14]。

渭河流域地形總特點是西北高東南低，陡坡自西向東逐漸變緩，河谷自上游至下游逐漸變寬。渭河流域北岸支流多發(fā)源于黃土丘陵和黃土高原，相對源遠流長，比降較小，含沙量大；南岸支流均發(fā)源于秦嶺山區(qū)，源短流急，谷狹坡陡，徑流較豐，含沙量小[15]。

1.1.3無定河流域

無定河是榆林地區(qū)最大的河流，它發(fā)源于定邊縣白于山北麓，流經(jīng)米脂、綏德到清澗縣川口以南20 km處注入黃河。全長491.0 km，流域面積30 260 km2，陜西境內(nèi)河長442.8 km，流域面積21 049.3 km2，流域降水量從東南向西北遞減，多年平均降雨量為409 mm[10]。

按地形地貌及水土流失特點，流域可分為風積沙丘沙蓋區(qū)、黃土斜梁高山區(qū)和黃土丘陵溝壑梁峁區(qū)三大區(qū)域。流域西北部為風積沙丘沙蓋區(qū)，約占全面積的54.3%，氣候干燥，地面起伏平緩，水網(wǎng)稀少；西南部為黃土斜梁高山區(qū)，約占全面積的11.4%，地面呈梁峁狀丘陵；東南部為黃土丘陵溝壑梁峁區(qū)，約占全面積的34.3%，地表支離破碎，梁峁起伏，溝壑縱橫。3個流域的自然地理條件對比見表1。

注：表1相關統(tǒng)計數(shù)據(jù)出自于陜西省2016年水資源公報

1.2 數(shù)據(jù)來源

本研究采用的數(shù)據(jù)主要有漢江流域石泉、安康水文站，渭河流域華縣、臨潼、咸陽、魏家堡、林家村、北道水文站以及無定河流域韓家峁、橫山、趙石窯、白家川、殿市、馬湖峪、青陽岔、李家河、曹坪、丁家溝、綏德水文站等共19個水文站點的實測年徑流及控制流域的面降雨數(shù)據(jù)。其中，年徑流數(shù)據(jù)來源于黃委水文資料匯編，降水數(shù)據(jù)來源中國氣象資料共享網(wǎng)，站點數(shù)據(jù)經(jīng)加權平均后獲得面降雨量數(shù)據(jù)。各流域站點分布見圖1。

圖1 各流域站點分布

2 研究方法

2.1 分段去趨勢

受水文序列非一致性的影響，其本身表現(xiàn)出時域上的差異性特征，而非單調(diào)性的趨勢變化。因此，在進行自相關性分析前，需要對水文序列進行去趨勢處理。本文使用了一種基于回歸變異的分段去趨勢方法來去除徑流序列的趨勢成分。其具體步驟如下：①采用Tomé和Miranda[18]提出的分段線性擬合技術，對水文時間序列進行回歸變異識別，得到該序列的多個回歸變異點；②擬合各分段回歸系數(shù)與回歸方程，分段去除趨勢，得到去趨勢序列，即平穩(wěn)序列；③對去趨勢序列進行自相關性分析。

2.2 BP神經(jīng)網(wǎng)絡

BP(back propagation)神經(jīng)網(wǎng)絡是1986年由Rumelhart和McClelland為首的科學家提出的一種按照誤差逆向傳播算法訓練的多層前饋神經(jīng)網(wǎng)絡，是目前應用最廣泛的神經(jīng)網(wǎng)絡模型之一。BP網(wǎng)絡能學習和存貯大量的輸入-輸出模式映射關系，而無需事前揭示描述這種映射關系的數(shù)學方程。它的學習規(guī)則是使用最速下降法，通過反向傳播來不斷調(diào)整網(wǎng)絡的權值和閾值，使網(wǎng)絡的誤差平方和最小[19-20]。

本文中，BP神經(jīng)網(wǎng)絡主要基于MATLAB軟件中的NEWFF函數(shù)開展，構建多輸入單輸出的神經(jīng)網(wǎng)絡模型，其中訓練誤差GOAL設置為0.000 1；其他參數(shù)均采用默認值。經(jīng)過反復迭代后，可確定各輸入變量的BP神經(jīng)網(wǎng)絡模型，得到預測結果。

3 計算結果

3.1 高階自相關識別

3.1.1漢江流域

通過分段去趨勢方法對漢江流域的石泉和安康水文站的實測年徑流序列進行處理，得到其去趨勢序列，并采用偏自相關(PACF)函數(shù)對其進行自相關性檢驗。結果表明石泉、安康兩站點的去趨勢年徑流序列的偏自相關系數(shù)值均未超出臨界值，即不存在顯著的高階自相關性，見圖2。

3.1.2渭河流域

利用分段去趨勢方法對渭河干流6個水文站的年徑流序列進行處理，得到其去趨勢序列，并采用偏自相關(PACF)函數(shù)對其進行自相關性檢驗。結果與漢江流域類似，華縣、臨潼、咸陽、魏家堡、林家村、北道6個站點的去趨勢年徑流序列的偏自相關系數(shù)值均未超出臨界值，即均不存在高階自相關性，見圖3。

a) 石泉

b) 安康圖2 去趨勢年徑流序列PACF檢驗

a) 華縣

b) 臨潼

c) 咸陽

d) 魏家堡

e) 林家村

f) 北道圖3 去趨勢年徑流序列PACF檢驗

3.1.3無定河流域

將分段去趨勢方法應用于無定河流域11個水文站的年徑流序列，得到去趨勢序列，并采用偏自相關(PACF)函數(shù)對其進行自相關性檢驗。結果表明韓家峁站存在6階顯著高階自相關性；橫山站存在7階顯著高階自相關性，6階偏自相關系數(shù)值未超出臨界值，但數(shù)值較大，在后續(xù)計算中應考慮6階的影響；趙石窯站存在11階顯著高階自相關性；白家川站存在2階和5階顯著高階自相關性；殿市存在一階自相關；馬湖峪、青陽岔、李家河、曹坪、丁家溝、綏德此6個站點偏自相關系數(shù)值均未超出臨界值，即均不存在顯著自相關性，具體結果見圖4。

a) 韓家峁

b) 橫山

c) 趙石窯

d) 白家川

e) 殿市

f) 馬湖峪圖4 去趨勢年徑流序列的PACF檢驗

g) 青陽岔

h) 李家河

i) 曹坪

3.2 高階自相關的物理解析

從上節(jié)的高階自相關識別結果來看，不同的流域檢驗的結果差異較大。漢江流域和渭河流域研究站點均未檢測到高階自相關，而位于陜北的無定河流域則在多站點出現(xiàn)了高階自相關成分。一般高階自相關性主要受水循環(huán)因子的周期變化以及地下水入滲到排泄的時差的影響。由于并不是3個流域的水文序列都表現(xiàn)為高階自相關，因此，本節(jié)主要從地下水文循環(huán)的物理機制角度，分析產(chǎn)生這種空間差異的原因。

j) 丁家溝

k) 綏德續(xù)圖4 去趨勢年徑流序列的PACF檢驗

3.2.1漢江流域

在漢江流域選取的石泉和安康2個站點均位于安康市，區(qū)域地貌單元為安康盆地漢階地地層，其中淺層含水巖組的含水層巖性為粗砂、礫石、卵石等，深度在18～33 m以下，含水巖組受大氣降水及區(qū)域地下水補給，承壓性一般，水量較豐富；淺層承壓水含水巖組的含水層巖性為粗砂礫石、卵礫石、漂卵石等，深度在70～80 m以下，含水巖組受地下水補給，承壓性好，水量較豐富[21]。同時該區(qū)域雨量充沛，地表水豐富，為地下水形成提供了良好的補給條件，但地形破碎，河溝縱橫，儲水條件差，不利于地下水貯存[27]。因此，地下水補給地表徑流量的部分較少，且由于地處山區(qū)補給過程所需的時間跨度也較短，故徑流序列中產(chǎn)生高階自相關性的物理機制并不顯著存在，所以未能在序列中檢測到高階自相關性。

3.2.2渭河流域

渭河流域干流上、中、下游河水與地下水轉(zhuǎn)化關系有所不同。在渭河干流，上游地區(qū)地下水補給河水；中游地區(qū)河水一側(cè)補給地下水、另一側(cè)地下水補給河水，豐、枯季節(jié)地下水與河水的轉(zhuǎn)化關系也時常發(fā)生變化；下游地區(qū)如渭南至潼關段，河床淤積抬升，一般河水補給地下水[15]。但由于渭河徑流量中地下水補給量有限，加之渭河流域地處關中平原區(qū)，受當?shù)厮牡刭|(zhì)條件影響，流域地下水循環(huán)較快，從入滲到排泄至河流的時滯較短，因此，其反映到年際尺度的水文序列變化上并不顯著。同時，近50 a來渭河流域降水量呈顯著減少趨勢[26]，導致以降水為地下水主要補給源的渭河流域的地下水蓄水量補充有所削減，加之區(qū)內(nèi)大量開采地下水，導致渭河基流量顯著減少，河川徑流中原有的地下水循環(huán)所蘊含的周期規(guī)律和相關關系在愈發(fā)弱化。由此可見，多因素的綜合影響導致了渭河流域各水文站的年徑流序列中未檢測到高階自相關性。

3.2.3無定河流域

無定河流域與漢江、渭河流域有較大不同，其河川徑流中30%～80%來自地下水補給[29]，因此，地下水循環(huán)中的典型規(guī)律會較為顯著地體現(xiàn)在河川徑流過程中。

從地貌類型來看，無定河流域風沙區(qū)的韓家峁、橫山、趙石窯序列和黃土丘陵溝壑梁峁區(qū)的白家川序列檢驗出高階自相關性；其余站序列均未檢驗出高階自相關性。

無定河流域風沙區(qū)地勢平坦，林草覆蓋較好，土質(zhì)為顆粒較大的沙土，地面滲漏強烈，下滲能力強[22]，地下水埋深較淺，降水可補給到淺層地下水，而后再經(jīng)由地下水徑流途徑，轉(zhuǎn)化為基流匯入河道。風沙區(qū)地面滲漏強烈，地下水補給徑流所占比重較大，一般達80%～90%以上[29]，因此地下水循環(huán)特征可得到較好的表征。而黃土區(qū)植被條件差，地形起伏大、溝谷深切，黃土區(qū)地面下3.5 m處的恒濕層阻斷了地表水與地下水之間的水力聯(lián)系[25]，使得地下水難以得到降水的有效補給，滲入土壤的水分通常無法到達飽和含水層，而以壤中流的方式補給到河流，因此在黃土區(qū)地下水循環(huán)路徑一般較短，短時滯在年際尺度的河川徑流變化中無法得到有效體現(xiàn)。黃土區(qū)地下水補給只占年徑流的30%左右[28-29]。白家川水文站是無定河流域的出口水文站，年徑流量遠大于其余水文站，且基徑比達71%[27]；由于黃土區(qū)地下水補給量相對較小，因此其更趨向于表現(xiàn)上游主導地下水過程的相關或耦合特征。因此，上述因素可較為清楚地解釋風沙區(qū)年徑流序列呈現(xiàn)出高階自相關性，而黃土區(qū)的年徑流序列不存在高階自相關性的現(xiàn)象。

除此之外，筆者在研究中發(fā)現(xiàn)，出現(xiàn)高階自相關的站點所控制的流域都存在與其他流域存在側(cè)向交換的情況，即在某一含水層與其他流域存在地下水連通或補給關系。因此筆者認為無定河流域存在水文序列的高階自相關性可能與流域地下水分水嶺不閉合存在一定的關系。但限于目前無相關成果，無法佐證。

3.3 驗證分析

上文對韓家峁、橫山、趙石窯、白家川4個站點顯著高階自相關性的物理成因進行了簡單的定性分析，物理機制定量論證受觀測資料限制目前尚無法開展。本文擬以BP神經(jīng)網(wǎng)絡為基礎構建降水-徑流模型，通過分析加入高階降雨序列對徑流預測效果的影響，從側(cè)面論證高階自相關的客觀存在性以及影響。

以橫山站實測年徑流和降雨序列(1957—2010年)為例，檢測結果顯示去趨勢年徑流序列存在7階顯著高階自相關性，同時6階偏自相關系數(shù)較大，不可忽視，故選取當年降雨序列P0、1—7階降雨序列(表示為P1…P7)，當年去趨勢年徑流序列R0、進行BP神經(jīng)網(wǎng)絡計算分析，模型的構建過程見圖5(P7+P0表示將兩組序列進行組合，其余類似)。

橫山站設置序列的訓練期為1964—1998年，輸入變量有5組，輸出變量有一組，將訓練期的輸入和輸出變量分別代入BP神經(jīng)網(wǎng)絡，構建多輸入單輸出的BP神經(jīng)網(wǎng)絡模型并得到5組訓練期預測值；設置序列的驗證期為1999—2010年，以構建的BP神經(jīng)網(wǎng)絡模型為基礎，將5組輸入變量代入其中，得到5組驗證期預測值，并將其與5組訓練期預測值進行組合，得到的5組預測序列分別記為R1—R5，其與當年去趨勢年徑流序列R0的對比，見圖6。

圖5 BP神經(jīng)網(wǎng)絡模型的構建過程(以橫山站為例)

由圖6和表2可以看出，當輸入變量依次為P0、P7+P0、P6+P7+P0、P1+P6+P7+P0即加入了不同階數(shù)的降雨序列時，預測序列R1—R4的預測效果逐漸變好，與R0的擬合度R2由0.336逐漸增加到0.769，表明當把不同階數(shù)的降雨序列加入當年降水序列中進行預測時能有效提高預測的效果，從側(cè)面驗證了徑流序列高階自相關的存在。但當加入1—6階的全部降雨序列時，預測序列R5的預測效果反而有所下降，與R0的擬合度降低至0.375。由此可見，在加入不同階數(shù)計算時應該首要考慮的是顯著階數(shù)、偏自相關系數(shù)較大值對應的階數(shù)，而把所有階數(shù)均加入計算時則有可能因為各階序列的重復疊加，使預測出現(xiàn)疊加誤差，進而導致預測效果有所下降。

a)對比R1

b) 對比R2

c) 對比R3

d) 對比R4

表2 橫山站5組預測序列與R0的擬合度

上文中所述的5組預測序列和R0均是不包含趨勢成分的殘余序列，并不存在實際的物理意義，遂將分段去趨勢步驟中所剔除的趨勢成分加回到預測序列與R0中，即得到重組序列RR1—RR5和實測年徑流序列R。此時重組序列則對應于實測年徑流序列的預測序列，重組序列與實測年徑流序列的相關關系，見圖7。與表2比較可知，將趨勢成分加回后，5組重組序列與實測年徑流序列的擬合度均有較大程度的提高。

a) RR1

b) RR2

c) RR3

d) RR4

e) RR5圖7 橫山站重組序列與實測年徑流序列的相關關系

將同樣的方法應用于韓家峁、白家川和趙石窯站，在構建BP神經(jīng)網(wǎng)絡模型時所涉及的變量見表3。 計算得到3個站的預測序列與當年去趨勢年徑流序列、重組序列與實測年徑流序列的擬合度見表4。

由表4可知，當把顯著階數(shù)的降雨序列納入預報因子并實施預測時，可有效提高相關的徑流序列擬合度的大??；但當加入過多階降雨序列時，擬合度反而有所下降?；谝陨戏治隹芍?，徑流序列中的高階自相關性是序列蘊含的規(guī)律性成分，在預測模型中考慮該特征，可有效提高徑流預測的效果。由此，在僅有短期實測徑流資料情況下，需要展延年徑流資料時，可利用序列的顯著低階與高階自相關性，優(yōu)化輸入變量，構建多輸入單輸出的BP神經(jīng)網(wǎng)絡模型，可在一定程度上有效提高預測的效果，并不失為一種提高徑流預測效果的新思路。

表3 BP神經(jīng)網(wǎng)絡模型中的變量

表4 韓家峁、白家川、趙石窯各序列擬合度R2對比

4 結論

在以往對水文時間序列自相關性的研究中，多聚焦于序列的一階自相關性，認為一階自相關成分剔除后，高階自相關即隨之消除或不再顯著。然而目前已有高階自相關的相關研究成果表明，其影響客觀存在且不容忽視。由于國內(nèi)外類似研究較少，水文序列中的高階自相關性對水文分析是否有影響、有何種影響尚不明確，故本文對上述問題做了一些有益的嘗試與分析。具體為使用分段去趨勢方法對不同流域?qū)崪y年徑流序列進行處理，并對去趨勢序列的高階自相關性進行識別；以BP神經(jīng)網(wǎng)絡為基礎，降水-徑流模型為載體，探討了加入低階或(與)高階降雨序列對徑流預測效果的影響，得到了以下結論。

a) 高階自相關性不同于一階自相關性，其并不是一個普遍問題。不同流域徑流序列中高階自相關性的存在與否以及結構特征均有所差異。如本次研究中，漢江流域的石泉和安康站的實測年徑流序列無高階自相關性；渭河流域7個水文站的實測年徑流序列亦無高階自相關性；而無定河流域的韓家峁、橫山、趙石窯、白家川4個水文站的實測年徑流序列被檢測到高階自相關成分，其余7站則未發(fā)現(xiàn)高階自相關性。

b) 徑流高階自相關性是否存在，受地貌類型、土壤特性、降水規(guī)律、地下水循環(huán)條件等多種因素的共同影響。被檢測到高階自相關性的無定河流域的韓家峁、橫山、趙石窯水文站具有類似的影響特征，如：①均位于風沙區(qū)；②地下水補給徑流所占比重較大，一般達80%～90%以上；③陜北地區(qū)的包氣帶厚度一般較大，地下水由入滲到補給河流這一過程所需的時間較長，因此在實測徑流序列中呈現(xiàn)出了高階自相關性。漢江與渭河流域河川徑流主要來源于地表水，地下水基流比重較小，因而相關特征在水文站的實測年徑流序列表現(xiàn)不強烈，即不能檢測到顯著高階自相關性。

c) 水文序列中存在高階自相關性可能會影響統(tǒng)計檢測結果的準確性，進而得到不科學的結論。但高價自相關性并非只有劣勢，在徑流預測中如考慮水文序列高階自相關性影響，實施預測模型輸入變量優(yōu)化，其可有效提高徑流預測預報的精度和效果。除此之外，在僅有短期實測徑流資料的情況下，高階自相關性的有效利用，也可為展延年徑流資料提供一種新思路。