Accumulated snow distribution prediction method and system based on double-branch model

Abstract

The invention discloses an accumulated snow distribution prediction method and system based on a double-branch model, and belongs to the technical field of remote sensing image processing, and the system comprises an accumulated snow data preprocessing module, an accumulated snow distribution information extraction module, a frequency domain interaction module, a multi-modal information fusion module, and an accumulated snow distribution condition prediction module. The accumulated snow data preprocessing module is used for processing remote sensing data, and the accumulated snow distribution information extraction module is used for extracting feature information related to accumulated snow in multi-source remote sensing data; the frequency domain interaction module is used for performing frequency domain analysis on the fused remote sensing data; the multi-modal information fusion module is used for integrating multi-source remote sensing data from different types of sensors and carrying out deep fusion on the multi-source remote sensing data; and the accumulated snow distribution condition prediction module is used for predicting the accumulated snow distribution conditions in different time periods by means of the multi-source remote sensing data fused by the multi-modal information module.

Description

Technical Field

[0001] The present invention belongs to the technical field of remote sensing image processing, and particularly relates to a snow cover distribution prediction method and system based on a dual-branch model. Background Art

[0002] Snow cover, as an essential part of the cryosphere, is not only a core variable for studying climate dynamics and the hydrological cycle but also a crucial link in the climate system. Its unique heat insulation effect and high albedo characteristics directly and profoundly affect the energy balance of the Earth's surface. Meanwhile, snow cover plays a pivotal role in the water cycle, and snowmelt runoff provides a stable and reliable water source for many water-scarce regions. However, the distribution of snow cover is not static. It is affected by multiple climate and geographical factors and shows complex and variable characteristics. Therefore, accurately obtaining and predicting the time series of snow cover distribution is of great significance for deeply understanding the climate system, optimizing water resource management, and preventing natural disasters.

[0003] Since the 1970s, satellite remote sensing with large-area synchronous observation has become an important tool for snow cover observation. Based on the characteristics that snow cover strongly reflects in the visible / near-infrared band and strongly absorbs in the short-wave infrared band, researchers widely use methods such as the ratio of each band of sensors such as Landsat and MODIS, and the normalized difference snow index (NDSI) to study the change of snow cover extent. However, because it is difficult to distinguish wet snow and dry snow from optical remote sensing data, and the radiation signal received by the satellite is also affected by weather factors such as clouds and fog, the identification and prediction related to snow cover are restricted.

[0004] Snow cover has a wide range and strong spatial heterogeneity in its distribution, which makes it difficult for a small number of stations to fully display the spatio-temporal variation characteristics of snow cover on a large spatial scale. The application of satellite observation technology effectively makes up for the deficiency of the insufficient number of stations in traditional in-situ snow measurement methods and can provide long-term and large-scale snow detection data. So far, it has brought great convenience to the research in fields such as hydrology, meteorology, earthquake prevention and disaster reduction, and climate change. Currently, optical remote sensing sensors are mainly used for remote sensing monitoring of snow cover area, snow reflectivity, etc., but it is difficult to effectively estimate parameters such as snow depth and snow water equivalent. Moreover, optical remote sensing is greatly affected by weather conditions, especially cloud cover, which seriously affects the accuracy and effectiveness of detection. In contrast, microwaves have the ability to penetrate clouds and snow cover layers and can detect information on the underlying surface of snow cover, featuring all-weather and large-scale snow detection.

[0005] Currently, for the realization of snow cover time series prediction, there are mainly two key engineering problems:

[0006] (1) The occlusion of clouds has a significant impact on the observation of snow cover distribution. When making snow cover predictions, in order to maintain temporal coherence, it is necessary to include image data that is partially or completely occluded by clouds. So, how to make full use of the complementary information of synthetic aperture radar (SAR) and multi-spectral optical (opt) to provide a richer and more accurate information basis for time series prediction;

[0007] (2) In the process of making temporal predictions of snow cover, it is crucial to fully integrate and utilize the snow cover-related information from two different modalities of synthetic aperture radar (SAR) and multi-spectral optical (opt), which will lay a solid foundation for the snow cover prediction task. The key lies in how to effectively maintain the continuity information of snow cover in the temporal and spatial dimensions, so as to ensure more accurate and reliable prediction results. Summary of the Invention

[0008] The technical problem to be solved by the present invention is to effectively overcome the influence of cloud occlusion on optical remote sensing data and maintain the continuity of snow cover in the temporal and spatial dimensions; the present invention proposes a method and system for predicting snow cover distribution, aiming to effectively solve the deficiencies in the prior art and improve the accuracy and reliability of snow cover distribution prediction. A method for predicting snow cover distribution according to the present invention includes the following steps:

[0009] Step 1, obtain remote sensing data of synthetic aperture radar (SAR) images and multi-spectral optical (opt) images, and preprocess the remote sensing data of the synthetic aperture radar (SAR) images and multi-spectral optical (opt) images to obtain preprocessed synthetic aperture radar (SAR) image data and multi-spectral optical (opt) image data;

[0010] Step 2, extract snow cover area image data based on the preprocessed synthetic aperture radar (SAR) image data ;

[0011] Extract snow cover area image data based on the preprocessed multi-spectral optical (opt) image data ;

[0012] Step 3, perform feature fusion on the extracted snow cover area image data and the snow cover area image data to obtain fused image features;

[0013] Step 4, predict the snow cover distribution based on the fused image features.

[0014] Further, the preprocessing of the remote sensing data of the synthetic aperture radar (SAR) images and multi-spectral optical (opt) images in Step 1 specifically includes:

[0015] Step 1.1, perform orbital correction, radiometric calibration, speckle noise filtering, geometric correction, and conversion of the backscattering coefficient to decibel format on the SAR image data in sequence to obtain the SAR image data after the first preprocessing;

[0016] Step 1.2, process the opt image data, specifically:

[0017] First, perform band selection on the opt image data to obtain the selected bands , expressed as:

[0018]

[0019] where, is the blue light band, is the green light band, is the red light band, is the near-infrared band, is the short-wave infrared band;

[0020] Perform radiometric calibration on the selected bands , and then perform weighted average fusion on the selected and to obtain the opt image data after the first preprocessing, expressed as:

[0021] where, is the image after band fusion, is the weight of the th band, is the band , , and after radiometric calibration;

[0022] Step 1.3, perform image registration on the SAR image data after the first preprocessing in Step 1.1 and the opt image data after the first preprocessing in Step 1.2, and resample the image after image registration using the Cubic Convolution method to obtain the synthetic aperture radar SAR image after preprocessing in Step 1 and the preprocessed multi-spectral optical opt image.

[0023] Furthermore, the extraction of snow-covered area image data from the multi-spectral optical opt image data described in Step 2 , specifically:

[0024] Use the green light band and the short-wave infrared band to calculate the Normalized Difference Snow Index NDSI:

[0025]

[0026] Among them, and are the reflectance in the green light band and the short-wave infrared band respectively.

[0027] Set a threshold to extract the snow cover T1, and extract the snow cover area according to the NDSI value to obtain the image data of the snow cover area :

[0028] .

[0029] Furthermore, step 3 includes the following steps:

[0030] Step 3.1, extract the shallow features of the snow cover area image data and the snow cover area image data respectively, which are expressed as:

[0031]

[0032] Among them, is the encoder SFE, is the shallow feature extracted from the synthetic aperture radar (SAR) image , is the shallow feature extracted from the multi-spectral optical (opt) image .

[0033] Step 3.2, extract the low-frequency basic features from the shallow features, which are expressed as:

[0034]

[0035] Among them represents the BTE encoder, represents the low-frequency basic feature extracted from the shallow feature , represents the low-frequency basic feature extracted from the shallow feature ;

[0036] Step 3.3, extract the high-frequency detailed features from the shallow features, which are expressed as:

[0037]

[0038] Among them, represents the encoder DCE, represents the high-frequency detailed feature extracted from the shallow feature , Represents that from shallow features The high-frequency detail features extracted from

[0039] Step 3.4, fusing the low-frequency basic features with the high-frequency detail features to obtain fused image features.

[0040] Furthermore, the step 3.4 is as follows:

[0041] Step 3.4.1: The snow area image data The low-frequency basic characteristics and high frequency detail features Convert the time domain to the frequency domain and convert the snow area image data The low-frequency basic characteristics and high frequency detail features Converting from the time domain to the frequency domain, it is represented as follows:

[0042]

[0043]

[0044] in, represents a two-dimensional Fourier transform operation, and They are the basic frequency domain features of optical and synthetic aperture radar SAR images, and They are the detailed frequency domain features of optical and synthetic aperture radar SAR images;

[0045] Step 3.4.2: In the frequency domain, the basic features of the multispectral optical opt image are transformed into Basic characteristics of synthetic aperture radar (SAR) images Fusion and detailed features of multi-spectral optical opt images and SAR image details The fusion is expressed as:

[0046]

[0047]

[0048] in, represents the fusion operation with channel attention mechanism, represents the basic features fused in the frequency domain, Represents the detailed features fused in the frequency domain;

[0049] Step 3.4.3, fusion features in the frequency domain and Convert it back to the time domain again, so as to obtain the fused basic features and the fused detailed features in the time domain.

[0050]

[0051] Among them, represents the two-dimensional inverse Fourier transform operation, represents the fused basic features in the time domain, represents the fused detailed features in the time domain;

[0052] Step 3.4.4, input the fused basic features and the fused detailed features in the time domain into the decoder to generate the finally fused image features , which is expressed as:

[0053]

[0054] Among them, is the decoded fused feature, represents the decoder.

[0055] Furthermore, in step 4, based on the fused image features, predict the snow cover distribution;

[0056] Specifically, a dual-branch model is used to predict the snow cover distribution; the dual-branch model includes a ConvLSTM module branch and a 3D-CNN module branch;

[0057] Use the sliding window method to extract fixed-length input windows and label windows from the fused image features and input them into the ConvLSTM module branch and the 3D-CNN module branch respectively for prediction;

[0058] The ConvLSTM module branch includes three layers of stacked ConvLSTM units. Each layer of ConvLSTM units uses a 3×3 convolutional kernel to extract spatio-temporal features and update the gated states in the ConvLSTM units. The output of the last layer of ConvLSTM units is used as the prediction result of the ConvLSTM module branch.

[0059] The 3D-CNN module branch includes a 3D-CNN module and a fully connected layer; the 3D-CNN module uses a three-dimensional convolutional kernel to extract the joint features of the time and space dimensions, and reduces the dimension of the joint features through three-dimensional pooling operations. The features of the last time step output by the 3D-CNN module are input into the fully connected layer for prediction to obtain the prediction result of the 3D-CNN module branch.

[0060] The prediction results of the ConvLSTM module branch and the prediction results of the 3D-CNN module branch are weighted and averaged to generate the final snow cover distribution prediction result 。

[0061] The present invention also proposes a snow cover distribution prediction system, including the following modules:

[0062] A snow data preprocessing module, which is used to perform orbital correction, radiometric calibration, speckle noise filtering, geocorrection on the synthetic aperture radar (SAR) image, and convert the backscattering coefficient into decibel format to obtain the SAR image data after the first preprocessing, and is used to perform band selection on the multi-spectral optical (opt) image, select bands B2, B3, B4, B8, and B11, where is the blue light band, is the green light band, is the red light band, is the near-infrared band, is the short-wave infrared band. After weighted fusion of the bands B2, B3, B4, and B8, they are registered with the SAR image data after the first preprocessing to obtain the preprocessed synthetic aperture radar (SAR) image data and multi-spectral optical (opt) image data;

[0063] A snow cover distribution information extraction module, which is used to extract the snow cover area image and snow cover area image data ;

[0064] A multi-modal information fusion module, which fuses the snow cover area image and snow cover area image data to obtain the fused image features;

[0065] A snow cover distribution prediction module, which predicts the snow cover distribution based on the fused image features.

[0066] Furthermore, the multi-modal information fusion module includes:

[0067] A shared feature encoder (SFE), which is used to extract the shallow features of the snow cover area image and snow cover area image data respectively to obtain the shallow feature and shallow feature ;

[0068] A basic Transformer encoder (BTE), which is used to extract the low-frequency basic features from the shallow feature , and extract from the shallow feature Low-frequency basic features extracted from ;

[0069] Detail CNN encoder DCE, used to extract high-frequency detail features from shallow features High-frequency detail features extracted from , and high-frequency detail features extracted from shallow features High-frequency detail features extracted from ;

[0070] Frequency domain interaction module, used to transform low-frequency basic features and high-frequency detail features into the frequency domain for fusion, and then back to the time domain to obtain the fused basic features and fused detail features;

[0071] Encoder, which decodes the fused basic features and fused detail features to obtain the finally fused image features .

[0072] Furthermore, the frequency domain interaction module includes a Fourier transform unit, a fusion unit, and an inverse Fourier transform unit;

[0073] The Fourier transform unit is used to transform the low-frequency basic features , low-frequency basic features , high-frequency detail features , high-frequency detail features into the frequency domain respectively;

[0074] Fusion unit, which fuses the low-frequency basic features transformed into the frequency domain to obtain the fused basic features, and fuses the high-frequency detail features transformed into the frequency domain to obtain the fused high-frequency detail features;

[0075] Inverse Fourier transform unit, used to perform inverse Fourier transform on the fused basic features and fused high-frequency detail features in the frequency domain to obtain the fused basic features and fused high-frequency detail features in the time domain.

[0076] Furthermore, the snow cover distribution prediction module includes a ConvLSTM module branch and a 3D-CNN module branch;

[0077] Using the sliding window method, fixed-length input windows and label windows are extracted from the fused image features and input into the ConvLSTM module branch and the 3D-CNN module branch for prediction respectively;

[0078] The ConvLSTM module branch includes three stacked ConvLSTM units. Each layer of ConvLSTM units uses a 3×3 convolutional kernel to extract spatio-temporal features and update the gated states in the ConvLSTM units. The output of the last layer of ConvLSTM units serves as the prediction result of the ConvLSTM module branch.

[0079] The 3D-CNN module branch includes a 3D-CNN module and a fully connected layer. The 3D-CNN module uses a three-dimensional convolutional kernel to extract joint features in the time and space dimensions and reduces the dimension of the joint features through three-dimensional pooling operations. The features of the last time step output by the 3D-CNN module are input into the fully connected layer for prediction to obtain the prediction result of the 3D-CNN module branch.

[0080] The prediction results of the ConvLSTM module branch and the 3D-CNN module branch are weighted and averaged to generate the final snow cover distribution prediction result. 。

[0081] Advantages: Compared with the prior art, the present invention has the following advantages:

[0082] 1. The present invention combines frequency domain interaction and multi-modal fusion to make full use of the complementary information of SAR and multi-spectral optical opt data, improve the quality of feature extraction, and effectively remove the influence of clouds.

[0083] 2. Different from traditional snow depth recognition, the present invention introduces time series information in snow cover prediction, so as to be able to more accurately predict the future distribution of snow cover, which is of great significance in practical applications. Previous prediction tasks often ignored the influence of time information, resulting in prediction results being limited to static analysis and failing to fully consider the dynamic changes in snow cover distribution.

[0084] 3. By analyzing the time series characteristics of snow cover distribution, the present invention identifies its trends and periodicities over time, aiming to improve the prediction accuracy and capture the complex dynamic relationships of snow cover distribution with seasonal changes, climate conditions and other influencing factors. Through this time series analysis, the formation and ablation processes of snow cover can be better understood, providing a scientific basis for the future changes of snow cover. Brief Description of the Drawings

[0085] Figure 1 is the overall idea framework diagram of the present invention;

[0086] Figure 2 is the dataset preprocessing process diagram of the present invention;

[0087] Figure 3 is the multi-modal data fusion flow chart of the present invention;

[0088] Figure 4 It is the flow chart of the dual-branch time series prediction model of the present invention. Specific implementation mode

[0089] The content of the present invention will be further explained below in conjunction with the attached drawings and specific embodiments.

[0090] The present invention provides a snow cover distribution prediction method and system based on a dual-branch model, so as to realize the prediction of the snow cover distribution in Xinjiang region. Specifically, by giving full play to the advantages of the dual-branch model, a rich, accurate and comprehensive information basis is provided for time series prediction. At the same time, the system of the present invention also aims to reduce the interference of weather factors such as clouds and fog on the satellite receiving radiation signal, so as to more accurately identify and predict the snow cover distribution.

[0091] Technical solution: A snow cover distribution prediction system based on a dual-branch model of the present invention includes a snow data preprocessing module, a snow cover distribution information extraction module, a frequency domain interaction module, a multi-modal information fusion module, and a snow cover distribution prediction module connected in sequence:

[0092] The snow data preprocessing module is used to preprocess multi-source remote sensing data to provide a good basis for subsequent prediction; the multi-source remote sensing data includes synthetic aperture radar (SAR) and multi-spectral optics (opt);

[0093] The snow cover distribution information extraction module is used to extract snow cover-related spatio-temporal features and snow cover distribution state features from multi-source remote sensing data;

[0094] The multi-modal fusion module is used to effectively integrate the data of SAR (synthetic aperture radar image) and opt (multi-spectral optics) images obtained from Sentinel-1 and Sentinel-2 satellites to obtain fusion features. The frequency domain interaction module is included in the multi-modal fusion module. This module preprocesses and interacts the spatio-temporal features related to snow cover and the snow cover distribution state features at the frequency domain level and outputs them to the multi-modal fusion module, that is, in the multi-modal fusion process, the fusion features are first frequency domain interacted, the fusion features are converted between time domain and frequency domain, and then multi-modal fusion outputs the fusion features to give full play to the advantages of each modal data;

[0095] The snow cover distribution prediction module simulates and predicts the snow cover distribution based on the fusion features obtained by the multi-modal fusion module to obtain the final snow cover distribution prediction;

[0096] As Figure 1 shown, the operation steps of the above snow data preprocessing module, snow cover distribution information extraction module, frequency domain interaction module, multi-modal information fusion module, and snow cover distribution prediction module are as follows:

[0097] (1)As shown Figure 2 below, the operation steps of the snow cover data preprocessing module are as follows:

[0098] Obtain SAR (Synthetic Aperture Radar) images and opt (multi-spectral optical) images from Sentinel-1 and Sentinel-2 satellites respectively. These two types of images have the advantages of high resolution, high revisit rate, and wide coverage. Moreover, the SAR (Synthetic Aperture Radar) image and opt (multi-spectral optical) image data can also be obtained for free, which enables the present invention to be widely applied in many different regions. All SAR (Synthetic Aperture Radar) image and opt (multi-spectral optical) image data are resampled using an alignment grid with 10m GSD.

[0099] (1.1)Synthetic Aperture Radar SAR data preprocessing steps:

[0100] (1.1.1)Unzip the downloaded remote sensing data SAR GRD file of the Synthetic Aperture Radar SAR, perform orbit correction Apply Orbit File, change the satellite state data in the xml file, and correcting this data will make the positioning accuracy higher, obtaining the image after orbit correction. Its core is to use an accurate orbit file to correct the orbit parameters in the original data, and its core can be simplified into the following form:

[0101]

[0102] Where represents the corrected orbit parameters, represents the original orbit parameters, is the orbit parameter, describing the position and motion state of the satellite in the orbit, represents the corrected one, that is, the orbit parameters after being adjusted by an accurate orbit file, represents the original one, that is, the uncorrected orbit parameters; represents the correction amount calculated through an accurate orbit file, represents the correction amount, that is, the difference between the orbit parameters in the accurate orbit file and the original orbit parameters.

[0103] (1.1.2)Perform radiometric calibration Calibrate on the image after orbit correction, convert the backscatter coefficient into a physical quantity, and eliminate the systematic error of the sensor itself, obtaining the image after error elimination. Its formula is:

[0104]

[0105] Where is the normalized backscatter coefficient, represents the backscatter coefficient, Indicates normalization. The normalized backscattering coefficient is used to eliminate the systematic error of the sensor itself; is the original backscattering signal intensity, represents power and is the original signal intensity received from the Synthetic Aperture Radar (SAR) sensor, is the unprocessed raw data, the signal intensity directly obtained from the sensor, which contains the systematic error of the sensor itself; is the calibration factor, usually determined by the gain and offset of the sensor.

[0106] (1.1.3) Perform Single Product Speckle Filtering on the image after error elimination to reduce and eliminate the speckle noise in the image and obtain the filtered image; the commonly used method is RefinedLee filtering. Its formula can be expressed as:

[0107]

[0108] where is the pixel value after filtering, is the image, is the one after filtering processing, is the coordinate position in the image; represents the average of the pixel values in the window ; represents the sum of all pixel values in the window ; represents the offset in the window ; is the horizontal offset, is the vertical offset; is the pixel value within the neighborhood centered on the pixel ; is the filtering window, is the number of pixels within the window.

[0109] (1.1.4)Perform geographic correction (Range-Doppler Terrain Correction) on the filtered image. Project the image from slant range to ground range through the terrain model and orbit information to obtain the geographically corrected image, which reflects the true geometric shape of the ground surface. Its core formula is:

[0110]

[0111] where is the image coordinate after geographic correction, represents the coordinate, Represents a geographically corrected coordinate image, which reflects the true geometry of the earth's surface. Is the slant range coordinate. Represents the slant range. Is the digital elevation model. Is the orbital parameter. Is the Range-Doppler terrain correction function. Represents the function. Represents Range-Doppler terrain correction, and this function is used to convert the slant range coordinate to the geographic coordinate.

[0112] (1.1.5)Convert the backscatter coefficient of the geographically corrected image to decibel format through Data Conversion to enhance the image contrast. The formula is:

[0113]

[0114] Where Is the backscatter coefficient in decibels. Is the normalized backscatter coefficient. Represents the backscatter coefficient. Represents normalization.

[0115] Finally, a single-channel synthetic aperture radar (SAR) image with a resolution of is obtained. That is, the shape of the synthetic aperture radar (SAR) image is .

[0116] (1.2)Optical image opt data preprocessing steps:

[0117] (1.2.1)For Sentinel-2 Level 2A data, perform data selection and band extraction (Band Select), and extract the 10m resolution bands: blue band B2, green band B3, red band B4, and near-infrared band B8. The formula is as follows:

[0118]

[0119] (1.2.2) For data, perform radiometric calibration to convert the digital value to reflectance. The formula is as follows:

[0120]

[0121] Where Is the surface reflectance. Is the radiance after atmospheric correction. Is the calibration coefficient.

[0122] (1.2.3) To make full use of high-resolution multispectral information, for the in bands are fused. The weighted average method is used for band fusion Layer Stacking, and its formula is as follows:

[0123]

[0124] where is the image after band fusion, is the weight of the th band, are the selected bands B2, B3, B4, and B8.

[0125] (1.3) Image registration

[0126] Image registration is to align images taken at different times, by different sensors, or with different resolutions into the same coordinate system. ENVI provides an automated registration tool (Image Registration Workflow), which can efficiently complete the image registration task of synthetic aperture radar SAR and multispectral optical opt. The image registration process of ENVI can be described by the following formula:

[0127] (1.3.1) Tie point generation: Based on the feature matching algorithm, Tie points are automatically generated. Tie points are the basis of image registration and are used to establish corresponding relationships between different images. The purpose of generating Tie points is to use these Tie points as references to ensure that the two images can be aligned spatially. Tie points can be expressed as:

[0128]

[0129] where, represents the coordinates of the th matching point on the reference image,

[0130] represents the coordinates of the matching point on the image to be registered. In this specific embodiment, the image to be registered is the synthetic aperture radar image SAR and the multispectral optical opt image in the same area at the same time after preprocessing.

[0131]

[0132] where, are the coordinates of the image to be registered. In this specific embodiment, the image to be registered is the synthetic aperture radar (SAR) image and the multi-spectral optical (opt) image of the same area at the same time after preprocessing. are the transformed coordinates. ~ 、 ~ are the polynomial coefficients.

[0133] (1.3.3) Resampling method: Use the cubic convolution interpolation method to resample the image after image registration to generate the final registered image.

[0134] Finally, a 4-channel multi-spectral optical (opt) image with a resolution of is obtained. That is, the shape of the multi-spectral optical (opt) image is .

[0135] (2) The operation steps of the snow cover distribution information extraction module are as follows:

[0136] Extract the snow cover distribution information from the synthetic aperture radar (SAR) and multi-spectral optical (opt) images obtained after being processed by the snow cover data preprocessing module.

[0137] (2.1) Snow cover information extraction based on the preprocessed synthetic aperture radar (SAR) data;

[0138] (2.1.1) Calculate the backscattering coefficient : Extract the backscattering coefficient from the synthetic aperture radar (SAR) image to distinguish wet snow and dry snow.

[0139] (2.1.2) Set the wet snow threshold: Judge the wet snow area according to the range of the backscattering coefficient:

[0140]

[0141] (2.2) Snow cover information extraction based on the preprocessed multi-spectral optical (opt) data;

[0142] (2.2.1) Calculate the Normalized Difference Snow Index (NDSI): Use the green band and the short-wave infrared band of Sentinel-2 to calculate the Normalized Difference Snow Index NDSI:

[0143]

[0144] where, and are the green band and the short-wave infrared band Reflectivity.

[0145] (2.2.2)Set a threshold to extract snow cover: Extract the snow cover area according to the NDSI value, and set the threshold to 0.4

[0146]

[0147] Generally, the spectral characteristics of snow cover are as follows: in the near-infrared band The reflectivity is low because snow cover absorbs more in the near-infrared region, and in the green band The reflectivity is high because snow cover reflects strongly in the visible light region. Therefore, after screening out the possible snow cover areas, use the near-infrared band to further exclude cloud interference, and use the green band to further exclude water body interference. By combining the near-infrared band and the green band reflectivity, the accuracy of snow cover identification is further improved.

[0148] (3)As Figure 3 shown, the operation steps of the multi-modal fusion module are as follows:

[0149] (3.1)Encoder part encoder

[0150] This encoder has a total of three components: a shared feature encoder (SFE) based on the Restormer block, a basic Transformer encoder (BTE) based on the LiteTransformer block, and a detailed CNN encoder (DCE) based on the invertible neural network (INN) block.

[0151] (3.1.1)Extract shallow features:

[0152] The Restormer block-based share feature encoder (SFE) extracts the common features in the two-modal images obtained through the snow cover distribution information processing module. The SFE module consists of multiple Restormer blocks and can extract global features from high-resolution images by applying self-attention in the feature dimension. The global features of synthetic aperture radar SAR include backscattering intensity, polarization features, and texture features, and the global features of multi-spectral optics opt include the normalized difference snow index (NDSI), color, and spectral features. Combining the two can more accurately identify and monitor snow cover.

[0153] Due to the different physical characteristics of the two modalities, directly extracting shared features may lead to information loss and confusion. Therefore, it is necessary to extract features from SAR and opt images separately to retain the unique information of each modality, that is, directly extracting shared features will result in the loss of modality-specific information. Therefore, it is necessary to extract and process them separately before extracting shared features.

[0154] The formula for extracting their respective shallow features is:

[0155]

[0156] Where is the encoder SFE, and are the synthetic aperture radar SAR and multispectral optical opt images obtained through the snow cover distribution information processing module, means inputting the synthetic aperture radar SAR into the SFE encoder to extract the features of the synthetic aperture radar SAR, means inputting the multispectral optical opt image into the same SFE encoder to extract the features of the multispectral optical opt image, is the shallow feature extracted from the synthetic aperture radar SAR image in, is the multispectral image opt image in the shallow feature extracted from.

[0157] (3.1.2)Low-frequency global feature extraction:

[0158] Base Transformer Encoder (BTE) is used to extract the low-frequency basic features of synthetic aperture radar SAR and multispectral optical opt from the shallow features , that is: , namely:

[0159]

[0160] Where represents the BTE encoder represents the low-frequency basic feature.

[0161] (3.1.3)High-frequency detail feature extraction:

[0162] Detail CNN Encoder (DCE) extracts high-frequency detail features from the shallow features , that is: , namely:

[0163]

[0164] Where, Denote the encoder DCE, represent the high-frequency detail features.

[0165] (3.2)Feature fusion

[0166] Perform a frequency-domain transformation on the low-frequency basic features and and the high-frequency detail features and , perform a frequency-domain transformation, perform feature fusion in the frequency domain, form a feature map containing rich information, and finally convert it back from the frequency domain to the time domain. The operation steps of the frequency-domain interaction module are as follows:

[0167] (3.2.1)Perform a two-dimensional Fourier transform (FFT) on the basic features and detail features of the synthetic aperture radar (SAR) image respectively, and perform a two-dimensional Fourier transform (FFT) on the basic features and detail features of the multi-spectral optical (opt) image, and convert them from the time domain to the frequency domain.

[0168]

[0169]

[0170] Among them, denote the two-dimensional Fourier transform operation, and are the basic frequency-domain features of the optical and synthetic aperture radar (SAR) images respectively, and are the detail frequency-domain features of the optical and synthetic aperture radar (SAR) images respectively.

[0171] (3.2.2)In the frequency domain, fuse the basic features of the multi-spectral optical (opt) image and the synthetic aperture radar (SAR) image, and fuse the detail features of the multi-spectral optical (opt) image and the synthetic aperture radar (SAR) image. Use a fusion layer with a channel attention mechanism to perform weighted sum and interaction on each frequency-domain feature to capture richer frequency-domain information.

[0172]

[0173]

[0174] Among them, denote the fusion operation with a channel attention mechanism, which can perform weighting according to the importance of each channel in the frequency domain, denote the basic features fused in the frequency domain, denote the detail features fused in the frequency domain.

[0175] (3.2.3) After the fusion operation, the fused features in the frequency domain are reconverted back to the time domain using the inverse Fourier transform (IFFT) to obtain the final fused features.

[0176]

[0177] Among them, represents the two-dimensional inverse Fourier transform operation, represents the fused base features in the time domain, represents the fused detail features in the time domain.

[0178] (3.3) Decoding operation

[0179] (3.3.1) Fuse the base feature and the detail feature to obtain the fused feature , which is expressed by the formula:

[0180]

[0181] Where is the fused feature, containing low-frequency and high-frequency information, represents the base feature and the detail feature; represents the concatenation operation, which concatenates the two feature maps in the channel dimension; represents the two-dimensional convolution operation, which is used to perform convolution on the feature map, and the convolution kernel size is set to 1. (3.3.2) Reduce the number of channels of the feature map through the convolution layer to reduce the computational complexity while maintaining important feature information, and input the feature map after channel reduction into the encoder layer to further extract and strengthen the features. The formula is as follows:

[0182]

[0183] Among them, is the finally obtained fused feature map; represents the two-dimensional convolution operation, and the convolution kernel size is set to 3; represents the activation function, which is used to map the values of the feature map to the range of (0, 1). The shape of the finally obtained feature map is [ h eig h t,widt h ,c h annels] .

[0184] (4)The snow accumulation distribution prediction module includes three stacked ConvLSTM layers and a 3D-CNN module;

[0185] The features obtained from the previous section As the input of the snow accumulation distribution prediction module, the running steps of the snow accumulation distribution prediction module are as follows:

[0186] As Figure 4 shown, the feature sequence after concatenating the features obtained from the previous section is used as the input of the snow accumulation prediction module, with the shape of [batch size, time steps, channels, height, width]. The features at each time step incorporate the optical image and multi-channel information, as well as the single-channel image of the synthetic aperture radar image.

[0187] Use the sliding window method to extract fixed-length input windows and label windows from the feature sequence. The input window contains the features of several time steps and serves as the input for ConvLSTM and 3D-CNN; the label window contains the target snow accumulation distribution for the next time step and is used for supervised learning. The formulas for selecting the input and label using the sliding window method are expressed as:

[0188] I n Input Sequence: X input [t:t+time_steps]

[0189] Label Sequence: Y label [t:t+tim e steps :t+tim e steps +label_steps]

[0190] Among them, X input [t:t+time_steps] represents the input sequence, starting from time step and being an input sequence segment with a length of ; Y label [t:t+tim e steps :t+tim e steps +label_steps] represents the label sequence, starting from time step and being a label sequence segment with a length of ; is the starting position of the current time step and represents the starting point of sequence segmentation, is the length of the input window, is the length of the label window.

[0191] (4.1) Running steps of the ConvLSTM module

[0192] Step 4.1.1, at each time t, the ConvLSTM receives the input sequence vector , the hidden state at the previous moment and the cell state , and calculates the values of the forget gate, input gate, and output gate through convolutional operations. The specific calculation formulas are as follows:

[0193] f t = σ ( ω f • X input [t] , h t -1 + b f )

[0194] i t =σ( ω i • X input [t], h t-1 + b i )

[0195] o t =σ( ω o • X input [t], h t-1 + b o )

[0196] where , and represent the forget gate, input gate, and output gate respectively; represents the sigmoid function; , , and are the weight matrices of the LSTM model respectively, is the hidden state at time t - 1, X input [t] is the feature of the input time series at time step t, , and represent the bias vectors.

[0197] Step 4.1.2, State Update and Hidden State Calculation: Update the cell state and the hidden state according to the output of the gating mechanism. The current candidate cell state The formula for

[0198] c t ̃ = tanh ( ω c • X input [t], h t-1 + b c )

[0199] where is the candidate cell state, represents the time step; represents the weight matrix, represents the candidate cell state; represents the hyperbolic tangent function, which maps the input to the range (-1, 1) and is used to represent the value of the cell state; represents the concatenation operation, which concatenates vectors together to form a larger vector, represents the hidden state at the previous time step; represents the bias vector, represents the candidate cell state.

[0200] The input gate and the forget gate respectively determine the and information proportion occupied in the current cell state , and the current cell state

[0201]

[0202] where represents the cell state, represents the time step cell state; represents the forget gate, represents the time step; represents the previous time point, represents the candidate cell state, which is the candidate memory cell at the current time step and is used to update the cell state; represents element-wise multiplication, which multiplies the corresponding elements of two matrices.

[0203] The output formula of the hidden layer is:

[0204]

[0205] where Denote the output of the hidden layer at time step t; Denote the activation value of the output gate; Represent the state of the current cell. Represent the multiplication of two matrix elements, Denote the cell state Apply the hyperbolic tangent function to scale the value of the cell state to the range (-1, 1), making the output more stable.

[0206] Adopt three stacked ConvLSTM layers, use multiple ConvLSTM Cells for hierarchical calculation, the input of each layer is the output of the previous layer, gradually capture spatial and temporal features at different levels, the number of hidden units in each layer is 256, use The convolutional kernel of is used to extract spatial features and update the gated state in the LSTM cell, supporting any value of the input time series length. The output of each layer includes the hidden state h and the updated memory cell c at each time step. Select the hidden state Output at the last time step of the last layer as the representation of temporal features to capture the temporal information at the current moment. Specifically, the final hidden state Can be expressed as:

[0207]

[0208] Among them, T represents the last time step of the time series, and L represents the last layer of ConvLSTM.

[0209] Step 4.1.3, obtain the prediction result of this branch of ConvLSTM , by Input into the fully connected layer for processing. The weight matrix of the fully connected layer is , and the bias vector is . The calculation formula of the prediction result is:

[0210]

[0211] Among them, Is the final prediction result obtained from this branch of ConvLSTM.

[0212] (4.2) Operating steps of the 3D-CNN module

[0213] Step 4.2.1, three-dimensional convolution operation: Define the three-dimensional convolutional kernel , where Is the number of channels of the input data, Is the size of the convolutional kernel in the time dimension, And Are the sizes of the convolutional kernel in the spatial dimension; The input tensor is , where is the number of time steps, is the number of channels, and is the spatial resolution. In 3D convolution, the th element of the output feature map can be calculated by the following formula:

[0214]

[0215] where represents the output of the th convolutional kernel at the position at time step ; is the weight of the th convolutional kernel, is the local region of the input data, is the bias term, is the activation function.

[0216] Step 4.2.2, 3D pooling operation: Define the 3D pooling window size as , then the formula for the pooling operation is:

[0217] Y p ool t,k,i,j =Pooling( Y conv [t :t+ D p ,i :i+ H p ,j :j+ W p ])

[0218] where, represents the pooling result of the th channel at the position at time step ; represents downsampling the local region of the input data through the pooling operation to reduce the spatial dimension of the feature map while retaining important information; Y conv [t :t+ D p ,i :i+ H p ,j :j+ W p ] represents the local region of the input data for the pooling operation, represents the current time step being processed, represents the channel index, i.e., the current channel number being processed, and respectively represent the positions of the currently processed feature map in the height and width directions; represents the coverage range of the pooling operation in the time dimension, represents the coverage range of the pooling operation in the height direction, represents the coverage range of the pooling operation in the width direction, represents the pooling operation.

[0219] Step 4.2.3, extract the features of the last time step from the output of the 3D-CNN module as the spatio-temporal feature representation at the current moment. The output feature map is , and the extracted features are:

[0220] F final = Y p ool [:,:,T,:,:]

[0221] Among them, is the final extracted feature, with a shape of [batch size, output channels, height, width], and T is the last time step of the number of time steps.

[0222] Step 4.2.4, input the extracted into the fully connected layer to generate the final prediction result. The weight matrix of the fully connected layer is , and the bias vector is . The calculation formula for the prediction result is:

[0223]

[0224] Among them, is the prediction result of the 3D-CNN branch, represents flattening the feature map into a one-dimensional vector.

[0225] (4.3) Snow distribution prediction

[0226] Integrate the prediction results of the ConvLSTM module and the 3D-CNN module to generate the final prediction result . The integration method is weighted average, and the formula is expressed as:

[0227]

[0228] Among them is a weight coefficient used to balance the contributions of the two branches, is the final prediction result.

Claims

1. A snow distribution prediction method, characterized in that: The following steps are involved: Step 1, acquiring remote sensing data of a synthetic aperture radar SAR image and a multispectral optical opt image, preprocessing the remote sensing data of the synthetic aperture radar SAR image and the multispectral optical opt image to obtain preprocessed synthetic aperture radar SAR image data and multispectral optical opt image data; Step 2: Extract snow area image data based on the preprocessed synthetic aperture radar SAR image data ; Extracting snow area image data based on preprocessed multispectral optical opt image data ; Step 3: extracting the snow area image data Image data of the snow area Perform feature fusion to obtain fused image features; Step 4: predicting snow distribution based on the fused image features.

2. A snow distribution prediction method according to claim 1, characterized in that: The step 1 of preprocessing the remote sensing data of the synthetic aperture radar SAR image and the multispectral optical opt image specifically includes: Step 1.1, performing orbit correction, radiation calibration, speckle noise filtering, geographic correction, and converting the backscatter coefficient into decibel format on the SAR image data in sequence to obtain the SAR image data after the first preprocessing; Step 1.2, processing the opt image data, specifically: First, perform band selection on the opt image data to obtain the selected band , expressed as: ； in, For the blue light band, For the green light band, For the red light band, The near-infrared band, It is the short-wave infrared band; For the selected band Perform radiation calibration and then select and Perform weighted average fusion to obtain the opt image data after the first preprocessing, expressed as: ； in, is the image after band fusion, It is The weight of the band, It is the bands B2, B3, B4 and B8 after radiometric calibration; Step 1.3, perform image registration on the SAR image data after the first preprocessing in step 1.1 and the opt image data after the first preprocessing in step 1.2, and resample the image after image registration using the cubic convolution method to obtain the synthetic aperture radar SAR image preprocessed in step 1 and the preprocessed multispectral optical opt image.

3. A snow distribution prediction method according to claim 2, characterized in that: Extracting snow area image data based on the multispectral optical opt image data in step 2 , specifically: Use green light band and shortwave infrared bands Calculate the Normalized Difference Snow Index NDSI: ； in, and Green band and shortwave infrared bands Reflectivity; Set the threshold to extract snow T1, and extract the snow area according to the NDSI value , get the snow area image data : 。 4. A snow distribution prediction method according to claim 1, characterized in that: The step 3 comprises the following steps: Step 3.1: Extract the snow area image data Image data of the snow area The shallow features of are expressed as: ； in, For encoder SFE, From synthetic aperture radar SAR images The shallow features extracted from For multispectral optical opt images The shallow features extracted from Step 3.2, extracting low-frequency basic features from the shallow features, expressed as: ； in Represents the BTE encoder, Represents shallow features The low-frequency basic features extracted from Represents shallow features The low-frequency basic features extracted from Step 3.3, extracting high-frequency detail features from the shallow features, expressed as: ； in, Indicates the encoder DCE, Represents that from shallow features The high-frequency detail features extracted from Represents that from shallow features The high-frequency detail features extracted from Step 3.4, fusing the low-frequency basic features with the high-frequency detail features to obtain fused image features.

5. A snow distribution prediction method according to claim 4, characterized in that: The step 3.4 is as follows: Step 3.4.1: The snow area image data The low-frequency basic characteristics and high frequency detail features Convert the time domain to the frequency domain and convert the snow area image data The low-frequency basic characteristics and high frequency detail features Converting from the time domain to the frequency domain, it is represented as follows: ； ； in, represents a two-dimensional Fourier transform operation, and They are the basic frequency domain features of optical and synthetic aperture radar SAR images, and They are the detailed frequency domain features of optical and synthetic aperture radar SAR images; Step 3.4.2: In the frequency domain, the basic features of the multispectral optical opt image are transformed into and basic features of synthetic aperture radar SAR images Fusion and detailed features of multi-spectral optical opt images and SAR image details The fusion is expressed as: ； ； in, represents the fusion operation with channel attention mechanism, represents the basic features fused in the frequency domain, Represents the detailed features fused in the frequency domain; Step 3.4.3, fusion features in the frequency domain and Reconvert back to the time domain to obtain the fused basic features and fused detail features in the time domain; ； in, represents the two-dimensional inverse Fourier transform operation, represents the basic features after fusion in the time domain, Represents the detailed features after fusion in the time domain; Step 3.4.4, the fused basic features in the time domain And the fused detail features Input into the decoder to generate the final fused image features , expressed as: ； in, is the fused feature after decoding, Represents a decoder.

6. A snow distribution prediction method according to claim 1, characterized in that: In step 4, the snow distribution is predicted based on the fused image features; Specifically, a dual-branch model is used to predict snow distribution; the dual-branch model includes a ConvLSTM module branch and a 3D-CNN module branch; Using a sliding window method, extracting an input window and a label window of fixed length from the fused image features, and inputting them into a ConvLSTM module branch and a 3D-CNN module branch for prediction respectively; The ConvLSTM module branch includes three layers of stacked ConvLSTM units, each layer of ConvLSTM units uses a 3×3 convolution kernel to extract spatiotemporal features and update the gating state in the ConvLSTM unit, and the output of the last layer of ConvLSTM units is used as the prediction result of the ConvLSTM module branch; The 3D-CNN module branch includes a 3D-CNN module and a fully connected layer; the 3D-CNN module uses a three-dimensional convolution kernel to extract joint features of time and space dimensions, and reduces the dimension of the joint features through a three-dimensional pooling operation, and the features of the last time step output by the 3D-CNN module are input to the fully connected layer for prediction, thereby obtaining a prediction result of the 3D-CNN module branch; The prediction results of the ConvLSTM module branch and the 3D-CNN module branch are weighted averaged to generate the final snow distribution prediction result. .

7. A snow distribution prediction system, characterized in that: Includes the following modules: Snow data preprocessing module is used to perform orbit correction, radiation calibration, speckle noise filtering, geographic correction, and convert the backscatter coefficient into decibel format on the synthetic aperture radar SAR image to obtain the SAR image data after the first preprocessing, and to perform band selection on the multispectral optical opt image to select the band , , , and ,in For the blue light band, For the green light band, For the red light band, The near-infrared band, The short-wave infrared band is a short-wave infrared band, and the bands B2, B3, B4 and B8 are weightedly fused and then registered with the SAR image data after the first preprocessing to obtain preprocessed synthetic aperture radar SAR image data and multispectral optical opt image data; The snow distribution information extraction module is used to extract the snow area image from the pre-processed synthetic aperture radar SAR image data and the multi-spectral optical opt image data. and snow area image data ; The multimodal information fusion module combines the snow area image and snow area image data Fusion, to obtain fused image features; a snow distribution prediction module, to predict snow distribution based on the fused image features.

8. A snow distribution prediction system according to claim 7, characterized in that: The multimodal information fusion module includes: A shared feature encoder SFE is used to extract the snow area image respectively and snow area image data The shallow features of and shallow features ; Basic Transformer encoder BTE, used to extract features from shallow layers Extract low-frequency basic features , and from the shallow features The low-frequency basic features extracted from ; Detailed CNN encoder DCE, used to extract features from shallow layers The high-frequency detail features extracted from , and from the shallow features The high-frequency detail features extracted from ; The frequency domain interaction module is used to transform the low-frequency basic features and high-frequency detail features into the frequency domain for fusion, and then transform them back into the time domain to obtain the fused basic features and fused detail features; The encoder decodes the fused basic features and the fused detail features to obtain the final fused image features. .

9. A snow distribution prediction system according to claim 8, characterized in that: The frequency domain interaction module includes a Fourier transform unit, a fusion unit and an inverse Fourier transform unit; The Fourier transform unit is used to transform the low-frequency basic features , low frequency basic characteristics , high frequency detail features , high frequency detail features Transform to frequency domain respectively; A fusion unit, which fuses the low-frequency basic features transformed into the frequency domain to obtain fused basic features, and fuses the high-frequency detail features transformed into the frequency domain to obtain fused high-frequency detail features; The inverse Fourier transform unit is used to perform inverse Fourier transform on the fused basic features and the fused high-frequency detail features in the frequency domain to obtain the fused basic features and the fused high-frequency detail features in the time domain.

10. A snow distribution prediction system according to claim 8, characterized in that: The snow distribution prediction module includes a ConvLSTM module branch and a 3D-CNN module branch; The sliding window method is used to extract the fused image features The fixed-length input window and label window are extracted and input into the ConvLSTM module branch and the 3D-CNN module branch for prediction respectively; The ConvLSTM module branch includes three layers of stacked ConvLSTM units, each layer of ConvLSTM units uses a 3×3 convolution kernel to extract spatiotemporal features and update the gating state in the ConvLSTM unit, and the output of the last layer of ConvLSTM units is used as the prediction result of the ConvLSTM module branch; The 3D-CNN module branch includes a 3D-CNN module and a fully connected layer; the 3D-CNN module uses a three-dimensional convolution kernel to extract joint features of time and space dimensions, and reduces the dimension of the joint features through a three-dimensional pooling operation, and the features of the last time step output by the 3D-CNN module are input to the fully connected layer for prediction, thereby obtaining a prediction result of the 3D-CNN module branch; The prediction results of the ConvLSTM module branch and the 3D-CNN module branch are weighted averaged to generate the final snow distribution prediction result. .

Images