基于減沙效益和經(jīng)濟(jì)效益的流域退耕還林方案優(yōu)化

喻 丹，董曉華※，彭 濤，劉 冀，劉 超，劉 芳

（1.三峽大學(xué)水利與環(huán)境學(xué)院，宜昌 443000；2.三峽庫區(qū)生態(tài)環(huán)境教育部工程研究中心，宜昌 443000；3.宜昌市夷陵區(qū)水資源利用服務(wù)中心，宜昌 443000；4.宜昌市水文局，宜昌 443000）

0 引言

土壤是人類賴以生存和發(fā)展的基礎(chǔ)自然資源，由于特殊的地理和歷史條件，中國水土流失量大且面積廣，是突出的生態(tài)環(huán)境問題之一。在眾多導(dǎo)致水土流失的因素之中，耕地又被認(rèn)為是水土流失的主要策源地[1]。為了遏制水土流失、改善生態(tài)環(huán)境，中國于1999 年啟動(dòng)實(shí)施了退耕還林工程，在長江、黃河和淮河等流域的長期實(shí)踐看，其控制水土流失的效果十分顯著。2014 年國家林業(yè)局啟動(dòng)了新一輪的退耕還林工程，旨在進(jìn)一步推進(jìn)中國的生態(tài)文明建設(shè)。

退耕還林是一項(xiàng)涉及社會(huì)、經(jīng)濟(jì)和生態(tài)多個(gè)方面的復(fù)雜系統(tǒng)性工程，工程的實(shí)施既要重視生態(tài)效益，又要兼顧社會(huì)經(jīng)濟(jì)效益，因此客觀評(píng)價(jià)退耕還林工程的成本有效性已成為衡量其是否成功的關(guān)鍵指標(biāo)。目前已有學(xué)者圍繞退耕還林實(shí)施后的生態(tài)效益評(píng)估[2-5]、經(jīng)濟(jì)效益評(píng)估[6]以及綜合效益評(píng)估[7]進(jìn)行了大量的研究。然而由于受到氣候、下墊面特征以及區(qū)域經(jīng)濟(jì)產(chǎn)業(yè)結(jié)構(gòu)的影響，退耕還林工程的各項(xiàng)效益存在著一定的空間差異[8-10]。

僅根據(jù)退耕還林工程實(shí)施后所進(jìn)行的“事后評(píng)價(jià)”難以直接給與管理決策者用于規(guī)劃退耕還林工程空間布置所需要的“事先評(píng)價(jià)”的信息。如何預(yù)先制定最具成本-效益的退耕還林方案，以平衡生態(tài)和經(jīng)濟(jì)之間的矛盾，已成為保證退耕還林工程可持續(xù)發(fā)展的關(guān)鍵。多數(shù)研究采用基于情景分析（scenario-based）的方法對(duì)退耕還林工程進(jìn)行優(yōu)化配置，即依據(jù)一定的準(zhǔn)則條件形成幾個(gè)退耕還林空間配置方案，再根據(jù)評(píng)價(jià)結(jié)果進(jìn)行方案比選、擇優(yōu)。黃康等[11-13]均根據(jù)不同的坡度設(shè)置了多個(gè)退耕還林情景，李越強(qiáng)等[14-15]則通過設(shè)置極端的情景，即假設(shè)全流域的耕地全都進(jìn)行退耕還林。于航等[16]以陡坡和河岸帶的耕地為目標(biāo)區(qū)域，設(shè)置了4 種不同的退耕還林情景。顯然這種基于情景分析的方法對(duì)于情景的假設(shè)是有限的，無法全局覆蓋退耕還林的所有情景，從而丟失潛在的更優(yōu)的配置方案。因此，需要探索基于全局搜索尋優(yōu)的方法實(shí)現(xiàn)退耕還林空間配置方案的優(yōu)化。

淮河流域是中國中東部地區(qū)重要的生態(tài)安全屏障，水土流失一直都是流域生態(tài)保護(hù)所面對(duì)的主要問題。2021 年度全國水土流失動(dòng)態(tài)監(jiān)測成果顯示淮河流域的水土流失表現(xiàn)為水力侵蝕為主，且主要分布在沂沭泗河、沙潁河以及淮河干流的上游區(qū)域[17]。在淮河典型流域開展退耕還林的空間配置研究對(duì)淮河流域水土流失綜合治理以及生態(tài)環(huán)境保護(hù)具有重要的意義。本文選擇淮河上游的息縣流域?yàn)檠芯繀^(qū)，先運(yùn)用生態(tài)水文模型模擬和GIS 空間疊置方法分別構(gòu)造退耕還林的生態(tài)效益目標(biāo)和經(jīng)濟(jì)效益目標(biāo)，再以子流域?yàn)樗亮魇е卫韱卧獦?gòu)成多目標(biāo)優(yōu)化問題，采用具有全局搜索能力的智能算法求解子流域尺度的退耕還林優(yōu)化方案，為淮河流域水土保持、退耕還林工程的科學(xué)規(guī)劃提供借鑒和指導(dǎo)。

1 研究區(qū)概況與數(shù)據(jù)來源

1.1 研究區(qū)概況

淮河流域（東經(jīng)111°55’～121°25’，北緯30°55’～36°36’）位于中國東部、長江和黃河兩大流域之間?；春痈闪靼l(fā)源于河南省南陽市桐柏縣和湖北省隨州市隨縣交界處的桐柏山太白頂，流經(jīng)河南、湖北、安徽和江蘇，全長1 000 km，流域總面積為27 萬km2。本文選取淮河干流上游的息縣水文站以上流域作為研究區(qū)（下文簡稱息縣流域），研究區(qū)概況如圖1 所示。

息縣流域面積為10 219.67 km2，流域內(nèi)以山區(qū)和丘陵區(qū)為主，地勢由西南和北部的山區(qū)逐漸向東降低。根據(jù)中科院資源環(huán)境科學(xué)數(shù)據(jù)中心提供的中國多期土地利用遙感監(jiān)測數(shù)據(jù)集（CNLUCC）以及沈玲怡等[18]對(duì)息縣流域2000—2020 年的土地利用變化分析可知，2000 年之后息縣流域并未發(fā)生明顯的土地利用變化，且耕地始終是流域主要的土地利用方式。以淮河干流為分界線，干流以北的耕地種植著玉麥和玉米為主的旱地作物，干流以南的耕地多為水田，主要種植水稻。其次是林地，主要分布在流域河源區(qū)域的桐柏山一帶。林地垂直分帶明顯，可分為中山落葉闊葉林，低山常綠落葉闊葉混交林，丘陵常綠針葉林及經(jīng)濟(jì)林[19]。然而，人口的增長和社會(huì)經(jīng)濟(jì)的發(fā)展導(dǎo)致了息縣流域土地資源的不合理利用，流域內(nèi)許多地區(qū)的自然植被遭到了破壞，加之耕地降低了土壤的抗侵蝕能力，因此引發(fā)了水土流失這一生態(tài)環(huán)境問題[20]。息縣流域已被水利部劃分為國家級(jí)水土流失重點(diǎn)預(yù)防區(qū)，此外還被納入第四期和第五期的全國水土保持重點(diǎn)工程建設(shè)項(xiàng)目規(guī)劃[21]。

1.2 數(shù)據(jù)來源

本研究使用的數(shù)據(jù)包括：研究區(qū)地形數(shù)據(jù)DEM（digital elevation model），來源于美國太空總署和國防部國家測繪局聯(lián)合測量的SRTM-DEM 數(shù)據(jù)（http://srtm.csi.cgiar.org/），其空間分辨率為90 m；土地利用空間分布情況從中國科學(xué)院資源環(huán)境科學(xué)數(shù)據(jù)中心（https://www.resdc.cn/Default.aspx）提供的中國多時(shí)期土地利用遙感監(jiān)測數(shù)據(jù)集中提取，空間分辨率為30 m；土壤數(shù)據(jù)來源于聯(lián)合國糧農(nóng)組織和維也納國際應(yīng)用系統(tǒng)研究所所構(gòu)建的世界土壤數(shù)據(jù)庫Harmonized World Soil Database -HWSD（https://www.fao.org/soils-portal/soil-survey/soilmaps-and-databases/harmonized-world-soil-databasev12/en/），空間分辨率為1 km；本研究采用中國GDP（gross domestic product）空間分布公里網(wǎng)格數(shù)據(jù)來表征研究區(qū)的社會(huì)經(jīng)濟(jì)狀況，該數(shù)據(jù)來源于中國科學(xué)院資源環(huán)境科學(xué)數(shù)據(jù)中心（https://www.resdc.cn/Default.aspx），空間分辨率為1 km；研究流域內(nèi)共計(jì)42 個(gè)雨量站的日降雨數(shù)據(jù)以及流域出口息縣水文站的流量觀測數(shù)據(jù)均由淮河水利委員會(huì)提供，時(shí)間范圍為1995—2010 年；其他氣象參數(shù)通過美國國家氣象預(yù)報(bào)中心提供的全球氣候再分析系統(tǒng)CFSR 計(jì)算得到。

2 研究方法

2.1 退耕還林的減沙效益

在降雨的驅(qū)動(dòng)作用下，土壤侵蝕剝離的泥沙伴隨著流域水文循環(huán)過程而遷移。流域生態(tài)水文模型不僅能對(duì)流域的水沙過程進(jìn)行有效模擬，還能定量評(píng)價(jià)水土保持措施對(duì)流域水沙的影響。SWAT（soil and water assessment tool）模型被廣泛應(yīng)用于預(yù)測不同土壤類別、土地利用類型以及管理措施下流域的水量、泥沙及水質(zhì)變化，具有較強(qiáng)的適用性[22]。根據(jù)研究區(qū)DEM 數(shù)據(jù)，SWAT 模型將全流域共劃分為49 個(gè)子流域，如圖1 所示。本研究選擇的土地利用數(shù)據(jù)時(shí)期為2010 年，土地利用分類依據(jù)的是中國土地資源分類系統(tǒng)，為了與SWAT 模型所含的土地利用數(shù)據(jù)庫對(duì)應(yīng)，如圖2a 所示將其重分類為耕地、林地、灌木、草地和水域，其中耕地和林地為研究區(qū)的主要土地利用類型，分別占全流域的55.41%和22.08%?；贖WSD 土壤提取本研究區(qū)的土壤類型空間分布如圖2b 所示，土壤類型的劃分依據(jù)FAO -90 土壤分類系統(tǒng)，其中堆積土和不飽和始成土為研究區(qū)的兩種主要的土壤類型，所占比例分別為28.66%和21.95%。由于SWAT 模型自帶的土壤數(shù)據(jù)庫沒有覆蓋到本文的研究區(qū)，故根據(jù)魏懷斌等[23]提出的方法計(jì)算土壤參數(shù)。根據(jù)相關(guān)研究[24]推薦的HRU（hydrologic response unit）劃分閾值（坡度10%、土地利用0、土壤10%），最終在研究流域劃分了528 個(gè)HRU。

圖2 研究區(qū)土地利用空間分布和土壤類型空間分布Fig.2 Land use and soil types throughout the study area

采用本研究收集到的研究區(qū)雨量站點(diǎn)的降雨數(shù)據(jù)和天氣發(fā)生器的氣象參數(shù)驅(qū)動(dòng)SWAT 模型，以研究流域出口息縣水文站實(shí)測的徑流和泥沙數(shù)據(jù)作為參照，以1995年為預(yù)熱期，1996—2003 年為率定期，2004—2010 年為驗(yàn)證期，采用SWAT-CUP 的SUFI-2 算法分別進(jìn)行月尺度的徑流參數(shù)率定以及月尺度的泥沙參數(shù)率定，選取統(tǒng)計(jì)指標(biāo)Nash 系數(shù)（Ens）、決定系數(shù)（R2）和百分比偏差（PBIAS）對(duì)SWAT 模型率定期和驗(yàn)證期的模擬精度進(jìn)行評(píng)定。

流域管理措施產(chǎn)生的效果會(huì)隨著配置地點(diǎn)地形、水文氣候、土壤等特征的不同而不同[25]。本研究利用SWAT模型來量化空間不同位置進(jìn)行退耕還林起到的水土保持效果。將建立SWAT 模型采用的2010 年土地利用為基準(zhǔn)，以子流域?yàn)橥烁€林的基本單元，每次僅選擇一個(gè)子流域（如子流域i）進(jìn)行退耕還林操作，再利用率定好的SWAT 模型進(jìn)行模擬，得到該子流域退耕還林后流域的產(chǎn)沙值Si，依次遍歷所有子流域，根據(jù)式（1）計(jì)算子流域退耕還林的泥沙削減系數(shù)Sci。

式中S0為基準(zhǔn)期土地利用條件下SWAT 模型模擬的流域產(chǎn)沙值，t；Si為子流域退耕還林后流域的產(chǎn)沙值，t；A耕i為 基準(zhǔn)期土地利用條件下第i個(gè)子流域中退耕還林的面積，即是基準(zhǔn)期該子流域耕地的面積，km2。

2.2 退耕還林的經(jīng)濟(jì)效益

土地是區(qū)域經(jīng)濟(jì)發(fā)展的最基本物質(zhì)載體，土地利用和經(jīng)濟(jì)發(fā)展有著密不可分的關(guān)系，社會(huì)經(jīng)濟(jì)結(jié)構(gòu)的變化直接決定了土地利用方式的變化?；春恿饔蜃鳛橹袊匾募Z食生產(chǎn)基地，其農(nóng)業(yè)產(chǎn)值和人均農(nóng)業(yè)產(chǎn)值均高于全國同期水平，流域內(nèi)的耕地資源有著極其重大的經(jīng)濟(jì)價(jià)值。因此可認(rèn)為在息縣流域?qū)嵤┩烁€林會(huì)給區(qū)域經(jīng)濟(jì)造成一定的影響。本研究采用資源環(huán)境科學(xué)與數(shù)據(jù)中心提供的中國GDP 公里網(wǎng)格分布數(shù)據(jù)來表征研究流域的經(jīng)濟(jì)發(fā)展水平，該數(shù)據(jù)是在綜合分析了人類活動(dòng)形成的土地利用格局與GDP 的空間互動(dòng)規(guī)律的基礎(chǔ)上，建立一二三產(chǎn)業(yè)GDP 數(shù)量與土地利用類型的空間相關(guān)性模型來實(shí)現(xiàn)GDP 的空間化，且認(rèn)為耕地、林地、草地、水體這些土地利用類型只貢獻(xiàn)第一產(chǎn)業(yè)GDP，第二三產(chǎn)業(yè)GDP 值僅與城鄉(xiāng)工礦居民地有關(guān)[26]。因此本研究針對(duì)的是退耕還林后第一產(chǎn)業(yè)的GDP 損失。

將2010 年的中國GDP 公里網(wǎng)格分布數(shù)據(jù)與基準(zhǔn)期2010 年的土地利用數(shù)據(jù)進(jìn)行空間疊置分析，分別得到基準(zhǔn)期土地利用條件下耕地和林地覆蓋的GDP 產(chǎn)值空間分布，再結(jié)合子流域空間分布可統(tǒng)計(jì)第i個(gè)子流域內(nèi)耕地和林地分別對(duì)應(yīng)的第一產(chǎn)業(yè)GDP，分別記為G耕i和G林i，并假設(shè)在同一個(gè)子流域內(nèi)，單位耕地或者單位林地產(chǎn)生的GDP 是均勻的。那么在基準(zhǔn)期土地利用條件下，第i個(gè)子流域的第一產(chǎn)業(yè)GDP 值G0i可表示為

選擇第i個(gè)子流域進(jìn)行退耕還林，退耕面積為A耕i，同時(shí)保持其他子流域的土地利用不變。退耕還林后該子流域原有耕地產(chǎn)生的GDP 產(chǎn)值就全部損失了，林地面積在原有的基礎(chǔ)上增加了。因此第i個(gè)子流域的第一產(chǎn)業(yè)GDP 值Gi可表示為

式中A林i為基準(zhǔn)期土地利用條件下第i個(gè)子流域內(nèi)林地的面積，km2。

由式（2）和式（3）即可推求退耕還林過程中第i個(gè)子流域第一產(chǎn)業(yè)GDP 產(chǎn)值的變化，記為△Gi：

為了消除退耕還林面積的影響，定義GDP 損失系數(shù)Gci為單位退耕還林面積的GDP 變化，具體表達(dá)式為

重復(fù)上述操作依次遍歷所有的子流域，即可得到空間不同子流域的GDP 損失系數(shù)。Gci＞0 表示GDP 在基準(zhǔn)期的基礎(chǔ)上增大，Gci＜0 表示GDP 相較于基準(zhǔn)期減小。

2.3 基于多目標(biāo)遺傳算法的退耕還林方案優(yōu)化

以小流域?yàn)閱卧乃亮魇ЬC合治理是中國水土保持工作的重要途徑，因此設(shè)定第i個(gè)子流域參與退耕還林面積xi為自變量。結(jié)合前述的泥沙削減系數(shù)Sci和GDP 損失系數(shù)Gci構(gòu)造研究區(qū)退耕還林的多目標(biāo)優(yōu)化問題，一方面生態(tài)效益目標(biāo)F1可由退耕還林后流域的泥沙產(chǎn)量來表示，如式（6）所示；另一方面用退耕還林導(dǎo)致的第一產(chǎn)業(yè)GDP 損失值來表征經(jīng)濟(jì)效益目標(biāo)F2，如式（7）所示。約束條件表示為每個(gè)子流域的退耕面積約束，如式（8）所示。

式中n為退耕還林的子流域個(gè)數(shù)。

本研究采用多目標(biāo)遺傳算法NSGA-II[27]對(duì)式（6）～（8）所描述的優(yōu)化問題進(jìn)行求解。NSGA-II 是基于Pareto 最優(yōu)概念的遺傳算法，用以解決多個(gè)沖突目標(biāo)之間的最優(yōu)化問題，其Pareto 前沿是由多個(gè)解構(gòu)成的最優(yōu)解集，且各Pareto 最優(yōu)解之間沒有優(yōu)劣之分。針對(duì)本研究的優(yōu)化問題，一個(gè)空間退耕還林方案xi(i=1,2,…,n)就是一個(gè)個(gè)體，NSGA-II 算法首先在約束條件下以隨機(jī)的方式生成初始化種群，即是一定數(shù)量的空間退耕還林方案；再根據(jù)生態(tài)效益目標(biāo)F1和經(jīng)濟(jì)效益目標(biāo)F2對(duì)每一個(gè)方案進(jìn)行多目標(biāo)評(píng)價(jià)；運(yùn)用算法的非支配排序、交叉、變異等操作產(chǎn)生新的子代種群；經(jīng)過迭代進(jìn)化最后得到最終的退耕還林Pareto 方案集。

優(yōu)化出的Pareto 前沿上存在一個(gè)向理想點(diǎn)方向的“凸起”，為Knee 解。Knee 解是Pareto 前沿上具有極大邊際效用的點(diǎn)，若沒有制定特定的偏好，則Knee 解被認(rèn)為是最具成本-效益的點(diǎn)，因此是最有可能被決策者優(yōu)先選擇的點(diǎn)。采用基于距離的方法確定Knee 解：根據(jù)Pareto前沿的極端點(diǎn)形成一條直線，計(jì)算Pareto 前沿上所有點(diǎn)到該直線的距離，距離最大的點(diǎn)解即為所求[28]。

3 結(jié)果與分析

3.1 SWAT 模型模擬結(jié)果

根據(jù)相關(guān)研究[29]確定如表1 所示的對(duì)徑流和泥沙模擬影響較大的敏感參數(shù)參與率定。SWAT 模型參數(shù)率定結(jié)果如表1 所示，徑流和泥沙模擬精度評(píng)定結(jié)果如表2所示。可以看出徑流模擬率定期和驗(yàn)證期的Ens值均在0.90 以上，決定系數(shù)R2的值接近于1，且總量誤差PBIAS控制在-10%以內(nèi)，說明月尺度徑流模擬的精度較高。相比之下，月尺度的泥沙模擬結(jié)果要略差于徑流。泥沙模擬率定期和驗(yàn)證期的Ens值均能夠保持在0.70 以上，R2值大于0.80，PBIAS的值在10%至20%之間，表現(xiàn)出一定的總量誤差，但依據(jù)精度指標(biāo)判定標(biāo)準(zhǔn)來看仍處于可接受的范圍。通過SWAT 模型模擬的月徑流和月泥沙與實(shí)測過程的對(duì)比線（圖3）可以看出，SWAT 模型能夠較好地重現(xiàn)研究流域徑流和泥沙隨時(shí)間的變化過程。綜合分析認(rèn)為本研究建立的SWAT 模型對(duì)研究區(qū)的月徑流和泥沙模擬具有較好的適用性。

表1 SWAT 模型徑流參數(shù)和泥沙參數(shù)率定結(jié)果Table 1 Optimal values for runoff and sediment parameters in SWAT model

表2 SWAT 模型精度評(píng)定結(jié)果Table 2 Accuracy evaluation results of the SWAT model

圖3 SWAT 模型模擬的月徑流和月泥沙與實(shí)測過程的對(duì)比Fig.3 Comparison of monthly runoff and sediment simulated by SWAT model with measured processes

3.2 子流域泥沙削減系數(shù)

根據(jù)現(xiàn)狀土地利用條件下的SWAT 模型模擬結(jié)果，統(tǒng)計(jì)得到基準(zhǔn)期流域多年平均泥沙產(chǎn)量為210.7 萬t，與鄒宏榮等[30]對(duì)息縣流域產(chǎn)沙模擬的結(jié)果較為一致。子流域?qū)嵤┩烁€林操作后，各子流域泥沙削減系數(shù)的空間分布如圖4 所示，空白部分表示該子流域沒有耕地，故不參與退耕還林。根據(jù)式(1)可知，計(jì)算的Sci數(shù)值越大，說明在該子流域?qū)嵤┩烁€林達(dá)到的減沙效果越好。由圖4 可以看出，Sci在子流域空間上呈現(xiàn)出一定的空間過渡趨勢，表現(xiàn)為上游河源區(qū)的值較大，中游到下游區(qū)域的值逐漸變小。說明對(duì)于單位面積的退耕還林操作，上游的河源區(qū)域能夠取得更加明顯的減沙效果。流域水系的河源區(qū)一般山高坡陡且降雨量較大，已有研究證明坡度和降雨對(duì)土壤侵蝕有著很大的促進(jìn)作用[31-33]，因此該區(qū)域?qū)儆谒亮魇舾袇^(qū)，不適合于發(fā)展耕地。Sci數(shù)值最小的3 個(gè)子流域依次為子流域18（26.70 t/km2）、子流域19（28.49 t/km2）和子流域12（29.22 t/km2），從空間分布上看這3 個(gè)子流域集中在流域下游干流以北的區(qū)域。

圖4 泥沙削減系數(shù)Sci 在子流域空間的分布Fig.4 Spatial distribution of sediment reduction coefficient Sci at sub-basin scale

3.3 子流域GDP 損失系數(shù)

結(jié)合現(xiàn)狀土地利用和GDP 網(wǎng)格數(shù)據(jù)疊置分析得到基準(zhǔn)期條件下研究流域總GDP 值為515.61 萬元，且各子流域內(nèi)林地和耕地覆蓋的GDP 值如表3 所示，對(duì)比可知對(duì)大部分子流域來說，耕地比林地產(chǎn)生了更高的GDP，這也與淮河流域農(nóng)業(yè)經(jīng)濟(jì)占主導(dǎo)地位這一事實(shí)相符。僅有部分子流域，如34、41 和44，其林地產(chǎn)生的GDP 略大于耕地。此外，子流域GDP 損失系數(shù)計(jì)算結(jié)果如表3所示，可以看出實(shí)施退耕還林后子流域GDP 的變化有著較大的空間差異，不僅出現(xiàn)了GDP 的降低，同時(shí)還出現(xiàn)了GDP 的增加。根據(jù)式(5)可知，從數(shù)值上看計(jì)算的Gci越大，說明在該子流域?qū)嵤┩烁€林操作后GDP 的損失越小，從經(jīng)濟(jì)角度上看就越適合進(jìn)行退耕還林。表3中Gci值大于0 的子流域有34 和41，說明這些區(qū)域原有的耕地產(chǎn)生的經(jīng)濟(jì)效益較低，在這些子流域進(jìn)行退耕還林能夠在不降低區(qū)域GDP 水平的前提下降低流域出口的泥沙含量，那么這些子流域內(nèi)的耕地應(yīng)盡可能多地參與退耕還林。退耕還林導(dǎo)致的GDP 損失最嚴(yán)重的為子流域30，其Gci值為-5 756.83 元/km2，表明該子流域內(nèi)每1 km2的退耕還林會(huì)使GDP 產(chǎn)值降低5 756.83 元。從地理位置上來看子流域30 位于信陽市市轄區(qū)附近，說明該區(qū)域內(nèi)的耕地為當(dāng)?shù)鼐用竦闹饕?jīng)濟(jì)來源，具有極其重要的經(jīng)濟(jì)價(jià)值，盲目地實(shí)施退耕還林會(huì)帶來較大的經(jīng)濟(jì)損失。除此之外，子流域4、11、15 和24 的Gci值均小于-1 000 元/km2，從地理位置上看這些子流域內(nèi)聚集著桐柏縣（子流域11）、信陽市（子流域4 和24）以及息縣（子流域15）的主要居民點(diǎn)，因此在這些區(qū)域進(jìn)行退耕還林時(shí)需要權(quán)衡經(jīng)濟(jì)損失和泥沙控制之間的平衡關(guān)系。

3.4 子流域退耕還林方案優(yōu)化結(jié)果

遺傳算法參數(shù)的選擇對(duì)優(yōu)化的最終結(jié)果有著一定的影響，為了確定合適的算法參數(shù)，對(duì)種群大小、迭代次數(shù)、交叉概率和變異概率這4 個(gè)參數(shù)分別進(jìn)行敏感性分析。如表4 所示設(shè)定不同的參數(shù)變化情景，采用OAT（one-at-a-time）方法，即每次更改一個(gè)參數(shù)的值同時(shí)將其他所有參數(shù)值固定在初始值。各參數(shù)情景的優(yōu)化結(jié)果如圖5 所示。由圖5a 可以看出當(dāng)?shù)螖?shù)從100 增加到600 時(shí)，優(yōu)化出的Pareto 前沿逐漸向著兩個(gè)目標(biāo)函數(shù)的最小值逼近，但是當(dāng)?shù)螖?shù)增加到1 000 時(shí)，Pareto前沿并沒有進(jìn)一步被優(yōu)化。圖5b 中顯示的種群大小大于以100 為分界線，種群大小超過100 得到的 Pareto 前沿明顯優(yōu)于種群大小為20 和50，當(dāng)種群大小從200 增加到500 時(shí)Pareto 前沿在兩個(gè)目標(biāo)函數(shù)的尋優(yōu)表現(xiàn)上并沒有顯著的提高。當(dāng)交叉概率從0.1 增加到0.7，圖5c 所示的Pareto 前沿也越來越優(yōu)，但是當(dāng)交叉概率進(jìn)一步增加到0.9 時(shí)，Pareto 前沿出現(xiàn)了反向移動(dòng)而變得更差。圖5 d 中所示的變異概率設(shè)置為0.000 5 得到的Pareto 前沿與初始值0.000 1 較為接近，且都劣于其余設(shè)定值得到的Pareto 前沿。綜合判斷最終確定遺傳算法的參數(shù)設(shè)定為：迭代次數(shù)600，種群大小400，交叉概率0.7，變異概率0.01。

圖5 遺傳算法參數(shù)敏感性分析結(jié)果Fig.5 Sensitivity analysis results of the genetic algorithm parameters

根據(jù)設(shè)定的遺傳算法參數(shù)進(jìn)行多目標(biāo)優(yōu)化，得到的Pareto 前沿圖6 所示，構(gòu)成了子流域尺度退耕還林多目標(biāo)優(yōu)化問題的最優(yōu)方案集，該方案集一方面將GDP 損失控制在105 萬元至196 萬元，損失率為30.13%至37.67%；另一方面將流域的產(chǎn)沙量控制在63.5 萬t 至97.9 萬t，相較于基準(zhǔn)期的泥沙產(chǎn)量分別削減了69.86%和53.54%。PIMENTEL[34]將土壤的形成速率0.5～1.0 t/hm2定義為區(qū)域可持續(xù)發(fā)展的土壤侵蝕水平，認(rèn)為土壤侵蝕的速度若大于這一臨界水平時(shí)，土壤有機(jī)體的多樣性和豐富度會(huì)大幅度降低從而不利于該區(qū)域農(nóng)業(yè)的可持續(xù)發(fā)展。據(jù)此標(biāo)準(zhǔn)換算出研究區(qū)息縣流域可持續(xù)發(fā)展水平下的土壤侵蝕量為51.2～102.5 萬t，不難看出對(duì)于優(yōu)化出的Pareto前沿上的所有解均沒有該水平的超過上限值，說明即使是在經(jīng)濟(jì)效益目標(biāo)的牽制下，優(yōu)化得到的退耕還林方案集的生態(tài)效益目標(biāo)仍然能夠達(dá)到較為滿意的程度。

圖6 多目標(biāo)優(yōu)化結(jié)果的Pareto 前沿及Knee 解Fig.6 Pareto frontier and Knee solution of multi-objective optimization results

退耕還林直接減少了糧食生產(chǎn)的基本物質(zhì)基礎(chǔ)，對(duì)區(qū)域糧食安全產(chǎn)生的影響不容忽視。人均耕地面積這一指標(biāo)常被用來分析地區(qū)或區(qū)域的糧食安全問題。經(jīng)統(tǒng)計(jì)得到Pareto 前沿解集中全流域的退耕面積為2 508.5～3 889.0 km2，根據(jù)研究流域所轄行政區(qū)統(tǒng)計(jì)年鑒統(tǒng)計(jì)研究流域的人口數(shù)量（2 753 960 人），由此計(jì)算退耕還林后剩余的人均耕地面積為1.04×10-3～1.54×10-3km2，均高于陳百明等[35]預(yù)測的近期人均耕地面積閾值0.74×10-3km2和遠(yuǎn)期人均耕地閾值0.64×10-3km2，說明本研究所求的退耕還林方案解集即使是減少了研究區(qū)耕地的面積也仍然能夠滿足人均糧食供給。另外，基于本研究構(gòu)造的優(yōu)化問題，即流域產(chǎn)沙量最小和GDP 損失最小，在優(yōu)化算法求解的過程中會(huì)引導(dǎo)著退耕還林面積更多地偏向于土壤易侵蝕且耕地產(chǎn)值較低的區(qū)域配置。時(shí)亞坤等[36-37]的研究均強(qiáng)調(diào)，一方面退耕還林實(shí)施區(qū)的耕地質(zhì)量差且糧食產(chǎn)量本就較低，另一方面耕地面積的減少促使農(nóng)戶從粗放式的生產(chǎn)模式轉(zhuǎn)變?yōu)榧s式，進(jìn)而有效改善了區(qū)域環(huán)境和土壤肥力，因此可認(rèn)為退耕還林的實(shí)施有助于提高糧食產(chǎn)量。綜合分析來看，本研究求解的退耕還林方案集能夠發(fā)揮水土保持和糧食供給的協(xié)同作用。

針對(duì)Pareto 前沿上的Knee 解（如圖6 所示），在土地利用現(xiàn)狀上確定全流域的退耕還林方案，如表5 所示按照子流域退耕還林面積百分比從大到小排列?？梢钥闯錾嫌巫恿饔蛲烁€林面積比例一般較大，基本上都超過了50%。退耕還林面積占比最小的幾個(gè)子流域（子流域19、30、4 和24）主要集中在信陽市轄區(qū)及周邊。根據(jù)表5 統(tǒng)計(jì)全流域退耕還林總面積為2 988.87 km2，占現(xiàn)狀土地利用耕地面積的44.19%。該退耕還林方案下對(duì)應(yīng)的兩個(gè)目標(biāo)函數(shù)分別為泥沙產(chǎn)量75.1 萬t 和GDP 損失130 萬元，經(jīng)計(jì)算可知該方案以犧牲約25%的經(jīng)濟(jì)損失換取了約65%的泥沙削減量。盛莉等[38]從土壤保持價(jià)值和涵養(yǎng)水源價(jià)值兩個(gè)方面估算了中國水土保持的生態(tài)功能服務(wù)價(jià)值，在全國平均水平得出單位重量（t）的土壤保持量可提供58.96 元的生態(tài)服務(wù)價(jià)值。據(jù)此可以推斷本研究中削減的泥沙量能夠提供約7 994 萬元的生態(tài)服務(wù)價(jià)值，遠(yuǎn)大于估算的GDP 損失值，因此可認(rèn)為優(yōu)化出的退耕還林方案對(duì)保護(hù)研究區(qū)生態(tài)環(huán)境有著巨大的潛力，具有較高的成本效益。

表5 Knee 解對(duì)應(yīng)的子流域尺度的退耕還林面積Table 5 GGP area at sub-basin scale corresponding to the Knee solution

考慮到退耕還林工程需要面對(duì)退耕耕地種類和還林林地種類的選擇，本研究遵循研究區(qū)自然規(guī)律，堅(jiān)持因地制宜的原則：子流域內(nèi)面積不占絕對(duì)主導(dǎo)的耕地類型優(yōu)先參與退耕；若子流域內(nèi)水田旱地面積相當(dāng)，則根據(jù)淮河干流南北的氣候差異，干流以北的耕地優(yōu)先考慮退水田，干流以南的區(qū)域優(yōu)先考慮退旱地；根據(jù)退耕區(qū)高程差異，高程800 m 以上種植落葉闊葉林，高程500～800 m 種植常綠闊葉林及落葉闊葉林，高程500 m 以下種植常綠針葉林及經(jīng)濟(jì)林。據(jù)此原則各子流域具體的退耕還林方式可分為五類：河源子流域11 內(nèi)的水田零散分布且海拔較高，因此該子流域內(nèi)優(yōu)先退水田還落葉闊葉林；相鄰的河源子流域26 內(nèi)面積較少的水田優(yōu)先退耕，退常綠闊葉林及落葉闊葉林；子流域10 和子流域36 分別位于桐柏縣和信陽浉河區(qū)，面積不占主導(dǎo)的旱地優(yōu)先退耕，還常綠闊葉林及落葉闊葉林；淮河干流以北的大部分子流域（1～9、12、13、16、18、19、21 和24）以及干流以南的部分子流域（25、27、35 和38）退水田還針葉林及經(jīng)濟(jì)林；位于淮河干流以南的剩余子流域則退旱地還常綠針葉林及經(jīng)濟(jì)林。

4 結(jié)論

本研究結(jié)合分布式水文模型模擬、空間疊置分析以及多目標(biāo)優(yōu)化方法，以流域泥沙產(chǎn)量和經(jīng)濟(jì)損失最小構(gòu)造多目標(biāo)優(yōu)化問題，求解得到子流域退耕還林最優(yōu)解集，并從中選取了最具成本效益的退耕還林方案。本文的主要研究結(jié)論有：

1) SWAT 模型徑流模擬率定期（1996—2003 年）和驗(yàn)證期（2004—2010 年）的Nash 系數(shù)Ens均大于0.90，決定系數(shù)R2均大于0.90，百分比偏差PBIAS分別為-4.26%和-5.73%；對(duì)于泥沙模擬率定期和驗(yàn)證期的Nash 系數(shù)Ens均大于0.70，決定系數(shù)R2均大于0.80，百分比偏差PBIAS均控制在20%之內(nèi)。由于本文在數(shù)據(jù)收集過程中受到一定的限制，只針對(duì)歷史時(shí)期的徑流和泥沙過程進(jìn)行了模擬。但是考慮到研究區(qū)的氣象和下墊面條件沒有發(fā)生較大的突變，因此認(rèn)為基于歷史時(shí)期模型率定和驗(yàn)證得到的參數(shù)仍能夠沿用到近期。因此，精度評(píng)定結(jié)果能夠說明SWAT 模型對(duì)研究流域徑流和泥沙的模擬具有較好的適用性。

2) 泥沙削減系數(shù)在空間上呈現(xiàn)出上游向下游逐漸減小的過渡趨勢，說明在上游進(jìn)行退耕還林操作能取得更好的泥沙控制效果。

3) GDP 損失系數(shù)其在空間上呈現(xiàn)出較大的差異性，在居民點(diǎn)聚集的區(qū)域進(jìn)行退耕還林操作會(huì)帶來較大的經(jīng)濟(jì)損失。

4) 子流域尺度退耕還林最優(yōu)方案集能夠在經(jīng)濟(jì)損失了30.13%～37.67%的同時(shí)將流域泥沙產(chǎn)量削減53.54%～69.86%，且能夠平衡水土保持、經(jīng)濟(jì)損失和糧食供給三者之間的矛盾。同時(shí)，最具成本效益的退耕還林方案對(duì)應(yīng)的退耕面積占總耕地面積的44.19%且主要集中于上游子流域，由此帶來的泥沙削減率約為65%，經(jīng)濟(jì)損失約為25%。