Method and system for simulating ocean phenomenon processes based on particle systems

Abstract

The application discloses a marine phenomenon process simulation method and system based on a particle system, and relates to the technical field of marine numerical simulation, which comprises the following steps: acquiring multi-source marine environment data, solving optimal weights to generate an equivalent stiffness field based on sparse regression and energy consistency functional minimization, calculating a local coherence time based on a logarithmic growth rate of a two-dimensional optimal transport distance of a particle group, and identifying a coherence dominant mode and an inertia dominant mode based on local linear stability analysis, constructing a particle velocity correction equation based on a local dissipation response approximation, wherein the particle velocity correction equation comprises a memory Stokes drift velocity and an effective transport velocity of the particles, constructing an energy balance determination rule to judge a gathering state of the particles based on an energy balance relationship between a coherent field work and relative kinetic energy, updating particle positions according to the effective transport velocity and the gathering state, and outputting a final particle transport trajectory for presenting a marine phenomenon process.

Description

Technical Field

[0001] This invention relates to the technical field of marine numerical simulation, and in particular to a method and system for simulating marine phenomena and processes based on particle systems. Background Technology

[0002] In the fields of marine pollution tracking, oil spill diffusion prediction, marine plastic waste enrichment analysis, and phytoplankton red tide evolution simulation, accurate prediction of the movement trajectory and aggregation pattern of marine surface materials is crucial. Traditional marine phenomenon simulation methods mainly rely on the Lagrange particle system model and introduce coherent structure theory to describe the aggregation and fragmentation effects of eddies and fronts on particles in the ocean, thereby macroscopically presenting specific marine phenomena. However, existing particle system-based simulation techniques have significant drawbacks when dealing with complex ocean dynamics: First, when calculating the equivalent stiffness or diffusion coefficient driving particle motion, multiple empirical weighting coefficients (such as the weights of wave gradient, density gradient, and convergence terms) need to be manually set. These artificial coefficients cannot adapt to the dynamic environment of different sea areas and spatiotemporal variations. Second, when defining the coherence time (memory time) that controls particle aggregation lifetime, it usually relies on the finite-time Lyapunov exponent FTLE, which requires setting a fixed integration time window and lacks physical consistency. Third, when identifying the switching between coherent and inertial dominant modes, as well as the triggering of particle aggregation and fragmentation, existing methods mostly use hard empirical thresholds, which makes particle systems prone to non-physical abrupt changes under complex sea conditions, resulting in poor generalization ability and an inability to truly reproduce the evolution of ocean phenomena.

[0003] Currently, Chinese invention patent CN118797896A discloses a method for constructing a digital marine environment simulation system, including: Step 1, collecting marine environmental data and constructing a basic marine scene; Step 2, collecting marine topographic data, constructing a seabed land model and adding it to the basic marine scene; Step 3, collecting submersible data and marine organism data, constructing a submersible model and a marine organism model respectively, and adding them to the basic marine scene; Step 4, simulating formation control, pursuit and combat, and moving wake waves of the submersible based on its motion attitude, and adding it to the submersible model; Step 5, simulating the motion patterns and behaviors of marine organisms, as well as the impact of water currents on marine organisms, and adding it to the marine organism model, thus completing the construction of the marine environment simulation system. This invention can accurately simulate physical phenomena and biological behaviors in the ocean, achieving a more realistic and detailed simulation of the marine environment. However, existing methods for simulating ocean phenomena typically employ globally fixed empirical diffusion coefficients, which are prone to over-diffusion at sub-mesoscale fronts, failing to accurately reproduce the banded enrichment of matter. The FTLE method, which relies on a fixed integration time window, is prone to ghost vortices and time lag artifacts, making it unable to accurately capture short-lived transient flow field structures. Ignoring the inertial hysteresis of macroscopic floating objects leads to severe inertial overshoot trajectory errors under complex wave conditions. Using geometric distance thresholds lacking physical basis to forcibly determine particle aggregation and fragmentation results in non-physical rigid aggregation or spurious fragmentation in the simulation results. Summary of the Invention

[0004] The technical problem solved by this invention is that existing particle system-based methods for simulating ocean phenomena typically use globally fixed empirical diffusion coefficients, which easily lead to over-diffusion at sub-mesoscale fronts and fail to accurately reproduce the banded enrichment of matter; the FTLE method, which relies on a fixed integration time window, is prone to producing ghost vortices and time lag artifacts, and cannot accurately capture short-lived transient flow field structures; ignoring the inertial hysteresis of macroscopic floating objects leads to severe inertial overshoot trajectory errors under complex wave conditions; and using geometric distance thresholds lacking physical basis to forcibly determine particle aggregation and fragmentation results in non-physical rigid aggregation or spurious fragmentation in the simulation results.

[0005] To solve the above-mentioned technical problems, the present invention provides the following technical solution: a method for simulating ocean phenomena based on particle systems, comprising the following steps:

[0006] Step S1: Acquire multi-source marine environmental data, and solve for the optimal weights based on sparse regression and energy consistency functional minimization to generate an equivalent stiffness field;

[0007] Step S2: Based on the logarithmic growth rate of the optimal transport distance of the particle cluster in two dimensions, calculate the local coherence time, and identify the coherence-dominant mode and the inertia-dominant mode based on local linear stability analysis.

[0008] Step S3: Construct the particle velocity correction equation based on the local dissipation response approximation. The particle velocity correction equation includes the particle's memory Stokes drift velocity and effective transport velocity.

[0009] Step S4: Based on the energy balance relationship between the work done by the coherent field and the relative kinetic energy, construct an energy balance judgment rule to determine the aggregation state of the particles, update the particle position according to the effective transport velocity and aggregation state, and output the final particle transport trajectory to present the ocean phenomenon process.

[0010] Preferably, step S1 includes the following sub-steps:

[0011] Step S11: Obtain multi-source marine environmental data of the target sea area. The multi-source marine environmental data includes surface flow velocity, seawater density, significant wave height, wave period data and coastline geometry data. Based on the surface flow velocity, calculate the mass derivative of the fluid micro-particles.

[0012] Step S12: Release virtual particles in the computational domain, use the lead particle tracking algorithm to perform trajectory integration on the virtual particles within a preset short time window, count the number of virtual particles on each grid node in the computational domain, and calculate the preliminary particle concentration field.

[0013] Step S13: Extract the original physical features that affect ocean particle transport. The original physical features include wave gradient features, density gradient features and Stokes convergence features. Introduce a unified feature force density scale parameter to perform dimensionless processing on the original physical features and perform standardization processing on columns to obtain a standardized feature matrix.

[0014] Step S14: Define the equivalent stiffness field as a linear combination of the features in each column of the standardized feature matrix, construct an energy-consistent functional, and solve for the optimal weights and generate the equivalent stiffness field based on sparse regression.

[0015] Preferably, step S12 includes the following sub-steps:

[0016] Step S121: Obtain the two-dimensional spatial range of the target sea area, discretize the two-dimensional spatial range into an Eulerian grid, and obtain the set of spatial coordinates of all grid nodes;

[0017] Virtual particles are uniformly released within the computational domain at a preset initial spatial density. Each virtual particle is assigned a unique index number, and the initial spatial coordinates and initial mass of each virtual particle are recorded to obtain a set of virtual particles.

[0018] Step S122: Set a preset short time window for pilot tracking. Based on the surface flow field velocity, use the fourth-order Runge-Kutta numerical integration algorithm to integrate the motion trajectory of all virtual particles within the preset short time window and obtain the end spatial coordinates of each virtual particle at the end of the integration.

[0019] Step S123: Traverse all grid nodes of the Euler grid, and use the smoothing kernel function to map the virtual particle mass at the end of the integration to the adjacent grid nodes to calculate the preliminary particle concentration field on each grid node.

[0020] Preferably, step S14 includes:

[0021] The mathematical expression for the equivalent stiffness field is:

[0022] ;

[0023] in, For the equivalent stiffness field, Let be the non-negative weight vector to be determined. , and Non-negative weights , and These are the dimensionless features in the standardized feature matrix;

[0024] Construct an energy consistency functional, the mathematical expression of which is:

[0025] ;

[0026] in, It is an energy-consistent functional. For the computational domain, The momentum residual term represents the squared error between the particle concentration diffusion flux divergence driven by the equivalent stiffness field and the actual inertial force of the fluid element. For L1 regularization terms, The regularization coefficient is . Let L1 be the norm of the weight vector;

[0027] Solving for optimal weights and generating an equivalent stiffness field based on sparse regression specifically includes:

[0028] The minimization problem of the energy consistency functional is transformed into a convex optimization problem with non-negativity constraints, and its mathematical expression is:

[0029] ;

[0030] ;

[0031] in, The optimal non-negative weight vector. Function representation of taking and Minimum The value, , and These are the optimal non-negative weights that need to be substituted into the equivalent stiffness field formula;

[0032] The above equations are solved using a convex optimization algorithm to obtain the optimal weight vector. The optimal weight vector is then substituted into the equivalent stiffness field formula to generate the final equivalent stiffness field.

[0033] Preferably, step S2 includes the following sub-steps:

[0034] Step S21: Initial position of the computational domain A first particle cluster and a second particle cluster are initialized at a certain point. The first particle cluster and the second particle cluster are spatially adjacent. The first particle cluster and the second particle cluster are subjected to time evolution based on the surface flow field velocity to obtain the probability measure distribution of the first particle cluster and the second particle cluster at time t.

[0035] Step S22: Calculate the two-dimensional optimal transport distance between the first particle cluster and the second particle cluster during the evolution process, and calculate the logarithmic growth rate of the two-dimensional optimal transport distance within a local time window;

[0036] Step S23: Calculate the inverse scale of the logarithmic growth rate and define it as the local coherence time. Construct the coherent tensor evolution equation. Linearize the coherent tensor evolution equation in the local equilibrium neighborhood to obtain a local linear operator. Perform eigenvalue decomposition on the local linear operator and extract the largest eigenvalue with the largest real part among all eigenvalues.

[0037] Step S24: When the real part of the largest eigenvalue is less than zero, it indicates that the disturbance decays over time and the structure is stable, and the current region is determined to be a coherent dominant mode;

[0038] When the real part of the largest eigenvalue is greater than or equal to zero, it indicates that the disturbance increases over time and the structure is easily broken, and the current region is determined to be in an inertial-dominated mode.

[0039] Preferably, step S23 specifically includes:

[0040] The mathematical expression for the local coherence time is:

[0041] ;

[0042] in, Initial position Local coherence time at that point This is the preset minimum growth rate lower limit. The growth rate is logarithmic.

[0043] A coherent tensor evolution equation is constructed, and then linearized in the local equilibrium neighborhood to obtain a linear perturbation equation:

[0044] ;

[0045] in, For partial derivative operators, This represents the perturbation of the coherence tensor. To be in the initial position The local linear operator obtained after freezing the background field coefficients;

[0046] For the local linear operator Perform eigenvalue decomposition, extract all eigenvalues, and obtain the largest eigenvalue with the largest real part. .

[0047] Preferably, step S3 includes the following sub-steps:

[0048] Step S31: Obtain the physical property parameters of the target particle and calculate the response time of the particle, wherein the response time is the ratio of the particle mass to the effective damping coefficient;

[0049] Step S32: Calculate the memory Stokes drift velocity of the particle based on the equivalent stiffness field and local coherence time, using the local dissipative response approximation theory.

[0050] Step S33: Construct a near-shore activation function based on spatial distance. Based on the pattern recognition results output in step S24, calculate the coherent induced velocity. Vector-add the surface flow velocity, the memory Stokes drift velocity, and the coherent induced velocity modulated by the near-shore activation function to obtain the effective transport velocity of the particles.

[0051] Preferably, step S4 includes the following sub-steps:

[0052] Step S41: Traverse any two adjacent first and second particles within the computational domain, and calculate the reduced mass and relative kinetic energy of the first and second particles;

[0053] Step S42: Perform line integration along the spatial line connecting the first particle and the second particle to calculate the work done by the coherent field on the first particle and the second particle, and obtain the work done by the coherent field.

[0054] Step S43: Obtain the local dissipated energy of the current sea area grid, add the relative kinetic energy and the local dissipated energy to obtain the separation energy threshold, determine the aggregation state of the first and second particles according to the energy balance judgment rule, update the particle position according to the effective transport velocity and aggregation state, and output the final particle transport trajectory to present the ocean phenomenon process.

[0055] Preferably, step S43 specifically includes:

[0056] The local dissipated energy is calculated using the equivalent stiffness field and particle spacing, and its mathematical expression is as follows:

[0057] ;

[0058] in, This represents the energy dissipation barrier caused by environmental turbulent diffusion on particle aggregation. This represents the equivalent stiffness field at the midpoint of the spatial connection. Spatial coordinates of the first particle Spatial coordinates of the second particle The Euclidean distance;

[0059] The relative kinetic energy is added to the local dissipated energy to obtain the separation energy threshold;

[0060] The specific rules for determining energy balance include:

[0061] When the work done by the coherent field is strictly greater than the separation energy threshold, it is determined that the first particle and the second particle have aggregated.

[0062] When the work done by the coherent field is not greater than the separation energy threshold, it is determined that the first particle and the second particle remain in a dispersed state. If the first particle and the second particle were originally in an aggregated state in the previous time step, it is determined that a split has occurred at the current time.

[0063] Based on the effective transport velocity, the spatial coordinates of all particles are updated in time step using a numerical integration algorithm. During the update process, coordinate merging and mass accumulation operations are performed on particle pairs that are determined to have aggregated, and coordinate separation and mass splitting operations are performed on particle clusters that are determined to have fragmented.

[0064] The simulation is executed in a loop with time increments until the preset total simulation duration is reached, and the final particle transport trajectory data is output to represent the ocean phenomenon process.

[0065] A particle-based ocean phenomenon process simulation system, used to execute a particle-based ocean phenomenon process simulation method, includes a basic field reconstruction module, a pattern recognition module, a dynamic velocity correction module, and an energy balance evolution module.

[0066] The basic field reconstruction module is used to acquire multi-source marine environmental data and solve for the optimal weights to generate an equivalent stiffness field based on sparse regression and energy consistency functional minimization.

[0067] The pattern recognition module is used to calculate the local coherence time based on the logarithmic growth rate of the two-dimensional optimal transport distance of the particle cluster, and to identify the coherent dominant mode and the inertial dominant mode based on local linear stability analysis.

[0068] The dynamic velocity correction module is used to construct a particle velocity correction equation based on the local dissipation response approximation. The particle velocity correction equation includes the particle's memory Stokes drift velocity and effective transport velocity.

[0069] The energy balance evolution module is used to construct energy balance judgment rules based on the energy balance relationship between the work done by the coherent field and the relative kinetic energy to determine the aggregation state of particles, update the particle position according to the effective transport velocity and aggregation state, and output the final particle transport trajectory to present the process of ocean phenomena.

[0070] The beneficial effects of this invention are as follows: By constructing an energy-consistent functional and solving it using sparse regression (LASSO), this invention achieves variational self-consistent estimation of the equivalent stiffness field. It completely eliminates the practice of manually setting gradient weights such as waves and density in traditional particle system simulations, enabling the model to automatically learn the dominant physical mechanisms from actual ocean dynamic residuals, significantly improving the model's adaptability and physical interpretability in different sea areas.

[0071] This invention introduces a two-dimensional optimal transport (OT) distance to define local coherence time and utilizes local linear stability analysis (calculating the maximum eigenvalue of the local linear operator) to identify the dominant coherent mode. This improvement eliminates the dependence of the traditional FTLE method on a fixed time window and the non-physical abrupt changes caused by the hard threshold switching of the Weissenberg number, making the identification of coherent structures more rigorous, smooth, and consistent with hydrodynamic stability theory.

[0072] This invention, during the particle dynamics evolution stage, approximates the memory Stokes drift velocity based on local dissipative response and constructs aggregation criteria based on the balance between work done by the coherent field and relative kinetic energy. This physically sufficient criterion based on energy and work replaces the traditional free energy difference threshold and effective temperature parameter, which lack physical basis, enabling spontaneous and self-consistent triggering of particle aggregation and fragmentation. This significantly improves the accuracy and robustness of predicting the evolution of marine phenomena (such as oil spill aggregation and red tide spread) under extreme sea conditions. Attached Figure Description

[0073] Figure 1 This is a basic flowchart of a particle system-based ocean phenomenon process simulation method provided in one embodiment of the present invention. Detailed Implementation

[0074] To make the above-mentioned objects, features and advantages of the present invention more apparent and understandable, the specific embodiments of the present invention will be described in detail below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments.

[0075] Reference Figure 1 As an embodiment of the present invention, a method for simulating ocean phenomena based on particle systems is provided, comprising the following steps:

[0076] Step S1: Acquire multi-source marine environmental data, and solve for the optimal weights based on sparse regression and energy consistency functional minimization to generate an equivalent stiffness field;

[0077] Step S2: Based on the logarithmic growth rate of the optimal transport distance of the particle cluster in two dimensions, calculate the local coherence time, and identify the coherence-dominant mode and the inertia-dominant mode based on local linear stability analysis.

[0078] Step S3: Construct the particle velocity correction equation based on the local dissipation response approximation. The particle velocity correction equation includes the particle's memory Stokes drift velocity and effective transport velocity.

[0079] Step S4: Based on the energy balance relationship between the work done by the coherent field and the relative kinetic energy, construct an energy balance judgment rule to determine the aggregation state of the particles, update the particle position according to the effective transport velocity and aggregation state, and output the final particle transport trajectory to present the ocean phenomenon process.

[0080] This application constructs a data-driven particle tracking framework that eliminates manual thresholding. This framework achieves physical consistency through four steps: variational minimization allows data to automatically learn equivalent stiffness; optimal transport distance (OT) and linear stability analysis spontaneously extract coherent time and patterns from the flow field topology; it abandons fixed time windows; derives velocity correction based on local dissipative response; solves the inertial overshoot problem of macroscopic floating objects; and reverts to classical mechanics, using work and energy balance to determine particle aggregation and fragmentation. This method establishes the adaptive, high-fidelity, zero-empirical-threshold physical thresholding framework of this invention in the field of ocean numerical simulation.

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

[0082] Step S11: Obtain multi-source marine environmental data for the target sea area. The multi-source marine environmental data includes surface flow velocity, seawater density, significant wave height, wave period data, and coastline geometry data. Based on the surface flow velocity, calculate the mass derivative of the fluid micro-particles.

[0083] Step S12: Release virtual particles in the computational domain, use the lead particle tracking algorithm to perform trajectory integration on the virtual particles within a preset short time window, count the number of virtual particles on each grid node in the computational domain, and calculate the preliminary particle concentration field.

[0084] Step S13: Extract the original physical features that affect ocean particle transport. The original physical features include wave gradient features, density gradient features and Stokes convergence features. Introduce a unified feature force density scale parameter to perform dimensionless processing on the original physical features, and perform standardization processing on the columns to obtain a standardized feature matrix.

[0085] Step S14: Define the equivalent stiffness field as a linear combination of the features in each column of the standardized feature matrix, construct an energy-consistent functional, and solve for the optimal weights and generate the equivalent stiffness field based on sparse regression.

[0086] In this embodiment, the mathematical expression for calculating the mass derivative of the fluid micro-particle based on the surface flow field velocity in step S11 is as follows:

[0087] ;

[0088] in, The mass derivative of the surface flow velocity. Let t be the surface flow velocity vector, and t be time. For local acceleration, For convective acceleration, This is the gradient operator.

[0089] In step S13, based on multi-source marine environmental data, the original physical features affecting marine particle transport are extracted. These original physical features include wave gradient features, density gradient features, and Stokes convergence features.

[0090] The mathematical expression for wave gradient characteristics is:

[0091] ;

[0092] The mathematical expression for the density gradient feature is:

[0093] ;

[0094] The mathematical expression for the Stokes convergence characteristic is:

[0095] ;

[0096] in, , and These are wave gradient characteristics, density gradient characteristics, and Stokes convergence characteristics, respectively. Let g be the wave-induced Stokes stress tensor, and g be the acceleration due to gravity. The density of seawater, For wave group speed, For free-running Stokes drift speed, The norm of a matrix or vector. For divergence operators;

[0097] Introducing a unified characteristic force density scale parameter, the mathematical expression for dimensionless processing of the original physical characteristics is as follows:

[0098] ;

[0099] in, For the k-th dimensionless feature, For the k-th original physical feature, For a unified characteristic force density scale parameter;

[0100] The dimensionless features are standardized column-wise so that the mean of each feature is zero and the variance is one, resulting in a standardized feature matrix. ,in , and This refers to the dimensionless characteristics after standardization.

[0101] It should be noted that after feature extraction, all physical quantities must be dimensionless and standardized. The logical connection of this step is that the original physical dimensions of wave gradient, density gradient, and convergence term are completely different. If regression calculation is performed directly, the optimization algorithm will be dominated by the feature with the largest absolute value. Dimensionless and standardized processing ensures that all features compete on the same mathematical scale, guaranteeing the numerical stability of subsequent regressions.

[0102] It should be noted that when simulating the enrichment of marine oil spills or microplastics into sub-mesoscale fronts, traditional Lagrange particle models typically use empirical random walks to simulate turbulent diffusion. Since the diffusion coefficient is often globally fixed or based on large-scale empirical formulas, the algorithm produces strong over-diffusion at the front. That is, particles that should be captured by the strong density gradient of the front and form elongated bands are forcibly dispersed by the excessively large random diffusion term, resulting in blurred boundaries and severely undervalued concentrations in the simulated pollution zone.

[0103] Step S1 of this method completely abandons the empirical diffusion coefficient by constructing an energy-consistent functional and utilizing LASSO sparse regression. When particles approach the front, the algorithm automatically identifies that the density gradient dominates the local residuals, thus consistently assigning a very high weight to the density gradient and dynamically generating an anisotropic equivalent stiffness field. This allows the random diffusion of particles in the front normal direction to be automatically suppressed by the physical mechanism, perfectly reproducing the banded enrichment phenomenon of matter at the front. Existing technologies require manual subjective setting of the weights of wave gradients, density gradients, and other factors when calculating equivalent stiffness. This application requires that the diffusion driven by the equivalent stiffness field must be physically rigorously explainable the actual motion residuals of the leading particles. By constructing an energy-consistent functional and introducing L1 regularization, the system can automatically realize the sparse selection of physical mechanisms, achieving complete data-driven and physically consistent results.

[0104] Step S12 includes the following sub-steps:

[0105] Step S121: Obtain the two-dimensional spatial range of the target sea area, discretize the two-dimensional spatial range into an Eulerian grid, and obtain the set of spatial coordinates of all grid nodes;

[0106] Virtual particles are uniformly released within the computational domain at a preset initial spatial density. Each virtual particle is assigned a unique index number, and the initial spatial coordinates and initial mass of each virtual particle are recorded to obtain a set of virtual particles.

[0107] Step S122: Set a preset short time window for pilot tracking. Based on the surface flow field velocity, use the fourth-order Runge-Kutta numerical integration algorithm to integrate the motion trajectory of all virtual particles within the preset short time window and obtain the end spatial coordinates of each virtual particle at the end of the integration.

[0108] Step S123: Traverse all grid nodes of the Euler grid, and use the smoothing kernel function to map the virtual particle mass at the end of the integration to the adjacent grid nodes to calculate the preliminary particle concentration field on each grid node.

[0109] In this embodiment, the mathematical expression for the virtual particle set is:

[0110] ;

[0111] in, A collection of virtual particles. For the p-th virtual particle at the initial time... The initial spatial coordinates, Let N be the initial mass of the p-th virtual particle, and N be the total number of virtual particles released.

[0112] The mathematical expression for obtaining the final spatial coordinates of each virtual particle at the end of the integration is:

[0113] ;

[0114] in, Let be the spatial coordinates of the p-th virtual particle at the end of the integration. To preset a short time window, For a moment In spatial coordinates The surface flow velocity vector at that location;

[0115] The mathematical expression for the mapping process is:

[0116] ;

[0117] in, For grid nodes The initial particle concentration field at that location, For grid nodes The corresponding control element area, For Gaussian kernel function, Let the terminal spatial coordinates of the p-th virtual particle be the coordinates of the grid node. The Euclidean distance between them This is the smooth length parameter.

[0118] The preset initial spatial density is dynamically allocated based on the total mass of the target ocean phenomenon and the computer's memory and computing power limits, and is the ratio of the desired number of particles per grid cell to the area of ​​the control surface cell. To ensure the continuity and statistical significance of subsequent kernel function mapping, the system requires each grid control surface cell to... The initial content must contain at least a preset number (4 to 10) of virtual particles.

[0119] The initial mass of the virtual particles follows the law of conservation of mass. The total mass of the simulated marine material (such as the total tonnage of oil spills or the total biomass of algae) is obtained, and the total number of virtual particles generated in the release area based on the initial spatial density is obtained. The initial mass of a single particle is the ratio of the total mass to the total number of particles.

[0120] A short time window is preset for pilot tracking to obtain the initial concentration field. The time must be long enough to show the background convergence trend of the flow field, but short enough to prevent particles from crossing multiple grids, which could lead to over-smoothing or distortion of the concentration field. The time window should be a small integer multiple of the system simulation time step, with a reasonable range of 10 minutes to 6 hours.

[0121] The area of ​​the control element is calculated based on the Earth's spherical projection coordinate system, using the latitude and longitude step size of the grid nodes to represent the actual physical area. The mathematical expression is as follows:

[0122] ;

[0123] in, For the Earth's radius, The latitude of the current node. and These are the grid radian steps for longitude and latitude, respectively.

[0124] The smoothing length parameter, based on the Smooth Particle Hydrodynamics (SPH) theory, must cover a sufficient number of neighboring particles to ensure the smoothness of the kernel density estimation. Its mathematical expression is:

[0125] ;

[0126] in, The smoothing coefficient is dynamically determined by the system using cross-validation. Values ​​to minimize the spatial variance of the reconstructed concentration field. and It refers to the spatial spacing of a two-dimensional grid in the x and y directions, which is the physical distance between adjacent grid nodes; its value range is usually 1.2 to 2.0 times the grid resolution.

[0127] It should be noted that the core purpose of this step is to provide a real physical observation target for subsequent variational self-consistent optimization. In the energy consistency functional of step S14, it is necessary to calculate the particle concentration diffusion flux divergence driven by the equivalent stiffness field. To solve this functional, the system must first know what the actual aggregation trend of particles is like in a very short time under the pure flow field advection. Step S12 is achieved by executing a very short time window. The leader trace is performed and a smooth kernel function is used. Transforming discrete Lagrange particle coordinates into a continuous Eulerian concentration field This preliminary concentration field It directly reflects the inherent convergence characteristics of the flow field itself, thus serving as the reference data input into the subsequent LASSO sparse regression. This forces the system to automatically learn the equivalent stiffness field weights that can perfectly explain this concentration change, ensuring the data-driven and physical self-consistency of the entire algorithm chain.

[0128] Step S14 includes:

[0129] The mathematical expression for the equivalent stiffness field is:

[0130] ;

[0131] in, For the equivalent stiffness field, Let be the non-negative weight vector to be determined. , and Non-negative weights , and These are the dimensionless features in the standardized feature matrix;

[0132] Construct an energy consistency functional, the mathematical expression of which is:

[0133] ;

[0134] in, It is an energy-consistent functional. For the computational domain, The momentum residual term represents the squared error between the particle concentration diffusion flux divergence driven by the equivalent stiffness field and the actual inertial force of the fluid element. For L1 regularization terms, The regularization coefficient is . It is the L1 norm of the weight vector, used to make the weight vector sparse;

[0135] Solving for optimal weights and generating an equivalent stiffness field based on sparse regression specifically includes:

[0136] The minimization problem of the energy-consistent functional is transformed into a convex optimization problem with non-negativity constraints, and its mathematical expression is:

[0137] ;

[0138] ;

[0139] in, The optimal non-negative weight vector. Function representation of taking and Minimum The value, , and These are the optimal non-negative weights that need to be substituted into the equivalent stiffness field formula;

[0140] The above equations are solved using a convex optimization algorithm to obtain the optimal weight vector. The optimal weight vector is then substituted into the equivalent stiffness field formula to generate the final equivalent stiffness field.

[0141] The specific process of solving the equation using the coordinate descent method with non-negativity constraints as a convex optimization algorithm includes:

[0142] The computational domain is discretized into a spatial set containing M grid nodes, and each grid node is obtained. The spatial coordinates of the location, where ;

[0143] The computation logic for M grid nodes is as follows: the original data grid based on the input ocean circulation model (ROMS, HYCOM) is spatially discretized, the number of longitude nodes and latitude nodes in the computation domain are extracted, the invalid nodes corresponding to the land mask are removed, and the total number of remaining valid water nodes is the product of the number of longitude nodes and latitude nodes and the difference between them.

[0144] In each grid node At each location, the particle diffusion flux divergence driven by three normalized dimensionless features is calculated, and the feature matrix is ​​constructed. Its mathematical expression is:

[0145] ;

[0146] in, Characteristic matrix The element in the i-th row and k-th column, These correspond to wave gradient characteristics, density gradient characteristics, and Stokes convergence characteristics, respectively. For grid nodes The k-th standardized dimensionless feature at position k, For grid nodes The initial particle concentration field at the location;

[0147] In each grid node At this point, the actual inertial force density of the fluid micro-element is calculated, and the target vector is constructed. Its mathematical expression is:

[0148] ;

[0149] in, target vector The i-th element in For grid nodes Seawater density at that location For grid nodes The mass derivative of the surface flow velocity at a given location.

[0150] The continuous energy consistency functional is transformed into a discrete LASSO regression objective function, the mathematical expression of which is:

[0151] ;

[0152] in, Let M be the discrete optimization objective function, and M be the total number of grid nodes. Let k be the weight coefficient to be determined. The regularization coefficient is used.

[0153] The regularization coefficient is the core penalty parameter controlling the sparsity of LASSO regression. Before solving the equivalent stiffness field, the system uses the L-curve method to calculate different... Using the sum of squared residuals and the L1 norm of the weights, we find the inflection point with the largest curvature (Pareto optimal solution) as the best solution for the current grid. The value is used to balance the fitting accuracy with the sparsity of the physical mechanism.

[0154] Physical constraints are set: Since the equivalent stiffness field must be positive to ensure the physical stability of the diffusion process, this is a mandatory requirement. .

[0155] Iterative updates based on coordinate descent include:

[0156] Initialize weight vector Set the current iteration number n=0;

[0157] For the k-th weight coefficient, the local bias residual vector after excluding this feature is calculated, and its mathematical expression is:

[0158] ;

[0159] in, For the grid nodes in the nth iteration The local partial residual vector at that location, This represents the current value of the j-th weight coefficient in the n-th iteration;

[0160] The mathematical expression for calculating the unregularized relevance and normalization factor of the k-th feature is as follows:

[0161] ;

[0162] ;

[0163] in, For unregularized correlation, Normalization factor;

[0164] The soft threshold operator, combining nonnegativity constraints and L1 regularization, updates the k-th weight coefficient, and its mathematical expression is:

[0165] ;

[0166] in, For the updated k-th weight coefficient, Used to strictly enforce nonnegativity constraints and achieve sparse truncation of features;

[0167] The above update formula is executed sequentially in the order of k=1, 2, 3 to complete one round of coordinate descent iteration, and then n=n+1.

[0168] Calculate the L2 norm of the change in the weight vector between two consecutive iterations. When the L2 norm of the change is strictly less than the preset convergence tolerance, the iteration is considered converged, the calculation is stopped, and the optimal weight vector is output. Substituting the optimal weight vector into the equivalent stiffness field formula, the final equivalent stiffness field is generated:

[0169] .

[0170] The preset convergence tolerance is the stopping criterion for the coordinate descent method. It is set based on the precision of computer floating-point numbers (such as the truncation error of Float64) and the relative error allowed in engineering. The L2 norm of the weight vector change between two adjacent iterations is also considered. When the value is less than the convergence tolerance, the objective function is considered to have reached its minimum point.

[0171] This step details how to translate abstract variational minimization into a computer-executable numerical algorithm. The core innovation lies in employing a coordinate descent method with non-negativity constraints. In ocean dynamics, a negative equivalent stiffness (diffusion coefficient) can lead to non-physical reverse diffusion, causing instantaneous computational crashes. This step cleverly combines L1 regularized soft thresholding with non-negativity physical constraints in the iterative formula. When a physical mechanism (such as wave gradient) contributes to the interpretation of local residuals... Unable to overcome regularization penalty At that time, the operator will directly assign the corresponding weights. Forced zeroing not only ensures extremely high computational efficiency in the solution process, but also enables automatic filtering of weak and irrelevant dynamic noise in the ocean at the mathematical level, accurately extracting the core physical mechanism that dominates particle transport in the current sea area.

[0172] Step S2 includes the following sub-steps:

[0173] Step S21: Initial position of the computational domain The first and second particle clusters are initialized at a certain point. The first and second particle clusters are spatially adjacent. The first and second particle clusters are subjected to time evolution based on the surface flow field velocity to obtain the probability measure distribution of the first and second particle clusters at time t.

[0174] In this embodiment, step S21 includes the following sub-steps:

[0175] Step S211: Select the initial spatial coordinates of the target sea area within the computational domain as the first anchor point. Define a small spatial perturbation vector The first anchor point is translated along the spatial perturbation vector to obtain the second anchor point. This ensures that the first anchor point and the second anchor point are strictly adjacent in space;

[0176] tiny spatial perturbation vector Used to initialize reference particle clusters at adjacent locations to calculate OT distances, its modulus is [not specified]. The infinitesimal assumption of the first-order Taylor expansion must be satisfied. The computational logic is to set its modulus to be strictly equal to a minimal proportion of the current grid space step size, i.e. The proportionality coefficient .

[0177] Centered on the first anchor point, the number of generators generated according to an isotropic Gaussian distribution is: Discrete particles, that is, randomly sampled particles in a two-dimensional plane according to a Gaussian distribution with the same variance and the same direction. Using the coordinates of each particle, we construct the initial state set for the first particle cluster:

[0178] ;

[0179] in, For the first particle cluster at the initial moment The set of spatial coordinates, Let be the spatial coordinates of the i-th particle in the first particle cluster at the initial moment. The total number of particles contained in a single particle cluster;

[0180] Gaussian distribution random sampling The coordinates of each particle are used to construct an empirical probability measure. The quantity is determined by the convergence of the Monte Carlo integral. Due to statistical error and Proportional, based on the maximum permissible statistical error (e.g., 5%) The required number of particles can be calculated by working backwards. .

[0181] Centered on the second anchor point, the number of generators generated according to the same isotropic Gaussian distribution is... The discrete particles are used to construct the initial state set of the second particle cluster:

[0182] ;

[0183] in, For the second particle cluster at the initial moment The set of spatial coordinates, Let be the spatial coordinates of the i-th particle in the second particle cluster at the initial moment.

[0184] Step S212: Extract the surface flow field velocity obtained in step S1 Using the fourth-order Runge-Kutta numerical integration algorithm, kinematic ordinary differential equations are established for each particle in the first and second particle clusters, and synchronous trajectory integration is performed in the time interval [0,t]. The mathematical expression is as follows:

[0185] ;

[0186] in, For time derivative operators, They represent the first and second particle clusters, respectively. For the first The i-th particle in the particle cluster at time i spatial coordinates, For a moment In spatial coordinates The surface flow velocity vector at that location;

[0187] Integrate to the target time t, obtaining the sets of spatial coordinates of the first and second particle clusters after their evolution at time t, denoted as follows: and .

[0188] Step S213: Transform the evolved discrete particle coordinates into a mathematical probability measure space. For the first particle cluster, assign each particle an equal mass weight. The empirical probability measure distribution of the first particle cluster at time t is constructed using the Dirac function, and its mathematical expression is:

[0189] ;

[0190] in, Let be the probability measure distribution of the first particle cluster at time t. Centered on spatial coordinates Dirac function at;

[0191] Similarly, for the second particle cluster, construct the empirical probability measure distribution of the second particle cluster at time t:

[0192] ;

[0193] in, Let be the probability measure distribution of the second particle cluster at time t. Centered on spatial coordinates The Dirac function at that location.

[0194] This step details how to prepare rigorous mathematical inputs for subsequent two-dimensional optimal transport (OT) distance calculations. Traditional finite-time Lyapunov exponential FTLE calculations only track the distance stretching rate between two single particles (point masses). However, ocean currents exhibit strong multi-scale turbulence, and the trajectory of a single particle is easily disturbed by localized micro-eddies, leading to violent oscillations in the calculated coherence time and a lack of statistical representativeness.

[0195] The innovative logic of step S21 lies in upgrading from single-particle tracking to particle cluster (probabilistic measure) tracking. This is achieved by initializing two clusters containing a large number of particles at adjacent anchor points. The Gaussian particle cluster, and its evolution result is strictly defined as an empirical probability measure composed of the Dirac function ( and The system can capture the overall statistical characteristics of fluid micro-elements during deformation, shearing, and stretching processes. This allows for the application of optimal transport theory (calculating the distribution) in step S22. Transport and remodeling into distribution This lays a solid data foundation at the lowest cost, fundamentally improving the noise resistance and physical robustness of marine coherent structure identification.

[0196] Step S22: Calculate the two-dimensional optimal transport distance between the first and second particle clusters during the evolution process, and calculate the logarithmic growth rate of the two-dimensional optimal transport distance within a local time window;

[0197] The mathematical expression for the optimal two-dimensional transport distance is:

[0198] ;

[0199] in, Let time t be the initial position The optimal two-dimensional transport distance at a given location. It is the optimal two-dimensional input operator. and These are the probability measure distributions of the first and second particle clusters at time t, respectively.

[0200] Within a local time window, the logarithmic growth rate of the optimal two-dimensional transport distance is calculated, and its mathematical expression is as follows:

[0201] ;

[0202] in, The growth rate is logarithmic. The length of the local time window. For time derivative, To prevent the logarithmic function's independent variable from being zero, the first numerical regularization term is used.

[0203] The local time window length is used to calculate the average logarithmic growth rate of the optimal transport (OT) distance. Its calculation is based on the integral time scale of local turbulence, obtained by integrating the historical velocity time series using the autocorrelation function. The timescale is set to 1 to 3 times the integral timescale to filter high-frequency noise and preserve coherent evolution trends.

[0204] The first numerical regularization term is a purely mathematical anti-singularity term used to prevent a fatal error caused by the complete overlap of two particle clusters when calculating the logarithmic growth rate, which would lead to an undefined logarithmic function. Its value is set to computer machine precision.

[0205] Step S23: Calculate the inverse scale of the logarithmic growth rate and define it as the local coherence time. Construct the coherent tensor evolution equation. Linearize the coherent tensor evolution equation in the local equilibrium neighborhood to obtain the local linear operator. Perform eigenvalue decomposition on the local linear operator and extract the largest eigenvalue with the largest real part among all eigenvalues.

[0206] Step S24: When the real part of the largest eigenvalue is less than zero, it indicates that the disturbance decays over time and the structure is stable, and the current region is determined to be a coherent dominant mode;

[0207] When the real part of the largest eigenvalue is greater than or equal to zero, it indicates that the disturbance increases over time and the structure is prone to breakage, thus determining that the current region is inertial-dominated.

[0208] It should be noted that this step replaces the hard physical threshold with optimal transport and linear stability analysis. This application spontaneously extracts the coherent time from the actual dispersion rate of particle clusters by calculating the growth rate of the two-dimensional optimal transport distance; simultaneously, through rigorous local linear stability analysis, the system automatically distinguishes between coherent and inertial modes based on the dynamic stability of the flow field itself. This approach eliminates the non-physical abrupt changes caused by empirical thresholds, making the identification of coherent structures smoother and based on rigorous fluid dynamics.

[0209] When identifying the encapsulation and transport of plankton (such as red tides) by ocean eddies, traditional models rely heavily on the Finite-Time Lyapunov Exponent (FTLE) to calculate coherent structures. However, the FTLE requires a fixed integration time window (e.g., 15 days). If the actual lifetime of a high-frequency transient eddy is only 3 days, the forced integration over 15 days will impose future flow field changes onto the current moment, leading to the appearance of ghost eddies or elongated artifacts in the calculation results, which seriously mislead the prediction of red tide outbreak paths.

[0210] Step S2 of this method introduces the growth rate of the two-dimensional optimal transport (OT) distance and local linear stability analysis. No longer requiring any pre-defined integration time window, it determines in real-time, frame-by-frame whether the current fluid particle is coherently dominated (real part of eigenvalue < 0, structurally stable) or inertially dominated (real part of eigenvalue ≥ 0, structurally broken) by calculating the maximum eigenvalue of the local linear operator. This transient response mechanism completely eliminates time window artifacts and can accurately capture the birth and death processes of short-lived vortices.

[0211] Step S23 specifically includes:

[0212] The mathematical expression for the local coherence time is:

[0213] ;

[0214] in, Initial position Local coherence time at that point This is the preset minimum growth rate lower limit. The growth rate is logarithmic.

[0215] The preset minimum growth rate lower limit is used to prevent the local coherence time from tending to infinity (i.e., dividing by zero). The calculation logic is to obtain the physically permissible lifespan of the largest ocean mesoscale eddy (e.g., half a year), and use its reciprocal as the absolute lower limit of the growth rate.

[0216] A coherent tensor evolution equation is constructed, and then linearized in the neighborhood of the local equilibrium state to obtain a linear perturbation equation:

[0217] ;

[0218] in, For partial derivative operators, This represents the perturbation of the coherence tensor. To be in the initial position The local linear operator obtained after freezing the background field coefficients;

[0219] For local linear operators Perform eigenvalue decomposition, extract all eigenvalues, and obtain the largest eigenvalue with the largest real part. .

[0220] In this embodiment, constructing the coherent tensor evolution equation specifically includes:

[0221] The stretching, shearing, and folding effects of ocean currents on particle clusters are abstracted as a second-order symmetric positive definite tensor, defined as the coherence tensor. The coherence tensor is directly obtained by calculating the spatial covariance matrix of the particle cluster and performing dimensionless transformation. The specific process includes:

[0222] Step A1: Extract the set of spatial coordinates of the first particle cluster after time t evolution in step S212. The mathematical expression for calculating the centroid coordinates of the first particle cluster at time t is as follows: ;

[0223] in, Let be the coordinates of the center of mass of the first particle cluster at time t. This represents the total number of particles contained in the first particle cluster. Let be the spatial coordinates of the i-th particle in the first particle cluster at time t.

[0224] Step A2: Using the center of mass coordinates as the origin, calculate the relative displacement vector of each particle in the first particle cluster relative to the center of mass. Its mathematical expression is:

[0225] ;

[0226] in, Let be the relative displacement vector of the i-th particle at time t.

[0227] Step A3: Summate the tensor outer products of the relative displacement vectors of all particles to calculate the spatial covariance matrix of the first particle cluster. Its mathematical expression is:

[0228] ;

[0229] in, Let be the spatial covariance matrix of the first particle cluster at time t. It is the outer product operator (tensor product) for vectors. Generate a second-order symmetric positive semi-definite matrix;

[0230] Since the particle clusters have a non-degenerate distribution in space, the spatial covariance matrix is ​​strictly a second-order symmetric positive definite matrix.

[0231] Step A4: Obtain the initial spatial variance when initializing the isotropic Gaussian distribution particle cluster in step S211. The spatial covariance matrix is ​​dimensionless by using the initial spatial variance to obtain the coherent tensor corresponding to the first particle cluster, whose mathematical expression is:

[0232] ;

[0233] in, Initial position The coherence tensor of the fluid element at time t. This is the initial spatial variance scalar;

[0234] When t=0, since the initial particle cluster has an isotropic Gaussian distribution, its covariance matrix is... This strictly guarantees the initial coherence tensor. ( (where is a second-order unit tensor) that satisfies the initial conditions for local equilibrium set in step S242.

[0235] It should be noted that the coherence tensor is constructed based on the covariance matrix of the first particle cluster at the initial position of the target, where the second particle cluster is only used as a reference perturbation particle cluster in the two-dimensional optimal transport OT distance calculation to characterize the sensitivity of local distribution evolution.

[0236] This step details how to rigorously quantify the abstract effects of stretching, shearing, and folding mathematically.

[0237] In traditional continuum mechanics, describing the deformation of fluid micro-clusters usually relies on the extremely complex calculation of the Cauchygreen deformation tensor. This requires high-precision spatial differentiation and long-term integration of the velocity gradient of the flow field. In discrete ocean grid data, this is prone to truncation errors and causes the tensor to lose its positive definiteness (physical distortion).

[0238] This application utilizes the statistical covariance of microscopic particle swarms to directly and equivalently replace the deformation tensor of macroscopic fluid microclusters. By calculating the second moment of the particle swarm relative to the center of mass (i.e., the spatial covariance matrix), the system can perfectly capture the topological effects of the flow field on matter: stretching is manifested as an increase in the largest eigenvalue of the covariance matrix (the particle cluster becomes longer in one direction); compression and folding are manifested as a decrease in the smallest eigenvalue (the particle cluster becomes flattened in another direction); shearing and rotation are manifested as a deflection of the eigenvector direction of the covariance matrix.

[0239] By dividing by the initial variance Dimensionless transformation yields the coherent tensor. It perfectly aligns with the mathematical definition of a second-order symmetric positive definite tensor and naturally satisfies the upper convection evolution equation given in step S241. This purely data-driven statistical computation method completely avoids complex grid space differentiation, greatly improving the numerical stability of complex ocean dynamics simulations while ensuring physical rigor.

[0240] Based on fluid kinematics and structural relaxation theory, the upper convection derivative is introduced to construct a nonlinear partial differential equation describing the evolution of the coherence tensor with respect to the surface flow field. Its mathematical expression is:

[0241] ;

[0242] in, The local time derivative of the coherence tensor. Let V be the surface flow field velocity vector. For spatial gradient operators, For convection terms, Let be the surface flow field velocity gradient tensor. This is the transpose of the velocity gradient tensor. The local coherence time is calculated in step S23. It is a second-order unit tensor. The relaxation dissipation term is used to restore the coherent tensor to an isotropic state.

[0243] At a given initial position At this point, the isotropic dispersion state of the particle cluster when it is not continuously stretched by the flow field is defined as a local equilibrium state, that is, the basic state tensor is defined. ;

[0244] Within the neighborhood of the local equilibrium state, a small second-order tensor perturbation is introduced. The coherence tensor is decomposed into the sum of the fundamental state and the perturbation:

[0245] ;

[0246] in, Spatial coordinates The coherence tensor at time t This represents the corresponding coherent tensor perturbation.

[0247] Second-order tensor perturbation In local linear stability analysis, this is introduced as a small perturbation in the theoretical derivation. When numerically solving for the eigenvalues ​​of the linear operator, it is not necessary to assign them actual physical values; instead, the Frobenius norm of its initial matrix is ​​normalized to 1, i.e. The maximum eigenvalue is directly extracted by solving the eigenvalue equation of the operator matrix, thus obtaining pure growth rate information. In theoretical derivation, it is an infinitesimal, and in numerical calculations, its norm is normalized to 1.

[0248] Transform the system to a Lagrangian servo system that translates with the fluid element and eliminate the convection term. The influence of local time derivatives Equivalent to the material derivative , recorded as ;

[0249] At the initial position At this point, the spatially dependent background field coefficients are frozen, and the velocity gradient tensor and local coherence time at that point are extracted as constant matrices and constant scalars:

[0250] ;

[0251] ;

[0252] in, To freeze the background field velocity gradient tensor, To freeze the local coherence time;

[0253] Substitute the decomposed coherent tensor into the coherent tensor evolution equation, and ignore second-order and higher-order infinitesimal terms containing perturbation products. The homogeneous linearized expansion is obtained as follows:

[0254] ;

[0255] Separate from the homogeneous linearized expansion the expression only related to the perturbation. The relevant homogeneous evolution part constructs a local linear operator acting on the second-order tensor space. Its mathematical expression is:

[0256] ;

[0257] in, To be in the initial position The local linear operator obtained after freezing the background field coefficients;

[0258] Substituting the local linear operator back into the expansion, we finally obtain the linear perturbation equation for the evolution of the control perturbation:

[0259] ;

[0260] By solving this linear perturbation equation, i.e., for the operator By performing eigenvalue decomposition, we can obtain the largest eigenvalue that determines the generation and destruction of local coherent structures.

[0261] This step details how to rigorously extract local linear operators for pattern recognition from complex nonlinear fluid dynamics equations.

[0262] In existing techniques for identifying coherent ocean structures (such as vortex boundaries and fronts), a dimensionless number, the Weissenberg number, is typically used. As a hard, empirical threshold, such as a mandatory provision For coherent clustering areas, This is the inertial dispersion region. This approach is extremely crude because it only considers the absolute magnitude (norm) of the velocity gradient, while completely ignoring the topological structure of the flow field (such as the difference in tensor directionality between rotation and shear).

[0263] This step introduces a local linear stability analysis of the entire tensor space. This is achieved by constructing a coherent tensor evolution equation that includes the upper convection derivative, and then analyzing the local equilibrium state (isotropic state). By performing rigorous Taylor expansion and coefficient freezing near the vicinity, a local linear operator in the form of a fourth-order tensor mapping was derived. This operator not only includes the relaxation dissipation rate. It also fully preserves the velocity gradient tensor. All orientation and rotation information. Subsequently, only the real part of the operator's largest eigenvalue needs to be calculated to mathematically rigorously determine whether a small disturbance will be amplified by the flow field (forming a coherent structure) or dissipated and smoothed out (regressing inertial dispersion), thus completely eliminating... Non-physical classification errors caused by hard thresholding.

[0264] Step S3 includes the following sub-steps:

[0265] Step S31: Obtain the physical property parameters of the target particle and calculate the response time of the particle. The response time is the ratio of the particle mass to the effective damping coefficient.

[0266] Step S32: Calculate the memory Stokes drift velocity of the particle based on the equivalent stiffness field and local coherence time, using the local dissipative response approximation theory.

[0267] Step S33: Construct a near-shore activation function based on spatial distance. Based on the pattern recognition results output in step S24, calculate the coherent induced velocity. Vector-add the surface flow velocity, the memory Stokes drift velocity, and the coherent induced velocity modulated by the near-shore activation function to obtain the effective transport velocity of the particles.

[0268] In this embodiment, in step S31, the physical property parameters of the target particle are obtained, including the particle's equivalent radius and particle density. The mass of the particle is calculated based on the sphere volume formula, and its mathematical expression is:

[0269] ;

[0270] in, For particle mass, For the equivalent radius, Particle density;

[0271] Obtain the dynamic viscosity coefficient of seawater in the current sea area. The effective damping coefficient of a particle is calculated based on Stokes' law of drag, and its mathematical expression is as follows:

[0272] ;

[0273] in, The effective damping coefficient, The dynamic viscosity coefficient of seawater;

[0274] The response time of the particle is calculated as the ratio of the particle mass to the effective damping coefficient, and its mathematical expression is:

[0275] ;

[0276] in, The response time of the particle.

[0277] In step S32, the effective wave height and wave period data extracted in step S1 are obtained, and the free Stokes drift velocity vector of the sea surface is calculated based on linear wave theory. ;

[0278] The Frobenius norm of the coherent tensor is calculated based on the equivalent stiffness field and local coherence time. ;

[0279] Based on the local dissipative response approximation theory, the inertial hysteresis of particles is coupled with the coherent structure of the flow field to calculate the memory Stokes drift velocity of particles. The mathematical expression for the memory Stokes drift velocity of particles is:

[0280] ;

[0281] in, To memorize the speed of Stokes Drift, For free-running Stokes drift speed, The response time of the particle. For local coherence time, Let Frobenius norm be the coherence tensor. For the equivalent stiffness field, To prevent the second numerical regularization term from having a denominator of zero;

[0282] The second numerical regularization term is a purely mathematical anti-singularity term, used to prevent the equivalent stiffness field from causing errors when calculating the Stokes drift velocity. The error of a denominator of zero caused by a value of 0 (such as in perfectly uniform still water) is addressed. Its value is set as a minimal scale of the characteristic scale of the stiffness field. The range of values ​​is... ~ .

[0283] In step S33, the vertical distance from the target particle's current spatial coordinates to the nearest coastline is obtained. ;

[0284] A near-shore activation function based on spatial distance is constructed to smoothly attenuate coherence effects in shallow near-shore waters. Its mathematical expression is as follows:

[0285] ;

[0286] in, For nearshore activation functions, This represents the vertical distance to the coastline. These are preset reference nearshore scale parameters;

[0287] The preset reference nearshore scale parameter determines the spatial decay rate of the nearshore activation function. By extracting isobath data from the coastal topography, the average distance from the coastline to the point where the water depth reaches the critical wave-breaking depth (usually 10 meters) is calculated. This distance is used as the reference nearshore scale parameter to ensure that the coherent effect is smoothly shut off within the wave-breaking zone. A smaller value is used for steep coasts, and a larger value is used for gentle tidal flats.

[0288] Based on the pattern recognition results output in step S24, calculate the coherent induction velocity:

[0289] When the current region is determined to be the coherent dominant mode, the spatial divergence of the coherent tensor is calculated as the coherent induced velocity, and its mathematical expression is:

[0290] ;

[0291] in, The coherent induced velocity coefficient;

[0292] When the current region is determined to be in the inertial-dominated mode, the coherent induced velocity is set to zero vector. .

[0293] The coherent induced velocity coefficient, as a mobility coefficient for dimension transformation, converts the equivalent force density (divergence) of the coherent field into induced velocity. Its calculation logic is based on Stokes' law of drag in fluid dynamics, namely, velocity equals force divided by the damping coefficient.

[0294] ;

[0295] in The dynamic viscosity coefficient of seawater, This is the reference particle size for the current simulated material.

[0296] The effective transport velocity of the particle is obtained by vector summing the surface flow velocity, the memory Stokes drift velocity, and the coherent induced velocity modulated by the near-shore activation function. Its mathematical expression is:

[0297] ;

[0298] in, For effective transport speed, For the surface flow velocity, For nearshore activation functions, The coherent induced velocity is activated only in the region determined to be the coherent dominant mode in step S2, and takes the value of zero vector in the region of the inertial dominant mode.

[0299] This step details how to overcome the physical limitations of traditional ocean simulations that equate macroscopic floating objects with massless water particles, achieving a unified understanding of multi-scale particle dynamics. In traditional Lagrange particle tracking models (LPTs), the calculation of wave-driven particle propulsion typically directly incorporates the full Stokes drift. This is superimposed on the flow velocity. This assumption holds true for extremely small particles such as microplastics; however, for macroscopic objects with significant mass and volume (such as containers, large seaweed rafts, and oil spill clumps), due to their enormous inertia, when the wave field changes abruptly (such as a sudden change in wind direction causing wave reversal), the object cannot respond instantaneously to the wave's drag, inevitably resulting in a memory lag. Traditional models ignore this, leading to severe inertial overshoot errors in predicted trajectories.

[0300] This application uses a local dissipative response approximation equation to self-consistently adjust the weight of wave drift. In the correction equation of step S32, the denominator of the correction coefficient includes a dimensionless term. The ratio of particle response time to flow field coherence time.

[0301] When the simulated object is a tiny particle with extremely small mass, the response time is... As the denominator approaches 0, the memory of Stokes Drift perfectly degenerates into the free Stokes Drift.

[0302] When the simulated object is a macroscopic heavy object with extremely large mass, the response time is... Significantly increased. At this point, if the flow field is in a strongly coherent state ( Larger and stiffer (If the denominator is smaller), the denominator will be significantly greater than 1, thus automatically weakening the instantaneous drag of the wave and producing a damping hysteresis effect that conforms to the laws of physics.

[0303] In addition, the near-shore activation function introduced in step S33 It cleverly solves the numerical instability problem in areas with complex nearshore topography. When particles are extremely close to the coast (…), When it approaches 0, The coherent induced velocity approaches 0, automatically shutting off the coherent induced velocity. This effectively mitigates the common numerical particle grounding phenomenon in traditional algorithms, where particles are pushed ashore by abnormal coherent forces and cannot return. This greatly improves the computational robustness of the model at complex coastline boundaries.

[0304] This step achieves self-consistency in velocity correction through the approximation of local dissipative response. Starting from local momentum balance and dissipative response, this application derives a rigorous correction equation that includes particle response time, coherence time, and equivalent stiffness field. When the particle mass or coherence stiffness is extremely large, the memory correction effect is significantly enhanced; while when the particle mass is extremely small, the memory correction naturally disappears, degenerating into free Stokes drift. This method requires no artificial coefficients and achieves the unification and self-consistency of multi-scale particle dynamics.

[0305] Step S4 includes the following sub-steps:

[0306] Step S41: Traverse any two adjacent first and second particles within the computational domain, and calculate the reduced mass and relative kinetic energy of the first and second particles;

[0307] Step S42: Perform line integration along the spatial line connecting the first particle and the second particle to calculate the work done by the coherent field on the first particle and the second particle, and obtain the work done by the coherent field.

[0308] Step S43: Obtain the local dissipated energy of the current sea area grid, add the relative kinetic energy and the local dissipated energy to obtain the separation energy threshold, determine the aggregation state of the first and second particles according to the energy balance judgment rule, update the particle position according to the effective transport velocity and aggregation state, and output the final particle transport trajectory to present the ocean phenomenon process.

[0309] In this embodiment, in step S41, the reduced mass and relative kinetic energy of any two adjacent first and second particles are calculated within the computational domain. The mathematical expression for this is:

[0310] ;

[0311] ;

[0312] in, To convert the mass, The mass of the first particle. The mass of the second particle. As relative kinetic energy, Let be the velocity vector of the first particle. The velocity vector of the second particle. Let L2 be the relative velocity between the two.

[0313] In step S42, a line integral is performed along the spatial line connecting the first particle and the second particle to calculate the work done by the coherent field on the particle pair. The mathematical expression for this work is:

[0314] Obtain the spatial coordinates of the first particle Spatial coordinates of the second particle Construct a spatial integration path from the first particle to the second particle. ;

[0315] Extract the coherent tensor generated in step S2, and calculate the divergence field of the coherent tensor in space. The divergence field represents the equivalent force density generated by the coherent field;

[0316] Along the spatial integration path, the continuous line integral of the work done by the coherent field on the particle pair is calculated, and its theoretical mathematical expression is as follows:

[0317] ;

[0318] in, Doing work for the coherent field, Let be the integral path from the position of the first particle to the position of the second particle. Let be the equivalent volume of the particle pair. Let be the divergence of the coherence tensor, and represent the equivalent power density generated by the coherent field. This is the equivalent force vector of the coherent field. Let be the path element vector.

[0319] During the computer numerical solution stage, since the first particle and the second particle are within the small spatial search radius defined in step S41, a first-order numerical discretization approximation is performed on the continuous line integral, and the midpoint of the spatial connection between the first particle and the second particle is taken. As a reference point;

[0320] The coherent tensor divergence value at the reference point is extracted and its dot product with the relative displacement vector of the two particles is performed to obtain the work done by the discretized coherent field. The specific calculation expression is as follows:

[0321] ;

[0322] in, For reference point The coherent tensor divergence vector calculated at that point Let be the relative displacement vector between the two particles.

[0323] Step S43 specifically includes:

[0324] The local dissipated energy is calculated using the equivalent stiffness field and particle spacing, and its mathematical expression is as follows:

[0325] ;

[0326] in, This represents the energy dissipation barrier caused by environmental turbulent diffusion on particle aggregation. This represents the equivalent stiffness field at the midpoint of the spatial connection. Spatial coordinates of the first particle Spatial coordinates of the second particle The Euclidean distance;

[0327] The relative kinetic energy is added to the local dissipated energy to obtain the separation energy threshold;

[0328] The specific rules for determining energy balance include:

[0329] When the work done by the coherent field is strictly greater than the separation energy threshold, it is determined that the first particle and the second particle have aggregated.

[0330] When the work done by the coherent field is not greater than the separation energy threshold, it is determined that the first particle and the second particle remain in a dispersed state. If the first particle and the second particle were originally in an aggregated state in the previous time step, it is determined that a split has occurred at the current time.

[0331] Based on the effective transport velocity, the spatial coordinates of all particles are updated in time step using a numerical integration algorithm. During the update process, coordinate merging and mass accumulation operations are performed on particle pairs that are determined to have aggregated, and coordinate separation and mass splitting operations are performed on particle clusters that are determined to have fragmented.

[0332] The simulation is executed in a loop with time increments until the preset total simulation duration is reached, and the final particle transport trajectory data is output to represent the ocean phenomenon process.

[0333] In this embodiment, within each time step, candidate particle pairs are sorted according to the difference between the work done by the coherent field and the separation energy, and the final interacting particle pair is determined using a non-overlapping matching rule; the same particle is only allowed to participate in the aggregation or fragmentation determination once within the same time step.

[0334] Based on the effective transport velocity calculated in step S3, the spatial coordinates of all particles are updated in time step using a numerical integration algorithm. The update process specifically includes:

[0335] Obtain the system's preset simulation time step Extract the effective transport velocity calculated in step S3, and construct a continuous velocity interpolation function with respect to spatial coordinates and time. The velocity interpolation function is used to obtain arbitrary spatial coordinates within the computational domain. and at any time The effective transport velocity vector.

[0336] The system's preset simulation time step The Courant-Friedrich-Lyuvis (CFL) stability condition for numerical fluid dynamics calculations must be strictly satisfied to prevent numerical divergence. The calculation logic is as follows:

[0337] ;

[0338] in, For spatial step size, To achieve the global maximum effective flow rate, This is the upper limit of the Courant number (for the fourth-order Runge-Kutta algorithm, it is usually taken as 0.5~0.8).

[0339] For each particle within the computational domain, including independent particles and aggregate particles generated by merging in step S43, obtain its position at the current time. spatial coordinates ;

[0340] The fourth-order Runge-Kutta numerical integration algorithm was used to calculate the effective transport velocity increments at four different spatial locations and time points:

[0341] Calculate the first velocity increment between the current time and the current position:

[0342] ;

[0343] Based on the first velocity increment, calculate the second velocity increment after half a time step:

[0344] ;

[0345] Based on the second velocity increment, calculate the third velocity increment after half a time step:

[0346] ;

[0347] Based on the third velocity increment, calculate the fourth velocity increment after one complete time step.

[0348] ;

[0349] in, , , and These are the first, second, third, and fourth velocity increment vectors, respectively. For time step, For the current moment, These are the current spatial coordinates.

[0350] Based on the weighted average principle, the four velocity increments are linearly combined to calculate the particle's velocity at the next moment. The spatial coordinates of are expressed mathematically as follows:

[0351] ;

[0352] in, The updated spatial coordinates for the next time step;

[0353] Traverse all particles within the computational domain to complete the synchronous update of the spatial coordinates of the global particle system.

[0354] The preset total simulation duration is determined by the physical life cycle of the specific marine phenomenon. It is automatically matched with big data based on the input phenomenon type (such as the red tide extinction cycle, oil spill weathering half-life, and golden time for search and rescue) or manually set by the user, serving as the global termination condition for time step updates.

[0355] It should be noted that if only the simple first-order Euler method is used, in areas with extremely large velocity gradients, such as the edge of an ocean eddy or a front, the particle trajectory will produce a serious truncation error, causing the particles to be thrown out of the eddy in a non-physical way.

[0356] In this application, the coherent structure itself is a highly nonlinear topological feature in the flow field. The RK4 algorithm, by performing four spatial and temporal trial interpolations within a single time step, can accurately capture the bending and shearing changes of the flow field within a small time scale, ensuring that particles strictly follow the effective transport velocity calculated in step S3 when updating their positions. This guarantees the computational accuracy of all the aforementioned physically consistent theories (such as energy balance and memory Stokes drift) at the engineering implementation level.

[0357] During the update process, particle state reconstruction based on the conservation of mass and momentum is performed:

[0358] For the first and second particles determined to have aggregated, coordinate merging and mass accumulation operations are performed to generate a new aggregate particle, the mathematical expression of which is:

[0359] ;

[0360] ;

[0361] in, The mass of the aggregate particles. The coordinates of the center of mass of the aggregate particles are given.

[0362] For aggregate particles that are determined to have split, coordinate separation and mass splitting operations are performed to divide their mass into two equal parts, and a small random spatial perturbation vector is applied to the centroid coordinates to generate two independent split particles.

[0363] Repeat the above steps until the preset total simulation time is reached, and output the final particle transport trajectory data to present the ocean phenomenon process.

[0364] It should be noted that this step uses coherent work and energy balance criteria to replace the thermodynamic free energy threshold. This application returns to the essence of classical mechanics, proposing a purely mechanical energy-based physical sufficiency criterion: quantifying the tendency of particle separation as the sum of relative kinetic energy and local dissipation, and converting the driving force of particle aggregation into work done by the coherent field. Aggregation will only occur spontaneously when the work done by the coherent field is sufficient to overcome the relative escape kinetic energy of the particles and the dissipation caused by environmental turbulence. This method completely eliminates artificial parameters such as effective temperature, allowing the aggregation and disintegration processes to be entirely determined by the energy distribution of the flow field itself, significantly improving the accuracy of predicting the evolution of ocean phenomena under extreme sea conditions.

[0365] It should be noted that in existing marine particle simulations (such as oil spill or Sargassum aggregation simulations), determining whether particles aggregate typically relies on two extremely coarse conditions: first, the geometric distance is less than a certain fixed threshold (e.g., 50 meters); second, an artificially defined thermodynamically effective temperature is used to compare the difference in free energy. This method not only lacks a hydrodynamic basis but also leads to non-physical rigid aggregation (oil films that cannot be torn apart by strong currents like steel plates) or spurious fragmentation (spreading out inexplicably in calm waters).

[0366] This application constructs a physically sufficient criterion from a purely mechanical energy perspective and achieves parameter self-consistency throughout the entire process. The resistance to aggregation is rigorously quantified into two parts: one is the particle's own escape tendency (relative kinetic energy). Secondly, there is the tearing effect of environmental turbulence. Ingeniously, this step directly utilizes the equivalent stiffness field learned through data-driven learning in step S1. To calculate local energy dissipation This allows the calculation of the energy barrier to be completely independent of manual settings, achieving deep self-consistency between the preceding and following steps. The power of accumulation is rigorously quantified as the work done by the coherent field. Aggregation only occurs spontaneously when the driving force strictly overcomes the resistance. This determining mechanism means that the tearing and merging of large-scale floating debris clumps is entirely determined by the energy distribution of the flow field itself. At strong shear fronts, due to the enormous work done, particles rapidly aggregate into bands; while in the region of high-dissipation breaking waves, due to... The aggregates are extremely large, and naturally disintegrate. This method completely eliminates artificial parameters such as effective temperature, significantly improving the accuracy and physical realism of predictions of the evolution of marine phenomena under extreme sea conditions.

[0367] Referring to Table 1, for different marine environmental conditions defined in the embodiments (submesoscale fronts, strong wind and wave breaking zones, transient eddy formation and extinction, abrupt changes in nearshore shallow waters, etc.), this system introduces virtual substances with different mass properties (such as oil spills, microplastics, and large seaweed rafts) for tracking and simulation verification. By applying the pattern recognition and energy balance evolution determination logic described in this invention, the key characteristic parameters and the final phenomenon reconstruction results are shown in Table 1. The values ​​of coherent field work and separation threshold in the table are normalized energy reference values.

[0368]

[0369] Table 1

[0370] As can be seen from the experimental data in Table 1, the multi-parameter joint judgment logic established in this invention forms a rigorous physical classification boundary. In particular, for the two working conditions that are prone to errors in traditional ocean numerical simulations—sub-mesoscale frontal over-diffusion and non-physical rigid accumulation of oil films under strong winds and waves—this invention accurately identifies the dominant flow field mode by using the real part of the maximum eigenvalue (the significant difference between <0 and ≥0), and strictly follows the energy balance law (whether the work done by the coherent field overcomes the relative kinetic energy and dissipation barrier), achieving a physically self-consistent switching between accumulation and disintegration states, thus possessing extremely high practical marine engineering forecasting value.

[0371] A particle-based ocean phenomenon process simulation system, used to execute a particle-based ocean phenomenon process simulation method, includes a basic field reconstruction module, a pattern recognition module, a dynamic velocity correction module, and an energy balance evolution module.

[0372] The basic field reconstruction module is used to acquire multi-source marine environmental data. Based on sparse regression and energy consistency functional minimization, it solves for the optimal weights to generate an equivalent stiffness field.

[0373] The pattern recognition module is used to calculate the local coherence time based on the logarithmic growth rate of the two-dimensional optimal transport distance of the particle cluster, and to identify the coherent dominant mode and the inertial dominant mode based on local linear stability analysis.

[0374] The dynamic velocity correction module is used to construct the particle velocity correction equation based on the local dissipation response approximation. The particle velocity correction equation includes the particle's memory Stokes drift velocity and effective transport velocity.

[0375] The energy balance evolution module is used to construct energy balance judgment rules based on the energy balance relationship between the work done by the coherent field and the relative kinetic energy to determine the aggregation state of particles, update the particle position according to the effective transport velocity and aggregation state, and output the final particle transport trajectory to present the process of ocean phenomena.

[0376] Those skilled in the art will understand that embodiments of the present invention can be provided as methods, systems, or computer program products. Therefore, the present invention can take the form of a completely hardware embodiment, a completely software embodiment, or an embodiment combining software and hardware aspects. Furthermore, the present invention can take the form of a computer program product implemented on one or more computer-usable storage media containing computer-usable program code. The storage medium can be implemented by any type of volatile or non-volatile storage device or a combination thereof, such as Static Random Access Memory (SRAM), Electrically Erasable Programmable Read-Only Memory (EEPROM), Erasable Programmable Read-Only Memory (EPROM), Programmable Red-Only Memory (PROM), Read-Only Memory (ROM), magnetic storage, flash memory, magnetic disk, or optical disk. These computer program instructions may also be stored in a computer-readable storage medium that can direct a computer or other programmable data processing device to function in a particular manner, such that the instructions stored in the computer-readable storage medium produce an article of manufacture including instruction means, which are implemented in a process Figure 1 One or more processes and / or boxes Figure 1 The function specified in one or more boxes.

[0377] It should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and are not intended to limit it. Although the present invention has been described in detail with reference to preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can be made to the technical solutions of the present invention without departing from the spirit and scope of the technical solutions of the present invention, and all such modifications or substitutions should be covered within the protection scope of the present invention.

Claims

1. A method for simulating ocean phenomena based on particle systems, characterized in that, Includes the following steps: Step S1: Acquire multi-source marine environmental data, and solve for the optimal weights based on sparse regression and energy consistency functional minimization to generate an equivalent stiffness field; Step S2: Based on the logarithmic growth rate of the optimal transport distance of the particle cluster in two dimensions, calculate the local coherence time, and identify the coherence-dominant mode and the inertia-dominant mode based on local linear stability analysis. Step S3: Construct the particle velocity correction equation based on the local dissipation response approximation. The particle velocity correction equation includes the particle's memory Stokes drift velocity and effective transport velocity. Step S4: Based on the energy balance relationship between the work done by the coherent field and the relative kinetic energy, construct an energy balance judgment rule to determine the aggregation state of the particles, update the particle position according to the effective transport velocity and aggregation state, and output the final particle transport trajectory to present the ocean phenomenon process.

2. The method for simulating ocean phenomena based on particle systems as described in claim 1, characterized in that: Step S1 includes the following sub-steps: Step S11: Obtain multi-source marine environmental data of the target sea area. The multi-source marine environmental data includes surface flow velocity, seawater density, significant wave height, wave period data and coastline geometry data. Based on the surface flow velocity, calculate the mass derivative of the fluid micro-particles. Step S12: Release virtual particles in the computational domain, use the lead particle tracking algorithm to perform trajectory integration on the virtual particles within a preset short time window, count the number of virtual particles on each grid node in the computational domain, and calculate the preliminary particle concentration field. Step S13: Extract the original physical features that affect ocean particle transport. The original physical features include wave gradient features, density gradient features and Stokes convergence features. Introduce a unified feature force density scale parameter to perform dimensionless processing on the original physical features and perform standardization processing on columns to obtain a standardized feature matrix. Step S14: Define the equivalent stiffness field as a linear combination of the features in each column of the standardized feature matrix, construct an energy-consistent functional, and solve for the optimal weights and generate the equivalent stiffness field based on sparse regression.

3. The method for simulating ocean phenomena based on particle systems as described in claim 2, characterized in that: Step S12 includes the following sub-steps: Step S121: Obtain the two-dimensional spatial range of the target sea area, discretize the two-dimensional spatial range into an Eulerian grid, and obtain the set of spatial coordinates of all grid nodes; Virtual particles are uniformly released within the computational domain at a preset initial spatial density. Each virtual particle is assigned a unique index number, and the initial spatial coordinates and initial mass of each virtual particle are recorded to obtain a set of virtual particles. Step S122: Set a preset short time window for pilot tracking. Based on the surface flow field velocity, use the fourth-order Runge-Kutta numerical integration algorithm to integrate the motion trajectory of all virtual particles within the preset short time window and obtain the end spatial coordinates of each virtual particle at the end of the integration. Step S123: Traverse all grid nodes of the Euler grid, and use the smoothing kernel function to map the virtual particle mass at the end of the integration to the adjacent grid nodes to calculate the preliminary particle concentration field on each grid node.

4. The method for simulating ocean phenomena based on particle systems as described in claim 3, characterized in that: Step S14 includes: The mathematical expression for the equivalent stiffness field is: ； in, For the equivalent stiffness field, Let be the non-negative weight vector to be determined. , and Non-negative weights , and These are the dimensionless features in the standardized feature matrix; Construct an energy consistency functional, the mathematical expression of which is: ； in, It is an energy-consistent functional. For the computational domain, The momentum residual term represents the squared error between the particle concentration diffusion flux divergence driven by the equivalent stiffness field and the actual inertial force of the fluid element. For L1 regularization terms, The regularization coefficient is . Let L1 be the norm of the weight vector; Solving for optimal weights and generating an equivalent stiffness field based on sparse regression specifically includes: The minimization problem of the energy consistency functional is transformed into a convex optimization problem with non-negativity constraints, and its mathematical expression is: ； ； in, The optimal non-negative weight vector. Function representation of taking and Minimum The value, , and These are the optimal non-negative weights that need to be substituted into the equivalent stiffness field formula; The above equations are solved using a convex optimization algorithm to obtain the optimal weight vector. The optimal weight vector is then substituted into the equivalent stiffness field formula to generate the final equivalent stiffness field.

5. The method for simulating ocean phenomena based on particle systems as described in claim 4, characterized in that: Step S2 includes the following sub-steps: Step S21: Initial position of the computational domain A first particle cluster and a second particle cluster are initialized at a certain point. The first particle cluster and the second particle cluster are spatially adjacent. The first particle cluster and the second particle cluster are subjected to time evolution based on the surface flow field velocity to obtain the probability measure distribution of the first particle cluster and the second particle cluster at time t. Step S22: Calculate the two-dimensional optimal transport distance between the first particle cluster and the second particle cluster during the evolution process, and calculate the logarithmic growth rate of the two-dimensional optimal transport distance within a local time window; Step S23: Calculate the inverse scale of the logarithmic growth rate and define it as the local coherence time. Construct the coherent tensor evolution equation. Linearize the coherent tensor evolution equation in the local equilibrium neighborhood to obtain a local linear operator. Perform eigenvalue decomposition on the local linear operator and extract the largest eigenvalue with the largest real part among all eigenvalues. Step S24: When the real part of the largest eigenvalue is less than zero, it indicates that the disturbance decays over time and the structure is stable, and the current region is determined to be a coherent dominant mode; When the real part of the largest eigenvalue is greater than or equal to zero, it indicates that the disturbance increases over time and the structure is easily broken, and the current region is determined to be in an inertial-dominated mode.

6. The method for simulating ocean phenomena based on particle systems as described in claim 5, characterized in that: Step S23 specifically includes: The mathematical expression for the local coherence time is: ； in, Initial position Local coherence time at that point This is the preset minimum growth rate lower limit. The growth rate is logarithmic. A coherent tensor evolution equation is constructed, and then linearized in the local equilibrium neighborhood to obtain a linear perturbation equation: ； in, For partial derivative operators, This represents the perturbation of the coherence tensor. To be in the initial position The local linear operator obtained after freezing the background field coefficients; For the local linear operator Perform eigenvalue decomposition, extract all eigenvalues, and obtain the largest eigenvalue with the largest real part. .

7. The method for simulating ocean phenomena based on particle systems as described in claim 6, characterized in that: Step S3 includes the following sub-steps: Step S31: Obtain the physical property parameters of the target particle and calculate the response time of the particle, wherein the response time is the ratio of the particle mass to the effective damping coefficient; Step S32: Calculate the memory Stokes drift velocity of the particle based on the equivalent stiffness field and local coherence time, using the local dissipative response approximation theory. Step S33: Construct a near-shore activation function based on spatial distance. Based on the pattern recognition results output in step S24, calculate the coherent induced velocity. Vector-add the surface flow velocity, the memory Stokes drift velocity, and the coherent induced velocity modulated by the near-shore activation function to obtain the effective transport velocity of the particles.

8. The method for simulating ocean phenomena based on particle systems as described in claim 7, characterized in that: Step S4 includes the following sub-steps: Step S41: Traverse any two adjacent first and second particles within the computational domain, and calculate the reduced mass and relative kinetic energy of the first and second particles; Step S42: Perform line integration along the spatial line connecting the first particle and the second particle to calculate the work done by the coherent field on the first particle and the second particle, and obtain the work done by the coherent field. Step S43: Obtain the local dissipated energy of the current sea area grid, add the relative kinetic energy and the local dissipated energy to obtain the separation energy threshold, determine the aggregation state of the first and second particles according to the energy balance judgment rule, update the particle position according to the effective transport velocity and aggregation state, and output the final particle transport trajectory to present the ocean phenomenon process.

9. The method for simulating ocean phenomena based on particle systems as described in claim 8, characterized in that: Step S43 specifically includes: The local dissipated energy is calculated using the equivalent stiffness field and particle spacing, and its mathematical expression is as follows: ； in, This represents the energy dissipation barrier caused by environmental turbulent diffusion on particle aggregation. This represents the equivalent stiffness field at the midpoint of the spatial connection. Spatial coordinates of the first particle Spatial coordinates of the second particle The Euclidean distance; The relative kinetic energy is added to the local dissipated energy to obtain the separation energy threshold; The specific rules for determining energy balance include: When the work done by the coherent field is strictly greater than the separation energy threshold, it is determined that the first particle and the second particle have aggregated. When the work done by the coherent field is not greater than the separation energy threshold, it is determined that the first particle and the second particle remain in a dispersed state. If the first particle and the second particle were originally in an aggregated state in the previous time step, it is determined that a split has occurred at the current time. Based on the effective transport velocity, the spatial coordinates of all particles are updated in time step using a numerical integration algorithm. During the update process, coordinate merging and mass accumulation operations are performed on particle pairs that are determined to have aggregated, and coordinate separation and mass splitting operations are performed on particle clusters that are determined to have fragmented. The simulation is executed in a loop with time increments until the preset total simulation duration is reached, and the final particle transport trajectory data is output to represent the ocean phenomenon process.

10. A particle system-based ocean phenomenon process simulation system, the system being used to execute a particle system-based ocean phenomenon process simulation method, characterized in that, It includes a fundamental field reconstruction module, a pattern recognition module, a dynamic velocity correction module, and an energy balance evolution module: The basic field reconstruction module is used to acquire multi-source marine environmental data and solve for the optimal weights to generate an equivalent stiffness field based on sparse regression and energy consistency functional minimization. The pattern recognition module is used to calculate the local coherence time based on the logarithmic growth rate of the two-dimensional optimal transport distance of the particle cluster, and to identify the coherent dominant mode and the inertial dominant mode based on local linear stability analysis. The dynamic velocity correction module is used to construct a particle velocity correction equation based on the local dissipation response approximation. The particle velocity correction equation includes the particle's memory Stokes drift velocity and effective transport velocity. The energy balance evolution module is used to construct energy balance judgment rules based on the energy balance relationship between the work done by the coherent field and the relative kinetic energy to determine the aggregation state of particles, update the particle position according to the effective transport velocity and aggregation state, and output the final particle transport trajectory to present the process of ocean phenomena.

Images