CN115879836B - 一种耦合机理模型的大豆作物遥感大区域快速估产方法 - Google Patents
一种耦合机理模型的大豆作物遥感大区域快速估产方法 Download PDFInfo
- Publication number
- CN115879836B CN115879836B CN202310213131.7A CN202310213131A CN115879836B CN 115879836 B CN115879836 B CN 115879836B CN 202310213131 A CN202310213131 A CN 202310213131A CN 115879836 B CN115879836 B CN 115879836B
- Authority
- CN
- China
- Prior art keywords
- remote sensing
- model
- data
- soybean
- ndvi
- 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.)
- Active
Links
- 235000010469 Glycine max Nutrition 0.000 title claims abstract description 79
- 244000068988 Glycine max Species 0.000 title claims abstract description 79
- 238000000034 method Methods 0.000 title claims abstract description 32
- 230000008878 coupling Effects 0.000 title claims abstract description 22
- 238000010168 coupling process Methods 0.000 title claims abstract description 22
- 238000005859 coupling reaction Methods 0.000 title claims abstract description 22
- 230000007246 mechanism Effects 0.000 title claims abstract description 20
- 238000011160 research Methods 0.000 claims abstract description 15
- 239000002689 soil Substances 0.000 claims abstract description 12
- 230000005540 biological transmission Effects 0.000 claims abstract description 5
- 238000004364 calculation method Methods 0.000 claims abstract description 5
- 230000005855 radiation Effects 0.000 claims abstract description 5
- 230000001351 cycling effect Effects 0.000 claims abstract description 4
- 238000007781 pre-processing Methods 0.000 claims abstract description 4
- 238000002310 reflectometry Methods 0.000 claims description 19
- 238000005070 sampling Methods 0.000 claims description 10
- 239000000203 mixture Substances 0.000 claims description 9
- 238000005457 optimization Methods 0.000 claims description 6
- 239000002245 particle Substances 0.000 claims description 6
- 230000003595 spectral effect Effects 0.000 claims description 6
- 238000001228 spectrum Methods 0.000 claims description 4
- 238000012937 correction Methods 0.000 claims description 3
- 238000009826 distribution Methods 0.000 claims description 3
- 239000011159 matrix material Substances 0.000 claims description 3
- 230000008569 process Effects 0.000 claims description 3
- 238000005316 response function Methods 0.000 claims description 3
- 238000013507 mapping Methods 0.000 abstract description 3
- 238000007726 management method Methods 0.000 description 8
- 241000196324 Embryophyta Species 0.000 description 3
- 230000035945 sensitivity Effects 0.000 description 3
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 description 2
- 238000010276 construction Methods 0.000 description 2
- 238000011161 development Methods 0.000 description 2
- 230000018109 developmental process Effects 0.000 description 2
- 238000005516 engineering process Methods 0.000 description 2
- 238000011156 evaluation Methods 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 230000003287 optical effect Effects 0.000 description 2
- 239000000049 pigment Substances 0.000 description 2
- 238000012545 processing Methods 0.000 description 2
- 230000004044 response Effects 0.000 description 2
- 238000012360 testing method Methods 0.000 description 2
- 238000012546 transfer Methods 0.000 description 2
- 238000012271 agricultural production Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000002457 bidirectional effect Effects 0.000 description 1
- 238000013135 deep learning Methods 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 238000004836 empirical method Methods 0.000 description 1
- 239000003337 fertilizer Substances 0.000 description 1
- 230000006870 function Effects 0.000 description 1
- 238000003306 harvesting Methods 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
- 230000029553 photosynthesis Effects 0.000 description 1
- 238000010672 photosynthesis Methods 0.000 description 1
- 238000011002 quantification Methods 0.000 description 1
- 239000002881 soil fertilizer Substances 0.000 description 1
- 230000017260 vegetative to reproductive phase transition of meristem Effects 0.000 description 1
Images
Classifications
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
- Y02A40/00—Adaptation technologies in agriculture, forestry, livestock or agroalimentary production
- Y02A40/10—Adaptation technologies in agriculture, forestry, livestock or agroalimentary production in agriculture
Landscapes
- Image Processing (AREA)
- Image Analysis (AREA)
Abstract
本发明涉及农业遥感领域，具体提供了一种耦合机理模型的大豆作物遥感大区域快速估产方法，包括如下步骤：S1：对研究区内的多源遥感卫星数据进行数据预处理；S2：通过对遥感影像中的NDVI数据自动分级来获取研究区内大豆作物不同生长状态的空间点位；S3：构建大豆作物生长模型的品质参数、土壤参数和管理参数；S4：将大豆作物的生长模式与辐射传输模型进行耦合；S5：对所有同化点，通过步骤S4进行依次循环并得到所有同化点的产量数据，并通过同化点的产量数据与时序遥感影像的指数数据构建多源回归模型，完成大区域大豆作物产量的计算。本发明估算效率和精度高，地域应用范围广，可快速实现大区域范围下的大豆作物产量制图。
Description
技术领域
本发明涉及农业遥感技术领域，尤其涉及一种耦合机理模型的大豆作物遥感大区域快速估产方法。
背景技术
大豆作为世界上最重要的油料作物之一，在农业生产中发挥着至关重要的作用，大豆收获前的产量预测与农户、农业企业和国家战略规划密切相关，准确的产量预测可以帮助个体农民改善田间管理，及时提高产量，并向保险公司提供相关的或有价值的信息，因此，需要对大豆估产方法进行设计。
目前，国内外遥感作物估产模型主要分为遥感植被指数统计估产模型、光能利用率模型、遥感和作物生长耦合模型等。1）遥感植被指数统计估产模型：主要基于遥感机理和经验方法来建立获取卫星遥感数据与作物产量之间的定量关系，未考虑作物本身的生长模式以及气象条件因素，通常只能得到大致的产量趋势而无法保证较高的精度，只适用于局部地区，时空可移植性较低；2）光能利用率模型：是基于遥感资料估算植被生产力的主要方法，它对光合作用做了理论上的简化和抽象，虽然对作物生长机理有一定的涉及，但部分参数的量化需要耗费较长的时间并进行野外测量来获取，且光能利用率模型很难考虑到人为因素的影响；3）遥感和作物生长耦合模型：通过将卫星遥感数据集成到植物作物生理生长模型中，主要是对影像中每个像元进行同化，具有时空可拓展性，但仅仅依靠这种逐像元同化的方式，计算时间过长，难以大面积推广。
综上所述，如何设计一种新的技大豆估产模型，以解决作物估产中存在的局限性问题并提高估算效率、估算精度及地域应用范围，是当下亟需解决的问题。
发明内容
本发明为解决上述问题，提供了一种耦合机理模型的大豆作物遥感大区域快速估产方法，通过结合时序遥感影像，气象数据参数，并将遥感影像的时序指数数据与同化产量数据进行多元回归模型构建，可快速实现大区域范围下的大豆作物产量制图。
为达到上述目的，本发明提出如下技术方案：一种耦合机理模型的大豆作物遥感大区域快速估产方法，包括如下步骤：
S1：对研究区内的多源遥感卫星数据进行数据预处理；
S11：获取多源遥感卫星数据的时序遥感影像；
S12：获取时序遥感影像的反射率数据和结构性参数；
S2：通过对遥感影像中的NDVI数据自动分级来获取研究区内大豆作物不同生长状态的空间点位；
S21：计算遥感影像中的NDVI数据；
S22：引入混合高斯模型GMM对NDVI进行聚类并自动化采集不同长势下大豆的发射率同化点；
S3：构建大豆作物生长模型的品质参数、土壤参数和管理参数；
S4：将大豆作物的生长模式与辐射传输模型进行耦合；
S41：引入DSSAT模型，利用粒子群算法根据步骤S3中的品质参数、土壤参数和管理参数的步长完成DSSAT模型的初始化；
S42：引入PROSAIL模型，根据初始化参数运行DSSAT模型，并得到大豆作物在整个生长周期内的LAI值和相应时间状态下的理论LAI数据序列，再将理论LAI数据序列及步骤S12中获取的结构性参数作为输入驱动PROSAIL模型，并得到该大豆同化点的光谱反射率曲线；
S43：计算当前大豆同化点红波段与近红外波段的反射率值，对DSSAT模型的参数进行优化并得到DSSAT模型的最佳参数，此时DSSAT模型运行结果中的最终产量为当前大豆同化点的最优产量；
S5：对步骤S2中得到的所有同化点，通过步骤S4进行依次循环并得到所有同化点的产量数据，并通过同化点的产量数据与时序遥感影像的指数数据构建多源回归模型，完成大区域范围内大豆作物产量的计算；
S51：计算遥感影像中的EVI数据、SIPI数据；并通过NDVI、EVI和SIPI数据计算该同化点的综合植被指数；
S52：构建多源回归模型并回归同化点的产量与综合植被指数的关系；
S53：构建研究区内时序遥感影像的NDVI、EVI与SIPI的指数集合，将步骤S52中的回归系数应用在研究区内的时序遥感影像中，完成大区域大豆作物的产量制图。
优选的，步骤S21中通过公式（1）计算遥感影像中的NDVI数据：
(1)
其中，
R nir 与
R r 分别为近红外波段与红波段的反射率值。
优选的，步骤S22中引入高斯混合模型GMM对NDVI进行聚类的过程如下：
S221：设置初始化高斯混合模型GMM的成分个数
K，并聚类不同长势下的NDVI；此处的
K为3，即大豆的长势分为长势好、长势中和长势差3个等级；
S222：计算大豆作物分布栅格影像中每个点的NDVI值属于各个高斯模型的概率，即通过公式（2）计算后验概率：
(2)
其中，
K为模型的个数，此处
K=3，
δ k 为属于第
k个高斯的概率，
δ k
＞0；为第
k个高斯的概率密度，
μ k 为均值向量, Σ k 为协方差矩阵；
S223：为每个NDVI点分别计算由该混合模型内每个高斯分量生成的概率，更新
δ k 、
μ k 、Σ k 参数，使得数据点属于长势好、长势中和长势差中某一类别的概率最大化；
S224：通过EM算法重复迭代步骤S222和步骤S223直到模型收敛；此时得到不同长势下NDVI的分级间隔点
θ 1 ，
θ 2 ；
S225：设NDVI1、NDVI2、NDVI3分别为长势差的点、长势中的点和长势好的点的NDVI值；根据需要同化点的数量，并根据阈值随机选取不同类别的同化采样点；NDVI1、NDVI2、NDVI3分别通过公式（3）、公式（4）和公式（5）进行判断：
。
优选的，步骤S3中的品质参数包括CSDL、PPSEN和EM-FL，土壤参数包括SLPF，管理参数包括种植时间、种植密度和施肥量。
优选的，步骤S43中当前大豆同化点红波段与近红外波段的反射率值通过公式（6）和公式（7）得到：
(6)
(7)
其中，与分别为PROSAIL模型模拟的当前大豆同化点的近红外波段反射率值与红波段反射率值；与为遥感卫星数据源中近红外波段与红波段的光谱响应函数；ρo为PROSAIL模型模拟的当前大豆同化点的光谱曲线。
优选的，设
R nir 与
R r 分别为采样点遥感影像中的近红外波段与红波段的反射率值；将
R nir 、
R r 、和代入公式（8）中并求出误差损失
L：
(8)
步骤S43中对DSSAT模型的参数进行优化具体为：根据粒子群优化算法重复公式（6）-公式（8）直至损失收敛，此时得到的DSSAT模型参数为最佳参数。
优选的，步骤S51中的EVI数据和SIPI数据分别通过公式（9）和公式（10）得到：
(9)
(10)
其中，
R b 为采样点遥感影像中的蓝波段的反射率值。
优选的，步骤S51中的综合植被指数通过公式（11）得到：
(11)
其中，
F t1 为
t1时相下的综合植被指数，
η和
λ分别为在构建综合植被指数特征时NDVI和EVI指数的权重。
优选的，通过公式（12）构建步骤S52中的多源回归模型并回归同化点的产量与综合植被指数的关系：
(12)
其中，
y s 为指数-时序回归模型得到的产量，{
α t1 ⋯
α tn }为
t1时相至
tn时相下回归模型的系数，
β为回归模型的偏置。
优选的，步骤S1中多源遥感卫星数据的数据源为哨兵2号、Landsat8、高分一号或Planet卫星；步骤S12中时序遥感影像的反射率数据通过大气校正方法获取；步骤S12中的结构性参数为太阳天顶角、卫星天顶角和相对方位角。
本发明有益效果是：
1、本发明一方面考虑了作物的生长机理与遥感获取的多时序、多阶段作物的光谱信息，并将两者有机的耦合起来，再综合利用气象资料，所得到的预测模型更适合于对大豆作物产量的预测；另一方面，为了可以将模型在大区域上进行应用，本发明中将遥感影像的时序指数数据与同化产量数据进行多元回归模型构建，可快速实现大区域范围下的作物产量制图。
2、本发明通过计算遥感影像中的NDVI数据并引入混合高斯模型GMM对NDVI进行聚类，可自动化采集不同长势下大豆的发射率同化点并获取研究区内大豆作物不同生长状态的空间点位，以此避免人工采集时的采集不均衡现象，避免人为因素的影响，提高估算效率和估算精度。
3、本发明可支持国内外多种光学载荷的卫星遥感数据，相较于当前流行的深度学习评估方法，本专利在同化方法构建中具备丰富的理论依据，模型可解释性强。
附图说明
图1是本发明实施例提供的估产方法的流程图。
图2是本发明实施例提供的大豆同化点产量与区域产量的精度散点图。
具体实施方式
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图1-2及具体实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，而不构成对本发明的限制。
一种耦合机理模型的大豆作物遥感大区域快速估产方法，包括如下步骤：
S1：对研究区内的多源遥感卫星数据进行数据预处理；多源遥感卫星数据的数据源为哨兵2号、Landsat8、高分一号或Planet卫星。
S11：获取多源遥感卫星数据的时序遥感影像。
S12：获取时序遥感影像的反射率数据和结构性参数，构建多源卫星数据的元数据查询表；时序遥感影像的反射率数据通过大气校正方法获取，结构性参数为太阳天顶角、卫星天顶角和相对方位角。
S2：通过对遥感影像中的NDVI数据自动分级来获取研究区内大豆作物不同生长状态的空间点位；
S21：通过公式（1）计算遥感影像中的NDVI数据：
(1)
其中，
R nir 与
R r 分别为近红外波段与红波段的反射率值；采集的这些大豆采样点的红波段与近红外波段反射率用于作物模型同化的损失函数的构建。
S22：引入混合高斯模型GMM对NDVI进行聚类并自动化采集不同长势下大豆的发射率同化点；高斯混合模型GMM对NDVI进行聚类的过程如下：
S221：设置初始化高斯混合模型GMM的成分个数
K，并聚类不同长势下的NDVI；此处的
K为3，即大豆的长势分为长势好、长势中和长势差3个等级；
S222：计算大豆作物分布栅格影像中每个点的NDVI值属于各个高斯模型的概率，即通过公式（2）计算后验概率：
(2)
其中，
K为模型的个数，此处
K=3，
δ k 为属于第
k个高斯的概率，
δ k
＞0；为第
k个高斯的概率密度，
μ k 为均值向量, Σ k 为协方差矩阵；
S223：为每个NDVI点分别计算由该混合模型内每个高斯分量生成的概率，更新
δ k 、
μ k 、Σ k 参数，使得数据点属于长势好、长势中和长势差中某一类别的概率最大化；
S224：通过EM算法重复迭代步骤S222和步骤S223直到模型收敛；此时得到不同长势下NDVI的分级间隔点
θ 1 ，
θ 2 ；
S225：设NDVI1、NDVI2、NDVI3分别为长势差的点、长势中的点和长势好的点的NDVI值；根据需要同化点的数量，并根据阈值随机选取不同类别的同化采样点；NDVI1、NDVI2、NDVI3分别通过公式（3）、公式（4）和公式（5）进行判断：
。
S3：构建大豆作物生长模型的品质参数、土壤参数和管理参数。根据专家知识与参考文献从大豆作物生长模型所需要的土壤参数、管理参数及品质参数三个大类中，选取候选的一般参数，通过设置候选参数不同的范围步长来运行DSSAT模型，并使用EFAST来分析不同参数范围的设置对输出结果中LAI与产量这两个值的波动敏感性。最终确定的品质参数包括CSDL、PPSEN和EM-FL，土壤参数包括SLPF，管理参数包括种植时间、种植密度和施肥量。
S4：将大豆作物的生长模式与辐射传输模型进行耦合；
S41：引入DSSAT模型，利用粒子群算法根据步骤S3中的品质参数、土壤参数和管理参数的步长完成DSSAT模型的初始化；其中气象数据来源于NASA提供的气象格网数据产品，根据研究区的位置下载整年的气象数据。
S42：引入PROSAIL模型，根据初始化参数运行DSSAT模型，并得到大豆作物在整个生长周期内的LAI值和相应时间状态下的理论LAI数据序列，再将理论LAI数据序列及步骤S12中获取的太阳天顶角、卫星天顶角和相对方位角参数作为输入驱动PROSAIL模型，并得到该大豆同化点的光谱反射率曲线；
S43：计算当前大豆同化点红波段与近红外波段的反射率值，当前大豆同化点红波段与近红外波段的反射率值通过公式（6）和公式（7）得到：
(6)
(7)
其中，与分别为PROSAIL模型模拟的当前大豆同化点的近红外波段反射率值与红波段反射率值；与为遥感卫星数据源中近红外波段与红波段的光谱响应函数；ρo为PROSAIL模型模拟的当前大豆同化点的光谱曲线。
设
R nir 与
R r 分别为采样点遥感影像中的近红外波段与红波段的反射率值；将
R nir 、
R r 、和代入公式（8）中并求出误差损失
L：
(8)
对DSSAT模型的参数进行优化并得到DSSAT模型的最佳参数，此时DSSAT模型运行结果中的最终产量为当前大豆同化点的最优产量；对DSSAT模型的参数进行优化具体为：根据粒子群优化算法重复公式（6）-公式（8）直至损失收敛，此时得到的DSSAT模型参数为最佳参数。
S5：对步骤S2中得到的所有同化点，通过步骤S4进行依次循环并得到所有同化点的产量数据，并通过同化点的产量数据与时序遥感影像的指数数据构建多源回归模型，完成大区域范围内大豆作物产量的计算。
S51：计算遥感影像中的EVI数据、SIPI数据；并通过NDVI、EVI和SIPI数据计算该同化点的综合植被指数，以较好的代表作物当前的生长状态。
EVI数据和SIPI数据分别通过公式（9）和公式（10）得到：
(9)
(10)
其中，
R b 为采样点遥感影像中的蓝波段的反射率值。
综合植被指数通过公式（11）得到：
(11)
其中，
F t1 为
t1时相下的综合植被指数，
η和
λ分别为在构建综合植被指数特征时NDVI和EVI指数的权重。
S52：构建多源回归模型并回归同化点的产量与综合植被指数的关系，以保证大区域范围下产量的评估能力。
S53：构建研究区内时序遥感影像的NDVI、EVI与SIPI的指数集合，将步骤S52中的回归系数应用在研究区内的时序遥感影像中，完成大区域大豆作物的产量制图。
通过公式（12）构建步骤S52中的多源回归模型并回归同化点的产量与综合植被指数的关系：
(12)
其中，
y s 为指数-时序回归模型得到的产量，{
α t1 ⋯
α tn }为
t1时相至
tn时相下回归模型的系数，
β为回归模型的偏置。
如图2所示，本实施例中的方法对于大豆作物产量估算的精度较高，产量均方根误差 (RMSE)可达 216.3 kg/ha，在处理效率中，经计算为每小时处理644平方公里的影像数据。
同时，本实施例中所提的NDVI（Normalized Difference Vegetation Index）为归一化植被指数；EVI（Enhanced Vegetation Index）为增强植被指数；SIPI（StructureInsensitive Pigment Index）为结构不敏感色素指数；LAI（Leaf Area Index）为叶面积指数；CSDL（Critical ShortDay Length）为临界短日照时长；PPSEN（Slope of the relativeresponse of development to photoperiod with time）为发育对光周期的相对响应随时间的斜率；EM-FL（Timebetween plant emergence and flower appearance）为植物出苗和开花之间的时间；SLPF（Soil and fertilizer parameters）为土壤肥料参数；EM（Expectation-Maximization Algorithm）为期望最大化算法；DSSAT（DecisionSupportSystem For Agro-technology Transfer）为农业技术转移决策支持系统；PROSAIL：PROSPECT叶片光学模型和植被冠层二向反射率SAIL模型耦合而成，是目前应用最广泛的植被冠层辐射传输模型；EFAST（Extended Fourier Amplitude Sensitivity Test）为扩展傅里叶幅度敏感性检验；RMSE：均方根误差。
尽管上面已经示出和描述了本发明的实施例，可以理解的是，上述实施例是示例性的，不能理解为对本发明的限制。本领域的普通技术人员在本发明的范围内可以对上述实施例进行变化、修改、替换和变型。
以上本发明的具体实施方式，并不构成对本发明保护范围的限定。任何根据本发明的技术构思所做出的各种其他相应的改变与变形，均应包含在本发明权利要求的保护范围内。
Claims (10)
1.一种耦合机理模型的大豆作物遥感大区域快速估产方法，其特征在于，包括如下步骤：
S1：对研究区内的多源遥感卫星数据进行数据预处理；
S11：获取多源遥感卫星数据的时序遥感影像；
S12：获取时序遥感影像的反射率数据和结构性参数；
S2：通过对遥感影像中的NDVI数据自动分级来获取研究区内大豆作物不同生长状态的空间点位；
S21：计算遥感影像中的NDVI数据；
S22：引入混合高斯模型GMM对NDVI进行聚类并自动化采集不同长势下大豆的发射率同化点；
S3：构建大豆作物生长模型的品质参数、土壤参数和管理参数；
S4：将大豆作物的生长模式与辐射传输模型进行耦合；
S41：引入DSSAT模型，利用粒子群算法根据步骤S3中的品质参数、土壤参数和管理参数的步长完成DSSAT模型的初始化；
S42：引入PROSAIL模型，根据初始化参数运行DSSAT模型，并得到大豆作物在整个生长周期内的LAI值和相应时间状态下的理论LAI数据序列，再将理论LAI数据序列及步骤S12中获取的结构性参数作为输入驱动PROSAIL模型，并得到该大豆同化点的光谱反射率曲线；
S43：计算当前大豆同化点红波段与近红外波段的反射率值，对DSSAT模型的参数进行优化并得到DSSAT模型的最佳参数，此时DSSAT模型运行结果中的最终产量为当前大豆同化点的最优产量；
S5：对步骤S2中得到的所有同化点，通过步骤S4进行依次循环并得到所有同化点的产量数据，并通过同化点的产量数据与时序遥感影像的指数数据构建多源回归模型，完成大区域范围内大豆作物产量的计算；
S51：计算遥感影像中的EVI数据、SIPI数据；并通过NDVI、EVI和SIPI数据计算该同化点的综合植被指数；
S52：构建多源回归模型并回归同化点的产量与综合植被指数的关系；
S53：构建研究区内时序遥感影像的NDVI、EVI与SIPI的指数集合，将步骤S52中的回归系数应用在研究区内的时序遥感影像中，完成大区域大豆作物的产量制图。
2.根据权利要求1所述的耦合机理模型的大豆作物遥感大区域快速估产方法，其特征在于，步骤S21中通过公式（1）计算遥感影像中的NDVI数据：
(1)
其中，R nir 与R r 分别为近红外波段与红波段的反射率值。
3.根据权利要求2所述的耦合机理模型的大豆作物遥感大区域快速估产方法，其特征在于，步骤S22中引入高斯混合模型GMM对NDVI进行聚类的过程如下：
S221：设置初始化高斯混合模型GMM的成分个数K，并聚类不同长势下的NDVI；此处的K为3，即大豆的长势分为长势好、长势中和长势差3个等级；
S222：计算大豆作物分布栅格影像中每个点的NDVI值属于各个高斯模型的概率，即通过公式（2）计算后验概率：
(2)
其中，K为模型的个数，此处K=3，δ k 为属于第k个高斯的概率，δ k ＞0；为第k个高斯的概率密度，μ k 为均值向量, Σ k 为协方差矩阵；
S223：为每个NDVI点分别计算由该混合模型内每个高斯分量生成的概率，更新δ k 、μ k 、Σ k 参数，使得数据点属于长势好、长势中和长势差中某一类别的概率最大化；
S224：通过EM算法重复迭代步骤S222和步骤S223直到模型收敛；此时得到不同长势下NDVI的分级间隔点θ 1 ，θ 2 ；
S225：设NDVI1、NDVI2、NDVI3分别为长势差的点、长势中的点和长势好的点的NDVI值；根据需要同化点的数量，并根据阈值随机选取不同类别的同化采样点；NDVI1、NDVI2、NDVI3分别通过公式（3）、公式（4）和公式（5）进行判断：
。
4.根据权利要求3所述的耦合机理模型的大豆作物遥感大区域快速估产方法，其特征在于，步骤S3中的品质参数包括CSDL、PPSEN和EM-FL，土壤参数包括SLPF，管理参数包括种植时间、种植密度和施肥量。
5.根据权利要求4所述的耦合机理模型的大豆作物遥感大区域快速估产方法，其特征在于，步骤S43中当前大豆同化点红波段与近红外波段的反射率值通过公式（6）和公式（7）得到：
(6)
(7)
其中，与分别为PROSAIL模型模拟的当前大豆同化点的近红外波段反射率值与红波段反射率值；与为遥感卫星数据源中近红外波段与红波段的光谱响应函数；ρo为PROSAIL模型模拟的当前大豆同化点的光谱曲线。
6.根据权利要求5所述的耦合机理模型的大豆作物遥感大区域快速估产方法，其特征在于，设R nir 与R r 分别为采样点遥感影像中的近红外波段与红波段的反射率值；将R nir 、R r 、和代入公式（8）中并求出误差损失L：
(8)
步骤S43中对DSSAT模型的参数进行优化具体为：根据粒子群优化算法重复公式（6）-公式（8）直至损失收敛，此时得到的DSSAT模型参数为最佳参数。
7.根据权利要求6所述的耦合机理模型的大豆作物遥感大区域快速估产方法，其特征在于，步骤S51中的EVI数据和SIPI数据分别通过公式（9）和公式（10）得到：
(9)
(10)
其中，R b 为采样点遥感影像中的蓝波段的反射率值。
8.根据权利要求7所述的耦合机理模型的大豆作物遥感大区域快速估产方法，其特征在于，步骤S51中的综合植被指数通过公式（11）得到：
(11)
其中，F t1 为t1时相下的综合植被指数，η和λ分别为在构建综合植被指数特征时NDVI和EVI指数的权重。
9.根据权利要求8所述的耦合机理模型的大豆作物遥感大区域快速估产方法，其特征在于，通过公式（12）构建步骤S52中的多源回归模型并回归同化点的产量与综合植被指数的关系：
(12)
其中，y s 为指数-时序回归模型得到的产量，{α t1 ⋯α tn }为t1时相至tn时相下回归模型的系数，β为回归模型的偏置。
10.根据权利要求1-9中任一项所述的耦合机理模型的大豆作物遥感大区域快速估产方法，其特征在于，步骤S1中多源遥感卫星数据的数据源为哨兵2号、Landsat8、高分一号或Planet卫星；步骤S12中时序遥感影像的反射率数据通过大气校正方法获取；步骤S12中的结构性参数为太阳天顶角、卫星天顶角和相对方位角。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310213131.7A CN115879836B (zh) | 2023-03-08 | 2023-03-08 | 一种耦合机理模型的大豆作物遥感大区域快速估产方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310213131.7A CN115879836B (zh) | 2023-03-08 | 2023-03-08 | 一种耦合机理模型的大豆作物遥感大区域快速估产方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN115879836A CN115879836A (zh) | 2023-03-31 |
| CN115879836B true CN115879836B (zh) | 2023-05-12 |
Family
ID=85761949
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202310213131.7A Active CN115879836B (zh) | 2023-03-08 | 2023-03-08 | 一种耦合机理模型的大豆作物遥感大区域快速估产方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN115879836B (zh) |
Families Citing this family (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116472833B (zh) * | 2023-04-07 | 2024-06-18 | 中国长江三峡集团有限公司 | 一种作物灌水施肥方案的模拟方法、装置、电子设备 |
| CN117195609B (zh) * | 2023-11-08 | 2024-01-26 | 吉林大学 | 一种区域玉米植株氮营养匮缺的定量诊断方法 |
| CN117876870B (zh) * | 2024-01-12 | 2024-10-18 | 中国农业科学院农业资源与农业区划研究所 | 一种基于多源遥感数据的农作物估产的方法和系统 |
Citations (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102162850A (zh) * | 2010-04-12 | 2011-08-24 | 江苏省农业科学院 | 基于模型的小麦产量遥感监测预报方法 |
| CA3007202A1 (en) * | 2015-12-02 | 2017-06-08 | The Climate Corporation | Forecasting field level crop yield during a growing season |
| CN107122739A (zh) * | 2017-01-23 | 2017-09-01 | 东北农业大学 | 基于Extreme数学模型重构VI时间序列曲线的作物估产模型 |
| CN107941713A (zh) * | 2017-10-17 | 2018-04-20 | 河海大学 | 一种基于耦合作物模型同化光谱反射率的水稻估产方法 |
| CN109829234A (zh) * | 2019-01-30 | 2019-05-31 | 北京师范大学 | 一种基于高分辨率遥感数据和作物模型的跨尺度高精度动态作物长势监测和估产方法 |
| CN110032939A (zh) * | 2019-03-13 | 2019-07-19 | 浙江工业大学 | 一种基于高斯混合模型的遥感时序数据拟合方法 |
| CN115439402A (zh) * | 2022-07-05 | 2022-12-06 | 中国科学院空天信息创新研究院 | 基于SAFY-Sugar模型的糖料蔗单产遥感估算方法和装置 |
| CN115759524A (zh) * | 2022-10-20 | 2023-03-07 | 中国农业大学 | 一种基于遥感影像植被指数的土壤生产力等级识别方法 |
-
2023
- 2023-03-08 CN CN202310213131.7A patent/CN115879836B/zh active Active
Patent Citations (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102162850A (zh) * | 2010-04-12 | 2011-08-24 | 江苏省农业科学院 | 基于模型的小麦产量遥感监测预报方法 |
| CA3007202A1 (en) * | 2015-12-02 | 2017-06-08 | The Climate Corporation | Forecasting field level crop yield during a growing season |
| CN107122739A (zh) * | 2017-01-23 | 2017-09-01 | 东北农业大学 | 基于Extreme数学模型重构VI时间序列曲线的作物估产模型 |
| CN107941713A (zh) * | 2017-10-17 | 2018-04-20 | 河海大学 | 一种基于耦合作物模型同化光谱反射率的水稻估产方法 |
| CN109829234A (zh) * | 2019-01-30 | 2019-05-31 | 北京师范大学 | 一种基于高分辨率遥感数据和作物模型的跨尺度高精度动态作物长势监测和估产方法 |
| CN110032939A (zh) * | 2019-03-13 | 2019-07-19 | 浙江工业大学 | 一种基于高斯混合模型的遥感时序数据拟合方法 |
| CN115439402A (zh) * | 2022-07-05 | 2022-12-06 | 中国科学院空天信息创新研究院 | 基于SAFY-Sugar模型的糖料蔗单产遥感估算方法和装置 |
| CN115759524A (zh) * | 2022-10-20 | 2023-03-07 | 中国农业大学 | 一种基于遥感影像植被指数的土壤生产力等级识别方法 |
Non-Patent Citations (2)
| Title |
|---|
| 作物模型应用与遥感信息集成技术研究进展;彭慧文;中国农业气象;第43卷(第8期);全文 * |
| 基于遗传算法优化BP 神经网络的玉米遥感估产方法;于海洋;世界地质;第39卷(第1期);全文 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN115879836A (zh) | 2023-03-31 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN115879836B (zh) | 一种耦合机理模型的大豆作物遥感大区域快速估产方法 | |
| CN113268923B (zh) | 一种基于模拟多光谱的夏玉米产量估算方法 | |
| CN105445229B (zh) | 基于氮平衡光谱法的小麦春季氮肥施用方法及其追氮量模型的构建方法 | |
| CN109187441B (zh) | 基于冠层光谱信息的夏玉米含氮量监测模型的构建方法 | |
| Inoue et al. | Analysis of spectral measurements in paddy field for predicting rice growth and yield based on a simple crop simulation model | |
| CN107505271B (zh) | 基于氮素组分辐射传输模型的植株氮素估算方法和系统 | |
| CN110222475A (zh) | 一种基于无人机多光谱遥感反演冬小麦植株含水率的方法 | |
| CN111044516B (zh) | 一种水稻叶绿素含量遥感估测方法 | |
| CN111898922B (zh) | 一种多尺度农作物产量评估方法及系统 | |
| Liu et al. | Hierarchical optimization control based on crop growth model for greenhouse light environment | |
| CN112380497B (zh) | 用于区域玉米的秸秆系数估算方法及系统 | |
| CN110189793A (zh) | 基于高光谱的小麦氮肥生理利用率估测模型构建及不同氮效率小麦品种分类 | |
| CN112287886A (zh) | 基于高光谱图像融合图谱特征的小麦植株氮含量估测方法 | |
| CN118072178B (zh) | 基于分类后百分比数据同化的玉米产量估算方法及系统 | |
| Zou et al. | Combining spectral and texture feature of UAV image with plant height to improve LAI estimation of winter wheat at jointing stage | |
| CN107274002A (zh) | 小麦氮素丰缺诊断方法及其诊断模型的构建方法 | |
| Zhangzhong et al. | Development of an evapotranspiration estimation method for lettuce via mobile phones using machine vision: Proof of concept | |
| CN116452359A (zh) | 冬小麦长势的综合评价方法 | |
| Xie et al. | Estimating stomatal conductance of citrus under water stress based on multispectral imagery and machine learning methods | |
| Gao et al. | Regional winter wheat maturity date prediction using remote sensing-crop model data assimilation and numerical weather prediction | |
| CN117457066B (zh) | 一种省域尺度的冬小麦籽粒蛋白质含量预测方法 | |
| CN115344997B (zh) | 夏玉米植株叶片-冠层-像元尺度氮浓度协同预测方法 | |
| Zha | Design of rice regional test information collection system based on cloud computing. | |
| Wang et al. | Study on remote sensing feature selection of green manure crop Astragalus sinicus based on multitemporal Sentinel-2 imagery | |
| Sheng et al. | Agricultural Mechanization Services, Adverse Selection and Multi-Stage Productivity of Small Farms: Evidence from Wheat Production in North China |
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 |