宣化區(qū)農(nóng)業(yè)水資源系統(tǒng)仿真與優(yōu)化配置

張思俊，何 莉，吳 霜，杜 煜，張照壟

(湖北工業(yè)大學電氣與電子工程學院，湖北 武漢 430068)

宣化區(qū)位于河北省張家口市區(qū)東南，為冀西北壩下雜糧主產(chǎn)區(qū)。境內(nèi)干旱少雨，“缺水”是制約全區(qū)旱作區(qū)農(nóng)業(yè)生產(chǎn)發(fā)展的“瓶頸”[1]。近年來區(qū)內(nèi)最大常年性河流洋河由于生產(chǎn)、生活污水不斷排入導致水質(zhì)污染，河下游是冬奧會期間滑雪場館的重要后備供水水源地[2]，對宣化區(qū)進行農(nóng)業(yè)水資源優(yōu)化配置意義重大。張珊等[3]構(gòu)建以水土資源限制為約束的多目標模型進行區(qū)域農(nóng)業(yè)水資源配置。王學全等[4]構(gòu)建水資源優(yōu)化配置系統(tǒng)動力學模型，對青海湖流域水資源進行優(yōu)化配置。結(jié)合模擬技術(shù)(SD)與多目標規(guī)劃(MOP)方法，以宣化區(qū)為研究對象，建立農(nóng)業(yè)水資源系統(tǒng)動力學模型，進行仿真預測并確定出敏感參數(shù)。在此基礎(chǔ)上，考慮水質(zhì)水量為約束，農(nóng)業(yè)水資源供需差額最小與農(nóng)業(yè)總產(chǎn)值最大為目標，建立多目標規(guī)劃模型，得出宣化區(qū)農(nóng)業(yè)水資源優(yōu)化配置方案，為宣化區(qū)的水資源合理利用與農(nóng)業(yè)發(fā)展提供參考依據(jù)。

1 研究方法

系統(tǒng)動力學(SD)是通過系統(tǒng)分析利用變量之間反饋關(guān)系對未來趨勢進行預測的方法[5]，但難以優(yōu)化敏感參數(shù)。多目標規(guī)劃(MOP)能夠兼顧多個目標優(yōu)化決策變量。將SD與MOP相結(jié)合，可功能互補。具體步驟如下：

1)確定系統(tǒng)邊界，建立SD模型。根據(jù)研究對象及目標來確定系統(tǒng)邊界，分析影響宣化區(qū)農(nóng)業(yè)水資源系統(tǒng)的影響因素及其反饋關(guān)系，利用仿真軟件Vensim構(gòu)建農(nóng)業(yè)水資源系統(tǒng)流程圖。

2)參數(shù)估計并分析敏感參數(shù)。采用回歸分析、表函數(shù)、取平均值等多種方法確定系統(tǒng)參數(shù)取值。運行模型并進行歷史檢驗，直到模型具有良好的魯棒性與可信度，再分析挑選出符合建模目的及背景的敏感參數(shù)。

3)構(gòu)建MOP模型并求解。以敏感參數(shù)為決策變量，引入水質(zhì)水量約束，綜合考慮農(nóng)業(yè)水資源供需平衡與農(nóng)業(yè)經(jīng)濟效益最大化等目標，建立農(nóng)業(yè)水資源MOP模型，并利用Matlab非線性規(guī)劃工具箱求解。

4)優(yōu)化方案仿真。將MOP模型優(yōu)化解代入SD模型，得到優(yōu)化后的農(nóng)業(yè)水資源配置方案。

2 SD仿真模型的構(gòu)建

2.1 系統(tǒng)結(jié)構(gòu)與主要變量分析

以張家口市宣化區(qū)的行政界線為空間邊界；以2006～2030年作為時間邊界，2016年為現(xiàn)狀年，目標年設(shè)為2020年，仿真預測的年限設(shè)為2017-2030年，仿真時間步長設(shè)定為1年。模型數(shù)據(jù)來源主要為《張家口市統(tǒng)計年鑒》(2006-2016)，《官廳水庫生態(tài)安全調(diào)查與評估報告》(2014)，2006-2016年張家口市水資源公報以及2006-2016年張家口市國民經(jīng)濟和社會發(fā)展統(tǒng)計公報、《張家口市農(nóng)業(yè)發(fā)展“十三五”規(guī)劃》等。

根據(jù)對地區(qū)農(nóng)業(yè)供用水結(jié)構(gòu)的分析，將農(nóng)業(yè)水資源系統(tǒng)劃分為3個子系統(tǒng)，即水資源供需子系統(tǒng)、水污染子系統(tǒng)和農(nóng)業(yè)經(jīng)濟子系統(tǒng)，水資源供需子系統(tǒng)考慮農(nóng)業(yè)可供水量及農(nóng)業(yè)用水量對系統(tǒng)的影響，總種植面積與牲畜存欄總量為狀態(tài)變量；水污染子系統(tǒng)考慮農(nóng)業(yè)污染即種植業(yè)、畜牧業(yè)所產(chǎn)生的COD入河量，其他如工業(yè)、生活所產(chǎn)生的COD入河量處理為常數(shù)或表函數(shù)，洋河COD負荷量為狀態(tài)變量；農(nóng)業(yè)經(jīng)濟子系統(tǒng)考慮種植業(yè)、畜牧業(yè)，而漁業(yè)、林業(yè)產(chǎn)值份額小不予考慮，各作物及牲畜單位用水效益處理成表函數(shù)。依據(jù)各子系統(tǒng)因果關(guān)系，運用Vensim軟件構(gòu)建農(nóng)業(yè)水資源系統(tǒng)動力學流程圖(圖1)。

(a)水資源供需子系統(tǒng)

(b)水污染子系統(tǒng)

(c)農(nóng)業(yè)經(jīng)濟子系統(tǒng)圖1 農(nóng)業(yè)水資源系統(tǒng)動力學模型流程圖

2.2 參數(shù)設(shè)定及模型檢驗

圖1中，模型狀態(tài)變量如總種植面積、牲畜存欄總量、洋河COD負荷量等的初始值采用2006年的真實值。常量如糧食作物用水排污系數(shù)、大牲畜用水排污系數(shù)、種植業(yè)污染物入河系數(shù)等由歷史數(shù)據(jù)取平均值得到。其他隨時間變化的變量如總種植面積變化率、糧食作物單位用水效益(糧食作物單位用水量、糧食種植面積占比)等變量運用回歸分析法或表函數(shù)法確定其取值。

以2006年作仿真起始年，選擇種植業(yè)用水量、畜牧業(yè)用水量、種植業(yè)總產(chǎn)值、畜牧業(yè)總產(chǎn)值等指標作為觀測變量進行歷史性檢驗，以考察模型的有效性與可信度。運行模型后，將觀測變量2006年、2009年、2012年、2016年的仿真值與實際值進行比較見表1，可以看到，各變量的仿真值與實際值誤差均在10%以內(nèi)，模型具有較好的魯棒性和可信度，可對水資源系統(tǒng)進行分析以及相關(guān)的預測。

2.3 敏感性分析

模型的敏感性分析是指改變模型參數(shù)取值，通過比較改變前后運行結(jié)果，分析參數(shù)對模型的影響。確定敏感參數(shù)需要以政策分析為指導，通過分析建模目的、研究背景、系統(tǒng)的因果關(guān)系等方面內(nèi)容來選定[6]。

表1 歷史檢驗

本文旨在通過優(yōu)化種植業(yè)、畜牧業(yè)水資源配置來緩解地區(qū)水資源供需矛盾與水質(zhì)污染。通過排除政策調(diào)控范疇外與對目標影響較小的參數(shù)，最終確定模型的敏感參數(shù)為糧食、油料、蔬菜、水果作物的單位用水量及種植面積占比，大牲畜、小牲畜的單位用水量及養(yǎng)殖數(shù)量占比。

3 MOP模型的構(gòu)建

將上述分析得到的敏感參數(shù)作為決策變量，以洋河COD濃度達標及供水能力為約束，農(nóng)業(yè)水資源供需差額最小、農(nóng)業(yè)總產(chǎn)值最大為目標，建立農(nóng)業(yè)水資源MOP模型。

目標函數(shù)：

1)農(nóng)業(yè)總產(chǎn)值最大化目標：

(1)

2)農(nóng)業(yè)水資源供需差額最小化目標：

(2)

約束條件：

1)等式約束

(3)

2)非線性約束

洋河COD濃度約束：

(4)

供水能力約束：

f2≥0

(5)

3)上下限約束

分析決策變量的歷史值，結(jié)合實際進行估計得到其取值范圍(表2)。

表2 決策變量取值范圍

式(1)到式(5)中：Pi表示第i種作物的單位用水效益(元/m3)；Ci表示第i種牲畜的單位用水效益(元/m3)；xi為第i種作物單位用水量(萬m3/hm2)；yi為第i種作物種植面積占比；ai為第i種牲畜單位用水量(m3/頭)；bi為第i種牲畜養(yǎng)殖數(shù)量占比；S為總種植面積(萬hm2)；L為總牲畜數(shù)量(萬頭)；γ1為種植業(yè)COD入河系數(shù)、γ2為畜牧業(yè)COD入河系數(shù)；αi分別為第i種作物COD排放系數(shù)(t/萬m3)、βi第i種牲畜COD排放系數(shù)(t/萬m3)；R1、R2分別為工業(yè)與生活COD處理入河量、工業(yè)與生活COD直排入河量(t)；W5洋河水量(萬m3)；W表示農(nóng)業(yè)可供水量(萬m3)。模型部分系數(shù)的取值見表3，由歷史數(shù)據(jù)計算取值得到。

表3 MOP模型系數(shù)取值

4 結(jié)果與分析

4.1 MOP模型有效性驗證

模型的目標農(nóng)業(yè)水資源供需差額、農(nóng)業(yè)總產(chǎn)值計量單位不同，其中農(nóng)業(yè)總產(chǎn)值為效益型屬性值，農(nóng)業(yè)水資源供需差額為偏離型屬性值，所以對各目標進行歸一化處理[7]。使用線性加權(quán)法將多目標問題轉(zhuǎn)化為單目標問題進行MOP模型的求解。構(gòu)造總體評價函數(shù)：

(6)

式中：f1、f2為歸一化后的農(nóng)業(yè)總產(chǎn)值、農(nóng)業(yè)水資源供需差額；λ1為農(nóng)業(yè)總產(chǎn)值權(quán)重，λ2為農(nóng)業(yè)水資源供需差額權(quán)重。

以2016年數(shù)據(jù)為例進行MOP模型有效性驗證，在Matlab軟件中利用非線性規(guī)劃算法工具箱求解MOP模型，得到最優(yōu)解集(即SD模型最優(yōu)敏感參數(shù))。表4為MOP優(yōu)化前后敏感參數(shù)及目標值對比。

表4中，2016年優(yōu)化前后糧食作物種植占比增加10.47%，蔬菜作物單位用水量增加20.55%；大、小牲畜的單位用水量分別減少22.24%、17.06%，根據(jù)表4中各作物及牲畜的用水效益及用水排污系數(shù)，各決策變量朝著高效益低排污方向調(diào)整。通過MOP優(yōu)化后2016年的農(nóng)業(yè)水資源供需差額降低了273.38萬m3，水供需矛盾得到緩解，洋河COD濃度降到20 mg/L達到水質(zhì)標準，農(nóng)業(yè)總產(chǎn)值降低了6.63%，經(jīng)濟效益略有犧牲。

表4 MOP優(yōu)化下2016年敏感參數(shù)及目標值對比

4.2 目標年水資源優(yōu)化配置

在Vensim軟件中運行原SD模型，得到2020年的SD模型各項參數(shù)取值后，運用MOP模型求得2020年的最優(yōu)敏感參數(shù)，進而得到2020年的最優(yōu)農(nóng)業(yè)水資源配置方案。2020年宣化區(qū)的農(nóng)業(yè)水資源配置優(yōu)化前后對比見圖2。各目標量的變化見表5。

圖2 2020年宣化區(qū)農(nóng)業(yè)水資源配置優(yōu)化前后對比

表5 2020年各目標變化

目標單位SDSD-MOP水資源供需差額萬m3-698.22-1.04洋河COD濃度mg/L20.5519.99農(nóng)業(yè)總產(chǎn)值萬元693736610574

從圖2可以看到，對比兩種方案，2020年SD-MOP方案中油料作物、大牲畜配水量分別減少了39.79%、24.87%，蔬菜配水量增加35.58%。從表5可看到，優(yōu)化后2020年農(nóng)業(yè)水資源供需差額降低了697.18萬m3，洋河COD濃度降低了2.75%，農(nóng)業(yè)總產(chǎn)值與2019年相比仍保持10.77%的增長率。通過SD-MOP水資源優(yōu)化配置，宣化區(qū)水資源供需矛盾得到有效緩解，水質(zhì)污染得到有效控制。

5 結(jié)論

針對農(nóng)業(yè)水資源系統(tǒng)，分析系統(tǒng)要素及其反饋關(guān)系，以宣化區(qū)為研究對象，建立SD模型進行地區(qū)農(nóng)業(yè)水資源系統(tǒng)變化仿真及趨勢預測，并建立了多目標規(guī)劃模型,求解得出農(nóng)業(yè)水資源優(yōu)化配置方案，優(yōu)化方案在滿足水質(zhì)達標及水量約束的前提下，合理的降低高排污低效益的油料、牲畜的水資源供給量，加大蔬菜等經(jīng)濟作物的水資源供給，實現(xiàn)水資源與農(nóng)業(yè)經(jīng)濟協(xié)調(diào)發(fā)展。同時將仿真模擬技術(shù)與多目標優(yōu)化方法相結(jié)合，發(fā)揮優(yōu)化搜索能力的同時，發(fā)揮仿真模型的可靠性與模擬性優(yōu)勢，使決策方案更具可行性。為后續(xù)農(nóng)業(yè)水資源優(yōu)化配置研究提供參考借鑒。