Remote sensing data fusion method based on bidirectional iteration

Abstract

The present application relates to a remote sensing data fusion method based on bidirectional iteration, time series matching of high temporal resolution remote sensing images and high spatial resolution remote sensing images is carried out through dynamic time warping to construct a reference image pair; based on the change trend of high temporal resolution remote sensing images, linear regression is used to predict high spatial resolution remote sensing images at the target time, and the high temporal resolution remote sensing images at the reference time and the target time are corrected through weighted regression; based on the change trend of the corrected high temporal resolution remote sensing images, linear regression is used to inversely deduce the high spatial resolution remote sensing images at the reference time, and difference comparison is carried out with the original high spatial resolution remote sensing images at the reference time; the difference is continuously reduced through iterative correction until the convergence condition is met, and the fusion sequence is obtained. Through bidirectional iterative correction, the accuracy of the predicted image is continuously optimized, and more accurate high temporal and spatial resolution images can be provided in a high dynamic environment.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of remote sensing data processing and fusion, and particularly relates to a remote sensing data fusion method based on bidirectional iteration. BACKGROUND

[0002] The demand 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 require monitoring of dynamic changes in the ground surface at high temporal and spatial scales, such as changes in vegetation cover, soil moisture, and hydrological processes, the data from a single sensor often cannot meet the needs, and therefore it is necessary to fuse different types of remote sensing data to obtain higher temporal and spatial resolution images.

[0003] Existing remote sensing data fusion methods can be broadly divided into direct and indirect strategies. Direct strategies introduce auxiliary data, such as land cover and land use databases, to help perform pixel decomposition on NDVI data to achieve the purpose of downscaling. Representative direct strategies include WLM (Busetto et al., 2008), STAVFM (Meng et al., 2013), NDVI-LMGM (Rao et al., 2015), and other models. In addition, machine learning-based methods have also been gradually 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 high temporal and spatial resolution images. Representative methods include STARFM (Gao et al., 2006), semi-physical fusion methods (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), and other models. These models usually assume that the ground object type does not change significantly during the prediction period, which is difficult to apply in areas where the ground object type changes significantly, such as wetlands and agricultural areas, resulting in insufficient accuracy of the fused images. SUMMARY

[0005] The present application aims to overcome the shortcomings of the prior art and provide a remote sensing data fusion method based on bidirectional iteration.

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

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

[0008] S1 respectively acquires a first remote sensing image and a second remote sensing image, the first remote sensing image is a high-time-resolution remote sensing image, and the second remote sensing image is a high-space-resolution remote sensing image;

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

[0010] S3 combines the reference data pair and linear regression, and predicts the second remote sensing image data at a target moment by using the trend of the time sequence of the first remote sensing image; then, the second remote sensing image data at the target moment is used to correct the first remote sensing image data at the target moment through weighted regression, and the second remote sensing image data at the reference moment is used to correct the first remote sensing image data at the reference moment through weighted regression;

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

[0012] S4 is repeated until a convergence condition is met, and a fusion sequence is acquired.

[0013] As a preferred embodiment, the method further comprises deleting the matching results with a water level difference greater than a critical value in the matching data, and retaining the data as the reference data pair; the critical value is a preset value. Preferably, the critical value is 0.5 m.

[0014] As a preferred embodiment, the manner of combining the reference data pair and the linear regression to predict the second remote sensing image data at the target moment by using the trend of the time sequence of the first remote sensing image is as follows:

[0015] The first remote sensing image data at the reference moment and the first remote sensing image data at the target moment are acquired, and a ratio of the first remote sensing image data at the reference moment to the first remote sensing image data at the target moment is calculated as a trend value; the second remote sensing image data matched with the first remote sensing image data at the reference moment is acquired, and a product of the second remote sensing image data and the trend value is calculated, that is, the second remote sensing image data at the target moment.

[0016] As a preferred embodiment, the manner of correcting the data through weighted regression is as follows:

[0017] a sliding window is set, the spectral difference and Euclidean distance between the center pixel of the first remote sensing image data and all pixels in the sliding window of the second remote sensing image data are calculated, and the pixels similar to the center pixel in time and space are selected;

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

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

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

[0021] Further, the image difference is represented by the root mean square error, standard deviation or coefficient of variation of the pixel value.

[0022] As a preferred embodiment, 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 application matches the time series of high temporal resolution and high spatial resolution remote sensing data based on dynamic time warping (DTW), so that the time characteristics of the two can be better aligned, thereby reducing the uncertainty caused by changes in ground object types. Through bidirectional iterative correction, the accuracy of the predicted image is continuously optimized, especially in areas with missing data or drastic changes. The present application effectively compensates for the limitations of existing fusion methods and models, and can provide more accurate high temporal and spatial resolution images in a high dynamic environment, providing strong data support for environmental monitoring, resource management and other applications. BRIEF DESCRIPTION OF DRAWINGS

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

[0025] Figure 2 is a flowchart of the method of the present application.

[0026] Figure 3 is a comparison of images before and after fusion.

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

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

[0029] The method of the present application is further described in the embodiments of the present application, taking a sub-basin of Poyang Lake Basin as an example.

[0030] As shown in Figure 1 , the study area is Poyang Lake and its floodplain wetlands, with a total area of about 3618 km 2 . The study area is distributed with two national nature reserves and 77 disc-shaped sub-lakes. Strong spatial heterogeneity and severe wet-dry changes bring unprecedented challenges to water body identification and water resource estimation.

[0031] As shown in Figure 2 , it is a technical flowchart of the present application, and the present embodiment includes the following steps:

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

[0033] The cloud and noise of the MODIS data product are removed by RMMEH filtering, and the maximum value calculation formula of the original value, the mean value and the median value is run: ; the calculation formula of the mean value : ; the calculation formula of the median value : ; finally, the entire time series is further smoothed by weighted moving average (Hanning filter) to remove the remaining noise, and the formula is: .

[0034] Step 2: Dynamic time warping matching DTW (dynamic time warping). 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 the target period of Landsat.

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

[0036] ;

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

[0038] Finally, the matching results with water level differences > 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 MODIS image observation date and the observed water level value corresponding to the matching Landsat image observation date.

[0039] Step 3: Linear regression prediction: Use the trend of the 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 baseline period, respectively; The base period Landsat NDVI is registered by dynamic time warping in step 2. The 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 present invention can maximize the matching of “contemporaneous” images 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 Corrected base period MODIS data 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 within 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] Where, represents the size (side length) of the moving window; b represents the band; is the combination weight determined by the spectrum and distance difference; Indicates spectral differences; represents the Euclidean distance; is a set of similar pixels represented by a binary matrix; 、 、 Represents 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: Iterate in both directions.

[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 baseline period were calculated by linear regression ( ),if and The larger the difference, the The Landsat of the target period was re-forecasted using linear regression, and the MODIS of the target period was corrected by weighted regression to obtain , repeat step 5, the molecules are iterated until the error (E) andthe difference) 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 is gradually reduced and terminate the iteration when the convergence condition is met.

[0051] In the bidirectional iteration process, regularization and convergence judgment are introduced by step 6. Calculate the difference between the current error and the error of the previous iteration after each iteration. If the relative change amplitude is less than the set threshold , i.e. , it is considered that the error has stabilized and reached convergence. Set a maximum number of iterations , if it has not converged after reaching the number, force stop the iteration. Finally, obtain the fusion data combined with MODIS time series and Landsat spatial information.

[0052] As shown in Figure 3 , the fusion results of 2004, 2008, 2012 and 2016 are extracted at four-year intervals for display. To illustrate the universality and robustness of the present application under changing hydrological conditions, the fusion results of the middle period of the dry season, the rising water period, the wet season and the recession period are selected respectively. Figure 3 Among them, the first and third columns respectively show the MOD13Q1 original data product before fusion in Poyang Lake National Nature Reserve and Nanji Wetland National Nature Reserve. The second and fourth columns respectively show the NDVI products of the two reserves after fusion at the same time. The fused NDVI can distinguish smaller water area; it has significantly improved in lake waterline delineation and river continuity preservation.

[0053] Figure 4 shows the correlation between the Poyang Lake water area extracted based on Landsat normalized water index (NDWI) and the Poyang Lake water area extracted based on the present application. The linear fitting coefficient R 2 of the two reaches 0.92, the root mean square error (RMSE) is 207.76 km 2 , and the average absolute error is 135.6 km 2 , which shows that the fusion data has high consistency with Landsat in spatial information. In addition to spatial information, the present application also fuses the high-frequency time information of MODIS. The NDVI sequence reconstructed by the present application can realize 8-day, 30-meter near real-time earth observation, making up for the shortcomings of single-source remote sensing images in time 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 using 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, using the second remote sensing image data at the target time, corrects the first remote sensing image data at the target time through weighted regression, and simultaneously 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; the method of correcting data through 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; 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. 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-predict the second remote sensing image data for the target period using the corrected time trend and the linear regression, and corrects the first remote sensing image data for the target period using the weighted regression. Repeat S4 until the convergence condition is met and obtain the fused sequence.

2. The method according to claim 1, characterized in that The method further includes deleting matching results in which the water level difference in the matching data is greater than a critical value, and retaining the data 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, wherein 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: Obtain the first remote sensing image data at the reference moment and the first remote sensing image data at the target moment, and calculate their ratio as the trend value; obtain the second remote sensing image data that matches the first remote sensing image data at the reference moment, and calculate the product of the second remote sensing image data and the trend value, that is, the predicted value of the second remote sensing image data at the target moment.

5. The method according to claim 1, wherein When the number of similar pixels is insufficient, the number of similar pixels is increased by expanding the size of the sliding window.

6. 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.

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

8. 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