CN112052988A - 耦合多目标优化和集合同化的作物产量估测方法及应用 - Google Patents
耦合多目标优化和集合同化的作物产量估测方法及应用 Download PDFInfo
- Publication number
- CN112052988A CN112052988A CN202010833836.5A CN202010833836A CN112052988A CN 112052988 A CN112052988 A CN 112052988A CN 202010833836 A CN202010833836 A CN 202010833836A CN 112052988 A CN112052988 A CN 112052988A
- Authority
- CN
- China
- Prior art keywords
- gpp
- remote sensing
- wofost
- model
- lai
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 238000000034 method Methods 0.000 title claims abstract description 34
- 238000005457 optimization Methods 0.000 title claims abstract description 13
- GVVPGTZRZFNKDS-JXMROGBWSA-N geranyl diphosphate Chemical compound CC(C)=CCC\C(C)=C\CO[P@](O)(=O)OP(O)(O)=O GVVPGTZRZFNKDS-JXMROGBWSA-N 0.000 claims abstract description 67
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 16
- 238000004088 simulation Methods 0.000 claims abstract description 11
- 238000005070 sampling Methods 0.000 claims abstract description 10
- 238000007781 pre-processing Methods 0.000 claims abstract description 8
- 208000036855 Left sided atrial isomerism Diseases 0.000 claims description 43
- 241000209140 Triticum Species 0.000 claims description 25
- 235000021307 Triticum Nutrition 0.000 claims description 25
- 238000012950 reanalysis Methods 0.000 claims description 15
- 238000002310 reflectometry Methods 0.000 claims description 10
- 238000011160 research Methods 0.000 claims description 7
- 238000013507 mapping Methods 0.000 claims description 3
- 244000068988 Glycine max Species 0.000 claims description 2
- 235000010469 Glycine max Nutrition 0.000 claims description 2
- 240000007594 Oryza sativa Species 0.000 claims description 2
- 235000007164 Oryza sativa Nutrition 0.000 claims description 2
- 240000008042 Zea mays Species 0.000 claims description 2
- 235000005824 Zea mays ssp. parviglumis Nutrition 0.000 claims description 2
- 235000002017 Zea mays subsp mays Nutrition 0.000 claims description 2
- 230000002776 aggregation Effects 0.000 claims description 2
- 238000004220 aggregation Methods 0.000 claims description 2
- 235000005822 corn Nutrition 0.000 claims description 2
- 235000009566 rice Nutrition 0.000 claims description 2
- 230000000052 comparative effect Effects 0.000 description 31
- 230000005855 radiation Effects 0.000 description 6
- 238000012795 verification Methods 0.000 description 5
- 238000004364 calculation method Methods 0.000 description 4
- 235000013339 cereals Nutrition 0.000 description 4
- 229930002875 chlorophyll Natural products 0.000 description 4
- 235000019804 chlorophyll Nutrition 0.000 description 4
- ATNHDLDRLWWWCB-AENOIHSZSA-M chlorophyll a Chemical compound C1([C@@H](C(=O)OC)C(=O)C2=C3C)=C2N2C3=CC(C(CC)=C3C)=[N+]4C3=CC3=C(C=C)C(C)=C5N3[Mg-2]42[N+]2=C1[C@@H](CCC(=O)OC\C=C(/C)CCC[C@H](C)CCC[C@H](C)CCCC(C)C)[C@H](C)C2=C5 ATNHDLDRLWWWCB-AENOIHSZSA-M 0.000 description 4
- 238000004519 manufacturing process Methods 0.000 description 3
- 238000004458 analytical method Methods 0.000 description 2
- 230000008859 change Effects 0.000 description 2
- 150000001875 compounds Chemical class 0.000 description 2
- 238000009826 distribution Methods 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000008878 coupling Effects 0.000 description 1
- 238000010168 coupling process Methods 0.000 description 1
- 238000005859 coupling reaction Methods 0.000 description 1
- 238000005336 cracking Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 230000003203 everyday effect Effects 0.000 description 1
- 238000001914 filtration Methods 0.000 description 1
- 238000003306 harvesting Methods 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
- 230000007246 mechanism Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000029553 photosynthesis Effects 0.000 description 1
- 238000010672 photosynthesis Methods 0.000 description 1
- 230000008569 process Effects 0.000 description 1
- 230000002123 temporal effect Effects 0.000 description 1
- 230000003685 thermal hair damage Effects 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
- G06Q10/00—Administration; Management
- G06Q10/04—Forecasting or optimisation specially adapted for administrative or management purposes, e.g. linear programming or "cutting stock problem"
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
- G06N3/00—Computing arrangements based on biological models
- G06N3/12—Computing arrangements based on biological models using genetic models
- G06N3/126—Evolutionary algorithms, e.g. genetic algorithms or genetic programming
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
- G06Q50/00—Information and communication technology [ICT] specially adapted for implementation of business processes of specific business sectors, e.g. utilities or tourism
- G06Q50/02—Agriculture; Fishing; Forestry; Mining
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Business, Economics & Management (AREA)
- Life Sciences & Earth Sciences (AREA)
- Theoretical Computer Science (AREA)
- Health & Medical Sciences (AREA)
- Economics (AREA)
- Human Resources & Organizations (AREA)
- Biophysics (AREA)
- Strategic Management (AREA)
- General Physics & Mathematics (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Bioinformatics & Computational Biology (AREA)
- General Health & Medical Sciences (AREA)
- Marketing (AREA)
- Tourism & Hospitality (AREA)
- Evolutionary Biology (AREA)
- General Business, Economics & Management (AREA)
- Game Theory and Decision Science (AREA)
- Genetics & Genomics (AREA)
- Primary Health Care (AREA)
- Entrepreneurship & Innovation (AREA)
- Operations Research (AREA)
- Quality & Reliability (AREA)
- Mining & Mineral Resources (AREA)
- Marine Sciences & Fisheries (AREA)
- Animal Husbandry (AREA)
- Agronomy & Crop Science (AREA)
- Physiology (AREA)
- Development Economics (AREA)
- Artificial Intelligence (AREA)
- Biomedical Technology (AREA)
- Computational Linguistics (AREA)
- Data Mining & Analysis (AREA)
- Evolutionary Computation (AREA)
- Molecular Biology (AREA)
- Computing Systems (AREA)
- General Engineering & Computer Science (AREA)
- Mathematical Physics (AREA)
- Software Systems (AREA)
Abstract
本发明属于农业遥感领域，具体公开了耦合多目标优化和集合同化的作物产量估测方法及应用。该方法具体步骤如下：S1、收集数据并预处理；S2、计算LSWI、EVI，再通过VPM模型计算得到遥感GPP；S3、利用MCMC采样方法标定WOFOST模型，之后计算WOFOST模型参数的四分位数；S4、设置WOFOST模型参数的区间范围；基于遥感LAI建立多目标函数，并将其最小化得到优化后的WOFOST模型参数的Pareto解集；S5、将Pareto解集输入WOFOST模型生成模拟GPP集，利用集合同化算法对模拟GPP集与遥感GPP添加高斯噪声生成的观测GPP集进行同化；S6、以同化后的GPP驱动WOFOST模型模拟输出作物产量。本发明估产方式准确性高。
Description
技术领域
本发明属于农业遥感领域，具体涉及一种耦合多目标优化和集合同化的作物产量估测方法及应用。
背景技术
全球气候变化背景下，日益频发的农业气象灾害(干旱、涝渍、高温热害等)，会对粮食产量造成巨大的影响。因此，开展区域尺度的粮食作物产量估算可用于评价气象变化对粮食产量的影响，指导生产改进，对粮食生产具有重要意义。传统农学、气象、统计的产量估测方法，受到理论和技术制约，还存在普适性差、缺乏可解释性等瓶颈问题。遥感与作物模型数据同化兼具遥感动态观测和机理过程明晰的优势，为破解产量精准估测提供了重要方法和解决途径。
但目前作物模型的标定多采用优化单目标函数的方式完成，这种单目标函数难以充分地刻画模拟值与观测值之间的差异，从而导致标定时未能充分地利用观测到的信息，考虑到田间数据的稀缺性，这无疑也是一种信息浪费。另外，现有的集合同化方法主要是对标定所得的一套参数添加高斯误差扰动来生成参数集合，在实际中常凭经验确定所添加的误差的大小，这使得标定的结果将因为操作人员的主观性而具有更大的不确定性。
因此，有必要提供一种新的作物产量估测方法及应用以解决现有技术中的问题。
发明内容
为解决现有技术中存在的问题，本发明提供一种采用耦合多目标优化和集合同化的产量估测方法，以综合利用观测数据，减少估测误差。
为了实现该目的，本发明的技术方案如下：
一种作物产量估测方法，其通过包含耦合多目标优化和集合同化步骤的方法来实现，具体步骤如下：
S1、收集研究区的作物分类数据；并将收集到的作物生育期内的遥感LAI、遥感反射率数据以及气象再分析数据分别进行预处理；
S2、利用步骤S1获得的预处理后的遥感反射率数据计算LSWI、EVI，再通过VPM模型计算得到遥感GPP；
S3、基于田间观测的LAI时间序列以及作物产量，利用MCMC采样方法标定WOFOST模型，之后依据MCMC采样方法获得的采样结果计算WOFOST模型参数的四分位数；
S4、以步骤S3获得的WOFOST模型参数的第一四分位数为下界、第三四分位数为上界设置WOFOST模型参数的区间范围；基于步骤S1获得的预处理后的遥感LAI建立多目标函数，采用NSGA-II方法最小化所述多目标函数得到优化后的WOFOST模型参数的Pareto解集；
S5、将步骤S4得到的WOFOST模型参数的Pareto解集输入WOFOST模型生成模拟GPP集，将步骤S2得到的遥感GPP添加高斯噪声生成观测GPP集，利用集合同化算法对模拟GPP集与观测GPP集进行同化，得到同化后的GPP；
S6、基于步骤S5得到的同化后的GPP驱动WOFOST模型模拟输出作物产量，并完成空间制图。
其中，步骤S1中，所述作物分类数据是指研究区内待估测产量的作物(如冬小麦)的空间分布图，该数据将用于过滤遥感影像中非此种作物的像元，使后续的研究可以针对研究区内的这种特定待估产作物。
所述气象再分析数据是指对历史的多种气象数据进行再分析形成的更完善的气象数据，比如全球NCEP再分析数据，ECMWF再分析数据，MERRA再分析数据或JRA再分析数据。
步骤S2中，EVI和LSWI的计算方法分别见公式(5)、(6)：
式中，ρnir，ρred，ρblue，以及ρswir分别表示遥感观测到的近红外，红光，蓝光以及短波红外波段的反射率。
VPM模型的计算方法见公式(7)到(12)：
GPP＝εg×APARchl (7)
FPARchl＝(EVI-0.1)×1.25 (8)
APARchl＝FPARchl×PAR (9)
εg＝ε0×Tscalar×Wscalar (10)
式中，GPP为总初级生产力，εg为待估产作物的光能利用率(C3作物设为0.42，C4作物设为0.63)，ε0为待估产作物的最大光能利用率，APARchl为叶绿素吸收的光合有效辐射总量，PAR为光合有效辐射(由气象再分析数据获得),FPARchl为叶绿素吸收的光合有效辐射的比例。Tscalar与Wscalar分别为温度和水分影响因子。T是白天平均气温(由气象再分析数据获得)，Tmax、Tmin以及Topt分别为光合作用最大、最小及最适温度(具体应用时其值分别设为-1，48以及30)。LSWImax为非降雪时期的年最大LSWI。
步骤S3中，LAI时间序列是指在作物不同的生育期观测的LAI。作物产量是指收割作物时的作物产量。
步骤S4中，所述的区间范围为最小化多目标函数时参数的可行域。
步骤S5中，添加高斯噪声是指对于遥感GPP添加一个均值为0，方差为σ2的高斯噪声，使得观测GPP集的均值为遥感GPP值，方差为σ2。如果有地面实测GPP可以验证遥感GPP，那么可设σ为遥感GPP与地面实测GPP之间的RMSE值；如果没有地面实测GPP，则可设置σ为遥感GPP值的10％。
本发明中，步骤S1中所述预处理的方式包括降尺度、插值、剔除异常值、聚合或数据填充中的一种或多种，以使得处理后的数据为时空连续的数据。
本发明中，步骤S4中所述基于步骤S1获得的预处理后的遥感LAI建立多目标函数的步骤具体采用公式(1)至公式(3)：
式中，θ表示待优化的WOFOST模型的参数，K表示获取遥感LAI的总次数，
分别表示第k次以及第k-1次观测的遥感LAI，
分别表示与第k次以及第k-1次观测的遥感LAI相对应的WOFOST模型模拟的LAI。
本发明中f1(θ)、f2(θ)、f3(θ)共同构成NSGA-II方法最小化的多目标函数，其中f1(θ)、f2(θ)分别适用于LAI较大时、LAI较小时遥感观测的LAI与模型模拟的LAI的比较，而f3(θ)则用于二者对全生育期的LAI曲线形态的比较。在实际应用时，多目标函数不限于以上三个，允许以不同的统计量刻画模拟值与观测值的差异。
本发明中，步骤S5中所述集合同化算法见公式(4)：
式中，
表示WOFOST生成的t时刻的模拟GPP集，
表示步骤S2获得的t时刻的遥感GPP,
表示步骤S5中t时刻的模拟GPP集的第i个成员，
为集合同化后的t时刻的GPP集的第i个成员。Var表示求取变量的方差。
本发明中，所述作物为小麦、水稻、玉米或大豆中一种，优选为冬小麦。
本发明还提供一种上述方法在估测作物产量中的应用。
本发明的有益效果至少在于：
本发明采用了先依据田间实测数据利用MCMC算法对作物模型完成初次标定，再利用多目标优化算法结合遥感观测LAI数据再次标定作物模型，得到模型参数的Pareto解集，最后利用Pareto解集生成模拟的GPP集合并通过集合同化的方式与遥感观测的GPP进行同化，这种方法在标定模型时综合利用了田间与遥感观测数据，且基于多目标优化所得的Pareto解集使得同化时生成模拟的GPP集合更具有科学性，减少了主观经验导致的同化误差。
附图说明
图1为本发明实施例1和对比例1、对比例2、对比例3的冬小麦产量估测结果比较图，其中，图1(a)为实施例1的产量估测结果图，图1(b)为对比例1的产量估测结果图，图1(c)为对比例2的产量估测结果图，图1(d)为对比例3的产量估测结果图。
图2为本发明实施例1和对比例1对比例2、对比例3的产量估测结果验证图，其中图2(a)为实施例1的产量估测结果验证图，图2(b)为对比例1的产量估测结果验证图，图2(c)为对比例2的产量估测结果验证图，图2(d)为对比例3的产量估测结果验证图。
具体实施方式
下面结合附图和实施例，对本发明的具体实施方式作进一步详细描述。以下实施例用于说明本发明，但不用来限制本发明的范围。
实施例1
本实施例中，以本发明的技术方案对2017年河北省冀州市的冬小麦产量进行估测。具体包括：
S1、收集2017年研究区的作物分类数据，并收集2016、2017年冬小麦生育期内的遥感LAI、遥感反射率数据以及气象再分析数据分别进行预处理；
其中，所述作物分类数据是指2017年河北省冀州市的冬小麦空间分布图，该数据将用于过滤遥感影像中非冬小麦的像元，后续的研究将针对冀州市内的冬小麦作物。
其中，所述遥感LAI、遥感反射率数据分别采用MOD15A2H、MOD09A1，其时间分辨率分别皆为8天，空间分辨率皆为500米。
其中，所述气象再分析数据是指对历史的多种气象数据进行再分析形成的更完善的气象数据，本实施例选用NCEP再分析数据，其时间分辨率为每天。
其中，所述预处理指将NCEP再分析数据插值为每8天的数据，使得插值后的NCEP再分析数据与MOD15A2H、MOD09A1的时间分辨率统一。
S2、利用S1所得预处理后的遥感反射率数据计算LSWI、EVI，得到时间分辨率为8天，空间分辨率为500米的EVI、LSWI数据；再通过VPM模型计算得到遥感GPP；
其中，遥感反射率数据指预处理后的MOD09A1。
其中，EVI和LSWI的计算方法分别为公式(5)、(6)：
式中，ρnir，ρred，ρblue，以及ρswir分别表示遥感观测到的近红外，红光，蓝光以及短波红外波段的反射率。
其中，VPM模型的计算方法为公式(7)到(12)：
GPP＝εg×APARchl (7)
FPARchl＝(EVI-0.1)×1.25 (8)
APARchl＝FPARchl×PAR (9)
εg＝ε0×Tscalar×Wscalar (10)
式中，GPP为总初级生产力，εg为冬小麦的光能利用率(冬小麦是C3作物，所以设为0.42)，ε0为冬小麦的最大光能利用率，APARchl是叶绿素吸收的光合有效辐射总量，PAR是光合有效辐射，由气象再分析数据获得(本实施例指预处理后的NCEP再分析数据中的短波辐射数据),FPARchl是叶绿素吸收的光合有效辐射的比例。Tscalar与Wscalar是温度和水分影响因子。T是白天平均气温，由气象再分析数据获得，Tmax、Tmin以及Topt是光合作用最大、最小及最适温度，其值分别设为-1，48以及30。LSWImax是非降雪时期的年最大LSWI。
S3、基于田间观测的LAI时间序列以及作物产量，利用MCMC采样方法标定WOFOST模型，之后依据采样结果计算WOFOST模型参数的四分位数；
其中，LAI时间序列是指在2016年10月至2017年6月观测的冬小麦在不同的生育期的LAI。作物产量是指2017年6月初收割时冬小麦的产量。
本实施例对河北省冀州市选取了两个样方，每个样方的大小为100米×100米，每个样方选取了5个样点，每个样点的大小为1米×1米。
S4、以S3获得的WOFOST模型参数的第一四分位数为下界、第三四分位数为上界设置WOFOST模型参数的区间范围；基于S1获得的预处理后的遥感LAI建立多目标函数，采用NSGA-II方法最小化上述多目标函数得到优化后的WOFOST模型参数的Pareto解集；
其中，遥感LAI指预处理后的MOD15A2H。
其中，所述的区间范围为最小化多目标函数时参数的可行域。
其中，所述基于S1获得的预处理后的遥感LAI建立多目标函数采用公式(1)至公式(3)：
式中，θ表示待优化的WOFOST模型的参数，K表示获取遥感LAI的总次数，
分别表示第k次以及第k-1次观测的遥感LAI，
分别表示与第k次以及第k-1次观测的遥感LAI相对应的WOFOST模型模拟的LAI。f1(θ)、f2(θ)、f3(θ)共同构成NSGA-II方法最小化的多目标函数，其中f1(θ)、f2(θ)分别适用于LAI较大时、LAI较小时遥感观测的LAI与模型模拟的LAI的比较，而f3(θ)则用于二者对全生育期的LAI曲线形态的比较。
S5、将S4得到的WOFOST模型参数的Pareto解集(解集包含50组参数)输入WOFOST模型生成模拟GPP集，将S2得到的遥感GPP添加高斯噪声生成观测GPP集，利用集合同化算法对模拟GPP集与观测GPP集进行同化，得到同化后的GPP；
其中，添加高斯噪声是指对于遥感GPP添加一个均值为0，方差为σ2的高斯噪声，使得观测GPP集的均值为遥感GPP值，方差为σ2。本实施例设置为遥感GPP值的10％。
其中，所述集合同化算法见公式(4)：
式中，
表示WOFOST生成的t时刻的模拟GPP集，
表示步骤S2获得的t时刻的遥感GPP,
表示步骤S5中t时刻的模拟GPP集的第i个成员，
为集合同化后的t时刻的GPP集的第i个成员。Var表示求取变量的方差。
S6、基于S5得到的同化后的GPP驱动WOFOST模型模拟输出冬小麦产量，并完成空间制图。
以该耦合多目标优化和集合同化GPP的方式获得的产量估测结果的平均值做图，得到的结果见图1中的图1(a)实施例1的产量估测结果图。将该估测结果与实测产量进行比较，得到产量估测结果验证图，见图2(a)。
对比例1
本对比例中，与实施例1一样对2017年河北省冀州市的冬小麦产量进行估测。依据田间实测数据利用MCMC算法对WOFOST模型完成标定后，再同化GPP估测产量结果。具体采用实施例1的步骤S1-S3后，直接执行步骤S5、步骤S6，但步骤S5中“将S4得到的WOFOST模型参数的Pareto解集输入WOFOST模型生成模拟GPP集”替换为“从S3得到的WOFOST模型参数的四分位数限定的参数范围中随机采样50组参数并输入WOFOST模型生成模拟GPP集”(即省略实施例1的步骤S4，变化步骤S5中输入WOFOST模型的参数，获得区域冬小麦产量同化后的结果)。以该方式获得的估测结果的平均值做图，得到的结果见图1中的图1(b)对比例1的产量估测结果图。将该估测结果与实测产量进行比较，得到产量估测结果验证图，见图2(b)。
对比例2
本对比例中，与实施例1一样对2017年河北省冀州市的冬小麦产量进行估测。先依据田间实测数据利用MCMC算法对WOFOST模型完成初次标定，再利用多目标优化算法结合遥感观测LAI数据再次标定WOFOST模型进而估测产量。具体采用实施例1的步骤S1-S4后，直接将步骤S4得到的WOFOST模型参数的Pareto解集输入WOFOST模型模拟冬小麦的产量(即省略实施例1的步骤S5与步骤S6，获得区域冬小麦产量同化前的结果)。以该方式获得的产量估测结果的平均值做图，得到的结果见图1中的图1(c)对比例2的产量估测结果图。将该估测结果与实测产量进行比较，得到产量估测结果验证图，见图2(c)。
对比例3
本对比例中，与实施例1一样对2017年河北省冀州市的冬小麦产量进行估测。依据田间实测数据利用MCMC算法对WOFOST模型完成标定后直接估测产量。具体采用实施例1的步骤S1-S3后，直接从步骤S3得到的WOFOST模型参数的四分位数限定的参数范围中随机采样50组参数并输入WOFOST模型模拟冬小麦的产量(即省略实施例1的步骤S4、步骤S5与步骤S6，获得区域冬小麦产量同化前的结果)。以该方式获得的产量估测结果的平均值做图，得到的结果见图1中的图1(d)对比例3的产量估测结果图。将该估测结果与实测产量进行比较，得到产量估测结果验证图，见图2(d)。
从图1和图2本发明实施例1和对比例1、对比例2、对比例3的冬小麦产量估测结果比较图和产量估测结果验证图中可以看出，实施例1的空间差异性(图1(a))最明显，其产量范围集中在4500kg/ha-7500kg/ha，且与实测产量的相关性最强(R2为0.3388)、精度最高(RMSE为628kg/ha)(图2(a))；对比例1的空间差异性(图1(b))以及与实测产量的相关性、精度(图2(b))弱于实施例1，其产量范围集中在4500kg/ha-7000kg/ha，R2为0.2149，RMSE为972kg/ha；对比例2的空间差异性(图1(c))弱于对比例1，其产量范围集中在4500kg/ha-6500kg/ha，虽然其精度稍高于对比例1(RMSE为920kg/ha)，但其与实测产量的相关性弱于实施例1(R2为0.1566)(图2(c))；对比例3的空间差异性(图1(d))以及与实测产量的相关性、精度(图2(d))最弱，其产量范围集中在4500kg/ha-6000kg/ha且大部分低于5500kg/ha，R2为0.1189，RMSE为1238kg/ha。
通过对比例2与对比例3的对比说明先依据田间实测数据利用MCMC算法对作物模型完成初次标定，再利用多目标优化算法结合遥感观测LAI数据再次标定作物模型进而估测产量较仅依据有限数目的样方(对比例中为2个样方)的田间实测数据利用MCMC算法对作物模型完成标定后直接估测产量的精度更高，相关性更强。
实施例1与对比例1估测产量的结果优于对比例2和对比例3，这说明了集合同化GPP的方式较直接利用标定的模型进行产量估测的精度更高，相关性更强。
实施例1与对比例1的对比说明，耦合多目标优化和集合同化的估测产量比不耦合多目标优化的集合同化的估测产量的精度更高，相关性更强。
以上实施方式仅用于说明本发明，而并非对本发明的限制，有关技术领域的普通技术人员，在不脱离本发明的精神和范围的情况下，还可以做出各种变化和变型，因此所有等同的技术方案也属于本发明的范畴，本发明的专利保护范围应由权利要求限定。
Claims (6)
1.一种作物产量估测方法，其特征在于，通过包含耦合多目标优化和集合同化步骤的方法来实现，具体步骤如下：
S1、收集研究区的作物分类数据；并将收集到的作物生育期内的遥感LAI、遥感反射率数据以及气象再分析数据分别进行预处理；
S2、利用步骤S1获得的预处理后的遥感反射率数据计算LSWI、EVI，再通过VPM模型计算得到遥感GPP；
S3、基于田间观测的LAI时间序列以及作物产量，利用MCMC采样方法标定WOFOST模型，之后依据MCMC采样方法获得的采样结果计算WOFOST模型参数的四分位数；
S4、以步骤S3获得的WOFOST模型参数的第一四分位数为下界、第三四分位数为上界设置WOFOST模型参数的区间范围；基于步骤S1获得的预处理后的遥感LAI建立多目标函数，采用NSGA-II方法最小化所述多目标函数得到优化后的WOFOST模型参数的Pareto解集；
S5、将步骤S4得到的WOFOST模型参数的Pareto解集输入WOFOST模型生成模拟GPP集，将步骤S2得到的遥感GPP添加高斯噪声生成观测GPP集，利用集合同化算法对模拟GPP集与观测GPP集进行同化，得到同化后的GPP；
S6、基于步骤S5得到的同化后的GPP驱动WOFOST模型模拟输出作物产量，并完成空间制图。
2.如权利要求1所述的方法，其特征在于，步骤S1中所述预处理的方式包括降尺度、插值、剔除异常值、聚合或数据填充中的一种或多种。
3.如权利要求1所述的方法，其特征在于，步骤S4中所述基于步骤S1获得的预处理后的遥感LAI建立多目标函数的步骤具体采用公式(1)至公式(3)：
式中，θ表示待优化的WOFOST模型的参数，K表示获取遥感LAI的总次数，
分别表示第k次以及第k-1次观测的遥感LAI，
分别表示与第k次以及第k-1次观测的遥感LAI相对应的WOFOST模型模拟的LAI。
4.如权利要求1所述的方法，其特征在于，步骤S5中所述集合同化算法见公式(4)：
式中，GPPt sim表示WOFOST生成的t时刻的模拟GPP集，GPPt obs表示步骤S2获得的t时刻的遥感GPP,
表示步骤S5中t时刻的模拟GPP集的第i个成员，
为集合同化后的t时刻的GPP集的第i个成员。
5.如权利要求1所述的方法，其特征在于，所述作物为小麦、水稻、玉米或大豆中的一种。
6.权利要求1-5任一项所述的方法在估测作物产量中的应用。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010833836.5A CN112052988B (zh) | 2020-08-18 | 2020-08-18 | 耦合多目标优化和集合同化的作物产量估测方法及应用 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010833836.5A CN112052988B (zh) | 2020-08-18 | 2020-08-18 | 耦合多目标优化和集合同化的作物产量估测方法及应用 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN112052988A true CN112052988A (zh) | 2020-12-08 |
| CN112052988B CN112052988B (zh) | 2024-03-22 |
Family
ID=73599634
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202010833836.5A Active CN112052988B (zh) | 2020-08-18 | 2020-08-18 | 耦合多目标优化和集合同化的作物产量估测方法及应用 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN112052988B (zh) |
Cited By (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112991247A (zh) * | 2021-03-04 | 2021-06-18 | 河南省气象科学研究所 | 冬小麦蒸散遥感反演及与作物模型同化方法 |
| CN115753625A (zh) * | 2022-11-02 | 2023-03-07 | 中国农业大学 | 区域作物的产量估测方法、装置、电子设备及存储介质 |
| CN116976516A (zh) * | 2023-08-01 | 2023-10-31 | 中国科学院空天信息创新研究院 | 一种干旱地区作物单产早期预测方法 |
| CN117314184A (zh) * | 2023-09-19 | 2023-12-29 | 沈阳农业大学 | 年总初级生产力时空变化结果的确定方法、装置及设备 |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2010166851A (ja) * | 2009-01-22 | 2010-08-05 | Chiharu Hongo | 作物の収量予測方法及びその装置 |
| CN104134095A (zh) * | 2014-04-17 | 2014-11-05 | 中国农业大学 | 一种基于尺度转换和数据同化的农作物产量估测方法 |
| CN108509836A (zh) * | 2018-01-29 | 2018-09-07 | 中国农业大学 | 双极化合成孔径雷达与作物模型数据同化的作物估产方法 |
| CN110222870A (zh) * | 2019-05-05 | 2019-09-10 | 中国农业大学 | 同化卫星荧光数据与作物生长模型的区域冬小麦估产方法 |
| CN110633841A (zh) * | 2019-08-13 | 2019-12-31 | 中国农业大学 | 基于集合采样的省级范围地块尺度快速数据同化产量预测方法 |
-
2020
- 2020-08-18 CN CN202010833836.5A patent/CN112052988B/zh active Active
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2010166851A (ja) * | 2009-01-22 | 2010-08-05 | Chiharu Hongo | 作物の収量予測方法及びその装置 |
| CN104134095A (zh) * | 2014-04-17 | 2014-11-05 | 中国农业大学 | 一种基于尺度转换和数据同化的农作物产量估测方法 |
| CN108509836A (zh) * | 2018-01-29 | 2018-09-07 | 中国农业大学 | 双极化合成孔径雷达与作物模型数据同化的作物估产方法 |
| CN110222870A (zh) * | 2019-05-05 | 2019-09-10 | 中国农业大学 | 同化卫星荧光数据与作物生长模型的区域冬小麦估产方法 |
| CN110633841A (zh) * | 2019-08-13 | 2019-12-31 | 中国农业大学 | 基于集合采样的省级范围地块尺度快速数据同化产量预测方法 |
Non-Patent Citations (1)
| Title |
|---|
| 黄健熙;黄海;马鸿元;李颖;侯英雨;何亮;朱德海;: "基于MCMC方法的WOFOST模型参数标定与不确定性分析", 农业工程学报, no. 16, 22 August 2018 (2018-08-22) * |
Cited By (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112991247A (zh) * | 2021-03-04 | 2021-06-18 | 河南省气象科学研究所 | 冬小麦蒸散遥感反演及与作物模型同化方法 |
| CN115753625A (zh) * | 2022-11-02 | 2023-03-07 | 中国农业大学 | 区域作物的产量估测方法、装置、电子设备及存储介质 |
| CN115753625B (zh) * | 2022-11-02 | 2024-04-19 | 中国农业大学 | 区域作物的产量估测方法、装置、电子设备及存储介质 |
| CN116976516A (zh) * | 2023-08-01 | 2023-10-31 | 中国科学院空天信息创新研究院 | 一种干旱地区作物单产早期预测方法 |
| CN116976516B (zh) * | 2023-08-01 | 2024-03-15 | 中国科学院空天信息创新研究院 | 一种干旱地区作物单产早期预测方法 |
| CN117314184A (zh) * | 2023-09-19 | 2023-12-29 | 沈阳农业大学 | 年总初级生产力时空变化结果的确定方法、装置及设备 |
| CN117314184B (zh) * | 2023-09-19 | 2024-04-30 | 沈阳农业大学 | 年总初级生产力时空变化结果的确定方法、装置及设备 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN112052988B (zh) | 2024-03-22 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN112052988B (zh) | 耦合多目标优化和集合同化的作物产量估测方法及应用 | |
| CN109829234A (zh) | 一种基于高分辨率遥感数据和作物模型的跨尺度高精度动态作物长势监测和估产方法 | |
| Meng et al. | Corn yield forecasting in northeast china using remotely sensed spectral indices and crop phenology metrics | |
| CN110309985A (zh) | 一种农作物产量预测方法及系统 | |
| CN110633841B (zh) | 基于集合采样的省级范围地块尺度数据同化产量预测方法 | |
| CN106485002B (zh) | 在复杂地形气候区域估算甘蔗潜在产量的方法 | |
| Wang et al. | Monitoring maize growth on the North China Plain using a hybrid genetic algorithm-based back-propagation neural network model | |
| CN107122739B (zh) | 基于Extreme数学模型重构VI时间序列曲线的作物估产模型 | |
| Wang et al. | Increased uncertainty in simulated maize phenology with more frequent supra-optimal temperature under climate warming | |
| CN116308866B (zh) | 一种基于冠层反射光谱的水稻穗生物量估算方法及系统 | |
| CN110245444B (zh) | 一种基于作物对环境的响应与适应机制的发育期模拟方法 | |
| CN110321774B (zh) | 农作物灾情评估方法、装置、设备及计算机可读存储介质 | |
| CN110793921A (zh) | 一种抽雄授粉期玉米洪涝灾害遥感监测与评估方法与系统 | |
| CN116108318A (zh) | 基于无人机多光谱影像的油菜氮肥推荐追施量计算的方法 | |
| Bai et al. | Growth simulation and yield prediction for perennial jujube fruit tree by integrating age into the WOFOST model | |
| CN108537679B (zh) | 遥感与作物模型融合的区域尺度作物出苗日期估算方法 | |
| CN115730523A (zh) | 一种基于深度学习的区域尺度作物产量近实时预测方法 | |
| CN110516943B (zh) | 一种基于地表温度的春灌期灌溉面积动态监测遥感方法 | |
| CN116485757A (zh) | 一种冬小麦全氮含量预测方法 | |
| CN115952421A (zh) | 耦合生态过程模型与机器学习算法的高精度时空模拟方法 | |
| CN116629453B (zh) | 一种适用于作物全生育时期的遥感估产方法 | |
| Viola et al. | Olive yield and future climate forcings | |
| CN111915096A (zh) | 基于作物模型、遥感数据和气候预测信息的作物产量早期预报技术 | |
| Cappelli et al. | GLORIFY: A new forecasting system for rice grain quality in Northern Italy | |
| CN112380497B (zh) | 用于区域玉米的秸秆系数估算方法及系统 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |