李鎮(zhèn)洋，李國芳，高倩雨，尹俊凱

(河海大學(xué)水文水資源學(xué)院，江蘇 南京 210098)

隨著經(jīng)濟(jì)社會的發(fā)展，形成徑流的環(huán)境在不斷發(fā)生變化。一方面，溫室氣體排放增多使溫室效應(yīng)加劇，影響了降雨；另一方面，人類大規(guī)模改變下墊面，不斷興建水庫、堤防等工程，影響了產(chǎn)流和匯流。變化的環(huán)境已經(jīng)影響到水文循環(huán)過程，破壞了部分地區(qū)徑流序列的一致性[1]，使水旱災(zāi)害防御和水資源開發(fā)利用面臨諸多新問題[2]。因此，研究變化環(huán)境對徑流的影響，對區(qū)域水資源管理、水問題治理，乃至新時(shí)期水文循環(huán)的研究具有重要意義。

在人類活動(dòng)較頻繁的地區(qū)，相比于氣候變化，短時(shí)間尺度上土地利用變化是區(qū)域水資源變化的主要驅(qū)動(dòng)因素[3]。傳統(tǒng)對比流域試驗(yàn)方法是一種可靠有效的方法，通過對比相似流域受不同人類活動(dòng)影響后的變化反映下墊面變化的水文響應(yīng)，但其著重于單一因子變化對徑流的作用，且不同時(shí)期的變量不易控制[4]；隨著水文模型的發(fā)展，基于物理過程模擬的分布式水文模型應(yīng)用越來越廣泛，如郭軍庭等[5]將SWAT模型應(yīng)用于潮河流域，榮艷淑等[6]將GAM模型應(yīng)用于洪家渡地區(qū)，但分布式模型的參數(shù)通常很多且難以測定和校準(zhǔn)[7]。

水利工程的興建在很大程度上改變了徑流原有的自然規(guī)律[8]，因而水利工程對徑流的影響也備受關(guān)注。Richter等[9-11]提出的水文變異指標(biāo)法(IHA)及變化范圍法(RVA)均可有效評價(jià)水利工程對河流的影響，如陳偉東等[12]采用IHA法分析錦江水庫對徑流的影響，黃草等[13]運(yùn)用IHA-RVA法分析水庫群對拉薩河徑流的影響；統(tǒng)計(jì)方法則是檢測大流域長期水文資料變化趨勢的常用方法之一[14]，如徐成漢[15]基于傾向分析和趨勢分析等方法評價(jià)三峽水庫蓄水前后烏江流域徑流的變化，趙萍等[16]通過離散系數(shù)和年內(nèi)不均勻系數(shù)等指標(biāo)評價(jià)卡普恰蓋水庫對伊犁河徑流的影響。但多數(shù)方法無法揭示水庫對流域水文循環(huán)的作用機(jī)制，也不能有效分離出除水庫外其他因素對徑流的作用。

綜上，土地利用變化和水利工程運(yùn)行綜合影響下區(qū)域徑流響應(yīng)研究是當(dāng)前熱點(diǎn)問題，但相關(guān)研究成果尚顯不足，對嫩江流域徑流演變的研究更是稀缺。因此，本文選取嫩江干流富拉爾基斷面為典型，分別采用統(tǒng)計(jì)指標(biāo)和SWAT模型，分析典型斷面以上土地利用變化和水庫運(yùn)行對流域徑流變化的影響，以期為嫩江流域水旱災(zāi)害防治和水資源開發(fā)利用提供參考依據(jù)。

1 研究區(qū)概況及數(shù)據(jù)資料

1.1 研究區(qū)概況

嫩江富拉爾基斷面以上流域位于我國東北部，流域面積12.39萬km2，其中尼爾基水庫以上地區(qū)主要為山區(qū)，尼爾基至富拉爾基斷面區(qū)間為山區(qū)與平原過渡段，研究區(qū)地形及水系如圖1所示。流域受溫帶大陸性季風(fēng)氣候影響，冬季干冷，夏季濕熱，年平均氣溫為2～4℃[17]，年降水量約500 mm，降水主要集中在夏季。

圖1 研究區(qū)示意圖

尼爾基水庫壩址以上集水面積6.64萬km2，總庫容86.1億m3，于2001年6月正式開工，2006年12月底主體工程全部完工，是以防洪、城鎮(zhèn)生活和工農(nóng)業(yè)供水為主的大型控制性工程。

1.2 數(shù)據(jù)資料

DEM數(shù)據(jù)來自中國科學(xué)院計(jì)算機(jī)網(wǎng)絡(luò)信息中心地理空間數(shù)據(jù)云平臺，分辨率為90 m。土地利用數(shù)據(jù)來自中國科學(xué)院資源環(huán)境科學(xué)數(shù)據(jù)中心，比例尺為1∶10萬。土壤數(shù)據(jù)來自國家青藏高原科學(xué)數(shù)據(jù)中心，比例尺為1∶100萬。氣象數(shù)據(jù)來自中國氣象數(shù)據(jù)網(wǎng)，時(shí)間為1981—2015年。水文數(shù)據(jù)選用富拉爾基站與氣象數(shù)據(jù)同時(shí)期的實(shí)測月徑流資料。

2 研究方法

2.1 研究思路

本文采用統(tǒng)計(jì)分析和模型模擬兩種方法，研究土地利用變化和水庫運(yùn)行對流域水文情勢的影響，分析流域徑流變化的趨勢。研究思路如下：①采用統(tǒng)計(jì)分析方法，評價(jià)1990—2010年期間流域土地利用變化，以及尼爾基建庫前后下游徑流變化；②驗(yàn)證SWAT模型對流域徑流模擬的適用性，并設(shè)置不同情景建立模型進(jìn)行對比，定量分析土地利用變化和水庫運(yùn)行對流域徑流變化的影響；③利用CA-Markov模型預(yù)測流域未來土地利用類型，基于預(yù)測結(jié)果設(shè)置情景，分析土地利用變化和水庫運(yùn)行綜合影響下流域徑流的變化趨勢。

2.2 土地利用變化分析方法

2.2.1土地利用動(dòng)態(tài)度

該指標(biāo)通過反映區(qū)域土地利用方式變化的速度和幅度，來描述土地利用類型變化的劇烈程度。

單一土地利用動(dòng)態(tài)度[18]的表達(dá)式為

(1)

式中：K為研究時(shí)段內(nèi)某一土地利用類型的動(dòng)態(tài)度；Ua、Ub分別為研究時(shí)段初和時(shí)段末某一土地利用類型的面積；T為時(shí)段長。

綜合土地利用動(dòng)態(tài)度[19]的表達(dá)式為

(2)

式中：LC為研究時(shí)段內(nèi)綜合土地利用動(dòng)態(tài)度；LUi為研究時(shí)段初第i類土地利用類型的面積；ΔLUij為研究時(shí)段內(nèi)第i類土地利用類型轉(zhuǎn)換為第j類土地利用類型的面積絕對值。

2.2.2土地利用轉(zhuǎn)移矩陣

該指標(biāo)可反映區(qū)域內(nèi)轉(zhuǎn)入和轉(zhuǎn)出的主要土地利用類型和土地利用類型間的轉(zhuǎn)化趨勢，表達(dá)式為

(3)

式中：P為土地利用轉(zhuǎn)移矩陣；Pij為第i類土地利用轉(zhuǎn)化為第j類土地利用的面積絕對值或占土地總面積的百分比。

2.3 SWAT模型

SWAT(soil and water assessment tool)模型是大中尺度的分布式水文模型，其在徑流模擬、泥沙模擬等方面得到廣泛應(yīng)用。模型基于物理機(jī)制，將流域的水文過程分為陸面部分和水面部分。流域水循環(huán)過程遵循水量平衡原理[20]，表達(dá)式為

式中：Wt為第i天的土壤最終含水量；W0為第i天的土壤初始含水量；Pday為第i天的降水量；Rsurf為第i天的地表徑流量；Ea為第i天的蒸散發(fā)量；Rseep為第i天土壤剖面地層中的側(cè)流總量；Rgw為第i天的地下徑流總量。

SWAT模型基于數(shù)字高程模型，根據(jù)控制點(diǎn)先將流域劃分為子流域，再結(jié)合土地利用和土壤類型來劃分水文響應(yīng)單元(HRUs)，最后對HRUs作匯流演算，在流域出口斷面得到總徑流。

2.4 土地利用類型變化預(yù)測方法

CA-Markov模型可用于預(yù)測流域未來土地利用類型。該模型由Markov模型和CA模型兩部分組成：Markov模型利用概率預(yù)測未來趨勢，可作為自發(fā)性下流域土地利用類型變化的預(yù)測模型；CA模型具有時(shí)空耦合特征，可用來預(yù)測自組織下流域土地利用類型變化的趨勢。

Markov模型的表達(dá)式為

Xk+1=PXk

(5)

式中：Xk、Xk+1分別為趨勢分析與預(yù)測對象在T=k和k+1時(shí)刻的狀態(tài)向量；P為土地利用轉(zhuǎn)移矩陣。

CA模型的表達(dá)式為

St+1=f(St,N)

(6)

式中：S為元胞的狀態(tài)集合；t和t+1為元胞所處的前后時(shí)刻；f為元胞轉(zhuǎn)化規(guī)則；N為元胞鄰域。

CA-Markov模型利用土地利用數(shù)據(jù)計(jì)算Markov轉(zhuǎn)移概率矩陣，并建立轉(zhuǎn)變適宜性圖像集，在確定CA濾波器和循環(huán)次數(shù)后，可得到預(yù)測的流域未來土地利用類型圖。

3 土地利用變化和水庫運(yùn)行對下游徑流影響分析

3.1 統(tǒng)計(jì)分析方法

3.1.1土地利用變化分析

研究區(qū)1990—2010年間的土地利用動(dòng)態(tài)度如表1所示，從中可以看出：水域和耕地的單一土地利用動(dòng)態(tài)度最大，分別為1.27%和1.20%；但流域的綜合土地利用動(dòng)態(tài)度為0.25%。說明1990—2010年間流域土地利用類型的變化幅度并不大。

表1 1990—2010年土地利用動(dòng)態(tài)度

1990—2010年間土地利用轉(zhuǎn)移矩陣如表2所示，從中可以看出：土地利用的轉(zhuǎn)移趨勢主要是由林地、草地和未利用土地向耕地、水域和建設(shè)用地轉(zhuǎn)化。但1990—2010年期間未變化的土地利用面積占全流域面積的84.9%，說明流域土地利用類型間的轉(zhuǎn)化并不明顯。

表2 1990—2010年研究區(qū)土地利用轉(zhuǎn)移矩陣 單位：km2

綜上，可以認(rèn)為流域的土地利用變化并不明顯，對徑流變化的影響較小。

3.1.2水庫對下游徑流影響分析

表3 尼爾基建庫前后富拉爾基站降水量及徑流量對比

圖2 尼爾基建庫前后富拉爾基站月平均流量對比

研究區(qū)水文站位置如圖1所示，本文將1983—2006年作為建庫前、2007—2015年作為建庫后，利用流域內(nèi)氣象站和富拉爾基水文站的實(shí)測資料，統(tǒng)計(jì)尼爾基建庫前后流域降水量和徑流量的變化，結(jié)果如表3和圖2所示。從中可以看出：水庫建成前后，流域的汛期降水量、非汛期降水量及多年平均降水量均差別不大；而水庫建成后，水庫下游多年平均徑流量有較大程度的減少，其中汛期(6月至9月)洪峰流量削減明顯，非汛期(10月至次年5月)枯季徑流顯著增加。據(jù)了解，這是由于尼爾基水庫作為北水南調(diào)的水源地，是北部引嫩總干渠的源頭、黑龍江省引嫩骨干工程的重要組成部分，年引水量超過23億m3[21]；同時(shí)，作為具有防洪、供水等多種功能的大型控制性工程，尼爾基水庫也在年內(nèi)發(fā)揮削減洪峰和調(diào)節(jié)徑流的作用。

與建庫前相比，建庫后流域降水量變化不大但下游徑流變化明顯，可以認(rèn)為尼爾基水庫運(yùn)行對下游徑流變化產(chǎn)生了較大的影響。

3.2 模型分析方法

3.2.1SWAT模型適用性分析

選取水庫建成前的1981—2000年資料做SWAT模型適用性研究，其中，1981—1982年為預(yù)熱期，1983—1994年為率定期，1995—2000年為驗(yàn)證期。進(jìn)行模型模擬效果評價(jià)時(shí)，采用納什效率系數(shù)NSE和決定系數(shù)R2。根據(jù)相關(guān)研究與應(yīng)用實(shí)例[22]，當(dāng)NSE>0.5、R2>0.6時(shí)，模型模擬的結(jié)果可被接受。

富拉爾基站1981—2000年的月徑流模擬結(jié)果如圖3和圖4所示，從中可以看出，實(shí)測值在每年主汛期前均有一個(gè)小峰，分析是春汛，由于模型模擬過程中未考慮流域內(nèi)冬季固態(tài)降水及春季融雪徑流，所以模擬值中未出現(xiàn)春汛。率定期NSE=0.81、R2=0.82，驗(yàn)證期NSE=0.84、R2=0.85，均達(dá)到了精度的要求，故認(rèn)為SWAT模型適用于嫩江富拉爾基斷面以上流域的徑流模擬。

圖3 富拉爾基站率定期(1983—1994年)月徑流模擬值與實(shí)測值對比

圖4 富拉爾基站驗(yàn)證期(1995—2000年)月徑流模擬值與實(shí)測值對比

3.2.2土地利用變化和水庫運(yùn)行影響模擬結(jié)果分析

為定量評估土地利用變化和尼爾基水庫運(yùn)行對徑流變化的影響，設(shè)置了3種情景。情景1:應(yīng)用原SWAT模型；情景2：應(yīng)用只考慮土地利用變化的SWAT模型；情景3：應(yīng)用考慮土地利用變化和尼爾基水庫運(yùn)行的SWAT模型。選取水庫建成后的2008—2015年資料進(jìn)行模擬研究，其中2008—2009年為預(yù)熱期，結(jié)果情景1NSE=0.71、R2=0.80，情景2NSE=0.71、R2=0.81，情景3NSE=0.93、R2=0.94。以2014年為例，3種情景下模擬流量與實(shí)測流量的對比如圖5所示。

圖5 富拉爾基站(2014年)月徑流模擬值與實(shí)測值對比

情景2和情景1相比，NSE沒有變化，R2從0.80提高到0.81，且兩者的模擬結(jié)果均與實(shí)測徑流過程差距明顯；兩者的月徑流模擬值十分相近，徑流量主要集中在汛期，非汛期徑流量很小，年內(nèi)分配極不均勻?？梢?，僅考慮土地利用變化以使得模擬值接近實(shí)測值的作用有限，土地利用變化對流域徑流變化的影響并不大。

情景3和情景2相比，NSE和R2分別從0.71和0.81提高到了0.93和0.94，且情景3的模擬結(jié)果與實(shí)測徑流過程十分接近；考慮水庫運(yùn)行后，水庫發(fā)揮削峰滯洪和調(diào)節(jié)徑流的作用，月徑流模擬值在非汛期有所增大，在汛期普遍減小，徑流量在年內(nèi)的分配更為均勻。可見，考慮水庫運(yùn)行以使得模擬值接近實(shí)測值的作用顯著，水庫運(yùn)行對流域徑流變化的影響較大。

綜上，選取R2作為定量分析指標(biāo)，計(jì)算所得到的流域徑流變化，7.14%是土地利用變化引起的，92.86%是尼爾基水庫運(yùn)行引起的。這與孫新國等[23]在研究豐滿Ⅱ區(qū)子流域時(shí)得出的結(jié)論一致：水利工程運(yùn)行對徑流變化占主導(dǎo)作用。

4 不同土地利用和水庫運(yùn)行情景下的徑流變化

4.1 流域未來土地利用類型預(yù)測

為進(jìn)一步研究變化環(huán)境下流域徑流的變化趨勢，利用CA-Markov模型模擬預(yù)測流域未來土地利用類型。先利用1990年和2000年的土地利用類型圖預(yù)測得到2010年的土地利用類型圖，并與實(shí)際圖進(jìn)行對比，利用Kappa系數(shù)衡量模型預(yù)測的準(zhǔn)確性。通常，當(dāng)Kappa大于0.6時(shí)，認(rèn)為具有高度的一致性。經(jīng)計(jì)算，2010年土地利用類型模擬結(jié)果的Kappa系數(shù)為0.921 7，即模擬結(jié)果與實(shí)際結(jié)果相近，表明CA-Markov模型可用于研究區(qū)土地利用類型的預(yù)測。因此基于該模型，利用2000年和2010年的土地利用類型圖，預(yù)測得到2020年的土地利用類型圖，如圖6所示。

圖6 研究區(qū)2020年土地利用模擬圖

4.2 情景建立與模擬結(jié)果分析

為研究在相同條件下流域徑流受土地利用變化和水庫運(yùn)行綜合影響的變化趨勢，本文在1990年、2000年、2010年和2020年4種土地利用方式下設(shè)置有/無水庫運(yùn)行情景，采用2008—2015年的氣象資料，選取流域出口處的富拉爾基站作為監(jiān)測站，基于SWAT模型進(jìn)行徑流模擬，結(jié)果如表4和圖7、圖8所示。

由表4和圖7、圖8，可以得出以下結(jié)論：①在相同的降水和水庫運(yùn)行條件下，隨著土地利用的變化，研究區(qū)年均徑流量呈先下降后上升趨勢。其中，1990—2010年為下降趨勢，2010—2020年為上升趨勢。分析原因可能是，1990—2010年流域的土地利用變化趨勢是林地、草地轉(zhuǎn)化為耕地，在國家政策的影響下2010—2020年土地利用變化趨勢是耕地、林地轉(zhuǎn)化為草地、水域；而耕地、林地面積和徑流量是負(fù)相關(guān)，草地面積和徑流量是正相關(guān)[24]。②相比無水庫運(yùn)行的情景，有水庫運(yùn)行的情景下，流域多年平均月徑流年內(nèi)變化有所減小。其中，汛期的洪峰流量有所削減，非汛期的枯季徑流有所增加。說明水庫發(fā)揮了削峰滯洪和調(diào)節(jié)徑流的作用，對流域徑流的年內(nèi)分配有明顯的影響。

表4 富拉爾基站不同情景多年平均徑流量模擬結(jié)果

圖7 富拉爾基站1990年和2000年不同情景下月徑流模擬結(jié)果

圖8 富拉爾基站2010年和2020年不同情景下月徑流模擬結(jié)果

5 結(jié) 論

a.1990—2010年流域綜合土地利用動(dòng)態(tài)度為0.25%，土地利用類型未變化的面積為84.9%，流域整體的土地利用變化不明顯，對徑流變化影響較?。荒釥柣畮爝\(yùn)行后下游徑流量變化顯著，其中汛期洪峰流量削減，非汛期枯季徑流增加，對徑流變化影響較大。

b.利用1981—2000年資料驗(yàn)證了SWAT模型適用于研究區(qū)徑流模擬；利用2008—2015年資料定量評估了土地利用變化和水庫運(yùn)行對徑流變化的影響，水庫運(yùn)行占92.86%，土地利用變化占7.14%。相比土地利用變化，水庫運(yùn)行對流域徑流變化的影響占主導(dǎo)作用。

c.基于CA-Markov模型預(yù)測的土地利用，設(shè)置情景分析得到，在相同的降水和水庫運(yùn)行條件下，1990—2020年的土地利用變化具有使流域年均徑流量先減小后增大的效應(yīng)；相比土地利用變化，水庫運(yùn)行對徑流變化的影響更明顯，對其年內(nèi)分配有顯著影響。