引入地形特征的田块尺度玉米遥感估产与空间格局分析

叶 强, 杨凤海*, 刘焕军,2, 孟祥添, 官海翔, 崔 杨

(1.东北农业大学公共管理与法学院，哈尔滨 150030；2.中国科学院东北地理与农业生态研究所，长春 130012)

玉米是中国三大主要粮食作物之一，种植面积及年产量达到中国粮食总产量的33.6%和36.1%。玉米作物长势动态监测及产量的准确预测，为规模化农业经营管理提供科技支撑，是精准农业发展的迫切需求，对国家粮食宏观价格调控，三农振兴以及粮食外贸等具有重要意义。

玉米产量的预测方法多采用传统估产法、农学估产法及气象统计法等。其中，农学预报法主要根据产量器官参数如有效穗数、百粒重、平均穗粒数等估算产量，受作物种类、品种、耕作制度等限制[1]；气象统计法以天气要素预测产量，发展较为成熟，但存在单一气象因子估产精度不足且各要素空间插值方法不同，难以划定气象站数据范围的问题[2]；传统估产法中地面采样工作采集样本数据规模较小、时间周期长、气象条件限制，样本实验室理化性质检测结果存在滞后性的问题，难以满足及时、准确的作物估产和田间作物生态保护的管理需求[3]。遥感技术突破地面采样的局限，以快速、无损、获取成本低、多尺度、多时间序列的优点成为精准农业研究热点[4]。Chen等[5]采用作物物候和叶面积指数引入MCWLA-Wheat同化模型提升相对估产精度16%。朱婉雪等[6]结合少量产量样本点与植被指数建立线性模型，决定系数R2最高达到0.7；Muoz等[7]应用产量和遥感参量建立指数模型和Richards非线性回归模型，偏差信息准则(deviance information criterion，DIC)最低为62.8。但是，对田块尺度作物遥感估产，测产多采用对角线、十字、“S”形，同级重复采样以保证单产精度。测产过程耗时耗力，破坏采样区间作物样本，存在主观因素影响，不同规模测产工作造成大量粮食损耗。

谷物联合收割机搭载测产系统可以实现同步收获、同步测产，从而节省人工，提高测产效率，降低作物粮食损耗量。该系统装载有差分校正的全球定位系统(global positioning system,GPS)能够收集有空间地理信息的产量数据，测产精度可达97%以上[8]。为提高遥感估产精度，以遥感数据作为主要信息源的统计模型不断优化，结合作物生长模型估产，或加入遥感辅助信息，如土地分类图、作物氮含量、作物株高[9]、叶面积指数、雷达数据、倾斜测量构建3D点云等[10]，以提高统计经验模型的可靠性和可延展性。但作物生理参数以及土壤化学性质随作物生长周期变化，多时期采集样本费时耗力。高精度地形数据空间分辨率高，地形演变在一定时期内稳定性强，影响作物水肥运移，间接作用于土壤理化性质和植被生长，造成作物产量空间异质性。地形因子作为辅助信息引入遥感估产模型具有合理性。大尺度区域研究中高程影响稻田空间格局和生产资源，如光照、土壤类型等[11]；李宇宇等[12]研究了赤峰草场发现低湿区牧草鲜重分别高于坡地和平地67.94%，45.7%；土壤属性与甘蔗产量在凹地较高，在平坦地形处空间变异性较弱[13]；坡度约为2°区域大豆长势最佳[14]；自然降雨条件下坡耕地坡度为15°设立高反坡阶玉米产量最高[15]。目前，田块尺度结合地形因子提高变量测产数据遥感反演估产精度的研究尚未完善。

以黑龙江省北安市规模化种植的春玉米为研究对象，利用作物关键生长期SPOT 6影像，提取7种典型植被指数，结合变量测产数据建立经验统计模型，确定适合春玉米遥感估产的最佳生长时期及最优植被指数；基于无人机获取的地形数据提取6种地形因子，应用多元逐步回归评价引入地形因子的春玉米遥感估产模型；应用聚类分析和空间统计分析探索产量空间分布格局。研究旨在筛选春玉米监测最优植被指数和最佳时期，建立以遥感信息结合地形因子的估产模型，确定玉米最佳生长区域，为农场精细生产管理、耕地高效利用、玉米单产精细尺度制图提供科学依据。

1 材料与方法

1.1 研究区概况

研究区位于黑龙江省北安市，地处小兴安岭南麓边缘，属于典型黑土区垄作耕地，年平均气温0.5 ℃，年平均降雨量为570 mm。变量测产数据采集区域位于该市赵光农场，农场经营管理条件一致，具有相同的机械化耕作方式与种植制度。主要种植作物为春玉米，春播晚熟，一年一熟。播种时间为每年4月底，成熟收获时间为同年9月底10月初。实验区田块面积5 hm2，中心经纬度为48°2′42.45″ N、126°45′19.45″ E，平均海拔为307.41 m，平均垄长为553.09 m，最大高程差为9.37 m，采用玉米大垄双行种植模式，垄宽为1.1 m。

1.2 数据获取与处理

编程订购2018年7月27日、8月19日、9月13日的SPOT 6多光谱卫星影像，分别为玉米生命周期中吐丝期(S1)、灌浆期(S2)和成熟期(S3)。卫星影像空间分辨率为6 m，包含蓝、绿、红和近红4个波段。使用ENVI 5.1对影像进行辐射定标、大气校正、正射校正等预处理。高精度地形数据获取，采用DJI Phantom 4 Pro RTK无人机飞行平台，搭载1 inch(1 inch=2.54 cm)2 000万像素CMOS传感器高清影像相机获取照片，飞行高度110 m，旁向重叠率为80%，航向重叠率为80%，气象条件晴朗无风，应用Pix4D Mapper软件对含有经纬度和高程信息的照片数据进行三维重建，经过软件自带的滤波功能处理生产空间分辨率为0.074 m栅格数据，基于卫星影像空间分辨率重采样为6 m。变量测产数据获取时间为2018年9月30日，利用搭载美国Raven Envizio Pro智能测产系统的联合收割机测定含空间地理信息的产量数据，收割机割幅6 600 mm，采集量为6大垄玉米，采样频率为1 Hz，平均步进距离1.27 m。因为田间复杂情况联合收割机存在轨迹弯曲、停滞不前、信号不稳定等问题，造成产量数据点缺失及分布不均。选择产量数据点相对完整，轨迹较直，行进速度相对稳定的田块变量测产数据。产量校正基于三倍标准差法将超出值域范围的样本数据视为异常点删除。基于影像数据计算栅格像元产量均值，最终保留578个产量点。作物产量属性特征如表1所示。

1.3 研究方法

1.3.1 植被指数与地形因子

为确定最适合春玉米估产的植被指数，根据国内外作物产量预测的相关文献和玉米自身光谱特性，选择7种光谱植被指数，构建作物估产模型，即归一化植被指数(normalized difference vegetation index，NDVI)、差值植被指数(difference vegetation index，DVI)、比值植被指数(ratio vegetation index，RVI)、无蓝色波段增强型植被指数(enhanced vegetation index without a blue band，EVI2EVI2)、次生修正土壤调节植被指数(modified secondary soil adjusted vegetation index，MSAVI2)、绿度归一化植被指数(green normalized difference vegetation index，GNDVI)、修正比值植被指数(modified simple ratio，MSR)[6]。地形因子中高程、坡度、坡向、曲率等影响作物水肥运移方向和速度，使用ArcGIS 10.3提取6种地形因子高程(elevation)、坡度(slope)、坡向(aspect)、曲率(curvature)、平面曲率(plane curvature)、剖面曲率(profile curvature)。

1.3.2 春玉米遥感估产模型

为快速、准确地估算玉米产量，本研究按照2:1的比例将数据集分为训练集和验证集，各含390个和188个样本点，采用线性回归模型，利用最小二乘法，构建基于植被指数与产量的经验统计模型，如式(1)所示。应用多元逐步回归分析，评价引入地形因子的遥感反演估产模型。精度评价指标选择决定系数R2、均方根误差(root mean squared error，RMSE)，检验遥感估产模型的拟合效果，其计算公式分别为

Yield=aVIt+b

(1)

(2)

(3)

1.3.3 空间格局分析

地理空间相关性在自然界中广泛存在，距离越近的事物之间空间相关性越强，地统计学分析基于此特性在地理空间研究中应用广泛。基于变量测产数据，使用ArcGIS 10.3中Geostatistics模块分析其趋势面特征及相关地形因子趋势面特征，应用K-mean聚类和空间统计分析探索春玉米产量的空间分布格局及地形因子对实际产量的影响。IBM SPSS Statistics 22中K-mean聚类分析是空间类型划分的主要工具，该算法快速、简单具有可伸缩性，能够高效处理大数据集，实现多维度的“自然”聚类，是应用范围最广泛的聚类算法。其算法过程是在数据集中随机筛选k个样本，每个数据作为簇中心，计算其于其他数据中心距离，根据距离最近原则划分簇类。所有数据划分后，计算每个簇均值，并以均值重新聚类直至满足收敛函数的要求，可表示为

(4)

为进一步分析田块尺度中基于地形因子的实测产量空间分布规律，本研究采用自然断点法，将不同地形因子分为高、中、低3类，使用分区中产量的均值和标准差评价产量组间变化。自然断点法具有组内间距最小，组间间距最大的特点，应用其对地理空间栅格数据分类，可以更加直观反映产量的空间分布格局。

2 结果与分析

2.1 春玉米不同生育期与植被指数的遥感估产评价

探究不同生育期春玉米卫星遥感估产精度，分别将春玉米吐丝期、灌浆期和成熟期的植被指数与实测单产数据建立经验统计模型。如表2所示，春玉米3个关键生长期的R2和RMSE中，灌浆期的R2普遍在0.6以上，最高达到0.73。估产模型精度由高到低的顺序是灌浆期、吐丝期、成熟期。灌浆期模型的RMSE比吐丝期和成熟期低，即灌浆期的产量估算与实际产量的偏差最小。主要是由于作物出苗、拔节至吐丝期前，营养物质主要集中在根茎叶等生长部位。吐丝期至结穗灌浆期作物光合作用形成的淀粉、蛋白质等大量营养物质转移至产量器官，是决定玉米产量的重要时期。而成熟期田间玉米叶片水分、叶绿素下降，存在倒伏现象，估产精度较差。综上，春玉米灌浆期决定产量器官营养物质积累，是遥感反演估产的最佳时期。

表2 春玉米3个生长期植被指数估产评价模型

不同时期、不同植被指数的产量预测精度存在显著差异。在吐丝期植被指数EVI2的R2最高，可达0.24。至灌浆期，RVI的预测精度达到0.73，其余植被指数依次为MSR、EVI2、NDVI、MSAVI2、GNDVI，R2均大于0.61，对应RMSE也比DVI的估产误差小。成熟期DVI的R2达到0.25，比吐丝期大，其余植被指数R2普遍较低，约为0.19。综上，估产模型R2最高可达0.73，RMSE为388.87 kg/hm2，最优植被指数为灌浆期的RVI，遥感预测模型为

Yield=1 482.261 RVI-4 771.111

(5)

2.2 引入地形辅助信息的遥感估产评价

分析不同地形因子与春玉米产量的敏感性作用。如表3所示，农作物实测产量与高程、曲率、剖面曲率、坡度和坡向均显著相关(P<0.05)，但不同地形因子对实测产量敏感性存在显著差异，其中高程和坡度与产量的相关系数达到0.25以上，P值均小于0.01，坡度与产量的相关系数比高程大，说明高程和坡度是决定水肥运移，影响作物产量的重要因素。

表3 不同植被指数和地形因子与农作物实测产量的相关系数

如表4所示，春玉米灌浆期植被指数结合地形因子的产量估算模型及对应的R2和RMSE中，RVI的R2最高达到0.79，其余植被指数依次为MSR、EVI2、NDVI、MSAVI2、GNDVI，R2均大于0.7，相应RMSE比DVI的误差小。对比灌浆期植被指数与实测产量建立的经验统计模型，R2提升最大的为DVI估产模型，达到10.17%，其余植被指数相对提升精度均大于5.5%以上，说明地形因子作为遥感辅助信息能够减小春玉米预测产量与实际产量之间的偏差。

表4 引入地形因子的玉米预测模型结果

2.3 春玉米产量空间分布格局

如图1所示，产量与地形因子趋势面分析图中x轴表示地理空间中东西方向，y轴表示地理空间中南北方向，z轴表示采集样本数据属性值在空间分布中的大小。图1(a)表示产量在东西方向和南北方向均呈现高-低-高的显著变化趋势，与图1(b)坡度因子趋势面的规律相反，而图1(c)高程趋势面在南北方向和东西方向上均逐渐降低。说明地形对农作物实测产量存在影响。使用聚类分析工具，设置3个聚类种类，加入实测产量、高程和坡度因子。K-mean聚类分析结果显示随高程增大坡度减小，春玉米产量逐渐增高，产量最佳的聚类中心为高程307.93 m与坡度1.39°。

图1 产量与地形趋势面分布

为进一步分析产量空间分布规律，基于自然断点法将高程与坡度因子分为3类生成9种不同地形因子组合。图2为地形因子对产量影响的空间分布结果。结果显示在高海拔-高坡度区域春玉米实测产量均值最低，为7 502.64 kg/hm2，标准差为802.62 kg/hm2，该区域海拔和坡度大，是水肥流失最严重的区域，难以满足玉米生长需求。其他地形组合除了低海拔-低坡度区域没有对应产量数据外，春玉米产量平均值均大于8 000 kg/hm2，在中海拔-低坡度区域春玉米实测产量均值最高，达到9 157.63 kg/hm2，标准差为357.48 kg/hm2(表5)。原因是岗地易出现缺肥而洼地富营养化，中海拔且坡度较小的地方水肥运移速度较慢，易于与田间施肥方案确定的肥量保持一致，满足作物生长所需的营养。同时在中海拔-低坡度、中海拔-中坡度与高海拔-低坡度区域春玉米实测产量均值大于田块产量平均水平，面积占比达到46.88%。高坡度区域产量均小于田块平均水平，面积占比达到23.22%，说明不同坡度通过水肥运移过程造成的作物长势差异显著。

图2 实测产量的空间分布

表5 不同地形因子对产量空间分布的影响

3 讨论

田块尺度的春玉米产量估算，对于农场规模化经营管理，粮食产量供给侧调整，储备粮更新，大型合作社生产等均具有实际需求。在制定科学管理决策时，当年农作物实际产量及其空间分布是重要参考依据，因此需要不断提高产量的预测精度与分析其空间变异性。本研究中使用变量测产设备采集的春玉米单产数据，尽管按照测产标准设置了智能测产系统的参数，但在实际应用中仍会存在产量数据缺失，因此后期需对产量数据做精细校正。遥感影像数据空间分辨率为6 m，基于地理空间相关性原则，相同像元内作物长势较为均匀，去除异常采样点数据后，按照影像像元大小计算实测产量算数平均值的方法具有一定的合理性。

3.1 不同生育期与植被指数对春玉米估产的影响

地物空间异质性随农作物生命周期而改变，基于不同时期遥感监测的估产模型精度存在显著差异。玉米从播种至吐丝期主要生长部位为根茎叶，反映植株的干生物量积累，估产精度较低。灌浆期氮素吸收率高，至乳熟期是玉米一生中光合作用最强的时期。该时期由光合作用产生的蛋白质、脂肪、淀粉等有机质转移至玉米籽粒中，故此时估产精度高。成熟期冠层水分与叶绿素含量下降，叶片发黄导致冠层光谱植被指数的估产精度降低。本研究中灌浆期估产精度R2大于0.6，最高达到0.73，是春玉米估产的最佳监测时期。RVI对绿色作物的植被覆盖度反应敏感，从出苗至灌浆期叶面积指数增加，估产精度高。成熟期逐渐下降，且易出现玉米倒伏现象，此时估产精度最差。当前研究卫星影像空间分辨率较大难以识别倒伏阴影，在今后的研究中可以考虑结合无人机影像滤波去除作物阴影对估产光谱指数的影响。

3.2 地形辅助信息对春玉米遥感估产的影响

作物估产引入不同辅助因子优化模型能够提高精度。地形获取途径多样，在短时间内获得的多尺度空间分辨率影像提升估产精度。高程将田块尺度内耕地划分为不同坡位，决定水肥运移的方向，影响农作物生长时期的太阳辐射量和土壤温度，阳坡面日照时间长，太阳辐射量高，土壤温度相对较高，而阴坡面相反；洼地汇聚田块降水，土壤水分含量高，岗地相对较低。坡度决定农作物生长过程中水肥运移速度。水土流失现象一般出现在坡度较大的区域，基于田间降水径流方向呈条带状分布，土壤养分随耕地表层土流失造成玉米产量下降，而低洼处易富营养化，造成土壤板结。坡面平缓的区域保持适合作物生长的土壤养分和水分，区域产量高。本研究中玉米产量与高程呈正相关，与坡度为负相关。今后研究中可以结合实测土壤理化性质数据构建作物生长模型探究地形对农作物生长要素的影响。

3.3 春玉米产量空间变异分析

地理空间统计能够反映春玉米产量数据在三维空间不同方向的变化趋势，难以定量分析田块尺度中春玉米产量的空间格局变化。K-mean聚类分析确定实际产量数据与地形因子的空间类型变化，通过地形因子重分类组合分析田块尺度中不同区域春玉米产量的均值与标准差，进一步明确田块尺度中精确的空间格局变化。实际生产管理中，田间地形主要通过水肥运移与太阳辐射量和地温等因素间接影响农作物生长，在农作物整个生命周期中，它对作物生长的作用具有连续性、长期性，时间跨度越大，地形对作物产量的影响更加显著。高海拔-高坡度区域春玉米产量均值最低，为7 502.64 kg/hm2，标准差为802.62 kg/hm2。在中海拔-低坡度区域产量最高，达到9 157.63 kg/hm2，标准差为357.48 kg/hm2。其次为中海拔-中坡度区域与高海拔-低坡度区域，产量均值大于田块尺度中产量均值，面积占比总计达到46.88%。说明坡度平缓且海拔较高处春玉米产量相对较高，与K-mean聚类分析结果基本一致。今后研究中可以进一步探究多年间相同地形条件下产量时空分异规律。

4 结论

结合田块尺度地形特征、作物不同生长期SPOT 6影像与变量测产数据，使用经验统计模型确定春玉米估产最佳生长期与最优植被指数，评价引入地形辅助信息的估产模型，分析田块尺度中产量空间格局分布。得出如下结论。

(1)春玉米遥感估产的最佳时期是灌浆期，其次是吐丝期和成熟期，RVI是最优植被指数，R2达到0.73，RMSE为388.87 kg/hm2。

(2)最佳估产模型引入地形辅助信息后精度提升5.6%，R2达到0.79，RMSE为346.03 kg/hm2。

(3)空间格局分析显示高海拔且高坡度区域产量最低，中海拔与低坡度区域产量最高。

(4)下一步将研究无人机影像滤除植被阴影估产与分析作物时空分异规律。本文研究建立了一种高精度的春玉米估产方法，确定玉米最佳生长区域的地形特征，为规模化农业精细管理、耕地整理与作物种植结构调整提供了科学依据。