随着各地数字烟草3S技术的兴起, 烟草冠层和叶片信息的定量无损精确提取, 对烟草生长的科学管理和监测具有重要意义。罗静等^[1]指出虽烟草的生长受多种因素的影响, 其过程极其复杂, 但其生长状况仍可以用一些与该过程密切相关的因子进行表征, 如叶面积指数 (LAI)。LAI是表征植被生长状况的关键物理量之一, 其大小与作物生长状态直接密切相关。且LAI的变化是一个随时间变化的动态过程, 不同时间处的LAI值之间具有一定的相关性, 形成依赖关系^[2]。因此可通过观测LAI的变化获取和判断烟草生长状态的变化。

目前针对于烟草LAI计算的研究相对较少。烟草LAI测算代表性方法主要如下：(1) 因农作物LAI随生长期变化的过程基本符合LOGISTIC曲线或其修正形式, 汪耀富等^[3]建立烤烟叶面积指数LAI随栽种时间t的普适模型LOGISTIC模型, 发现相关系数达到97%以上；(2) 刘国顺等^[4]分析了17种光谱参数与烟草叶面积指数的关系, 通过建立逐步回归模型对烟草LAI进行估测；(3) 张正杨^[5]、王建伟^[6]等探究了RVI、NDVI等常用的植被指数和LAI的相关关系, 发现植被指数可精确反演烟草LAI, 其中与LAI相关性最强的常用植被指数为NDVI。

LOGISTIC模型是一个对观测LAI进行统计拟合的经验模型, 用数学模型表征LAI的变化趋势, 但其更多的是对生长过程的数学语言描述, 虽然这种描述是有价值的, 但得到的结果与实际存在一定的偏差, 很难精确模拟诸如突发性灾害等造成的作物生长状态剧烈不规律变化。另外作物模型的实际应用在初始值获取和参数区域化方面遇到很多问题；遥感观测是通过作物反射光谱特征获取作物群体瞬间生长状况, 遥感观测手段能够很好地及时反映作物生长的区域分布状况, 但易受遥感观测平台稳定性、地表异质性和遥感反演算法等因素的影响, 结果存在较大的不确定性。另外遥感技术是瞬时的, 易受大气条件、可行观测周期等因素限制而不能对作物进行无间断连续监测, 更无法对作物生长的相关状态进行预报。刘国顺^[4]用多光谱参数回归分析方法拟合烟草作物LAI时使用多达17种光谱参数, 多光谱参数参与可有效增加计算的精度, 但参数获取的可行性则限制了该方法的实际应用；另外LAI与光谱反射率存在着复杂的非线性关系, 张正杨^[5]则发现主成分分析和神经网络等方法反演精度相较于植被指数反演方法更高。

在作物长势的监测上, 遥感信息的实时性、宏观性与作物生长模型的连续性、机理性形成良好的互补性关系^[7], 可利用作物生长模型的连续性模拟弥补遥感观测的瞬时性, 利用遥感观测获取数据的及时性为作物生长模型提供参数, 而数据同化技术则可将拟合模型和观测数据两种信息有效结合起来。数据同化将新的观测数据引入到过程模型中, 不断减少或者滤掉过程模型的噪声, 使得模型轨迹模拟轨迹更加贴近自然界的真实状态^[8-9]。遥感获取的观测数据常常是瞬时的物理特征, 作物生长模型是一种面向过程、机理性的动态模型^[10], 利用数据同化技术将烟草LAI的多时相观测数据和作物生长动态变化模型二者信息相融合, 具有改善和提高烟草LAI数据精度的能力。

综合这两种信息的成熟数据同化方法通常有滤波算法或变分方法两种^[11], 而滤波算法因其实用性, 易于实现, 可移植性强的特点而被广泛采用。目前滤波方法尤其是集合卡尔曼滤波算法在农作物的应用方面已经趋于成熟。王东伟^[12]以2004年北京冬小麦为例, 对冠层观测反射率进行卡尔曼滤波同化, 初步实现了集合卡尔曼滤波算法在冬小麦LAI反演中的应用, 结果显示集合卡尔曼滤波算法能够使LAI同化结果接近真实值；黄健熙等^[13]选择河北衡水地区冬小麦为研究对象, 通过集合卡尔曼算法获取时间序列最优的LAI, 并依此估算区域内冬小麦产量, 结果表明同化后的冬小麦产量比未同化的产量预测精度有显著提高；陈思宁^[14]基于集合卡尔曼滤波对东北玉米产量及LAI进行估计, 对同化前后模拟LAI对比发现同化后其轨迹更接近实测值, 更符合玉米的生长发育趋势。以上研究均发现, 应用集合卡尔曼滤波同化方法于农作物LAI估算能够显著提高估算的精度。

遥感观测是一种快速获取烟草信息的手段, 模拟模型能够对烟草生长过程相对较好的整体拟合, 卡尔曼滤波同化技术则能够综合二者各自的优势, 有利于提高数据的预测精度。文中尝试以广东省南雄市烟草为例, 构建南雄市烟草LOGISTIC生长模型, 利用集合卡尔曼滤波数据同化技术综合遥感观测LAI数据和作物生长模型二者信息进行同化分析, 以获取精度相对较高的烟草LAI数据。

1 研究区域

文中以广东省韶关市南雄地区为实验区。南雄市位于广东省东北部, 113°55′30″—114°44′38″E, 24°56′59″—25°25′20″N, 该地区是广东省烟草主要种植县市之一。研究区内烟草类型主要为烤烟。烟草生长期分苗床和大田, 广东省烟草生长周期约为180d。苗床期指播种到移栽时间段, 即1—2月底 (一般为60d), 大田期指烟苗移栽到采收完毕时间段, 即3—6月 (一般为100—120d)。而4—6月广东省处于梅雨季节, 这3个月雨水天气居多。

2 数据与材料

实验中光谱数据主要是使用光谱仪在野外实测得到。实验使用荷兰Avantes公司AvaSpec-ULS2048光纤光谱仪烟草野外实地测量光谱反射率, 光谱范围取300—1100nm, 光谱采样区间为0.6 nm, 光谱分辨率为1.4 nm。选择晴朗无云无风天气, 于10:00—14:00时间段测定光谱反射率。测量时传感器探头垂直向下, 光谱仪视场角为15°, 合理控制光谱仪距离烟草植被冠层顶部垂直高度。测量时合理挑选样株, 严格使用标准白板校正, 每个样株进行10次光谱测量, 以其平均值作为该样株的光谱反射值。烟草光谱采样时间为2014年2—6月, 间隔基本为每半个月1次, 得到烟草的连续时间序列光谱数据。

LAI实测数据主要是在野外对烟草实地测量得到。实验组中实测LAI数据分为两组, 建模数据和实验数据, 建模LAI数据来源于多期对南雄市多个采样点的烟草实地测量, 用于建立适合于南雄的LOGISTIC模型；实验LAI数据主要位于南雄市黄坑镇许村 (25°13′N, 114°28′E) 附近多个采样点, 与烟草光谱测量工作同时进行, 对每个烟草样本测量光谱的同时对其测量叶面积, 作为验证数据。因为烟草的叶片为长椭圆形, 根据统计分析得出每个叶片大致为整个叶片最长与最宽之积面积的63.45%, 故本实验中叶面积的计算公式为：叶片长×叶片宽×0.6345。

实验中的天气温度数据取自2011年到2014年韶关市气象局发布的逐日温度统计数据, 日积温数据则取日平均温度。

3 研究方法 3.1 集合卡尔曼滤波 (EnKF) 方法

卡尔曼滤波是一种将模型模拟和外部观测结合起来估计目标参数的优化方法, 是一种用Monte Carlo的集合预报方法估计预报误差协方差的顺序同化算法, 它能够有效降低估算过程中的误差, 提高预测精度。最大的优点是不需要预报算子的切线性模式和伴随模式, 对于非线性很强、不连续的动态模型也有很好的模拟。

集合卡尔曼滤波包括预报和更新两个部分。在预报部分, 初始状态向量的集合通过过程模型模拟来获得预报场的集合, 用预报集合计算预报误差协方差矩阵；在更新部分, 利用观测向量和状态向量的误差协方差矩阵更新每个集合, 得到分析场的集合, 最后将分析场集合的均值作为模型状态的后验估计值^[11]。当前计算结果只与前一状态估计值、当前状态测量值有关, 并且易于实时处理, 无需存储大量的数据。

集合卡尔曼滤波通过不断迭代预测和更新两个过程完成数据的同化, 图 2表示了集合卡尔曼滤波算法一个完整的迭代过程。

在预测步骤, 将t时刻分析阶段得到的分析集合X_{i, t}^a(初始化集合输入X_{i, 0}^a) 作为初始场, 通过模型算子M(.) 向前预报至t+1观测资料时刻, 得到t+1时刻的预报集合X_{i, t+1}^f。更新步骤则是通过增益矩阵K_t+1对预报集合X_{i, t+1}^f和观测数据集合Y_{i, t+1}进行同化计算, 同化得到的状态值集合平均值X_{i, t+1}^a作为同化值储存, 而X_{i, t+1}^a则继续参与下一个预测与更新过程的计算。

3.2 LOGISTIC模型

烟草生长需要适宜的温度, 积温对烤烟大田生长发育有一定的影响, 烟草完成自己的生长周期, 需要一定的温度积累。如果生长期间的昼夜平均温度较低, 植株为满足自己所需要的温度总和, 会使生育期延长^[15]。在LOGISTIC模型中, 为消除烟草生长发育随品种、播期、地域的不同, 对模型的时间进行统一处理。用活动积温表示生育期长度, 并作归一化处理^[16], 使用归一化积温来统一时间尺度。积温AT与归一化积温DS的计算方法分别如下：

(1)

式中, i表示时间序列数, n为总时间总数, T_i为逐日平均温度, AT表示总积温值, DS_i为i时刻的归一化积温值。本实验只对大田期的烟草进行试验, 对DS取值范围为0—1。以DS为自变量, LAI为因变量建立简化的LOGISTIC模型^[3]得到：

(2)

LAI是叶面积指数, 这里LAI_max是生育期中LAI最大值, DS为归一化积温, A, B, C是常量参数。

4 研究结果 4.1 南雄市烟草LAI数据计算

本实验采用直接建立起南雄市烟草归一化植被指数 (NDVI) 与LAI之间的回归关系的方法计算烟草LAI, 该方法简单灵活, 应用广泛, 是一种最为常用的方法。通过光谱仪实测数据计算NDVI由以下公式得到：

(3)

式中, R₈₀₀和R₆₇₀分别表示在光谱值为800和670时对应烟草植被冠层光谱反射率值。

LAI与NDVI的关系模型主要有指数关系和对数关系2种类型^[17-18]：指数模型和对数模型, 为简化LAI的获取函数, 把获得的NDVI数据和实测LAI数据进行了回归分析, 并结合已有的烟草实验^[1]以及其它参考^[19-20], 经计算得到的LAI与NDVI对数模型如下：

在95%的置信区域时相关系数和调整后决定系数值分别为0.8953和0.8804, 标准误差RMSE为0.6375。利用上述模型, 通过平均光谱反射率计算得到各个时间南雄市烟草LAI理论计算值 (图 3)。LAI总体呈现先上升后下降趋势, 在前30d左右缓慢上升, 75d左右达到最大值4.72, 75d以后LAI值则急剧下降。

4.2 南雄烟草LOGISTIC模型构建

在对建模LAI样本数据初步筛选后, 对烟草实测LAI和南雄市归一化积温数据DS进行的简化后的LOGISTIC模型构建, 通过计算, 得到方程为：

由图 4中可以看出, 实测烟草LAI整体呈现先上升后下降的趋势, LAI值前期缓慢上升, 生长后期则快速下降, 在归一化积温0.6左右 (生长期大约80d左右) 时达到最大值。LOGISTIC模型对烟草实测LAI的整体拟合精度良好, 和烟草实测LAI的变化趋势基本保持一致。

4.3 数据同化

同化过程采用集合卡尔曼滤波同化算法进行, 运行的状态变量是LAI, 将实验区NDVI计算得到LAI和简化的LOGISTIC模型拟合LAI数据通过集合卡尔曼滤波算法进行同化。其中, 模型算子M是简化的作物生长模型LOGISTIC模型, 观测变量LAI是由实测光谱反射率数据计算NDVI进而获得的LAI值, 计算中初值X_{i, 0}^a用初始时刻LAI值平均值表示。南雄市烟草LAI数据同化结果如图 5。

假设实际测量LAI的平均值为真实LAI值, 从图 5可以得到, 单独使用LOGISTIC能够较好的整体拟合烟草LAI的变化趋势, 计算简单, 能够直接使用数学模型表示, 但是计算误差相对较大, 得到的LAI值会出现前期偏小, 后期偏大的情况, 分别在移栽后15—45d和90—120d内非常明显。NDVI计算LAI方法和LAI同化方法结果则相比较为精确, 基本与LAI实测均值曲线吻合。与NDVI计算LAI方法相比, 数据同化综合NDIV计算LAI结果和LOGISTIC模型拟合结果, 将其按权重分配, 得到的LAI值介于两者之间, 尽可能的消除二者的误差, 得到的LAI同化数据则更接近实测平均LAI, 精度相比较于单独NDVI计算LAI方法或模型拟合方法有所提高。同化LAI值呈先缓慢上升, 后急剧下降的趋势, 在大田期75d左右达到最大值, 最高LAI值达到4.7, 与烟草南雄市烟草LAI生长走势基本吻合。

与实际测量LAI数据作回归分析, 发现3种方法的相关系数都相对较高, 决定系数均在0.8以上, 表明3种方法均可一定程度上反应烟草LAI的生长情况。由图 6发现：同化得到的LAI值与观测数据基本一致, 同化LAI值与LAI实际观测样本的回归分析关系式为y=0.9904x+0.0126, 在95%的置信区域时相关系数为0.9456, 相比较高于NDVI计算LAI方法的0.8859。

根据y=x散点图对LAI的分布情况比较发现：数据同化方法、NDVI计算LAI方法和LOGISTIC模型拟合3种方法均可一定程度上表征烟草LAI的变化状态, 其中LOGISTIC模型拟合方法效果最差, NDIV计算LAI方法次之, 集合同化方法效果最好。3种方法在烟草作物成熟时期LAI拟合均取得很好的效果, LAI>4时尤为明显, 3种方法LAI值均均匀分布在y=x曲线两侧；而在LAI＜4时, NDVI计算LAI方法和LOGISTIC拟合LAI方法则会发生一定程度的离散, NDVI计算LAI方法在LAI＜2时尤为明显, LOGISTIC则会出现在每个时间节点拟合LAI相同的情况, 虽然此时实际的LAI值不同, 数据同化的总体效果则相对较好, 基本均匀分布于y=x的两侧。

5 结论与讨论

数据同化技术, 可将拟合模型和观测数据两者信息有效的结合起来, 观测数据依次被引入到过程模型中, 不断减少或者滤掉过程模型的噪声, 使得同化后结果更加贴近自然界的真实状态。文中通过实测光谱遥感数据计算了广东省南雄市2014年2—6月的烟草叶面积指数LAI, 模拟了南雄市烟草简化作物生长模型LOGISTIC模型, 并采用LOGISTIC模型作为集合卡尔曼滤波算法的动态模型, 对烟草遥感计算得到的LAI数据进行卡尔曼滤波数据同化, 初步实现了集合卡尔曼滤波同化算法在烟草LAI反演中的应用。结果显示同化得到的LAI结果基本和真实值吻合, 精度有所提高, LAI曲线更符合烟草的实际生长状况。

相较于传统选择多种微观指标进行回归分析方法, 该方法以LAI数据为对象, 直接对LAI进行计算处理, 无需测量和筛选与LAI高度相关的微观指标, 操作方法相对较为简单, 易于实现。数据同化方法、NDVI计算LAI方法和LOGISTIC模型拟合3种方法对比发现：3种方法均能一定程度的表征烟草LAI的生长变化, 但相比较与单一的LOGISTIC模型拟合方法或NDVI计算LAI方法, 集合卡尔曼滤波同化算法综合遥感观测数据和LOGISITC模型二者信息共同拟合, 得到的LAI同化结果更为接近真实值, 结果的总体误差更小。研究结果表明, 同化遥感信息与作物生长模型对提高烟草作物LAI信息监测是有效的, 可用于获取精度相对更高的烟草LAI信息。

文中发现：LOGISTIC模型不能有效的指示烟草的生长状况, 主要是因为LOGISTIC模型更多的是拟合得到积温和LAI的回归关系, LOGISTIC能够有效用数学模型表现烟草LAI变化趋势, 但得到的结果更多与积温数据相关, 而同一时间点积温值相同, 导致LAI计算结果在同一时间点全部相同, 不能表现烟草不同植株生长的差异性, 更不能有效表征出诸如天气、病虫害等外界因素导致烟草LAI迅速变化的状况, 导致拟合结果与实际情况稍有偏差。NDVI计算LAI方法在烟草长势最好时候计算效果较佳, 在初生长期和枯萎期则会发生LAI值过大过小的现象, 究其原因可能是生长前期和后期烟草的冠层面积过小导致野外测量的不稳定, 如烟草早期尚未成熟和后期叶片成熟枯萎, 植被冠层叶面积过小, 会加入周边地物的光谱信息或者光谱测量结果为叶面光谱信息。数据同化方法则表现最优, 数据同化得到的LAI基本均匀分布于y=x曲线两侧, 总体拟合效果良好。主要是集合卡尔曼滤波同化算法具有动态调节的优势, 它能够综合作物生长模型和观测数据两者的优势, 使得整体的拟合效果达到最优。该过程可有效的减小测量值的误差敏感性, 有利于减小同化后值的误差。

总体来说, 目前烟草研究相对于小麦、玉米等作物较为落后。烟草研究尚未出现专门的研究模型, 研究中对烟草模拟简单使用基于统计原理的LOGISTIC模型, 该模型更多的是一种统计性表达, 探究并应用一种如水稻、小麦等普适机理性模型可进一步增加数据同化算法面向烟草应用的专业性。另外, 实验应用光谱数据初步验证了集合卡尔曼滤波数据同化算法在烟草光谱应用的可能性和效果, 但相关实验应用和验证的范围相对狭窄, 进一步应用如卫星数据更有利于烟草LAI同化工作的大面积和自动化的进行。