Sentinel-1 VH polarization signal time sequence processing method based on years of spatio-temporal information

Abstract

The invention belongs to the technical field of remote sensing monitoring, and particularly provides a Sentinel-1VH polarized signal time sequence processing method based on years of spatio-temporal information, which is used for solving the problem that time sequence analysis is interfered by noise in a phenological detection process in the existing radar time sequence processing method. And the defect that a traditional time domain filtering method is directly used on a radar time sequence curve is overcome. According to the method, autocorrelation in a VH time sequence and regularity of agricultural activities are fully considered, multiple sampling points in a space window are utilized to obtain an annual time sequence curve averaged for many years, a regression relation is established between the annual time sequence curve and a real observed current-year curve after macroscopic stretching translation, and then local stretching of the curve is carried out to obtain a VH time sequence curve. The current year-round VH polarization signal time sequence observation curve of the target pixel of the detected crop is obtained and is used for subsequent phenology detection, and the accuracy of phenology detection is effectively improved.

Description

Technical Field

[0001] The invention belongs to the technical field of remote sensing monitoring, and specifically provides a Sentinel-1VH polarization signal time series processing method based on multi-year spatiotemporal information. Background Art

[0002] Sentinel-1 synthetic aperture radar data is a series of earth observation satellites launched by the European Space Agency (ESA). It uses the C-band (wavelength of about 5.6 cm) for imaging and can obtain high-resolution radar image data of the earth's surface in all weather and all day conditions. At present, Sentinel-1 synthetic aperture radar data has been widely used in various surface monitoring, including: research on land cover types, surface sedimentation, and vegetation growth phenology; 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 phenological information; although Sentinel-1 radar data can be unaffected by conditions such as cloud cover, its scattering intensity data is more susceptible to environmental noise and crop growth status when conducting time series analysis research, 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 noise as much as possible, researchers have proposed a variety of methods. The early methods were based on function fitting. This method used different mathematical functions to simulate the annual trajectory of radar time series, smoothed the original time series, and defined specific thresholds to mark the key time nodes of crop growth. In recent years, in order to make full use of synthetic aperture radar time series data for phenological detection, geometric features based on various types of polarization data have been proposed to find the relationship between the actual phenological stage of crops (such as sowing, germination, jointing, etc.) and the geometric properties of the time series (such as the mutation point or local extreme value of the slope). After smoothing, the geometric feature points can be directly extracted and correspond to the known phenological stages. This method can directly capture the phenological turning points of crops and reduce the errors caused by improper function fitting.

[0004] It can be seen that most existing methods focus on directly obtaining phenological periods through the geometric characteristics of the curve itself, and only perform simple time domain filtering or function fitting smoothing for related noise; however, these methods can only handle some high-frequency noise, and in time series data analysis, there is a lack of good solutions for low-frequency noise caused by environmental influences; such as soil characteristics, extreme environmental influences, weed growth, etc. These noises will have a significant impact on the geometric characteristics of radar time series data, thereby interfering with phenological detection results, especially when conducting phenological detection in a large area, due to different crop planting cycles and growth environments, such effects are further prominent. Therefore, it is urgent to further improve the Sentinel-1 VH polarization signal time series processing method. Summary of the invention

[0005] The purpose of the present invention is to propose a Sentinel-1 VH polarization signal time series processing method based on multi-year spatiotemporal information, in order to solve the problem that the time series analysis of the existing radar time series processing method is interfered by noise during phenological detection, and the deficiency of the traditional time domain filtering method directly used on the radar time series curve. The present invention fully considers the autocorrelation in the VH time series and the regularity of agricultural activities, and uses multiple sampling points in the spatial window to obtain the annual time series curve after multi-year averaging. After macroscopic stretching and translation, a regression relationship is established with the actual observed curve of the year, and then the optimal curve is obtained after local stretching of the curve.

[0006] To achieve the above object, the technical solution adopted by the present invention is:

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

[0008] S1. Obtain the annual VH time series curve of the target pixel of the detected crop in the current year, set a spatial window with the target pixel as the center, sample other pixels of the detected crop within the spatial window, and obtain the average VH time series curve of each sampled 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, the average VH time series curve of the filtered sampling pixel is translated and stretched to match, and the regression relationship between the annual VH time series curve of the target pixel and each average VH time series curve after translation and stretching is established one by one through the GMFR algorithm. The predicted VH time series curve of the target pixel in the current year is generated according to the regression relationship, and the MSD value is calculated with the annual VH time series curve of the filtered target pixel, and the predicted VH time series curve with the smallest MSD value is selected as the preferred predicted time series curve; among all sampling pixels, the preferred predicted time series curve with the smallest MSD value is selected as the optimal predicted time series curve;

[0011] S4. Divide the optimal prediction curve into segments with the maximum value as the center, and stretch and match the front segment and the back segment respectively; combine the stretch matching results of the front segment and the back segment to form the final prediction time series curve of the target pixel for detecting crops in the current year.

[0012] Furthermore, 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 detection crop, calculate the average value, and obtain the average VH time series curve:

[0013]

[0014] Among them, t represents the time variable, N i Indicates the sampling years of the ith sampling pixel, VH j (t) represents the VH time series curve of the i-th sampling pixel in the j-th sampling year, 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] Furthermore, in step S3, the calculation process of translation and stretching is:

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

[0017] Among them, 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].

[0018] Furthermore, in step S3, the generation process of the predicted VH timing curve is:

[0019] The average VH time series curve of the i-th sampling pixel after filtering is expressed as VH i ″(t), according to the time scale translation parameter and time scale stretch 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 translation and stretching processing;

[0020] The annual VH time series curve VH′ of the target pixel is established by GMFR algorithm tar (t) and the average VH time series curve of the sampled pixels VH i ″(T), the regression coefficients α and β are obtained, and the predicted VH time series curve VH of the target pixel of the detected crop in the current year is generated pre (t): VH pre (t) = α × VH″(T) + β.

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

[0022] The annual VH time series curve of the filtered target pixel is expressed as VH′ tar (t), the predicted VH timing curve VH pre (t) and the annual VH time series curve VH′ tar (t) Compare and calculate the MSD value:

[0023]

[0024] Wherein, Year represents the total duration of the VH timing curve.

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

[0026] T′=xscale×|tt max |

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

[0028] Furthermore, in step S4, the specific process of stretching matching is as follows:

[0029] The annual VH time series curve of the filtered target pixel is expressed as VH′ tar (t), according to the maximum value of the optimal prediction curve, the annual VH time series curve VH′ tar (t) perform the same segmentation;

[0030] For the front or back section of the optimal prediction curve, the stretched segmented curve is compared with the annual VH time series curve VH′ tar (t) Calculate the MSD value of the corresponding segment and take the stretching result with the smallest MSD value as the segmented stretching matching result.

[0031] Based on the above technical solution, the beneficial effects of the present invention are:

[0032] The present invention provides a Sentinel-1 VH polarization signal time series processing method based on multi-year spatiotemporal information. The method comprises the following steps: setting a spatial window around a target pixel and sampling other pixels of a detected crop within the spatial window to obtain a multi-year VH polarization signal time series of the sampled pixel, and fully integrating the time and space information within the spatial window around the target pixel; adopting an average calculation method for multi-year VH polarization signal time series of each sampled pixel, generating a multi-year average VH polarization signal time series for each sampled pixel, so that the complex environmental noise contained in the VH polarization signal time series is smoothed; adopting an overall translation and stretching calculation of the multi-year average VH polarization signal time series to avoid the macroscopic difference between the multi-year average VH polarization signal time series and the VH polarization signal time series actually observed in the current year; adopting a GMFR algorithm to establish a regression relationship between the multi-year average VH polarization signal time series after macroscopic stretching and translation and the VH polarization signal time series actually observed in the current year , based on the regression relationship, the predicted VH polarization signal time series actually observed in the current year is obtained, and the influence of some abnormal outlier VH values ​​is excluded; the regression prediction of the VH polarization signal time series actually observed in the current year and the MSD value calculation method between the VH polarization signal time series actually observed in the current year are adopted to obtain the optimal regression prediction of the VH polarization signal time series actually observed in the current year with the minimum MSD value; the method of stretching on both sides with the maximum value of VH as the center point is adopted to further make the optimal regression prediction of the VH polarization signal time series actually observed in the current year more able to reflect the regularity between crop growth and VH polarization signal changes locally; finally, the present invention has the ability to accurately restore the change law of the VH polarization signal time series based on many years of spatiotemporal information, exclude the noise influence of complex environmental factors on the VH polarization signal, and effectively identify the VH polarization signal pattern related to the crop phenological cycle. The VH polarization signal time series processed based on the present invention provides a better data quality foundation for the monitoring and prediction of crop phenology, and can further improve the phenological detection accuracy based on the VH polarization signal time series. BRIEF DESCRIPTION OF THE DRAWINGS

[0033] Figure 1 It is an algorithm flow chart of the Sentinel-1 VH polarization signal time series processing method based on multi-year spatiotemporal information;

[0034] Figure 2 It is a box plot of all samples between Phenocam ground phenology observation values ​​and corn phenology estimation values ​​in the embodiments of the present invention and the comparative examples.

[0035] Figure 3It is a box plot of all samples between the Phenocam ground phenology observation values ​​and the soybean phenology estimation values ​​in the embodiments of the present invention and the comparative examples.

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

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

[0038] Figure 6 The graph is a curve showing the change of the ratio of phenological phases of corn in Nebraska over time in the embodiments of the present invention and the comparative example.

[0039] Figure 7 The graph is a curve showing the change of the ratio of soybean phenological phases in Nebraska over time in the embodiments of the present invention and the comparative example.

[0040] Figure 8 It is a scatter plot comparing the rapeseed phenology detection results and phenology labels in Sichuan Province in the embodiments of the present invention and the comparative examples. DETAILED DESCRIPTION

[0041] In order to make the purpose, technical solution and beneficial effects of the present invention more clearly understood, the present invention is further described in detail below with reference to the accompanying drawings and embodiments.

[0042] The present invention provides a Sentinel-1 VH polarization signal time series processing method based on multi-year spatiotemporal information. The method uses multi-year averaged VH polarization time series data to generate a processed time series through macro and local stretching and translation calculations, which is used to reduce the noise generated by environmental factors and restore the corresponding law between vegetation growth phenology and time series, and can be used for phenology detection. The specific process is as follows: Figure 1 As shown, the following steps are included:

[0043] S1. Obtain the annual VH time series curve VH of the target pixel for detecting crops in the current year tar (t), set a spatial window with the target pixel as the center, and sample other pixels of the detection crop within the spatial window; for each sampled pixel: obtain the annual VH time series curve of the sampled pixel in each year of planting the detection crop, calculate the average value, and generate an average VH time series curve based on multi-year VH data;

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

[0045]

[0046] Among them, t represents the time variable, N i Indicates the sampling years of the ith sampling pixel, VH j (t) represents the VH time series curve of the i-th sampling pixel in the j-th sampling year, 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. Perform the same mean filtering on the target pixel's annual VH time series curve and the average VH time series curve of all sampled pixels. The target pixel's annual VH time series curve after filtering is expressed as VH′ tar (t), the average VH time series curve of the sampled pixels after filtering is expressed as VH i ″(t);

[0048] S3, the average VH timing curve of the sampled pixels after filtering is translated and stretched to match, and the average VH timing curve of the sampled pixels after translation and stretching is expressed as VH i ″(T); Establish the annual VH time series curve VH′ of the target pixel one by one through the GMFR algorithm tar (t) and each average VH time series curve VH i ″(T), excluding the influence of extreme outliers, and generating the predicted VH time series curve VH of the target pixel of the detected crop in the current year based on the regression relationship pre (t), and the VH time series curve VH′ of the target pixel after filtering throughout the year tar (t) Calculate the MSD value; for each sampling pixel, select the predicted VH timing curve with the smallest MSD value as the preferred predicted timing curve, that is, complete the translation and stretching matching of the current sampling pixel; further, among all sampling pixels, select the preferred predicted timing curve with the smallest MSD value as the optimal predicted timing curve;

[0049] Specifically, the translation and stretching process is expressed as:

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

[0051] Among them, 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 size is 6 days; xscale represents the time scale stretching parameter, the value range of xscale is [0.9, 1.1], and the step size is 0.05;

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

[0053] For each sampled pixel: establish the annual VH time series curve VH′ of the target pixel through the GMFR algorithm tar (t) and the average VH time series curve of the sampled pixels VH i The regression relationship of ″(T) is used to obtain the regression coefficients, including the slope parameter α and the intercept parameter β. According to α and β, the predicted VH time series curve VH of the target pixel of the detected crop in the current year is generated. pre (t):

[0054] VH pre (t) = α × VH″(T) + β

[0055] After the regression relationship is established, the predicted VH time series data of the target pixel for detecting crops in the current year are obtained by predicting the average time series data;

[0056] The predicted VH timing curve VH pre (t) and the annual VH time series curve VH′ of the filtered target pixel tar (t) Compare and calculate the MSD value:

[0057]

[0058] Among them, Year represents the total duration of the VH timing curve;

[0059] S4, segmenting the optimal prediction curve with the maximum value as the center, and stretching and matching the front and rear segments respectively; combining the stretching and matching results of the front and rear segments to form the final prediction time series curve of the target pixel of the detected crop in the current year, which is used as the data basis for subsequent phenological detection;

[0060] Specifically, for phenological detection, the optimal prediction curve is segmented with the maximum value as the center, into the front segment and the back segment. The front segment can be regarded as the time series before the crop growing season, and the back segment can be regarded as the time series after the crop growing season. The front segment and the back segment are stretched and matched respectively, which is specifically expressed as:

[0061] T′=xscale×|tt max |

[0062] Where T′ represents the time variable after stretching treatment, t max Indicates the maximum moment of the optimal prediction curve; the value range of the time scale stretching parameter xscale is [0.9, 1.1], and the step size is 0.05;

[0063] Similar to the formula for translation and stretching matching in S2, in S2, stretching is performed with the starting time as the starting point, and in S4, the optimal prediction curve is stretched with the maximum value as the center;

[0064] According to the maximum time t of the optimal prediction curve max The annual VH time series curve VH′ of the filtered target pixel tar (t) The same segmentation is performed. On this basis, for the front or back segment of the optimal prediction curve, the stretched segmented curve is compared with the annual VH time series curve VH′ tar (t) Calculate the MSD value of the corresponding segment, and take the stretching result with the smallest MSD value as the segmented stretching matching result; Combine the stretching matching results of the front segment and the back segment to obtain the annual VH polarization signal time series observation curve of the target pixel for detecting crops in the current year after the Sentinel-1 VH polarization signal time series processing based on multi-year spatiotemporal information in the present invention.

[0065] Taking American corn, soybean and Chinese rapeseed as examples, the beneficial effects of the Sentinel-1 VH polarization signal time series processing method based on multi-year spatiotemporal information in the present invention are described in detail below.

[0066] The Sentinel-1 VH polarization signal time series processing method based on multi-year spatiotemporal information was used to conduct phenological detection experiments in corn and soybean planting areas in the United States and rapeseed planting areas in Sichuan Province, China. The LOESS filter group was used as a comparative example to compare and analyze the time series data processed by the LOESS filter. The three test scenarios have different characteristics:

[0067] The first test scenario used PhenoCam phenological data of corn and soybean in the United States. Its characteristics were to obtain radar time series data based on the positioning of a single phenological site and then conduct phenological detection. The digital photos provided by the Phenocam website were used to interpret the phenological period as label data.

[0068] The second test scenario is corn and soybeans in 12 states in the U.S. Corn Belt. It is characterized by regional phenology detection over a large spatial range. The crop area that reaches a certain phenological period in the survey provided by the U.S. Department of Agriculture in the Crop Progress Report of each state is used as label data. Corn and soybeans in the United States are generally planted around May.

[0069] The third test scenario was Sichuan Province, China. As a cloudy and foggy area, it is very difficult to obtain optical data. This can verify the availability of Sentinel-1 radar data for phenological detection in this area and the effectiveness of the processing algorithm. This test scenario uses farmers to take pictures on the spot and agricultural experts to interpret the photos to obtain phenological period labels.

[0070] In order to conduct quantitative evaluation, two methods, SMF-S algorithm and DTW algorithm, were used for phenological detection. The phenological detection results were compared with the label results, and the MAE value was calculated.

[0071] A quantitative evaluation was conducted on the corn and soybean Phenocam sites in the United States. The test method was based on the time series data processed by the method of the present invention and the time series data filtered by the LOESS algorithm. The SMF-S algorithm and the DTW two phenological detection algorithms were used to compare the phenological detection accuracy. As shown in Tables 1 and 2, the phenological detection accuracy of the three phenological detection algorithms on corn and soybean crops based on the time series data processed by the method of the present invention and the time series data filtered by the LOESS algorithm are shown respectively.

[0072] Table 1

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

[0074] Table 2

[0075] Phenological period VC V2 R7 SMF-S (LOESS filter) 16.13 14.29 16.97 SMF-S (the present invention) 11.03 8.63 15.59 DTW (LOESS filter) 19.50 19.37 17.31 DTW (present invention) 16.69 11.24 14.62

[0076] like Figure 2 , Figure 3As shown, box plots of all samples between Phenocam ground phenological observations and phenological estimates detected by the two algorithms based on two types of data are shown; by averaging the MAEs of all phenological periods, for corn, the MAEs of phenological detection by the two algorithms after processing by the method of the present invention are 12.23 and 21.77, respectively, while the corresponding MAEs for soybean are 11.75 and 14.19; after phenological detection of time series data simply after LOESS filtering, for corn, the MAEs of phenological detection by the two algorithms are 18.40 and 24.92, respectively, while the MAEs for phenological detection by soybean are 15.80 and 18.73, respectively; for the two algorithms, the data processed by the method of the present invention obviously have more advantages in the field of phenological detection, especially when the SMF-S algorithm performs phenological detection on the data processed by the method of the present invention, the final accuracy is improved by 4 to 6 days compared with the original data (corn: 18.40 days vs. 12.23 days; soybean: 15.80 days vs. 11.75 days). Judging from the results, in the phenological detection of the two algorithms, the data processed by the method of the present invention has stronger stability. Compared with the original data, the algorithm reduces the impact of low-frequency environmental noise and restores the regularity of the VH time series and the crop phenological growth cycle. At the same time, with Phenocam stations spread across the United States, the processing of all sample time series is relatively stable within a large and complex spatial range, suppressing the obstacles to VH time series analysis caused by multiple environmental factors.

[0077] The phenology of corn and soybean in the 12 states of the Corn Belt of the United States was detected in a large spatial range. Based on the Cropland Data Layer data released by the United States Department of Agriculture, 200 crop pixels were selected in each state and used as sample data. All sample data were processed by LOESS filtering and the method of the present invention to perform phenology detection on the data. Based on the Crop Progress Report (CPR) provided by the United States Department of Agriculture, the 50% area ratio was used as the threshold, and the date when the phenological area in the state's agricultural survey reached 50% was obtained as the label of the state's crop phenological period. The median of the final result of the phenological detection in each state was compared with the phenological period result. Figure 4 , Figure 5As shown, a scatter plot is shown comparing the phenological detection results of corn and soybean in each state with the phenological labels provided by CPR. It can be seen from the figure that the performance of the data processed by the method of the present invention is better than that of the data processed only by the LOESS filtering algorithm in both algorithms; in addition, the estimation of all samples by the method of the present invention is relatively stable, without outliers, that is, without significantly large errors, indicating that the present invention is robust; and, for some points that are obviously outliers in the original data, after being processed by the method of the present invention, the phenological detection results return to the normal range; the phenological detection MAE is greatly improved. As far as the SMF-S algorithm is concerned, the average MAE of phenological detection of soybeans and corn decreases by about 3 to 4 days before and after being processed by this method (soybeans: 11.32 days vs. 8.56 days; corn: 11.10 days vs. 7.61 days). Figure 6 , Figure 7 As shown, the curves of the change over time of the proportion of crop area reaching the phenological period recorded by the CPR data of corn and soybean crops in Nebraska, the curve of the change over time of the proportion of sample data reaching the phenological period after being processed by the method of the present invention detected by the SMF-S algorithm, and the curve of the change over time of the proportion of sample data reaching the phenological period after being processed by the LOESS filter are compared. By comparing the changing trends of the three curves, it is found that the proportion of sample data reaching the phenological period after being processed by the method of the present invention detected by the SMF-S algorithm is more consistent with the changing trend of the corresponding phenological period area proportion investigated by the CPR data, which can reflect the growth phenological cycle of crops in Nebraska.

[0078] The phenological detection of some rapeseed in Sichuan Province, China, was carried out based on the interpretation of multiple field photos taken by farmers to obtain the rapeseed growth phenological labels of the located fields; the phenological detection MAE of the sample data processed by the LOESS filtering algorithm and the method of the present invention under the two phenological detection algorithms SMF-S and DTW was compared. The LOESS algorithm filtered out some low-frequency noise and effectively smoothed the time series data, but did not effectively remove the impact of environmental factors on the crop phenological growth signal. On the basis of effectively removing noise, the method of the present invention retained the regularity between the crop phenological growth law and the SAR scattering intensity signal, as shown in Table 3, which shows the accuracy of the two algorithms in phenological detection under two different processing methods.

[0079] Table 3

[0080] Phenological period Bud stage Early flowering period Blooming period Pod-setting period Maturity Harvest period SMF-S (LOESS filter) 19.50 24.33 19.38 18.63 26.25 32.11 SMF-S (the present invention) 11.00 11.78 9.44 8.56 8.44 14.00 DTW (LOESS filter) 34.39 25.39 17.28 14.83 17.11 7.78 DTW (present 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 impact on the accuracy of phenological detection, and the DTW overall matching algorithm has a large deviation, but the method of the present invention still improves the final accuracy of phenological detection to a certain extent; the local matching method of SMF-S makes it less affected by environmental noise, but the MAE of phenological detection of data after LOESS filtering is still generally greater than 18 days, and the MAE of time series data processed by the method of the present invention is generally around 10 days, and the accuracy of pod-setting and maturity stages is the best (pod-setting stage: 18.63 days vs. 8.56 days; maturity stage: 26.25 days vs. 8.44 days). Figure 8 As shown, a scatter plot is shown which compares the rapeseed phenology detection results in Sichuan Province with the field monitoring phenology labels. It can be seen that the data processed by the method of the present invention can significantly improve the accuracy in the field of phenology detection, and there is still an obvious denoising effect even in mountainous and complex areas such as Sichuan Province where the small farmer economy is dominant.

[0082] The above description is only a specific implementation mode of the present invention. Any feature disclosed in this specification, unless otherwise stated, can be replaced by other alternative features that are equivalent or have similar purposes; all the disclosed features, or all the steps in the 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 polarization signals based on multi-year spatiotemporal information, characterized in that: The following steps are involved: S1. Obtain the annual VH time series curve of the target pixel of the detected crop in the current year, set a spatial window with the target pixel as the center, sample other pixels of the detected crop within the spatial window, and obtain the average VH time series curve of each sampled pixel; 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; S3. For each sampling pixel, the average VH time series curve of the filtered sampling pixel is translated and stretched to match, and the regression relationship between the annual VH time series curve of the target pixel and each average VH time series curve after translation and stretching is established one by one through the GMFR algorithm. The predicted VH time series curve of the target pixel in the current year is generated according to the regression relationship, and the MSD value is calculated with the annual VH time series curve of the filtered target pixel, and the predicted VH time series curve with the smallest MSD value is selected as the preferred predicted time series curve; among all sampling pixels, the preferred predicted time series curve with the smallest MSD value is selected as the optimal predicted time series curve; S4. Divide the optimal prediction curve into segments with the maximum value as the center, and stretch and match the front segment and the back segment respectively; combine the stretch matching results of the front segment and the back segment to form the final prediction time series curve of the target pixel for detecting crops in the current year.

2. According to claim 1, the method for processing the time series of Sentinel-1 VH polarization signals based on multi-year spatiotemporal information is characterized in that: 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 detection crop, calculate the average value, and obtain the average VH time series curve: Among them, t represents the time variable, N i Indicates the sampling years of the ith sampling pixel, VH j (t) represents the VH time series curve of the i-th sampling pixel in the j-th sampling year, VH′ i (t) represents the average VH time series curve of the i-th sampling pixel, and I represents the number of sampling pixels.

3. According to claim 1, the method for processing the time series of Sentinel-1 VH polarization signals based on multi-year spatiotemporal information is characterized in that: In step S3, the calculation process of translation and stretching is: T=xscale×(t+tshift) Among them, 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. According to claim 1, the method for processing the time series of Sentinel-1 VH polarization signals based on multi-year spatiotemporal information is characterized in that: In step S3, the generation process of the predicted VH timing curve is: The average VH time series curve of the i-th sampling pixel after filtering is expressed as VH″ i (t), according to the time scale translation parameter and time scale stretch parameter, the average VH time series curve VH″ i (t) Perform translation and stretching calculations to obtain the corresponding average VH timing curve VH″ i (T), T represents the time variable after translation and stretching processing; The annual VH time series curve VH′ of the target pixel is established by GMFR algorithm tar (t) and the average VH time series curve VH″ of the sampling pixels i (T), the regression coefficients α and β are obtained, and the predicted VH time series curve VH of the target pixel of the detected crop in the current year is generated pre (t): VH pre (t) = α × VH″(T) + β.

5. According to claim 1, the method for processing the time series of Sentinel-1 VH polarization signals based on multi-year spatiotemporal information is characterized in that: In step S3, the calculation process of the MSD value is: The annual VH time series curve of the filtered target pixel is expressed as VH′ tar (t), the predicted VH timing curve VH pre (t) and the annual VH time series curve VH′ tar (t) Compare and calculate the MSD value: Wherein, Year represents the total duration of the VH timing curve.

6. The method for processing the time series of Sentinel-1 VH polarization signals based on multi-year spatiotemporal information according to claim 1, characterized in that: In step S4, the stretching calculation process is: T′=xscale×|t-t max | Where T′ represents the time variable after stretching treatment, t max Indicates the maximum moment of the optimal prediction curve; xscale represents the time scale stretching parameter, and its value range is [0.9,1.1].

7. The method for processing the time series of Sentinel-1 VH polarization signals based on multi-year spatiotemporal information according to claim 1, characterized in that: In step S4, the specific process of stretch matching is as follows: The annual VH time series curve of the filtered target pixel is expressed as VH′ tar (t), according to the maximum value of the optimal prediction curve, the annual VH time series curve VH′ tar (t) perform the same segmentation; For the front or back section of the optimal prediction curve, the stretched segmented curve is compared with the annual VH time series curve VH′ tar (t) Calculate the MSD value of the corresponding segment and take the stretching result with the smallest MSD value as the segmented stretching matching result.

Images