李 敏， 朱清科

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

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

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

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

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

1 研究區(qū)域與時段

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

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

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

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

黃河年輸沙量承載了黃土高原生態(tài)環(huán)境狀況、水土流失與水土保持的信息。

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

2 研究方法

2.1 原理與方法

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

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

圖3 1980—2015年黃河中游黃土高原地區(qū)水土流失治理面積變化(資料來源：黃河上中游管理局)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 降雨因子選擇

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

2.3 水土保持因子選擇與預(yù)處理

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

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

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

2.4 研究時段劃分

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

3 結(jié)果與分析

3.1 1954—1966年計算結(jié)果分析

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

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

統(tǒng)計分析表明，由于1954—1966年研究區(qū)域黃河年輸沙量較大，年際間的波動劇烈；同時該時段水土保持工程面積較小，治理薄弱，水土流失治理度偏低，黃河4站年輸沙量對水土流失治理度沒有統(tǒng)計學(xué)意義上的響應(yīng)(圖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

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

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

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

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

3.2 1967—1987年計算結(jié)果分析

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

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

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

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

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

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

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

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

A為1967年水土流失治理度水平下年輸沙量響應(yīng)曲面,B為1987年水土流失治理度水平下年輸沙量響應(yīng)曲面。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 年輸沙量對不同水土流失治理度的響應(yīng)曲面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年計算結(jié)果分析

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

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

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

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

3.3.2 1988—2015年時段減沙分析 進(jìn)一步分析1988—2015年回歸結(jié)果，計算在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。

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

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

4 結(jié)論與討論

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

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

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