A snow accumulation distribution prediction method and system based on a dual-branch model

Abstract

The present invention discloses a snow cover distribution prediction method and system based on a dual-branch model, belonging to the technical field of remote sensing image processing. The system includes a snow cover 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. The snow cover data preprocessing module is used to process remote sensing data. The snow cover distribution information extraction module is used to extract feature information related to snow cover from multi-source remote sensing data. The frequency domain interaction module is used to perform frequency domain analysis on the fused remote sensing data. The multi-modal information fusion module is used to integrate multi-source remote sensing data from different types of sensors and perform deep fusion on them. The snow cover distribution prediction module is used to predict the snow cover distribution 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 indispensable part of the cryosphere, is not only a core variable for studying climate dynamics and hydrological cycles, 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. At the same time, snow cover plays a pivotal role in the water cycle, and snowmelt runoff provides a stable and reliable water source for many water-scarce areas. However, the distribution of snow cover is not constant. It is affected by a variety of climatic and geographical factors and shows complex and changeable 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 heterogeneity in spatial 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 defect of insufficient number of stations in traditional snow cover field measurement methods and can provide long-term and large-scale snow cover 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 cover 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, can detect the information of the underlying surface under the snow cover, and have the characteristics of all-weather and large-scale snow cover 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. Then, 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, 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 both 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 to 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 orbit correction, radiometric calibration, speckle noise filtering, geocorrection, 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 after radiometric calibration , , and ;

[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 perform resampling on the registered image using the cubic convolution interpolation method to obtain the synthetic aperture radar (SAR) image after preprocessing in Step 1 and the multi-spectral optical (opt) image after preprocessing.

[0023] Further, the snow-covered area image data is extracted based on the multi-spectral optical (opt) image data 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 image data of the snow cover area and the image data of the snow cover area 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 features extracted from the shallow features , represents the low-frequency basic features extracted from the shallow features ;

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

[0037]

[0038] Among them, represents the encoder DCE, represents the high-frequency detail features extracted from the shallow features , 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 specifically 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 Detailed Features of Synthetic Aperture Radar (SAR) Images 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 to obtain the fused base features and fused detail features in the time domain.

[0050]

[0051] 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;

[0052] Step 3.4.4, input the fused base features and the fused detail 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, use a dual-branch model 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 a fixed-length input window and a label window 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 stacked ConvLSTM units. Each layer of ConvLSTM unit uses a 3×3 convolutional kernel to extract spatio-temporal features and update the gated state in the ConvLSTM unit. The output of the last layer of ConvLSTM unit 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 multispectral 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 multispectral 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 multimodal 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 multimodal 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 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 to 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 ; And high-frequency detail features extracted from shallow features ; ;

[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 low-frequency basic features , low-frequency basic features , high-frequency detail features , high-frequency detail features Into the frequency domain respectively;

[0074] The fusion unit 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] The inverse Fourier transform unit is 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 is used 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 three-dimensional convolutional kernels 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 to 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] Beneficial effects: 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 identification, the present invention introduces time series information in snow cover prediction, so as 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 of 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 process of snow formation and ablation can be better understood, providing a scientific basis for the future changes of snow cover. 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. Detailed implementation manners

[0089] The content of the present invention will be further explained below in conjunction with the accompanying 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 received 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 cover 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 cover 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 the spatio-temporal features related to snow cover and the 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 subjected to frequency domain interaction, and 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 running steps of the above snow cover 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 in Figure 2 , the running steps of the snow cover data preprocessing module are as follows:

[0098] Obtain SAR (Synthetic Aperture Radar) images and opt (multispectral 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 (multispectral 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) images and opt (multispectral optical) image data are resampled using an alignment grid of 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 status data in the xml file, and correcting this data will result in higher positioning accuracy to obtain 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, which describes 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 the 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 to obtain 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 that has undergone 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 Geometric 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 geometrically corrected image, which reflects the true geometric shape of the ground surface. Its core formula is:

[0110]

[0111] where is the image coordinate after geometric correction, represents the coordinate, Represents a geocorrected coordinate image that 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 a function. Represents Range-Doppler terrain correction, which is used to convert the slant range coordinate to a geographic coordinate.

[0112] (1.1.5) Convert the backscatter coefficient of the geocorrected 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) to 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) Perform radiometric calibration on the data 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, the bands in are subjected to fusion processing. 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) that 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, and

[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 SAR image and the multispectral optical (opt) image of 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 a synthetic aperture radar (SAR) image and a multispectral 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 multispectral optical (opt) image with a resolution of is obtained. That is, the shape of the multispectral 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 multispectral 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 multispectral optical (opt) data;

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

[0143]

[0144] where and are the green band and the shortwave infrared band respectively.Reflectance.

[0145] (2.2.2) Set the 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 reflectance is low because snow cover absorbs more in the near-infrared region, and in the green light band The reflectance 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 light band to further exclude water body interference. By combining the near-infrared band and the green light band reflectance, 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, they need to be extracted separately and then the shared features are extracted.

[0154] The formulas for extracting their respective shallow features are as follows:

[0155]

[0156] Where is the encoder SFE, and are the synthetic aperture radar (SAR) and multispectral optical (opt) images obtained after being processed by the snow 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 ; is the shallow feature extracted from the multispectral image (opt) image ;

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

[0158] The Base Transformer Encoder (BTE) is used to extract the low-frequency basic features of the 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 features.

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

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

[0163]

[0164] Where Denote the encoder DCE, represent 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 to 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 two-dimensional Fourier transforms (FFT) on the basic features and detail features of the synthetic aperture radar (SAR) image respectively, and perform two-dimensional Fourier transforms (FFT) on the basic features and detail features of the multi-spectral optical (opt) image to 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 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 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] Among them 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 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 (0, 1), and 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 content of 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 splicing the features obtained from the content of the previous section is used as the input of the snow accumulation prediction module, with a 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 convolution 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 is:

[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 moment ; 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 time the cell state at the moment; 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 is used to multiply the corresponding elements of two matrices.

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

[0204]

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

[0206] Adopts three stacked ConvLSTM layers, uses multiple ConvLSTM Cells for hierarchical calculation, the input of each layer is the output of the previous layer, gradually capturing spatial and temporal features at different levels, the number of hidden units in each layer is 256, and uses 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 of the last time step of the last layer Output as the representation of the time feature to capture the temporal information at the current moment. Specifically, the final hidden state

[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 the ConvLSTM branch , by Input to the fully connected layer for processing. The weight matrix of the fully connected layer is , and the bias vector is . The calculation formula for the prediction result is:

[0210]

[0211] Among them, is the final prediction result obtained by the ConvLSTM branch.

[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, reducing 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 accumulation 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 method for predicting snow accumulation distribution, characterized in that, It includes the following steps: Step 1: Obtain remote sensing data of a synthetic aperture radar (SAR) image and a multispectral optical (opt) image, and preprocess the remote sensing data of the SAR image and the multispectral optical (opt) image to obtain preprocessed SAR image data and multispectral optical (opt) image data; Step 2: Extract snow-covered area image data based on the preprocessed synthetic aperture radar (SAR) image data ; Extract snow-covered area image data based on preprocessed multi-spectral optical opt image data ; Step 3, perform feature fusion on the extracted snow-covered area image data and the snow-covered area image data to obtain the fused image features; It includes the following steps: Step 3.1, extract the image data of the snow-covered area respectively and the image data of the snow-covered area of the shallow features, expressed as: ； Among them, is the encoder SFE, is the shallow feature extracted from the synthetic aperture radar (SAR) image and is the shallow feature extracted from the multi-spectral optical (opt) image ;​ Step 3.2: Extract low-frequency basic features from the shallow features, expressed as: ； Among them represents the BTE encoder, represents the low-frequency basic features extracted from the shallow features and represents the low-frequency basic features extracted from the shallow features as well;​ Step 3.3: Extract high-frequency detail features from the shallow features, expressed as: ； Among them, represents the encoder DCE, represents the high-frequency detail features extracted from the shallow features and represents the high-frequency detail features extracted from the shallow features as well; Step 3.4: Fuse the low-frequency basic features and the high-frequency detail features to obtain fused image features, specifically as follows: Step 3.4.1, convert the low-frequency basic features and high-frequency detail features of the snow-covered area image data from the time domain to the frequency domain, and convert the low-frequency basic features and high-frequency detail features of the snow-covered area image data from the time domain to the frequency domain, which is expressed as follows: ； ； Among them, represents a two-dimensional Fourier transform operation, and are the basic frequency domain features of optical and synthetic aperture radar (SAR) images respectively, and are the detailed frequency domain features of optical and synthetic aperture radar (SAR) images respectively; Step 3.4.2, in the frequency domain, fuse the basic features of the multi-spectral optical opt image and the basic features of the synthetic aperture radar SAR image and fuse the detailed features of the multi-spectral optical opt image and the detailed features of the synthetic aperture radar SAR image as follows: ； ； Among them, represents the fusion operation with channel attention mechanism, represents the basic features obtained by fusion in the frequency domain, represents the detailed features obtained by fusion in the frequency domain; Step 3.4.3, convert the fused features in the frequency domain and back to the time domain again, so as to obtain the fused basic features and the fused detail features in the time domain; ； Among them, represents a two-dimensional inverse Fourier transform operation, represents the fused basic features in the time domain, represents the fused detailed features in the time domain; Step 3.4.4, input the fused base features in the time domain and the fused detail features into the decoder to generate the finally fused image features , which is expressed as: ； Among them, is the decoded fused feature, represents the decoder; Step 4, based on the fused image features , predict the snow cover distribution.

2. The method for predicting snow cover distribution according to claim 1, wherein The preprocessing of the remote sensing data of the SAR image and the multispectral optical (opt) image in Step 1 specifically includes: Step 1.1: Perform orbital correction, radiometric calibration, speckle noise filtering, geocorrection, 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; Step 1.2: Process the opt image data, specifically: First, perform band selection on the opt image data to obtain the selected bands , which is expressed as: ； Among them, 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; For the selected bands perform radiometric calibration, and then perform weighted average fusion on the selected and to obtain the opt image data after the first preprocessing, expressed as: ； Among them, is the image after band fusion, is the weight of the th band, and are the radiometrically calibrated bands B2, B3, B4, and B8; Step 1.3: Register the SAR image data after the first preprocessing in Step 1.1 with the opt image data after the first preprocessing in Step 1.2, and resample the registered image using the Cubic Convolution method to obtain the preprocessed SAR image and the preprocessed multispectral optical (opt) image in Step 1.

3. The method for predicting snow cover distribution according to claim 2, wherein Extracting snow-covered area image data based on the multi-spectral optical opt image data described in step 2 , specifically as follows: Using the green light band and the shortwave infrared band Calculate the Normalized Difference Snow Index (NDSI): ； Among them, and are the reflectivities in the green light band and the short-wave infrared band respectively; Set a threshold to extract snow cover T1, and extract the snow cover area based on the NDSI value , and obtain the image data of the snow cover area : 。 4. The snow accumulation distribution prediction method according to claim 1, wherein In step 4, based on the fused image features , predict the snow cover distribution; 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; Extract a fixed-length input window and a label window from the fused image features using the sliding window method and input them into the ConvLSTM module branch and the 3D-CNN module branch respectively for prediction; The ConvLSTM module branch includes three stacked ConvLSTM units. Each layer of ConvLSTM unit uses a 3×3 convolutional kernel to extract spatio-temporal features and update the gated state in the ConvLSTM unit. The output of the last layer of ConvLSTM unit 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 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; 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 .

5. A snow accumulation distribution prediction system, characterized in that, It includes the following modules: Snow data preprocessing module, which is used to perform orbit correction, radiometric calibration, speckle noise filtering, geocorrection, and conversion of the backscattering coefficient to 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, and select bands , , , and , 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 , , and , register with the SAR image data after the first preprocessing to obtain the preprocessed synthetic aperture radar (SAR) image data and multispectral optical (opt) image data; A snow accumulation distribution information extraction module, configured to extract snow-covered area images from the preprocessed synthetic aperture radar (SAR) image data and multispectral optical (opt) image data respectively and snow-covered area image data ; The multi-modal information fusion module fuses the snow-covered area image and the snow-covered area image data to obtain the fused image features; A snow cover distribution prediction module that predicts the snow cover distribution based on the fused image features; The multi-modal information fusion module includes a shared feature encoder (SFE), a basic Transformer encoder (BTE), a detail CNN encoder (DCE), a frequency domain interaction module, and an encoder; Shared Feature Encoder (SFE) for separately extracting the snow-covered area image and the snow-covered area image data to obtain the shallow features, resulting in shallow features and shallow features ; The basic Transformer encoder BTE is used to extract low-frequency basic features from shallow features and the low-frequency basic features extracted from shallow features ; and the low-frequency basic features extracted from shallow features ; The detailed CNN encoder DCE is used to extract high-frequency detailed features from shallow features and high-frequency detailed features extracted from shallow features ; and high-frequency detailed features extracted from shallow features ; A frequency-domain interaction module, which 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 to the time domain to obtain the fused basic features and the fused detail features; An encoder decodes the fused base features and the fused detail features to obtain the finally fused image features ; 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 , the low-frequency basic features , the high-frequency detail features , the high-frequency detail features to the frequency domain respectively; The fusion unit 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; The inverse Fourier transform unit is used to perform an 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.

6. The snow accumulation distribution prediction system according to claim 5, characterized in that, The snow accumulation distribution prediction module includes a ConvLSTM module branch and a 3D-CNN module branch; Extract a fixed-length input window and label window from the fused image features using the sliding window method and input them into the ConvLSTM module branch and the 3D-CNN module branch for prediction respectively; The ConvLSTM module branch includes three stacked ConvLSTM units. Each layer of the ConvLSTM unit uses a 3×3 convolution kernel to extract spatio-temporal features and update the gating states in the ConvLSTM unit. The output of the last layer of the ConvLSTM unit 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 the joint features of the time and space dimensions, and reduces the dimension of the joint features through a three-dimensional pooling operation. 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; Perform weighted averaging on the prediction results of the ConvLSTM module branch and the prediction results of the 3D-CNN module branch to generate the final snow accumulation distribution prediction result .

Images