基于Pareto解集的治沙區非線性分式多目標水土資源優化配置模型

郭文賢 郭 萍 潘 琦 羅 彪

(中國農業大學 水利與土木工程學院，北京 100083)

土地沙漠化是干旱、半干旱及部分半濕潤地區由于人地關系不協調所造成的土地退化，以風沙活動為主要標志[1]。全球土地沙漠化面積為3 600×104km2，對全球24%的陸地面積和14%的人口存在潛在影響[2]。我國是世界上土地沙漠化危害最嚴重的國家之一[3]。截至2014年，我國土地沙漠化總面積達172.12×104km2，占國土總面積的17.93%[4]。沙漠化地區往往具有水資源匱乏、土壤保水能力差、地表植被稀疏等特征，其脆弱的生態環境嚴重制約了區域經濟發展和人民生活質量。沙漠化治理是指由水(文)、土(壤)、氣(候)和生(物)4個主要要素組成的生態系統退化之后在人工干預下的整體恢復過程[5]。常見的沙漠化治理措施主要有機械固沙、化學固沙和生物固沙，其中，機械固沙和化學固沙分別通過設置障礙物和改變沙地土壤孔隙結構來穩固流沙，生物固沙則通過土壤-植物-大氣連續體改善沙漠化地區生態環境，從而保證沙漠化治理的可持續性，在實際工程中普遍采用。國內外對于生物固沙的研究多從固沙植物特性[6-9]、固沙林效益[10-11]、固沙林體系建設[12]、固沙過程中水資源的利用[13-14]等方面展開探索，取得了一定成果。但現有的研究大多只聚焦于治沙過程中的單一方面，較少從生態系統整體上對防沙治沙進行綜合分析和治理。因此，如何充分考慮治沙過程中水、土地、生物等要素之間的關系和作用，從生態系統可持續發展的角度進行沙漠化防治，是當前亟待解決的問題。

依據決策主體的意愿，農業水土資源優化配置模型可分為提高灌溉效率的模型[15-16]、水土資源決策系統開發模型[17]和調整農業結構的模型[18]等三類。Ren等[19]構建的考慮地下水和地表水聯合運用的模糊分式規劃模型在優化灌水的同時提高了用水效率；Guo等[20]構建了考慮節水和質量保證的不確定規劃模型來提高河西走廊制種玉米的灌溉效率；Zhang等[21]開發了集決策支持、信息管理和檢索于一體的灌區水資源優化配置與灌溉決策支持系統；張帆等[22]構建了分式兩階段隨機規劃模型來對甘肅省民勤縣作物進行種植結構優化。模型的求解方法也從最小偏差法[23]、模糊協調法[24]、理想點法[25]等傳統方法逐漸擴展到智能算法[26-27]。對于復雜的多目標問題，智能算法具有更加強大的全局尋優能力。隨著對智能算法的不斷改進，蟻群算法[28]、粒子群算法[29]、遺傳算法[30]等智能算法在求解多目標優化問題方面得到廣泛應用，Guo等[31]將遺傳算法和逼近理想解排序法相結合求解了非線性多層次多目標優化模型；Wang等[32]采用并列選擇遺傳算法求解了分布式AquaCrop—非線性多目標模糊相關機會規劃模型。沙漠化治理問題與農業水土資源優化配置問題有著極高的相似性，都包含了種植結構、系統效益、用水效率和植物耗水量等指標的優化問題，然而當前水土資源優化配置在治沙過程中的研究還存在一定的空白。

基于以上問題，本研究構建了非線性分式多目標水土資源優化配置模型，模型考慮了治沙過程中水資源、土地資源和生物資源之間的作用與制約關系，將治沙過程中水、土地、生物要素聯合配置，以期在提高沙漠化治理效益的同時兼顧水資源利用效率，采用多目標遺傳算法求解模型，以討論不同目標之間的博弈過程，最終得到一系列優化配置方案，為沙漠化治理的管理者和決策者提供豐富有效的理論依據和科學指導。

1 模型建立與求解

1.1 非線性分式規劃原理

管理者進行決策時，為衡量系統的效率，通常會進行比率優化，當構建的模型中目標為物理或經濟數量的分數，我們將這種模型形式稱為分式規劃問題，分式規劃可以很好的反映研究對象的效率問題，同時也體現了分子分母2個模型之間的均衡[22]，當分式規劃的子函數存在非線性關系時，該分式規劃為非線性分式規劃，典型的非線性分式規劃表達式如下：

(1)

式中：f(x)為目標函數；g(x)、h(x)為關于x的函數，g(x)和h(x)至少存在一個是關于x的非線性函數；A為m×n的矩陣；x和b分別是n維和m維的列向量。

1.2 非線性分式多目標水土資源優化配置模型

人工治沙林的生態效益體現在固碳制氧、改善小氣候、保育土壤、降低風速、保護生物多樣性、改變沙地水分狀況等方面。治沙初期，管理和決策者要考慮治沙林種植植物存活的問題，植物存活數量占種植總數量的比例即為植物存活率。謝菲等[33]通過對沙生植物不同梯度的水分處理，得到植物存活率與植物用水量呈二次函數關系的結論，本研究引用此二次函數模型來表征植物存活率。植物存活率和累積生物量結合體現為治沙區固碳制氧總價值，這一指標在治沙初期能夠很好地描述治沙林生態效益，故本研究選擇治沙初期固碳制氧總價值反映治沙林的生態效益。然而，水資源是制約沙漠化地區生態發展的主要因素，干旱區沙漠治理不僅要考慮生態效益的最大化，還要考慮有限水資源的利用效率，本研究以植物存活數量和總用水量的比值表征治沙區水資源利用效率。國內外關于水土資源優化配置的研究成果豐富，但大多將水資源和土地資源分開研究，水土資源聯合優化配置較少，治沙區植物配水量和種植面積相互影響，如梭梭種植比例的增大會導致配水總量的增大，配水量減小時管理者會傾向于種植耗水量更小的花棒，因此本研究將植物配水量和種植面積同時設定為決策變量。

為對治沙區水土資源進行優化配置，采用多目標非線性分式規劃方法構建水土資源優化模型，目標函數為最大化固碳制氧總價值和最大化用水效率，約束條件有可供水量約束、最小需水量約束、面積約束和非負約束。模型涉及到的相關參數和變量意義見表1。

本模型中的目標函數如下：

1)最大化固碳制氧總價值

(2)

Vi=1.63×Bi×0.272 9×10-6×250+ 1.2×Bi×10-6×369.7

(3)

2)最大化治沙區用水效率

(4)

約束條件：

1)可供水量約束

(5)

2)最小需水量約束

0.1mij+P0j≥Wij,min

(6)

3)面積約束

(7)

Ai,min≤Ai

(8)

4)非負約束

處女作《回到明朝當王爺》橫掃網絡，囊括多項年終大獎；《大爭之世》《步步生蓮》入圍臺灣圖書借閱榜前二十名；《錦衣夜行》榮獲臺灣金石堂銷量第一名。

mij≥0

(9)

Ai≥0

(10)

(11)

1.3 模型求解

本研究建立的非線性分式多目標優化模型的2個目標相互沖突，不存在能使2個目標同時達到最優的解，但存在這樣的一個解：進一步優化一個目標時，必然存在另一個目標的劣化，這樣的解稱為非劣最優解。對于一個多目標優化問題可能存在多個非劣最優解，即多目標優化問題的非劣解集。遺傳算法作為一種魯棒的全局優化算法，提供了一種求解非線性、多模型、多目標等復雜優化問題的通用框架，并且已應用于許多水量和水質管理問題[34]。模型采用帶精英策略的并列選擇法進行求解，采用實數編碼的遺傳算法，相較于二進制編碼的遺傳算法，實數編碼的遺傳算法具有收斂速度更快、精度更高等優勢。求解過程中的每一步算子結束后都進行驗證，確保個體在可行域內，加快了模型的收斂速度。模型求解思路如圖1所示，求解步驟的具體描述如下：

表1 本研究中相關參數意義說明Table 1 Description of related parameters in this studynni

圖1 遺傳算法流程Fig.1 Flow chart of genetic algorithm

1)產生初始種群。初始化運行參數，采用英國謝菲爾德大學遺傳算法工具箱在解空間內生成隨機分布的種群，篩選出滿足約束條件的個體，種群個體數為100，將其作為初始種群。

2)根據適應度函數選擇。本研究選擇并列選擇法[35]求解算法，構建的水土資源優化配置模型有2個目標，故將父代種群中的個體均等的劃分為2個種群，種群1以固碳制氧總價值為目標函數，種群2以用水效率為目標函數，2個目標函數在各自種群獨立的進行運算，以目標函數值作為適應度，根據得到的2個種群的適應度，以0.9為代溝選擇出適應度高的個體，合并2個種群并去除不滿足約束條件的個體，得到選擇后的種群。

3)交叉和變異。將選擇后的種群內的個體隨機搭配成對作為2個父代，對2個父代進行整體算術雜交，得到重組種群，將交叉重組得到的種群以0.05的概率進行變異，本研究采用實值變異算子，保證了種群的多樣性，將交叉變異之后得到的種群作為新一代的種群。

4)終止條件。設置最大迭代次數為100，當迭代次數超過最大迭代次數時，終止運行，記錄得到的最終種群。判斷最終種群是否使Pareto前沿面穩定，若前沿面穩定，結束程序，反之重新進入循環。

2 案例研究

2.1 研究區概況

治沙區S共有面積500 hm2，地處干旱半干旱過渡區，屬于溫帶荒漠干旱氣候區，干燥少雨，蒸發強烈，年蒸發量為1 800～2 800 mm，年蒸發量為年降雨量的15倍以上，有效降雨為0[36]，水資源十分匱乏，存在較為嚴重的荒漠化現象，風沙流動強烈，生態環境極其脆弱，僅僅依靠天然沙生植物遠不足以緩解生態危機，人工種植沙生植物是遏制沙化面積擴大、改善當地生態環境的重要途徑。區域內治沙林建設主要存在以下問題：1)水資源極度緊缺，區域內惡劣的自然環境導致水資源總量不足，在滿足生活生產用水后可用于治沙林建設的水資源進一步減少；2)需要人工灌水保證治沙植物的存活，區域內土壤為風沙土，養分特征為平氮、少磷、缺鉀、有機質含量低，土壤貧瘠，持水能力差，不利于植物生長，無人工灌水時治沙植物無法存活，無法達到理想的生態恢復效果。面對以上問題，該區域治沙管理者開始重點關注如何兼顧治沙林生態效益和水資源利用效率，因此，上節中提出的模型在治沙區水土資源優化配置研究中意義重大。研究表明花棒，梭梭，沙拐棗均具有根系發達、耐旱保水、生長迅速的特點，為治沙優勢植物[37]，故治沙區選擇花棒，梭梭，沙拐棗為治沙植物。

2.2 相關參數的確定

花棒、梭梭、沙拐棗根據相關試驗，獲得其灌水量與存活率之間的函數關系[33]。3種植物的主要生長季為4—10月，種植密度為2 000株/hm2，生長季內的每個月均需灌溉，生長季內可供水量為3 500 m3/hm2。3種植物的騰發量根據相關試驗獲得[38-39]，考慮到沙漠區植物的抗旱性，植物最小需水量由植物騰發量乘以折減系數0.5確定，模型所需的植物生物量數據、最小種植面積、存活率相關關系等參數已有相關研究[40-43]，具體數據見表2和表3。

3 結果與分析

3.1 水土資源優化結果

圖3為非劣集的植物配水量。可見：隨著目標傾向于提高固碳制氧總價值(方案編號減小)，花棒、梭梭、沙拐棗全生長季配水量逐漸增大，這是因為配水量直接影響植物的存活率，從而影響植物的固碳制氧價值。圖4為3種植物固碳制氧價值對配水量的響應曲線。可見：3種植物的單位面積固碳制氧價值都隨著配水量的增加而增大，但固碳制氧價值大小差異明顯，從大到小依次為梭梭、沙拐棗、花棒；曲線的斜率代表響應速率，3種植物的響應速率都隨著配水量的增加而減小,模型中3種植物灌水量與存活率之間都為二次函數關系，因此3種植物對配水量的響應存在相似的變化規律；植物單株生物量只影響植物固碳制氧價值的大小，對響應趨勢和響應速率沒有影響，3種植物的生物量關系為梭梭>沙拐棗>>花棒，說明生物量會使同一配水量下的植物固碳制氧價值存在差異。因此當決策者傾向于追求固碳制氧總價值時，建議增加系統配水量，同時可以通過調整植物種植比例來提高梭梭的配水比例，限制花棒配水比例，以保證治沙林生態效益的提升。

表2 最小種植面積、單株生物量及灌溉水量和存活率相關關系的參數Table 2 Minimum planting area, biomass and correlation coefficient between irrigation amount and survival rate

表3 4—10月植物最小需水量Table 3 Minimum water requirement of plants from April to October mm

圖2 用水效率和固碳制氧總價值博弈下初始種群、最終種群與Pareto解集的表現Fig.2 Initial population, final population and Pareto solution set under the trade-off between water use efficiency and total value of carbon sequestration and oxygen production

圖3 非劣方案集配水量Fig.3 Water distribution of noninferior scheme set

圖4 植物配水量與固碳制氧價值相關性Fig.4 Correlation between plant water distribution and income from carbon sequestration and oxygen

本研究將用水效率定義為植物總成活數量和總配水量的比值，即灌水量與存活率函數關系的導數。可見：隨著單位面積配水量的增加，3種植物的用水效率逐漸減小。花棒用水效率從0.84株/m3下降至0.54株/m3，梭梭用水效率從0.69株/m3下降至0.45株/m3，沙拐棗用水效率從0.60株/m3下降至0.42株/m3；同時，在配水量相同的情況下，花棒的用水效率顯著高于梭梭和沙拐棗。隨著配水量增加，函數導數逐漸減小，用水效率逐漸降低，模型若想獲得較大的用水效率，將傾向于減少植物單位面積配水量。因此，當決策者傾向于追求用水效率時，應減小配水量，同時提高花棒的配水比例，以保障對水資源的高度節約。

圖5 植物配水量與用水效率相關性Fig.5 Correlation between plant water distribution and water use efficiency

非劣集植物種植面積配置結果見圖6。可見：當模型傾向固碳制氧總價值最大時，花棒和沙拐棗配置面積接近最小配置面積，這是因為梭梭的單株生物量遠大于花棒和沙棗拐，生物量越大，固碳制氧總價值越大；當模型傾向用水效率最大時，花棒為模型面積配置中的優勢植物，由謝菲等[33]的研究可知，在相同灌水量處理下，花棒的存活率大于梭梭和沙拐棗，3種植物種植密度相同，相同面積下花棒存活數量最多，存活數量和配水量的比值最大，即花棒的用水效率最大；為保證用水效率最大化，模型更傾向于減小梭梭、沙拐棗種植面積，增加花棒種植面積；無論優先考慮固碳制氧總價值目標還是優先考慮用水效率目標，沙拐棗在面積配置上均不占優勢，這是因為沙棗拐的最小需水量遠高于花棒和梭梭，分別是花棒的1.40倍和梭梭的1.46倍，且與梭梭相比，沙拐棗單株生物量不占優勢。在對有限的土地資源進行配置時，模型更傾向于降低沙拐棗的配置面積，以保證花棒、梭梭的種植面積處于較高水平，從而獲得更大的固碳制氧總價值以及更高的用水效率。然而，治沙林植物越單一，抗病蟲害能力越弱，越容易退化。梭梭、花棒和沙拐棗從形態特征來看分別屬于小喬木、半灌木和灌木，種植多種形態的植物增加了治沙林物種豐富度和組成多樣性，同時為其他動植物提供了多元的生存空間，有利于增強生態系統的穩定性，降低治沙林退化風險，能更好的滿足治沙區不同階段的植被演替，廖空太等[44]通過對治沙林適應性和抗逆性的研究也表明梭梭、花棒和沙拐棗混交種植優于梭梭和花棒混交種植。因此建議決策者在對植物種植面積進行規劃配置時僅將沙拐棗作為豐富物種多樣性的植物，給予其允許的最小配置面積。

圖6 非劣方案集種植面積配置Fig.6 Planting area configuration of non-inferior scheme

3.2 Pareto解集和目標博弈過程

本研究將水資源和土地資源同時優化，進一步擴大了尋優范圍，單一優化水資源或者土地資源時方案選擇范圍大大縮小，得到的非劣解數量也低于水土資源同時優化時的非劣解數量，同時優化水土資源能夠保證治沙區“地盡其力，水盡其用”。治沙區固碳制氧總價值和用水效率2個目標均為越大越優，在Pareto前沿面上，非劣配置方案的一個目標值增大必然會引起另一個目標值的減小，從第1號方案至111號方案，模型優化結果由傾向固碳制氧總價值最大逐步過渡到傾向用水效率最大，灌溉用水效率從0.49株/m3增至0.80株/m3，而固碳制氧價值從4.98萬元下降至1.38萬元。對比發現，不存在一種方案使得優化的2個目標值同時達到最大，即2個目標之間存在較為激烈的博弈關系。為量化固碳制氧總價值目標和用水效率目標之間的博弈關系，本研究引入變化率CR。

(12)

式中：CR為用水效率和固碳制氧總價值博弈的變化率;f2，m+1和f2，m分別為為第m+1和第m個方案時用水效率目標值;f1，m+1和f1，m分別為為第m+1和第m個方案時固碳制氧總價值目標值。

圖7為變化率實際值和變化率的擬合曲線。可見：方案1～21中，變化率擬合值大于1，邊際效應大，即固碳制氧總價值的減小對于用水效率的提升貢獻很大，但此時用水效率小，不適合水資源匱乏的治沙區；方案104～111中，變化率擬合值小于0.1，邊際效應小，即固碳制氧總價值的減小對于用水效率的提升貢獻很小，若選擇該范圍內的方案，治沙區生態恢復速度無法保障；方案22～103在保障生態恢復速度的同時兼顧了用水效率，是較為理想的選擇范圍。

圖7 用水效率和固碳制氧總價值博弈的變化率及其擬合曲線Fig.7 Change rate and fitting curve of the trade-off between water use efficiency and total value of carbon sequestration and oxygen production

3.3 典型方案的對比

為具體分析不同目標傾向下植物的水土資源配置情況，選取第1、55、111號方案為代表進行分析，3組方案分別代表固碳制氧總價值目標權重最大時、2個目標權重處于平衡時和用水效率目標權重最大時的最優解，3個典型方案的固碳制氧總價值依次為4.93、4.68和1.38萬元，用水效率依次為0.49、0.63 和0.79株/m3，圖8為典型方案配水量和植物種植面積。當追求固碳制氧總價值最大時(1號方案)，花棒、梭梭和沙拐棗的配水量得到充分滿足，面積配置優先滿足單位面積固碳制氧價值更高的梭梭，花棒和沙拐棗種植面積均為最小種植面積；當追求用水效率最大時(111號方案)，可供水量不足，面積配置優先滿足用水效率更高的花棒，梭梭和沙拐棗種植面積接近各自最小種植面積，總配水量相較于55號方案減少了27.08%，總配水量和梭梭種植面積同時減小導致固碳制氧總價值減小了70.29%，用水效率增加了26.11%。對比以上3個典型方案，方案1雖然能得到最大的固碳制氧價值，但配水量遠高于另外2種方案，對水資源造成一定程度的浪費，難以適用于水資源匱乏、取水不便的治沙區；方案111用水效率最大、配水量最小，但固碳制氧總價值遠小于方案1和方案55的固碳制氧總價值，很有可能造成前期投入巨大但生態恢復成果不佳的情況，難以完成治沙區生態建設的目標；55號方案的固碳制氧總價值略小于1號方案，用水效率在3種方案中居中，即在少量減少固碳制氧總價值的情況下用水效率提升明顯；考慮到治沙區需要保證一定的治沙效益，又存在水資源極度匱乏的問題，建議決策者在55號方案附近選擇決策結果。

圖8 不同典型方案下的配水量和種植面積Fig.8 Water distribution and acreage under different typical schemes

3種典型方案下植物對固碳制氧總價值的貢獻和用水效率的貢獻見表4。可以看出，3種方案下植物對固碳制氧總價值的貢獻從大到小排序均為梭梭>>沙拐棗>花棒，這是因為植物單株生物量關系為梭梭>沙拐棗>>花棒，方案111雖然增加了花棒的種植面積，但梭梭最小種植面積產生的固碳制氧價值仍大于花棒最大種植面積產生的固碳制氧價值；3種方案下植物對用水效率的貢獻為花棒略大于梭梭和沙拐棗，灌水量和成活率的函數關系決定了植物用水效率，從2個目標的角度看，治沙優勢物種分別為梭梭和花棒，與3.1小節的結論一致。因此建議決策者優先保證梭梭的種植面積，以獲得更大的固碳制氧總價值，提升花棒的成活率，當配水量對花棒成活率提升不明顯時可以考慮對花棒增加施肥、病蟲害防治等農業措施，以保證治沙區用水效率。

3.4 敏感性分析

沙漠化地區水資源問題突出，水是沙漠化地區修復植被、恢復生態的關鍵因素。為測試模型對水資源的敏感性，選擇9種供水量情景對模型進行求解，求解得到的Pareto前沿面如圖9所示，可以看出，供水量對Pareto前沿面存在明顯的影響，供水量越小，Pareto前沿面越接近原點，這意味著隨著可供水量減小，固碳制氧總價值和用水效率都相應減小；當供水量超過3 500 m3/hm2時，Pareto前沿面不再向遠離原點方向移動，即超過3 500 m3/hm2的供水量存在水資源浪費的問題，因為本研究基于灌水量和存活率的二次函數關系，當灌水量增加到一定程度時，植物存活率不再增大，表現在Pareto前沿面上即為前沿面不再移動。

表4 不同典型方案下植物對固碳制氧總價值(f1)和用水效率(f2)的貢獻Table 4 Contribution of plants to the total value of carbon sequestration and oxygen production (f1) and water efficiency (f2) under different typical scheme

圖9 不同供水量情形下Pareto前沿面Fig.9 Pareto fronts under different water supplies

不同用水量情景下最大固碳制氧總價值存在較大差異，而最大用水效率相差較小，都在0.78～0.80株/m3之間，說明可供水量對固碳制氧總價值的影響大于對用水效率的影響。為量化供水量的變化對固碳制氧總價值的影響，選取上文55號典型方案(目標處于平衡傾向)的用水效率附近不同供水量的最優解，得到目標傾向處于平衡時固碳制氧總價值對不同供水量的響應關系，當供水量小于3 000 m3/hm2時，單位面積上(1 hm2)增加1 m3供水量，固碳制氧總價值增加20.76元，供水量的增大對固碳制氧總價值提升明顯，供水量超過3 000 m3/hm2時，單位面積上(1 hm2)增加1 m3供水量，固碳制氧總價值只增加0.84元，供水量的增大對固碳制氧總價值的提升影響微弱，即當目標傾向處于平衡時，最佳供水量為3 000 m3/hm2，最佳供水量能在獲得理想的固碳制氧總價值的同時充分利用水資源。

4 結 論

本研究將水土資源優化配置方法應用于防沙治沙，構建了治沙區非線性分式多目標水土資源優化配置模型，采用遺傳算法求解該模型，主要結論如下：

1)當決策偏好傾向于固碳制氧總價值目標時，植物配水量會向灌水上限靠近，當決策偏好從固碳制氧價值最大向用水效率最大轉變時，植物單位面積配水量呈減少趨勢。無論目標傾向如何變化，沙拐棗面積配置均不占優勢。

2)水資源和土地資源同時優化擴大了尋優范圍，可為決策者提供更加豐富的Pareto解集，Pareto解集可以體現目標之間豐富的博弈信息，通過量化2個目標之間的博弈信息可以篩選出能同時兼顧生態效益和用水效率的方案選擇范圍。

3)55號方案同時兼顧了治沙區的生態恢復和水資源利用效率，是更合理的選擇。從固碳制氧總價值角度看，梭梭為治沙優勢物種，從用水效率角度分析，花棒為優勢物種。

4)超過3 500 m3/hm2的供水量存在水資源浪費現象；當目標處于平衡傾向時，治沙區最佳供水量為3 000 m3/hm2。

綜上所述，本研究的優點在于將水土資源優化配置與防沙治沙相結合，可以為治沙區水資源和土地資源配置提供決策依據。但在實際治沙過程中，基礎數據可能存在較大的不確定性，不同生長季的植物對水分的敏感程度也存在差異，在今后的研究中將會繼續對模型進行改進，增加基礎數據的不確定性和不同生長季配水的優先序，使之適用于不確定情形下治沙區水土資源優化配置。