環(huán)境變化對石棧道站徑流量的影響研究

孟鵬飛，任 政，石 彬，賈騰飛，王鈺娟

(河北工程大學 水利水電學院，河北 邯鄲 056038)

1 研究背景

濁漳河流域是我國水事糾紛沖突最頻繁的地區(qū)之一[1]，也是最易受環(huán)境變化影響的地區(qū)之一[2]。石棧道站是濁漳北源極其重要的徑流觀測站，至20世紀60年代以來，徑流呈現(xiàn)顯著的減少趨勢，給流域水資源管理帶來了巨大挑戰(zhàn)。有效的水資源管理對區(qū)域社會經(jīng)濟健康發(fā)展至關重要，未來區(qū)域水資源規(guī)劃需要全面了解環(huán)境變化對徑流過程的影響[3]，定量確定徑流變化的主導因素即歸因分析已成為制定未來有效水資源管理的前提之一。

目前，氣候變化和人類活動是影響徑流變化的主要環(huán)境因素。彈性系數(shù)法和水文模擬法為應用于定量評估環(huán)境變化對徑流影響的兩種主要方法。彈性系數(shù)法對原始數(shù)據(jù)的要求較低，能單獨評估各影響因素的影響值[4]，如PENG等[5]采用雙累積曲線法和彈性系數(shù)法研究太湖流域，結果表明徑流變化主要受人類活動的影響。然而，彈性系數(shù)法僅能定量分析年尺度上的環(huán)境變化的影響[6]。

水文模擬法基于物理分布的水文模型，可應用在月尺度或日尺度等多種尺度上[7]，該模型模擬結果精度較高，更能合理地模擬降雨徑流過程，還可與多種模型耦合嵌套使用[8]。ZENG等[9]采用SIMHYD模型對漳河流域不同時間尺度的徑流進行研究，結果表明氣候變化和人類活動對月徑流量變化的影響在不同月份是不同的。黃亞等[10]將 CMIP5和VIC模型嵌套，分析了紅水河上游流域未來徑流量的變化情況。

本研究以石棧道水文站為研究對象，以水文氣象資料為基礎，使用VIC水文模型從不同年代來定量評估環(huán)境變化對石棧道水文站年徑流的影響，以期為石棧道站洪水預報、濁漳河流域水資源規(guī)劃提供相關依據(jù)。

2 研究區(qū)域及數(shù)據(jù)來源

濁漳河是海河流域漳衛(wèi)河水系漳河上游的一條支流，分為北源、南源和西源。濁漳河北源發(fā)源于山西省榆社縣柳樹溝。石棧道站是濁漳北源上游重要水文站，在濁漳河上游石棧道村的榆社河上，位于112°47′～113°9′ E，37°4′～37°24′ N之間，其控制流域面積為702 km2，研究區(qū)域水系及氣象、水文站點分布如圖1所示。

圖1 研究區(qū)域水系及氣象、水文站點分布圖

氣象數(shù)據(jù)包括：(1)研究區(qū)內石棧道、新莊、社城、水澤、官寨、郜村、銀郊和雙峰8個降水觀測站，時段為1959-2015年的實測日數(shù)據(jù)，數(shù)據(jù)來源為《水文年鑒》；(2)研究區(qū)周邊邢臺、陽城、介休、平定、太原、安澤、襄垣、長治和榆社共9個國家氣象站，時段為1953-2015年的平均氣溫、最高氣溫、最低氣溫和平均風速的日數(shù)據(jù)，數(shù)據(jù)來源為《中國氣象數(shù)據(jù)共享服務網(wǎng)》。石棧道站日平均徑流數(shù)據(jù)來源為《水文年鑒》，時段為1959-1991年和2006-2015年。其他數(shù)據(jù)，如DEM數(shù)據(jù)來源于地理空間數(shù)據(jù)云，植被數(shù)據(jù)采用馬里蘭大學發(fā)展的全球1 km的土地覆蓋數(shù)據(jù)，土壤數(shù)據(jù)采用Harmonized World Soil Database。

3 研究方法

3.1 Mann-Kendall趨勢檢驗

給定水文時間序列X{xi,i=1,2,…,n}進行趨勢檢驗，定義統(tǒng)計量S的計算公式為[11]：

(1)

其中函數(shù)sgn(·)定義如下：

(2)

Mann(1945)和Kendall(1975)指出：統(tǒng)計量S近似服從正態(tài)分布，統(tǒng)計量的均值和方差分別為E(S)=0和Var(S)=n(n-1)(2n+5)/18，則水文時間序列標準化后的檢驗統(tǒng)計量Z定義為[11]

(3)

式中：n為水文時間序列的長度。

根據(jù)Mann-Kendall法原理，以統(tǒng)計量Z值判斷徑流序列變化趨勢[12]。統(tǒng)計量Z為正表示序列呈現(xiàn)增加趨勢，Z為負表示序列具有遞減趨勢[11]。若|Z|>1.96，序列變化趨勢通過5%顯著水平檢驗。

3.2 Mann-Kendall突變檢驗

對水文時間序列X{xi,i=1,2,…,n}進行突變檢驗，構建秩序列：

(4)

其中rij定義如下：

(5)

產(chǎn)生的秩序列Sk近似服從正態(tài)分布，Sk的均值E(Sk)和方差Var(Sk)的計算公式為：

(6)

(7)

定義統(tǒng)計變量：

(8)

(k=1,2,…,n),UF1=0

將所有的UFk值繪制成曲線UF。按時間順序逆序排序時間序列X，重復上述步驟，得到反序列的UBk，然后乘以(-1)，得到一條曲線UB。在同一個圖中，繪制這兩條曲線。在臨界值內，兩條曲線的交點為水文時間序列的突變點。

3.3 Pettitt突變檢驗

假設水文時間序列X在時刻t發(fā)生了突變，將時刻t作為序列X的分割點，將序列分為兩部分，定義統(tǒng)計量Ut,n[11]：

(9)

(t=2,3,4,…,n)

式中：sgn(·)為符號函數(shù)，見公式(2)所示，則統(tǒng)計量Ut,n的極小值點為序列的突變點。

突變點可能的顯著性水平P定義為:

(10)

若P≤0.05，則認為該點通過5%顯著水平檢驗。

3.4 構建基于VIC模型的徑流變化歸因分析模型

根據(jù)石棧道站年徑流序列，進行趨勢和突變分析，將徑流序列劃分為基準期和變化期。基于基準期的水文氣象資料，構建基于VIC徑流模擬模型，并率定和檢驗模型的各項參數(shù)，模型的各項參數(shù)見表1?；谧兓诘臍庀筚Y料，模擬石棧道站變化期逐月徑流過程，按照以下公式定量評估變化期環(huán)境變化對徑流變化的貢獻率[13]。

表1 VIC模型參數(shù)及取值范圍

ΔWT=WHR-WB

(11)

ΔWH=WHR-WHN

(12)

ΔWC=WHN-WB

(13)

ηH=ΔWH/ΔWT

(14)

ηC=ΔWC/ΔWT

(15)

式中：ΔWT為石棧道站徑流變化量,m3；WHR為石棧道站變化期多年平均徑流量,m3；WB為石棧道站基準期多年平均徑流量,m3；ΔWH為人類活動對徑流的貢獻量,m3；WHN為石棧道站變化期水文模型模擬的多年平均徑流量,m3；ΔWC為氣候變化對徑流的貢獻量,m3；ηH為人類活動對徑流的貢獻百分比；ηC為氣候變化對徑流的貢獻百分比。

4 結果與分析

4.1 徑流趨勢分析

根據(jù)石棧道站1959-1991年連續(xù)實測徑流資料，繪制年徑流與年份相關圖和5 a滑動平均年徑流與年份相關圖，如圖2所示。從圖2可以看出，無論是年徑流過程和5 a滑動平均年徑流過程，石棧道站徑流過程均呈現(xiàn)減少趨勢，多年平均遞減率為188×104m3/a。

圖2 石棧道水文站年徑流趨勢分析

為進一步驗證石棧道徑流變化趨勢，采用Mann-Kendall進行趨勢檢驗，見表2所示。

表2 石棧道站實測徑流量趨勢性檢驗表

從表2中可以看出，統(tǒng)計量Z<0，年徑流序列呈減少趨勢；且|Z|>1.96，年徑流序列變化趨勢通過5%顯著水平檢驗。因此，石棧道站1959-1991年徑流序列，年平均遞減率188×104m3/a，且遞減趨勢顯著，這與相關文獻分析結果是一致的[14]。

4.2 徑流基準期劃分

根據(jù)石棧道站1959-1991年連續(xù)實測徑流資料，分別采用Mann-Kendall和Pettitt方法進行突變檢驗，結果見圖3和4所示。

圖3 基于Mann-Kendall方法的石棧道水文站1959-1991年徑流量突變檢驗

由圖3可以看出，Mann-Kendall檢驗的UF、UB曲線的相交點位于1979-1980年，石棧道水文站年徑流量的突變點在1979年前后，突變檢驗結果見表3。由圖4可以看出，Pettitt檢驗值的極小值出現(xiàn)在1979年且P=0.0467(結果見表3)，石棧道水文站年徑流量的突變點為1979年，且突變通過5%顯著水平檢驗，突變顯著。

圖4 基于Pettitt方法的石棧道水文站1959-1991年徑流量突變檢驗

表3 石棧道站年徑流量時間序列突變檢驗結果

根據(jù)Mann-Kendall和Pettitt法對石棧道站1959-1991年徑流量序列的突變檢驗結果，綜合選擇1979年作為年徑流量序列的突變點。突變點將年徑流量時間序列劃分為徑流量基準期(1959-1979年)和變化期(1980-1991年、2006-2015年)，這與相關文獻分析結果是一致的[14]。

4.3 基于VIC模型的徑流量變化歸因分析

選取基準期的年徑流量序列，構建基于VIC模型，模擬石棧道站徑流量過程。在VIC模型參數(shù)率定過程中,考慮到模型的狀態(tài)變量對模擬值的影響,將1959-1962年定為模型的預熱期,1963-1972年定為模型的參數(shù)率定期,1973-1979年為模型參數(shù)的驗證期。

選取9個國家氣象站氣溫數(shù)據(jù)，率定研究區(qū)的日平均氣溫、日最高氣溫和日最低氣溫梯度，以研究區(qū)最近的榆社站氣溫為基礎，插值研究區(qū)1 km×1 km的網(wǎng)格日氣溫值，網(wǎng)格風速采用榆社站同期風速值，網(wǎng)格日降水為研究區(qū)8個降水觀測站的插值結果，構建石棧道站基準期徑流量過程模擬的VIC模型，調整VIC模型參數(shù)，具體參數(shù)值見表4所示。

表4 VIC模型各參數(shù)率定結果

對比模型模擬值和實測值，采用Nsc系數(shù)作為率定和驗證參數(shù)的目標函數(shù)，并計算其相對誤差Re。若Nsc系數(shù)越高，則意味著模擬值與實測值吻合度越強。1963-1979年石棧道站逐月徑流量實測值和模型模擬值對比結果如圖5所示。

從圖5中可以看到，VIC模型模擬過程與石棧道站實測徑流量過程擬合較好，洪峰出現(xiàn)的時間基本一致。VIC模型模擬的月徑流量結果中，率定期的Nsc系數(shù)和相對誤差Re分別為0.88和3.77%，驗證期的Nsc系數(shù)和相對誤差Re分別為0.92和-7.54%，其中相對誤差Re絕對值均控制在10%以內。分析年平均徑流量時，發(fā)現(xiàn)年徑流量的模擬值與實測值較為接近，率定期的Nsc系數(shù)和相對誤差Re分別為0.96和3.24%，驗證期的Nsc系數(shù)和相對誤差Re分別為0.81和-1.80%，其中相對誤差Re絕對值均控制在5%以內，具體統(tǒng)計結果見表5。

圖5 1963-1979年石棧道水文站流域逐月徑流量實測與模擬結果

表5 石棧道站VIC模型目標函數(shù)和相對誤差統(tǒng)計表

VIC模型目標函數(shù)和相對誤差均達到模型要求的0.75以上及10%以下。因此，VIC模型還原石棧道站的徑流量變化期過程較為合理。模型在研究區(qū)域內具有較高的適用性和可信度。

根據(jù)基準期構建的VIC模型，結合相關數(shù)據(jù)資料，模擬出石棧道站徑流量變化期的月徑流量數(shù)據(jù)。依據(jù)實測數(shù)據(jù)資料，并將徑流量變化期分為2個時段，即1980-1991年和2006-2015年。對比分析變化期實測徑流量和模型模擬徑流量結果，如圖6所示。評估石棧道站徑流量對環(huán)境變化的響應值，具體計算結果見表6。

圖6 石棧道水文站徑流量變化期2個時段逐月徑流量實測域模擬結果

表6 氣候變化與人類活動對流域徑流量變化的影響 mm，%

由表6中可知：(1)與基準期相比，1980-1991年實測降水量有所下降，2006-2015年實測降水量基本持平。而實測徑流量過程與降水的變化趨勢有所不同，實測徑流量在不同時段內均有所下降。(2)徑流量的變化量在兩個時段均為負值，這說明由人類活動引起的下墊面的變化減小了流域的產(chǎn)流能力。人類活動在這兩個時段內對區(qū)間產(chǎn)流減小的貢獻率分別為85.6%和84.6%。氣候變化在這兩個時段內對區(qū)間產(chǎn)流減小的貢獻率分別為14.4%和15.4%，貢獻率隨著區(qū)間降水量的增加而增大。(3)環(huán)境變化對徑流量的貢獻率在各個時段是不盡相同的，在1980-1991年和2006-2015年兩個時段內，均是人類活動占主導地位。

5 結 論

(1)石棧道站年徑流量呈現(xiàn)明顯下降趨勢，根據(jù)Mann-Kendall突變檢驗1979-1980年為徑流量突變點，Pettitt檢驗突變點為1979年，綜合判斷分析，確定突變年份為1979年。將1959-1979年作為石棧道徑流量時間序列的基準期。

(2)VIC模型模擬的基準期月徑流量結果中，率定期的NSC系數(shù)和相對誤差Re分別為0.88和3.77%，驗證期的Nsc系數(shù)和相對誤差Re分別為0.92和-7.54%。分析基準期的年平均徑流量，發(fā)現(xiàn)年徑流量的模擬值與實測值較為接近，率定期的Nsc系數(shù)和相對誤差Re分別為0.96和3.24%，驗證期的Nsc系數(shù)和相對誤差Re分別為0.81和-1.80%，均達到模型要求的0.75以上及10%以下的要求。因此，VIC模型還原石棧道站變化期的徑流量過程較為合理。模型在研究區(qū)域內具有較高的適用性和可信度。

(3)根據(jù)VIC模型模擬變化期的徑流量過程，人類活動和氣候變化的影響值在不同時段是不盡相同的，在1980-1991年和2006-2015年兩個時段內，濁漳河石棧道站徑流量變化受人類活動影響較大。人類活動對徑流量減少的影響分別為85.6%和84.6%，遠大于氣候變化14.4%和15.4%的影響。