A multi-source remote sensing spatio-temporal fusion method for Bohai coastal sea ice parameters, medium and system

Abstract

This invention provides a multi-source remote sensing spatiotemporal fusion method, medium, and system for sea ice parameters in the Bohai Sea nearshore area, belonging to the field of sea ice analysis technology. This invention utilizes an iterative optimal transmission spatiotemporal registration algorithm constrained by ice drift physical constraints to achieve sub-pixel-level spatiotemporal alignment of three-source data. The results are input into an artificial intelligence model that guides a multi-scale deformation map convolutional fusion based on ice type priors, simultaneously outputting a sea ice distribution probability map and a preliminary sea ice concentration estimate map. Then, an iterative posterior expectation map and a confidence interval map are output through a sparse Bayesian multi-source ice concentration joint inversion algorithm. Finally, after logical consistency verification, median filtering, and shoreline clipping, sea ice distribution raster products and sea ice concentration raster products are generated. This invention solves the technical problem that multi-source remote sensing data fusion cannot effectively model nonlinear spatiotemporal correlations.

Description

Technical Field

[0001] This invention belongs to the field of sea ice analysis technology, specifically, it relates to a multi-source remote sensing spatiotemporal fusion method, medium and system for sea ice parameters in the Bohai Sea nearshore area. Background Technology

[0002] In the field of remote sensing monitoring of nearshore sea ice in the Bohai Sea, existing technologies rely on three types of platforms for sea ice parameter inversion: optical satellites, UAV multispectral systems, and shore-based X-band radars. Optical satellites provide macroscopic coverage and rich spectral information, UAVs can acquire detailed nearshore ice conditions at small scales, and shore-based radars have all-weather, high-frequency, continuous observation capabilities. All three are widely used in operational sea ice monitoring. At the fusion level, existing methods mainly employ statistical strategies such as linear interpolation, weighted averaging, and Kalman filtering to combine multi-source observation results in a deterministic linear manner.

[0003] However, the aforementioned linear statistical fusion methods have inherent limitations. First, different sensors differ significantly in imaging mechanisms, spatial resolution, and noise characteristics, making simple weighting unable to adaptively handle the complementary relationships between heterogeneous data. Second, there is a complex nonlinear response between sea ice concentration and multi-source observation signals, making it difficult for linear models to accurately characterize the mapping relationships under complex scenarios such as varying ice types, ice-water mixing, and cloud obstruction. Third, existing fusion results lack spatiotemporal smoothness under physical constraints, frequently exhibiting fusion artifacts and temporal jumps in nearshore broken ice and thin ice regions.

[0004] In other words, existing technologies have a technical problem: multi-source remote sensing data fusion cannot effectively model nonlinear spatiotemporal correlations. Summary of the Invention

[0005] In view of this, the present invention provides a method, medium and system for multi-source remote sensing spatiotemporal fusion of sea ice parameters for the Bohai Sea nearshore area, which can solve the technical problem in the prior art that multi-source remote sensing data fusion cannot effectively model nonlinear spatiotemporal correlations.

[0006] The present invention is implemented as follows: The first aspect of the present invention provides a multi-source remote sensing spatiotemporal fusion method for sea ice parameters in the Bohai Sea nearshore area, comprising the following steps:

[0007] We acquired optical images from the Haiyang-1 C / D satellite, multispectral images from UAVs, and backscattering maps from shore-based X-band radar. We performed radiometric calibration, atmospheric correction, cloud masking, and clutter suppression on the three-source data, then uniformly projected and resampled them to obtain the three-source preprocessed feature tensors.

[0008] Using the transit time of Haiyang-1C / D satellite as the time anchor point, and based on the ice drift physical constraint iterative optimal transmission spatiotemporal registration algorithm, the spatiotemporal alignment transformation field of multi-source data is calculated for the three-source preprocessed feature tensor, and the three-source preprocessed feature tensor is transformed into a three-source aligned feature tensor.

[0009] The normalized difference ice and snow index, near-infrared brightness ratio index, gray-level co-occurrence matrix texture features, and Gabor multi-directional filtering response are calculated for the three-source aligned feature tensor to obtain the three-source enhanced feature tensor. Based on the thresholds of cloud coverage index and ice drift amplitude index, low-light enhancement processing is performed on the optical data part or low-light enhancement processing is skipped.

[0010] The three-source enhanced feature tensor is input into the ice type prior-guided multi-scale deformation map convolutional fusion model. After superpixel map construction, deformation map convolutional aggregation and coarse-fine two-level map iterative refinement, the sea ice distribution probability map and the initial sea ice concentration map are output simultaneously.

[0011] Based on the sparse Bayesian multi-source ice concentration joint inversion algorithm, the initial sea ice concentration map is used as the initial value. The three-source aligned feature tensor is fused to iteratively output the posterior expectation map of sea ice concentration and the confidence interval map of sea ice concentration.

[0012] A distribution threshold is applied to the sea ice distribution probability map to generate a binary sea ice mask. Median filtering is applied to the posterior expectation map of sea ice concentration. Logical consistency check is performed. Land areas are clipped based on shoreline vectors. Finally, sea ice distribution raster products and sea ice concentration raster products are output.

[0013] Specifically, the unified projection and resampling involves projecting onto the WGS84UTMZone51N coordinate system and resampling to 50m pixels.

[0014] Specifically, the time difference between the three data sources is controlled within a time difference threshold of 1 hour, and the spatial registration error is less than 0.5 pixels.

[0015] Specifically, the ice drift physical constraint iterative optimal transport spatiotemporal registration algorithm is based on the optimal transport theory framework. It uses the cross-correlation method to estimate the initial ice drift velocity field, constructs a transport cost matrix that integrates Euclidean distance and flow velocity divergence minimization constraints, uses the Sinkhorn-Knopp iterative algorithm to solve the regularized optimal transport mapping, introduces anisotropic diffusion constraints, and stops iterating when the change in the transport mapping is lower than the convergence threshold.

[0016] The clutter suppression process includes ice wave separation processing, which specifically involves extracting pixel-level Doppler spectral features through short-time Fourier transform, using the Doppler frequency shift range of ocean waves and the upper limit of ice drift frequency as separation thresholds, and inputting them into a one-dimensional lightweight convolution classifier for fine separation of ice, water, and wave mixtures.

[0017] The formula for calculating the near-infrared luminance ratio index is as follows: ,in For near-infrared reflectivity, For red band reflectivity, This refers to the reflectivity in the green band.

[0018] The cloud coverage rate index is calculated from the pixel ratio of the cloud mask in the optical image. When the cloud coverage rate index is greater than the cloud coverage rate threshold, the low-light enhancement process is skipped. When the cloud coverage rate index is not greater than the cloud coverage rate threshold, the low-light enhancement process based on the radiative transfer model is performed. The cloud coverage rate threshold is 0.6.

[0019] The input layer of the ice-type prior-guided multi-scale deformation map convolutional fusion model receives a three-source enhanced feature tensor with a spatial size of 512×512 pixels. The superpixel map construction module divides superpixel nodes at a coarse resolution of 50m, and the edge weights between nodes are defined by spatial adjacency and ice dynamics correlation.

[0020] Specifically, the deformation map convolutional aggregation is achieved by dynamically adjusting the sampling positions of neighboring nodes using a learnable offset field. This learnable offset field is predicted from local features by a lightweight convolutional network and optimized end-to-end during training.

[0021] The ice type prior guidance layer utilizes the historical ice type distribution prior probability map and injects context vectors of four types of ice types—gray ice, white ice, fixed ice, and broken ice—into each deformation map convolutional layer through conditional normalization. The scaling and offset parameters of the conditional normalization are generated by the historical ice type distribution prior probability map through a fully connected layer.

[0022] The training of the ice-type prior-guided multi-scale deformation map convolutional fusion model adopts a composite loss function. The sea ice distribution head adopts Dice loss, and the density regression head adopts smooth L1 loss. The weighting coefficients of the two are automatically updated by the dynamic loss weight adjustment function. The sum of the two weighting coefficients is always 1, and the value of each weighting coefficient is limited to the lower limit of the weighting coefficient to the upper limit of the weighting coefficient.

[0023] The dynamic loss weight adjustment function is based on the fusion imbalance index. ,when Increase the Dice loss weighting coefficient when it is not less than the upper imbalance threshold. When the value is less than the lower imbalance threshold, the weighting coefficient of the smoothing L1 loss is increased; otherwise, it remains unchanged. The adjustment amount each time is the adjustment step size.

[0024] Specifically, the sparse Bayesian multi-source ice concentration joint inversion algorithm expands the sea ice concentration field into a linear superposition of an overcomplete dictionary composed of historical typical ice condition samples, uses evidence maximization to iteratively estimate the hyperparameters of each basis function, automatically drives the weights of unrelated basis functions to approach zero, and updates the posterior estimates of the sea ice concentration field and the variance of the observation noise of each source at each iteration.

[0025] The sea ice concentration confidence interval map is the posterior standard deviation map of sea ice concentration output by the sparse Bayesian learning framework. The larger the posterior standard deviation, the lower the prediction confidence of the corresponding pixel.

[0026] The distribution threshold is 0.5, the median filter window is 3×3, the lower limit of the weighting coefficient is 0.3, the upper limit of the weighting coefficient is 0.7, the adjustment step size is 0.05, the upper imbalance threshold is 1.5, the lower imbalance threshold is 0.67, the wave Doppler frequency shift range is 0.1 to 0.3 Hz, and the upper limit of the ice drift frequency is 0.01 Hz.

[0027] A second aspect of the present invention provides a computer-readable storage medium storing program instructions, which, when executed in a computer, are used to perform the above-described method for multi-source remote sensing spatiotemporal fusion of sea ice parameters for the Bohai Sea nearshore area.

[0028] A third aspect of the present invention provides a multi-source remote sensing spatiotemporal fusion system for sea ice parameters in the Bohai Sea nearshore area, comprising the aforementioned computer-readable storage medium, wherein the system is a computer, the computer-readable storage medium is disposed within the system, and the system is provided with a microprocessor for executing program instructions stored in the computer-readable storage medium.

[0029] This invention employs an iterative optimal transmission spatiotemporal registration algorithm constrained by ice drift physical constraints to achieve geometrically consistent alignment of multi-source data at the sub-pixel level. It utilizes an ice-type prior-guided multi-scale deformation map convolutional fusion model to model the fusion problem as a graph reasoning task. Furthermore, it outputs a posterior expectation map and a confidence interval map of sea ice concentration through a sparse Bayesian multi-source ice concentration joint inversion algorithm, thus solving the technical problem that multi-source remote sensing data fusion cannot effectively model nonlinear spatiotemporal correlations.

[0030] This invention leverages a graph structure that naturally adapts to the topological relationships of sea ice spatial distribution. Deformation graph convolution dynamically changes the receptive field according to the ice block shape, overcoming the limitation of traditional linear methods in describing complex nonlinear mappings. An ice-type prior guidance layer embeds domain-specific prior knowledge into the feature propagation process, enabling differentiated aggregation strategies for different ice-type regions. Physical constraints ensure the spatiotemporal rationality of the fusion results. The sparse Bayesian framework, through an automatic correlation determination mechanism, automatically reduces the corresponding source weights in areas of low observation quality and concentrates weights in areas of high observation quality, thus outputting physically interpretable inversion results even under degradation conditions such as cloud cover and radar clutter.

[0031] In summary, this invention solves the technical problem mentioned in the background art that multi-source remote sensing data fusion cannot effectively model nonlinear spatiotemporal correlations. Attached Figure Description

[0032] Figure 1 This is a flowchart of the method of the present invention.

[0033] Figure 2 This is a schematic diagram of the spatial distribution of sea ice probability map and preliminary sea ice concentration estimate map in the nearshore study area of ​​Liaodong Bay, Bohai Sea.

[0034] Figure 3 This is a time series curve showing the changes in sea ice concentration in the study area from January 5 to January 10, 2026.

[0035] Figure 4 This is a spatial distribution map of the sparse Bayesian inversion weights of the three source data in each pixel of the study area. Detailed Implementation

[0036] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions in the embodiments of the present invention will be clearly and completely described below.

[0037] like Figure 1 The diagram shown is a flowchart of a multi-source remote sensing spatiotemporal fusion method for sea ice parameters in the Bohai Sea nearshore area, provided by the first aspect of this invention. This method includes the following steps:

[0038] S01. Acquire optical images from Haiyang-1 C / D satellite, multispectral images from UAVs, and backscattering maps from shore-based X-band radar. Perform radiometric calibration, atmospheric correction, cloud masking, and clutter suppression on the three-source data respectively. Project them uniformly to the WGS84UTMZone51N coordinate system and resample them to 50m pixels to obtain the three-source preprocessed feature tensors.

[0039] S02. Using the transit time of Haiyang-1 C / D satellite as the time anchor point, the time difference of the three-source data is controlled within 1 hour. Based on the ice drift physical constraint iterative optimal transmission spatiotemporal registration algorithm, the spatiotemporal alignment transformation field of the three-source preprocessed feature tensor is calculated, and the three-source preprocessed feature tensor is transformed into a three-source aligned feature tensor. The spatial registration error is less than 0.5 pixels.

[0040] S03. Calculate the normalized difference ice and snow index, near-infrared brightness ratio index, gray-level co-occurrence matrix texture features and Gabor multi-directional filtering response for the three-source aligned feature tensor to obtain the three-source enhanced feature tensor; based on the thresholds of cloud coverage index and ice drift amplitude index, perform low-light enhancement processing or skip low-light enhancement processing on the optical data part of the three-source enhanced feature tensor.

[0041] S04. Input the three-source enhanced feature tensor into the ice type prior-guided multi-scale deformation map convolutional fusion model. After superpixel map construction, deformation map convolutional aggregation and coarse-fine two-level map iterative refinement, the sea ice distribution probability map and the initial sea ice concentration estimate map are output simultaneously.

[0042] S05. Based on the sparse Bayesian multi-source ice concentration joint inversion algorithm, using the initial sea ice concentration map as the initial value, the optical reflectivity, radar backscattering and UAV multispectral index in the three-source aligned feature tensor are fused to iteratively output the posterior expectation map of sea ice concentration and the confidence interval map of sea ice concentration.

[0043] S06. Apply a threshold of 0.5 to the sea ice distribution probability map to generate a sea ice binary mask, apply a 3×3 window mid-range filter to the sea ice concentration posterior expectation map, perform a logical consistency check, clip the land area based on the shoreline vector, and finally output the sea ice distribution raster product and the sea ice concentration raster product.

[0044] The principle and implementation of the ice drift physical constraint iterative optimal transmission spatiotemporal registration algorithm are as follows: Based on the optimal transmission theory, multi-temporal ice field registration is modeled as a metric space mapping problem with physical constraints. First, based on the time-series backscatter map of shore-based X-band radar, the initial ice drift velocity field is estimated using the cross-correlation method, thereby constructing a transmission cost matrix. The cost function simultaneously integrates Euclidean distance and ice dynamic constraints, where the ice dynamic constraint is the minimization of current velocity divergence. The Sinkhorn-Knopp iterative algorithm is used to solve the regularized optimal transmission mapping. In each iteration, the pixel correspondence is updated and the ice drift velocity field estimation is corrected in reverse. Anisotropic diffusion constraints are introduced to force the ice drift velocity in the near-shore fixed ice area to zero, while allowing large displacement transmission in open water. Iteration stops when the change in the transmission mapping is lower than the convergence threshold, and the spatiotemporal alignment transformation field of multi-source data is output. The convergence threshold was obtained through multiple iterative experiments on historical shore-based X-band radar time-series data of the Bohai Sea. Specifically, for 100 sets of ice field sequences with known true values, transmission mapping changes of 0.01 pixels, 0.05 pixels, and 0.1 pixels were set respectively. The trade-off between registration error and the number of iterations was statistically analyzed, and the transmission mapping change corresponding to when the registration error no longer decreased significantly was taken as the convergence threshold. The technical effect of the algorithm is that traditional feature point matching methods fail in broken ice scenarios due to texture repetition and feature sparsity. However, the optimal transmission framework transforms registration into a global probability allocation problem. Physical constraints enable adaptive switching of alignment strategies between nearshore fixed ice and drifting ice regions, thereby obtaining physically reasonable sub-pixel-level registration results even in low-texture regions such as broken ice and fractured ice edges, providing geometrically consistent input for subsequent multi-source fusion.

[0045] The formula for calculating the normalized difference snow index is as follows: ;in For green band reflectivity, Both are dimensionless reflectance values, representing the shortwave infrared band reflectance.

[0046] The formula for calculating the near-infrared luminance ratio index is as follows: ;in For near-infrared reflectivity, For red band reflectivity, The reflectance is measured in the green band, and all three values ​​are dimensionless reflectance values. The near-infrared brightness ratio index utilizes the difference between the strong absorption characteristics of suspended sediment in the 865nm band and the high reflectance of thin ice in the 865nm band to construct a distinguishing ability between turbid water and thin ice, compensating for the problem of missed detection of thin ice due to spectral overlap caused by the normalized difference ice and snow index in the Liaodong Bay environment with high suspended sediment. The classification threshold of the near-infrared brightness ratio index was obtained through the following experiment: Ten sets of on-site synchronous spectral measurement samples were collected in Liaodong Bay during winter, and the near-infrared brightness ratio index values ​​were measured for thin ice, turbid water, and open water respectively. Distribution histograms were plotted, and the near-infrared brightness ratio index value at the point of least overlap between the two distributions was taken as the classification threshold. Its stability was confirmed by 5-fold cross-validation.

[0047] The cloud coverage rate index is calculated from the pixel ratio of the cloud mask in the optical image and is a dimensionless proportional value. The ice drift amplitude index is estimated by cross-correlation of the backscattering maps of adjacent time points from the shore-based X-band radar and is in m / s. When the cloud coverage rate index is greater than 0.6, low-light enhancement processing is skipped for the optical data portion of the three-source enhancement feature tensor; when the cloud coverage rate index is not greater than 0.6, low-light enhancement processing based on the radiative transfer model is performed on the optical data portion of the three-source enhancement feature tensor. The cloud coverage rate index threshold of 0.6 was obtained through statistical analysis of historical satellite archived data from 10 winters in the Bohai Sea from 2015 to 2024. Specifically, the change curves of sea ice identification accuracy under different cloud coverage rate index conditions were statistically analyzed, and the cloud coverage rate index value at the inflection point where the accuracy significantly decreased was taken as the threshold.

[0048] The specific structure of the ice-type prior-guided multi-scale deformation map convolutional fusion model is as follows: The input layer receives three-source enhancement feature tensors, corresponding to the optical enhancement features of the Haiyang-1 C / D satellite, the multispectral enhancement features of the UAV, and the backscattering enhancement features of the shore-based X-band radar, respectively. The spatial size of each tensor is 512×512 pixels, and the number of channels is determined according to the number of each source band. The superpixel map construction module divides the study area into several superpixel nodes with a coarse resolution of 50m. The node features are formed by multi-scale convolutional encoding and concatenation of the three-source enhancement feature tensors at corresponding spatial positions. The edge weights between nodes are defined by the spatial adjacency relationship and the ice dynamic correlation. The ice dynamic correlation is quantized by the ice drift velocity field estimated by the shore-based X-band radar. The deformation map convolution module dynamically adjusts the sampling positions of neighboring nodes at each superpixel node using a learnable offset field, so that the map receptive field adapts to the shape and degree of ice fragmentation. The learnable offset field is predicted from local features by a lightweight convolutional network, and the learnable offset is optimized end-to-end during training. The ice type prior guidance layer utilizes a prior probability map of historical ice type distribution. Through conditional normalization, it injects contextual vectors of four ice types—gray ice, white ice, fixed ice, and broken ice—into each deformation map convolutional layer, guiding the network to adopt differentiated feature propagation strategies in different ice type regions. The scaling and offset parameters of the conditional normalization are generated from the prior probability map of historical ice type distribution via a fully connected layer. The coarse-to-fine two-level graph iterative refinement module first completes the first round of graph inference at a coarse resolution of 50m, obtaining a coarse-scale sea ice distribution probability map and a preliminary estimate of sea ice concentration. Subsequently, it dynamically updates the superpixel node features and reconstructs the edge connections between nodes using UAV multispectral enhancement features. Within the coverage area of ​​the UAV multispectral image, it refines the graph resolution to the original resolution of the UAV multispectral image, completing the second round of graph inference and forming a multi-resolution fused representation. The decoder uses a feature pyramid network structure for feature upsampling and multi-scale fusion, restoring the output resolution to 50m. The output layer features a dual-head design: the sea ice distribution head outputs a single-channel sea ice distribution probability map, activated by a Sigmoid function; the density regression head outputs a single-channel preliminary sea ice density estimate map, linearly mapped to 0% to 100%. The ice-type prior-guided multi-scale deformation graph convolutional fusion model models the fusion problem as a graph reasoning task. The graph structure naturally adapts to the topological relationships of sea ice spatial distribution. Deformation graph convolution allows the receptive field to dynamically change with the shape of the ice block. The ice-type prior-guided layer embeds domain knowledge into the feature propagation process. Iterative refinement of the coarse and fine graphs enhances the macroscopic distribution and local detail information, thus obtaining physically reasonable and spatially detailed fusion results in nearshore broken ice, thin ice, and ice-water mixed areas. This overcomes the limitation of traditional weighted average methods in modeling nonlinear spatiotemporal correlations.

[0049] The specific steps for establishing the training dataset for the ice-type prior-guided multi-scale deformation map convolutional fusion model include: collecting 12 observation records from icebreaker voyages in the Bohai Sea during winter from 2015 to 2024, acquiring concurrent optical images from the Haiyang-1 C / D satellite, multispectral images from UAVs, and backscatter maps from shore-based X-band radar, constructing no fewer than 200 sets of three-source synchronous observation samples; labeling data according to the World Meteorological Organization's sea ice concentration level 7 standard, with experts manually delineating ice zone boundaries on high-resolution images, combining the normalized difference ice and snow index with morphological semi-automatic segmentation results, and obtaining a true sea ice concentration map after manual correction; specifically expanding the sample of mixed thin ice and turbid water scenes to 30% of the total sample size; dividing the training set, validation set, and test set in a 7:2:1 ratio, covering all stages of the initial glacial period, the peak glacial period, and the melting glacial period.

[0050] The training steps of the ice-type prior-guided multi-scale deformation map convolutional fusion model specifically include: The first stage involves pre-training on 10,000 sets of multi-source synthetic data generated by a physical simulation platform. This multi-source synthetic data covers different ice types, varying cloud cover, and radar noise scenarios with different signal-to-noise ratios. The second stage involves fine-tuning on the Bohai Sea measured dataset. During training, geometric enhancement (±30° rotation and horizontal / vertical flipping), photometric perturbations (±20% each for brightness and contrast), and random block masks to simulate cloud occlusion are applied in real time. A composite loss function is used: Dice loss for the sea ice distribution head and smoothed L1 loss for the density regression head. The weighting coefficients of both are automatically updated during training by a dynamic loss weight adjustment function. The optimizer is AdamW, with an initial learning rate of... The weight decays to The batch size is 8, the input block size is 512×512 pixels, the maximum training is 200 rounds, and early stopping is triggered if the validation set loss does not decrease within 15 rounds.

[0051] The dynamic loss weight adjustment function is used to automatically adjust the weighting coefficients of the Dice loss of the sea ice distribution head and the smoothed L1 loss of the concentration regression head during training. The dynamic loss weight adjustment function calculates the fusion imbalance index based on the Dice loss value of the sea ice distribution head and the smoothed L1 loss value of the concentration regression head on the validation set of the current training round. The calculation formula of the fusion imbalance index is expressed as follows: ;in This represents the Dice loss value for the current sea ice distribution. The initial baseline value for the Dice loss of the sea ice distribution head was set at the end of the first round of training. This represents the L1 loss value of the head smoothing regression for the current round's density. The initial baseline value for the dense regression head smoothing L1 loss is set at the end of the first round of training. This is a dimensionless index for the degree of fusion imbalance. When When the sea ice distribution head Dice loss weighting coefficient is increased by 0.05, the concentration regression head smoothing L1 loss weighting coefficient is correspondingly decreased by 0.05; when At that time, the weighted coefficients of the Dice loss in the sea ice distribution head and the L1 loss in the concentration regression head remain unchanged; when At that time, the weighting coefficient of the density regression head smoothing L1 loss is increased by 0.05, and the weighting coefficient of the sea ice distribution head Dice loss is correspondingly decreased by 0.05; the sum of the two weighting coefficients is always kept at 1, and the value of each weighting coefficient is limited to between 0.3 and 0.7. The initial baseline values ​​of the sea ice distribution head Dice loss and the density regression head smoothing L1 loss are automatically recorded after the first round of training. The fusion imbalance index thresholds of 1.5 and 0.67 are obtained by analyzing 10 training experiments with different initial weighting coefficient configurations. Specifically, the convergence speed and final accuracy of the validation set accuracy are statistically analyzed under each configuration, and the boundary of the fusion imbalance index interval that makes the convergence speed of the two tasks most balanced is taken as the threshold.

[0052] The principle and implementation of the sparse Bayesian multi-source ice concentration joint inversion algorithm are as follows: A sparse Bayesian learning framework is used to probabilistically model the mapping relationship between multi-source observations and sea ice concentration in the three-source aligned feature tensor. The sea ice concentration field is expanded into a linear superposition of an overcomplete dictionary, which consists of historical typical ice condition samples. Optical reflectivity, radar backscattering, and UAV multispectral index are used as multi-view measurement vectors. Evidence maximization is used to iteratively estimate the hyperparameters of each basis function, automatically driving the weights of unrelated basis functions to approach zero to achieve automatic sparsification. Each iteration simultaneously updates the posterior estimates of the sea ice concentration field and the variance of noise from each source observation, ultimately outputting a posterior expectation map and a confidence interval map of sea ice concentration. The technical advantage of the algorithm is that traditional regression methods treat multi-source observations as equally weighted deterministic mappings, which cannot quantify the differences in observation quality and the uncertainty of sea ice concentration estimation. In contrast, the sparse Bayesian learning framework uses an automatic correlation determination mechanism to automatically reduce the weight of the corresponding source and expand the confidence interval of sea ice concentration in areas with low observation quality, while concentrating the weights in areas with high observation quality to improve accuracy. Thus, even under observation degradation conditions such as cloud cover and radar clutter, it can still output physically interpretable posterior expectation maps and confidence interval maps of sea ice concentration.

[0053] The clutter suppression processing of the shore-based X-band radar backscatter map includes ice wave separation processing. This ice wave separation processing utilizes the time-varying texture characteristics of the shore-based X-band radar time-series backscatter map for differentiation. The wave Doppler frequency shift ranges from 0.1 to 0.3 Hz, and the ice drift frequency is less than 0.01 Hz. Pixel-level Doppler spectral features are extracted using short-time Fourier transform and input into a one-dimensional lightweight convolution classifier for fine separation of the ice, water, and wave mixture, resulting in a shore-based X-band radar backscatter map free from wave Bragg scattering interference. The upper limit of the ice drift frequency (0.01 Hz) and the wave Doppler frequency shift range (0.1 to 0.3 Hz) are obtained through spectral statistical analysis of historical time-series backscatter maps of the Bohai Sea during winter. Specifically, 50 time-series samples each from known ice areas and known open water areas are selected, and short-time Fourier transform power spectra are extracted for each. The peak distribution intervals of the two target spectra are statistically analyzed, and the frequency boundaries of the two distributions without overlap are taken as the separation threshold.

[0054] The superpixel node refers to a local homogeneous region unit formed by aggregating the study area grid according to spatial proximity and spectral similarity, serving as the basic computational node in graph inference. The learnable offset field refers to a two-dimensional offset field predicted from local features by a lightweight convolutional network, allowing the sampling position of the deformation map convolution to dynamically float on the feature map, thus enabling the receptive field shape to adapt to the geometric changes of the target. Conditional normalization refers to standardizing the input features of the deformation map convolutional layer by calculating the mean and variance per channel, and then performing an affine transformation on the standardized result using scaling and offset parameters generated from the historical ice type distribution prior probability map, aligning the feature distribution of different ice type regions to the semantic space of the corresponding ice type. The feature pyramid network structure refers to a decoder structure that progressively fuses deep low-resolution semantic features with shallow high-resolution detail features through top-down paths and lateral connections, outputting a multi-scale feature map. The sea ice concentration confidence interval map is a posterior standard deviation map of sea ice concentration output by a sparse Bayesian learning framework, reflecting the uncertainty of the posterior expected value of sea ice concentration for each pixel. A larger posterior standard deviation indicates a lower prediction confidence for that pixel. The historical ice type distribution prior probability map is obtained from statistical analysis products of historical sea ice in the Bohai Sea from 2015 to 2024, reflecting the historical frequency of each pixel in the study area belonging to one of the four ice types: gray ice, white ice, fixed ice, and broken ice. The shoreline vector is extracted from the land-water boundary of the Haiyang-1 C / D satellite optical imagery and is used to crop the land area in the sea ice distribution raster product and the sea ice concentration raster product.

[0055] The specific implementation of step S01 is as follows: When acquiring optical images from the Haiyang-1 C / D satellite, the original digital quantization values ​​are first radiometrically calibrated according to the sensor calibration coefficients to convert the pixel values ​​into radiance. Then, an atmospheric correction model based on the radiative transfer equation is used, with input auxiliary parameters such as aerosol optical thickness and water vapor content, to convert the radiance into surface reflectance. Cloud masking uses a combination of reflectance thresholds from the blue band and shortwave infrared band to identify and mark cloud-covered pixels, generating an effective observation mask. For UAV multispectral images, radiometric calibration is first performed on each band using a randomly mounted radiometric calibration board to eliminate sensor response inconsistencies. Then, orthorectification is performed using flight attitude and position data to eliminate geometric distortions caused by terrain and flight attitude. The backscattering map from the onshore X-band radar is converted into a backscattering intensity map using a range Doppler imaging algorithm. Multi-look processing is employed to reduce speckle noise, followed by short-time Fourier transform to extract pixel-level Doppler spectral features. Using the upper limit of ice drift frequency (0.01 Hz) and the wave Doppler frequency shift range (0.1 to 0.3 Hz) as separation boundaries, a one-dimensional lightweight convolution classifier is input to separate ice and waves, removing wave Bragg scattering interference. After processing the three-source data, they are uniformly projected onto the WGS84UTMZone51N coordinate system and resampled to 50m pixels using bilinear interpolation, forming a spatially consistent three-source preprocessed feature tensor. This provides geometrically normalized input for subsequent spatiotemporal registration.

[0056] The specific implementation of step S02 is as follows: using the transit time of the Haiyang-1 C / D satellite as the time anchor point, the flight window of the UAV and the data collection period of the shore-based radar are both controlled within one hour before and after the satellite transit to ensure that the ice condition described by the three sources of data is approximately synchronized. The core principle of the ice drift physical constraint iterative optimal transmission spatiotemporal registration algorithm is the optimal transmission theory, which models the multi-temporal ice field registration as a probability allocation problem in the measurement space. The specific implementation process is as follows: First, based on the time-series backscatter map of the shore-based X-band radar, the correlation coefficient peak is calculated within a local window using the cross-correlation method to estimate the initial ice drift velocity field. A transmission cost matrix is ​​constructed using the initial velocity field. The cost function integrates the Euclidean spatial distance and the ice dynamics constraint of minimizing the current velocity divergence. The latter penalizes physically unreasonable large divergence drift patterns. The Sinkhorn-Knopp iterative algorithm is used to solve the regularized optimal transmission problem. In each iteration, the pixel correspondence is updated and the ice drift velocity field estimation is corrected in reverse. Anisotropic diffusion constraints are introduced to force the drift velocity of the nearshore fixed ice area to zero. Large displacement transmission is allowed in open water, enabling the algorithm to adaptively switch the registration strategy in different ice types. When the change in transmission mapping is lower than the convergence threshold, the iteration stops, and the spatiotemporal alignment transformation field of the multi-source data is output. The three-source preprocessed feature tensors are spatially remapped according to this transformation field to obtain a three-source aligned feature tensor with a spatial registration error of less than 0.5 pixels.

[0057] The specific implementation of step S03 is as follows: For the optical and UAV multispectral components in the three-source aligned feature tensor, first calculate the normalized difference ice and snow index, using the formula: Ice-water separation is achieved by utilizing the spectral characteristics of ice and snow—high reflectance in the green band and low reflectance in the short-wave infrared band; the near-infrared brightness ratio index is calculated using the formula: This study utilizes the difference between the strong absorption characteristics of suspended sediment in the near-infrared band and the high near-infrared reflectivity of thin ice to compensate for the under-detection problem of the normalized difference ice and snow index in environments with high suspended sediment. It calculates the gray-level co-occurrence matrix texture features, extracting statistics such as contrast, energy, and correlation to describe the spatial gray-level relationships within the pixel neighborhood. Gabor multi-directional filtering is applied to multi-band images to extract edge responses in four directions: 0°, 45°, 90°, and 135°, enhancing the discernibility of ice fragment boundaries. The above indices and features are concatenated to form a three-source enhanced feature tensor. Subsequently, the cloud coverage index, i.e., the proportion of cloud mask pixels in the optical image, is calculated; and the ice drift amplitude index, i.e., the ice drift velocity estimated by the cross-correlation of backscatter maps from adjacent time points using shore-based X-band radar, is calculated. When the cloud coverage index is greater than 0.6, the optical data portion skips low-light enhancement processing to avoid introducing false information due to enhancement under cloud cover; when the cloud coverage index is not greater than 0.6, low-light enhancement processing based on the radiative transfer model is performed on the optical data portion to compensate for the reduced signal-to-noise ratio caused by the low solar altitude angle in the Bohai Sea during winter.

[0058] The specific implementation of step S04 is as follows: After inputting the three-source enhanced feature tensor into the ice-type prior-guided multi-scale deformation map convolutional fusion model, the superpixel map construction module divides the study area into several superpixel nodes based on spatial proximity and spectral similarity. The node features are constructed by concatenating the multi-scale convolutional codes of the three-source enhanced feature tensor at corresponding positions. The edge weights between nodes comprehensively consider the spatial adjacency relationship and the correlation with ice dynamics. The correlation with ice dynamics is quantified by the ice drift velocity field estimated in step S02, thereby making the graph structure reflect the actual physical topology of the ice field. The deformation map convolutional aggregation module uses a lightweight convolutional network at each superpixel node to predict a two-dimensional learnable offset field from local features, so that the sampling positions of neighboring nodes dynamically float on the feature map, making the effective receptive field of the graph adapt to the geometry and fragmentation of the ice block, overcoming the problem that fixed-neighborhood convolution cannot adapt to irregular ice block boundaries. The ice type prior guidance layer reads the historical ice type distribution prior probability map obtained from the analysis of historical sea ice in the Bohai Sea from 2015 to 2024. It generates context vectors for four ice types: gray ice, white ice, fixed ice, and broken ice. These vectors are then injected into each deformation map convolutional layer through conditionally normalized affine transformation parameters, enabling the network to employ differentiated feature propagation strategies for different ice type regions. The coarse-to-fine two-level graph iterative refinement module first completes the first round of graph inference at a coarse resolution of 50m, obtaining a coarse-scale sea ice distribution probability map and a preliminary estimate of sea ice concentration. Then, within the coverage area of ​​UAV multispectral imagery, it dynamically updates node features and reconstructs edge connections using UAV multispectral enhancement features, refining the graph resolution to the original resolution of the UAV imagery, completing the second round of graph inference. The decoder adopts a feature pyramid network structure, fusing multi-scale features step-by-step through top-down paths and lateral connections to restore the output resolution to 50m. The dual-head output layer outputs the sea ice distribution probability map and the preliminary estimate of sea ice concentration through Sigmoid activation and linear mapping, respectively.

[0059] The specific implementation of step S05 is as follows: The sparse Bayesian multi-source ice concentration joint inversion algorithm uses a sparse Bayesian learning framework to probabilistically model the mapping relationship between multi-source observations and sea ice concentration. The algorithm expands the sea ice concentration field into a linear superposition of an overcomplete dictionary composed of historical typical ice condition samples. Each basis function in the overcomplete dictionary corresponds to a historical typical ice condition pattern. Optical reflectivity, radar backscattering, and UAV multispectral index are used as multi-view measurement vectors, respectively. The hyperparameters of each basis function are iteratively estimated using the principle of maximizing evidence. The magnitude of the hyperparameters automatically controls the weight of the corresponding basis function in the final solution. When the hyperparameters approach infinity, the weight of the corresponding basis function approaches zero, thereby achieving automatic sparsification and retaining only the basis functions most relevant to the current ice field. Each iteration simultaneously updates the posterior mean of the sea ice concentration field and the posterior estimate of the variance of the noise of each source observation, enabling the algorithm to adaptively adjust the contribution weights of each source in different observation quality regions. After the iteration converges, the posterior expectation map of sea ice concentration is output as the final concentration inversion result. At the same time, the posterior standard deviation map is output as the confidence interval map of sea ice concentration. The larger the posterior standard deviation, the lower the prediction confidence of the corresponding pixel, which provides a quantitative basis for subsequent quality control.

[0060] The specific implementation of step S06 is as follows: Apply a distribution threshold of 0.5 to the sea ice distribution probability map output in step S04. Pixels with a probability value greater than 0.5 are determined to be ice-covered, and a sea ice binary mask is generated. Apply a 3×3 window mid-range filter to the sea ice concentration posterior expectation map output in step S05 to eliminate isolated salt-and-pepper noise points while preserving the spatial sharpness of the ice edge line. Perform a logical consistency check and force the density values ​​corresponding to pixels determined to be non-ice areas in the binary mask to zero to ensure the physical consistency between the distribution map and the density map. Based on the shoreline vector extracted from the land-water boundary of the Haiyang-1 C / D satellite optical image, crop the sea ice distribution raster product and the sea ice concentration raster product to remove land pixels, and finally output the sea ice distribution raster product and the sea ice concentration raster product facing the nearshore area of ​​the Bohai Sea.

[0061] It should be noted that the key technologies of this invention include the synergistic effect of three technologies: ice drift physical constraint iterative optimal transmission spatiotemporal registration technology, ice type prior guidance multi-scale deformation graph convolutional fusion technology, and sparse Bayesian multi-source ice density joint inversion technology. The optimal transmission registration technology transforms multi-source data registration into a global probability allocation problem. Physical constraints enable adaptive switching of alignment strategies between nearshore fixed ice and drifting ice regions, providing geometrically consistent input for subsequent fusion and overcoming the failure problem of traditional feature point matching in low-texture scenarios with broken ice. The graph convolutional fusion technology models the spatial topology of sea ice with a graph structure. Deformation convolution dynamically changes the receptive field with the shape of the ice block, and ice type prior guidance embeds domain knowledge into feature propagation, establishing a nonlinear mapping between multi-source observations and sea ice parameters. The sparse Bayesian inversion technology, through an automatic correlation determination mechanism, adaptively allocates the weights of each source and quantifies the inversion uncertainty under degraded observation conditions. The three technologies are progressively enhanced from three levels: geometric alignment, nonlinear fusion, and uncertainty quantification. They work together to ensure the physical rationality and accuracy consistency of the fusion results in complex nearshore scenarios, and have significant advantages over the traditional weighted average method in both ice-bearing and cloud-covered areas.

[0062] It should be noted that this invention also solves the following technical problem: In the high suspended sediment environment of Liaodong Bay in the Bohai Sea, the traditional normalized difference ice and snow index suffers from a high rate of missed detection of thin ice due to the overlap of the spectral responses of sediment and thin ice in the green band and shortwave infrared band. This invention introduces the near-infrared brightness ratio index, which utilizes the strong absorption characteristics of suspended sediment in the 865nm band and the relatively high reflectivity difference of thin ice in this band to construct a spectral index with stronger distinguishing ability between turbid water and thin ice, thus overcoming the limitations of the normalized difference ice and snow index in the special spectral environment of Liaodong Bay. The classification threshold of the near-infrared brightness ratio index is obtained by plotting the distribution histogram of the synchronous spectral measurement samples of thin ice, turbid water, and open water, and confirming the index value at the point of least overlap between the two distributions through 5-fold cross-validation, ensuring the stable applicability of the threshold under the typical optical conditions of the Bohai Sea, thereby effectively reducing the problem of underestimated sea ice area caused by missed detection of thin ice in the high suspended sediment environment.

[0063] A second aspect of the present invention provides a computer-readable storage medium storing program instructions, which, when executed in a computer, are used to perform the above-described method for multi-source remote sensing spatiotemporal fusion of sea ice parameters for the Bohai Sea nearshore area.

[0064] A third aspect of the present invention provides a multi-source remote sensing spatiotemporal fusion system for sea ice parameters in the Bohai Sea nearshore area, comprising the aforementioned computer-readable storage medium. The system can be any one of a computer, a server, or a microcontroller. The computer-readable storage medium is disposed within the system, and the system is provided with a microprocessor that executes the program instructions stored in the computer-readable storage medium.

[0065] Specifically, the principle of this invention is as follows: The reason why this invention can solve the above-mentioned technical problems lies in the following points. First, the ice drift physical constraint iterative optimal transmission spatiotemporal registration algorithm models multi-temporal ice field registration as a metric space mapping problem with physical constraints, introduces ice dynamic constraints and anisotropic diffusion constraints, and enables the alignment strategy of nearshore fixed ice and drifting ice regions to adaptively switch, thus ensuring the geometric consistency of multi-source input data from the root. Second, the ice type prior-guided multi-scale deformation graph convolutional fusion model uses graph reasoning as a framework, fuses spatial adjacency and ice dynamic correlation between node edge weights, and the learnable offset field of deformation graph convolution makes the receptive field adaptive to the ice block geometry. The ice type prior-guided layer embeds the historical ice type prior probability graph into feature propagation through conditional normalization, enabling the network to establish differentiated nonlinear mapping relationships in different ice type regions. Finally, the sparse Bayesian learning framework iteratively estimates the hyperparameters of each basis function through evidence maximization, automatically driving the weights of unrelated basis functions to zero, achieving adaptive weight allocation of the observation quality of each source, and outputting statistically significant posterior expected values ​​and confidence intervals, thus enabling the fusion results to have interpretable uncertainty quantification capabilities. These three progressively advanced steps collectively ensure the physical rationality and nonlinear modeling capabilities of the method in complex nearshore scenarios.

[0066] The following provides a specific embodiment 1 of the present invention, and the specific implementation of each step in this embodiment 1 is described in detail below.

[0067] The specific implementation of step S01 is as follows: acquire optical images from the Haiyang-1 C / D satellite, multispectral images from UAVs, and backscatter maps from shore-based X-band radar. Perform radiometric calibration, atmospheric correction, cloud masking, and clutter suppression processing on the three sources of data respectively. Project them uniformly onto the WGS84UTMZone51N coordinate system and resample them to 50m pixels to obtain the three-source preprocessed feature tensors. Among them, clutter suppression of the shore-based X-band radar backscatter map includes ice wave separation processing. The texture time-varying characteristics of the time-series backscatter map are used for differentiation. The Doppler frequency shift range of the sea waves is 0.1 to 0.3 Hz, and the ice drift frequency is less than 0.01 Hz. Pixel-level Doppler spectral features are extracted by short-time Fourier transform and input into a one-dimensional lightweight convolution classifier for fine separation of ice, water, and wave mixtures to obtain a backscatter map with wave Bragg scattering interference removed.

[0068] The specific implementation of step S02 is as follows: using the transit time of the Haiyang-1 C / D satellite as the time anchor point, the time difference of the three-source data is controlled within 1 hour. An iterative optimal transmission spatiotemporal registration algorithm with ice drift physical constraints is used to complete the spatiotemporal alignment of the multi-source data, with a spatial registration error of less than 0.5 pixels. This algorithm models multi-temporal ice field registration as a metric space mapping problem with physical constraints, using optimal transmission theory as a framework. First, based on the time-series backscattering map of the shore-based X-band radar, the initial ice drift velocity field is estimated using the cross-correlation method. The estimation formula is expressed as follows:

[0069] ;

[0070] In the formula, For pixels The initial ice drift velocity at the location, in pixels per second, can be converted by multiplying by the pixel size of 50m. , These are the spatial coordinates of the current pixel, in pixels. These are the cell coordinates within the search window, in pixels. for Backscattering plot at time 1 The dimensionless grayscale value at that location for The dimensionless grayscale value of the backscattering image at the offset position at time step. The spatial offset that maximizes the cross-correlation value, expressed in pixels. The time interval between adjacent frames, in seconds. For reference only, the empirical value is 1 second. This represents the offset that maximizes the cross-correlation function. The transport cost matrix is ​​constructed using the initial ice drift velocity field, and the cost function is expressed as follows:

[0071] ;

[0072] In the formula, A pixel with dimensionless to pixel The transmission cost, , Each pixel , The spatial coordinate vector, in meters. for and The Euclidean distance between them, in meters. The reference distance normalization factor is empirically set at 50m. For pixels The ice drift velocity vector at the location, in units of , For pixels The divergence of the ice drift velocity field, in units of , This is the absolute value of the divergence, in units of . , As a reference divergence normalization factor, the empirical value is , Let be the weighting coefficient for ice dynamics constraints, dimensionless, with an empirical value of 0.5. The Sinkhorn-Knoop iterative algorithm is used to solve for the regularized optimal transport map. The objective function of the regularized optimal transport problem is expressed as follows:

[0073] ;

[0074] In the formula, The elements of the dimensionless transport mapping matrix represent pixels. to pixel The quality distribution ratio, The normalization factor is a one-dimensional normalizer with an empirical value of 1. is the entropy regularization coefficient, dimensionless, with an empirical value of 0.1. The Shannon entropy of the transmission mapping matrix, dimensionless. This is used to smooth the transport mapping. The Sinkhorn-Knoop iterative algorithm works by alternately applying... To solve the above problem, normalize the rows and columns of the given information. The iteration update formula is expressed as follows:

[0075] ;

[0076] In the formula, For the first The elements of the transport mapping matrix in the round of iterations are of dimensionless quality. The elements of the kernel matrix are dimensionless. and The first The row scaling factor and column scaling factor in each iteration are both dimensionless and are obtained through alternating updates. The row scaling factor update formula is: The formula for updating the column scaling factor is: initial value The default setting is an all-one vector. Each iteration updates the pixel correspondence and reverses the ice drift velocity field estimation. The convergence criterion formula is as follows:

[0077] ;

[0078] In the formula, The Frobenius norm of the difference between two adjacent rounds of transport mapping matrices, dimensionless. The total number of pixels, dimensionless. To achieve a convergence threshold with a dimensionless value, it was determined by setting transmission mapping changes of 0.01, 0.05, and 0.1 pixels for 100 sets of known true ice field sequences, and then statistically analyzing the trade-off between registration error and the number of iterations. Anisotropic diffusion constraints were introduced, and the ice drift velocity in the nearshore fixed ice region was also considered. Forced zeroing, open water allows for large displacement transmission, and finally outputs a multi-source data spatiotemporal aligned transform field and transforms the three-source preprocessed feature tensor into a three-source aligned feature tensor.

[0079] The specific implementation of step S03 is as follows: The normalized difference ice and snow index, near-infrared brightness ratio index, gray-level co-occurrence matrix texture features, and Gabor multi-directional filtering response are calculated on the three-source aligned feature tensor to obtain the three-source enhanced feature tensor. The calculation formula for the normalized difference ice and snow index is expressed as follows:

[0080] ;

[0081] In the formula, The normalized difference ice and snow index is a dimensionless value. The green band reflectance is dimensionless. Here, the reflectance is in the short-wave infrared band, dimensionless. The formula for calculating the near-infrared luminance ratio index is as follows:

[0082] ;

[0083] In the formula, The near-infrared luminance ratio index is a dimensionless quantity. Near-infrared reflectance, dimensionless. The reflectivity is in the red band, dimensionless. The reflectance is in the green band, dimensionless. By utilizing the difference between the strong absorption characteristics of suspended sediment in the 865nm wavelength band and the high reflectivity of thin ice in the 865nm wavelength band, a method is constructed to distinguish between turbid water and thin ice, thus compensating for the shortcomings of traditional methods. The problem of missed detection of thin ice due to spectral overlap in the high suspended sediment environment of Liaodong Bay. The classification threshold was obtained through the following experiment: Ten sets of on-site synchronous spectral measurement samples were collected in Liaodong Bay during winter, and measurements were taken on thin ice, turbid water, and open water. Values, plot distribution histograms, and take the values ​​where the overlap between the two distributions is minimal. The value was used as a classification threshold, and its stability was confirmed through 5-fold cross-validation. Based on the cloud coverage metric... ice drift amplitude index ,when Skip low-light enhancement processing when Low-light enhancement processing based on a radiative transfer model is performed on the optical data portion. Among other things, Calculated from the pixel ratio of the optical image cloud mask, dimensionless Estimated by cross-correlation of backscatter maps at adjacent time points from shore-based X-band radar, in units of The cloud coverage threshold of 0.6 was obtained through statistical analysis of historical satellite archive data from 10 winters in the Bohai Sea from 2015 to 2024, and was determined by taking the inflection point value at which sea ice identification accuracy significantly decreased.

[0084] The specific implementation of step S04 is as follows: The three-source enhanced feature tensor is input into the ice-type prior-guided multi-scale deformation map convolutional fusion model. The input layer spatial size is 512×512 pixels. The superpixel map construction module divides the study area into several superpixel nodes with a coarse resolution of 50m, and the node feature vectors... It is composed of multi-scale convolutional encodings of three-source enhanced feature tensors at corresponding spatial locations, and has one dimension. The edge weights between nodes... Defined jointly by spatial adjacency and ice dynamic correlation, the ice dynamic correlation is quantified by the ice drift velocity field estimated by shore-based X-band radar. The calculation formula is expressed as follows:

[0085] ;

[0086] In the formula, For nodes With nodes The boundary weights between them are dimensionless. , They are nodes , The spatial coordinate vector, in meters. for and The Euclidean distance between them, in meters. This is a spatial distance attenuation scale, measured in meters (m), with an empirical value of 150 m. This is the ice dynamics correlation weighting coefficient, dimensionless, with an empirical value of 0.5. For nodes With nodes The angle between the direction of the line connecting the two points and the direction of the ice drift velocity, in rad. For nodes The corresponding normalization factor, dimensionless, For nodes The set of neighboring nodes. The deformation map convolution module is composed of a learnable offset field at each superpixel node. Dynamically adjust the sampling positions of neighboring nodes. The dimensionless two-dimensional offset, in pixels, is predicted from local features and optimized end-to-end by a lightweight convolutional network. The deformation map convolution aggregation formula is expressed as follows:

[0087] ;

[0088] In the formula, For the first Layer The feature vector of each node, with dimensionless... For graph convolutional layer index, For the first Layer learnable weight matrix, dimensionless. For the first Layer nodes Offset The feature vector obtained by post-sampling has one dimension. is the activation function for the linear rectifier unit. The ice type prior guidance layer utilizes the prior probability map of historical ice type distribution, and injects context vectors of four ice types—gray ice, white ice, fixed ice, and broken ice—into each deformation map convolutional layer through conditional normalization. The conditional normalization formula is expressed as follows:

[0089] ;

[0090] In the formula, For the first The node The eigenvalues ​​after channel condition normalization are dimensionless. For the first The node The eigenvalues ​​before channel normalization have one dimension. For the first The mean of the channel within the batch, dimensionless. For the first The standard deviation of the channel within the batch, dimensionless. This is a numerically stable term, with one dimension, and is assumed to be [value missing]. , and These are the scaling and offset parameters generated from the prior probability map of historical ice type distribution through a fully connected layer, respectively, both of dimensionless. This is for feature channel indexing. The coarse-to-fine two-level graph iterative refinement module completes the first round of graph inference at a coarse resolution of 50m. Then, it dynamically updates node features and reconstructs edge connections using UAV multispectral enhancement features, refining the graph resolution to the original UAV resolution within the UAV image coverage area, completing the second round of graph inference. The decoder uses a feature pyramid network structure for feature upsampling and multi-scale fusion, restoring the output resolution to 50m. The output layer sea ice distribution head outputs a single-channel sea ice distribution probability map. Activate, and the concentration regression head outputs a single-channel preliminary sea ice concentration map, which is then linearly mapped to 0%–100%.

[0091] The specific implementation of step S05 is as follows: Based on the sparse Bayesian multi-source ice concentration joint inversion algorithm, using the initial sea ice concentration estimate map as the initial value, the algorithm integrates optical reflectivity, radar backscattering, and UAV multispectral index from the three-source aligned feature tensor, and iteratively outputs the posterior expectation map and the sea ice concentration confidence interval map. This algorithm expands the sea ice concentration field into a linear superposition of overcomplete dictionaries. The multi-source observation model formula is expressed as follows:

[0092] ;

[0093] In the formula, For the first The source's observation vector, dimensionless. The three categories are optical, radar, and UAV multispectral. For the first The source corresponds to an overcomplete dictionary matrix, with one dimension, composed of historical typical ice condition samples. Let be the sea ice concentration coefficient vector, dimensionless. For the first The observed noise vector of the source, dimensionless, follows a pattern with zero mean and variance. Gaussian distribution, Let be the observation noise variance with dimensionless variance. The hyperparameters of each basis function are estimated iteratively using evidence maximization. Automatic sparsity is achieved by automatically driving the weights of unrelated basis functions to approach zero, and the hyperparameters... The update formula is expressed as follows:

[0094] ;

[0095] In the formula, For the first The update values ​​of the hyperparameters of the basis functions, with dimensions one. For the current number The hyperparameter values ​​of the basis functions are dimensionless. The posterior covariance matrix is ​​the first... One diagonal element, dimensionless. A vector of sea ice concentration coefficients No. The posterior mean of each component, dimensionless. This is a reference value for sea ice concentration, dimensionless, and taken as 1 (corresponding to 100% concentration). Observation noise variance. The update formula is expressed as follows:

[0096] ;

[0097] In the formula, For the first The square norm of the source observation residuals, dimensionless. For the first Reference value for the variance of source observation noise, dimensionless, default value is 1. A vector of sea ice concentration coefficients The posterior mean vector, with one dimension. For the first Source observation vector The dimension is one, and the summation term in the denominator is the effective degrees of freedom correction term, also with one dimension. Each iteration simultaneously updates... and The posterior estimate is used to finally output the posterior expected value map of sea ice concentration. The confidence interval plot for sea ice concentration is calculated using the posterior standard deviation. Characterization, dimensionless The larger the value, the lower the prediction confidence of the corresponding pixel.

[0098] The specific implementation of step S06 is as follows: apply a threshold of 0.5 to the sea ice distribution probability map to generate a sea ice binary mask, apply a 3×3 window mid-range filter to the sea ice concentration posterior expectation map, perform logical consistency verification, crop the land area based on the shoreline vector, the shoreline vector is extracted from the land-water boundary of the Haiyang-1 C / D satellite optical image, and finally output the sea ice distribution raster product and the sea ice concentration raster product.

[0099] The specific implementation of the dynamic loss weighting function is as follows: Let the Dice loss weighting coefficient of the sea ice distribution head be... The weighted coefficient of the L1 loss for the head smoothing of the density regression is: Both satisfy Each value is limited to a range of 0.3 to 0.7. (This refers to the fusion imbalance index.) The calculation formula is expressed as follows:

[0100] ;

[0101] In the formula, The fusion imbalance index is a dimensionless quantity. The Dice loss value for the current sea ice distribution, dimensionless. The initial baseline value for the Dice loss of sea ice distribution head at the end of the first round of training, dimensionless. This represents the L1 loss value of the head smoothing regression for the current round, with one dimension. The initial baseline value for the dense regression head smoothed L1 loss at the end of the first round of training is given, with a dimension of one. When hour, Increase by 0.05. The corresponding reduction is 0.05; when hour, and Remain unchanged; when hour, Increase by 0.05. The corresponding reduction is 0.05. The fusion imbalance index thresholds of 1.5 and 0.67 were obtained by analyzing 10 training experiments with different initial weighting coefficient configurations. Specifically, the convergence speed and final accuracy of the validation set were statistically analyzed under each configuration, and the boundary of the fusion imbalance index interval that makes the convergence speed of the two tasks most balanced was taken as the threshold.

[0102] To better understand and implement this invention, the following is a specific application scenario of the invention, Example 2: To illustrate the effect of the invention, technicians set up a test environment and, by acquiring typical winter ice condition observation data of a certain area near the coast of Liaodong Bay in the Bohai Sea, comprehensively applied the multi-source remote sensing spatiotemporal fusion method proposed in this invention to carry out a complete inversion experiment on the distribution range and density of sea ice.

[0103] This experiment uses a typical ice-affected area near the Liaodong Bay in the Bohai Sea on January 4, 2026 as the research object, selecting the transit time of the Haiyang-1C satellite on that day as the time anchor point. The acquisition of the three data sources is as follows: After radiometric calibration and atmospheric correction, the effective observation coverage of the Haiyang-1C satellite optical image is approximately 74%, with a cloud coverage index of 0.26, below the threshold of 0.6. Therefore, low-light enhancement processing based on the radiative transfer model was performed on the optical data. The UAV multispectral image covers key nearshore areas with a spatial resolution of 0.4m. After orthorectification, it was resampled to 50m to match the satellite image. The shore-based X-band radar provides a continuous backscattering intensity sequence, which, after ice wave separation processing, removes Bragg scattering interference from sea waves to generate an effective backscattering map. The time difference between the three data sources is controlled within 1 hour. After processing using the ice drift physical constraint iterative optimal transmission spatiotemporal registration algorithm, the spatial registration error verification result is better than 0.5 pixels, meeting the geometric consistency requirements for subsequent fusion.

[0104] In the feature enhancement step, the normalized difference ice and snow index distribution ranged from -0.12 to 0.78, and the near-infrared brightness ratio index distribution ranged from 0.31 to 2.15 in the study area. The study area contains typical water bodies with high suspended sediment. When using only the normalized difference ice and snow index for discrimination, thin ice pixels showed significant missed detections. Introducing the near-infrared brightness ratio index significantly improved the ability to distinguish between thin ice and turbid water. Gray-level co-occurrence matrix texture features and Gabor multi-directional filtering responses further enhanced the discriminability of ice fragment boundaries. The concatenated features formed a three-source enhanced feature tensor, which served as the input to the graph convolutional fusion model.

[0105] In the inference stage of the multi-scale deformation graph convolutional fusion model guided by ice type priors, the study area was divided into several superpixel nodes. The learnable offset field of the deformation graph convolution exhibited a significant irregular sampling pattern in areas with dense ice fragmentation, and the receptive field shape dynamically changed with the degree of ice fragmentation. Historical ice type distribution prior probability maps showed that the nearshore area of ​​the study area was dominated by fixed ice and broken ice. Based on this, the ice type prior guidance layer injected differentiated contextual vectors into the corresponding areas, making graph inference tend to smoothly aggregate large receptive fields in fixed ice areas and finely capture small receptive fields in broken ice areas. After two levels of coarse and fine iterative refinement, the details of sea ice boundaries in the UAV multispectral coverage area were significantly improved. The overall spatial distribution of the sea ice distribution probability map and the preliminary sea ice concentration estimate output by this inference is as follows: Figure 2 As shown.

[0106] The sparse Bayesian multi-source ice concentration joint inversion step uses the aforementioned preliminary estimate as the initial value, integrates optical reflectivity, radar backscattering, and UAV multispectral index from the three-source aligned feature tensor, and outputs a posterior expectation map and a posterior standard deviation map of sea ice concentration after iterative evidence maximization convergence. The posterior standard deviation map shows that the posterior standard deviation is generally low in the nearshore fixed ice region and high in the ice-water mixing transition zone, consistent with physical expectations, reflecting the adaptive weight allocation capability of the sparse Bayesian framework for the observation quality of different regions.

[0107] The post-processing steps involve applying a distribution threshold of 0.5 to the sea ice distribution probability map to generate a binary mask, applying a 3×3 median filter to the concentration posterior expectation map, performing a logical consistency check, and cropping the land area based on the shoreline vector. This yields the final sea ice distribution raster product and sea ice concentration raster product for the study area. To quantitatively evaluate the inversion accuracy, the fused results are compared with synchronous field observation data; the statistical results are shown in Table 1.

[0108] Table 1. Evaluation Results of Sea Ice Concentration Inversion Accuracy

[0109]

[0110] As shown in Table 1, this method outperforms the traditional weighted average method in both sea ice distribution identification and concentration inversion, especially in the thin ice detection rate, demonstrating the synergistic effect of near-infrared brightness ratio index and graph convolutional nonlinear modeling.

[0111] This experiment also conducted a spatiotemporal plausibility check on continuous time-series sea ice products from January 5th to January 10th, 2026, such as... Figure 3 As shown, the ice edge line is continuous and smooth, without obvious jagged edges or isolated patches. The ice condition evolution in the time series conforms to physical laws and there are no abnormal jumps. This indicates that the physical constraints of ice drift spatiotemporal registration and the sparse Bayesian joint inversion jointly ensure the temporal consistency of the product.

[0112] like Figure 4 As shown, the spatial distribution of the contribution weights of each source observation indicates that in cloud-covered areas, the weight of optical data is automatically reduced to near zero, while the weight of shore-based radar is correspondingly increased; in open water areas with low radar backscattering contrast, the weight of optical data is relatively high; the weight contribution of UAV multispectral data to the thin ice boundary is most prominent in the coverage area, which is consistent with the physical characteristics of each sensor.

[0113] The technological advancements of this invention compared to traditional methods are reflected in the following aspects. First, traditional linear weighted fusion methods treat each source data as an equal or fixed weight combination, failing to dynamically adjust the credibility of each source based on local scenarios. This invention, however, utilizes a sparse Bayesian automatic correlation determination mechanism to achieve adaptive allocation of source weights from a probabilistic perspective, outputting physically interpretable results even under degraded observation conditions. Second, traditional methods cannot model the nonlinear response relationship between sea ice concentration and multi-source observations. This invention's ice-type prior-guided multi-scale deformation graph convolutional fusion model establishes a nonlinear mapping using a graph inference framework, allowing the receptive field to dynamically change with ice block shape. Prior knowledge is embedded in feature propagation, fundamentally breaking through the expression limits of linear methods. Furthermore, traditional feature point matching methods fail in low-texture scenarios with broken ice due to feature sparsity. This invention's optimal transmission registration framework transforms registration into a global probability allocation problem, with physical constraints ensuring sub-pixel-level geometric consistency in complex nearshore scenarios, laying a reliable foundation for subsequent fusion.

[0114] It should be noted that the variables involved in this invention are explained in detail in Tables 2 and 3.

[0115] Table 2. Variable Explanation Table (Part 1)

[0116]

[0117] Table 3. Variable Explanation Table (Part Two)

[0118]

[0119] The above description is merely a specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any changes or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in the present invention should be included within the scope of protection of the present invention.

Claims

1. A multi-source remote sensing spatiotemporal fusion method for sea ice parameters in the Bohai Sea nearshore area, characterized in that, Includes the following steps: We acquired optical images from the Haiyang-1 C / D satellite, multispectral images from UAVs, and backscattering maps from shore-based X-band radar. We performed radiometric calibration, atmospheric correction, cloud masking, and clutter suppression on the three-source data, then uniformly projected and resampled them to obtain the three-source preprocessed feature tensors. Using the transit time of Haiyang-1C / D satellite as the time anchor point, and based on the ice drift physical constraint iterative optimal transmission spatiotemporal registration algorithm, the spatiotemporal alignment transformation field of multi-source data is calculated for the three-source preprocessed feature tensor, and the three-source preprocessed feature tensor is transformed into a three-source aligned feature tensor. The normalized difference ice and snow index, near-infrared brightness ratio index, gray-level co-occurrence matrix texture features, and Gabor multi-directional filtering response are calculated for the three-source aligned feature tensor to obtain the three-source enhanced feature tensor. Based on the thresholds of cloud coverage index and ice drift amplitude index, low-light enhancement processing is performed on the optical data part or low-light enhancement processing is skipped. The three-source enhanced feature tensor is input into the ice type prior-guided multi-scale deformation map convolutional fusion model. After superpixel map construction, deformation map convolutional aggregation and coarse-fine two-level map iterative refinement, the sea ice distribution probability map and the initial sea ice concentration map are output simultaneously. Based on the sparse Bayesian multi-source ice concentration joint inversion algorithm, the initial sea ice concentration map is used as the initial value. The three-source aligned feature tensor is fused to iteratively output the posterior expectation map of sea ice concentration and the confidence interval map of sea ice concentration. A distribution threshold is applied to the sea ice distribution probability map to generate a binary sea ice mask. Median filtering is applied to the posterior expectation map of sea ice concentration. Logical consistency check is performed. Land areas are clipped based on shoreline vectors. Finally, sea ice distribution raster products and sea ice concentration raster products are output.

2. The multi-source remote sensing spatiotemporal fusion method for sea ice parameters in the Bohai Sea nearshore area according to claim 1, characterized in that, The unified projection and resampling specifically involves projecting onto the WGS84UTMZone51N coordinate system and resampling to 50m pixels.

3. The multi-source remote sensing spatiotemporal fusion method for sea ice parameters in the Bohai Sea nearshore area according to claim 2, characterized in that, The sparse Bayesian multi-source ice concentration joint inversion algorithm specifically unfolds the sea ice concentration field into a linear superposition of an overcomplete dictionary composed of historical typical ice condition samples, uses evidence maximization to iteratively estimate the hyperparameters of each basis function, automatically drives the weights of unrelated basis functions to approach zero, and updates the posterior estimates of the sea ice concentration field and the variance of the observation noise of each source at each iteration.

4. The multi-source remote sensing spatiotemporal fusion method for sea ice parameters in the Bohai Sea nearshore area according to claim 3, characterized in that, The ice drift physical constraint iterative optimal transport spatiotemporal registration algorithm is specifically based on the optimal transport theory framework. It uses the cross-correlation method to estimate the initial ice drift velocity field, constructs a transport cost matrix that integrates Euclidean distance and flow velocity divergence minimization constraints, uses the Sinkhorn-Knopp iterative algorithm to solve the regularized optimal transport mapping, introduces anisotropic diffusion constraints, and stops iterating when the change in the transport mapping is lower than the convergence threshold.

5. The multi-source remote sensing spatiotemporal fusion method for sea ice parameters in the Bohai Sea nearshore area according to claim 4, characterized in that, The clutter suppression process includes ice wave separation processing, which specifically involves extracting pixel-level Doppler spectral features through short-time Fourier transform, using the wave Doppler frequency shift range and the upper limit of ice drift frequency as separation thresholds, and inputting them into a one-dimensional lightweight convolution classifier for fine separation of ice, water, and wave mixtures.

6. The multi-source remote sensing spatiotemporal fusion method for sea ice parameters in the Bohai Sea nearshore area according to claim 5, characterized in that, The formula for calculating the near-infrared luminance ratio index is as follows: ,in For near-infrared reflectivity, For red band reflectivity, This refers to the reflectivity in the green band.

7. The multi-source remote sensing spatiotemporal fusion method for sea ice parameters in the Bohai Sea nearshore area according to claim 6, characterized in that, The cloud coverage rate index is calculated from the pixel ratio of the cloud mask in the optical image. When the cloud coverage rate index is greater than the cloud coverage rate threshold, the low-light enhancement process is skipped. When the cloud coverage rate index is not greater than the cloud coverage rate threshold, the low-light enhancement process based on the radiative transfer model is performed. The cloud coverage rate threshold is 0.

6.

8. The multi-source remote sensing spatiotemporal fusion method for sea ice parameters in the Bohai Sea nearshore area according to claim 7, characterized in that, The input layer of the ice-type prior-guided multi-scale deformation map convolutional fusion model receives a three-source enhanced feature tensor with a spatial size of 512×512 pixels. The superpixel map construction module divides superpixel nodes at a coarse resolution of 50m, and the edge weights between nodes are defined by spatial adjacency and ice dynamics correlation.

9. A computer-readable storage medium, characterized in that, The computer-readable storage medium stores program instructions, which, when executed in a computer, are used to perform the spatiotemporal fusion method for multi-source remote sensing of sea ice parameters for the Bohai Sea nearshore area as described in any one of claims 1-8.

10. A multi-source remote sensing spatiotemporal fusion system for sea ice parameters in the Bohai Sea nearshore area, characterized in that, The system comprises the computer-readable storage medium of claim 9, wherein the system is a computer, the computer-readable storage medium is disposed within the system, and the system is provided with a microprocessor that executes program instructions stored in the computer-readable storage medium.

Images