金沙江下游區(qū)間來(lái)沙驅(qū)動(dòng)因子分析及產(chǎn)沙預(yù)測(cè)模型

譚寓寧，劉懷湘，陸永軍，3

(1. 南京水利科學(xué)研究院水文水資源及水利工程科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，江蘇 南京 210029；2. 四川大學(xué)水力學(xué)與山區(qū)河流開(kāi)發(fā)保護(hù)國(guó)家重點(diǎn)實(shí)驗(yàn)室，四川 成都 610065；3. 長(zhǎng)江保護(hù)與綠色發(fā)展研究院，江蘇 南京 210098)

隨著近年來(lái)長(zhǎng)江上游水電開(kāi)發(fā)的快速推進(jìn)，金沙江流域水庫(kù)群的蓄水?dāng)r沙導(dǎo)致劇烈的水沙調(diào)整，其中下游區(qū)間烏東德、白鶴灘、溪洛渡和向家壩4級(jí)巨型水庫(kù)的影響尤為顯著。向家壩、溪洛渡水庫(kù)分別于2012年10月和2013年5月蓄水運(yùn)用，蓄水初期向家壩下游年輸沙量在2 a內(nèi)即銳減至不足1998年前輸沙量均值的1%[1]。2011年以來(lái)金沙江流域水庫(kù)群年均總攔沙量達(dá)1.92億t，其中2013—2016年溪洛渡和向家壩兩庫(kù)年均攔沙量即占1.05億t[2]。

水庫(kù)攔沙淤積將影響長(zhǎng)期運(yùn)行效益，其中入庫(kù)沙量是決定淤積總量的控制性因素。金沙江流域來(lái)水來(lái)沙不均衡，水沙異源等現(xiàn)象十分突出[2-4]。攀枝花和屏山水文站分別為金沙江下游入口、出口控制站，兩者區(qū)間內(nèi)支流來(lái)沙量可占入庫(kù)總沙量的近80%。以區(qū)間為主的來(lái)沙特征將顯著影響庫(kù)區(qū)輸沙過(guò)程、淤積形態(tài)及優(yōu)化調(diào)控思路?，F(xiàn)有研究對(duì)該區(qū)間來(lái)沙往往基于少量支流水文站的觀測(cè)資料進(jìn)行簡(jiǎn)單估算[2，4-5]，較為精細(xì)化的研究非常缺乏。金沙江下游流域面積達(dá)500 km2以上的支流有32條，只有少數(shù)建有水文站，部分站點(diǎn)還存在時(shí)間序列短、位置遠(yuǎn)離河口等問(wèn)題。已有數(shù)據(jù)表明，各支流輸沙模數(shù)在3個(gè)量級(jí)內(nèi)變化，單位面積產(chǎn)沙能力差異極大，無(wú)法簡(jiǎn)單推算無(wú)資料支流的來(lái)沙。因此現(xiàn)有研究相對(duì)不足，難以認(rèn)識(shí)區(qū)間來(lái)沙全貌，從而也將影響對(duì)干流庫(kù)群泥沙淤積的預(yù)測(cè)能力。

區(qū)間來(lái)沙問(wèn)題實(shí)際上主要為流域產(chǎn)沙的問(wèn)題。山區(qū)河流泥沙有多種來(lái)源，包括降水導(dǎo)致的片蝕、細(xì)溝侵蝕和滑坡、泥石流、崩岸等重力侵蝕等[6-8]。影響山區(qū)流域產(chǎn)沙的因素如氣候、植被、地形等會(huì)由于當(dāng)?shù)氐耐寥狼治g和泥沙輸移機(jī)理而呈現(xiàn)出不同的驅(qū)動(dòng)因子組合[9-10]。再加上這些地區(qū)氣候、下墊面條件復(fù)雜且存在相互影響[11]，特別是支流小流域?qū)Ω饕蜃訑_動(dòng)的響應(yīng)比大流域更為直接、迅速[12]，對(duì)于該類(lèi)地區(qū)實(shí)際泥沙特性的研究更有難度。利用數(shù)學(xué)模型方法可有效評(píng)估流域產(chǎn)沙能力和預(yù)測(cè)擾動(dòng)影響下的泥沙響應(yīng)[13-14]。為此，獲取驅(qū)動(dòng)因子與產(chǎn)沙之間的定量關(guān)系并擇優(yōu)進(jìn)入模型搭建是該方法中的核心問(wèn)題。過(guò)去研究中通常使用直接線(xiàn)性回歸逐一確定因子間關(guān)系，以及通過(guò)偏最小二乘回歸(Partial Least Squares Regression，PLSR)法對(duì)全部因子回歸得到經(jīng)驗(yàn)?zāi)Ｐ蚚15-17]。但該類(lèi)回歸方法可能無(wú)法揭示因子間的實(shí)質(zhì)聯(lián)系，而且輸入過(guò)多的因子變量還會(huì)導(dǎo)致對(duì)因變量的重復(fù)解釋。此類(lèi)模型體量臃腫，實(shí)際應(yīng)用中需要大量的參數(shù)化和校準(zhǔn)，在資料匱乏的山區(qū)流域通常難以適用。因此，只有消除冗余因素、降低參數(shù)維度、識(shí)別和篩選出目標(biāo)主控因子，才能合理認(rèn)識(shí)目標(biāo)流域的產(chǎn)沙機(jī)制并搭建產(chǎn)沙預(yù)測(cè)模型。

明確金沙江下游區(qū)間泥沙全貌時(shí)空特征，辨析泥沙來(lái)源格局成因與主控因子，對(duì)合理利用管理與調(diào)控流域泥沙具有重要意義，可為金沙江下游巨型梯級(jí)水庫(kù)運(yùn)行方式的進(jìn)一步優(yōu)化提供依據(jù)，也可為高強(qiáng)度人類(lèi)活動(dòng)影響下長(zhǎng)江流域水沙通量變化研究奠定基礎(chǔ)。本研究基于金沙江18個(gè)支流水文站觀測(cè)數(shù)據(jù)，采用Spearman秩相關(guān)方法揭示因子集內(nèi)部以及與產(chǎn)沙之間的矩陣關(guān)系，通過(guò)偏相關(guān)方法對(duì)因子集進(jìn)一步降維，識(shí)別能夠反映區(qū)間來(lái)沙特性的多因子組合，采用PLSR-多元雙重回歸法構(gòu)建山區(qū)流域輸沙模數(shù)預(yù)測(cè)模型。

1 研究區(qū)域

金沙江多年平均懸移質(zhì)輸沙量為2.47億t，約占長(zhǎng)江上游輸沙量的47%，為長(zhǎng)江流域重要沙源。雅礱江匯口至宜賓為金沙江下游(圖1)，全長(zhǎng)約768 km，流域面積約為8.6萬(wàn)km2。降水主要集中在5—10月，可占年降水的80%以上，汛期平均降水量為855 mm。氣溫差異明顯，中部山區(qū)5—10月平均溫度可低至8 ℃，而干熱河谷地區(qū)溫度可達(dá)25 ℃。山高谷深的地勢(shì)使得區(qū)域內(nèi)熱量和水分條件都明顯垂向分化，形成了垂直氣候帶。區(qū)間內(nèi)如美姑河、黑水河等(圖2)眾多支流均具有典型山區(qū)河流特征，河道劇烈下切，河谷最大深度可達(dá)約3 100 m，形成了典型的“V”型侵蝕河谷地貌。流域內(nèi)崩塌、滑坡、泥石流等重力侵蝕頻繁，產(chǎn)沙豐富。

圖1 研究區(qū)域概況Fig.1 Distribution features of the study area

圖2 金沙江下游典型支流現(xiàn)場(chǎng)照片F(xiàn)ig.2 Typical tributary site photos

2 研究方法

2.1 數(shù)據(jù)來(lái)源

本研究中，支流水文站(16個(gè)位于金沙江下游，2個(gè)位于中游段出口區(qū)，見(jiàn)圖1和表1)泥沙觀測(cè)數(shù)據(jù)來(lái)源于《長(zhǎng)江流域水文年鑒》，選擇其中水電開(kāi)發(fā)影響較小時(shí)間段(表1)的天然泥沙序列作為訓(xùn)練集數(shù)據(jù)。地形數(shù)據(jù)來(lái)源于30 m×30 m分辨率ASTER GDEM(https：∥urs.earthdata.nasa.gov)，通過(guò)處理得到空間尺度、河流水系等。土地覆被數(shù)據(jù)來(lái)源于資源環(huán)境科學(xué)與數(shù)據(jù)中心(http：∥www.resdc.cn)1 km分辨率土地利用遙感監(jiān)測(cè)數(shù)據(jù)集(1980年，2010年)，根據(jù)第二次全國(guó)土地調(diào)查土地利用/覆蓋分類(lèi)體系，按照一級(jí)類(lèi)型進(jìn)行重新合并得到。氣溫、降水?dāng)?shù)據(jù)來(lái)源于國(guó)家地球系統(tǒng)科學(xué)數(shù)據(jù)中心(http：∥www.geodata.cn)利用原始?xì)庀笳军c(diǎn)資料插值生成的1 km×1 km分辨率逐月平均數(shù)據(jù)集(1970—2018年，取支流泥沙資料對(duì)應(yīng)的5—10月數(shù)據(jù))[18]，在該數(shù)據(jù)庫(kù)生成5 km×5 km分辨率逐月歸一化植被指數(shù)(INDV)遙感數(shù)據(jù)集(1982—2018年)。所有數(shù)據(jù)產(chǎn)品均經(jīng)過(guò)嚴(yán)格復(fù)合審查，精度得到保證。

表1 支流水文站數(shù)據(jù)序列

2.2 因子選擇

選取地形、氣候、土地覆被和空間尺度幾類(lèi)因素[6，19-21]的原始數(shù)據(jù)作為因子集(表2)，有助于明晰自然環(huán)境條件與懸移質(zhì)來(lái)沙的直接因果關(guān)系，同時(shí)也降低了后續(xù)建模參數(shù)要求。

表2 選用因子基本概況

2.3 分析方法原理

2.3.1 Spearman秩相關(guān)分析

Spearman秩相關(guān)是一種有效非參數(shù)統(tǒng)計(jì)方法，適用于不是正態(tài)雙變量或總體分布未知的資料[3，11]。計(jì)算公式如下：

(1)

式中：R為Spearman秩相關(guān)系數(shù)；di為變量X第i個(gè)觀測(cè)值xi的秩和變量Y第i個(gè)觀測(cè)值yi的秩之間的差值；n為樣本容量。R越接近1，說(shuō)明兩者間相關(guān)性越大。相關(guān)系數(shù)檢驗(yàn)p<0.05時(shí)相關(guān)性顯著，p<0.01時(shí)相關(guān)性極顯著，其余表示不顯著。

偏相關(guān)方法量化去除控制變量影響后某2個(gè)選定變量之間的關(guān)聯(lián)程度[19，21]，原理是分別對(duì)控制變量和選定變量進(jìn)行回歸，然后計(jì)算這2個(gè)回歸結(jié)果殘差之間的相關(guān)性。

2.3.2 PLSR回歸分析

Q2=1.0-SE/SS

(2)

(3)

(4)

式中：SE為預(yù)測(cè)誤差平方和；SS為殘差平方和；N為成分ti的個(gè)數(shù)。

3 產(chǎn)沙異性歸因分析

3.1 泥沙不平衡分布格局

懸移質(zhì)輸沙量(W)和懸移質(zhì)輸沙模數(shù)(M)是比較不同尺度流域產(chǎn)沙能力的常用指標(biāo)[10，22]。對(duì)于小流域山區(qū)河流，即使環(huán)境條件穩(wěn)定時(shí)產(chǎn)沙隨機(jī)性也較大，即年際變化較大，因此，從多年平均尺度進(jìn)行分析更為準(zhǔn)確。金沙江下游泥沙分配極不均勻，如牛欄江大沙店站(5-DSD)年輸沙量達(dá)1 112萬(wàn)t，遠(yuǎn)高于其他水文站(圖3(a))，但一定程度上是由于其集水面積較大；與之相比圖3(b)的輸沙模數(shù)分布更能反映研究區(qū)域內(nèi)單位面積產(chǎn)沙能力，各地區(qū)之間也更具有可比性。如昭覺(jué)河(1-ZJ)、黑水河(4-NN)的輸沙量遠(yuǎn)小于牛欄江(5-DSD)，但輸沙模數(shù)反而略大于后者。

圖3 支流水文站集水區(qū)實(shí)測(cè)泥沙指標(biāo)分布Fig.3 Distribution of observed sediment index in catchments of tributary stations

3.2 各因子類(lèi)別內(nèi)部降維

單個(gè)因子無(wú)法完整預(yù)報(bào)流域產(chǎn)沙，但使用過(guò)多因子又會(huì)導(dǎo)致重復(fù)解釋的問(wèn)題。因此，在完整地定量化研究因子-產(chǎn)沙關(guān)系前需先進(jìn)行降維處理，即去除代表性弱、關(guān)系不大或是可被替代的因子。開(kāi)展Spearman秩相關(guān)分析得到表2中因子與W、M以及各因子之間的矩陣關(guān)系(圖4,*代表p≤0.05,**代表p≤0.01)，然后從侵蝕輸沙的物理意義出發(fā)，在地形、氣候、土地覆被和空間尺度幾大類(lèi)別中篩選各類(lèi)別的代表性指標(biāo)。

(1) 地形類(lèi)別。各坡度因子(Smean，S8，S18，S25)內(nèi)部之間的相關(guān)系數(shù)高達(dá)0.91～0.97，即研究區(qū)域內(nèi)各坡度占比具有較高的一致性，因此過(guò)多的坡度因子在指標(biāo)信息上是冗余的。S8與M的相關(guān)性最好(R=0.64)，S18次之(R=0.61)，而S25和Smean對(duì)于M的解釋能力較為一般(R=0.53)。由于各集水區(qū)內(nèi)部坡度的分布差異性較大，單一的平均坡度Smean影響產(chǎn)沙的物理機(jī)制相對(duì)不明確是可以預(yù)見(jiàn)的，而中坡度占比S8則擁有預(yù)報(bào)M變化的相對(duì)優(yōu)勢(shì)。

高程因子中Hmax與M相關(guān)性最顯著(R=0.83)，其次是相關(guān)性明顯弱一級(jí)的Hmean(R=0.58)及Hmin(R=-0.12)，且Hmean與Hmax之間本身自相關(guān)(R=0.72)，因此，Hmax可替代Hmean。地形對(duì)于流域產(chǎn)沙的影響由多過(guò)程耦合而成，從侵蝕源頭來(lái)看，高海拔(高Hmax)山區(qū)流域內(nèi)通常也具有極大的落差，對(duì)應(yīng)的是流域內(nèi)坡面更大的潛在勢(shì)能、更低的地表入滲率和更大的產(chǎn)流量，產(chǎn)沙能力相應(yīng)增強(qiáng)[17]。從泥沙輸移的角度來(lái)看，高山地區(qū)眾多的坡面(高S8)減小了侵蝕泥沙顆粒隨后沉積的概率，所以坡度因子直接控制的是流域泥沙的沉積匯過(guò)程。從降維意義而言，可選取S8與Hmax的參數(shù)組合來(lái)替代其他地形因子。

(2) 氣候類(lèi)別。T與M具有極顯著的負(fù)相關(guān)性(R=-0.81)，而P與M的關(guān)系則相對(duì)一般(R=0.59)。氣溫與產(chǎn)沙之間的強(qiáng)關(guān)系存在多方面原因：首先氣溫關(guān)系到降雨概率、雨時(shí)、落區(qū)以及范圍和強(qiáng)度大小等直接影響侵蝕源[21]，如圖4中T與P的負(fù)相關(guān)性(R=-0.70)量化反映了研究區(qū)域內(nèi)明顯的干-熱與冷-濕分布差異(圖1)；其次，T與Hmax的負(fù)相關(guān)性(R=-0.64)也說(shuō)明氣溫通過(guò)地形地理?xiàng)l件關(guān)系到產(chǎn)沙能力。由于P很大程度上能被T反映，因此，選取T代表氣候因子。

(3) 土地覆被類(lèi)別。fa與M相關(guān)性極差，而ff、fg與M的相關(guān)性近似處于同一等級(jí)(ff：R=-0.45，fg：R=0.57)。本質(zhì)上草地和林地都抑制了土壤侵蝕，但草地水土保持能力遠(yuǎn)弱于林地[23]。fg與M的正相關(guān)性很大程度上歸結(jié)于研究區(qū)域特性：土地覆被以林地和草地為主，有14個(gè)集水區(qū)的林地、草地占比總和超過(guò)70%，ff與fg存在一定此消彼長(zhǎng)的現(xiàn)象(R=-0.52)。從土地覆被與其他因子之間聯(lián)系也可看出，Hmax、T與fg分別為正相關(guān)(R=0.42)、負(fù)相關(guān)(R=-0.53)，與ff則相反，即形成了可見(jiàn)的二元空間分布格局：高海拔、低溫情況下草地優(yōu)勢(shì)，林地劣勢(shì)，反之亦然。也就是說(shuō)，fg與M的正相關(guān)性是由抗侵蝕能力更強(qiáng)的林地減少所導(dǎo)致，這也從側(cè)面證明了單因子回歸分析的局限性。因此，選取更符合物理意義的ff代表土地覆被因子。

(4) 空間尺度類(lèi)別。FF與M相關(guān)性極差，因此，只選取相關(guān)性更好的A作為該類(lèi)別的代表。

圖4 水文站集水區(qū)多因子Spearman秩相關(guān)分析Fig.4 Spearman rank correlations between factors

3.3 因子類(lèi)別間降維

基于篩選因子開(kāi)展偏秩相關(guān)性分析，最終分離出能夠定量反映產(chǎn)沙多過(guò)程的驅(qū)動(dòng)因子組合。如圖5所示，輸入S8、Hmax、T、ff、A共計(jì)5個(gè)因子，然后將顯著性檢驗(yàn)結(jié)果p作為選擇控制變量閾值(優(yōu)先p<0.01)，按照偏相關(guān)系數(shù)大小順序逐步篩選出余下最能解釋M的因子變量，直到偏相關(guān)系數(shù)均不通過(guò)顯著性檢驗(yàn)為止(即圖中所有變量p>0.05，為空心圓)，此時(shí)表明余下自變量已不足以解釋因變量變化，輸出6個(gè)驅(qū)動(dòng)因子組合。

Hmax作為單因子對(duì)M變化已有一定解釋能力。相對(duì)地控制住T后，此時(shí)A與M表現(xiàn)出了明顯的正相關(guān)性。類(lèi)似地控制住S8后M與ff之間負(fù)相關(guān)性明顯加強(qiáng)，說(shuō)明流域內(nèi)中等坡地與產(chǎn)沙的強(qiáng)關(guān)系掩蓋了林地對(duì)M的抑制作用，去除前者影響時(shí)后者的相關(guān)性有了大幅提高。在S8組里得到了4個(gè)因子組合S8-ff，S8-T-A，S8-T-ff-A和S8-Hmax-ff。不難看出，圖5一階控制和二階控制中M—A的正相關(guān)性均在控制T后顯現(xiàn)，說(shuō)明當(dāng)山區(qū)流域總體上侵蝕潛能接近時(shí)，流域間存在明顯的M—A尺度效應(yīng)。以上結(jié)論都是傳統(tǒng)的單因子回歸方法所無(wú)法反映的。

一些研究認(rèn)為M—A相關(guān)是因?yàn)殡S著集水面積的擴(kuò)大，類(lèi)似于谷底這種泥沙易沉積的緩坡比例相應(yīng)增加，因此M會(huì)隨著A的增大而減少[10]。但本研究中各級(jí)坡地占比與A沒(méi)有表現(xiàn)出明顯相關(guān)性(R=0.06～0.25，圖4)，由于河道下切劇烈、遍布“V”型河谷，緩坡寬谷段未隨流域面積增大。究其原因是由于A與Hmin負(fù)相關(guān)(R=-0.76)，因此，在控制T的前提下，A增大對(duì)應(yīng)于Hmin減小與流域內(nèi)高差增大，提高侵蝕勢(shì)能與產(chǎn)沙，故此時(shí)M—A呈正相關(guān)。

圖5 偏秩相關(guān)分析及驅(qū)動(dòng)因子組合結(jié)果Fig.5 Partial rank correlation analysis and factor combinations

4 產(chǎn)沙預(yù)測(cè)模型

4.1 PLSR回歸模型

表3 PLSR回歸分析結(jié)果

4.2 多元回歸模型

考慮到M與自然環(huán)境因素之間的非線(xiàn)性關(guān)系，因此，對(duì)比線(xiàn)性模型4#進(jìn)一步建立相似環(huán)境下的多元回歸模型。按照偏相關(guān)分析(圖5)的控制順序，首先測(cè)試S8-M的多種函數(shù)擬合結(jié)果發(fā)現(xiàn)指數(shù)形式擬合優(yōu)度(R2=0.43)最好，以此為基礎(chǔ)進(jìn)一步評(píng)估結(jié)合T、A多種形式的多元回歸效果(圖6(a)—圖6(c))。經(jīng)過(guò)比對(duì)發(fā)現(xiàn)T和A分別取指數(shù)和乘冪形式的擬合效果最佳(圖6(b))，此時(shí)模型足夠解釋約86.4%的M變化性，但各模型在靠近預(yù)測(cè)結(jié)果(M*)極值時(shí)失穩(wěn)(圖6(d))。這是由于泥沙訓(xùn)練集主體含有相當(dāng)部分的重力侵蝕占比，未修正前的模型會(huì)高估少數(shù)支流泥石流、滑坡水平。如當(dāng)輸沙模數(shù)很大時(shí)，重力侵蝕可能已達(dá)到上限，其余侵蝕產(chǎn)沙比例上升，此時(shí)模型預(yù)測(cè)值偏高需向下修正。假定修正參數(shù)k為實(shí)測(cè)值與預(yù)測(cè)值之比，對(duì)k和M*擬合得到修正系數(shù)公式(圖6(e))。模型形式如下：

Mpre=kM*

(5)

式中：Mpre為修正后的預(yù)測(cè)輸沙模數(shù)；k為修正系數(shù)；M*為未修正的預(yù)測(cè)輸沙模數(shù)。其公式如下：

M*=exp(0.7S8-1.56T)A0.28

(6)

(7)

過(guò)去類(lèi)似模型大多止步于統(tǒng)一的線(xiàn)性或冪函數(shù)形式，并沒(méi)有對(duì)模型進(jìn)一步優(yōu)化[15-16]。模型7#R2與ERMS明顯優(yōu)于模型4#結(jié)果。其中，T可反映流域總體潛在侵蝕源大小，S8則主要描述過(guò)程中沉積匯可能，A則與流域內(nèi)泥沙來(lái)源、侵蝕類(lèi)型組成等部分特性有間接聯(lián)系，因此，本模型也能反映流域的主要侵蝕-輸沙特征。

圖6 多元回歸模型擬合修正Fig.6 Scatters of multivariate regression model results

4.3 模型驗(yàn)證與對(duì)比

對(duì)于金沙江下游干流河段，在向家壩、溪洛渡蓄水前(2012—2013年)可近似看作多年沖淤平衡的自然狀態(tài)，即區(qū)間各支流來(lái)沙總和按區(qū)間出、入口水文站輸沙量差值來(lái)估計(jì)。分別選取攀枝花站(及雅礱江桐子林站)、華彈站和屏山站多年平均輸沙量資料作為實(shí)測(cè)值對(duì)以上模型進(jìn)行驗(yàn)證，并通過(guò)2時(shí)段對(duì)比檢驗(yàn)?zāi)Ｐ驮跁r(shí)程上的預(yù)測(cè)精度。表4顯示，華彈—屏山區(qū)間模擬值與實(shí)測(cè)值總體吻合較好。攀枝花—華彈區(qū)間存在輸沙量異常高的小江等支流，考慮到該類(lèi)支流無(wú)水文站泥沙資料，本研究根據(jù)東川泥石流觀測(cè)研究站估算小江年輸沙量來(lái)得到區(qū)間來(lái)沙總量，模擬結(jié)果也與攀枝花—華彈實(shí)測(cè)值接近(表4)。除了區(qū)間來(lái)沙的模型計(jì)算總量符合干流實(shí)測(cè)值外，圖7顯示，模型對(duì)不同時(shí)段劃分后的產(chǎn)沙計(jì)算值也與各支流對(duì)應(yīng)時(shí)段的實(shí)測(cè)值吻合較好。

表4 金沙江下游區(qū)間來(lái)沙量模型驗(yàn)證

金沙江下游目前還沒(méi)有系統(tǒng)性的產(chǎn)沙模型研究，僅有部分針對(duì)土壤侵蝕和產(chǎn)流的模擬[24-25]，不易實(shí)現(xiàn)例如本模型對(duì)區(qū)間來(lái)沙的估算。相比于對(duì)金沙江中游流域無(wú)資料地區(qū)進(jìn)行的類(lèi)似產(chǎn)沙模擬研究[26]，模型7#所需輸入數(shù)據(jù)更易于獲得，輸出結(jié)果表現(xiàn)良好，避免了采用其他模型時(shí)遇到的諸多限制。

現(xiàn)有產(chǎn)沙預(yù)測(cè)模型也難以直接應(yīng)用到本研究區(qū)域。如圖8所示，使用ART模型[27]或更新后的BQART模型[28]往往會(huì)高估中小型流域M值[29]。該類(lèi)模型側(cè)重于大中型流域(流域面積為104～106km2)泥沙通量的估算，而研究區(qū)域18個(gè)集水區(qū)平均面積僅2 780 km2，引入時(shí)還需降尺度處理。相比之下，結(jié)合土壤侵蝕模型(RUSLE)與陳治諫等[30]所得泥沙輸移比(RSD)推算出來(lái)的結(jié)果稍好。但RSD在各個(gè)子流域之間各不相同，由于低產(chǎn)沙帶輸沙難度更大，在計(jì)算中套用目前僅有的流域均值會(huì)高估該區(qū)域M。對(duì)比發(fā)現(xiàn)本研究得到的模型7#表現(xiàn)最好，因此，能夠有效地模擬流域范圍內(nèi)無(wú)資料地區(qū)的M。

圖7 劃分時(shí)段后的模型計(jì)算值驗(yàn)證Fig.7 Model validation for two periods

4.4 泥沙時(shí)空變化格局

根據(jù)模型對(duì)金沙江下游整個(gè)區(qū)間多年平均M的空間分布進(jìn)行推算(圖9)。參照?qǐng)D3(b)發(fā)現(xiàn)，除了已有泥沙資料的少量支流外仍有大量未監(jiān)測(cè)高產(chǎn)沙地帶，主要集中在白鶴灘和溪洛渡兩大巨型水庫(kù)庫(kù)區(qū)，其M明顯高于其余地區(qū)特別是入口近攀枝花和出口近宜賓一帶。

圖9 金沙江下游區(qū)間多年平均輸沙模數(shù)分布Fig.9 Predicted M distribution of the study area

圖10 2007—2018年相對(duì)1970—1986年輸沙模數(shù)變化Fig.10 Variations of M from 1970—1986 to 2007—2018

時(shí)間變化上，由圖10可以看出，大部分地區(qū)M近年比20世紀(jì)七八十年代明顯減小，減小幅度一般在50～300 t/(km2·a)。有學(xué)者統(tǒng)計(jì)得出1991—2005年金沙江下游“長(zhǎng)治”工程(包含下墊面治理及塘堰、攔沙壩等水保工程)水土流失治理減沙效益僅為4.9%，尤其對(duì)重力侵蝕效果不佳，不是支流來(lái)沙減少的主因[31]，即這一趨勢(shì)很大程度上是自然條件變化導(dǎo)致。Li等[9]在地理位置相鄰的西江流域中同樣得到了長(zhǎng)序列下的M減小趨勢(shì)，并且指出升溫刺激了植被長(zhǎng)勢(shì)，INDV增加。經(jīng)本研究統(tǒng)計(jì)，金沙江下游INDV在2個(gè)時(shí)間段間也升高了9.8%～20.7%，研究區(qū)域溫度的升高可能提供了更好的植被發(fā)育條件，對(duì)產(chǎn)沙的抑制作用也因此增強(qiáng)。同時(shí)，中部高產(chǎn)沙區(qū)的減沙程度總體上高于出、入口地區(qū)，即研究區(qū)域產(chǎn)沙的空間不均衡性近年來(lái)有所削弱，若自然條件變化趨勢(shì)不變則有可能進(jìn)一步均衡化，減輕白鶴灘、溪洛渡水庫(kù)淤積壓力。

相關(guān)研究結(jié)果有利于區(qū)間內(nèi)泥沙溯源，進(jìn)而對(duì)干流烏東德-白鶴灘-溪洛渡-向家壩梯級(jí)庫(kù)區(qū)支流來(lái)沙與泥沙淤積的不均勻分布進(jìn)行精細(xì)化研究，如支流河口攔門(mén)沙等，也有利于為梯級(jí)水庫(kù)的調(diào)度運(yùn)行提供參考。

5 結(jié) 論

本文采用金沙江18個(gè)支流水文站多年泥沙資料，分析流域泥沙和地形、氣候、下墊面、空間尺度及人類(lèi)活動(dòng)等因素的關(guān)系，對(duì)金沙江下游區(qū)間多年來(lái)泥沙分布格局、驅(qū)動(dòng)因子以及變化趨勢(shì)進(jìn)行探究。主要結(jié)論如下：

(1) 明晰了能夠解釋金沙江下游來(lái)沙異性格局的多個(gè)自然環(huán)境因子，揭示了各因子對(duì)流域產(chǎn)沙的潛在影響以及各因子自身的變化獨(dú)立性與信息冗余程度，彌補(bǔ)了單因子回歸分析等傳統(tǒng)研究方法的不足。

(2) 分離出多個(gè)關(guān)鍵產(chǎn)沙驅(qū)動(dòng)因子組合，繼而獲得了與類(lèi)似研究相比擬合優(yōu)度更高的因子和產(chǎn)沙關(guān)系模型。其中，以氣溫、8°以上坡度占比、集水面積為變量構(gòu)建的預(yù)測(cè)模型能夠解釋約92%的輸沙模數(shù)變化性，模型參數(shù)需求較低且經(jīng)驗(yàn)證預(yù)測(cè)精度良好，在時(shí)空尺度上都具有穩(wěn)健表現(xiàn)。

(3) 本文模型計(jì)算結(jié)果表明，金沙江下游無(wú)實(shí)測(cè)資料區(qū)間輸沙模數(shù)空間分布極不均衡，在87～1 189 t/(km2·a)間變化，其中高產(chǎn)沙地帶主要集中在白鶴灘和溪洛渡兩大巨型水庫(kù)區(qū)間；同時(shí)，研究區(qū)域輸沙模數(shù)近50 a來(lái)減少約50～300 t/(km2·a)，空間不均衡性有所削弱。本文可為解決類(lèi)似的山區(qū)支流泥沙與水庫(kù)淤積計(jì)算問(wèn)題提供參考。