Inversion method for complex spatially distributed eddy viscosity coefficients in intratidal ocean models

Abstract

This invention discloses an inversion method for complex spatially distributed eddy viscosity coefficients in marine internal tide models, belonging to the field of marine numerical simulation technology. Key improvements of this invention include: the introduction of "independent point interpolation parameterization," a decisive solution to the computational efficiency bottleneck, transforming the problem from an "extremely difficult to compute" high-dimensional state to a "computable" low-to-medium-dimensional state, providing a feasible basis for subsequent operations; and the adoption of a "spectral projection gradient algorithm" as the optimization engine, further improving efficiency on the basis of dimensionality reduction, ensuring that the optimal solution is approximated with as few iterations as possible in each iteration step, thereby translating the advantages of dimensionality reduction into significant savings in actual computation time. This invention achieves a synergistic leap in accuracy, efficiency, and applicability, not only providing a powerful new tool for internal tide research but also potentially promoting the development of high-precision marine environmental forecasting and climate models, possessing significant scientific value and application potential.

Description

Technical Field

[0001] This invention belongs to the field of marine numerical simulation technology, specifically relating to an inversion method for the eddy viscosity coefficient of complex spatial distribution in an internal ocean tide model. Background Technology

[0002] I. The Importance of Internal Moisture Simulation and Eddy Viscosity Coefficient

[0003] Internal tides, tidal fluctuations occurring on density gradients within the ocean, are a crucial link in global ocean energy transport and dissipation. They arise when barotropic tides encounter seafloor topography (such as seamounts and mid-ocean ridges), transferring enormous tidal energy from large-scale to mesoscale levels, and ultimately dissipating it to smaller scales through turbulent mixing. This process plays a vital role in maintaining the meridional overturning circulation of the global ocean, influencing the climate system, and regulating nutrient transport in marine ecosystems. Therefore, accurately simulating the generation, propagation, and dissipation of internal tides is one of the core topics of modern physical oceanography.

[0004] In numerical models, the eddy viscosity coefficient is a core parameter characterizing turbulent momentum exchange and energy dissipation, and its value and spatial distribution directly determine the accuracy of endorheic simulations. Specifically:

[0005] (1) Key parameters of energy dissipation: During propagation, the internal tides experience energy attenuation due to mechanisms such as shear instability and topographic friction. The eddy viscosity coefficient is the core physical quantity characterizing this momentum exchange and energy dissipation rate. It quantifies the momentum transport effect of turbulence or vortices on the mean flow, and directly determines the energy attenuation rate, mixing efficiency, and the final propagation distance of the internal tides. An inaccurate eddy viscosity coefficient will lead to the simulated internal tide amplitude being too strong or too weak, and the energy distribution will deviate significantly from the actual situation.

[0006] (2) Controlling the propagation and breakup of internal tides: Internal tides, especially baroclinic internal tides, may break up along their propagation path, triggering intense turbulent mixing. The magnitude and spatial distribution of the eddy viscosity coefficient directly affect the stability of internal tides. In the model, it acts as a "valve" for the transfer of energy from an resolvable wave scale to an unresolvable turbulent scale. If this coefficient is not set properly, the location and intensity of the breakup of internal tides cannot be accurately reproduced, leading to misjudgments of internal ocean mixing.

[0007] (3) Model results are highly sensitive to coefficients: Ocean numerical models, especially three-dimensional non-hydrostatic models or specific parameterization schemes used for endorheic studies, are extremely sensitive to the values ​​of eddy viscosity coefficients in their simulation results (such as vertical velocity structure, energy flux, mixing rate, etc.). Small changes in coefficients can lead to orders of magnitude differences in simulation results. Therefore, the accuracy of this coefficient is one of the keys to the success or failure of endorheic simulations.

[0008] II. Existing Inversion Algorithms for Eddy Viscosity Coefficient and Research

[0009] For this key parameter, researchers primarily employ the adjoint assimilation method for inversion studies. This method optimizes unknown parameters under the constraints of a numerical model by assimilating observed data. Its core idea is to use an ocean dynamic model as a constraint to invert unobservable ocean elements (such as eddy viscosity coefficients) by assimilating observable ocean elements. For example:

[0010] (1) Gao Xiumin (2010) applied the accompanying assimilation method to simulate the tidal wave in the Bohai Sea. By assimilating the harmonic constant data of 19 tide gauge stations along the Bohai Sea coast, she optimized the bottom friction coefficient and horizontal eddy viscosity coefficient (assumed to be constant) in the tidal wave model, thus improving the simulation results of each tidal constituent.

[0011] (2) Chen Haibo (2012) established an isodense coordinate end-tidal assimilation model and conducted an inversion experiment on the horizontal eddy viscosity coefficient on a narrow, ideal terrain (such as a parabolic ridge), successfully reconstructing the preset simple spatial distribution. This shows that the assimilation method has a certain inversion capability for the eddy viscosity coefficient under well-controlled ideal conditions. At the same time, this study also applied the model to end-tidal simulation, and by assimilating T / P altimeter data, optimized the open boundary conditions, bottom friction coefficient, interlayer friction coefficient, and horizontal eddy viscosity coefficient of each layer, thus improving the accuracy of end-tidal simulation.

[0012] III. Major Problems and Shortcomings of Existing Technologies

[0013] Despite initial progress in previous studies, existing techniques still have significant limitations in characterizing the complex spatial distribution of eddy viscosity coefficients, mainly in the following aspects:

[0014] (1) Inversion capability is limited to simple scenarios: Most current studies, especially in the field of internal tide simulation (e.g., Chen Haibo, 2012), have been successfully verified on highly idealized terrains (e.g., parabolic ridges, linear slopes). These terrains are far removed from the complex terrains such as seamounts, canyons, and ridges that exist in the real ocean. At the same time, the distribution of inverted coefficients is mostly constant or simple linear, which cannot effectively describe the highly complex, non-uniform, three-dimensional heterogeneous eddy viscosity coefficient field formed by the interaction of complex terrain, stratification, and multi-scale dynamic processes in the real ocean.

[0015] (2) Insufficient application and verification for complex real terrain: Although Gao Xiumin (2010) applied it in the Bohai Sea, the applicability and effectiveness of the existing method in real sea areas with extremely complex seabed topography and multi-scale dynamic processes still need to be verified. Extending the inversion algorithm from ideal terrain to real complex terrain and accurately inverting the corresponding complexity coefficient distribution is the main challenge we are currently facing.

[0016] In summary, existing research based on methods such as adjoint assimilation has made valuable explorations in the inversion of eddy viscosity coefficients, and some studies have been able to invert spatially distributed coefficients under ideal topography and simple distribution assumptions. However, they have generally failed to effectively achieve high-precision and stable inversion of eddy viscosity coefficients with fine three-dimensional spatial distribution characteristics under the complex environment of the real ocean. This lack of ability to describe key parameters has become a bottleneck restricting further improvement in the accuracy of endocytic simulations. Summary of the Invention

[0017] The purpose of this invention is to provide an inversion method for the complex spatial distribution of eddy viscosity coefficient in an internal oceanic tide model. This method can efficiently and accurately invert the fine three-dimensional spatial distribution of eddy viscosity coefficient under complex real oceanic environments, thereby fundamentally breaking through the current accuracy bottleneck of internal tide simulation.

[0018] To achieve the above objectives, the present invention specifically needs to solve the following key technical problems:

[0019] (1) How to improve the ability to characterize complex spatial distributions: Develop a new inversion framework that can break through the assumption of "constant" or "simple distribution" and directly reconstruct the eddy viscosity coefficient field with complex spatial structure from the observation data, so as to more realistically reflect the turbulent dissipation process inside the ocean.

[0020] (2) How to ensure the effectiveness and robustness of the method in real complex terrain: make the algorithm not only applicable to ideal terrain, but also directly and stably applied to sea areas with real complex terrain, and obtain physically reasonable inversion results with significantly improved accuracy.

[0021] (3) How to achieve synergistic optimization of computational efficiency and inversion accuracy: Design efficient numerical optimization strategies to ensure high-precision inversion of complex distributions while keeping computational costs within an acceptable range, thus solving the dilemma of existing methods being "fast but inaccurate" or "accurate but difficult to calculate".

[0022] In summary, the core motivation of this invention lies in bridging the significant gap between existing inversion technology capabilities and the real-world requirements of ocean physics. By providing an advanced algorithm capable of efficiently and accurately inverting complex spatially distributed eddy viscosity coefficients, this invention will provide indispensable key technical support for the accurate tracking of internal tidal energy dissipation paths, reliable estimation of ocean mixing rates, and improvement of global climate models.

[0023] To achieve the above objectives, the present invention adopts the following technical solution:

[0024] The core concept of this invention lies in constructing an efficient and high-precision inversion framework to solve the problem of inverting spatially complex eddy viscosity coefficients by assimilating ocean observation data (such as sea surface height, water temperature, salinity, or current fields). To overcome the high computational cost of traditional gradient definition methods in calculating the gradient of the objective function, this invention innovatively employs a dual strategy of "dimensionality reduction" and "accelerated convergence." Specifically, the independent point scheme significantly reduces the dimensionality of the optimization problem while maintaining the spatial complexity of the inversion parameters; simultaneously, the spectral projection gradient (SPG) algorithm is used to efficiently utilize gradient information and greatly accelerate the convergence speed of the optimization process. These two techniques work synergistically to ensure the computational feasibility and inversion efficiency of this invention when dealing with complex real-world ocean applications.

[0025] The eddy viscosity coefficient inversion method provided by this invention consists of the following steps:

[0026] S1: Construct a forward numerical model of the ocean that includes internal tidal dynamic processes and acquire observational data to be assimilated, such as satellite altimeter data and anchored current velocities;

[0027] S2: Parameterization and Dimensionality Reduction Based on Independent Point Scheme: When defining the spatial distribution of the eddy viscosity coefficient across the entire computational domain, instead of directly inverting the value of each grid point, a set of independent control points—far fewer than the total number of model grids—is arranged within the computational domain. The coefficient values ​​at these independent points are the core parameters to be inverted. Through a pre-defined spatial interpolation method (such as bilinear interpolation, Cressman interpolation, etc.), the continuous and complex spatial distribution across the entire field is reconstructed from the coefficient values ​​of the independent points.

[0028] S3: Calculate the gradient of the objective function using the gradient definition method: This is performed in the reduced-dimensional parameter space established in S2, defining an objective function that measures the difference between the model simulation results and the observed data; then, using the gradient definition method, the gradient of the objective function with respect to the coefficients of all independent points is calculated by perturbing the coefficient values ​​at each independent point and running the forward model.

[0029] S4: Apply the SPG algorithm for efficient iterative optimization: Input the gradient calculated in S3 into the SPG algorithm; The SPG algorithm uses the current iteration point and its gradient information, combined with historical iteration information, to intelligently determine the search direction and step size, and applies physical constraints on the range of values ​​of independent point coefficients in the process.

[0030] S4 is the solver for the inversion problem. It receives the gradient provided by S3 and performs an iterative search in the reduced-dimensional parameter space defined by S2. Finally, through the coordinated execution of the above four steps, an eddy viscosity coefficient field with a complex spatial distribution that can fit the observation data with high accuracy is output.

[0031] Furthermore, in S1, the forward model is used to simulate ocean conditions and parameterize the eddy viscosity coefficient; the observational data provides constraints for the inversion, ensuring that the inversion results conform to actual ocean phenomena. The construction of the ocean numerical forward model is as follows: This model is a nonlinear primitive equation model, using isodense coordinates in the vertical direction and spherical coordinates in the horizontal direction, considering the hydrostatic approximation and the free sea surface. The model employs internal and external mode separation techniques. The horizontal gradient term of sea surface atmospheric pressure is ignored, and it is assumed that the layered structure is horizontally uniform, i.e., the initial thickness of each layer except the bottom layer is... The horizontal gradients of all layers are zero. Let the density of each layer be... Dynamic thickness is Therefore, the dynamic mass is Taking tidal force as and Thus, the internal modal governing equations are obtained:

[0032] Layer 1 (Surface):

[0033] ;

[0034] ;

[0035] ;

[0036] No. layer( ):

[0037] ;

[0038] ;

[0039] ;

[0040] No. Layer (bottom layer):

[0041] ;

[0042] ;

[0043] ;

[0044] in, For time, and These are geographical longitude and latitude, respectively, and are unit vectors. , and East, north, and upward are respectively considered positive. and They are respectively and The velocity component in the direction; The radius of the Earth; For Coriolis parameters, where This represents the magnitude of the Earth's angular velocity of rotation. It is the acceleration due to gravity; The horizontal eddy viscosity coefficient; the Laplace operator is specifically expressed as:

[0045] ,

[0046] .

[0047] The vertical eddy viscosity term can be written in the form of interlayer friction:

[0048] , ,

[0049] in, For the first Layer and first The average density of the layer, For the first Layer and first The average initial thickness of the layer, is the interlayer friction coefficient. The bottom friction term in the momentum equation for the bottom layer (bottom layer) takes the following form:

[0050] ,

[0051] in The coefficient of friction is the lowest.

[0052] The external modal equations are obtained by vertically integrating the internal modal equations; firstly, the barotropic flow is defined as the mass-weighted average of the velocities in each layer, i.e.

[0053] , , ;

[0054] The momentum equations of each layer of the internal mode are... By weighted averaging and directly adding the continuity equations of each layer, the external modal equations can be obtained:

[0055] ;

[0056] ;

[0057] ;

[0058] in,

[0059] ,

[0060] ,

[0061] ,

[0062] ,

[0063] For the general still water depth, The vertical average density of the water column. and The first term on the right is external pressure, caused by the undulations of the sea surface; the second term is internal pressure, caused by the horizontal density gradient of the fluid in each layer. and The first term on the right is the mass-weighted average of the convection terms for each layer, the second term is the mass-weighted average of the turbulent viscous forces for each layer, and the third term is the difference between the local acceleration of the barotropic flow and the mass-weighted average of the local acceleration of the velocity for each layer.

[0064] Furthermore, step S2 serves as a prerequisite and preparation for subsequent optimization problems. By reducing the number of parameters to be optimized, it directly alters the dimensionality of the optimization problem in S4. Its core function is to "reduce the dimensionality" of the inversion problem. It transforms a high-dimensional problem requiring the solution of tens of thousands or even millions of grid parameters into a medium- to low-dimensional problem requiring the optimization of only tens to hundreds of independent point parameters, thereby significantly reducing the total computational cost required for subsequent gradient calculations and optimization iterations. Simultaneously, it maintains the continuity and complexity of the coefficient space distribution through interpolation.

[0065] Furthermore, in S2, considering that in practical applications, available observational data may be scarce while the number of control variables to be determined is large, which could lead to the ill-conditioned nature of the inverse problem, we assume that the control variables are spatially distributed. To achieve this spatial distribution, an "independent point scheme" is adopted in the optimization process. This involves selecting some points in the computational domain as independent points, and the spatial distribution of the control variables is obtained by linear interpolation of the values ​​at these independent points. This reduces the number of control variables to ensure the well-posedness of the inverse problem and physically maintains the continuity and smoothness of the spatial distribution of the control variables.

[0066] Furthermore, the specific idea of ​​S2 is as follows: select some grid points in the computational domain (denoted as...). As an independent point, let For the control variable values ​​at independent points, calculate the values ​​at each grid point in the calculation region (denoted as...). The value of ) Depend on Obtained through linear interpolation, and The following relationship must be satisfied:

[0067] ,

[0068] in, These are weighting coefficients, for example, they can be expressed in the Cressman form, i.e.:

[0069] ,

[0070] in, For grid and The center distance between them The radius of influence. The cost function can be derived. about The gradient is:

[0071] .

[0072] Furthermore, S3 relies heavily on the parameterization system established in S2 and provides crucial gradient information input for the optimization algorithm in S4. Although the gradient definition method itself has a high computational cost, it becomes acceptable based on the dimensionality reduction in S2. This step ensures the universality and numerical accuracy of gradient calculation, pointing the optimization algorithm in the correct direction of objective function descent.

[0073] Furthermore, the specific implementation of the gradient definition method in S3 involves estimating the spatially distributed horizontal eddy viscosity coefficient by assimilating satellite altimeter data. First, a cost function is constructed based on simulated and observed water levels:

[0074] ;

[0075] in This refers to the index of the isodense water layer where the internal modal variables reside. For the spatiotemporal computation region of the model, The control variables representing the model are (open boundary conditions, bottom friction coefficient, interlayer friction coefficient, or horizontal eddy viscosity coefficient). and They represent the first Simulated and observed values ​​of water level at the upper interface. The weight coefficient is set to 1 at the water point and 0 elsewhere. The gradient of the cost function with respect to the control parameters is calculated using the gradient definition method, i.e.:

[0076] ,

[0077] in, This represents the disturbance value of the control parameters.

[0078] Compared with the prior art, the present invention has at least the following beneficial effects:

[0079] This invention, by combining the "independent point dimensionality reduction scheme" with the "SPG optimization algorithm," successfully overcomes the bottlenecks of traditional methods in terms of characterization capability, computational efficiency, and practical application, achieving the following significant beneficial effects:

[0080] (1) Significantly improves the accuracy of internal tidal simulation and realizes the true reproduction of energy dissipation process.

[0081] Existing methods assume the eddy viscosity coefficient to be constant or have a simple spatial distribution, failing to reflect the complex heterogeneity of this parameter in the real ocean. This leads to distortion in the simulated internal tide energy dissipation process, which is a core bottleneck restricting simulation accuracy. This invention, for the first time within a computationally acceptable range, achieves high-precision inversion of the complex spatial distribution of the eddy viscosity coefficient. This directly brings about a qualitative leap in the accuracy of internal tide simulation.

[0082] (2) Significantly improves the efficiency of inversion calculation, making the inversion of complex distributions feasible for practical applications.

[0083] Traditional gradient definition methods, when directly applied to each grid point, result in computational costs that increase exponentially with the number of grid points, making the inversion task computationally infeasible. This invention achieves an order-of-magnitude improvement in computational efficiency through a dual approach of "dimensionality reduction" and "efficient convergence."

[0084] (3) Enhancing the robustness and applicability of the method in real complex sea areas

[0085] Existing methods are mostly limited to ideal terrain and are prone to failure in real, complex terrain due to instability in the optimization process, or may reproduce physically unreasonable coefficient distributions. This invention exhibits excellent robustness and can be stably applied to the complex environment of the real ocean. Attached Figure Description

[0086] Figure 1 This is a top-down view of the experimental area. White represents land, red dots indicate the locations of T / P altimeter data nadir observation points in waters deeper than 200m, green circles indicate independent point locations, and the black solid line represents the 500m isobath.

[0087] Figure 2 The spatial distribution diagrams of the given horizontal eddy viscosity coefficient are shown; (a) is spatial distribution 1, and (b) is spatial distribution 2.

[0088] Figure 3 The cost function changes during the assimilation process of the two inversion examples.

[0089] Figure 4 The inversion results are shown in Figure 1 for the spatial distribution of the horizontal eddy viscosity coefficient; (a) is the inversion result for spatial distribution 1, and (b) is the inversion result for spatial distribution 2.

[0090] Figure 5 This is a flowchart illustrating the inversion process of the spatial distribution of the horizontal eddy viscosity coefficient. Detailed Implementation

[0091] To make the objectives, technical solutions, and advantages of this invention clearer, the following detailed description, in conjunction with specific embodiments and accompanying drawings, further illustrates the invention. Obviously, the described embodiments are only a portion, not all, of the embodiments disclosed in this invention. All other embodiments obtained by those skilled in the art based on the embodiments disclosed in this invention without inventive effort are within the scope of protection of this invention.

[0092] Although preliminary progress has been made in retrieving eddy viscosity coefficients based on methods such as adjoint assimilation, these existing studies have significant limitations in dealing with the complexity of real marine environments, mainly reflected in their severely insufficient ability to describe the complex spatial distribution of eddy viscosity coefficients.

[0093] Example 1

[0094] The eddy viscosity coefficient inversion method provided in this embodiment, such as Figure 5 As shown, it consists of the following steps:

[0095] S1: Construct a forward numerical model of the ocean incorporating internal tidal dynamics and acquire observational data to be assimilated, such as satellite altimeter data and anchored current velocities. The forward model is used to simulate ocean conditions and parameterize the eddy viscosity coefficient. The observational data provides constraints for the inversion, ensuring that the inversion results conform to actual ocean phenomena. Specifically, the forward numerical model of the ocean is constructed as follows: This model is a nonlinear primitive equation model, using isodense coordinates in the vertical direction and spherical coordinates in the horizontal direction, considering the hydrostatic approximation and the free sea surface. The model employs internal and external mode separation techniques. The horizontal gradient term of sea surface atmospheric pressure is ignored, and it is assumed that the layered structure is horizontally uniform, i.e., the initial thickness of each layer except the bottom layer is uniform. The horizontal gradients of all layers are zero. Let the density of each layer be... Dynamic thickness is Therefore, the dynamic mass is Taking tidal force as and Thus, the internal modal governing equations are obtained:

[0096] Layer 1 (Surface):

[0097] ;

[0098] ;

[0099] ;

[0100] No. layer( ):

[0101] ;

[0102] ;

[0103] ;

[0104] No. Layer (bottom layer):

[0105] ;

[0106] ;

[0107] ;

[0108] in, For time, and These are geographical longitude and latitude, respectively, and are unit vectors. , and East, north, and upward are respectively considered positive. and They are respectively and The velocity component in the direction; The radius of the Earth; For Coriolis parameters, where This represents the magnitude of the Earth's angular velocity of rotation. It is the acceleration due to gravity; The horizontal eddy viscosity coefficient; the Laplace operator is specifically expressed as:

[0109] ,

[0110] .

[0111] The vertical eddy viscosity term can be written in the form of interlayer friction:

[0112] , ,

[0113] in, For the first Layer and first The average density of the layer, For the first Layer and first The average initial thickness of the layer, is the interlayer friction coefficient. The bottom friction term in the momentum equation for the bottom layer (bottom layer) takes the following form:

[0114] ,

[0115] in The coefficient of friction is the lowest.

[0116] The external modal equations are obtained by vertically integrating the internal modal equations; firstly, the barotropic flow is defined as the mass-weighted average of the velocities in each layer, i.e.

[0117] , , ;

[0118] The momentum equations of each layer of the internal mode are... By weighted averaging and directly adding the continuity equations of each layer, the external modal equations can be obtained:

[0119] ;

[0120] ;

[0121] ;

[0122] in:

[0123] ,

[0124] ,

[0125] ,

[0126] ,

[0127] For the general still water depth, The vertical average density of the water column. and The first term on the right is external pressure, caused by the undulations of the sea surface; the second term is internal pressure, caused by the horizontal density gradient of the fluid in each layer. and The first term on the right is the mass-weighted average of the convection terms for each layer, the second term is the mass-weighted average of the turbulent viscous forces for each layer, and the third term is the difference between the local acceleration of the barotropic flow and the mass-weighted average of the local acceleration of the velocity for each layer.

[0128] S2: Parameterization and Dimensionality Reduction Based on Independent Point Scheme: When defining the spatial distribution of the eddy viscosity coefficient across the entire computational domain, instead of directly inverting the value of each grid point, a set of independent control points—far fewer than the total number of model grids—is arranged within the computational domain. The coefficient values ​​at these independent points are the core parameters to be inverted. Using a pre-defined spatial interpolation method (such as bilinear interpolation, Cressman interpolation, etc.), the continuous and complex spatial distribution of the entire field is reconstructed from the coefficient values ​​of the independent points. A selection of grid points (denoted as...) are then chosen within the computational domain. As an independent point, let For the control variable values ​​at independent points, calculate the values ​​at each grid point in the calculation region (denoted as...). The value of ) Depend on Obtained through linear interpolation, and The following relationship must be satisfied:

[0129] ,

[0130] in, These are weighting coefficients, for example, they can be expressed in the Cressman form, i.e.:

[0131] ,

[0132] in, For grid and The center distance between them The radius of influence. The cost function can be derived. about The gradient is:

[0133] .

[0134] S3: Calculating the gradient of the objective function using the gradient definition method: This is performed in the reduced-dimensional parameter space established in S2, defining an objective function that measures the difference between the model simulation results and the observed data; subsequently, using the gradient definition method, the gradient of the objective function with respect to the coefficients of all independent points is calculated by perturbing the coefficient values ​​at each independent point and running the forward model; the specific implementation of the gradient definition method: estimating the spatially distributed horizontal eddy viscosity coefficient by assimilating satellite altimeter data. First, a cost function is constructed based on simulated and observed water levels:

[0135] ;

[0136] in This refers to the index of the isodense water layer where the internal modal variables reside. For the spatiotemporal computation region of the model, The control variables representing the model are (open boundary conditions, bottom friction coefficient, interlayer friction coefficient, or horizontal eddy viscosity coefficient). and They represent the first Simulated and observed values ​​of water level at the upper interface. The weight coefficient is set to 1 at the water point and 0 elsewhere. The gradient of the cost function with respect to the control parameters is calculated using the gradient definition method, i.e.:

[0137] ,

[0138] in, This represents the disturbance value of the control parameters.

[0139] S4: Applying the SPG algorithm for efficient iterative optimization: The gradient calculated in S3 is input into the SPG algorithm. The SPG algorithm intelligently determines the search direction and step size by utilizing the current iteration point and its gradient information, combined with historical iteration information, and applies physical constraints on the range of independent point coefficients during this process. Detailed steps of the SPG algorithm can be found in the literature (Birgin EG, Martinez JM, Raydan M (2000) Nonmonotone spectral projected gradient methods on convex sets. SIAM Journal on Optimization, 10:1196–1211. https: / / doi.org / 10.1137 / S1052623497330963). The SPG algorithm is renowned for its superlinear convergence characteristics. Compared to traditional optimization algorithms (such as the steepest descent method), it significantly reduces the number of iterations required to reach convergence. This, combined with the dimensionality reduction effect of S2, achieves a doubling of computational efficiency, ensuring that the entire inversion process obtains a high-precision solution within an acceptable timeframe.

[0140] S4 is the solver for the inversion problem. It receives the gradient provided by S3 and performs an iterative search in the reduced-dimensional parameter space defined by S2. Finally, through the coordinated execution of the above four steps, an eddy viscosity coefficient field with a complex spatial distribution that can fit the observation data with high accuracy is output.

[0141] Example 2:

[0142] Twin experiments are a way to verify the effectiveness of inversion algorithms. The specific steps are as follows: run a forward model using the given spatially distributed horizontal eddy viscosity coefficient, and use the sea surface water level value at the observation point as the simulation result as the "observation data". Then, assign an initial guess value to the horizontal eddy viscosity coefficient, optimize the horizontal eddy viscosity coefficient of independent points by assimilating the observation data, and thus invert the spatially distributed horizontal eddy viscosity coefficient.

[0143] In this example, the simulation area for the numerical experiment is the sea area near the Luzon Strait. Figure 1The geographic range is 115.5°E~127°E and 15.5°N~26°N, with 139×127 horizontal grids and a horizontal resolution of 5′×5′. Vertically, it is divided into two layers, with the initial layer interface set at a depth of 500m. Coriolis parameters are taken as local values, and the angular frequency of the M2 tidal wave is 1.4050789025E⁻⁴s⁻¹. The time step is 4.969 min (1 / 150 of the M2 tidal cycle). The harmonic constant of the M2 tidal cycle from T / P altimeter satellite data at a depth of 200m or more is used as the observation data. Since the T / P altimeter's nadir observation point is not necessarily at the center of the model's computational grid, and there may be multiple observation points within the same grid, we take the average of all observation data within the same grid as the observation value at the center of that grid. The location of the observation point obtained from this is shown below. Figure 1 As shown, there are a total of 780 observation points used for assimilation.

[0144] The spatial distribution of the horizontal eddy viscosity coefficient is optimized by assimilating observational data. We assume that the value of the horizontal eddy viscosity coefficient depends only on the horizontal position and remains unchanged in the vertical direction. Two spatial distributions of the horizontal eddy viscosity coefficient are given, both formed by parabolic surfaces of revolution, one opening downwards (distribution 1) and the other opening upwards (distribution 2). The values ​​of both distributions are between 100 and 300. The given spatial distributions are as follows: Figure 2 As shown. Accordingly, we use two inversion examples to invert these two distributions respectively. When inverting the horizontal eddy viscosity coefficient, to realize its spatial distribution, we adopt the independent point scheme, such as... Figure 1 As shown, some independent points are evenly distributed in the calculation area, with the shortest distance between independent points being 1.5°, and the influence radius in the independent point scheme is set to 2.25°.

[0145] The assimilation process involved 100 iterations, and the changes in the cost function during the assimilation process were as follows: Figure 3 As shown, after assimilation, the cost functions of both experiments decreased significantly compared to their values ​​before assimilation, by about four orders of magnitude. This indicates that the difference between the simulation results and the observation results has been greatly reduced, and therefore the model's assimilation effect on the observation data is significant.

[0146] The optimized spatial distribution of eddy viscosity coefficient is as follows: Figure 4 As shown. It can be intuitively seen that in the two examples, the inverted spatial distribution ( Figure 4 ) and a given spatial distribution ( Figure 2 The results show good consistency, with correlation coefficients all above 0.9. In both examples, the mean square error of the inverted distribution after assimilation and the given distribution decreased by about an order of magnitude compared with that before assimilation, indicating that the horizontal eddy viscosity coefficient was effectively optimized by assimilating the observation data.

[0147] In summary, the inversion of the spatially distributed horizontal eddy viscosity coefficient using this adjoint assimilation model is highly successful, and the independent point scheme is an effective means of realizing the spatial distribution of the horizontal eddy viscosity coefficient. The two twin experimental examples above verify that the eddy viscosity coefficient inversion algorithm of this invention is very effective. In practical applications, there are many choices for setting the independent point scheme for the spatial distribution of the horizontal eddy viscosity coefficient, which can be flexible and varied. For example, near closed boundaries and steep terrain, the spatial variation of the horizontal eddy viscosity force is relatively complex, and the independent points can be set more densely, with a corresponding smaller radius of influence. In waters far from the lateral boundary, the spatial variation of the horizontal eddy viscosity force is relatively gentle, and the independent points can be set relatively sparsely, with a corresponding larger radius of influence, or even the horizontal eddy viscosity coefficient can be directly treated as a constant.

[0148] In one embodiment, in an ideal experiment targeting a typical complex topographic region, the complex coefficient distribution inverted by this invention was used to simulate internal tides. The results showed a reduction in root mean square error of approximately 20% compared to internal tide amplitudes observed by satellite altimeters. Simultaneously, the spatial distribution of the simulated internal tide energy flux divergence exhibited higher consistency with the observed ocean mixing "hotspots." This effect is directly attributable to the independent point scheme. This scheme ensures the spatial continuity of the coefficient field through interpolation functions, and its flexible control points allow it to fully fit complex spatial variations driven by topography, stratification, etc., thus providing a more physically realistic description of the energy dissipation process.

[0149] In one embodiment, in a real-world terrain case with over 100,000 grid points, directly inverting the coefficients of all grid points would require an estimated 1000 CPU hours for a single gradient calculation. However, by employing the independent point scheme of this invention (setting only 500 control points), the parameters to be optimized are reduced by three orders of magnitude. Combined with the superlinear convergence characteristics of the SPG algorithm, the total inversion time is shortened to approximately 1 / 50th of the direct method, enabling inversion to be completed in days rather than years, thus possessing the capability to handle practical problems. The independent point scheme fundamentally reduces the dimensionality of the optimization problem. The SPG algorithm, by intelligently determining the search direction using historical gradient information, reduces the number of iterations required for convergence by approximately 60% compared to the traditional steepest descent method. Together, they solve the core problems of being "too difficult to compute" and "too slow to compute."

[0150] In one embodiment, the present invention was applied to the Luzon Strait, a complex sea area with the strongest endorheic energy globally. It successfully inverted the non-uniform, high-intensity patchy distribution of the eddy viscosity coefficient in this region, formed by the combined effects of strong tidal currents, steep topography, and mesoscale eddies. The inversion process was stable, without any non-convergence or numerical oscillations. The SPG algorithm itself has good numerical stability, and its built-in projection step effectively constrains the solution within a physically reasonable range (e.g., ensuring the viscosity coefficient is non-negative). Simultaneously, the ill-posedness of the optimization problem is reduced after dimensionality reduction, further improving the stability of the solution.

[0151] The specific embodiments described above further illustrate the purpose, technical solution, and beneficial effects of the present invention. It should be understood that the above descriptions are merely specific embodiments of the present invention and are not intended to limit the scope of the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. A method for inverting complex spatially distributed eddy viscosity coefficients in an oceanic internal tide model, characterized in that, Comprise the following steps: S1: Constructing a numerical forward model of ocean including internal tidal dynamic process, and obtaining observation data to be assimilated; S2: Parameterizing and dimensionality reduction based on independent point scheme: arranging a group of independent control points in the calculation domain, the number of which is less than the total number of model grid, the coefficient value of these independent points is the core parameter to be inverted; through the preset spatial interpolation method, the continuous and complex spatial distribution of the whole field is reconstructed from the coefficient value of the independent points; S3: Calculate the gradient of the objective function by gradient definition method: in the dimensionality reduction parameter space established in S2, define an objective function to measure the difference between the model simulation results and the observation data; then, by perturbing the coefficient value of each independent point and running the forward model, the gradient of the objective function with respect to all independent point coefficients is calculated by gradient definition method; S4: Apply SPG algorithm for efficient iterative optimization: input the gradient calculated in S3 into SPG algorithm; SPG algorithm uses the current iteration point and its gradient information, combined with historical iteration information, to intelligently determine the search direction and step length, and in this process, the value range of the independent point coefficient is physically constrained; Wherein, the S4 is the solver of the inversion problem, which receives the gradient provided by S3 and performs iterative search in the dimensionality reduction parameter space defined in S2; finally, through the cooperative execution of the above four steps, a complex spatial distribution of the eddy viscosity coefficient field with high precision fitting observation data is finally output.

2. The inversion method of claim 1, wherein, S1, the ocean numerical forward model, is a nonlinear primitive equation model that uses isodense coordinates in the vertical direction and spherical coordinates in the horizontal direction, considering hydrostatic approximation and free sea surface. The model employs internal and external mode separation techniques; it neglects the horizontal gradient term of sea surface atmospheric pressure and assumes that the layered structure is horizontally uniform, i.e., the initial thickness of each layer except the bottom layer is uniform. The horizontal gradients of all layers are 0, and the density of each layer is assumed to be... Dynamic thickness is Therefore, the dynamic mass is Taking tidal force as and Thus, the internal modal governing equations are obtained: The first layer is the surface layer: ； ； ； No. layer( ): ； ； ； First layer: ； ； ； where is time, and are the geographical longitude and latitude, respectively, and , and are unit vectors in the east, north and up directions, respectively; and are the velocity components in the and directions, respectively; is the Earth radius; is the Coriolis parameter, where is the magnitude of the Earth's rotation angular velocity; is the gravitational acceleration; is the horizontal eddy viscosity coefficient; the Laplace operator is explicitly given by: ， ； The vertical eddy viscosity term is written as the form of interlayer friction: ， ， wherein, is the first layer and the second average density of the first layer and the second average initial thickness of the first layer and the second is the interfacial friction coefficient, and the bottom friction term in the momentum equation for the first layer takes the form: ， wherein is the bottom friction coefficient; The external modal equation is obtained by vertically integrating the internal modal equation; first, define the barotropic flow as the mass weighted average of the flow velocity of each layer, that is: ， ， ； The momentum equations of the inner mode are weighted and averaged with the weight of 1 / 2, and the continuous equations of the layers are directly added to obtain the equation of the outer mode. ； ； ； Wherein, ， ， ， ， is the total hydrostatic depth, is the vertically averaged density of the water column, and The first term on the right is the external pressure, due to the sea surface elevation, and the second term is the internal pressure, due to the horizontal density gradient of the layers, and The first term on the right is the mass-weighted average of the advection terms in the layers, the second term is the mass-weighted average of the turbulent viscous forces in the layers, and the third term is the difference between the mass-weighted average of the local acceleration of the barotropic flow and the local acceleration of the flow velocity in the layers.

3. The inversion method of claim 1, wherein, The core role of S2 is to reduce the dimension of the inversion problem, which converts a high-dimensional problem with tens of thousands or even millions of grid parameters to be solved into a low-dimensional problem with only dozens to hundreds of independent point parameters to be optimized; assuming that the control variable is spatially distributed, and to achieve the spatial distribution of the control variable, the "independent point scheme" is adopted in the optimization process of the control variable, that is, some points in the calculation region are selected as independent points, and the spatial distribution of the control variable is obtained by linear interpolation of the values of the independent points, thereby reducing the number of control variables to ensure the well-posedness of the inverse problem, and physically maintaining the continuity and smoothness of the spatial distribution of the control variable.

4. The inversion method of claim 3, wherein, In S2, some grid points in the calculation region are selected, denoted as As an independent point, let The control variable value at the independent point is calculated, and the grid points in the calculation region are denoted as The value of the grid point is The value of the grid point is Linear interpolation is performed to obtain And The following relationship is satisfied: ， wherein, is a weight coefficient, in the form of Cressman, i.e.: ， where, is the center distance between the grid and , is the influence radius, the cost function can be derived as The gradient of is: 。 5. The inversion method of claim 3, wherein, The gradient definition method in S3 is specifically: estimating the spatial distribution of horizontal eddy viscosity coefficient by assimilating satellite altimeter data First, a cost function is constructed based on simulated and observed water levels: ； where is the index of the water layer where the internal mode variable is located, is the spatio-temporal calculation region of the model, represents the control variable of the model, and respectively represent the simulation value and the observed value of the water level on the upper interface of the layer, is the weight coefficient, which is 1 at the water point and 0 elsewhere, and the gradient of the cost function with respect to the control parameter is calculated using the gradient definition method, i.e.: ， wherein, is the perturbation value to the control parameter.

Images