Remote sensing data fusion method based on bidirectional iteration

Abstract

The invention relates to a remote sensing data fusion method based on bidirectional iteration, and the method comprises the steps: carrying out the time sequence matching of a high-time-resolution remote sensing image and a high-spatial-resolution remote sensing image through dynamic time warping, and constructing a reference image pair; based on the variation trend of the high-temporal-resolution remote sensing image, predicting the high-temporal-resolution remote sensing image at the target moment by using linear regression, and correcting the high-temporal-resolution remote sensing images at the reference moment and the target moment through weighted regression; based on the change trend of the corrected high-temporal-resolution remote sensing image, reversely deducing the high-spatial-resolution remote sensing image at the reference moment by using linear regression, and performing difference comparison with the high-spatial-resolution remote sensing image at the original reference moment; and continuously reducing the difference through an iterative correction mode until a convergence condition is met, and obtaining a fusion sequence. According to the method, the precision of the predicted image is continuously optimized through bidirectional iterative correction, and a more accurate high-temporal-spatial-resolution image can be provided in a high dynamic environment.

Description

Technical Field

[0001] The invention relates to the field of remote sensing data processing and fusion, and in particular to a remote sensing data fusion method based on bidirectional iteration. Background Art

[0002] The need for remote sensing data fusion stems from the trade-off between the temporal and spatial resolutions of single remote sensing data. MODIS data has high temporal resolution but low spatial resolution, while Landsat data has high spatial resolution but low temporal resolution. For applications that need to monitor dynamic changes in the surface at high temporal and spatial scales, such as changes in vegetation cover, soil moisture, and hydrological processes, data from a single sensor is often difficult to meet the needs, so it is necessary to obtain images with higher temporal and spatial resolution by fusing different types of remote sensing data.

[0003] Existing remote sensing data fusion methods can be roughly divided into two categories: direct strategies and indirect strategies. The direct strategy introduces auxiliary data, such as land cover and land use databases, to help decompose NDVI data into pixels to achieve the purpose of downscaling. Representative direct strategies include models such as WLM (Busetto et al., 2008), STAVFM (Meng et al., 2013), and NDVI-LMGM (Rao et al., 2015). In addition, machine learning-based methods are gradually being applied to remote sensing data fusion, such as random forests and deep learning techniques, which can capture complex nonlinear relationships, but often require a large amount of labeled data, which is challenging in highly dynamic areas.

[0004] Indirect strategies use higher-resolution remote sensing data (such as Landsat) to decompose lower-resolution remote sensing data (such as MODIS) to generate images with high temporal and spatial resolution. Representative methods include STARFM (Gao et al., 2006), semi-physical fusion method (Roy et al., 2008), STARCH (Hilker et al., 2009), ESTARFM (Zhu et al., 2010), STDAF (Wu et al., 2012) and ISTDFA (Wu et al., 2018). These models usually assume that the type of land object does not change much during the prediction period, which is difficult to apply in areas where the type of land object changes significantly, such as wetlands and agricultural areas, resulting in insufficient accuracy of the fused image. Summary of the invention

[0005] The purpose of the present invention is to overcome the deficiencies of the above-mentioned prior art and provide a remote sensing data fusion method based on bidirectional iteration.

[0006] To achieve the above object, the present invention adopts the following technical solutions:

[0007] The remote sensing data fusion method based on bidirectional iteration includes:

[0008] S1 acquires a first remote sensing image and a second remote sensing image respectively, wherein the first remote sensing image is a remote sensing image with high temporal resolution, and the second remote sensing image is a remote sensing image with high spatial resolution;

[0009] S2 resamples the first remote sensing image to the spatial resolution of the second remote sensing image, and matches the time series of the first remote sensing image and the second remote sensing image by a dynamic time warping matching method to obtain matching data as a reference data pair;

[0010] S3 combines the reference data pair and linear regression to predict the second remote sensing image data at the target time using the trend of the first remote sensing image time series; then uses the second remote sensing image data at the target time to correct the first remote sensing image data at the target time through weighted regression, and uses the second remote sensing image data at the reference time to correct the first remote sensing image data at the reference time through weighted regression;

[0011] S4 reversely calculates the second remote sensing image data at the reference time based on the time trend of the corrected first remote sensing image and the second remote sensing image data predicted for the target period, and if the difference between the reversely calculated data and the second remote sensing image data at the original reference time exceeds a preset threshold, re-predicts the second remote sensing image data for the target period by using the corrected time trend and the linear regression, and corrects the first remote sensing image data for the target period by using the weighted regression;

[0012] Repeat S4 until the convergence condition is met to obtain the fused sequence.

[0013] As a preferred implementation, the method further includes deleting the matching results in the matching data where the water level difference is greater than a critical value, and the retained data is used as a reference data pair; the critical value is a preset value. Preferably, the critical value is 0.5 m.

[0014] As a preferred implementation, the method of combining the reference data pair and linear regression to predict the second remote sensing image data at the target time using the trend of the first remote sensing image time series is:

[0015] The first remote sensing image data at the reference moment and the first remote sensing image data at the target moment are obtained, and their ratio is calculated as the trend value; the second remote sensing image data matching the first remote sensing image data at the reference moment is obtained, and the product of the second remote sensing image data and the trend value is calculated, i.e., the second remote sensing image data at the target moment.

[0016] As a preferred implementation, the method of correcting data by weighted regression is:

[0017] Setting a sliding window, calculating the spectral difference and Euclidean distance from the central pixel of the first remote sensing image data to all pixels of the second remote sensing image data within the sliding window, and selecting pixels that are similar to the central pixel in time and space;

[0018] For each of the central pixels, a final prediction value is generated using a filtering operator according to the selected similar pixels as a correction value of the central pixel.

[0019] Furthermore, when the number of similar pixels is insufficient, the number of similar pixels is increased by expanding the size of the sliding window.

[0020] Furthermore, in S4, the difference change is calculated after each iteration to ensure that the image difference gradually decreases and the iteration is terminated when a convergence condition is met or a maximum number of iterations is exceeded; the convergence condition is: the image difference is less than a preset threshold.

[0021] Furthermore, the image difference is characterized by the root mean square error, standard deviation or coefficient of variation of pixel values.

[0022] As a preferred implementation, the monitoring frequency of the first remote sensing image is ≤8d, and the spatial resolution of the second remote sensing image is less than or equal to 30m.

[0023] The present invention matches the time series of remote sensing data with high temporal resolution and high spatial resolution based on dynamic time warping (DTW), so that the temporal characteristics of the two can be better aligned, thereby reducing the uncertainty caused by changes in the types of objects. Through bidirectional iterative correction, the accuracy of the predicted image is continuously optimized, especially in areas where data is missing or changes dramatically. The present invention effectively makes up for the limitations of existing fusion methods and models, and can provide more accurate high temporal and spatial resolution images in high dynamic environments, providing strong data support for applications such as environmental monitoring and resource management. BRIEF DESCRIPTION OF THE DRAWINGS

[0024] Figure 1 It is the location and geographical overview of the study area in the embodiment.

[0025] Figure 2 is a flow chart of the method of the present invention.

[0026] Figure 3 It is the image comparison before and after fusion.

[0027] Figure 4 It is a comparison of water area extracted from fused data and Landsat NDWI. DETAILED DESCRIPTION

[0028] The specific implementation of the present invention is further described in detail below in conjunction with the accompanying drawings and examples. The following implementation cases are used to illustrate the present invention, but are not used to limit the scope of the present invention.

[0029] The embodiment of the present invention takes a sub-basin of the Poyang Lake Basin as an example to further illustrate the method of the present invention.

[0030] like Figure 1 As shown in Figure 2, the study area is Poyang Lake and its floodplain wetlands, with a total area of ​​approximately 3618 km 2 There are two national nature reserves and 77 saucer-shaped sub-lakes scattered across the study area. The strong spatial heterogeneity and dramatic dry-wet changes have brought unprecedented challenges to water body identification and water resource estimation.

[0031] like Figure 2 FIG. 1 is a technical flow chart of the present invention, and this embodiment includes the following steps:

[0032] Step 1: Data preparation. Download MODIS MOD13Q1 product data and Landsat data from 2000 to 2016, use RMMEH (Running Median, Mean, Maximum, Endpoint handling, and Hanning) to filter MODIS data, then resample MODIS data to the resolution of Landsat data to obtain reconstructed MODIS NDVI data. Calculate Landsat NDVI for Landsat data.

[0033] The MODIS data products are declouded and de-noised by RMMEH filtering, and the maximum value calculation formula of the original value, mean and median is run: ; Mean The calculation formula is: ; Median The calculation formula is: ; Finally, the entire time series is further smoothed by weighted moving average (Hanning filter) to remove the remaining noise. The formula is: .

[0034] Step 2: Dynamic time warping (DTW): Match the time series of MODIS and Landsat data to determine the optimal correspondence between them and provide a “reference image pair” for the next step of predicting Landsat in the target period.

[0035] Through dynamic programming technology, the time series is "stretched" or "compressed" to find an optimal matching path to maximize the similarity between them and reduce the matching error. The recursive formula is:

[0036] ;

[0037] Is to reach the position The cumulative minimum distance of Indicates the first and The distance between elements; Represents the minimum cumulative distance among adjacent locations.

[0038] Finally, the matching results with water level difference > 0.5 m were deleted to ensure the relative consistency of the dry-wet distribution of the “reference image pair”; the water level difference here refers to the difference between the observed water level value corresponding to the observation date of the MODIS image and the observed water level value corresponding to the observation date of the matching Landsat image.

[0039] Step 3: Linear regression prediction: Use the trend of MODIS time series to predict the Landsat data at the target time.

[0040] The resampled MOD13Q1 NDVI data and Landsat NDVI theoretically have the same change trend. Step 3 is to integrate the time information of MODIS data. ;trend Where: is the Landsat NDVI value at the predicted target moment; and are the MODIS NDVI values ​​at the target time and the base period, respectively; The base period Landsat NDVI registered by dynamic time warping in step 2. Long-term changes in land features will lead to trends The previous method used the "image pair" at a certain historical moment as the only benchmark. The change of ground objects will cause the trend R of MODIS and the trend R' of Landsat in the previous and next periods to be inconsistent, resulting in The prediction result is inaccurate. The present invention can match the "contemporary" images to the maximum extent through step 2, thereby reducing or avoiding trend distortion caused by changes in ground objects.

[0041] Step 4: Weighted regression prediction. Using Landsat data at the target time Correct the MODIS data at the target time get , and using the baseline Landsat data Correction of MODIS data in the baseline period get The steps of weighted regression prediction are as follows:

[0042] Set a sliding window, calculate the spectral difference (NDVI value difference can be evaluated based on conventional indicators such as root mean square error, standard deviation or coefficient of variation) and Euclidean distance from the MODIS central pixel to all Landsat pixels in the sliding window, and select pixels that are similar in time and space to the target pixel (i.e., the MODIS central pixel) through a hierarchical strategy. When the number of similar pixels is insufficient, the number of similar pixels can be increased by expanding the size of the sliding window to ensure that the model has enough data for fusion and prediction. For each target pixel, the model uses a filtering operator to generate the final predicted NDVI value based on the selected set of similar pixels. The filtering operator combines spatial and spectral information, and weights the spectral difference and distance of similar pixels to calculate the corrected value of the target pixel:

[0043]

[0044]

[0045]

[0046]

[0047] In the formula, represents the size (side length) of the moving window; b represents the band; is the combination weight determined by the spectrum and the distance difference; Indicates spectral differences; represents the Euclidean distance; is a set of similar pixels represented by a binary matrix; , , Respectively represent the coordinates ( ) in the reference period, the NDVI value of the MODIS pixel in the reference period, and the NDVI value of the MODIS pixel at the target moment (the next moment in the reference period).

[0048] Step 5: Bidirectional iteration.

[0049] Based on the time trend of the modified MODIS ( ) and the Landsat data predicted for the target period ( ), and the Landsat data of the base period were calculated by linear regression ( ),if and The larger difference is based on The Landsat of the target period is re-predicted by linear regression, and the MODIS of the target period is corrected by weighted regression to obtain , repeat step 5, The numerator is iterated until the error ( and The difference of ) converges to an acceptable range.

[0050] Step 6: Introduce regularization and convergence judgment. Calculate the error change after each iteration to ensure that the error gradually decreases and terminate the iteration when the convergence condition is met.

[0051] In the bidirectional iteration process, regularization and convergence judgment are introduced through step 6. After each iteration, the difference between the current error and the error of the previous iteration is calculated. If the relative change is less than the set threshold ,Right now: , the error is considered to be stable and converged. Set a maximum number of iterations If convergence is not achieved after reaching this number, the iteration is forced to stop. Finally, the fused data combining MODIS time series and Landsat spatial information is obtained.

[0052] like Figure 3 As shown, the fusion results of 2004, 2008, 2012 and 2016 are extracted at intervals of four years for display. In order to illustrate that the present invention has universality and robustness under changing hydrological conditions, the fusion results of the dry season, flood season, flood season and the middle period of the receding season are selected for four representative years. Figure 3 In the figure, the first and third columns show the MOD13Q1 original data products of Poyang Lake National Nature Reserve and Nanji Wetland National Nature Reserve before fusion. The second and fourth columns show the fused NDVI products of the two reserves for the same period. The fused NDVI can distinguish smaller water surfaces; it has significantly improved the depiction of lake waterside lines and the continuity of rivers.

[0053] Figure 4 The correlation between the water area of ​​Poyang Lake extracted based on Landsat Normalized Water Index (NDWI) and the water area of ​​Poyang Lake extracted based on the present invention is shown. The linear fitting coefficient R 2 Reached 0.92, the root mean square error (RMSE) was 207.76 km 2 , the mean absolute error is 135.6 km 2 , indicating that the fused data has a high consistency with Landsat in spatial information. In addition to spatial information, the present invention also integrates the high-frequency temporal information of MODIS. The NDVI sequence reconstructed by the present invention can achieve 8-day, 30-meter near-real-time ground observation, making up for the shortcomings of single-source remote sensing images in temporal resolution or spatial resolution.

Claims

1. A remote sensing data fusion method based on bidirectional iteration, characterized in that: include: S1 acquires a first remote sensing image and a second remote sensing image respectively, wherein the first remote sensing image is a remote sensing image with high temporal resolution, and the second remote sensing image is a remote sensing image with high spatial resolution; S2 resamples the first remote sensing image to the spatial resolution of the second remote sensing image, and matches the time series of the first remote sensing image and the second remote sensing image by a dynamic time warping matching method to obtain matching data as a reference data pair; S3 combines the reference data pair and linear regression to predict the second remote sensing image data at the target time using the trend of the first remote sensing image time series; then uses the second remote sensing image data at the target time to correct the first remote sensing image data at the target time through weighted regression, and uses the second remote sensing image data at the reference time to correct the first remote sensing image data at the reference time through weighted regression; S4 reversely calculates the second remote sensing image data at the reference time based on the time trend of the corrected first remote sensing image and the second remote sensing image data predicted for the target period, and if the difference between the reversely calculated data and the second remote sensing image data at the original reference time exceeds a preset threshold, re-predicts the second remote sensing image data for the target period by using the corrected time trend and the linear regression, and corrects the first remote sensing image data for the target period by using the weighted regression; Repeat S4 until the convergence condition is met to obtain the fused sequence.

2. The method according to claim 1, characterized in that: It also includes deleting matching results in which the water level difference in the matching data is greater than a critical value, and the retained data is used as a reference data pair; the critical value is a preset value.

3. The method according to claim 2, characterized in that The critical value is 0.5m.

4. The method according to claim 1, characterized in that: The method of combining the reference data pair and linear regression to predict the second remote sensing image data at the target time using the trend of the first remote sensing image time series is: The first remote sensing image data at the reference moment and the first remote sensing image data at the target moment are obtained, and their ratio is calculated as the trend value; the second remote sensing image data matching the first remote sensing image data at the reference moment is obtained, and the product of the second remote sensing image data and the trend value is calculated, i.e., the predicted value of the second remote sensing image data at the target moment is calculated.

5. The method according to claim 1, characterized in that The method of correcting data by weighted regression is: Setting a sliding window, calculating the spectral difference and Euclidean distance from the central pixel of the first remote sensing image data to all pixels of the second remote sensing image data within the sliding window, and selecting pixels that are similar to the central pixel in time and space; For each of the central pixels, a final prediction value is generated using a filtering operator according to the selected similar pixels as a correction value of the central pixel.

6. The method according to claim 5, characterized in that When the number of similar pixels is insufficient, the number of similar pixels is increased by expanding the size of the sliding window.

7. The method according to claim 1, characterized in that In S4, the difference change is calculated after each iteration to ensure that the difference gradually decreases and the iteration is terminated when a convergence condition is met or a maximum number of iterations is exceeded; the convergence condition is: the difference is less than a preset threshold.

8. The method according to claim 7, characterized in that The differences are characterized by the root mean square error, standard deviation or coefficient of variation of the pixel values.

9. The method according to claim 1, characterized in that: The monitoring frequency of the first remote sensing image is ≤8d, and the spatial resolution of the second remote sensing image is ≤30m.

Images