Multi-temporal remote sensing image fusion method, device and equipment and storage medium

Abstract

The present disclosure provides a multi-temporal remote sensing image fusion method, device, equipment and storage medium. The method of the present disclosure comprises: reading and preprocessing multi-temporal remote sensing image data; performing linear normalization for each band of the multi-temporal remote sensing image data and saving global normalization parameters; initializing an adaptive Kalman filter for each pixel point in the multi-temporal remote sensing image data, the initialization of each adaptive Kalman filter comprising setting a forgetting factor of the adaptive Kalman filter; performing processing using the adaptive Kalman filter to determine the pixel value of each pixel point in the fused image data to obtain fused image data, the processing of the adaptive Kalman filter comprising a prediction operation and an update operation, the update operation comprising adjusting an observation noise covariance matrix using the forgetting factor; de-normalizing the fused image data according to the global normalization parameters; saving the fused image data and setting its geographic information. The present disclosure can improve the multi-temporal remote sensing image fusion quality.

Description

Technical Field

[0001] This disclosure relates to the field of remote sensing image processing technology, specifically to a method, apparatus, device, and storage medium for multi-temporal remote sensing image fusion. Background Technology

[0002] In the field of remote sensing technology, multi-temporal remote sensing imagery contains rich geographic information and is of great significance for geographic monitoring and environmental change research. However, remote sensing images acquired at different times may be affected by sensor errors, illumination variations, and other factors, easily introducing noise. Directly using these images for analysis may affect the accuracy of the results. Among related technologies, image fusion methods suffer from problems such as information loss and lack of dynamic adjustment capabilities when processing multi-temporal remote sensing images. They also cannot adaptively optimize for time-series data and cannot fully leverage the advantages of multi-temporal imagery. Therefore, there is a need for an effective method to fuse multi-temporal remote sensing images to improve image quality. Summary of the Invention

[0003] To address the aforementioned technical problems, embodiments of this disclosure provide a method, apparatus, device, and storage medium for multi-temporal remote sensing image fusion.

[0004] According to one aspect of this disclosure, a multi-temporal remote sensing image fusion method is provided, comprising:

[0005] Read and preprocess multi-temporal remote sensing image data to ensure that all images in the multi-temporal remote sensing image data have the same size and the same band;

[0006] The multi-temporal remote sensing image data is subjected to band-by-band linear normalization and the global normalization parameters are saved.

[0007] An adaptive Kalman filter is initialized for each pixel in the multi-temporal remote sensing image data. The initialization of each adaptive Kalman filter includes setting the forgetting factor of the adaptive Kalman filter.

[0008] The adaptive Kalman filter is used to perform the following processing to determine the pixel value of each pixel in the fused image data, thereby obtaining the fused image data of the multi-temporal remote sensing image data: The measured value of the current pixel in the first image of the multi-temporal remote sensing image data is taken as the initial state estimate of the adaptive Kalman filter; the prediction operation of the adaptive Kalman filter is performed to obtain the state prediction value of the current pixel in the second image; an update operation is performed using the state prediction value of the current pixel in the second image and the measured value of the current pixel in the second image to fuse the state prediction value and the measured value to obtain the state estimate value of the current pixel in the second image; this process is iterated until all images of the multi-temporal remote sensing image data have undergone prediction and update operations, and the final state estimate value is taken as the pixel value of the current pixel in the fused image data. The update operation includes adjusting the observation noise covariance matrix and the process noise covariance matrix using the forgetting factor.

[0009] The fused image data is denormalized according to the global normalization parameters;

[0010] The fused image data is saved and its geographic information is set, including geographic transformation information and projection information.

[0011] In some embodiments of the first aspect of this disclosure, the prediction operation includes: determining the predicted state value at the current time based on the state transition matrix of the adaptive Kalman filter and the state estimate value at the previous time; and determining the predicted state covariance matrix at the current time based on the state transition matrix, the state estimate covariance matrix at the previous time, and the process noise covariance matrix updated at the previous time.

[0012] The update operation includes: calculating the residual and residual covariance matrix at the current time based on the current state prediction value and the current observation value; determining the current observation noise covariance matrix based on the previous observation noise covariance matrix, the forgetting factor, the current residual, the observation matrix of the adaptive Kalman filter, and the current state prediction covariance matrix; determining the current Kalman gain matrix based on the current state prediction covariance matrix, the observation matrix of the adaptive Kalman filter, and the current residual covariance matrix; determining the current state estimate value based on the current state prediction value, the current Kalman gain matrix, and the current residual; determining the current state estimate covariance matrix based on the current Kalman gain matrix and the current state prediction covariance matrix; and determining the updated process noise covariance matrix at the current time based on the forgetting factor, the previous process noise covariance matrix, the current Kalman gain matrix, the current residual, the state transition matrix of the adaptive Kalman filter, the current state estimate covariance matrix, and the previous state estimate covariance matrix.

[0013] In some embodiments of the first aspect of this disclosure, the reading and preprocessing of multi-temporal remote sensing image data includes: when reading the multi-temporal remote sensing image data, using the GDAL library to open the file of each image in the multi-temporal remote sensing image data one by one, and throwing an exception if it cannot be opened; checking the number of bands of each image in the multi-temporal remote sensing image data, using the number of bands of the first image in the multi-temporal remote sensing image data as a standard, and throwing an exception if the number of bands of other images is the same as the standard; cropping each image in the multi-temporal remote sensing image data, using the number of rows and columns of the first image as a standard, cropping all images to the corresponding size of the upper left corner to ensure that all images in the multi-temporal remote sensing image data have the same size.

[0014] In some embodiments of the first aspect of this disclosure, the global normalization parameters are stored using a dictionary data structure.

[0015] In some embodiments of the first aspect of this disclosure, initializing an adaptive Kalman filter for each pixel in the multi-temporal remote sensing image data further includes: when initializing the adaptive Kalman filter for each pixel in the multi-temporal remote sensing image data, setting the state transition matrix, observation matrix, process noise covariance matrix, and measurement noise covariance matrix of the adaptive Kalman filter.

[0016] In some embodiments of the first aspect of this disclosure, saving the fused image data and setting its geographic information includes: when saving the fused image data, using a driver in the GDAL library to create a new image file, setting the number of rows, columns, bands and data type of the image, writing the fused image data into the corresponding band and setting the geographic transformation information and projection information of the image.

[0017] In some embodiments of the first aspect of this disclosure, the formula for calculating the inverse normalization is:

[0018] ,

[0019] in, It is the inverse normalized value of the b-th band. It is the normalized pixel value of the fused image data of the b-th band. It is the global maximum value of the b-th band in the global normalized parameters. It is the global minimum value of the b-th band in the global normalized parameters.

[0020] According to a second aspect of this disclosure, a multi-temporal remote sensing image fusion apparatus is provided, comprising:

[0021] The data reading module is used to read and preprocess multi-temporal remote sensing image data so that all images in the multi-temporal remote sensing image data have the same size and the same band.

[0022] The normalization module is used to perform band-by-band linear normalization on the multi-temporal remote sensing image data and save the global normalization parameters;

[0023] A filtering module is used to initialize an adaptive Kalman filter for each pixel in the multi-temporal remote sensing image data. The initialization of each adaptive Kalman filter includes setting a forgetting factor for the adaptive Kalman filter. The module also performs the following processing using the adaptive Kalman filter to determine the pixel value of each pixel in the fused image data, thereby obtaining fused image data of the multi-temporal remote sensing image data: taking the measured value of the current pixel in the first image of the multi-temporal remote sensing image data as the initial state estimate of the adaptive Kalman filter; performing a prediction operation of the adaptive Kalman filter to obtain the state prediction value of the current pixel in the second image; performing an update operation using the state prediction value of the current pixel in the second image and the measured value of the current pixel in the second image to fuse the state prediction value and the measured value to obtain the state estimate value of the current pixel in the second image; iterating in this way until all images of the multi-temporal remote sensing image data have completed the prediction and update operations, and taking the final state estimate value as the pixel value of the current pixel in the fused image data. The update operation includes adjusting the observation noise covariance matrix and the process noise covariance matrix using the forgetting factor.

[0024] The denormalization module is used to denormalize the fused image data according to the global normalization parameters;

[0025] The image storage module is used to store the denormalized fused image data and set its geographic information, which includes geographic transformation information and projection information.

[0026] According to a third aspect of this disclosure, an electronic device is provided, the electronic device including a processor and a memory, the memory storing a computer program that, when executed by the processor, causes the processor to perform the methods described above.

[0027] According to a fourth aspect of this disclosure, a computer storage medium is provided, on which a computer program is stored, which, when executed by a processor, causes the processor to perform the methods described above.

[0028] This disclosure eliminates the impact of differences in dimensions and numerical ranges between different bands through band-by-band global normalization, ensuring that image data is unified to the same scale, making data processing more stable, and preserving the unique radiometric characteristics of each band. It introduces an adaptive Kalman filter to reduce reliance on prior information from the noise model, introduces a forgetting factor to balance the impact of old and new data on system estimation, and combines prediction residuals to adaptively update state estimation and noise covariance matrix, effectively adapting to real-time noise changes, preserving the temporal variation patterns of the images, effectively utilizing information from multi-temporal images, better reflecting changes in geographic information, and preserving the spatial reference information of the original images. This enables effective fusion of multi-temporal remote sensing images, better preserving useful information in the images, and improving the fusion quality of multi-temporal remote sensing images. This disclosure can be applied to, but is not limited to, remote sensing image processing and geographic monitoring. Attached Figure Description

[0029] Figure 1 This is a flowchart illustrating the multi-temporal remote sensing image fusion method provided in the embodiments of this disclosure;

[0030] Figure 2 This is a schematic diagram of Sentinel-2 multi-temporal multispectral remote sensing images from 2019 and 2020.

[0031] Figure 3 This is a schematic diagram of the 2021 and 2022 images from the Sentinel-2 multi-temporal multispectral remote sensing imagery.

[0032] Figure 4 This is a schematic diagram of the 2023 Sentinel-2 multi-temporal multispectral remote sensing image and a schematic diagram of the fused image obtained by fusing the Sentinel-2 multi-temporal multispectral remote sensing image using the method of this embodiment.

[0033] Figure 5 This is a schematic diagram of the spectral curves of the reference image and the fused image of the Sentinel-2 multi-temporal multispectral remote sensing image plotted on the same pixel, where (a) is a schematic diagram of the spectral curve of the reference image of the Sentinel-2 multi-temporal multispectral remote sensing image.

[0034] (b) is a schematic diagram of the spectral curve of the fused image obtained by fusing Sentinel-2 multi-temporal multispectral remote sensing images using the method of the embodiments of this disclosure;

[0035] Figure 6 This is a schematic diagram of the structure of the multi-temporal remote sensing image fusion device provided in the embodiments of this disclosure;

[0036] Figure 7 This is a schematic structural block diagram of an electronic device provided in an embodiment of this disclosure. Detailed Implementation

[0037] Hereinafter, exemplary embodiments according to the present disclosure will be described in detail with reference to the accompanying drawings. Obviously, the described embodiments are merely some embodiments of the present disclosure, and not all embodiments of the present disclosure, and it should be understood that the present disclosure is not limited to the exemplary embodiments described herein.

[0038] Figure 1 A schematic flowchart of the multi-temporal remote sensing image fusion method provided in this disclosure is shown. See also Figure 1 The multi-temporal remote sensing image fusion method of this disclosure includes the following steps:

[0039] Step 101: Read and preprocess multi-temporal remote sensing image data to ensure that all images in the multi-temporal remote sensing image data have the same size and the same band.

[0040] Step 102: Perform band-by-band linear normalization on the multi-temporal remote sensing image data and save the global normalization parameters;

[0041] Step 103: Initialize an adaptive Kalman filter for each pixel in the multi-temporal remote sensing image data. The initialization of each adaptive Kalman filter includes setting the forgetting factor of the adaptive Kalman filter.

[0042] Step 104: Using an adaptive Kalman filter, perform the following processing to determine the pixel value of each pixel in the fused image data, obtaining fused image data of multi-temporal remote sensing image data: Take the measured value of the current pixel in the first image of the multi-temporal remote sensing image data as the initial state estimate of the adaptive Kalman filter; perform the prediction operation of the adaptive Kalman filter to obtain the state prediction value of the current pixel in the second image; use the state prediction value of the current pixel in the second image and the measured value of the current pixel in the second image to perform an update operation to fuse the state prediction value and the measured value to obtain the state estimate value of the current pixel in the second image; iterate in this way until all images of the multi-temporal remote sensing image data have completed the prediction operation and update operation, and take the final state estimate value as the pixel value of the current pixel in the fused image data, wherein the update operation includes adjusting the observation noise covariance matrix and the process noise covariance matrix using the forgetting factor;

[0043] Step 105: Denormalize and fuse the image data according to the global normalization parameters;

[0044] Step 106: Save the fused image data and set its geographic information, which includes geographic transformation information and projection information.

[0045] In this embodiment of the disclosure, the multi-temporal remote sensing image data includes data from multiple images of the same area acquired at different times. That is, all images in the multi-temporal remote sensing image data cover the same geographical area, and each image was captured at a different time.

[0046] The following uses Sentinel-2 multitemporal multispectral remote sensing imagery (containing 12 bands) as an example to illustrate the specific implementation of this disclosure. The Sentinel-2 multitemporal multispectral remote sensing imagery spans five years from 2019 to 2023, with a spatial resolution of 10 meters and an image size of 520*480. Figure 2 The imagery from 2019 and 2020 is shown in the Sentinel-2 multi-temporal multispectral remote sensing imagery. Figure 3 The imagery shown is from Sentinel-2 multi-temporal multispectral remote sensing, including images from 2021 and 2022. Figure 4 The image shown is from Sentinel-2 multitemporal multispectral remote sensing imagery for 2023, and the fused image obtained by fusing Sentinel-2 multitemporal multispectral remote sensing imagery using the method of embodiments of this disclosure.

[0047] In some examples, in step 101, when reading multi-temporal remote sensing image data, the Geospatial Data Abstraction Library (GDAL) is used to open the file of each image in the multi-temporal remote sensing image data one by one. If the file cannot be opened, an exception is thrown to ensure data availability. The number of bands in each image in the multi-temporal remote sensing image data is checked, using the number of bands in the first image in the multi-temporal remote sensing image data as a standard. If the number of bands in other images does not match the standard, an exception is thrown to ensure the consistency of the image data. Each image in the multi-temporal remote sensing image data is cropped, using the number of rows and columns of the first image as a standard, and all images are cropped to the corresponding size from the upper left corner to ensure that all images in the multi-temporal remote sensing image data have the same size, so that subsequent processing can be performed uniformly. Thus, by reading and preprocessing multi-temporal remote sensing image data, a unified data format can be provided for subsequent fusion processing.

[0048] Taking Sentinel-2 multi-temporal multispectral remote sensing imagery (containing 12 bands) as an example, in step 101, the GDAL library can be used to read remote sensing image data from five periods. When checking the number of bands, the number of bands in the first image is used as the standard. If the number of bands in other images differs from this, the entire system throws a ValueError exception, prompting "Image band count mismatch, please re-enter remote sensing imagery," and terminates the system operation, preventing subsequent fusion based on an adaptive Kalman filter to ensure the consistency of the image data involved in the fusion. When cropping images, the number of rows and columns in the first image is used as the standard, and the upper left corner area of ​​520*480 pixels is cropped from all images to ensure that subsequent processing can be performed uniformly and that the spatial range is consistent.

[0049] In step 102, for each band, the global maximum and minimum values ​​of all images to be fused (i.e., all images in the preprocessed multi-temporal remote sensing image data) in that band are linearly normalized and the global normalization parameters are saved to ensure that all data are unified to the same scale so that the data can be processed stably.

[0050] In some examples, the formula for calculating band-by-band global normalization is as follows:

[0051] ,

[0052] in, It is the globally normalized value of the b-th band. These are all pixel values ​​in the b-th band. It is the global minimum value of the b-th band. It is the global maximum value of the b-th band. It is a very small value, in order to avoid the case where the denominator is 0.

[0053] in, It represents all pixel values ​​in band b of the image, extracted directly from the remote sensing image data. It is the global minimum value of band b. We count all pixel values ​​and take the smallest value. It is the global maximum value of band b. It is calculated by counting all pixel values ​​in band_data_b and taking the largest value. To avoid the minimum value where the denominator is 0, a preset fixed value can be used, with a range between 10⁻⁶ and 10⁻¹⁰.

[0054] 'b' represents the number of bands, which depends on the sensor type of the remote sensing image. For example, RGB imagery has 3 bands. Sentinel-2 imagery typically has 12 bands, and LandSat imagery typically has 7 bands. The aforementioned global normalization formula allows for the same normalization operation to be performed on each band. The parameters corresponding to different bands (i.e.,...) , , The values ​​of ) are completely different, but the constant ϵ can remain consistent. That is, regardless of the band, All pixel values ​​in this band are taken. Each represents the global minimum value of that band. Each represents the global maximum value of its respective band, while ϵ represents the minimum value to avoid a denominator of 0. The original pixel value ranges differ across bands (e.g., infrared pixel values ​​are typically higher than visible light band values), so their maximum / minimum value statistics are completely independent.

[0055] When performing global normalization on all images in multi-temporal remote sensing image data, this embodiment of the disclosure adds an extremely small value. This can effectively avoid situations where the denominator is 0.

[0056] Furthermore, in step 102, the global normalization parameters can be stored using a dictionary data structure for easy use during denormalization in subsequent step 105.

[0057] In step 103, for each pixel, i.e. each geographic location, an adaptive Kalman filter is initialized for each pixel in the multi-temporal remote sensing image data. That is, an independent adaptive Kalman filter is initialized for each geographic coordinate location in the multi-temporal remote sensing image data that spans all temporal phases of the complete time series.

[0058] In step 103, when initializing the adaptive Kalman filter for each pixel in the multi-temporal remote sensing image data, in addition to setting the forgetting factor of the adaptive Kalman filter, the state transition matrix, observation matrix, process noise covariance matrix and measurement noise covariance matrix of the adaptive Kalman filter are also set.

[0059] Specifically, the state transition matrix, observation matrix, process noise covariance matrix, measurement noise covariance matrix, and forgetting factor of the adaptive Kalman filter can be set according to the actual situation. For example, the state transition matrix of the adaptive Kalman filter can be assumed to be a simple uniform velocity model, the observation matrix can map the state to the measurement space, and the process noise covariance matrix, measurement noise covariance matrix, and forgetting factor can be adjusted according to the actual noise situation. The initial state estimate can be set as a zero vector. Such settings can reasonably initialize the adaptive Kalman filter according to the characteristics and noise situation of the image, providing a foundation for subsequent filtering processing.

[0060] Taking Sentinel-2 multi-temporal multispectral remote sensing imagery (containing 12 bands) as an example, the state vector dimension is set to 2, and the state vector contains pixel values ​​and their rates of change, and the state transition matrix... Assuming a uniform velocity model with a time interval dt = 1 (year), the observation matrix... The state is mapped to the measurement space. The process noise covariance matrix and the measurement noise covariance matrix can be adjusted according to the actual noise situation. The forgetting factor can be set to 0.99, and the initial state estimate can be set to the zero vector.

[0061] Adaptive Kalman filtering is an improvement on Kalman filtering, mainly in that it adaptively updates the noise covariance matrix to adapt to real-time noise changes. Adaptive Kalman filtering describes the relationship between the system state and the observed values ​​based on state equations and observation equations. Its state equation is shown below:

[0062] ,

[0063] in, is the state vector at the current moment; F is the state transition matrix of the adaptive Kalman filter, which describes the dynamic change of the system state from the previous moment to the current moment; It is the state vector from the previous moment; It is the process noise vector, which follows a Gaussian distribution.

[0064] The observation equation for adaptive Kalman filtering is shown below:

[0065] ,

[0066] in, is the observation vector at the current moment, i.e., the sensor's measurement value; H is the observation matrix of the adaptive Kalman filter, which maps the state to the observation space; It is the observation noise vector, which follows a Gaussian distribution.

[0067] In step 104, the prediction operation may include: determining the predicted state value at the current time based on the state transition matrix of the adaptive Kalman filter and the state estimate value at the previous time step; and determining the predicted state covariance matrix at the current time step based on the state transition matrix, the state estimate covariance matrix at the previous time step, and the updated process noise covariance matrix at the previous time step. That is, the prediction process of the adaptive Kalman filter in step 104 predicts the state at the current time step based on the state at the previous time step, and includes two steps: state prediction and state covariance prediction.

[0068] In some examples, the formula for calculating the state prediction value is as follows:

[0069] ,

[0070] in, It is the predicted state value at the current moment; It is the state transition matrix of the adaptive Kalman filter, which describes the dynamic change of the system state from the previous time step to the current time step; It is the state estimate from the previous moment.

[0071] In some examples, the formula for calculating the state prediction covariance matrix is ​​as follows:

[0072] ,

[0073] in, It is the state prediction covariance matrix at the current moment, representing the uncertainty in estimating the current state; It is the state transition matrix of the adaptive Kalman filter; It is the state estimation covariance matrix of the previous time step; It is the transpose of the state transition matrix; It is the adaptively updated process noise covariance matrix of the previous time step, representing the uncertainty caused by noise factors during the system's state transition process.

[0074] In step 104, the update operation of the adaptive Kalman filter may include: calculating the residual and residual covariance, calculating the Kalman gain matrix at the current time, adaptively adjusting the noise matrix (i.e., the observation noise covariance matrix and the process noise covariance matrix) in combination with the forgetting factor, and updating the state estimate and the state estimate covariance matrix.

[0075] Specifically, the update operation may include: calculating the residual and residual covariance matrix at the current time based on the current state prediction value and the current observation value; determining the current observation noise covariance matrix based on the previous observation noise covariance matrix, the forgetting factor, the current residual, the observation matrix of the adaptive Kalman filter, and the current state prediction covariance matrix; determining the current Kalman gain matrix based on the current state prediction covariance matrix, the observation matrix of the adaptive Kalman filter, and the current residual covariance matrix; determining the current state estimate value based on the current state prediction value, the current Kalman gain matrix, and the current residual; determining the current state estimate covariance matrix based on the current Kalman gain matrix and the current state prediction covariance matrix; and determining the updated process noise covariance matrix at the current time based on the forgetting factor, the previous process noise covariance matrix, the current Kalman gain matrix, the current residual, the state transition matrix of the adaptive Kalman filter, the current state estimate covariance matrix, and the previous state estimate covariance matrix.

[0076] In some examples, the residuals are calculated using the following formulas:

[0077] ,

[0078] in, It is the residual at the current moment, representing the difference between the actual measured value and the predicted value; It is the observation value at the current moment; It is the observation matrix of the adaptive Kalman filter; It is the predicted state value at the current moment.

[0079] In some examples, the formula for calculating residual covariance is as follows:

[0080] ,

[0081] in, It is the residual covariance matrix at the current moment; It is the observation matrix of the adaptive Kalman filter; It is the state prediction covariance matrix at the current moment; It is the observation noise covariance matrix updated at the previous time step.

[0082] In some examples, the formula for calculating the observation noise covariance matrix is ​​as follows:

[0083] ,

[0084] in, It is the observation noise covariance matrix at the current moment; It is the forgetting factor of the adaptive Kalman filter, used to adjust the weights of new and old data; It is the residual at the current moment; It is the observation matrix; It is the state prediction covariance matrix at the current moment.

[0085] In some examples, the Kalman gain matrix is ​​calculated using the following formula:

[0086] ,

[0087] in, It is the Kalman gain matrix, which can balance the uncertainty of the predicted and observed values ​​and assign weights to them to obtain a more accurate state estimate; It is the state prediction covariance matrix at the current moment; It is the transpose of the observation matrix; It is the inverse of the residual covariance matrix at the current moment.

[0088] In some examples, the update formula for the state estimate is as follows:

[0089] ,

[0090] in, It is the state estimate at the current moment; It is the predicted state value at the current moment; It is the Kalman gain matrix at the current moment, used to weigh the predicted and measured values; It is the residual at the current moment.

[0091] In some examples, the updated covariance of the state estimate is calculated as follows:

[0092] ,

[0093] in, It is the covariance matrix of the current state estimate; It is the identity matrix; It is the Kalman gain matrix; It is the state prediction covariance matrix at the current moment.

[0094] In some examples, the update formula for the process noise covariance matrix is ​​as follows:

[0095] ,

[0096] in, It is the process noise covariance matrix updated at the current moment; It is a forgetting factor; It is the process noise covariance matrix updated in the previous time step; It is the Kalman gain matrix; It is the residual at the current moment; It is the state transition matrix; It is the covariance matrix of the current state estimate; It is the state estimation covariance matrix of the previous time step.

[0097] Taking the Sentinel-2 multi-temporal multispectral remote sensing image (containing 12 bands) as an example, in step 104, the adaptive Kalman filter of each pixel is used to perform the following processing to determine the pixel value of each pixel in the fused image data, thereby completing its fusion: the pixel value of the current pixel in the first image, i.e., the 2019 image, is used as the initial state estimate of the current pixel, and the rate of change is initialized to 0. The 12 bands of each pixel are processed separately. Taking the 2020 image as an example, the prediction stage includes: predicting the state observation value and state observation of 2020 based on the 2019 state estimate. The covariance matrix update phase includes: inputting the observed value of the current pixel in 2020, calculating the prediction residual by combining it with the state prediction value, calculating the residual covariance matrix, updating the observation noise covariance matrix, calculating the Kalman gain matrix, performing state update to obtain the state estimate of the current pixel in 2020, updating the state estimate covariance matrix to obtain the state estimate covariance matrix of the current pixel in 2020, and updating the process noise covariance matrix; the aforementioned prediction and update phases are executed iteratively until iterates to 2023, and the final state estimate is taken as the pixel value of the fused image at the current pixel.

[0098] As can be seen above, the pixel value of each pixel in the fused image data is determined by iteratively performing prediction and update using an adaptive Kalman filter for each pixel, fully utilizing the information in the multi-temporal images and improving the accuracy of fusion. Thus, the fusion of multi-temporal, multispectral remote sensing images is achieved through adaptive Kalman filtering.

[0099] In this embodiment, the measured value of the first image is used as the initial state estimate. For each pixel in the image, the measured value of that pixel in each image is sequentially input into an adaptive Kalman filter for prediction. Then, the prediction residual is calculated based on the predicted value and the actual observed value. A forgetting factor is introduced to balance the influence of new and old data on the current estimate. The process noise covariance matrix and the observation noise covariance matrix are dynamically adjusted in combination to adapt to the real-time changes in noise.

[0100] In step 105, inverse normalization is performed based on the global normalization parameters during normalization, which can restore the fused image data to the original data range.

[0101] In some examples, the formula for calculating inverse normalization is as follows:

[0102] ,

[0103] in, It is the inverse normalized value of the b-th band. It is the normalized pixel value of the fused image data of the b-th band. It is the global maximum value of the b-th band in the global normalized parameters. It is the global minimum value of the b-th band in the global normalized parameters.

[0104] In step 106, when saving the fused image data, a new image file can be created using the driver in the GDAL library. The number of rows, columns, bands, and data type of the image can be set. The fused image data is then written to the corresponding bands, and the geographic transformation and projection information of the image is set. This ensures the accuracy and usability of the fused image in geospatial context, while also ensuring that the fused image is correctly saved and has accurate geographic information, facilitating subsequent applications and analysis.

[0105] Taking Sentinel-2 multi-temporal multispectral remote sensing imagery (containing 12 bands) as an example, the driver in the GDAL library can be used to create new image files, set the number of rows, columns, bands, and data types of the image, write the fused image data to the corresponding bands, and set the geographic transformation information and WGS84 projection information of the image to ensure the accuracy and usability of the fused imagery in geospatial context.

[0106] Taking the Sentinel-2 multi-temporal multispectral remote sensing image (containing 12 bands) as an example, in order to verify the effectiveness of the multi-temporal remote sensing image fusion method provided in this disclosure, the peak signal-to-noise ratio (PSNR), root mean square error (RMSE), structural similarity index (SSIM), spectral angle (SAM), and correlation coefficient (CC) are introduced as verification indicators to quantitatively evaluate the fused image data obtained through this disclosure. At the same time, the spectral curves of the reference image and the fused image are plotted on the same pixel to compare and analyze the effectiveness of the fused image data.

[0107] Peak signal-to-noise ratio (PSNR) measures the noise level between a fused image and a reference image. A higher PSNR indicates that the fused image retains more signal information, i.e., useful information in the fused image, and also means that the noise level in the fused image is lower, making the image look clearer.

[0108] The formula for calculating PSNR is as follows:

[0109] ,

[0110] in, The maximum pixel value. Mean squared error is used to assess the difference between predicted and observed values.

[0111] in, The calculation formula is as follows:

[0112] ,

[0113] in, Indicates the number of rows in the image. Indicates the number of columns in the image; Indicates the location of the reference image pixel values, Indicates the location of the fused image The pixel value.

[0114] The root mean square error (RMSE) reflects the pixel-level deviation between the fused image and the reference image. The smaller the value, the smaller the error between the fused image and the reference image, indicating that the fusion accuracy is higher.

[0115] The formula for calculating RMSE is as follows:

[0116] ,

[0117] in, This represents the mean square error.

[0118] The Structural Similarity Index (SSIM) can be used to evaluate the similarity of two images (i.e., the fused image and the reference image) in terms of brightness, contrast, and structure, providing a more comprehensive assessment of the image quality of the fused image. A SSIM value closer to 1 indicates a higher structural similarity between the fused image and the reference image.

[0119] The formula for calculating SSIM is as follows:

[0120] ,

[0121] in, , These are the average brightness values ​​(i.e., average pixel values) of the reference image and the fused image, respectively. , These are the variances of the reference image and the fused image, respectively. It is the covariance between the brightness values ​​of two images, representing the correlation between their brightness changes; , It is a small constant introduced to avoid the denominator being zero.

[0122] Spectral angle (SAM) measures the similarity of spectral information between two images (i.e., a fused image and a reference image), and can be used to assess spectral fidelity. The smaller the spectral angle (SAM), the lower the spectral distortion.

[0123] The formula for calculating the spectral angle (SAM) is as follows:

[0124] ,

[0125] in, This indicates the angle between two images (i.e., the merged image and the reference image); This represents the spectral reflectance of the reference image in the i-th band; The spectral reflectance of the fused image in the i-th band is represented; n is the total number of bands; arccos() is the inverse cosine function used to calculate the included angle.

[0126] The correlation coefficient (CC) measures the linear correlation between a reference image and a fused image. Its value ranges from -1 to 1, where 1 indicates a perfect positive correlation, -1 indicates a perfect negative correlation, and 0 indicates no linear correlation between the two variables. A CC closer to 1 indicates a more consistent trend in pixel value changes between the two images. The formula for calculating the correlation coefficient (CC) is as follows:

[0127] ,

[0128] in, This represents the i-th observation of the reference image; express The average value of the variable; This represents the i-th observation in the fused image; express The average value of the variable; n represents the total number of observations.

[0129] Taking the Sentinel-2 multi-temporal multispectral remote sensing image (containing 12 bands) as an example, the above indicators are calculated and shown in Table 1 below.

[0130] Table 1:

[0131] Reference Image Fusion Image PSNR RMSE SSIM SAM CC 2023 fused 32.57 0.085 0.91 0.079 0.94

[0132] As can be seen from Table 1, the fused image obtained by fusing Sentinel-2 multi-temporal multispectral remote sensing images using the method of the present invention has high quality.

[0133] Figure 5 This diagram illustrates the spectral curves of the reference image and the fused image of the Sentinel-2 multi-temporal multispectral remote sensing image plotted on the same pixel block. Figure 5(a) shows a schematic diagram of the spectral curves of the reference image for the Sentinel-2 multi-temporal multispectral remote sensing image. Figure 5 (b) is a schematic diagram of the spectral curve of the fused image obtained by fusing Sentinel-2 multi-temporal multispectral remote sensing images using the method of this embodiment. See also Figure 5 (a) and Figure 5 In (b), the spectral curve of the fused image shows a more complex shape, revealing more details and subtle variations, with more peaks and troughs. This indicates that the fused image combines the spectral information of multiple images at this pixel, and can reflect a more detailed material composition and richer spectral features.

[0134] The beneficial effects of the embodiments disclosed herein include at least the following:

[0135] 1) By using band-by-band global normalization, it is possible to ensure that image data are unified to the same scale, eliminate the influence of differences in dimensions and numerical ranges between different bands, make data processing more stable, and preserve the unique radiometric characteristics of each band, thereby improving the quality of fused images.

[0136] 2) The introduction of a forgetting factor balances the impact of old and new data on system estimation. Combined with the prediction residual, the state estimation and noise covariance matrix are adaptively updated. This effectively adapts to the real-time changes in noise, preserves the temporal variation of the images, effectively utilizes the information in multi-temporal images, better reflects the changes in geographic information, and preserves the spatial reference information of the original images, which is convenient for subsequent analysis and application.

[0137] 3) It has good versatility and can be widely used in remote sensing image processing, geographic monitoring, environmental research and other fields, providing high-quality image data support for related research and applications.

[0138] Figure 6 A schematic diagram of the device provided in an embodiment of this disclosure is shown. See also Figure 6 The multi-temporal remote sensing image fusion apparatus provided in this disclosure embodiment may include:

[0139] The data reading module 601 is used to read and preprocess multi-temporal remote sensing image data so that all images in the multi-temporal remote sensing image data have the same size and the same band.

[0140] The normalization module 602 is used to perform band-by-band linear normalization on multi-temporal remote sensing image data and save the global normalization parameters.

[0141] The filtering module 603 is used to initialize an adaptive Kalman filter for each pixel in the multi-temporal remote sensing image data. The initialization of each adaptive Kalman filter includes setting a forgetting factor for the adaptive Kalman filter. The module then uses the adaptive Kalman filter to perform the following processing to determine the pixel value of each pixel in the fused image data, thereby obtaining fused image data of the multi-temporal remote sensing image data: taking the measured value of the current pixel in the first image of the multi-temporal remote sensing image data as the initial state estimate of the adaptive Kalman filter; performing a prediction operation of the adaptive Kalman filter to obtain the state prediction value of the current pixel in the second image; using the state prediction value of the current pixel in the second image and the measured value of the current pixel in the second image to perform an update operation to fuse the state prediction value and the measured value to obtain the state estimate value of the current pixel in the second image; this process is iterated until all images in the multi-temporal remote sensing image data have completed the prediction and update operations, and the final state estimate value is taken as the pixel value of the current pixel in the fused image data. The update operation includes adjusting the observation noise covariance matrix and the process noise covariance matrix using the forgetting factor.

[0142] The denormalization module 604 is used to denormalize the fused image data according to the global normalization parameters.

[0143] The image storage module 605 is used to store the denormalized fused image data and set its geographic information, which includes geographic transformation information and projection information.

[0144] Furthermore, the filtering module 603 can be specifically used to perform prediction and update operations in the manner described in the aforementioned method.

[0145] Furthermore, the data reading module 601 can be specifically used for: when reading multi-temporal remote sensing image data, using the GDAL library to open the file of each image in the multi-temporal remote sensing image data one by one, and throwing an exception if it cannot be opened; checking the number of bands of each image in the multi-temporal remote sensing image data, using the number of bands of the first image in the multi-temporal remote sensing image data as the standard, and throwing an exception if the number of bands of other images is not the standard; cropping each image in the multi-temporal remote sensing image data, using the number of rows and columns of the first image as the standard, cropping all images to the corresponding size of the upper left corner, to ensure that all images in the multi-temporal remote sensing image data have the same size.

[0146] Furthermore, the normalization module 602 can be specifically used to: store global normalization parameters using a dictionary data structure.

[0147] Furthermore, the filtering module 603 can be specifically used to: set the state transition matrix, observation matrix, process noise covariance matrix and measurement noise covariance matrix of the adaptive Kalman filter when initializing the adaptive Kalman filter for each pixel in the multi-temporal remote sensing image data.

[0148] Furthermore, the image storage module 605 can be used to create new image files using the driver in the GDAL library when storing fused image data, set the number of rows, columns, bands and data types of the image, write the fused image data to the corresponding bands and set the geographic transformation information and projection information of the image.

[0149] In practical applications, multi-temporal remote sensing image fusion devices can be implemented as software, hardware, or a combination of both. In some examples, multi-temporal remote sensing image fusion devices can be implemented as follows: Figure 7 The electronic device shown or the software deployed in the electronic device.

[0150] Figure 7 A schematic structural diagram of an electronic device provided according to an embodiment of this disclosure is shown. See also... Figure 7 The electronic device 700 provided in this embodiment may include a processor 701 and a memory 702. The memory 702 stores a computer program, which, when run by the processor 701, causes the processor 701 to execute the aforementioned multi-temporal remote sensing image fusion method.

[0151] Electronic device 700 may include one or more processors 701 and memory 702.

[0152] The processor 701 may be a central processing unit (CPU) or other form of processing unit with data processing capabilities and / or instruction execution capabilities, and may control other components in the electronic device to perform desired functions.

[0153] The memory 702 may include one or more computer program products, which may include various forms of computer-readable storage media, such as volatile memory and / or non-volatile memory. The volatile memory may include, for example, random access memory (RAM) and / or cache memory. The non-volatile memory may include, for example, read-only memory (ROM), hard disk, flash memory, etc. One or more computer program instructions may be stored on the computer-readable storage medium, and a processor may execute the program instructions to implement the methods of the various embodiments of this disclosure described above and / or other desired functions.

[0154] Depending on the specific application, the electronic device 700 may also include any other suitable components.

[0155] In addition to the methods and devices described above, embodiments of this disclosure may also be computer program products, including computer program instructions that, when executed by a processor, cause the processor to perform the steps in the multi-temporal remote sensing image fusion method of embodiments of this disclosure.

[0156] The computer program product can be written in any combination of one or more programming languages ​​to perform the operations of the embodiments of this disclosure. The programming languages ​​include object-oriented programming languages ​​such as Java and C++, as well as conventional procedural programming languages ​​such as C or similar languages. The program code can be executed entirely on a user's computing device, partially on a user's computing device, as a standalone software package, partially on a user's computing device and partially on a remote computing device, or entirely on a remote computing device or server.

[0157] Furthermore, embodiments of this disclosure may also be computer-readable storage media storing a computer program thereon, which, when run by a processor, causes the processor to execute the steps in the multi-temporal remote sensing image fusion method of embodiments of this disclosure.

[0158] The computer-readable storage medium may be any combination of one or more readable media. A readable medium may be a readable signal medium or a readable storage medium. A readable storage medium may be, for example, an electrical, magnetic, optical, electromagnetic, infrared, or semiconductor system, apparatus, or device, or any combination thereof. More specific examples of readable storage media (a non-exhaustive list) include: an electrical connection having one or more wires, a portable disk, a hard disk, random access memory (RAM), read-only memory (ROM), erasable programmable read-only memory (EPROM or flash memory), optical fiber, portable compact disk read-only memory (CD-ROM), optical storage device, magnetic storage device, or any suitable combination thereof.

[0159] The above description has been given for purposes of illustration and description. Furthermore, this description is not intended to limit the embodiments of this disclosure to the forms disclosed herein. Although numerous exemplary aspects and embodiments have been discussed above, those skilled in the art will recognize certain variations, modifications, alterations, additions, and sub-combinations therein.

Claims

1. A method for fusing multi-temporal remote sensing images, characterized in that, include: Read and preprocess multi-temporal remote sensing image data to ensure that all images in the multi-temporal remote sensing image data have the same size and the same band; The multi-temporal remote sensing image data is subjected to band-by-band linear normalization and the global normalization parameters are saved. An adaptive Kalman filter is initialized for each pixel in the multi-temporal remote sensing image data. The initialization of each adaptive Kalman filter includes setting the forgetting factor of the adaptive Kalman filter. The adaptive Kalman filter is used to perform the following processing to determine the pixel value of each pixel in the fused image data, thereby obtaining the fused image data of the multi-temporal remote sensing image data: The measured value of the current pixel in the first image of the multi-temporal remote sensing image data is taken as the initial state estimate of the adaptive Kalman filter; the prediction operation of the adaptive Kalman filter is performed to obtain the state prediction value of the current pixel in the second image; an update operation is performed using the state prediction value of the current pixel in the second image and the measured value of the current pixel in the second image to fuse the state prediction value and the measured value to obtain the state estimate value of the current pixel in the second image; this process is iterated until all images of the multi-temporal remote sensing image data have undergone prediction and update operations, and the final state estimate value is taken as the pixel value of the current pixel in the fused image data. The update operation includes adjusting the observation noise covariance matrix and the process noise covariance matrix using the forgetting factor. The fused image data is denormalized according to the global normalization parameters; The fused image data is saved and its geographic information is set, including geographic transformation information and projection information.

2. The method according to claim 1, characterized in that, The prediction operation includes: determining the predicted state value at the current time based on the state transition matrix of the adaptive Kalman filter and the state estimate value at the previous time step; and determining the predicted state covariance matrix at the current time step based on the state transition matrix, the state estimate covariance matrix at the previous time step, and the process noise covariance matrix updated at the previous time step. The update operation includes: Calculate the residuals and residual covariance matrix at the current time based on the predicted state values ​​and the observed values ​​at the current time. The observation noise covariance matrix at the current time is determined based on the observation noise covariance matrix at the previous time step, the forgetting factor, the residual at the current time step, the observation matrix of the adaptive Kalman filter, and the state prediction covariance matrix at the current time step. The Kalman gain matrix at the current time is determined based on the state prediction covariance matrix at the current time, the observation matrix of the adaptive Kalman filter, and the residual covariance matrix at the current time. The state estimate at the current time is determined based on the predicted state value at the current time, the Kalman gain matrix at the current time, and the residual at the current time. The state estimation covariance matrix at the current time is determined based on the Kalman gain matrix and the state prediction covariance matrix at the current time. The updated process noise covariance matrix at the current time is determined based on the forgetting factor, the process noise covariance matrix at the previous time, the Kalman gain matrix at the current time, the residual at the current time, the state transition matrix of the adaptive Kalman filter, the state estimation covariance matrix at the current time, and the state estimation covariance matrix at the previous time.

3. The method according to claim 1, characterized in that, The reading and preprocessing of multi-temporal remote sensing image data includes: when reading the multi-temporal remote sensing image data, using the GDAL library to open the file of each image in the multi-temporal remote sensing image data one by one, and throwing an exception if it cannot be opened; checking the number of bands of each image in the multi-temporal remote sensing image data, using the number of bands of the first image in the multi-temporal remote sensing image data as the standard, and throwing an exception if the number of bands of other images does not match the standard; cropping each image in the multi-temporal remote sensing image data, using the number of rows and columns of the first image as the standard, cropping all images to the corresponding size of the upper left corner to ensure that all images in the multi-temporal remote sensing image data have the same size.

4. The method according to claim 1, characterized in that, The global normalization parameters are stored using a dictionary data structure.

5. The method according to claim 1, characterized in that, The step of initializing an adaptive Kalman filter for each pixel in the multi-temporal remote sensing image data further includes: when initializing the adaptive Kalman filter for each pixel in the multi-temporal remote sensing image data, setting the state transition matrix, observation matrix, process noise covariance matrix, and measurement noise covariance matrix of the adaptive Kalman filter.

6. The method according to claim 1, characterized in that, Saving the fused image data and setting its geographic information includes: when saving the fused image data, using the driver in the GDAL library to create a new image file, setting the number of rows, columns, bands and data types of the image, writing the fused image data into the corresponding bands and setting the geographic transformation information and projection information of the image.

7. The method according to claim 1, characterized in that, The formula for calculating the inverse normalization is: ， in, It is the inverse normalized value of the b-th band. It is the normalized pixel value of the fused image data of the b-th band. It is the global maximum value of the b-th band in the global normalized parameters. It is the global minimum value of the b-th band in the global normalized parameters.

8. A multi-temporal remote sensing image fusion device, comprising: The data reading module is used to read and preprocess multi-temporal remote sensing image data so that all images in the multi-temporal remote sensing image data have the same size and the same band. The normalization module is used to perform band-by-band linear normalization on the multi-temporal remote sensing image data and save the global normalization parameters; A filtering module is used to initialize an adaptive Kalman filter for each pixel in the multi-temporal remote sensing image data. The initialization of each adaptive Kalman filter includes setting a forgetting factor for the adaptive Kalman filter. The module also performs the following processing using the adaptive Kalman filter to determine the pixel value of each pixel in the fused image data, thereby obtaining fused image data of the multi-temporal remote sensing image data: taking the measured value of the current pixel in the first image of the multi-temporal remote sensing image data as the initial state estimate of the adaptive Kalman filter; performing a prediction operation of the adaptive Kalman filter to obtain the state prediction value of the current pixel in the second image; performing an update operation using the state prediction value of the current pixel in the second image and the measured value of the current pixel in the second image to fuse the state prediction value and the measured value to obtain the state estimate value of the current pixel in the second image; iterating in this way until all images of the multi-temporal remote sensing image data have completed the prediction and update operations, and taking the final state estimate value as the pixel value of the current pixel in the fused image data. The update operation includes adjusting the observation noise covariance matrix and the process noise covariance matrix using the forgetting factor. The denormalization module is used to denormalize the fused image data according to the global normalization parameters; The image storage module is used to store the denormalized fused image data and set its geographic information, which includes geographic transformation information and projection information.

9. An electronic device, characterized in that, The electronic device includes a processor and a memory, the memory storing a computer program that, when executed by the processor, causes the processor to perform the method as described in any one of claims 1-7.

10. A computer storage medium, characterized in that, The computer storage medium stores a computer program, which, when executed by a processor, causes the processor to perform the method as described in any one of claims 1-7.

Images