本研究数据采集于国家精准农业示范研究基地(北京市昌平区小汤山镇, 地理坐标40° 00'— 40° 21'N, 116° 34'— 117° 00'E, 平均高程36 m), 该区域具有典型半湿润温带季风气候特征。 其东西轴向延伸40 m(保护行不计), 南北跨度15 m, 共布设48个小区, 单小区占地面积32.5 m²。 实验区气候条件与空间配置均满足马铃薯生育期观测需求。

本实验数据来源于2017年(实验1, Exp.1)和2019年(实验2, Exp.2)两个年份的实验。 不同马铃薯品种分别在三种密度和四种氮肥梯度下进行处理, 还包括六个不同钾肥和灌溉处理的小区。 马铃薯分别于2017年3月12日和2019年3月28日播种, 2017年6月29日和2019年7月3日收获。 具体的实验区位置和实验设计见图1和表1。

图1

Fig.1

	Figure Option ViewDownloadNew Window
	图1 实验区和实验设计Fig.1 Study area and experimental design

表1

Table 1

表1(Table 1)

表1 不同地点和年份的马铃薯详细种植资料 Table 1 Detailed potato planting information in different locations and years Year No Cultivars Treatments Plots 2017 Exp.1 Zhongshu 143, Favorita Density rate: 60 000, 75 000 and 90 000 tubers· hm-2 N rate: 0, 75, 150, and 225 kg· ha-1 Pure N 18 24 Zhongshu 143 Irrigation rate: rainfed/manual irrigation 6 2019 Exp.2 Zhongshu 3, Zhongshu 5 Density rate: 60 000, 72 000 and 84 000 tubers· hm-2 N rate: 0, 112, 225, and 337 kg· ha-1 Pure N 18 24 Zhongshu 3 K rate: 0, 495, and 990 kg· hm-2 K2O 6

表1 不同地点和年份的马铃薯详细种植资料 Table 1 Detailed potato planting information in different locations and years

1.2.1 无人机多光谱数据影像获取及预处理

本研究以八旋翼无人机为遥感平台, 搭载Parrot Sequoia多光谱成像系统开展田间观测。 于2017年5月18日、 6月2日、 6月15日、 6月29日, 以及2019年5月13日、 5月28日、 6月10日、 6月20日和7月3日进行飞行作业, 获取马铃薯从现蕾期到成熟期的多光谱影像。 飞行器作业时严格保持距离冠层20米的相对高度, 航向与旁向影像重叠率分别设置为80%和85%, 各次飞行任务的起降点位与航线规划均保持高度一致。 为避免光照对实验结果的影响, 选择太阳高度角稳定的正午时段(12:00— 14:00)开展航摄工作, 同时排除云层干扰以确保光照条件均一化。 每次数据采集前均执行白板校正流程, 消除大气辐射对光谱数据的影响。 多光谱传感器的工作参数详见表2。

表2

Table 2

表2(Table 2)

表2 Parrot Sequoia多光谱相机技术参数 Table 2 Technical specifications of parrot sequoia multispectral camera 传感器 参数 多光谱传感器 单波段像素: 120万, 照片尺寸: 1 280× 960 多光谱波段 绿(550 nm± 40 nm) 红(660 nm± 40 nm) 红边(735 nm± 10 nm) 近红外(790 nm± 40 nm)

表2 Parrot Sequoia多光谱相机技术参数 Table 2 Technical specifications of parrot sequoia multispectral camera

多光谱影像处理流程基于Agisoft PhotoScan Professional平台实施。 具体步骤为: 运用地面基准点完成影像空间校准与畸变校正, 通过空间坐标解算建立初步点云模型, 经密度优化后输出高精度数字正射影像。 进一步将标准化影像数据导入ENVI 5.6中, 按实验样方边界截取目标区域, 计算不同波段反射率均值, 以此表征马铃薯冠层光谱响应特征。

1.2.2 农学参数的获取

地上生物量通过收获法获取。 选取具有代表性的3棵植株样本, 进行取样并迅速带回实验室进行茎叶分离。 样本用流水清洗之后, 在105 ℃下杀青, 随后在80 ℃下烘干48 h以上, 直到质量恒定。 最后, 称量样本的干质量, 将茎叶干质量相加, 计算出每个小区的马铃薯生物量, 最终通过群体密度和样本干质量得出每个小区的总生物量。 具体计算公式如式(1)— 式(3)

AGBl=ml3cd(1)

AGBv=mv3cd(2)

AGBt=ml+mv3cd(3)

式(1)— 式(3)中, m_l和m_v分别为所有叶片和茎秆部分的干质量; 3代表马铃薯取样的3棵植株; c_d为种植密度。

1.2.3 气候数据的获取

本研究使用的气象数据来自欧洲中期天气预报中心(ECMWF)发布的ERA5数据集, 该数据集提供全球范围内的历史气象栅格数据, 能够有效弥补实际观测数据中的缺失或不足。 从ERA5数据集中获取了对应研究区和时间段的逐日最高气温和最低气温, 空间分辨率为0.125° , 能够满足对马铃薯发育进程研究的量化需求。

1.2.4 有效积温的计算

马铃薯是一种耐寒的作物, 生长需要较凉的气候。 马铃薯幼芽生长的最佳温度是13~18 ℃。 如果温度低于4 ℃, 种薯就无法发芽。 如果气温超过21 ℃, 马铃薯的生长会受到影响, 生长速度明显变慢。 当土壤温度超过25 ℃时, 块茎的生长几乎会停止。 在中原地区, 由于7月的平均气温超过25 ℃, 因此只在早春和秋季种植马铃薯。 以实验区为例, 北京地处华北平原, 3月中旬气温超过4 ℃, 而7月初气温可达25 ℃。 因此, 本研究将4 ℃定为马铃薯的生物学最低温度, 用于计算有效积温。 有效积温的计算公式如式(4)和式(5)

GDD=∑DOY=1DOY(Tavg-Tbase)(4)

Tavg=Tmax+Tmin2, Tbase< Tavg< TupperTavg=Tbase, Tavg≤TbaseTavg=Tupper, Tavg≥Tupper(5)

式(4)和式(5)中, GDD为有效积温; DOY是儒略日指标, 代表一年中的某一天; T_avg为日平均气温; T_base为马铃薯生长发育所需要的基本温度, 本研究使用4 ℃作为基本温度; T_upper为马铃薯生长发育所需要的上限温度; T_max和T_min分别为日最高气温和日最低气温, 上述温度的单位均为℃。

1.3.1 植被指数的选取

以往研究中, 许多植被指数已广泛应用于马铃薯地上生物量的估算, 但与马铃薯叶生物量相关的植被指数研究较少。 本研究结合前人的研究成果, 选取了与AGB密切相关的植被指数。 初步筛选了14个植被指数, 包括NDVI、 GNDVI、 MSR和MTVI1等, 具体的计算公式列在表3中。

表3

Table 3

表3(Table 3)

表3 植被指数 Table 3 Vegetation index 缩写 植被指数名称 公式 参考文献 DVI 差异植被指数 (difference vegetation index) DVI=NIR-RED [16] GNDVI 绿色归一化差异植被指数 (green normalized difference vegetation index) GNDVI=(NIR-GRE)/(NIR+GRE) [14] MSR 改进型简单比值指数 (modified simple ratio) MSR=(NIR/RED-1)/ (NIR/RED+1) [14] MSAVI 修正土壤调整植被指数 (modified soil-adjusted vegetation index) MSAVI=NIR-0.5× (2×NIR+1)2-8×(NIR-RED)+0.5 [15] NDVI 归一化差异植被指数 (normalized difference vegetation index) NDVI=(NIR-RED)/(NIR+RED) [16] NLI 非线性植被指数 (non-linear index) NLI=(NIR2-RED)/(NIR2+RED) [16] OSAVI 优化土壤调整植被指数 (optimized soil-adjusted vegetation index) OSAVI=1.16× (NIR-RED)/(OIR+RED+0.16) [14] RDVI 再归一化差异植被指数 (renormalized difference vegetation index) RDVI=(NIR-RED)/ NIR+RED [14] RVI 比值植被指数 (ratio vegetation index) RVI=NIR/RED [16] SAVI 土壤调整植被指数 (soil-adjusted vegetation index) SAVI=1.5× (NIR-RED)/(NIR+RED+0.5) [16] MTVI1 改进的三角植被指数1 (modified triangular vegetation index 1) MTVI1=1.2× (1.2× (NIR-GRE)-2.5× (RED-GRE)) [15] MTVI2 改进的三角植被指数2 (modified triangular vegetation index 2) MTVI2= 1.5×(1.2×(NIR-GRE)-2.5×(RED-GRE))((2×NIR+1)2-(6×NIR-5×RED))-0.5 [15] NDRE 归一化红边差异指数 (normalized difference red edge index) NDRE=(NIR-REG)/(NIR+REG) [14] TVI 三角植被指数 (triangular vegetation index) TVI=60× (NIR-GRE)-100× (RED-GRE) [15] 注: GRE、 RED、 REG和NIR分别表示绿光、 红光、 红边光和近红外波段的反射率值 Note: GRE, RED, REG, and NIR represent the reflective values of the green, red, red edge, and near-infrared bands respectively

表3 植被指数 Table 3 Vegetation index

1.3.2 植被指数的筛选方法

植被指数的组合筛选采用了三种评价指标: 相关性分析(R)、 变量投影重要性分析(VIP)和随机森林特征重要性评估(FIM)。 通过这些方法, 从不同角度衡量植被指数与LGB全生育期之间的关系。 综合分析评价结果, 能够全面评估每个植被指数对LGB反演的贡献, 从而筛选出最具反演价值的指数或指数组合, 优化模型结构, 减少冗余变量, 提高LGB反演模型的精度和效率。

通过VIP分析筛选最具有解释意义的变量, 优化模型结构并减少冗余变量。 VIP值越大, 自变量对因变量的解释力越强。 VIP值大于1表示对叶片生物量有显著解释作用, 0.8~1之间表示中等程度的解释意义, 小于0.8则无明显解释意义, 可考虑剔除。 计算公式如式(6)所示

VIP=p∑h=1mRd(Y, th)whj2∑h=1mRd(Y, th)(6)

式(6)中, p是参与建模的自变量个数, m是原变量中所提取的成分个数, t_h代表第h个成分, Rd(Y, t_h)代表成分t_h对Y的解释能力, w_hj是特征向量w_h的第j个分量, 表示该变量在主成分上的权重。

随机森林特征重要性评估(FIM)通过分析特征对模型估算能力的贡献, 帮助选择重要特征并优化模型。 特征重要性由袋外错误率均值和标准差决定, FIM结果的取值包括正、 零和负。 正数表示特征扰动使模型估算能力减弱, 负数则表示增强, 0表示无影响。 为了直观比较, 本研究将FIM结果进行标准化。

1.4.1 AGB动态分配模型的构建

根据马铃薯植株生物量分配规律, 本研究从整体与部分之间的关系出发, 分析了LGB与AGB比值在整个生育期内的变化规律。 基于此, 构建了以有效积温(GDD)为自变量、 LGB/AGB比值为因变量的生物量动态分配模型, 旨在更好地理解作物生物量的动态变化过程, 模型如式(7)所示

LGBiAGBi=f(GDD)(7)

式(7)中, LGB_i为叶片生物量(kg· ha^-1), AGB_i为地上生物量(kg· ha^-1), f(GDD)为LGB/AGB随生育时期的函数。

不同作物的LGB/AGB与DOY的关系可能不同, 这取决于实际数据的特点和作物的生物学特性。 拟合LGB/AGB与GDD的关系可以采用多种数学模型, 如线性、 幂函数或对数方程。 为了建立准确的模型, 需要在不同生育阶段详细测定作物生物量, 并记录生长和收获日期以计算GDD, 结合环境因素分析, 最终通过实验数据确定最合适的模型。

1.4.2 基于AGB动态分配模型的AGB估算

本研究通过结合AGB动态分配模型与叶片生物量估算模型, 构建了一个以植被指数(VIs)和有效积温(GDD)为输入变量的马铃薯AGB估算模型。 该模型的核心思想是将植被指数估算的LGB作为中间量, 再结合AGB动态分配模型确定的LGB/AGB与GDD的关系, 构建AGB估算模型。 通过这种方式, 模型能够考虑到作物在不同生长阶段的生物量比例变化, 从而有效克服“ 饱和现象” 带来的限制, 提高AGB估算的精度和可靠性。 估算模型公式如式(9)所示

LGB=f(VIs)(8)

AGB=f(VIs)f(GDD)(9)

有效积温通过欧洲中期天气预报中心(ECMWF)发布的ERA5数据集计算, 获取实验区每日的日平均气温, 并减去马铃薯的最低生物学温度4 ℃, 得到每日有效积温, 然后, 将每日有效积温累加, 得到全生育期的有效积温数据; 植被指数的筛选采用VIP分析, 通过SIMCA-P软件计算植被指数的VIP值, 量化评估植被指数对LGB和AGB的解释能力; 随机森林特征重要性评价则利用Python 3.11计算各植被指数的特征重要性。

在构建马铃薯AGB动态分配模型时, 选取2017年和2019年实验区全生育期的重复一和重复二共计320个样本数据用于AGB动态分配模型的构建。 重复三共计160个样本数据, 其中划分为70%的训练集和30%的测试集, 用于机器学习方法的AGB估算。 为了估算马铃薯地上生物量, 本研究采用三种机器学习方法: 其中两种非线性回归方法高斯过程回归(Gaussian process regression, GPR)和偏最小二乘回归(partial least square regression, PLSR)以及一种线性回归方法随机森林(random forest, RF)。 模型的精度评估采用决定系数(R²)、 均方根误差(RMSE)和归一化均方根误差(NRMSE)作为评价指标。

马铃薯不同生育期的LGB和AGB变化如图2所示。 2017年和2019年的数据变化趋势如图3所示: 在五个生育时期内, 叶生物量(LGB)和地上生物量(AGB)变化趋势基本一致。 从现蕾期开始, 叶生物量逐渐上升, 到中后期的淀粉积累期达到顶峰, 随后开始下降。 与叶生物量变化相似, 地上生物量也在后期下降。 五个时期之间, LGB和AGB的变化规律相似, 表明LGB在AGB中占有很大比例。 因此, 通过估算LGB的变化, 可以间接估算AGB, 这为全生育期内马铃薯AGB的估算提供了可行性。

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 马铃薯各生育期LGB和AGB变化Fig.2 LGB and AGB dynamics in potato growth stages

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图3 实验区每日平均气温(℃)Fig.3 Daily mean temperature (℃)

从ERA5数据集中下载了2017年和2019年每日地表2 m的平均气温数据, 用以计算实验区马铃薯从1月到7月初的全生育期气温数据(DOY=年内天数, 1月1日为第1天), 如图3所示。 马铃薯的生长温度要求因地区纬度不同而有所差异。 以实验区为例, 北京地处华北平原, 3月中旬气温超过4 ℃, 7月初气温可达25 ℃。 因此, 本研究以4 ℃为基准进行有效积温(GDD)的计算, 首先将每日平均气温减去4 ℃, 然后逐日累加得到有效积温数据。 2017年, 马铃薯播种日期为3月12日, 收获日期为6月29日, DOY为180天, 累计有效积温为1 611.75 ℃; 2019年, 马铃薯播种日期为3月28日, 收获日期为7月3日, DOY为184天, 累计有效积温为1 697.62 ℃。

本研究根据三种评价指标进行分析, 选择了6种在马铃薯全生育期被广泛应用且表现优秀的植被指数(MTVI2、 NDVI、 RVI、 MTVI1、 GNDVI和TVI)作为模型输入变量。 植被指数的选择遵循VIP值大于1、 相关性系数和FIM尽可能影响大的原则。 这6个变量在FIM、 R²和VIP三项指标中表现出色, 如图4和表4所示。 具体而言, MTVI2和MTVI1在捕捉高密度植被光谱特性方面具有优势, 能够有效反映生物量较大的植被区域; NDVI和GNDVI是通用性较强的植被指数, 适用于评估植被覆盖度和低密度植被的生物量特征; RVI直接反映植被密度, 与叶面积指数(LAI)和生物量

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图4 三种评价指标结果图Fig.4 Results of three evaluation metrics

表4

Table 4

表4(Table 4)

表4 所选指数的评价指标计算结果 Table 4 Evaluation metrics results for selected indices 植被指数 评价指标 R VIP FIM MTVI1 0.40 1.38 0.64 MTVI2 0.41 1.31 0.69 NDVI 0.35 0.88 0.73 RVI 0.36 0.89 0.68 GNDVI 0.24 1.10 0.74 TVI 0.40 1.19 0.53

表4 所选指数的评价指标计算结果 Table 4 Evaluation metrics results for selected indices

具有显著相关性; TVI则对植被生长变化敏感, 能够捕捉生长旺季的生物量变化。 这些指数具有较强的互补性, 能够从多个维度(光谱反射、 植被覆盖和生长变化)提供可靠的输入, 显著提升了反演叶片生物量模型的精度和稳定性。

2.4.1 AGB动态分配模型的构建

基于马铃薯全生育期生物量变化规律的分析, 本研究引入了有效积温(GDD)作为生育阶段的量化指标, 成功构建了AGB动态分配模型, 如式(10)所示。 根据图5的结果, 马铃薯的LGB/AGB与GDD之间的关系可以通过多项式线性函数最佳拟合, 拟合精度为R²=0.88, 表明模型具有较高的准确性, 能够有效描述叶片生物量与整体植株生物量之间的关系。

图5

Fig.5

	Figure Option ViewDownloadNew Window
	图5 AGB动态分配模型构建Fig.5 AGB dynamic allocation model construction

模型显示, LGB/AGB随GDD的增加呈现先上升后下降的趋势。 在生长初期, 随着GDD增加, LGB/AGB比值上升, 并在GDD为920 ℃时达到最高值。 此时, 叶片生物量占地上生物量的比例较高, 表明光合作用的积累优先支持叶片的生长。 随着生长季节的推进, 叶片生物量增长到顶峰并达到饱和, 比值为0.69, 对应的GDD约为920 ℃, 此时处于块茎增长期。 随着植株进入生殖生长阶段, 地上生物量继续增加, 但叶片生物量几乎不变, 其占比开始下降, LGB/AGB比值随之下降。 最终, 在成熟期植株的营养积累放缓, 地下块茎完成积累, 地上植株开始枯萎, 叶片生物量占比达到最低值。

LGBAGB=-0.4009(GDD2)+0.7399(GDD)+0.3437(10)

2.4.2 AGB估算结果与分析

为验证构建的AGB动态分配模型的效果, 本研究将全生育期筛选出的6种植被指数与随机森林(RF)、 高斯过程回归(GPR)和偏最小二乘回归(PLSR)方法结合, 建立AGB估算模型, 具体结果见表5。 PLSR方法与AGB动态分配模型结合后的精度最低, 建模集和验证集的R²分别为0.34和0.20, RMSE分别为460.28和412.98 kg· ha^-1, NRMSE分别为17.77%和22.33%。 相比之下, GPR和RF方法的精度更高, 其中GPR稍差一些(建模: R²=0.55, RMSE=333.14 kg· ha^-1, NRMSE=12.8%; 验证: R²=0.50, RMSE=455.67 kg· ha^-1, NRMSE=18.9%), 而RF方法精度最高(建模: R²=0.80, RMSE=256.73 kg· ha^-1, NRMSE=9.91%; 验证: R²=0.76, RMSE=211.91 kg· ha^-1, NRMSE=11.46%)。 在三种方法中, RF表现最好, 建模集和验证集的NRMSE分别为9.91%和11.46%, 表明该模型在训练数据和验证数据上的估算误差较小。 尽管验证集的R²和RMSE略逊于建模集, 但两者差异较小, 表明该模型具有较好的泛化能力, 能够在新的数据集上保持稳定的估算性能。 综上所述, RF-AGB动态分配模型在建模集和验证集上的表现最好, 估算结果与实测值之间的相关性较高, 证明该模型在实际应用中具有较大潜力, 可用于准确估算马铃薯全生育期的地上生物量。

表5

Table 5

表5(Table 5)

表5 基于AGB动态分配模型的马铃薯AGB估算 Table 5 Potato AGB estimation based on the AGB dynamic allocation model Methods Calibration Validation R2 RMSE/(kg· ha-1) NRMSE/% R2 RMSE/(kg· ha-1) NRMSE/% RF 0.80 256.73 9.91 0.76 211.91 11.46 GPR 0.55 333.14 12.80 0.50 455.67 18.90 PLSR 0.34 460.28 17.77 0.20 412.96 22.33

表5 基于AGB动态分配模型的马铃薯AGB估算 Table 5 Potato AGB estimation based on the AGB dynamic allocation model

为了进一步评估AGB动态分配模型的可靠性, 本研究设计了基于机器学习方法的对比实验, 对比了三种机器学习方法在直接估算与结合模型时的表现(表6)。 通过方法特性与数据结果的交叉分析, 揭示动态模型的作用机制。 方法与机器学习反演LGB相同, 具体而言, 通过计算植被指数与AGB之间的三种相关性评价方法, 筛选出最佳的植被指数组合, 并使用这些指数作为输入变量, 直接估算AGB。 最终, 将这些结果与基于AGB动态分配模型的估算结果进行精度对比, 以全面评估本研究模型在马铃薯AGB估算中的表现。 三种机器学习方法直接估算AGB结果如表6所示。

表6

Table 6

表6(Table 6)

表6 三种机器学习方法直接估算AGB结果 Table 6 Direct AGB estimation via 3 machine learning methods Methods Calibration Validation R2 RMSE/(kg· ha-1) NRMSE/% R2 RMSE/(kg· ha-1) NRMSE/% RF 0.52 350.67 13.51 0.35 507.93 21.10 GPR 0.75 301.12 11.60 0.43 496.21 20.61 PLSR 0.29 475.91 18.37 0.15 392.51 22.14

表6 三种机器学习方法直接估算AGB结果 Table 6 Direct AGB estimation via 3 machine learning methods

根据表6和图6的结果, RF在直接估算时, 验证集R²(0.35)显著低于校准集(0.52), 表明其决策树结构易受训练数据噪声干扰。 然而, 结合动态模型后, 验证集R²提升至0.76, NRMSE降低至11.46%。 通过引入马铃薯AGB动态分配模型, 为RF提供了生理学先验知识, 从而有效抑制了过拟合。 GPR在直接估算时, 校准集R²高达0.75, 但验证集R²骤降至0.43, 反映其核函数在光谱特征上过度拟合训练数据。 结合动态模型后, 验证集R²提升至0.50, 表明GDD的时序信息部分弥补了核函数在高维数据中的不稳定性, 尽管其非线性建模优势仍未被充分挖掘。 PLSR在直接估算时, 验证集R²仅为0.15, 结合动态模型后略有提升至0.20, 两者表现最差。 尽管动态模型引入了非线性生物量分配规律, 但PLSR的线性潜在变量无法有效融合此类信息, 凸显了其在复杂农学场景中的局限性。

图6

Fig.6

	Figure Option ViewDownloadNew Window
	图6 (a)— (c)为基于AGB动态分配模型的机器学习方法估算AGB结果; (d)— (f)为机器学习方法直接估算AGB结果Fig.6 (a)— (c) Machine Learning-based AGB estimation with dynamic allocation models; (d)— (f) Machine Learning-driven AGB retrieval via direct inversion

综合来看, 动态模型在改善RF和GPR的估算精度方面表现显著, 但对于PLSR来说, 其线性假设与动态模型的非线性特征协同失效, 未能有效提升估算性能。 基于AGB动态分配模型的AGB估算精度明显优于直接估算AGB, 三种机器学习方法的精度均有所提升。 其中, RF模型的精度提升最大(建模: R²=0.80, NRMSE=9.91%; 验证: R²=0.76, NRMSE=11.46%)。 从图6(d)— (f)可见, 直接使用机器学习方法估算马铃薯AGB时, 会出现过拟合现象, 导致全生育期的AGB估算效果较差且可靠性较低。 相比之下, 基于AGB动态分配模型的估算结果在R²和NRMSE上均有所提升, 能够更准确地将建模集上的规律应用于全生育期的验证数据, 提供更可靠的AGB估算。