馮 釗，張凱旋

（1.清華大學(xué)地球系統(tǒng)科學(xué)系/地球系統(tǒng)數(shù)值模擬教育部重點(diǎn)實(shí)驗(yàn)室，北京 100084；2.全球變化與中國綠色發(fā)展協(xié)同創(chuàng)新中心，北京 100875；3.西北農(nóng)林科技大學(xué)林學(xué)院，陜西 咸陽 712100）

【研究意義】農(nóng)業(yè)健康發(fā)展是我國“十三五”規(guī)劃的明確要求。Logistic模型方程的產(chǎn)生之初是用于描述人口及動植物、微生物種群的數(shù)量增長。經(jīng)近兩百年的發(fā)展，現(xiàn)被用于經(jīng)濟(jì)增長、農(nóng)業(yè)發(fā)展、生態(tài)系統(tǒng)評估等多個方面。該曲線形如“S”型，能基于統(tǒng)計觀測數(shù)據(jù)預(yù)估研究對象在未來的數(shù)量變化趨勢及環(huán)境最大容納量（K值）?；谠撃Ｐ头治龅霓r(nóng)業(yè)水土資源等指標(biāo)變化，可為當(dāng)?shù)剞r(nóng)業(yè)水資源部門、土地管理部門提供水土資源利用的參考信息，以促進(jìn)糧食安全生產(chǎn)?！厩叭搜芯窟M(jìn)展】在農(nóng)業(yè)生態(tài)系統(tǒng)中，產(chǎn)生“S”型的限制因子有多種［1］，如管理水平、水土資源等［2］。在曲線的農(nóng)業(yè)利用方面，喬瑞波等［3］基于Logistic模型評估浙江省耕種面積的極值變化；陳淇等［4］利用Logistic模型測算農(nóng)業(yè)土地面積流轉(zhuǎn)；徐建英等［5］利用Logistic模型分析臥龍地區(qū)土地資源的流轉(zhuǎn)；呂德東等［6］利用Logistic模型分析農(nóng)作物大豆中豆莢結(jié)構(gòu)與食心蟲之間的抗性關(guān)系；Shabani等［7］利用雙Logistic模型估算在水脅迫條件下，模擬農(nóng)作物冬油菜生長與干物質(zhì)的變化量；陳其清等［8］利用Logistic方程計算農(nóng)業(yè)經(jīng)濟(jì)變化。水稻是我國重要的糧食作物之一，從東北到海南均有分布。近年來，市場對東北大米偏好增加，加速了東北水稻的擴(kuò)大種植［9］，導(dǎo)致近20年種植面積翻兩番［10］，并且總產(chǎn)量持續(xù)增長。黑龍江省較遼寧省和吉林省增長更為突出。水稻對農(nóng)業(yè)水土資源的需求導(dǎo)致黑龍江農(nóng)業(yè)生態(tài)系統(tǒng)不斷調(diào)整種植結(jié)構(gòu)、部分地市將旱地變成水澆地，出現(xiàn)“旱改水”［11］的現(xiàn)象。Liu等［2］利用數(shù)學(xué)模型評價了黑龍江三江平原地區(qū)水土資源恢復(fù)力。

【本研究切入點(diǎn)】利用Logistic模型（數(shù)學(xué)模型）和RMSE模型（檢驗(yàn)?zāi)Ｐ停┑慕y(tǒng)計方法，結(jié)合黑龍江農(nóng)業(yè)生態(tài)系統(tǒng)近20年水土資源利用的相關(guān)指標(biāo)。深入研究黑龍江省水土資源的利用情況和變化過程?！緮M解決的關(guān)鍵問題】在當(dāng)前農(nóng)業(yè)科技管理水平下，分析糧食產(chǎn)量、種植面積、農(nóng)業(yè)水資源供應(yīng)等指標(biāo)的趨勢變化，并擬合、計算、評估出單產(chǎn)、總產(chǎn)、種植面積等指標(biāo)極值與時間的相關(guān)關(guān)系［12］，驗(yàn)證Logistic模型可應(yīng)用于農(nóng)業(yè)生態(tài)系統(tǒng)中大尺度、多指標(biāo)的變化，同時探究國家農(nóng)業(yè)政策、農(nóng)業(yè)管理信息、農(nóng)業(yè)科技水平和氣候條件等內(nèi)外限制因子對農(nóng)業(yè)生態(tài)系統(tǒng)的影響。

1 研究區(qū)域概況

1.1 糧食作物種植概況

黑龍江境內(nèi)有松嫩平原、三江平原等平原地區(qū)，盛產(chǎn)大豆、玉米、水稻等糧食作物。該地區(qū)年降雨量為400～650 mm，年均溫為-5～5℃，初霜凍一般出現(xiàn)在9月份，終霜凍一般在4月底至5月初，全年平均積溫的波動范圍為1 800～2 800℃（在≥10℃的條件下）［13］。由于農(nóng)作物對溫度和水資源需求，水稻種植區(qū)域多集中在該省東南部［14］。1985—2016年《國家統(tǒng)計年鑒》數(shù)據(jù)表明，1985年水稻種植總面積約占全國的5%，2015年占10%。2016年水稻種植總面積約為3.20×106hm2，約占全國水稻種植總面積（3.02×107hm2）的10.6%；玉米種植總面積為（5.22×106hm2），約占全國玉米總種植面積（3.68×107hm2）的14.2%；大豆種植總面積為3.05×106hm2，約占全國大豆種植總面積（9.70×106hm2）的31.4%［14］。2017年，黑龍江省糧食總產(chǎn)量在全國省份排名第一，其中水稻種植總面積同比增幅最大。

1.2 農(nóng)業(yè)水資源使用概況

農(nóng)業(yè)供水來源包括地表水和地下水。降雨影響該地區(qū)的地表水和地下水安全。王振芬等［15］研究表明，1961—2011年黑龍江降雨量出現(xiàn)逐漸減少的趨勢。由于該地區(qū)耗水作物種植面積擴(kuò)大，農(nóng)業(yè)水資源的供應(yīng)量在近20年連年攀升，但地表水的供應(yīng)量未出現(xiàn)明顯增加，導(dǎo)致了地下水的加速開采，以彌補(bǔ)農(nóng)業(yè)對水資源的需求。截至2014年，地下水的供應(yīng)量約為2004年的1.5倍。據(jù)預(yù)測，該地區(qū)農(nóng)業(yè)“旱改水”的持續(xù)增加將加快水資源消耗，在未來甚至?xí)霈F(xiàn)水資源匱乏、極端干旱等現(xiàn)象［16］。

2 材料與方法

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

農(nóng)業(yè)生態(tài)系統(tǒng)中作物產(chǎn)量、面積等數(shù)據(jù)來源于1995—2017年《黑龍江省統(tǒng)計年鑒》和《國家統(tǒng)計年鑒2017》。農(nóng)業(yè)用水?dāng)?shù)據(jù)來源于1995—2017年《黑龍江統(tǒng)計年鑒》和2005—2017年《黑龍江水利年鑒》。

2.2 研究方法

Logistic模型可將研究對象的變化范圍分成開始期、加速期、轉(zhuǎn)折期、減速期和飽和期5個時期。本研究中，開始期研究對象數(shù)量發(fā)展較緩慢，不斷適應(yīng)外界條件的變化；加速期是基于開始期發(fā)生數(shù)量上的急劇增加，增長率逐漸達(dá)到最大；轉(zhuǎn)折期是介于加速期和減速期之間的分界點(diǎn)；減速期是研究對象增加量逐漸減小的時期；而在飽和期則研究對象數(shù)量趨于穩(wěn)定，出現(xiàn)最大環(huán)境容納量（K值），發(fā)展達(dá)到最大限度［17］。Logistic模型的具體表達(dá)式為：式中，N表示在生態(tài)系統(tǒng)中研究對象的數(shù)量大??；K值表示環(huán)境容納量極大值；t表示時間；r為不同Logistic曲線下的常數(shù)，能夠影響“S”型曲線坡度。

3 結(jié)果與分析

3.1 土地資源利用變化分析

3.1.1 水稻和糧食作物種植面積的變化 從圖1A可見，水稻種植總面積（Rice total Plating Area，RPA）的增長期大約開始于2002年。2002—2015年出現(xiàn)明顯的加速期、轉(zhuǎn)折期和減速期。水稻種植面積在2015年開始出現(xiàn)拐點(diǎn)，極值在2020年出現(xiàn)，大小約為4.05×106hm2。對比水稻種植面積變化，從圖1B可知，糧食作物種植總面積（Crop total Planting Area，CPA）增長趨勢與水稻的增長趨勢相近，增長期也開始于2002年，但減速期拐點(diǎn)在2014年產(chǎn)生，而飽和期在2016年出現(xiàn)，產(chǎn)生極值的大小約為1.43×107hm2?？梢姡Z食作物的種植總面積K值出現(xiàn)的年份稍早于水稻種植總面積出現(xiàn)的年份。

圖1 1996—2016年黑龍江省水稻（A）及糧食作物（B）種植總面積變化Fig. 1 Changes in RPA（A） and CPA（B） in Heilongjiang Province from 1996 to 2016

3.1.2 水稻和糧食作物單產(chǎn)的變化 單產(chǎn)是農(nóng)田系統(tǒng)管理水平、物候條件、土壤肥力的綜合體現(xiàn)。由圖2可知，加速期開始于2008年，減速期止于2014年，糧食作物單產(chǎn)（Crops Per unit Yield，CPY）的極值約為5 275 kg/hm2。由圖3可知，水稻單產(chǎn)（Rice Per Unit Yield，RPY）在1980年以前就進(jìn)入加速期，進(jìn)入飽和期年份約為2016年，預(yù)期得到水稻單產(chǎn)極值為7 030 kg/hm2。對比糧食作物單產(chǎn)的變化，水稻單產(chǎn)的增長時間長，增長期持續(xù)時間長，最終進(jìn)入穩(wěn)定期的年份滯后，“S”型曲線的坡度比作物單產(chǎn)的坡度小。

圖2 1996—2016年黑龍江省糧食作物單產(chǎn)變化Fig. 2 Changes in CPY in Heilongjiang Province from 1996 to 2016

圖3 1980—2016年黑龍江省水稻單產(chǎn)變化Fig. 3 Changes in RPY in Heilongjiang Province from 1980 to 2016

圖4 1996—2016年黑龍江糧食作物總產(chǎn)量（A）與2004—2016年農(nóng)業(yè)總供水量（B）變化Fig. 4 Changes in CTY（A） in Heilongjiang Province from 1996 to 2016 and ATW（B） from 2004 to 2016

3.2 農(nóng)業(yè)水資源利用和糧食總產(chǎn)量變化結(jié)果

從圖4A可以看出，糧食作物總產(chǎn)量（Crop Total Yield，CTY）的變化趨勢與圖1A相同，在2004年開始進(jìn)入加速期，在2014年進(jìn)入拐點(diǎn)，達(dá)到穩(wěn)定期。此外，糧食作物總產(chǎn)量極值在2020年產(chǎn)生，約為6.39×107t。從圖4B可知，水資源供水極值與水稻的變化趨勢基本一致。供水極值在2016年左右產(chǎn)生，約為3.16×108m3。2016年農(nóng)業(yè)水資源總供給量（Agricultural Total Water supply，ATW）已進(jìn)入曲線的飽和期，表明黑龍江省農(nóng)業(yè)供水水平已達(dá)上限。

3.3 水土資源指標(biāo)模型驗(yàn)證及相關(guān)分析

為進(jìn)一步驗(yàn)證Logistic模型的可靠性，通過國家和黑龍江省長期統(tǒng)計結(jié)果與模型模擬數(shù)據(jù)的比較分析，計算Logistic模型在黑龍江農(nóng)業(yè)生態(tài)系統(tǒng)水土資源和部分農(nóng)業(yè)指標(biāo)模擬的符合度。同時采用相對均方根誤差（Root Mean Square Error，RMSE） 模型檢驗(yàn)方法評估Logistic模擬值與統(tǒng)計值的符合度， 其中R2為擬合優(yōu)度（決定系數(shù)），n為統(tǒng)計數(shù)量，r為Pearson相關(guān)系數(shù)。當(dāng) RMSE 值小于10%時，模擬值與統(tǒng)計值一致性非常好；當(dāng) RMSE值在 10%～20%之間時，表明二者一致性較好；當(dāng) RMSE 值在 20%～30% 時，二者一致性一般；而當(dāng) RMSE 值大于 30%時，表明二者存在較大偏差，一致性較差［18］。

不同類別的Logistic模型模擬檢驗(yàn)結(jié)果如表1 所示。由表1可知，Logistic模型模擬的6種類別的統(tǒng)計值與模擬值間相關(guān)性極好。6種類別的RMSE 值均小于 10%，說明二者間存在非常好的一致性。6種類別的擬合優(yōu)度分別為0.97、0.98、0.91、0.92、0.97和0.89，都達(dá)到了0.01水平下的極顯著相關(guān)關(guān)系?？傮w而言，Logistic模型的模擬效果較好，利用該模型模擬黑龍江省水土資源和部分農(nóng)業(yè)指標(biāo)的變化趨勢是可行的。

表1 不同類別模型與統(tǒng)計檢驗(yàn)結(jié)果Table1 Different categories and statistical test results

3.4 不同積溫帶地市土地資源利用變化分析

黑龍江農(nóng)業(yè)信息網(wǎng)公布信息表明，黑龍江地區(qū)有多種積溫帶分布。其中，第一、第二積溫帶在西南和東南均有分布，第三、第四積溫帶在黑龍江省中間分布，第五、第六積溫帶大體位于西北方，部分位于東北方。溫度差異決定了積溫帶的分布不同，也決定了種植作物分布差異。該地區(qū)水稻、玉米的生長積溫均值約為2 500℃，大豆、小麥的生長積溫均值約為1 500℃。根據(jù)《黑龍江統(tǒng)計年鑒2000—2017》的數(shù)據(jù)和不同積溫帶分布，可得典型積溫帶不同地市的土地利用數(shù)值與Logistic模型之間的關(guān)系。從表2可見，六地市中不同類型（面積和單產(chǎn)）指標(biāo)大體符合Logistic模型。其中，種植面積的擬合優(yōu)度明顯高于產(chǎn)量的，且前者模擬值和統(tǒng)計值二者間具有較強(qiáng)的相關(guān)性。而糧食作物單產(chǎn)、水稻單產(chǎn)的擬合優(yōu)度范圍跨度較大，存在部分指標(biāo)的擬合優(yōu)度偏低，如伊春糧食作物單產(chǎn)。均方根誤差RMSE也反映出不同積溫帶地市的統(tǒng)計值與Logistic模型之間的差值，也是重要的參考指標(biāo)。從R2和r、RMSE等3種統(tǒng)計指標(biāo)分析可得，種植面積變化的擬合程度高于單產(chǎn)、總產(chǎn)的變化。伊春的水稻單產(chǎn)和大興安嶺的糧食單產(chǎn)變化不遵循Logistic變化。

表2 典型積溫地市帶不同類別模型與統(tǒng)計檢驗(yàn)結(jié)果Table2 Different categories and statistical test results of typical Accumulated Temperature Zones（ATZ）

4 討論

4.1 水土資源利用與產(chǎn)生“S”型變化的主要原因

糧食全要素生產(chǎn)率和技術(shù)效率［19］影響農(nóng)業(yè)的發(fā)展，是農(nóng)業(yè)生產(chǎn)水平的反映。糧食單產(chǎn)、總產(chǎn)的提高和種植面積的優(yōu)化也是衡量生產(chǎn)效率和技術(shù)的重要參考指標(biāo)。

黑龍江農(nóng)業(yè)水土資源產(chǎn)生的“S”型變化與該地區(qū)農(nóng)業(yè)發(fā)展的限制因子有關(guān)。國家政策（經(jīng)濟(jì)政策、水資源政策等）、農(nóng)業(yè)科技水平（機(jī)械化、基因工程育種等）和農(nóng)業(yè)管理（施肥管理、灌溉等）是形成該曲線變化的主要限制因子。以水稻為例，國家從2004年恢復(fù)對黑龍江地區(qū)的最優(yōu)保底政策，并同時實(shí)行“一免兩補(bǔ)”，有效刺激水稻種植總面積的增長，該現(xiàn)象的出現(xiàn)與圖1A中的變化趨勢保持一致。此外，氣候變暖、抗凍品種（如富士光和龍稻5號等稻米）的培育和農(nóng)業(yè)災(zāi)害防治水平的提高也有利于水稻單產(chǎn)和總產(chǎn)的持續(xù)增加。從水稻單產(chǎn)變化的時間分析可知，水稻單產(chǎn)增長始于20世紀(jì)80年代，近40年穩(wěn)步提高，這種變化趨勢與農(nóng)業(yè)技術(shù)發(fā)展、變革的時間基本一致。另外，水稻種植面積的極值，既受外部因子的影響，也受地形、溫度、水資源供應(yīng)等內(nèi)部限制因子的制約，Logistic模型中K值的出現(xiàn)是該農(nóng)業(yè)生態(tài)系統(tǒng)發(fā)展的必然結(jié)果。

對農(nóng)業(yè)水資源利用分析可知，水稻是對水資源需求最高的糧食作物。水稻種植面積影響農(nóng)業(yè)水資源總量的供給變化，水稻種植總面積的“S”型變化也決定了水資源供應(yīng)的“S”型變化。若未來水資源灌溉技術(shù)出現(xiàn)革新，那么農(nóng)業(yè)水資源供應(yīng)量將在現(xiàn)有的“S”型變化曲線的基礎(chǔ)上發(fā)生明顯下降。而其變化依然是形成一個新的“S”型的遞減規(guī)律，最終產(chǎn)生新的K值。

黑龍江近20年主要糧食作物的單產(chǎn)，為水稻＞玉米＞小麥＞大豆。當(dāng)前糧食作物單產(chǎn)的“S”型變化是基于當(dāng)前糧食作物的種植比例下統(tǒng)計學(xué)指標(biāo)均值的體現(xiàn)，是降雨、溫度和土壤質(zhì)地等內(nèi)部因子和外部國家政策、管理水平和農(nóng)業(yè)科技水平等因子綜合作用的結(jié)果。當(dāng)未來國家政策、農(nóng)業(yè)科技水平和農(nóng)業(yè)管理中有任何一個限制因子在當(dāng)前基礎(chǔ)上發(fā)生變化，農(nóng)業(yè)生態(tài)系統(tǒng)的“S”型曲線將在現(xiàn)有基礎(chǔ)上發(fā)生變化，進(jìn)入下一個Logistic周期。

4.2 Logistic模型與RMSE的不確定性

Logistic方程是基于現(xiàn)有數(shù)據(jù)而產(chǎn)生的經(jīng)驗(yàn)?zāi)Ｐ?。對于預(yù)測農(nóng)業(yè)生態(tài)系統(tǒng)指標(biāo)變化，其前提條件是，該農(nóng)業(yè)生態(tài)系統(tǒng)中農(nóng)業(yè)管理、農(nóng)業(yè)科技水平、氣候條件等不會發(fā)生明顯變化，即依據(jù)當(dāng)前氣候和農(nóng)業(yè)管理、資源水平而產(chǎn)生的預(yù)測極值。不確定性主要來源于：（1）統(tǒng)計學(xué)誤差，即已獲得統(tǒng)計學(xué)數(shù)據(jù)存在人為統(tǒng)計偏差及系統(tǒng)誤差；（2）未來的農(nóng)業(yè)生態(tài)系統(tǒng)科技管理水平提高和氣候變化等限制因子改變。均方根誤差RMSE常用于衡量統(tǒng)計值與模擬值之間的誤差大小，通常與R2等指標(biāo)同時出現(xiàn)。在研究典型積溫帶中，佳木斯的作物種植面積出現(xiàn)RMSE值小、而R2也偏小的情況。這可能與統(tǒng)計值數(shù)量以及模型的時間尺度有關(guān)。因此，衡量模型的擬合優(yōu)劣，需要評估多種指標(biāo)，否則也會出現(xiàn)類似于RMSE失真的可能性。

4.3 未來氣候變化的水土資源利用與K值的關(guān)系

據(jù)IPCC第五次研究報告表明，到2100年，地球平均溫度將比當(dāng)前提高約0.3～4.8℃［20］，但東北地區(qū)降雨量在未來無明顯變化。根據(jù)Thomthwait方程［21］，當(dāng)溫度升高，蒸騰作用加劇，作物生長周期的需水量也隨之加劇。在現(xiàn)有灌溉水量條件下，水稻、大豆等耗水量偏大的糧食作物產(chǎn)量發(fā)生水脅迫［22］，在未來產(chǎn)值可能低于預(yù)期K值。若保證糧食產(chǎn)量不變，則種植面積K值會相應(yīng)減小。而市場對東北水稻的需求、國家政策對大豆總產(chǎn)的需求都將增大作物種植面積的K值。若未來氣候變化條件不明顯，農(nóng)業(yè)生態(tài)系統(tǒng)的最大值也依然會受到突發(fā)因素（土壤質(zhì)地、洪澇、部分干旱［12］等）的干擾，其真實(shí)產(chǎn)量N將圍繞Logistic模型K值上下波動［17］。

4.4 Logistic模型周期與K值的關(guān)系

Logistic模型K值的產(chǎn)生是生態(tài)系統(tǒng)限定因子，包含國家政策、農(nóng)業(yè)科技水平和農(nóng)業(yè)管理等綜合作用的結(jié)果。當(dāng)農(nóng)業(yè)系統(tǒng)中限定因子發(fā)生改變，K值也將相應(yīng)地發(fā)生變化?；?、水利設(shè)施、農(nóng)業(yè)管理和種植品種均為該農(nóng)業(yè)生態(tài)系統(tǒng)外部的限制因子，而土壤質(zhì)地、氣候、積溫帶等則屬于該系統(tǒng)內(nèi)部的限制因子。因此，當(dāng)某種限制因子不再制約農(nóng)業(yè)的發(fā)展時，農(nóng)業(yè)生態(tài)系統(tǒng)各變量將擬合一個新的Logistic方程，并以當(dāng)前時期為開始期，進(jìn)入下一個“S”型變化周期，依據(jù)外部變化可增可減，從而產(chǎn)生新的K值。

4.5 黑龍江農(nóng)業(yè)系統(tǒng)水土資源的潛在風(fēng)險

根據(jù)Logistic方程分析，水稻種植面積極值的產(chǎn)生時間晚于農(nóng)業(yè)水資源用量極值的產(chǎn)生時間，因此，該系統(tǒng)可能面臨水資源虧損的風(fēng)險。在未來氣候變化條件下，溫度升高和降雨的減少也會帶來水資源虧損的風(fēng)險。另外，長時間的“旱改水”導(dǎo)致旱地土壤細(xì)菌、真菌和放線菌比值與總量的變化，將影響土壤的特性。國家化肥零增長［23］政策一方面有助于降低化土壤性質(zhì)變化的風(fēng)險，但另一方面也給維持糧食高產(chǎn)量帶來了巨大的挑戰(zhàn)。

5 結(jié)論

（1）本研究采用 RMSE模型檢驗(yàn)的方法，對統(tǒng)計值與Logistic模擬值之間進(jìn)行了擬合相關(guān)性分析，從模型的模擬值與統(tǒng)計值間的 RMSE 值可知，Logistic模型可以較好地模擬黑龍江水土資源的動態(tài)變化。

（2）根據(jù)模型評估分析表明，黑龍江農(nóng)業(yè)生態(tài)系統(tǒng)農(nóng)業(yè)水資源供應(yīng)量在2016年產(chǎn)生極值，約為3.16×108m3。糧食作物種植總面積在2016年達(dá)到極值，約為1.43×107hm2。同時，在2020年，水稻種植面積極值約為4.05×106hm2，水稻單產(chǎn)極值約為7 030 kg/hm2。此外，糧食作物總產(chǎn)極值為6.39×107t，潛在增量約為1.14×107t。未來還可從優(yōu)化管理、提高水資源利用率、增加土壤肥力、基因工程培育高產(chǎn)作物等方面來提高黑龍江農(nóng)業(yè)生態(tài)系統(tǒng)總生產(chǎn)力。