A method for processing time series of Sentinel-1 VH polarized signals based on multi-year spatio-temporal information

Abstract

The application belongs to the technical field of remote sensing monitoring, and particularly provides a Sentinel-1 VH polarization signal time sequence processing method based on multi-year space-time information, which is used to solve the problems that the existing radar time sequence processing method is interfered by noise in the process of performing phenology detection, and the traditional time domain filtering method is directly used on the radar time sequence curve. The application fully considers the autocorrelation in the VH time sequence and the regularity of agricultural activities, uses multiple sampling points in the spatial window to obtain the annual time sequence curve after being averaged for multiple years, establishes a regression relationship with the real observed curve of the current year after macroscopic stretching and translation, and then obtains the annual VH polarization signal time sequence observation curve of the target pixel of the detected crop in the current year through local stretching of the curve, which is used for subsequent phenology detection and effectively improves the accuracy of phenology detection.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of remote sensing monitoring, and particularly provides a Sentinel-1 VH polarization signal time series processing method based on multi-year spatio-temporal information. BACKGROUND

[0002] Sentinel-1 synthetic aperture radar data is a series of earth observation satellites launched by the European Space Agency (ESA), which uses C-band (wavelength about 5.6 cm) for imaging, and can obtain high-resolution radar image data of the earth's surface under all-weather and all-day conditions; at present, Sentinel-1 synthetic aperture radar data has been widely used in various ground monitoring, including: land cover type, land subsidence, and vegetation growth phenology research; especially in the field of crop phenology detection, based on its characteristics of not being affected by clouds and fog, it can detect crop growth conditions multiple times and obtain phenology information; although Sentinel-1 radar data can not be affected by cloud cover and other conditions, but based on its scattering intensity data, when conducting time series analysis research, it is more easily affected by environmental noise and crop growth state, resulting in abnormal environmental noise.

[0003] At present, in order to detect crop phenology based on Sentinel-1 VH polarization time series and eliminate environmental factor noise as much as possible, researchers have proposed a variety of methods. Early methods are based on function fitting, such methods use different mathematical functions to simulate the annual trajectory of radar time series, smooth the original time series, and define a specific threshold to mark the key time nodes of crop growth. In recent years, in order to fully utilize synthetic aperture radar time series data for phenology detection, geometric features based on multiple different types of polarization data have been proposed to find the correlation between the actual phenology stage of crops (such as sowing, germination, and jointing) and the geometric properties of time series (such as the mutation point of slope or local extreme value). After smoothing, the geometric feature points can be directly extracted and correspond to the known phenology stage. Such methods can directly capture the turning points of crop phenology and reduce errors caused by improper function fitting.

[0004] Therefore, it can be seen that most of the existing methods focus on directly obtaining the phenological phase through the geometric characteristics of the curve itself, and only do simple time domain filtering or function fitting smoothing for related noise; however, these methods can only process part of the high-frequency noise, and in the analysis of time series data, there is no good solution for low-frequency noise caused by environmental influences; for example, soil characteristics, extreme environmental influences, weed growth, etc., which will have a significant impact on the geometric characteristics of radar time series data, thereby interfering with the phenological detection results, especially when performing phenological detection in a large area, due to the different crop planting periods and growth environments, such influences are further highlighted. Therefore, it is urgent to further improve the Sentinel-1 VH polarization signal time series processing method. SUMMARY

[0005] The present application aims to provide a Sentinel-1 VH polarization signal time series processing method based on multi-year spatiotemporal information, to solve the problem that the existing radar time series processing method is disturbed by noise in the process of phenological detection, and the deficiency of traditional time domain filtering method directly used on radar time series curve. The present application fully considers the autocorrelation in the VH time series and the regularity of agricultural activities, uses multiple sampling points in the spatial window to obtain the annual time series curve after multi-year averaging, establishes a regression relationship with the real observed curve after macroscopic stretching and translation, and obtains the optimal curve after local stretching of the curve.

[0006] To achieve the above object, the technical scheme adopted by the present application is:

[0007] A Sentinel-1 VH polarization signal time series processing method based on multi-year spatiotemporal information, characterized in that it comprises the following steps:

[0008] S1, obtaining the annual VH time series curve of the target pixel of the detected crop in the current year, setting a spatial window centered on the target pixel, sampling other pixels of the detected crop in the spatial window, and obtaining the average VH time series curve of each sampling pixel;

[0009] S2, performing the same mean filtering on the annual VH time series curve of the target pixel and the average VH time series curve of all sampling pixels;

[0010] S3, for each sampling pixel, shift and stretch the average VH time series curve of the filtered sampling pixel, establish the regression relationship between the annual VH time series curve of the target pixel and each shifted and stretched average VH time series curve by GMFR algorithm one by one, generate the predicted VH time series curve of the target pixel in the current year according to the regression relationship, calculate the MSD value with the filtered annual VH time series curve of the target pixel, and select the predicted VH time series curve with the minimum MSD value as the optimal predicted time series curve; in all sampling pixels, select the optimal predicted time series curve with the minimum MSD value as the optimal predicted time series curve;

[0011] S4, segment the optimal predicted curve with the maximum value as the center, stretch and match the front segment and the rear segment respectively, and combine the stretching and matching results of the front segment and the rear segment to form the final predicted time series curve of the target pixel of the detected crop in the current year.

[0012] Further, in step S1, the sampling process is: for each sampling pixel, obtain the annual VH time series curve of the sampling pixel in each year of planting the detected crop, calculate the average value to obtain the average VH time series curve:

[0013]

[0014] Wherein, t represents the time variable, N i represents the sampling year number of the i th sampling pixel, VH j (t) represents the VH time series curve of the j th sampling year of the i th sampling pixel, VH i (t) represents the average VH time series curve of the i th sampling pixel, and I represents the number of sampling pixels.

[0015] Further, in step S3, the calculation process of shift and stretch is:

[0016] T=xscale×(t+tshift)

[0017] Wherein, T represents the time variable after shift and stretch processing, t represents the original time variable; tshift represents the time scale shift parameter, the value range of tshift is [-30, 30]; xscale represents the time scale stretch parameter, the value range of xscale is [0.9, 1.1].

[0018] Further, in step S3, the generation process of the predicted VH time series curve is:

[0019] The filtered average VH time series curve of the i th sampling pixel is represented as VH i (t), according to the time scale shift parameter and the time scale stretch parameter, the average VH time series curve VH i"(t) and the average VH time series curve VH i "(T), T represents the time variable after translation and stretching processing;

[0020] The annual VH time series curve VH' of the target pixel is established by the GMFR algorithm tar "(t) and the average VH time series curve VH i "(T) of the target pixel in the current year, and the regression coefficients a and b are obtained, and the predicted VH time series curve VH pre "(t) of the target pixel in the current year, and the regression coefficients a and b are obtained, and the predicted VH time series curve VH pre "(t) = a x VH"(T) + b.

[0021] Further, in step S3, the calculation process of MSD value is:

[0022] The annual VH time series curve of the filtered target pixel is represented as VH' tar "(t), and the predicted VH time series curve VH pre "(t) and the annual VH time series curve VH' tar "(t) is calculated:

[0023]

[0024] Where Year represents the total length of the VH time series curve.

[0025] Further, in step S4, the calculation process of stretching is:

[0026] T' = xscale x |t-t max |

[0027] Where T' represents the time variable after stretching processing, t max represents the maximum value time of the optimal prediction curve; xscale represents the time scale stretching parameter, and the value range is [0.9, 1.1].

[0028] Further, in step S4, the specific process of stretching matching is:

[0029] The annual VH time series curve of the filtered target pixel is represented as VH' tar "(t), and the same segmentation is performed on the annual VH time series curve VH' tar "(t) according to the maximum value time of the optimal prediction curve;

[0030] For the front segment or the rear segment of the optimal prediction curve, the MSD value is calculated for the corresponding segment of the stretched segmented curve and the annual VH time series curve VH' tar "(t), and the stretching result with the smallest MSD value is taken as the segmented stretching matching result.

[0031] Based on the above technical scheme, the application has the beneficial effects of:

[0032] The application provides a Sentinel-1 VH polarization signal time sequence processing method based on multi-year space-time information, which sets a space window around a target pixel, samples other pixels of the detected crop in the space window, acquires a multi-year VH polarization signal time sequence of the sampled pixels, and fully integrates the time and space information in the space window around the target pixel; the application adopts an average calculation method of the multi-year VH polarization signal time sequence of each sampled pixel, generates a multi-year average VH polarization signal time sequence for each sampled pixel, and smoothes the complex environmental noise contained in the VH polarization signal time sequence; the application adopts overall translation and stretching calculation of the multi-year average VH polarization signal time sequence to avoid macroscopic differences between the multi-year average VH polarization signal time sequence and an actual observed VH polarization signal time sequence of the current year; the application adopts a GMFR algorithm to establish a regression relationship between the macroscopically stretched and translated multi-year average VH polarization signal time sequence and the actual observed VH polarization signal time sequence of the current year, obtains a predicted actual observed VH polarization signal time sequence of the current year according to the regression relationship, and excludes the influence of part of abnormal outlier VH values; the application adopts a method of calculating MSD values between the regression predicted actual observed VH polarization signal time sequence of the current year and the actual observed VH polarization signal time sequence of the current year to obtain an optimal regression predicted actual observed VH polarization signal time sequence with the minimum MSD value; the application adopts a method of stretching the left and right sides with the VH maximum value as the center point to further make the optimal regression predicted actual observed VH polarization signal time sequence of the current year more capable of reflecting the regularity between crop growth and VH polarization signal change in the local; finally, the application can accurately restore the change rule of the VH polarization signal time sequence based on multi-year space-time information, exclude the noise influence of complex environmental factors on the VH polarization signal, and effectively identify the VH polarization signal mode related to the crop phenology cycle. The VH polarization signal time sequence processed based on the application provides a higher data quality basis for crop phenology monitoring and prediction, and can further improve the phenology detection accuracy based on the VH polarization signal time sequence. BRIEF DESCRIPTION OF DRAWINGS

[0033] Figure 1 The algorithm flowchart of the Sentinel-1 VH polarization signal time sequence processing method based on multi-year space-time information;

[0034] Figure 2 The box plot of all samples between the Phenocam ground phenology observation values and the corn phenology estimated values in the embodiment and the comparative example of the application.

[0035] Figure 3This is a box plot of all samples between Phenocam ground phenological observations and soybean phenological estimates in the embodiments and comparative examples of the present invention.

[0036] Figure 4 This is a scatter plot comparing the corn phenological detection results and phenological labels of various states in the Corn Belt in the embodiments and comparative examples of the present invention.

[0037] Figure 5 This is a scatter plot comparing soybean phenological detection results with phenological labels in various states of the Corn Belt in the embodiments and comparative examples of the present invention.

[0038] Figure 6 The figures show the changes in the proportion of phenological stages of maize in Nebraska over time in the embodiments and comparative examples of this invention.

[0039] Figure 7 The figures show the changes in the proportion of soybean phenological stages in Nebraska over time in the embodiments and comparative examples of the present invention.

[0040] Figure 8 This is a scatter plot comparing the rapeseed phenological detection results and phenological labels in Sichuan Province in the embodiments and comparative examples of the present invention. Detailed Implementation

[0041] To make the objectives, technical solutions, and beneficial effects of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings and embodiments.

[0042] This invention provides a time series processing method for Sentinel-1 VH polarization signals based on multi-year spatiotemporal information. Using multi-year averaged VH polarization time series data, and through macroscopic and local stretching and translation calculations, a processed time series is generated. This method reduces noise from environmental factors, restores the correspondence between vegetation growth phenology and the time series, and can be used for phenological detection. The specific process is as follows: Figure 1 As shown, it includes the following steps:

[0043] S1. Obtain the annual VH time-series curve of the target pixel of the current year for the detected crop. tar (t), set a spatial window centered on the target pixel, and sample other pixels of the detected crop within the spatial window; for each sampled pixel: obtain the annual VH time series curve of the sampled pixel in each year when the detected crop is planted, calculate the average value, and generate an average VH time series curve based on VH data over many years;

[0044] Specifically, the VH time series curve of the target pixel of the detected crop in the current year is obtained, a spatial window with a side length of 3 km is generated outwardly from the target pixel, based on the crop classification result (such as the Cropland DataLayer data of the United States Department of Agriculture), other pixels of the detected crop in the spatial window are sampled to obtain the average VH time series curve of each sampling pixel, which is represented as:

[0045]

[0046] wherein t represents a time variable, N i represents the sampling years of the i th sampling pixel, VH j (t) represents the VH time series curve of the j th sampling year of the i th sampling pixel, VH' i (t) represents the average VH time series curve of the i th sampling pixel, and I represents the number of sampling pixels.

[0047] S2, the same mean filtering is performed on the VH time series curve of the target pixel and the average VH time series curve of all sampling pixels, the filtered VH time series curve of the target pixel is represented as VH' tar (t), and the filtered average VH time series curve of the sampling pixel is represented as VH i "(t)";

[0048] S3, the average VH time series curve of the sampling pixel after filtering is matched by translation and stretching, the average VH time series curve of the sampling pixel after translation and stretching is represented as VH i "(T); the regression relationship between the VH time series curve VH' tar (t) of the target pixel and each average VH time series curve VH i "(T) is established one by one through the GMFR algorithm, the influence of extreme outliers is excluded, the predicted VH time series curve VH pre (t) of the target pixel of the detected crop in the current year is generated according to the regression relationship, and the MSD value is calculated with the VH time series curve VH' tar (t) of the target pixel after filtering; for each sampling pixel, the predicted VH time series curve with the smallest MSD value is selected as the preferred predicted time series curve, that is, the translation and stretching matching of the current sampling pixel is completed; further, in all sampling pixels, the preferred predicted time series curve with the smallest MSD value is selected as the optimal predicted time series curve.

[0049] Specifically, the translation and stretching processing is represented as:

[0050] T = xscale × (t + tshift)

[0051] Wherein, T represents the time variable after translation and stretching processing, t represents the original time variable; tshift represents the time scale translation parameter, the value range of tshift is [-30, 30], and the step is 6 days; xscale represents the time scale stretching parameter, the value range of xscale is [0.9, 1.1], and the step is 0.05;

[0052] According to the time scale translation parameter tshift and the time scale stretching parameter xscale, the average VH time sequence curve VH i (t) of each sampling pixel is translated and stretched to obtain the corresponding average VH time sequence curve VH i ″(T);

[0053] For each sampling pixel: the annual VH time sequence curve VH′ tar (t) of the target pixel is established by the GMFR algorithm, and the regression relationship between VH i ″(T) of the sampling pixel is obtained, and the regression coefficient is obtained, including the slope parameter a and the intercept parameter b, and the predicted VH time sequence curve VH pre (t) of the target pixel of the detected crop in the current year is generated according to a and b;

[0054] VH pre (t)=a×VH″(T)+b

[0055] After the regression relationship is established, the predicted VH time sequence data of the target pixel of the detected crop in the current year is obtained by the average time sequence data prediction;

[0056] The predicted VH time sequence curve VH pre (t) is compared with the annual VH time sequence curve VH′ tar (t) of the filtered target pixel, and the MSD value is calculated:

[0057]

[0058] Wherein, Year represents the total length of the VH time sequence curve;

[0059] S4, the optimal prediction curve is segmented with the maximum value as the center, and the front segment and the rear segment are stretched and matched respectively; the stretched and matched results of the front segment and the rear segment are combined to form the final prediction time sequence curve of the target pixel of the detected crop in the current year, which is used as the data basis for subsequent phenology detection;

[0060] Specifically, for the phenology detection, the optimal prediction curve is segmented with the maximum value as the center, into a front segment and a rear segment, the front segment can be regarded as the time sequence before the crop growing season, and the rear segment can be regarded as the time sequence after the crop growing season, and the front segment and the rear segment are stretched and matched respectively, which is specifically represented as:

[0061] T' = xscale * |t - t max |

[0062] Wherein, T' represents time variable after stretching processing, t max represents the maximum value time of the optimal prediction curve; the value range of the time scale stretching parameter xscale is [0.9, 1.1], and the step is 0.05;

[0063] Similar to the formula matching translation and stretching in S2, the stretching in S2 takes the starting time as the starting point, and the optimal prediction curve in S4 takes the maximum value as the center to stretch;

[0064] According to the maximum value time t max The annual VH time series curve VH' of the filtered target pixel tar (t) is segmented in the same way, and on this basis, for the front segment or the rear segment of the optimal prediction curve, the MSD value is calculated for the corresponding segment of the stretched segmented curve and the annual VH time series curve VH' of the filtered target pixel tar (t), and the stretching result with the minimum MSD value is taken as the segmented stretching matching result; the stretching matching results of the front segment and the rear segment are combined to obtain the annual VH polarized signal time series observation curve of the target pixel in the current year of the crop detected after the processing of the Sentinel-1 VH polarized signal time series based on multi-year spatiotemporal information in the application.

[0065] The following takes American corn, soybean and Chinese rape as examples to explain the beneficial effects of the Sentinel-1 VH polarized signal time series processing method based on multi-year spatiotemporal information in the application in detail.

[0066] The Sentinel-1 VH polarized signal time series processing method based on multi-year spatiotemporal information in the application is used for phenological detection experiment test in the corn and soybean planting areas in the United States and the rape planting areas in Sichuan Province, China, and the LOESS filtering group is taken as a comparative example, and the time series data processed by LOESS filtering is compared and analyzed; the three test scenarios have different characteristics respectively.

[0067] The first test scenario uses American corn and soybean PhenoCam phenological data, and the characteristic is to obtain radar time series data based on single phenological site positioning, and then perform phenological detection; the digital photos provided by the Phenocam website are used for phenological period interpretation as label data;

[0068] The second test scenario is corn and soybean in 12 states of the US Corn Belt, characterized by regional phenology detection in a large spatial range; the labeled data is the crop area reaching a certain phenological stage in the survey provided in the Crop Progress Report of each state provided by the US Department of Agriculture, and the US corn and soybean are generally planted around May;

[0069] The third test scenario is Sichuan Province, China, which is a cloudy and foggy area, and it is very difficult to obtain optical data in the area. The availability of Sentinel-1 radar data for phenology detection in the area and the effectiveness of the processing algorithm can be verified. The test scenario uses a method of taking photos by farmers and interpreting the photos by agricultural experts to obtain phenological stage labels.

[0070] For quantitative evaluation, two methods, SMF-S algorithm and DTW algorithm, are used for phenology detection, and the phenology detection results are compared with the labeled results to calculate the MAE value.

[0071] The quantitative evaluation of the US corn and soybean Phenocam sites is performed, and the test method is based on the time series data processed by the method and the time series data filtered by the LOESS algorithm. The SMF-S algorithm and the DTW algorithm are used for phenology detection, and the phenology detection accuracy is compared. As shown in Tables 1 and 2, the phenology detection accuracy of the three phenology detection algorithms based on the time series data processed by the method and the time series data filtered by the LOESS algorithm on corn and soybean crops is shown.

[0072] Table 1

[0073] phenophase V2.5 V7 VT R6 SMF-S (LOESS filtered) 20.19 15.58 18.91 18.90 SMF-S (invention) 10.98 12.36 11.87 13.71 DTW (LOESS filtered) 19.82 19.36 28.54 32.94 DTW (invention) 19.05 22.02 18.76 27.26

[0074] Table 2

[0075] phenophase VC V2 R7 SMF-S (LOESS filtered) 16.13 14.29 16.97 SMF-S (invention) 11.03 8.63 15.59 DTW (LOESS filtered) 19.50 19.37 17.31 DTW (invention) 16.69 11.24 14.62

[0076] As Figure 2 , Figure 3As shown, the box plot of all samples between the ground phenology observation value of Phenocam and the phenology estimation value detected by the two algorithms based on two types of data is shown; by averaging the MAE of all phenological periods, for corn, the MAE of the phenology detection of the two algorithms after the processing of the method is 12.23, 21.77, and for soybeans, the corresponding MAE is 11.75, 14.19; after the phenology detection of the time series data processed only by LOESS filtering, for corn, the MAE of the phenology detection of the two algorithms is 18.40, 24.92, and for soybeans, the MAE of the phenology detection is 15.80, 18.73; for the two algorithms, the data processed by the method has obvious advantages in the field of phenology detection, especially when the SMF-S algorithm detects the phenology of the data processed by the method, the final accuracy is improved by 4 to 6 days (corn: 18.40 days vs. 12.23 days; soybeans: 15.80 days vs. 11.75 days) compared with the original data. From the results, in the phenology detection of the two algorithms, the data processed by the method has stronger stability, compared with the original data, the algorithm reduces the influence of low-frequency environmental noise, and restores the regularity of the VH time series and the crop phenological growth cycle; at the same time, under the condition that the Phenocam sites are distributed throughout the United States, the processing of all sample time series is relatively stable in a larger complex spatial range, and the obstacles caused by multiple environmental factors to VH time series analysis are suppressed.

[0077] For the phenology detection of corn and soybeans in a large spatial range of 12 states in the United States corn belt, based on the Cropland Data Layer data released by the United States Department of Agriculture, 200 crop pixels in each state range are selected as sample data; LOESS filtering and the method are used to process all sample data, and the data are subjected to phenology detection. Based on the Crop Progress Report phenology report (CPR) provided by the United States Department of Agriculture, the date when the area of the crop phenology reaches 50% in the state agricultural survey is taken as the label of the crop phenological period in the state, and the final result of the phenology detection in each state is compared with the result of the phenological period. As shown in Table 1, the MAE of the phenology detection of the two algorithms after the processing of the method is 12.23, 21.77 for corn, and 11.75, 14.19 for soybeans, respectively, which is significantly lower than the MAE of the phenology detection of the time series data processed only by LOESS filtering, which is 18.40, 24.92 for corn and 15.80, 18.73 for soybeans, respectively. The method has obvious advantages in the field of phenology detection, especially when the SMF-S algorithm detects the phenology of the data processed by the method, the final accuracy is improved by 4 to 6 days (corn: 18.40 days vs. 12.23 days; soybeans: 15.80 days vs. 11.75 days) compared with the original data. Figure 4 、 Figure 5As shown in the figure, the scatter plot of the corn and soybean phenology detection results of each state compared with the phenology label provided by the CPR is shown, and it can be seen from the figure that the data processed by the method of the application is better than the data processed only by the LOESS filtering algorithm in both algorithms; in addition, the estimation of all samples by the method of the application is relatively stable, without outliers, i.e. without obvious larger errors, indicating that the application has robustness; and some points that are obviously outliers in the original data are returned to the normal range after being processed by the method of the application; the MAE of the phenology detection is greatly improved, and for the SMF-S algorithm, the average MAE of the phenology detection of soybeans and corn is reduced by about 3-4 days before and after processing by the method (soybeans: 11.32 days vs. 8.56 days; corn: 11.10 days vs. 7.61 days). As shown in Figure 6 、 Figure 7 As shown in the figure, the change curves of the proportion of the crop area reaching the phenological stage recorded by the CPR data of the corn and soybean crops in Nebraska, the proportion of the sample data reaching the phenological stage detected by the SMF-S algorithm processed by the method of the application and the proportion of the sample data reaching the phenological stage processed by the LOESS filtering algorithm are compared, and the change trends of the three curves are compared. It is found that the proportion of the sample data reaching the phenological stage detected by the SMF-S algorithm processed by the method of the application is more consistent with the change trend of the proportion of the area reaching the corresponding phenological stage investigated by the CPR data, which can reflect the growth phenological cycle of the crops in Nebraska.

[0078] The phenology detection of part of the rape in Sichuan Province of China is based on the interpretation of the photos taken by farmers in the field to obtain the growth phenology label of the rape in the located field. The MAE of the phenology detection of the sample data processed by the method of the application and the data processed only by the LOESS filtering algorithm under the SMF-S and DTW two phenology detection algorithms is compared. The LOESS algorithm filters out part of the low-frequency noise, effectively smoothing the time series data, but does not better remove the influence of environmental factors on the crop phenology growth signal, and the method of the application retains the regularity between the crop phenology growth rule and the SAR scattering intensity signal on the basis of effectively removing the noise. As shown in Table 3, the accuracy of the phenology detection of the two algorithms under two types of data processed in different ways is shown.

[0079] Table 3

[0080] phenophase rosette stage first flowering stage full flowering stage podding stage mature stage harvest stage SMF-S (LOESS filtered) 19.50 24.33 19.38 18.63 26.25 32.11 SMF-S (invention) 11.00 11.78 9.44 8.56 8.44 14.00 DTW (LOESS filtered) 34.39 25.39 17.28 14.83 17.11 7.78 DTW (invention) 48.44 19.39 8.56 12.11 7.61 6.33

[0081] The results show that the complex terrain and climate environment have a great influence on the accuracy of phenology detection, and the DTW overall matching algorithm has a large deviation, but the method still has a certain improvement on the final accuracy of phenology detection; the local matching mode of SMF-S makes it less affected by environmental noise, but the MAE of phenology detection of the data processed by LOESS filtering is still generally greater than 18 days, while the time series data processed by the method has an MAE of about 10 days, and the accuracy of the podding period and the mature period is optimal (podding period: 18.63 days vs. 8.56 days; mature period: 26.25 days vs. 8.44 days). As shown in FIG. 8, the scatter plot of the rapeseed phenology detection results in Sichuan Province compared with the field monitoring phenology label is shown, and it can be seen that the data processed by the method can significantly improve the accuracy in the field of phenology detection, and even in the complex mountainous area of Sichuan Province dominated by smallholder economy, there is still a significant denoising effect. Figure 8

[0082] The above is only a specific implementation of the present application, and any feature disclosed in the specification can be replaced by other equivalent or similar purpose alternative features unless specifically stated; all features disclosed, or steps in all methods or processes, except for mutually exclusive features and / or steps, can be combined in any way.​

Claims

1.A method for processing time series of Sentinel-1 VH polarized signals based on multi-year spatio-temporal information, characterized in that, The method comprises the following steps: S1, obtaining a full-year VH time series curve of a target pixel of a detected crop in a current year, setting a spatial window centered on the target pixel, sampling other pixels of the detected crop in the spatial window, and obtaining an average VH time series curve of each sampling pixel; S2, performing the same mean filtering on the full-year VH time series curve of the target pixel and the average VH time series curves of all sampling pixels; S3, for each sampling pixel, performing translation and stretching matching on the filtered average VH time series curve of the sampling pixel, establishing a regression relationship between the full-year VH time series curve of the target pixel and each average VH time series curve after translation and stretching through the GMFR algorithm, generating a predicted VH time series curve of the target pixel in the current year according to the regression relationship, calculating an MSD value of the filtered full-year VH time series curve of the target pixel, and selecting a predicted VH time series curve with the smallest MSD value as an optimal predicted time series curve; S4, segmenting the optimal predicted curve with the maximum value as the center, performing stretching matching on the front segment and the rear segment respectively, and combining the stretching matching results of the front segment and the rear segment to obtain a final predicted time series curve of the target pixel of the detected crop in the current year. 2.The method according to claim 1, wherein, In step S1, the sampling process is as follows: for each sampling pixel, obtaining a full-year VH time series curve of the sampling pixel in each year of the detected crop, calculating the average value, and obtaining an average VH time series curve. where t represents a time variable, N i denotes the number of sampling years of the i-th sampling pixel, VH j (t) represents the VH time series curve of the j-th sampling year of the i-th sampling pixel, VH' i (t) represents the average VH time series curve of the i-th sampling pixel, I represents the number of sampling pixels. 3.The method according to claim 1, wherein, In step S3, the calculation process of translation and stretching is as follows: T = xscale × (t + tshift) Wherein, T represents the time variable after translation and stretching processing, t represents the original time variable; tshift represents the time scale translation parameter, and the value range of tshift is [-30, 30]; xscale represents the time scale stretching parameter, and the value range of xscale is [0.9, 1.1]. 4.The method according to claim 1, wherein, In step S3, the generation process of the predicted VH time series curve is as follows: The average VH time series curve of the filtered i-th sampling pixel is expressed as VH" i (t), according to the time scale translation parameter and the time scale stretching parameter, the average VH time series curve VH" i (t) is translated and stretched to obtain the corresponding average VH time series curve VH" i (T), T represents the time variable after the translation and stretching processing. The GMFR algorithm is used to establish the VH time series curve VH' of the target pixel tar (t) of the average VH time series curve VH" of the sampling pixel i The regression relationship of (T) is obtained to obtain the regression coefficients α and β, and a predicted VH time series curve VH of the target pixel in the current year of the detected crop is generated pre (t) of the VH time series curve VH pre (t) = α × VH"(T) + β. 5.The method of claim 1, wherein, In step S3, the calculation process of the MSD value is as follows: The annual VH time series curve of the filtered target pixel is expressed as VH' tar (t), the predicted VH time series curve VH pre (t) is compared with the annual VH time series curve VH' tar (t) to calculate the MSD value: Wherein, Year represents the total length of the VH time series curve. 6.The method of claim 1, wherein, In step S4, the calculation process of stretching is as follows: T' = xscale x |t - t max | where T' represents the time variable after the stretching treatment, t max represents the maximum value of the optimal prediction curve, xscale represents the time scale stretching parameter, and the value range is [0.9, 1.1]. 7.The method according to claim 1, wherein, In step S4, the specific process of stretching matching is as follows: The annual VH time series curve of the filtered target pixel is expressed as VH' tar (t), the annual VH time series curve VH'(t) is segmented according to the maximum value time of the optimal prediction curve. tar (t). For the front section or the rear section of the optimal prediction curve, the stretched section curve is matched with the annual VH time series curve VH' tar (t) Calculate the MSD value for the corresponding section, and take the stretching result with the minimum MSD value as the section stretching matching result.

Images