基于SWAT模型的清水河流域年徑流侵蝕功率空間分布

楊 光, 魯克新, 李 鵬,, 柳海亮, 劉桂華, 曹峰華, 馬天文, 王得軍, 文妙霞, 孫景梅

(1.西安理工大學 省部共建西北旱區生態水利國家重點實驗室, 西安 710048; 2.西安理工大學 旱區生態水文與災害防治國家林業局重點實驗室, 西安 710048)

SWAT模型是美國農業部農業研究中心于20世紀90年代研發的以SWRRB為框架，同時融入了CREAMS、GLEAMS、EPIC和ROTO等模型優點的流域尺度長時段分布式水文模型[8]。國內外學者多使用SWAT模型進行流域水沙研究。Zhang等[9]有關土地利用變化對于渾河流域產水產沙的影響的模擬結果表明，林地全年產沙量減少，水分入滲增加；草地和林地產流產沙對土地利用變化的響應相似，但前者的水土保持效果弱于后者。邢立文等[10]研究結果表明，氣候變化對產沙的影響比產流大，產流、產沙與降雨呈正相關關系而與氣溫呈負相關關系，降雨量的變化對產流和產沙的影響比氣溫大。劉世梁等[11]有關瀾滄江中游小流域水土流失變化與NDVI時空變化的研究結果表明，瀾滄江中游小流域的徑流產沙量在空間尺度上隨著流域植被NDVI的增大而減小，在時間尺度上與流域植被NDVI間的關系不顯著。朱楠等[12]研究了羅玉溝土地利用結構對流域水沙的影響，結果表明林地的減水、減沙效果最好，坡耕地的產流、產沙量最大；當林地分布在流域上游時，到達流域中游出口的水沙量最小。然而，以往研究學者主要利用SWAT模型研究氣候變化和土地利用變化對流域水沙的影響等，而對于黃河流域上游徑流侵蝕功率的空間分布規律的研究還不夠深入。

本研究以寧夏回族自治區水土流失嚴重的清水河流域為研究對象，基于SWAT模型開展清水河流域徑流模擬，并使用程圣東[4]提出的能更好表示流域年徑流過程的年徑流侵蝕功率公式計算各子流域出口斷面的年徑流侵蝕功率，進而進行年徑流侵蝕功率的空間分布特征和尺度效應研究，以期從徑流侵蝕功率角度解釋流域水力侵蝕的空間分布情況，為后續在清水河流域不同區域開展針對性的水土保持措施布設以及生態治理工作提供理論支撐。

1 研究區概況

清水河是黃河寧夏段最大的支流，地理位置為東經105°09′—116°43′，北緯35°49′—37°30′，發源于六盤山東麓固原市，向北流經原州區、海原縣、同心縣、中寧縣等縣(區)，在中寧縣的泉眼山西側注入黃河，全長320 km，流域面積14 481 km2，河道平均縱比降1.49‰。流域地貌以黃土丘陵溝壑區為主，地勢南高北低，河源海拔2 489 m，入黃口海拔1 190 m，相對高差1 299 m。

清水河流域多年平均降水量349 mm，年降水量空間分布不均，由南到北遞減，流域上、下游的年降水量分別為600，200 mm。流域內年平均氣溫由南向北遞減，北部中衛市年平均氣溫8℃以上，中部同心縣年平均氣溫7～8℃，南部固原市年平均氣溫5～6℃。清水河流域日照多、濕度小、風大，水面蒸發強烈，多年平均水面蒸發量1 272 mm，年水面蒸發量介于900～1 300 mm[13]。清水河流域多年平均入黃泥沙量約占寧夏黃河流域入黃泥沙總量的49%[14]。

2 材料與方法

2.1 數據來源

SWAT模型的輸入數據包括流域數字高程模型(DEM)數據、土地利用數據、土壤類型數據及逐日降水量等氣象數據。根據土壤類型數據使用SWAT模型的SPAW工具計算得到土壤下滲率、電導度、可持水量等土壤參數數據，進而制備SWAT模型運行所需的土壤參數數據庫。根據逐日降水量氣象數據，使用SWAT Weather工具計算各氣象站點的各月月降水量的均值標準差、最高和最低氣溫均值和標準差、月均濕度、月均太陽輻射以及月均風速數值，制備天氣發生器，為后續流域徑流模擬提供參考。

本研究使用的清水河流域空間數據包括中科院地理空間數據云30 m分辨率的DEM數據、中國科學院遙感所解譯的1980年流域土地利用數據以及世界土壤數據庫(Harmonized World Soil Database，HWSD)中的1∶100萬比例尺流域土壤類型數據。

本研究使用的水文氣象數據包括：1976—1986年《黃河流域水文年鑒》摘錄歷年的清水河流域22個雨量站的逐日降水量以及清水河干流韓府灣水文站的逐日平均流量；來自國家氣象科學數據中心的中衛、中寧、同心、海原、固原和西吉6個氣象站點的逐日降水量、氣溫、風速、太陽輻射以及相對濕度等氣象數據。清水河流域水文站、氣象站和雨量站空間分布圖見圖1。

圖1 清水河流域DEM數據、站點分布

2.2 模型建立

本研究使用ARCSWAT 2012版本對清水河流域進行徑流模擬。在進行流域劃分時，本研究把集水區面積閾值設置為10 000 hm2并以清水河入黃口作為流域出口，將整個清水河流域劃分為93個子流域，詳見圖1B。在建立水文響應單元(HRU)時，需要輸入土地利用數據、土壤類型數據和坡度數據。

土地利用類型對流域產匯流具有重要影響，將預先準備好的1980年清水河流域的土地利用數據重分類為耕地、林地、草地、水域、城鎮居民用地以及未利用地6類土地利用類型。統計結果表明，清水河流域最主要的土地利用類型為草地，面積占流域面積的56%；其次為耕地，面積占流域面積的37%，而其他4種土地利用類型面積合計占流域面積的比例不足10%。清水河流域的耕地主要沿河道兩岸分布，而流域北部分布有大范圍的草地以及未利用地。清水河流域1980年土地利用類型空間分布圖見圖2。

圖2 清水河流域1980年土地利用類型空間分布

輸入土壤數據，從HSWD的土壤數據庫中提取清水河流域的土壤數據，重新分類為23種不同的土壤類型；將流域內的坡度劃分為5級(0～15%，15%～25%，25%～35%，35%～45%，45%～99.99%)。

在HRU定義時，將土地利用數據、土壤類型數據以及流域坡度的閾值都設置為5%，共劃分出2 058個HRU。通過閾值的設定，可以減少面積過小的土地利用、土壤類型以及坡度對于模型模擬結果精度的影響[15]。將整理好的氣象數據如逐日的降水、氣溫、風速、太陽輻射以及相對濕度數據輸入SWAT模型，隨后將其寫入數據庫后即可運行SWAT模型。

為了確保SWAT模型模擬結果的準確性，本研究以清水河干流韓府灣水文站1976年、1977年作為預熱期，以1978—1982年歷年各月實測月平均流量數據進行模型參數率定，以1983—1986年歷年各月實測月平均流量數據進行模型驗證。

2.3 模型參數敏感性分析

影響SWAT模型模擬結果的參數眾多，為了減少模型運行的不確定性和提高模型的運行效率，需要對模型的相關參數進行敏感性分析。

本研究使用LH-OAT算法對SWAT模型的參數進行敏感性分析，以p和t檢驗值作為參數敏感程度的評判標準，p值越接近于0且t值越大，則說明參數越敏感[16]。本研究選擇對模型模擬結果敏感性較高的參數參與后續的模型參數率定和驗證。

2.4 模型率定與驗證

SWAT模型運行成功后，需要驗證模型模擬結果的可靠性。本研究選取清水河干流的韓府灣水文站的實測徑流資料對模型模擬結果進行驗證。選取參數敏感性分析中的敏感性較高的參數參與模型率定和驗證。使用SUFI-2算法對率定期參數進行率定，得到最優參數值，進而將最優參數值帶入驗證期進行參數驗證。

本研究使用決定系數R2以及納什效率系數NS對清水河SWAT模型的適用性進行評估。R2表示兩個變量間的相關程度，值越接近于1說明兩個變量間的相關關系越好；NS一般用于評價觀測值與模擬值的擬合程度，值越接近于1說明擬合情況越好，當NS=1時，觀測值和模擬值完全相等。一般研究[17-21]認為，當R2和NS均大于0.5時，模型模擬結果可信。如果模型率定期和驗證期的模擬結果均可信，則認為最優參數值適用于研究流域，模型模擬結果能較好地表示研究流域的實際徑流情況。

2.5 徑流侵蝕功率

程圣東[4]將徑流侵蝕功率由場次尺度推廣到年尺度的，考慮年徑流深為Hy，年內最大月平均流量為Qm，得到年尺度徑流侵蝕功率計算公式如下：

(1)

(2)

(3)

3 結果與分析

3.1 模型率定與驗證

使用清水河干流韓府灣水文站1978—1986年的實測月平均流量數據對清水河流域SWAT模型進行參數率定和驗證。以1979—1982年的月平均流量為率定期數據，采用LH-OAT方法進行模型參數敏感性分析，選擇18個敏感性強的參數進行率定，結果見表1。

從表1中可以看出，對清水河流域徑流模擬結果影響最大的參數為決定流域徑流量大小的SCS徑流曲線數CN2，CN2越大則地表徑流量就越大；其次是主要影響下滲情況的表層土壤容重SOL_BD以及主要表示河岸地下徑流對于河道流量的補給快慢的主河道調蓄系數ALPHA_BNK。另外，其余參數如降雪氣溫SFTMP、淺層地下水徑流系數GWQMN、基流消退系數ALPHA_BF、主河道曼寧系數CH_N2都對流域徑流模擬結果的影響較為敏感。

表1 清水河流域SWAT模型參數敏感性排名

將收集到的韓府灣水文站實測月平均流量數據劃分為預熱期(1976年和1977年)、率定期(1978—1982年)和驗證期(1983—1986年)3個階段。在SWAT-CUP軟件中，使用率定期數據對上文選取的18個參數對進行參數率定，將參數進行優化迭代，直至決定系數R2>0.6且納什效率系數NS>0.5，則可認為模型的模擬結果可信。將率定期得到的最優參數值回代到驗證期后再次運行SWAT模型，并根據驗證期的月流量模擬值與實測值再次計算決定系數R2和納什效率系數NS。若驗證期的決定系數R2>0.6且納什效率系數NS>0.5，則認為SWAT模型的模擬結果可信，確定的模型最優參數值可用于后續的月流量模擬計算。

韓府灣水文站率定期(1978—1982年)和驗證期(1983—1986年)月平均流量實測值與模擬值對比結果見圖3。

圖3 1978-1986年韓府灣站月平均流量過程模擬值與實測值對比

經計算，模型率定期的決定系數R2為0.76、納什效率系數NS為0.74；驗證期的決定系數R2為0.88，納什效率系數NS為0.85。因此認為清水河流域月流量過程的模擬結果可信，SWAT模型在清水河流域具有良好的適用性。

利用SWAT模型可以得到1978—1986年清水河流域入黃口斷面以及流域93個子流域出口斷面的月流量過程數據以及歷年年徑流深等數據。

3.2 徑流侵蝕功率空間分布

清水河流域SWAT模型輸出文件中的Sub表中的子流域產水量(WYLD)字段包含1978—1986年93個子流域的逐月月產水量數據，以此可以推求出各子流域的多年平均徑流深H；經統計分析，從Sub表中查得各子流域歷年的最大月徑流深hm。將徑流侵蝕功率公式進一步變換得到公式(4)—(5)。

hmA′=QmΔt

(4)

(5)

式中：hm為年內最大月徑流深(m)；A′為子流域面積(km2)；Δt為時段長(s),此處平均按每月30.4 d計算各月總秒數,Δt=2.627×106s。

由公式(5)得到1978—1986年清水河各子流域歷年年內最大月平均流量模數Q′m，最后由公式(1)計算得到相應年份各子流域的年徑流侵蝕功率。利用1978—1986年清水河各子流域歷年的年徑流侵蝕功率結果求得多年平均徑流侵蝕功率，并基于ArcGIS平臺制作清水河流域多年平均徑流侵蝕功率空間分布圖(圖4)。

從圖4可以看出，清水河流域多年平均徑流侵蝕功率呈現“支流大，干流?。粬|部大，西部小”的空間分布規律，即各支流區域的多年平均徑流侵蝕功率大于清水河干流，且流域西部地區的多年平均徑流侵蝕功率明顯小于東部地區。

圖4 1976-1986年清水河流域多年平均徑流侵蝕功率空間分布 圖5 1976-1986年清水河流域各子流域出口斷面多年平均徑流侵蝕功率空間分布

3.3 年徑流侵蝕功率的空間尺度效應

為了探究流域徑流匯聚效應與流域徑流侵蝕功率的關系，本節進一步研究了子流域年徑流侵蝕功率與子流域出口斷面集水面積間的相關關系。

SWAT模型輸出的Rch表包含各子流域出口斷面的集水面積和月流量(FLOW-OUT)。子流域出口斷面的集水面積指子流域出口斷面以上的集水區總面積。子流域出口斷面流量(FLOW-OUT)指子流域出口斷面以上匯聚到子流域出口處的流量。以各子流域出口斷面的集水面積A′和多年平均月平均流量(FLOW-OUT)，依據公式(1)—(3)計算各子流域出口斷面多年平均徑流侵蝕功率。

由圖5可知，在各級支流的徑流匯聚過程中，子流域出口斷面的多年平均徑流侵蝕功率都有減少的趨勢，尤其以流域上游地區較為明顯，但此趨勢在清水河下游集水面積較大的子流域間不太明顯。因此，本研究將出口斷面集水面積小于4 000 km2的子流域的集水面積與子流域出口斷面多年平均徑流侵蝕功率進行擬合分析，試圖找到子流域出口斷面集水面積與子流域出口斷面多年平均徑流侵蝕功率間的相關關系，二者關系散點圖見圖6。由圖6可知，子流域出口斷面多年平均徑流侵蝕功率隨著子流域出口斷面集水面積的增大而迅速減小，直到出口斷面集水面積達到某一數值時子流域出口斷面多年平均徑流侵蝕功率對其敏感性降低以至于趨于穩定。

圖6 子流域出口斷面集水面積與多年平均徑流侵蝕功率擬合分析

經回歸分析，當子流域出口斷面集水面積小于4 000 km2時，清水河子流域出口斷面多年平均徑流侵蝕功率E與子流域出口斷面集水面積Ac間呈如下顯著的冪函數關系

(6)

經統計分析，當子流域出口斷面集水面積大于4 000 km2時，出口斷面多年平均徑流侵蝕功率數值穩定在1.56×10-5m4/(s·km2)左右。分析其主要原因可能是清水河流域右岸支流的匯水量過大，導致在下游地區干流與其他支流間的多年平均徑流侵蝕功率空間分布規律不明顯。

為確定子流域出口斷面集水面積的臨界值，將公式(6)進行求導得到如下導數方程：

(7)

|E′|<0.001說明子流域出口斷面多年平均徑流侵蝕功率變化速率減小，其值趨于穩定[6]。當出口斷面集水面積Ac大于84.85 km2時，|E′|<0.001，子流域出口斷面多年平均徑流侵蝕功率不再隨著子流域集水面積的增大而明顯減少，因此影響清水河流域多年平均徑流侵蝕功率大小的子流域出口斷面集水面積空間閾值為84.85 km2。

依據魯克新等[3]提出的流域次暴雨輸沙模型研究結果，流域次暴雨徑流侵蝕功率與流域輸沙模數呈冪函數關系，徑流侵蝕功率越大則流域輸沙模數越大。因此，在進行清水河流域生態治理時，優先選擇處于流域上中游區且出口斷面集水面積小于84.85 km2的小流域，可取得良好的水土流失治理效果。

4 討 論

4.1 年徑流侵蝕功率空間分布規律

通過本研究以及龔珺夫等[5-6]在無定河及延河的徑流侵蝕功率空間分布的研究發現，不同流域的年徑流侵蝕功率存在區域性差異的同時，也有共同的特點。從空間角度來看，流域年徑流侵蝕功率的空間分布情況與流域下墊面條件、氣候情況相關，不同流域會具有不同的年徑流侵蝕功率空間分布特點。從流域水系的角度來看，不同流域的多年平均徑流侵蝕功率都存在“支流大，干流小”的空間分布規律。

4.2 年徑流侵蝕功率的空間尺度效應

為了更好研究子流域出口斷面徑流侵蝕功率與子流域出口斷面集水面積間的相關關系，本研究相較于龔珺夫[5-6]的研究選取了程圣東[4]提出的年徑流侵蝕功率，使用年徑流深作為徑流侵蝕能量的參數，以期更好地表達流域的徑流侵蝕功率情況。研究結果表明，清水河子流域多年平均徑流侵蝕功率與子流域出口斷面集水面積間呈顯著的冪函數關系，與龔珺夫等[5-6]的研究結果相似。不同流域的徑流侵蝕功率的子流域集水面積閾值并不相同。對于清水河流域而言，小于流域集水面積閾值的子流域均位為各支流的河源處，即圖5中徑流侵蝕功率較大的區域，這對于清水河流域生態治理具有一定的實際指導意義。

4.3 不足與展望

本研究基于SWAT模型模擬研究了清水河流域的年徑流侵蝕功率的空間分布規律及其空間尺度效應，但并未闡明清水河流域的多年平均徑流侵蝕功率存在“支流大，干流小”空間分布規律的機理。另外，本研究沒有進行清水河流域產沙輸沙量模擬，因此無法將模擬得到的清水河流域徑流侵蝕功率用于流域輸沙模數計算。以上兩個方面的不足還有待今后深入開展相關研究，以期為流域侵蝕產沙預報與流域生態治理提供理論指導。

5 結 論

(1) 本文構建的清水河流域SWAT模型的率定期(1978—1982年)決定系數R2為0.76，納什效率系數NS為0.74；驗證期(1983—1986年)決定系數R2為0.88，納什效率系數NS為0.85，因此SWAT模型在清水河流域徑流模擬中具有良好的適用性，模擬的徑流結果可以較好地表示流域的實際徑流情況。

(2) 清水河流域的多年平均年徑流侵蝕功率具有“支流大，干流?。粬|部大，西部小”的空間分布規律。

(3) 清水河子流域多年平均徑流侵蝕功率隨著出口斷面集水面積的增大而迅速減小，而當子流域出口斷面集水面積達到臨界閾值時出口斷面徑流侵蝕功率對其敏感性降低并趨于穩定。當子流域的出口斷面集水面積小于4 000 km2時，子流域多年平均徑流侵蝕功率與出口斷面集水面積間呈顯著的冪函數關系。清水河上中游流域多年平均徑流侵蝕功率的流域出口斷面集水面積臨界閾值約為84.85 km2。優先選擇處于清水河上中游區且出口斷面集水面積小于84.85 km2的小流域進行水土保持措施布設及生態治理，可取得良好的水土流失治理效果。