CN115063610B - 基于Sentinel-1、2影像的大豆种植区识别方法 - Google Patents
基于Sentinel-1、2影像的大豆种植区识别方法 Download PDFInfo
- Publication number
- CN115063610B CN115063610B CN202210596186.6A CN202210596186A CN115063610B CN 115063610 B CN115063610 B CN 115063610B CN 202210596186 A CN202210596186 A CN 202210596186A CN 115063610 B CN115063610 B CN 115063610B
- Authority
- CN
- China
- Prior art keywords
- sentinel
- feature
- soybean
- features
- time sequence
- 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
- 244000068988 Glycine max Species 0.000 title claims abstract description 102
- 235000010469 Glycine max Nutrition 0.000 title claims abstract description 100
- 238000000034 method Methods 0.000 title claims abstract description 48
- 238000007637 random forest analysis Methods 0.000 claims abstract description 26
- 238000012706 support-vector machine Methods 0.000 claims abstract description 24
- 238000007781 pre-processing Methods 0.000 claims abstract description 12
- 230000006870 function Effects 0.000 claims description 51
- 238000004364 calculation method Methods 0.000 claims description 33
- 238000003066 decision tree Methods 0.000 claims description 24
- 238000012549 training Methods 0.000 claims description 20
- 238000002310 reflectometry Methods 0.000 claims description 14
- 238000000605 extraction Methods 0.000 claims description 13
- 238000013145 classification model Methods 0.000 claims description 11
- 230000008569 process Effects 0.000 claims description 9
- 238000013507 mapping Methods 0.000 claims description 8
- 238000003971 tillage Methods 0.000 claims description 6
- 238000004422 calculation algorithm Methods 0.000 claims description 5
- 238000005192 partition Methods 0.000 claims description 5
- 238000012545 processing Methods 0.000 claims description 4
- 230000001174 ascending effect Effects 0.000 claims description 3
- 230000008859 change Effects 0.000 claims description 3
- 239000003638 chemical reducing agent Substances 0.000 claims description 3
- 238000010276 construction Methods 0.000 claims description 3
- 238000006073 displacement reaction Methods 0.000 claims description 3
- 230000009977 dual effect Effects 0.000 claims description 3
- 238000003384 imaging method Methods 0.000 claims description 3
- 238000005305 interferometry Methods 0.000 claims description 3
- 230000000873 masking effect Effects 0.000 claims description 3
- 238000003908 quality control method Methods 0.000 claims description 3
- 238000005070 sampling Methods 0.000 claims description 3
- 230000002123 temporal effect Effects 0.000 claims description 3
- 238000001228 spectrum Methods 0.000 abstract description 4
- 230000007547 defect Effects 0.000 abstract description 2
- 238000000691 measurement method Methods 0.000 abstract 1
- 230000000875 corresponding effect Effects 0.000 description 6
- 238000010586 diagram Methods 0.000 description 5
- 238000011160 research Methods 0.000 description 5
- 238000012216 screening Methods 0.000 description 4
- 230000003595 spectral effect Effects 0.000 description 4
- 241000196324 Embryophyta Species 0.000 description 3
- 238000005516 engineering process Methods 0.000 description 3
- 238000010801 machine learning Methods 0.000 description 3
- 238000004088 simulation Methods 0.000 description 3
- 238000004458 analytical method Methods 0.000 description 2
- 230000015572 biosynthetic process Effects 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 238000011835 investigation Methods 0.000 description 2
- 238000012821 model calculation Methods 0.000 description 2
- 230000009467 reduction Effects 0.000 description 2
- 238000003786 synthesis reaction Methods 0.000 description 2
- 240000008042 Zea mays Species 0.000 description 1
- 235000005824 Zea mays ssp. parviglumis Nutrition 0.000 description 1
- 235000002017 Zea mays subsp mays Nutrition 0.000 description 1
- 238000009825 accumulation Methods 0.000 description 1
- 238000012512 characterization method Methods 0.000 description 1
- 238000007635 classification algorithm Methods 0.000 description 1
- 235000005822 corn Nutrition 0.000 description 1
- 230000002596 correlated effect Effects 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 238000012417 linear regression Methods 0.000 description 1
- 239000011159 matrix material Substances 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000012544 monitoring process Methods 0.000 description 1
- 238000005457 optimization Methods 0.000 description 1
- WPPDXAHGCGPUPK-UHFFFAOYSA-N red 2 Chemical compound C1=CC=CC=C1C(C1=CC=CC=C11)=C(C=2C=3C4=CC=C5C6=CC=C7C8=C(C=9C=CC=CC=9)C9=CC=CC=C9C(C=9C=CC=CC=9)=C8C8=CC=C(C6=C87)C(C=35)=CC=2)C4=C1C1=CC=CC=C1 WPPDXAHGCGPUPK-UHFFFAOYSA-N 0.000 description 1
- 238000007619 statistical method Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
- G06V10/00—Arrangements for image or video recognition or understanding
- G06V10/40—Extraction of image or video features
- G06V10/62—Extraction of image or video features relating to a temporal dimension, e.g. time-based feature extraction; Pattern tracking
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
- G06N20/00—Machine learning
- G06N20/10—Machine learning using kernel methods, e.g. support vector machines [SVM]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
- G06N20/00—Machine learning
- G06N20/20—Ensemble learning
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/10—Segmentation; Edge detection
- G06T7/11—Region-based segmentation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/60—Analysis of geometric attributes
- G06T7/62—Analysis of geometric attributes of area, perimeter, diameter or volume
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
- G06V10/00—Arrangements for image or video recognition or understanding
- G06V10/70—Arrangements for image or video recognition or understanding using pattern recognition or machine learning
- G06V10/764—Arrangements for image or video recognition or understanding using pattern recognition or machine learning using classification, e.g. of video objects
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
- G06V10/00—Arrangements for image or video recognition or understanding
- G06V10/70—Arrangements for image or video recognition or understanding using pattern recognition or machine learning
- G06V10/77—Processing image or video features in feature spaces; using data integration or data reduction, e.g. principal component analysis [PCA] or independent component analysis [ICA] or self-organising maps [SOM]; Blind source separation
- G06V10/774—Generating sets of training patterns; Bootstrap methods, e.g. bagging or boosting
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/10—Image acquisition modality
- G06T2207/10024—Color image
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/10—Image acquisition modality
- G06T2207/10032—Satellite or aerial image; Remote sensing
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/10—Image acquisition modality
- G06T2207/10048—Infrared image
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/20—Special algorithmic details
- G06T2207/20112—Image segmentation details
- G06T2207/20132—Image cropping
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/30—Subject of image; Context of image processing
- G06T2207/30181—Earth observation
- G06T2207/30188—Vegetation; Agriculture
Landscapes
- Engineering & Computer Science (AREA)
- Theoretical Computer Science (AREA)
- Physics & Mathematics (AREA)
- Computer Vision & Pattern Recognition (AREA)
- General Physics & Mathematics (AREA)
- Software Systems (AREA)
- Evolutionary Computation (AREA)
- Medical Informatics (AREA)
- Artificial Intelligence (AREA)
- Computing Systems (AREA)
- Multimedia (AREA)
- Databases & Information Systems (AREA)
- General Health & Medical Sciences (AREA)
- Health & Medical Sciences (AREA)
- Data Mining & Analysis (AREA)
- General Engineering & Computer Science (AREA)
- Mathematical Physics (AREA)
- Geometry (AREA)
- Image Processing (AREA)
Abstract
本发明涉及基于Sentinel‑1、2影像的大豆种植区识别方法及其面积测算方法，与现有技术相比解决了大豆与其他作物光谱相似度高导致其依靠高维特征难以实现种植区识别的缺陷。本发明包括以下步骤：Sentinel‑1、2影像的获取和预处理；时间序列特征提取；支持向量机模型的构建；优选特征子集确定；大豆种植区识别。本发明借助GEE云计算平台，利用线性谐波模型提取大豆生长季内Sentinel‑1、2影像的时间序列特征，然后构造支持向量机模型，同时结合随机森林分类模型及斯皮尔曼相关系数探究大豆识别优选特征子集，最终利用支持向量机模型提取大豆种植区并测算面积。
Description
技术领域
本发明涉及遥感数据处理技术领域，具体来说是基于Sentinel-1、2影像的大豆种植区识别方法。
背景技术
传统方法获取大豆种植面积依赖于田间调查与统计上报，这种调查方式高成本低效率，且不具备详细刻画种植区域空间分布的能力。遥感技术具有覆盖范围广、采集数据快、动态对地观测等优点，能够通过低成本高效率的方式实现农作物种植面积提取与监测。
目前基于时间序列数据尤其是植被指数的时序数据已经成为农作物遥感分类识别研究的新机遇。植被指数时间序列数据能够精确地反映作物物候信息(如大豆的出现、结荚、灌浆、成熟)，有效削弱“同物异谱，同谱异物”现象，在农作物遥感分类识别与信息提取研究中发挥了重要作用。增强植被指数和归一化差异植被指数是目前利用遥感技术进行农作物信息识别与提取最常用的植被指数指标。线性谐波模型被广泛的应用于时间序列数据特征提取。特征筛选能够降低冗余信息，缩减特征数量，减少计算运行时间。因此，研究基于优选的时间序列特征提取大豆种植区的技术迫在眉睫。
现有技术中，一种大豆生长季空间分布图的生成方法和系统，其构造大豆生长季内Sentinel-2光谱波段的时间序列合成影像，然后结合随机森林分类模型探究大豆最早识别的时间窗口，其次通过评估时间窗口内所有特征的重要性进一步筛选特征子集，最终绘制出大豆空间分布图。另有，基于植被指数时序谱特征的大豆种植区提取方法和系统，通过构建中高时空分辨率的时间序列植被指数，利用线性谐波模型提取研究区植被指数时序谱特征。其次，基于关键生育期遥感原始影像波段反射率，通过百分位数、最大值、最小值、平均值、标准差方法提取光谱特征。
然而，在种植结构复杂或者田块破碎的地区，仅仅利用光谱特征进行大豆识别时，很难将与大豆光谱相似度高的作物如玉米进行区分。另外，机器学习在实际应用时，输入特征个数越多，不仅会增加模型计算开销，同时还可能引起“维度灾难”，即随着高维特征的引入，分类效率和分类精度均会受到不同程度的影响。最后，仅仅利用单一机器学习分类模型进行特征分析及分类时，结果容易受到模型本身“偏好”而产生一定的误差。
因此，如何设计出一种高精确度的大豆种植区快速准确识别方法已经成为急需解决的技术问题。
发明内容
本发明的目的是为了解决现有技术中大豆与其他作物光谱相似度高导致其依靠高维特征难以实现种植区识别的缺陷，提供一种基于Sentinel-1、2影像的大豆种植区识别方法来解决上述问题。
为了实现上述目的，本发明的技术方案如下：
一种基于Sentinel-1、2影像的大豆种植区识别方法，包括以下步骤：
Sentinel-1、2影像的获取和预处理：获取大豆生长季内Sentinel-2卫星影像，进行预处理后计算5个植被指数并选择5个原始波段；获取大豆生长季内Sentinel-1影像，并进行预处理后选择2个原始波段；分别对所选的植被指数及原始波段进行耕地掩膜处理；
时间序列特征提取：利用2阶线性谐波模型提取大豆生长季内掩膜后的Sentinel-2卫星影像的5个取植被指数及5个原始波段的时间序列特征，作为Sentinel-2候选时间序列特征；利用2阶线性谐波模型提取大豆生长季内掩膜后的Sentinel-1卫星影像的2个原始波段的时间序列特征，作为Sentinel-1候选时间序列特征；汇总Sentinel-2和Sentinel-1所有候选时间序列特征，得到时间序列特征；
支持向量机模型的构建：构建支持向量机分类模型；
优选特征子集确定：构建随机森林分类模型，利用随机森林分类模型计算所有时间序列特征的重要性得分；分别计算Sentinel-2和Sentinel-1候选时间序列特征中每两个候选特征之间的斯皮尔曼相关系数，按照时间序列特征的重要性降序去除斯皮尔曼相关系数大于0.5或小于-0.5的候选特征，得到优选特征子集；
大豆种植区识别：将优选特征子集输入支持向量机分类模型训练，得到分类结果并提取出大豆种植区。
所述Sentinel-1、2影像的获取和预处理包括以下步骤：
在GlobeLand30网站下载30m空间分辨率的土地覆盖数据，在ArcGIS软件中进行镶嵌及剪裁，然后利用ArcGIS软件中提取工具获得耕地数据，并上传至GEE云计算平台；
在GEE云计算平台中调用大豆生长季内所有的Sentinel-2卫星大气底端反射率数据，利用质量控制QA60波段将云量限制在10％以内，并将所有波段值除以10000得到各个波段反射率值；
计算大豆生长季内Sentinel-2卫星影像的增强型植被指数EVI、归一化植被指数NDVI、归一化耕作指数NDTI、归一化差值衰变指数NDSVI、红边位置指数REPI，共5个植被指数；
并选择大豆生长季内Sentinel-2卫星影像的红边1、红边2、红边3、短波红外1和短波红外2，共5个原始波段；
利用耕地数据对大豆生长季内Sentinel-2卫星影像的5个植被指数和5个原始波段进行耕地掩膜，
其中Sentinel-2卫星影像的植被指数计算公式如下所示：
其中，Red、Blue、NIR分别是Sentinel-2的红光、蓝光和近红外波段反射率值；
其中，NIR和Red分别是Sentinel-2的近红外和红光波段反射率值；
其中，SWIR1和SWIR2分别是Sentinel-2的短波红外1和短波红外2波段反射率值；
其中，SWIR1和Red分别是Sentinel-2的短波红外1和红光波段反射率值；
其中，Red是Sentinel-2的红光波段反射率值，RE1、RE2和RE3是Sentinel-2的三个红边波段即红边1、红边2、红边3波段反射率值；
在GEE云计算平台中调用大豆生长季内所有的Sentinel-1卫星Level-1地距影像数据，成像方式为干涉测量宽幅模式，利用Sentinel-1工具箱进行预处理；选择Sentinel-1卫星影像的VV、VH波段，共2个原始波段；利用耕地数据对大豆生长季内Sentinel-1卫星影像的2个原始波段进行掩膜处理。
所述时间序列特征提取包括以下步骤：
利用2阶线性谐波模型提取大豆生长季内掩膜后的Sentinel-2卫星影像的5个植被指数和5个原始波段的时间序列特征；利用2阶线性谐波模型提取大豆生长季内掩膜后的Sentinel-1卫星影像的2个原始波段的时间序列特征；其中，2阶线性谐波模型计算公式如下：
f(t)＝b+C1cos(3πt)+D1sin(3πt)+C2cos(6πt)+D2sin(6πt)，
其中，f(t)是第t时刻分别拟合的Sentinel-2卫星影像的5个植被指数和5个原始波段值以及Sentinel-1卫星影像的2个原始波段值，b是常数项，C1、C2、D1、D2是余弦函数和正弦函数的系数，t是一年中的一天，
取相位和振幅作为各个时间序列特征输入，其中振幅1定义为二维矢量[C1，D1]的长度，振幅2定义为二维矢量[C2，D2]的长度，相位1定义为二维矢量[C1，D1]的角度，相位2定义为二维矢量[C2，D2]的角度；
汇总Sentinel-2和Sentinel-1所有候选时间序列特征，得到时间序列特征。其中汇总的计算公式如下所示：
VI＝VI1.addBands(VI2).addBands(VI3)…addBands(VIx)
其中，VI1表示第一个候选时间序列特征，VIx表示为最后一个候选时间特征，.addBands()为添加特征的函数。
所述支持向量机模型的构建包括以下步骤：
在GEE云计算平台中调用支持向量机分类模型，选择核函数类型及其参数；
设定核函数类型选择线性核函数，cost设置为1，其他参数保持默认值。
所述优选特征子集确定包括以下步骤：
在GEE云计算平台中调用随机森林分类模型，设置的参数包括：决策树的数量、最大特征数、决策树深度、叶子节点最少样本数、节点划分的最小样本数，其中将决策树的数量设置为300，其他参数保持默认值；
将时间序列特征输入随机森林分类模型进行训练，其步骤如下：
基于时间序列特征集，采用随机且有放回地抽样Bagging方法产生每棵决策树的训练子集；
设定随机森林由多棵CART决策树组成，CART决策树采用基尼Gini系数选择最优特征；其具体步骤如下：
对于特征X，根据某一属性值将其划分为两个子集，计算该属性值进行节点划分得到的基尼值，基尼值的计算公式如下所示：
其中，k表示类别个数，pk样本属于第k类的概率；
计算每个特征的基尼系数，选择基尼系数最小的特征X进行节点划分，基尼系数Gini(D,X)的计算公式如下所示：
其中，D表示给定样本集合，D1和D2是将D划分成的两个子集；
遍历特征X的所有属性值，选择基尼系数最小的属性值作为特征X的最优划分节点值；
不断遍历这棵树的特征子集，重复步骤5221)和步骤5222)直至所有的特征都被选择完毕或子数据集都属于同一类；
利用explain函数计算随机森林训练过程中所有时间序列特征的重要性得分，其中，重要性得分的计算步骤如下所示：
每个特征基尼指数GIm的计算公式为：
其中，k表示类别个数，pmk表示节点m中类别k所占的比例；
计算特征Xj在节点m的重要性即节点m分枝前后的基尼指数变化量的计算公式如下：
其中，GIl和GIr分别表示分枝后两个新节点的基尼指数；
计算特征Xj在第i颗树的重要性其计算公式如下：
其中，n表示为决策树的棵树，集合M表示特征Xj在决策树i中出现的节点，m为M中的某个节点；
计算特征Xj最终的重要性得分，其计算公式如下：
其中，c表示特征的数量；
在GEE云计算平台中调用ee.Reducer.spearmansCorrelation()算法，分别计算Sentinel-2和Sentinel-1候选时间序列特征中每两个候选特征之间的斯皮尔曼相关系数，斯皮尔曼相关系数的计算公式如下所示：
其中，分子为两个序列之间的误差之和，反映两个候选特征之间的差异；分母为与序列长度相关的一个常数；
计算过程为：对于维度为n的两个候选特征X和Y，Xi、Yi分别表示它们对应的第i(1≤i≤n)个特征值；对X和Y按照升序或降序的相同方式进行排列，得到新的序列x、y，其中，元素xi为Xi在X中的排行、yi为Yi在Y中的排行，获得差分集合di＝xi-yi，最后代入公式求解；
令优选特征子集从空集开始，按时间序列特征重要性分数降序的顺序，首先将最重要的候选特征添加进特征子集中，从第二重要的候选特征开始，查看其与优选特征子集中已含特征之间的斯皮尔曼相关系数，若与优选特征子集中已含特征的相关系数大于0.5或小于-0.5，则舍弃当前候选特征，以此类推直至最后一个候选特征，最终得到优选特征子集。
所述大豆种植区识别包括以下步骤：
将优选特征子集输入支持向量机模型进行训练，其中训练具体步骤如下：
令φ(x)表示将样本x映射后的特征向量，构造在特征空间中划分超平面所对应的模型f(x)表示为：
其中，wT为法向量，b为位移项，；
计算超平面最大间隔，将||w||2最小化：
其中，为最小化函数，s.t.为该最小化函数中参数w、b的约束条件，m为样本数量，yi为某一样本对应的类别；
该公式表示在约束条件下求||w||2最小的参数对w、b；
构造线性核函数：
其中，k(xi，xj)为线性核函数，xi，xj为所选取的样本值，是第i个样本映射后的特征向量，/>是第j个样本映射后的特征向量；
构建拉格朗日函数：
其中，w和b为模型的参数，α＝(α1；α2；...；αm)为拉格朗日乘子，m为样本数量，yi为某第i个样本对应的类别，是第i个样本映射后的特征向量；
利用拉格朗日乘子法，求对偶问题：
其中，α＝(α1；α2；...；αm)为拉格朗日乘子，m为样本数量，yi，yi为样本类别，k(xi，xj)为线性核函数，为最小化函数，s.t./>αi≥0(i＝1，2，…，m)为该最小化函数中α的约束条件；
得到最优的超平面F(x)
其中，αi为拉格朗日乘子，k(xi，xj)为线性核函数，m为样本数量，b是参数；
利用最优的超平面将分类样本分布在这个平面的两侧，得到最终分类结果；
利用eq()函数提取分类结果中的大豆种植区。
基于Sentinel-1、2影像的大豆种植区的面积测算方法，在GEE云计算平台中调用ee.Image.pixelArea()算法计算大豆种植区面积，计算每个大豆像元面积并累加得到所有大豆像元总面积；测算出大豆种植面积。
有益效果
本发明的基于Sentinel-1、2影像的大豆种植区识别方法，与现有技术相比借助GEE云计算平台，利用线性谐波模型提取大豆生长季内Sentinel-1、2影像的时间序列特征，然后构造支持向量机模型，同时结合随机森林分类模型及斯皮尔曼相关系数探究大豆识别优选特征子集，最终利用支持向量机模型提取大豆种植区并测算面积。
本发明经研究区结果分析如下：(1)分类总体精度为93.75％，Kappa系数为0.74，大豆F1指数为0.96，大豆种植区测算面积为923.32千公顷；(2)利用特征筛选后的优选特征子集，可以提升分类精度的同时减少65％的输入数据量。
附图说明
图1为本发明的方法顺序图；
图2为耕地数据以及集样本集分布图；
图3为2020.10.1-2021.05.31期间Sentinel-2影像(RGB波段合成)图；
图4为2020.10.1-2021.05.31期间Sentinel-1影像(VV波段)图；
图5为本发明所示所述特征提取方法提取的EVI时间序列特征图；
图6为利用本发明所述方法所生成的候选时间序列特征重要性降序排列图；
图7为利用本发明所述方法所生成的候选时间序列特征之间斯皮尔曼相关系数图，其中：(a)为Sentinel-2候选时间序列特征图，(b)为Sentinel-1候选时间序列特征图；
图8为利用本发明所述方法所生成的2020/2021生长季巴西圣卡塔琳娜州大豆种植区图。
具体实施方式
为使对本发明的结构特征及所达成的功效有更进一步的了解与认识，用以较佳的实施例及附图配合详细的说明，说明如下：
随机森林分类模型能够处理高维遥感数据，并评估每个特征重要性；然而，随机森林分类模型本身对数据噪声不敏感从而可能导致过拟合。支持向量机模型在小样本训练集上优于其他分类模型，通过优化目标使结构化风险最小，避免了过拟合问题。借助随机森林分类模型基于基尼指数进行高维特征重要性评估，能够准确地得到每个特征的重要性分数，再结合斯皮尔曼相关系数进行特征优选，能够减少冗余信息获得重要且低相关的低维特征子集，最终利用支持向量机模型实现准确、快捷地分类。
如图1所示，本发明所述的基于Sentinel-1、2影像的大豆种植区识别方法，其包括以下步骤：
第一步，Sentinel-1、2影像的获取和预处理。获取大豆生长季内Sentinel-2卫星影像，进行预处理后计算5个植被指数并选择5个原始波段；获取大豆生长季内Sentinel-1影像，并进行预处理后选择2个原始波段；分别对所选的植被指数及原始波段进行耕地掩膜处理。
(1)在GlobeLand30网站下载30m空间分辨率的土地覆盖数据，在ArcGIS软件中进行镶嵌及剪裁，然后利用ArcGIS软件中提取工具获得耕地数据，并上传至GEE云计算平台。
(2)在GEE云计算平台中调用大豆生长季内所有的Sentinel-2卫星大气底端反射率数据，利用质量控制QA60波段将云量限制在10％以内，并将所有波段值除以10000得到各个波段反射率值。
(3)计算大豆生长季内Sentinel-2卫星影像的增强型植被指数EVI、归一化植被指数NDVI、归一化耕作指数NDTI、归一化差值衰变指数NDSVI、红边位置指数REPI，共5个植被指数；
并选择大豆生长季内Sentinel-2卫星影像的红边1、红边2、红边3、短波红外1和短波红外2，共5个原始波段；
利用耕地数据对大豆生长季内Sentinel-2卫星影像的5个植被指数和5个原始波段进行耕地掩膜，
其中Sentinel-2卫星影像的植被指数计算公式如下所示：
其中，Red、Blue、NIR分别是Sentinel-2的红光、蓝光和近红外波段反射率值；
其中，NIR和Red分别是Sentinel-2的近红外和红光波段反射率值；
其中，SWIR1和SWIR2分别是Sentinel-2的短波红外1和短波红外2波段反射率值；
其中，SWIR1和Red分别是Sentinel-2的短波红外1和红光波段反射率值；
其中，Red是Sentinel-2的红光波段反射率值，RE1、RE2和RE3是Sentinel-2的三个红边波段即红边1、红边2、红边3波段反射率值。
(4)在GEE云计算平台中调用大豆生长季内所有的Sentinel-1卫星Level-1地距影像数据，成像方式为干涉测量宽幅模式，利用Sentinel-1工具箱进行预处理；选择Sentinel-1卫星影像的VV、VH波段，共2个原始波段；利用耕地数据对大豆生长季内Sentinel-1卫星影像的2个原始波段进行掩膜处理。
第二步，时间序列特征提取。利用2阶线性谐波模型提取大豆生长季内掩膜后的Sentinel-2卫星影像的5个取植被指数及5个原始波段的时间序列特征，作为Sentinel-2候选时间序列特征；利用2阶线性谐波模型提取大豆生长季内掩膜后的Sentinel-1卫星影像的2个原始波段的时间序列特征，作为Sentinel-1候选时间序列特征；汇总Sentinel-2和Sentinel-1所有候选时间序列特征，得到时间序列特征。
线性谐波模型可填补由于噪声等导致的陆地卫星缺失的观测值，已被广泛应用于提取时间序列特征。其优点是可以在GEE云计算平台中直接调用内置的线性回归函数实现，特别是在大区域尺度上能够方便快捷的应用。其中2阶线性谐波模型较为常用。
时间序列特征提取具体步骤如下：
(1)利用2阶线性谐波模型提取大豆生长季内掩膜后的Sentinel-2卫星影像的5个植被指数和5个原始波段的时间序列特征；利用2阶线性谐波模型提取大豆生长季内掩膜后的Sentinel-1卫星影像的2个原始波段的时间序列特征；其中，2阶线性谐波模型计算公式如下：
f(t)＝b+C1cos(3πt)+D1sin(3πt)+C2cos(6πt)+D2sin(6πt)，
其中，f(t)是第t时刻分别拟合的Sentinel-2卫星影像的5个植被指数和5个原始波段值以及Sentinel-1卫星影像的2个原始波段值，b是常数项，C1、C2、D1、D2是余弦函数和正弦函数的系数，t是一年中的一天，
取相位和振幅作为各个时间序列特征输入，其中振幅1定义为二维矢量[C1，D1]的长度，振幅2定义为二维矢量[C2，D2]的长度，相位1定义为二维矢量[C1，D1]的角度，相位2定义为二维矢量[C2，D2]的角度。
(2)汇总Sentinel-2和Sentinel-1所有候选时间序列特征，得到时间序列特征。其中汇总的计算公式如下所示：
VI＝VI1.addBands(VI2).addBands(VI3)…addBands(VIx)
其中，VI1表示第一个候选时间序列特征，VIx表示为最后一个候选时间特征，.addBands()为添加特征的函数。
第三步，支持向量机模型的构建：构建支持向量机分类模型。
(1)在GEE云计算平台中调用支持向量机分类模型，选择核函数类型及其参数。
(2)设定核函数类型选择线性核函数，cost设置为1，其他参数保持默认值。
第四步，优选特征子集确定。构建随机森林分类模型，利用随机森林分类模型计算所有时间序列特征的重要性得分；分别计算Sentinel-2和Sentinel-1候选时间序列特征中每两个候选特征之间的斯皮尔曼相关系数，按照时间序列特征的重要性降序去除斯皮尔曼相关系数大于0.5或小于-0.5的候选特征，得到优选特征子集。
优选特征子集确定包括以下步骤：
(1)在GEE云计算平台中调用随机森林分类模型，设置的参数包括：决策树的数量、最大特征数、决策树深度、叶子节点最少样本数、节点划分的最小样本数，其中将决策树的数量设置为300，其他参数保持默认值。
(2)将时间序列特征输入随机森林分类模型进行训练，其步骤如下：
A1)基于时间序列特征集，采用随机且有放回地抽样Bagging方法产生每棵决策树的训练子集；
A2)设定随机森林由多棵CART决策树组成，CART决策树采用基尼Gini系数选择最优特征；其具体步骤如下：
A21)对于特征X，根据某一属性值将其划分为两个子集，计算该属性值进行节点划分得到的基尼值，基尼值的计算公式如下所示：
其中，k表示类别个数，pk样本属于第k类的概率；
A22)计算每个特征的基尼系数，选择基尼系数最小的特征X进行节点划分，基尼系数Gini(D，X)的计算公式如下所示：
其中，D表示给定样本集合，D1和D2是将D划分成的两个子集；
A23)遍历特征X的所有属性值，选择基尼系数最小的属性值作为特征X的最优划分节点值；
A24)不断遍历这棵树的特征子集，重复步骤5221)和步骤5222)直至所有的特征都被选择完毕或子数据集都属于同一类；
A3)利用explain函数计算随机森林训练过程中所有时间序列特征的重要性得分，其中，重要性得分的计算步骤如下所示：
A31)每个特征基尼指数GIm的计算公式为：
其中，k表示类别个数，pmk表示节点m中类别k所占的比例；
A32)计算特征Xj在节点m的重要性即节点m分枝前后的基尼指数变化量的计算公式如下：
其中，GIl和GIr分别表示分枝后两个新节点的基尼指数；
A33)计算特征Xj在第i颗树的重要性其计算公式如下：
其中，n表示为决策树的棵树，集合M表示特征Xj在决策树i中出现的节点，m为M中的某个节点；
A34)计算特征Xj最终的重要性得分，其计算公式如下：
其中，c表示特征的数量；
(4)在GEE云计算平台中调用ee.Reducer.spearmansCorrelation()算法，分别计算Sentinel-2和Sentinel-1候选时间序列特征中每两个候选特征之间的斯皮尔曼相关系数，斯皮尔曼相关系数的计算公式如下所示：
其中，分子为两个序列之间的误差之和，反映两个候选特征之间的差异；分母为与序列长度相关的一个常数；
计算过程为：对于维度为n的两个候选特征X和Y，Xi、Yi分别表示它们对应的第i(1≤i≤n)个特征值；对X和Y按照升序或降序的相同方式进行排列，得到新的序列x、y，其中，元素xi为Xi在X中的排行、yi为Yi在Y中的排行，获得差分集合di＝xi－yi，最后代入公式求解。
(5)令优选特征子集从空集开始，按时间序列特征重要性分数降序的顺序，首先将最重要的候选特征添加进特征子集中，从第二重要的候选特征开始，查看其与优选特征子集中已含特征之间的斯皮尔曼相关系数，若与优选特征子集中已含特征的相关系数大于0.5或小于-0.5，则舍弃当前候选特征，以此类推直至最后一个候选特征，最终得到优选特征子集。
利用随机森林分类模型基于平均不纯度下降的方法评估每个候选特征的重要性。平均不纯度下降是一种基于基尼系数进行特征重要性评估的方法，在随机森林中，当训练决策树的时候，可以计算出每个特征减少了多少树的不纯度，对于一个决策树森林来说，可以计算出每个特征平均减少的不纯度，并把它平均减少的不纯度作为特征选择的值。其值越大则认为该特征的分类能力越强，在模型中的重要性越大，反之亦然。然后结合斯皮尔曼相关系数，减少特征之间的信息冗余。斯皮尔曼相关系数是用于评价两个变量之间相关性的一种统计学方法，它利用单调方程评价两个统计变量的相关性。如果数据中没有重复值，并且当两个变量完全单调相关时，斯皮尔曼相关系数则为+1或-1，其最为显著的特点是无需考察变量的样本规模或总体分布特性，具有快捷、稳健的特点.
机器学习在实际应用时，输入特征个数越多，不仅会增加模型计算开销，同时还可能引起“维度灾难”。特征筛选能够降低冗余信息，缩减特征数量，减少计算运行时间。本发明采用结合特征重要性及斯皮尔曼相关系数共同确定优选特征子集。
第五步，大豆种植区识别：将优选特征子集输入支持向量机分类模型训练，得到分类结果并提取出大豆种植区。
(1)将优选特征子集输入支持向量机模型进行训练，其中训练具体步骤如下：
B1)令φ(x)表示将样本x映射后的特征向量，构造在特征空间中划分超平面所对应的模型f(x)表示为：
其中，wT为法向量，b为位移项，；
B2)计算超平面最大间隔，将||w||2最小化：
其中，为最小化函数，S.t.为该最小化函数中参数w、b的约束条件，m为样本数量，yi为某一样本对应的类别；
该公式表示在约束条件下求||w||2最小的参数对w、b；
B3)构造线性核函数：
其中，k(xi，xj)为线性核函数，xi，xj为所选取的样本值，是第i个样本映射后的特征向量，/>是第j个样本映射后的特征向量；
B4)构建拉格朗日函数：
其中，w和b为模型的参数，α＝(α1；α2；...；αm)为拉格朗日乘子，m为样本数量，yi为某第i个样本对应的类别，是第i个样本映射后的特征向量；
B5)利用拉格朗日乘子法，求对偶问题：
其中，α＝(α1；α2；...；αm)为拉格朗日乘子，m为样本数量，yi，yi为样本类别，k(xi，xj)为线性核函数，为最小化函数，s.t./> 为该最小化函数中α的约束条件；
B6)得到最优的超平面F(x)
其中，αi为拉格朗日乘子，k(xi，xj)为线性核函数，m为样本数量，b是参数。
(2)利用最优的超平面将分类样本分布在这个平面的两侧，得到最终分类结果。
(3)利用eq()函数提取分类结果中的大豆种植区。
在此，还提供一种基于Sentinel-1、2影像的大豆种植区的面积测算方法，在GEE云计算平台中调用ee.Image.pixelArea()算法计算大豆种植区面积，计算每个大豆像元面积并累加得到所有大豆像元总面积；测算出大豆种植面积。
下面结合仿真实验对本发明的效果做进一步的说明：
1.仿真实验条件：
本发明借助Google Earth Engine(GEE)平台，GEE作为云数据平台拥有强大的计算能力，不仅可以方便的调用、分析和处理各种卫星影像、地理空间数据集，而且还可以提供多种分类算法接口。用户免费注册账号后利用JavaScript编程语言即可使用。另外，选择巴西圣卡塔琳娜州作为研究区。
为了减少非耕地对分类结果的影响，在GlobeLand30网站下载30m空间分辨率的土地覆盖数据，利用ArcGIS软件处理后，得到耕地数据，如图2所示。
为了训练及验证模型分类结果的有效性，利用地面调查的样本集，此处样本类型分为2类，即大豆及非大豆。样本集分布如图2所示。
另外，图3、4分别为2020.10.1-2021.05.31期间Sentinel-2影像(RGB波段合成)、Sentinel-1影像(VV波段)。本发明采用的精度指标有基于混淆矩阵的总体精度(OA)、Kappa系数和F1指数。
2.仿真实验内容及结果分析：
表1为本发明所示所述特征提取方法提取的所有时间序列特征，共48个。图5为本发明所示所述特征提取方法提取的EVI时间序列特征。
表1时间序列特征表
图6为利用本发明所述方法所生成的候选时间序列特征重要性降序排列图，从图6可知，整体来看前6个特征相对重要，且均为线性谐波模型提取的2阶振幅特征，可见2阶振幅特征对大豆识别更加敏感。另外，除了NDTI及VH的2阶振幅特征分别与其相邻特征之间特征重要性差距较大，其他特征之间差距不显著，说明特征之间存在一定程度的信息冗余。
图7为利用本发明所述方法所生成的候选时间序列特征之间斯皮尔曼相关系数图，其中(a)为Sentinel-2候选时间序列特征，(b)为Sentinel-1候选时间序列特征；由图7(a)可见，只有1个特征(REPI的1阶相位)与其他特征之间呈弱相关性，剩余的47个特征之间均存在不同程度的相关性。充分说明了Sentinel-2时间序列特征之间存在大量的冗余。由图7(b)可见8个特征之间均存在不同程度的相关性，充分说明了Sentinel-1时间序列特征之间也存在大量的冗余。
表2为利用本发明所述方法所生成的优选特征子集，由表2可知，降维后的优选特征子集中仅仅含17个特征，相对于原始48个候选特征，数据量下降了65％。且分类总体精度为93.75％，Kappa系数为0.74，大豆F1指数为0.96。
为了验证本发明所提出的优选特征子集的性能，将原始48个候选特征直接输入支向量机进行分类，结果显示分类总体精度为90.62％，Kappa系数为0.61，大豆F1指数为0.94。由此可知优选特征子集能够将输入数据量缩减65％，同时还能提高分类精度(总体精度、Kappa系数和大豆F1指数分别提了3.45％、21.31％和2.13％)。
表2优选特征子集表
| 序号 | 优选特征子集 |
| 1 | NDVI_amplitude2 |
| 2 | VV_amplitude2 |
| 3 | VV_phase2 |
| 4 | NDSVI_amplitude1 |
| 5 | REPI_amplitude1 |
| 6 | VV_phase1 |
| 7 | B11_amplitude2 |
| 8 | B6_amplitude2 |
| 9 | NDTI_phase2 |
| 10 | EVI_phase1 |
| 11 | B11_phase1 |
| 12 | B6_phase1 |
| 13 | REPI_phase1 |
| 14 | B11_phase2 |
| 15 | B11_amplitude1 |
| 16 | NDSVI_phase1 |
| 17 | B6_amplitude1 |
图8为利用本发明所述方法所生成的2020/2021生长季巴西圣卡塔琳娜州大豆种植区图，如图8所示，利用本发明所述方法经测算大豆种植区面积为923.32千公顷。
以上显示和描述了本发明的基本原理、主要特征和本发明的优点。本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实施例和说明书中描述的只是本发明的原理，在不脱离本发明精神和范围的前提下本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明的范围内。本发明要求的保护范围由所附的权利要求书及其等同物界定。
Claims (6)
1.一种基于Sentinel-1、2影像的大豆种植区识别方法，其特征在于，包括以下步骤：
11)Sentinel-1、2影像的获取和预处理：获取大豆生长季内Sentinel-2卫星影像，进行预处理后计算5个植被指数并选择5个原始波段；获取大豆生长季内Sentinel-1影像，并进行预处理后选择2个原始波段；分别对所选的植被指数及原始波段进行耕地掩膜处理；
12)时间序列特征提取：利用2阶线性谐波模型提取大豆生长季内掩膜后的Sentinel-2卫星影像的5个取植被指数及5个原始波段的时间序列特征，作为Sentinel-2候选时间序列特征；利用2阶线性谐波模型提取大豆生长季内掩膜后的Sentinel-1卫星影像的2个原始波段的时间序列特征，作为Sentinel-1候选时间序列特征；汇总Sentinel-2和Sentinel-1所有候选时间序列特征，得到时间序列特征；
13)支持向量机模型的构建：构建支持向量机分类模型；
14)优选特征子集确定：构建随机森林分类模型，利用随机森林分类模型计算所有时间序列特征的重要性得分；分别计算Sentinel-2和Sentinel-1候选时间序列特征中每两个候选特征之间的斯皮尔曼相关系数，按照时间序列特征的重要性降序去除斯皮尔曼相关系数大于0.5或小于-0.5的候选特征，得到优选特征子集；
15)大豆种植区识别：将优选特征子集输入支持向量机分类模型训练，得到分类结果并提取出大豆种植区。
2.根据权利要求1所述的基于Sentinel-1、2影像的大豆种植区识别方法，其特征在于，所述Sentinel-1、2影像的获取和预处理包括以下步骤：
21)在GlobeLand30网站下载30m空间分辨率的土地覆盖数据，在ArcGIS软件中进行镶嵌及剪裁，然后利用ArcGIS软件中提取工具获得耕地数据，并上传至GEE云计算平台；
22)在GEE云计算平台中调用大豆生长季内所有的Sentinel-2卫星大气底端反射率数据，利用质量控制QA60波段将云量限制在10％以内，并将所有波段值除以10000得到各个波段反射率值；
23)计算大豆生长季内Sentinel-2卫星影像的增强型植被指数EVI、归一化植被指数NDVI、归一化耕作指数NDTI、归一化差值衰变指数NDSVI、红边位置指数REPI，共5个植被指数；
并选择大豆生长季内Sentinel-2卫星影像的红边1、红边2、红边3、短波红外1和短波红外2，共5个原始波段；
利用耕地数据对大豆生长季内Sentinel-2卫星影像的5个植被指数和5个原始波段进行耕地掩膜，
其中Sentinel-2卫星影像的植被指数计算公式如下所示：
其中，Red、Blue、NIR分别是Sentinel-2的红光、蓝光和近红外波段反射率值；
其中，NIR和Red分别是Sentinel-2的近红外和红光波段反射率值；
其中，SWIR1和SWIR2分别是Sentinel-2的短波红外1和短波红外2波段反射率值；
其中，SWIR1和Red分别是Sentinel-2的短波红外1和红光波段反射率值；
其中，Red是Sentinel-2的红光波段反射率值，RE1、RE2和RE3是Sentinel-2的三个红边波段即红边1、红边2、红边3波段反射率值；
24)在GEE云计算平台中调用大豆生长季内所有的Sentinel-1卫星Level-1地距影像数据，成像方式为干涉测量宽幅模式，利用Sentinel-1工具箱进行预处理；选择Sentinel-1卫星影像的VV、VH波段，共2个原始波段；利用耕地数据对大豆生长季内Sentinel-1卫星影像的2个原始波段进行掩膜处理。
3.根据权利要求1所述的基于Sentinel-1、2影像的大豆种植区识别方法，其特征在于，所述时间序列特征提取包括以下步骤：
31)利用2阶线性谐波模型提取大豆生长季内掩膜后的Sentinel-2卫星影像的5个植被指数和5个原始波段的时间序列特征；利用2阶线性谐波模型提取大豆生长季内掩膜后的Sentinel-1卫星影像的2个原始波段的时间序列特征；其中，2阶线性谐波模型计算公式如下：
f(t)＝b+C1cos(3πt)+D1sin(3πt)+C2cos(6πt)+D2sin(6πt)，
其中，f(t)是第t时刻分别拟合的Sentinel-2卫星影像的5个植被指数和5个原始波段值以及Sentinel-1卫星影像的2个原始波段值，b是常数项，C1、C2、D1、D2是余弦函数和正弦函数的系数，t是一年中的一天，
取相位和振幅作为各个时间序列特征输入，其中振幅1定义为二维矢量[C1，D1]的长度，振幅2定义为二维矢量[C2，D2]的长度，相位1定义为二维矢量[C1，D1]的角度，相位2定义为二维矢量[C2，D2]的角度；
32)汇总Sentinel-2和Sentinel-1所有候选时间序列特征，得到时间序列特征，其中汇总的计算公式如下所示：
VI＝VI1.addBands(VI2).addBands(VI3)…addBands(VIx)
其中，VI1表示第一个候选时间序列特征，VIx表示为最后一个候选时间特征，.addBands()为添加特征的函数。
4.根据权利要求1所述的基于Sentinel-1、2影像的大豆种植区识别方法，其特征在于，所述支持向量机模型的构建包括以下步骤：
41)在GEE云计算平台中调用支持向量机分类模型，选择核函数类型及其参数；
42)设定核函数类型选择线性核函数，cost设置为1，其他参数保持默认值。
5.根据权利要求1所述的基于Sentinel-1、2影像的大豆种植区识别方法，其特征在于，所述优选特征子集确定包括以下步骤：
51)在GEE云计算平台中调用随机森林分类模型，设置的参数包括：决策树的数量、最大特征数、决策树深度、叶子节点最少样本数、节点划分的最小样本数，其中将决策树的数量设置为300，其他参数保持默认值；
52)将时间序列特征输入随机森林分类模型进行训练，其步骤如下：
521)基于时间序列特征集，采用随机且有放回地抽样Bagging方法产生每棵决策树的训练子集；
522)设定随机森林由多棵CART决策树组成，CART决策树采用基尼Gini系数选择最优特征；其具体步骤如下：
5221)对于特征X，根据某一属性值将其划分为两个子集，计算该属性值进行节点划分得到的基尼值，基尼值的计算公式如下所示：
其中，k表示类别个数，pk样本属于第k类的概率；
5222)计算每个特征的基尼系数，选择基尼系数最小的特征X进行节点划分，基尼系数Gini(D,X)的计算公式如下所示：
其中，D表示给定样本集合，D1和D2是将D划分成的两个子集；
5223)遍历特征X的所有属性值，选择基尼系数最小的属性值作为特征X的最优划分节点值；
5224)不断遍历这棵树的特征子集，重复步骤5221)和步骤5222)直至所有的特征都被选择完毕或子数据集都属于同一类；
53)利用explain函数计算随机森林训练过程中所有时间序列特征的重要性得分，其中，重要性得分的计算步骤如下所示：
531)每个特征基尼指数GIm的计算公式为：
其中，k表示类别个数，pmk表示节点m中类别k所占的比例；
532)计算特征Xj在节点m的重要性即节点m分枝前后的基尼指数变化量的计算公式如下：
其中，GIl和GIr分别表示分枝后两个新节点的基尼指数；
533)计算特征Xj在第i颗树的重要性其计算公式如下：
其中，n表示为决策树的棵树，集合M表示特征Xj在决策树i中出现的节点，m为M中的某个节点；
534)计算特征Xj最终的重要性得分，其计算公式如下：
其中，c表示特征的数量；
54)在GEE云计算平台中调用ee.Reducer.spearmansCorrelation()算法，分别计算Sentinel-2和Sentinel-1候选时间序列特征中每两个候选特征之间的斯皮尔曼相关系数，斯皮尔曼相关系数的计算公式如下所示：
其中，分子为两个序列之间的误差之和，反映两个候选特征之间的差异；分母为与序列长度相关的一个常数；
计算过程为：对于维度为n的两个候选特征X和Y，Xi、Yi分别表示它们对应的第i(1≤i≤n)个特征值；对X和Y按照升序或降序的相同方式进行排列，得到新的序列x、y，其中，元素xi为Xi在X中的排行、yi为Yi在Y中的排行，获得差分集合di＝xi－yi，最后代入公式求解；
55)令优选特征子集从空集开始，按时间序列特征重要性分数降序的顺序，首先将最重要的候选特征添加进特征子集中，从第二重要的候选特征开始，查看其与优选特征子集中已含特征之间的斯皮尔曼相关系数，若与优选特征子集中已含特征的相关系数大于0.5或小于-0.5，则舍弃当前候选特征，以此类推直至最后一个候选特征，最终得到优选特征子集。
6.根据权利要求1所述的基于Sentinel-1、2影像的大豆种植区识别方法，其特征在于，所述大豆种植区识别包括以下步骤：
61)将优选特征子集输入支持向量机模型进行训练，其中训练具体步骤如下：
611)令表示将样本x映射后的特征向量，构造在特征空间中划分超平面所对应的模型f(x)表示为：
其中，wT为法向量，b为位移项；
612)计算超平面最大间隔，将||w||2最小化：
其中，为最小化函数，s.t.为该最小化函数中参数w、b的约束条件，m为样本数量，yi为某一样本对应的类别；
该公式表示在约束条件下求||w||2最小的参数对w、b；
613)构造线性核函数：
其中，k(xi,xj)为线性核函数，xi,xj为所选取的样本值，是第i个样本映射后的特征向量，/>是第j个样本映射后的特征向量；
614)构建拉格朗日函数：
其中，w和b为模型的参数，α＝(α1；α2；…；αm)为拉格朗日乘子，m为样本数量，yi为某第i个样本对应的类别，是第i个样本映射后的特征向量；
615)利用拉格朗日乘子法，求对偶问题：
其中，α＝(α1；α2；…；αm)为拉格朗日乘子，m为样本数量，yi,yi为样本类别，k(xi,xj)为线性核函数，为最小化函数，/>为该最小化函数中α的约束条件；
616)得到最优的超平面F(x)
其中，αi为拉格朗日乘子，k(xi,xj)为线性核函数，m为样本数量，b是参数；
62)利用最优的超平面将分类样本分布在这个平面的两侧，得到最终分类结果；
63)利用eq()函数提取分类结果中的大豆种植区。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210596186.6A CN115063610B (zh) | 2022-05-30 | 2022-05-30 | 基于Sentinel-1、2影像的大豆种植区识别方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210596186.6A CN115063610B (zh) | 2022-05-30 | 2022-05-30 | 基于Sentinel-1、2影像的大豆种植区识别方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN115063610A CN115063610A (zh) | 2022-09-16 |
| CN115063610B true CN115063610B (zh) | 2024-03-12 |
Family
ID=83199173
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202210596186.6A Active CN115063610B (zh) | 2022-05-30 | 2022-05-30 | 基于Sentinel-1、2影像的大豆种植区识别方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN115063610B (zh) |
Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN114494909A (zh) * | 2022-02-16 | 2022-05-13 | 中国科学院空天信息创新研究院 | 一种大豆生长季空间分布图的生成方法和系统 |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US8108324B2 (en) * | 2008-05-15 | 2012-01-31 | Intel Corporation | Forward feature selection for support vector machines |
-
2022
- 2022-05-30 CN CN202210596186.6A patent/CN115063610B/zh active Active
Patent Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN114494909A (zh) * | 2022-02-16 | 2022-05-13 | 中国科学院空天信息创新研究院 | 一种大豆生长季空间分布图的生成方法和系统 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN115063610A (zh) | 2022-09-16 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Solano et al. | A methodology based on GEOBIA and WorldView-3 imagery to derive vegetation indices at tree crown detail in olive orchards | |
| Pott et al. | Satellite-based data fusion crop type classification and mapping in Rio Grande do Sul, Brazil | |
| Peña-Barragán et al. | Object-based crop identification using multiple vegetation indices, textural features and crop phenology | |
| CN105678281B (zh) | 基于光谱和纹理特征的地膜覆盖农田遥感监测方法 | |
| CN108458978B (zh) | 基于敏感波段和波段组合最优的树种多光谱遥感识别方法 | |
| CN111598045A (zh) | 一种基于对象图谱和混合光谱的遥感耕地变化检测方法 | |
| CN112861435B (zh) | 一种红树林质量遥感反演方法及智能终端 | |
| Shi et al. | Evaluation of wavelet spectral features in pathological detection and discrimination of yellow rust and powdery mildew in winter wheat with hyperspectral reflectance data | |
| CN114821362A (zh) | 一种基于多源数据的水稻种植面积提取方法 | |
| Liang et al. | Improved estimation of aboveground biomass in rubber plantations by fusing spectral and textural information from UAV-based RGB imagery | |
| CN116645603A (zh) | 一种大豆种植区识别和面积测算方法 | |
| CN107688777A (zh) | 一种协同多源遥感影像的城市绿地提取方法 | |
| CN112836725A (zh) | 基于时序遥感数据的弱监督lstm循环神经网络稻田识别方法 | |
| Xu et al. | A 21-year time series of global leaf chlorophyll content maps from MODIS imagery | |
| Xu et al. | Retrieving global leaf chlorophyll content from MERIS data using a neural network method | |
| Yuan et al. | Research on rice leaf area index estimation based on fusion of texture and spectral information | |
| CN116883853B (zh) | 基于迁移学习的作物时空信息遥感分类方法 | |
| CN116665073A (zh) | 一种基于多源数据的玉米产量遥感估测方法 | |
| CN115063610B (zh) | 基于Sentinel-1、2影像的大豆种植区识别方法 | |
| CN115953685A (zh) | 多层多尺度分割的农用大棚类型信息提取方法及系统 | |
| CN112577954B (zh) | 城市绿地地上生物量估测方法 | |
| Lottering et al. | Spatially optimizing vegetation indices integrated with sparse partial least squares regression to detect and map the effects of Gonipterus scutellatus on the chlorophyll content of eucalyptus plantations | |
| You et al. | Crop Mapping of Complex Agricultural Landscapes Based on Discriminant Space | |
| Dutta | Assessment of tea bush health and yield using geospatial techniques | |
| Zhang et al. | Estimation of grassland height using optical and SAR remote sensing data |
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 |