油菜是中国传统优势油料作物，总种植面积约占全国油料作物总面积的50%，是重要的食用油原料和油粕原料^[1]。对油菜进行快速准确的监测，不仅关系到保障食用油稳定供应，也关系到中国油料产品贸易与油料产业的发展。

传统的农作物产量预测方法主要有抽样、气象估产、农学估产和统计模型估产等，但在监测过程中，这些估产方式都表现出明显的局限性^[2]。迅速发展的遥感技术可提供对地宏观、动态、快速的观测，为作物长势的动态监测和产量估算提供一种全新的科学手段。作物遥感监测研究多使用卫星遥感数据，借助光谱植被指数等对作物长势进行监测。Gitelson等^[3]基于小麦冠层光谱特性，研究了不同植被指数与植被覆盖率的关系；Viña等^[4]利用植被指数对玉米与大豆的叶面积指数进行遥感估算；Becker-Reshef等^[5]通过时间序列归一化植被指数（normalized difference vegetation index, NDVI）数据对堪萨斯州小麦产量进行估算，并在乌克兰地区实现提前预测；秦占飞等^[6]分析了不同波段组合，构建了4种植被指数与水稻叶面积指数反演模型；王蕾等^[7]利用条件植被温度指数预测了县域尺度的小麦单产；江东等^[8]以中国河南省为例，分析了不同地区冬小麦高、中、低产县NDVI时间线曲线和不同生长阶段的关系。有研究^[9]认为，NDVI与产量相关，进行遥感估产最普遍的方式是在NDVI与作物产量之间建立经验模型。广泛应用的产量与遥感光谱指数的单纯统计模式操作简单，便于推广。然而，传统植被遥感中的NDVI无法体现农作物的具体植株器官形态及发育，忽略了不同器官在光合作用、淀粉累计等产量形成环节中的差异，无法满足农作物精准分析的客观需求。另外，从产量构成角度来看，物理意义不够明确也使植被指数模型在跨地区与跨年的应用效果不佳。因此，从作物成像和产量形成机理上对油菜进行定量研究尤为迫切。

遥感影像中的光谱信息可以反映作物生长情况，不同时期各端元丰度数据可反映作物的物候发育、物质生产与分配，与作物最终产量紧密相联。由于混合像元的广泛存在，丰度数据可以改善植被指数单纯基于像元反射率的估产方式，能从亚像元层次与植株器官发育的角度对NDVI估产模型进行修正。本文基于无人机遥感的方式，从单一时期丰度数据与产量构成关系分析出发，紧密结合油菜自身生长的农学特性，提出了新型的利用多时相丰度信息与NDVI相结合的分析方法，建立了适合油菜种植实际特点的估产模型。

1 实验数据与采集方法

本文实验数据包括无人机影像数据、地物端元光谱数据和油菜产量数据。实验田内共有48块实验小区，其中24块小区直播种植，另24块小区移栽种植。无人机影像数据包括40 m航高的油菜叶期影像、50 m航高的花期影像以及50 m航高的角果期影像，使用的Mini-MCA（multiple camera array）6波段相机的波段范围分别为485~495 nm、545~555 nm、665~675 nm、715~725 nm、780~820 nm、890~910 nm。对3个时期无人机影像的6个波段分别进行几何纠正和辐射定标。几何纠正时，以第5波段作为基准波段进行配准。经过几何纠正后，使用与Mini-MCA传感器配套的标准定标毯进行场地定标，得到地表反射率数据。地物端元光谱使用ASD（analytical spectral devices）光谱仪采集，探测波段范围为350~2 500 nm。采集时，光谱仪探头垂直向下，距地物0.2 m左右，每次记录5个光谱作为一采样光谱，选定位置分布均匀的3~5个点，将其光谱平均值作为该地物端元的光谱反射值。由于采集到的ASD光谱数据在对应MCA相机6个波段范围内的波动较小，对每个端元参照无人机影像的6个波段范围，分别将3个时期的地物光谱数据进行均值处理，得到平均反射率数据。

2 混合光谱分解

油菜作物结构复杂，传感器成像时，单一像元内往往包含不同的地物覆盖成分，混合光谱现象广泛存在。为准确获取目标地物的光谱信息，提高遥感影像的应用精度，必须进行混合光谱分解^[10]。具体步骤分为端元光谱获取和混合像元分解。端元获取的常见方式有图像端元提取法、物理端元提取法、光谱库获取法等。本文使用物理端元提取法，通过野外采集地物光谱曲线的方式获取端元光谱。常见的混合光谱分解模型有线性模型、概率模型、随机几何模型等。其中，线性模型使用最为广泛，原理也较为简单，一般形式如下：

$\mathit{\boldsymbol{P}} = \sum\limits_{i = 1}^N {{c_i}{e_i} + e = } Ec + n$

(1)

式中，P为像元的波谱向量；N为地物的端元数量；E为在每个波段范围内端元的反射向量；c为端元所占比例系数向量；c_i为第i个端元所占像元的面积比；e_i为第i个端元的反射率；e和n为误差项。在实际光谱解混中，线性模型会受到两个约束条件限制——非负性约束和归一化约束^[11]，即：

$\sum\limits_{i = 1}^N {{c_i} = 1,0 \le {c_i} \le 1} $

(2)

3 多时相分析模型

油菜产量的形成需要经过漫长的生育期，不同时期植株器官发育均影响着作物最终经济产量。本文通过探讨不同时期各端元丰度与产量的相关性，将NDVI估产模式置于产量形成过程的具体生物分析基础上，建立结合油菜生育期特征的估产模型。

3.1 模型基本流程

为了探究油菜在不同生长期植株器官发育对产量的影响，分别将叶期、花期、角果期经解混得到的端元丰度数据与产量进行相关分析，探寻对产量最为敏感的要素。在此基础上，进行多时相相关分析并建立多时相估产模型，如图 1所示。首先对无人机影像进行几何纠正和辐射定标处理，得到作物地表反射率影像；然后根据无人机波段范围，对采集到的地物端元光谱数据进行均值处理，得到端元光谱矩阵；接着基于全约束混合光谱线性解混模型对地表反射率影像进行解混，得到端元丰度数据；依据丰度与产量的相关性分析，得到产量最为敏感的端元；最后对比结合NDVI估产模式，通过多时相相关性分析，建立多时相估产模型。

决定油菜籽粒产量的主要因素有单位面积角果数、角粒数和千粒重。最终产量是干物质的积累，约90%~95%产量来源于光合产物。产量的形成主要分为3个阶段：从出苗到花芽分化前，叶片同化作用产生的干物质供植株建成和后续花芽分化，该阶段为角粒数和单位面积角果数的积淀过程；初花期到终花期是角粒数和单位面积角果数的定型过程；在角果期，角果层的成熟过程决定了千粒重大小^[12]。由于作物每一生长时段周围环境的可变性，只通过作物某一生长阶段的瞬时信息预测产量，往往考虑得不够全面。为较为准确地表征这种积累过程，使用多元线性回归的方式结合3个生长期的数据建立多时相估产模型。多元线性回归模型的一般形式为：

$ y = {\beta _0} + {\beta _1}{x_1} + {\beta _2}{x_2} + \ldots + {\beta _i}{x_i} $

(3)

式中，y为因变量；x为自变量；i为自变量个数；β₀与β_i分别是常数项与回归系数。建立多时相线性模型时，自变量的选取需考虑油菜不同生长期中有助于最终产量形成的部分。

3.2 评价指标

皮尔逊相关系数常用于度量两个变量X和Y之间的相关性（线性相关），对单一时期敏感要素与产量进行相关分析时，本文使用皮尔逊相关系数（ρ_X，Y）进行衡量：

$ {\rho _{X,Y}} = \frac{{ \sum \nolimits^ ({X_i} - \bar X)\left( {{Y_i} - \bar Y} \right)}}{{\sqrt[{}]{{ \sum \nolimits^ {{({X_i} - \bar X)}^2}}}\sqrt[{}]{{ \sum \nolimits^ {{({Y_i} - \bar Y)}^2}}}}} $

(4)

式中，X_i、Y_i分别为第i个样本的实际值；X与Y分别为X和Y的平均值。通常情况下，通过表 1所示的取值范围来判断相关强度。

由于实际种植条件的限制，一般难以获取大量的产量数据。针对本文实验样本数量较少的情况（一组24个），使用留一交叉验证的方式来建立最终估产模型。留一交叉验证是最为常见的交叉验证方法之一，广泛应用于模型建立和验证，能够充分利用实验数据^[13]。其基本思想是每次从n个样本中选取1个作为验证集，剩下的n－1个样本作为训练集进行训练。训练过程重复n次，将n次训练过程的结果取平均值作为泛化误差的估计，得到最终模型的系数和精度，具体计算公式如下：

$C = \frac{1}{n} \sum\limits_{i = 1} ^n {C_i}$

(5)

${\overline {R{\rm{}}} ^2} = 1 - \left( {1 - {R^2}} \right)\frac{{n - 1}}{{n - p - 1}}$

(6)

${\rm{RMSE}} = \sqrt[{}]{{\frac{1}{n} \sum\limits_{i = 1} ^n {{({X_{{\rm{obs}},j}} - {X_{{\rm{model}},j}})}^2}}}$

(7)

式中，n为样本个数；C与$C_i$分别为最终模型系数与第i次训练集训练得到的模型系数；$R^2$与R²分别为判定系数与调整后的判定系数；$p$为自变量个数；RMSE（root mean squared error）表示均方根误差；X_{obs, j}与X_{model, j}分别为第j个样本的样本实际值和模型预测值。

4 油菜长势监测实验分析

由于移栽油菜排列较整齐，植株间距较大，而直播油菜生长较为密集，故两种种植方式下油菜区域结构有所差别，需分别进行敏感要素分析与建模，在此基础上分析二者的共性规律，进而得到一般性估产方法。

4.1 丰度与长势相关性分析 4.1.1 叶期端元丰度

在油菜叶期，实验田内主要为叶子和土壤。从油菜植株的生长发育角度来看，叶片作为油菜光合作用和蒸腾作用的主要场所，其同化作用过程中干物质的积累对油菜植株、角果以及籽粒的生长发育有很大影响。叶期叶片分为长柄叶和短柄叶，其中长柄叶主要供应根系养分；短柄叶养分供应兼顾上下，是春壮稳长所调控的一组叶片^[14]。随着油菜的生长，不同类型叶片净光合速率的变化存在差异。

考虑到本文实验土壤环境与叶片的发育情况，选取5种地物端元进行光谱采集，分别为长柄叶、短柄叶、干土、半干土、湿土。基于全约束混合像元线性解混的方法计算得到叶期移栽和叶期直播各端元的丰度图（见图 2、图 3），最终叶丰度为长柄叶与短柄叶两种端元丰度之和。对叶端元丰度与产量进行相关分析可知，移栽油菜与产量相关系数为0.889，直播油菜与产量相关系数为0.818。两种种植方式下油菜叶丰度与产量均为极强相关，相关系数均高于0.8，能很好地反映叶丰度与产量之间的联系。

4.1.2 花期端元丰度

在开花后，油菜花序会逐渐遮挡植株茎秆，叶片光合作用会逐渐减弱，光合产物主要提供给茎干生长、花序分化。从产量积累过程来看，叶片仍然是光合作用的主要场所。花期叶片包括无柄叶、短柄叶和长柄叶，其中短柄叶的光合贡献更显著。由于干物质逐渐向花转移，油菜叶子会开始变黄，此时，实验田内地物组成较为复杂，不同地物的光谱反射曲线如图 4所示。

从图 4可以看到，花端元的光谱曲线走势较为接近；不同类型的叶片在绿波段的反射率有着较大差异；土壤和油菜植株的反射率差别很大。由于无人机相机只有6个波段，满足线性解混模型条件下求解丰度时最多只能选取6种端元的要求^[15]。选取上层花、无柄叶、短柄叶、干土、半干土、湿土6种端元进行解混，得到各端元丰度图（见图 5、图 6）。分别将花丰度和叶片总丰度与产量进行皮尔逊相关分析，统计相关系数，见表 2。

从表 2可以看出，不同种植方式下，不同端元与产量相关性不同。对移栽油菜而言，叶丰度与产量的相关性极强，花丰度与产量的相关性中等，而直播油菜则相反。值得注意的是，直播油菜叶丰度与产量是负相关的，这是由于直播油菜植株的单株产能不高，但最终亩产密度大。

4.1.3 角果期端元丰度

角果期是油菜生殖生长阶段进行养分积累和物质合成的关键时期。终花后，角果和叶片一同进行光合作用；角果层形成后，由于对叶片的严重遮盖，叶片开始衰老和掉落，短柄叶和无柄叶表面积减少。角果皮产生的同化物转移分配给籽粒，供籽粒充分灌浆^[16]。随着角果鲜重逐渐加大，果序轴被压弯，导致角果下坠或重叠在一起^[17]。根据角果期实验田地物组成，选取无柄叶、短柄叶、上层角果、下层角果、干土、湿土6种端元，并解算丰度图（见图 7、图 8）。分别将角果丰度与叶片丰度与产量进行回归分析，统计相关系数，如表 3所示。可以看到，移栽油菜的叶丰度与产量仍保持极强的相关性，角果丰度与产量相关性较弱；而直播油菜的情况较为复杂，叶丰度和产量相关系数为0.531，角果丰度和产量相关系数为0.766，两者差异不明显。

4.2 多时相模型

油菜籽粒产量的高低和品质好坏主要受油菜植株光合作用的调控，大约90%~95%的产量来自光合产物^[18]。根据产量形成过程与单一时期丰度数据分析可知，在油菜生殖生长过程中，不同器官发育状况都影响着最终产量。相比仅用一个自变量估测产量，采用多个自变量的最优组合共同估测因变量的方式更符合实际，也更有效^[19]。大量估产研究表明，NDVI是反映植被生长状态及植被覆盖度的最佳因子^[20]。为更好地评价多时相丰度数据综合估产的效果，分别选取以下两种不同的自变量方案进行综合回归分析。

1）方案1：将3个生长期NDVI数据作为自变量建立估产模型，模型精度如表 4所示。

2）方案2：分别在3个时期内各选取一种端元丰度数据作为自变量建立估产模型。根据单时期端元丰度与产量相关分析可知，对于移栽油菜，3个生长期的叶丰度与产量均极强相关，且明显强于其他端元丰度；而直播油菜情况较为复杂，没有在多个时期保持与产量相关性突出的端元丰度。因此，针对移栽油菜，方案2分别选择3个时期的叶丰度作为自变量；而对于直播油菜，尝试4种自变量组合（见表 5）。使用不同的自变量组合时，模型精度如表 6所示。

从表 6可以看到，方案2中不同丰度组合下模型调整后测定的系数均在0.74以上（A组合在0.8以上），说明从植株发育与干物质的积累方面预测最终产量是可行的。移栽油菜的方案2相比方案1拥有更好的拟合效果；而直播油菜4种端元丰度组合方式的效果均略逊于方案1，其中A组合表现最佳，这是由于直播植株生长密集、冠层结构复杂造成的。方案1主要通过增强近红外与红光通道反射率的对比度来监测植株的生长情况；而方案2从植株器官发育的角度出发，利用不同端元丰度数据反映干物质的积累与转移过程，进而估测产量，解混得到的丰度数据也可在亚像元层次反映植株的生长情况。

3）为结合光谱指数和丰度信息的优势，本文提出方案3，结合方案2中提出的不同端元丰度组合方式，将3个生长期的3种端元丰度与对应时期的NDVI相乘，乘积作为自变量与产量进行建模。模型最终评价指标统计结果如表 7所示。

对比方案1中的R²和RMSE，方案3中的A丰度组合方式可明显提高估产精度，这说明丰度数据可以改善植被指数单纯基于像元反射率的估产效果。由于移栽油菜冠层结构较直播油菜简单，方案3的估产精度提高效果更为明显，虽然不同种植方式的精度提高程度有所差异，但两种种植条件下最优丰度自变量组合选取一致，均为A组合。这说明在选取合适的端元丰度组合时，丰度对NDVI均有一定的修正效果，且在不同种植方式下均适用。在实际生产中, 中国油菜多采用移栽种植，而国外部分地区则采用直播种植。表 8列出了选取A组合时不同种植方式下方案3估产模型的回归公式。

在实际应用中，产量往往需要提前预测，以便采取及时有效的增产措施。仅选取叶期和花期数据分别按照上述3种自变量方案（其中方案2、方案3的自变量均使用A丰度组合）建立多时相估产模型，模型回归统计结果如表 9所示。

由表 9可知，仅用前两个生育期数据（叶期和花期）的模型效果略逊于三时期综合估产，但方案3的R²仍在0.85以上，能较准确地反映最终产量。对比方案1，不同种植方式下A丰度组合的方案3估产效果仍保持提高，说明丰度对NDVI的修正效果是稳定且可靠的。表 10给出了提前估产时方案3的估产模型公式（A丰度组合）。

5 结语

本文从油菜产量形成的农学意义出发，对油菜3个典型生长期进行长势监测分析，综合考虑反映油菜整体生长状态的NDVI和体现植株器官构成的丰度信息，得到了能进行产量预估的无人机遥感模型。实验结果表明：

1）在不同的种植方式下，不同器官端元丰度与产量的相关规律有所差异，移栽油菜叶丰度在3个生长期与产量均有极强的相关性，而直播油菜由于冠层结构复杂，植株密度大，没有明显突出的强相关端元丰度。

2）从干物质的生产分配的角度出发，利用多时相丰度数据建立的估产模型（方案2）也是可行和有效的，但精度略低于NDVI建立的估产模型。

3）在选取合适的端元丰度组合时，丰度数据与NDVI结合进行综合估产的模型精度最高，说明丰度可从产量构成及亚像元层次改善NDVI的估产效果。虽然不同种植方式的提高程度有所差异，但两种种植条件下最优丰度自变量组合的选取是一致的（均为A组合），说明丰度对NDVI的修正效果是稳定且可靠的，适用于产量的提前预测。

由于油菜植株是立体结构，使用线性混合光谱模型进行分解时不能完全描述地物情况。后续研究还需考虑不同种植条件下植株的生长差异，从更严密的角度进行混合光谱分析，建立更为一般化的估产模型。