黃河流域水沙變化趨勢多模型預(yù)測及其集合評(píng)估

胡春宏，張曉明，于坤霞，徐夢(mèng)珍，趙 陽

(1.中國水利水電科學(xué)研究院 流域水循環(huán)模擬與調(diào)控國家重點(diǎn)實(shí)驗(yàn)室，北京 100048；2.西安理工大學(xué) 省部共建西北旱區(qū)生態(tài)水利國家重點(diǎn)實(shí)驗(yàn)室，陜西 西安 710048；3.清華大學(xué) 水沙科學(xué)與水利水電工程國家重點(diǎn)實(shí)驗(yàn)室，北京 100084)

1 前言

水沙變化情勢深刻影響著水資源的開發(fā)利用與社會(huì)可持續(xù)發(fā)展。近幾十年來，隨著水利水保工程的興建并發(fā)揮作用、資源開發(fā)和城鎮(zhèn)化進(jìn)程的持續(xù)推進(jìn)，以及氣候變化影響，全球主要大江大河的水、沙量均發(fā)生了明顯變化[1-3]。據(jù)統(tǒng)計(jì)，全球主要江河1990年代以來126億t的入海沙量減少了30%，其中，47%的河流輸沙量減少、22%的河流徑流量減少、19%的河流兩者都減少[4]。水沙變化如此之大、如此之快，主導(dǎo)因素是什么，未來趨勢如何，一直是科學(xué)界和工程實(shí)踐關(guān)注的重點(diǎn)。水文模型是研究流域水文過程的重要工具，其主要是基于數(shù)學(xué)或物理方程對(duì)復(fù)雜水文過程進(jìn)行概化表達(dá)[5]，具有化繁為簡的優(yōu)點(diǎn)[6]，但水文過程的簡化處理不可避免地使得模擬和預(yù)測結(jié)果存在不確定性。因此，盡管關(guān)于水文模型預(yù)測的研究成果很多，但由于各類模型結(jié)構(gòu)的差異而導(dǎo)致模型輸入數(shù)據(jù)的類型、格式與分辨率以及結(jié)構(gòu)參數(shù)各異，模擬結(jié)果“百家爭鳴”[7-8]。如何提高水沙變化歸因分析和趨勢預(yù)測的可靠性，降低結(jié)果的不確定性，是當(dāng)前研究亟需解決的問題。

目前，利用多模式集合技術(shù)提高評(píng)估結(jié)果的精度并降低其不確定性，受到國內(nèi)外研究者高度重視。1999年，Krishnamurti等[9]首次提出了超級(jí)集合預(yù)報(bào)模式，該方法利用集合預(yù)報(bào)的多個(gè)參與模式的表現(xiàn)效果來確定其權(quán)重系數(shù)，進(jìn)而對(duì)預(yù)報(bào)結(jié)果進(jìn)行修正，以獲得最佳的預(yù)報(bào)性能。Cartwright等[10]基于9個(gè)模式參與降水量集合預(yù)報(bào)，發(fā)現(xiàn)超級(jí)集合預(yù)報(bào)方法的預(yù)報(bào)效果優(yōu)于任意單一模式。Krishnamurti等[11]對(duì)降水日變化的研究也印證了多模式超級(jí)集合方法能有效降低評(píng)估結(jié)果的不確定性，并且預(yù)報(bào)效果較單個(gè)模型更好。伍清等[12]基于11個(gè)氣候模式對(duì)中國西南地區(qū)的氣溫進(jìn)行了集合預(yù)報(bào)，并證明集合預(yù)報(bào)可以有效降低模擬誤差。Li等[13]的研究表明多模型集合可以有效提高徑流月變化的估測精度，具有良好的開發(fā)前景?？梢?，多模型集合預(yù)報(bào)在提高模擬結(jié)果科學(xué)性和可信度方面具有顯著優(yōu)勢。近年來，不斷發(fā)展和完善的區(qū)間預(yù)報(bào)方法[14-15]，成為水文學(xué)領(lǐng)域評(píng)估水文事件不確定性的重要手段。

借鑒以往研究者的集合預(yù)報(bào)思想，結(jié)合區(qū)間預(yù)報(bào)方法，本文試圖在識(shí)別各類模擬方法適用性和剖析研究成果不確定性來源的基礎(chǔ)上，提出針對(duì)流域水沙變化歸因分析和趨勢預(yù)測的集合評(píng)估方法，并圍繞黃河流域水沙預(yù)測現(xiàn)有方法的預(yù)測結(jié)果開展集合評(píng)估，以期科學(xué)評(píng)價(jià)未來氣候和下墊面變化下黃河流域水沙量，進(jìn)而推動(dòng)水沙過程集合預(yù)報(bào)技術(shù)的發(fā)展。

2 研究區(qū)概況與研究方法

2.1 研究區(qū)概況黃河干流全長5464 km，流域面積79.5萬km2(包括內(nèi)流區(qū))，內(nèi)蒙古河口鎮(zhèn)(水文控制站為頭道拐)以上為上游，河口鎮(zhèn)至鄭州桃花峪(水文控制站為花園口)為中游，桃花峪以下為下游。黃土高原水土流失面積45.4萬km2，占黃河流域水土流失總面積46.5萬km2的98%。黃土高原向黃河貢獻(xiàn)了97%的泥沙。黃土高原與黃河流域位置圖及干流主要水文站分布如圖1所示。黃河中游流域以半干旱大陸性季風(fēng)氣候?yàn)橹?，年均降水量由西北部?50 mm遞增至東南部的650 mm。由于黃河中游大部分支流流經(jīng)水土流失嚴(yán)重的黃土高原地區(qū)，是黃河粗泥沙的主要來源。自1980年代以來，國內(nèi)外學(xué)者圍繞黃河中游開展了大量水沙變化機(jī)理與水沙模型構(gòu)建及預(yù)測的研究[8，16]，本研究以黃河流域潼關(guān)以上區(qū)域?yàn)閷?duì)象，開展典型流域水沙變化歸因不確定性分析與水沙預(yù)測集合評(píng)估研究。

圖1 研究區(qū)域

2.2 黃河流域水沙預(yù)測方法黃河水沙變化預(yù)測成果直接關(guān)系到黃河治理開發(fā)方略制定，因此，1980年代至今一直有國家計(jì)劃課題資助開展相關(guān)研究，并取得了相應(yīng)條件下的預(yù)測結(jié)果。其中，由水利部主導(dǎo)的黃河水沙變化研究基金第一期提出“水文法”和“水保法”[17]，這類方法選取的變量較少，是一種典型的基于歷史觀測數(shù)據(jù)的回歸模型。1950年代以來，隨著對(duì)下滲理論、土壤水運(yùn)動(dòng)理論以及河道水力學(xué)理論等方面認(rèn)識(shí)的不斷加深，開始發(fā)展了概念性模型[18]和分布式水文模型[19-20]。隨著信息技術(shù)的發(fā)展，機(jī)器學(xué)習(xí)方法成為挖掘不同特征之間復(fù)雜關(guān)系的有效方法，可真實(shí)的對(duì)自然規(guī)律做出預(yù)測，也成為流域水沙預(yù)測的重要方法[21]。本研究選用流域水沙變化歸因和預(yù)測應(yīng)用較廣的方法作為集合評(píng)估對(duì)象開展研究。

(1)經(jīng)驗(yàn)?zāi)Ｐ?。包括水文?LR)和雙累積曲線法(DMC)等，通過基準(zhǔn)期和人類活動(dòng)影響期建立降雨-徑流(輸沙)關(guān)系或降雨-徑流(輸沙)連續(xù)累計(jì)值關(guān)系，來分離人類活動(dòng)與氣候變化對(duì)水文過程的影響，預(yù)測流域水沙量。該方法計(jì)算簡單，但物理意義缺失，在數(shù)據(jù)精度受限下是大尺度流域水沙變化預(yù)測的有效方法。本研究采用1954—2018年(日歷年)共65 a的樣本作為數(shù)據(jù)集，其中1954—1989年為基準(zhǔn)期，1990—2018年為驗(yàn)證期。

(2)彈性系數(shù)法。徑流(輸沙)變化率對(duì)氣候因子變化率的比值[22]，以基于Budyko假設(shè)的水熱耦合平衡方程為核心，能實(shí)現(xiàn)不同影響因素的定量分析，相較于經(jīng)驗(yàn)?zāi)Ｐ臀锢頇C(jī)制更強(qiáng)，但該方法不能解決各輸入項(xiàng)之間有潛在聯(lián)系的問題，且流域特征概化，無法實(shí)現(xiàn)精細(xì)化的模擬和預(yù)測。

(3)機(jī)器學(xué)習(xí)模型(ML)?；趯?duì)徑流輸沙演變規(guī)律的大數(shù)據(jù)分析，挖掘徑流形成的驅(qū)動(dòng)要素，提取徑流輸沙預(yù)測因子，借助機(jī)器學(xué)技術(shù)，建立徑流輸沙預(yù)測模型。該模型屬于考慮多因子同徑流輸沙的關(guān)系經(jīng)驗(yàn)?zāi)Ｐ秃秃谙淠Ｐ?，并?duì)流域特征進(jìn)行了概化。本研究中機(jī)器學(xué)習(xí)方法使用水土保持措施面積數(shù)據(jù)和歸一化植被指數(shù)(NDVI)來反映人類活動(dòng)對(duì)水沙的影響。

(4)GAMLSS(Generalized Additive Models for Location，Scale and Shape)。GAMLSS是位置、尺度、形狀的廣義可加模型，模型能夠非常靈活地建立概率分布參數(shù)與協(xié)變量之間的函數(shù)關(guān)系，屬于黑箱概化模型。本文將水土保持措施面積作為影響徑流量和輸沙量的協(xié)變量。

(5)SWAT(Soil and Water Assessment Tool)。SWAT是一種以日為時(shí)間步長運(yùn)行的連續(xù)時(shí)間分布式流域水文模型，通過修正的徑流曲線數(shù)法(SCS法)和Green &Ampt下滲法估算地表徑流，運(yùn)用修正的通用土壤流失方程計(jì)算每個(gè)水文響應(yīng)單元的水力侵蝕[23]。該模型可以長時(shí)間連續(xù)計(jì)算且效率高，適用多種土壤類型，但模型計(jì)算方程多為經(jīng)驗(yàn)公式。

(6)流域水沙動(dòng)力過程模型(DWN)。DWN是由降雨產(chǎn)流模型、坡面產(chǎn)沙模型、溝坡重力侵蝕模型及溝道水沙演進(jìn)與沖淤模型四個(gè)基礎(chǔ)模型組成，并在流域單元的基礎(chǔ)上結(jié)合河網(wǎng)水系進(jìn)行水沙動(dòng)力學(xué)過程模擬的分布式模型[20]。模型動(dòng)力學(xué)機(jī)理完善，雙層率定技術(shù)解決多參數(shù)率定的計(jì)算效率問題，模型復(fù)雜度高；但模型參數(shù)率定困難，計(jì)算量大，數(shù)據(jù)獲取難度也大。

(7)多因子驅(qū)動(dòng)的黃河流域分布式水循環(huán)模型(Multi-factors Driven Water Processes Model in Yellow River Basin，MFD-WESP)。MFD-WESP以流域二元水循環(huán)WEP- L模型為基礎(chǔ)，增加了黃河源區(qū)凍土水熱耦合模塊，并將匯流過程由“坡面-河道”系統(tǒng)改進(jìn)為“坡面-溝壑-河道”系統(tǒng)。該模型可考慮源區(qū)凍土、黃土區(qū)溝壑產(chǎn)匯流和平原區(qū)取用水等的綜合影響[24]。

(8)HydroTrend。由Syvitski等開發(fā)，是基于水量平衡原理和流域長期產(chǎn)沙經(jīng)驗(yàn)?zāi)Ｐ偷募偸剿哪Ｐ?，最適于在年代尺度以上研究氣候變化和人類活動(dòng)(尤其是土地利用和水庫)對(duì)河流流量和輸沙量的影響[25]。模型依賴的流域信息少，能夠在宏觀的角度提供多年平均的變化趨勢。本研究直接使用由陳蘊(yùn)真[26]改進(jìn)的HydroTrend4Yellow模型開展沙量預(yù)測。

(9)產(chǎn)沙指數(shù)模型?；谶b感技術(shù)提取黃土高原不同侵蝕類型區(qū)林草、梯田等空間分布信息，構(gòu)建適用于各類型區(qū)林草梯田有效覆被率-產(chǎn)沙指數(shù)關(guān)系，來推算黃河流域主要產(chǎn)沙區(qū)沙量[27]。

(10)人工智能模型?；诩蓪W(xué)習(xí)中的極端梯度提升樹(XGBoost)算法設(shè)計(jì)，構(gòu)建大數(shù)據(jù)驅(qū)動(dòng)的黃河流域徑流與侵蝕產(chǎn)沙智能預(yù)測模型[28]，屬于機(jī)器學(xué)習(xí)模型中的一種。

(11)BP神經(jīng)網(wǎng)絡(luò)模型。基于黃河流域降水和干流水文站水沙資料，構(gòu)建的以年降雨量預(yù)測值的累計(jì)值和年徑流量累計(jì)值作為輸入層，輸沙量年累計(jì)值作為輸出層的BP網(wǎng)絡(luò)神經(jīng)模型，預(yù)測未來入黃沙量。

2.3 數(shù)據(jù)來源本研究構(gòu)建了流域水沙變化歸因分析和模擬預(yù)測的多方法的標(biāo)準(zhǔn)化數(shù)據(jù)庫。如表1所示，降水?dāng)?shù)據(jù)采用了水文年鑒、國家級(jí)氣象站數(shù)據(jù)、CMORPH衛(wèi)星融合數(shù)據(jù)和氣候模式數(shù)據(jù)。氣候模式數(shù)據(jù)：GCM(Global Climate Model)數(shù)據(jù)使用融合模式比較計(jì)劃第五階段CMIP5模式中的9個(gè)不同模型CMCC-CMS 、GFDL-ESM2M、IPSL-CM5A-LR、CNRM-CM5、CSIRO-Mk3-6-0、EC-EARTH_QM、EC-EARTH_RCM、MIROC-ESM-CHEM和NorESM1-M，獲取2021—2070年的日降水、月氣溫和月蒸發(fā)數(shù)據(jù)。周文翀等[29]利用格點(diǎn)化觀測資料，對(duì)CMIP5提供的18個(gè)全球氣候模式在黃河流域的模擬能力進(jìn)行評(píng)估，結(jié)果表明MIROC-ESM-CHEM、CSIRO-Mk3-6-0、NorESM1-M、CNRM-CM5和EC-EARTH，5個(gè)全球氣候模式綜合評(píng)分較優(yōu)，其中EC-EARTH模式包括統(tǒng)計(jì)降尺度EC-EARTH_QM和動(dòng)力降尺度EC-EARTH_RCM兩種模式?？紤]到GFDL-ESM2M、IPSL-CM5A-LR、CMCC-CMS能捕捉到極端天氣[30]，因此本研究預(yù)測未來水沙時(shí)選取了這9種氣候模式。未來植被條件基于2000—2018年MODIS NDVI數(shù)據(jù)，其時(shí)間分辨率為16 d，空間分辨率為500 m，采取最大值合成法(MVC)獲得逐年NDVI值，以2016—2018三年的平均NDVI值為基準(zhǔn)，根據(jù)計(jì)算的Theil Sen斜率值，后續(xù)預(yù)測中假定增長率不變，每隔五年預(yù)測一次年及年內(nèi)各個(gè)月份的NDVI值。未來淤地壩布局根據(jù)水利普查數(shù)據(jù)及未來淤積庫容公式推算。未來梯田分布根據(jù)《黃河流域生態(tài)保護(hù)和高質(zhì)量發(fā)展規(guī)劃綱要》新要求，統(tǒng)籌各侵蝕分區(qū)水土流失類型、地形、降雨等因素，為滿足黃土高原農(nóng)業(yè)人口發(fā)展、耕地需求、建設(shè)潛力和相關(guān)規(guī)劃要求推算。

表1 研究數(shù)據(jù)

3 研究結(jié)果與分析

3.1 流域水沙變化歸因-預(yù)測不確定性解析

3.1.1 流域水沙變化歸因中的不確定性 黃河流域年均輸沙量從1919—1959年16億t/a銳減至2010—2020年1.83億t/a，減幅89%。為探析黃河水沙變化機(jī)理及影響因素貢獻(xiàn)率，自1980年代以來，國內(nèi)學(xué)者運(yùn)用多種方法開展大量研究。如圖2所示，為2012年以來針對(duì)皇甫川流域開展的水沙變化影響因素貢獻(xiàn)率研究結(jié)果[31]。其中，研究成果中40%采用經(jīng)驗(yàn)?zāi)Ｐ头?，而基于物理過程模型的貢獻(xiàn)率研究較少。相對(duì)于徑流，經(jīng)驗(yàn)?zāi)Ｐ偷难芯拷Y(jié)果變化區(qū)間較大(氣候變化貢獻(xiàn)率為10%～50%)，基于Budyko假設(shè)的彈性系數(shù)法得到的氣候變化平均貢獻(xiàn)率最小，為9%。相對(duì)于輸沙，6種水文法求得的氣候變化平均貢獻(xiàn)率為35%，5種雙累積曲線法為29%，水文法的結(jié)果略高于雙累積曲線法。顯然，造成研究結(jié)果不一致的主要原因在于研究時(shí)期、基準(zhǔn)期、研究方法等不同。因此，黃河流域水沙變化影響因素的貢獻(xiàn)率仍未形成統(tǒng)一認(rèn)識(shí)，影響了區(qū)域治理措施合理布局的科學(xué)指導(dǎo)，亟需集合評(píng)估方法來評(píng)價(jià)各類成果的可信度。

注：LR表示水文法，DMC表示雙累積曲線法，SCR表示SCRAQ法，NP-E表示非參數(shù)式彈性系數(shù)法，B-E表示基于Budyko假設(shè)的彈性系數(shù)法，DM表示Decomposition Method，P-P表示基于過程的物理模型法。圖2 皇甫川流域徑流和輸沙變化歸因既有成果對(duì)比解析

3.1.2 流域水沙模擬精度評(píng)價(jià)的不確定性 流域內(nèi)降雨產(chǎn)流產(chǎn)沙是一個(gè)復(fù)雜的過程，歷來備受水文學(xué)者的關(guān)注?；诂F(xiàn)有降雨-水沙經(jīng)驗(yàn)公式分析表明，降雨-水沙經(jīng)驗(yàn)公式在形式、變量構(gòu)成、數(shù)據(jù)要求、模擬精度等方面均具有多樣性，并無固定的輸入要求，由此導(dǎo)致水沙模擬的輸出亦具有差異性。本研究收集整理了黃河流域水沙模擬研究在使用不同精度評(píng)價(jià)指標(biāo)下的預(yù)測結(jié)果誤差范圍及可信度數(shù)據(jù)[32]，如圖3所示，繪制箱線圖進(jìn)行現(xiàn)有方法的精度分析。

圖3 現(xiàn)有水沙預(yù)測模型模擬精度的可信度分析

評(píng)價(jià)模型模擬精度的常用指標(biāo)是相關(guān)系數(shù)r、確定性系數(shù)R2、納什效率系數(shù)NSE和相對(duì)誤差RE，其中r、R2和NSE越接近1、RE越小表明模型模擬精度越高。圖3顯示不同模型在黃河流域不同子流域水沙模擬中精度判別指標(biāo)對(duì)比情況，其中NSE和RE統(tǒng)計(jì)結(jié)果顯示徑流和輸沙模擬率定期精度整體高于驗(yàn)證期，但是r和R2統(tǒng)計(jì)結(jié)果顯示徑流和輸沙在驗(yàn)證期的模擬精度優(yōu)于率定期。出現(xiàn)這種情況的主要原因：一是不同模型評(píng)價(jià)指標(biāo)由于評(píng)價(jià)側(cè)重點(diǎn)不同，對(duì)模型精度的評(píng)定結(jié)果不完全一致；二是絕大部分研究僅用1個(gè)或2個(gè)評(píng)價(jià)指標(biāo)來評(píng)價(jià)模型的模擬精度，圖3中不同評(píng)價(jià)指標(biāo)的樣本容量差異很大，評(píng)價(jià)率定期徑流模擬精度的NSE樣本容量最大48，評(píng)價(jià)驗(yàn)證期輸沙模擬精度的R2樣本容量最小，僅為5，因此不同評(píng)價(jià)指標(biāo)結(jié)果的代表性不同；三是部分研究未劃分率定期和驗(yàn)證期，驗(yàn)證期樣本容量遠(yuǎn)小于率定期。

因此，不同模型的水沙模擬精度差別較大，且在不同模型評(píng)價(jià)指標(biāo)下得到的評(píng)價(jià)結(jié)果精度之間存在差異，僅采用一種評(píng)價(jià)指標(biāo)評(píng)價(jià)模型精度不全面。因此，在黃河流域水沙變化趨勢預(yù)測時(shí)，不僅要采用多個(gè)指標(biāo)評(píng)價(jià)不同水沙預(yù)測模型在率定期模擬精度，更要全面評(píng)價(jià)模型驗(yàn)證期模擬精度。

3.2 流域水沙變化預(yù)測集合評(píng)估

3.2.1 集合評(píng)估技術(shù)的提出 黃河流域是水沙-地貌-生態(tài)多過程耦合的動(dòng)態(tài)系統(tǒng)?；谏瞎?jié)分析表明，在預(yù)測流域水沙變化時(shí)，由于復(fù)雜系統(tǒng)內(nèi)多因素的耦合驅(qū)動(dòng)及其時(shí)空的非均質(zhì)性變化，各類模型在表征這種群體性、時(shí)滯性和非線性關(guān)系時(shí)，可能存在某種假定、某些算法的忽視或某些參數(shù)的概化表達(dá)，導(dǎo)致模型輸出結(jié)果存在不確定性，同時(shí)觀測誤差也會(huì)引發(fā)模型精度評(píng)價(jià)的不確定性。如圖4所示，所有模型模擬結(jié)果不確定均由數(shù)據(jù)輸入、模型參數(shù)、模型結(jié)構(gòu)以及觀測等不確定性引發(fā)，這也導(dǎo)致各類方法和模型在應(yīng)用時(shí)均存在不同程度的局限性。

圖4 水文模型模擬/預(yù)測-診斷中的不確定性診斷

為確保流域水沙模擬和預(yù)測的結(jié)果可信，需構(gòu)建基于各類方法輸入、參數(shù)、結(jié)構(gòu)和輸出等綜合評(píng)價(jià)的集合評(píng)估技術(shù)，科學(xué)確定輸出值及其置信區(qū)間。如圖5所示，為流域水沙變化趨勢預(yù)測集合評(píng)估技術(shù)路線。即，基于模型不確定性傳遞的評(píng)價(jià)，篩選影響模擬結(jié)果的可表征輸入、參數(shù)和模型結(jié)構(gòu)等的不確定性指標(biāo)與精度指標(biāo)，構(gòu)建綜合評(píng)價(jià)矩陣，基于TOPSIS模型或模糊決策理論等分析不同評(píng)估指標(biāo)的敏感性，并確定各類方法模擬精度與不確定性指標(biāo)權(quán)重。圍繞年尺度、不同降雨水平年、多評(píng)價(jià)指標(biāo)、多評(píng)估方法的多維度集合評(píng)估，提出既有方法模擬精度的概率函數(shù)與不確定性區(qū)間分布；基于極大似然法思想對(duì)各類方法加權(quán)平均，利用BMA模型中反距離加權(quán)與貝葉斯加權(quán)，并根據(jù)既有方法模擬精度和不確定性區(qū)間權(quán)重分布，搭建多元耦合的非線性水沙變化最優(yōu)評(píng)價(jià)模型，集合預(yù)報(bào)典型流域未來水沙變化，并針對(duì)水沙變化預(yù)測置信度識(shí)別體系，確定預(yù)測結(jié)果的置信度及置信區(qū)間。

圖5 水沙變化趨勢預(yù)測集合評(píng)估技術(shù)路線

集合評(píng)估不同于集合預(yù)報(bào)，前者是在對(duì)各類預(yù)測結(jié)果綜合評(píng)價(jià)基礎(chǔ)上提出最接近真值的值，后者是對(duì)同一系統(tǒng)分段模擬來集合輸出一個(gè)結(jié)果，后者是前者的對(duì)象。水沙變化不同模擬和預(yù)測方法的集合評(píng)估，有助于突破黃河水沙趨勢預(yù)測難于達(dá)成一致共識(shí)的桎梏，為科學(xué)預(yù)測未來常態(tài)水沙情勢提供重要依據(jù)[8]。

3.2.2 黃河流域典型水沙預(yù)測方法集合評(píng)估 現(xiàn)有流域水沙預(yù)測方法種類較多，如2.2節(jié)所列，MFD-WESP只模擬預(yù)測黃河流域徑流量，產(chǎn)沙指數(shù)模型只預(yù)測流域產(chǎn)沙量，HydroTrend輸出結(jié)果為多年均值。因此，考慮各類方法輸入的數(shù)據(jù)格式標(biāo)準(zhǔn)一致，輸出的結(jié)果為長序列年尺度徑流量和輸沙量，選用LR、ML、GAMLSS和SWAT作為代表性模型，基于黃河流域潼關(guān)站水沙預(yù)測值和實(shí)測值(Observation，簡稱OBS)比較，來集合評(píng)估各類方法的適用性?？紤]到模擬流域的空間尺度限制和模擬結(jié)果的精度，上述模型分別對(duì)皇甫川、延河以及無定河三支流的歷史徑流輸沙給出預(yù)測。其中對(duì)于黃河全流域推算，則通過非線性擬合建立1960—2018年潼關(guān)水沙量與三支流水沙量的冪函數(shù)關(guān)系來求算。如圖6所示，由于黃河流域有“水沙異源”的特點(diǎn)，黃河上游為主要產(chǎn)流區(qū)，黃河中游為主要產(chǎn)沙區(qū)，皇甫川、延河和無定河三支流流域均屬于黃河中游主要產(chǎn)沙區(qū)，故三支流流域的年輸沙量之和與潼關(guān)年輸沙量相關(guān)性顯著，而徑流量的相關(guān)關(guān)系則相對(duì)較弱，進(jìn)而使得推算的徑流量準(zhǔn)確性低于輸沙量。

圖6 黃河流域潼關(guān)站-三支流之和的徑流/輸沙相關(guān)關(guān)系

如圖7所示，為黃河潼關(guān)站不同模型水、沙模擬值與實(shí)測值比較。結(jié)果表明：各模型對(duì)年徑流、年輸沙的模擬在2000年前基準(zhǔn)期表現(xiàn)較好，2000年之后的變化期則差異分化，特別是年輸沙的模擬值存在較大差異；GAMLSS相對(duì)其他模型，年徑流和年輸沙的模擬效果均表現(xiàn)優(yōu)異。

圖7 各模型黃河潼關(guān)站年徑流和年輸沙模擬值與實(shí)測值比較

在上述多方法模擬結(jié)果集合評(píng)估的基礎(chǔ)上，搭建了基于貝葉斯模型加權(quán)平均方法的水文模型，即BMA模型，集合預(yù)報(bào)了黃河流域徑流量和輸沙量，并開展對(duì)以上四類方法與BMA模型進(jìn)行綜合評(píng)價(jià)。鑒于各方法對(duì)年徑流、年輸沙模擬在2000年前后的差異，采用確定性模型精度評(píng)價(jià)指標(biāo)納什效率系數(shù)NSE、相對(duì)誤差RE、一致性指數(shù)D、相對(duì)均方根誤差RMSE、相對(duì)平方均方誤差MSESQ、相對(duì)對(duì)數(shù)均方誤差MSELN、校正決定性系數(shù)adjustedR2共7個(gè)，對(duì)各模型模擬序列所劃分的基準(zhǔn)期(2000年之前)和變化期(2001—2018年)進(jìn)行了相應(yīng)評(píng)價(jià)指標(biāo)的計(jì)算?？紤]到評(píng)價(jià)指標(biāo)的量綱不統(tǒng)一，將評(píng)價(jià)指標(biāo)進(jìn)行了無量綱化處理，即將評(píng)價(jià)指標(biāo)的實(shí)際量值轉(zhuǎn)化為[0，1]區(qū)間上的無量綱數(shù)，結(jié)果如圖8所示。年徑流模擬方面，GAMLSS相對(duì)表現(xiàn)最好；年輸沙模擬方面，GAMLSS和ML優(yōu)于LR和SWAT，BMA模型的模擬精度介于單個(gè)模型之間。

圖8 各模型模擬評(píng)價(jià)指標(biāo)隸屬度計(jì)算值統(tǒng)計(jì)圖

根據(jù)上述結(jié)果，按照制定的模型集合評(píng)估準(zhǔn)則對(duì)5種方法在輸入、模型結(jié)構(gòu)和輸出三方面，進(jìn)行數(shù)據(jù)需求、物理基礎(chǔ)、應(yīng)用效率、輸出尺度和預(yù)測精度五指標(biāo)的評(píng)價(jià)。雖然DWN在模擬水沙過程時(shí)對(duì)輸入數(shù)據(jù)的標(biāo)準(zhǔn)要求高、參數(shù)率定困難，但作為物理意義明確、動(dòng)力學(xué)機(jī)理完善的分布式模型，流域水沙預(yù)測仍是必要手段。通過模型各方面排序賦分，模型評(píng)分結(jié)果如圖9所示。不同的模型在不同的需求方面表現(xiàn)不一，在實(shí)際應(yīng)用中，應(yīng)基于模擬條件和輸出成果精度等多維度權(quán)衡來篩選模型。

圖9 模型集合評(píng)估評(píng)分雷達(dá)圖

3.3 黃河流域未來50 a水沙變化趨勢預(yù)測集合評(píng)估黃河水沙受自然過程與人工措施疊加影響，加之氣候與下墊面等邊界條件的不斷變化，黃河水沙預(yù)測結(jié)果存在較大不確定性。2000年以前水沙預(yù)測方法主要為水文法和水保法，受措施統(tǒng)計(jì)類型、系列長短與數(shù)據(jù)標(biāo)準(zhǔn)化程度等影響，且流域下墊面變化相對(duì)劇烈，潼關(guān)站沙量預(yù)測成果之間的差異很大，介于5億t/a～21億t/a之間[8]。2000年以來，隨著黃河流域生態(tài)建設(shè)與水土保持大力推進(jìn)、對(duì)水沙變化機(jī)理認(rèn)識(shí)不斷深入、水沙變化預(yù)測的數(shù)據(jù)處理技術(shù)和模擬技術(shù)等日臻成熟，黃河未來水沙變化預(yù)測成果間的差異逐步縮小，潼關(guān)站預(yù)測沙量介于3億t/a～8億t/a之間，但社會(huì)爭議仍較大[8]。

基于以上評(píng)價(jià)分析，篩選典型模型作為集合評(píng)估對(duì)象來預(yù)測黃河潼關(guān)站未來50 a水沙量。由于水文統(tǒng)計(jì)模型在對(duì)水保措施未來持續(xù)開展所發(fā)揮的減水減沙效益考慮不足，在2000年以后估算得到徑流量和輸沙量均偏大；考慮了水土保持措施的GAMLSS和ML在基準(zhǔn)期和變化期模擬精度均較高；SWAT描述了坡面上的降水產(chǎn)流、徑流產(chǎn)沙和溝道水沙輸移等物理過程，因此在黃土高原典型流域的水沙過程模擬運(yùn)用比較廣；DWN通過子流域單元?jiǎng)澐郑紤]了氣候、地形地貌、土壤植被等的空間差異性，并通過河網(wǎng)水沙演進(jìn)與河道沖淤計(jì)算，重現(xiàn)流域內(nèi)任何河道斷面和流域出口的水沙過程，物理基礎(chǔ)堅(jiān)實(shí)。因此，在對(duì)未來黃河水沙預(yù)測時(shí)，水文統(tǒng)計(jì)模型的預(yù)測結(jié)果不納入集合評(píng)估樣本，其余均用于集合評(píng)估對(duì)象。

如圖10所示，在下墊面和水土保持措施等標(biāo)準(zhǔn)化數(shù)據(jù)及不同氣候模式輸入下相同模型年徑流和年輸沙預(yù)測的結(jié)果差異較大，相同氣候模式不同模型預(yù)測的變化趨勢一致，如GAMLSS和DWN均顯示CNRM-CM5_QM氣候模式條件下年輸沙量的預(yù)測結(jié)果較大，CSIRO-Mk3-6-0_QM氣候模式條件下年輸沙量的預(yù)測結(jié)果小，ML和DWN在IPSL-CM5A-LR氣候模式下年均輸沙量模擬結(jié)果相當(dāng)。另外，不同模型因?yàn)樽陨斫Y(jié)構(gòu)的局限性導(dǎo)致個(gè)別年份預(yù)測值異常，而開展多模型預(yù)測結(jié)果的集合評(píng)估，其目的也正是用來消減因模型結(jié)構(gòu)引起的模擬結(jié)果不確定性。為此，針對(duì)以上SWAT、GAMLSS、ML、DWN模型對(duì)潼關(guān)站未來50 a(2021—2070年)水沙變化趨勢預(yù)測，并結(jié)合“十三五”重點(diǎn)研發(fā)專項(xiàng)“黃河水沙變化機(jī)理與趨勢預(yù)測”項(xiàng)目組采用的MFD-WESP、HydroTrend、產(chǎn)沙指數(shù)模型、人工智能模型、BP神經(jīng)網(wǎng)絡(luò)模型等五類模型關(guān)于黃河潼關(guān)站未來50 a水沙變化趨勢預(yù)測[33]，最后對(duì)上述9類模型預(yù)測結(jié)果集合評(píng)估，如圖11所示。潼關(guān)站未來50 a年均徑流量為240億m3，年徑流量90%的置信區(qū)間范圍為[164，328]億m3。未來50 a年均輸沙量為和2.45億t，年輸沙量90%的置信區(qū)間范圍是[0.79，5.12]億t。近期王光謙等在假定現(xiàn)有水利水保措施攔沙效率不變的條件下，輸入CMIP5-RCP4.5情景預(yù)測未來50 a潼關(guān)站年均徑流量為247億m3，未來50 a均輸沙量在3億t/a左右[34]；“十三五”重點(diǎn)研發(fā)計(jì)劃專項(xiàng)“黃河下游河道與灘區(qū)治理”項(xiàng)目組[35]預(yù)測潼關(guān)站未來50 a的年均輸沙量為2.2億t以下，年均徑流量為300億m3～350億m3；李雅娟等預(yù)測潼關(guān)站 RCP4.5情景下未來50 a的徑流量為213.3億m3，輸沙量為3.42億t[36]。以上學(xué)者的研究進(jìn)一步驗(yàn)證了黃河流域水沙多模型預(yù)測集合評(píng)估的成果。綜合以上分析，作者認(rèn)為黃河流域未來50 a的徑流量年均約為300億m3左右，輸沙量約為3億t左右。

圖10 潼關(guān)2021—2070年不同方法的年輸沙量預(yù)測值

圖11 潼關(guān)未來年均流量和年均輸沙量預(yù)測

4 結(jié)論

自1980年代中期以來，黃河水沙情景發(fā)生顯著變化，直接影響著黃河重大水沙調(diào)控工程布局與流域生態(tài)保護(hù)工程規(guī)劃設(shè)計(jì)。為此，眾多學(xué)者圍繞黃河流域水沙變化歸因與黃河水沙變化趨勢等關(guān)鍵科學(xué)問題開展了卓有成效的研究。然而，黃河水沙變化受氣候、生態(tài)、地質(zhì)地貌和強(qiáng)人類活動(dòng)等多界面、多系統(tǒng)和多因素耦合影響，現(xiàn)有流域水沙變化歸因和預(yù)測方法及模型因輸入、結(jié)構(gòu)、輸出等差異而導(dǎo)致研究結(jié)果有巨大不確定性，難以形成統(tǒng)一認(rèn)識(shí)。本研究通過探討流域水沙變化歸因和預(yù)測結(jié)果的不確定性來源，提出了多模型集合評(píng)估技術(shù)，并基于統(tǒng)一標(biāo)準(zhǔn)化的流域輸入數(shù)據(jù)，辨識(shí)現(xiàn)有黃河流域水沙預(yù)測模型的適用性，進(jìn)而集合評(píng)估黃河流域水沙變化趨勢。

(1)系統(tǒng)梳理了黃河流域水沙變化歸因和預(yù)測既有方法，目前關(guān)于典型流域水沙變化歸因方法應(yīng)用中水文法、雙累積曲線法和彈性系數(shù)法等經(jīng)驗(yàn)?zāi)Ｐ驼急雀?，物理過程模型應(yīng)用的比較少，且因數(shù)據(jù)需求、構(gòu)成變量及模擬精度評(píng)價(jià)方法等差異性等引起方法或模型測算結(jié)果巨大的不確定性，如經(jīng)驗(yàn)?zāi)Ｐ陀?jì)算得到的皇甫川流域水沙變化貢獻(xiàn)率中氣候變化占比從10%至50%不等，差異顯著。同時(shí)基于不同評(píng)價(jià)指標(biāo)得到的模型精度評(píng)價(jià)結(jié)果存在差異，僅采用一種評(píng)價(jià)指標(biāo)評(píng)價(jià)模型精度不全面。

(2)明確了模型在對(duì)水沙模擬時(shí)因數(shù)據(jù)輸入、模型參數(shù)、模型結(jié)構(gòu)以及觀測等不確定性引發(fā)的其在應(yīng)用中的局限性，為此構(gòu)建了模型適用性判別準(zhǔn)則和評(píng)價(jià)技術(shù)以及水沙變化多模型集合評(píng)價(jià)-多結(jié)果加權(quán)融合-BMA集合預(yù)測的集合評(píng)估技術(shù)，并基于數(shù)據(jù)需求、物理基礎(chǔ)、應(yīng)用效率、輸出尺度和預(yù)測精度等指標(biāo)集合評(píng)估了LR、SWAT、ML和GAMLSS典型代表性模型的適用性。

(3)基于下墊面和水土保持措施等標(biāo)準(zhǔn)化數(shù)據(jù)及不同氣候模式輸入，篩選了SWAT、GAMLSS、ML、DWN、MFD-WESP、HydroTrend、產(chǎn)沙指數(shù)模型、人工智能模型、BP神經(jīng)網(wǎng)絡(luò)模型9類模型關(guān)于黃河潼關(guān)站未來50 a水沙變化趨勢預(yù)測結(jié)果，集合評(píng)估并預(yù)測了潼關(guān)站未來50 a年均徑流量約為240億m3、年均輸沙量為2.45億t，年徑流量和年輸沙量90%的置信區(qū)間范圍分別為[164，328]億m3和[0.79，5.12]億t。近期許多學(xué)者關(guān)于黃河水沙變化趨勢預(yù)測成果與本研究基于多模型預(yù)測集合評(píng)估成果一致。