Marine space-time multi-scale multi-source multi-variable variational data assimilation method

Abstract

The invention belongs to the technical field of ocean data assimilation, and relates to an ocean space-time multi-scale multi-source variable variation data assimilation method. A space-time four-dimensional multi-scale variational analysis method is adopted, natural constraint conditions are formed through a seawater state equation, static balance and dynamic height integration, no human assumption is needed, and a more objective ocean three-dimensional gridding thermohaline objective analysis product can be constructed only by relying on mutual complementation among different observations. According to the method, spatial-temporal multi-scale information in different observation data is effectively extracted, and meanwhile, the correlation among temperature, salt and flow fields is established, so that high-quality, high-resolution and long-time-sequence objective analysis of the ocean three-dimensional gridding temperature-salt flow field can be realized, and technical support is provided for developing high-precision ocean reanalysis products.

Description

Technical Field

[0001] This invention belongs to the field of marine data assimilation technology and relates to a method for marine spatiotemporal multi-scale, multi-source, and multivariate variational data assimilation. Background Technology

[0002] Today, the continuous development of satellite remote sensing technology provides a wealth of high-resolution, wide-ranging, and continuous remote sensing data on sea surface temperature (SST), salinity, and sea surface height anomalies for ocean research. However, this data only contains sea surface information and lacks information on underwater structures. In contrast, past in-situ observation data acquired by instruments such as CTD (Conductivity-Temperature-Depth), MBT (Mechanical Bathythermograph), XBT (Expendable Bathythermograph), and Argo (Array for Real-time Geostrophic Oceanography) buoys constitute a global ocean temperature and salinity profile dataset. In-situ ocean observation data can reveal underwater temperature and salinity structures, but its distribution is sparse in horizontal space and discontinuous in time. How to combine and assimilate sea surface data obtained from satellite remote sensing with in-situ observation data to complement each other and obtain a more accurate three-dimensional ocean temperature and salinity field, thereby acquiring information on the internal structural characteristics and changes of the ocean, is a cutting-edge issue that urgently needs to be addressed in the field of ocean reanalysis.

[0003] Data assimilation refers to a method that, based on mathematical theory, finds an optimal solution between the background field and actual observations, taking into account the spatiotemporal distribution of data and errors between the observation field and the background field. Currently, various data assimilation methods have been widely applied in operational reanalysis systems both domestically and internationally to conduct reanalysis of oceanographic fields and correct ocean forecasting models. However, in addition to the basic assumption that the dynamic height integral has a zero flow layer, which is commonly used in data assimilation, many assimilation systems have added other assumptions. The following are some of the artificial assumptions used in more mature assimilation systems: (1) It is assumed that there is a linear relationship between sea surface elements and underwater elements. For example, the temperature and salinity relationship between the sea surface and underwater is calculated based on the model output results. The temperature and salinity profile "pseudo-observation" data are obtained by inverting the satellite altimetry data. Then, the temperature is assimilated using the multigrid three-dimensional variational method. Then, the salinity is calculated from the temperature analysis field according to the temperature-salinity relationship to obtain the salinity background field used for assimilation. Finally, the salinity is assimilated to obtain the three-dimensional temperature and salinity analysis field. (2) It is assumed that the vertical background field error covariance matrix of the elements is defined. Data assimilation methods such as 3D-Var (Three Dimensional Variation), OI (Optimal Interpolation), and EnKF (Ensemble Kalman Filter) used in assimilation all require the background field error covariance matrix to be defined manually. (3) It is assumed that the water body rise and fall process is an adiabatic process. These assumptions undoubtedly increase human error in the assimilation process. How to establish a more objective relationship between temperature, salinity, and flow is an urgent problem to be solved to further improve the accuracy of assimilation. Summary of the Invention

[0004] To address the problems existing in the prior art, this invention provides a method for assimilating multi-scale, multi-source, and multi-variable variational data in the ocean. Without making any relational assumptions, it effectively extracts multi-scale information from different observational data while establishing correlations between temperature, salinity, and current fields. This enables high-quality, high-resolution, long-term objective analysis of the ocean's three-dimensional gridded temperature, salinity, and current fields, providing technical support for the development of high-precision ocean reanalysis products. This method is mainly applied to the reanalysis of ocean observation data, enabling the development of high-quality, high-resolution, long-term objective analysis products of ocean three-dimensional gridded temperature, salinity, and current fields.

[0005] The technical solution adopted by the present invention to achieve the above objectives is as follows:

[0006] A method for assimilating multi-scale, multi-source, and multivariate variational data from the ocean in a spatiotemporal manner, specifically including the following steps:

[0007] S1. In the observation increment calculation module, the FGAT (First Guess at Appropriate Time) method is first used to calculate the daily average temperature, salinity and sea surface height anomalies as the background field using the model output results, and then projected onto the observation points to calculate the observation increment. Finally, the observation increments of temperature, salinity and SSHa (Sea Surface Height Anomaly, SSHa) and the anomaly field of the model output results relative to the background field are obtained.

[0008] S2. In the multigrid 3D variational data assimilation module, a preliminary analysis grid background field is first generated using the daily mean field. Then, the analysis grid is sequentially refined based on the bisection method, and 3D variational data assimilation is performed on the refined grid. The analysis field obtained on the coarser grid in the previous step is used as the background field, and the analysis increment is also the increment relative to the previous analysis result. It is determined whether the maximum grid multiplicity has been reached. If not, the iterative calculation continues; if it has been reached, the final temperature and salinity analysis field is obtained.

[0009] S3. Based on the dynamic height anomaly calculation in the ocean current geostrophic adjustment module, generate an ocean current SSHa analysis field that matches the temperature and salinity analysis field, and use the assimilated temperature and salinity analysis field to perform geostrophic adjustment on the flow field of the above average model to obtain an ocean current analysis field that matches it.

[0010] S4, the post-processing module of the analysis field uses the anomaly field calculated in S1 to superimpose the tidal signal on the analysis field to obtain the final analysis field, and then enters S1 to start the analysis calculation for the next day.

[0011] Step S1 specifically includes:

[0012] Calculation of FGAT increments from temperature and salinity observations:

[0013] The hourly output results are averaged over a total of 25 hours before and after the analysis field time or the start time of the model integration, and the model tide information is filtered out. The hourly output results are as follows:

[0014]

[0015] in, The representative model outputs hourly sea surface height results. The representative mode outputs hourly temperature results. The representative mode outputs hourly salinity results. The representative mode outputs hourly results of the east-west velocity components of the sea surface. The representative model outputs hourly results of the north-south velocity components of the sea surface.

[0016] Record the differences between the analysis field time-time model results and the average results of this model. and The daily mean fields for temperature and salinity are as follows:

[0017]

[0018] The spatiotemporal points of the field observations are calculated, and spatiotemporal interpolation is performed on the model results. The observation results are then subtracted from the model spatiotemporal interpolation results to obtain the observation increment:

[0019]

[0020] Step 2 specifically includes:

[0021] S21, Multigrid 3D variational form:

[0022] Let column vector X b Representing the background field, column vector Y obs Represents the observation field, column vector X a Let R represent the analysis field, H() be the observation field error covariance matrix, H() be the observation projection operator, and S be the smoothing matrix, which is the integral of the square of the Laplace operator of the control variables over the entire space. In the multigrid three-dimensional variational data assimilation method, the objective functional should take the following form:

[0023]

[0024] Wherein, the superscript (n) represents the nth grid. P is the analysis increment for the nth grid. (n) This is the projection matrix from the nth mesh to the original mesh; coarse meshes correspond to long-wavelength modes, and fine meshes correspond to short-wavelength modes. The wavelength or correlation scale is expressed by the mesh size. The final analysis result is expressed as:

[0025]

[0026] The coarseness of the grid is used to describe the correlation scale in the background field error covariance matrix. Three-dimensional variational analysis is performed on the increment of the observed field relative to the background field on a set of grids from coarse to fine. In each analysis, the analysis field obtained from the previous analysis on the coarser grid is used as the new background field and substituted into the analysis on the next analysis on the finer grid. The increment of each analysis is also the increment relative to the new background field obtained from the previous analysis on the coarser grid. Finally, the analysis results of each grid are superimposed to obtain the final analysis result.

[0027] S22. Based on this, we extend it to a spatiotemporal four-dimensional multigrid variational analysis framework:

[0028] Within the spatiotemporal four-dimensional variational analysis framework, multigrid technology is used to perform variational analysis on spatiotemporal coarse grids to fine grids, thereby optimally combining large-scale, mesoscale to small-scale information from ocean numerical models and observational data.

[0029] The objective function of the spatiotemporal multi-scale fast data assimilation framework is as follows:

[0030]

[0031] in, and The control variables are temperature increment and salinity increment, both including four dimensions: longitude, zonation, vertical, and time. The subscripts T, S, SST, SSS, and SSHa represent temperature, salinity, sea surface temperature, sea surface salinity, and sea surface height anomaly, respectively. ΔY is the observation increment; h(·) is the sea surface height anomaly increment obtained by the dynamic height integration of the temperature and salinity analysis results. The gradient of the above objective function is:

[0032]

[0033] For the analysis grid setup, a bisection method is used to refine the analysis grid layer by layer. The number of grid points in each layer of the three spatial dimensions and one temporal dimension satisfies 2-1 with the number of layers. (n-1) The relationship is +1, where n is the grid multiplicity; the horizontal grid and the time grid are uniformly densed, and the vertical grid adopts a variable grid form with denser top and sparser bottom.

[0034] Step 3 specifically includes:

[0035] S31. Seawater state equation, static equilibrium, and dynamic height integral constraints:

[0036] Changes in sea surface height are determined by variations in temperature and salinity distribution below the sea surface. The dynamic height anomaly is calculated as follows:

[0037]

[0038] Where ρ(P, T, S) is the UNESCO 1981 seawater density equation of state, and T m S m These are the multi-year average temperature and salinity values ​​calculated from other historical temperature and salinity data; T and S are respectively the multi-year average temperature and salinity values ​​calculated from other historical temperature and salinity data. T and X S ;z K The reference depth is z; the vertical water depth is P; the pressure is P.

[0039] In numerical experiments, it is necessary to calculate the dynamic height on a discrete vertical grid. Therefore, the integral of the above equation is rewritten as a summation as follows:

[0040]

[0041] Where k is the vertical depth layer number; K is the depth layer number of the reference depth; T b S b The background temperature and salinity values ​​are denoted as Δz; the depth difference between two adjacent depth layers is Δz; the actual satellite remote sensing sea level anomaly is the abnormal value of the daily average sea level height on a certain day relative to the annual average sea level height over many years. The selection of the annual average sea level height value over many years has an impact on the analysis results.

[0042] S32. Static stability constraints:

[0043] Introduce the following function α(T,S) into the objective function:

[0044]

[0045] in:

[0046]

[0047] As a weighting coefficient, when unstable stratification occurs, i.e. when the upper layer density is greater than the lower layer density, the α(T,S) function will increase. Introducing this function into the objective function can force the optimization of the state field to proceed in the direction of stable stratification.

[0048] S33, Ocean current geostrophic adjustment:

[0049] By applying geostrophic adjustment to the averaged model flow field using the assimilated temperature-salinity analysis field, a matching ocean current analysis field is obtained. The thermal wind equation is:

[0050]

[0051] The geostrophic increment of the analytical field relative to the background field is calculated as follows:

[0052]

[0053] By superimposing the geostrophic flow increment onto the ocean current background field, we obtain the average ocean current analysis field:

[0054]

[0055] The present invention has the following beneficial effects:

[0056] 1. Fusion of multi-source observation data

[0057] This invention can effectively extract spatiotemporal multi-scale information from different observation data, establish the correlation between temperature, salinity, and flow field, and make full use of the complementary advantages of on-site observation and satellite remote sensing data in the spatiotemporal dimension, so that the two data sources complement each other in time and space, and finally obtain analysis results that match underwater and sea surface.

[0058] 2. Multivariate analysis and natural constraints

[0059] This invention differs from univariate assimilation by simultaneously incorporating temperature, salinity, and SSHa satellite remote sensing data. Through a dynamic height integration process that includes the seawater density equation of state, it establishes a more objective and natural connection between various observational data. Without requiring artificial assumptions, salinity analysis is supplemented and constrained by temperature and sea surface height data, resulting in a positive adjustment and achieving a "1+1>2" effect. Therefore, this scheme can generate mutually matched reanalysis data on temperature, salinity, current, and sea surface height anomalies, thus providing strong data support for ocean forecasting and prediction, as well as research on ocean dynamics and thermodynamics. Attached Figure Description

[0060] Figure 1 This is a schematic diagram of a spatiotemporal multi-scale rapid data assimilation framework.

[0061] Figure 2 The implementation process of spatiotemporal multi-scale, multi-source, and multivariate variational data assimilation method.

[0062] Figure 3 This is a scatter plot of the assimilation analysis results based on the present invention.

[0063] Figure 4 This is a comparison chart of the SST spatial distribution results of the assimilation analysis based on this invention.

[0064] Figure 5 This is a comparison chart of the spatial distribution of SSHa based on the assimilation analysis results performed according to this invention. Detailed Implementation

[0065] The technical solution of the present invention will be described in detail below with reference to the accompanying drawings and embodiments:

[0066] Example

[0067] like Figure 2 The method for assimilating multi-scale, multi-source, and multivariate variational data from the ocean, as shown, specifically includes the following steps:

[0068] S1. In the observation increment calculation module, the FGAT (First Guess at Appropriate Time) method is first used to calculate the daily average temperature, salinity and sea surface height anomalies as the background field using the model output results, and then projected onto the observation points to calculate the observation increment. Finally, the observation increments of temperature, salinity and SSHa (Sea Surface Height Anomaly, SSHa) and the anomaly field of the model output results relative to the background field are obtained.

[0069] Calculation of FGAT increments from temperature and salinity observations:

[0070] The hourly output results are averaged over a total of 25 hours before and after the analysis field time or the start time of the model integration, and the model tide information is filtered out. The hourly output results are as follows:

[0071]

[0072] in, The representative model outputs hourly sea surface height results. The representative mode outputs hourly temperature results. The representative mode outputs hourly salinity results. The representative mode outputs hourly results of the east-west velocity components of the sea surface. The representative model outputs hourly results of the north-south velocity components of the sea surface.

[0073] Record the differences between the analysis field time-time model results and the average results of this model. and The daily mean fields for temperature and salinity are as follows:

[0074]

[0075] The spatiotemporal points of the field observations are calculated, and spatiotemporal interpolation is performed on the model results. The observation results are then subtracted from the model spatiotemporal interpolation results to obtain the observation increment:

[0076]

[0077] S2. In the multigrid 3D variational data assimilation module, a preliminary analysis grid background field is first generated using the daily mean field. Then, the analysis grid is sequentially refined based on the bisection method, and 3D variational data assimilation is performed on the refined grid. The analysis field obtained on the coarser grid in the previous step is used as the background field, and the analysis increment is also the increment relative to the previous analysis result. It is determined whether the maximum grid multiplicity has been reached. If not, the iterative calculation continues; if it has been reached, the final temperature and salinity analysis field is obtained.

[0078] S21, Multigrid 3D variational form:

[0079] Let column vector X b Representing the background field, column vector Y obs Represents the observation field, column vector X a Let R represent the analysis field, H() be the observation field error covariance matrix, H() be the observation projection operator, and S be the smoothing matrix, which is the integral of the square of the Laplace operator of the control variables over the entire space. In the multigrid three-dimensional variational data assimilation method, the objective functional should take the following form:

[0080]

[0081] Wherein, the superscript (n) represents the nth grid. P is the analysis increment for the nth grid. (n) This is the projection matrix from the nth mesh to the original mesh; coarse meshes correspond to long-wavelength modes, and fine meshes correspond to short-wavelength modes. The wavelength or correlation scale is expressed by the mesh size. The final analysis result is expressed as:

[0082]

[0083] The coarseness of the grid is used to describe the correlation scale in the background field error covariance matrix. Three-dimensional variational analysis is performed on the increment of the observed field relative to the background field on a set of grids from coarse to fine. In each analysis, the analysis field obtained from the previous analysis on the coarser grid is used as the new background field and substituted into the analysis on the next analysis on the finer grid. The increment of each analysis is also the increment relative to the new background field obtained from the previous analysis on the coarser grid. Finally, the analysis results of each grid are superimposed to obtain the final analysis result.

[0084] S22. Based on this, we extend it to a spatiotemporal four-dimensional multigrid variational analysis framework:

[0085] Within the framework of spatiotemporal four-dimensional variational analysis, multigrid techniques are used to perform variational analysis on spatiotemporal coarse to fine grids, such as... Figure 1 Large-scale, meso-scale, and micro-scale information from ocean numerical models and observational data are successively and optimally combined.

[0086] The objective function of the spatiotemporal multi-scale fast data assimilation framework is as follows:

[0087]

[0088] in, and The control variables are temperature increment and salinity increment, both including four dimensions: longitude, zonation, vertical, and time. The subscripts T, S, SST, SSS, and SSHa represent temperature, salinity, sea surface temperature, sea surface salinity, and sea surface height anomaly, respectively. ΔY is the observation increment; h(·) is the sea surface height anomaly increment obtained by the dynamic height integration of the temperature and salinity analysis results. The gradient of the above objective function is:

[0089]

[0090] For the analysis grid setup, a bisection method is used to refine the analysis grid layer by layer. The number of grid points in each layer of the three spatial dimensions and one temporal dimension satisfies 2-1 with the number of layers. (n-1) The relationship is +1, where n is the grid multiplicity; the horizontal grid and the time grid are uniformly densed, and the vertical grid adopts a variable grid form with denser top and sparser bottom.

[0091] S3. Based on the dynamic height anomaly calculation in the ocean current geostrophic adjustment module, generate an ocean current SSHa analysis field that matches the temperature and salinity analysis field, and use the assimilated temperature and salinity analysis field to perform geostrophic adjustment on the flow field of the above average model to obtain an ocean current analysis field that matches it.

[0092] S31. Seawater state equation, static equilibrium, and dynamic height integral constraints:

[0093] Changes in sea surface height are determined by variations in temperature and salinity distribution below the sea surface. The dynamic height anomaly is calculated as follows:

[0094]

[0095] Where ρ(P, T, S) is the UNESCO 1981 seawater density equation of state, and T m S m These are the multi-year average temperature and salinity values ​​calculated from other historical temperature and salinity data; T and S are respectively the multi-year average temperature and salinity values ​​calculated from other historical temperature and salinity data. T and X S ;z K The reference depth is z; the vertical water depth is P; the pressure is P.

[0096] In numerical experiments, it is necessary to calculate the dynamic height on a discrete vertical grid. Therefore, the integral of the above equation is rewritten as a summation as follows:

[0097]

[0098] Where k is the vertical depth layer number; K is the depth layer number of the reference depth; T b S bThe background temperature and salinity values ​​are denoted as Δz; the depth difference between two adjacent depth layers is Δz; the actual satellite remote sensing sea level anomaly is the abnormal value of the daily average sea level height on a certain day relative to the annual average sea level height over many years. The selection of the annual average sea level height value over many years has an impact on the analysis results.

[0099] S32. Static stability constraints:

[0100] Introduce the following function α(T, S) into the objective function:

[0101]

[0102] in:

[0103]

[0104] As a weighting coefficient, when an unstable stratification occurs in the density, i.e. when the density of the upper layer is greater than that of the lower layer, the α(T,S) function will become larger. Introducing this function into the objective function can force the optimization of the state field to proceed in the direction of stable stratification.

[0105] S33, Ocean current geostrophic adjustment:

[0106] By applying geostrophic adjustment to the averaged model flow field using the assimilated temperature-salinity analysis field, a matching ocean current analysis field is obtained. The thermal wind equation is:

[0107]

[0108] The geostrophic increment of the analytical field relative to the background field is calculated as follows:

[0109]

[0110] By superimposing the geostrophic flow increment onto the ocean current background field, we obtain the average ocean current analysis field:

[0111]

[0112] S4, the post-processing module of the analysis field uses the anomaly field calculated in S1 to superimpose the tidal signal on the analysis field to obtain the final analysis field, and then enters S1 to start the analysis calculation for the next day.

[0113] Figure 3 The scatter plot shows the assimilation analysis results based on this invention. The closer the scatter points are to y = x, the higher the analysis accuracy.

[0114] Figure 4 The image shows a comparison of the spatial distribution of SST based on the assimilation analysis results of this invention (unit: °C). From left to right, the images represent the background field, the observation field, and the analysis results.

[0115] Figure 5 The image shows a comparison of the spatial distribution of SSHa in the assimilation analysis results based on this invention. From left to right, the images represent the background field, the observation field, and the analysis results.

Claims

1. A method for assimilating multi-scale, multi-source, and multi-variable variational data from the ocean in a spatiotemporal manner, characterized in that: Specifically, the steps include the following: S1. In the observation increment calculation module, the FGAT method is first used to calculate the daily average temperature, salinity and sea surface height anomalies as the background field using the model output results, and then projected onto the observation points to calculate the observation increments. Finally, the observation increments of temperature, salinity and SSHa and the anomaly field of the model output results relative to the background field are obtained. S2. In the multigrid 3D variational data assimilation module, a preliminary analysis grid background field is first generated using the average field of the day. Then, the analysis grid is successively refined based on the bisection method, and 3D variational data assimilation is performed on the refined grid. The analysis field obtained on the coarser grid in the previous step is used as the background field, and the analysis increment is also the increment relative to the previous analysis result. It is determined whether the maximum grid multiplicity has been reached. If it has not been reached, the iterative calculation continues. If it has been reached, the final temperature and salinity analysis field is obtained. S3. Based on the dynamic height anomaly calculation in the ocean current geostrophic adjustment module, generate an ocean current SSHa analysis field that matches the temperature and salinity analysis field, and use the assimilated temperature and salinity analysis field to perform geostrophic adjustment on the flow field of the above average model to obtain an ocean current analysis field that matches it. S4, the post-processing module of the analysis field uses the anomaly field calculated in S1 to superimpose the tidal signal on the analysis field to obtain the final analysis field, and then enters S1 to start the analysis calculation for the next day.

2. The ocean spatiotemporal multi-scale multi-source multivariate variational data assimilation method according to claim 1, characterized in that, Step S1 specifically includes: Calculation of FGAT increments from temperature and salinity observations: The hourly output results are averaged over a total of 25 hours before and after the analysis field time or the start time of the model integration, and the model tide information is filtered out. The hourly output results are as follows: in, The representative model outputs hourly sea surface height results. The representative mode outputs hourly temperature results. The representative mode outputs hourly salinity results. The representative mode outputs hourly results of the east-west velocity components of the sea surface. The representative model outputs hourly results of the north-south velocity components of the sea surface. Record the differences between the analysis field time-time model results and the average results of this model. and The daily mean fields for temperature and salinity are as follows: The spatiotemporal points of the field observations are calculated, and spatiotemporal interpolation is performed on the model results. The observation results are then subtracted from the model spatiotemporal interpolation results to obtain the observation increment:

3. The ocean spatiotemporal multi-scale multi-source multivariate variational data assimilation method according to claim 1, characterized in that, Step 2 specifically includes: S21, Multigrid 3D variational form: Let column vector x b Representing the background field, column vector Y obs Represents the observation field, column vector X a Let R represent the analysis field, H() be the observation field error covariance matrix, H() be the observation projection operator, and S be the smoothing matrix, which is the integral of the square of the Laplace operator of the control variables over the entire space. In the multigrid three-dimensional variational data assimilation method, the objective functional should take the following form: Wherein, the superscript (n) represents the nth grid. P is the analysis increment for the nth grid. (n) This is the projection matrix from the nth mesh to the original mesh; coarse meshes correspond to long-wavelength modes, and fine meshes correspond to short-wavelength modes. The wavelength or correlation scale is expressed by the mesh size. The final analysis result is expressed as: The coarseness of the grid is used to describe the correlation scale in the background field error covariance matrix. Three-dimensional variational analysis is performed on the increment of the observed field relative to the background field on a set of grids from coarse to fine. In each analysis, the analysis field obtained from the previous analysis on the coarser grid is used as the new background field and substituted into the analysis on the next analysis on the finer grid. The increment of each analysis is also the increment relative to the new background field obtained from the previous analysis on the coarser grid. Finally, the analysis results of each grid are superimposed to obtain the final analysis result. S22. Based on this, we extend it to a spatiotemporal four-dimensional multigrid variational analysis framework: Within the spatiotemporal four-dimensional variational analysis framework, multigrid technology is used to perform variational analysis on spatiotemporal coarse grids to fine grids, thereby optimally combining large-scale, mesoscale to small-scale information from ocean numerical models and observational data. The objective function of the spatiotemporal multi-scale fast data assimilation framework is as follows: in, and The control variables are temperature increment and salinity increment, both including four dimensions: longitude, zonation, vertical, and time. The subscripts T, S, SST, SSS, and SSHa represent temperature, salinity, sea surface temperature, sea surface salinity, and sea surface height anomaly, respectively. ΔY is the observation increment; h(·) is the sea surface height anomaly increment obtained by the dynamic height integration of the temperature and salinity analysis results. The gradient of the above objective function is: For the analysis grid setup, a bisection method is used to refine the analysis grid layer by layer. The number of grid points in each layer of the three spatial dimensions and one temporal dimension satisfies 2-1 with the number of layers. (n-1) The relationship is +1, where n is the grid multiplicity; the horizontal grid and the time grid are uniformly densed, and the vertical grid adopts a variable grid form with denser top and sparser bottom.

4. The ocean spatiotemporal multi-scale multi-source multivariate variational data assimilation method according to claim 1, characterized in that, Step 3 specifically includes: S31. Seawater state equation, static equilibrium, and dynamic height integral constraints: Changes in sea surface height are determined by variations in temperature and salinity distribution below the sea surface. The dynamic height anomaly is calculated as follows: Where ρ(P, T, S) is the UNESCO 1981 seawater density equation of state, and T m S m These are the multi-year average temperature and salinity values ​​calculated from other historical temperature and salinity data; T and S are respectively the multi-year average temperature and salinity values ​​calculated from other historical temperature and salinity data. T and X S ;z K The reference depth is z; the vertical water depth is P; the pressure is P. In numerical experiments, it is necessary to calculate the dynamic height on a discrete vertical grid. Therefore, the integral of the above equation is rewritten as a summation as follows: Where k is the vertical depth layer number; K is the depth layer number of the reference depth; T b S b The background temperature and salinity values ​​are denoted as Δz; the depth difference between two adjacent depth layers is Δz; the actual satellite remote sensing sea level anomaly is the abnormal value of the daily average sea level height on a certain day relative to the annual average sea level height over many years. The selection of the annual average sea level height value over many years has an impact on the analysis results. S32. Static stability constraints: Introduce the following function α(T,S) into the objective function: in: As a weighting coefficient, when unstable stratification occurs, i.e. when the upper layer density is greater than the lower layer density, the α(T,S) function will increase. Introducing this function into the objective function can force the optimization of the state field to proceed in the direction of stable stratification. S33, Ocean current geostrophic adjustment: By applying geostrophic adjustment to the averaged model flow field using the assimilated temperature-salinity analysis field, a matching ocean current analysis field is obtained, and the thermal wind equation is: The geostrophic increment of the analytical field relative to the background field is calculated as follows: By superimposing the geostrophic flow increment onto the ocean current background field, we obtain the ocean current average analysis field:

Images