李 敏， 朱清科

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

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

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

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

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

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

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

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

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

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

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

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

2 研究方法

2.1 原理與方法

根據(jù)水文學(xué)和水土保持學(xué)基本原理，河流年輸沙量對(duì)降雨和水土保持的響應(yīng)遵從統(tǒng)計(jì)學(xué)規(guī)律。因此，筆者將研究區(qū)域黃河年輸沙量的產(chǎn)生過程視為一個(gè)灰箱或者黑箱，采用逐步回歸分析方法，篩選具有物理成因概念和含義的解釋因子(自變量)，獲得符合水文學(xué)和水土保持學(xué)基本原理、仿真模擬效果好的計(jì)算結(jié)果，以此研究黃河年輸沙量對(duì)水土保持響應(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ū)間(以下簡(jiǎn)稱“河龍區(qū)間”)、源自汛期降雨(暴雨)；因此，對(duì)于降雨因子，初步選取4站區(qū)域年降雨量，即河龍區(qū)間(包括河龍區(qū)間5個(gè)亞區(qū))年降雨量、7—8月降雨量、最大3日降雨量、最大一日降雨量等。以上共初步選擇22個(gè)降雨因子。

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

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

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

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

2.4 研究時(shí)段劃分

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

3 結(jié)果與分析

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

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

圖5 1954—1966年實(shí)測(cè)年輸沙量與計(jì)算年輸沙量比較Fig.5 Comparison of measured annual sediment transport with calculated annual sediment transport in 1954 to 1966

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

圖4 1954—1966年黃河4站年輸沙量與水土流失治理度散點(diǎn)圖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)過對(duì)22個(gè)降雨因子和水土流失治理度因子的逐步回歸計(jì)算結(jié)果表明，僅“河吳區(qū)間7—8月降雨量”1個(gè)因子可以“解釋”該時(shí)段年輸沙量近78%(R=0.892 330,調(diào)整的R2=0.777 731)的變化。線性回歸顯著性水平α=0.001。

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

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

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

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

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

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

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

圖6 1967—1987年實(shí)測(cè)年輸沙量與計(jì)算年輸沙量比較Fig.6 Comparison of measured annual sediment transport with calculated annual sediment load in 1967 to 1987

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

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

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

這一現(xiàn)象說明水土流失治理不僅大幅度減少了黃河年輸沙量，而且使年輸沙量對(duì)降雨變化響應(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 年輸沙量對(duì)不同水土流失治理度的響應(yīng)曲面Fig.7 Response surface of annual sediment transport with different controlled percentage of soil erosion

圖8 1988—2015年實(shí)測(cè)年輸沙量與計(jì)算年輸沙量比較Fig.8 Comparison of measured annual sediment transport and calculated annual sediment transport in 1988 to 2015

3.3 1988—2015年計(jì)算結(jié)果分析

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

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

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

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

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

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

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

4 結(jié)論與討論

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

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

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