一种农作物叶面积指数同化方法 【技术领域】
本发明涉及叶面积指数技术领域， 特别涉及一种农作物叶面积指数同化方法。背景技术 叶面积指数是一个重要的植物学参数， 在农作物长势监测、 产量估算等领域应用 广泛。 在农业定量遥感研究中， 叶面积指数的提取和反演估算主要是基于遥感观测数据、 作 物生长模型、 辐射传输模型等开展的。
目前， 叶面积指数估算方法主要有经验公式法、 非参数法、 物理模型反演法和数据 同化算法等。 其中， 经验公式法从数据自身统计特性出发， 通过分析光谱数据或各类植被指 数与叶面积指数的定量关系实现叶面积指数的反演估算 ； 非参数法以其所描述的数据关系 为对象进行数值拟合， 数据量大小在一定程度上影响叶面积指数估算结果的可靠性 ； 物理 模型反演法以物理机制为基础， 选用合适的数学工具反演估算叶面积指数。
数据同化是近年来发展较为迅速的一类用于地表参数反演估算的研究思路， 其在 考虑数据时空分布、 观测场误差和背景场误差的基础上， 将观测数据融入模型的动态运行 过程， 以实现目标参数的提取。现有的经典同化算法有滤波法、 变分法和启发式优化算法 等。 其中， 集合卡尔曼滤波同化算法主要用于单状态变量优化， 无法直接应用到基于多状态 变量优化求解的叶面积指数同化反演过程 ； 四维变分同化算法适于解决多状态变量优化问 题， 但线性模式和伴随模式的引入增加了计算复杂度 ； 非常快速模拟退火算法独立于目标 泛函， 优于处理非线性及不连续性系统。 三种算法在各自的适用范围和使用条件的基础上， 各具优势 ； 同时， 每种算法又无法兼具另外两种算法的优势。因此， 如果能够提供一种综合 三种算法优势的同化方法， 将对农作物长势监测、 产量估算等领域具有重要意义。
发明内容 ( 一 ) 要解决的技术问题
本发明要解决的技术问题是 ： 如何提供一种农作物叶面积指数同化方法， 以提高 农作物叶面积指数估算的精度。
( 二 ) 技术方案
为解决上述技术问题， 本发明提供一种农作物叶面积指数同化方法， 其包括步 骤：
B： 依据四维变分的代价函数构建目标泛函 ；
C： 在当前同化过程下， 依据非常快速模拟退火算法求取所述目标泛函的最优解 ；
D： 判断当前同化过程的次数是否大于总观测次数， 如果是， 将所述最优解作为待 调整参数的最终值， 同化过程结束 ； 否则， 将所述最优解作为下一次同化过程中所述待调整 参数的当前值， 执行所述步骤 C。
优选地， 在所述步骤 B 之前还包括步骤 A ： 基于农作物时序观测数据和地面先验知 识设定待调整参数的初值， 根据所述待调整参数的初值和所述作物生长模型生成符合高斯
分布的背景场数据集合初值， 以及所述背景场数据集合的误差协方差矩阵初值。
优选地， 所述步骤 A 具体包括步骤 ：
A1 ： 选取对叶面积指数变动敏感的参数作为作物生长模型的待调整参数 X ；
A2 ： 基于农作物时序观测数据和地面先验知识设定所述待调整参数的初值 X(t0) ； 其中， t0 表示同化开始时刻 ；
A3 ： 根据所述待调整参数的初值 X(t0) 和所述作物生长模型， 生成符合高斯分布的 a 背景场数据集合初值 X (t0)， 以及所述背景场数据集合的误差协方差矩阵初值 Pa(t0)。
优选地， 所述步骤 B 中的目标泛函 J(X) 的表达式如下 ：
其中， R(ti) 表示第 i 个观测时刻由观测数据测量误差和模型误差组成的观测误 差； Y (ti) 表示第 i 个观测时刻的光谱数据集 ； Xa(tk) 表示第 k 次同化过程时的背景场数据 集合的当前值， k 的初值为 0 ； Pa(tk) 表示 Xa(tk) 的误差协方差矩阵的当前值 ； Num 表示所述 总观测次数 ； M 表示作物生长模型算子 ； H 表示 PROSAIL 辐射传输模型算子。
优选地， 所述步骤 C 具体包括步骤 ：
C1 ： 在当前同化过程下， 设定标准温度 T0， 设定优化过程的状态变量 X0(tk，j) 的初
0值并且设定初始温度 T1 ＝ T0exp(-c) ； 其中，表示 Xa(tk) 的算术平均值， c 为常数， j 的初值为 0 ； 0
C2 ： 在当前温度 TI 下， 根据当前状态变量 X0(tk， j) 计算下一个状态变量 X (tk， j+1) ； I 表示迭代次数， 初值为 1 ； 0
C3 ： 计算 ΔJ ＝ J(X0(tk， 判断 ΔJ 是否大于 0， 如果是， 接受 X0(tk， j))-J(X (tk， j+1))， 否则， 以概率 exp(-ΔJ/cTI) 接受 X0(tk， j+1) ； j+1)，
C4 ： 判断是否达到平衡， 如果是， 执行 C5 ； 否则， 计算下一温度 TI+1 ＝ T0exp(-c×I1/ n )， I 自加 1， j 自加 1， 然后执行步骤 C2 ； 其中， n 为待调整参数的维数 ；
C5 ： 判断当前温度 TI 是否低于预设最低温度， 或者已经连续出现了预设阈值个状 0 态变量均未被接受， 如果是， 将 X (tk， 执行步骤 D ； 否则， TI+1 j+1) 作为所述目标泛函的最优解， ＝ TI， I 自加 1， j 自加 1， 然后执行步骤 C2。
优选地， 所述步骤 C2 中通过下面两个公式计算下一个状态变量 X0(tk， j+1) ： 0 |2u-1|
δX0(tk， -1] ； j) ＝ sgn(u-0.5)×X (tk， j)×[(1+1/TI) 0 0 0
X (tk， j+1) ＝ X (tk， j)+δX (tk， j)(Rh-Rl) ；
其中， sgn 表示符号函数 ； u 表示随机生成的均匀分布变量， 并且 u ∈ [0， 1] ； Rh 表 示所述待调整参数 X 的物理上限 ； Rl 表示所述待调整参数 X 的物理下限。
优选地， 所述步骤 C4 具体包括步骤 ：
C41 ： 通过下面公式计算标准值 Stand ：
C42 ： 判断所述标准值 Stand 是否小于预设极小值 ε， 如果是， 认为达到平衡， 执行 1/n C5 ； 否则， 认为未达到平衡， 计算下一温度 TI+1 ＝ T0exp(-c×I )， I 自加 1， j 自加 1， 然后执
行步骤 C2。
优选地， 所述预设极小值 ε 为 0.01。
优选地， 所述标准温度 T0 为 100℃， 所述预设最低温度为 10℃。
优选地， 所述步骤 D 具体包括步骤 ：
D1 ： 判断当前同化次数 k 是否大于所述总观测次数 Num， 如果是， 将所述最优解作 为所述待调整参数的最终值， 同化过程结束 ； 否则， 执行步骤 D2 ；
D2 ： k 自加 1， 将所述最优解作为第 k 次同化过程中所述待调整参数的当前值 X(tk)， 根据所述待调整参数的当前值 X(tk) 和所述作物生长模型生成符合高斯分布的背景 场数据集合当前值 Xa(tk)， 以及所述背景场数据集合的误差协方差矩阵当前值 Pa(tk)， 执行 步骤 C。
( 三 ) 有益效果
本发明所述农作物叶面积指数同化方法， 融合了现有多种估算方法的优势， 通过 将时序观测数据和作物生长模型模拟数据融入到背景场数据集合的迭代更新模式和代价 函数的优化求解模式中， 在获取全局最优解的同时， 降低了计算复杂度， 提高了叶面积指数 的估算精度和抗数据饱和性。 附图说明
图 1 是本发明实施例所述的农作物叶面积指数同化方法流程图 ； 图 2 是本发明实施例所述方法与现有技术方法的叶面积指数估算效果对比示意图。 具体实施方式
下面结合附图和实施例， 对本发明的具体实施方式作进一步详细描述。以下实施 例用于说明本发明， 但不用来限制本发明的范围。
图 1 是本发明实施例所述的农作物叶面积指数同化方法流程图。如图 1 所示， 所 述方法包括 ：
步骤 A ： 基于农作物时序观测数据和地面先验知识设定待调整参数的初值， 根据 所述待调整参数的初值和所述作物生长模型生成符合高斯分布的背景场数据集合初值， 以 及所述背景场数据集合的误差协方差矩阵初值。
所述步骤 A 具体包括 ：
步骤 A1 ： 选取对叶面积指数变动敏感的参数作为作物生长模型的待调整参数 X。 在本发明实施例中， 以冬小麦作为示例农作物， 因此， 选取灌溉日期、 灌溉量、 施肥日期和施 肥量作为 CERES-WHEAT 作物生长模型的待调整参数 X。 但是， 本发明并不具体限定所述待调 整参数 X 和作物生长模型的类别， 比如当以水稻作为叶面积指数同化方法的客体时， 还可 以选择 CERES-RICE 作物生长模型。
步骤 A2 ： 基于农作物时序观测数据和地面先验知识设定所述待调整参数的初值 X(t0) ； 其中， t0 表示同化开始时刻 ；
步骤 A3 ： 根据所述待调整参数的初值 X(t0) 和所述作物生长模型， 生成符合高 a 斯分布的背景场数据集合初值 X (t0)， 以及所述背景场数据集合的误差协方差矩阵初值Pa(t0)。
步骤 B ： 依据四维变分的代价函数构建目标泛函。 所述目标泛函 J(X) 的表达式如下 ：其中， R(ti) 表示第 i 个观测时刻由观测数据测量误差和模型误差组成的观测误 差， i 的取值从 1 至 Num ； Num 表示所述总观测次数 ； Y0(ti) 表示第 i 个观测时刻的光谱数据 集； Xa(tk) 表示第 k 次同化过程时的背景场数据集合的当前值， k 的初值为 0 ； Pa(tk) 表示 Xa(tk) 的误差协方差矩阵的当前值 ； M 表示作物生长模型算子 ； H 表示 PROSAIL 辐射传输模 型算子。
步骤 C ： 在当前同化过程下， 依据非常快速模拟退火算法求取所述目标泛函的最 优解。
所述步骤 C 具体包括 ：
步骤 C1 ： 在当前同化过程下， 设定标准温度 T0， 设定优化过程的状态变量 X0(tk，j)
的初值并且设定初始温度 T1 ＝ T0exp(-c) ； 其中，表示 Xa(tk) 的算术平均值， c 为常数， j 的初值为 0。所述标准温度 T0 一般为 100℃。 0
步骤 C2 ： 在当前温度 TI 下， 根据当前状态变量 X0(tk， j) 计算下一个状态变量 X (tk， 初值为 1。 j+1)。I 表示迭代次数， 0
步骤 C3 ： 计算 ΔJ ＝ J(X0(tk， 判断 ΔJ 是否大于 0， 如果是， 接受 j))-J(X (tk， j+1))， 0 0 X (tk， 否则， 以概率 exp(-ΔJ/cTI) 接受 X (tk， j+1)， j+1)。
步 骤 C4 ： 判 断 是 否 达 到 平 衡， 如 果 是， 执 行 C5 ； 否 则， 计 算 下 一 温 度 TI+1 ＝ 1/n T0exp(-c×I )， I 自加 1， j 自加 1， 然后执行步骤 C2。其中， n 为待调整参数的维数。
所述步骤 C4 具体包括 ：
步骤 C41 ： 通过下面公式计算标准值 Stand ：
步骤 C42 ： 判断所述标准值 Stand 是否小于预设极小值 ε， 如果是， 认为达到平衡， 1/n 执行 C5 ； 否则， 认为未达到平衡， 计算下一温度 TI+1 ＝ T0exp(-c×I )， I 自加 1， j 自加 1， 然后执行步骤 C2。所述预设极小值 ε 一般为 0.01。
步骤 C5 ： 判断当前温度 TI 是否低于预设最低温度， 或者已经连续出现了预设阈值 0 个状态变量均未被接受， 如果是， 将 X (tk，j+1) 作为所述目标泛函的最优解， 执行步骤 D ； 否 则， TI+1 ＝ TI， I 自加 1， j 自加 1， 然后执行步骤 C2。所述预设最低温度一般为 10℃， 所述预 设阈值一般为 5。
步骤 D ： 判断当前同化过程的次数是否大于总观测次数， 如果是， 将所述最优解作 为待调整参数的最终值， 同化过程结束 ； 否则， 将所述最优解作为下一次同化过程中所述待 调整参数的当前值， 执行所述步骤 C。
所述步骤 D 具体包括 ：
步骤 D1 ： 判断当前同化次数 k 是否大于所述总观测次数 Num， 如果是， 将所述最优
解作为所述待调整参数的最终值， 同化过程结束 ； 否则， 执行步骤 D2 ；
步骤 D2 ： k 自加 1， 将所述最优解作为第 k 次同化过程中所述待调整参数的当前值 X(tk)， 根据所述待调整参数的当前值 X(tk) 和所述作物生长模型生成符合高斯分布的背景 场数据集合当前值 Xa(tk)， 以及所述背景场数据集合的误差协方差矩阵当前值 Pa(tk)， 执行 步骤 C。
图 2 是本发明实施例所述方法与现有技术方法的叶面积指数 (LAI) 估算效果对比 示意图。图 2 中左侧子图是基于集合卡尔曼滤波同化算法得到的模拟 LAI 与实测 LAI 的散 点图， 中间子图是基于四维变分同化算法得到的模拟 LAI 与实测 LAI 的散点图， 右侧子图是 基于本发明方法得到的模拟 LAI 与实测 LAI 的散点图。通过图 2 可以看到， 本发明方法相 比现有技术方法具有以下优势 ：
(1) 在整体估算精度方面， 本发明方法的拟合结果最佳。 集合卡尔曼滤波顺序同化 算法、 四维变分同化算法和本发明方法的拟合结果为 ： 实测叶面积指数和模拟叶面积指数 的决定系数分别为 0.42、 0.60 和 0.80， 均方根误差分别为 0.84、 0.69 和 0.49， 准确度分别 为 0.92、 1.00 和 1.00。
(2) 在抗数据饱和性方面， 通过对叶面积指数大于等于 3.00( 图中圆圈对应的散 点 ) 的饱和数据进行拟合效果分析可知， 本发明方法的抗数据饱和性效果最优。集合卡尔 曼滤波顺序同化算法、 四维变分同化算法和本发明的同化方法的拟合结果为 ： 实测叶面积 指数和模拟叶面积指数的决定系数分别为 0.01、 0.19 和 0.40， 均方根误差分别为 0.67、 0.63 和 0.51， 准确度分别为 0.17、 0.44 和 0.64。 本发明实施例所述农作物叶面积指数同化方法， 融合了现有多种估算方法的优 势， 通过将时序观测数据和作物生长模型模拟数据融入到背景场数据集合的迭代更新模式 和代价函数的优化求解模式中， 在获取全局最优解的同时， 降低了计算复杂度， 提高了叶面 积指数的估算精度和抗数据饱和性。
以上实施方式仅用于说明本发明， 而并非对本发明的限制， 有关技术领域的普通 技术人员， 在不脱离本发明的精神和范围的情况下， 还可以做出各种变化和变型， 因此所有 等同的技术方案也属于本发明的范畴， 本发明的专利保护范围应由权利要求限定。