CN110909821B - A method for high spatiotemporal resolution vegetation index data fusion based on crop reference curve - Google Patents

Abstract

The invention discloses a method for reconstructing a vegetation index data set with high temporal and spatial resolution based on a crop reference curve. The method includes the following steps: Step 1: constructing a sample library of crop reference curves; Multiple MODIS reference curves of different crops are collected to form a MODIS crop reference curve sample library. At the same time, combined with MODIS local curves, they are used as MODIS candidate curves to participate in curve matching; Step 2: Determine the initial fitting curve; Step 3: Determine the final fitting curve According to the local difference between the initial Landsat fitting curve and the original Landsat time series NDVI images, the final Landsat fitting curve is obtained after quadratic fitting and smoothing, thereby obtaining the MODIS temporal resolution and Landsat spatial resolution at the same time. high spatial and temporal resolution vegetation index data.

Description

High spatiotemporal resolution vegetation index data fusion based on crop reference curve method

technical field

The invention relates to the technical field of agricultural remote sensing, in particular to a method for reconstructing a vegetation index with high temporal and spatial resolution based on a crop reference curve, which is suitable for agricultural remote sensing monitoring research of different farmland systems.

Background technique

When monitoring the dynamic changes of the Earth system at the regional or global scale, such as changes in vegetation cover and land use, Earth observation satellites can provide abundant data and monitoring methods to obtain surface information, thereby providing scientific basis for resource management and policy formulation. However, due to the mutual constraints of high temporal resolution and high spatial resolution of remote sensing data obtained by earth observation satellites, remote sensing monitoring represented by field-scale research areas has relatively high temporal and spatial resolutions. High requirements and therefore specific methods are required to generate remote sensing data with high spatial and temporal resolution. Taking Landsat and MODIS images as examples, Landsat images have a medium-high spatial resolution of 30 m, which can better meet the requirements of field-scale remote sensing monitoring, but its 16-day reentry cycle and serious cloud coverage greatly limit Landsat. On the contrary, MODIS images have the ability to obtain data every day, but its low spatial resolution (250-1000m) makes it difficult to achieve accurate monitoring of complex farmland systems.

In order to meet the above-mentioned requirements of remote sensing monitoring for high temporal and spatial resolution data, if only from the perspective of data sources, it is possible to obtain enough high temporal and spatial resolution images by launching more Earth observation satellites with different characteristics. This method requires more advanced software and hardware technology and huge R&D and production costs, and is not practical. Therefore, based on the existing multi-source remote sensing images, the method of using spatiotemporal fusion technology to obtain high spatiotemporal resolution images has attracted attention since it was proposed, and has achieved great development and application. This technology combines the time, space and spectral characteristics of different sensors well, and through the algorithm, there is a certain correspondence between the high-spatial and low-temporal resolution data and the high-temporal and low-spatial resolution data obtained by different satellites in the same area. Simulate and fuse to generate images with high temporal and spatial resolution, which well solves the limitation of inconsistent temporal and spatial resolutions on the efficiency and accuracy of remote sensing monitoring. In the development process in recent decades, the classical spatiotemporal fusion algorithm STARFM has achieved good application results in many fields such as vegetation dynamic change monitoring, crop phenology information extraction, water resources monitoring, and surface temperature monitoring, and has also produced a series of Improved algorithms, such as STAARCH, ESTARFM, FSDAF, etc. The theoretical basis of this type of spatiotemporal fusion algorithm is usually that for high and low spatial resolution images of the same region and the same time period, the image features of the two in the same time phase are in a one-to-one correspondence, so the continuous images on the low spatial resolution images can be combined. The temporal phase change feature is mapped to the high spatial resolution image, which makes up for the lack of temporal phase change information caused by the long acquisition period of the high spatial resolution image, and thus well combines the spatial structure information and the The temporal dimension information provided by low spatial resolution imagery generates image data with high spatial and temporal resolution. However, this type of spatiotemporal fusion method also has some limitations. First, the low spatial resolution images that provide temporal dimension information have serious mixed pixel phenomenon, resulting in the obtained spatiotemporal fusion product with high temporal and spatial resolution, but the ground The accuracy of object recognition and classification is often not high, and misjudgment is prone to occur, especially for the classification and recognition of linear objects, small objects or irregular objects. Therefore, this type of algorithm is more suitable for remote sensing in spatially homogeneous areas. Monitoring identifies, but is less suitable for areas with greater spatial heterogeneity. Secondly, this type of spatiotemporal fusion algorithm requires multiple pairs of high spatial low temporal resolution and low spatial high temporal resolution images acquired under clear sky conditions and matching and fitting. Due to the high temporal resolution of low spatial resolution images, It can meet the requirements for the number of available images under clear sky conditions, but high spatial resolution images are easily affected by factors such as cloud coverage due to the long acquisition period, which limits the number of images that can be used for matching. Especially in areas with cloudy and rainy weather, there are often fewer high spatial resolution images available, and the error of fusion results is large.

Through existing research, it is found that 84% of the global agricultural production units have an area of less than 0.02km ² (approximately 140*140m), and there are certain regional differences in different regions. For example, the area of agricultural production units in the United States is relatively large. , the average and median values can reach 0.193km2 and 0.278km2, and more than half of the fields in some African countries are less than ^0.004km2 (approximately 63*63m), and such smaller fields usually have larger high spatial heterogeneity. Therefore, when using remote sensing means to monitor the agricultural production status at the field scale, the commonly used spatiotemporal fusion algorithm has great limitations, and the application effect is not ideal. At the same time, considering that field-scale agricultural production units are widely distributed around the world, the accurate monitoring of field-scale agricultural production activities through remote sensing methods is of great significance for reducing global poverty and ensuring food security.

SUMMARY OF THE INVENTION

The technical problem to be solved by the present invention is to provide a method for fusion of vegetation index data with high spatial and temporal resolution based on crop reference curves, and to perform spatial and temporal fusion reconstruction based on multi-source remote sensing images, so as to obtain a method suitable for spatial heterogeneity. High temporal and spatial resolution vegetation index data for remote sensing monitoring of higher complex farmland systems.

In order to solve the above technical problems, the present invention proposes a method for fusing vegetation index data with high spatial and temporal resolution based on crop reference curves.

A method for high temporal and spatial resolution vegetation index data fusion based on crop reference curves, comprising: step 1: constructing a sample library of crop reference curves; obtaining different crops based on pure pixels by using crop category information and spatial positions of a crop classification map CDL Multiple MODIS reference curves were obtained to form a MODIS crop reference curve sample library, and combined with MODIS local curves as MODIS candidate curves to participate in curve matching; Step 2: Determine the initial fitting curve; Fit and match the image pixel values at the corresponding positions of all the curves in the MODIS candidate curve, and find the candidate curve with the highest fitting degree as the best reference curve based on the MODIS crop reference curve; then use the fitting parameters of the best reference curve Fitting the Landsat time series NDVI images to obtain the Landsat initial fitting curve; Step 3: Determine the final fitting curve; for the local differences between the Landsat initial fitting curve and the original Landsat time series NDVI images, the second fitting is performed. The final fitting curve of Landsat is obtained after combined and smoothed, thereby obtaining vegetation index data with high spatial and temporal resolution with MODIS temporal resolution and Landsat spatial resolution at the same time.

Said method, said step 1 includes the following steps:

Step A: Use the crop category information in the crop classification map CDL to set a moving window with a size of 25*25 pixels. By moving the pixel window, when more than 95% of the pixels in the window belong to the crop classification map. When there is a single crop category, the window pixel is considered to be a pure pixel, and the position information of the CDL pure pixel window is obtained;

Step B: Preprocess the low spatial resolution MODIS data by resampling and reprojection so that it has the same spatial resolution and projection information as the high spatial resolution Landsat data and the crop classification map CDL data. The data is geometrically registered, and then the MODIS time series vegetation index NDVI data is obtained by calculation; at the same time, some Landsat7 images in the high spatial resolution Landsat data are interpolated to obtain the data covering the study area completely, and then the Landsat time series NDVI data is obtained by calculation. data;

Step C: Using the position information of the CDL pure pixel window obtained in Step A, select the position corresponding to the part of the pixel area in the middle of the window on the MODIS time series vegetation index NDVI data, and extract the position of the pixel contained in the position. The average NDVI value is obtained to obtain multiple MODIS time series NDVI curves of different types of crops; for the outliers that still exist in some curves, linear interpolation method is used to eliminate the outliers, and then these curves are used as MODIS reference curves of different ground objects to construct samples. library.

Said method, said step 2 includes the following steps:

Step D: Use the original value of the MODIS time series vegetation index NDVI data calculated in step B to obtain the time series vegetation index NDVI curve of each pixel of the MODIS image, which is called the MODIS local curve of the pixel, and then comprehensively consider The MODIS local curve and the MODIS reference curve sample library obtained in step C are combined together to participate in the curve matching as a MODIS candidate curve;

Step E: Using the Landsat time series NDVI data calculated in step B, analyze the correlation index R ² of all image pixels of the data and all curves in the MODIS crop reference curve sample library obtained in step C; and the Landsat time The correlation index R ^{2 of all image pixels of the sequence NDVI data and the MODIS local curve obtained in step D, the curve with the largest correlation index R 2 is} used as the best reference curve of MODIS for the image pixels of the Landsat time series NDVI data, and then through Formula (1) is used to further fit the Landsat time series NDVI data and the MODIS best reference curve in step B to obtain the Landsat fitting curve.

L(x)=a*M(x+x0)+b (1)

where x represents a certain day of the year, L(x) is the Landasat time series NDVI curve function, M(x) is the MODIS optimal reference curve function, and x0 represents the possible time offset between the two curves The range is between ±30 days, a and b are the fitting parameters calculated by the least squares method;

Step F: Calculate the correlation index R 2 between the pixel values of the Landsat time series NDVI image in Step B and the pixel values of all candidate curves of the MODIS at the corresponding position, analyze the correlation index R ² between the ^two , and select the most MODIS value. best reference curve;

Step G: Use steps E and F to confirm the best reference curve. At this time, the best time offset x0 is 0, and the corresponding fitting parameters a and b are obtained at the same time. Then use the MODIS best reference curve and the corresponding Fit the parameters to get the Landsat initial fit curve.

In the method, the third step includes the following steps: Step H: Calculate the difference between the NDVI pixel value (clear sky pixel) of the original Landsat time series and the NDVI fitting value of the corresponding pixel on the Landast initial fitting curve. The difference between the two values, the linear interpolation method is used to fit the linear interpolation result to the initial fitting curve to reduce the difference between the two, and finally the Gaussian function in the IDL programming platform is used to further smooth the optimized fitting curve. , so that the curve contains all the available original Landsat time series NDVI values as much as possible, and the curve is called the final fitting curve.

In the described method, in step F, for each MODIS crop reference curve, there are 60 correlation indices R ² , the value range of which is concentrated between 0.11 and 0.92, and then the candidate curve with the largest correlation index is selected as the MODIS crop. The curve to be matched of the reference curve, and then compare the correlation index between the curve to be matched and the Landsat time series NDVI data corresponding to all MODIS crop reference curves, as well as the correlation index between the MODIS local curve and the Landsat time series NDVI data. The value range is centered between 0.78 and 0.92. If the MODIS local curve correlation index R ² value is the largest, the MODIS local curve is directly used as the best reference curve; if the correlation index R ² maximum value appears in the curve to be matched in the MODIS crop reference curve, the correlation index R of all candidate curves is calculated. ² , select the candidate curve higher than the average value of R ² and average it again, and take the final average curve as the best reference curve.

In the described method, the MODIS reference curves of different ground objects are obtained by using the pure pixel area selected by the results of the known crop classification map CDL data with higher classification accuracy, and the specific setting parameters in the step A are based on: Empirical parameters of multiple experiments.

Beneficial effects:

According to the feature of linear relationship between vegetation index data of different spatial resolutions of different sensor images, the invention proposes a data spatiotemporal fusion method based on crop reference curve, which is suitable for remote sensing monitoring of complex farmland systems with high spatial heterogeneity. Using the CDL data of the known crop classification layer, the moving window was set to determine the pure pixel area of different types of crops, and the MODIS crop reference curve sample library was constructed based on the pure pixel area. By comparing the time series NDVI of Landsat image pixels The correlation between the data and the MODIS crop reference curve and the local curve of the original NDVI data of MODIS corresponding pixels, the selection method of the best reference curve is constructed, and based on the fitting parameters of the best reference curve and the original Landsat time series NDVI data The initial fitting curve is obtained; the role of the original effective value of Landsat in the curve fitting is preserved to the greatest extent by methods such as linear interpolation and Gaussian function fitting, which further reduces the phenomenon of initial overfitting, and the final fitting curve is obtained, This results in a vegetation index dataset with both high temporal and high spatial resolution. This method comprehensively considers the two major problems of serious mixed pixels in low spatial resolution images and the lack of matching data in traditional fusion methods, and uses known reliable classification data to reduce the impact of mixed pixels. The time dimension information of the optimal reference curve obtained from the high-resolution image is "transferred", which avoids the problem of limited matching number of high-spatial-resolution images caused by cloud coverage, and can maximize the use of all available original image information to generate The vegetation index data with high spatial and temporal resolution can be well applied in remote sensing monitoring of complex farmland systems.

Description of drawings

Figure 1 shows the MODIS reference curves of four typical features; the upper left is forest, the upper right is double crop, the lower left is corn, and the lower right is soybean;

Figure 2 shows the MODIS reference curve, the best reference curve, the initial fitting curve, the final fitting curve and the Landsat time series NDVI values of four typical ground objects; the upper left is corn, the upper right is soybean, the lower left is double crop, the lower right is the forest; (1) MODIS reference curve, (2) best reference curve, (3) initial fitting curve, (4) final fitting curve, dots represent Landsat time series NDVI values;

Figure 3 shows the comparison of NDVI fitting curves based on crop reference curve and STARFM on three crop types in 2013; Figure (a) is the crop classification diagram, Figures (b), (c), and (d) are respectively ( a) Curves of corn, soybean, and double-crop crops corresponding to the dot positions. The small dots in the figures (b), (c), and (d) are the time-series Landsat images, and the large dots are the STARFM paired images. (1) is the fitting curve based on the crop reference curve method, (2) is the STARFM method fitting curve, and (3) is the MODIS reference curve;

Detailed ways

The present invention will be described in detail below with reference to specific embodiments.

The present invention proposes a method for high temporal and spatial resolution vegetation index data fusion based on a crop reference curve, including the following steps:

S1: The study area of this example is located in the Choptank River Basin on the east coast of the United States, including parts of Maryland, Delaware, and New Jersey. The land use pattern in this area is mainly 58% agricultural land, 33% forest and With only 9% of urban built-up land, the most common cropping systems include corn-soybean rotation and double-cropping soybean-wheat. Taking the 2013 and 2014 crop classification map CDL data as reliable crop distribution data, and then setting a moving window with a size of 25*25 pixels, by moving the pixel window, more than 95% of the pixels in the window belong to crop classification When there is a single crop category in the CDL, the window is considered to be a pure pixel window, and the location information of the CDL pure pixel window is obtained. In this research example, corn, soybean, double crop (winter wheat and soybean) and Pure pixel windows of four types of objects such as forests;

S2: Use the MRT (MODIS Reprojection Tool) tool to analyze the low spatial resolution MODIS data (MCD43A4 collection 6, the resolution is 500m, the time range is from January 1, 2013 to December 31, 2014, and the image row is h12v05) Perform preprocessing such as resampling and reprojection to make it have the same spatial resolution and projection information as the high spatial resolution Landsat data and crop classification map CDL data, and the resampling method is bilinear interpolation. , the reprojection is UTM (UNIVERSAL TRANSVERSE MERCARTOR GRIDSYSTEM) projection, and the MODIS data and Landsat data are subjected to automatic geometric registration (automatic image matching method), and then the MODIS time series vegetation is obtained by calculation Exponential NDVI data; meanwhile, the high spatial resolution Landsat data used in this study contains 16 periods of Landsat7 data and 18 periods of Landsat8 data available in 2013 and 2014, with row and column numbers 33 and 14; There is an obvious banding problem, so the GNSPI (Geostatistical Neighborhood Similar PixelInterpolator) interpolation method is used to interpolate the Landsat7 images with band gaps to obtain data that completely covers the study area, and finally the Landsat time series NDVI data is obtained by calculation;

S3: Using the position information of the CDL pure pixel window obtained in step S1, select the position corresponding to the partial pixel area (16*16 pixels) in the middle of the CDL pure pixel window on the MODIS time series vegetation index NDVI data , and extract the average NDVI value of the pixels contained in the position. Therefore, according to the number of pure pixel windows of different types of crops, multiple MODIS time series NDVI curves of different types of crops can be obtained. In this research example, 23 corn curves, 1 soybean curve, 5 double crop (winter wheat and soybean) curves, and 26 forest curves were obtained in 2013, and 24 corn curves, 3 soybean curves, 10 There are two double crop (winter wheat and soybean) curves and 23 forest curves; since some of the curves still have outliers, this example uses linear interpolation to eliminate the outliers, and then uses these curves as MODIS crop reference curves for different crops. A MODIS crop reference curve sample library, Figure 1 is the MODIS reference curve diagram of four typical features in 2014;

S4: Using the original value of the MODIS time series vegetation index NDVI data calculated in step S2, the time series vegetation index NDVI curve of each pixel of the MODIS image is obtained, which is called the MODIS local curve of the pixel. Since the Landsat image pixels and the corresponding MDOSI image pixels may both be pure pixels, the NDVI curves of the two have a good correlation, or both are mixed pixels but may still have a good correlation. , so this study comprehensively considered the MODIS local curve and the MODIS crop reference curve sample library obtained in step S3, and combined the two as a candidate curve to participate in the curve matching;

S5: Using the Landsat time series NDVI data calculated in step S2, analyze the correlation index between all image pixels of the data and all the curves in the MODIS crop reference curve sample library obtained in step S3 and the MODIS local curve obtained in step S4 R ² size, take the curve with the largest correlation index R ² as the best MODIS reference curve of the image pixel of Landsat time series NDVI data, and then use formula (1) to further fit the Landsat time series NDVI data and MODIS in step S2 The best reference curve is obtained to obtain the Landsat fitting curve.

L(x)=a*M(x+x0)+b (1)

where x represents a certain day of the year, L(x) is the Landasat time series NDVI curve function, M(x) is the MODIS optimal reference curve function, and x0 represents the possible time offset between the two curves The range is between ±30 days, a and b are the fitting parameters calculated by the least squares method;

S6: By calculating the correlation index R 2 between the pixel values of the Landsat time series NDVI image in step S2 and the pixel values of all candidate curves of the corresponding MODIS, analyze the correlation index R ² between the ^two and select the best MODIS Reference curve. For each MODIS crop reference curve, there are 60 correlation indices R ² , whose value range is concentrated between 0.11 and 0.92, and then select the candidate curve with the largest correlation index as the curve to be matched for the MODIS crop reference curve, and then compare The magnitude of the correlation index between the curves to be matched and Landsat time series NDVI data corresponding to all MODIS crop reference curves, as well as the magnitude of the correlation index between the MODIS local curve and Landsat time series NDVI data, the values range from 0.78 to 0.92 between. If the MODIS local curve correlation index R ² value is the largest, the MODIS local curve is directly used as the best reference curve; if the correlation index R ² maximum value appears in the curve to be matched in the MODIS crop reference curve, the correlation index R of all candidate curves is calculated. ² , select the candidate curve higher than the average value of R ² and average it again, and take the final average curve as the best reference curve. As shown in Figure 2, the black dots are the Landsat time series NDVI values calculated in step S2, the curve (1) is the MODIS reference curve obtained in step S3, and the curve (2) is the MODIS maximum obtained by fitting in steps S5 and S6. best reference curve;

S7: Use steps S5 and S6 to confirm the best reference curve. At this time, the best time offset x0 is 0, and the corresponding fitting parameters a and b are obtained at the same time, and then the MODIS best reference curve and the corresponding fitting parameters are obtained. Fit the parameters to get the Landsat initial fitting curve. As shown in Figure 2, curve (3) is the initial fitting curve. At this time, the initial fitting curve has a good fitting relationship with the Landsat time series NDVI values of the target pixel area, but there are still some original Landsat NDVI values that cannot be Fits well to the initial fitting curve;

S8: Calculate the difference between the NDVI pixel value (clear sky pixel) of the original Landsat time series and the NDVI fitting value of the corresponding pixel on the Landast initial fitting curve obtained in step S7, and use the linear interpolation method to fit the linear interpolation result. Combined to the initial fitting curve to reduce the difference between the two, and finally use the Gaussian function in the IDL programming platform to further smooth the optimized fitting curve, so that the curve contains all the available original Landsat time as much as possible The sequence NDVI value is called the final fitting curve. As shown in Figure 2, curve (4) is the Landsat final fitting curve. At this time, the Landsat time series NDVI value and the Landsat final fitting curve have a good fit relation.

In order to compare the effect of the data reconstruction method based on the crop reference curve in this study, this paper also uses the classic spatiotemporal fusion algorithm STARFM to carry out comparative experiments, which requires multiple pairs of paired images acquired under clear sky conditions to be fitted. Figure 3 shows the NDVI fitting curves of three different crops based on the crop reference curve and the STARFM algorithm in the study area in 2013, in which curve (3) is the MODIS reference curve, and the black big dots are the available data for the STARFM algorithm. Cloud paired images, curve (2) is the fitting curve based on the STARFM algorithm, black dots are all available Landsat images based on the crop reference curve reconstruction method used in this study, and curve (1) is the crop reference curve reconstruction method The resulting fitted curve. It can be found that the coarse-resolution MODIS image is affected by resolution and mixed pixels, while the classical algorithm STARFM is affected by cloud coverage, which leads to the lack of available paired images. phenological change characteristics. The data reconstruction method used in this study can greatly avoid the influence of cloud cover, and the obtained high temporal and spatial resolution data can well reflect the key phenological changes of different crops.

To sum up, this study uses the known crop classification data to set the pure pixel area, and uses the concept of pure pixel to obtain the crop reference curve based on low spatial resolution, and then fits the high spatial resolution through translation stretching and other fitting methods. The image and crop reference curve are fitted, converted and adjusted by secondary optimization, and the temporal dimension information of the low spatial resolution crop reference curve is well converted to the high spatial resolution image, and the high spatial and temporal resolution based on the crop reference curve is realized. The fusion and reconstruction of vegetation index data can well solve the interference and limitation of the original data spatiotemporal fusion algorithm for mixed pixels and cloud coverage and other factors for data fusion in areas with high spatial heterogeneity. System remote sensing monitoring is of great significance.

It should be understood that, for those skilled in the art, improvements or changes can be made according to the above description, and all these improvements and changes should fall within the protection scope of the appended claims of the present invention.

Claims (4)

1. A method for carrying out high spatial-temporal resolution vegetation index data fusion based on a crop reference curve is characterized by comprising the first step of constructing a crop reference curve sample library, obtaining a plurality of MODIS reference curves of different crops based on pure pixels by utilizing crop category information and spatial positions of a crop classification diagram, forming an MODIS crop reference curve sample library, and simultaneously combining MODIS local curves to be used as MODIS candidate curves to participate in curve matching, the second step of determining an initial fitting curve, fitting and matching pixel values of L andsat time sequence NDVI images and image pixel values of all corresponding positions of all curves in the MODIS candidate curves, finding a candidate curve with the highest fitting degree as an optimal reference curve based on the MODIS crop reference curve, fitting L andsat time sequence NDVI images by utilizing fitting parameters of the optimal reference curve to obtain L andsat initial fitting curves, the third step of determining a final fitting curve, fitting and smoothing the final fitting curve to obtain MODisat time resolution and space resolution after fitting of a secondary vegetation index with high spatial resolution, so that the MODIS L and space resolution are obtained;

the first step comprises the following steps:

step A: setting a moving window with the size of 25 × 25 pixels by utilizing the crop category information in the crop classification map, and when more than 95% of pixels in the window belong to a single crop category in the crop classification map through the moving pixel window, considering the window pixel as a pure pixel and acquiring the position information of the pure pixel window of the crop classification map;

b, resampling and reprojection preprocessing are carried out on the MODIS data with low spatial resolution, so that the MODIS data has the same spatial resolution and projection information as the L andsat data with high spatial resolution and the crop classification map data, geometric registration is carried out on the MODIS data and the L andsat data, and then the data are calculated to obtain the NDVI data of the MODIS time sequence vegetation index;

and C: b, selecting positions corresponding to part of pixel areas in the middle of the window on the MODIS time sequence vegetation index NDVI data by using the position information of the pure pixel windows of the crop classification images obtained in the step A, and extracting to obtain NDVI average values of pixels contained in the positions to obtain a plurality of MODIS time sequence NDVI curves of different types of crops; for the abnormal values of part of curves, a linear interpolation method is adopted to eliminate the abnormal values, and the curves are used as MODIS reference curves of different ground objects to construct a sample library;

the second step comprises the following steps:

step D: obtaining a time-series vegetation index NDVI curve of each pixel of the MODIS image by using the original value of the MODIS time-series vegetation index NDVI data obtained in the step B, referring to the MODIS local curve of the pixel, then comprehensively considering the MODIS local curve and the MODIS reference curve sample library obtained in the step C, and combining the MODIS local curve and the MODIS reference curve sample library together to be used as an MODIS candidate curve to participate in curve matching;

step E, utilizing L andsat time sequence NDVI data obtained by calculation in the step B to analyze the correlation indexes R of all image pixels of the data and all curves in the MODIS crop reference curve sample library obtained in the step C²The size of the data and the correlation index R of all image pixels of the L andsat time sequence NDVI data and the MODIS local curve obtained in the step D²Size, correlation index R²The maximum curve is used as an MODIS optimal reference curve of an image pixel of L andsat time series NDVI data, and then L andsat time series NDVI data in the step B and the MODIS optimal reference curve are further fitted through a formula (1) to obtain a L andsat fitting curve;

L(x)＝a*M(x+x0)+b (1)

wherein x represents a day of the year, L (x) is a L andasat time series NDVI curve function, M (x) is an MODIS optimal reference curve function, x0 represents the offset in time which may exist between two curves, the offset is within +/-30 days, and a and b are fitting parameters calculated by using a least square method;

step F, calculating a correlation index R between the L andsat time sequence NDVI image pixel value in the step B and all candidate curve pixel values of the corresponding position MODIS²Analysis of the correlation index R between the two²Selecting an optimal MODIS reference curve;

and G, confirming an optimal reference curve by using the steps E and F, wherein the optimal time offset x0 is 0, obtaining corresponding fitting parameters a and b simultaneously, and then obtaining a L andsat initial fitting curve by using the MODIS optimal reference curve and the corresponding fitting parameters.

2. The method of claim 1, wherein step three comprises the steps of step H of calculating the difference between the original L andsat time series NDVI pixel values and the corresponding pixel NDVI fitted values on the L andast initial fitted curve obtained in step G, fitting the linear interpolation result to the initial fitted curve by linear interpolation to reduce the difference between the two, and finally smoothing the optimized fitted curve by using a Gaussian function in the ID L programming platform to make the curve contain all available original L andsat time series NDVI values as much as possible, and then calling the curve as a final fitted curve.

3. The method as claimed in claim 2, wherein in step F, there are 60 correlation indices R for each MODIS crop reference curve²The value range is concentrated between 0.11 and 0.92, then the candidate curve with the maximum correlation index is selected as the curve to be matched of the MODIS crop reference curve, the correlation index between the curve to be matched corresponding to all the MODIS crop reference curves and the L andsat time sequence NDVI data is compared, the value range of the correlation index between the MODIS local curve and the L andsat time sequence NDVI data is concentrated between 0.78 and 0.92, if the correlation index R of the MODIS local curve is concentrated between 0.11 and 0.92²If the value is maximum, directly taking the MODIS local curve as the optimal reference curve; if the correlation index R²If the maximum value appears in the curve to be matched of the MODIS crop reference curve, calculating all the candidate curve correlation indexes R²Is selected to be higher than R²And averaging the candidate curves of the average value again, and taking the finally obtained average curve as an optimal reference curve.

4. The method according to claim 1, wherein the MODIS reference curves of different land features are obtained by using the selected pure pixel areas with higher classification precision results in the known crop classification map data, and the specific setting parameters in the step A are empirical parameters according to a plurality of experiments.

Images