Sparse ais ocean current field reconstruction method based on physical information neural network

Abstract

The present application belongs to the field of artificial intelligence and marine environment observation technology, and relates to a sparse AIS sea current field reconstruction method based on a physical information neural network. The method constructs a PI-NODE model, takes sparse AIS inversion velocity as input, and introduces CMEMS background priori and marine physical soft constraint in the training stage to realize spatiotemporal continuous sea current field reconstruction. The model includes a confidence weighted observation model, a multi-scale parallel encoder, a gated data-bounded guided double branch neural ODE, a flow function auxiliary decoder and a physical soft constraint system. The present application provides a sea surface current field reconstruction method which cooperatively fuses sparse observation, continuous background priori and physical constraint. The flow function auxiliary decoder is introduced to fuse the prediction result of the direct velocity branch and the auxiliary velocity field derived from the flow function according to a learnable weight, so as to reduce unreasonable structures in the output flow field.

Description

Technical Field

[0001] This invention belongs to the field of artificial intelligence and marine environmental observation technology, and relates to a sparse AIS ocean current field reconstruction method based on physical information neural network. Background Technology

[0002] High-resolution sea surface currents influence ship drift, pollutant diffusion, material transport, and ocean dynamics. Therefore, obtaining continuous current fields with high spatial resolution is a crucial task in ocean observation. However, due to limitations in the physical characteristics and engineering conditions of observation technologies, existing methods struggle to simultaneously achieve both coverage and observational accuracy: drifting buoys provide reliable point velocities, but their spatial distribution is limited; high-frequency ground-wave radar can continuously observe nearshore surface current fields, but its coverage is limited by shore-based stations; satellite remote sensing offers wide-area observation capabilities, but optical remote sensing is susceptible to cloud conditions and lighting conditions; synthetic aperture radar typically only acquires the velocity component along the radar line of sight, and the results are affected by wind, waves, non-geophysical Doppler signals, and revisit periods.

[0003] To address the limitations of existing observations in terms of coverage, the Automatic Identification System (AIS) provides an alternative source of sea surface current field estimation data. Originally used for ship identification and navigation safety, AIS continuously records information such as ship position, speed over land, heading over land, bow direction, and time. Since ship motion over land is influenced by ocean currents, these records can also be used to estimate ocean current speeds. This invention refers to the process of estimating ocean current speeds using AIS ship motion information as AIS current inversion, and the resulting discrete two-dimensional velocity field, composed of east-west and north-south velocities arranged according to spatial location, is simply called the AIS inversion field.

[0004] A series of methodological explorations have been conducted both domestically and internationally regarding the AIS current inversion problem. However, due to the objective limitations of ship route distribution, the AIS inverted current field still suffers from significant spatial sparsity and large-scale gaps in open sea areas. Daily inversion results are insufficient to cover the entire study area, and local anomalies are prone to appear in low-density areas due to insufficient course diversity. Overall, it exhibits a non-uniform distribution characteristic of "dense along shipping routes and sparse along distant routes." Therefore, how to further reconstruct a spatially continuous, temporally consistent, and physically self-consistent high-resolution sea surface current field based on the sparse current field inverted by AIS remains a key scientific problem that needs to be solved. Summary of the Invention

[0005] To address the challenge of reconstructing continuous sea surface current fields from sparse AIS inversion fields, this invention provides a method for reconstructing sparse AIS ocean current fields based on a Physical Information Neural Network (PI-NODE). This method utilizes the Physical Information Neural ODE (PI-NODE) ​​framework, combining AIS inversion observations, CMEMS background priors, and soft constraints from ocean physics during the training phase to reconstruct a spatiotemporally continuous, high-resolution sea surface current field and improve the physical plausibility of the reconstruction results.

[0006] The technical solution provided by this invention is as follows: A sparse AIS ocean current field reconstruction method based on physical information neural network, comprising the following steps: (1) Construct the PI-NODE model, which includes a multi-scale parallel encoder, a gated data-bounded guided dual-branch constant differential equation, and a stream function-assisted decoder; (2) The multi-scale parallel encoder extracts features from the input three-channel observation field, which includes the east-west and north-south velocity fields and the corresponding confidence fields; the outputs of each branch are aligned to the same spatial size, and feature fusion and latent variable mapping are performed in the channel dimension to obtain the initial latent state; (3) The gated data-bounded guided dual-branch ordinary differential equations take the initial latent state as the initial value, calculate the latent state change rate under the combined action of the data branch and the Tanh bounded guided branch, and obtain a discrete latent state sequence after solving; (4) Each latent state is restored to spatial features through a fully connected layer, and then enters the stream function-assisted decoder for decoding to obtain the sea surface current field prediction sequence within the observation window; (5) The prediction results of the overlapping windows are weighted and synthesized to reconstruct the prediction sequence and obtain the daily continuous sea surface flow field.

[0007] Preferably, the PI-NODE model further includes a confidence-weighted observation model, which constructs grid confidence based on the number of samples and the consistency of neighborhood flow direction, in order to adjust the AIS reconstruction error.

[0008] Preferably, the multi-scale parallel encoder includes three parallel convolutional branches, with the kernel scale of the first two convolutional layers respectively adopting... , and Features are extracted from different spatial ranges.

[0009] Preferably, the stream function-assisted decoder includes two parallel decoding branches: a direct speed branch and a stream function branch. The direct speed branch retains local expressive power, while the stream function branch provides an auxiliary structure with low divergence tendency. The two branches are fused through learnable global coefficients.

[0010] Preferably, during the training phase, the PI-NODE model compares the predicted velocity field for each day with the AIS inversion velocity for that day to calculate the reconstruction error.

[0011] Preferably, during the training phase, the PI-NODE model constructs a physical soft-constraint loss based on the changes and spatial structure of the predicted velocity fields over multiple consecutive days. ; in, This is the loss due to physical soft constraints; , , and These are the spatial-temporal consistency losses of divergence and vorticity, temporal continuity, and spatial Laplace smoothing, respectively. , , and These are the weighting coefficients for the corresponding loss terms.

[0012] Preferably, the total training loss of the PI-NODE model during the training phase is: ; Where L is the total loss; Huber loss weighted by confidence level; For latent variable regularization, For a physically soft-constrained system, For velocity modulus constraints, For CMEMS background priors; , and These are the fixed weights for latent variable regularization, velocity modulus constraint, and background prior term, respectively. To follow the training rounds The changing physical scaling factor is adjusted only. The intensity of the effect. The various loss terms and physical scaling methods are given in the specific implementation details.

[0013] Preferably, the training samples consist of a multi-day window; the normalized AIS velocity field and confidence field of the first day are fed into a multi-scale parallel encoder to generate an initial latent state; the gated data-bounded guided dual-branch constant differential equation evolves the latent state at multiple time steps, including the first day, starting from this initial latent state.

[0014] Compared with the prior art, the present invention has the following beneficial effects: (1) A PI-NODE model framework for sparse AIS observations is proposed: the AIS observation field is encoded by a variational autoencoder and the state evolution within the time window is described by a neural ordinary differential equation in the latent space; CMEMS is used as the background prior in the training phase and physical constraints are applied to the output flow field.

[0015] (2) Design a multi-scale parallel encoder and variational mapping module. The encoder uses 3×3, 5×5 and 7×7 convolutional branches to extract flow field features of different spatial ranges, completes feature fusion through channel attention, and maps the fusion result to the initial latent state required by the neural ODE.

[0016] (3) Construct a gated dual-branch neural ODE, calculate the latent state change rate by the data branch and the Tanh bounded guiding branch respectively, and then fuse the two by gating weights to complete the continuous evolution within the time window.

[0017] (4) Introduce a stream function-assisted decoder to fuse the prediction results of the direct velocity branch with the auxiliary velocity field derived from the stream function according to learnable weights, thereby reducing unreasonable structures in the output flow field.

[0018] (5) The model adjusts the reconstruction error between the predicted velocity and the AIS inversion velocity on the same day based on the grid confidence, and combines latent variable regularization, CMEMS background prior and constraints such as divergence, vortex evolution, time continuity, spatial smoothness and velocity modulus; among which the weight of physical constraints is gradually increased during the training process. Attached Figure Description

[0019] Figure 1 This is the overall framework of the PI-NODE model in the embodiments of the present invention; Figure 2 This is a schematic diagram of the processing flow of a multi-scale parallel encoder in an embodiment of the present invention; Figure 3 This is a schematic diagram of the processing flow of the stream function-assisted decoder in an embodiment of the present invention; Figure 4 This is a schematic diagram illustrating the relationship between the reconstruction sequence, each constraint term, and parameter updates in an embodiment of the present invention; Figure 5 For the same batch of buoy samples in this embodiment of the invention, PI-NODE and CMEMS are... , Component density scatter plot verification; where (a) and (b) are respectively PI-NODE... , Components; (c) and (d) are CMEMS components respectively. , Components; the dot color represents the kernel density in the scatter space, and the dashed line is the 1:1 reference line. Detailed Implementation

[0020] The present invention will now be described in further detail with reference to the accompanying drawings and specific embodiments. It should be understood that the specific embodiments described herein are merely illustrative of the invention and are not intended to limit the invention.

[0021] Example 1: AIS inversion observations are spatially sparse and unevenly distributed, and cannot be directly converted into stable, continuous sea surface current fields. To address this issue, this example provides a method for reconstructing sparse AIS ocean current fields based on a physical information neural network. For example... Figure 1 As shown, this method mainly includes three stages: PI-NODE model construction and training, reconstruction / inference, and application. The construction of the PI-NODE model includes building an observation confidence model, a multi-scale variational encoder, a gated data-bounded guided two-branch Neural ODE, a stream function-assisted decoder, and physical soft constraints. During the model training stage, a large-scale reference is supplemented with background priors from CMEMS (Copernicus Marine Environment Monitoring Service).

[0022] Figure 1 The training and reconstruction process of the PI-NODE model is demonstrated: each sample contains data from 8 consecutive days. The normalized AIS velocity field and confidence field of the first day are used to generate the initial latent state. The gated data-bounded guided dual-branch neural ODE uses this initial latent state as a starting point and evolves the latent state over 8 time steps, including the first day. The stream function-assisted decoder recovers the two-dimensional velocity field step by step. During training, the model compares the predicted velocity field for each day with the AIS inverted velocity for that day to calculate the reconstruction error. Simultaneously, based on the changes and spatial structure of the predicted velocity fields over 8 consecutive days, physical soft constraints are constructed, with CMEMS only used as background prior in the loss calculation. During reconstruction, the prediction results of overlapping windows are weighted and synthesized to obtain the daily continuous sea surface current field.

[0023] The components of the PI-NODE model are explained below.

[0024] 1. Confidence-weighted observation model The reliability of observations retrieved through AIS clustering inversion is not spatially uniform. There are typically more samples near the main shipping lanes and fewer samples in areas far from the lanes. Furthermore, when the current directions in the inversion results of adjacent grids are inconsistent, even with similar sample numbers, local inversion results are more prone to instability. Therefore, this invention constructs grid confidence based on the number of samples and the consistency of neighboring current directions, and uses this confidence to adjust for AIS reconstruction errors.

[0025] (1) Velocity normalization The east-west and north-south velocity components obtained from AIS inversion are mapped to a 1 / 12° regular grid through linear interpolation and then normalized. ; ; in, and The first Daily grid The east-west and north-south AIS inversion speed at the location; , These are the grid row and column indexes, respectively. and For the corresponding normalized velocity component; and To maintain a fixed normalization scale, all values ​​are taken as 3 m / s; This means that values ​​exceeding [-1, 1] will be truncated to the corresponding endpoint.

[0026] Normalization aligns the input velocity with the numerical range of the decoder's Tanh output, reducing the impact of outliers on training. The reconstruction result is then multiplied by the same normalization scale to recover the physical unit. The tilde above the velocity symbol below indicates normalization.

[0027] (2) Observation density factor The number of effective samples participating in AIS inversion within a grid reflects the richness of local observations. To reduce the impact of differences in sample numbers across different grids on the confidence level, an observation density factor is defined: ; in, For the first Daily grid The observation density factor at that location The effective number of samples for this grid participating in AIS inversion is represented by log, where log represents the natural logarithm. Adding 1 ensures that the calculation is still possible even when the sample number is 0; the logarithmic transformation is used to reduce the dominance of extremely high sample numbers on the confidence level.

[0028] (3) Neighborhood directional coherence factor To measure the consistency of flow direction between the central grid and the surrounding grids, the flow angle of each grid is first calculated: .

[0029] in, For the first Daily grid The flow direction angle at a point is expressed in radians; atan2 determines the quadrant of the flow direction based on the two velocity components and their signs.

[0030] Based on the flow direction difference between the central grid and its 3×3 neighborhood, the neighborhood directional coherence factor is defined as follows: ; in, For grid The 3×3 neighborhood centered on the grid includes 9 locations, including the central grid. This represents the number of neighboring grid cells. , Used as the row and column index for the neighborhood grid. This is the cosine mean of the difference between the neighborhood flow direction and the central flow direction. The closer the value is to 1, the more consistent the directions; a value close to 0 indicates that the cosines of the direction differences generally cancel each other out; a negative value indicates that terms with a direction difference greater than 90° dominate the cosine summation. Calculations are only performed on internal grids with complete 3×3 neighborhoods; the boundary grid is set to 0.

[0031] (4) Confidence fusion The observation density factor and directional coherence factor are combined linearly and mapped using a sigmoid method to form the final confidence score. ; ; in, For grid confidence; and These are fixed weights for observation density and directional coherence, respectively. In this embodiment, we take... , Sigmoid takes real numbers as input. Mapped to (0,1), exp represents the natural exponential function.

[0032] (5) Confidence-weighted reconstruction loss The confidence level is used to adjust the weight of the reconstruction error. The residual between the normalized prediction field and the AIS target field is defined as: ; in, To reconstruct the residuals, For the model's predicted field, For AIS target fields, both contain normalized east-west and north-south velocity components; the cap indicates the prediction. For effective sea area masking, values ​​are 1 within the sea area and 0 elsewhere; ⊙ indicates element-wise multiplication. This indicates that the residual outside the sea area is set to zero.

[0033] To mitigate the impact of abnormal residuals on training, Huber loss is employed: ; ; in, To reconstruct residuals The corresponding Huber loss is δ, which is the transition parameter between the quadratic segment and the linear segment. In this embodiment, δ = 0.05.

[0034] The confidence-weighted reconstruction loss is: ; in, This indicates that the arithmetic mean of all confidence-weighted loss values ​​obtained within the parentheses is used as the reconstruction loss. The weight 0.3 + 0.7C ensures that high-confidence observations contribute more, while retaining the basic weight of low-confidence observations.

[0035] 2. Multi-scale parallel encoder The flow field retrieved by AIS simultaneously includes local velocity variations, tortuous flow bands, and large-scale background structures. To avoid bias towards a particular spatial pattern by a single convolutional scale, the multi-scale parallel encoder of this invention employs three parallel convolutional branches of 3×3, 5×5, and 7×7 to sequentially complete feature extraction, fusion, channel attention recalibration, and latent variable mapping.

[0036] (1) Multi-scale parallel feature extraction For an 8-day observation window, the encoder only receives the three-channel observation field from the first day: ; in, For the first day's input, Indicates the first day; , and These are the normalized east-west velocity field, north-south velocity field, and confidence field for that day; , Let represent the number of rows and columns of the target grid, and ℝ represent the real number field.

[0037] The encoder was used only for the three-channel observation field on the first day of the window. As input, the remaining 7 days of AIS observations within the window are not input to the encoder, but are compared with the model predictions at the corresponding time steps during the training phase to calculate the reconstruction error. The complete 8-day prediction sequence generated by the model is further used to calculate the physical soft constraints.

[0038] The encoder is configured with three parallel convolutional branches. All three branches receive... However, the kernel scale of the first two convolutional layers adopts different methods. , and Features are extracted from different spatial ranges. The outputs of the three branches are denoted as follows: .

[0039] in, , and These represent feature extraction branches using 3×3, 5×5, and 7×7 first two convolutional kernels, respectively. , and For the corresponding output, the number of channels in each branch increases sequentially from 3 to 24, 48, 96, and 128, with convolution strides of 1, 2, 2, and 2 respectively, and the LeakyReLU activation function is used. The outputs of each branch are cropped to the same spatial size and then concatenated.

[0040] (2) Feature fusion and latent variable mapping The aligned three-scale features are concatenated along the channel dimension as follows: ; in, For the concatenated feature tensor, Concat indicates concatenation along the channel dimension; , 384 represents the number of rows and columns of the feature map after each branch is aligned, and 384 represents the total number of channels for the three 128-channel branches.

[0041] Flatten: ; in, This represents a 384-dimensional channel response vector; GAP indicates global average pooling, applied to all channels. The average value is taken from each spatial location; Flatten indicates the flattening operation, which arranges the pooling results into a vector.

[0042] Then, channel weights are generated through two layers of fully connected mapping: ; in, It is a 384-dimensional channel weight vector; and The learnable weight matrices have shapes of 96×384 and 384×96, respectively. and The corresponding 96-dimensional and 384-dimensional biases are used. ReLU sets the negative inputs to zero and retains the non-negative inputs, while Sigmoid maps the channel weights to (0,1).

[0043] Apply channel weights to the concatenated features: ; in, Let be the feature tensor after channel attention weighting, and ⊙ denote the element-wise multiplication after spatial broadcasting. Each attention weight applies to all channels within the corresponding channel. The channel has a spatial location; when the weight is close to 1, the response of the channel is basically preserved, and when it is close to 0, the response is suppressed.

[0044] Subsequently, the 1×1 convolution reduced the number of channels from 384 to 256: ; in, is the compressed 256-channel feature tensor; Conv represents the convolution mapping, with its subscript 1×1 indicating the kernel size and the superscript 384→256 indicating the number of input and output channels; LeakyReLU preserves a small slope in the negative input region.

[0045] Flatten the feature tensor: ; in, for The result after flattening is 256× h × w 3D encoded vector.

[0046] The two fully connected headers map the encoding vector to the mean vector and log-variance vector of the latent variable distribution, respectively: ; ; in, It is a 60-dimensional mean vector. It is a 60-dimensional log-variance vector; and These are the standard deviation and variance of each latent variable, respectively. This indicates a fully connected mapping.

[0047] (3) Obtain the initial latent state by reparameterization: ; in, It is a 60-dimensional standard normal random vector. It is a 60×60 identity matrix; This represents element-wise multiplication. The exponential transformation yields the stepwise standard deviation. This represents the 60-dimensional initial latent state after reparameterization. The complete encoder process from the first frame's three-channel observation field to the latent variable distribution is shown below. Figure 2 .

[0048] 3. Gated data - Bounded guided dual-branch neural ODE (ordinary differential equation) The encoder encodes the three-channel observation field on the first day of the window into a 60-dimensional initial latent state. .by As an initial value, the neural ODE calculates the latent state change rate under the combined action of the data branch and the Tanh bounded guiding branch, and adjusts the contribution of the two branches dimension by dimension through the gating vector.

[0049] After being solved by ODE, the model obtains the corresponding latent states at 8 discrete time steps within the window, and decodes them into velocity fields one by one according to the time step.

[0050] It should be noted that the Tanh-guided branch is a multilayer perceptron that learns latent state change patterns and is not directly embedded in the ocean dynamics equations. The model first decodes the latent state at each time step into a velocity field, and then calculates the physical soft constraint loss based on the predicted velocity field. These losses, along with the reconstruction error, participate in the optimization and update the encoder, neural ODE, and decoder parameters through backpropagation.

[0051] The latent state evolution equation of the neural ODE is: ; in, For the first The 60-dimensional potential state of the day, abbreviated on the right as follows: ; A gated vector for assigning weights to the two branches one dimension at a time. For the data branch, a multilayer perceptron of 60→128→256→128→60 is used, and the hidden layer is activated using Softplus; The Tanh-bounded guided branch is used, employing a 60→128→128→60 multilayer perceptron. The two hidden layers use Tanh activation, and the output layer is a linear mapping.

[0052] The gating vector is defined as: ; in, The gating weight matrix is ​​60×60. The bias vector is 60-dimensional, and both are learnable parameters. The 60 components are located in (0,1), and the contributions of the data branch and the Tanh bounded guiding branch are adjusted dimension by dimension.

[0053] Using the encoder output as the initial value, the neural ODE is solved using the fourth-order Runge-Kutta method, and the gradient is calculated during training using the adjoint method to obtain the latent states at 8 time steps within the window: ; in, Let Θ be the latent state change function of the gated ODE mentioned above, and let Θ be the set of learnable parameters of the ODE module. This is the initial potential state. For the integration time node, take The time coordinate uses the day number within the window, with an interval of 1 between adjacent integration time nodes. This indicates the ODE integral process that supports the adjoint method of inverse differentiation. This represents the fourth-order Runge–Kutta method.

[0054] The eight latent states are input into the decoder to form the predicted velocity field sequence for each time step. During the training phase, the confidence-weighted reconstruction error and physical soft constraint loss are calculated based on this predicted sequence, and the model parameters are updated through backpropagation using the adjoint method.

[0055] 4. Stream function-assisted decoder The neural ODE outputs a sequence of latent variables, which still needs to be reconstructed into a two-dimensional velocity field. This invention sets up two parallel decoding branches: a direct velocity branch retains local expressive power, and a stream function branch provides an auxiliary structure with low divergence tendency, achieving fusion through learnable global coefficients. Spatial decoding, stream function differencing, and branch fusion are described below.

[0056] (1) Dual-branch decoder structure The total dimension of the latent state variables is 60, with the first 50 dimensions entering the decoder: ; in, Indicates the first The latent state of the day, Its decoded subvector consisting of dimensions 1 to 50, with superscript... This is used to label the subvector. The remaining 10 dimensions participate in the full latent state update through the ODE fully connected layer, indirectly affecting the decoded subvector in subsequent time steps.

[0057] Fully connected layers restore the decoded subvectors to spatial feature tensors: ; in, This is the initial feature tensor for decoding with 256 channels, and Reshape represents shape rearrangement; and Let the number of rows and columns of this feature map be taken as follows: , ; This indicates rounding down. This indicates taking the larger value.

[0058] The direct speed branch consists of three layers of transposed convolutions, LeakyReLU activation, and one layer of two-dimensional convolutions, with the number of channels increasing proportionally. Changes; its output is recovered to [value] by bilinear interpolation. And mapped to the normalized velocity range via Tanh: ; in, For the direct velocity branch in the 1st Normalized prediction field for the day, and These are its east-west and north-south velocity components, respectively, with subscripts. Identify the direct speed branch.

[0059] The stream function branches use the same upsampling channel structure, but ultimately output a single-channel scalar field: ; in, From Learnable decoding mappings for generating stream function fields. For the goal A stream function on the grid is used to assist the scalar field. The same mapping is repeated at each time step. The time index is omitted from the notation.

[0060] (2) Stream function constraints and branch fusion In the stream function-assisted field, central differences with unit intervals are used along the grid row and column indices. The internal grid is defined as follows: ; ; in, and These are the central difference operators along the row and column indices, respectively. For grid The stream function auxiliary scalar value at the location; the row and column directions correspond to the north-south and east-west directions respectively, and the denominator 2 represents the total length of the interval between two adjacent indices.

[0061] Two auxiliary velocity components are constructed from the above difference: ; in, and These are the east-west and north-south auxiliary components before Tanh mapping, with superscripts. Identify unmapped values. Copy the padding boundary for both components and align them using bilinear interpolation. After Tanh mapping, we get: ; in, For the first The normalized auxiliary prediction field of the day, and It has its east-west and north-south components. This branch is used to characterize the relative spatial structure.

[0062] The predicted fields of the two branches are fused using the following formula: ; ; in, For the first The fusion normalization prediction field of the day, The fusion ratio of the direct velocity branch, The fusion ratio of auxiliary branches for the stream function; This is a global scalar learned through backpropagation. It is shared across batches, time, speed channels, and spatial locations. See the dual-branch decoding and fusion process. Figure 3 .

[0063] 5. Marine physical soft confinement system The stream function auxiliary branch only constitutes a part of the final output; therefore, during the training phase, physical soft constraints need to be applied to the fused normalized prediction sequence to suppress local anomalous convergence, vorticity abrupt changes, time jumps, and spatial noise. The two components of the fused velocity field are written as: ; in, and They are respectively fusion prediction fields Within the window Normalized east-west and north-south velocity components of the sun .

[0064] To calculate the spatial variation of the predicted field, a first-order spatial difference is defined: ; ; in, For the scalar field on the target grid, For it in the grid The value at; and They are respectively along the east-west direction and north-south direction The first-order central difference operator; , This represents the spatial grid spacing in the corresponding direction.

[0065] The second-order spatial difference is defined as: ; ; in, and respectively along , The second-order central difference operator for the direction. The above difference is computed on an internal grid with the desired neighborhood.

[0066] The horizontal Laplace operator is defined as: ; in, It is a horizontal Laplace operator, obtained by adding the second-order central differences in two directions.

[0067] (1) Spatial-temporal consistency constraints of divergence and vorticity The divergence constraint is defined as: ; in, For divergence loss; The divergence outside the sea area is set to zero based on the spatial location corresponding to the difference result. This represents the arithmetic mean of the values ​​over the corresponding batch, time, and internal grid, including the zero-valued locations. This loss is used to suppress anomalous local convergence and divergence in the fused velocity field.

[0068] Relative vorticity is written as: ; in, The discrete relative vorticity field is obtained by spatial difference and sea area masking for normalized predicted velocity.

[0069] For adjacent discrete time steps, construct spatiotemporal consistent regular residuals of vorticity: ; in, For the first in the window The day arrives Daily vorticity spatiotemporal consistency normalized residual The difference in vorticity between adjacent time steps. This is the local discrete advection term; This corresponds to 7 pairs of adjacent dates within an 8-day window. The residual is constructed using the normalized velocity and the day number within the window, and is used to regularize the spatiotemporal variations of the vorticity field.

[0070] vorticity loss is defined as: ; in, This is the mean square loss of the vorticity residuals, used to reduce discontinuous vorticity evolution in the prediction sequence.

[0071] (2) Time continuity constraint To limit anomalous jumps in the velocity field between adjacent time steps, the following definition is made: ; in, The vector representing the change in normalized predicted velocity over two consecutive days, after... The change outside the sea area after processing is 0.

[0072] Write the velocity change vector as two components and calculate the magnitude after adding the numerical stabilization term: ; in, and The first Each sample window, adjacent day pairs Grid The normalized changes in east-west and north-south velocity at the location; subscript , Identify the two velocity components. The magnitude of the velocity change vector after adding a numerical stabilization term, within the square root. Used to ensure numerical stability.

[0073] The time continuity loss is written as: ; in, For speed channel index; The threshold for abrupt changes in the normalized velocity change modulus; coefficient 2 is the fixed weight of the second term relative to the first term. The first term penalizes the square of the velocity changes over all adjacent days, and the second term further penalizes... Exceed The part.

[0074] (3) Spatial Laplace smoothing constraint Spatial Laplace smoothing loss is defined as: ; in, The Laplace mean square loss of the two velocity components after sea area masking is given. The value is taken based on the spatial location corresponding to the difference result. This loss is used to reduce isolated grid noise and local spatial oscillations.

[0075] (4) Velocity modulus constraint To penalize excessively large normalized prediction speed modulus, we define: ; in, The velocity modulus constraint loss. The square root is the value of the numerical stability constant. The normalized velocity modulus is then used; ReLU is used to extract the portion exceeding 1, and then the average is taken for batch, time, and spatial location.

[0076] (5) Physical soft constraints The physical soft constraint loss involved in progressive modulation is: ; in, This is the loss due to physical soft constraints; , , and These are the spatial-temporal consistency losses of divergence and vorticity, temporal continuity, and spatial Laplace smoothing, respectively. , , and These are the weighting coefficients for the corresponding loss terms.

[0077] The equivalent weights of the divergence term are combined according to the following formula: =0.18; in, For divergence penalty weights, The weights are for stream function consistency penalties. Both terms use the divergence loss of the fused velocity field, therefore the weights are summed to obtain... The asterisk indicates the equivalent weight after merging.

[0078] The remaining weights are: ; ; .

[0079] The aforementioned fixed weights are used to balance the relative contributions of each physical constraint during training.

[0080] 6. Model Training Objectives and Incremental Optimization Strategies (1) Overall training objectives Incorporating comprehensive data reconstruction, latent variable regularization, physical soft constraints, and CMEMS background priors, the training objective (total loss function) is written as: ; Where L is the total loss; The Huber reconstruction loss is weighted by confidence level. For latent variables, KL regularization term, For physical soft constraint combination loss, For velocity modulus constraint loss, For CMEMS background prior loss; , and These are the fixed weights for the corresponding loss terms; To follow the training rounds A varying scaling factor is used to adjust... The intensity of the effect.

[0081] This embodiment takes , , All three are fixed during training, controlling the relative strength of latent variable regularization, background prior, and velocity modulus penalty, respectively.

[0082] (2) Latent variable regularization and background prior The KL regularization term of a variational encoder is defined as: ; in, This represents the KL regularization term for the posterior distribution of the latent variables relative to the standard normal distribution. The square, logarithmic, and exponential operations within the square brackets are calculated for each latent variable. The mean is calculated for the batch sample and the 60 latent variable dimensions. This term is used to limit the extent to which the posterior distribution deviates from the standard normal distribution.

[0083] CMEMS background velocity uses the same fixed-scale normalization as AIS, and its background prior loss is: ; in, To normalize the CMEMS background velocity field, For the corresponding valid background indication field; For the sample window index in the batch, The day number within the window. For speed channel index, , For grid row and column indexes; and These represent the predicted and background values ​​at the same time step, channel, and grid location, respectively. The value is 1 when the corresponding background component is valid and located in the sea area, and 0 otherwise. The denominator is the number of valid background components. This loss is only calculated when the denominator is greater than 0. CMEMS provides a background reference through this loss, and no additional input to CMEMS is required during the reconstruction stage.

[0084] (3) Progressive training and overlapping window reconstruction Training progress is represented by the ratio of the current training round to the total number of rounds: ; in, To normalize the training progress, For the current training round, To determine the total number of training rounds, this embodiment uses E=200; This indicates taking the smaller value so that the training progress does not exceed 1.

[0085] Based on the training progress, the physical constraint scaling factor is set in three stages: ; in, The scaling factor for physical constraints is: it increases linearly from 0.1 to 1.0 when the training progress is below 0.3; it remains at 1.0 when the progress is between 0.3 and 0.6; thereafter, it increases to 1.5 when the progress reaches 1 as the training progress increases.

[0086] The training samples consist of an 8-day window; only the observations from the first day are entered into the encoder, the 8-day AIS inversion fields are used as the corresponding day's reconstruction targets, and the 8-day prediction sequences are used for physical loss. The CMEMS same-day field is only used for background prior loss. This embodiment uses Adam for 200 training rounds with an initial learning rate of... Weight decay Cosine annealing hot restart is used, with an initial restart cycle of 50 rounds, a cycle multiplication factor of 2, and a minimum learning rate. The evaluation uses the checkpoints of the models with the lowest total training loss.

[0087] During the reconstruction phase, model parameters are frozen and predictions are made window by window. The weights of each day's position within a window are defined as follows: ; in, For the first in the window The weight of each position, ; The window length is T=8 in this embodiment, i.e., 8 consecutive days. The center of the weighting function is T / 2=4, and the scale is T / 4=2.

[0088] Multi-window forecasts for the same date are synthesized using the following formula: ; in, For date The multi-window weighted normalized prediction field, with the horizontal line indicating weighted synthesis; The starting date index is The window for dates Normalized prediction; For the set of date indices covered by this window, take ; For date Within this window, For the corresponding weights; sum and iterate through all covered dates. The window is the sum of weights in the denominator.

[0089] After weighted synthesis, the two normalized velocity components are multiplied by 3 m / s to restore physical units; Gaussian smoothing is used to reduce isolated mesh noise, and the final flow velocity modulus is limited to within 3 m / s. The relationship between the reconstruction sequence, constraints, and parameter updates is detailed in [link to documentation]. Figure 4 .

[0090] Example 2: To evaluate the technical effect of the method of the present invention, the following simulation experiment was conducted.

[0091] 1. Data and Research Area (1) Study area The research area of ​​this invention is located in the core area of ​​the Agulhas Current in the southeastern sea off the coast of South Africa, with a spatial range of 23°E–29°E and 36°S–32°S. The Agulhas Current is one of the strongest western boundary currents in the Southern Hemisphere, transporting warm water southwestward along the east coast of South Africa. In the study area, it forms a clear main axis, a significant lateral velocity gradient, and is accompanied by current band bending and mesoscale disturbances (Le Goff C, Boussidi B, Mironov A, et al. 2021. Monitoring the greater Agulhas Current with AIS data information. Journal of Geophysical Research: Oceans, 126(5): e2021JC017228, doi: 10.1029 / 2021JC017228; Benaïchouche Simon, et al. 2021. Unsupervised reconstruction of sea surface currents from AISmaritime traffic data using trainable variational models. Remote Sensing, 13(16): 3162, doi: 10.3390 / rs13163162).

[0092] The selection of this region is based on two conditions directly related to this invention. Firstly, the strong currents and their complex spatial structure present a clear challenge to reconstructing continuous flow fields under sparse observation conditions. Secondly, the southeastern coast of South Africa is an important international shipping route, with ship tracks concentrated along major shipping lanes, while significant observation gaps remain in the open sea areas far from these lanes. Therefore, the simultaneous inclusion of both densely populated and sparsely populated AIS observation areas within the same study area reflects the non-uniform spatial support characteristics of ship observations in real-world sea areas.

[0093] (2) AIS inverted sea surface current field data The AIS data covers June 2022 and contains 612,558 dynamic records after data cleaning. Each record contains information such as UTC time, ship latitude and longitude, speed over land (SOG), heading over land (COG), and heading. These variables are the basis for subsequent sea surface current speed inversion.

[0094] Based on the established direction-aware adaptive clustering and ship motion vector decomposition method, AIS records that are spatially adjacent and whose heading information meets the solution requirements are organized into effective sample clusters, and the eastward speed is obtained from the ship motion information within each cluster. and northbound speed After outlier filtering, the data was converted to daily mean sea surface currents on a 1 / 12° grid. Each valid grid point contains longitude, latitude, and longitude. , The `count` field indicates the number of valid AIS samples that participated in the grid point inversion and aggregation.

[0095] (3) CMEMS data CMEMS data is derived from the Global Ocean Physics Analysis and Forecast (product ID: GLOBAL_ANALYSISFORECAST_PHY_001_024; https: / / data.marine.copernicus.eu / products) product published by the Copernicus Marine Service. This invention utilizes the daily sea surface velocity data with a spatial resolution of 1 / 12° to extract the eastward velocity. and northbound speed And extract the effective ocean grid within the study area in June 2022.

[0096] (4) Drifting buoy data The drifting buoy data were obtained from the Global Drifter Program (GDP) of the Atlantic Oceanic and Meteorological Laboratory (NOAA / AOML) of the National Oceanic and Atmospheric Administration (https: / / www.aoml.noaa.gov / phod / gdp). Four surface drifting buoys that entered the study area in June 2022 were selected, with WMO numbers 1401817, 1401818, 1601639, and 1601724. The original data contained 86 valid records in the study area. Each record included the WMO number, observation time, latitude and longitude, and eastward and northward drift velocities (ElipotShane, Lumpkin R, Perez RC, et al. 2016. A global surface drifter data set at hourly resolution. Journal of Geophysical Research: Oceans, 121(5): 2937-2966, doi: 10.1002 / 2016JC011716).

[0097] Buoy records are Lagrange observations obtained as the buoy moves. This invention preserves the original time, position, and velocity components of each record.

[0098] (5) Sentinel-1 OCN radial velocity data Sentinel-1 data was obtained from the CopernicusBrowser (https: / / browser.dataspace.copernicus.eu) provided by the Copernicus Data Space Ecosystem. This invention acquired the Sentinel-1A interferometric wide-swath (IW) Level-2 oceanographic products (OCN) intersecting the study area in June 2022 and extracted the radial velocity (RVL) and related quality control variables. RVL was estimated by SAR Doppler shift and represents the signed velocity component in the radar line-of-sight direction. Eleven valid transits were obtained during the study period, on June 2, 5, 7, 12, 14, 17, 19, 21, 24, 26, and 29.Due to the different satellite orbit positions and image slices corresponding to different transits, the coverage of each SAR observation in the study area varies, and not every transit covers the entire study area (Moiseev A, Johannessen JA, Johnsen H. 2022. Towards retrieving reliable ocean surfacecurrents in the coastal zone from Sentinel-1 Doppler shift observations. Journal of Geophysical Research: Oceans, 127: e2021JC018201, doi: 10.1029 / 2021JC018201; Hauser Danièle, et al. 2023. Satellite remote sensing of surfacewinds, waves, and currents: Where are we now? Surveys in Geophysics, 44(5):1357-1446, doi: 10.1007 / s10712-023-09771-2; Martin ACH, Gommenginger CP, Jacob B, Staneva J. 2022. First multi-year assessment of Sentinel-1 radialvelocity products using HF radar currents in a coastal environment. RemoteSensing of Environment, 268: 112758, doi: 10.1016 / j.rse.2021.112758).

[0099] The OCN RVL record includes radial velocity standard deviation, satellite heading angle, radar incident angle, landmarks, land coverage, and Doppler estimation confidence flags. For each transit, land pixels with land coverage of at least 99%, valid confidence flags, and radial velocity standard deviations not exceeding 0.1 m / s are first selected. The systematic bias is estimated using the median radial velocity of these land pixels, and then subtracted from the ocean pixels in the same transit. Ocean pixels used for comparison must meet the following criteria: landmarks are oceanic, land coverage is 0%, confidence flags are valid, and radial velocity standard deviations do not exceed 0.1 m / s. After correction, the data is aggregated to a 1 / 12° grid based on the reciprocal of the error variance, with each grid containing at least 5 pixels. The two-dimensional flow field is projected onto the radar line-of-sight direction based on the satellite heading angle and incident angle of the same transit and compared with the RVL. A single transit only provides one line-of-sight direction component and cannot replace the two-dimensional data. , Component verification.

[0100] 2. Experimental Results and Analysis (1) Evaluation scheme and control settings The PI-NODE model of this invention is evaluated from three aspects: spatial structure, independent buoy accuracy, and ablation of key components. First, the AIS inversion field, CMEMS daily field, and PI-NODE reconstructed field are compared on the same day. The CMEMS and PI-NODE flow fields are then projected onto the Sentinel-1 radar line-of-sight direction, and their consistency with SAR radial velocity is compared. Second, the two methods are compared on the same batch of buoy samples. , Component error.

[0101] To verify the effectiveness of each key module of the method, ablation experiments were further constructed. Each model was trained for 200 epochs using the same training samples, optimizer, number of training epochs, and fixed training random seeds. The checkpoint with the lowest training loss was used for evaluation. The complete model and its variants were computed on the same 20 inference random seeds and the same batch of 86 buoy records. , RMSE (Meaning).

[0102] (2) Results of global flow field reconstruction

[0103] To further verify whether the above spatial corrections are consistent with independent remote sensing observations, Sentinel-1 OCN radial velocity products from 16:53 UTC on the same day were selected for spatial comparison. SAR pixels after quality screening were aggregated to a 1 / 12° grid based on the inverse of the error variance; PI-NODE and CMEMS... , The components are projected onto the same radar line-of-sight direction and form 802 common samples at the same date and grid location.

[0104] To compare the radial errors of the two types of flow fields relative to SAR, the radial error advantage is defined as follows: ΔE = |V r CMEMS - V r ,SAR| - |V r ,PI-NODE - V r ,SAR| ; in, , and These are the SAR radial velocity, PI-NODE projected radial velocity, and CMEMS projected radial velocity at the same grid location, all three using the same radar line-of-sight direction. This represents the difference in absolute error between the two models relative to SAR. This indicates that PI-NODE is closer to the SAR radial velocity. This indicates that CMEMS is closer.

[0105] Across all 802 common grid cells, the radial RMSEs of PI-NODE and CMEMS relative to SAR were 0.347 m / s and 0.399 m / s, respectively, with PI-NODE showing a 13.1% reduction. It should be noted that the Sentinel-1 OCN product provides the radar line-of-sight velocity component, and this comparison only covers a single pass. Therefore, this result is used to supplement the explanation of the improved consistency of PI-NODE in this pass and observation direction, and does not replace the buoy pair. , Independent verification of components is not equivalent to the true value of the complete two-dimensional flow field.

[0106] (3) Independent buoy verification Drifting buoys were paired with PI-NODE and CMEMS results according to the observation date: PI-NODE read the reconstructed value of the target grid cell at 1 / 12° of the buoy's location, and CMEMS read the value of the grid point closest to the buoy's position from the same day's data. Only records with valid values ​​for both flow fields were retained, resulting in 86 common samples. Table 1 lists the verification results for both components. PI-NODE's... , The component RMSE values ​​were 0.377 m / s and 0.295 m / s, respectively, lower than those of CMEMS (0.453 m / s and 0.381 m / s), representing decreases of 16.8% and 22.6%, respectively; the corresponding R... 2 They increased from 0.508 and 0.621 to 0.659 and 0.773, respectively.

[0107] Table 1. Samples from the same batch of buoys , Component verification results .

[0108] Figure 5 Provide the density scatter distribution for the same sample. PI-NODE in , The scatter points on both components are more concentrated near the 1:1 reference line, especially The high-density sample bands of the components are closer to the reference line; the CMEMS scatter plots are more dispersed within the same value range. This result is consistent with the error indices in Table 1; MSE is the mean square error, RMSE is the root mean square error, and R0 is the root mean square error. 2 The coefficient of determination.

[0109] (4) Ablation test To examine the contribution of each key component to the reconstruction accuracy, this invention sets up a complete model and five retraining ablation variants. In Table 2, number 1 represents the complete model; number 2 removes spatial confidence information, that is, sets the encoder confidence channel to zero and replaces the grid-by-grid reconstruction weights with the global average confidence; number 3 replaces the multi-scale encoder with a single 3×3 convolution branch; number 4 replaces the gated dual-branch ODE with a single data branch; number 5 removes the stream function-assisted decoding branch; and number 6 removes physical soft constraints.

[0110] For each inference random seed, the RMSE rate of change of numbers 2–6 relative to number 1 was first calculated, and then the mean and standard deviation of the 20 paired results were calculated. The full model and the corresponding variant used the same random seed, thereby reducing the interference of latent variable sampling differences on component comparisons. The results are shown in Table 2.

[0111] Table 2. Paired ablation results of key PI-NODE components ; Note: √ indicates that the component is retained, × indicates that the component is removed or replaced; the values ​​are the mean ± standard deviation of the rate of change of RMSE relative to the corresponding full model under 20 paired inference random seeds, with positive values ​​indicating an increase in error. Each model is trained for 200 epochs and the checkpoint with the lowest training loss is used.

[0112] Table 2 shows the five simplified settings. , All component average RMSE change rates were positive. After replacing the multi-scale encoder with a single 3×3 branch in number 3, , The component RMSE increased by 12.8%±1.7% and 7.8%±2.6% respectively, showing the largest increase among the tested settings, indicating that parallel extraction of features at different spatial scales helps improve the accuracy of sparse AIS flow field reconstruction. After removing spatial confidence information in setting 2, the RMSE of the two components increased by 7.1%±1.7% and 2.8%±1.6% respectively, indicating that the reliability of grid-by-grid observations helps the model utilize spatially unevenly distributed AIS information.

[0113] After replacing the gated dual-branch ODE in item 4 with a single data branch, , The RMSE of the components increased by 8.2%±1.5% and 3.3%±1.9%, respectively; after removing the stream function-assisted decoding branch in component 5, the two components increased by 7.9%±1.7% and 3.9%±1.7%, respectively; after removing the physical soft constraint in component 6, the two components increased by 5.2%±0.9% and 4.0%±2.3%, respectively. Therefore, spatial confidence, multi-scale coding, gated dual-branch latent state evolution, stream function-assisted decoding, and physical soft constraints together constitute the effective functional design of the complete model.

[0114] This invention addresses the problem that sparse AIS inversion observations cannot directly form a continuous sea surface current field, and proposes the PI-NODE model framework. This model uses AIS inversion data... , and confidence level As encoder input, a continuous flow field is generated through multi-scale parallel encoding, gated dual-branch neural ODE, and stream function-assisted decoding, and oceanic physical soft constraints are applied during the training phase. CMEMS provides a large-scale background reference only through prior loss and is not used as encoder input or involved in the reconstruction phase.

[0115] In June 2022, experiments in the Agulhas Current region evaluated the method of this invention from three aspects: spatial structure, accuracy of independent buoys, and ablation of key components. Global results show that the PI-NODE model provided by this invention can express local spatial variations in the AIS support region while forming a continuous flow field. On 86 drifting buoy records paired on the same day, the PI-NODE model... , The component RMSE values ​​were 0.377 m / s and 0.295 m / s, respectively, lower than those of CMEMS (0.453 m / s and 0.381 m / s); the corresponding R... 2The values ​​improved from 0.508 and 0.621 to 0.659 and 0.773, respectively. On June 2nd, across 802 paired SAR grids, the radial RMSEs for PI-NODE and CMEMS were 0.347 m / s and 0.399 m / s, respectively, demonstrating that the PI-NODE model had lower errors in the transit and radar line-of-sight directions. In the retraining ablation, the two-component average RMSE increased for all five simplified settings, further supporting the effectiveness of the model's key functional design.

Claims

1. A sparse AIS ocean current field reconstruction method based on physical information neural networks, characterized in that, Includes the following steps: (1) Construct the PI-NODE model, which includes a multi-scale parallel encoder, a gated data-bounded guided dual-branch constant differential equation, and a stream function-assisted decoder; (2) The multi-scale parallel encoder extracts features from the input three-channel observation field, which includes the east-west and north-south velocity fields and the corresponding confidence fields; the outputs of each branch are aligned to the same spatial size, and feature fusion and latent variable mapping are performed in the channel dimension to obtain the initial latent state; (3) The gated data-bounded guided dual-branch constant differential equations take the initial latent state as the initial value, calculate the latent state change rate under the combined action of the data branch and the Tanh bounded guided branch, and obtain a discrete latent state sequence after solving. (4) Each latent state is restored to spatial features through a fully connected layer, and then enters the stream function-assisted decoder for decoding to obtain the sea surface current field prediction sequence within the observation window; (5) The prediction results of the overlapping windows are weighted and synthesized to obtain the daily continuous sea surface flow field.

2. The sparse AIS ocean current field reconstruction method based on physical information neural network according to claim 1, characterized in that: The multi-scale parallel encoder includes three parallel convolutional branches, with the kernel scale of the first two convolutional layers respectively adopting... , and Features are extracted from different spatial ranges.

3. The sparse AIS ocean current field reconstruction method based on physical information neural network according to claim 1, characterized in that: The stream function-assisted decoder includes two decoding branches set in parallel: a direct speed branch and a stream function branch. The direct velocity branch preserves local expressive power, while the stream function branch provides an auxiliary structure with low divergence tendency; the two branches are fused through learnable global coefficients.

4. The sparse AIS ocean current field reconstruction method based on physical information neural network according to claim 1, characterized in that: The PI-NODE model also includes a confidence-weighted observation model, which constructs grid confidence based on the number of samples and the consistency of neighborhood flow direction to adjust the AIS reconstruction error.

5. The sparse AIS ocean current field reconstruction method based on physical information neural network according to claim 1, characterized in that: During the training phase, the PI-NODE model compares the predicted velocity field for each day with the AIS inversion velocity for that day to calculate the reconstruction error.

6. The sparse AIS ocean current field reconstruction method based on physical information neural network according to claim 1, characterized in that: During the training phase, the PI-NODE model constructs a physical soft-constraint loss based on the changes and spatial structure of the predicted velocity fields over multiple consecutive days. ； in, This is the loss due to physical soft constraints; , , and These are the spatial-temporal consistency losses of divergence and vorticity, temporal continuity, and spatial Laplace smoothing, respectively. , , and These are the weighting coefficients for the corresponding loss terms.

7. The sparse AIS ocean current field reconstruction method based on physical information neural network according to claim 1, characterized in that: The total training loss of the PI-NODE model during the training phase is: ； Where L is the total loss; Huber loss weighted by confidence level; For latent variable regularization, For physical soft constraint loss, For velocity modulus constraints, For CMEMS background priors; , and These are the fixed weights for latent variable regularization, velocity modulus constraint, and background prior term, respectively. The physical scaling factor varies with the training epoch e, only adjusting... The intensity of the effect.

8. The sparse AIS ocean current field reconstruction method based on physical information neural network according to claim 1, characterized in that: The training samples consist of consecutive multi-day windows; the normalized AIS velocity field and confidence field of the first day are fed into a multi-scale parallel encoder to generate the initial latent state; the gated data-bounded guided dual-branch constant differential equations evolve the latent state at multiple time steps, including the first day, starting from this initial latent state.

Images