CN112861810B - 一种基于时序遥感观测数据的人工林种植时间自动检测方法 - Google Patents
一种基于时序遥感观测数据的人工林种植时间自动检测方法 Download PDFInfo
- Publication number
- CN112861810B CN112861810B CN202110309180.1A CN202110309180A CN112861810B CN 112861810 B CN112861810 B CN 112861810B CN 202110309180 A CN202110309180 A CN 202110309180A CN 112861810 B CN112861810 B CN 112861810B
- Authority
- CN
- China
- Prior art keywords
- time
- index
- planting
- remote sensing
- artificial forest
- 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
- 238000001514 detection method Methods 0.000 title claims abstract description 40
- 238000000034 method Methods 0.000 claims abstract description 30
- 238000009499 grossing Methods 0.000 claims abstract description 17
- 238000007781 pre-processing Methods 0.000 claims abstract description 9
- 230000003595 spectral effect Effects 0.000 claims abstract description 7
- 238000012545 processing Methods 0.000 claims abstract description 5
- 230000035772 mutation Effects 0.000 claims description 70
- 230000008859 change Effects 0.000 claims description 18
- 238000000605 extraction Methods 0.000 claims description 17
- 239000002689 soil Substances 0.000 claims description 17
- 238000012549 training Methods 0.000 claims description 12
- 238000012937 correction Methods 0.000 claims description 11
- 238000002310 reflectometry Methods 0.000 claims description 11
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 claims description 11
- 238000004364 calculation method Methods 0.000 claims description 10
- 239000004576 sand Substances 0.000 claims description 10
- 238000007637 random forest analysis Methods 0.000 claims description 8
- 208000027066 STING-associated vasculopathy with onset in infancy Diseases 0.000 claims description 7
- 230000001404 mediated effect Effects 0.000 claims description 5
- 238000001228 spectrum Methods 0.000 claims description 4
- 230000015572 biosynthetic process Effects 0.000 claims description 3
- 230000003750 conditioning effect Effects 0.000 claims description 3
- 230000006872 improvement Effects 0.000 claims description 3
- 238000012163 sequencing technique Methods 0.000 claims description 3
- 238000003786 synthesis reaction Methods 0.000 claims description 3
- 238000012795 verification Methods 0.000 claims description 3
- 230000002159 abnormal effect Effects 0.000 abstract description 10
- 238000004422 calculation algorithm Methods 0.000 abstract description 9
- 238000011160 research Methods 0.000 description 8
- 230000008569 process Effects 0.000 description 6
- 238000010586 diagram Methods 0.000 description 5
- 238000007635 classification algorithm Methods 0.000 description 4
- 238000005516 engineering process Methods 0.000 description 4
- 239000000284 extract Substances 0.000 description 4
- 238000003066 decision tree Methods 0.000 description 3
- 241000238631 Hexapoda Species 0.000 description 2
- 241000607479 Yersinia pestis Species 0.000 description 2
- 238000004458 analytical method Methods 0.000 description 2
- 201000010099 disease Diseases 0.000 description 2
- 208000037265 diseases, disorders, signs and symptoms Diseases 0.000 description 2
- 238000011156 evaluation Methods 0.000 description 2
- 230000007774 longterm Effects 0.000 description 2
- 238000012216 screening Methods 0.000 description 2
- 238000012360 testing method Methods 0.000 description 2
- 239000002028 Biomass Substances 0.000 description 1
- 238000012300 Sequence Analysis Methods 0.000 description 1
- 238000010521 absorption reaction Methods 0.000 description 1
- 238000009825 accumulation Methods 0.000 description 1
- 230000006978 adaptation Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 229930002875 chlorophyll Natural products 0.000 description 1
- 235000019804 chlorophyll Nutrition 0.000 description 1
- 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 1
- 238000013145 classification model Methods 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 125000004122 cyclic group Chemical group 0.000 description 1
- 230000003247 decreasing effect Effects 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000001914 filtration Methods 0.000 description 1
- 238000003306 harvesting Methods 0.000 description 1
- 238000010801 machine learning Methods 0.000 description 1
- 238000013507 mapping Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000005457 optimization Methods 0.000 description 1
- 238000000638 solvent extraction Methods 0.000 description 1
- 238000009331 sowing 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
- 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
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
- G06N3/00—Computing arrangements based on biological models
- G06N3/004—Artificial life, i.e. computing arrangements simulating life
- G06N3/006—Artificial life, i.e. computing arrangements simulating life based on simulated virtual individual or collective life forms, e.g. social simulations or particle swarm optimisation [PSO]
-
- 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/10—Office automation; Time management
- G06Q10/109—Time management, e.g. calendars, reminders, meetings or time accounting
-
- 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)
- Theoretical Computer Science (AREA)
- Physics & Mathematics (AREA)
- Business, Economics & Management (AREA)
- Data Mining & Analysis (AREA)
- General Physics & Mathematics (AREA)
- Human Resources & Organizations (AREA)
- Life Sciences & Earth Sciences (AREA)
- General Health & Medical Sciences (AREA)
- Strategic Management (AREA)
- Evolutionary Computation (AREA)
- Artificial Intelligence (AREA)
- General Engineering & Computer Science (AREA)
- Health & Medical Sciences (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Tourism & Hospitality (AREA)
- Entrepreneurship & Innovation (AREA)
- General Business, Economics & Management (AREA)
- Bioinformatics & Computational Biology (AREA)
- Marketing (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Economics (AREA)
- Evolutionary Biology (AREA)
- Mathematical Physics (AREA)
- Computational Linguistics (AREA)
- Marine Sciences & Fisheries (AREA)
- Mining & Mineral Resources (AREA)
- Software Systems (AREA)
- Primary Health Care (AREA)
- Animal Husbandry (AREA)
- Molecular Biology (AREA)
- Agronomy & Crop Science (AREA)
- Computing Systems (AREA)
- Multimedia (AREA)
- Biophysics (AREA)
- Biomedical Technology (AREA)
- Operations Research (AREA)
- Quality & Reliability (AREA)
- Image Processing (AREA)
Abstract
本发明公开了一种基于时序遥感观测数据的人工林种植时间自动检测方法，包括1)遥感影像获取；2)影像预处理；3)人工林种植范围提取；4)人工林种植时间信息提取.本发明依据种植人工林时植被指数变化特征，建立综合考虑植被光谱变化特征且通过时序数据处理与子空间变化探测的人工林种植时间的方法，提升人工林种植时间探测的科学性和可靠性，该方法结合平滑、条件设置等操作来减少异常值和外界因素的影响，提升算法针对干扰的鲁棒性。
Description
技术领域
本发明属于遥感数据检测领域，具体涉及一种基于时序遥感观测数据的人工林种植时间自动检测方法的。
背景技术
人工林是指采用人工播种、栽植或扦插等方法和技术措施营造培育而成的森林，通常使用培育好的树种幼苗进行移植，按用途分为人工用材林、人工薪炭林、人工经济林、人工防护林等，种植人工林在治理水土流失、防治沙漠化、木材采伐良性循环等方面具有重要意义。人工林种植区域在植被对环境影响的研究中需要考核植被生长前后的周边环境变化，因此，精确的人工林种植时空分布信息是人工林种植情况多维表达的基础，也是人工林对环境影响研究的关键因素。
我国的人工林种植面积居于世界首位，但种植人工林的记录手段比较缺乏，大部分仅依靠“行政单位+面积+年份”的记录模式。目前，随着遥感技术长时间的发展已经积累了大量的时序影像资料，然而，利用遥感影像对大面积人工林种植时空分布的探测方法研究相对较少，目前已有的技术通过对植被的生物量进行年龄估计，或针对种植人工林时的时序植被指数特征的研究需要先验数据且易受环境等因素的干扰，不具有适用性且精度较低。
传统方法中针对经济型人工林重构时间序列归一化植被指数NDVI；通过高分辨率影像分类识别小区域尺度人工林获取其空间分布；再利用初始种植时间和初始砍伐时间来获取其间差值，通过差值来迭代获取每一次砍伐后种植的时间以获取时间分布；融合空间分布与时间分布来获得人工林的时空分布。其缺点在于：①仅适用于经济林的种植时间判断，不适用于所有类型的人工林；②判断方法需要砍伐间隔的先验数据；③小区域砍伐的间隔可能有不同，不具有科学性。
或者通过分类算法获取人工林的空间分布，利用已有的森林扰动算法进行时序植被指数的突变点探测，以探测人工林区域的时序植被指数突变时间；其缺点在于：①需要对算法设计一系列的控制参数和滤波过程来防止过拟合现象，同时还需要对不同地区的植被生长情况进行人为的参数调整，鲁棒性较差；②森林扰动算法探测所有的扰动点，不单纯地针对人工林的种植时间时序植被指数特征，因此需要一种基于时序遥感观测数据的人工林种植时间自动检测方法。
发明内容
为了解决上述问题，本发明提出一种基于时序遥感观测数据的人工林种植时间自动检测方法，针对现有人工林种植时间的反演技术的局限性，提出了一种新的基于时序遥感数据的突变探测方法来对人工林种植时间进行精确自动提取。
本发明包括以下步骤：
(1)获取遥感影像；
(2)对所述遥感影像进行预处理；
(3)对所述遥感影像进行人工林种植范围提取；
(4)在所述遥感影像的人工林种植范围进一步进行人工林种植时间信息提取；
(5)在所述时间信息提取中选择突变点并优化所述突变点的精度以输出种植时间点。
进一步地，所述人工林种植范围提取包括：
①选择训练样本及配置分类指数；
利用经过预处理后的遥感影像进行假彩色合成以选择用于人工林分布提取的样本，除了样本自身光谱信息外，选择以下指数建立样本特征空间：归一化植被指数NDVI(Normalized difference vegetation index)、增强植被指数EVI(Enhanced vegetationindex)、土壤调解植被指数SAVI(Soil adjusted vegetation index)、归一化沙地指数NDSI(Normalized difference sand index)、新型水体指数NWI(New water index)、归一化建筑指数NDBI(Normalized difference building index)。
②训练分类器并分类最新影像；将所述样本点输入随机森林分类器进行分类器训练，根据所选样本的在遥感传感器反射光谱以及特征指数中的不同特征取值空间来进行类别判断，建立分类规则，在不同的表征指数中具有特定的取值空间，其中，输入指标包括：蓝光波段-Blue、绿光波段-Green、红光波段-Red；近红外波段-NIR以及两个短波红外波段-SWIR1、SWIR2的反射率；各项指数包括归一化植被指数NDVI、增强植被指数EVI、土壤调解植被指数SAVI、归一化沙地指数NDSI、新型水体指数NWI、归一化建筑指数NDBI。
③对所述最新影像进行精度验证与改进，将所选样本集合按照7:3的比例进行随机划分，70％用于分类获得一个分类结果，30％用于验证该分类结果的精度。
进一步地，所述人工林种植时间的信息提取包括以下步骤：
a通过多卫星平台地表反射率交叉定标，选择沙地、裸土以及植被典型地物区域作为交叉标定区域；
b窗口平滑时序影像并提取时序植被指数，以平滑目标影像作为时间窗口中心，联合时间窗口内的影像数据进行平均处理作为平滑目标影像的值，通过窗口平滑后的时序影像集提取一种植被指数波段作为突变探测的实施目标，可选植被指数包含：差值植被指数DVI、比值植被指数RVI、归一化植被指数NDVI、增强植被指数EVI、土壤调节植被指数SAVI等；
c时序曲线子空间划分与判断指标计算，在获得平滑处理的长时序植被指数曲线的基础上，将长时序曲线划分成不同的子空间，对各时间点前后子空间内的时序植被指数分别进行线性拟合，获取前后子空间各自的拟合斜率，采用前后子空间的斜率差作为突变探测的判断指标；
d突变点确定,对所述斜率差的最大值作为该突变点发生时间。
进一步地，优化所述突变点的精度的方法包括查找不同平滑窗口尺度对应的适宜子空间尺度以输出种植时间，包括以下步骤：
(1)在阈值空间内对平滑窗口尺度及子空间尺度进行循环设置，所述阈值空间范围为1-9；
(2)若根据时序突变指标检测得到人工林植被指数合理突变点，则判定为该平滑窗口尺度对应的适宜子空间尺度，输出种植时间；
(3)若判定结果不合理，则判定为非适宜子空间尺度，并对子空间尺度进行扩大，再根据突变指标序列计算结果确定是否检测合理突变点；如果检测成功，停止调整，输出种植时间，否则，继续调整子空间尺度并进行类似判断；
(4)通过调整平滑窗口仍不能找到合理的突变点，则根据多子空间尺度和多平滑窗口大小下时序突变指标排序结果，选择最大值发生时间最为频繁点为突变点，确定种植时间。
进一步地，所述时序遥感影像预处理包括辐射定标、大气校正、几何校正和正射校正。
进一步地，所述平滑目标影像的公式为
其中，x0为平滑中心数据，n为窗口大小，平滑目标影像两侧影像数量为w。
本发明的有益效果在于：
1)基于遥感影像光谱信息、光谱指数信息等实现人工林种植范围自动提取；
2)克服传统时序突变探测易受干扰的缺点，提升人工林种植时间探测的科学性和精度；
3)适用于防护林以及经济林等不同类型的人工林，排除对先验数据的依赖，能够实现自动检测；
4)能够精确探测人工林种植时间。
附图说明
图1本发明的技术流程图；
图2前后子空间拟合线对比示意图；
图3突变点前后Sdiff对比示意图
具体实施方式
以下由特定的具体实施例以及附图说明本发明的实施方式，熟悉此技术的人士可由本说明书所揭露的内容轻易地了解本发明的其他优点及功效。
如图1所示本发明的具体流程如下：
1)遥感数据获取
本发明选用美国Landsat数据为基础遥感影像数据集，其主要特点是目前已具有较长时序观测，大部分地区可以获得1980以后至今的遥感数据，而且其相对较高的空间分辨率能够相对准确地表征地表人工林与其他地类之间的差异，为了后续针对人工林开展相关分析，本发明主要采用生长季的遥感影像数据。
2)遥感影像预处理
在获取生长季的时序遥感影像之后，需要对时序遥感影像集进行必要的预处理，包含辐射定标、大气校正、几何校正、正射校正等，辐射定标是为了减少传感器本身的误差，确保获得准确的辐射值；大气校正是为了减少由于大气散射、吸收、反射等引起的误差；几何校正减少由于系统因素和非系统因素引起的几何畸变；正射校正减少由于地形起伏和传感器方位角引起的畸变，在开展基础影像预处理的基础上，针对生长季遥感影像可能有云雾覆盖的特征，首先需采用云-阴影识别方法，结合地表反射率数据的质量文件，获取无云时序遥感影像数据集，去云操作主要为：通过云-阴影通用算法(CFmask)为每个像素计算云、云阴影分数，以判断每个像素的云、云阴影的不同分布情况，该算法主要通过分类进行云、云阴影的提取，具体为标记云、云阴影训练样本来训练分类器，再通过决策树分类器对影像进行分类，迭代估计云高后投影到地面来创建云阴影蒙版，在使用CFmask算法识别每年生长季影像云-阴影条件的基础上，根据每个像素分配云、云阴影分数，在每个像元位置筛选云、云阴影指数较低的像素，通过这些像素来合成每年生长季的无云影像，基于上述算法循环操作便可以得到每年生长季的无云影像，构成时序无云遥感数据集。
3)人工林种植范围信息提取
开展人工林种植时间信息提取，其重要前提条件是要准确获得人工林的空间分布信息，本发明选择以上时序影像，基于最新时序影像开展人工林种植范围信息提取，计算机自动分类算法是实现这一过程的重要方式，考虑到随机森林监督分类算法能够利用很多种资料产生高准确度的分类器，依靠多种参数对不同类型样本点建立多颗分类决策树形成随机森林分类器，再对分类影像中所有点的各类参数依靠分类器中的决策树进行分类。
因此，本发明在人工种植范围信息提取部分选择随机森林分类器作为分类算法，开展人工林种植范围提取。主要步骤包括：
①选择训练样本及配置分类指数
利用经过预处理后的遥感影像进行假彩色合成以选择用于人工林分布提取的样本，人工林提取精度一定程度取决于选择样本的准确性、样本数量的充足性以及各类样本之间数量的均衡性，因此，在人工林样本选择的同时，同步提取各类非人工林地物样本。
与此同时，为了基于样本数据开展后续人工林空间分布信息的提取，除了样本自身光谱信息外，选择以下指数建立样本特征空间：归一化植被指数NDVI(Normalizeddifference vegetation index)、增强植被指数EVI(Enhanced vegetation index)、土壤调解植被指数SAVI(Soil adjust vegetation index)、归一化沙地指数NDSI(Normalizeddifference sand index)、新型水体指数NWI(New water index)、归一化建筑指数NDBI(Normalized difference building index)。
在最新影像中计算上述指数作为分类器的分类指标，计算公式如下表：
表1典型光谱指数计算公式，其中，NIR为近红外波段反射率；Red为红光反射率；Blue为蓝光反射率；SWIR1与SWIR2分别是两个宽度范围的短波红外反射率；Ultra-blue为深蓝色和紫色反射率。
②训练分类器并分类最新影像
将所选的样本点输入随机森林分类器进行训练，根据所选样本的在反射光谱以及特征指数中的不同特征取值空间来进行类别判断，建立分类规则，不同的地物在不同的波段中具有不同反射率大小，在不同的表征指数中具有特定的取值空间。
其中，输入指标包括：蓝光波段-Blue、绿光波段-Green、红光波段-Red；近红外波段-NIR以及两个短波红外波段-SWIR1、SWIR2的反射率；各项指数包括归一化植被指数NDVI、增强植被指数EVI、土壤调解植被指数SAVI、归一化沙地指数NDSI、新型水体指数NWI、归一化建筑指数NDBI。
采用上述训练样本及指标构建的分类模型来分类最新影像，以获得整个区域各类地物空间分布并提取相应的人工林空间分布区域。
③精度验证与改进
为了检验分类的精度，将所选样本集合按照7:3的比例进行随机划分，70％用于分类获得一个分类结果，30％用于验证该分类结果的精度，若分类结果的精度评估指数未达到精度要求，则返回训练样本选择进行训练样本的检验和选择，以优化样本的准确性、充足性、均衡性，优化样本选择之后再进行随机森林分类器分类，循环处理至满足精度要求，从分类结果中提取人工林种植范围。
4)人工林种植时间信息提取
在人工林的种植范围内进行突变探测获取人工林种植时间，首先需要了解人工林种植时的植被指数特征。防护林多在沙地、荒漠、稀疏草地等植被覆盖度较低的区域移植培育好的树种幼苗，种植前植被盖度较小种植后植被盖度增大，植被指数呈现平缓后突增的过程；而经济林、薪炭林等是包括开荒初期以及砍伐后进行种植幼苗，砍伐后植被盖度急剧下降种植后植被盖度瞬间提升，植被指数呈现先下降再升高的变化。所以，种植各类人工林都将导致植被指数的突增，人工林种植区植被指数的突增时间即为其种植时间。
为准确提取其突变时间，本发明以长时序的植被指数为数据基础，针对其时序变化曲线提出了相应的突变时间探测方法，提取人工林种植时间，主要步骤包括：①多卫星平台地表反射率交叉定标；②窗口平滑时序影像并提取时序植被指数；③时序曲线子空间划分与判断指标计算；④突变点确定；⑤适宜平滑窗口尺度及子空间尺度判断。
①多卫星平台地表反射率交叉定标
为了降低不同传感器之间的系统误差，对时序遥感影像进行交叉定标，选择沙地、裸土、植被等典型地物区域作为交叉标定区域；由于Landsat 7的运行周期与Landsat 5及Landsat 8存在交集，故选择Landsat 7-ETM+传感器来标定Landsat 5-TM传感器以及Landsat 8-OLI传感器；将ETM+分别与TM、OLI在交叉标定区域同一时间的各波段反射率分别进行线性拟合，获得拟合直线：
R标定＝a*R待标定+b (1)
其中：a、b分别为交叉标定系数，R标定为标定传感器在标定区域的反射率，R待标定为待标定传感器在标定区域的反射率；利用各波段的交叉标定系数a、b来对TM、OLI传感器的各波段反射率通过式1进行标定计算，得到标定后的R标定。
②窗口平滑时序影像并提取时序植被指数
受原始遥感数据自身精度、云雾覆盖、地形起伏等因素的影响，各像元点提取的长时序植被指数数据中必然会存在一定的异常值，这些数据的直接引入会对突变探测造成突变特征混淆，而可能存在的数据遗漏点也会影响后续的时序分析。
为了减弱异常值和缺失值的影响，本发明首先对时序影像进行窗口平滑，种植人工林时的植被指数突增是一个长期的变化，而异常值常出现于特定时刻，因此选择一定大小的时间窗口平滑能够压缩异常信号，具体平滑操作如下：以平滑影像作为时间窗口中心，联合时间窗口内的影像数据进行平均处理作为目标窗口影像的值。
具体操作可参考公式2：
其中，x0为平滑中心数据，n为窗口大小，平滑目标影像两侧影像数量为w。
由于可探测时间序列初试时间之前及当前时间之后缺少可用数据，为了满足初始年份影像以及当前年份影像作为窗口中心点进行平滑的需求，将初始影像复制w份作为初始时间前w年的数据，将当前影像复制w张作为当前时间后w年的数据，复制影像只用作平滑使用不参与突变探测。
通过窗口平滑后的时序影像集提取一种植被指数波段作为突变探测的实施目标，可选植被指数包含：差值植被指数DVI、比值植被指数RVI、归一化植被指数NDVI、增强植被指数EVI、土壤调节植被指数SAVI等，本发明选取归一化植被指数NDVI作为时序突变探测的实施目标，NDVI对低地植被覆盖区域植被敏感的特性适合用于人工林种植初期的突变探测。
③时序曲线子空间划分与判断指标计算
在获得平滑处理的长时序植被指数曲线的基础上，将长时序曲线划分成不同的子空间，以带探测时间点Y为例，假设子空间的宽度为W，时间点Y-W至时间点Y为时序植被指数的前子空间，时间点Y时间至时间点Y+W为时序植被指数的后子空间。
根据此类划分标准，在时间序列上依次对每年的数据进行子空间划分。其中，为了满足研究时间序列初始和当前年份的子空间拟合需求，将初始植被指数和当前植被指数复制W并向外依次填充，但这些复制数据仅作拟合使用不作为突变探测点。
在子空间划分的基础上，对各时间点前后子空间内的时序植被指数进行线性拟合，获取前后子空间的拟合斜率Sbef和Saft，及对应的截距b1、b2。
如图2所示为，前后子空间拟合线示意图为了将探测点前后的变化趋势融合到一个指标之内，本发明提出采用后子空间拟合线斜率Saft与前子空间拟合线斜率Sbef差值为突变探测的判断指标，即前后子空间的斜率差Sdiff(式3)。
Sdiff＝Saft-Sbef (3)
Sdiff可以通过其正负性来表示探测点的变化趋势(正值则表示增加的变化趋势，负值则表示下降的变化趋势)；其绝对值可以表示变化程度的高低(绝对值越大，变化程度越大)。依据以上公式，每个像元点在时间序列的每个时间点都可以得到对应的Sdiff，获得突变指标的时间序列。
④突变点确定
为了去除外界因素对时序植被指数的影响，需要对判断指标筛选进行条件设置。遥感影像获取的是地表的反射率，水体覆盖会导致植被指数降低至负数，植被指数呈现下降至一个负值后再上升的规律；病虫害、干旱等导致植被叶绿素含量下降产生异常值，时序植被指数在较高指数附近出现小幅的凹陷。两者皆会产生较高的判断指标Sdiff数值，与真实的种植时间产生混淆。为了消除水体产生的影响，在筛选判断指标Sdiff需设置突变点NDVI值大于0的条件；为了避免在生长过程中的病虫害、干旱等自然因素引起的植被指数曲线变化与人工林种植时植被指数曲线变化类似进而导致突变信息提取有误，本发明设置突变点NDVI值在研究区裸土NDVI值附近，具体数值大小可根据研究区裸土NDVI大小来定。
在以上条件限定的前提下，人工林时序植被指数一般呈现先平滑后增加的趋势，体现在前子空间内时序NDVI平缓，斜率趋于0；后子空间内时序NDVI上升，斜率为正值。因此，突增点处的Sdiff为一个正值减去一个趋于0的值，结果为一个正值。而在突变点之前，后子空间的斜率由于放缓导致Sdiff相对突变点较小；而在突变点之后，前子空间的斜率会提升导致Sdiff相对突变点也较小。Sdiff在突变点处达到最大值。因此，在这样的变化过程中，人工林种植时间即为该突变点发生时间，也就是筛选Sdiff的最大值。如图3所示为，突变点前后Sdiff对比示意图。
⑤种植时间判断优化
在基于上述方法开展种植时间判断的时候，可能由于部分时间植被干扰或者遥感信号不确定引起探测突变点存在一定的偏差，进而干扰了种植时间的探测精度。因此，为了消除异常值的干扰而不影响探测精度，找出准确的种植时间，本发明还提出了平滑窗口和子空间尺度优化方法。具体如下：
根据突变指标计算方法获得时序突变指标，并依据大小顺序提取突变指标的前3个最大值。如果最大值与其他两个数值差异显著，则认为该突变点即为种植时间点；如果最大值与其他两个数值差异不大，则说明三个时间存在一定相似性，其突变点判断有异常值的影响。在这种情况下，本发明对平滑窗口尺度和子空间尺度分别进行嵌套循环设定，查找不同平滑窗口尺度对应的适宜子空间尺度：
(1)平滑窗口大小初始设置为3，子空间尺度初始设定为3；
(2)若根据时序突变指标检测得到人工林植被指数合理突变点，则判定为该平滑窗口尺度对应的适宜子空间尺度，输出种植时间；
(3)若判定结果不合理，则判定为非适宜子空间尺度，并对子空间尺度进行扩大，再根据突变指标序列计算结果确定是否检测合理突变点；如果检测成功，停止调整，输出种植时间，否则，继续调整子空间尺度并进行类似判断；
(4)如果判定结果一直不合适，且子空间尺度调整至最大阈值5仍不能得到合理结果，则调整平滑窗口大小，平滑窗口大小增加2至5，子空间尺度仍从3开始；
(5)重复步骤2-4操作，在不能找到合适突变点的情况下，再次调整平滑窗口大小(窗口最大阈值为7)，并继续重复步骤2-4操作；
(6)如果步骤5-6通过调整平滑窗口仍不能找到合理的突变点，则根据多子空间尺度和多平滑窗口大小下时序突变指标排序结果，选择最大值发生时间最为频繁点为突变点，确定种植时间。
最后，根据人工林空间范围探测结果和种植时间反演结果，联合完成人工林精细制图，为后续分析人工林长势和评估提供数据支撑。
本发明提出了一种新的探测线性突变点的方法来探测人工林的种植时间。利用植被指数作为表征地面植被情况的指标，针对种植人工林时该区域植被指数突增的特征，设计出突变判断指数代表探测点前后的变化趋势，对时序植被指数进行突变探测以探测人工林种植时间。基于人工林光谱、指数等信息，采用机器学习算法建立人工林种植空间范围提取方法。提出了利用时序植被指数突增点作为人工林的种植时间，为表征人工林的种植时间提出了一个具体的指代对象。传统的突变探测依靠累积值、者单纯的阈值或回归等方式来进行突变判断指数的构造，极易受到异常值的影响。该方法使用判断指标Sdiff将探测点前后的趋势变化差异结合到一个指标之中，利用窗口平滑和子空间尺度的大小调整减小异常值和缺失值的影响，使突变判断指标的最值控制在拐点处。本发明只需要时序遥感影像作为研究数据，排除了对先验数据的依赖性。
上述实施例仅例示性说明本发明的原理及其功效，而非用于限制本发明。任何熟悉此技术的人士皆可在不违背本发明的精神及范畴下，对上述实施例进行修饰或改变。因此，举凡所属技术领域中具有通常知识者在未脱离本发明所揭示的精神与技术思想下所完成的一切等效修饰或改变，仍应由本发明的权利要求所涵盖。
Claims (5)
1.一种基于时序遥感观测数据的人工林种植时间自动检测方法，其特征在于通过时序遥感数据自动获取其种植时间，不需要其他先验知识，包括以下步骤
1)获取遥感影像；
2)对所述遥感影像进行预处理；
3)对所述遥感影像进行人工林种植范围提取；
4)在所述遥感影像的人工林种植范围进行人工林种植时间信息提取，所述人工林种植时间的信息提取包括以下步骤：
a通过多卫星平台地表反射率交叉定标，选择沙地、裸土以及植被典型地物区域作为交叉标定区域；
b窗口平滑时序影像并提取时序植被指数，以平滑目标影像作为时间窗口中心，联合时间窗口内的影像数据进行平均处理作为平滑目标影像的值，通过窗口平滑后的时序影像集提取一种植被指数波段作为突变探测的实施目标，可选植被指数包含：差值植被指数DVI、比值植被指数RVI、归一化植被指数NDVI、增强植被指数EVI、土壤调节植被指数SAVI；
c时序曲线子空间划分与判断指标计算，在获得平滑处理的长时序植被指数曲线的基础上，将长时序曲线划分成不同的子空间，对各时间点前后子空间内的时序植被指数分别进行线性拟合，获取前后子空间各自的拟合线斜率，采用前后子空间的斜率差作为突变探测的判断指标；
d突变点确定,对所述斜率差的最大值作为该突变点发生时间；
5)在所述时间信息提取中选择突变点并优化所述突变点的精度以输出种植时间点。
2.根据权利要求1所述一种基于时序遥感观测数据的人工林种植时间自动检测方法，其特征在于，所述人工林种植范围提取包括：
①选择训练样本及配置分类指数；
利用经过预处理后的遥感影像进行假彩色合成以选择用于人工林分布提取的样本，除了遥感卫星自身的光谱反射信息外，选择以下指数建立样本特征空间：归一化植被指数NDVI、增强植被指数EVI、土壤调解植被指数SAVI、归一化沙地指数NDSI、新型水体指数NWI、归一化建筑指数NDBI；
②训练分类器并分类最新影像；将样本点输入随机森林分类器进行训练，根据所选样本的在反射光谱以及特征指数中的不同特征取值空间来进行类别判断，建立分类规则，在不同的表征指数中具有特定的取值空间，其中，输入指标包括：蓝光波段-Blue、绿光波段-Green、红光波段-Red；近红外波段-NIR以及两个短波红外波段-SWIR1、SWIR2的反射率；各项指数包括归一化植被指数NDVI、增强植被指数EVI、土壤调解植被指数SAVI、归一化沙地指数NDSI、新型水体指数NWI、归一化建筑指数NDBI；
③对所述最新影像进行精度验证与改进，将所选样本集合按照7:3的比例进行随机划分，70％用于分类获得一个分类结果，30％用于验证该分类结果的精度。
3.根据权利要求1所述一种基于时序遥感观测数据的人工林种植时间自动检测方法，其特征在于，优化所述突变点的精度的方法包括查找不同平滑窗口尺度对应的适宜子空间尺度以输出种植时间，包括以下步骤：
(1)在阈值空间内对平滑窗口尺度及子空间尺度进行循环设置，所述阈值空间为1-9；
(2)若根据时序突变指标检测得到人工林植被指数合理突变点，则判定为该平滑窗口尺度对应的适宜子空间尺度，输出种植时间；
(3)若判定结果不合理，则判定为非适宜子空间尺度，并对子空间尺度进行扩大，再根据突变指标序列计算结果确定是否检测合理突变点；如果检测成功，停止调整，输出种植时间，否则，继续调整子空间尺度并进行类似判断；
(4)通过调整平滑窗口仍不能找到合理的突变点，则根据多子空间尺度和多平滑窗口大小下时序突变指标排序结果，选择最大值发生时间最为频繁点为突变点，确定种植时间。
4.根据权利要求1所述一种基于时序遥感观测数据的人工林种植时间自动检测方法，其特征在于，所述遥感影像预处理包括辐射定标、大气校正、几何校正和正射校正。
5.根据权利要求1所述一种基于时序遥感观测数据的人工林种植时间自动检测方法，其特征在于，所述平滑目标影像的公式为
其中，x0为平滑中心数据，n为窗口大小，平滑目标影像两侧影像数量为w。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110309180.1A CN112861810B (zh) | 2021-03-23 | 2021-03-23 | 一种基于时序遥感观测数据的人工林种植时间自动检测方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110309180.1A CN112861810B (zh) | 2021-03-23 | 2021-03-23 | 一种基于时序遥感观测数据的人工林种植时间自动检测方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN112861810A CN112861810A (zh) | 2021-05-28 |
| CN112861810B true CN112861810B (zh) | 2021-11-02 |
Family
ID=75992430
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202110309180.1A Active CN112861810B (zh) | 2021-03-23 | 2021-03-23 | 一种基于时序遥感观测数据的人工林种植时间自动检测方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN112861810B (zh) |
Families Citing this family (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115346120B (zh) * | 2022-08-16 | 2023-06-20 | 中国科学院空天信息创新研究院 | 一种草地地上生物量及其固碳量遥感估算方法 |
| CN115471760B (zh) * | 2022-10-31 | 2023-02-28 | 吉林高分遥感应用研究院有限公司 | 玉米长势监测方法、系统、电子设备及计算机存储介质 |
| CN116306819B (zh) * | 2023-03-22 | 2024-05-03 | 大连海事大学 | 基于光谱重构的高光谱交叉定标方法、装置及电子设备 |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2002360070A (ja) * | 2001-06-12 | 2002-12-17 | Kansai Electric Power Co Inc:The | 植物の活力度評価法 |
| CN109840516A (zh) * | 2019-03-06 | 2019-06-04 | 福州大学 | 一种基于时序遥感影像的水体变化自动识别方法 |
| CN111832506A (zh) * | 2020-07-20 | 2020-10-27 | 大同煤矿集团有限责任公司 | 一种基于长时序植被指数的重建植被遥感判别方法 |
-
2021
- 2021-03-23 CN CN202110309180.1A patent/CN112861810B/zh active Active
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2002360070A (ja) * | 2001-06-12 | 2002-12-17 | Kansai Electric Power Co Inc:The | 植物の活力度評価法 |
| CN109840516A (zh) * | 2019-03-06 | 2019-06-04 | 福州大学 | 一种基于时序遥感影像的水体变化自动识别方法 |
| CN111832506A (zh) * | 2020-07-20 | 2020-10-27 | 大同煤矿集团有限责任公司 | 一种基于长时序植被指数的重建植被遥感判别方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN112861810A (zh) | 2021-05-28 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN112861810B (zh) | 一种基于时序遥感观测数据的人工林种植时间自动检测方法 | |
| CN111461052A (zh) | 基于迁移学习的多个生育期小麦倒伏区域识别方法 | |
| Liu et al. | Estimating potato above-ground biomass by using integrated unmanned aerial system-based optical, structural, and textural canopy measurements | |
| CN105678281B (zh) | 基于光谱和纹理特征的地膜覆盖农田遥感监测方法 | |
| CN106372592B (zh) | 一种基于冬小麦面积指数的冬小麦种植面积计算方法 | |
| CN112800973B (zh) | 一种基于植被物候特征决策的互花米草提取方法 | |
| Jin et al. | High-throughput measurements of stem characteristics to estimate ear density and above-ground biomass | |
| CN111461053A (zh) | 基于迁移学习的多个生育期小麦倒伏区域识别系统 | |
| CN113029971B (zh) | 一种作物冠层氮素监测方法及系统 | |
| CN110398466A (zh) | 基于遥感反演的农作物生长状态监测方法 | |
| CN108710864B (zh) | 基于多维度识别及图像降噪处理的冬小麦遥感提取方法 | |
| CN116665073A (zh) | 一种基于多源数据的玉米产量遥感估测方法 | |
| CN112287886B (zh) | 基于高光谱图像融合图谱特征的小麦植株氮含量估测方法 | |
| Liu et al. | Estimating wheat fractional vegetation cover using a density peak k-means algorithm based on hyperspectral image data | |
| CN112711989A (zh) | 基于雷达遥感与光学遥感的玉米秸秆覆盖度估算方法 | |
| CN108764284B (zh) | 一种对松树病死木的高分辨率影像的分类去噪方法及系统 | |
| CN114387516A (zh) | 一种针对复杂地形环境下中小田块的单季稻sar识别方法 | |
| CN116543316B (zh) | 一种利用多时相高分辨率卫星影像识别稻田内草皮的方法 | |
| CN112215090A (zh) | 融合物候知识的遥感水稻制图方法及其应用 | |
| CN105678280A (zh) | 基于纹理特征的地膜覆盖农田遥感监测方法 | |
| CN116563721B (zh) | 基于分层分类思想的烟田提取方法 | |
| GASTELLU-ETCHEGORRY et al. | Computer-assisted land cover mapping with SPOT in Indonesia | |
| CN115063610A (zh) | 基于Sentinel-1、2影像的大豆种植区识别方法及其面积测算方法 | |
| Santana et al. | Counting of Coffee Trees Based on Convolutional Neural Network Applied to RGB Images Obtained by RPA. Sustainability 2023, 15, 820 | |
| CN118470441B (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 |