基于系統(tǒng)動力學(xué)的民勤綠洲農(nóng)業(yè)系統(tǒng)碳排放仿真模擬研究

陳軍娟，燕振剛，李 薇，李 廣

(1.甘肅農(nóng)業(yè)大學(xué)信息科學(xué)技術(shù)學(xué)院，蘭州 730070；2.甘肅農(nóng)業(yè)大學(xué)財經(jīng)學(xué)院，蘭州 730070；3.甘肅農(nóng)業(yè)大學(xué)林學(xué)院，蘭州 730070)

【研究意義】近年來，溫室氣體所引起的氣候變化問題受到世界各國的廣泛關(guān)注。其中，農(nóng)業(yè)是產(chǎn)生溫室氣體的原因之一[1]。據(jù)已有的研究顯示，全球人為導(dǎo)致產(chǎn)生的溫室氣體里面農(nóng)業(yè)產(chǎn)生的溫室氣體占將近30%[2]，想要實現(xiàn)減少我國溫室氣體排放目標，還應(yīng)從治理農(nóng)業(yè)碳排放角度入手[2]。為促進我國經(jīng)濟的快速發(fā)展，減少碳排放對我國趁早實現(xiàn)2030年前實現(xiàn)碳達峰，2060年前實現(xiàn)碳中和提供重要意義和一定的理論參考價值?！厩叭搜芯窟M展】近年來，國內(nèi)外學(xué)者對農(nóng)業(yè)碳排放進行了廣泛的研究，相關(guān)研究主要在于農(nóng)業(yè)生產(chǎn)活動對碳排放的影響、土地利用產(chǎn)生碳排放的測算以及影響農(nóng)業(yè)碳排放因素的分析等方面，例如：李穎[4]以黃土高原退耕區(qū)縣南溝為例對農(nóng)業(yè)活動中碳產(chǎn)品顯化及市場對接進行了研究；舒璜[5]以南昌市為研究區(qū),基于2000—2017年土地利用動態(tài)變化，運用遙感解譯手段提取土地利用變化數(shù)據(jù),分析土地利用數(shù)量、結(jié)構(gòu)和利用程度的變化及趨勢；Edwards等[6]研究表明灌溉方式不同，也會影響農(nóng)業(yè)碳排放量不同；楊龍等[7]從基于2019年楊凌區(qū)及周邊的圍觀調(diào)查數(shù)據(jù)分析了種植業(yè)投入要素中土地投入、農(nóng)藥投入以及勞動力投入對農(nóng)業(yè)碳排放的影響，劉立平[8]利用Kaya恒等式分析了影響河南省農(nóng)業(yè)碳排放的因素，發(fā)現(xiàn)導(dǎo)致農(nóng)業(yè)碳排放量增加的主導(dǎo)因素是經(jīng)濟的發(fā)展。趙先超等[9]運用LMDI模型分析了湖南省農(nóng)業(yè)碳排放的影響因素，結(jié)果表明農(nóng)業(yè)經(jīng)濟水平和產(chǎn)業(yè)結(jié)構(gòu)是促進農(nóng)業(yè)碳排放量增加的主要因素。李彩弟等[10]對張掖市農(nóng)業(yè)碳排放影響因素進行了研究，研究結(jié)果表明經(jīng)濟因素會使農(nóng)業(yè)碳排放增加，而效率、結(jié)構(gòu)和勞動力因素能減少碳排放，陳林等[11]利用IPCC清單估算法對宜賓市2001—2015年農(nóng)業(yè)碳排放進行了估算，發(fā)現(xiàn)近15年來碳排放量呈增加趨勢，并且導(dǎo)致碳排放的主導(dǎo)因素是農(nóng)用物資和農(nóng)用能源因素?！颈狙芯壳腥朦c】目前在基于系統(tǒng)動力學(xué)農(nóng)業(yè)系統(tǒng)碳排放預(yù)測方面的研究較少，基于此，本文以基于民勤綠洲2004—2019年農(nóng)業(yè)生產(chǎn)等數(shù)據(jù)作為基礎(chǔ)，運用系統(tǒng)動力學(xué)仿真模擬來計算民勤綠洲2004—2019年耕地—農(nóng)業(yè)和家畜的碳排放量并預(yù)測未來11年碳排放變化趨勢?！緮M解決的關(guān)鍵問題】以期為民勤綠洲農(nóng)業(yè)生產(chǎn)引導(dǎo)低碳方式做出理論參考。

1 材料與方法

1.1 研究區(qū)域概況

民勤縣是甘肅省武威市下轄縣，位于甘肅省河西走廊的東北部(101°49′41″～104°12′10″E，38°3′45″～39°27′37″N)，在石羊河流域下流，南靠武威，西鄰鎳都金昌，東北和西北面與內(nèi)蒙古的左、右旗相接。東西長206 km，南北寬156 km，全縣總面積1.59萬km2，耕地面積63745.19萬km2，主要農(nóng)作物為春小麥、玉米和大麥等；農(nóng)作物播種面積達58.01千hm2，平均海拔1400 m，屬溫帶大陸性干旱氣候區(qū)，年均氣溫為9.2 ℃，全年無霜期210 d[12]。年平均降水量113.2 mm，降雨主要集中在每年6—9月，8月最多，占全年降水量的66%。民勤縣日照時間長，特別適宜農(nóng)作物生長，是甘肅省重要的商品糧食基地縣。

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

本實驗中所涉及的數(shù)據(jù)均來源于2004—2019年《甘肅統(tǒng)計年鑒》《甘肅發(fā)展年鑒》《甘肅農(nóng)村年鑒》，主要數(shù)據(jù)有農(nóng)作物的化肥、農(nóng)藥、農(nóng)膜、灌溉、翻耕、機械總動力等生產(chǎn)投入數(shù)據(jù)，其中化肥為實用量；翻耕面積為農(nóng)作物播種面積；灌溉為有效灌溉面積。

1.3 研究方法

系統(tǒng)動力學(xué)(System dynamics,簡稱SD)是一門分析研究信息反饋系統(tǒng)的學(xué)科，也是一門認識并解決系統(tǒng)問題的交叉綜合性學(xué)科，通過定性與定量相結(jié)合的方法研究系統(tǒng)復(fù)雜問題并相應(yīng)的建立系統(tǒng)模型[13]。農(nóng)業(yè)碳排放是指農(nóng)業(yè)生產(chǎn)過程中所產(chǎn)生的碳排放，對于民勤綠洲農(nóng)業(yè)系統(tǒng)碳排放的測算，農(nóng)業(yè)系統(tǒng)的碳排放主要由農(nóng)業(yè)原料的投入產(chǎn)生。其中農(nóng)業(yè)原料的投入包含了增加土壤肥力的化肥、用于害蟲防治的農(nóng)藥、農(nóng)膜、農(nóng)業(yè)灌溉以及灌溉消耗的農(nóng)機動力的使用等[14]。

ECcrop=Cf+Cm+Gi+Ca+Cp

(1)

式中，CEcrop是指農(nóng)田農(nóng)作物生產(chǎn)過程的總碳排放量；Cf是指農(nóng)田化肥施用所產(chǎn)生的碳排放；Cm是農(nóng)業(yè)機械使用過程中燃燒燃料所產(chǎn)生的碳排放；Ci是灌溉農(nóng)田過程中所產(chǎn)生的碳排放；Ca是農(nóng)膜使用過程中所產(chǎn)生的碳排放量，Cp是農(nóng)藥使用過程中所產(chǎn)的碳排放量。

Cf=Gf×Qf

(2)

Cm=Sp×Qp+Wm×Qm

(3)

Ci=Si×Qi

(4)

Ca=Ga×Qa

(5)

Cp=Gp×Ap

(6)

式(2)中，Gf指農(nóng)作物生產(chǎn)過程中所施用的化肥數(shù)量,Qf是指化肥的碳排放系數(shù)；式(3)中，農(nóng)業(yè)機械使用過程中產(chǎn)生的總碳排放量是Cm由農(nóng)作物種植面積翻耕和農(nóng)業(yè)機械動力兩部分組成，Sp是農(nóng)作物的種植面積，Wm表示農(nóng)業(yè)機械使用過程的總動力，Qp的碳排放轉(zhuǎn)換系數(shù)為16.47 kg/hm2；式(4)中，灌溉過程中所產(chǎn)生的碳排放量Ci由農(nóng)作物的灌溉面積Si和灌溉系數(shù)Qi相乘所得；式(5)中，農(nóng)膜使用所產(chǎn)生的碳排放量Ca是由農(nóng)膜使用量Ga與農(nóng)膜碳排放系數(shù)Qa相乘所得；式(6)中，農(nóng)藥使用所產(chǎn)生的碳排放量Cp是由農(nóng)藥的使用量Gp與農(nóng)藥碳排放系數(shù)Ap的乘積[14]。其中碳排放系數(shù)見表1。

表1 碳排放系數(shù)Table 1 Carbon emission coefficient

家畜的養(yǎng)殖過程中在農(nóng)業(yè)的生產(chǎn)中也產(chǎn)生碳排放，主要來源于飼料生產(chǎn)時的碳排放，家畜胃腸道碳排放，糞便碳排放以及飼料種植和加工時產(chǎn)生的碳排放[14]。民勤綠洲農(nóng)戶畜禽養(yǎng)殖有牛、羊和豬三類，主要以養(yǎng)羊為主。碳排放量以年末的存欄數(shù)量來計算，計算公式為：

C牲畜=∑iSi牲畜×Q牲畜

(7)

式中，C牲畜為牲畜的碳排放量，Si牲畜表示第i中牲畜的年末存欄頭數(shù)，Qi牲畜表示第i中牲畜的碳排放系，Q羊=0.075 t/a,Q牛=0.796 t/a，Q豬=0.082 t/a。

2 結(jié)果與分析

2.1 系統(tǒng)邊界的確定

空間邊界設(shè)定為民勤縣，時間邊界為：2004—2030年，模擬基期年為2004年，主要歷史數(shù)據(jù)時段為2004—2019年。模擬步長定為1年。

2.2 模型的構(gòu)建

2.2.1 農(nóng)業(yè)生產(chǎn)投入消耗碳排放的SD模擬模型構(gòu)建 利用系統(tǒng)動力學(xué)原理和方法構(gòu)建農(nóng)業(yè)生產(chǎn)投入碳排放的因果關(guān)系圖和存量流量圖，是農(nóng)業(yè)投入消耗碳排放動態(tài)模擬研究的重要組成部分。Vensim-PLE模型的因果樹可以簡潔明了地表達農(nóng)業(yè)投入消耗碳排放系統(tǒng)中重要元素之間的相互影響,而存量流量圖是在因果關(guān)系圖的基礎(chǔ)上進而明確以及設(shè)定系統(tǒng)變量和參數(shù)，農(nóng)業(yè)生產(chǎn)系統(tǒng)的因果樹和存量流量圖如圖1～2所示。農(nóng)業(yè)生產(chǎn)系統(tǒng)存量流量圖模塊中共有26個量，其中6個狀態(tài)變量，7個速率變量和13個輔助變量，模塊中還有一個隱藏變量‘Time’，其他為常量(表2)。

圖1 農(nóng)業(yè)生產(chǎn)系統(tǒng)的因果樹Fig.1 The cause tree of agricultural production system

圖2 農(nóng)業(yè)生產(chǎn)系統(tǒng)的存量流量Fig.2 Inventory flow of agricultural production system

表2 模型變量Table 2 Model variables

2.2.2 家畜碳排放的SD模擬模型構(gòu)建 家畜碳排放系統(tǒng)的因果樹和存量流量圖如圖3～4所示。家畜生產(chǎn)系統(tǒng)存量流量圖模塊中共有13個量，其中3個狀態(tài)變量，3個速率變量和7個輔助變量，模塊中還有一個隱藏變量‘Time’，其他為常量(表3)。

圖3 家畜產(chǎn)生碳排放的原因樹Fig.3 Causative tree of carbon emissions produced by livestock

圖4 家畜產(chǎn)生碳排放的存量流量Fig.4 Inventory flow of carbon emissions produced by livestock

表3 模型變量Table 3 Model variables

2.3 模型檢驗

模型建立后應(yīng)對模型進行直觀與運行檢驗、歷史仿真檢驗，以判斷模型的合理性、真實性、有效性以及穩(wěn)定性[14]。使用Vensim-PLE軟件的內(nèi)置方程式檢查功能可以直觀地檢查模型的合理性并運行。通過檢查和測試結(jié)果表明該模型沒有錯誤提示。因此，本文構(gòu)建的農(nóng)業(yè)系統(tǒng)碳排放模擬模型是合理的。

模型構(gòu)建完成之后往往達不到理想的精度，通常需經(jīng)過歷史仿真檢驗，精度在±10%內(nèi)可認為模型有效[14]。將2004—2019年的歷史統(tǒng)計數(shù)據(jù)代入模型進行模擬驗證，農(nóng)作物系統(tǒng)中選取農(nóng)業(yè)機械動力、農(nóng)膜和化肥作為檢驗變量，家畜系統(tǒng)中選取牛、羊、豬三類作為檢驗變量，將模型測算出的預(yù)測值與歷史實際值進行誤差計算，結(jié)果如表4～5所示。這些指標的相對誤差值都在±1%以內(nèi)，說明模型通過了有效性檢驗。這表明系統(tǒng)動力學(xué)模型的模擬結(jié)果是可靠的，不但可以通過調(diào)整關(guān)鍵參數(shù)進行模擬實驗，還能模擬碳排放的狀態(tài)和農(nóng)業(yè)系統(tǒng)的變化趨勢。

表4 2004—2019年農(nóng)作物投入數(shù)據(jù)歷史檢驗結(jié)果Table 4 Results of historical inspection of crop input data from 2004 to 2019

2.4 系統(tǒng)仿真結(jié)果分析

2.4.1 農(nóng)作物生產(chǎn)投入所產(chǎn)生的碳排放 通過對甘肅省民勤縣農(nóng)作物碳排放系統(tǒng)的仿真，得到民勤縣2004—2030年農(nóng)作物生產(chǎn)投入碳排放模擬仿真結(jié)果(表6，圖5)。將2004—2030年民勤農(nóng)作物生產(chǎn)投入所產(chǎn)生的碳排放分為3個階段，2004—2008年為第一階段，該階段的碳排放量增長較不穩(wěn)定。2009年由于農(nóng)作物播種面積和有效灌溉面積的減少，農(nóng)作物碳排放量相比2008年碳排放減少了238.09×104t，2009—2013年為第二階段，相比上一階段，該階段的總碳排放量比上一階段總碳排放量減少了1824×104t，第三階段2014—2030年碳排放呈一定趨勢緩慢增長，到2030年碳排放量達3277.93×104t，年均增長率為2.9%。

圖5 2004—2030年農(nóng)作物生產(chǎn)投入碳排放仿真結(jié)果Fig.5 Simulation results of carbon emissions from crop production inputs from 2004 to 2030

表5 2004—2019年家畜存欄數(shù)據(jù)歷史檢驗結(jié)果Table 5 Results of historical examination of livestock data from 2004 to 2019 (萬頭)

表6 2004—2030年農(nóng)作物生產(chǎn)投入碳排放Table 6 Carbon emissions from crop production inputs from 2004 to 2030 (萬t)

對農(nóng)作物生產(chǎn)投入所產(chǎn)生的碳排放進行仿真，在其它投入原料增加速率不變的情況下將農(nóng)業(yè)機械總動力的增加速率總體減少2.5%，則從仿真結(jié)果來看，2004—2025年碳排放變化趨勢相似，在2025—2030年碳排放變化趨勢相對平緩，到2030年碳排放量為2930.56×104t，比仿真前減少了347.37×104t，減排效果顯著。

2.4.2 家畜所產(chǎn)生的碳排放 通過對甘肅省民勤縣家畜碳排放系統(tǒng)的仿真，得到民勤縣2004—2030年家畜所產(chǎn)生的碳排放模擬仿真結(jié)果(表7，圖6)。

圖6 2004—2030年家畜所產(chǎn)生的碳排放仿真結(jié)果Fig.6 Simulation results of carbon emissions produced by livestock from 2004 to 2030

表7 2004—2030年家畜所產(chǎn)生的碳排放Table 7 Carbon emissions produced by livestock from 2004 to 2030 (萬t)

家畜的存欄數(shù)直接影響著家畜產(chǎn)生碳排放的變化趨勢，按照目前家畜碳排放增長趨勢看，家畜碳排放量從2018—2030年持續(xù)增長，這與我國提出的2030年前達到碳排放峰值的目標產(chǎn)生偏差，因此應(yīng)制定低碳發(fā)展規(guī)劃，以有效抑制家畜碳排放量的增長態(tài)勢[16]。

從Vensim-PLE軟件模擬家畜2004—2030年年末存欄數(shù)得知，家畜的存欄數(shù)在一定程度上都在增長，其中羊的增長速度最快，增長速率達85.6%。為了有效制約這種碳排放量的增長態(tài)勢，相對減少羊的存欄數(shù)改變增加速率，調(diào)整羊的增加速率之后2004—2025年家畜碳排放相比自然增長碳排放變化趨勢相似，但碳排放量總體稍高。綜合來看，到2025年2種情景的碳排放量值相等，且未來幾年調(diào)整羊的增加速率后家畜碳排放量增長緩慢且平穩(wěn)，到2030年調(diào)整后的碳排放量比自然增長碳排放量減少了約13.82×104t，年均增長率降低了1.9%，總的來說效果較顯著。

3 討 論

農(nóng)業(yè)碳排放研究一直是學(xué)術(shù)界碳排放研究的熱點話題之一，本文通過對民勤綠洲2004—2019年農(nóng)業(yè)生產(chǎn)碳排放值評估及變化趨勢進行分析，構(gòu)建了民勤綠洲農(nóng)業(yè)系統(tǒng)碳排放SD模型，通過設(shè)定作物原料投入使用增加速率和家畜增加速率進行動態(tài)仿真模擬分析，得到家畜的碳排放量年均增長速率比作物碳排放量年均增長速率高波動大，這與彭敏[14]研究的重慶市土地利用耕地—農(nóng)業(yè)生產(chǎn)系統(tǒng)是碳源系統(tǒng)結(jié)果一致。通過對比作物自然增長和調(diào)整作物農(nóng)業(yè)機械增加速率產(chǎn)生的碳排放量發(fā)現(xiàn)，到2030年調(diào)整后的作物碳排放量比調(diào)整前碳排放量減少了0.12倍，同樣通過調(diào)整家畜羊的增加速率調(diào)整后家畜的碳排放量比調(diào)整前減少了0.6倍。從統(tǒng)計年鑒數(shù)據(jù)來看，2004—2019年民勤綠洲農(nóng)作物播種面積有所減少，從62.66千hm2減少到58.01千hm2，出現(xiàn)這種趨勢的變化，原因有兩個：一是由于人口城鎮(zhèn)化，一部分農(nóng)田可能會被棄置和耕種，所以播種面積會略有下降；二是耕地面積減少，農(nóng)業(yè)產(chǎn)業(yè)結(jié)構(gòu)調(diào)整及國家推行的“退耕還林”導(dǎo)致了耕地面積的減少。通過對農(nóng)業(yè)生產(chǎn)系統(tǒng)所產(chǎn)生的碳排放仿真與分析，2004—2030年農(nóng)業(yè)碳排放總量呈穩(wěn)定增長的趨勢，到2030年農(nóng)作物生產(chǎn)投入碳排放量和家畜碳排放量分別達到2930.56萬t和23.33萬t，這與李波等[15]測算我國1993—2008年農(nóng)業(yè)碳排放量，發(fā)現(xiàn)自1993年以來我國農(nóng)業(yè)碳排放處于階段性的上升態(tài)勢一致。也與吳萌等[16]研究武漢市土地利用碳排放系統(tǒng)的仿真與分析結(jié)果是一致的。農(nóng)業(yè)生產(chǎn)投入的農(nóng)資數(shù)量與農(nóng)業(yè)產(chǎn)業(yè)結(jié)構(gòu)和農(nóng)業(yè)技術(shù)水平直接影響著農(nóng)業(yè)碳排放的增加，文中影響農(nóng)業(yè)碳排放增加的因素有農(nóng)用化肥、農(nóng)膜以及農(nóng)業(yè)機械的使用，其中主導(dǎo)因素是農(nóng)業(yè)機械。這與何艷秋等[17]研究的影響農(nóng)業(yè)碳排放主導(dǎo)因素的為農(nóng)地利用和農(nóng)業(yè)生產(chǎn)技術(shù)包括農(nóng)用機械，且種植業(yè)和畜牧業(yè)同時發(fā)展一致；同時，姚延婷等[18]發(fā)現(xiàn)農(nóng)業(yè)機械總動力、農(nóng)業(yè)柴油量、化肥施用量是農(nóng)業(yè)溫室氣體排放的主要原因。通過調(diào)整農(nóng)業(yè)機械總動力的增加速率，到2030年碳排放量為2930.56×104t，比仿真前減少了347.37×104t，減排效果顯著。目前運用系統(tǒng)動力學(xué)方法的研究內(nèi)容大多數(shù)集中在研究土地利用碳排放方面，例如：陳紫君等[19]用系統(tǒng)動力學(xué)方法研究了湖南省土地利用碳排放情景方案設(shè)定與仿真模擬，通過綜合情景分析，最終目標年土地利用碳排放量較趨勢發(fā)展情景減少了560.6萬t，減排效果顯著。陳亞麗[20]運用系統(tǒng)動力學(xué)方法對山西省土地利用碳排放進行了研究。此方法也被廣泛應(yīng)用于土地利用變化特征分析[21]、低碳城市發(fā)展研究[22-24]、碳排放預(yù)測研究[25]、能源消耗碳排放[26]以及旅游業(yè)碳排放模擬仿真等領(lǐng)域[27]。

4 結(jié) 論

基于系統(tǒng)動力學(xué)的原理和方法，本文利用Vensim-PLE軟件對民勤的耕地—農(nóng)業(yè)和家畜系統(tǒng)的碳排放進行了建模和模擬，并模擬了民勤2004—2030年耕地—農(nóng)業(yè)和家畜系統(tǒng)碳排放的發(fā)展趨勢。得出主要結(jié)論：農(nóng)業(yè)生產(chǎn)的原料投入使用量和家畜的年末存欄數(shù)直接影響著農(nóng)業(yè)碳排放量，農(nóng)業(yè)碳排放系統(tǒng)是一個復(fù)雜的系統(tǒng)，系統(tǒng)仿真過程是一個不斷調(diào)整模型參數(shù)和變量關(guān)系的過程，因此，仿真結(jié)果只代表一種可能性，并不是全部結(jié)果，目的是直觀地表達不同情景下農(nóng)業(yè)碳排放的變化趨勢，為民勤縣制定低碳農(nóng)業(yè)政策提供參考。通過基于系統(tǒng)動力學(xué)模型檢驗結(jié)果可知，將系統(tǒng)動力學(xué)模型用于農(nóng)業(yè)碳排放動態(tài)仿真模擬是可行的。本研究也存在一定的局限性。模型中變量的選擇和參數(shù)的設(shè)置存在一定的主觀性，如何與其他優(yōu)化模型相結(jié)合，更好地模擬農(nóng)業(yè)碳排放系統(tǒng)，還需要進一步深入研究。