技术路线如图1所示。 首先使用差分法结合光谱突变法提取了研究区内冬小麦每年的种植面积分布; 其次使用CASA模型在两种模拟情景下模拟每年的植被NPP; 最后, 结合收获指数和提取的冬小麦面积, 获取2000年— 2016冬小麦单产的时空分布并进行分析。

1.2.1 冬小麦面积提取方法

在研究时段内, 山东省逐年冬小麦面积分布的提取方法使用了先前研究工作提出的提取流程, 以进一步验证该提取方法的普适性。 提取过程中, 使用农业气象站观测播种期与成熟期和气象再分析资料获取的站点气温、 降水和辐射建立多元线性回归模型, 以使用气象再分析资料模拟冬小麦播种期和成熟期的空间分布; 采用Savitzky-Golay(S-G)滤波重构时序EVI序列, 根据模拟的播种期和成熟期的空间分布, 截取每个像元在播种期至成熟期之间的重构EVI序列, 采用差分法计算截取后得到的EVI序列的峰值频数; 获取每个像元成熟期前后的EVI, 采用光谱突变法计算其在冬小麦成熟后相较于成熟前的EVI突变情况(Slope); 提取波峰个数为2和光谱下降斜率低于负0.02的像元, 并与旱地分布取交集, 获取冬小麦的空间分布。 其中, 光谱突变法计算Slope的公式如式(1)

Slope=EVIm+16-EVIm-16EVIm-16(1)

式(1)中, EVIm-16和EVIm+16分别为冬小麦成熟期前、 后16 d的EVI值。

1.2.2 CASA模型简介

首先使用光能利用率模型CASA模拟山东省植被NPP, 模型的总体设计如图2所示。 CASA模型最早由Potter[2]提出, 模型表述如式(2)— 式(4)

NPP(x, t)=APAR(x, t)×ξ(x, t)(2)

APAR(x, t)=RSG(x, t)×FPAR(x, t)×0.5(3)

ξ(x, t)=Tξ1(x, t)×Tξ2(x, t)×Wξ(x, t)×ξmax(4)

其中, APAR(x, t)为像元x在t月份吸收的光合有效辐射, 单位为MJ· m-2· mon-1; ξ (x, t)为像元x在t月份的实际光能利用率, 单位为g· C· mJ-1。 RSG(x, t)为太阳总辐射, 单位为MJ· m-2· mon-1; FPAR(x, t)表示植被对入射光合有效辐射的吸收比例, 无量纲; 0.5表示植被所能利用的光合有效辐射占太阳总辐射的比例。 Tξ 1(x, t)和Tξ 2(x, t)表示温度胁迫因子; Wξ (x, t)为水分胁迫因子; ξ max为植被最大光能利用率。 FPAR(x, t), Tξ 1(x, t), Tξ 2(x, t)和Wξ (x, t)的计算可见参考文献[3]。

为探究最大光能利用率对作物单产模拟的影响, 在驱动CASA模型时, 设置了Case1和Case2两种模拟情景, 其中Case1设置ξ max为不随时间变化的固定值2.8 g· C· mJ-1, Case2设置ξ max为随时间变化的序列值。 参照前人研究工作, 在Case2中设定2014年冬小麦ξ max为2.8 g· C· mJ-1, 按ξ max每年增加0.03 g· C· mJ-1· a-1的速率[9]推算2000年— 2016年间其他年份的冬小麦最大光能利用率。

1.2.3 冬小麦单产估算、 验证与分析

光能利用率模型模拟作物单产的原理如式(5)和式(6)

Yield(x)=NPPsum(x)×Tc×HI(5)

NPPsum(x)=∑t=35NPP(x, t)(6)

其中, Yield(x)为像元x的模拟单产(kg· hm-2); NPPsum为冬小麦生长季内NPP的和, 根据山东省冬小麦物候特征, 取3月— 5月NPP的和[5]; Tc为植物体内碳素转换为干物质的转换系数, 取值为2.22[5]; HI为收获指数, 取其值为0.48, 与已有研究[5]中的收获指数取值接近; NPP(x, t)为像元x在t月份的NPP模拟值, 单位为g· C· m-2· mon-1。

使用提取的冬小麦面积分布掩膜模拟结果, 获取Case1和Case2的冬小麦模拟单产, 采用年鉴统计单产数据分别对两种模拟结果进行验证, 对比验证结果, 选择精度较高的模拟结果作为研究区冬小麦单产, 分析冬小麦单产的时空特征。