趙雪巖，張 鑫，孫 媛

（西北農林科技大學 水利與建筑工程學院，陜西 楊凌 712100）

隨著社會發(fā)展，土地利用變化通過調節(jié)產匯流和降水蒸發(fā)過程來影響水文循環(huán)和水資源量的變化［1－3］。自1960 年以來，無定河流域受人類活動（如退耕還林還草、防風治沙、開荒種田等）的影響，其水文循環(huán)發(fā)生了顯著變化［4－7］。前人大多聚焦于無定河流域植被恢復前后水文通量的大幅度變化進行研究，如王計平等［4］認為1995 年前后無定河流域土地流轉頻繁；任宗萍等［5］認為水土治理措施使無定河支流大理河徑流量在1996 年出現(xiàn)拐點；張守紅等［6］認為水保措施是無定河流域徑流量減少的主要原因。然而，水保措施實施后徑流、生態(tài)的大幅變化趨勢是不可持續(xù)的，國家山水林田湖草沙等自然資源共同發(fā)展演化的趨勢是從大修大補的退耕還林（草）轉向小微型、質量型發(fā)展［8］。

SWAT（Soil and Water Assessment Tool）模型是一個基于物理機制的、復雜的半分布式水文模型［9－16］，該模型能夠模擬流域水文循環(huán)過程和土地利用變化過程等地表自然物理化學過程，并在國內外得到了廣泛應用［17－20］。Yang 等［9］認為灤河流域退耕還林使夏季流量增大，其他季節(jié)流量減??；傅春等［10］研究表明撫河流域林地增加、耕地減少使徑流量減少；楊倩等［11］研究表明漢江上游退耕還林使年徑流量減少；Ma 等［12］研究表明岷江流域耕地增加、草地減少使產水量減少；李爍陽等［13－16］在渭河流域和汾河上游均開展不同角度的研究。

本文以無定河流域為研究對象，運用SWAT 模型，利用五期土地利用數(shù)據(jù)（2005 年、2015 年、2020年、2025 年、2030 年，其中2025 年和2030 年為CAMarkov 模型預測結果）驅動模型，分析4 種土地利用變化及其對徑流的影響。

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

1.1 研究區(qū)概況

無定河流域位于北緯37°2′—39°1′、東經107°47′—110°34′（見圖1），多年平均徑流量為11.51 億m3。流域處在東部季風氣候和西北干旱氣候的過渡帶，多年平均氣溫為7.8～9.6 ℃，多年平均降水量為300～500 mm，降水量由東南向西北遞減，年際變化大，降水多集中在夏秋兩季，年內分配不均［4－7］。

圖1 無定河流域高程及水文站和氣象站分布

1.2 數(shù)據(jù)來源與處理

1.2.1 空間數(shù)據(jù)處理

空間數(shù)據(jù)介紹詳細見表1。其中，2005 年、2015年、2020 年三期土地利用數(shù)據(jù)用于搭建SWAT 模型的土地利用數(shù)據(jù)庫，原始數(shù)據(jù)中的二級分類需要重分類為草地、林地、耕地、未利用地、居民用地和水域6 種類型（見表2，其中土地利用分類屬性前的數(shù)字為代號）；2010 年土地利用數(shù)據(jù)用于構建CA－Markov 模型中的土地轉移概率矩陣。HWSD 用于搭建土壤屬性數(shù)據(jù)庫，同時還需輸入默認參數(shù)值，如上層鹽度默認0.1等；依托HWSD，利用SPAW 軟件進行計算后獲取其他參數(shù)。另外，石灰性雛形土和石灰性沖積土根據(jù)屬性的不同進一步細分為1 類和2 類。土壤類型和面積統(tǒng)計詳見表3 和圖2（空間參考坐標系統(tǒng)一設置為WGS 1984－utm－48n）。

表1 空間數(shù)據(jù)介紹

表2 土地利用數(shù)據(jù)重分類

表3 不同類型土壤的面積統(tǒng)計

圖2 無定河流域土壤分類

1.2.2 氣象水文數(shù)據(jù)處理

氣象數(shù)據(jù)包括逐日降水量、最高／最低氣溫、太陽輻射、風速、相對濕度等，來源于中國氣象數(shù)據(jù)網，氣象站分布見圖1。水文數(shù)據(jù)為逐日平均流量，來源于黃河流域水文年鑒（中游區(qū)上段），水文站分布見圖1。兩者時間段均為2006—2018 年。另外，各類氣象數(shù)據(jù)需要按照站編號、站名、經緯度、高程做成索引表（txt文件），缺測數(shù)據(jù)使用天氣發(fā)生器模擬進行插補。

2 研究方法

2.1 土地利用數(shù)據(jù)模擬

CA－Markov 模型結合了CA 模型模擬復雜系統(tǒng)空間變化的能力和Markov 模型的定量預測優(yōu)勢，可作為土地利用模擬工具，并在IDRISI 軟件中實現(xiàn)，其實質是建立土地轉移概率矩陣并利用卡帕系數(shù)Kappa評估概率矩陣模擬出的土地利用數(shù)據(jù)和實測數(shù)據(jù)的差異，從而建立模型用于輸出預測數(shù)據(jù)［21－23］。

式中：PO為恰當被模擬的部分；PC為由偶然因素導致的被希望糾正的部分；PP為觀測和模擬完全匹配的理想比例，取為1；Lij為i類土地利用類型轉變?yōu)閖類土地利用類型的概率；n為土地利用類型總數(shù)。

如果Kappa ＝1，說明二者是一致的；如果0.75 ≤Kappa ＜1，說明高度一致；如果0.5 ≤Kappa ＜0.75，說明中度一致；如果Kappa＜0.5，說明不一致［21］。

本文輸出2025 年和2030 年的無定河流域土地利用數(shù)據(jù)的主要步驟［21－22］：①利用Markov 模型建立時間間隔為5 a（2005—2010 年）的過渡概率矩陣；②設定預測的起始年份為2010 年，并利用CA 模型設置循環(huán)迭代次數(shù)為5 次，分別用來模擬2015 年土地利用數(shù)據(jù)。③模擬2015 年土地利用數(shù)據(jù)與實際2015 年土地利用數(shù)據(jù)進行比較，以Kappa來評估模擬效果。

2.2 SWAT 模型

在SWAT 模型中，一個流域被劃分為由河流網絡相互連接的許多子流域，每個子流域根據(jù)土地利用類型、土壤類型和坡度進一步劃分為水文響應單元（HRUs），模擬水的流動軌跡從HRUs 到子流域再到流域出口斷面。水文過程基于如下水量平衡方程［18－20］：

式中：SWt為最終土壤含水量，mm；SW為初始土壤含水量，mm；R、Qsurf、ET、P、QR分別為滲漏量、徑流量、蒸散發(fā)量、日降雨量、地下水含量，mm。

常用Nash－Suttcliffe 系數(shù)（Ens）、相關系數(shù)（r）［18］和百分比偏差（PBIAS）［24］來評估SWAT 模型的模擬效果。

式中：Qoi為i時的觀測徑流量；Qsi為i時的模擬徑流量；為觀測徑流量平均值；為模擬徑流量平均值；n為觀測數(shù)據(jù)總數(shù)［24］。

當r ＝1 時，結果非常吻合；r≥0.5 時，模擬結果可以接受。當Ens＝1 時，模擬結果和實測值完全吻合；0.5＜Ens＜1 時，模擬結果可以接受。當PBIAS ＝0 時，模擬結果非常好；PBIAS ＜0 時，模擬結果偏大；PBIAS ＞0 時，模擬結果偏小［24］。

年內分配完全調節(jié)系數(shù)Cr可表示徑流量年內分配不均勻程度，值越大不均勻程度越大［25］。

式中：t為月份；R（t）為t月的徑流量；為多年平均徑流量。

3 結果與分析

3.1 無定河流域土地利用變化分析

CA－Markov 模型運行結果表明：模擬的2015 年土地利用數(shù)據(jù)和實測數(shù)據(jù)的Kappa系數(shù)為0.930 7，說明二者高度一致。所以用2005—2010 年土地轉移概率矩陣和2010 年土地利用數(shù)據(jù)分別迭代15、20 次來模擬2025 年、2030 年土地利用數(shù)據(jù)。

本文用于土地利用分析的數(shù)據(jù)分辨率為1 km。由表4 可知：①無定河流域的土地利用類型主要是草地、耕地和未利用地，其中草地占流域總面積的40%左右，耕地占30%左右，未利用地占22%左右，其余為林地、水域和居民用地，共占8%左右；②草地面積先增加后減小，耕地先增加后減小而后保持維定，未利用地先減少后增加，居民用地持續(xù)擴張且偶有波動；③水域面積變化很小，說明流域25 a 內的水庫、河流、灘地基本保持不變。由圖3 可知：①草地主要集中在流域中上游地區(qū)，并在下游的耕地中有少量零星草地；②耕地主要集中在流域中下游；③未利用地大多分布在中上游地區(qū)，呈點狀分布，并且2030 年未利用地在東北部大塊狀出現(xiàn)，應引起關注；④居民用地增加主要集中在河道附近和流域東北部。

表4 不同時期土地類型面積統(tǒng)計

圖3 各土地利用類型分布情況

本文涉及4 種土地利用變化方式，見表5（記為C1、C2、C3、C4）。為探究C1～C4 的規(guī)律，建立4 個土地轉移矩陣（見表6～表9）。綜合圖4、表6～表9 可知：①居民用地雖然在流域內占比小，但是增長幅度最大，C1～C4 分別增長110%、220%、120%、210%；②C1～C4的轉移矩陣中，互相轉化的土地利用類型的面積遠小于沒有發(fā)生轉化的面積，說明2005—2030 年土地流轉僅在小范圍內進行；③C1，未利用地面積的1.2%和草地面積的0.6%共同轉向耕地，使耕地增加1.3%；④C2，4.7%未利用地和4.3%耕地轉向草地，使草地增加3%，同時1.27%的草地轉出為居民用地；⑤C3，未利用地減少10.4%，轉出為耕地、居民用地和草地；⑥C4，1.1%草地和2.1%未利用地轉向耕地，使耕地增加2.9%，與C1 相比變化幅度更大，同時4.3%的草地轉為耕地、未利用地和居民用地。

表5 不同土地利用變化方式

表6 C1 土地利用轉移矩陣 km2

表7 C2 土地利用轉移矩陣 km2

表8 C3 土地利用轉移矩陣 km2

表9 C4 土地利用轉移矩陣 km2

3.2 SWAT 模型校準和驗證

根據(jù)本文研究需要，并結合前人研究經驗［26－30］，采用1 km 的土壤數(shù)據(jù)集和土地利用數(shù)據(jù)集構建模型，共劃分18 個子流域（見圖5）和372 個HRUs。把2006—2013 年設定為模型率定期，把2014—2018 年設定為模型驗證期。2005 年土地利用數(shù)據(jù)用于模型初始構建，利用白家川水文站還原后的逐月平均流量（以下稱為還原值）進行校準。同時使用SWAT－CUP自動校準和驗證參數(shù)，見表10。由圖6 可知：①率定期和驗證期的多年月平均流量模擬值和還原值的變化趨勢有較好的一致性，擬合的線率分別為0.913、0.853 3；②率定期r、Ens、PBIAS分別為0.85、0.84、5.4%，驗證期r、Ens、PBIAS分別為0.89、0.88、5.0%。綜上，在無定河流域建立的SWAT 模型有良好適用性，可以用來研究C1～C4 對徑流的影響。

圖5 子流域分布

表10 參數(shù)率定結果

圖6 模型率定期和驗證期流量的模擬值和還原值

3.3 分析土地利用變化對徑流影響的試驗設計

氣象數(shù)據(jù)不變，分析C1～C4 對徑流的影響，具體試驗設計見表11。其中，定義Q1～Q5 為情景1～5 下得到的徑流模擬值的時間序列，R1～R5 為由流量時間序列得到的徑流深，ΔR1～ΔR4 為C1～C4 引起的徑流深變化。

表11 試驗設計

3.4 不同土地利用變化方式對徑流的影響分析

由圖7 可知：①ΔR1～ΔR4 的變化趨勢基本一致，且C1～C4 都會導致流域徑流量增加，徑流深多年總增加量從大到小排序為ΔR2＞ΔR4＞ΔR3＞ΔR1，說明C1～C4 均對流域徑流起促進作用，雖然C1～C4 中沒有進行流轉的土地占各自土地利用類型面積的88%以上（見圖4、表6～表9），且草地、耕地和未利用地面積各有增減，對徑流有不同程度的促進或抑制作用，但是居民用地面積均增加，說明居民用地面積增加對徑流的促進作用占主導地位，城市化程度高的流域，不透水地面的面積大，下滲量和下滲強度降低，屬于超滲產流，且居民用地主要分布在河道兩側，從而使匯流時間變短，徑流量增大，這與楊滿根［31］在淮河流域得出的結論一致。②ΔR2 多年變化總量達9.54 mm，共計2.74 億m3，占流域內現(xiàn)有水量的23.81%，說明C2 對徑流的促進作用最大。變化量最小的是ΔR1，多年變化總量達3.97 mm，C1 中居民用地面積增加1.1 倍，說明C1 對徑流的促進作用最小。③ΔR3 整體上略大于ΔR1，變化總量為4.40 mm；ΔR4 整體上略小于ΔR2，變化總量為8.40 mm，針對C4 使ΔR4 增加較大的情況，應引起關注。

圖7 土地利用變化引起的徑流深變化

為進一步探究ΔR的年內變化情況，本文從多年月平均徑流深和年內分配完全調節(jié)系數(shù)（Cr）兩方面入手。對R1～R5 中每年相同月份求平均值，而后分別和情景1 做差比較，得ΔR的年內變化，由圖8 可知：①ΔR1～ΔR3 中，多年月平均徑流深的年內增加趨勢一致，呈先減小后增大而后減小的趨勢。多年月平均徑流深的增大量排序為ΔR2＞ΔR3＞ΔR1＞ΔR4，ΔR2 的年內變化總量達0.92 mm，說明C1～C3 對ΔR的各月流量均起促進作用，且C2 對各月流量的促進作用最顯著，原因是C2 草地增加最多，由耕地和未利用地轉入，使地表粗糙度增大，這與孫鐵軍等［17］的結論一致。且經過多年水保措施改造后土壤含水量高［17］，雨水落在草地上形成超滲產流；居民用地面積增加最多，形成不透水面積最大［31］，兩方面原因共同促使各月徑流深增加量最多。②ΔR4 中，1—6 月徑流深減小總量共計0.26 mm，7—12 月徑流深增大總量共計0.68 mm，說明C4 對上半年徑流起抑制作用，對下半年徑流起促進作用，且促進程度比抑制程度更高，原因可能是C4 未利用地和草地向耕地轉移，使耕地面積大量增加，農作物在春天需要大量使用河水灌溉［7］，河水從河道進入土壤和地下，從而起到抑制徑流的作用；下半年降水量增加，灌溉用水量少，從而起到促進徑流的作用，針對草地和未利用地大量向耕地轉移，導致上半年各月徑流減少的情況應引起關注。③從年內分布的整體來看，ΔR1～ΔR4 的最大月份均在8 月，最枯月份均在3 月；C1～C4 均使7—12 月ΔR遠大于1—6 月的，降水是本流域水資源的主要補給方式，且多集中在夏季和秋季，河流進入汛期，C1～C4 對汛期徑流的影響更大。

圖8 C1～C4 引起的徑流深變化（ΔR）的年內分布

依據(jù)Q1～Q5 中每月Q值按式（7）和式（8）計算每年的Cr值，而后分別和情景1 的Cr值做差得ΔCr，繪制2006—2018 年ΔCr分布圖。由圖9 可知：①2006—2018 年的ΔCr1～ΔCr4 均大于0，變化幅度在0.14 之間，說明一定程度上C1～C4 均使徑流深的年內分配不均勻，原因是C1～C4 各種土地類型面積的12%左右發(fā)生流轉，造成下墊面在小范圍內變化。②2006—2018 年的ΔCr1、ΔCr3 值均小于ΔCr2、ΔCr4，說明C1、C3 使徑流深年內分配相對來說更均勻些，雖然一年內大部分降水集中在夏秋兩季，使河水上漲，但是C1、C3 使耕地面積不同程度增加，無定河流域地處陜北，多種玉米和馬鈴薯［7］，一年內均需消耗不同程度的灌溉用水量，因此使徑流深年內分配不均勻程度較小。③ΔCr1～ΔCr4 中，2014—2018 年的差值整體均呈現(xiàn)上升趨勢，ΔCr2 的上升趨勢最為顯著，說明C2 對2014—2018 年的徑流深年內分配不均勻程度影響最大。

圖9 C1～C4 引起的年內分配不均勻系數(shù)的分布

4 結 論

為研究無定河流域2005 年后的土地利用的小幅度變化對徑流的影響，采用2005 年、2015 年、2020 年三期土地利用數(shù)據(jù)和基于CA－Markov 模型預測的2025 年、2030 年兩期數(shù)據(jù)，分析4 種土地利用變化方式，并構建SWAT 模型，研究不同土地利用變化方式對徑流的影響。

（1）2005 年后無定河流域草地、耕地和未利用地占總流域面積的90%以上。2005—2030 年耕地、草地在小范圍內互相轉化；未利用地持續(xù)減少，主要轉化為耕地和草地；居民用地占流域面積較少，但擴張極其迅速。

（2）基于2005 年土地利用數(shù)據(jù)構建的SWAT 模型，經過參數(shù)率定后，在率定期和驗證期的r和Ens均在0.8 以上，說明對無定河流域水文過程的模擬效果良好，適用于研究本流域土地利用變化對徑流的影響。

（3）4 種土地利用變化均使本流域2006—2018 年徑流量有不同程度增加，且居民用地面積增加在促進徑流增加中占主導作用。

（4）4 種土地利用變化方式對7—12 月徑流的促進作用均遠大于1—6 月；草地增加地表粗糙度，居民用地增加不透水面積，二者是使ΔR2 最大的主要原因；第4 種土地利用變化對上半年徑流起抑制作用，對下半年徑流起促進作用。

（5）4 種土地利用變化均使2006—2018 年的徑流年內分配不均勻，第1 和第3 種土地利用變化使徑流年內分配相較于其他兩種變化來說均勻些。

（6）第4 種土地利用變化雖然未利用地總面積減少，卻在東北部大塊出現(xiàn)，并且引起徑流年際增加量大、上半年徑流減少下半年徑流增加、歷年徑流深的年內分配極不均勻，該情況應引起關注。

本文利用土地轉移概率矩陣進行未來土地利用預測，沒有考慮政策導向和氣候變化的影響，同時所采用土地利用數(shù)據(jù)的分辨率也有待于進一步提高，進而提高成果的精度，這將是下一步研究的重點。