融合无人机影像光谱与纹理特征的冬小麦氮营养指数估算

杨福芹，冯海宽，肖天豪，李天驰，郭向前

（1. 河南工程学院土木工程学院，河南 郑州 451191；2. 国家农业信息化工程技术研究中心，北京 100097；3. 河南省测绘工程院，河南 郑州 450003；4. 河南省地质矿产勘查开发局测绘地理信息院，河南 郑州450006）

“有收无收在于水，收多收少在于肥”表明了矿质营养在农业生产中的重要性，而在作物生长中所需的七种主要矿质元素（氮、钾、钙、镁、磷、硫、硅）中，氮素所占比例超过40%，氮素是合成小麦子粒中氨基酸和蛋白质的主要成分，也是叶绿素、植物激素及维生素的主要成分[1]。当氮肥供应不足时，植株表现为矮小、瘦弱、直立，分蘖少或无，穗小粒少；当氮肥供应过量时，植株表现为营养生长旺盛，植株高大细长，开花少，秕粒多，易倒伏，贪青晚熟等现象，导致小麦的质量及产量显著下降[2]。此外，过多使用氮肥还会造成农业资源的巨大浪费，直接或间接地危害人畜健康[3]。因此，及时掌握作物的氮素营养状况，根据其需求进行氮肥管理，对节约成本，减少环境污染具有重要意义[1-4]。

传统小麦氮素测定的方法主要有田间调查、取样与室内样品生化分析，但这些方法具有滞后性、破坏性和非动态性，且成本较高；而且由于受理化分析条件限制，只能选取植株的叶茎进行分析，“以点带面”与实际作物氮素营养状况存在很大差异，难以精确的描述植株整体氮素水平。近年来，随着无人机遥感技术的不断发展，利用无人机搭载高光谱相机、多光谱相机、数码相机等传感器获取遥感影像，利用光谱信息或影像纹理特征在作物生态参数反演、农田生态环境信息检测等方面得到广泛应用[5-9]，如基于雷达数据、图像数据及实测数据估算草地生物量[10]、获取多光谱影像数据用于反演大豆叶面积[11]。研究发现，影像纹理特征估算作物参数（生物量、叶绿素、氮含量）的能力高于单一的植被指数，但在作物氮营养指数方面缺乏探讨。本研究将图像指数和纹理特征进行融合形成“图-谱”融合指标，探究“图-谱”融合指标估算氮素营养诊断的能力，以实现融合光谱信息和纹理特征的冬小麦氮素营养诊断，从而更好的监测作物长势，为精准氮肥管理提供一种技术手段。

1 材料与方法

1.1 试验设计

试验于2014—2015年在北京市昌平区小汤山镇东北部国家精准农业研究示范基地进行，基地处于40°00′N～40°21′N，116°34′E～117°00′E，海拔高度36 m。前茬作物为玉米，试验土壤类型为潮土，土壤中0～30 cm土层中硝态氮含量为3.16～14.82 mg/kg，全氮含量1.0～1.2 g/kg，有机质含量为15.8～20.0 g/kg，有效磷含量为3.14～21.18 mg/kg，速效钾含量为86.83～

120.62 mg/kg。灌溉方式采用喷灌灌溉，传统种植方式。采用品种、氮肥和水分的正交试验。试验设计为2个品种、4个氮素水平和3个水分水平。2个品种分别为京9843（J9843）和中麦175（ZM175）；4个氮素水平分别为未施尿素（0 kg/hm2，N1）、195 kg/hm2尿素（1/2正常，N2）、390 kg/hm2尿素（正常，N3）、585 kg/hm2尿素（3/2正常，N4）；3个水分水平分别为雨养（W1）、正常水（W2）、2倍正常水（W3）。试验田东西向长度为84 m，南北向长度为32 m，共48个小区，每个小区面积为48 m2，16个小区为一组，重复3次，重复1和重复2用于建模，重复3用于验证，如图1所示。

1.2 数据获取

1.2.1 地上部生物量及氮含量测定 试验于冬小麦开花期（2015年5月13日）进行，主要获取冬小麦植株氮含量及地上部生物量等数据。在室外调查固定样方0.3 m×1 m的每个小区内，破坏性选取具有代表性的20株小麦作为样本，器官分离后分别置于纸袋中，置于烘箱105 ℃杀青30 min，然后设置温度到75 ℃烘干至恒重（48 h以上），将烘干后的各器官进行称重。将测定过的地上部生物量样品粉碎，采用凯式定氮法测定植株不同器官的氮含量。

1.2.2 氮营养指数 氮营养指数（nitrogen nutrition index，NNI）描述为作物实测植株氮浓度与临界氮浓度的比值[12]，计算公式为：

式中：N为实测植株氮浓度，g/100 g，Nct为临界氮浓度。根据Lemaire等[13]的定义：临界氮浓度（Nct）为作物地上生物量达到最佳生长速度所需要的最低氮浓度，计算公式为：

图1 试验设计Fig. 1 Experimental design

式中：a表示植株地上部生物量为1 t/hm2时的临界氮浓度，取值为5.35，b为决定临界氮浓度稀释曲线斜率的参数，取值为-0.442，DM为地上部生物量，t/hm2。

1.2.3 无人机数码影像获取及处理 无人机影像采用八旋翼电动无人机，搭载DSC-QX100数码相机获得。DSC-QX100数码相机包含红光波段（700 nm）、绿光波段（550 nm）和蓝光波段（470 nm）。无人机选择天空晴朗无云，中午12:00-13:00对冬小麦试验田进行航拍，飞行高度为50 m，地面分辨率为0.013 dpi。将获取的数码影像数据采用俄罗斯Agisoft LLC公 司 开 发 的Agisoft PhotoScan Professional软件处理，主要步骤为：1）数码影像和POS点数据输入，2）点云生成，3）网格构建，4）纹理生成，5）生成无人机高清数字高程模型（digital elevation model，DEM）和数字正射影像（digital orthophoto map，DOM）。最后获取的无人机影像如图2所示。

图2 无人机数码影像Fig. 2 Digital image of UAV

1.3 研究方法

1.3.1 图像指数的选取 利用无人机数码影像可见光波段构建的图像指数可以较好的反映作物的氮素营养状况[14]。本研究通过获得的DOM高清正射影像，利用ENVI软件获取每个试验小区图像的平均红光值R，平均绿光值G和平均蓝光值B。根据图像R、G、B灰度平均值计算了3个归一化特征参数，分别为归一化红光r、归一化绿光g和归一化蓝光b，其公式分别为r=R/R+G+B、g=G/R+G+B和b=B/R+G+B。根据3个归一化特征参数r、g和b，在借鉴前人研究的基础上，选取了能反映作物氮素营养状况的16个可见光图像指数：红蓝比值指数[15]、绿蓝比值指数[15]、红蓝差值指数[15]、红蓝和指数[15]、绿蓝差值指数[15]、红蓝植被指数[15]、三波段植被指数[15]、超绿指数EXG[16]、沃贝克指数WI[16]、红绿植被指数GRVI[17]、修正红绿植被指数MGRVI[17]、红绿蓝植被指数RGBVI[17]、超红指数EXR[18]、归一化差异植被指数NDI[19]、大气阻抗植被指数VARI[20]和超绿超红差分指数EXGR[20]。

1.3.2 影像纹理特征的提取 纹理特征是图像灰度等级的变化，不仅可以表示图像的均匀、细致、粗糙等现象，而且可以揭示图像中地物与其周围环境的关系，是遥感影像的重要特征[21]。研究表明，影像纹理特征在作物长势监测中取得了较好的试验结果[22]。本研究利用ENVI中的灰度共生矩阵（Gray Level Cooccurrence Matrix，GLCM）对可见光波段进行0°、45°、90°、135°4个方向8个纹理特征的提取，对不同方向的纹理特征进行平均，得到各波段的8个纹理特征。每个波段的8个纹理特征值分别为：均值（mean）、方差（variance，var）、同质性（homogenetity，hom）、对比度（contrast，con）、差异性（dissimilarity，dis）、熵（entropy，ent）、二阶距（second moment，sm）和相关性（correlation，cor）[22]。其中，mean_R、var_R、hom_R、con_R、dis_R、ent_R、sm_R、cor_R表示红波段对应的纹理特征；mean_G、var_G、hom_G、con_G、dis_G、ent_G、sm_G、cor_G表示绿波段对应的纹理特征；mean_B、var_B、hom_B、con_B、dis_B、ent_B、sm_B、cor_B表示蓝波段对应的纹理特征。

1.3.3 灰色关联分析 灰色关联分析（grey relation analysis，GRA）作为一种灰色系统分析方法，适用于研究因变量受到其他多个因素影响的强弱关系[23]。分析方法包括以下步骤：首先将氮营养指数视为参考序列，图像指数视为比较序列；随后对参考序列和比较序列进行无量纲化处理；最后计算灰色关联度。

1.3.4 多重共线性分析 对于图像指数之间存在的多重共线性，采用方差膨胀因子（variance inflation factor，VIF）进行衡量。当图像指数之间的方差膨胀因子值较大时，会使回归模型出现较大估算误差[24]。其计算公式如下：

式中：Ri2表示第i个图像指数与其他图像指数之间的决定系数。通常情况下当VIF＜10时可以视为图像指数之间不存在多重共线性；当10≤VIF≤20时，图像指数之间存在一定的多重共线性；当VIF＞20时，图像指数之间存在严重的多重共线性。

1.3.5 偏最小二乘回归 偏最小二乘回归法（partial least squares regression，PLSR）与传统多元线性回归分析法相比，能在自变量之间存在较严重的多重共线性条件下进行回归计算，同时该方法还包含了典型相关分析与主成分分析法[25]。

1.3.6 模型精度评定 为了评价模型的效果，选取决定系数（coeffcient of determination,R2）和均方根误差（root mean square error, RMSE）作为结果精度评价指标。计算公式为：

式中：n表示为样本总数，Xi、Yi分别表示第i个样本氮营养指数的实测值与估算值表示样本氮营养指数实测平均值与预测平均值。

2 结果与分析

2.1 图像指数与氮营养指数的相关性

表1为图像指数与冬小麦氮营养指数相关性的分析结果。从表中可以看出，除了g/b，g-b和WI图像指数外，其他图像指数与氮营养指数均达到0.01显著水平（P＜0.01），与冬小麦氮营养指数相关性的绝对值在0.371～0.773之间。对于3个归一化特征参数，图像指数与冬小麦氮营养指数相关性绝对值最大的是归一化蓝光b，其相关性为0.773。对于其他图像指数，与氮营养指数相关性绝对值最大的是归一化差值植被指数(r-b)/(r+b)，值为0.794，与氮营养指数相关性最差的是沃贝克指数WI，其相关性仅为0.036。为了更好的与冬小麦影像纹理特征融合，选择相关性较好的(r-b)/(r+b)、r/b、r-b、b、(r-g-b）/(r+g)、r+b、r、VARI、EXR、NDI、MGRVI图像指数与影像纹理特征进行融合。

2.2 纹理特征与氮营养指数的相关性

表2为纹理特征与冬小麦氮营养指数相关性的分析结果。从表中可以看出，除了contrast和mean纹理特征与氮营养指数相关性未达到0.01显著水平（P＜0.01）外，其他纹理特征与氮营养指数的相关性都达到了0.05显著水平（P＜0.05）。在红波段和绿波段，8个纹理特征与氮营养指数相关性绝对值的大小顺序都相同，依次为cor_R＞var_R＞dis_R＞hom_R＞sm_R＞ent_R＞mean_R和cor_G＞var_G＞dis_G＞hom_G＞sm_G＞ent_G＞mean_G；在 蓝波段，8个纹理特征与氮营养指数相关性绝对值的大小顺序依次为cor_B＞var_B＞dis_B＞hom_B＞ent_B＞sm_B ＞mean_B，仅红波段和绿波段与蓝波段纹理特征的熵和二阶矩有很小的变动，其他相关性顺序保持不变。根据纹理特征与氮营养指数相关性的优劣，筛选出var_R、cor_R、var_G、dis_G、cor_G、var_B、cor_B与图像指数进行融合。

表1 图像指数与氮营养指数的相关性Table 1 Correlation coeffcient between image index and nitrogen nutrition index

表2 纹理特征与氮营养指数的相关性Table 2 Correlation coeffcient between texture features and nitrogen nutrition index

2.3 “图-谱”融合指标与氮营养指数的相关性

将图像指数与纹理特征进行相乘或相除，构建既有光谱信息又有纹理特征的“图-谱”融合指标，探究“图-谱”融合指标反演氮营养指数的能力。将“图-谱”融合指标与氮营养指数进行相关性分析，根据相关性的优劣筛选出14个“图-谱”融合指标（表3）。从表中可以看出，“图-谱”融合指标与氮营养指数的相关性均达到0.01显著水平（P＜0.01），其相关性绝对值取值在0.5～0.9之间。其中相关性最好的“图-谱”融合指标是cor_R*((r-b)/(r+b))，其绝对值为0.819，相关性最差的“图-谱”融合指标是var_B/(1.4r-g)，其值为0.537。

表3 融合指标与氮营养指数的相关性Table 3 Correlation coeffcient between fusion index and nitrogen nutrition index

2.4 氮营养指数估算模型变量选择

为进一步提高模型精度，对14种“图-谱”融合指标进行灰色关联分析，探究融合指标与氮营养指数的贴近程度，结果如表4所示。可以看出，“图-谱”融合指标中灰色关联度的大小顺序为cor_B/(r-b)＞ cor_R/(r-b)＞ dis_G/(1.4r-g)＞ cor_G/((r-g-b）/(r+g))＞ var_G/(1.4r-g)＞ var_R/1.4r-g＞ var_B/(1.4r-g)＞ dis_G×(r-b)＞ var_G×((g-r)/(g+r-b))＞ var_B×((g-r)/(g+r-b))＞ var_R×(g-r)/(g+r-b)＞ cor_R×((r-b)/(r+b))＞ cor_B×((r-b)/(r+b))＞ cor_G×((r-b)/(r+b))，关联度最好的“图-谱”融合指标是cor_B/(r-b)，其值为0.865。为防止偏最小二乘回归模型的入选参量之间的高相关，采用方差膨胀因子对选取的“图-谱”融合指标进行多重共线性分析，结果如图3所示。可以看出，当0＜VIF＜10时，“图-谱”融合指标间不存在多重共线性，如var_R/1.4r-g与cor_R×((r-b)/(r+b))，cor_R/(r-b)与var_R/1.4r-g，cor_R/(r-b)与cor_R×((r-b)/(r+b))；当10≤VIF≤20时，“图-谱”融合指标间存在一定的多重共线性，如dis_G/(1.4r-g)与var_G×((g-r)/(g+r-b))，cor_G/((r-g-b）/(r+g))与cor_B×((r-b)/(r+b))；当VIF＞20时，“图-谱”融合指标间存在严重的多重共线性，如var_R/1.4r-g与var_R×(g-r)/(g+r-b)，var_G×((g-r)/(g+r-b))与var_R×(g-r)/(g+r-b)，var_G/(1.4r-g)与var_R×(g-r)/(g+r-b)。结合表4和图3，选取cor_B/(r-b)、dis_G/(1.4r-g)、cor_G/((r-g-b)/(r+g))、var_G/(1.4r-g)和dis_G×(r-b)作为多元变量入选偏最小二乘回归估算模型。

表4 融合指标与氮营养指数的灰色关联分析Table 4 Gray Correlation degree and rank of evaluation index for nitrogen nutrition index

图3 融合指标间的方差膨胀因子Fig. 3 Variance inflation factor among fusion index

2.5 氮营养指数反演建模与验证

2.5.1 基于图像指数的氮营养指数反演 采用偏最小二乘回归算法，基于5个图像指数建立的氮营养指数回归模型见式（6）。

式中：x1表示修正红绿植被指数MGRVI，x2表示超红指数EXR，x3表示归一化差异植被指数NDI，x4表示大气阻抗植被指数VARI，x5表示归一化蓝光b。冬小麦氮营养指数的反演结果如图4所示。从图中可以看出，冬小麦图像指数反演氮营养指数的建模结果为决定系数0.593 8，均方根误差为0.104 1（图4-A）。实测值与预测值分布在1∶1线附近，建模效果较好。验证结果的均方根误差为0.227 0（图4-B）。实测值与预测值大多分布在1∶1线之下，一部分预测值被严重低估，实测值与预测值之间存在较大误差。结果表明用图像指数反演氮营养指数虽然建模效果较好，但验证效果不佳。

2.5.2 基于纹理特征的氮营养指数反演 采用偏最小

二乘回归算法，基于5个纹理特征建立的氮营养指数回归模型见式（7）。

式中：y1表示var_R，y2表示cor_R，y3表示var_G，y4表示dis_G，y5表示cor_B。

冬小麦氮营养指数的反演结果如图5所示。可以看出，冬小麦纹理特征反演氮营养指数的建模结果为决定系数0.584 8，均方根误差为0.105 0（图5-A），实测值与预测值分布在1∶1线附近，建模效果较好。验证结果的均方根误差为0.115 3（图5-B），实测值与预测值同样分布在1∶1线附近，实测值与预测值之间的误差较小，验证效果较用图像指数反演氮营养指数效果好。

2.5.3 基于“图-谱”融合指标的氮营养指数反演采用偏最小二乘回归算法，基于5个“图-谱”融合指标建立的氮营养指数回归模型见式（8）。

式 中：V14表 示cor_B/(r-b)，V8表 示dis_G/(1.4rg)，V10表示cor_G/((r-g-b）/(r+g))，V6表示var_G/(1.4r-g)，V7表示dis_G×(r-b)。

冬小麦氮营养指数的反演结果如图6所示。可以看出，“图-谱”融合指标反演氮营养指数的建模结果决定系数为0.644 3，均方根误差为0.097 0（图6-A），实测值与预测值分布在1∶1线附近，建模效果较好。验证结果的均方根误差为0.114 0（图6-B），实测值与预测值同样分布在1∶1线附近，实测值与预测值之间的误差较小。研究发现，基于融合指标cor_B/(r-b)、dis_G/(1.4r-g)、cor_G/((r-g-b）/(r+g))、var_G/(1.4r-g)和dis_G×(r-b)所构建的氮营养指数模型，相比单一图像指数、单一纹理特征，两者融合反演的氮营养指数模型精度明显提高，这是因为氮营养指数反演模型中既含有丰富的纹理信息，又含有一定的光谱信息，该模型综合考虑了纹理特征及光谱信息对氮营养指数的贡献，相比图像指数、纹理特征构建的氮营养指数模型，融合图像指数与纹理特征共同反演氮营养指数的效果最好。

图4 基于图像指数的氮营养指数建模与验证Fig. 4 Relationship between calibration and validation results of nitrogen nutrition index

图5 基于纹理特征的氮营养指数建模与验证Fig. 5 Relationship between calibration and validation results of nitrogen nutrition index

图6 基于融合指标的氮营养指数建模与验证结果Fig. 6 Relationship between calibration and validation results of nitrogen nutrition index

3 讨论

氮营养指数可以直观地反映植株体内氮素的营养状况，若氮营养指数NNI＞1，表明作物植株体内氮素含量过高；若NNI=1，表明作物体内植株氮素含量达到最佳；若NNI＜1，表明植株体内氮素含量供应不足[26]。本研究利用图像指数预测冬小麦氮营养指数NNI的变化范围为0.06～1.03，在相同施氮水平下预测的氮营养指数最小值远远小于实测的氮营养指数最小值，预测的氮营养指数最大值较接近实测最大值，但依然偏小；利用影像纹理特征预测冬小麦氮营养指数NNI的变化范围为0.23～1.00，在相同施氮水平下预测的氮营养指数最小值仍小于实测的氮营养指数最小值，但效果好一点，预测的氮营养指数的最大值依然小于实测的氮营养指数实测值；利用“图-谱”融合指标预测的冬小麦氮营养指数NNI的变化范围为0.32～1.04，在相同施氮水平时比实测的氮营养指数的最大值和最小值都偏小，但都接近实测值。Liu等[27]基于无人机高光谱影像，预测的北京地区开花期冬小麦氮营养指数的范围在0.51～1.30之间，预测值与实测值的RMSE为0.074。王仁红等[28]利用地面非成像光谱仪获取高光谱数据，利用已有光谱指数构建北京地区冬小麦氮营养指数估算模型，其预测的氮营养指数范围在0.3～1.3之间，预测值与实测值的RMSE在

0.138～2.859之间。相对于无人机高光谱传感器和地面非成像高光谱传感器，本研究使用的是数码相机，只提供红绿蓝三个波段的光谱信息，用图像指数或纹理特征反演氮营养指数的效果远远不如前面所述。将图像指数与纹理特征融合后，氮营养指数反演模型的估测效果有了很大的提高，其效果与用高光谱数据估测的效果接近。此外，前人在“图-谱”信息综合反演作物参数时，多是利用多光谱影像或高光谱影像，而利用数码影像的比较少见。

当然，本研究也存在一些不足之处。本研究仅利用了冬小麦开花期的数码影像及地面实测数据，所使用的样本量较少，对构建的模型的通用性有一定的影响，未来的研究应增加不同生育期、不同年限、不同区域，拓展模型应用范围，提高模型的稳定性和普适性。另外，该模型对于不同作物能否进行氮素营养诊断并没有做探讨，这是将来研究方向。

4 结论

现有研究多是直接选取对氮营养指数敏感的光谱特征和纹理特征，基于多元回归模型反演氮营养指数。本研究从数码影像DN值、纹理特征以及两者融合入手，利用数码影像形成的图像指数与纹理特征采用相乘或相除的方式融合形成“图-谱”融合指标，探讨“图-谱”融合指标估算氮营养指数的能力。研究发现，基于融合指标所构建的氮营养指数模型，相比图像指数和影像纹理特征，其反演的模型精度明显提高，主要是因为融合指标同时考虑了数码影像光谱信息与纹理特征对氮营养指数的贡献，与前人[7,9,29]的相关研究结果一致。本文主要结论如下：

1）整合灰色关联度-方差膨胀因子，筛选出与氮营养指数灰色关联度较好，指标与指标间多重共线性较小的五个“图-谱”融合指标，分别为cor_B/(r-b)、dis_G/(1.4r-g)、cor_G/((r-g-b）/(r+g))、var_G/(1.4r-g)和dis_G×(r-b)。

2）相比单一图像指数、纹理特征，融合图像指数与纹理特征构建的“图-谱”融合指标反演冬小麦氮营养指数模型估算精度较高（R2=0.644 3），高于分别基于图像指数（R2=0.593 8）和纹理特征（R2=0.584 5）构建的氮营养指数模型。研究表明采用纹理特征和图像光谱信息相融合的方式进行建模，可以显著提高氮素营养状况的估算精度。