0.   引　言

准确估测农作物产量对于保障粮食安全至关重要。现有估产方法可以归纳为传统估产和基于遥感手段估产两类。其中，传统估产方法，如抽样调查和实地观测等，虽然估测精度高，但是估测过程中需要耗费大量人力和物力，估测成本相对较高。此外，传统方法只能对少数样本进行实测，无法满足大范围的估产工作^[1]。随着遥感技术的不断发展，遥感手段逐渐成为农作物估产的主流方法^[2]。目前，基于遥感技术的农作物估产方法可归纳为统计回归法、作物生长模型法、生物量和收获指数法、机器学习及深度学习法等。其中统计回归法原理简单、效率较高，但缺乏作物生长机理的融合^[3]；作物生长模型法虽然与作物生长机理紧密结合，但需要大量大田作物参数及各种影响因子对模型进行修正^[4]；机器学习及深度学习法估产精度较高，但其“黑箱”模型的特点使其在机理方面难以解释^[5]；生物量和收获指数法既可以很好地与作物生长机理相结合，也具备效率高、成本低的特点，在作物估产研究中被广泛应用^[6]。

在生物量和收获指数法中，陆地生态系统总初级生产力（gross primary production，GPP）和陆地生态系统净初级生产力（net primary production，NPP）具有关键作用。近些年来，采用以遥感方法为主的光能利用率模型，如卡内基-艾姆斯-斯坦福方法（Carnegie Ames Stanford approach，CASA）^[7]、净光合生产力模型（photosynthesis， PSN）^[8]和植被光合作用模型 （vegetation photosynthesis model，VPM）^[9]等，进行GPP和NPP估算，成为陆地生态系统生产力模拟研究的一个重要发展方向。其中，相较于CASA和PSN模型，VPM模型结合了增强植被指数（enhanced vegetation index，EVI）和陆表水指数（land surface water index，LSWI）^[9]，在农田C3和C4作物GPP模拟中应用更加广泛^[10]，并表现出误差低、估产精度高^[11]的特点。

基于传统VPM模型的农作物估产工作中，需要确定最大光能利用率、收获指数和植物碳素含量等关键参数。各项关键参数因作物类型、作物品种和研究地区的差异而有所不同，如玉米的最大光能利用率在美国为3.31 g/MJ（以C计，下同），在中国范围内则介于2.17～2.29 g/MJ之间^[12]；中国不同地区的收获指数波动范围为其平均值的20%左右^[13]。因此，现有研究多通过实地测量的方式获取这些关键参数，以实现准确的农作物估产^[14]。然而，这种数据获取方式成本高、效率低，与利用遥感手段进行低成本、高效率的农作物估产目标相违背，严重制约了基于遥感技术的VPM等模型在地块尺度上快速、准确实现估产工作的能力^[15]。

为解决大范围估产参数确定难、成本高的问题，本研究在传统VPM估产模型基础上，提出了VPM与转产系数相结合的估产模型；并以吉林省玉米为研究对象，基于Sentinel-2遥感数据完成了10m分辨率的玉米估产研究。

1.   研究方法与数据来源

1.1   研究区概况

本研究选取中国吉林省内玉米年播种总面积超过4×10⁴ hm²的30个区县作为研究区（图1）。吉林省（40°52′～46°18′N，121°38′～131°19′E）是中国重要的粮食生产基地，2021年全省玉米产量达3.20×10⁷ t，位列全国第二。吉林省地形表现为东南高、西北低的趋势，并以东北–西南向纵贯全省的大黑山山脉为界，分为东部长白山地和西部松辽平原两大区域。吉林省属于温带大陆性气候，冬季长而严寒，夏季短而炎热，年平均降水量500～600 mm，年均温在−2～6 ℃之间。研究区内粮食作物主要以玉米为主，据吉林省统计局统计，2016—2021年间研究区的年平均玉米播种面积和产量分别为3.88×10⁶ hm²和2.81×10⁷ t，分别占研究区粮食作物总播种面积和总产量的77.32%和80.01%。

图  1  研究区概况

Figure  1.  Study area

下载: 全尺寸图片 幻灯片

1.2   数据来源

本研究所使用的数据包括遥感数据、气象数据、玉米种植分布数据以及农业统计数据。遥感数据包括欧洲航天局的Sentinel-2数据和美国宇航局的MOD09A1数据，其中Sentinel-2数据空间分辨率为10 m，时间分辨率为5 d；MOD09A1数据空间分辨率为500 m，时间分辨率为8 d。气象数据主要包括光合有效辐射（photosynthetically active radiation，PAR）和气温，其中PAR数据来源于全球陆表特征参量数据产品数据集 （global land surface satellite，GLASS），空间分辨率为5 km，时间分辨率为1 d；气温数据采用国家青藏高原数据中心（https://data.tpdc.ac.cn/）发布的“中国1 km分辨率逐月平均气温数据集”^[16]。玉米种植分布数据来源于国家生态科学数据中心（http://www.nesdc.org.cn/）发布的“2016—2021年中国玉米种植分布数据集（crop distribution，CD）”，空间分辨率为30 m，时间分辨率为1 a^[17]。此外，为探讨玉米种植分布数据的空间分辨率对研究结果的影响，本研究选择该平台发布的“全国三大粮食作物1 km种植分布数据集”中玉米种植分布数据用于对比研究^[18]。2016—2021年研究区内30个区县的玉米总产量和播种面积来源于吉林省统计局^[19]（http://tjj.jl.gov.cn/）。最后，本研究将不同空间分辨率的栅格类型数据重采样为10 m空间分辨率数据，进而展开后续估产工作。

1.3   研究方法

1.3.1   VPM模型基本原理

VPM模型作为光能利用率模型^[20]的典型代表之一，被广泛应用于全球植被生产力的模拟^[21]。VPM模型中将植被冠层吸收光合有效辐射的比例分为叶绿素吸收部分和非光合植被吸收部分。由于估产所需变量为NPP，而传统VPM模型模拟对象为GPP，因此结合GPP与NPP的转换系数进行调整。主要算法如下：

式中${V}_{\mathrm{N}\mathrm{P}\mathrm{P}}$为净初级生产力；${V}_{\mathrm{G}\mathrm{P}\mathrm{P}}$为总初级生产力；${V}_{\mathrm{C}\mathrm{U}\mathrm{E}}$表示碳利用率，即GPP转换成NPP的转化系数；${\varepsilon }_{g}$表示光能利用率，g/MJ；${V}_{{FPAR}_{chl}}$为叶绿素吸收部分的植被冠层吸收光合有效辐射比例；$Q_{PAR}$表示光合有效辐射，MJ/m²；${\varepsilon }_{0}$表示最大光能利用率，g/MJ； ${T}_{scalar}$、${P}_{scalar}$和${W}_{scalar}$分别表示温度、水分和叶龄在冠层尺度上对最大光能利用率的胁迫系数^[22]。

进一步，NPP转化为单产的计算公式如下^[23]：

式中$Y$为农作物单产，t/hm²；$b$为作物地下部分生物量较之整株植物所占的比例；$V_{\mathrm{HI}}$ 为收获指数；$c$为植物碳素含量；$\omega$为作物收获后储藏期的含水率系数。依据前人研究^[10]，$b$、$V_{\mathrm{HI}}$ 、$c$、$\omega$分别取值为0.1、0.49、0.45、14%。

1.3.2   转产系数的构建与计算

根据式（1）～（4），传统VPM模型中农作物单产的计算公式可以综合表达为

式中${\varepsilon }_{0}$、$V_{\mathrm{HI}}$、${V}_{\mathrm{C}\mathrm{U}\mathrm{E}}$ 、$c$、$b$ 和$\omega$等关键参数通常为固定值，忽略了不同作物长势、环境、气候条件、作物品种和农业管理条件下参数的差异性，进而导致模拟结果存在较大偏差；$Q_{PAR}$、${V}_{{FPAR}_{chl}}$、${T}_{scalar}$、${P}_{scalar}$和${W}_{scalar}$为动态观测值。为解决关键参数值获取难、成本高等问题，本研究提出VPM动态观测指数${\mathrm{A}\mathrm{P}\mathrm{A}\mathrm{R}}_{{\varepsilon }_{g}}$和转产系数$a$以代替多个区域化固定参数和动态观测变量，以对原模型参数进行整合。

本文所定义的VPM动态观测指数$\mathrm{A}\mathrm{P}\mathrm{A}\mathrm{R}_{\varepsilon_g}$（式中用以$I_{{\mathrm{A}\mathrm{P}\mathrm{A}\mathrm{R}}_{{\varepsilon }_{g}}}$表示）计算如下：

本文依照VPM原理提出将众多相对固定的关键参数综合成为一个固定的转产系数a。结合以后可降低估产成本，提高估产精度。如果具备更为精细和丰富产量数据则可以得到更小范围的转产系数，以代替原模型中相对固定的众多关键参数。本研究所定义的转产系数$a$计算公式如下：

转产系数$a$受研究区环境、作物类型、作物品种和管理条件等影响，且需依托大量实测数据，获取不易。本研究通过作物统计单产数据和${\mathrm{A}\mathrm{P}\mathrm{A}\mathrm{R}}_{{\varepsilon }_{g}}$之间的线性关系计算得到。计算如下：

最后结合各地级市不同区县的${\mathrm{A}\mathrm{P}\mathrm{A}\mathrm{R}}_{{\varepsilon }_{g}}$与对应区县的统计单产建立截距为0的线性回归模型，进而模拟得到各地级市的转产系数$a$。

1.4   模型验证

本研究将2016—2021年间研究区各区县${\mathrm{A}\mathrm{P}\mathrm{A}\mathrm{R}}_{{\varepsilon }_{g}}$与其对应的统计单产数据集以7∶3的比例随机分为训练集和测试集。首先通过训练集得到地级市的转产系数$a$，计算出研究区各区县的估产结果后，进而在验证集上获取最终的模型估产精度结果。本研究所采用的精度评价指标包括决定系数^[24]（coefficient of determination，R²）、均方根误差（root mean square error，RMSE）、平均相对误差（mean relative error，MRE）和归一化均方根误差^[25]（normalized root mean squared error，NRMSE）。

2.   结果与分析

2.1   精度验证

相比传统VPM模型，本研究方法估产精度相对更高。对比结果显示（图2a、2b），本研究模型估测单产和统计单产拟合程度更高，R²（0.53）高于传统VPM模型（0.35），提升0.18；估产结果误差（RMSE=0.81 t/hm²、MRE=9.40%、NRMSE =11.73%）也较传统VPM估产模型（RMSE=1.03 t/hm²、MRE=13.19%、NRMSE =14.97%）要更小，RMSE降低3.24%。

图  2  模型精度对比

注：R²代表决定系数，RMSE代表均方根误差，MRE代表平均相对误差，NRMSE代表归一化均方根误差。

Figure  2.  Comparison of model accuracy

Note: R² stands for coefficient of determination, RMSE stands for root mean square error, MRE stands for mean relative error, NRMSE stands for normalized root mean square error.

下载: 全尺寸图片 幻灯片

在空间上（图3a、3b），本研究区内有12个县区（40%）的RMSE有所下降。主要原因是由于本研究模型可以得到更为本地化的转产系数，并且转产系数比原固定的关键参数更加接近于真实值，使估产结果因地区品种、气候等因素的差异引起的误差有所减小。

图  3  模型均方根误差RMSE空间分布对比

Figure  3.  Comparison of model RMSE spatial distribution

下载: 全尺寸图片 幻灯片

2.2   吉林省玉米单产估测

基于VPM模型原理与转产系数，本研究估测了2016—2021年吉林省30个区县玉米单产（图4）。研究区内各地级市转产系数如表1所示，其中四平市转产系数最高，通化和白城市的转产系数则较低。估产结果表明，研究区的高产区主要集中在吉林省中部，并呈现由中部向西北和东南地区递减的趋势。这主要受到地形、气候因素的影响，如吉林省东南地区多为山地，玉米单产偏低；而吉林省西部和中部多为平原地区，并且河流分布广泛，有利于作物生长，玉米单产偏高。

图  4  2016—2021年研究区玉米估测单产空间分布

Figure  4.  Spatial distribution of estimated maize yield from 2016 to 2021 in the study area

下载: 全尺寸图片 幻灯片

表  1  吉林省各地级市转产系数α

Table  1.  Yield conversion coefficient α of each city in Jilin Province

城市 City	转产系数 Yield conversion coefficient/(×10−2)		城市 City	转产系数 Yield conversion coefficient/(×10−2)
长春	1.21		四平	1.33
通化	1.08		吉林	1.10
松原	1.19		白城	1.05

下载: 导出CSV 

| 显示表格

从时间上看（图5），研究区玉米单产高值主要出现在2016年和2017年，平均单产分别达到8.60和9.02 t/hm²；2018—2021年的平均单产则为6.62 t/hm²。与统计结果对比显示，估测单产结果与统计单产在数值分布以及年际变化趋势上基本保持一致，但个别年份略有偏差。如统计单产最低值（6.40 t/hm²）出现在2018年，而估测单产最低值出现在2020年（6.15 t/hm²）。出现该问题的原因可能是未考虑农业技术进步所带来的影响，例如，新的种植品种、育种技术或农业机械的引入可能会在某些年份提高产量，而在其他年份可能没有同样的效果；偶然因素，如在农业生产中，偶然因素（如意外事件、天灾等）可能会导致个别年份的产量偏差，这些因素可能无法完全通过统计或估测方法考虑到。

图  5  2016—2021年研究区玉米估测单产与统计单产对比

Figure  5.  Comparison of estimated and statistical maize yield from 2016 to 2021 in the study area

下载: 全尺寸图片 幻灯片

2.3   敏感性分析

本研究通过对比调整转产系数的精细程度（市级与省级转产系数）、玉米种植分布的空间分辨率（30 m与1 km）以及遥感数据的空间分辨率（10 m与500 m），对本研究中估产模型进行敏感性分析（表2）。

表  2  敏感性分析表

Table  2.  Sensitivity analysis table t·hm⁻²

项目Item		转产系数α		1 km CD	MOD09A1
		市级 City-level	省级 Province-level
城市 City	白城	1.29	1.36	1.24	2.13
	四平	0.94	1.30	1.34	1.74
	长春	1.07	1.13	1.08	1.11
	松原	0.93	0.96	1.03	0.99
	通化	0.71	0.84	1.09	0.79
	吉林	0.59	0.71	0.96	0.59
年份 Year	2016	1.38	1.39	1.61	1.84
	2017	0.91	0.98	0.97	1.38
	2018	0.80	1.04	1.13	1.27
	2019	0.67	0.86	1.01	1.22
	2020	0.87	0.94	0.98	0.94
	2021	0.77	0.94	0.98	1.31
注：从左到右4列数据分别是使用不同尺度（市级、省级）的转产系数、不同空间分辨率的玉米分布数据（1 km）和遥感数据（500 m）估产结果的RMSE值。
Note:Four columns from left to right were the RMSE of estimated maize yield by using different scale conversion coefficients (city level, provincial level), maize distribution data (1 km), and remote sensing data (500 m) with different spatial resolutions, respectively.

下载: 导出CSV 

| 显示表格

具体来说，在时间上和空间上与使用地级市级转产系数、省级转产系数、1 km分辨率玉米种植分布数据和500 m分辨率MOD09A1遥感数据的估产结果进行精度对比。结果显示，结合地级市级转产系数、10 m空间分辨率的遥感数据和30 m空间分辨率作物种植分布数据，模拟效果更好。具体而言，在市级层面上面，使用MOD09A1遥感数据的估产结果中，个别地级市的RMSE出现高值，并且在逐年的对比中也总是表现出较大的RMSE；使用省级转产系数和使用1 km玉米种植分布数据的估产结果中，RMSE在时间和空间上基本相同，但从整体上看，前者的RMSE略低。综合来看，各数据对估产结果影响从大到小依次是遥感数据空间分辨率、作物分布数据空间分辨率和所使用的转产系数精细程度。

3.   讨　论

为解决传统VPM模型中区域化相对固定的参数多且较难获取的问题，本研究基于传统VPM模型估产原理，分别提出转产系数 $a$ 和VPM动态观测指数$\mathrm{A}\mathrm{P}\mathrm{A}\mathrm{R}_{\varepsilon_g}$以代替多个区域化固定参数以及多个动态观测变量，简易化了获取区域性的关键参数过程。本研究模型通过$\mathrm{A}\mathrm{P}\mathrm{A}\mathrm{R}_{\varepsilon_g}$和单产的线性关系得到各地市的转产系数，以替代传统的一个省一个参数的问题。研究发现，各地级市的转产系数存在较大差异（表1），并且能够有效减少省级尺度转产系数带来的估产误差（图5）。与传统VPM模型模拟结果对比发现，本研究模型R²提升0.23，RMSE、MRE和NRMSE分别下降了0.23 t/hm²、3.66%和3.53%，模拟精度更高。此外，与其他研究方法相比，本研究模型也具有一定的优势。如王鹏新等^[26]使条件植被温度指数和叶面积指数结合机器学习进行华北地区玉米估产时，R²为0.30；周西嘉等^[27]使用Sentinel-2数据进行LAI重建进而得到冬小麦的估产结果，其R²为0.52。除模拟精度外，本研究基于Sentinel-2遥感数据实现了10 m空间分辨率的玉米估产，在空间分辨率上具有一定的优势^[28]。本研究所使用数据的开源性相较于部分使用无人机和非公开数据集具有成本优势和容易获取的优势^[29]；最后，本文模型的估产机理性相较于深度学习等人工智能模型相比，具有较强的可解释性等特点^[30]。

本研究也存在一定的局限性和不确定性。首先，从数据角度来看，1）光学遥感数据、气象数据作为该模型的重要输入变量，受时间分辨率、多数据空间分辨率不一致和云等因素的影响较大，可能导致模拟精度降低；2）本模型使用的作物种植分布范围数据存在误差因素，限制了其应用范围；3）地块级作物实际单产和统计数据缺乏导致本研究未能实现获取更为精细尺度的县级转产系数；4）本文模型采用气象数据和遥感数据对研究区的玉米进行估产，但未考虑土壤、地形、管理措施等多种因素影响。其次，从方法和试验过程来看，由于数据和GEE平台计算限制，未采取更高时空分辨率数据和时空融合方法估测值。此外，本研究模型对吉林省30个县区玉米单产估产效果较好，但对其他地区或多作物估产研究的适用性并不清楚。因此，可以从提升输入数据分辨率、降低作物分布数据误差、增补地块级实际单产数据和增加更多的估产因子等方面着手改进方案，以提高估产模型的数据丰富度和作物适用性，提高估产结果精度。

4.   结　论

本研究基于VPM模型估产原理，整合原模型的各动态变量形成动态观测指数${\mathrm{A}\mathrm{P}\mathrm{A}\mathrm{R}}_{{\varepsilon }_{g}}$，将原模型中的固定关键参数整合为转产系数，提出了一种新的思路来解决传统模型中固定参数多且难以获取的问题，对大规模、低成本的精确作物估产工作具有一定的现实意义。

1）相较于传统综合植被光合作用模型（vegetation photosynthesis model，VPM），本研究提出的结合转产系数的估产模型估产精度更高（改进前和改进后模型决定系数分别为0.35和0.53，均方根误差分别为1.03和0.81 t/hm²、平均相对误差分别为13.19%和9.40%，归一化均方根误差分别为14.97%和11.73%。改进后的模型，既保持了良好的机理性，又可以在缺少实测参数的情况下，基于县级统计数据完成低成本、大面积的地块尺度作物估产。

2）本研究结合Sentinel-2遥感数据，得到2016—2021年吉林省10 m分辨率玉米单产分布图，显示中部高产（7～13 t/hm²），东西部较低（5～9 t/hm²）。年际变化显示，2016—2017年单产较高，而2018—2021年偏低。

3）使用精细转产系数、高分辨率遥感及玉米种植分布数据能提升估产精度。模型精度的敏感性程度依次为：遥感数据空间分辨率、作物分布数据空间分辨率、转产系数精细度。