LARGE-SCALE CROP PHENOLOGY EXTRACTION METHOD BASED ON SHAPE MODEL FITTING METHOD

Abstract

Disclosed is a large-scale crop phenology extraction method based on a shape model fitting method. The method comprises: acquiring a multi-year vegetation index time sequence curve in a localized geographic region; performing smooth fitting on the vegetation index time sequence curve by using a dual logistic function fitting means; establishing shape models by using reference curves and reference points of agrometeorological stations; performing shape model fitting by means of transformation; and obtaining a phenological period extraction value of the localized geographic region by means of calculation using the optimal scaling parameter. According to the present invention, macroscopic features of the curve are used, such that the influence of localized fluctuation and noise of the curve can be reduced, and a better extraction precision is obtained; and each phenological period of a crop can be extracted at the same time.

Claims

1. A large-scale crop phenology extraction method based on a shape model fitting method, comprising: step S1: obtaining a vegetation index time sequence curve, wherein all pixels where a crop whose phenology is to be extracted is located in a localized geographic region are selected according to crop surface classification data, a vegetation index of each of the pixels is calculated based on satellite remote sensing image data to obtain the vegetation index time sequence curve of each of the pixels, and a multi-year vegetation index time sequence curve of pixels where agrometeorological stations are located is obtained according to geographic coordinate data of the agrometeorological stations; step S2: performing smooth fitting on the vegetation index time sequence curve, wherein the vegetation index time sequence curve of each of the pixels where the crop whose phenology is to be extracted is located and each of the pixels where the agrometeorological stations are located are smoothly fitted by using a dual logistic function fitting means, and a smooth vegetation index time sequence curve is obtained; step S3: establishing shape models, wherein a shape model reference point is established by using a ground phenological observation value recorded by the agrometeorological stations, a shape model reference curve, which is corresponded, is established by using the vegetation index time sequence curve of each of the pixels where the agrometeorological stations are located, the shape models is formed by using the shape model reference point and the shape model reference curve, one of the shape models is established for each of the agrometeorological stations, and a plurality of the shape models are established for a plurality of station sites; step S4: performing a shape model fitting, wherein a most suitable shape model corresponding to each of the pixels where the crop whose phenology is to be extracted is located is selected according to a spatial relationship between each of the pixels and the agrometeorological stations in the localized geographic region, the shape model reference curve is fitted to the smooth vegetation index time sequence curve of the pixels, a root mean square error (RMSE) between two curves, which are fitted, is calculated, and a scaling parameter corresponding to a minimum root mean square error is treated as an optimal scaling parameter of the pixels; and step S5: extracting phenological periods from a region, wherein the shape model reference point is transformed according to the optimal scaling parameter to obtain a phenology extraction result of each of the pixels where the crop whose phenology is to be extracted is located.

2. The large-scale crop phenology extraction method based on the shape model fitting method according to claim 1, wherein a wide dynamic range vegetation index (WDRVI) is calculated as follows: $\mathrm{WDRVI}=\frac{\alpha .Math.{\rho }_{\mathrm{NIR}}-{\rho }_{RED}}{\alpha .Math.{\rho }_{\mathrm{NIR}}+{\rho }_{RED}}$ in the formula, ρ.sub.NIR is a reflectance in a near-infrared band, ρ.sub.RED is a reflectance in a red band, and α is a weighting coefficient, whose value varies according to different crops, generally is 0.1 to 0.2.

3. The large-scale crop phenology extraction method based on the shape model fitting method according to claim 1, wherein the dual logistic function fitting means in the step S2 specifically comprises: fitting the vegetation index time sequence curve of an entire growth period of the crop to a double logistic function curve, wherein a fitted formula is expressed as follows: $\mathrm{VI}\left(t\right)=V{I}_1+\left(V{I}_m-V{I}_1\right)\times \left(\frac{1}{1+\exp \left(-S^{\prime }\times \left(t-S\right)\right)}+\frac{1}{1+\exp \left(A^{\prime }\times \left(t-A\right)\right)}-1\right)$ where VI.sub.l represents a mean vegetation index in a crop non-growing period (i.e., a first half) in the vegetation index time sequence curve, VI.sub.m represents a peak value of the vegetation index time sequence curve, t represents time, S and A respectively represent time corresponding to a rising inflection point and a falling inflection point in the vegetation index time sequence curve, S′ and A′ respectively represent curvatures at the rising inflection point and the falling inflection point, and VI(t) represents the vegetation index at time t.

4. The large-scale crop phenology extraction method based on the shape model fitting method according to claim 1, wherein in the step S3, the shape model reference point uses a mean value of multi-year ground observation phenological periods of the agrometeorological stations, and the shape model reference curve uses a mean value of multi-year historical vegetation index time sequence curves of each of the pixels where the agrometeorological stations are located.

5. The large-scale crop phenology extraction method based on the shape model fitting method according to claim 1, wherein in the step S4, a selection of the most suitable shape model for each of the pixels where the crop whose phenology is to be extracted is located according to the spatial relationship between the pixels and the agrometeorological stations are comprehensively made based on a spatial distance and a geographical environment factor and is processed and obtained by weighted calculation.

6. The large-scale crop phenology extraction method based on the shape model fitting method according to claim 1, wherein in the shape model fitting performed in the step S4, the following formula is used to fit and transform the shape model reference curve to the smooth vegetation index time sequence curve of the crop whose phenology is to be extracted, and a specific formula is as follows:
g(x)=yscale×h(xscale×(x.sub.0+tshift)) where g(x) represents the shape model reference curve after fitting and transformation, h(x) represents the shape model reference curve, xscale and yscale respectively represent the scaling parameter in x-axis and y-axis directions, tshift represents an offset in the x-axis direction, and x.sub.0 represents the shape model reference point, and the following formula is then used to process and obtain the root mean square error (RMSE) between the two curves of the shape model reference curve after fitting and transformation and the smooth vegetation index time sequence curve of the crop whose phenology is to be extracted: $\begin{array}{cc}\mathrm{RMSE}=\sqrt{\frac{1}{365/n}\underset{t=n,2n,.Math.365}{.Math.}{\left(f\left(x\right)-g\left(x\right)\right)}^2}& \end{array}$ where x represents a daily date, n represents a time interval of a curve, and f(x) represents the smooth vegetation index time sequence curve of the crop whose phenology is to be extracted.

7. The large-scale crop phenology extraction method based on the shape model fitting method according to claim 1, wherein in the step S5, the following formula is used to transform the shape model reference point according to the optimal scaling parameter to obtain the phenological periods to be extracted of each of the pixels where the crop whose phenology is to be extracted is located, and extracted values of the phenological periods of all of the pixels in the localized geographic region are integrated to obtain a phenological distribution in time and space of the region:
X.sub.est=xscale×(x.sub.0+tshift) where x.sub.0 represents a phenological period of the shape model reference point, and X.sub.est is an extracted and obtained phenological period of the crop whose phenology is to be extracted.

8. The large-scale crop phenology extraction method based on the shape model fitting method according to claim 1, wherein a size of the localized geographic region is not limited and at least comprises one of the pixels.

9. The large-scale crop phenology extraction method based on the shape model fitting method according to claim 1, wherein the phenological periods to be extracted comprises each of phenological periods of an entire growth period of the crop.

10. The large-scale crop phenology extraction method based on the shape model fitting method according to claim 1, wherein the pixels are a smallest pixel unit in the satellite remote sensing image data.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

[0029] FIG. 1 is a map of a localized geographic region and distribution of agrometeorological station sites according to a specific embodiment of the invention.

[0030] FIG. 2 is a schematic chart of a large-scale crop phenology extraction method based on a shape model fitting method according to a specific embodiment of the invention.

[0031] FIG. 3 is a schematic graph of smooth fitting of a double logistic function curve according to a specific embodiment of the invention.

[0032] FIG. 4 is a schematic flow chart of establishment of a shape model according to a specific embodiment of the invention.

[0033] FIG. 5 is a schematic graph of the shape model according to a specific embodiment of the invention.

[0034] FIG. 6 is a schematic graph of a transformation process of the shape model according to a specific embodiment of the invention.

DESCRIPTION OF THE EMBODIMENTS

[0035] The principles and features of the invention are described below together with the accompanying drawings, and the examples are only used to explain the invention and is not used to limit the scope of the invention. Based on the embodiments of the invention, all other embodiments obtained by a person having ordinary skill in the art without making any inventive effort fall within the scope that the invention seeks to protect.

[0036] The effect of large-scale phenology extraction is affected by factors such as the geographical environment of the localized geographic region and the growth characteristics of the research crop. A localized geographic region provided in this embodiment is the northeastern region of China, and the research crop is spring maize. Four key phenological periods of the three-leaf period, heading period, silking period, and heading period of spring maize in Northeast China are extracted by using a shape model fitting method in this embodiment. In a map of a localized geographic region and distribution of agrometeorological station sites as shown in FIG. 1, three agrometeorological station sites are selected in the localized geographic region, namely 1-Shuangcheng, 2-Lishu, and 3-Haicheng.

[0037] As shown in FIG. 2, the embodiments of the invention and the implementation process thereof include the following step.

[0038] In step S1, a vegetation index time sequence curve is obtained. According to crop surface classification data, MODIS satellite remote sensing data corresponding to pixels of a planted crop of spring maize in Northeast China is obtained, and a vegetation index of each pixel is calculated to obtain the corresponding vegetation index time sequence curve. According to geographic attribute data of the three agrometeorological stations, the MODIS satellite remote sensing data and vegetation index time sequence 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). The wide dynamic range vegetation index has high sensitivity in high biomass areas, and its calculation formula is shown in formula (1):

[00004]$\begin{array}{cc}\mathrm{WDRVI}=\frac{\alpha .Math.{\rho }_{\mathrm{NIR}}-{\rho }_{RED}}{\alpha .Math.{\rho }_{\mathrm{NIR}}+{\rho }_{RED}}& \left(1\right)\end{array}$

in the formula, ρ.sub.NIR is reflectance in a near-infrared band, ρ.sub.RED is reflectance in a red band, and α is a weighting coefficient, which is 0.2 in this embodiment.

[0039] In step S2, smooth fitting is performed on the time sequence curve. As shown in FIG. 3, the vegetation index time sequence curve of each pixel is smoothly fitted by using a double logistic function fitting means, the vegetation index curve of an entire growth period of the crop is fitted to a double logistic function by using an iterative nonlinear least squares method, and a smooth vegetation index time sequence curve is obtained. The double logistic function is expressed as formula (2):

[00005]$\begin{array}{cc}\mathrm{VI}\left(t\right)={\mathrm{VI}}_1+\left({\mathrm{VI}}_m-{\mathrm{VI}}_1\right)\times \left(\frac{1}{1+\exp \left(-S^{\prime }\times \left(t-S\right)\right)}+\frac{1}{1+\exp \left(A^{\prime }\times \left(t-A\right)\right)}-1\right)& \left(2\right)\end{array}$

[0040] In the step S3, a shape model is established. In the schematic flow chart of establishment of a shape model as shown in FIG. 4, a shape model reference point uses a mean value of multi-year ground observation phenological periods of one agrometeorological station, and a shape model reference curve uses a mean value of multi-year historical vegetation index time sequence curves of the agrometeorological station after abnormal data is removed. The shape model corresponds with the agrometeorological station one-to-one. In this embodiment, three shape models are established for the three agrometeorological stations. FIG. 5 is a schematic graph of one shape model provided by this embodiment.

[0041] In step S4, shape model fitting method is performed. According to a spatial relationship of the three agrometeorological stations of the pixels in Northeast China, a most suitable shape model is selected for each pixel according to a minimum weighted geographic distance. In a schematic graph of a transformation process of the shape model as shown in FIG. 6, the shape model reference curve is fitted to the smooth vegetation index time sequence curve of the pixel, as shown in formula (3):

g(x)=yscale×h(xscale×(x.sub.0+tshift))  (3)

[0042] A root mean square error (RMSE) between the two curves is calculated to be treated as a scaling parameter exhibiting an optimal fitting effect, and acting as an optimal scaling parameter for the pixel, the calculation of the minimum root mean square error is shown in formula (4):

[00006]$\begin{array}{cc}RMSE=\sqrt{\frac{1}{73}{.Math.}_{t=5,10,.Math.365}{\left(f\left(x\right)-g\left(x\right)\right)}^2}& \left(4\right)\end{array}$

[0043] The results of the optimal scaling parameters for different provinces calculated in this embodiment are as follows:

TABLE-US-00001 TABLE 1 Optimal scaling parameters for different provinces in the Northeast region (2003-2013) xscale yscale tshift Province Mean Max Min Mean Max Min Mean Max Min Heilongjiang 1.01 1.17 0.9 1.01 1.17 0.9 1.01 1.17 0.9 Jilin 0.99 1.14 0.9 0.99 1.14 0.9 0.99 1.14 0.9 Liaoning 1.03 1.24 0.9 1.03 1.24 0.9 1.03 1.24 0.9

[0044] In step S5, a phenological period is extracted from the region. According to the optimal scaling parameter of each pixel, the shape model reference point is transformed to obtain the phenological period of each pixel to be extracted. The calculation formula of an extraction process of the phenological period is provided as follows:

X.sub.est=xscale×(x.sub.0+tshift)  (5)

[0045] The extraction precision RMSE of the four key phenological periods in this embodiment is 5.80 to 10.28 days. The extraction precision of most phenological dates (91.8%) of the studied pixels is within 15 days. Specific phenological extraction results are shown in the following table:

TABLE-US-00002 TABLE 2 Phenology extraction results of different provinces in Northeast China (2003-2013) Phenological period Region RMSE (days) Three-leaf period Heilongjiang 5.90 Jilin 5.80 Liaoning 6.77 All 6.17 Heading period Heilongjiang 8.07 Jilin 7.02 Liaoning 8.46 All 7.87 Silking period Heilongjiang 10.49 Jilin 7.29 Liaoning 10.24 All 9.45 Maturity period Heilongjiang 12.85 Jilin 10.25 Liaoning 10.58 All 11.28

[0046] According to the method of the present invention, macroscopic features of the vegetation index curve in the growth period are used, such that the influence of localized fluctuation and noise of the curve can be reduced, and precise extraction results are obtained. In addition, in the method, each phenological period of the crop may be extracted at the same time, and there is no need to respectively set different measurement standards, such that the method can be applied to large-scale crop phenology extraction.