20世紀中期以來不同時段黃河年輸沙量對水土保持的響應

李 敏， 朱清科

(1.水利部 黃河水利委員會 黃河上中游管理局， 710021，西安； 2.北京林業大學林業生態工程教育部工程研究中心，100083，北京)

黃土高原是世界上水土流失最嚴重的地區，流經黃土高原的黃河是世界上年輸沙量最大的河流之一。20世紀中期黃土高原地區的水土保持就納入國家重點治理。經過近70年的大規模治理，水土保持工程對遏制黃土高原水土流失產生了明顯的作用，進而對減少黃河年輸沙量發揮了巨大的“治本”效果。

近期，國家主要領導人視察黃河，提出：治理黃河，重在保護，要在治理。要堅持山水林田湖草綜合治理、系統治理、源頭治理。深入研究黃河年輸沙量對水土保持的響應特性與效應，對于評判黃河治理效益，確定黃河治理重點具有重要意義。

黃河輸沙量的變化與降雨和水土保持緊密相關。以往在研究水土保持與黃河年輸沙量的關系中，多采用水文法和水保法、按照年代劃分時段進行計算[1-5]。

理論和實踐表明，在自然和人類活動的“干擾”下，河流輸沙量對水土保持的響應遵從統計學原理。筆者在以往分析50多年整體系列的基礎上[11]，進一步研究了不同時段黃河年輸沙量對水土保持的響應特性。

1 研究區域與時段

黃河泥沙基本上來自流經黃土高原地區的河口鎮(頭道拐)到龍門區間的各支流和渭河、北洛河、汾河，其中主要產沙區域為黃河的河口鎮(頭道拐)到龍門區間(以下簡稱“河龍區間”)。本研究區域為渭河(華縣站)、北洛河(狀頭站)、汾河(河津站)控制區域和河龍區間(以下簡稱“4站區域”)。4站區域總面積30.05萬km2(圖1)。

圖1 研究區域示意圖Fig.1 Schematic diagram of the study area

黃河自1919年設站觀測，距今不到100年。本研究時段為1954—2015年，共62年。

1954—2015年(不計頭道拐以上來沙)黃河年輸沙量波動變化，進入21世紀平均值僅2.6億t(圖2)[6-10]，其中2015年僅0.55億t。

黃河年輸沙量承載了黃土高原生態環境狀況、水土流失與水土保持的信息。

從20世紀50年代初期開始，黃河中游黃土高原地區的水土保持就開始了大規模的實施。20世紀70年代以前，由于受技術和資金條件的限制，工程質量普遍較低，水毀嚴重，期間還經歷過毀林開荒等人為破壞，水土保持設施的面積和數量增長緩慢。進入20世紀80年代，黃河中游水土保持工程質量得到提升，資金投入得到保障，法規監督得以落實，水土流失治理面積基本以線性趨勢持續增加(圖3)。

2 研究方法

2.1 原理與方法

根據水文學和水土保持學基本原理，河流年輸沙量對降雨和水土保持的響應遵從統計學規律。因此，筆者將研究區域黃河年輸沙量的產生過程視為一個灰箱或者黑箱，采用逐步回歸分析方法，篩選具有物理成因概念和含義的解釋因子(自變量)，獲得符合水文學和水土保持學基本原理、仿真模擬效果好的計算結果，以此研究黃河年輸沙量對水土保持響應的特性。

圖2 1954—2015年黃河4站區域年輸沙量變化圖Fig.2 Change of annual sediment transport in the area of 4 stations in the Yellow River in 1954 to 2015

圖3 1980—2015年黃河中游黃土高原地區水土流失治理面積變化(資料來源：黃河上中游管理局)Fig.3 Change of the controlled area of soil erosion in the Loess Plateau of the middle reaches of the Yellow River in 1980 to 2015(Source: Upper and Middle Yellow River Bureau)

2.2 降雨因子選擇

已知黃河泥沙主要源自河口鎮到龍門區間(以下簡稱“河龍區間”)、源自汛期降雨(暴雨)；因此，對于降雨因子，初步選取4站區域年降雨量，即河龍區間(包括河龍區間5個亞區)年降雨量、7—8月降雨量、最大3日降雨量、最大一日降雨量等。以上共初步選擇22個降雨因子。

2.3 水土保持因子選擇與預處理

計算中對水土保持坡面措施面積、大型壩(骨干壩)壩控面積進行“耦合”處理[11-12]，以使坡面措施與淤地壩相互協調，既不重復計算，又不產生矛盾。

為了使今后對不同時間和空間尺度的計算結果具有可比性，筆者將“水土流失治理度”作為參選自變量。水土流失治理度的定義為

式中：D為水土流失治理度，%；A1為水土保持綜合治理面積，萬km2；A2為水土流失面積，萬km2 [13-15]。

2.4 研究時段劃分

對于黃河年輸沙量，一般認為：在20世紀60年代前后為“天然狀態”的輸沙量，即流域下墊面破壞嚴重，水土保持等治理工程措施薄弱階段；在70—80年代為水土保持措施減沙效果顯現時期，該時期減沙效果波動變化；80年代后期以來水土保持工程項目持續實施，減沙效果越亦顯著。據此，本文將1954—2015年劃分為年輸沙量“高位運行，劇烈波動” “波動遞減” “低位運行”3個時段分別進行計算分析。

3 結果與分析

3.1 1954—1966年計算結果分析

觀察1954—1966年黃河4站年輸沙量變化(圖2中1954—1966年時段)可以看出：該時段黃河年輸沙量均值達到15.7億t；波動劇烈，年際間波動幅度最大超過20億t。簡單回歸分析表明，盡管該時段年輸沙量波動幅度很大，但是時段變化趨勢微弱，基本沒有統計學意義上的增減變化趨勢。

圖5 1954—1966年實測年輸沙量與計算年輸沙量比較Fig.5 Comparison of measured annual sediment transport with calculated annual sediment transport in 1954 to 1966

統計分析表明，由于1954—1966年研究區域黃河年輸沙量較大，年際間的波動劇烈；同時該時段水土保持工程面積較小，治理薄弱，水土流失治理度偏低，黃河4站年輸沙量對水土流失治理度沒有統計學意義上的響應(圖4)。

圖4 1954—1966年黃河4站年輸沙量與水土流失治理度散點圖Fig.4 Scatter plot of annual sediment load with controlled percentage of soil erosion from 1954 to 1966 in the area of four stations in the Yellow River

經過對22個降雨因子和水土流失治理度因子的逐步回歸計算結果表明，僅“河吳區間7—8月降雨量”1個因子可以“解釋”該時段年輸沙量近78%(R=0.892 330,調整的R2=0.777 731)的變化。線性回歸顯著性水平α=0.001。

比較回歸模型計算值與實測值(圖5)看出，計算值與實測值擬合較好，反映了黃河年輸沙量年際變化來自于對河吳區間7—8月降雨量年際變化的響應。

其中，1966年回歸標準化殘差為2.184 6，該年的計算值明顯小于實測值。初步分析，造成這一現象的原因是該年黃河流域發生較大的局部暴雨，沖毀了部分水土保持工程措施，造成以往多年攔淤減少的泥沙在當年被沖入黃河，使“河吳區間7—8月降雨量”不能夠“解釋”這些輸沙量，因而計算年輸沙量明顯小于實測年輸沙量[16-18]。

對1954—1966年時段的回歸計算結果說明，這一時段黃河年輸沙量的變化主要來自汛期降雨的變化，由于該時段水土流失治理措施規模較小，輸沙量對其沒有響應。

3.2 1967—1987年計算結果分析

3.2.1 計算結果的總體分析 觀察1967—1987年黃河4站年輸沙量變化(圖2 1967—1987年時段)可以看出：該時段黃河年輸沙量波動下降，簡單回歸分析表明，該時段年輸沙量呈指數減少趨勢。

經過對1967—1987年22個降雨因子和水土流失治理度因子的逐步回歸篩選計算，建立了由“河吳區間7—8月降雨量”“水土流失治理度”“河龍區間最大一日降雨量”“ln區域年降雨量”等4個自變量組成的非線性回歸方程。回歸方程R=0.954 689，顯著性水平α=0.001。

根據回歸方程，對計算年輸沙量與實測年輸沙量比較(圖6)看出，計算值與實測值擬合很好，不但描述了輸沙量年際間的波動變化，而且反映了長期遞減的趨勢。

圖6 1967—1987年實測年輸沙量與計算年輸沙量比較Fig.6 Comparison of measured annual sediment transport with calculated annual sediment load in 1967 to 1987

回歸分析表明，1967—1987時段，黃河年輸沙量的主要影響因子仍然是汛期降雨，年際間輸沙量的波動變化響應來自汛期降雨，“河吳區間7—8月降雨量”因子“解釋”了該時段黃河年輸沙量接近63%的變化，其次“水土流失治理度”1個因子可以“解釋”黃河年輸沙量接近50%的變化，黃河年輸沙量波動下降的現象主要是對該因子的響應結果。

3.2.2 回歸方程的可視性分析 由于是非線性關系，不同水土流失治理度對年輸沙量的影響程度不同，影響的過程(趨勢)也不同。根據回歸方程做出在相同降雨條件下不同水土流失治理度不同視角的三維曲面圖形，直觀對比分析不同水土流失治理度對年輸沙量的影響程度和影響趨勢。為了做出三維圖形，對貢獻率最小的“ln區域年降雨量”固定為時段的平均值，選擇時段起始年(1967年)和時段末尾年(1987年)，分別做出該2個年份年輸沙量對水土流失治理度的響應曲面進行比較分析(圖7)。

從圖7總體上看出，1967年和1987年的水土流失治理度對年輸沙量影響顯著不同。在曲面的最高點，即2個降雨因子均處于最大值時，1987年輸沙量響應曲面比1967年輸沙量響應曲面降低超過20億t；在曲面的最低點，即2個降雨因子均處于最小值時，2個年輸沙量響應曲面相差幅度約5億t。對比分析說明，隨著水土流失治理度的增加，一方面，總體上年輸沙量顯著減少，另一方面，年輸沙量響應曲面在高降雨條件下的減少幅度大于低降雨條件[11]。

這一現象說明水土流失治理不僅大幅度減少了黃河年輸沙量，而且使年輸沙量對降雨變化響應的“敏感度”降低。

A為1967年水土流失治理度水平下年輸沙量響應曲面,B為1987年水土流失治理度水平下年輸沙量響應曲面。A is the annual sediment transport response surface for the controlled percentage of soil erosion in 1967, and B is the annual sediment transport response surface for the controlled percentage of soil erosion in 1987.圖7 年輸沙量對不同水土流失治理度的響應曲面Fig.7 Response surface of annual sediment transport with different controlled percentage of soil erosion

圖8 1988—2015年實測年輸沙量與計算年輸沙量比較Fig.8 Comparison of measured annual sediment transport and calculated annual sediment transport in 1988 to 2015

3.3 1988—2015年計算結果分析

3.3.1 回歸計算結果分析 在1988—2015年時段，黃河4站年輸沙量仍然以指數變化趨勢急劇波動下降，顯著減少(圖2 1988—2015年時段)。進入21世紀，年輸沙量在6億t以下繼續波動減少，其中2015年僅0.55億t。

1988—2015年時段，黃河年輸沙量對水土保持的響應十分顯著。回歸分析表明，僅“水土流失治理度”1個因子就能夠“解釋”黃河年輸沙量接近80%的變化(圖8)，回歸方程的相關系數R達到0.881 389，顯著性水平α=0.001。

如前分析，黃河年輸沙量的年際間波動來自于降雨量的年際變化，而長期遞減趨勢則主要來自于對水土流失治理度的響應，所以“水土流失治理度”的變化特性較好地擬合了黃河4站年輸沙量的變化趨勢。

圖9 1988—2015年回歸結果預測分析Fig.9 Prediction and analysis of regression results in 1988 to 2015

3.3.2 1988—2015年時段減沙分析 進一步分析1988—2015年回歸結果，計算在80%和95%置信度時的回歸計算值上限和下限，分析黃河4站年輸沙量變化趨勢(圖9)。2008年時水土流失治理度超過50%，計算年輸沙量不到2億t，其80%置信度上限5.9億t；2015年水土流失治理度超過60%，計算年輸沙量僅0.7億t，其80%置信度上限2.7億t，95%置信度上限5.6億t。

這一結果表明，根據1988—2015年數據建立的回歸模型分析，在當前水土流失治理程度條件下，以95%的置信度預測，黃河年輸沙量將小于6億t。

與研究區域人為破壞嚴重、水土流失治理薄弱狀態下的年輸沙量比較，當前黃河中游黃土高原地區的水土保持治理措施年減沙量超過10億t。

4 結論與討論

1)依據水土保持和降雨與黃河年輸沙量間的統計學關系，劃分3個時段，采用回歸分析方法對觀測數據進行計算。結果表明：1954—1966時段流域下墊面破壞嚴重、水土流失治理度偏低，黃河四站年輸沙量對水土保持沒有響應；1967—1987時段水土流失治理度逐漸增加，黃河四站年輸沙量開始對其有所響應；1988—2015時段水土流失治理度持續增加，黃河四站年輸沙量對其響應強烈，僅水土流失治理度1個因子就較好地“解釋”了20世紀80年代以來黃河年輸沙量銳減的趨勢。

2)在計算時段內(以1967—1987年時段為例)，相同降雨條件下，水土流失治理度的增加一方面顯著減少黃河年輸沙量，另一方面，隨水土流失治理度的提高，年輸沙量對降雨的響應趨于“遲緩”。由此，在1988—2015年時段，僅由水土流失治理度1個因子就能夠很好“解釋” 黃河年輸沙量的變化。

3)根據現有數據分析，在水土流失治理度超過60%(2015年治理度水平)時，以95%置信度預測黃河年輸沙量小于6億t。當前由于實施水土保持治理工程，黃河中游黃土高原地區的年度水土流失量減少了10億t以上。