Denoising reconstruction method and device for irregular multi-year time sequence NDVI data

Abstract

The invention discloses a denoising reconstruction method and device for irregular multi-year time sequence NDVI (normalized difference vegetation index) data, and the method comprises the steps: reading NDVI time sequence data of each pixel of a multi-year fusion image, sorting according to the number of years and days, and carrying out the first denoising reconstruction of the NDVI time sequence data of each year; and constructing a multi-year synchronous set of each pixel according to a first denoising reconstruction result, performing secondary denoising reconstruction on the multi-year synchronous set of each pixel, storing the denoising reconstruction result into a set LMN, writing data in the set LMN into the image, and obtaining a denoised NDVI time sequence image. According to the method, the dependence of a traditional algorithm on continuous time series data is improved, and the availability and continuity of vegetation remote sensing monitoring data in cloudy and rainy areas are improved.

Description

Technical Field

[0001] The present invention relates to a method for reconstructing NDVI data, and in particular, to a denoising and reconstruction method and device for irregular multi-year time-series NDVI data. Background Art

[0002] Remote sensing is the most effective means for large-scale monitoring of the dynamic changes of ground objects. In recent years, with the rapid development of remote sensing science and technology, earth observation satellites with different bands and different spatio-temporal spectral resolutions have been launched one after another, and earth observation data resources with better time continuity, higher spatial resolution, and more information can be obtained. The time-series data of vegetation parameters obtained based on remote sensing are widely used in the monitoring of surface vegetation dynamics. The Normalized Difference Vegetation Index (NDVI) is one of the most commonly used remote sensing vegetation parameters, which reflects the unique spectral reflection characteristics of surface vegetation in the visible and near-infrared bands and is the best indicator factor for vegetation growth status and vegetation coverage. Since the satellite remote sensing NDVI time-series data show cycles and changes related to the biological characteristics of vegetation, they are widely used in research fields such as regional and global-scale ecological environment monitoring, climate change analysis, and agricultural production assessment.

[0003] China's tropical and subtropical regions are affected by the monsoon climate, with frequent cloudy and rainy weather, resulting in a reduction in the number of effective optical images that can be used. In addition, the region is mountainous and hilly, with a complex surface environment and a high degree of landscape fragmentation. Therefore, it is difficult to achieve high-precision monitoring of regional vegetation dynamics using a single satellite remote sensing data source, and it is necessary to coordinate the comprehensive advantages of multiple satellite sensors for earth observation. Among the current mainstream optical remote sensing data sources, the Landsat series of satellites can finely depict surface details with a spatial resolution of 30 meters, but its 16-day revisit cycle is difficult to meet the needs of high-frequency monitoring of vegetation growth; although MODIS satellite data has the advantage of high temporal resolution with daily coverage, its maximum spatial resolution is 250 meters, and the low spatial resolution limits the ability to resolve local vegetation characteristics. To this end, many scholars have conducted vegetation change monitoring based on the NDVI time series dataset constructed by the spatiotemporal fusion of Landsat and MODIS, aiming to build a high spatiotemporal resolution NDVI continuous observation dataset by coordinating the high spatial resolution of Landsat and the high temporal resolution of MODIS, thereby improving the efficiency of vegetation monitoring. However, such fused images still face some technical bottlenecks in practical applications: first, natural factors such as cloud cover and atmospheric disturbances lead to spatiotemporal discontinuity in image sequences, affecting data integrity; second, complex noise sources such as satellite sensor errors, terrain irradiance differences, surface cover heterogeneity, and human activity interference make the fused NDVI time series data mixed with high-frequency random noise and systematic deviations. Such noise will mask the true signal of vegetation characteristics, resulting in a significant decrease in the accuracy of ecosystem parameter inversion.

[0004] Current research on noise removal and reconstruction of time series NDVI data focuses on noise suppression of NDVI time series with regular step lengths (such as 8-day / 16-day intervals). Stable denoising can be achieved by relying on classical filtering or interpolation algorithms. However, in multi-source data fusion scenarios, due to factors such as asynchronous acquisition time of the original images and cloud pollution removal, NDVI series often present irregular step length characteristics, which are manifested as irregular time intervals or partial missing of data points. Such irregular time series will destroy the stationarity assumption of traditional algorithms, resulting in problems such as time series distortion or over-smoothing in the denoising results. Summary of the invention

[0005] In view of the problems existing in the prior art, the purpose of the present invention is to provide a denoising and reconstruction method and device for irregular NDVI time series. The method improves the dependence of traditional algorithms on continuous time series data, and the reconstruction results are more accurate.

[0006] In order to achieve the above-mentioned object of the invention, the present invention provides the following technical solutions:

[0007] A denoising and reconstruction method for irregular multi-year time series NDVI data includes the following steps:

[0008] (1) Read the NDVI time series data of each pixel in the multi-year fused image and store it in the NDVI data set LM;

[0009] (2) Sort the NDVI data in the NDVI data set LM by the number of days in a year to obtain the NDVI time series set LMY;

[0010] (3) Read the NDVI data of any year from the NDVI time series set LMY, store it in the set LMY1, and perform the first denoising and reconstruction on the set LMY1;

[0011] (4) Loop through step (3) to traverse the NDVI data of each year in the set LMY, complete the first denoising and reconstruction of the NDVI data of all years, and store the denoising and reconstruction results in the set LM';

[0012] (5) Based on the set LM', construct the multi-year synchronous set of each pixel and store it in the NDVI time series set LM";

[0013] (6) Read the NDVI data of any pixel from the NDVI time series set LM" and perform secondary denoising and reconstruction;

[0014] (7) Loop through step (6) to traverse the NDVI data of all pixels in the set LM", complete the secondary denoising and reconstruction of the NDVI data of all pixels, store the denoising and reconstruction results in the set LMN, and write the data in the set LMN into the image to obtain the denoised NDVI time series image.

[0015] Further, step (3) includes:

[0016] (3-1) Read the NDVI data of any year from the NDVI time series set LMY and store it in the set where, represents the NDVI data of the pixel in the y0th year, the dth day, the ith row, and the jth column in the set LMY1, and k, g, f represent the number of days, the number of rows, and the number of columns respectively;

[0017] (3-2) Arbitrarily select an element in the set LMY1 Record the indexes of all missing value positions in, and check the lengths of all consecutive missing value segments in; if there is a consecutive missing value segment with a length greater than or equal to 10, then use the head of the consecutive missing value segment as the segmentation point to segment it from the front valid data segment, and use the tail as the segmentation point to segment it from the rear valid data segment, so as to segment the consecutive missing value segment into an independent missing value data segment; if the length is less than 10, no processing is performed, and finally is divided into several valid data segments and missing value data segments and stored in the set LMY2;

[0018] (3 - 3) Arbitrarily select an element lmy2 from the set LMY2 u . If all the data in the element lmy2 u are missing values, skip the processing; otherwise, count the number n of NDVI data of the element lmy2 u . If n ≤ 3, skip the processing; if 3 < n ≤ 30, execute the first outlier test method, if n > 30, execute the second outlier test method, and correct the detected outliers to the mean of the two adjacent NDVI data; if the outlier is at the head or tail position, correct the outlier to the NDVI data at the adjacent position, thus completing the denoising and reconstruction of lmy2 u .

[0019] (3 - 4) Loop and execute step (3 - 3) until all elements of the set LMY2 are traversed, completing the denoising and reconstruction of .

[0020] (3 - 5) Loop and execute steps (3 - 2) to (3 - 4) until all elements in the set LMY1 are traversed, completing the first - round denoising and reconstruction of all NDVI data in the set LMY1, and storing the processed data segments in the set LM'.

[0021] Furthermore, the specific steps of the first outlier test method described in step (3.3) include:

[0022] (3 - 3 - 1 - 1) Arrange the NDVI data in lmy2 u in ascending order to obtain the NDVI data sequence {ndvi1, …, ndvi v , …, ndvi n}, where ndvi1, ndvi v , ndvi n represent the NDVI data on the t1, t v , t n - th days respectively;

[0023] (3 - 3 - 1 - 2) Set the significance level α and the critical value Q;

[0024] (3 - 3 - 1 - 3) Calculate the difference d1 between ndvi n and ndvi1: d1=(ndvi n - ndvi1);

[0025] (3 - 3 - 1 - 4) If there is high - end noise in lmy2 u , calculate the difference d2 between ndvi n and ndvi n-1 : d2=(ndvi n - ndvin-1 ) and calculate the ratio Q of d2 to d1 c ;

[0026] (3 - 3 - 1 - 5) If lmy2 u has low - end noise, calculate the difference d3 = (ndvi2 - ndvi1) between ndvi2 and ndvi1, and calculate the ratio Q of d3 to d1 d ;

[0027] (3 - 3 - 1 - 6) If Q c > Q, determine that the NDVI data ndvi u in lmy2 n is an outlier; if Q d > Q, determine that the NDVI data ndvi1 in lmy2 u is an outlier.

[0028] Furthermore, the specific steps of the second outlier test method in step (3.3) include:

[0029] (3 - 3 - 2 - 1) For lmy2 u , use the Savitzky - Golay polynomial function method for fitting to obtain the NDVI fitting data;

[0030] (3 - 3 - 2 - 2) Based on the day - of - year data T corresponding to the NDVI data in lmy2 u , construct the hat matrix H according to the following formula:

[0031]

[0032] H = X(X ′ X) -1 X ′

[0033] where ndvi1, ndvi v , ndvi n represent the NDVI data on the t1, t v , t n days respectively, X is the design matrix, the first column of X is a vector of all 1s, and the second column is the day - of - year data corresponding to the NDVI data in lmy2 u ; X ′ represents the transpose of X; the diagonal element h vv of H represents the weight of the ndvi v data during fitting;

[0034] (3 - 3 - 2 - 3) Based on the NDVI data and the NDVI fitting data in lmy2 u , calculate the studentized deleted residual according to the following formula;

[0035]

[0036] In the formula, t v represents the studentized deleted residual, p represents the order selected by the Savitzky-Golay polynomial function method, r v represents the residual, and ndvi v represents the v-th NDVI data in lmy2 u and represents the corresponding value in the NDVI fitting data corresponding to ndvi v , and SSE represents the sum of the squares of the residuals;

[0037] (3-3-2-4) Calculate the threshold BC according to the following formula:

[0038]

[0039] In the formula, b is the significance value, represents the critical value of the t-distribution with degrees of freedom n-p-1 and significance level ;

[0040] (3-3-2-5) Determine whether each NDVI data ndvi u in lmy2 v satisfies -e*BC < ndvi v < e*BC, where e is the correction factor. If so, determine that ndvi v is an outlier.

[0041] Furthermore, step (5) includes:

[0042] (5-1) Extract the NDVI data of each pixel every ten days over the years from the set LM′ represents the multi-year synchronous data set of the pixel in the r0-th row and l0-th column, represents the NDVI data of the pixel in the r0-th row and l0-th column on the m-th and m+9-th days in the y-th year, k represents the total number of days, and t represents the total number of years;

[0043] (5-2) Store the multi-year synchronous data sets of all pixels into the set where g and f represent the total number of rows and columns of the fused image, respectively.

[0044] Furthermore, step (6) includes:

[0045] (6-1) Read the multi-year synchronous data set of any pixel from LM″

[0046] (6-2) Elimination All missing values ​​in the data are stored in the collection

[0047] (6-3) Statistics The number of NDVI data in the table is N. If N≤3, the processing is skipped. <N≤30,则执行第一异常值检验方法,若N≥30,则执行第二异常值检验方法;采用除异常值外其它NDVI数据的平均值替换检验到的异常值,完成 Denoising and reconstruction.

[0048] A computer device comprises a memory, a processor and a computer program stored in the memory and executable on the processor, wherein the processor executes the computer program to implement the above method.

[0049] A computer-readable storage medium having a computer program / instruction stored thereon, wherein the computer program / instruction implements the above method when executed by a processor.

[0050] A computer program product comprises a computer program / instruction, wherein the computer program / instruction implements the above method when executed by a processor.

[0051] Compared with the prior art, the present invention has the following beneficial effects: the irregular multi-year time series NDVI data is subjected to two denoising processes. First, the local outliers are subjected to Dixon test and Studentized residual test within the year to correct and smooth the data; second, based on the periodicity of the NDVI values ​​over the same period of the year, the data quality is further optimized between years, the outliers in the data are accurately identified and removed, and more reliable NDVI time series data is restored, eliminating the influence of outliers and making the reconstruction results more accurate. BRIEF DESCRIPTION OF THE DRAWINGS

[0052] Figure 1 It is a flow chart of a denoising and reconstruction method for irregular multi-year time series NDVI data provided by the present invention;

[0053] Figure 2 It is a flow chart of the first denoising and reconstruction of the fused image provided by the present invention;

[0054] Figure 3 It is a flow chart of secondary denoising and reconstruction of fused images provided by the present invention;

[0055] Figure 4 is a representative fusion image data diagram used in an embodiment of the present invention, wherein the framed portion is the experimental area of ​​this embodiment;

[0056] Figure 5It is a display diagram of Dixon outlier detection for NDVI data of a certain pixel in a certain year and a certain time series (continuous length less than 30 days) in an embodiment of the present invention;

[0057] Figure 6 It is a display diagram of studentized residual outlier detection for NDVI data of a certain pixel in a certain year and a certain time series (continuous length greater than 30 days) in an embodiment of the present invention;

[0058] Figure 7 It is a display diagram of studentized residual outlier detection for NDVI data of a certain pixel in the same period of multiple years (day of the year 61 - 70) in an embodiment of the present invention;

[0059] Figure 8 It is a display diagram of consecutive time series images of the same year for the first noise detection in the experimental area of an embodiment of the present invention; where (a) is the fused image of the experimental area on the 271st day of 2013, (a′) is the position of noise points in the experimental area on the 271st day, (b) is the fused image of the experimental area on the 272nd day of 2013, (b′) is the position of noise points in the experimental area on the 272nd day, (c) is the fused image of the experimental area on the 273rd day of 2013, (c′) is the position of noise points in the experimental area on the 273rd day, (d) is the fused image of the experimental area on the 274th day of 2013, (d′) is the position of noise points in the experimental area on the 274th day;

[0060] Figure 9 It is a display diagram of multi - year synchronous images for the secondary noise detection in the experimental area of an embodiment of the present invention; where (a) is the fused images of the experimental area on the 271st, 272nd, 273rd, and 274th days of 2013 after the first denoising and reconstruction, (a′) is the position of noise points in the experimental area of 2013 after the first denoising and reconstruction, (b) is the fused images of the experimental area on the 271st, 272nd, 273rd, and 274th days of 2015 after the first denoising and reconstruction, (b′) is the position of noise points in the experimental area of 2015 after the first denoising and reconstruction, (c) is the fused images of the experimental area on the 271st, 272nd, 273rd, and 274th days of 2018 after the first denoising and reconstruction, (c′) is the position of noise points in the experimental area of 2018 after the first denoising and reconstruction, (d) is the fused images of the experimental area on the 271st, 272nd, 273rd, and 274th days of 2019 after the first denoising and reconstruction, (d′) is the position of noise points in the experimental area of 2019 after the first denoising and reconstruction;

[0061] Figure 10 It is a result display diagram of the experimental area after denoising and reconstruction in an embodiment of the present invention; where (a) is the reconstructed image of the experimental area on the 273rd day of 2013, (b) is the reconstructed image of the experimental area on the 273rd day of 2015, (c) is the reconstructed image of the experimental area on the 273rd day of 2018, (d) is the reconstructed image of the experimental area on the 273rd day of 2019. Detailed implementation mode

[0062] Next, the technical solutions in the embodiments of the present invention will be clearly and completely described in conjunction with the accompanying drawings in the embodiments of the present invention. The experimental data of this embodiment is based on the GEE platform and uses the ESTARFM (Enhanced Spatial and Temporal Adaptive Reflectance Fusion Model) algorithm to fuse Landsat 8 images and MODIS product data. The projection coordinate system used for this image is WGS_1984_UTM_Zone_50N. This will be further illustrated below in conjunction with the accompanying drawings and by describing a specific embodiment.

[0063] The embodiment of the present invention provides a denoising and reconstruction method for irregular multi-year time-series NDVI data, as Figure 1 shown, including the following steps:

[0064] (1) Read the NDVI time-series data of each pixel in the multi-year fused image and store it in the NDVI data set LM = {lm (y,d,r,l) |y = 1, 2,..., t; d = 1, 2,..., k; r = 1, 2,..., g; l = 1, 2,..., f}; where, lm (y,d,r,l) represents the NDVI data of the y-th year, d-th day, r-th row, and l-th column of the Landsat and MODIS fused image; t, k, g, and f respectively represent the total number of years, total number of days, total number of rows, and total number of columns. In this embodiment, t = 10, k = 1747, g = 982, f = 557, and the experimental area is as Figure 4 shown.

[0065] (2) Sort the NDVI data in the NDVI data set LM according to the year and day to obtain the NDVI time-series set LMY = {lmy (y,d,1,1) , lmy (y,d,1,3) , …, lmy (y,d,i,j) , …, lmy (y,d,g,f)}, lmy (y,d,i,j) represents the NDVI data of the pixel in the i-th row and j-th column on the d-th day in the y-th year in the set LMY1.

[0066] (3) Read the NDVI data of any year from the NDVI time-series set LMY, store it in the set LMY1, and perform the first denoising and reconstruction on the set LMY1.

[0067] As Figure 2 shown, this step includes:

[0068] (3-1) Read the NDVI data of any year from the NDVI time-series set LMY and store it in the set where, Denote the NDVI data of the pixel at the \(i\)-th row and \(j\)-th column on the \(d\)-th day in the \(y_0\)-th year in the set \(LMY1\), where \(k\), \(g\), and \(f\) represent the number of days, the number of rows, and the number of columns respectively; in this embodiment, \(y_0 = 2013\), \(d = 225\), \(i = 200\), \(j = 200\).

[0069] (3 - 2) Randomly select an element from the set \(LMY1\). Record the indices of the positions of all missing values (Nodata) in, and check the lengths of all consecutive missing value segments in; if there exists a consecutive missing value segment with a length greater than or equal to 10, then use the head of this consecutive missing value segment as the splitting point to split it from the valid data segment in front, and use the tail as the splitting point to split it from the valid data segment behind, so as to split the consecutive missing value segment into an independent missing value data segment; if the length is less than 10, then no processing is done, and finally is divided into several valid data segments and missing value data segments, and stored in the set \(LMY2\).

[0070] (3 - 3) Randomly select an element \(lmy2\) from the set \(LMY2\). u , if all the data in the element \(lmy2\) u are missing values, then skip the processing; otherwise, count the number \(n\) of NDVI data in the element \(lmy2\) u , if \(n\leq3\), then skip the processing; if \(3 < n\leq30\), then execute the first outlier test method, if \(n > 30\), then execute the second outlier test method, and correct the detected outliers to the mean of the two adjacent NDVI data; if the outlier is at the head or tail position, then this outlier is corrected to the NDVI data at the adjacent position, thus completing the denoising and reconstruction of \(lmy2\) u .

[0071] (3 - 4) Loop and execute step (3 - 3) until all elements in the set \(LMY2\) are traversed, and complete the denoising and reconstruction of .

[0072] (3 - 5) Loop and execute steps (3 - 2) to (3 - 4) until all elements in the set \(LMY1\) are traversed, complete the first denoising and reconstruction of all NDVI data in the set \(LMY1\), and store the processed data segments in the set \(LM'\).

[0073] Among them, the specific first outlier test method is the Dixon test, including:

[0074] (3 - 3 - 1 - 1) Arrange the NDVI data in \(lmy2\) u in ascending order to obtain the NDVI data sequence \(\{ndvi1,\ldots,ndvi v ,\ldots,ndvi n}, ndvi1, ndvi v , ndvi n represent the NDVI data on the 1st, v-th, and n-th days respectively;

[0075] (3 - 3 - 1 - 2) Set the significance level α and set the critical value Q;

[0076] (3 - 3 - 1 - 3) Calculate the difference d1 between ndvi n and ndvi1: d1 = (ndvi n - ndvi1);

[0077] (3 - 3 - 1 - 4) If there is high - end noise in lmy2 u , calculate the difference d2 between ndvi n and ndvi n-1 : d2 = (ndvi n - ndvi n-1 ), and calculate the ratio of d2 to d1

[0078] (3 - 3 - 1 - 5) If there is low - end noise in lmy2 u , calculate the difference d3 between ndvi2 and ndvi1: d3 = (ndvi2 - ndvi1), and calculate the ratio of d3 to d1

[0079] (3 - 3 - 1 - 6) If Q c > Q, then determine that the NDVI data ndvi u in lmy2 n is an outlier; if Q d > Q, then determine that the NDVI data ndvi1 in lmy2 u is an outlier.

[0080] In this embodiment, α is 0.05. Table 1 is as follows, Q c = 0.56, Q d = 0.002, Q = 0.44. The Dixon test for outlier processing effect is as Figure 5 shown;

[0081] Table 1 Critical value table of Dixon test method

[0082] n Critical value Q for α = 0.01 Critical value Q for α = 0.05 5 0.850 0.739 6 0.796 0.681 7 0.756 0.639 8 0.726 0.608 9 0.790 0.700 10 0.767 0.676 … … …

[0083] Among them, the second outlier test method is specifically the studentized residual detection method, including:

[0084] (3 - 3 - 2 - 1) For lmy2 u , use the Savitzky - Golay polynomial function method for fitting to obtain the NDVI fitting data;

[0085] (3 - 3 - 2 - 2) Based on lmy2 u For the annual day - of - year data T corresponding to the NDVI data, construct the hat matrix H according to the following formula: (indicating the influence degree of each data point on the fitting result):

[0086]

[0087] H = X(X′X) -1 In the formula of X′, ndvi1, ndvi v , ndvi n respectively represent the NDVI data on the t1, t v , t n days. X is the design matrix. The first column of X is a vector of all 1s, and the second column is the annual day - of - year data corresponding to the NDVI data in lmy2 u ; X′ represents the transpose of X; the diagonal element h of H vv represents the weight of the ndvi v data during fitting;

[0088] (3 - 3 - 2 - 3) Based on lmy2 u For the NDVI data and NDVI fitting data in lmy2, calculate the studentized deleted residual according to the following formula;

[0089]

[0090] In the formula, t v represents the studentized deleted residual, p represents the order selected by the Savitzky - Golay polynomial function method, p = 2, indicating the use of quadratic polynomial regression. In this embodiment, p = 2, r v represents the residual, ndvi v represents the v - th NDVI data in lmy2 u , represents the corresponding value in the NDVI fitting data to ndvi v ; SSE represents the sum of squared residuals, h vv is the v - th diagonal element of the hat matrix H;

[0091] (3 - 3 - 2 - 4) Calculate the threshold BC according to the following formula:

[0092]

[0093] In the formula, b is the significance value, represents the critical value of the t - distribution with degrees of freedom n - p - 1 and significance level ;

[0094] (3-3-2-5) Determine for lmy2 u in each NDVI data ndvi v , whether it satisfies -e*BC < ndvi v < e*BC, where e is the correction factor. If so, then determine that ndvi v is an outlier. In this embodiment, b = 0.05, e = 0.86, and the effect of outlier processing by studentized residual detection is as Figure 6 shown.

[0095] (4) Loop through step (3) to traverse the NDVI data for each year in the set LMY, complete the first denoising and reconstruction of the NDVI data for all years, and store the denoising and reconstruction results in the set LM'. In this embodiment, the result of the first denoising and reconstruction of the experimental area of the fused image is as Figure 8 shown.

[0096] (5) Based on the set LM', construct the multi-year synchronous set for each pixel and store it in the NDVI time series set LM''.

[0097] This step includes:

[0098] (5-1) Extract the NDVI data for every ten days within each year for each pixel from the set LM' represents the multi-year synchronous data set of the pixel at row r0 and column l0, represents the NDVI data of the pixel at row r0 and column l0 on the mth and m + 9th days in the yth year;

[0099] (5-2) Store the multi-year synchronous data sets of all pixels in the set

[0100] (6) Read the NDVI data of any pixel from the NDVI time series set LM ″″ and perform secondary denoising and reconstruction.

[0101] As Figure 3 shown, this step specifically includes:

[0102] (6-1) Read the multi-year synchronous data set of any pixel from LM ″″ and

[0103] (6-2) Eliminate all missing values in and store the processed data in the set

[0104] (6-3) Statistic The number N of NDVI data. If N ≤ 3, the process is skipped; if 3 < N ≤ 30, the first outlier test method is executed; if N ≥ 30, the second outlier test method is executed. Replace the detected outliers with the average of other NDVI data except the outliers to complete denoising and reconstruction.

[0105] In this step, the first outlier test method and the second outlier test method are the same as those in step (3). However, since the inter-annual detection has nothing to do with the autocorrelation of the NDVI time series, in the second outlier test method, the set LM ″″″ is first sorted in ascending order so that the Savitzky-Golay function polynomial method can obtain the fitting data.

[0106] In this embodiment, for a pixel set of the multi-year synchronous fusion image, the denoising and reconstruction effect is as Figure 7 shown.

[0107] (7) Loop through step (6) to traverse the NDVI data of all pixels in the set LM ″″ to complete the secondary denoising and reconstruction of the NDVI data of all pixels, store the denoising and reconstruction results in the set LMN, write the data in the set LMN into the image, and obtain the denoised NDVI time series image.

[0108] In this embodiment, the first denoising and reconstruction result and the secondary noise detection result of the fusion image experimental area are as Figure 9 shown, and the final result of the experimental area after denoising and reconstruction for each year on the 273rd day is as Figure 10 shown.

[0109] The embodiment of the present invention also provides a computer device to provide services for the implementation of the above method of the present invention. The device may include: a memory storing computer-executable programs; a processor coupled to the memory; the processor calls the computer-executable programs stored in the memory to execute the steps in the method described in Embodiment 1.

[0110] The embodiment of the present invention also provides a storage medium containing computer-executable programs, and the computer-executable programs are used to execute the method of Embodiment 1 when executed by a computer processor.

[0111] The embodiment of the present invention also provides a computer product, such as an app on a mobile phone, a tablet, an installation program on a computer, etc. The product includes computer programs / instructions, and the computer programs / instructions implement the method described in Embodiment 1 when executed by a processor.

[0112] It should be understood that the above embodiments and the descriptions in the specification only illustrate the principles, main features and advantages of the present invention. Without departing from the spirit and scope of the present invention, the present invention will have various changes and improvements, and all such changes and improvements fall within the protection scope of the present invention.

Claims

1. A denoising and reconstruction method for irregular multi-year time series NDVI data, characterized in that: The steps include: (1) Read the NDVI time series data of each pixel in the multi-year fusion image and store it in the NDVI data set LM; (2) Sort the NDVI data in the NDVI data set LM according to the number of days in the year to obtain the NDVI time series set LMY; (3) Read the NDVI data of any year from the NDVI time series set LMY, store it in the set LMY1, and perform the first denoising and reconstruction of the set LMY1; (4) Loop through step (3), traverse the NDVI data of each year in the set LMY, complete the first denoising and reconstruction of the NDVI data of all years, and store the denoising and reconstruction results in the set LM′; (5) Based on the set LM′, construct a multi-year set for each pixel and store it in the NDVI time series set LM″; (6) Read the NDVI data of any pixel from the NDVI time series set LM″ and perform secondary denoising and reconstruction; (7) Loop through step (6), traverse the NDVI data of all pixels in the set LM″, complete the secondary denoising and reconstruction of the NDVI data of all pixels, store the denoising and reconstruction results in the set LMN, write the data in the set LMN into the image, and obtain the denoised NDVI time series image.

2. The denoising and reconstruction method for irregular multi-year time series NDVI data according to claim 1, characterized in that: Step (3) includes: (3-1) Read the NDVI data of any year from the NDVI time series set LMY and store it in the set in, represents the NDVI data of the pixel in the i-th row and j-th column of the d-th day in the y0-th year in the set LMY1, where k, g, and f represent the number of days, rows, and columns, respectively; (3-2) Take any element from the set LMY1 Record The indices of all missing values ​​in , check The length of all continuous missing value segments in ; if the length of a continuous missing value segment is greater than or equal to 10, the head of the continuous missing value segment is used as the segmentation point to separate it from the previous valid data segment, and the tail is used as the segmentation point to separate it from the subsequent valid data segment, so as to separate the continuous missing value segment into an independent missing value data segment; if the length is less than 10, no processing is done, and finally Divide into several valid data segments and missing value data segments, and store them into the set LMY2; Randomly select an element lmy2 from the set LMY2 u . If all the data in the element lmy2 u are missing values, skip the processing; otherwise, count the number n of NDVI data of the element lmy2 u . If n ≤ 3, skip the processing; if 3 < n ≤ 30, execute the first outlier test method, if n > 30, execute the second outlier test method, and correct the detected outliers to the mean of the two adjacent NDVI data. If the outlier is at the head or tail position, the outlier is corrected to the NDVI data at the adjacent position, thus completing the denoising and reconstruction of lmy2 u . (3-4) Loop through step (3-3) until all elements of the set LMY2 are traversed. Denoising and reconstruction; (3-5) Loop through steps (3-2) to (3-4) until all elements in the set LMY1 are traversed, the first denoising and reconstruction of all NDVI data in the set LMY1 is completed, and the processed data segments are stored in the set LM′.

3. The denoising and reconstruction method for irregular multi-year time series NDVI data according to claim 2 is characterized in that: The specific steps of the first outlier detection method in step (3.3) include: (3-3-1-1) for lmy2 u The NDVI data in are arranged in ascending order to obtain the NDVI data sequence {ndvi1,…,ndvi v ,…,ndvi n },ndvi1,ndvi v ,ndvi n Respectively represent the t1,t v ,t n Daily NDVI data; (3-3-1-2) Set the significance level α and the critical value Q; (3-3-1-3) Calculate NDVI n The difference between ndvi1 and d1 = (ndvi n -ndvi1); (3-3-1-4) If lmy2 u If there is high-end noise in the n With NDVI n-1 The difference d2 = (ndvi n -ndvi n-1 ), and calculate the ratio Q of d2 to d1 c ; (3-3-1-5) If lmy2 u If there is low-end noise in the circuit, the difference between ndvi2 and ndvi1 is calculated as d3 = (ndvi2-ndvi1), and the ratio Q of d3 to d1 is calculated. d ; (3-3-1-6) If Q c >Q, then determine lmy2 u NDVI data in China n is an outlier; if Q d >Q, then determine lmy2 u The NDVI data ndvi1 is an abnormal value.

4. The denoising and reconstruction method for irregular multi-year time series NDVI data according to claim 2, characterized in that: The specific steps of the second outlier detection method in step (3.3) include: (3-3-2-1) for lmy2 u , use the Savitzky-Golay polynomial function method to fit and get the NDVI fitting data; (3-3-2-2) Based on lmy2 u The annual daily data T corresponding to the NDVI data in the text is used to construct the hat matrix H according to the following formula: H=X(X′X) -1 X′ Where, ndvi1,ndvi v ,ndvi n Respectively represent the t1,t v ,t n NDVI data for the day, X is the design matrix, the first column of X is a vector of all 1s, and the second column is lmy2 u The annual daily data corresponding to the NDVI data in the figure; X′ represents the transpose of X; the diagonal element h of H vv Indicates ndvi v The weight of the data in the fitting process; (3-3-2-3) Based on lmy2 u The studentized deleted residuals were calculated from the NDVI data and the NDVI fitting data according to the following formula; Where, t v represents the studentized deleted residual, p represents the number of times the Savitzky-Golay polynomial function method is selected, and r v Residual, ndvi v lmy2 u The vth NDVI data in Indicates the NDVI fitting data with ndvi v The corresponding value of , SSE means the residual sum of squares; (3-3-2-4) Calculate the threshold BC according to the following formula: In the formula, b is the significance value, Indicates that the degrees of freedom are np-1 and the significance level is The critical value of the t-distribution when ; (3-3-2-5) Determine for lmy2 u each NDVI data ndvi in v , whether it satisfies -e*BC < ndvi v < e*BC, where e is the correction factor. If so, then determine that ndvi v is an outlier.

5. The denoising and reconstruction method for irregular multi-year time series NDVI data according to claim 2, characterized in that: After correction, lmy2 u Use Savitzky-Golay filtering.

6. The denoising and reconstruction method for irregular multi-year time series NDVI data according to claim 1, characterized in that: Step (5) comprises: (5-1) Extract the NDVI data of each pixel every ten days for many years from the set LM′ represents the multi-year dataset of the pixel in row r0 and column l0. It represents the NDVI data of the pixel in the r0th row and the l0th column on the mth and m+9th day in the yth year, k represents the total number of days, and t represents the total number of years; (5-2) Store all multi-year datasets of all pixels into a collection Among them, g and f represent the total number of rows and columns of the fused image, respectively.

7. The denoising and reconstruction method for irregular multi-year time series NDVI data according to claim 1, characterized in that: Step (6) comprises: (6-1) Read multi-year datasets of any pixel from LM″ (6-2) Elimination All missing values ​​in the data are stored in the collection (6-3) Statistics Count the number N of NDVI data. If N ≤ 3, skip the processing; if 3 < N ≤ 30, execute the first outlier test method; if N ≥ 30, execute the second outlier test method; replace the detected outliers with the average value of other NDVI data except the outliers to complete denoising and reconstruction.

8. A computer device comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein: The processor executes the computer program to implement the method according to any one of claims 1 to 7.

9. A computer-readable storage medium having a computer program / instruction stored thereon, characterized in that: The computer program / instructions, when executed by a processor, implement the method of any one of claims 1-7.

10. A computer program product comprising a computer program / instructions, characterized in that When the computer program / instructions are executed by a processor, the method according to any one of claims 1 to 7 is implemented.

Images