A numerical simulation-based multi-element coupling sandy coast ecological early warning method

Abstract

The present application provides a kind of sandy coast ecological early warning method based on numerical simulation multi-element coupling, belong to sandy coast ecological early warning technical field, the present application obtains boundary condition by arranging wave observation buoy array and sediment sampling point in sandy coast monitoring area, constructs three-dimensional hydrodynamic-sediment transport coupling calculation grid and carries out grid optimization using dynamic partition algorithm based on hilbert curve mapping, using multi-time step algorithm combined with operator splitting technique to establish conservation type coupling interface respectively solve wave propagation equation and geomorphology evolution equation, the physical field data obtained by calculation is input into multi-field coupling prediction model to output prediction value, according to the prediction value, start fluid-solid coupling calculation and calculate ecological risk index to divide early warning level, solve the technical problem that the accuracy and efficiency of multi-physical field coupling calculation in sandy coast ecological early warning are difficult to balance.

Description

Technical Field

[0001] This invention belongs to the field of ecological early warning technology for sandy coasts, and more specifically, it relates to an ecological early warning method for sandy coasts based on multi-element coupling of numerical simulation. Background Technology

[0002] Early warning systems for sandy coasts require comprehensive consideration of the coupled effects of wave dynamics, sediment transport, and geomorphological evolution. Traditional methods employ direct coupling calculations with a uniform time step to achieve multiphysics numerical simulations. However, since wave propagation occurs on a timescale of seconds while geomorphological evolution occurs on a timescale of hours, a uniform time step results in extremely low computational efficiency, while a coarse time step fails to accurately capture the rapid changes in wave propagation. In current sandy coast monitoring and early warning practices, the timescales of the wave field and the geomorphological field differ by more than three orders of magnitude, making it difficult for existing technologies to achieve long-term geomorphological evolution prediction while maintaining accuracy in wave propagation calculations. In other words, existing technologies face the challenge of balancing accuracy and efficiency in multiphysics coupling calculations for sandy coast ecological early warning. Summary of the Invention

[0003] In view of this, the present invention provides a method for ecological early warning of sandy coasts based on multi-element coupling of numerical simulation, which can solve the technical problem of difficulty in balancing the accuracy and efficiency of multi-physics field coupling calculation in ecological early warning of sandy coasts in the prior art.

[0004] This invention is implemented as follows: It provides a multi-element coupled early warning method for sandy coastlines based on numerical simulation. Wave observation buoy arrays and sediment sampling points are deployed in the sandy coastline monitoring area to collect wave parameter and sediment grain size distribution data. A two-dimensional frequency-direction spectrum is obtained as the far-field boundary condition through a wave energy spectrum direction decomposition and reconstruction algorithm. A three-dimensional hydrodynamic-sediment transport coupled computational grid is constructed, and a dynamic partitioning algorithm based on Hilbert curve mapping is used to locally refine the unstructured grid in complex terrain areas. A multi-time-step algorithm is used to solve the wave propagation equation and the geomorphic evolution equation. An explicit calculation format for the wave process time step is set for the wave propagation equation, and an implicit calculation format for the geomorphic evolution time step is set for the geomorphic evolution equation. A conserved coupling interface is established through operator splitting technology, and turbulence model parameters are input as random parameters into a polynomial chaotic expansion model. The model employs Legendre polynomial series expansion and Galerkin projection to transform it into a deterministic extended equation system for solution. An algebraic multigrid preconditioner combined with a generalized minimum residual iteration method is used to solve the sparse linear equation system formed after discretization of the three-dimensional hydrodynamic control equations. The calculated velocity field, wave height field, and suspended sediment concentration field data are input into a multi-field coupled prediction model, which outputs predicted values ​​for coastal erosion rate, shoreline retreat distance, and damaged area in ecologically sensitive zones. When the predicted values ​​exceed a threshold, fluid-structure interaction calculations of the sediment-ecology module are initiated. A hypergrid-based conservation interpolation algorithm is used to transfer bed shear stress and sediment flux, calculating the bed shear stress distribution within the ecologically sensitive zone. When the bed shear stress exceeds the critical initiation shear stress, it is determined to be in an erosion state. The proportion of erosion-state grid cells to the total area of ​​the ecologically sensitive zone is statistically analyzed as an ecological risk index, and warning levels are classified based on the ecological risk index.

[0005] The specific steps of the wave energy spectrum direction decomposition and reconstruction algorithm are as follows: using the finite-direction wave data collected by the wave observation buoy array as the observation constraint, using the maximum entropy method to invert the two-dimensional frequency-direction spectrum function, constructing the Lagrangian function to transform the constrained optimization problem into an unconstrained problem, and solving the nonlinear optimization equation through the conjugate gradient iteration method.

[0006] The maximum entropy method expresses the entropy function as follows: after taking the negative logarithm of the integral of the two-dimensional frequency-direction spectrum function in the frequency-direction space, the integral is applied over the entire domain. Under the conditions of satisfying observation constraints and non-negativity constraints, the entropy function reaches its maximum value. After introducing observation constraints through Lagrange multipliers, the optimization equation is obtained by taking the variational equation of the two-dimensional frequency-direction spectrum function.

[0007] Among them, the Tikhonov regularization term is introduced into the objective function to suppress numerical oscillations in high-frequency directions. The Tikhonov regularization term is obtained by adding an integral term of the second derivative of the two-dimensional frequency-direction spectrum function with respect to the direction angle to the objective function. The smoothness is controlled by multiplying the integral term by a regularization parameter. The regularization parameter is adaptively determined by the L-curve method based on the noise level of the finite direction wave data.

[0008] Specifically, the dynamic partitioning algorithm based on Hilbert curve mapping involves mapping the geometric center coordinates of the two-dimensional unstructured grid cells formed after local refinement of the unstructured grid to the one-dimensional sequence index of the Hilbert curve, and dividing the one-dimensional sequence index into continuous segments according to the number of parallel computing processes and allocating them to each process.

[0009] Specifically, after mesh refinement or coarsening, the one-dimensional sequence index of the Hilbert curve is recalculated and load balancing detection is triggered. When the ratio of the standard deviation of the number of mesh cells in each process to the average value exceeds 0.15, incremental repartitioning is initiated. A multi-level graph segmentation algorithm is used to identify the minimum cut set on the boundary of each partition and migrate the boundary cells to the adjacent process.

[0010] The multi-level graph segmentation algorithm includes a coarsening stage, an initial segmentation stage, and a refinement stage. In the coarsening stage, a layer-by-layer unstructured grid cell is merged using a multi-edge matching algorithm to construct a hierarchical structure. In the initial segmentation stage, a greedy algorithm or a spectral segmentation algorithm is used at the coarsest layer to obtain the initial partitions. In the refinement stage, the Kernighan-Lin algorithm is used to optimize the partition boundaries during the layer-by-layer refinement process.

[0011] The steps of establishing a conserved coupling interface using operator splitting technology specifically involve decomposing the wave propagation equation and the geomorphic evolution equation into convection operators, diffusion operators, and source term operators according to physical processes. Within the wave process time step of solving the wave propagation equation, only the convection operators and diffusion operators are solved. Within the geomorphic evolution time step of solving the geomorphic evolution equation, the wave radiation stress of multiple wave process time steps is accumulated as the input of the source term operator.

[0012] Among them, a two-way feedback mechanism between the wave field and the geomorphic field is established at the conservation coupling interface through flux conservation conditions. The flux conservation conditions require that the integral of wave radiation stress on the seabed surface is equal to the momentum change caused by sediment transport. The error introduced by the operator splitting technique is compensated by introducing a modified flux term at the conservation coupling interface. The modified flux term is explicitly estimated based on the solution of the previous time step.

[0013] The steps of performing Legendre polynomial series expansion on random parameters in the polynomial chaotic expansion model are as follows: the mixed length coefficient, bottom roughness coefficient, and suspended sediment settling velocity coefficient are expressed as probability distribution functions, the orthogonal polynomial family corresponding to the probability distribution function is selected as the basis function, and the mixed length coefficient, bottom roughness coefficient, and suspended sediment settling velocity coefficient are expanded into a finite series sum of orthogonal polynomial families.

[0014] This invention establishes a conservation-type coupling interface by employing a multi-time-step algorithm combined with operator splitting technology. It sets an explicit calculation format of 0.01 to 0.1 seconds for the wave propagation equation to capture rapidly changing wave processes, and an implicit calculation format of 1 to 24 hours for the geomorphic evolution equation to describe slowly evolving terrain changes. At the conservation-type coupling interface, a bidirectional feedback mechanism between the wave field and the geomorphic field is established through flux conservation conditions, overcoming the low computational efficiency of traditional unified time-step methods. Through operator splitting technology, the wave propagation equation and the geomorphic evolution equation are decomposed into convection operators, diffusion operators, and source term operators according to physical processes. These are solved separately at different time scales, and the error introduced by operator splitting is compensated by correcting the flux term. This ensures both the time accuracy of wave propagation calculations and the computational efficiency for long-term geomorphic evolution prediction. In summary, this invention solves the technical problem mentioned in the background art of balancing accuracy and efficiency in multi-physics coupling calculations for ecological early warning of sandy coastlines. Attached Figure Description

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

[0016] Figure 2 The flowchart shows the time progression process for the multi-time-step algorithm.

[0017] Figure 3 This is a map showing the topographic conditions of the water depth.

[0018] Figure 4 This is a simulation result of the nearshore wave field of a sandy coast using Spectral Waves FM.

[0019] Figure 5 This is a simulation result of storm surge on nearshore sandy coasts using Flow MODEL FM.

[0020] Figure 6 This is a simulation result of nearshore currents on sandy coasts using Flow MODEL FM.

[0021] Figure 7 This is a diagram showing the simulation results of nearshore biological oxygen demand in sandy coastal areas from ECO Lab.

[0022] Figure 8 This is a figure showing the results of the ECO Lab numerical simulation of nearshore chemical oxygen demand for sandy coastlines. Detailed Implementation

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

[0024] like Figure 1 The diagram shown is a flowchart of a multi-factor coupled ecological early warning method for sandy coastlines based on numerical simulation provided by this invention. This method includes the following steps: S01. Deploy wave observation buoy arrays and sediment sampling points in the sandy coastal monitoring area to collect data on wave effective wave height, wave peak period, wave main wave direction angle and sediment grain size distribution. Obtain the complete two-dimensional frequency-direction spectrum as far-field boundary conditions through the wave energy spectrum direction decomposition and reconstruction algorithm. S02. Construct a three-dimensional hydrodynamic-sediment transport coupled computational grid, and use a dynamic partitioning algorithm based on Hilbert curve mapping to locally refine the unstructured grid in the complex topographic area of ​​the sandy coast monitoring area. Set a semi-coarsening strategy in the water depth direction to maintain the vertical grid accuracy, and establish a real-time load monitoring mechanism. S03. A multi-time-step algorithm is used to solve the wave propagation equation and the landform evolution equation. An explicit calculation format with a wave process time step of 0.01 seconds to 0.1 seconds is set for the wave propagation equation, and an implicit calculation format with a landform evolution time step of 1 hour to 24 hours is set for the landform evolution equation. A conserved coupling interface is established through operator splitting technology. S04. Input the mixing length coefficient, bottom roughness coefficient and suspended sediment settling velocity coefficient in the turbulence model as random parameters into the polynomial chaotic expansion model, perform Legendre polynomial series expansion on the random parameters, and use the Galerkin projection method to transform the stochastic partial differential equations into a deterministic extended equation system and solve it. S05. The sparse linear equation system formed by discretizing the three-dimensional hydrodynamic control equations is solved by combining algebraic multigrid preconditioners with the generalized minimum residual iterative method. Multi-level iterative parameters of coarsening-interpolation-smoothing are set, and the sparse matrix-vector multiplication operation is accelerated by using a graphics processor. S06. Input the calculated velocity field, wave height field and suspended sediment concentration field data into the multi-field coupled prediction model. The multi-field coupled prediction model outputs the predicted values ​​of the coastal erosion rate, the predicted value of the shoreline retreat distance and the predicted value of the damaged area of ​​the ecologically sensitive area for the next 7 to 30 days. S07. When the predicted value of the coastal erosion rate exceeds 150% of the annual average erosion rate threshold or the predicted value of the shoreline retreat distance exceeds 5 meters, the fluid-structure interaction calculation of the sediment-ecology module is initiated, and a conservation interpolation algorithm based on supergrid is used to transfer the bed shear stress and sediment flux in the overlapping area of ​​the fluid grid and sediment grid. S08. Calculate the distribution of bed shear stress and sediment resuspension flux in the ecologically sensitive area. When the bed shear stress in the distribution of bed shear stress exceeds the critical starting shear stress, it is determined to be an erosion state. The proportion of the erosion state grid cells to the total area of ​​the ecologically sensitive area is used as the ecological risk index. S09. Based on the ecological risk index, the early warning level is divided. When the ecological risk index is less than 10%, it is determined to be in a normal state. When the ecological risk index is between 10% and 30%, it is determined to be in a state of concern and the monitoring frequency is increased. When the ecological risk index is between 30% and 60%, it is determined to be in an early warning state and the protection measures optimization calculation is initiated. When the ecological risk index exceeds 60%, it is determined to be in an alarm state and an emergency response plan is output.

[0025] The specific steps of the wave energy spectrum direction decomposition and reconstruction algorithm include: using the finite-direction wave data collected by the wave observation buoy array as observation constraints, using the maximum entropy method to invert the two-dimensional frequency-direction spectrum function, constructing a Lagrangian function to transform the constrained optimization problem into an unconstrained problem, solving the nonlinear optimization equation through the conjugate gradient iteration method, introducing a Tikhonov regularization term in the objective function to suppress numerical oscillations in high-frequency directions, and using a mixture of Gaussian distribution models to fit the peak direction and direction broadening parameters of each wave system for wave systems with multiple wave peaks.

[0026] The entropy function of the maximum entropy method is expressed as follows: after taking the negative logarithm of the integral of the two-dimensional frequency-direction spectrum function in the frequency-direction space, the integral is applied over the entire domain. Under the conditions of satisfying the observation constraints and non-negativity constraints, the entropy function reaches its maximum value. After introducing the observation constraints through Lagrange multipliers, the variational equation of the two-dimensional frequency-direction spectrum function is obtained to obtain the optimization equation.

[0027] The Tikhonov regularization term is achieved by adding an integral term of the second derivative of the two-dimensional frequency-direction spectrum function with respect to the direction angle to the objective function. The integral term is multiplied by a regularization parameter to control the smoothness. The regularization parameter is adaptively determined using the L-curve method based on the noise level of the finite-direction wave data.

[0028] The specific steps of the dynamic partitioning algorithm based on Hilbert curve mapping include: mapping the geometric center coordinates of the two-dimensional unstructured grid cells formed after local densification of the unstructured grid to the one-dimensional sequence index of the Hilbert curve; dividing the one-dimensional sequence index into continuous segments according to the number of parallel computing processes and allocating them to each process; recalculating the one-dimensional sequence index of the Hilbert curve after grid densification or coarsening and triggering load balancing detection; starting incremental repartitioning when the ratio of the standard deviation to the average value of the number of grid cells in each process exceeds 0.15; and using a multi-level graph segmentation algorithm to identify the minimum cut set on the boundary of each partition and migrate the boundary cells to adjacent processes.

[0029] The Hilbert curve mapping ensures that spatially adjacent two-dimensional unstructured grid cells also remain adjacent in the one-dimensional sequence index, reducing communication overhead in parallel computing.

[0030] The multi-level graph segmentation algorithm includes a coarsening stage, an initial segmentation stage, and a refinement stage. In the coarsening stage, the two-dimensional unstructured mesh cells are merged layer by layer to construct a hierarchical structure through a multiple edge matching algorithm. In the initial segmentation stage, a greedy algorithm or a spectral segmentation algorithm is used at the coarsest layer to obtain the initial partitions. In the refinement stage, the partition boundaries are optimized during the layer-by-layer refinement process using the Kernighan-Lin algorithm.

[0031] The steps of establishing a conservation coupling interface using the operator splitting technique specifically include: decomposing the wave propagation equation and the landform evolution equation into convection operators, diffusion operators, and source operators according to physical processes; solving only the convection operators and the diffusion operators within the wave process time step of solving the wave propagation equation; accumulating wave radiation stresses of multiple wave process time steps as input to the source operator within the landform evolution time step of solving the landform evolution equation; and establishing a bidirectional feedback mechanism between the wave field and the landform field at the conservation coupling interface through flux conservation conditions.

[0032] The flux conservation condition requires that the integral of the wave radiation stress on the seabed surface equals the momentum change caused by sediment transport. The error introduced by the operator splitting technique is compensated by introducing a modified flux term at the conservation coupling interface. The modified flux term is explicitly estimated based on the solution of the previous time step.

[0033] The bidirectional feedback mechanism includes topographic feedback and dynamic feedback. The topographic feedback changes the wave propagation path and wave height distribution by updating the water depth field, while the dynamic feedback changes the bed shear stress distribution and sediment transport direction by updating the wave field.

[0034] The specific steps of performing Legendre polynomial series expansion on the random parameters in the polynomial chaotic expansion model include: representing the mixing length coefficient, the bottom roughness coefficient, and the suspended sediment settling velocity coefficient as probability distribution functions; selecting an orthogonal polynomial family corresponding to the probability distribution function as a basis function; expanding the mixing length coefficient, the bottom roughness coefficient, and the suspended sediment settling velocity coefficient into a finite series sum of the orthogonal polynomial family; expanding the state variables in the stochastic partial differential equation into a series sum of the same orthogonal polynomial family; using the Galerkin projection method to perform inner product projection on the expanded stochastic partial differential equation in the orthogonal polynomial space to obtain and solve the deterministic extended equation system with respect to the expansion coefficients; and directly calculating the statistical moment information of the state variables based on the expansion coefficients.

[0035] The selection principle for the orthogonal polynomial family is as follows: when the mixing length coefficient, the bottom roughness coefficient, and the suspended sediment settling velocity coefficient follow a uniform distribution, the Legendre polynomial is selected; when the mixing length coefficient, the bottom roughness coefficient, and the suspended sediment settling velocity coefficient follow a normal distribution, the Hermit polynomial is selected; and when the mixing length coefficient, the bottom roughness coefficient, and the suspended sediment settling velocity coefficient follow a beta distribution, the Jacobi polynomial is selected.

[0036] The statistical moment information includes the expected value and the variance. The expected value is equal to the zero-order expansion coefficient, and the variance is equal to the sum of the squares of the higher-order expansion coefficients multiplied by the square of the norm of the corresponding orthogonal polynomial family.

[0037] The steps of the algebraic multigrid preconditioner specifically include: constructing a coarsening operator to form a coarse grid hierarchy by aggregating adjacent grid nodes; constructing an interpolation operator to define the mapping relationship between coarse grid variables and fine grid variables; constructing a constraint operator to define the projection relationship between fine grid residuals and coarse grids; performing smooth iteration in the finest grid layer to reduce high-frequency error components; solving the residual equation in the coarse grid layer to correct low-frequency error components; and passing the coarse grid correction value back to the fine grid through the interpolation operator for post-smoothing iteration.

[0038] The coarsening operator adopts the Ruge-Stuben coarsening strategy, which identifies coarse grid nodes by analyzing the strong connection relationships of the coefficient matrix. The strong connection determination criterion is: when the absolute value of the matrix element between nodes is greater than 0.25 times the maximum absolute value of the off-diagonal element in that row, it is determined to be a strong connection.

[0039] The smoothing iteration employs Gauss-Seidel iteration or incomplete LU decomposition iteration, and the number of iterations of the smoothing iteration is adaptively adjusted according to the residual descent rate.

[0040] The multi-field coupled prediction model has the following structure: the input layer receives the values ​​of six physical field variables at spatial grid nodes: the horizontal and vertical components of the velocity field, the effective wave height and spectral peak period of the wave height field, and the surface and bottom concentrations of the suspended sediment concentration field. Spatial features are extracted through a three-dimensional convolutional layer. The encoder part adopts a representation learning framework based on the information bottleneck principle, containing four residual convolutional blocks. Each residual convolutional block contains two convolutional layers and one skip connection. The kernel size is 3×3×3, and the number of channels is 64, 128, 256, and 512 respectively. Batch connections are made after each residual convolutional block. A normalization layer and a modified linear unit activation function are connected to a variational coding module at the output of the encoder section. The variational coding module includes a mean coding branch and a variance coding branch. The mean coding branch and the variance coding branch respectively map the feature vector to the mean vector and log-variance vector of the latent representation through a fully connected layer. The dimension of the latent representation is set to 128. The decoder section adopts a mirror-symmetric structure and gradually restores the spatial resolution through a transposed convolutional layer. The output layer generates the predicted values ​​of the coastal erosion rate, the predicted values ​​of the shoreline retreat distance, and the predicted values ​​of the damaged area of ​​the ecologically sensitive area through a 1×1×1 convolutional layer.

[0041] The representation learning framework based on the information bottleneck principle applies information compression constraints to the latent representation layer, forcing the encoder to extract the key features most relevant to the prediction target from the input data while filtering redundant information. The core of this framework is to minimize the mutual information between the input data and the latent representation while maximizing the mutual information between the latent representation and the prediction target. Since direct computation of the mutual information is infeasible in high-dimensional continuous space, a variational approximation technique is used to introduce a variational distribution to approximate the true posterior distribution, transforming the computation of the mutual information into an optimization problem of a variational lower bound. Reparameterization techniques are used to make the sampling process differentiable, thus supporting backpropagation training. A mutual information regularization term is introduced into the loss function to balance compression ratio and fidelity. The coefficients of the mutual information regularization term are dynamically adjusted as Lagrange multipliers. When the prediction error is large, the compression intensity is reduced to allow more information to be transmitted; when the prediction error is small, the compression intensity is increased to improve the model's generalization ability.

[0042] The information compression constraint enhances the robustness of the multi-field coupled prediction model against observation noise and parameter uncertainty because redundant information and noise components are effectively suppressed during encoding. The latent representation preserves stable structured patterns in the data rather than specific numerical details. For the multi-factor coupled problem of ecological early warning for sandy coastlines, complex nonlinear interactions exist between different physical fields. Some coupling relationships contribute significantly to the early warning results, while others have a weaker impact. The representation learning framework based on the information bottleneck principle can automatically identify dominant coupling patterns and encode them into compact representations, avoiding overfitting of the multi-field coupled prediction model to accidental correlations in the training data. This significantly improves the prediction accuracy for unseen conditions and the ability to capture long-term evolutionary trends.

[0043] Furthermore, the uncertainty quantification function provided by the representation learning framework based on the information bottleneck principle reflects the prediction confidence through the variance component of the latent representation. When the input data deviates from the training distribution, the increased variance automatically triggers an early warning reliability assessment, providing decision-makers with risk interval information in addition to point prediction values, thus enhancing the practicality and reliability of the early warning system. The representation learning framework based on the information bottleneck principle continuously monitors the information content of the latent representation during the optimization process through the mutual information regularization term, enabling the encoder part to compress the representation dimension as much as possible while ensuring prediction performance, forming a stable feature space insensitive to input disturbances. When wave propagation paths shift due to terrain changes or sediment particle size distribution fluctuates due to seasonal factors, the multi-field coupled prediction model can still maintain stable prediction output without drastic fluctuations. This stability stems from the representation learning framework based on the information bottleneck principle forcing the encoder part to learn the essential correlations of physical laws rather than the surface statistical characteristics of the data.

[0044] The steps for establishing the training dataset for the multi-field coupled prediction model specifically include: selecting 10 to 20 typical storm events from historical periods as training samples, with each typical storm event containing continuous observation data from 48 hours before the event to 72 hours after the event; performing numerical simulations on each training sample using the three-dimensional hydrodynamic-sediment transport coupled computational grid to obtain time-series data of the velocity field, wave height field, and suspended sediment concentration field at the spatial grid nodes; extracting physical field data from the time-series data for 24 hours before and after the storm peak as input features, and extracting measured coastal erosion rate, measured shoreline retreat distance, and measured damaged area of ​​ecologically sensitive areas 7 to 30 days after the storm as output labels; normalizing the input features and output labels to map the numerical range to the interval between 0 and 1, with normalization parameters including the minimum and maximum values ​​of each physical field variable; and dividing the training samples into a training set and a validation set in an 8:2 ratio, with the training set used for model parameter updates and the validation set used for hyperparameter selection and early termination strategy determination.

[0045] The training steps of the multi-field coupled prediction model specifically include: initializing the convolutional kernel weights of the encoder and decoder parts using the Kaiming He initialization method, initializing the bias term to 0, and randomly initializing the weights of the variational coding module using a normal distribution; setting the batch size to 8, the initial learning rate to 0.001, updating the model parameters using the Adam optimizer, and setting the momentum parameters to 0.9 and 0.999; constructing a composite loss function including a reconstruction loss term, a mutual information regularization term, and a prediction loss term, wherein the reconstruction loss term uses mean squared error to measure the difference between the decoded output and the input data, the mutual information regularization term uses KL divergence to measure the difference between the latent representation distribution and the standard normal distribution, and the prediction loss term uses mean squared error to measure the predicted coastal erosion rate. The differences between the predicted value of the shoreline retreat distance, the predicted value of the damaged area of ​​the ecologically sensitive area, and the output label are calculated. The weight coefficients of the reconstruction loss term, the mutual information regularization term, and the prediction loss term are set to 1.0, 0.5, and 2.0, respectively. The composite loss function value is calculated by forward propagation on the training set, and the gradient is calculated and the model parameters are updated by backpropagation algorithm. After each training epoch, the prediction error is evaluated on the validation set. When the prediction error on the validation set does not decrease for 10 consecutive epochs, an early stopping mechanism is triggered to terminate the training. The model parameters with the smallest prediction error on the validation set are saved as the final model. During the training process, the learning rate is adjusted every 5 epochs. The learning rate decays to 0.1 times the initial value of the learning rate according to the cosine annealing strategy.

[0046] The specific steps of the hypergrid-based conservation interpolation algorithm include: identifying overlapping regions between the fluid grid and the sediment grid in space; constructing a hypergrid unit for each pair of overlapping fluid grid units and sediment grid units, where the vertex of the hypergrid unit is the intersection of the boundaries of the two grid units; calculating the area of ​​each hypergrid unit and recording its affiliation with the source grid unit and the target grid unit; interpolating and transferring the bed shear stress distribution in the fluid grid; calculating the area weight of all hypergrid units within each sediment grid unit, where the area weight is equal to the area of ​​the hypergrid unit divided by the total area of ​​the sediment grid units; the bed shear stress value of the target grid unit is equal to the sum of the products of the bed shear stress values ​​of the source grid units corresponding to all hypergrid units and the area weight; interpolating and transferring the sediment flux in the sediment grid; calculating the sediment flux value received by the fluid grid unit using the same area weight method; and verifying conservation by checking whether the flux integrals of the source grid unit and the target grid unit are equal.

[0047] The hypergrid-based conservation interpolation algorithm employs a covariant component transformation method for conservation interpolation of vector fields. This method converts the components of the vector field in the source grid coordinate system into coordinate-independent geometric objects, which are then re-decomposed into coordinate components in the target grid coordinate system. This avoids physical inconsistencies caused by directly interpolating coordinate components.

[0048] The critical initiation shear stress is determined based on a Shields curve, which describes the relationship between the dimensionless critical shear stress and the particle Reynolds number. The dimensionless critical shear stress is equal to the critical initiation shear stress divided by the product of the underwater specific weight of the sediment particles and the median particle size. The particle Reynolds number is equal to the product of the frictional velocity and the median particle size divided by the kinematic viscosity of water. The frictional velocity is equal to the square root of the bottom shear stress divided by the water density. For sandy sediments with a median particle size of 0.1 mm to 0.5 mm, the dimensionless critical shear stress ranges from 0.03 to 0.06.

[0049] The steps for optimizing the protective measures specifically include: when the early warning state is triggered, adding an artificial sandbar or underwater submerged dike to the three-dimensional hydrodynamic-sediment transport coupled calculation grid, adjusting the position parameters, elevation parameters, and length parameters of the artificial sandbar or underwater submerged dike, rerunning the three-dimensional hydrodynamic-sediment transport coupled calculation grid to calculate the impact of the artificial sandbar or underwater submerged dike on wave propagation and sediment transport, evaluating the reduction of the ecological risk index under different protection schemes, and selecting the protection scheme with the largest reduction in the ecological risk index and an engineering investment cost lower than the budget constraint as the optimization result output.

[0050] The present invention also provides a sandy coast ecological early warning system implemented by a computer, wherein the computer is equipped with a storage medium, the storage medium stores program instructions, and the program instructions execute the above-mentioned sandy coast ecological early warning method based on numerical simulation and multi-element coupling when the computer is run.

[0051] The specific implementation methods of the above steps are described in detail below.

[0052] The specific implementation of step S01 involves deploying a wave observation buoy array along the shoreline and offshore directions in the sandy coastal monitoring area. The buoy spacing is determined based on the monitoring area and wavelength characteristics, with a typical spacing of 200 to 500 meters. The buoys are equipped with accelerometers and gyroscopes to collect real-time data on significant wave height, peak period, and main wave direction angle. The sampling frequency is set to 2 Hz to 10 Hz to capture the high-frequency components of wave motion. Simultaneously, representative cross-sections within the monitoring area are selected to deploy sediment sampling points at depths ranging from 0 to 0.5 meters above the surface. A laser grain size analyzer is used to determine the sediment grain size distribution data, obtaining the median grain size and cumulative grain size curve. The wave data collected by the buoy array in finite directions is used as observation constraints and input into a wave energy spectrum direction decomposition and reconstruction algorithm. This algorithm is based on the maximum entropy method to invert the two-dimensional frequency-direction spectrum function. The core of the maximum entropy method is to maximize the information entropy of the spectrum while satisfying the observation constraints, thereby obtaining a spectrum estimate with minimal assumptions. The entropy function is constructed by taking the negative logarithm of the integral of the two-dimensional frequency-direction spectrum function over the frequency-direction space and integrating it over the entire domain. Observation constraints and nonnegativity constraints are introduced into the objective function using Lagrange multipliers. Variational optimization of the two-dimensional frequency-direction spectrum function yields a nonlinear optimization equation. This equation is solved using the conjugate gradient iteration method until the objective function converges. The convergence criterion is that the relative change of the objective function in adjacent iteration steps is less than 1%. To suppress numerical oscillations in the high-frequency direction, a Tikhonov regularization term is introduced into the objective function. This regularization term is the integral of the two-dimensional frequency-direction spectrum function with respect to the second derivative of the direction angle multiplied by a regularization parameter. The regularization parameter is adaptively determined based on the signal-to-noise ratio of the observed data using the L-curve method, and its typical value range is [range missing]. to For complex wave systems with multiple wave crests, a Gaussian mixture model is used to fit the peak direction and directional broadening parameters of each wave system. The weighting coefficients, mean, and variance of the Gaussian components are iteratively optimized using the expectation-maximization algorithm, and finally a complete two-dimensional frequency-direction spectrum is obtained as the far-field boundary condition for numerical simulation.

[0053] The specific implementation of step S02 involves constructing a three-dimensional hydrodynamic-sediment transport coupled computational grid covering the sandy coastline monitoring area. Horizontally, an unstructured triangular grid is used to adapt to complex shoreline and topographic variations. Vertically, a Sigma coordinate system is used to map the irregular water depth in physical space into a regular hierarchical structure in computational space. For complex terrain areas with large topographic gradients and drastic water depth variations, a dynamic partitioning algorithm based on Hilbert curve mapping is used for local grid refinement. The refinement criterion is that the refinement operation is triggered when the water depth difference between adjacent grid cells exceeds 30% of the shallower cell's water depth. The Hilbert curve is a space-filling curve that can map two-dimensional grid cells to a one-dimensional sequence while maintaining spatial locality; that is, spatially adjacent grid cells remain adjacent in the one-dimensional sequence. This property significantly reduces the communication overhead between processes in parallel computing. Using the geometric center coordinates of the unstructured grid cells as input, a recursive subdivision algorithm is used to calculate the one-dimensional sequence index of each grid cell on the Hilbert curve. The one-dimensional sequence index is divided into continuous segments according to the number of parallel computing processes and assigned to each process. After mesh refinement or coarsening, the Hilbert curve mapping is recalculated, triggering load balancing detection. The load balancing evaluation metric is the ratio of the standard deviation to the average number of mesh cells in each process. When this ratio exceeds 0.15, a load imbalance is determined, and incremental repartitioning is initiated. Repartitioning employs a multi-level graph partitioning algorithm, treating the mesh as a graph structure with mesh cells as nodes and adjacency relationships as edges. A coarsened hierarchical structure is constructed by merging mesh cells layer by layer using a multi-edge matching algorithm. Initial partitions are obtained at the coarsest layer using a greedy algorithm or spectral partitioning algorithm. Then, the Kernighan-Lin algorithm optimizes the partition boundaries during the layer-by-layer refinement process to minimize boundary cut sets. A semi-coarsening strategy is adopted in the vertical mesh settings, where the number of vertical mesh layers decreases moderately with increasing water depth. The number of vertical layers is set to 10 to 15 in shallow water and 5 to 8 in deep water, maintaining vertical mesh accuracy while reducing computational load. A real-time load monitoring mechanism is established to record the number of mesh cells and computation time for each process, outputting load statistics every 100 time steps to optimize the partitioning strategy.

[0054] The specific implementation of step S03 involves using a multi-time-step algorithm to solve the wave propagation equation and the geomorphic evolution equation separately, taking into account the time scale differences. The wave propagation equation describes the propagation and evolution of waves in space, with a characteristic time scale on the order of the wave period, typically 1 to 20 seconds. The wave process time step is set to 0.01 to 0.1 seconds, and an explicit time-progression scheme, such as the fourth-order Runge-Kutta method or the Adams-Bashforth method, is used for solving it. The advantage of explicit schemes is their high computational efficiency and ease of parallel implementation. The geomorphic evolution equation describes the change in seabed elevation due to sediment transport, with a characteristic time scale of the tidal period or longer, typically several hours to several days. The geomorphic evolution time step is set to 1 to 24 hours, and an implicit time-progression scheme, such as the backward Eulerian method or the Crank-Nicolson method, is used for solving it. Implicit schemes have the characteristic of unconditional stability, allowing for larger time steps without inducing numerical instability. A conservation-type coupling interface between the wave field and the geomorphic field is established using operator splitting technology. This technique decomposes the wave propagation equation and the geomorphic evolution equation into convection operators, diffusion operators, and source operators based on physical processes. Within the wave process time step, only the convection and diffusion operators are solved to calculate the energy dissipation caused by wave propagation and wave breaking. Within the geomorphic evolution time step, wave radiation stress from multiple wave process time steps is accumulated as the input to the source operator. Wave radiation stress represents the momentum flux gradient caused by wave propagation, driving nearshore currents and sediment transport. A flux conservation condition is established at the coupling interface, requiring that the integral of wave radiation stress at the seabed surface equals the momentum change caused by sediment transport. A correction flux term is introduced to compensate for the error introduced by operator splitting; this correction flux term is explicitly estimated based on the solution from the previous time step. The two-way feedback mechanism includes topographic feedback and dynamic feedback. Topographic feedback changes the wave propagation path and wave height distribution by updating the water depth field, while dynamic feedback changes the distribution of bed shear stress and sediment transport direction by updating the wave field. This two-way coupling ensures that the interaction between wave dynamics and geomorphological evolution is accurately simulated.

[0055] The specific implementation of step S04 involves treating the mixing length coefficient, bottom roughness coefficient, and suspended sediment settling velocity coefficient in the turbulence model as random parameters. These parameters exhibit spatial variability and measurement uncertainty in the actual marine environment, with typical uncertainty ranging from 20% to 50% of the mean. A polynomial chaotic expansion model is used to probabilistically characterize the random parameters and perform uncertainty propagation analysis. The basic idea of ​​polynomial chaotic expansion is to represent random variables as a series sum of orthogonal polynomial families, with the series coefficients being deterministic quantities. The corresponding orthogonal polynomial family is selected based on the probability distribution type of the random parameters: Legendre polynomials are chosen when the parameters follow a uniform distribution, Hermitian polynomials when they follow a normal distribution, and Jacobi polynomials when they follow a beta distribution. This selection principle is based on the natural correspondence between orthogonal polynomials and probability measures, enabling the achievement of a given approximation accuracy with the fewest expansion terms. The mixing length coefficient, bottom roughness coefficient, and suspended sediment settling velocity coefficient are expanded into a finite series sum of orthogonal polynomial families. The expansion order is determined according to the required accuracy, typically ranging from order 3 to order 5. The state variables in the stochastic partial differential equations, such as flow velocity, wave height, and suspended sediment concentration, are also expanded into a series sum of the same family of orthogonal polynomials. The Galerkin projection method is then used to project the expanded stochastic partial differential equations onto the orthogonal polynomial space using an inner product. This inner product operation leverages the orthogonality of the polynomials to transform the stochastic equations into a deterministic extended system of equations with respect to the expansion coefficients. Solving this deterministic extended system yields the expansion coefficients. The statistical moments of the state variables can be directly calculated from these coefficients. The expected value equals the zeroth-order expansion coefficients, and the variance equals the sum of the squares of the higher-order expansion coefficients multiplied by the square of the norm of the corresponding orthogonal polynomial. Compared to the Monte Carlo method, polynomial chaotic expansion has a higher convergence speed, converging exponentially for smooth stochastic response functions, significantly reducing the computational cost of uncertainty quantification.

[0056] The specific implementation of step S05 involves using an algebraic multigrid preconditioner combined with a generalized minimum residual iterative method to solve the large-scale sparse linear equation system formed after discretizing the three-dimensional hydrodynamic governing equations. The hydrodynamic governing equations include continuity and momentum equations, which are discretized using the finite volume method or finite element method to obtain a linear equation system with a coefficient matrix dimension reaching [missing value]. to Direct solution methods, such as Gaussian elimination, have a computational complexity of cubed dimension, which is too costly for large-scale problems. Algebraic multigrid methods are efficient iterative solutions. Their core idea is to iterate across multiple grid levels. Fine grid layers reduce high-frequency error components, while coarse grid layers correct low-frequency error components. A coarsening operator is constructed using the Ruge-Stuben coarsening strategy to aggregate adjacent grid nodes into a coarse grid hierarchy. This strategy is based on strong connections in the coefficient matrix; a strong connection is defined when the absolute value of a matrix element between nodes is greater than 0.25 times the maximum absolute value of the off-diagonal elements in that row. An interpolation operator is constructed to define the mapping relationship from coarse to fine grid variables, with interpolation weights calculated based on the connection pattern of the coefficient matrix. A constraint operator is constructed to define the projection relationship from the fine grid residuals to the coarse grid, typically taken as the transpose of the interpolation operator. Smoothing iteration is performed at the finest grid layer to reduce high-frequency errors. This smoothing iteration uses Gaussian-Seidel iteration or incomplete LU decomposition iteration, with the number of iterations adaptively adjusted according to the residual descent rate, typically 2 to 5 times. Low-frequency errors are corrected by solving the residual equation at the coarse-grid layer. The coarse-grid correction values ​​are then passed back to the fine-grid layer via interpolation operators for post-smoothing iteration. Algebraic multigrid preconditioners improve the convergence speed of the iterative method from linear convergence dependent on grid size to geometric convergence independent of grid size, making the total computational complexity linearly related to the matrix dimension. Graphics processors (GPUs) are used to accelerate sparse matrix-vector multiplication operations. GPUs have numerous parallel computing cores, making them suitable for data-parallel operations like sparse matrix-vector multiplication, achieving speedups of 10 to 50 times.

[0057] The specific implementation of step S06 involves inputting the calculated velocity field, wave height field, and suspended sediment concentration field data into a multi-field coupled prediction model. This model, based on a deep learning architecture, can learn the complex nonlinear coupling relationships between multiple physical fields and predict the coastal evolution trend over the next 7 to 30 days. The input layer receives the values ​​of six physical field variables at spatial grid nodes: the horizontal and vertical components of the velocity field, the effective wave height and spectral peak period of the wave height field, and the surface and bottom concentrations of the suspended sediment concentration field. Spatial features are extracted through three-dimensional convolutional layers, and the convolutional kernels can capture the local spatial patterns of the physical fields. The encoder part adopts a representation learning framework based on the information bottleneck principle, containing four residual convolutional blocks. Each residual convolutional block contains two convolutional layers and one skip connection. The skip connection directly passes the input to the output, alleviating the gradient vanishing problem in deep networks. The convolutional kernel size is... The number of channels is 64, 128, 256, and 512, increasing layer by layer to extract more abstract high-level features. A batch normalization layer and a modified linear unit activation function (MLU) are connected after each residual convolutional block. The batch normalization layer standardizes the feature distribution to accelerate training convergence, while the MLU introduces nonlinear transformation capabilities. A variational coding module is connected to the encoder output, containing mean coding and variance coding branches. These branches map the feature vectors to the mean and log-variance vectors of the latent representation through fully connected layers, respectively. The latent representation dimension is set to 128. The representation learning framework based on the information bottleneck principle applies information compression constraints to the latent representation layer, forcing the encoder to extract the key features most relevant to the prediction target while filtering redundant information. Its core principle is to minimize the mutual information between the input data and the latent representation while maximizing the mutual information between the latent representation and the prediction target. Since direct computation of mutual information is infeasible in high-dimensional continuous space, variational approximation techniques are used to introduce a variational distribution to approximate the true posterior distribution. Reparameterization techniques are used to make the sampling process differentiable, thus supporting backpropagation training. The decoder section employs a mirror-symmetric structure, gradually restoring spatial resolution through transposed convolutional layers. The output layer... The convolutional layer generates predicted values ​​for coastal erosion rate, shoreline retreat distance, and damaged area in ecologically sensitive areas.

[0058] The specific implementation of step S07 involves determining whether to initiate fluid-structure interaction (FSI) calculations for the sediment-ecology module based on the predicted coastal erosion rate and shoreline retreat distance output by the multi-field coupling prediction model. The annual average erosion rate threshold is determined statistically based on historical observation data, typically ranging from 0.5 m / year to 2 m / year. When the predicted coastal erosion rate exceeds 150% of the annual average erosion rate threshold or the predicted shoreline retreat distance exceeds 5 meters, the risk of coastal erosion is deemed to have increased significantly, requiring the initiation of refined FSI calculations to assess the damage to ecologically sensitive areas. The FSI calculations involve two independent grid systems: a fluid grid and a sediment grid. The fluid grid is used to solve the hydrodynamic equations, while the sediment grid is used to solve the sediment transport equations. While there is spatial overlap, the grid division methods differ. A hypergrid-based conservation interpolation algorithm is used to transfer bed shear stress and sediment flux in the overlapping areas of the fluid and sediment grids. The core of this algorithm is the construction of hypergrid cells to achieve conservation data transfer. For each pair of overlapping fluid and sediment grid cells, the intersection of the two cell boundaries is calculated as the vertex of the supergrid cell. The area of ​​the supergrid cell is calculated using the coordinates of the polygon vertex, and the affiliation of each supergrid cell with the source and target grid cells is recorded. The distribution of bed shear stress in the fluid grid is interpolated and transferred. The area weight of all supergrid cells within each sediment grid cell is calculated, where the area weight equals the area of ​​the supergrid cell divided by the total area of ​​the sediment grid cells. The bed shear stress value of the target grid cell is equal to the sum of the products of the bed shear stress values ​​of the corresponding source grid cells and their area weights. The sediment flux in the sediment grid is interpolated and transferred, and the same area weight method is used to calculate the sediment flux received by the fluid grid cell. Conservation is verified by checking whether the flux integrals of the source and target grid cells are equal; the conservation error should be less than 1%.

[0059] The specific implementation of step S08 involves calculating the distribution of bottom shear stress and sediment resuspension flux within the ecologically sensitive area. Bottom shear stress, generated by the combined effects of water flow velocity and wave motion, is a key dynamic factor driving sediment initiation and transport. The critical initiation shear stress is determined using a Shields curve, which describes the relationship between the dimensionless critical shear stress and the particle Reynolds number. The dimensionless critical shear stress equals the critical initiation shear stress divided by the product of the underwater specific gravity of the sediment particles and the median particle size. The particle Reynolds number equals the product of the frictional velocity and the median particle size divided by the water kinematic viscosity. The frictional velocity equals the bottom shear stress divided by the square root of the water density. For sandy sediments with a median particle size of 0.1 mm to 0.5 mm, the dimensionless critical shear stress ranges from 0.03 to 0.06. When the bottom shear stress exceeds the critical initiation shear stress, it is considered an erosion state. Sediment particles overcome gravity and interparticle cohesion and begin to move, leading to a decrease in seabed elevation and destruction of ecological habitats. The number of grid cells in an eroded state within the ecologically sensitive area is counted, and their proportion of the total area of ​​the ecologically sensitive area is calculated as the ecological risk index. Sediment resuspension flux is calculated based on the difference between the bed shear stress and the critical initiation shear stress. The larger the difference, the higher the resuspension flux. The increase in resuspension flux leads to an increase in water turbidity, which affects the photosynthesis of underwater vegetation and the living environment of benthic organisms.

[0060] The specific implementation of step S09 involves classifying early warning levels based on the ecological risk index and outputting corresponding management measure recommendations. When the ecological risk index is less than 10%, it is considered a normal state, indicating that the impact of waves and currents on the ecologically sensitive area is within acceptable limits, requiring no special measures and maintaining the regular monitoring frequency. When the ecological risk index is between 10% and 30%, it is considered a state of concern, indicating that some areas show signs of erosion but have not yet posed a serious threat. It is recommended to increase the monitoring frequency from once a month to once a week to closely track the coastal evolution trend. When the ecological risk index is between 30% and 60%, it is considered an early warning state, indicating significant erosion and a further deterioration trend. It is necessary to initiate protective measure optimization calculations, adding artificial sandbars or underwater submerged dikes to the three-dimensional hydrodynamic-sediment transport coupled calculation grid, adjusting the position, elevation, and length parameters of the artificial sandbars or underwater submerged dikes, rerunning the coupled calculation to evaluate the reduction in the ecological risk index under different protection schemes, and selecting the protection scheme with the largest reduction in the ecological risk index and an engineering investment cost lower than the budget constraint. When the ecological risk index exceeds 60%, it is considered an alarm state, indicating that the ecologically sensitive area is facing a serious risk of damage. An emergency response plan needs to be issued immediately, including measures such as temporarily closing the ecologically sensitive area, implementing emergency bank protection projects, and allocating emergency supplies. At the same time, an inter-departmental coordination mechanism should be activated to ensure the rapid implementation of emergency measures.

[0061] Specifically, the principle of this invention is as follows: This invention solves the problem of balancing accuracy and efficiency in multiphysics coupling calculations because wave propagation and landform evolution have different physical timescale characteristics. Wave propagation involves short water movement response times, requiring fine time steps for capture; while landform evolution is a cumulative effect of wave action, requiring long integration times to observe significant changes. The operator splitting technique decomposes the coupling equations according to physical mechanisms, allowing different physical processes to be solved using their own suitable time steps. When transferring physical quantities across different timescales through a conservation-type coupling interface, mass and momentum conservation are ensured, avoiding the computational resource waste caused by traditional methods that are forced to use the smallest time step to meet numerical stability conditions. The corrected flux term is explicitly estimated based on the solution of the previous time step to compensate for the truncation error introduced by operator splitting, ensuring the computational accuracy of cross-timescale coupling. Therefore, the technical solution of this invention logically achieves a balance between accuracy and efficiency.

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

[0063] The specific implementation of step S01 is as follows: Wave observation buoy arrays are deployed at intervals of 200 to 500 meters in the sandy coastal monitoring area. Each buoy is equipped with a pressure-type wave height sensor and an inclination sensor, with a sampling frequency set to 2 to 4 Hz to continuously collect the effective wave height. Wave peak cycle and the angle of the main wave direction ,in The unit is meters. The unit is seconds. The unit is degrees, with a range of 0 to 360 degrees. Sediment sampling points were established in both the intertidal and subtidal zones, with sampling depths from 0 to 10 cm above the surface. Sediment grain size distribution was determined using a laser grain size analyzer, and the median grain size was calculated. The unit is millimeters. The wave energy spectrum direction decomposition and reconstruction algorithm uses the maximum entropy method to invert the two-dimensional frequency-direction spectrum function. ,in Wave frequency, measured in Hertz (Hz). The wave direction angle is expressed in radians. The unit is square meters per hertz per radian. The expression for the entropy function is: .

[0064] In the formula, It is a dimensionless entropy value; The maximum cutoff frequency is measured in Hertz (Hz) and is typically between 0.5 and 1.0 Hz. The reference spectral density, measured in square meters per hertz per radian, is set to 1 square meter per hertz per radian for dimensionless processing. This is done while satisfying observational constraints. Nonnegative constraints lower envoy Reaching the maximum value, where The observed one-dimensional spectrum is expressed in square meters per hertz. The Lagrangian function is constructed as follows: .

[0065] In the formula, Let be a Lagrange multiplier, and be a function of frequency, dimensionless; This represents the index of the discrete frequency point. Through... The variational equation is obtained to obtain the optimization equation. The nonlinear optimization equation is solved using the conjugate gradient iteration method. A Tikhonov regularization term is introduced into the objective function: .

[0066] In the formula, This is a dimensionless regularization parameter, determined using the L-curve method based on the noise level of the wave data. Its typical value range is... to ; As a reference angle, it is set to 1 radian for dimensionless processing; Let be the second partial derivative of the two-dimensional frequency-direction spectrum with respect to the direction angle, expressed in square meters per hertz per cubic radian. For a wave system with multiple wave crests, the two-dimensional frequency-direction spectrum is represented by a Gaussian mixture distribution model as follows: .

[0067] In the formula, The number of wave systems is typically taken as 2 to 4; For the first The peak amplitude of each wave system, expressed in square meters per hertz per radian. For the first The peak frequency of each wave series, in Hertz; For the first The peak direction of each wave system, in radians; For the first The frequency broadening parameter of the wave system is expressed in Hertz, with an empirical value of 0.01 to 0.05 Hz. For the first The directional broadening parameter of the wave system is expressed in radians, with an empirical value of 0.2 to 0.5 radians. For wave series numbering.

[0068] The specific implementation of step S02 is as follows: A three-dimensional hydrodynamic-sediment transport coupled computational grid is constructed, using an unstructured triangular grid in the horizontal direction and a vertical grid... The coordinate system is divided into 10 to 20 layers. A dynamic partitioning algorithm based on Hilbert curve mapping is used to determine the geometric center coordinates of the two-dimensional grid cells. Mapping to one-dimensional sequence index ,in and The unit is meters. Number the grid cells. For dimensionless integer indices. The Hilbert curve mapping formula is: ,in A recursively defined space-filling curve function ensures that adjacent spatial units remain adjacent in a one-dimensional sequence. This is based on the number of parallel computing processes. Divide the one-dimensional sequence index into consecutive segments, the first... The index range allocated to each process is ,in The total number of grid cells. The process number is a number ranging from 1 to 1. Calculate the load balancing metric after mesh refinement or coarsening: .

[0069] In the formula, It is a dimensionless load balance index; The standard deviation of the number of grid cells in each process is calculated using the following formula: ,in For the first The number of grid cells per process; The average number of grid cells for each process is calculated using the following formula: .when Incremental repartitioning is initiated at specific times. The coarsening stage of the multi-level graph partitioning algorithm employs edge matching, with the merging criterion being edge weight. ,in For grid cells and unit The edge weights between them are dimensionless and are defined as follows: , For unit and unit The length of the shared side, in meters. For unit and unit The center distance, in meters. To be with unit Number all adjacent cells. During the refinement stage, the partition boundaries are optimized using the Kernighan-Lin algorithm, with the objective function being to minimize the cut edge weights. .

[0070] In the formula, This represents the sum of dimensionless edge cut weights; and It consists of two partitions; For partitioning The grid cell number in the data; For partitioning The grid cell numbering in the data. The vertical grid uses a semi-coarsening strategy, and the surface grid spacing... and bottom grid spacing All are 0.5 to 2 meters apart, with the middle layer grid spacing being... It can be widened to 5 to 10 meters.

[0071] The specific implementation of step S03 is as follows: A multi-time-step algorithm is used to solve the wave propagation equation and the geomorphic evolution equation separately. The wave propagation equation adopts an explicit calculation format, and the wave process time step is... Setting the time to 0.01 to 0.1 seconds satisfies the CFL stability condition. ,in This is a grid scale, with the unit being meters. The group velocity, measured in meters per second, is calculated using the following formula: .

[0072] In the formula, Wave number, expressed per meter; Water depth, measured in meters; The phase velocity is expressed in meters per second, and the calculation formula is as follows: .

[0073] In the formula, The acceleration due to gravity is taken as 9.81 m / s². The landform evolution equation employs an implicit calculation scheme, with the landform evolution time step being... Set to 1 to 24 hours. Operator splitting technique applies wave propagation equations. Decomposed into convection operators Source term operator ,in This is the wave energy density vector, in joules per square meter. Time, in seconds. This is the gradient operator, expressed in meters. Wave dissipation rate, expressed in watts per square meter. Geomorphological evolution equation. Decomposed into input operators Source term operator ,in This refers to seabed elevation, in meters. The porosity of the sediment is dimensionless, with an empirical value of 0.3 to 0.5. This refers to sediment transport flux, expressed in kilograms per meter per second. This represents the wave radiation stress term, expressed in meters per second. The conservatism coupling interface is established based on the flux conservation condition: .

[0074] In the formula, is the wave radiation stress tensor, in Pascals; The density of the sediment is expressed in kilograms per cubic meter, and is taken as 2650 kilograms per cubic meter. The density of the water is expressed in kilograms per cubic meter, and is taken as 1025 kilograms per cubic meter. This represents the seabed surface area, in square meters. Corrected flux term: .

[0075] In the formula, To correct the flux term, the unit is Pascal; The number of wave time steps contained within a geomorphic evolution time step; For the first The radiation stress tensor for each wave time step, in Pascals. Number the wave time step.

[0076] The specific implementation of step S04 is: to adjust the mixing length coefficient. Bottom roughness coefficient and suspended sediment settling velocity coefficient Let be a uniformly distributed random parameter, where The unit is meters, and the value ranges from 0.01 to 0.1 meters. The unit is meters, and the range of values ​​is... to rice, The unit is meters per second, and the value ranges from 0.001 to 0.01 meters per second. A Legendre polynomial series expansion is used. For uniformly distributed random variables... Dimensionless, Legendre polynomial is defined as , The recursive relation is: .

[0077] In the formula, The order of the polynomial; for The Legendre polynomial is dimensionless. Expanding the random parameters as follows: , , ,in The order of the expansion is typically taken as 3 to 5. , , These are the expansion coefficients, with units of meters, meters per second, and meters per second, respectively. The sequence number is the expansion term number. The velocity, a state variable in the stochastic partial differential equation, is... Expanded to: .

[0078] In the formula, These are spatial coordinate vectors, with units of meters. Time, in seconds; This is the expansion coefficient vector of the flow velocity, in meters per second; The index is the number of the expansion term. The inner product projection of the expanded stochastic partial differential equation onto the orthogonal polynomial space is performed using the Galerkin projection method. ,in Pressure, measured in Pascals. This refers to the fluid density, expressed in kilograms per cubic meter. The coefficient of kinematic viscosity is expressed in square meters per second. This represents the inner product operator in the probability space. Let be the order of the projective polynomial. This yields the deterministic extended system of equations: .

[0079] In the formula, For the first The order expansion coefficient vector, in meters per second; The inner product coefficients of the product of three polynomials are dimensionless and are calculated using the following formula: ; The inner product coefficients of the binomial polynomial product are dimensionless and are calculated using the following formula: ; This is the expansion coefficient of the kinematic viscosity, expressed in square meters per second. is the expansion coefficient of pressure, in Pascals. The expected value of the state variable is... The unit is meters per second, and the variance is: .

[0080] In the formula, The unit is square meters per second squared; The norm square of the Legendre polynomial is dimensionless and is calculated using the following formula: .

[0081] The specific implementation of step S05 is as follows: The sparse linear equation system formed after discretization of the three-dimensional hydrodynamic control equations is solved using an algebraic multigrid preconditioner combined with the generalized minimum residual iteration method. ,in The coefficient matrix, The vector represents the unknown solution, and its unit is meters per second. The vector is the right-hand side term, in meters per second. Algebraic multigrid preconditioners construct coarsening operators. and interpolation operators The operator is restricted to ,in Number the grid levels. This represents the matrix transpose. The Ruge-Stuben coarsening strategy determines strong connections based on the following criteria: ,in Coefficient matrix The Line number Column elements, and Number the grid nodes. For nodes All connected nodes are numbered, and strongly connected node pairs are aggregated into coarse mesh nodes. The V-loop iteration process is as follows: at the finest mesh layer... Perform smooth iteration To reduce high-frequency errors, among which For the number of iterations, For the first The solution vector for the nth iteration, in meters per second. When using Gauss-Seidel iteration as a smoothing operator... , It is a diagonal matrix. It is a lower triangular matrix. Calculate the residuals: .

[0082] In the formula, For the first The residual vector of the layer mesh, in meters per second; For the first The coefficient matrix of the layer mesh; For the first The solution vector of the layered mesh, in meters per second. Projected onto a coarse mesh by a constraint operator. Solving the residual equation on a coarse mesh ,in This is the coarse grid correction value vector, in meters per second. For the first The coefficient matrix of the layer mesh, For the first The residual vector of the layer mesh, in meters per second. Correction values ​​are transferred through interpolation operators. ,in For the first The correction vector of the layer mesh, in meters per second, updates the solution. And perform post-smoothing iteration. The generalized minimum residual iteration method in the Krylov subspace... To find an approximate solution, we construct orthogonal basis vectors through Arnoldi iteration and minimize the residual norm. ,in Describes the Euclidean norm operator. For the first The approximate solution vector obtained through iterative steps, measured in meters per second. Graphics processors accelerate sparse matrix-vector multiplication. It is stored in CSR format, where and It is a vector, with units of meters per second. The CSR format contains three arrays: Store non-zero element values. Store column indexes, Store the starting position of each row. Each thread calculates the dot product of a row: .

[0083] In the formula, For output vector The Each component is measured in meters per second. For input vector The Each component is measured in meters per second. Coefficient matrix The One non-zero element value; and For the first Walking The start and end indices of the array.

[0084] The specific implementation of step S07 is as follows: when the predicted coastal erosion rate value Exceeding the annual average erosion rate threshold 150% of Predicted value of shoreline setback distance More than 5 meters At meter , initiate fluid-structure interaction calculations for the sediment-ecology module, where and All units are meters per year. The unit is meters. A conservation interpolation algorithm based on hypergrids is used in fluid mesh cells. and sediment grid cells Constructing supermesh cells in overlapping regions superscript Indicates the first A fluid mesh cell, superscript Indicates the first One sediment grid cell. Hypergrid cell area: .

[0085] In the formula, The area of ​​the hypergrid cell is in square meters. Subbed shear stress. The interpolation formula for transferring data from the fluid grid to the sediment grid is: .

[0086] In the formula, For the first The bed shear stress of a fluid grid cell, in Pascals; For the first The bed shear stress of a sediment grid cell, in Pascals; To be compatible with sediment grid cells The number of overlapping fluid mesh cells; The area weight is dimensionless, and the calculation formula is: .

[0087] In the formula, For the first Total area of ​​each sediment grid cell, in square meters. Sediment flux. The interpolation formula for transferring data from the sediment grid to the fluid grid is: .

[0088] In the formula, For the first Sediment flux per sediment grid cell, in kilograms per meter per second; For the first Sediment flux per fluid grid cell, in kilograms per meter per second; For use with fluid mesh elements The number of overlapping sediment grid cells; The area weight is dimensionless, and the calculation formula is: .

[0089] In the formula, For the first The total area of ​​the fluid mesh cells is expressed in square meters. The conservation verification condition is as follows: .

[0090] In the formula, This represents the total number of fluid mesh cells; This represents the total number of sediment grid cells.

[0091] The specific implementation of step S08 is as follows: Calculate the distribution of subgrade shear stress within the ecologically sensitive area. And sediment resuspension flux: .

[0092] In the formula, For spatial location The subbed shear stress at the point is expressed in Pascals. The resuspension flux is expressed in kilograms per square meter per second. The resuspension factor is expressed in kilograms per square meter per second per pascal, with an empirical value of [value missing]. to kilograms per square meter per second per pascal; Critical starting shear stress, measured in Pascals. The dimensionless critical shear stress is determined based on the Shields curve: .

[0093] In the formula, It is a dimensionless critical shear stress; This refers to the density of sediments, expressed in kilograms per cubic meter. This refers to the density of water, expressed in kilograms per cubic meter. This is the acceleration due to gravity, measured in meters per second squared. Median particle size, in meters. Particle Reynolds number: .

[0094] In the formula, The dimensionless particle Reynolds number; The frictional velocity is expressed in meters per second, and the calculation formula is as follows: .

[0095] In the formula, The viscosity coefficient of water kinematics, in square meters per second, is assumed to be... Square meters per second. For median particle size It consists of sandy sediments ranging from 0.1 to 0.5 mm in size. The value ranges from 0.03 to 0.06. When The site was determined to be in a state of erosion, with an ecological risk index of: .

[0096] In the formula, It is a dimensionless ecological risk index; This represents the sum of the areas of the eroded grid cells, in square meters. This represents the total area of ​​the ecologically sensitive zone, expressed in square meters.

[0097] The specific implementation method of step S09 is the same as described above, and will not be repeated in detail here.

[0098] When a warning state is triggered, the optimization calculation of protective measures adds artificial sandbars or underwater submerged dikes to the three-dimensional hydrodynamic-sediment transport coupled calculation grid. The impact of the protective structure on wave propagation is determined by modifying the wave energy equation. Implementation, in which Wave energy density, expressed in joules per square meter. The dissipation term for the protective structure, expressed in watts per square meter, is calculated using the following formula: .

[0099] In the formula, The dissipation factor, measured in seconds, is related to the elevation of the protective structure. and length The relevant calculation formula is as follows: .

[0100] In the formula, This is a dimensionless dissipation parameter, with an empirical value of 0.1 to 0.5; The elevation of the protected structure is shown in meters. The length of the protective structure is in meters; Water depth at the protected structure, in meters. Assess the reduction in ecological risk index under different protection schemes: .

[0101] In the formula, The reduction rate of the dimensionless ecological risk index; This is the ecological risk index without protective measures, and is dimensionless. To adopt the first The ecological risk index following the No. 1 protection plan is dimensionless. Number the protection scheme. Select the one that meets the requirements. And engineering investment costs The protection scheme is output as the optimization result, in which For the first The investment cost of the protection solution is in yuan. Due to budget constraints, the unit is yuan.

[0102] To better understand and implement this invention, the following is a specific application scenario example 2: To verify the effectiveness of this invention, technicians built a test environment combining actual data and simulation, selecting a sandy coastal monitoring area as the research object. This area has a coastline length of approximately 12 kilometers, with nearshore water depths ranging from 0 to 15 meters, and distinct intertidal and subtidal ecologically sensitive zones. Technicians deployed five wave observation buoys in the monitoring area, spaced approximately 2.5 kilometers apart, and simultaneously set up 15 sediment sampling points, covering the intertidal zone, subtidal zone, and offshore shallow waters. The observation period selected was a typhoon event in October 2024, lasting 72 hours, including continuous observation data from 48 hours before the event and 72 hours after its conclusion.

[0103] The data collected by technicians through the wave observation buoy array is shown in Table 1.

[0104] Table 1. Statistical Table of Wave Observation Data

[0105] Sediment sampling and analysis showed that the median grain size in the intertidal zone was 0.28 mm, belonging to the medium sand category, while the median grain size in the subtidal zone was 0.15 mm, belonging to the fine sand category. The sediments were well sorted, with a grain size distribution standard deviation of 0.12 mm. Technicians used a wave energy spectrum direction decomposition and reconstruction algorithm to process the wave data in a limited number of directions. A two-dimensional frequency-direction spectrum function was obtained through maximum entropy inversion. A Tikhonov regularization term was introduced into the objective function, and the regularization parameter was determined to be 0.08 using the L-curve method, successfully suppressing numerical oscillations in high-frequency directions.

[0106] Technicians constructed a three-dimensional coupled hydrodynamic-sediment transport computational grid. Horizontally, an unstructured triangular grid with a total of 45,000 cells was used, with localized refinement in shallow near-shore waters and areas with dramatic topographic changes. The minimum grid size was 8 meters, while the grid size in deep offshore waters was 50 meters. Vertically, the grid was divided into 10 layers using a Sigma coordinate system. The surface and bottom layers had grid thicknesses of 5% and 8% of the water depth, respectively, while the intermediate layers were proportionally distributed. Figure 3 As shown, the water depth topography of the monitoring area exhibits a distinct coastal sandbar-valley system. The water depth at the top of the sandbars is approximately 3 meters, and the water depth at the bottom of the valleys is approximately 8 meters. The topographic slope varies between 1:50 and 1:200. Technicians employed a dynamic partitioning algorithm based on Hilbert curve mapping to distribute the grid across 16 parallel computing processes. Hilbert curve mapping ensured the adjacency of spatially adjacent grid cells in the one-dimensional sequence index. The initial load balancing index was 0.12, meeting the load balancing requirements.

[0107] Technicians used a multi-time-step algorithm to solve the wave propagation equation and the geomorphic evolution equation, such as Figure 2 As shown, the wave propagation equation uses an explicit calculation format with a time step of 0.05 seconds, while the geomorphic evolution equation uses an implicit calculation format with a time step of 6 hours. A conserved coupling interface is established using operator splitting technology to decompose the wave propagation equation into convection and diffusion operators. Wave radiation stress for 120 wave process time steps is accumulated as the source term input within each geomorphic evolution time step. The mixing length coefficient in the turbulence model is set as a uniformly distributed random parameter, ranging from 0.35 to 0.45; the bottom roughness coefficient ranges from 0.018 to 0.028; and the suspended sediment settling velocity coefficient ranges from 0.8 to 1.2. A Legendre polynomial series expansion is performed on the random parameters using a polynomial chaotic expansion model, with the expansion order set to 3. The stochastic partial differential equations are then transformed into a deterministic extended equation set using the Galerkin projection method. The extended equation set is 20 times larger than the original equation set.

[0108] Technicians used an algebraic multigrid preconditioner combined with a generalized minimum residual iteration method to solve the sparse linear equation system formed after discretization of the three-dimensional hydrodynamic control equations. The coefficient matrix had 1.8 million non-zero elements and a matrix size of 450,000 x 450,000. A four-layer coarse-grid hierarchy was constructed using the Ruge-Stuben coarsening strategy, with the number of nodes in the coarse-grid layers being 110,000, 28,000, 7,000, and 1,800, respectively. The strong connectivity threshold was set to 0.25. Three Gauss-Seidel smoothing iterations were performed in the finest grid layer, and the residual equations were solved in the coarse-grid layer. The iteration convergence criterion was that the residual norm decreased to a certain percentage of the initial residual. The average number of iterations is 15, and the time per iteration is 0.8 seconds. Utilizing a graphics processing unit (GPU) to accelerate sparse matrix-vector multiplication operations results in approximately 8 times the computational efficiency compared to a CPU implementation.

[0109] like Figure 4 As shown, the nearshore wave field of a sandy coast simulated by the Spectral Waves FM module reveals that during a typhoon's passage, the significant wave height in shallow nearshore waters decreases from 4.8 meters to 2.5 meters due to shallow water deformation and bottom friction. The wave breaking zone is located in the area with a water depth of 2.5 to 3.5 meters, and the wave propagation direction is deflected due to topographic refraction. Figure 5 As shown, the Flow MODEL FM module simulation of nearshore storm surge on sandy coasts shows that the maximum water level increase occurred 6 hours after the typhoon passed, with a nearshore water level increase of 1.2 meters. The superposition of the storm surge and the astronomical tide caused the tide level to exceed the warning level by 0.8 meters. Figure 6As shown, the nearshore current results simulated by the Flow MODEL FM module for sandy coasts show that during storms, the coastal current velocity reaches 0.8 to 1.2 meters per second, the offshore current velocity reaches 0.6 meters per second in the trough area, and the bottom current velocity is about 30% to 40% lower than the surface current velocity.

[0110] Technicians input the calculated velocity field, wave height field, and suspended sediment concentration field data into a multi-field coupled prediction model. The input data includes the values ​​of six physical field variables—horizontal velocity component, vertical velocity component, significant wave height, spectral peak period, surface suspended sediment concentration, and bottom suspended sediment concentration—across 45,000 grid nodes. The multi-field coupled prediction model employs a representation learning framework based on the information bottleneck principle. The encoder extracts spatial features through four layers of residual convolutional blocks, with a latent representation dimension of 128. The mean vector and log-variance vector output by the variational coding module are used for uncertainty quantification. The model training dataset contains 15 historical typhoon events, with a total of 15 training samples, divided into a training set of 12 and a validation set of 3 at an 8:2 ratio. Model training uses a batch size of 8, an initial learning rate of 0.001, Adam optimizer momentum parameters of 0.9 and 0.999, and the weight coefficients for the reconstruction loss term, mutual information regularization term, and prediction loss term in the composite loss function are 1.0, 0.5, and 2.0, respectively. After 80 training rounds, the prediction error on the validation set converged to a minimum. The model output predicted a coastal erosion rate of 0.18 meters per day, a shoreline retreat distance of 5.4 meters, and a damaged area of ​​1.2 hectares for the next 30 days.

[0111] Because the predicted coastal erosion rate exceeded the annual average erosion rate threshold by 150%, technicians initiated fluid-structure interaction calculations for the sediment-ecological module. A sediment grid was constructed based on the fluid grid, with a total of 30,000 grid cells, slightly larger in scale than the fluid grid. A hypergrid-based conservation interpolation algorithm was used to transfer bed shear stress and sediment flux in the overlapping area of ​​the fluid and sediment grids. The hypergrid had 72,000 cells, and conservation interpolation was achieved using an area-weighted method, with flux integration error before and after interpolation less than 0.5%. Technicians calculated the bed shear stress distribution within the ecologically sensitive area and determined the critical initiation shear stress for sandy sediments with a median grain size of 0.28 mm to be 0.8 Pascals based on the Shields curve. Statistical results showed that grid cells with bed shear stress exceeding the critical initiation shear stress accounted for 42% of the total area of ​​the ecologically sensitive area, indicating an erosion state and an ecological risk index of 42%.

[0112] like Figure 7As shown, the ECO Lab module simulation results of nearshore biological oxygen demand (BOD) on sandy coasts indicate that storm disturbances caused resuspension of bottom sediments, releasing the carried organic matter into the water. BOD in the nearshore area increased from the normal level of 2.5 mg / L to 4.8 mg / L, indicating a significant increase in the rate of oxygen consumption in the water. Figure 8 As shown in the figure, the nearshore chemical oxygen demand (COD) results simulated by the ECO Lab module show that the COD increased from the normal level of 8 mg / L to 15 mg / L, exceeding the water quality standard limit, indicating that sediment resuspension has had a significant impact on the aquatic environment.

[0113] Based on an ecological risk index of 42%, falling within the 30% to 60% range, technicians determined it to be a warning state and initiated optimization calculations for protective measures. They added an underwater submerged breakwater to the three-dimensional hydrodynamic-sediment transport coupled calculation grid. The breakwater was located 800 meters offshore, at an elevation 1.5 meters below sea level, and 4 kilometers long. After rerunning the numerical simulation, the underwater submerged breakwater reduced nearshore wave height by 0.6 meters, significantly weakened the distribution of bottom shear stress, and lowered the ecological risk index to 22%, reaching a state of concern, demonstrating a significant protective effect.

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

Claims

1. A multi-factor coupled ecological early warning method for sandy coastlines based on numerical simulation, characterized in that, Wave observation buoy arrays and sediment sampling points were deployed in the sandy coastal monitoring area to collect wave parameter and sediment grain size distribution data. A two-dimensional frequency-directional spectrum was obtained using a wave energy spectrum direction decomposition and reconstruction algorithm as the far-field boundary condition. A three-dimensional hydrodynamic-sediment transport coupled computational grid was constructed, and a dynamic partitioning algorithm based on Hilbert curve mapping was used to locally refine the unstructured grid in complex terrain areas. A multi-time-step algorithm was used to solve the wave propagation equation and the geomorphic evolution equation. An explicit calculation format for the wave process time step was set for the wave propagation equation, and an implicit calculation format for the geomorphic evolution time step was set for the geomorphic evolution equation. A conserved coupling interface was established using operator splitting technology. Turbulence model parameters were input as random parameters into a polynomial chaotic expansion model for Legendre polynomial series expansion, and Galerkin projection was used. The method transforms the equations into a deterministic extended system for solving. An algebraic multigrid preconditioner combined with a generalized minimum residual iterative method is used to solve the sparse linear system of equations formed after discretization of the three-dimensional hydrodynamic control equations. The calculated velocity field, wave height field, and suspended sediment concentration field data are input into a multi-field coupled prediction model to output predicted values ​​of coastal erosion rate, shoreline retreat distance, and damaged area in ecologically sensitive areas. When the predicted values ​​exceed the threshold, the fluid-structure interaction calculation of the sediment-ecology module is initiated. A hypergrid-based conservation interpolation algorithm is used to transfer the bed shear stress and sediment flux, and the bed shear stress distribution in the ecologically sensitive area is calculated. When the bed shear stress exceeds the critical starting shear stress, it is determined to be in an erosion state. The proportion of erosion state grid cells to the total area of ​​the ecologically sensitive area is statistically analyzed as the ecological risk index, and the early warning level is classified according to the ecological risk index.

2. The method according to claim 1, characterized in that, The steps of the wave energy spectrum direction decomposition and reconstruction algorithm are as follows: using the finite-direction wave data collected by the wave observation buoy array as the observation constraint, using the maximum entropy method to invert the two-dimensional frequency-direction spectrum function, constructing the Lagrangian function to transform the constrained optimization problem into an unconstrained problem, and solving the nonlinear optimization equation through the conjugate gradient iteration method.

3. The method according to claim 2, characterized in that, The entropy function of the maximum entropy method is expressed as the integral of the two-dimensional frequency-direction spectrum function over the entire domain after taking the negative logarithm of the integral over the frequency-direction space.

4. The method according to claim 3, characterized in that, Tikhonov regularization terms are introduced into the objective function to suppress numerical oscillations in high-frequency directions. The Tikhonov regularization term is obtained by adding an integral term of the second derivative of the two-dimensional frequency-direction spectrum function with respect to the direction angle to the objective function. The smoothness is controlled by multiplying the integral term by a regularization parameter. The regularization parameter is adaptively determined by the L-curve method based on the noise level of the finite-direction wave data.

5. The method according to claim 4, characterized in that, The steps of the dynamic partitioning algorithm based on Hilbert curve mapping are as follows: the geometric center coordinates of the two-dimensional unstructured grid cells formed after local densification of the unstructured grid are mapped to the one-dimensional sequence index of the Hilbert curve, and the one-dimensional sequence index is divided into continuous segments according to the number of parallel computing processes and allocated to each process.

6. The method according to claim 5, characterized in that, After mesh refinement or coarsening, the one-dimensional sequence index of the Hilbert curve is recalculated and load balancing detection is triggered.

7. The method according to claim 6, characterized in that, Multi-level graph segmentation algorithms include a coarsening stage, an initial segmentation stage, and a refinement stage.

8. The method according to claim 7, characterized in that, The steps of establishing a conserved coupling interface using operator splitting technology are as follows: the wave propagation equation and the landform evolution equation are decomposed into convection operators, diffusion operators and source term operators according to the physical process. Within the wave process time step of solving the wave propagation equation, only the convection operators and diffusion operators are solved. Within the landform evolution time step of solving the landform evolution equation, the wave radiation stress of multiple wave process time steps is accumulated as the input of the source term operator.

9. The method according to claim 8, characterized in that, This includes establishing a two-way feedback mechanism between the wave field and the geomorphic field at the conserved coupling interface by using flux conservation conditions.

10. The method according to claim 9, characterized in that, The steps of performing Legendre polynomial series expansion on random parameters in the polynomial chaotic expansion model are as follows: the mixed length coefficient, bottom roughness coefficient, and suspended sediment settling velocity coefficient are expressed as probability distribution functions, and the orthogonal polynomial family corresponding to the probability distribution function is selected as the basis function.

Images