基于CASA-WOFOST耦合模型的大豆单产遥感估算研究

图1 双山农场位置及地块分布

Fig.1 Location of Shuangshan farm and the plot distribution

2.2　数据源及预处理

2.2.1　时间序列遥感数据

考虑研究区影像可获得性和数据空间分辨率，本文选择环境与灾害监测预报小卫星星座 A、B 星（HJ-1 A/B）数据作为主要遥感数据源，由中国资源卫星应用中心网站获取(http://www.cresda.com)。HJ-1A/B星搭载宽幅CCD相机，幅宽700 km，空间分辨率为30 m，红、绿、蓝及近红外等4个谱段推扫成像，本文选择2016年环境星数据8期，可覆盖研究区大豆各个生长发育阶段，以满足作物单产估算对数据源时间分辨率的要求。环境星数据成像时间及质量如表1所示。对遥感影像进行辐射定标、大气校正、几何精校正等预处理。

表1 研究区时间序列HJ-CCD数据列表

Table 1 Time-series HJ-CCD data list of study area

影像编号	获取时间	成像传感器	轨道编号	影像质量
N1	2016-05-22	HJ-1B CCD1	451-53	优
N2	2016-06-01	HJ-1A CCD2	450-56	优
N3	2016-06-27	HJ-1B CCD2	449-56	良（薄云）
N4	2016-07-06	HJ-1B CCD1	449-55	良（薄云）
N5	2016-07-19	HJ-1B CCD1	453-53	优
N6	2016-07-31	HJ-1B CCD1	450-53	优
N7	2016-08-20	HJ-1B CCD2	452-56	优
N8	2016-09-30	HJ-1B CCD2	450-56	良（薄云）

2.2.2　气象数据

气象数据是农作物单产遥感估算模型的重要驱动数据，本研究所使用的气象数据由美国国家海洋和大气管理局(ftp://ftp.ncdc.noaa.gov/pub/data/ gsod/2017/)免费提供。采集的气象数据主要包括平均温度、最高温度、最低温度、降雨量、相对湿度、平均风速、平均海平面气压。日照时数由欧盟官方网站免费提供，网址为(http://spirits.jrc.ec.europa.eu/files/ecmwf/ope/asia/rad/)，所获取的气象数据时间分辨率为1 d。分析气象数据得知，研究区在大豆生长发育阶段平均温度为18.455 ℃，平均降雨量为32.535 mm，平均日照时数为8.569 h，与研究区实际情况一致。

2.2.3　田间调查数据

本研究是针对大豆地块开展作物单产遥感估算研究，田间调查主要实测一些模型参数以及单产数据，以便于进行耦合模型模拟结果的精度验证。本文实地调查了20个具有代表性的大豆地块（图1），数据采集时间为2016年10月初前后，在农场收割完毕后按照田块统计其产量，并根据田块的大小可计算出实测地块的单产。

2.2.4　其他基础数据

本研究所用的基础数据主要包括研究区地块边界矢量数据以及作物分布图。其中地块边界矢量数据通过空间分辨率6.5 m的RapidEye数据绘制，作物分布图根据农场耕作计划以及遥感影像分类结果获取，并据此提取研究区大豆地块。

3 研究方法

3.1　总体技术路线

首先利用遥感数据，利用线性插值和S-G（Savitzky-Golay）迭代滤波求取作物时序NDVI，再根据作物分布数据、气象数据等，实现以日时间为步长运行CASA模型，逐日模拟作物生物量，再通过将WOFOST模型作物生长模块的干物质分配过程与CASA模型进行耦合，模拟出研究区大豆单产。并与单独利用CASA模型模拟大豆单产进行对比分析，其主要流程如图2所示。

图2

图2 研究方法流程图

Fig.2 Methods flowchart

3.2　时间序列NDVI

NDVI时间序列能够有效地反应作物从播种到收获NDVI变化的全过程，是求取光合有效辐射吸收比FPAR，支撑运行CASA模型的关键，而由表1可知，本研究所利用的HJ-CCD影像平均时间间隔为17 d，这样达不到本研究实现对研究区农作物连续监测的时间分辨率要求。为了填补两个时相影像之间空缺的数值，需要获取时间序列DNVI，这样能更精确地刻画作物生长细节，且WOFOST模型是以日时间步长计算作物，计算出每日NDVI为后续模型的耦合提供相同的时间尺度。

由于作物的生长是一个缓慢的过程，作物DNVI值不可能出现急剧上升或者下降的情况^[17]，因此本文首先根据相邻点进行线性插值来替代作物生长过剩中受到环境影响产生的噪声^[18]，然后利用迭代滤波方法（Savitzky-Golay（S-G）滤波）^[19]进一步对NDVI数据进行拟合，消除环境噪声的影响，得到如图3所示的NDVI拟合结果。

图3

图3 大豆NDVI时间序列曲线图

Fig.3 NDVI time-series curve of soybean

3.3　CASA模型及参数标定

3.3.1　CASA模型

CASA模型^[20]作为光能利用率模型应用最广泛的模型之一，其产量预估模式为生物量与收获指数相结合，公式为：

Y i e l d = \frac{B i o m a s s \times H I}{1 - θ_{g r a i n}}

(1)

其中：Biomass表示作物生长期内的生物量；HI表示收获指数；研究区收获指数值为0.47^[21]；θ_grain指作物籽粒含水量，研究区大豆籽粒含水量值为12.5%^[22]。生物量和单产之间的转换就是通过收获指数来实现的。作物生长期内的生物量公式为：

B i o m a s s = \sum_{t = t_{0}}^{t_{n}} N P P (x, t) = \sum_{t = t_{0}}^{t_{n}} (A P A R (x, t) \times ε (x, t))

(2)

其中：Biomass表示作物生长期内的生物量，NPP为作物生长期内的净初级生产力，t为作物的不同生长时段，以天为单位，x为作物生长的空间位置，t₀表示作物开始生长阶段，t_n是作物收获阶段，APAR(x,t)为光合有效辐射，单位为MJ/m²， $ε (x, t)$ 为实际光能利用率，单位为gC/MJ。

A P A R (x, t) = S O L (x, t) \times F P A R (x, t) \times 0.5

(3)

其中：SOL(x,t)为地表太阳总辐射量，单位为MJ/m²，0.5表示植被所能利用太阳辐射的比率，FPAR(x,t)为光合有效辐射吸收比率。

ε (x, t) = T_{ε 1} (x, t) \times T_{ε 2} (x, t) \times W_{ε} (x, t) \times ε_{m a x}

(4)

其中： $ε (x, t)$ 为实际光能利用率，单位为gC/MJ，T_ε1(x,t)、 T_ε1(x,t)为温度胁迫因子，表示低温、高温对作物光能利用率胁迫作用的大小；W_ε (x,t)表示水分胁迫因子；ε_max表示没有任何条件的限制下作物的最大光能利用率。

3.3.2　CASA模型参数标定

（1）光合有效辐射吸收比FPAR。Los和Field等^[23,24]通过将模拟的FPAR与实测FPAR进行对比分析发现由NDVI模拟出的FPAR比实际偏高，SR模拟出的FPAR比实际偏低。Los^[23]引入修正系数α，其取值范围为0.1~0.9，测试结果当α=0.5时，实测值与FPAR相关性最好，因此，本文求取作物光合有效辐射吸收比公式为：

\begin{array}{l} F P A R_{N D V I} (x, t) = \frac{(N D V I (x, t) - N D V I m i n) (F P A R m a x - F P A R m i n)}{(N D V I m a x - N D V I m i n)} + F P A R m i n \end{array}

(5)

F P A R_{S R} (x, t) = \frac{(S R (x, t) - S R m i n) (F P A R m a x - F P A R m i n)}{(S R m a x - S R m i n)} + F P A R m i n

(6)

S R (x, t) = \frac{1 + N D V I (x, t)}{1 - N D V I (x, t)}

(7)

F P A R (x, t) = \frac{F P A R_{N D V I} (x, t) + F P A R_{S R} (x, t)}{2}

(8)

其中：FPAR_NDVI(x,t)为利用归一化植被指数估算的FPAR结果，FPAR_SR(x,t)表示利用比值植被指数估算的FPAR结果。此外，计算FPAR时取NDVI的最大取值和最小取值为NDVI 95%和5%下侧分位数，分别取值为0.175、0.736，基于Perman等^[25]提出的利用纬度、日照时数、日序等参数估算太阳总辐射。

（2）实际光能利用率(ε)计算。光能利用率(ε)是指植被将吸收的光合有效辐射(APAR)通过冠层光合作用转化为有机碳的效率，反映了植被将太阳能转换为有机物的能力。Potter等^[20]认为植被在理想条件下具有最大的光能利用率，但随着生长季内的环境条件的改变，实际植被对光的利用效率也不断发生变化的，且主要受到温度和水分因子影响。本研究将大豆生长的最适温度定为20 ℃。

根据双山农场所处的气候环境以及农场管理条件，认为研究区大豆出现水分胁迫影响的可能极低，因此本研究认为在大豆生长发育过程中不受水分胁迫，将水分胁迫因子W_ε (x,t)设为1。大豆属于C3作物^[26]，其最大光能利用率上限可达到3.0 gC/MJ，结合研究区气候、植被类型等特点，将最大光能利用率设为2.6 gC/MJ^[27]。

3.4　CASA模型与WOFOST模型耦合

WOFOST(World Food Studies)模型是由荷兰瓦赫宁根大学(Wageningen)和世界粮食组织(CWFS)共同研发的一个作物生长模拟模型^[28]，它是一种以1 d为步长、综合作物、土壤、气候以及管理等因素模拟作物生长动态过程的面向过程、机理性的动态解释性模型^[29,30]。WOFOST模型主要包括作物生长模拟模块、土壤水分模拟模块、作物蒸散模块三部分，这3个模块协同作用，共同组成WOFOST模型对作物呼吸作用、光合作用、生长发育、干物质积累及分配等生长模拟的全过程。模型的运行原理是通过作物生长模拟模块和作物蒸散模块以CO₂为核心驱动因子求出潜在产量，然后根据水分、养分等限制算出胁迫产量，再减去作物生长呼吸作用以及蒸腾作用得到净同化量，最终根据存储器官分配系数计算出作物最终产量。

CASA模型建立初期大量应用于估算全球尺度下的陆地生态系统净初级生产力，应用CASA模型对小区域尺度甚至田块尺度的农作物监测研究相对较少。CASA模型的较为简单实用，模型运行速度较快，需要获取或标定的参数较少且较容易获得，但是CASA模型模拟作物的净初级生产力(NPP)后需要构建一个作物收获指数来实现与净初级生产力之间的转化，然而作物由光合作用的同化量到生物量再到作物最终的产量涉及很多复杂的生理生化过程，经过多步的转化，并不能由一个简单的公式求得的收获指数来实现这一过程。由于WOFOST模型能够很详细地模拟作物生长发育的各个过程，这导致WOFOST模型中参数非常多且有些参数获取代价大，参数标定困难，且模型运行速度慢。本文通过引入WOFOST模型作物生长模拟模块中的干物质分配过程代替收获指数模拟作物单产，求出由CASA模型模拟作物生物量扣除作物呼吸消耗再分配至存储器官的干物质重，从而模拟出作物产量。

WOFOST模型作物生长模型已在我国东北及华北地区取得比较理想的效果^[31,32,33]，但是对于不同的研究区域，由于生态环境、气候条件的差异以及不同土壤和管理条件差异，需要对WOFOST模型进行参数本地化，否则将会降低模型的通用性。本文主要通过文献资料调研^[33]、遥感数据和田间实测相结合的方式标定耦合模型中WOFOST模型模块的参数，标定结果如表2所示。

表2 WOFOST模型参数标定值

Table 2 The value of input crop parameters for the WOFOST model

参数名称	参数定义	数值	单位
TSUM1	出苗到开花所需积温	1 062.39	℃
TSUM2	开花到成熟所需积温	1 230.94	℃
Q10	温度升10 ℃呼吸消耗相对变化量	2	—
FO1(DVS=0.1)	DVS=0.1存储器官分配系数	0	kg/kg
FO2(DVS=0.5)	DVS=0.5存储器官分配系数	0	kg/kg
FO3(DVS=1.1)	DVS=1.1存储器官分配系数	0.23	kg/kg
FO4(DVS=1.7)	DVS=1.7存储器官分配系数	0.76	kg/kg
CVL	叶片干物质转换效率	0.76	kg/kg
CVS	茎干物质转换效率	0.72	kg/kg
CVO	存储器官干物质转换效率	0.45	kg/kg
CVR	根干物质转换效率	0.76	kg/kg

3.5　耦合模型模拟精度评价指标

本研究耦合模型模拟精度主要采用均方根误差(Root Mean Square Error, RMSE )、Nash-Sutcliffe效率系数(NSE)和线性回归判定系数(R²)。公式如下：

R M S E = \sqrt[]{\frac{1}{n} \sum_{i = 1}^{n} (X - X_{o b s})^{2}}

(9)

N S E = 1 - \frac{\sum_{i = 1}^{n} (X - X_{o b s})^{2}}{\sum_{i = 1}^{n} (X_{o b s} {- \bar{X})}^{2}}

(10)

其中：n表示观测值个数；X表示耦合模型模拟值；X_obs表示观测值。RMSE越低，模拟精度越高，NSE越接近于1，模拟精度越高。

4 结果与分析

4.1　单产结果分析

如图4所示为本研究利用CASA-WOFOST耦合估产模型模拟双山农场单产结果图，图5为研究区大豆单产像元频率分布直方图。由图5可知，双山农场大豆单产差异性不大，在1 500~1 900 kg/hm²之间分布，分析原因可能与研究区域过小，各区域自然状况、降雨、温度等变化不明显，造成研究区范围内单产结果不能够呈现一定的规律性变化。在像元尺度上分析大豆单产，CASA模型模拟单产的值比耦合模型的值更集中在一定范围内且CASA模型模拟结果总体比CASA-WOFOST耦合模型模拟的值更小，大豆单产大多分布在1 650~1 850 kg/hm²，CASA模型和CASA-WOFOST耦合模型单产小于1 650 kg/hm²的像元占总像元分别为7.48%、3.29%；1 650~18 50 kg/hm²像元占总像元分别为89.04%、81.22%；大于1 850 kg/hm²的分别占3.48%、15.48%。但从地块尺度分析两者模拟结果可知（图5），农场内部及田块内部均出现了一定的空间异质性，其中机务二队小五队由于地势较低，使得田块单产较机务一队、机务二队、机务三队整体偏低，分析田块内部可知地块内部的单产较高，地块边缘的单产较低，分析由于地块边缘受到人为因素干扰过大，地块内部单产由于耕作管理合理，水分以及肥料适宜且相对受到人为因素较少，使得单产较大。

图4

图4 大豆单产模拟结果分布图

Fig.4 Yield simulation results distribution of soybean

图5

图5 大豆单产统计像元值频率分布直方图

Fig.5 Yield statistical histogram of soybean

4.2　单产模拟精度分析

研究区共有40余个大豆地块，选取了其中20个具有代表性的地块，在大豆收获时统计每个地块的平均单产，作为本文的单产模拟的验证数据，验证地块名称及实测单产以及模拟单产如表3所示。其中一队4#下号和一队5#大区地块在CASA模型和耦合模型估产下都存在较大的误差，分析原因可能是该地块面积较小且作物种植稀疏，使得模型估产结果受土壤背景的影响较大，此外，可能由于人为因素在实地测量采样过程中采集数据存在一定误差。

表3 地块尺度下双山农场基地2016年大豆单产模拟精度分析

Table 3 Analysis results of Soybean Yield Simulation of Shuangshan farm at field scale in 2016

地块名称	实测单产/(kg/hm²)	CASA模型/(kg/hm²)	CASA-WOFOST模型/(kg/hm²)	CASA模型		CASA-WOFOST模型
地块名称	实测单产/(kg/hm²)	CASA模型/(kg/hm²)	CASA-WOFOST模型/(kg/hm²)	绝对误差(kg/hm²)	相对误差	绝对误差(kg/hm²)	相对误差
一队大头	1 713.62	1 644.47	1 704.25	69.14	4.20%	9.37	0.55%
一分部1#	1 704.57	1 730.42	1 696.50	25.84	1.49%	8.08	0.48%
一队副4#	1 823.64	1 807.28	1 763.26	16.36	0.91%	60.38	3.42%
一队4#下	1 662.37	1 779.45	1 705.93	117.07	6.58%	43.55	2.55%
一队5#3区	1 863.63	1 793.76	1 872.85	69.87	3.89%	9.22	0.49%
一队三条田	1 731.70	1 690.88	1 723.16	40.82	2.41%	8.54	0.50%
一队2#上	1 693.52	1 643.58	1 665.95	49.94	3.04%	27.57	1.65%
一队7#	1 773.40	1 788.88	1 818.46	15.48	0.87%	45.06	2.48%
一队8#	1 694.53	1 701.44	1 729.98	6.92	0.41%	35.45	2.05%
一分部10#	1 805.05	1 725.68	1 779.31	79.36	4.60%	25.74	1.45%
一分部9#	1 544.82	1 504.08	1 594.70	40.74	2.71%	49.88	3.13%
二分部2#	1 762.85	1 721.71	1 817.00	41.14	2.39%	54.15	2.98%
一队6#南	1 962.80	1 968.48	1 992.67	5.69	0.29%	29.88	1.50%
一分部2#	1 660.87	1 592.08	1 667.86	68.78	4.32%	7.00	0.42%
一队10#	1 783.95	1 653.80	1 727.90	130.14	7.87%	56.05	3.24%
一队11#	1 734.72	1 670.50	1 715.77	64.22	3.84%	18.95	1.10%
一队5#1区	1 835.19	1 793.54	1 840.15	41.66	2.32%	4.96	0.27%
一队12#	1 766.87	1 832.17	1 763.12	65.30	3.56%	3.75	0.21%
一队5#2区	1 807.06	1 815.55	1 814.28	8.49	0.47%	7.22	0.40%
一队5#大区	1 593.05	1 664.25	1 678.66	71.20	4.28%	85.62	5.10%

图6为精度对比分析，CASA模型线性回归判定系数(R²)为0.668 53，耦合模型R²为0.844 72，耦合模型的预测结果比较理想。分析不同的精度评价指标如表4所示，CASA模型和CASA-WOFOST耦合模型均方根误差(RMSE )分别为51.41 kg/hm²、29.52 kg/hm²；Nash-Sutcliffe效率系数(NSE)分别为0.614、0.813；平均绝对误差分别为51.41 kg/hm²、29.52 kg/hm²；相对误差分别为3.02%、1.70%。均表示耦合模型具有较高的模拟精度，说明本文基于CASA-WOFOST耦合估产模式，区域田块尺度上具有较高的模拟精度，可为该条件下的作物估产提供理论支持。

图6

图6 单产模拟结果精度分析

Fig.6 Accuracy analysis of yield simulation results

表4 双山农场基地2016年大豆单产模拟精度分析

Table 4 Analysis results of Soybean yield simulation of in Shuangshan farm in 2016

精度评价指标	CASA模型	CASA-WOFOST模型
均方根误差(RMSE)/(kg/hm²)	61.39	37.31
Nash-Sutcliffe效率系数	0.614	0.813
线性回归判定系数(R²)	0.668 53	0.844 72
绝对误差/(kg/hm²)	51.41	29.52
相对误差	3.02	1.70

5 结语

（1）通过对比分析目前主流估产模型的优势和局限性，寻找更高精度的作物估产模型，基于CASA-WOFOST耦合估产模式成功模拟出黑龙江双山农场2016大豆单产，本文主要是以CASA模型作为耦合模型的主体，通过引入WOFOST模型中的作物生长模拟模块中的干物质分配，代替CASA模型中的收获指数，规避由于收获指数而带来的不必要的误差，提升模型的精度，耦合得到CASA-WOFOST单产估算模式，不仅可以发挥CASA模型的优势，有较高的运行速度，还能发挥WOFOST模型的机理优势，此外还能有效发挥遥感数据的作用。CASA适合在大区域范围内进行作物估产，首先利用CASA模型在田块尺度上进行估产，结果并不理想，因此对CASA模型进行改造，使其与WOFOST模型进行耦合，耦合后的模型在不明显影响模型运行速度的前提下具有更高模拟精度，适合在小区域田块尺度上进行单产估算。此外，根据相邻点线性插值和S-G迭代滤波方法，成功对NDVI进行时间序列重构，实现以天为步长连续监测作物。

（2）采用CASA-WOFOST耦合模型，综合考虑作物外部光能利用等环境因素与作物生长发育全过程，模拟大豆单产，线性回归判定系数(R²)为0.844 72，RMSE为37.31 kg/hm²，Nash-Sutcliffe效率系数为0.813，平均绝对误差与相对误差分别为29.52 kg/hm²、1.70%，预测结果均比利用CASA模型与收获指数模拟单产效果好。说明基于CASA-WOFOST的耦合估产模式，在区域田块尺度上具有较高的模拟精度，可为该条件下的作物估产提供理论支持。

本文所应用的CASA-WOFOST耦合估产模式可以很大程度地提高作物估产精度，为区域农作物估产提供参考，下一步工作有：①在田块尺度上使用更高分辨率的遥感影像；②将模型应用于其他类型作物，探讨模型的适用性；③考虑水分胁迫、温度胁迫以及养分胁迫对模型模拟能力的影响，以提高耦合模型对极端天气的适应性。

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

Cheng

Zhiqiang

, Meng

Jihua

Research Advances and Perspectives on Crop Yield Estimation Models

[J]. Chinese Journal of Eco-Agriculture, 2015, 23(4):402-415.

程志强,蒙继华.

作物单产估算模型研究进展与展望

[J].中国生态农业学报,2015, 23(4):402-415.

[2]

Zhou

Lei

,Li

Gang

, Jia

Dewei

, et al.

Yield Estimation of Winter Wheat in Henan Province based on the Vegetation Photosynthesis Model

[J]. Chinese Journal of Agricultural Resources and Regional Planning, 2017, 38(6):108-115.

周磊,李刚,贾德伟,等.

基于光能利用率模型的河南省冬小麦单产估算研究

[J].中国农业资源与区划, 2017, 38(6): 108-115.