CN113205014B - 一种基于图像锐化的时序数据耕地提取方法 - Google Patents
一种基于图像锐化的时序数据耕地提取方法 Download PDFInfo
- Publication number
- CN113205014B CN113205014B CN202110421926.8A CN202110421926A CN113205014B CN 113205014 B CN113205014 B CN 113205014B CN 202110421926 A CN202110421926 A CN 202110421926A CN 113205014 B CN113205014 B CN 113205014B
- Authority
- CN
- China
- Prior art keywords
- image
- data
- band
- time sequence
- vegetation index
- 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
- 238000000605 extraction Methods 0.000 title claims abstract description 25
- 238000003707 image sharpening Methods 0.000 title claims abstract description 23
- 238000000034 method Methods 0.000 claims abstract description 24
- 239000002689 soil Substances 0.000 claims abstract description 19
- 238000007781 pre-processing Methods 0.000 claims abstract description 11
- 238000007637 random forest analysis Methods 0.000 claims abstract description 11
- 238000005070 sampling Methods 0.000 claims abstract description 9
- 230000009466 transformation Effects 0.000 claims description 20
- 238000012549 training Methods 0.000 claims description 16
- 238000012795 verification Methods 0.000 claims description 15
- 238000004364 calculation method Methods 0.000 claims description 14
- 208000027066 STING-associated vasculopathy with onset in infancy Diseases 0.000 claims description 11
- 238000011160 research Methods 0.000 claims description 9
- 238000012937 correction Methods 0.000 claims description 8
- 238000012545 processing Methods 0.000 claims description 8
- 238000012360 testing method Methods 0.000 claims description 8
- 238000004422 calculation algorithm Methods 0.000 claims description 6
- 238000011156 evaluation Methods 0.000 claims description 6
- 230000005855 radiation Effects 0.000 claims description 6
- 238000004088 simulation Methods 0.000 claims description 6
- 230000003595 spectral effect Effects 0.000 claims description 6
- 238000006243 chemical reaction Methods 0.000 claims description 5
- 238000010219 correlation analysis Methods 0.000 claims description 5
- 238000005520 cutting process Methods 0.000 claims description 4
- 238000009826 distribution Methods 0.000 claims description 4
- 238000012986 modification Methods 0.000 claims description 4
- 230000004048 modification Effects 0.000 claims description 4
- 238000005316 response function Methods 0.000 claims description 4
- DIWRORZWFLOCLC-UHFFFAOYSA-N Lorazepam Chemical compound C12=CC(Cl)=CC=C2NC(=O)C(O)N=C1C1=CC=CC=C1Cl DIWRORZWFLOCLC-UHFFFAOYSA-N 0.000 claims description 3
- 238000013075 data extraction Methods 0.000 abstract description 3
- 238000012544 monitoring process Methods 0.000 description 22
- 238000005516 engineering process Methods 0.000 description 6
- 230000004927 fusion Effects 0.000 description 6
- 238000010276 construction Methods 0.000 description 3
- 238000011161 development Methods 0.000 description 3
- 230000000694 effects Effects 0.000 description 3
- 238000010801 machine learning Methods 0.000 description 3
- 230000006872 improvement Effects 0.000 description 2
- 238000004519 manufacturing process Methods 0.000 description 2
- 238000012216 screening Methods 0.000 description 2
- 230000001932 seasonal effect Effects 0.000 description 2
- 238000007792 addition Methods 0.000 description 1
- 238000013473 artificial intelligence Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000003750 conditioning effect Effects 0.000 description 1
- 238000012790 confirmation Methods 0.000 description 1
- 238000003066 decision tree Methods 0.000 description 1
- 230000007812 deficiency Effects 0.000 description 1
- 230000002950 deficient Effects 0.000 description 1
- 238000013467 fragmentation Methods 0.000 description 1
- 238000006062 fragmentation reaction Methods 0.000 description 1
- 238000011835 investigation Methods 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
- 238000002310 reflectometry Methods 0.000 description 1
- 230000002123 temporal effect Effects 0.000 description 1
- 238000012876 topography Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
- G06V20/00—Scenes; Scene-specific elements
- G06V20/10—Terrestrial scenes
- G06V20/188—Vegetation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F16/00—Information retrieval; Database structures therefor; File system structures therefor
- G06F16/50—Information retrieval; Database structures therefor; File system structures therefor of still image data
- G06F16/58—Retrieval characterised by using metadata, e.g. metadata not derived from the content or metadata generated manually
- G06F16/583—Retrieval characterised by using metadata, e.g. metadata not derived from the content or metadata generated manually using metadata automatically derived from the content
- G06F16/5838—Retrieval characterised by using metadata, e.g. metadata not derived from the content or metadata generated manually using metadata automatically derived from the content using colour
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F16/00—Information retrieval; Database structures therefor; File system structures therefor
- G06F16/50—Information retrieval; Database structures therefor; File system structures therefor of still image data
- G06F16/58—Retrieval characterised by using metadata, e.g. metadata not derived from the content or metadata generated manually
- G06F16/5866—Retrieval characterised by using metadata, e.g. metadata not derived from the content or metadata generated manually using information manually generated, e.g. tags, keywords, comments, manually generated location and time information
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F18/00—Pattern recognition
- G06F18/20—Analysing
- G06F18/21—Design or setup of recognition systems or techniques; Extraction of features in feature space; Blind source separation
- G06F18/214—Generating training patterns; Bootstrap methods, e.g. bagging or boosting
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F18/00—Pattern recognition
- G06F18/20—Analysing
- G06F18/24—Classification techniques
- G06F18/243—Classification techniques relating to the number of classes
- G06F18/24323—Tree-organised classifiers
-
- G06T5/73—
-
- G06T5/80—
-
- G06T5/90—
-
- 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
- G06T2207/10036—Multispectral image; Hyperspectral 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
- G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
- G06V20/00—Scenes; Scene-specific elements
- G06V20/10—Terrestrial scenes
- G06V20/194—Terrestrial scenes using hyperspectral data, i.e. more or other wavelengths than RGB
-
- 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
Abstract
本发明公开了一种基于图像锐化的时序数据耕地提取方法，包括如下步骤：1)输入单时相高空间分辨率遥感影像全色波段并进行数据预处理；2)输入多时相Sentinel‑2遥感影像多光谱波段并进行数据预处理；3)基于单时相高空间分辨率全色波段和多时相Sentinel‑2多光谱波段进行Gram‑Schmidt图像锐化操作，生成高空间分辨率时序数据；4)基于高空间分辨率时序数据，计算归一化差异植被指数和土壤调节植被指数，生成植被指数时间序列数据集；5)构建分类样本库并进行分层抽样；6)利用随机森林的分类方法对作物类型进行识别，并对比原始时间序列数据提取结果与锐化后的时间序列提取结果；7)将提取的所有类型的耕地进行整合，生成耕地要素提取结果。
Description
技术领域
本发明涉及一种遥感农情监测方法，特别是指一种基于图像锐化的时序数据耕地提取方法。
背景技术
以土地资源为生产对象的农业作为国家的第一产业，支撑着国民经济的建设与发展，也是联合国可持续发展的重要组成部分。实时准确地监测耕地也是21世纪粮食安全的关键所在。耕地资源要素空间分布是耕地管理、精准农业以及土地确权的基本要求。然而，传统的耕地资源监测主要是以实地调查监测的方式为主，这种方式耗时耗力，满足不了大范围、快速准确监测的要求。随着遥感技术的发展和各种智能化技术的提高，越来越多的生产部门将遥感作为主要是监测方式进行耕地资源监测，其周期性、时效性以及经济性的特点受到很大的利用。因此，以遥感技术为基础的耕地资源监测是今后国情监测、土地权属分配等领域的重要支撑技术。
当前，基于遥感技术的耕地资源监测受到诸多因素的影响：1)复杂自然环境中，耕地要素在区域之间的差异明显，特别是不同区域的地形地貌的不同，导致耕地资源要素在遥感影像上呈现完全不一致的特点。例如：东北地区主要以平原为主，地块规整；而南方低山丘陵区地势连绵起伏，小农特征明显，耕地破碎化程度严重，边界模糊；2)耕地具有典型的季节性特征，相同区域内由于作物类型与种植模式的差异，导致区域内耕地类型异质性较大，复杂的物候特征增加了其监测的难度；3)高时空遥感数据可用性不足。尽管越来越多的数据源能够为农业监测带来契机，但是遥感数据却依然无法在时间和空间分辨率上达到一种平衡
多源遥感影像时空融合虽然可以在一定程度上解决遥感数据在时间分辨率和空间分辨率上难以统一的问题，但是多源数据间跨尺度和跨模态匹配问题仍是一大挑战，并且时空融合结果大多是构建中分辨率数据集，而高空间分辨率时序数据集仍然很匮乏。随着各种智能化技术的提高，以机器学习为主的人工智能方法逐渐适应了遥感领域，此类方法模型具有良好的鲁棒性，并且能产生相对较好的监测结果，但是依然受到农业监测最主要的一个问题：真实样本数据缺乏，特别是季节性描述样本的不足。因此，针对耕地监测以及当今农业监测领域所面临的耕地监测难、监测粒度不精确、数据的高时空特性不足等问题，通过融合多源不同尺度多源遥感数据，构建高空间率时序数据集，基于机器学习的方法实现高精度、低成本、大范围的耕地资源要素遥感监测研究，是当前农业遥感研究和应用重要一项内容。
发明内容
本发明的目的在于提供一种能够有效进行耕地资源监测的基于图像锐化的时序数据耕地要素识别方法。
为实现上述目的，本发明所提供的一种基于图像锐化的时序数据耕地提取方法，包括如下步骤：
1)输入单时相高空间分辨率遥感影像全色波段并进行数据预处理；预处理步骤包括正射校正、辐射校正、几何配准、影像裁剪等。在几何配准时，需要手工选择若干控制点与辅助数据(例如谷歌数据等)进行几何配准。
2)输入L2A级别的多时相时间序列Sentinel-2影像并对其进行批量预处理，即通过欧空局SNAP平台进行格式转换输出与批量裁剪处理；
3)对于步骤1)和步骤2)的结果，基于单时相高空间分辨率全色波段和多时相Sentinel-2多光谱波段进行Gram-Schmidt图像锐化操作，生成高空间分辨率时序数据；
4)对于步骤3)中的高空间分辨率时序数据，计算归一化差异植被指数和土壤调节植被指数，生成高空间分辨率植被指数时间序列数据集；并通过对比锐化前与锐化后的时序数据集进行植被指数相关性分析。
5)根据步骤4)得到的高空间分辨率植被指数时间序列数据集，结合野外采样各地物的样本点，构建分类样本库并进行分层抽样；
6)利用步骤5)中构建的分类样本库，先利用随机森林算法进行模型训练，通过测试集的表现确定最优模型参数；将训练好的模型应用至全区域，实现全区域作物类型的识别，并对比原始时间序列数据提取结果与锐化后的时间序列提取结果；
7)利用步骤6)得到的作物类型的结果，将提取的所有类型的作物进行整合，形成耕地要素数据，将所有的非作物类型进行整合，形成非耕地数据，并进行精度验证；最终生成耕地资源要素空间分布数据。
优选地，所述步骤1)中，预处理包括辐射校正、正射校正、几何校正、影像融合以及影像裁剪。
优选地，所述步骤1)中，高空间分辨率全色波段影像预处理包括辐射校正、正射校正、几何配准以及影像裁剪。在进行几何配准时，需要依托于高精度谷歌数据、天地图等开源数据，手动选择若干控制点进行几何纠正。
优选地，所述步骤2)中，输入时间序列L2A级别的Sentinel-2影像，并进行批量影像裁剪处理；其中，输入的Sentinel-2的影像必须经过严格的云量控制，例如：通过整景云量小于5％或者研究区云量小于2％来进行数据的筛选；所利用到的时序Sentinel-2是由蓝、绿、红以及近红外四个波段组成，空间分辨率为10米。
优选地，所述步骤3)中具体包括以下步骤：
3.1)使用多波段低分辨率影像对单波段高分辨率影像进行模拟。模拟的方法根据权重Wi进行模拟，即模拟的全色波段影像灰度值
(Bi为多光谱低空间分辨率影像第i波段灰度值,Wi为各波段的光谱响应函数对应的权重，k为影像的波段的数量)；模拟的高分辨率波段影像在后面的处理中被称作Gram-Schmidt第一分量进行Gram-Schmidt变换；
3.2)利用模拟的高分辨率波段影像作为Gram-Schmidt变换第一分量来对模拟的高分辨波段影像和低分辨率波段影像进行Gram-Schmidt变换。在进行Gram-Schmidt变换时，第T个变换分量由前T-1个变换分量构造，即：
其中，GST是GS变换后产生的第T个正交分量，BT是原始低空间分辨率多光谱影像的第T个波段影像，uT是第T个原始低空间分辨率多光谱波段影像灰度值的均值，i,j分别表示低空间分辨率影像的行数和列数，
为原始低空间分辨率遥感影像第T个波段与GS变换后产生的第l个分量间的协方差，均值uT公式如下：
其中，C是影像的列数，R是影像行数，协方差
计算公式如下：
其中，标准差σ的计算公式如下：
3.3)根据GS第一分量，即模拟波段的均值u和标准差σ对全色波段进行修改，修改公式如下：
P′＝(P×k1)+k2
k2＝uintensity-(k1×upan)
其中：P为原始全色波段影像的灰度值，eintensity和epan分别为GS第一分量和全色波段影像灰度值的标准差。
3.4)最后将修改后的全色波段作为第一分量，进行GS逆变换，输出N+1个波段，去除第一个波段，就是融合后的高分辨率多光谱影像结果，逆变换的计算公式如下：
优选地，所述步骤4)具体包括如下步骤：
4.1)基于高空间时序数据进行归一化差异植被指数和土壤调节植被指数计算，得到植被指数高空间分辨率时序数据集；所述的两种植被指数的计算公式如下：
归一化差异植被指数NDVI：
NDVI＝(Nir-Red)/(Nir+Red)
土壤调节植被指数SAVI：
SAVI＝(Nir-Red)*(1+L)/(Nir+Red+L)
式中，Nir为近红外波段地表反射率值，Red为红波段地表反射率值，L是随着植被密度变化的参数，取值范围是从0～1，当植被覆盖度很高时为0，很低时为1，通常取0.5时SAVI消除土壤反射率的效果最好。
4.2)将原始的时间序列数据进行归一化差异植被指数和土壤调节植被指数的计算，得到植被指数原始时序数据集；
4.3)随机选取区域内1～5万个点，对锐化前与锐化后不同植被指数进行相关性分析；随机选点是通过ArcGIS的随机产生点的工具进行，并通过多值提取至点，进行不同时相的点之间的相关性大小；
优选地，所述步骤5)，主要是通过野外采集的各土地利用类型(主要是作物类型)的样本点，建立该区域作物类型样本库，按照一定的比例进行分层抽样，获取训练数据集与验证数据集。
优选地，所述步骤6)具体包括如下步骤：
6.1)将训练数据集进行划分，抽取10％作为训练模型的测试数据，根据测试数据集的准确率，选择合适的训练参数，例如：随机树的数量，叶子节点的数量以及最大叶深度等；
6.2)将训练好的模型应用至整个研究区，完成作物类型分类，包括锐化后的高空间时序数据集和锐化前的时序数据集等。
优选地，所述步骤7)具体包括如下步骤：
7.1)将不同的作物类型进行整合，形成耕地数据；将非作物类型整合，形成非耕地数据。
7.2)将验证集中的不同作物类型的样本点合并成一类，即耕地样本点，并进行精度评价。
7.3)如精度验证结果如果与期望精度不符，则返回步骤5)重新进行分层抽样并重复步骤6)与步骤7)，直至与期望精度相符，最终完成耕地类型识别。
与现有的技术相比，本发明的有益效果在于：本发明基于单时相高空间分辨率全色波段影像和时间序列Sentinel-2数据，突破了单一数据源在空间分辨率或时间分辨率上的制约，大大提高了耕地监测的颗粒度，经过图像锐化产生的高时空时序数据集既满足了遥感耕地资源要素监测所需要的季节性信息，又实现了亚米级别的要素类型提取；同时结合鲁棒性好，泛化能力强的随机森林算法，解决了耕地要素监测所面临的高时空数据不足、提取精度较低以及监测颗粒度较粗等瓶颈问题。图像锐化时序数据集的构建从高时空数据融合的角度解决了遥感影像在空间和时间分辨率不平衡的问题，在保证精度的情况下，实现作物类型与耕地要素的精细化识别。
附图说明
图1为本发明所提供的一种基于图像锐化的时序数据耕地提取方法的流程示意图。
具体实施方式
下面结合附图对本发明作进一步的详细说明。
如图1所示，本发明所提供的一种基于图像锐化的时序数据耕地提取方法，包括如下步骤：
1)输入单时相高空间分辨率全色波段地表反射率遥感影像，对该数据进行辐射校正、正射校正、几何配准以及影像裁剪工作；其中，单时相高空间分辨率全色波段遥感影像为空间分辨率1米及1米以内的四波段影像，并且获取时间要在年周期内。辐射校正通过不同传感器的定标系数来完成；正射校正利用全球30米数字高程模型数据辅助完成；几何配准需要基于谷歌地球、天地图等高分辨率无偏移辅助数据，手动选择地面控制点完成。
2)输入时间序列L2A级别的Sentinel-2影像，并进行批量影像裁剪处理；其中，多时相时间序列影像是指在年周期内获取不同月份的遥感影像，时相范围应考虑研究区不同作物类型的生长物候时间。其中，输入的Sentinel-2的影像必须经过严格的云量控制，例如：通过整景云量小于5％或者研究区云量小于2％来进行数据的筛选；所利用到的时序Sentinel-2是由蓝、绿、红以及近红外四个波段组成，空间分辨率为10米；
3)基于单时相高空间分辨率全色波段和多时相Sentinel-2多光谱波段进行Gram-Schmidt图像锐化操作，生成高空间分辨率时序数据，Gram-Schmidt图像锐化改进主成分变换算法中信息过分集中的问题，不受波段限制，较好的保持空间纹理信息，尤其能高保真保持光谱特征。专为最新高空间分辨率影像设计，能较好保持影像的纹理和光谱信息。Gram-Schmidt图像锐化包括如下步骤：
3.1)使用多波段低分辨率影像对单波段高分辨率影像进行模拟。模拟的方法根据光谱响应函数Wi进行模拟，即模拟的全色波段影像灰度值
(Bi为多光谱地空间分辨率影像第i波段灰度值，Bi为多光谱低空间分辨率影像第i波段灰度值,Wi为各波段的光谱响应函数对应的权重，k为影像的波段的数量)；模拟的高分辨率波段影像在后面额处理中被称作Gram-Schmidt第一分量进行Gram-Schmidt变换；
3.2)利用模拟的高分辨率波段影像作为Gram-Schmidt变换第一分量来对模拟的高分辨波段影像和低分辨率波段影像进行Gram-Schmidt变换。在进行Gram-Schmidt变换时，第T个变换分量由前T-1个变换分量构造，即：
其中，GST是GS变换后产生的第T个正交分量，BT是原始低空间分辨率多光谱影像的第T个波段影像，uT是第T个原始低空间分辨率多光谱波段影像灰度值的均值，i,j分别表示低空间分辨率影像的行数和列数，
为原始低空间分辨率遥感影像第T个波段与GS变换后产生的第l个分量间的协方差，均值uT公式如下：
其中，C是影像的列数，R是影像行数，协方差
计算公式如下：
其中，标准差σ的计算公式如下：
3.3)根据GS第一分量，即模拟波段的均值u和标准差σ对全色波段进行修改，修改公式如下：
P′＝(P×k1)+k2
k2＝uintensity-(k1×upan)
其中：P为原始全色波段影像的灰度值，eintensity和epan分别为GS第一分量和全色波段影像灰度值的标准差。
3.4)最后将修改后的全色波段作为第一分量，进行GS逆变换，输出N+1个波段，去除第一个波段，就是融合后的高分辨率多光谱影像结果，逆变换的计算公式如下：
4)基于高空间分辨率时序数据，计算归一化差异植被指数和土壤调节植被指数，生成植被指数时间序列数据集。其中，归一化差异植被指数NDVI由于计算简单，较好的消除太阳高度角、地形、天气等辐射变化的影响，已经成为植被生长状态及植被覆盖度最佳指示因子；土壤调节植被指数作为NDVI的后续，是基于NDVI和大量观测数据提出的，主要用于减少土壤背景的影响。具体包括如下步骤：
4.1)基于高空间时序数据进行归一化差异植被指数和土壤调节植被指数计算，得到植被指数高空间分辨率时序数据集；所述的两种植被指数的计算公式如下：
归一化差异植被指数NDVI：
NDVI＝(Nir-Red)/(Nir+Red)
土壤调节植被指数SAVI：
SAVI＝(Nir-Red)*(1+L)/(Nir+Red+L)
式中，Nir为近红外波段地表反射率值，Red为红波段地表反射率值，L是随着植被密度变化的参数，取值范围是从0～1，当植被覆盖度很高时为0，很低时为1，通常取0.5时SAVI消除土壤反射率的效果最好。
4.2)将原始的时间序列数据进行归一化差异植被指数和土壤调节植被指数的计算，得到植被指数原始时序数据集；所述的两种植被指数的计算公式如下：
归一化差异植被指数NDVI：
NDVI＝(Nir-Red)/(Nir+Red)
土壤调节植被指数SAVI：
SAVI＝(Nir-Red)*(1+L)/(Nir+Red+L)
式中，Nir为近红外波段地表反射率值，Red为红波段地表反射率值，L是随着植被密度变化的参数，取值范围是从0～1，当植被覆盖度很高时为0，很低时为1，通常取0.5时SAVI消除土壤反射率的效果最好。
4.3)基于步骤4.1和4.2的结果，随机选取区域内1～5万个点，对锐化前与锐化后不同植被指数进行相关性分析；随机选点是通过ArcGIS的随机产生点的工具进行，并通过多值提取至点，进行不同时相的点之间的相关性大小；相关性大小通过相关性系数R2表示，相关系数越大，证明融合前后影像灰度值差异就越小，融合精度越准确，反之亦然。
5)结合野外采样各地物的样本点，构建分类样本库；综合考虑不同地物的物候曲线特征差异性较大，先将总体样本按其类型进行分层，然后在每层中随机抽取样本，并将抽取的样本按一定比例分为训练样本与验证样本。其中训练样本集与验证样本集的抽取比例优选为2:1。
6)利用随机森林的分类方法对作物类型进行识别，并对比原始时间序列数据提取结果与锐化后的时间序列提取结果。随机森林就是通过集成学习的思想将多棵树集成的一种算法，它的基本单元是决策树，而它的本质属于机器学习的一大分支—集成学习，随机森林最主要的优势在于不需要特征降维，适合容量较大的数据集并且具有良好的鲁棒性。主要包括以下步骤：
6.1)将训练数据集进行划分，抽取10％作为训练模型的测试数据，根据测试数据集的准确率，选择合适的训练参数，例如：随机树的数量，叶子节点的数量以及最大叶深度等；
6.2)将训练好的模型应用至整个研究区，完成作物类型分类，包括锐化后的高空间时序数据集和锐化前的时序数据集等；
6.3)对比原始时间序列数据提取结果与锐化后的时间序列提取结果。
7)基于步骤6)得到的不同作物类型的结果，将提取的所有类型的作物进行整合，形成耕地要素数据，将所有的非作物类型进行整合，形成非耕地数据并完成精度评价，主要包括步骤：
7.1)将不同的作物类型进行整合，形成耕地数据；将非作物类型整合，形成非耕地数据；
7.2)将验证集中的不同作物类型的样本点合并成一类，即耕地样本点，并进行精度评价；
7.3)如精度验证结果如果与期望精度不符，则返回步骤5)重新进行分层抽样并重复步骤6)与步骤7)，直至与期望精度相符，最终完成耕地类型识别
7.1)将不同的作物类型进行整合，形成耕地数据；将非作物类型整合，形成非耕地数据。
7.2)将验证集中的不同作物类型的样本点合并成一类，即耕地样本点，并进行精度评价。精度验证的指标选择总体准确率(Overall accuracy)、精确率(Precision)、召回率(Recall)以及F-score四个指标，其中计算公式如下：
accuracy＝(TP+TN)/(TP+FN+FP+TN)
precision＝TP/(TP+FP)
recall＝TP/(TP+FN)
F-score＝2*(precision*recall)/(precision+recall)
式中，TP为实际是正样本预测为正样本的样本数，FP为实际是负样本预测为正样本的样本数，TN为实际是负样本预测为负样本的样本数，FN为实际是正样本预测为负样本的样本数。
7.3)如精度验证结果如果与期望精度不符，则返回步骤5)重新进行分层抽样并重复步骤6与步骤7，直至与期望精度相符，最终完成耕地类型识别。
本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。
Claims (8)
1.一种基于图像锐化的时序数据耕地提取方法，其特征在于，包括如下步骤：
步骤1)，输入单时相高空间分辨率遥感影像全色波段并进行数据预处理；
步骤2)，输入L2A级别的多时相时间序列Sentinel-2影像并对其进行批量预处理；
步骤3)，对于步骤1)和步骤2)的结果，基于单时相高空间分辨率全色波段和多时相Sentinel-2多光谱波段进行Gram-Schmidt图像锐化操作，生成高空间分辨率时序数据；
所述步骤3)具体包括如下步骤：
3.1)使用多波段低分辨率影像对单波段高分辨率影像进行模拟，模拟的方法是根据权重Wi进行模拟，即模拟的全色波段影像灰度值
其中Bi为多光谱地空间分辨率影像第i波段灰度值，Wi为各波段的光谱响应函数对应的权重，k为影像的波段的数量；模拟的高分辨率波段影像在后面的处理中被称作Gram-Schmidt第一分量进行Gram-Schmidt变换；
3.2)利用模拟的高分辨率波段影像作为Gram-Schmidt变换第一分量来对模拟的高分辨波段影像和低分辨率波段影像进行Gram-Schmidt变换，在进行Gram-Schmidt变换时，第T个变换分量由前T-1个变换分量构造，即：
其中，GST是GS变换后产生的第T个正交分量，BT是原始的空间分辨率多光谱影像的第T个波段影像，uT是第T个原始低空间分辨率多光谱波段影像灰度值的均值，i,j分别表示低空间分辨率影像的行数和列数，
为原始低空间分辨率遥感影像第T个波段与GS变换后产生的第l个分量间的协方差，均值uT公式如下：
其中，C是影像的列数，R是影像行数，协方差
计算公式如下：
3.3)根据GS第一分量，即模拟波段的均值u和标准差σ对全色波段进行修改，修改公式如下：
P′＝(P×k1)+k2
k2＝uintensity-(k1×upan)
其中：P为原始全色波段影像的灰度值，eintensity和epan分别为GS第一分量和全色波段影像灰度值的标准差；
3.4)最后将修改后的全色波段作为第一分量，进行GS逆变换，输出N+1个波段，去除第一个波段，就是融合后的高分辨率多光谱影像结果，逆变换的计算公式如下：
步骤4)，对于步骤3)中的高空间分辨率时序数据，计算归一化差异植被指数和土壤调节植被指数，生成高空间分辨率植被指数时间序列数据集；并通过对比锐化前与锐化后的时序数据进行植被指数相关性分析；
步骤5)，根据步骤4)得到的高空间分辨率植被指数时间序列数据集，结合野外采样各地物的样本点，构建分类样本库并进行分层抽样；
步骤6)，利用步骤5)中构建的分类样本库，先利用随机森林算法进行模型训练，通过测试集的表现确定最优模型参数；将训练好的模型应用至全区域，实现全区域作物类型的识别，并对比原始时间序列数据提取结果与锐化后的时间序列提取结果；
步骤7)，利用步骤6)得到的作物类型的结果，将提取的所有类型的作物进行整合，形成耕地要素数据，将所有的非作物类型进行整合，形成非耕地数据并完成精度评价；最终生成耕地资源要素空间分布数据。
2.根据权利要求1所述的一种基于图像锐化的时序数据耕地提取方法，其特征在于：所述步骤1)中，高空间分辨率全色波段影像预处理包括辐射校正、正射校正、几何配准以及影像裁剪，在进行几何配准时，依托于高精度谷歌数据、天地图开源数据，手动选择若干控制点进行几何纠正；步骤2中所述预处理是通过欧空局SNAP平台进行格式转换输出与批量裁剪处理。
3.根据权利要求1所述的一种基于图像锐化的时序数据耕地提取方法，其特征在于：所述步骤2)中，输入时间序列L2A级别的Sentinel-2影像，并进行批量影像裁剪处理；其中，输入的Sentinel-2的影像经过严格的云量控制，即通过整景云量小于n1％或者研究区云量小于n2％来进行数据的筛选；所利用到的时序Sentinel-2是由蓝、绿、红以及近红外四个波段组成。
4.根据权利要求1所述的一种基于图像锐化的时序数据耕地提取方法，其特征在于：所述步骤4)具体包括如下步骤：
4.1)基于高空间时序数据进行归一化差异植被指数和土壤调节植被指数计算，得到植被指数高空间分辨率时序数据集；两种植被指数的计算公式如下：
归一化差异植被指数NDVI：
NDVI＝(Nir-Red)/(Nir+Red)
土壤调节植被指数SAVI：
SAVI＝(Nir-Red)*(1+L)/(Nir+Red+L)
式中，Nir为近红外波段地表反射率值，Red为红波段地表反射率值，L是随着植被密度变化的参数，取值范围是从0～1；
4.2)将原始的时间序列数据进行归一化差异植被指数和土壤调节植被指数的计算，得到植被指数原始时序数据集；
4.3)随机选取区域内若干个点，采用相关性系数R2对锐化前与锐化后不同植被指数进行相关性分析；随机选点是通过ArcGIS的随机产生点的工具进行，并通过多值提取至点，进行不同时相的点之间的相关性大小。
5.根据权利要求1所述的一种基于图像锐化的时序数据耕地提取方法，其特征在于：所述步骤5)是通过野外采集的各土地利用类型的样本点，建立该区域作物类型样本库，按照一定的比例进行分层抽样，获取训练数据集与验证数据集。
6.根据权利要求1所述的一种基于图像锐化的时序数据耕地提取方法，其特征在于：所述步骤6是通过随机森林的算法，实现不同作物类型的识别，并对比原始时间序列数据提取结果与锐化后的时间序列提取结果，包含以下步骤：
6.1)将训练数据集进行划分，抽取n3％作为训练模型的测试数据，根据测试数据集的准确率，选择合适的训练参数，包括随机树的数量，叶子节点的数量以及最大叶深度；
6.2)将训练好的模型应用至整个研究区，完成作物类型分类，包括锐化后的高空间时序数据集和锐化前的时序数据集；
6.3)对比原始时间序列数据提取结果与锐化后的时间序列提取结果。
7.根据权利要求1所述的一种基于图像锐化的时序数据耕地提取方法，其特征在于：所述步骤7具体包含以下步骤：
7.1)将不同的作物类型进行整合，形成耕地数据；将非作物类型整合，形成非耕地数据；
7.2)将验证集中的不同作物类型的样本点合并成一类，即耕地样本点，并进行精度评价；
7.3)若精度验证结果与期望精度不符，则返回步骤5)重新进行分层抽样并重复步骤6)与步骤7)，直至与期望精度相符，最终完成耕地类型识别。
8.根据权利要求7所述的一种基于图像锐化的时序数据耕地提取方法，其特征在于：所述步骤7.3)中，精度验证的指标选择总体准确率、精确率、召回率以及F-score四个指标。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110421926.8A CN113205014B (zh) | 2021-04-20 | 2021-04-20 | 一种基于图像锐化的时序数据耕地提取方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110421926.8A CN113205014B (zh) | 2021-04-20 | 2021-04-20 | 一种基于图像锐化的时序数据耕地提取方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN113205014A CN113205014A (zh) | 2021-08-03 |
| CN113205014B true CN113205014B (zh) | 2022-06-07 |
Family
ID=77027293
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202110421926.8A Active CN113205014B (zh) | 2021-04-20 | 2021-04-20 | 一种基于图像锐化的时序数据耕地提取方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN113205014B (zh) |
Families Citing this family (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113762226B (zh) * | 2021-11-09 | 2022-01-07 | 成都理工大学 | 基于高光谱分辨率调整提升树种识别精度的方法及系统 |
| CN115457356A (zh) * | 2022-08-16 | 2022-12-09 | 湖北省交通规划设计院股份有限公司 | 面向地质勘察的遥感图像融合方法、装置、设备及介质 |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105447494A (zh) * | 2015-12-01 | 2016-03-30 | 二十一世纪空间技术应用股份有限公司 | 基于多源多时相遥感影像数据的菜田监测方法 |
| CN109389049A (zh) * | 2018-09-19 | 2019-02-26 | 中国科学院东北地理与农业生态研究所 | 基于多时相sar数据与多光谱数据的作物遥感分类方法 |
| CN111598019A (zh) * | 2020-05-19 | 2020-08-28 | 华中农业大学 | 基于多源遥感数据的作物类型与种植模式识别方法 |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US11068737B2 (en) * | 2018-03-30 | 2021-07-20 | Regents Of The University Of Minnesota | Predicting land covers from satellite images using temporal and spatial contexts |
-
2021
- 2021-04-20 CN CN202110421926.8A patent/CN113205014B/zh active Active
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105447494A (zh) * | 2015-12-01 | 2016-03-30 | 二十一世纪空间技术应用股份有限公司 | 基于多源多时相遥感影像数据的菜田监测方法 |
| CN109389049A (zh) * | 2018-09-19 | 2019-02-26 | 中国科学院东北地理与农业生态研究所 | 基于多时相sar数据与多光谱数据的作物遥感分类方法 |
| CN111598019A (zh) * | 2020-05-19 | 2020-08-28 | 华中农业大学 | 基于多源遥感数据的作物类型与种植模式识别方法 |
Non-Patent Citations (1)
| Title |
|---|
| 基于随机森林模型的陆地卫星-8遥感影像森林植被分类;张晓羽等;《东北林业大学学报》;20160625(第06期);第53-57页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN113205014A (zh) | 2021-08-03 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN111598019B (zh) | 基于多源遥感数据的作物类型与种植模式识别方法 | |
| CN106780091B (zh) | 基于植被指数时空间统计特征的农业灾害信息遥感提取方法 | |
| Xiao et al. | Use of general regression neural networks for generating the GLASS leaf area index product from time-series MODIS surface reflectance | |
| CN105740759B (zh) | 基于多时相数据中特征提取的中稻信息决策树分类方法 | |
| CN107527014A (zh) | 县级作物种植面积遥感统计抽样调查方案设计方法 | |
| CN108458978B (zh) | 基于敏感波段和波段组合最优的树种多光谱遥感识别方法 | |
| CN112749627A (zh) | 一种基于多源遥感影像的烟草动态监测的方法和装置 | |
| CN113205014B (zh) | 一种基于图像锐化的时序数据耕地提取方法 | |
| Zhang et al. | Spatial domain bridge transfer: An automated paddy rice mapping method with no training data required and decreased image inputs for the large cloudy area | |
| Prasetyo et al. | Rice productivity prediction model design based on linear regression of spectral value using NDVI and LSWI combination on landsat-8 imagery | |
| Yu et al. | Deep convolutional neural networks for estimating maize above-ground biomass using multi-source UAV images: A comparison with traditional machine learning algorithms | |
| CN112800827A (zh) | 高光谱图像分类实验方法 | |
| Koranteng et al. | Remote sensing study of land use/cover change in West Africa | |
| Aguirre-Salado et al. | Modelling site selection for tree plantation establishment under different decision scenarios | |
| Liu et al. | Evaluating how lodging affects maize yield estimation based on UAV observations | |
| Ayub et al. | Wheat Crop Field and Yield Prediction using Remote Sensing and Machine Learning | |
| CN115830464A (zh) | 基于多源数据的高原山地农业大棚自动提取方法 | |
| WO2023131949A1 (en) | A versatile crop yield estimator | |
| CN115205688A (zh) | 茶树种植区域提取方法及系统 | |
| CN114611699A (zh) | 土壤水分降尺度方法、装置、电子设备及存储介质 | |
| Singh et al. | Assessment of digital image classification algorithms for forest and land-use classification in the eastern Himalayas using the IRS LISS III sensor | |
| Bao et al. | A fine digital soil mapping by integrating remote sensing-based process model and deep learning method in Northeast China | |
| do Amaral et al. | Geoprocessing Applied to Crop Management | |
| Niu et al. | Extraction of Winter Wheat Planting Area and Monitoring of Growth During key Fertility Periods Based on Time Series Sentinel-2 Images | |
| Driese et al. | Rule‐based Integration of Remotely‐sensed Data and GIS for Land Cover Mapping in NE Costa Rica |
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 |