马铃薯试验在北京市昌平区小汤山镇国家精准农业研究示范基地进行。 基地位于北纬40° 10'34″, 东经116° 26'39″, 平均海拔为36 m, 气候类型为暖温带半湿润大陆性季风气候, 年均降水量为640 mm, 年均温度为10 ℃, 年均无霜期180 d。 试验为小区随机试验设计, 采用了中薯5(P1)和中薯3(P2)2个均为早熟的马铃薯品种, 试验田共计48个小区, 每个小区面积为6.5 m× 5 m。 试验区共设密度试验(N区)、 氮素试验(S区)、 钾肥试验(K区)3个试验区, 每种试验3种重复, 每种重复进行不同程度的密度、 氮素和钾肥处理, 其中密度设置了3种水平(T1: 60 000株· hm^-2、 T2: 72 000株· hm^-2、 T3: 84 000株· hm^-2), 氮素设置了4种水平(N0: 0 kg尿素· hm^-2、 N1: 244.65 kg尿素· hm^-2、 N2: 489.15 kg尿素· hm^-2(正常处理, 15 kg纯氮)、 N3: 733.5 kg尿素· hm^-2), 钾肥设置了2种水平(K0: 0 kg钾肥· hm^-2、 K1: 970.5 kg钾肥· hm^-2(N区和S区均为K1处理)、 K2: 1 941 kg钾肥· hm^-2)。 为了更好的获取试验田的位置, 在试验小区周围均匀布控11个GCP(k1— k11由0.3 m× 0.3 m的木板和埋于地下的木桩组成, 其上有黑白标志的塑料板, 目的是准确确定木板的中心位置)并用差分GPS测定其三维空间位置, 详细的试验设计见图1所示。

图1

Fig.1

	Figure Option ViewDownloadNew Window
	图1 马铃薯试验位置和试验设计Fig.1 The experimental area and design of potato

选择天气晴朗无云的时间分别在现蕾期(2019年5月13日)、 块茎形成期(2019年5月28日)、 块茎增长期(2019年6月10日)和淀粉积累期(2019年6月20)进行无人机数码和高光谱遥感作业。 飞行时间为各时期的中午12:00, 此时太阳光照强度稳定, 飞行高度为50 m。 无人机搭载的高清数码相机型号为Cyber-shot DSC-QX100, 其主要参数为: 尺寸62.5 mm× 62.5 mm× 55.5 mm, 质量179 g; 2 090万像素CMOS传感器; 焦距10 mm(定焦拍摄)。 无人机搭载的高光谱传感器为德国Cubert公司生产的UHD185机载成像光谱仪(测量过程中保证光谱仪垂直向下), 该传感器大小为195 nm× 67 nm× 60 nm, 质量470 g, 波长范围为450~950 nm, 共有125个光谱通道, 采样间隔4 nm, 光谱分辨率8 nm, 数字分辨率12位。 遥感作业前, 在地面利用黑白板进行高光谱辐射校正, 并保持各时期的飞行航线一致。

无人机获取的数码影像需要进行拼接处理, 首先将带有经纬度、 高度、 偏航角、 俯仰角和旋转角共6种POS信息的每张影像导入PhotoScan软件中, 还原影像拍摄时的姿态和空间位置; 将GCP添加到软件, 利用GCP的三维坐标信息优化照片对齐, 生成飞行区域的稀疏点云; 基于稀疏点云构建密集点云, 进而生成飞行区域的空间格网并建立纹理信息; 生成飞行区域马铃薯试验田的数字正射影像(digital orthophoto map, DOM)。

无人机获取的高光谱影像处理主要包括辐射校正、 影像拼接、 影像融合和光谱提取4部分。 首先完成由高光谱影像DN值(digital number, DN)到地表反射率的辐射定标。 利用Agisoft PhotoScan 软件进行高光谱影像拼接, 变为cue格式数据, 再提取子波段为jpg格式。 最后再对各个拼接的子波段影像进行合并生成马铃薯试验田的高光谱影像数据。 基于Arcmap10.2软件, 绘制马铃薯每个小区矢量数据, 基于IDL语言统计每个小区平均光谱作为马铃薯冠层光谱, 得到高光谱反射率数据。

地面数据获取与无人机遥感作业同步进行, 分别获取马铃薯现蕾期到淀粉积累期4个关键时期的AGB数据。 马铃薯AGB通过收获法获取, 选取能够代表每个小区长势水平的3棵植株样本, 将其茎叶分离, 随后用清水洗净, 105 ℃杀青, 80 ℃烘干48 h以上, 直到质量恒定再称量。 将植株茎和叶的干质量求和得到样本干质量, 最后通过样本干重和种植密度得到每个小区的马铃薯AGB。

参考相关文献, 选取在AGB监测方面较好的高光谱植被指数, 如PBI, PRI, LCI, EVI, OSAVI和DVI这6种, 用于构建AGB模型, 从而监测AGB, 具体的表达式见表1。

表1

Table 1

表1(Table 1)

表1 植被指数计算公式 Table 1 Calculation formula of vegetation indexes 植被 指数 公式 参考 文献 PBI R810/R560 [8] PRI (R530-R570)/(R530+R570) [8] LCI (R850-R710)/(R850+R670) [8] EVI 2.5× (802-670)/(802+6× 670-7.5× 474+1) [7] OSAVI 1.16× (R800-R670)/(R800+R670+0.16) [7] DVI R890-R670 [7]

表1 植被指数计算公式 Table 1 Calculation formula of vegetation indexes

IWD基于DWT技术将一幅大小为m× n的影像利用包括低通(lowpass, LP)和高通(highpass, HP)的级联滤波器组进行二维影像分解, 最终得到4幅子影像, 分解过程如图2所示。 通过IWD获得的各个子影像大小为原始影像的一半, 其中包含了4种信息, 分别为低频信息LL、 水平方向的高频信息LH、 垂直方向的高频信息HL和对角线方向的高频信息HH。 3种高频信息代表了丰富的植被冠层光谱, 忽略了土壤背景的影响, 从而能够有效的估算作物不同生育期的AGB。 数码影像含有R, G和B 3个波段, 为了能够充分利用冠层影像的光谱信息, 将3个波段的同种高频信息按照均等权重进行平均, 最终每个小区得到上述3种高频信息。

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 影像小波分解过程Fig.2 Image wavelet decomposition process

CWT是一种集信号处理与分析为一体的强有力积分变换工具, 可对信号进行细化分解, 凸显出信号中隐藏的更多重要信息, 在高光谱数据降维处理中被广泛应用。 其原理是通过平移和缩放的小波母函数与马铃薯原始冠层光谱反射率数据进行卷积运算, 将光谱数据分解成不同尺度下和不同波长的小波能量系数。 本研究利用CWT分别获取每个生育期前10尺度下的小波能量系数, 其表达式为

Wf(a, b)=< f, Ψa, b> =∫-∞+∞fλΨa, bλdλ(1)

Ψa, b(λ)=1aΨλ-ba(2)

式(1)和式(2)中, Ψ _{a, b}(λ )为小波母函数(采用Gaussian二阶导函数), a为伸缩因子或尺度因子, b为平移因子, λ 为马铃薯冠层光谱反射率的波段(450~950 nm), f(λ )为马铃薯冠层光谱反射率数据, W_f(a, b)为小波能量系数, 是由分解尺度(1, 2, …, m)和波长(1, 2, …, n)组成的二维矩阵。

通过PLSR分析马铃薯每个生育期的高光谱遥感数据, 构建AGB估算模型。 PLSR是利用了多元线性回归(multiple linear regression, MLR)、 典型相关分析和主成分分析(principal components analysis, PCA)为一体, 提供一种以多个自变量构建线性回归建模方法。 当输入参数个数较多时, 且存在多重相关性, 而观测变量数据较少时, PLSR构建的模型具有传统的典型回归分析所没有的优点, 利用最小化误差的平方和构建最佳的模型, 以达到较好的估算效果。 因此得到的分析结果, 除了提供一个合理的模型外, 还可以同时完成一些类似于PCA和典型相关分析的流程。

对每个生育期分别挑选2/3样本数据(32个, 重复1和重复2)作为建模集, 1/3样本数据(16个, 重复3)作为验证集构建马铃薯AGB估算模型。 采用决定系数(coefficient of determination, R²)、 均方根误差(root mean square error, RMSE)、 标准均方根误差(normalized root mean square error, NRMSE)评价模型的精度。 R²越接近于1, RMSE和NRMSE越低, 其估测模型的精度越高。

将选取的6种植被指数分别与马铃薯每个生育期的AGB进行相关性分析, 得到相关性系数绝对值, 结果见表2所示。 由表可知, 每个生育期的植被指数与AGB均达到0.01显著水平, 相关性较高, 能够很好地估算AGB。 每个生育期与AGB相关性最高的植被指数不同, 现蕾期, EVI的相关性系数绝对值最高, 为0.673。 块茎形成期, DVI的相关性系数绝对值最高, 为0.734。 块茎增长期和淀粉积累期, 相关性系数绝对值最高的都为PBI, 相关系数分别为0.748和0.651。

表2

Table 2

表2(Table 2)

表2 植被指数与马铃薯AGB的相关性 Table 2 Correlation between vegetation indexes and AGB of potato 植被指数 每个生育期的相关性 现蕾期 块茎形成期 块茎增长期 淀粉积累期 PBI 0.657 0.730 0.748 0.651 PRI 0.643 0.733 0.745 0.638 LCI 0.639 0.656 0.731 0.624 EVI 0.673 0.703 0.735 0.642 OSAVI 0.641 0.716 0.732 0.622 DVI 0.617 0.734 0.742 0.650

表2 植被指数与马铃薯AGB的相关性 Table 2 Correlation between vegetation indexes and AGB of potato

利用1.5方法对马铃薯4个生育期的数码影像进行IWD, 将得到的3种高频信息与AGB进行相关性分析, 得到相关性系数绝对值, 结果见表3。 由表3可知, 不同生育期提取的3种高频信息均与马铃薯AGB达到0.01显著水平。 整个生育期, 同种高频信息和AGB的相关性变化趋势与植被指数相关性变化趋势的保持一致。 每个生育期得到的3种高频信息与AGB的相关性高于植被指数。

表3

Table 3

表3(Table 3)

表3 高频信息与马铃薯AGB的相关性 Table 3 Correlation between high frequency information and AGB of potato 生育期 高频信息相关性 LH HL HH 现蕾期 0.663 0.672 0.691 块茎形成期 0.745 0.736 0.750 块茎增长期 0.758 0.747 0.762 淀粉积累期 0.674 0.688 0.701

表3 高频信息与马铃薯AGB的相关性 Table 3 Correlation between high frequency information and AGB of potato

每个生育期前10尺度下的小波系数与AGB决定系数如图3所示。 由图可知, 每个生育期的小波系数与AGB的相关性差异明显, 整个生育期呈现先升高后降低的趋势, 这与通过原始光谱反射率构建的植被指数和利用IWD技术提取的高频信息与AGB的相关性变化趋势保持一致。 经过分析发现, 现蕾期和淀粉积累期, 当分解尺度在2⁷以上时, 分解得到的小波能量系数已无法表征马铃薯的AGB信息, 而块茎形成期和块茎增长期, 当分解尺度在2⁷以上时, 小波能量系数与AGB在所有波长处均呈正相关关系且相关系数均达到0.5以上。 现蕾期(相关系数绝对值达到0.6以上)、 块茎形成期(相关系数绝对值达到0.7以上)、 块茎增长期(相关系数绝对值达到0.7以上)和淀粉积累期(相关系数绝对值达到0.7以上)的主要波段区域见表4所示。

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图3 小波系数与马铃薯AGB的相关性Fig.3 Correlation between wavelet coefficient and AGB of potato

表4

Table 4

表4(Table 4)

表4 各生育期与AGB相关的主要波段区域(nm) Table 4 The main bands related to AGB in different growth periods (nm) 尺度 现蕾期 块茎形成期 块茎增长期 淀粉积累期 C1 466~478, 698~714, 734~766 698~718, 738~778 470~474, 694~706, 734~762 574, 698~710, 734~742 C2 466~490, 606~630, 694~714, 738~798 586~642, 694~714, 742~794, 886~894 650~658, 682~706, 738~774, 790~802, 858~862, 890~910 578~582, 610~622, 694~706, 738~762 C3 462~514, 602~618, 682~718, 746~854 602~618, 674~714, 750~838, 870~894 610~706, 738~898, 950 598~614, 750~778 C4 474~482, 570~718, 746~834 558~710, 758~902 566~706, 734~890, 950 C5 458~694 462~678, 770~950 454~670, 722~950 C6 454~586 454~562, 854~950 454~550, 726~950 C7 806~950

表4 各生育期与AGB相关的主要波段区域(nm) Table 4 The main bands related to AGB in different growth periods (nm)

基于表4每个生育期挑选的与AGB高度相关的敏感波段区间中, 筛选出绝对相关系数较高的前10个波段作为输入变量, 估算马铃薯AGB。 不同生育期筛选的不同尺度下最优波段的结果见表5所示。 由表与数据可知, 每个生育期选中的波段尺度主要集中在第1, 2, 3和4尺度, 每种尺度下筛选的敏感波段与AGB均达到0.01显著水平。 每个生育期筛选的不同尺度下的波段与AGB的相关性高于植被指数和高频信息。

表5

Table 5

表5(Table 5)

表5 前10敏感波段及对应的尺度 Table 5 Top 10 sensitive bands and corresponding scales 序号 现蕾期 块茎形成期 块茎增长期 淀粉积累期 尺度 波长/nm |R| 尺度 波长/nm |R| 尺度 波长/nm |R| 尺度 波长/nm |R| 1 C1 474 0.719 C3 606 0.763 C2 750 0.776 C1 702 0.732 2 C2 486 0.703 C2 618 0.759 C4 658 0.776 C1 706 0.732 3 C3 502 0.702 C3 602 0.755 C2 742 0.774 C2 582 0.730 4 C3 498 0.701 C2 614 0.755 C1 694 0.769 C2 618 0.730 5 C3 506 0.700 C2 642 0.754 C4 566 0.769 C1 710 0.728 6 C3 494 0.699 C1 738 0.753 C3 650 0.769 C2 750 0.727 7 C3 490 0.697 C1 702 0.751 C3 694 0.768 C3 602 0.724 8 C1 478 0.692 C1 706 0.751 C3 670 0.766 C3 750 0.724 9 C3 482 0.692 C3 610 0.751 C1 474 0.765 C2 614 0.721 10 C3 510 0.689 C1 762 0.750 C1 762 0.765 C3 766 0.721

表5 前10敏感波段及对应的尺度 Table 5 Top 10 sensitive bands and corresponding scales

基于VIs, VIs+IWD和VIs+CWT3种变量采用PLSR方法建立每个生育期的AGB估算模型, 并计算每种模型的精度评价指标, 其结果见表6所示。 以3种变量估算马铃薯AGB均在块茎增长期表现效果最佳, 整个生育期, 同种模型变量估算AGB的精度从现蕾期到块茎增长期逐渐变好, 随后开始变差。 每个生育期仅以VIs为变量估算马铃薯AGB的精度较低, 现蕾期到块茎增长期建模R²范围是0.54~0.70, RMSE范围是272.84~98.88 kg· hhm^-12, NRMSE的变化为19.22%~11.63%, 块茎增长期到淀粉积累期建模R²范围是0.70~0.62, RMSE范围是98.88~212.63 kg· hm^-2, NRMSE的变化为11.63%~17.62%, 验证效果同建模效果相一致, 表明基于VIs估算AGB效果较差。 以VIs+IWD为模型因子, 现蕾期到块茎增长期建模R²从0.65变化到0.78, RMSE从254.63 kg· hm^-2变化到86.45 kg· hm^-2, NRMSE从17.64%变化到10.17%, 验证R²(0.69~0.82)也不断增大, RMSE和NRMSE都不断减少, 块茎增长期到淀粉积累期, 建模R²范围从0.78下降到0.66, RMSE从86.45 kg· hm^-2增加到205.05 kg· hm^-2, NRMSE从10.17%增加到16.53%, 验证效果与建模效果保持一致, 验证R²(0.82~0.71)逐渐降低, RMSE和NRMSE都逐渐增加。 以VIs+CWT为模型输入变量, 现蕾期到块茎增长期建模R²范围是0.68~0.85, RMSE范围是244.57~74.25 kg· hm^-2, NRMSE的变化为16.83%~9.27%, 块茎增长期到淀粉积累期建模R²范围是0.85~0.72, RMSE范围是74.25~181.28 kg· hm^-2, NRMSE的变化为9.27%~14.27%, 验证R²也是先升高后降低以及RMSE和NRMSE先降低后升高。 对比每个生育期的3种AGB估算模型可知, VIs结合IWD或CWT明显提高了AGB估算精度, 且VIs+CWT为变量估算AGB效果更好。

表6

Table 6

表6(Table 6)

表6 基于不同变量的马铃薯AGB估算精度对比 Table 6 Comparison accuracy of potato AGB estimation based on different variables 生育期 变量 建模 验证 R2 RMSE/(kg· hm-2) NRMSE/% R2 RMSE/(kg· hm-2) NRMSE/% VIs 0.54 272.84 19.22 0.64 215.89 18.44 现蕾期 VIs+IWD 0.65 254.63 17.64 0.69 200.33 16.88 VIs+CWT 0.68 244.57 16.83 0.73 180.20 15.53 VIs 0.64 189.95 16.64 0.70 160.33 14.22 块茎形成期 VIs+IWD 0.72 162.22 14.37 0.78 141.89 12.18 VIs+CWT 0.75 150.58 13.45 0.82 133.46 11.93 VIs 0.70 98.88 11.63 0.75 88.69 9.74 块茎增长期 VIs+IWD 0.78 86.45 10.17 0.82 79.95 8.49 VIs+CWT 0.85 74.25 9.27 0.88 72.06 7.86 VIs 0.62 212.63 17.62 0.68 192.84 14.26 淀粉积累期 VIs+IWD 0.66 205.05 16.53 0.71 178.37 13.09 VIs+CWT 0.72 181.28 14.27 0.78 167.70 12.79

表6 基于不同变量的马铃薯AGB估算精度对比 Table 6 Comparison accuracy of potato AGB estimation based on different variables

将6种植被指数、 3种高频信息和10种小波系数分别与AGB进行相关性分析, 结果表明, 随着生育期推进, 3种类型的输入变量与AGB的相关性均呈现先升高后降低的变化趋势, 这是因为马铃薯从现蕾期刚刚开始生长, 主要为营养生长和生殖生长, 随着时间推移, 茎节伸长, 叶片扩大, 到了块茎增长期, 茎节和叶片等营养器官生长最为旺盛, 此时植被覆盖度最大, 随后由于地上部有机物需要不断向地下块茎输送, 外加连续多天大雨, 造成地上部叶片迅速枯黄脱落, 淀粉积累期的植被覆盖度明显降低, 使得上述3种变量与AGB的联系由好变差。 每个生育期提取的高频信息和小波系数与AGB相关性都高于植被指数, 说明通过IWD和CWT技术能够排除土壤背景的影响, 使得提取出的作物冠层光谱信息能够与AGB充分联系, 解决作物生长后期植被指数存在的饱和性问题, 提高光谱信息的敏感度; 然而小波系数与AGB相关性更高, 这是因为高光谱传感器所含波段数目较多, 利用CWT技术能够细化光谱信息, 深度挖掘潜在的敏感光谱信息, 凸显出更多与AGB相关的有效信息, 使得筛选出的敏感波段下的小波能量系数与AGB的相关性更高。

本研究每个生育期以VIs结合CWT技术估算效果最好、 VIs结合IWD技术次之、 仅以VIs估算效果最差, 这是因为作物器官垂直生长, 冠层结构包含了较少的像元, 这对作物冠层光谱的贡献不大, 导致利用作物冠层光谱信息组建植被指数估算AGB存在困难, 而IWD技术结合数码影像提取高频信息以及CWT技术结合马铃薯高光谱影像提取小波能量系数, 不受高覆盖冠层结构的限制, 使得提取的影像信息能充分反演AGB。 但是, 成像光谱仪相当昂贵, 而数码相机却十分便宜, 以VIs+IWD为输入变量估算AGB, 不仅模型预测效果较植被指数的好, 而且节约了成本。 因此, 未来可以使用廉价的数码和多光谱相机以及低成本的无人机去提取作物冠层的高频信息结合植被指数来估算作物多生育期的AGB。 未来的研究, 应该注重使用其他的技术结合多种类型的作物冠层影像对不同年限和不同作物的AGB及其他参数进行遥感估算。