A sea surface current field retrieval and reconstruction method based on stationary satellite high-frequency images

Abstract

The application discloses a kind of sea surface flow field inversion and reconstruction method based on stationary satellite high frequency image, belong to marine remote sensing technical field.It includes the following steps: satellite image is preprocessed and locally enhanced;Initial inversion is carried out using dense optical flow algorithm for single time phase flow field;The flow field data cube of space-time grid uniformization is constructed;Time series reconstruction is carried out using two-way Kalman filtering, through forward recursion, backward recursion and two-way fusion, using time series information to correct abnormal vector and fill in missing data;Finally, based on the physical consistency principle of seawater approximate incompressibility, the reconstructed flow field is subjected to iterative divergence constraint, and the non-physical disturbance is weakened.The application effectively solves the problem that the traditional optical remote sensing flow measurement method is easy to produce abnormal vector and poor spatiotemporal consistency under complex sea conditions such as cloud cover and weak texture by introducing time series correlation and physical constraint, significantly improves the accuracy, stability and physical reasonableness of sea surface flow field inversion.

Description

Technical Field

[0001] This invention belongs to the field of marine remote sensing and marine dynamic environment monitoring technology, specifically relating to a method for inverting and reconstructing sea surface flow fields based on high-frequency geostationary satellite images. Background Technology

[0002] Ocean surface currents are a key element in describing the marine dynamic environment and have significant applications in marine disaster prevention and mitigation, pollutant dispersion tracking, maritime search and rescue, and marine ecological research. Traditional ocean current observation methods mainly rely on buoy arrays (such as Argo buoys and drifting buoys) and shipborne ADCP (Advanced Digital Acceleration Profiler) observations. While these methods offer high measurement accuracy, they are limited by equipment deployment costs and maintenance difficulties. Furthermore, these methods suffer from sparse spatial coverage and low temporal resolution, making it difficult to achieve large-scale, high-time synchronous monitoring. With the development of satellite remote sensing technology, high-frequency observations based on geostationary ocean color satellites (such as GOCI-II) have provided a new approach to current field acquisition. The core idea is to process continuous temporal ocean color images (such as chlorophyll and suspended matter concentrations) to analyze the temporal movement patterns of marine tracers, thereby achieving continuous estimation of the surface velocity field. These methods can acquire high spatial resolution velocity information without direct contact with the water body, offering advantages over traditional methods such as wide coverage and all-weather operation.

[0003] Existing ocean remote sensing velocities based on satellite imagery mainly include the Maximum Cross-Correlation (MCC) method and optical flow estimation algorithms (such as dense optical flow methods). The MCC method determines pixel displacement by searching for the maximum correlation coefficient between sub-windows of two consecutive image frames; optical flow methods estimate pixel motion by modeling grayscale changes in local neighborhoods of the image. These methods can achieve certain velocity inversion results under conditions of clear tracer texture and high observation quality. However, under complex actual sea conditions, the image matching process is prone to ambiguity or complete failure due to factors such as clouds, solar flares, sensor noise, weak sea surface texture, and local non-rigid deformation. This results in poor computational stability of traditional methods, easily disrupting the continuity of spatiotemporal evolution. The inverted flow field often exhibits numerous spurious vectors with distorted directions and abnormal amplitudes (such as local or patchy flow direction deflections) in space, and severely lacks temporal continuity.

[0004] In recent years, to improve the accuracy of inversion results, existing studies often employ post-processing strategies such as smoothing, interpolation, or simple averaging to correct the inverted flow field. However, traditional linear smoothing or mathematical interpolation typically lacks characterization of the temporal correlation and physical consistency of fluid motion, easily leading to the loss of structural information during correction, or failing to effectively identify and correct pseudo-vectors that are "incorrect in direction but continuous in value." Especially in cases of fluctuating observation quality, cloud edges, and weakly textured regions, non-physical divergence and convergence, as well as temporal discontinuities, may still occur, making it difficult to meet the stability and reliability requirements of operational applications. Although data-driven methods such as deep learning have shown potential in some tasks, they usually rely on large-scale, high-quality ground truth data for training. However, high-density, high-precision ground truth flow field data in marine environments are difficult to obtain, hindering their widespread application.

[0005] Therefore, in ocean flow field inversion missions based on geostationary satellites, it is urgent to explore a method that can adaptively identify and correct anomalous vectors (including erroneous flow direction vectors, direction deflection vectors, etc.) in the inverted flow field under the influence of factors such as cloud cover, weak texture, and observation noise, and achieve spatiotemporal consistent reconstruction by utilizing the dynamic correlation of the time dimension. Summary of the Invention

[0006] This invention addresses the problems of existing technologies by providing a method for sea surface current field inversion and reconstruction based on high-frequency geostationary satellite imagery. It addresses issues such as the generation of anomalous vectors (including erroneous flow direction vectors and direction deflection vectors) in existing satellite imagery motion estimation under complex observation conditions such as cloud cover, poor spatiotemporal consistency, and high noise in single-frame inversion. This invention fully leverages the advantages of high-frequency continuous observation by geostationary satellites and introduces a bidirectional Kalman filter time-series reconstruction mechanism with physical consistency constraints. The method can adaptively fuse forward prediction and backward smoothing updates based on time series during low-quality observation stages, such as cloud cover changes and abrupt changes in observation quality, to improve the spatiotemporal continuity of the flow field. Furthermore, through physical consistency constraints and inertial smoothing in the time dimension, it enhances the evolution continuity of the flow field under complex observation conditions, correcting local erroneous flow directions and mitigating non-physical disturbances. Through these processes, this method improves the physical and dynamic consistency of the inverted flow field and provides a highly stable and physically plausible technical solution for long-term continuous monitoring of ocean dynamic processes.

[0007] To address the above technical problems, this invention provides the following technical solution: a method for inverting and reconstructing sea surface current fields based on high-frequency geostationary satellite imagery, comprising the following steps:

[0008] Step S1: Acquire geostationary satellite image data of the target sea area at multiple consecutive times, preprocess the image data to enhance the texture features of the tracer, and generate an initial validity mask for identifying invalid observation areas;

[0009] Step S2: Based on the enhanced images of adjacent time moments, a motion estimation algorithm is used to calculate the pixel displacement vector, obtain the initial flow field observation value at a single time moment, and smooth the initial flow field observation value. Invalid regions are marked according to the initial validity mask.

[0010] Step S3: Spatiotemporally align and uniformly rasterize the initial flow field observations at multiple times, map them to a standard geospatial grid, construct a multidimensional spatiotemporal flow field data cube containing spatial and temporal dimensions, and simultaneously integrate the validity masks of each time to form a spatiotemporal observation state matrix.

[0011] Step S4, reconstructing the flow field data based on a bidirectional temporal filtering method, including:

[0012] Forward recursion: Recursively recursively from front to back in time sequence, using the estimation results of the previous time moment to predict the current time moment, and updating it in combination with the valid observations of the current time moment. For invalid observation areas, the observation update weight is reduced or the update is skipped, so that the estimation results maintain the temporal continuity.

[0013] Backward recursion: Recursively recursively from back to front in time sequence, using the estimation results of the next time moment to predict the current time moment, and updating it by combining the effective observations of the current time moment;

[0014] Two-way fusion: The forward recursive estimation results and the backward recursive estimation results are weighted and fused according to their respective uncertainties to obtain the time-series reconstructed flow field;

[0015] Step S5: Based on the principle of physical consistency of fluid motion, a divergence constraint is applied to the temporal reconstructed flow field. The divergence gradient is calculated iteratively and the velocity vector is corrected so that the flow field divergence tends to meet the continuity requirement. The result of spatiotemporally continuous sea surface flow field inversion after correction by physical consistency constraint is output.

[0016] Furthermore, step 2 described above includes the following sub-steps:

[0017] S2.1 Using the Farneback dense optical flow algorithm, the signal intensity in the local neighborhood of the image is approximated by a quadratic polynomial, as follows:

[0018]

[0019]

[0020] In the formula, p= The spatial coordinate vector of the pixel. This represents the coefficient matrix of the quadratic terms, describing the shape of the quadratic surface in the local neighborhood. , Let b1 and b2 represent the constant term (scalar), b1 and b2 represent the coefficient vector of the first-order term, describing the gradient direction, and T represent the transpose. d is the global displacement, satisfying: ,

[0021] The displacement of each pixel is obtained by solving the quadratic polynomial using the least squares method. It employs a multi-level pyramid strategy to capture displacements at different scales, with parameter settings including: number of pyramid levels, average window size, polynomial order, and polynomial Gaussian standard deviation.

[0022] S2.2 Apply two-dimensional Gaussian convolution smoothing to the horizontal velocity component u and the vertical velocity component v respectively, and define the Gaussian kernel function:

[0023] ,

[0024] In the formula, The standard deviation of the Gaussian kernel is represented. This represents the weight value at spatial coordinates (x, y), used to perform a weighted average of the pixel to be smoothed and its neighborhood.

[0025] Calculate the smoothed flow field observations for:

[0026] ,

[0027] In the formula, This represents the original optical flow inversion result at pixel (x,y) in frame t.

[0028] S2.3 Outlier Removal: If the flow rate amplitude of a pixel... Greater than the maximum flow amplitude If it is not a divergence point in the algorithm, it will be forcibly marked as an invalid value and included in subsequent masks. middle.

[0029] Furthermore, in step S4 above, the bidirectional fusion specifically involves weighting and fusing the state estimate and uncertainty at time t obtained from forward recursion and the state estimate and uncertainty at time t obtained from backward recursion, based on the uncertainty information, to obtain the bidirectional reconstructed flow field:

[0030] ,

[0031] in, This represents the fusion weights of the forward estimation. This represents the uncertainty at time t during the forward recursion process. This represents the uncertainty at time t during the backward recursion process. This represents a tiny positive number introduced to ensure the stability of numerical calculations;

[0032] ,

[0033] in, This represents the final reconstructed ocean surface flow field state vector at time t, containing both meridional and latitudinal velocity components. This represents the forward recursive state estimation vector at time t, obtained through forward recursion. The backward recursive state estimation vector at time t is obtained by backward recursion, utilizing observation information from subsequent times.

[0034] Another aspect of the present invention provides an electronic device, including a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the computer program to implement the steps of any of the methods described in the present invention.

[0035] The present invention also provides a computer-readable storage medium having a computer program stored thereon, wherein the computer program, when executed by a processor, implements the steps of any of the methods described in the present invention.

[0036] Compared with the prior art, the beneficial technical effects of the present invention using the above technical solution are as follows:

[0037] 1. Improve the accuracy of sea surface current direction inversion and its consistency with physical structure.

[0038] To address the issue of traditional optical flow methods being susceptible to image noise interference and resulting in chaotic flow directions at cloud edges and in areas with weak texture, this invention introduces a spatiotemporal smoothing mechanism based on physical divergence constraints and bidirectional Kalman filtering to dynamically correct the flow field vectors. Experiments show that the flow field reconstructed by this invention exhibits a significantly reduced Directional Mean Absolute Error (MAE) compared to ocean numerical models (including but not limited to HYCOM). This indicates that this method can more accurately reconstruct the rotational morphology of sub-mesoscale structures such as vortices and fronts, and has significant practical value for applications highly sensitive to direction, such as marine oil spill tracking and search and rescue target drift prediction.

[0039] 2. Improve the stability of flow field inversion results under high-frequency remote sensing conditions.

[0040] By jointly processing high-frequency geostationary satellite image sequences, a continuity constraint is imposed on the flow field results in the time dimension, effectively overcoming the problem of unstable results that easily occur under conditions of image quality fluctuations in traditional single-phase or independent inversion methods between adjacent phases. Especially in the scenario of high-frequency geostationary satellite image observation, when some phases are obscured by clouds or the image quality drops sharply, this invention can utilize the backward smoothing strategy in the two-way mechanism to introduce high-quality observation results from adjacent phases, and perform reverse correction and modification of the flow field during the contaminated period, thereby significantly improving the consistency and global stability of the inverted flow field in the time series.

[0041] 3. Enhance the spatial continuity of the flow field and reduce noise interference.

[0042] In the multi-temporal fusion and reconstruction process, the consistency of the flow field in both spatial and temporal dimensions is comprehensively considered, and strict hydrodynamic divergence constraints are applied, resulting in a more continuous spatial distribution of the final output sea surface flow field. This mechanism significantly reduces discrete noise points common in traditional methods and adaptively identifies and eliminates directionally disordered local spurious vectors. Through spatiotemporal joint constraints, this invention can effectively suppress non-physical disturbances introduced by image noise, matching ambiguities, and other factors while maintaining the main flow direction and spatial structure characteristics. Attached Figure Description

[0043] Figure 1 This is an overall flowchart of the method of the present invention;

[0044] Figure 2 This is a schematic diagram of the timing reconstruction principle based on bidirectional Kalman filtering in the method of this invention;

[0045] Figure 3 This is a schematic diagram of the iterative divergence constraint logic steps based on physical consistency in the method of this invention;

[0046] Figure 4 This is a time evolution comparison chart of the flow field inversion accuracy index of the method of the present invention. Detailed Implementation

[0047] To better understand the technical content of the present invention, specific embodiments are described below in conjunction with the accompanying drawings.

[0048] In this invention, various aspects of the invention are described with reference to the accompanying drawings, in which numerous illustrative embodiments are shown. Embodiments of the invention are not limited to those depicted in the drawings. It should be understood that the invention is implemented through any of the various concepts and embodiments described above, as well as the concepts and embodiments described in detail below, because the concepts and embodiments disclosed herein are not limited to any particular implementation. Furthermore, some aspects of the invention disclosed may be used alone or in any suitable combination with other aspects of the invention disclosed.

[0049] To address the shortcomings of existing sea surface current field remote sensing inversion technologies, such as geostationary ocean color satellites (e.g., GOCI-II) which have the advantage of high-frequency observation, but are susceptible to interference from factors such as cloud cover and solar flares in complex sea conditions leading to spatiotemporal discontinuities in the inversion data, and the limited robustness of conventional motion estimation methods in weakly textured regions and the tendency to generate local spurious vectors, this embodiment proposes a sea surface current field inversion and reconstruction method based on high-frequency geostationary satellite imagery.

[0050] like Figure 1 As shown, the method of this invention mainly includes five core steps. First, the system mines weak ocean tracer textures through local adaptive enhancement processing; second, it obtains initial flow field observations using motion estimation methods such as dense optical flow; then, it constructs a bidirectional temporal recursive Kalman estimation model, combines an effectiveness mask to deweight or skip updates for invalid observations, and uses adaptive filling / smoothing to maintain the continuous evolution of occluded regions, while performing fusion updates on effective regions; finally, it introduces physical consistency iterative correction based on divergence constraints to weaken non-physical disturbances and improve the physical consistency of the reconstructed flow field.

[0051] The following is a detailed explanation of each step:

[0052] Step S1, Satellite Image Preprocessing and Local Texture Enhancement: Acquire geostationary satellite image data of the target sea area at multiple consecutive times, preprocess the image data to enhance the texture features of the tracer, and generate an initial validity mask for identifying invalid observation areas.

[0053] Due to the complex environment of ocean color remote sensing, the radiation signals of tracers such as chlorophyll concentration and total suspended matter are often submerged in strong atmospheric scattering, solar flares, and sea surface reflections, resulting in low signal-to-noise ratios in the original images. Directly inputting such images into optical flow algorithms can easily lead to false displacements or incorrect flow directions. Therefore, this embodiment first constructs an image preprocessing and local adaptive enhancement workflow based on quality control. The embodiment selects GOCI-II satellite L2 level hourly average products as the input data source, covering chlorophyll concentration (Chl). First, the pixel quality flags in the NetCDF file are read.

[0054] Step S1 includes the following sub-steps:

[0055] S1.1 Mask Generation Based on Quality Flag: To ensure the purity of the observation data, a strict mask generation logic is defined. Let the original flag value at time t and spatial coordinates (x, y) be... By constructing a mask All pixels affected by cloud cover, sea ice cover, land areas, or inversion anomalies are uniformly marked as invalid, as follows: ,

[0056] In the formula, Indicates the mask bits corresponding to cloud / cloud occlusion. Indicates the mask bits corresponding to sea ice. Indicates the mask bits corresponding to the land area. The mask bits represent the corresponding values ​​for quality failure, invalidity, or retrieval failure. This represents the original quality identifier value read at pixel (x, y) at time t. This indicates a bitwise AND operation.

[0057] Based on the above judgment logic, all pixels affected by cloud cover, sea ice coverage, land areas, or inversion anomalies can be uniformly marked as invalid.

[0058] S1.2 Local Normalization Processing: After removing obviously invalid pixels using a quality mask, large-scale brightness gradient problems still exist, such as overall brightness tilt caused by changes in the solar altitude angle, or regional radiation differences caused by changes in atmospheric conditions. If global enhancement is directly applied to the entire image, noise is often amplified simultaneously, while weak texture areas remain difficult to distinguish.

[0059] Therefore, this embodiment employs a sliding window-based local normalization method for adaptive image enhancement. The sliding window size is set to W×W pixels. In the code implementation of this embodiment, the parameter is set to 20, i.e., a 20×20 neighborhood window is used. For each valid pixel (x, y) in the image, the local mean of its neighborhood is calculated. and local standard deviation :

[0060] ,

[0061] ,

[0062] In the formula, the size of the sliding window is W×W, and k=(W-1) / 2. The number of valid pixels within the sliding window. Indicates the original image in pixels ( The pixel value (or grayscale value) at that location.

[0063] Then, a normalized image is obtained by performing a standardization mapping. As shown in the following formula:

[0064] ,

[0065] In the formula, Indicates the original image at the center pixel The pixel value (or grayscale value) at that location. It is a very small constant to prevent the denominator from being zero.

[0066] Through the above processing, the brightness value of each pixel is recalibrated with reference to its local statistical characteristics, thereby effectively eliminating large-scale background gradients and significantly amplifying local subtle texture changes. This processing method enables ocean fronts and vortex edges, which were originally in the weak gradient region, to obtain weights on a numerical scale comparable to those in the strong gradient region, providing more stable input features for subsequent optical flow algorithms.

[0067] S1.3, Linear stretching based on percentile truncation: To adapt to the optical flow algorithm's requirement for the input data range (typically 8-bit unsigned integers in the range [0, 255]), linear stretching based on percentile truncation is performed. Perform robust linear stretching. Calculate the cumulative distribution function (CDF) of all effective pixels in the image and determine the 2% quantile. and 98% quantile This strategy effectively eliminates extreme value interference caused by a very small number of noise points, resulting in a more enhanced image. The calculation formula is as follows:

[0068] ,

[0069] In the formula, This represents the lower bound of the normalization mapping. This represents the upper bound of the normalized mapping. clip Represents the cutoff / limiting function This indicates rounding to the nearest integer, converting a floating-point number to an integer.

[0070] This processing method effectively suppresses the impact of a very small number of outliers on the overall dynamic range, making the image grayscale distribution more concentrated and stable. After the above processing steps, the mesoscale eddies, frontal structures, and banded texture features in the original blurred ocean surface image become clearly distinguishable.

[0071] Step S2 involves calculating pixel displacement vectors using a motion estimation algorithm based on the enhanced images of adjacent time points (times t-1 and t) to obtain the initial flow field observations at a single time point. These initial flow field observations are then smoothed, and invalid regions are marked according to the initial validity mask. Specifically, this includes the following sub-steps:

[0072] S2.1 Initial Inversion: Utilizing the Farneback dense optical flow algorithm, this embodiment preferably employs the dense optical flow algorithm based on polynomial expansion proposed by Gunnar Farneback. Compared to sparse optical flow (such as Lucas-Kanade), this algorithm provides pixel-level dense vector fields, making it more suitable for representing continuous fluid motion. The algorithm approximates the signal intensity within the local neighborhood of the image using a quadratic polynomial:

[0073]

[0074]

[0075] In the formula, p= The spatial coordinate vector of the pixel. These are the x and y coordinates of the pixel, respectively. The quadratic coefficient matrix (symmetric matrix) describes the shape of the quadratic surface in the local neighborhood. , b1 and b2 represent the constant term (scalar), b1 and b2 represent the linear coefficient vector (row vector) describing the gradient direction, and T represents the transpose. d is the global displacement, satisfying: ,

[0076] The displacement of each pixel is obtained by solving the quadratic polynomial using the least squares method. In practice, in order to capture motion at different scales and prevent aliasing caused by large displacements, a multi-layered image pyramid strategy is adopted.

[0077] The specific parameter settings for this embodiment are as follows:

[0078] Pyramid levels: 5, allowing the algorithm to iterate gradually from coarse to fine scales, capturing motion with a maximum displacement not exceeding 1 / 8 of the image size.

[0079] Average window size (winsize): 50, which balances noise suppression and detail preservation.

[0080] Polynomial order (poly_n): 7, neighborhood fitting range.

[0081] Polynomial Gaussian standard deviation (poly_sigma): 2, used for the smoothness of the weighted fit.

[0082] S2.2 Multi-scale Gaussian smoothing and physical scale matching: The optical flow field obtained from the original inversion often contains a large amount of non-physical high-frequency noise, which mainly comes from: (a) sensor thermal noise; (b) sea surface wave flicker; (c) pseudo-motion at cloud edges. In addition, the pixel resolution of GOCI-II (about 250m) is much higher than that of the ocean dynamic background field (such as the HYCOM model, about 9km), and direct comparison lacks physical meaning.

[0083] In this embodiment, the horizontal (east-west) velocity component (u component) and vertical (north-south) velocity component (v component) obtained by inversion are smoothed by two-dimensional Gaussian convolution, as shown in the following formula:

[0084] ,

[0085] In the formula, The standard deviation of the Gaussian kernel is represented. This represents the weight value at spatial coordinates (x, y), used to perform a weighted average of the pixel to be smoothed and its neighborhood.

[0086] Calculate the smoothed flow field observations for:

[0087] ,

[0088] In the formula, This represents the original optical flow inversion result at pixel (x,y) in frame t.

[0089] S2.3 Outlier Removal: If the flow rate amplitude of a pixel... Greater than the maximum flow amplitude For example, a value of 3.0 m / s exceeds the general ocean current limit, so it is considered a divergence point in the algorithm, forcibly marked as an invalid value, and included in subsequent masking. middle.

[0090] Step S3 involves performing spatiotemporal alignment and unified rasterization on the initial flow field observations at multiple times, mapping them to a standard geospatial grid, constructing a multidimensional spatiotemporal flow field data cube containing spatial and temporal dimensions, and simultaneously integrating the validity masks of each time point to form a spatiotemporal observation state matrix.

[0091] S3.1 In this embodiment, the data is uniformly mapped onto a standard geographic reference grid. Specifically, a spatial index based on a k-dimensional tree is constructed for the observation data at each time step. For each target grid point on the standard grid, the radius R is determined according to the target grid resolution, and effective initial velocity observation points within this radius R are found. Spatial aggregation is performed using the inverse distance weighting (IDW) method, with spatial weights... And the initial values ​​of the horizontal and vertical velocity components assigned to the standard grid points after aggregation. The calculation is as follows:

[0092] ,

[0093] ,

[0094] in, Let represent the Euclidean distance from the i-th valid observation point in the neighborhood to the center of the target grid, and m represent the total number of valid observation points participating in the spatial aggregation calculation within the search radius (or local neighborhood). This represents a small positive number introduced to prevent the denominator from being zero; p is the distance decay exponent. , For the velocity components at the effective observation points, For the corresponding spatial weights, This is the initial flow velocity value assigned to the standard grid point after aggregation.

[0095] S3.2 Constructing the spatiotemporal flow field data cube and mask matrix: After spatial grid unification and aggregation, the aggregation results of N consecutive time steps are overlaid along the time dimension to form data with dimensions of... The horizontal (east-west) velocity component sequences u_maps and the vertical (north-south) velocity component sequences v_maps are generated, where H and W are the spatial dimensions of a unified standard grid. The effective generation dimension of the velocity components at each time step is also considered. The observation validity mask matrix is ​​used to determine whether, for any time t, a grid point has insufficient valid observation samples within the search radius. If such a grid point is not represented by enough valid observation samples, it will remain a missing value in u_maps and v_maps and will be marked as invalid in the mask. This serves as the basis for skipping observation updates and performing continuous reconstruction in the subsequent bidirectional temporal reconstruction process.

[0096] Step S4 involves reconstructing the flow field data based on a bidirectional temporal filtering method. This is the core step of the invention. Figure 2 As shown, to address the issues of observation gaps and temporal discrepancies caused by factors such as cloud cover, a linear dynamic system is constructed for the initial flow field observation sequence, and a bidirectional temporal recursion and fusion approach is adopted to achieve spatiotemporal continuous reconstruction of the flow field.

[0097] (1) State-space model and observation reliability assessment

[0098] For the flow field data cube constructed in step S3, the state vector of the system at time t is defined as:

[0099] ,

[0100] in, ) represent the east-west and north-south velocity components at position (x,y) at time t, respectively, and T represents the transpose.

[0101] To suppress outlier vectors, this embodiment introduces a dynamic observation uncertainty adjustment mechanism based on the observation amplitude, which is achieved by constructing the observation noise covariance matrix R. t To quantify the reliability of observations, this mechanism increases the observation noise weight R in low-velocity regions or anomalously low-value points. t This reduces the Kalman gain during the update process, thereby automatically suppressing interference from spurious vectors at the system level.

[0102] ,

[0103] in, Let be the observed flow velocity amplitude at time t. Represents the reference observation noise covariance. As a reference speed threshold, To prevent division by zero of tiny positive numbers, This represents the state transition matrix.

[0104] (2) Forward recursion: Recursively estimate from front to back according to time sequence, use the estimation results of the previous time to predict the current time, and update the results by combining the valid observations of the current time. Reduce the observation update weight or skip the update for invalid observation areas to keep the estimation results in time sequence continuity.

[0105] Specifically, the process is executed frame by frame from t=1 to N, and evolution prediction is performed based on the state transition matrix I and the process noise covariance Q.

[0106] The prediction step involves calculating the prior state estimate as follows:

[0107] = ,

[0108] ,

[0109] in, The prior state estimation vector at time t is composed of the east-west flow velocity u and the north-south flow velocity v. This represents the posterior state estimate vector at the previous time t-1; The prior error covariance matrix at time t is used to quantify the uncertainty of the predicted value. Let represent the posterior error covariance matrix of the previous time step t-1; It represents the process noise covariance matrix, characterizing the physical randomness of the flow field evolution over time;

[0110] Based on mask M t Update posterior state: If M t=1, meaning a valid observation, then the Kalman gain is calculated as follows, and then a state update is performed;

[0111] ,

[0112] in, The Kalman gain at time t determines the weight of the observation's contribution to the state update. The observation noise covariance matrix at time t is dynamically calculated based on the observed flow velocity amplitude.

[0113] If M t =0, meaning the observation is invalid or obscured by clouds, then the observation update is skipped or the observation weight is reduced accordingly, causing the state to degenerate into prediction / inertial propagation;

[0114] Update steps:

[0115] The state update formula is as follows:

[0116] ,

[0117] in, This represents the posterior state estimate vector at time t. This represents the prior state estimate vector at time t. This represents the Kalman gain at time t. This represents the observation vector at the current time t. This represents the observation residual.

[0118] The uncertainty update formula is as follows:

[0119] ,

[0120] in, This represents the posterior error covariance at time t. Let represent the prior error covariance at time t.

[0121] Spatial adaptive repair: After each time-instance update, adaptive smoothing is performed on the flow field. The distance d from the pixel to the reliable observation area is calculated using distance transformation. Based on the distance, weights are constructed, and multi-scale weighted fusion is performed on the Gaussian smoothing results of at least two different scales to eliminate isolated anomalous jump points and improve spatial continuity.

[0122] Forward recursion output: After completing the state update (or skip update) and spatial adaptive repair at time t, the forward state estimate and its uncertainty at that time are obtained and saved / recorded to form a forward estimation sequence.

[0123] (3) Backward recursion: Recursively recursively recursively from back to front, using the estimation results of the next time moment to predict the current time moment, and updating it in combination with the effective observations of the current time moment.

[0124] Specifically, the process is executed frame by frame in reverse from t=N to 1, and evolution prediction is made based on the state transition matrix I and the process noise covariance Q. However, the prediction is based on the posterior estimate of the next time step, so as to use the information of the subsequent time step to correct the estimation bias that may exist in the previous time step.

[0125] The prediction step involves calculating the prior state estimate as follows:

[0126] ,

[0127] ,

[0128] in, This represents the prior state estimate vector at time t, which is derived from the posterior estimate at the next time t+1. Backpropagation yielded, This represents the prior error covariance at time t. This represents the posterior error covariance at time t+1. Represents the process noise covariance matrix;

[0129] Based on mask M t =Update posterior state: If M t =1, meaning a valid observation, then calculate the Kalman gain and perform a state update;

[0130] If M t =0, meaning the observation is invalid or obscured by clouds, then the observation update is skipped or the observation weight is reduced accordingly, causing the state to degenerate into prediction / inertial propagation, thereby avoiding the introduction of abnormal displacement by invalid observations;

[0131] Update steps:

[0132] The state update formula is as follows:

[0133] ,

[0134] in, This represents the posterior state estimate vector at time t. This represents the prior state estimate vector at time t. This represents the observation vector at the current time t. This represents the observation residual.

[0135] The uncertainty update formula is as follows:

[0136] ,

[0137] in, This represents the posterior error covariance at time t. Let represent the prior state estimation vector at time t.

[0138] Spatial adaptive repair: After each time-instance update, adaptive smoothing is performed on the flow field. The distance d from the pixel to the reliable observation area is calculated using distance transformation. Based on the distance, weights are constructed, and multi-scale weighted fusion is performed on the Gaussian smoothing results of at least two different scales to eliminate isolated anomalous jump points and improve spatial continuity.

[0139] Backward recursion output: After completing the state update or skip update and spatial adaptive repair at time t, the backward state estimate and its uncertainty at that time are obtained and saved / recorded to form a backward estimation sequence.

[0140] (4) Two-way fusion: The state estimate and uncertainty at time t obtained by forward recursion and the state estimate and uncertainty at time t obtained by backward recursion are weighted and fused according to the uncertainty information to obtain the two-way reconstructed flow field:

[0141] ,

[0142] In the formula, This represents the fusion weights of the forward estimation. This represents the uncertainty at time t during the forward recursion process. This represents the uncertainty at time t during the backward recursion process. This represents a tiny positive number introduced to ensure the stability of numerical calculations.

[0143] The final reconstructed ocean surface flow field state vector at time t The calculation is as follows:

[0144] ,

[0145] in, This represents the final reconstructed ocean surface flow field state vector at time t, containing both meridional and latitudinal velocity components. This represents the forward recursive state estimation vector at time t, obtained through forward recursion. The backward recursive state estimation vector at time t is obtained by backward recursion, utilizing observation information from subsequent times.

[0146] Step S5: Although the flow field reconstructed in step S4 is continuous in time, it may contain non-physical noise (such as local strong-radial scattering points) introduced by independent pixel calculations in space. To ensure physical rationality, this step introduces divergence constraints based on physical consistency to perform final spatial dimension refinement of the flow field.

[0147] like Figure 3 As shown, based on the principle of physical consistency of fluid motion, divergence constraints are applied to the temporal reconstructed flow field. The divergence gradient is calculated iteratively and the velocity vector is corrected so that the flow field divergence tends to meet the continuity requirements. The result is a spatiotemporally continuous sea surface flow field inversion result after correction by physical consistency constraints.

[0148] (1) Explanation of constraint principle and triggering of iterative loop

[0149] Based on the approximately incompressible physical property of seawater, the horizontal divergence of the sea surface current field in a two-dimensional plane should approach zero. Its mathematical definition is:

[0150]

[0151] in, This represents the two-dimensional horizontal divergence of the sea surface current field. These represent the partial derivatives of the velocity components in their respective spatial coordinate directions.

[0152] Before entering the constraint process, set the maximum number of iterations N and the constraint strength coefficient. The system initializes the iteration counter i=1 and constructs an iterative correction closed loop. The system continuously determines whether the current counter i satisfies the condition. If satisfied, a single divergence gradient descent correction loop is triggered; if satisfied... If the error occurs, the correction will be terminated and the process will jump to the product packaging output.

[0153] (2) Calculation of divergence field and its spatial gradient based on central difference

[0154] Within a single cycle, the system uses the central difference method to calculate the sum of the gradients of each component of the flow field at the current moment, thus obtaining the discrete divergence field.

[0155]

[0156] in, Let represent the discrete horizontal divergence value at spatial coordinates (x, y) during the i-th iteration. Let x and y represent the east-west (meridian) and north-south (latitudinal) velocity components at position (x, y) during the i-th iteration, respectively. These represent the spatial step size (grid spacing) of the grid in the longitudinal and latitudinal directions, respectively.

[0157] Subsequently, spatial gradient calculations were performed on the divergence field again to extract gradient components that characterize the drasticness of divergence changes and their spatial orientation. and :

[0158]

[0159]

[0160] in, and These represent the spatial gradient components of the current divergence field in the meridional and latitudinal directions, respectively, pointing in the direction of the fastest increase in divergence value.

[0161] (3) Physical correction and iterative state update of flow velocity vector

[0162] The aforementioned spatial gradient is used as a physical penalty term, combined with the constraint strength coefficient. Perform a reverse approximation correction on the current flow velocity vector:

[0163]

[0164]

[0165] in, This indicates the velocity vector component used in the next iteration or final output after this correction. This represents the step size strength coefficient for divergence constraint correction.

[0166] Repeat the above steps until the preset maximum number of iterations is reached.

[0167] (4) Final product output

[0168] The spatiotemporally continuous flow field, corrected by the aforementioned physical consistency constraints, is the final high-precision ocean surface flow field product output. The product is saved as a grid containing latitude and longitude coordinates, corrected velocity components u and v, and flow direction, ensuring the continuity and physical rationality of the product across spatiotemporal scales.

[0169] Experimental Results and Quantitative Analysis

[0170] To fully verify the effectiveness and robustness of the method of this invention, this embodiment selects GOCI-II image data from six consecutive observation times (with a time interval of 1 hour) in a typical sea area as the experimental input set. The experiment uses the flow field state at the preceding time as prior information to perform bidirectional Kalman filtering spatiotemporal reconstruction on the continuous time period, and compares the inversion results with the HYCOM global ocean mixed coordinate model data of the same period.

[0171] (1) Time evolution analysis of accuracy index

[0172] like Figure 4 As shown, the trends of flow direction error and flow velocity error at these six consecutive time points are statistically analyzed.

[0173] Flow direction accuracy: The mean absolute error (MAE) of the initial flow field fluctuates significantly, with an average value of around 22.6°. However, the method of this invention (solid line) significantly reduces the flow direction error and stabilizes it at around 16° through spatiotemporal constraints, demonstrating excellent direction correction capability.

[0174] Flow rate accuracy: The root mean square error (RMSE) curve of the flow rate of the method of the present invention is always below that of the prior art and the trend is more stable, indicating that the algorithm has good noise resistance and robustness.

[0175] (2) Comprehensive performance evaluation of multi-dimensional indicators

[0176] To further quantitatively evaluate the performance of the inverted flow field in terms of vector decomposition accuracy and spatial structure similarity, this embodiment statistically analyzes the average values ​​of various indicators at the aforementioned six consecutive time points. In addition to the conventional flow direction and velocity errors, the root mean square error of the velocity component (RMSE-U, RMSE-V) and the velocity scalar correlation coefficient are also introduced for comprehensive evaluation. See Table 1 below:

[0177] Table 1: Comparison of Comprehensive Performance Indicators of Initial Flow Field and Reconstructed Flow Field

[0178] Evaluation indicators Initial flow field Reconstructing the flow field MAE_dir(°) 22.58 16.59 RMSE_speed (m / s) 0.227 0.184 RMSE-U (m / s) 0.216 0.165 RMSE-V (m / s) 0.198 0.145 Corr_U 0.748 0.773 Corr_V 0.766 0.828

[0179] Based on the statistical results, the method of this invention was used to compare and evaluate the "initial flow field" and "reconstructed flow field" within the same study area and time period. The results show that the reconstructed flow field achieves consistent improvement in both flow direction and velocity. Specifically, the mean absolute error of flow direction (MAE_dir) decreased from 22.58° to 16.59°, a reduction of approximately 26.5%, indicating that in cloud-occupied and weakly textured regions, temporal reconstruction can significantly suppress flow direction deviations introduced by missing measurements and noise. The root mean square error of velocity (RMSE_speed) decreased from 0.227 m / s to 0.184 m / s, a reduction of approximately 18.9%, demonstrating that the reconstruction process has a smoothing and corrective effect on random fluctuations in velocity amplitude. Furthermore, the errors of both the east-west and north-south components decreased simultaneously: RMSE_U decreased from 0.216 m / s to 0.165 m / s (a decrease of approximately 23.6%), and RMSE_V decreased from 0.198 m / s to 0.145 m / s (a decrease of approximately 26.8%). This indicates that the present invention not only improves the overall velocity magnitude but also enhances the quantitative accuracy and directional stability of the vector components. Simultaneously, the correlation indices were further enhanced: Corr_U increased from 0.748 to 0.773, and Corr_V increased from 0.766 to 0.828, indicating that the reconstructed flow field and the reference background field have higher consistency in spatial distribution and variation trends, especially in the V component, demonstrating the present invention's ability to recover vector field structure information under multi-temporal constraints.

[0180] While the present invention has been described above with reference to preferred embodiments, it is not intended to limit the invention. Those skilled in the art can make various modifications and refinements without departing from the spirit and scope of the invention. Therefore, the scope of protection of the present invention shall be determined by the claims.

Claims

1. A method for inverting and reconstructing sea surface current fields based on high-frequency geostationary satellite imagery, characterized in that, Includes the following steps: Step S1: Acquire geostationary satellite image data of the target sea area at multiple consecutive times, preprocess the image data to enhance the texture features of the tracer, and generate an initial validity mask for identifying invalid observation areas; Step S2: Based on the enhanced images of adjacent time moments, a motion estimation algorithm is used to calculate the pixel displacement vector, obtain the initial flow field observation value at a single time moment, and smooth the initial flow field observation value. Invalid regions are marked according to the initial validity mask. Step S3: Spatiotemporally align and uniformly rasterize the initial flow field observations at multiple times, map them to a standard geospatial grid, construct a multidimensional spatiotemporal flow field data cube containing spatial and temporal dimensions, and simultaneously integrate the validity masks of each time to form a spatiotemporal observation state matrix. Step S4, reconstructing the flow field data based on a bidirectional temporal filtering method, including: Forward recursion: Recursively recursively from front to back in time sequence, using the estimation results of the previous time moment to predict the current time moment, and updating it in combination with the valid observations of the current time moment. For invalid observation areas, the observation update weight is reduced or the update is skipped, so that the estimation results maintain the temporal continuity. Backward recursion: Recursively recursively recursively from back to front, using the estimation results of the next time step to predict the current time step, and updating it by combining the effective observations of the current time step; Two-way fusion: The forward recursive estimation results and the backward recursive estimation results are weighted and fused according to their respective uncertainties to obtain the time-series reconstructed flow field; Step S5: Based on the principle of physical consistency of fluid motion, a divergence constraint is applied to the temporal reconstructed flow field. The divergence gradient is calculated iteratively and the velocity vector is corrected so that the flow field divergence tends to meet the continuity requirement. The result of spatiotemporally continuous sea surface flow field inversion after correction by physical consistency constraint is output.

2. The method for sea surface current field inversion and reconstruction based on high-frequency geostationary satellite imagery according to claim 1, characterized in that, Step S1 includes the following sub-steps: S1.1, Generating a mask based on quality identifiers: By constructing a mask All pixels affected by cloud cover, sea ice cover, land areas, or inversion anomalies are uniformly marked as invalid, as follows: ， In the formula, Indicates the mask bits corresponding to cloud / cloud occlusion. Indicates the mask bits corresponding to sea ice. Indicates the mask bits corresponding to the land area. The mask bits indicating quality failure / invalidity / retrieval failure This represents the original quality identifier value read at pixel (x, y) at time t. Indicates a bitwise AND operation; S1.2 Local Normalization Processing: For each pixel (x, y) in the image, calculate the local mean of its neighborhood. and local standard deviation As shown in the following formula: ， ， In the formula, the size of the sliding window is W×W, and k=(W-1) / 2. The number of valid pixels within the sliding window. Indicates the original image in pixels ( The pixel value or grayscale value at (). Then, a normalized image is obtained by performing a standardization mapping. As shown in the following formula: ， In the formula, Indicates the original image at the center pixel Pixel value at that location, It is a very small constant to prevent the denominator from being zero; S1.3, Linear stretching based on percentile truncation: for normalized images Linear stretching is performed to remove extreme value interference caused by noise, resulting in an enhanced image. The calculation formula is as follows: ， In the formula, This represents the lower bound of the normalization mapping. This represents the upper bound of the normalized mapping. clip Represents the cutoff / limiting function This indicates rounding to the nearest integer, converting a floating-point number to an integer.

3. The method for sea surface current field inversion and reconstruction based on high-frequency geostationary satellite imagery according to claim 1, characterized in that, Step S2 includes the following sub-steps: S2.1 Using the Farneback dense optical flow algorithm, the signal intensity in the local neighborhood of the image is approximated by a quadratic polynomial, as follows: ， ， In the formula, p= The spatial coordinate vector of the pixel. This represents the coefficient matrix of the quadratic terms, describing the shape of the quadratic surface in the local neighborhood. , b1 and b2 represent the constant term, the local brightness offset, the coefficient vector of the first-order term, describing the gradient direction, and T represents the transpose. d is the global displacement, satisfying: , The displacement of each pixel is obtained by solving the quadratic polynomial using the least squares method. It employs a multi-level pyramid strategy to capture displacements at different scales, with parameter settings including: number of pyramid levels, average window size, polynomial order, and polynomial Gaussian standard deviation. S2.2 Apply two-dimensional Gaussian convolution smoothing to the horizontal velocity component u and the vertical velocity component v respectively, and define the Gaussian kernel function: ， In the formula, The standard deviation of the Gaussian kernel is represented. This represents the weight value at spatial coordinates (x, y), used to perform a weighted average of the cell to be smoothed and its neighborhood; Calculate the smoothed flow field observations for: ， In the formula, This represents the original optical flow inversion result at pixel (x,y) in frame t; S2.3 Outlier Removal: If the flow rate amplitude of a pixel... Greater than the maximum flow amplitude If it is not a divergence point in the algorithm, it will be forcibly marked as an invalid value and included in subsequent masks. middle.

4. The method for sea surface current field inversion and reconstruction based on high-frequency geostationary satellite imagery according to claim 1, characterized in that, Step S3 includes the following sub-steps: S3.1 Based on the initial flow field data of N consecutive time moments output in step S2, construct a spatial index based on kd-tree for the observation data at each time moment. For each target grid point on the standard grid, determine the radius R according to the target grid resolution, and find the effective initial flow velocity observation points within the radius R. Use the inverse distance weighting method for spatial aggregation. S3.2 After spatial grid unification aggregation, the aggregation results of N consecutive time steps are overlaid along the time dimension to form dimensions of... The east-west velocity component sequences u_maps and north-south velocity component sequences v_maps are generated, where H and W are the spatial dimensions of a unified standard grid. The effective generation dimension of the velocity components at each time step is also considered. The observation validity mask matrix is ​​used to determine the grid point. For any time t, if a grid point has no preset number of valid observation samples within the search radius, then the grid point remains a missing value in u_maps and v_maps and is marked as invalid in the mask. This serves as the basis for skipping observation updates and performing continuous reconstruction in the subsequent bidirectional time series reconstruction process.

5. The method for inverting and reconstructing sea surface current fields based on high-frequency geostationary satellite imagery according to claim 1, characterized in that, In step S4, the forward recursion serves as the initial temporal constraint, executed frame by frame from t=1 to N, and evolution prediction is performed based on the state transition matrix I and the process noise covariance Q. The prediction step involves calculating the prior state estimate as follows: ＝ ， ， in, The prior state estimation vector at time t is composed of the east-west flow velocity u and the north-south flow velocity v. This represents the posterior state estimate vector at the previous time t-1; The prior error covariance matrix at time t is used to quantify the uncertainty of the predicted value. Let represent the posterior error covariance matrix of the previous time step t-1; It represents the process noise covariance matrix, characterizing the physical randomness of the flow field evolution over time; According to the mask M t Update posterior state: if M t = 1, i.e. valid observation, compute Kalman gain: ， in, The Kalman gain at time t determines the weight of the observation's contribution to the state update. The observation noise covariance matrix at time t is dynamically calculated based on the observed flow velocity amplitude. Where speed is the observed flow velocity amplitude at time t; r0 represents the baseline observation noise covariance, and speed0 is the reference velocity threshold. To prevent division by zero of tiny positive numbers, Represents the state transition matrix; If M t =0, meaning the observation is invalid or obscured by clouds, then the observation update is skipped or the observation weight is reduced accordingly, causing the state to degenerate into prediction / inertial propagation; Update steps: The state update formula is as follows: ， in, This represents the posterior state estimate vector at time t. Let represent the prior state estimation vector at time t. This represents the Kalman gain at time t. This represents the observation vector at the current time t. Represents the observation residual; The uncertainty update formula is as follows: ， in, This represents the posterior error covariance at time t. This represents the prior error covariance at time t; Spatial adaptive repair: After each time-instance update, adaptive smoothing is performed on the flow field. The distance d from the pixel to the reliable observation area is calculated using distance transformation. Based on the distance, weights are constructed, and multi-scale weighted fusion is performed on the Gaussian smoothing results of at least two different scales to eliminate isolated anomalous jump points and improve spatial continuity. Forward recursion output: After completing the state update and spatial adaptive repair at time t, the forward state estimate and its uncertainty at that time are obtained and saved / recorded to form a forward estimation sequence.

6. The method for inverting and reconstructing sea surface current fields based on high-frequency geostationary satellite imagery according to claim 5, characterized in that, In step S4, the backward recursion is a reverse time constraint, which is executed frame by frame from t=N to 1. Evolution prediction is performed based on the state transition matrix I and the process noise covariance Q. The prediction is based on the posterior estimate of the next time step, and the estimation bias in the previous time step is corrected by using the information of the subsequent time step. The prediction step involves calculating the prior state estimate as follows: ， ， in, This represents the prior state estimate vector at time t, which is derived from the posterior estimate at the next time t+1. Backpropagation yielded, This represents the prior error covariance at time t. This represents the posterior error covariance at time t+1. Represents the process noise covariance matrix; Based on mask M t Update posterior state: If M t =1, meaning a valid observation, then calculate the Kalman gain and perform a state update; If M t =0, meaning the observation is invalid or obscured by clouds, then the observation update is skipped or the observation weight is reduced accordingly, causing the state to degenerate into prediction / inertial propagation, thus avoiding the abnormal displacement introduced by invalid observations; Update steps: The state update formula is as follows: ， in, This represents the posterior state estimate vector at time t. This represents the observation vector at the current time t; The uncertainty update formula is as follows: ， in, This represents the posterior error covariance at time t. Represents the prior state estimation vector at time t; Spatial adaptive repair: After each time-instance update, adaptive smoothing is performed on the flow field. The distance d from the pixel to the reliable observation area is calculated using distance transformation. Based on the distance, weights are constructed, and multi-scale weighted fusion is performed on the Gaussian smoothing results of at least two different scales to eliminate isolated anomalous jump points and improve spatial continuity. Backward recursion output: After completing the state update or skip update and spatial adaptive repair at time t, the backward state estimate and its uncertainty at that time are obtained and saved / recorded to form a backward estimation sequence.

7. The method for sea surface current field inversion and reconstruction based on high-frequency geostationary satellite imagery according to claim 1, characterized in that, In step S4, the bidirectional fusion specifically involves: weighting and fusing the state estimates and uncertainties at time t obtained from forward recursion and the state estimates and uncertainties at time t obtained from backward recursion according to the uncertainty information to obtain the bidirectional reconstructed flow field. ， in, This represents the fusion weights of the forward estimation. This represents the uncertainty at time t during the forward recursion process. This represents the uncertainty at time t during the backward recursion process. This represents a tiny positive number introduced to ensure the stability of numerical calculations; ， in, This represents the final reconstructed ocean surface flow field state vector at time t, containing east-west and north-south velocity components. The forward recursive state estimation vector at time t is obtained by forward recursion. Let represent the backward recursive state estimation vector at time t, obtained by backward recursion.

8. The method for inverting and reconstructing sea surface current fields based on high-frequency geostationary satellite imagery according to claim 1, characterized in that, Step S5 includes the following sub-steps: S5.1 Set the maximum number of iterations and constraint strength coefficient; S5.2 Calculate the discrete divergence field of the flow field at the current moment using the central difference method; S5.3 Calculate the spatial gradient of the discrete divergence field; S5.

4. The spatial gradient is used as a physical penalty term, and combined with the constraint strength coefficient, the current flow velocity vector is corrected. S5.5 Repeat steps S5.2 to S5.4 until the preset maximum number of iterations is reached.

9. An electronic device comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that, When the processor executes the computer program, it implements the steps of the method according to any one of claims 1 to 8.

10. A computer-readable storage medium having a computer program stored thereon, characterized in that, When the computer program is executed by a processor, it implements the steps of the method according to any one of claims 1 to 8.

Images