1952年以来我国大豆单产变异特征及其影响因素研究*

秦婷婷,曹鑫悦,周泽群,褚超群,方雨桐,曲乐安,支俊俊**,王 震,耿 涛

(1.安徽师范大学地理与旅游学院 芜湖 241002; 2.中国矿业大学(北京)地球科学与测绘工程学院 北京 100083)

随着居民消费结构升级,大豆()需求快速增加,国内产需缺口不断扩大。2018年我国大豆85%依靠国际市场进口,主要来自于巴西、美国和阿根廷。一旦国际市场出现影响大豆供给的不利因素,国内大豆及关联产业将会受到巨大影响。2019年农业农村部印发大豆振兴计划实施方案,提出了提高大豆种植面积、提高大豆单产、提升品质以及绿色发展等大豆生产振兴目标。目前,我国大豆单产水平仅是世界大豆平均单产水平的70%,是大豆主产国美国的60%。因此,在当前我国耕地总量难以明显增加的现状下,提升大豆的单产水平成为提高我国大豆总产量的首要可行举措。目前国内外学者对制约大豆单产的影响因素研究多为田间尺度对照试验或长时间序列区域尺度下某单一因素对大豆单产的影响,如CO浓度、套作种植、气候变化等因素。研究结果表明,CO浓度升高可提高30.9%的大豆单株荚数和百粒重,与玉米()间作较单作具有显著的间作优势,最高温度和降水量对大豆单产有积极的显著影响,1971—2016年间霜冻对大豆产量的影响在一定程度上以负效应为主。

然而,目前有关我国大豆单产制约因素的研究缺少多个因素贡献率之间的定量化计算与比较,大豆单产与各影响因素之间的量化关系尚不明确,且鲜有在全国尺度上开展大豆单产影响因素及其空间异质性的研究。在此背景下,本文以全国各省市1952—2017年大豆产量相关的统计年鉴数据为基础,从大豆种植的管理措施、自然因素、科技水平、社会因素、经济因素等方面选取影响大豆单产的因子,基于增强回归树模型,量化各因子的相对重要性及其与大豆单产之间的关系,揭示1952年以来长时间序列下全国尺度和4个大豆主产区之间的大豆单产驱动力时空分异特征,以及改革开放前后大豆单产的驱动因子。研究结果对切实提高我国大豆单产水平、指导和调控农业种植结构、规划国家和地区社会经济发展,以及确保国家粮食安全等具有重要意义。

1 材料与方法

1.1 试验区概况

为揭示我国地域的大豆单产驱动因子空间差异特征,本文综合气候类型、耕作制度和品种类型等因素,将我国分为4个大豆主产区: 1)北方春大豆区,包括黑龙江、吉林、辽宁、内蒙古、新疆、甘肃、陕西北部(榆林、延安); 2)黄淮海流域夏大豆区,包括山东、山西、河南、河北、安徽北部(宿州、阜阳、亳州)、江苏北部(宿迁、徐州、淮安、连云港); 3)长江流域春夏大豆区,包括湖南、湖北、江西、浙江、四川、安徽中南部、江苏中南部、陕西南部(汉中、安康、商洛); 4)南方多熟大豆区,包括广东、广西、云南、贵州、福建、海南。

1.2 研究数据

1.2.1 数据来源

本文所使用的相关农业统计数据主要来源于全国各省市(不含港、澳、台地区)统计部门正式出版的统计年鉴、农村统计年鉴、农业年鉴、农业统计资料、农村经济年鉴和城市年鉴等,并利用“国家农业科学数据共享中心”平台获取的数据对缺失数据进行了补充完善。具体搜集了全国各省市(不含港、澳、台地区)1952—2017年(1952年、1965年、1978年、1990年、2000年、2010年、2017年共7个年份)农作物总播种面积、大豆种植面积、大豆产量、农用化肥施用量、农药施用量、有效灌溉面积、农用机械总动力、第一产业生产总值、地区生产总值、文盲率、城市化率、农作物受灾面积、研发经费支出等数据。气象数据来源于国家农业科学数据共享中心网站(http://data.cma.cn/)发布的各气象观测站点(共2421个站点)数据,包括年均温、年均降水量和年均日照时间,对这些数据使用ArcGIS软件空间插值统计分析,获得了1952—2017年全国各省市500 m × 500 m的气象要素(年均值)空间插值栅格数据。对数据进行单位统一换算、指标计算、涉及价格数据基准转换(如将所有GDP数据统计换算成可比价)、异常值删除等处理,得到共2008个样本数据。

1.2.2 指标选择

本文从大豆种植的管理措施、自然因素、科技水平、社会因素、经济因素等方面选取了13个驱动因子(表1)作为自变量,以大豆单产作为目标变量,构建增强回归树模型。

表1 用于模型构建的大豆单产影响因子Table 1 Selected influencing factors using for soybean yield modeling

1.3 研究方法

1.3.1 增强回归树模型

增强回归树模型(boosted regression trees,BRT)也称为促进回归树模型。BRT采用无放回抽样技术从训练样本中抽取一定比例(一般≥50%)的样本用于建模,分析自变量对因变量的影响程度,未被抽取的样本用于精度验证,最后对生成的全部分类回归树取均值并输出,从而提高模型的稳定性和预测精度。BRT的最大优点是不需要考虑自变量之间的交互作用,数据允许存在缺省值并支持多种数据类型,可量化各自变量因子对目标变量的影响程度(相对贡献率)以及与目标变量之间的非线性关系,清晰直观、易于解释。

采用开源编程软件R (version 3.6.3),基于统计整理得到的1952—2017年(1952年、1965年、1978年、1990年、2000年、2010年、2017年共7个年份)全国各省市共2008个样本点数据,将大豆单产作为因变量,13种驱动因子作为自变量,分别建立全时期(1952—2017)、改革开放前(1952年、1965年、1978年)和改革开放后(1990年、2000年、2010年、2017年) 3个模型,使用gbm软件包(version 2.1.8)建立增强回归树模型。

1.3.2 模型参数设置

模型中需要调试的参数有4个,分别是: baggingfraction (抽样率),指的是随机抽取训练样本的比例,用于构建单个决策树; learning rate (学习率),决定了每棵决策树对最终模型的贡献大小; tree complexity(树复杂性),决定了模型中预测因子交互作用的水平;ntree (树数量),指的是模型中决策树的数量。本研究将4个参数分别设置为0.60、0.005、3和20 000。最佳参数配置由基于10-折交叉验证(10-fold crossvalidation)的最高预测精度来确定,模型的预测精度通过重复运行50次计算取平均值。

1.3.3 模型性能指标

模型性能通过4个指标进行评估: 决定系数(coefficient of determination,)、相关系数(Pearson’s correlation coefficient,)、平均绝对误差(mean absolute prediction error,MAE)和均方根误差(root mean square error,RMSE)。其中,和的值越大、MAE和RMSE的值越小代表模型的性能越优越。具体计算公式如下:

2 结果与分析

2.1 1952—2017年大豆单产变异特征

图1展示了1952—2017年期间7个年份大豆单产的变异特征。该时期内大豆单产总体呈增长趋势,1952年大豆单产均值为0.816 t·hm,2017年增长至1.853 t·hm,是1952年的2.27倍,年均增长率为1.27%。值得注意的是,1965年的大豆单产均值为0.715 t·hm,相比1952年下降12.4%。从1965年开始,大豆单产进入快速增长期,至2000年增长至1.656 t·hm,是1965年的2.32倍; 该期间年均增长率为2.43%,约为1952—2017年期间年均增长率的两倍。2000年之后,大豆单产增长趋于平缓,至2017年期间的年均增长率降至0.67%,仅为1952—2017年期间年均增长率的一半。图1的箱形图展示了各年份大豆单产的地区变异,总体而言,各年份的大豆单产变异系数为34.1%～73.2%,属于中等程度变异,离散程度相对较高,表明全国各地市大豆单产之间存在较大的差异。

图1 1952—2017年中国大豆单产的变异特征Fig.1 Variations of soybean yield in different years from 1952 to 2017 in China

2.2 大豆单产模型的性能分析

全国尺度及4个大豆主产区的各模型性能结果如表2所示。全国尺度大豆单产驱动力BRT模型的性能结果为0.433,为0.639,RMSE为0.097,MAE为0.467,表明模型可有效解释全国尺度下43.3%的大豆单产变异性。最高解释度出现在黄淮海流域夏大豆区(67.8%),南方多熟大豆区的解释度也在60.0%以上,长江流域春夏大豆区的解释度则接近60.0%,而北方春大豆区的解释度最低,为43.1%。北方春大豆区跨度大,气候、土壤和种植结构等方面差异明显,区域内大豆单产影响机制差异较大,导致该模型性能低于其他地区。4个大豆主产区模型性能指标中的相关系数均高于0.63,MAE和RMSE均保持较低值,表明所构建的BRT模型性能良好。将全国尺度数据划分改革开放前后两个时期分别建模,改革开放前为0.486,改革开放后为0.494,模型性能较为稳定。

表2 用于大豆单产驱动力分析的增强回归树模型（boosted regression trees,BRT）的性能Table 2 Performance of boosted regression trees (BRT) models for soybean yield analysis

2.3 大豆单产驱动因子的相对重要性

由模型得到各驱动因子的相对影响,即各因子对大豆单产的贡献率(图2)。1952年至2017年,对大豆单产影响最大的因子为大豆播种面积占农作物总播种面积的百分比(SAP),相对重要性为20.9%。我国文盲率(ILR)从78.6%下降至3.0%,其对大豆单产的影响程度为18.9%,是促使大豆单产提升的重要因素。每公顷化肥(折纯)施用量(FCP)和每公顷农药(折纯)施用量(PCP)影响程度分别为10.6%和8.5%。年均降水量(AVP)的值变化幅度不大,有效灌溉面积占农作物播种面积中百分比(EIAP)增加明显,两因子影响程度分别为7.6%与7.4%。2017年全国平均每公顷农业机械总动力(AMP)为5.9 kW·hm,是1952年0.0013 kW·hm的4538倍,但其影响程度仅为6.1%。年均温(AVT)因子紧随其后,影响程度为6.0%。其他因子影响程度较低,依次为研发支出所占地区生产总值的百分比(RDIG,4.2%)、年均日照时间(AST,3.8%)、人口城镇化率(PUR,2.8%)、第一产业占地区生产总值的百分比(PIG,2.0%)、受灾害面积占农作物总播种面积的百分比(DAP,1.3%)。

图2 1952—2017年（A）、改革开放前（1952—1978年,B）和改革开放后（1985—2017年,C）大豆单产各驱动因子的相对重要性（因子简写见表1,图中误差线表示模型运行50次计算求得的各因子重要性的标准差）Fig.2 Relative influence of each driving factor on soybean yield during the periods from 1952 to 2017(A),before the reform and opening up (B,1952—1978),and after the reform and opening up(C,1985—2017) (See Table 1 for factors abbreviations,error bars represent standard deviations of the variable importance averaged over 50 model runs)

大豆单产主要贡献因子,改革开放前为大豆播种面积占农作物总播种面积的百分比(SAP,23.1%)、文盲率(ILR,16.2%)和每公顷化肥(折纯)施用量(FCP,9.4%),改革开放后为大豆播种面积占农作物总播种面积的百分比(SAP,20.9%)、每公顷农业机械总动力(AMP,14.5%)和年均温(AVT,12.4%)。可以看出,无论改革前、改革后,还是全时期,SAP都是各个时期对大豆单产影响最大的因子,其重要性均超过20%。

2.4 各驱动因子与大豆单产之间的非线性交互关系

各驱动因子与大豆单产之间的关系均呈非线性,大豆单产与各驱动因子之间的定量关系如局部依赖图(图3)所示。1952—2017年,全国大豆播种面积占农作物总种植面积的百分比(SAP)由8.3%下降至6.0%,大豆单产随大豆SAP的升高而下降直至趋于平稳。文盲率(ILR)高于20%时对大豆单产影响微弱,但在0～20%范围内对随着ILR升高,大豆单产显著降低。随着每公顷化肥(折纯)施用量(FCP)提升,大豆单产明显提升,当施用量超0.5 t·hm后,大豆单产基本保持不变; 每公顷农药(折纯)施用量(PCP)因子、每公顷农业机械总动力(AMP)因子影响机制与其类似。

图3 各驱动因子与大豆单产之间的非线性交互关系（指标全称见表1,因子名后括号内分别为因子的单位及相对重要性）Fig.3 Interactive nonlinear relationships between soybean yield and its impact driving factors (See Table 1 for factors abbreviations,the unit and relative importance of each factor were showed in the bracket after factor’s names).

大豆单产随着年均降水量(AVP)的增加波动变化,年均降水量在500～1000 mm区间时,年均降水量与单产呈正相关; 大豆单产与有效灌溉面积占农作物总播种面积的百分比(EIAP)因子整体呈正相关;年均温(AVT)在0～5 ℃时,大豆单产最高,随着温度上升大豆单产总体呈下降趋势; 日照时间(AST)的增加亦不能直接提高大豆单产。研发支出占地区生产总值(RDIG)的增加对大豆单产提升具有积极影响,对大豆单产的贡献率为4.2%。城镇化(PUR)的提升对大豆单产总体呈负面影响。第一产业占地区总产值的百分比(PIG)对大豆单产不仅贡献率较低,对大豆单产提升也较为微弱。

2.5大豆各主产区的核心驱动因子

大豆各主产区主要驱动因子不同,但文盲率在各主产区均为主要驱动因子之一,主要驱动因子空间分异如图4所示。影响北方春大豆区大豆单产的最重要因素为每公顷农业机械总动力(AMP,13.1%)、文盲率(ILR,11.8%)、大豆播种面积占农作物总种植面积的百分比(SAP,11.7%),影响黄淮海流域夏大豆区大豆单产的最重要因素为每公顷化肥施用量(FCP,25.6%)、每公顷农药(折纯)施用量(PCP,18.4%)、文盲率(ILR,15.1%),影响长江流域春夏大豆区大豆单产的最重要因素依次为研发支出占地区生产总值的百分比(RDIG,21.5%)、有效灌溉面积占农作物播种面积的百分比(EIAP,14.3%)、文盲率(ILR,11.8%),影响南方多熟大豆区大豆单产的最重要因素依次为文盲率(ILR,35.4 %)、每公顷化肥(折纯)施用量(FCP,22.9%)、每公顷农药(折纯)施用量(PCP,7.8%)。

图4 1952—2017年中国大豆主产区的大豆单产主要驱动因子（指标全称见表1）Fig.4 Top three determinants of soybean yield in the four major soybean producing areas of China during the period from 1952 to 2017 (See Table 1 for factor abbreviations)

3 讨论

3.1 全国尺度下大豆单产的影响因子

研究结果表明,全时期、改革开放前、改革开放后3个时期大豆单产的最主要影响因子均为大豆播种面积占农作物总播种面积的百分比。改革开放前大豆单产的其他重要影响因子还包括文盲率与每公顷化肥(折纯)施用量,改革开放后还包括每公顷农业机械总动力与年均温。文盲率降低,农民知识水平提高,有助于提高单产。1949年我国文盲率高达80%,改革开放以后,国家出台一系列政策来支持农村义务教育事业发展,2000年我国已经实现了“基本普及九年义务教育、基本扫除青壮年文盲”的历史性任务,2017年我国文盲率仅3%。义务教育的普及提高了农村劳动力的知识水平,职业教育培养了一批农村实用型人才,显著降低了文盲率,使得农民在优质品种选择、作物生长管理以及农业科技成果应用等方面的知识水平大幅提升,并进一步提升农业生产效率和农业技术水平,从而提高大豆单产。

3.2 影响因子的相对重要性及非线性关系

扩大大豆种植面积、提高大豆在农作物播种总面积中的占比能够有效提高大豆单产,其机理在于提高大豆在当地农作物生产中的地位能够促使农民提高耕作管理水平(如减少重迎茬比例、增加轮作),增加农业生产资料的投入(如优良品种、化肥和农药施用量、机械化作业等),以及吸引政府加大耕作补贴、农业保险、贷款等方面的扶持力度,从而降低生产成本和提高自然灾害风险的应对水平,进一步提高大豆单产。因此,在提高大豆种植面积的同时,需增加其他生产要素的投入力度,以进一步形成规模化效益。本研究结果中大豆单产与大豆播种面积占农作物总种植面积的百分比之间的交互作用关系反映的主要是大豆播种面积占比在时间尺度上的变化,与1952—2017年期间大豆播种面积由1.1679万hm下降至0.8245万hm,但大豆单产从0.816 t·hm公顷增长至1.853 t·hm的趋势一致。为进一步揭示同一年份不同地区之间大豆单产存在差异性的关键影响因子,在今后研究中可考虑搜集县级统计资料以显著扩大样本数量,从而构建单个年份的大豆单产与各指标之间的相互作用关系来进行分析。

从大豆生长情况看,多年使用有机肥料土壤中钾、磷、氮等物质的含量已经接近饱和,仅增加氮肥、钾肥等有机型的肥料施用量而不提高其他微量元素如锌、锰等的投入,则很难对大豆生长起促进作用。农药对防病治虫、促进粮食和农业稳产高产至关重要,但区域中也存在滥用现象。因此,化肥与农药须合理施用,否则容易造成土壤养分失调、地力下降,导致大豆抗逆性差、产量降低。大豆单产与生育期的水分、光照、热量等条件关系密切,大豆单产可能受积温和种植制度等的影响较大,因此,有必要进一步分区进行统计分析。

3.3 大豆主产区驱动因子空间分异特征

由于大豆四大主产区所处地理环境、种植规模、经济发展水平等条件不同,大豆单产的核心驱动力存在空间差异。北方春大豆区,大豆种植规模大,地形辽阔平整,区域内各省份大豆全程机械化水平略有不同。其中,东北春大豆区中集中分布的松花江和辽河流域平原地带基本实现了大豆全程机械化生产,在耕整地、播种、收获等各个环节均有适合该地区生产的成熟机具,大豆生产综合机械化水平高。黄土高原春大豆种植区域气候寒凉,土质瘠薄,品种中耐瘠薄、耐旱的黑豆较多。西北春大豆种植区域年均降水量少,土壤蒸发量大,种植大豆须人工灌溉。

黄淮海流域夏大豆区是我国大豆的第二大主产区,是高蛋白大豆产区。每公顷化肥施用量是该地区的最重要驱动因子,土壤养分供应、氮磷钾养分吸收规律、田间施肥技术、肥料损失等的区域异质性高,不同区域应根据当地的气候状况、土壤特征、种植制度和大豆肥效等因素,制定专用复合(混)肥料配方。黄淮海流域夏大豆产区在耕种环节仍然存在上茬作物秸秆处理难的问题,大豆收获主要依靠稻麦联合收获机,机具适应性差、收获损失大、破碎程度高。

四大产区中长江流域最注重科技投入,1952—2017年期间研发支出占地区生产总值的百分比平均为0.8%,大豆单产的提高需要加大科技投入,实施大豆重大科研联合攻关,培育出适合地区种植的优质大豆种苗,并加快大豆产业领域新技术、新成果的应用与转化。

与北方及黄淮海大豆产区相比,南方基本以种植水稻(Oryza)为主,旱地则有玉米()、甘薯()、油菜()等作物,耕作制度比较复杂。大豆种植方式多样,如与旱地作物实行间套混作或利用田埂隙地栽培等。南方大豆产区机械化生产受品种、地形、种植模式等多方面因素制约,大豆机械化收获推进缓慢,机械化生产水平相对较低,尤其是收获环节,缺乏专用的大豆收获机具。可加大大豆耕种收获等环节机具购置补贴力度,推广适用的大豆全程机械化技术及装备,提升大豆生产科技水平。第一产业占地区生产总值百分比和年均降水量是该地区除化肥外的主要影响因子,南方多熟大豆产区雨量充沛,为大豆生长提供了丰富的水分,1952—2017年间年均降水量为1599.4 mm,是北方春大豆产区的3.5倍。

3.4 增强回归树模型的特点及结果不确定性

与前人研究通常使用的线性回归模型和决策树模型等相比,本文所构建的增强回归树模型可量化各驱动因子对大豆单产的影响(相对重要性)以及大豆单产与单个驱动因子之间的非线性关系(局部依赖图),其拟合多个树的过程能最大程度地弥补单一树模型预测能力弱的缺点,比大多传统模型方法有更强的预测能力,采用的无放回抽样技术可提升模型稳定性与模拟精度,并能有效评估模型结果的不确定性(由经过多次模型重复运行计算结果的标准差表达)。

本文结果的不确定性主要来源于数据的搜集质量,由于统计数据时间跨度大、覆盖范围广,早期的统计数据与实际数据可能存在一定偏差; 其次,尽管本文对全部数据进行了统一基准等处理,但各时期统计口径及标准等方面存在的差异可能增大模型结果的不确定性; 再次,由于本文使用的统计数据是以地级市为单位进行的汇总数据,模型拟合结果只能细化表达至地级市之间的空间差异特征。

4 结论

本研究构建的增强回归树模型不仅量化了各影响因子对于大豆单产的相对重要性,并且有效揭示了大豆单产与各影响因子之间的非线性关系。全国尺度下,各时期大豆单产的最重要影响因子均为大豆播种面积占农作物总播种面积的百分比,对大豆单产有重要影响的因子在改革开放前还包括文盲率和每公顷化肥(折纯)施用量,在改革开放后还包括每公顷农业机械总动力和年均温。不同大豆主产区的大豆单产驱动力存在空间异质性,但文盲率在各主产区均为主要驱动因子之一。研究结果可为国家及各省市今后提升大豆单产水平所应采取的有效措施以及调整农业种植结构等提供决策依据,各省市需合理施肥和播撒农药,努力提高机械化水平和农业生产者的知识水平。