WO2021098471A1 - 一种基于形态模型法的大范围作物物候提取方法 - Google Patents

一种基于形态模型法的大范围作物物候提取方法 Download PDF

Info

Publication number
WO2021098471A1
WO2021098471A1 PCT/CN2020/124834 CN2020124834W WO2021098471A1 WO 2021098471 A1 WO2021098471 A1 WO 2021098471A1 CN 2020124834 W CN2020124834 W CN 2020124834W WO 2021098471 A1 WO2021098471 A1 WO 2021098471A1
Authority
WO
WIPO (PCT)
Prior art keywords
morphological model
curve
vegetation index
crop
pixel
Prior art date
Application number
PCT/CN2020/124834
Other languages
English (en)
French (fr)
Inventor
林涛
林智贤
江昊
应义斌
Original Assignee
浙江大学
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by 浙江大学 filed Critical 浙江大学
Priority to US17/778,014 priority Critical patent/US11734925B2/en
Publication of WO2021098471A1 publication Critical patent/WO2021098471A1/zh

Links

Images

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06VIMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V20/00Scenes; Scene-specific elements
    • G06V20/10Terrestrial scenes
    • G06V20/188Vegetation
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06QINFORMATION 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/00Administration; Management
    • G06Q10/06Resources, workflows, human or project management; Enterprise or organisation planning; Enterprise or organisation modelling
    • G06Q10/063Operations research, analysis or management
    • G06Q10/0639Performance analysis of employees; Performance analysis of enterprise or organisation operations
    • G06Q10/06393Score-carding, benchmarking or key performance indicator [KPI] analysis
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T11/002D [Two Dimensional] image generation
    • G06T11/20Drawing from basic elements, e.g. lines or circles
    • G06T11/203Drawing of straight lines or curves
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/0002Inspection of images, e.g. flaw detection
    • G06T7/0004Industrial image inspection
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/10Segmentation; Edge detection
    • G06T7/149Segmentation; Edge detection involving deformable models, e.g. active contour models
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06VIMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V10/00Arrangements for image or video recognition or understanding
    • G06V10/40Extraction of image or video features
    • G06V10/60Extraction of image or video features relating to illumination properties, e.g. using a reflectance or lighting model
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06VIMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V10/00Arrangements for image or video recognition or understanding
    • G06V10/70Arrangements for image or video recognition or understanding using pattern recognition or machine learning
    • G06V10/764Arrangements for image or video recognition or understanding using pattern recognition or machine learning using classification, e.g. of video objects
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06VIMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V10/00Arrangements for image or video recognition or understanding
    • G06V10/70Arrangements for image or video recognition or understanding using pattern recognition or machine learning
    • G06V10/77Processing 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/7715Feature extraction, e.g. by transforming the feature space, e.g. multi-dimensional scaling [MDS]; Mappings, e.g. subspace methods
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06VIMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V10/00Arrangements for image or video recognition or understanding
    • G06V10/70Arrangements for image or video recognition or understanding using pattern recognition or machine learning
    • G06V10/77Processing 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/776Validation; Performance evaluation
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06VIMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V20/00Scenes; Scene-specific elements
    • G06V20/10Terrestrial scenes
    • G06V20/13Satellite images
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/10Image acquisition modality
    • G06T2207/10032Satellite or aerial image; Remote sensing
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/20Special algorithmic details
    • G06T2207/20112Image segmentation details
    • G06T2207/20124Active shape model [ASM]
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/30Subject of image; Context of image processing
    • G06T2207/30181Earth observation
    • G06T2207/30188Vegetation; Agriculture

Definitions

  • the invention relates to the field of agricultural remote sensing, in particular to a large-scale crop phenology extraction method based on a morphological model method.
  • Agricultural production is the cornerstone of human society, and extensive and reliable agricultural information is an important foundation for ensuring regional and global food security.
  • the phenology of crops is reflected in the seasonal changes in their physiological characteristics and their physiological responses to environmental conditions along with growth, which can reflect the impact of climate change on agricultural ecosystems.
  • the phenological information of crops can provide important decision-making information for crop irrigation scheduling and fertilization management; accurate phenological period is also an important input information for accurate crop yield prediction models; in addition, phenological period is also an important basis for remote sensing of crop identification.
  • the methods of phenological monitoring mainly include field observation method, accumulated temperature model method and remote sensing monitoring.
  • the field observation method is relatively simple and easy to implement, and the results are more accurate.
  • the quality of the obtained data is difficult to guarantee, and it requires a lot of manpower and material resources.
  • the observation period is long and the area is small, making it difficult to apply to a large range.
  • Accumulated temperature model method requires higher input information, including crop sowing date, temperature and accumulated temperature information corresponding to different phenological periods, which limits its application range.
  • the remote sensing technology does not have the above shortcomings, it has the advantages of wide monitoring range, high time efficiency, objective accuracy and so on.
  • Threshold method which uses fixed or dynamic thresholds to estimate the phenology stage
  • Moving window method which uses detection timing The change of the vegetation index curve in the defined moving time window is used to determine the phenological period
  • the function fitting method is to use mathematical functions such as logistic function to fit the time series vegetation index curve of the crop to be extracted, and The phenological period is determined according to the characteristic points existing on the fitted curve.
  • these methods mainly focus on the local morphological feature points or feature segments of the vegetation index time series curve.
  • the present invention proposes a large-scale crop phenology extraction method based on morphological model method, which uses satellite remote sensing data to establish vegetation index time series curve, uses morphological model method to perform curve fitting, and combines ground observation data , Calculate the phenology extraction value according to the optimal fit scaling parameter, which solves the limitations of the current remote sensing crop phenology extraction application such as small scale, low accuracy, and poor versatility.
  • Step S1 Vegetation index time series curve acquisition: According to the classification data of the crop surface, select all the pixels where the crops to be extracted are located in the local geographic area, and calculate the vegetation index of each pixel based on the satellite remote sensing image data, and then obtain each The vegetation index time series curve of the pixel; and, according to the geographic coordinate data of the agricultural gas station, the multi-year vegetation index time series curve of each pixel where the agricultural gas station is located;
  • the crop surface classification data is the existing data, and the spatial distribution information of different crops can be obtained by performing crop identification on the satellite remote sensing data.
  • the pixel is the pixel in the satellite remote sensing image.
  • Step S2 Time series curve smooth fitting: Use the double logistic function fitting method to smoothly fit the vegetation index time series curve of each pixel where the crops to be extracted and the agricultural gas station are located to obtain the smoothed vegetation index time series curve;
  • Step S3 Establishment of morphological model: use the ground phenological observations recorded by the agricultural gas station to establish the morphological model reference point, and use the vegetation index time series curve of the pixel where the agricultural gas station is located to establish the corresponding morphological model reference curve, from the morphological model reference point and morphology
  • the model reference curve constitutes a morphological model.
  • Each agricultural gas station establishes a morphological model, and multiple stations establish multiple morphological models;
  • Ground phenological observations refer to the phenological period of ground crops, such as the phenological period of wheat.
  • Step S4 Morphological model fitting: According to the spatial relationship between each pixel in the local geographic area and the agricultural gas station (spatial relationship refers to the distance, similar or dissimilar geographic environment), select each pixel where the crop to be extracted is located Corresponding to the most suitable morphological model of agricultural gas station, fit the morphological model reference curve to the smooth vegetation index time series curve of the pixel, calculate the root mean square error (RMSE) between the two curves after fitting, and calculate the mean square The zoom parameter corresponding to the minimum root error is used as the best zoom parameter of the pixel, and the zoom parameter has the best fitting effect;
  • RMSE root mean square error
  • Step S5 Regional phenological phase extraction: transform the morphological model reference point according to the optimal scaling parameter to obtain the phenological extraction result of each pixel where the crop to be extracted is located.
  • the vegetation index used in the present invention is specifically the Wide Dynamic Range Vegetation Index (WDRVI).
  • WDRVI Wide Dynamic Range Vegetation Index
  • the wide dynamic range vegetation index adds a weighting coefficient on the basis of NDVI, which can improve the sensitivity of NDVI in areas with high vegetation coverage.
  • the calculation formula of the Wide Dynamic Range Vegetation Index (WDRVI) is shown in the following formula:
  • ⁇ NIR is the reflectivity in the near-infrared band
  • ⁇ RED is the reflectivity in the red band
  • is the weighting coefficient. Its value varies with different crops, generally 0.1-0.2.
  • the method of fitting the dual logistic function in step S2 is specifically as follows: Fit the vegetation index time series curve of the entire crop growth period to the dual logistic function curve.
  • the dual logistic function includes six parameters, VI 1 , VI m , S, A , S′ and A′, the formula after fitting is expressed as:
  • VI 1 represents the average vegetation index in the non-growth period (that is, the first half) of the vegetation index time series curve
  • VI m represents the peak value of the vegetation index time series curve
  • t represents the time
  • S and A respectively represent the vegetation index time series curve.
  • the time corresponding to the rising inflection point and the falling inflection point, S′ and A′ represent the curvature at the rising inflection point and the falling inflection point, respectively
  • VI(t) represents the vegetation index at time t.
  • the rising inflection point is the inflection point with the largest slope
  • the falling inflection point is the inflection point with the largest slope
  • the morphological model corresponds to the agricultural gas station one-to-one, representing the typical seasonal vegetation index response of the crop in a certain spatial range of the corresponding agricultural gas station.
  • the morphological model reference point uses the average of the phenological period of the agricultural gas station's multi-year ground observations.
  • the morphological model reference curve uses the mean value of the time series curve of the multi-year historical vegetation index of all pixels where the agricultural gas station is located.
  • the most suitable morphological model for each pixel where the crop to be extracted is located is selected according to the spatial relationship between the pixel and the agricultural gas station, and the spatial distance and geographical environment factors are comprehensively considered and obtained by weighted calculation processing.
  • the specific implementation can be The similarity is calculated by weighting, and the one with the highest similarity is regarded as the most appropriate.
  • step S4 the following formula is used to fit and transform the morphological model reference curve into the smooth vegetation index time series curve of the crop to be extracted.
  • the specific formula is as follows:
  • g(x) represents the morphological model reference curve after fitting transformation
  • h(x) represents the morphological model reference curve
  • xscale and yscale respectively represent the scaling parameters in the x-axis and y-axis directions
  • tshift represents the x-axis direction Offset
  • x 0 represents the reference point of the morphological model
  • RMSE root mean square error
  • x represents the day-by-day date
  • n represents the time interval of the curve
  • f(x) represents the smooth vegetation index time series curve of the crop to be extracted.
  • the optimal scaling parameter is determined, and the fitting effect is the best at this time.
  • step S5 the following formula is used to transform the reference point of the morphological model according to the optimal scaling parameters to obtain the phenological period to be extracted for each pixel of the crop to be extracted, and to integrate the extraction of the phenological period of all pixels in the local geographic area Value to get the temporal and spatial distribution of phenology in the area:
  • x 0 represents the phenological period at the reference point of the morphological model
  • X est is the phenological period of the crop to be extracted.
  • the size of the local geographic area is not limited, and includes at least one pixel.
  • the phenological phases to be extracted include the various phenological phases of the entire growth period of the crop, taking corn crops as an example, including but not limited to the emergence phase, the three-leaf phase, the green phase, the heading phase, and the mature phase.
  • the pixel is the smallest pixel unit in the satellite remote sensing image data. According to different satellite remote sensing products, the pixel size can be 30m ⁇ 30m, 250m ⁇ 250m, 500m ⁇ 500m, etc.
  • the method of the invention solves the limitation of the existing method for extracting phenology of crops.
  • the morphological model is used to represent the typical growth curve of the crop, and then the morphological model and the crop growth curve are fitted through a transformation method, and finally the obtained optimal scaling parameters are used to calculate the phenological period of the extract.
  • This method utilizes the macro characteristics of the vegetation index time series curve, which can reduce the influence of local fluctuations and noise of the curve, and has better phenological extraction accuracy; in addition, this method can extract the various phenological periods of crops at the same time without setting different settings separately.
  • the metric standard can be applied to a wide range of crop phenology extraction.
  • the method of the present invention utilizes the macro characteristics of the curve, can reduce the influence of local fluctuations and noise of the curve, and has better extraction accuracy; and can extract various phenological periods of crops at the same time without setting different metric standards separately. It is applied to extract phenology of large-scale crops.
  • Figure 1 is a local geographic area and agrogas station distribution map according to a specific embodiment of the present invention
  • FIG. 2 is a schematic diagram of a large-scale crop phenology extraction method based on a morphological model method according to a specific embodiment of the present invention
  • Fig. 3 is a schematic diagram of smooth fitting of dual logic function curves according to a specific embodiment of the present invention
  • Figure 4 is a schematic diagram of a morphological model establishment process according to a specific embodiment of the present invention.
  • Figure 5 is a schematic diagram of a morphological model according to a specific embodiment of the present invention.
  • Figure 6 is a schematic diagram of a morphological model transformation process according to a specific embodiment of the present invention.
  • the effect of large-scale phenology extraction is affected by factors such as the geographical environment of the local geographical area and the growth characteristics of the researched crops.
  • the local geographic area of this example is the Northeast of my country, and the research crop is spring corn.
  • This example uses the morphological model method to analyze the four key phenology of three-leaf, heading, milking, and heading stages of spring corn in Northeast my country. Phenomenological extraction was carried out during the period.
  • the local geographic area and the distribution map of agricultural gas stations three agricultural gas stations are selected in the local geographic area, namely 1-Shuangcheng, 2-Lishu, and 3-Haicheng.
  • the embodiment of the present invention and its implementation process include the following:
  • Step S1 Vegetation index time series curve acquisition: According to the crop surface classification data, the MODIS satellite remote sensing data corresponding to the pixel where the crop is spring corn in Northeast my country is obtained, the vegetation index of each pixel is calculated, and the corresponding vegetation index time series curve is obtained . According to the geographic attribute data of the three agricultural gas stations, the MODIS satellite remote sensing data and vegetation index time series curves of the corresponding pixels from 2003 to 2013 are obtained.
  • the vegetation index used in this embodiment is the Wide Dynamic Range Vegetation Index (WDRVI).
  • WDRVI Wide Dynamic Range Vegetation Index
  • the wide dynamic range vegetation index has higher sensitivity in high biomass areas, and its calculation formula is shown in formula (1):
  • ⁇ NIR is the reflectivity in the near-infrared band
  • ⁇ RED is the reflectivity in the red band
  • is the weighting coefficient, which is 0.2 in this embodiment.
  • Step S2 Time series curve smooth fitting: As shown in Figure 3, the double logistic function fitting method is used to smoothly fit the vegetation index time series curve of each pixel, and the iterative nonlinear least square method is used to calculate the vegetation of the entire crop growth period. The exponential curve is fitted to the double logistic function, and the smoothed vegetation index time series curve is obtained.
  • the dual logic function expression is shown in formula (2):
  • Step S3 Morphological model establishment: The morphological model establishment process schematic diagram shown in Figure 4, the morphological model reference point uses the average value of the ground observation phenology period of the agricultural gas station for many years, and the morphological model reference curve uses the agricultural gas station to remove the abnormal data The mean value of the time series curve of vegetation index for many years of history.
  • the morphological model corresponds to the agricultural gas station one to one. In this embodiment, three morphological models are established for the three stations. As shown in FIG. 5, a schematic diagram of a morphological model established in this embodiment;
  • Step S4 morphological model fitting: According to the spatial relationship of the three agricultural gas stations of each pixel in the northeast region, the most suitable morphological model for each pixel is selected according to the minimum weighted geographic distance.
  • the morphological model transformation process as shown in Figure 6 is a schematic diagram of the morphological model reference curve fitted to the smooth vegetation index time series curve of the pixel, as shown in equation (3):
  • Table 1 The best scaling parameters of different provinces in the Northeast region (2003-2013)
  • Step S5 Regional phenological phase extraction: According to the optimal scaling parameters of each pixel, the morphological model reference point is transformed to obtain the phenological phase to be extracted for each pixel.
  • the calculation formula of the phenological phase extraction process is as formula (5):
  • the RMSE of the extraction accuracy of the four key phenological periods in this embodiment is 5.80 to 10.28 days.
  • the extraction accuracy of most phenological dates (91.8%) of the research pixels is within 15 days.
  • the specific phenological extraction results are shown in the following table:
  • the method of the present invention utilizes the macroscopic characteristics of the vegetation index curve during the growth period, can reduce the influence of local fluctuations and noise of the curve, and the phenological phase extraction result is accurate; in addition, the method can simultaneously extract the various phenological phases of the crop without requiring separate settings. Different metric standards can be applied to extract phenology of large-scale crops.

Landscapes

  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Multimedia (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Software Systems (AREA)
  • General Health & Medical Sciences (AREA)
  • Health & Medical Sciences (AREA)
  • Databases & Information Systems (AREA)
  • Artificial Intelligence (AREA)
  • Computing Systems (AREA)
  • Evolutionary Computation (AREA)
  • Medical Informatics (AREA)
  • Business, Economics & Management (AREA)
  • Human Resources & Organizations (AREA)
  • Astronomy & Astrophysics (AREA)
  • Remote Sensing (AREA)
  • Strategic Management (AREA)
  • Educational Administration (AREA)
  • Economics (AREA)
  • Entrepreneurship & Innovation (AREA)
  • Development Economics (AREA)
  • Quality & Reliability (AREA)
  • Game Theory and Decision Science (AREA)
  • Marketing (AREA)
  • Operations Research (AREA)
  • Tourism & Hospitality (AREA)
  • General Business, Economics & Management (AREA)
  • Image Processing (AREA)

Abstract

本发明公开了一种基于形态模型法的大范围作物物候提取方法。获取局部地理区域多年植被指数时序曲线;利用双重逻辑函数拟合方式对植被指数时序曲线进行平滑拟合;使用农气站的参考曲线和参考点建立形态模型;利用变换进行形态模型拟合;使用最优缩放参数计算获得局部地理区域物候期提取值。本发明的方法利用了曲线的宏观特征,可以减少曲线局部波动和噪声的影响,有更好的提取精度;并且可同时对作物的各个物候期进行提取,无需分别设定不同的度量标准,可应用于大范围作物物候提取。

Description

一种基于形态模型法的大范围作物物候提取方法 技术领域
本发明涉及农业遥感领域,具体涉及一种基于形态模型法的大范围作物物候提取方法。
背景技术
农业生产是人类社会的基石,广泛而可靠的农情信息是确保区域和全球粮食安全的重要基础。作物的物候体现在其生理特性及其对环境条件的生理响应随生长而发生的季节性变化,这种变化可以反映气候变化对农业生态系统的影响。作物的物候信息可以为作物灌溉调度和施肥管理提供重要的决策信息;准确的物候期也是作物精确产量预测模型的重要输入信息;另外,物候期也是遥感作物识别的重要基础。
目前物候监测的方法主要有野外观测法、积温模型法和遥感监测等。野外观测法相对简单易行,结果较为准确,但由于主观因素的影响,获得的数据质量难以保证,而且需要大量人力物力,观测周期长、面积小,难以应用在较大范围。积温模型法对输入信息要求较高,包括作物播种日期、气温和不同物候期所对应的积温信息等,这限制了它的应用范围。而遥感技术并无上述缺点,其具有监测范围广、时效高、客观准确等优势。
近年来,国内外学者对遥感物候提取展开了诸多研究,主要有以下几种方法:一、阈值法,即通过使用固定阈值或动态阈值来估计物候阶段;二、移动窗口法,即通过检测时序植被指数曲线在定义的移动时间窗口中的变化来确定物候时期;三、函数拟合法,即用如逻辑函数(Logistic Function)等数学函数来拟合待提取物候的作物的时序植被指数曲线,并根据拟合曲线上存在的特征点来确定物候期。然而这些方法主要关注植被指数时序曲线的局部形态特征点或特征段,这样一方面很容易受到时序曲线上的如噪声带来的局部波动的影响,从而给物候提取结果带来误差;另一方面不同物候期和不同作物需要应用不同的度量标准,方法的通用性较差,很难大范围应用。
发明内容
为了解决背景技术中存在的问题,本发明提出了一种基于形态模型法的大范围作物物候提取方法,利用卫星遥感数据建立植被指数时序曲线,使用形态模型法进行曲线拟合,结合地面观测数据,根据最优拟合缩放参数计算物候提取值,解决了目前遥感作物物候提取应用尺度小、准确性低、通用性差等局限性。
本发明所采用的技术方案是:
步骤S1:植被指数时序曲线获取:根据作物地表分类数据,在局部地理区域中选取待提取物候的作物所在的所有像元,基于卫星遥感图像数据,计算上述各个像元的植被指数,进而得到各个像元的植被指数时序曲线;并且,根据农气站的地理坐标数据,获取农气站所在各个像元的多年植被指数时序曲线;
所述作物地表分类数据是已有的数据,通过对卫星遥感数据进行作物识别,得到不同作物种植的空间分布信息。所述的像元为卫星遥感图像中的像元。
步骤S2:时序曲线平滑拟合:使用双重逻辑函数拟合方式对待提取物候的作物所在各个像元和农气站所在各个像元的植被指数时序曲线进行平滑拟合,得到平滑后的植被指数时序曲线;
步骤S3:形态模型建立:使用农气站记录的地面物候观测值建立形态模型参考点,并使用农气站所在像元的植被指数时序曲线建立对应形态模型参考曲线,由形态模型参考点和形态模型参考曲线构成形态模型,每个农气站建立一个形态模型,多个站点分别建立多个形态模型;
地面物候观测值是指地面作物的物候期,例如小麦的物候期。
步骤S4:形态模型拟合:根据局部地理区域各个像元和农气站的空间关系(空间关系是指距离远近、地理环境相似或者不相似),分别选取与待提取物候的作物所在各个像元对应的最合适的农气站的形态模型,将形态模型参考曲线拟合到像元的平滑植被指数时序曲线,计算拟合后的两曲线之间的均方根误差(RMSE),将均方根误差最小时对应的缩放参数作为该像元的最佳缩放参数,该缩放参数拟合效果最好;
步骤S5:区域物候期提取:根据最佳缩放参数将形态模型参考点进行变换,求得待提取物候的作物所在各个像元的物候提取结果。
本发明采用的植被指数具体为宽动态范围植被指数(WDRVI),所 述的宽动态范围植被指数在NDVI的基础上增加了一个加权系数,可以提高NDVI在高植被覆盖度区域的灵敏度。所述的宽动态范围植被指数(WDRVI)计算公式如下式所示:
Figure PCTCN2020124834-appb-000001
式中,ρ NIR为近红外波段反射率,ρ RED为红色波段反射率,α为加权系数,其值根据不同作物而变化,一般取0.1-0.2。
所述步骤S2中的双重逻辑函数拟合方式具体如下:将整个作物生长期的植被指数时序曲线拟合到双重逻辑函数曲线,双重逻辑函数包含六个参数,VI 1、VI m、S、A、S′和A′,拟合后的公式表示为:
Figure PCTCN2020124834-appb-000002
其中,VI 1表示植被指数时序曲线中处于作物非生长期一段(即前半段)的平均植被指数,VI m代表植被指数时序曲线的峰值,t表示时间,S和A分别代表植被指数时序曲线中上升拐点和下降拐点对应的时间,S′和A′分别代表上升拐点和下降拐点处的曲率;VI(t)表示时间t下的植被指数。
上升拐点是指斜率最大的拐点,下降拐点是指斜率最大的拐点。
所述步骤S3中,形态模型与农气站一一对应,代表作物在对应农气站一定空间范围内的典型季节性植被指数响应,形态模型参考点使用农气站多年地面观测物候期的平均值,形态模型参考曲线使用农气站所在的所有像元的多年历史植被指数时序曲线的均值。
所述步骤S4中,与待提取物候的作物所在各个像元最合适的形态模型根据像元和农气站的空间关系选取,综合考虑空间距离和地理环境因素并加权计算处理获得,具体实施可加权计算相似度,以相似度最高的作为最合适的。
所述步骤S4的形态模型拟合过程中,采用以下公式将形态模型参考曲线拟合变换到待提取物候的作物的平滑植被指数时序曲线,具体公式如下:
g(x)=yscale×h(xscale×(x 0+tshift))
其中,g(x)代表经过拟合变换后的形态模型参考曲线,h(x)代表形态模型参考曲线,xscale、yscale分别代表x轴、y轴方向上的缩放参数,tshift 代表x轴方向的偏移量,x 0表示在形态模型参考点;
再采用以下公式处理获得经过拟合变换后的形态模型参考曲线和待提取物候的作物的平滑植被指数时序曲线的两曲线间的均方根误差(RMSE):
Figure PCTCN2020124834-appb-000003
其中,x代表逐日日期,n代表曲线的时间间隔,f(x)代表到待提取物候的作物的平滑植被指数时序曲线。
根据均方根误差的最小值确定最佳缩放参数,此时的拟合效果最好。
所述步骤S5中,采用以下公式根据最佳缩放参数将形态模型参考点进行变换,求得待提取物候的作物所在各个像元待提取的物候期,综合局部地理区域所有像元物候期的提取值,得到该区域的物候时空分布情况:
X est=xscale×(x 0+tshift)
其中,x 0表示在形态模型参考点的物候期,X est是提取得到待提取物候的作物的物候期。
所述局部地理区域大小不限,至少包括一个像元。
所述待提取物候期包括作物整个生长期的各个物候期,以玉米作物为例,包括但不限于出苗期、三叶期、返青期、抽穗期、成熟期等。
所述像元为卫星遥感图像数据中的最小像素单元,根据不同卫星遥感产品,像元大小可以为30m×30m、250m×250m、500m×500m等。
本发明方法解决了现有方法进行作物提取物候的局限性。使用形态模型代表作物的典型生长曲线,再通过变换的方式拟合形态模型和作物生长曲线,最后使用得到的最优缩放参数计算待提取物候期。该方法利用到了植被指数时序曲线的宏观特征,可以减少曲线局部波动和噪声的影响,有更好的物候提取精度;另外,该方法可以同时对作物的各个物候期进行提取,无需分别设定不同的度量标准,可应用于大范围作物物候提取。
本发明的有益效果是:
本发明的方法利用了曲线的宏观特征,可以减少曲线局部波动和噪声的影响,有更好的提取精度;并且可同时对作物的各个物候期进行提取,无需分别设定不同的度量标准,可应用于大范围作物物候提取。
附图说明
图1为根据本发明一个具体实施例的局部地理区域及农气站点分布图;
图2为根据本发明一个具体实施例的一种基于形态模型法的大范围作物物候提取方法示意简图;
图3为根据本发明一个具体实施例的双重逻辑函数曲线平滑拟合示意
简图;
图4为根据本发明一个具体实施例的形态模型建立流程示意简图;
图5为根据本发明一个具体实施例的形态模型示意简图;
图6为根据本发明一个具体实施例的形态模型变换过程示意简图;
具体实施方式
以下结合附图对本发明的原理和特征进行描述,所举实例只用于解释本发明,并非用于限定本发明的范围。基于本发明中的实施例,本领域普通技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例,都属于本发明保护的范围。
大范围的物候提取效果受局部地理区域地理环境和研究作物生长特性等因素影响。本实施例的局部地理区域为我国的东北地区,研究作物为春玉米,本实施例使用形态模型法对我国东北地区春玉米的三叶期、抽穗期、乳熟期和抽穗期四个关键物候期进行了物候提取。如图1所示的局部地理区域及农气站点分布图,在局部地理区域选取三个农气站点,分别为1-双城、2-梨树、3-海城。
如图2所示,本发明的实施例及其实施过程包括以下:
步骤S1:植被指数时序曲线获取:根据作物地表分类数据,获取我国东北地区中种植作物为春玉米的像元对应的MODIS卫星遥感数据,计算各个像元的植被指数,得到对应的植被指数时序曲线。根据三个农气站的地理属性数据,获取对应像元的2003-2013年的MODIS卫星遥感数据和植被指数时序曲线。本实施例中使用的植被指数为宽动态范围植被指数(WDRVI),宽动态范围植被指数在高生物量区域的灵敏度较高,其计算公式如式(1)所示:
Figure PCTCN2020124834-appb-000004
式中,ρ NIR为近红外波段反射率,ρ RED为红色波段反射率,α为加权系数,在本实施例中取0.2。
步骤S2:时序曲线平滑拟合:如图3所示,使用双重逻辑函数拟合方式对各个像元的植被指数时序曲线进行平滑拟合,使用迭代非线性最小二乘法将整个作物生长期的植被指数曲线拟合到双重逻辑函数,得到平滑后的植被指数时序曲线。双重逻辑函数表示如式(2)所示:
Figure PCTCN2020124834-appb-000005
步骤S3:形态模型建立:如图4所示的形态模型建立流程示意简图,形态模型参考点使用农气站多年地面观测物候期的平均值,形态模型参考曲线使用农气站除去异常数据后的多年历史植被指数时序曲线的均值。形态模型与农气站一一对应,本实施例中三个站点分别建立三个形态模型。如图5所示,为本实施例建立的一个形态模型示意简图;
步骤S4:形态模型拟合:根据东北地区各个像元三个农气站的空间关系,依据最小加权地理距离分别选取各个像元最合适的形态模型。如图6所示的形态模型变换过程示意简图,将形态模型参考曲线拟合到像元的平滑植被指数时序曲线,如式(3)所示:
g(x)=yscale×h(xscale×(x 0+tshift))      (3)
计算两曲线之间的最小均方根误差(RMSE),作为拟合效果最好的缩放参数,且作为该像元的最佳缩放参数,最小均方根误差的计算如式(4)所示:
Figure PCTCN2020124834-appb-000006
本实施例计算得到的不同省份的最佳缩放参数结果如下表:
表1:东北区域不同省份的最佳缩放参数(2003-2013年)
Figure PCTCN2020124834-appb-000007
步骤S5:区域物候期提取:根据各个像元的最佳缩放参数,将形态模型参考点进行变换,求得各个像元待提取的物候期,物候期提取过程计算公式如式(5):
X est=xscale×(x 0+tshift)          (5)
本实施例中四个关键物候期的提取精度RMSE为5.80到10.28天。研 究像元大部分物候日期(91.8%)的提取精度在15天以内。具体的物候提取结果如下表所示:
表2:东北区域不同省份的物候提取结果(2003-2013年)
Figure PCTCN2020124834-appb-000008
本发明的方法利用了生长期内植被指数曲线的宏观特征,可以减少曲线局部波动和噪声的影响,物候期提取结果准确;另外,该方法可以同时对作物的各个物候期进行提取,无需分别设定不同的度量标准,可应用于大范围作物物候提取。

Claims (10)

  1. 一种基于形态模型法的大范围作物物候提取方法,包括如下步骤:
    步骤S1:植被指数时序曲线获取:根据作物地表分类数据,在局部地理区域中选取待提取物候的作物所在的所有像元,基于卫星遥感图像数据,计算上述各个像元的植被指数,进而得到各个像元的植被指数时序曲线;并且,根据农气站的地理坐标数据,获取农气站所在各个像元的多年植被指数时序曲线;
    步骤S2:时序曲线平滑拟合:使用双重逻辑函数拟合方式对待提取物候的作物所在各个像元和农气站所在各个像元的植被指数时序曲线进行平滑拟合,得到平滑后的植被指数时序曲线;
    步骤S3:形态模型建立:使用农气站记录的地面物候观测值建立形态模型参考点,并使用农气站所在像元的植被指数时序曲线建立对应形态模型参考曲线,由形态模型参考点和形态模型参考曲线构成形态模型,每个农气站建立一个形态模型,多个站点分别建立多个形态模型;
    步骤S4:形态模型拟合:根据局部地理区域各个像元和农气站的空间关系,分别选取与待提取物候的作物所在各个像元对应的最合适的形态模型,将形态模型参考曲线拟合到像元的平滑植被指数时序曲线,计算拟合后的两曲线之间的均方根误差(RMSE),将均方根误差最小时对应的缩放参数作为该像元的最佳缩放参数;
    步骤S5:区域物候期提取:根据最佳缩放参数将形态模型参考点进行变换,求得待提取物候的作物所在各个像元的物候提取结果。
  2. 根据权利要求1所述的一种基于形态模型法的大范围作物物候提取方法,其特征在于:所述的宽动态范围植被指数(WDRVI)计算如下:
    Figure PCTCN2020124834-appb-100001
    式中,ρ NIR为近红外波段反射率,ρ RED为红色波段反射率,α为加权系数,其值根据不同作物而变化,一般取0.1-0.2。
  3. 根据权利要求1所述的一种基于形态模型法的大范围作物物候提取方法,其特征在于:所述步骤S2中的双重逻辑函数拟合方式具体如下:将整个作物生长期的植被指数时序曲线拟合到双重逻辑函数曲线,拟合后的公式表示为:
    Figure PCTCN2020124834-appb-100002
    其中,VI 1表示植被指数时序曲线中处于作物非生长期一段(即前半段)的平均植被指数,VI m代表植被指数时序曲线的峰值,t表示时间,S和A分别代表植被指数时序曲线中上升拐点和下降拐点对应的时间,S′和A′分别代表上升拐点和下降拐点处的曲率;VI(t)表示时间t下的植被指数。
  4. 根据权利要求1所述的一种基于形态模型法的大范围作物物候提取方法,其特征在于:所述步骤S3中,形态模型参考点使用农气站多年地面观测物候期的平均值,形态模型参考曲线使用农气站所在的所有像元的多年历史植被指数时序曲线的均值。
  5. 根据权利要求1所述的一种基于形态模型法的大范围作物物候提取方法,其特征在于:所述步骤S4中,与待提取物候的作物所在各个像元最合适的形态模型根据像元和农气站的空间关系选取,综合考虑空间距离和地理环境因素并加权计算处理获得。
  6. 根据权利要求1所述的一种基于形态模型法的大范围作物物候提取方法,其特征在于:所述步骤S4的形态模型拟合过程中,采用以下公式将形态模型参考曲线拟合变换到待提取物候的作物的平滑植被指数时序曲线,具体公式如下:
    g(x)=yscale×h(xscale×(x 0+tshift))
    其中,g(x)代表经过拟合变换后的形态模型参考曲线,h(x)代表形态模型参考曲线,xscale、yscale分别代表x轴、y轴方向上的缩放参数,tshift代表x轴方向的偏移量,x 0表示在形态模型参考点;
    再采用以下公式处理获得经过拟合变换后的形态模型参考曲线和待提取物候的作物的平滑植被指数时序曲线的两曲线间的均方根误差(RMSE):
    Figure PCTCN2020124834-appb-100003
    其中,x代表逐日日期,n代表曲线的时间间隔,f(x)代表到待提取物候的作物的平滑植被指数时序曲线。
  7. 根据权利要求1所述的一种基于形态模型法的大范围作物物候提取方法,其特征在于:所述步骤S5中,采用以下公式根据最佳缩放参数将形态模型参考点进行变换,求得待提取物候的作物所在各个像元待提取的物候期,综合局部地理区域所有像元物候期的提取值,得到该区域的物候时空分布情况:
    X est=xscale×(x 0+tshift)
    其中,x 0表示在形态模型参考点的物候期,X est是提取得到待提取物候的作物的物候期。
  8. 根据权利要求1所述的一种基于形态模型法的大范围作物物候提取方法,其特征在于:所述局部地理区域大小不限,至少包括一个像元。
  9. 根据权利要求1所述的一种基于形态模型法的大范围作物物候提取方法,其特征在于:所述待提取物候期包括作物整个生长期的各个物候期。
  10. 根据权利要求1所述的一种基于形态模型法的大范围作物物候提取方法,其特征在于:所述像元为卫星遥感图像数据中的最小像素单元。
PCT/CN2020/124834 2019-11-19 2020-10-29 一种基于形态模型法的大范围作物物候提取方法 WO2021098471A1 (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
US17/778,014 US11734925B2 (en) 2019-11-19 2020-10-29 Large-scale crop phenology extraction method based on shape model fitting method

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
CN201911134130.3A CN111104858B (zh) 2019-11-19 2019-11-19 一种基于形态模型法的大范围作物物候提取方法
CN201911134130.3 2019-11-19

Publications (1)

Publication Number Publication Date
WO2021098471A1 true WO2021098471A1 (zh) 2021-05-27

Family

ID=70420758

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/CN2020/124834 WO2021098471A1 (zh) 2019-11-19 2020-10-29 一种基于形态模型法的大范围作物物候提取方法

Country Status (3)

Country Link
US (1) US11734925B2 (zh)
CN (1) CN111104858B (zh)
WO (1) WO2021098471A1 (zh)

Cited By (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113421255A (zh) * 2021-07-21 2021-09-21 中国科学院地理科学与资源研究所 一种基于栅格的耕地复种指数提取方法及系统
CN113642464A (zh) * 2021-08-13 2021-11-12 福州大学 结合twdtw算法和模糊集的时序遥感影像作物分类方法
CN113850139A (zh) * 2021-08-25 2021-12-28 南京林业大学 一种基于多源遥感的森林年际物候监测方法
CN114092831A (zh) * 2021-12-02 2022-02-25 中国科学院东北地理与农业生态研究所 一种草本沼泽植被物候信息提取方法
CN114332492A (zh) * 2021-12-15 2022-04-12 南京林业大学 一种基于长时间序列的植被特征提取及选择方法
CN114529097A (zh) * 2022-02-26 2022-05-24 黑龙江八一农垦大学 多尺度农作物物候期遥感降维预测方法
CN114544515A (zh) * 2022-02-23 2022-05-27 中国矿业大学 一种草原物候遥感监测方法及系统
CN114782838A (zh) * 2022-06-17 2022-07-22 中化现代农业有限公司 水稻识别方法、装置、电子设备和存储介质
CN117689959A (zh) * 2024-01-30 2024-03-12 中交第二公路勘察设计研究院有限公司 一种融合植被生命周期特征的遥感分类方法

Families Citing this family (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111104858B (zh) 2019-11-19 2022-07-15 浙江大学 一种基于形态模型法的大范围作物物候提取方法
CN112667708A (zh) * 2020-12-24 2021-04-16 河南中原光电测控技术有限公司 一种基于ndvi时间序列的冬小麦物候期监测方法及装置
CN113822360B (zh) * 2021-09-24 2022-04-22 中化现代农业有限公司 一种作物的物候期确定方法及装置
CN113673490B (zh) * 2021-10-21 2022-01-04 武汉大学 一种物候期自适应的作物生理参数遥感估测方法及系统
CN116579521B (zh) * 2023-05-12 2024-01-19 中山大学 产量预测时间窗口确定方法、装置、设备及可读存储介质
CN116597157B (zh) * 2023-07-11 2023-09-26 中国科学院地理科学与资源研究所 基于特征光谱变化的植物物候提取方法及装置
CN117423003B (zh) * 2023-12-19 2024-03-19 山东科技大学 一种冬小麦越冬期苗情分级遥感监测方法
CN117689120B (zh) * 2024-02-04 2024-04-26 北京工业大学 一种农业火灾的火辐射功率作物精细化分配方法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104850694A (zh) * 2015-05-13 2015-08-19 福州大学 基于生长期植被指数增量的冬小麦遥感监测方法
US20180330486A1 (en) * 2017-05-12 2018-11-15 Harris Lee Cohen Computer-implemented methods, computer readable medium and systems for a precision agriculture platform
CN111104858A (zh) * 2019-11-19 2020-05-05 浙江大学 一种基于形态模型法的大范围作物物候提取方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110175931B (zh) * 2019-05-10 2020-04-24 北京师范大学 一种大范围快速提取农作物种植面积和物候信息的方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104850694A (zh) * 2015-05-13 2015-08-19 福州大学 基于生长期植被指数增量的冬小麦遥感监测方法
US20180330486A1 (en) * 2017-05-12 2018-11-15 Harris Lee Cohen Computer-implemented methods, computer readable medium and systems for a precision agriculture platform
CN111104858A (zh) * 2019-11-19 2020-05-05 浙江大学 一种基于形态模型法的大范围作物物候提取方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
DING XIAO: "Study on Distribution of Crop's Structure in Heilongjiang", AGRICULTURE, CHINA MASTER’S THESES FULL-TEXT DATABASE, 1 June 2014 (2014-06-01), XP055814109 *
HU ZONGCHEN: "Mapping Major Crops Spatial and Temporal Distribution in China Using MODIS Time-Series Data", INFORMATION & TECHNOLOGY, CHINA MASTER'S THESES FULL-TEXT DATABASE, 23 June 2013 (2013-06-23), XP055814116 *

Cited By (13)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113421255A (zh) * 2021-07-21 2021-09-21 中国科学院地理科学与资源研究所 一种基于栅格的耕地复种指数提取方法及系统
CN113642464B (zh) * 2021-08-13 2023-09-15 福州大学 结合twdtw算法和模糊集的时序遥感影像作物分类方法
CN113642464A (zh) * 2021-08-13 2021-11-12 福州大学 结合twdtw算法和模糊集的时序遥感影像作物分类方法
CN113850139A (zh) * 2021-08-25 2021-12-28 南京林业大学 一种基于多源遥感的森林年际物候监测方法
CN113850139B (zh) * 2021-08-25 2022-04-08 南京林业大学 一种基于多源遥感的森林年际物候监测方法
CN114092831A (zh) * 2021-12-02 2022-02-25 中国科学院东北地理与农业生态研究所 一种草本沼泽植被物候信息提取方法
CN114332492A (zh) * 2021-12-15 2022-04-12 南京林业大学 一种基于长时间序列的植被特征提取及选择方法
CN114332492B (zh) * 2021-12-15 2024-03-19 南京林业大学 一种基于长时间序列的植被特征提取及选择方法
CN114544515A (zh) * 2022-02-23 2022-05-27 中国矿业大学 一种草原物候遥感监测方法及系统
CN114544515B (zh) * 2022-02-23 2024-05-14 中国矿业大学 一种草原物候遥感监测方法及系统
CN114529097A (zh) * 2022-02-26 2022-05-24 黑龙江八一农垦大学 多尺度农作物物候期遥感降维预测方法
CN114782838A (zh) * 2022-06-17 2022-07-22 中化现代农业有限公司 水稻识别方法、装置、电子设备和存储介质
CN117689959A (zh) * 2024-01-30 2024-03-12 中交第二公路勘察设计研究院有限公司 一种融合植被生命周期特征的遥感分类方法

Also Published As

Publication number Publication date
US11734925B2 (en) 2023-08-22
US20220406054A1 (en) 2022-12-22
CN111104858B (zh) 2022-07-15
CN111104858A (zh) 2020-05-05

Similar Documents

Publication Publication Date Title
WO2021098471A1 (zh) 一种基于形态模型法的大范围作物物候提取方法
CN104778451B (zh) 一种考虑草地高度因子的草地生物量遥感反演方法
CN111368736B (zh) 基于sar和光学遥感数据的水稻精细化估产方法
CN108446999B (zh) 基于冠气温差和遥感信息进行灌区不同作物et估算方法
CN106897707B (zh) 基于多源中分的特征影像时间序列合成方法及装置
CN111241912A (zh) 一种基于机器学习算法的多植被指数的水稻产量估测方法
CN101858971A (zh) 一种基于modis数据的水稻单产遥感估算方法
CN107122739B (zh) 基于Extreme数学模型重构VI时间序列曲线的作物估产模型
CN109284771A (zh) 一种番茄生长模型判定方法及装置
CN104850694A (zh) 基于生长期植被指数增量的冬小麦遥感监测方法
CN102663396B (zh) 一种水稻乳熟期自动检测的方法
CN102542560A (zh) 一种水稻移栽后密度自动检测的方法
CN103971199A (zh) 一种大范围农作物长势的遥感评级方法
CN105954207A (zh) 一种枣树冠层叶绿素b含量的检测方法
Hu et al. Phenotyping of poplar seedling leaves based on a 3D visualization method
CN114881127A (zh) 一种基于高分辨率遥感卫星影像的农作物精细分类方法
CN103413172A (zh) 玉米苗期缺苗数的测量方法及装置
CN110910266A (zh) 一种基于遥感的灌区农业灌溉用水总量估算方法
CN103471514B (zh) 一种基于机器视觉和一元线性回归分析的大蒜分级方法
CN102621075B (zh) 一种水稻抽穗期自动检测的方法
CN109035209A (zh) 甘蔗分蘖期自动观测方法
CN107491724A (zh) 一种冬小麦穗分化期识别方法及装置
TAO et al. Fusing multi-source data to map spatio-temporal dynamics of winter rape on the Jianghan Plain and Dongting Lake Plain, China
Wei et al. Application of remote sensing technology in crop estimation
CN115841615A (zh) 一种基于无人机多光谱数据的烟草产量预测方法及装置

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 20890173

Country of ref document: EP

Kind code of ref document: A1

NENP Non-entry into the national phase

Ref country code: DE

122 Ep: pct application non-entry in european phase

Ref document number: 20890173

Country of ref document: EP

Kind code of ref document: A1