一种基于无人机高光谱遥感的农作物叶绿素含量反演方法与流程

1.本发明涉及遥感信息的处理与应用技术领域，尤其是一种基于无人机高光谱遥感的农作物叶绿素含量反演方法。


背景技术：

2.传统农业中，农作物的生长参数监测主要通过人工破坏性取样结合手工测量来获取相关的农学指标。传统方法难以应用于大面积作物生长监测，且传统方法工作量大。然而，农作物生长及其进行光合作用和生理代谢的过程中，其长势包括一些生理生化参量，如叶片叶绿素含量、叶面积指数、植株含氮量等与冠层光谱特征密切相关。因此，可利用遥感手段实时获取农作物光谱信息，间接推演其生理生化参量等长势信息。卫星遥感虽然能进行大范围监测且成本较低，但容易受到云层的影响且运行周期长，缺乏精准农业所要求的时空分辨率。近年来，无人机在国内外被广泛应用于农业监测的各个领域，无人机操作灵活、数据获取速度较快、时空分辨率高，可应用于复杂的农田环境。


技术实现要素：

3.为了克服上述现有技术中的缺陷，本发明提供一种基于无人机高光谱遥感的农作物叶绿素含量反演方法，有利于大面积的农作物叶绿素含量监测，且不会对农作物产生破坏，为农作物叶片叶绿素含量的监测提供了一种新颖且可靠的方法，具有较高的实用价值。
4.为实现上述目的，本发明采用以下技术方案，包括：
5.一种基于无人机高光谱遥感的农作物叶绿素含量反演方法，包括以下步骤：
6.s1，对田地中的叶片叶绿素含量进行测量，田地划分为m个田块区，得到每个田块区的叶绿素含量测定值；
7.根据各个田块区的叶绿素含量测定值，构建标签集y，y＝{y
m
}，m＝1,2,
…
,m，y
m
为第m个田块区的叶绿素含量测定值；
8.s2，利用无人机获取田地的高光谱遥感影像，在n个波段上，得到每个田块区在各个波段上的反射率均值；
9.根据各个田块区在各个波段上的反射率均值，构建样本集x，根据各个田块区在各个波段上的反射率均值，构建样本集x，根据各个田块区在各个波段上的反射率均值，构建样本集x，为第m个田块区在第n个波段λ
n
上的反射率均值；
10.s3，对样本集x中田块区的向量进行聚类，得到j个聚类集；
11.每个聚类集中包括若干个田块区的向量；田块区在各个波段上的反射率均值所构成的向量即为该田块区的向量；
12.其中，第j个聚类集为即
13.j＝1,2,
…
,j；
14.k
j
为第j个聚类集c
j
中的田块区的向量的总数量；
15.为第j个聚类集c
j
中的第k个田块区的向量，即第j个聚类中心z
j
中的第k个田块
区在各个波段上的反射率均值所构成的向量，即即为第j个聚类集c
j
中第k个田块区在第n波段上的反射率均值；
16.s4，分别对各个聚类集进行特征波段的筛选，从各个聚类集中分别筛选出d个特征波段；根据各个聚类集中的各个田块区在各个特征波段上的反射率均值，分别构建各个聚类集的特征矩阵；
17.其中，第j个聚类集c
j
的特征矩阵为的特征矩阵为
18.为第j个聚类集c
j
中的第k个田块区在第d个特征波段上的反射率均值；
19.s5，将各个聚类集的特征矩阵分别与标签集y中对应的田块区的叶绿素含量测定值相对应，进行训练得到反演模型；
20.反演模型的输入为某田块区在各个波段上的反射率均值所构成的向量，即输入为某田块区的向量，反演模型的输出为该田块区的叶绿素含量的反演值。
21.步骤s2中，第m个田块区在第n个波段λ
n
上的反射率均值的计算方式如下所示：
[0022][0023]
其中，为第n个波段λ
n
上的第m个田块区内的第a个像素点的反射率值；n＝1,2,
…
,n；m＝1,2,
…
,m；a＝1,2,
…
,a
m
，a
m
为第m个田块区内的像素点总数量。
[0024]
步骤s3中，对样本集x中田块区的向量进行聚类，具体方式如下所示：
[0025]
s31，在样本集x中随机选取j个田块区的向量，j＜m；将此j个田块区的向量作为聚类算法中的聚类中心的初始值，构成聚类中心集z；将光谱角距离作为聚类算法中的聚类距离参量；
[0026]
s32，分别计算样本集x中各个田块区的向量与聚类中心集z中的各个聚类中心之间的光谱角距离，将样本集x中的各个田块区的向量分别聚类至与其光谱角距离最小的聚类中心上，得到各个聚类中心所聚类的田块区的向量，各个聚类中心所聚类的田块区的向量即构成该聚类中心所对应的聚类集，共得到j个聚类集；
[0027]
样本集x中的第m个田块区的向量x
m
与聚类中心集z中的第j个聚类中心z
j
之间的光谱角距离sad(x
m
,z
j
)，计算方式如下所示：
[0028][0029]
其中，sad(x
m
,z
j
)为样本集中的第m个田块区的向量x
m
与聚类中心集z中的第j个聚类中心z
j
之间的光谱角距离；
[0030]
x
m
为样本集x中的第m个田块区的向量，即为样本集x中的第m个田块区在各个波段上的反射率均值所构成的向量，即即为样本集x中的第m个田块区在第n个波段λ
n
上的反射率均值；m＝1,2,
…
,m；
[0031]
z
j
为聚类中心集z中的第j个聚类中心z
j
，即为聚类中心集z中的第j个田块区在各
个波段上的反射率均值所构成的向量即即为聚类中心集z中第j个田块区在第n个波段λ
n
上的反射率均值；j＝1,2,
…
,j；
[0032]
s33，分别计算各个聚类集的均向量；
[0033]
其中，第j个聚类集c
j
的均向量δ
j
的计算方式如下所示：
[0034][0035]
s34，分别根据各个聚类集的均向量对该聚类集的聚类中心进行更新，得到最优的聚类中心，以及得到最优的聚类中心所对应聚类集即最优的聚类集；
[0036]
其中，计算第j个聚类集c
j
的中心向量δ
j
与该第j个聚类中心z
j
的差值，若差值小于或等于设定的阈值t，则该第j个聚类中心z
j
不需要进行更新，该第j个聚类中心z
j
即为最优的聚类中心，该第j个聚类集c
j
即为最优的聚类集；若差值大于设定的阈值t，则对该第j个聚类中心z
j
进行更新，将该第j个聚类中心z
j
的值更新为第j个聚类集c
j
的中心向量δ
j
的值，重新进行聚类，即重新执行步骤s32～s34，直至得到最优的第j个聚类中心z
j
，该最优的第j个聚类中心z
j
所对应聚类集c
j
即为最优的聚类集。
[0037]
步骤s4中，对第j个聚类集进行特征波段筛选，具体方式如下所示：
[0038]
s41，初始化d＝1，在第j个聚类集中，随机选择一个波段λ
ω
，将波段λ
ω
的向量作为第d次投影的起始向量h
d
；第j个聚类集中的波段λ
ω
上的各个田块区的反射率均值所构成的向量即为该波段λ
ω
的向量的向量即即为第j个聚类集c
j
中第k个田块区在波段λ
ω
上的反射率均值；
[0039]
s42，对第j个聚类集进行第d次投影，分别计算第d次投影的起始向量h
d
对对第j个聚类集中各个波段的向量的投影向量；
[0040]
其中，第d次投影的起始向量h
d
对第j个聚类集中第n个波段λ
n
的向量的投影向量计算方式如下所示：
[0041][0042]
其中，n＝1,2,
…
,n且,n且表示的的转置；h
dt
表示的h
d
的转置；
[0043]
为第d次投影的起始向量h
d
对第j个聚类集c
j
中第n个波段λ
n
的向量的投影向量；
[0044]
第j个聚类集c
j
中第n个波段λ
n
的向量为：第j个聚类集c
j
中第n个波段λ
n
上的各个田块区反射率均值所构成的向量，即即为第j个聚类集c
j
中第k个田块区在第n波段上的反射率均值；n＝1,2,
…
,n；
[0045]
s43，记录最大的投影向量所对应的波段，该波段即为第j个聚类集的第d次投影的特征波段
[0046]
s44，对第j个聚类集进行下一次投影即第d 1次投影，将第d次投影的特征波段的反射率均值向量x
d
的值作为第d 1次投影的起始向量h
d 1
的值，按照步骤s42～s43的方式
进行第d 1次投影，得到第d 1次投影的特征波段
[0047]
s45，若d 1＝d，则不再进行下一次投影，即完成第j个聚类集的d个特征波段的筛选；若d 1＜d，则按照步骤s44的方式对第j个聚类集继续进行下一次投影；
[0048]
s46，从第j个聚类集c
j
中筛选出d个特征波段后，根据第j个聚类集c
j
中的各个田块区在各个特征波段上的反射率均值，构建第j个聚类集c
j
的特征矩阵
[0049]
步骤s5中，反演模型的训练方式如下所示：
[0050]
s51，各个聚类集的特征矩阵均相当于一个决策树，聚类集的特征波段即为对应的决策树的特征波段，将各个聚类集的特征矩阵分别与标签集y相对应；
[0051]
s52，分别基于各个聚类集的特征矩阵和标签集y中对应的田块区的叶绿素含量测定值，训练生成各个决策树，将多个决策树集成起来，建立随机森林回归模型即反演模型。
[0052]
利用反演模型进行预测时，输入测试样本，测试样本为某个田块区在各个波段上的反射率均值所构成的向量，随机森林回归模型中的每个决策树分别提取测试样本中对应的特征波段上的反射率均值，每个决策树根据对应的特征波段上的反射率均值进行预测，每个决策树均输出该田块区的叶绿素含量预测值，对各个决策树输出的该田块区的叶绿素含量预测值计算平均值，并将该平均值作为最终的该田块区的叶绿素含量的反演值。
[0053]
n＝176。
[0054]
本发明的优点在于：
[0055]
(1)传统的农作物叶绿素含量监测是通过破坏性取样，且只能实现单取样点的监测，难以应用于大面积的农作物监测。本发明将无人机高光谱遥感技术引入到农作物叶绿素含量监测中，实现了由点到面即由单取样点到田块区的突破，有利于大面积的农作物叶绿素含量监测，且不会对农作物产生破坏。
[0056]
(2)本发明利用光谱角距离对田块区进行光谱聚类，可以将光谱中所含有的特征以及光谱之间存在的差异之处，通过聚类后最大限度的表现出来，区别于以往的整体建模。本发明先聚类后再建模即建立反演模型，基于差异化建模的思想，可以使得随机森林回归模型即反演模型的精度更高。
[0057]
(3)本发明进行特征波段的筛选，降低建模维度，降低反演模型的冗余性和复杂度，避免过拟合现象。
附图说明
[0058]
图1为本发明的流程图。
[0059]
图2为获取的试验田地的无人机高光谱遥感影像。
[0060]
图3为试验田地中部分田块区的叶绿素含量的的反演结果示意图。
[0061]
图4为反演模型的验证精度图。
具体实施方式
[0062]
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他
实施例，都属于本发明保护的范围。
[0063]
由图1所示，本发明的一种基于无人机高光谱遥感的农作物叶绿素含量反演方法，包括以下步骤：
[0064]
s1，对试验田地中的叶片叶绿素含量进行测量，田地划分为m个田块区，得到每个田块区的叶绿素含量测定值。根据各个田块区的叶绿素含量测定值，构建标签集y，y＝{y
m
}，m＝1,2,
…
,m，y
m
为第m个田块区的叶绿素含量测定值。
[0065]
步骤s1中，对试验田地进行划分，划分为m个田块区，本实施例中，m＝126；使用日产spad
‑
502型手持式叶绿素仪对每个田块区均进行若干个采样点的叶片叶绿素含量测定，每个田块区筛选出一个采样点的叶片叶绿素含量作为该田块区的叶绿素含量测定值。
[0066]
s2，利用无人机获取试验田地的高光谱遥感影像，在n个波段上，得到每个田块区在各个波段上的反射率均值。根据各个田块区在各个波段上的反射率均值，构建样本集x，在各个波段上的反射率均值。根据各个田块区在各个波段上的反射率均值，构建样本集x，为第m个田块区在第n个波段λ
n
上的反射率均值。
[0067]
步骤s2中，使用大疆m600pro六悬翼无人机搭载gaiasky
‑
mini2高光谱成像仪，获取蒙城马店试验站长期定位试验田地(116
°
37
′
e，33
°
13
′
n)的拔节期小麦的高光谱遥感影像。无人机获取该试验田地的多幅高光谱遥感影像，使用hispectralstitcher软件对该试验田地的多幅高光谱遥感影像进行拼接，得到试验田地完整的高光谱遥感影像。使用specview软件对试验田地的高光谱遥感影像进行预处理，所述预处理包括：辐射定标和大气校正，预处理后得到高光谱遥感影像中各个像素点的反射率值。
[0068]
步骤s2中，第m个田块区在第n个波段λ
n
上的反射率均值的计算方式如下所示：
[0069][0070]
其中，为第n个波段λ
n
上的第m个田块区内的第a个像素点的反射率值；n＝1,2,
…
,n；m＝1,2,
…
,m；a＝1,2,
…
,a
m
，a
m
为第m个田块区内的像素点总数量。
[0071]
步骤s2中，还可使用envi软件中的roi工具从每个田块区中提取一个感兴趣区，在n个波段上，分别导出各个波段上的每个田块区的感兴趣区内的各个像素点的反射率值，计算每个田块区的感兴趣区在各个波段上的反射率均值，将每个田块区的感兴趣区在各个波段上的反射率均值作为该田块区的在各个波段上的反射率均值。即，第m个田块区在第n个波段λ
n
上的反射率均值为：
[0072][0073]
其中，为第n个波段λ
n
上的第m个田块区的感兴趣区内的第a
′
个像素点的反射率值；n＝1,2,
…
,n；m＝1,2,
…
,m；a
′
＝1,2,
…
,a
′
m
，a
′
m
为第m个田块区的感兴趣区内的像素点总数量。
[0074]
s3，对样本集x中田块区的向量进行聚类，得到j个聚类集。每个聚类集中包括若干个田块区的向量；田块区在各个波段上的反射率均值所构成的向量即为该田块区的向量。
[0075]
步骤s3中，对样本集x中田块区的向量进行聚类，具体方式如下所示：
[0076]
s31，在样本集x中随机选取j个田块区的向量，且j＜m，将此j个田块区的向量作为聚类算法中的聚类中心的初始值，构成聚类中心集z；将光谱角距离作为聚类算法中的聚类距离参量。
[0077]
s32，分别计算样本集x中各个田块区的向量与聚类中心集z中的各个聚类中心之间的光谱角距离，将样本集x中的各个田块区的向量分别聚类至光谱角距离最小的聚类中心上，得到各个聚类中心所聚类的田块区的向量，各个聚类中心所聚类的田块区的向量即构成该聚类中心所对应的聚类集。
[0078]
其中，样本集x中的第m个田块区的向量x
m
与聚类中心集z中的第j个聚类中心z
j
之间的光谱角距离sad(x
m
,z
j
)，计算方式如下所示：
[0079][0080]
其中，sad(x
m
,z
j
)为样本集中的第m个田块区的向量x
m
与聚类中心集z中的第j个聚类中心z
j
之间的光谱角距离。
[0081]
x
m
为样本集x中的第m个田块区的向量，即为样本集x中的第m个田块区在各个波段上的反射率均值所构成的向量，即即为样本集x中的第m个田块区在第n个波段λ
n
上的反射率均值；m＝1,2,
…
,m。
[0082]
z
j
为聚类中心集z中的第j个聚类中心z
j
，即为聚类中心集z中的第j个田块区在各个波段上的反射率均值所构成的向量，即即为聚类中心集z中第j个田块区在第n个波段λ
n
上的反射率均值；j＝1,2,
…
,j。
[0083]
聚类中心集z中的第j个聚类中心z
j
所对应的聚类集，即为第j个聚类集c
j
，即即为第j个聚类中心z
j
所聚类的第k个田块区的向量，即第j个聚类集c
j
中的第k个田块区的向量，即第j个聚类中心z
j
所聚类的第k个田块区在各个波段上的反射率均值所构成的向量，即k
j
为第j个聚类集c
j
中的田块区的向量的总数量，为第j个聚类集c
j
中第k个田块区在第n波段上的反射率均值。
[0084]
当某田块区的向量与聚类中心相等时，此时二者之间的光谱角距离为最小，因此聚类集中包括与聚类中心相等的某田块区的向量，即相当于包括聚类中心本身。
[0085]
s33，分别计算各个聚类集的均向量。
[0086]
其中，第j个聚类集c
j
的均向量δ
j
的计算方式如下所示：
[0087][0088]
s34，分别根据各个聚类集的均向量对该聚类集的聚类中心进行更新，得到各个最优聚类中心，以及各个最优聚类中心所对应的聚类集即最优聚类集。
[0089]
其中，计算第j个聚类集c
j
的中心向量δ
j
与该第j个聚类中心z
j
的差值，若差值小于或等于设定的阈值t，则该第j个聚类中心z
j
不需要进行更新，该第j个聚类中心z
j
即为最优
聚类中心，该第j个聚类集c
j
即为最优聚类集；若差值大于设定的阈值t，则对该第j个聚类中心z
j
进行更新，将该第j个聚类中心z
j
的值更新为第j个聚类集c
j
的中心向量δ
j
的值，重新进行聚类，即重新执行步骤s32～s34，直至得到最优的第j个聚类中心z
j
，该最优的第j个聚类中心z
j
所对应聚类集c
j
即为最优聚类集。
[0090]
s4，分别对各个聚类集进行特征波段的筛选，从各个聚类集中分别筛选出d个特征波段；根据各个聚类集中的各个田块区在各个特征波段上的反射率均值，分别构建各个聚类集的特征矩阵。
[0091]
步骤s4中，对步骤s3得到的最优的第j个聚类集c
j
进行特征波段筛选，具体方式如下所示：
[0092]
s41初始化d＝1，在第j个聚类集中，随机选择一个波段λ
ω
，将波段λ
ω
的向量作为第d次投影的起始向量h
d
；第j个聚类集中的波段λ
ω
上的各个田块区的反射率均值所构成的向量即为该波段λ
ω
的向量的向量即即为第j个聚类集c
j
中第k个田块区在波段λ
ω
上的反射率均值。
[0093]
s42，对第j个聚类集c
j
进行第d次投影，分别计算第d次投影的起始向量h
d
对第j个聚类集c
j
中各个波段的向量的投影向量。
[0094]
第d次投影的起始向量h
d
对第j个聚类集c
j
中第n个波段λ
n
的向量的投影向量计算方式如下所示：
[0095][0096]
其中，n＝1,2,
…
,n且,n且表示的的转置；h
dt
表示的h
d
的转置。
[0097]
为第d次投影的起始向量h
d
对第j个聚类集c
j
中第n个波段λ
n
的向量的投影向量。
[0098]
第j个聚类集c
j
中第n个波段λ
n
的向量为：第j个聚类集c
j
中第n个波段λ
n
上的各个田块区反射率均值所构成的向量，即即为第j个聚类集c
j
中第k个田块区在第n波段上的反射率均值；n＝1,2,
…
,n。
[0099]
s43，记录最大的投影向量所对应的波段，该波段即为第j个聚类集c
j
的第d次投影的特征波段
[0100]
s44，对第j个聚类集c
j
进行下一次投影即第d 1次投影，将第d次投影的特征波段的反射率均值向量x
d
的值作为第d 1次投影的起始向量h
d 1
的值，按照步骤s42～s43的方式进行第d 1次投影，得到第d 1次投影的特征波段
[0101]
s45，若d 1＝d，则不再进行下一次投影，即完成第j个聚类集c
j
的d个特征波段的筛选；若d 1＜d，则按照步骤s44的方式对第j个聚类集c
j
继续进行下一次投影。
[0102]
s46，从第j个聚类集c
j
中筛选出d个特征波段后，根据第j个聚类集c
j
中的各个田块区在各个特征波段上的反射率均值，构建第j个聚类集c
j
的特征矩阵
为第j个聚类集c
j
中的第k个田块区在第d个特征波段上的反射率均值。
[0103]
s5，将各个聚类集的特征矩阵分别与标签集y中对应的田块区的叶绿素含量测定值相对应，进行训练得到反演模型；
[0104]
反演模型的输入为某田块区在各个波段上的反射率均值所构成的向量，即输入为该田块区的向量，反演模型的输出为该田块区的叶绿素含量的反演值。
[0105]
步骤s5中，反演模型的训练方式如下所示：
[0106]
s51，各个聚类集的特征矩阵均相当于一个决策树，聚类集的特征波段即为对应的决策树的特征波段，将各个聚类集的特征矩阵分别与标签集y相对应；
[0107]
s52，分别基于各个最优聚类集的特征矩阵和标签集y中对应的田块区的叶绿素含量测定值，训练得到各个决策树，将多个决策树集成起来，建立随机森林回归模型即反演模型。
[0108]
利用反演模型进行预测时，输入测试样本，测试样本为某个田块区在各个波段上的反射率均值所构成的向量，随机森林回归模型中的每个决策树分别根据测试样本中对应的特征波段上的反射率均值进行预测，每个决策树均输出该田块区的叶绿素含量预测值，对各个决策树输出的该田块区的叶绿素含量预测值计算平均值，并将该平均值作为最终的该田块区的叶绿素含量的反演值。
[0109]
本实施例中，试验区为农业农村部蒙城砂姜黑土生态环境站内(116
°
37
′
e，33
°
13
′
n)，图2为获取的试验田地的无人机高光谱遥感影像，白色虚线框内为本次试验田地。图3为试验田地中部分田块区的叶绿素含量的模型预测结果，结果显示，整体上，拔节期冬小麦叶绿素含量值的范围为36～53，模型可以实现对冬小麦叶片叶绿素含量也就是相对叶绿素含量的由点到面的精确预测。图4为反演模型的验证精度图，由图4所示，模型预测的spad值与实测spad值之间存在较强的正相关性，决定系数r2即拟合优度达到0.77。
[0110]
本发明以无人机高光谱遥感影响为基础，充分利用小麦田地的丰富的高光谱信息，结合各田块区的小麦叶片叶绿素含量的测定值，采用pre
‑
clustered rf分步回归的高级机器学习算法，建立反演模型，对拔节期冬小麦叶片叶绿素含量进行反演，获得较好的反演结果，能直观、准确地展示长期定位试验田地冬小麦叶绿素含量的空间分布及其特征，为小麦叶片叶绿素信息监测及土壤肥力对叶绿素含量的影响判断提供了一种新颖且可靠的方法，具有较高的实用价值。
[0111]
以上仅为本发明创造的较佳实施例而已，并不用以限制本发明创造，凡在本发明创造的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明创造的保护范围之内。