A method and system for gridding marine observation data based on a climate model

Abstract

The application discloses a kind of ocean observation data gridding method and system based on climate model.The application provides more optimal initial field and flow-dependent ocean parameter space by introducing coupled model comparison plan model set as dynamic set sample, and is more in line with real ocean physical variability, and then more reasonably reconstructs irregular ocean subsurface observation data into data set with complete global ocean coverage.In addition, multi-scale cyclic reconstruction technology is used, while gradually reconstructing multiple typical scale variability, the problem that traditional grid assimilation method can introduce deviation to climate state system is avoided.This method promotes the construction of more optimal long-time complete reliable ocean subsurface grid data set, and provides more complete and accurate data set for climate change and ocean science research.The application is mainly used for but not limited to ocean subsurface observation data gridding and construction of global complete grid data set.

Description

TECHNICAL FIELD

[0001] The application belongs to the field of ocean observation, data processing and gridding, and particularly relates to a method and system for gridding ocean observation data based on a climate model. BACKGROUND

[0002] Long-term complete and reliable ocean temperature and salinity gridded data is the basis for climate change and ocean science research. Limited by observation means and observation cost, global ocean in-situ observation data is sparse and uneven in space-time distribution and some element observations are scarce, which restricts the in-depth study of ocean temperature and salinity changes and their mechanisms.

[0003] In order to calculate global physical parameters (e.g. temperature, salinity, heat content) using in-situ observation data which cannot cover the global ocean, internationally, a certain "Mapping method" (gridding method) is generally used: that is, according to a certain principle, the variables in the unobserved area are "guessed". There are currently five mainstream Mapping methods internationally, but different methods have different advantages and disadvantages: (1) "Simple average" method: using the assumption that "the change in the observed area is consistent with the change in the unobserved area", which is not entirely correct, has greater uncertainty; (2) "Objective analysis" method: this method depends on the selection of prior field, the construction of background covariance (transferring information from the observed area to the unobserved area), and the general background field is the pre-defined static climatological background covariance and the initial guess field (both based on World Ocean Atlas 1994), which has a systematic bias towards the climatological mean in the areas with sparse observations; (3) "Dimension reduction interpolation" method: using short-term variability of the ocean to represent long-term variability, ignoring the characteristics of ocean variability changing over time; (4) "Sea level anomaly (SLA) reconstruction" method: using the linear correlation between sea level anomaly and subsurface data, projecting the sea surface height field to the subsurface for subsurface reconstruction. This method can only be used in years with satellite observations (after 1993); in addition, the correlation between sea surface height and upper layer variables is small in many areas: for example, the correlation between sea surface height and temperature is low in high latitudes, and if there are no in-situ temperature profile observations in these areas, relying only on satellite altimeter and sea surface temperature data will bring bias to the estimation of ocean variables; (5) "Regional calculation" method: dividing the global into two regions: ship sampling area (areas that can be covered by traditional ship observations) and Argo supplemented area (areas that cannot be covered by traditional ship observations but can be better sampled by Argo system), ignoring the interannual and interdecadal variability in the Argo supplemented area, so it can only be used for long-term trend evaluation. (6) Machine learning method: using machine learning methods to fuse satellite observation data at sea surface and profile observation data in the ocean. The defect of this method is that the amount of in-situ observation data is small, which restricts the accuracy of the reconstruction of the trained machine learning model. SUMMARY

[0004] The purpose of the present application is to overcome the shortcomings of the prior art, and to provide a method and system for gridding ocean observation data based on a climate model. The present application reconstructs ocean subsurface observation data into a global grid data set. By introducing a coupled model comparison plan as a dynamic set of samples, the present application provides a more optimal initial field and a flow-dependent ocean parameter space, and a more realistic ocean physical variability. Furthermore, the present application gradually reconstructs multiple typical scale variabilities while avoiding the problem of introducing a bias towards a climate state system deviation in traditional grid assimilation methods. The present application promotes the construction of a more optimal, long-term, complete and reliable ocean subsurface grid data set, and provides a more complete and accurate data set for climate change and ocean science research. The present application is mainly used for, but not limited to, ocean subsurface observation data gridding and the construction of a global complete grid data set. The present application specifically includes the following contents:

[0005] S1: Read the ocean data of the tth month with a resolution of Lat°xLon°, divide the global into a number of grids (number of grid points in the longitude direction*number of grid points in the latitude direction), Lat° is the resolution in the latitude direction, the number of grid points in the latitude direction is 180° / Lat°, Lon° is the resolution in the longitude direction, and the number of grid points in the latitude direction is 360° / Lat°; as shown in the following table: Figure 1

[0006] S1-1: Read the ocean state observation data (such as ocean heat content data, dissolved oxygen data, salinity data, etc.) of the tth month: first, divide the ocean depth into layers, and define a time window vector Window1=[1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,4,4,8,8,8,8,10,10,12,12,12,12,12,12,18] for different depths, with a size of 41, corresponding to depths from the 1st layer to the 27th layer, a time window of 1 month, and depths from the 28th layer to the 41st layer, a time window of 4, 4, 8, 8, 8, 8, 10, 10, 12, 12, 12, 12, 12, 12, 18 months, respectively. Then, read all the observation data in the specified time window, and the specified time window is the time window containing the observation data of the tth month; the time window of different depth layers is different, and the observation time window of the deepest 41st layer is 18 months.

[0007] S1-2: Read the climate mean state observation error data of the tth month: read the climate mean state observation error data as the uncertainty data of the observation, denoted as ​d2 is the number of depth in the climate mean state observation error data, d3 is the number of latitude grid points, d4 is the number of longitude grid points, Lat°x d3 = 180°, Lon°x d4 = 360°.

[0008] S1-3: Read the CMIP (Coupled Model Intercomparison Project) mode data of the tth month as an independent reliable data source, denoted as wherein d1 is the number of CMIP modes, d2 is the number of depths in each mode;

[0009] S1-4: Read the ocean and land boundary topography data of the tth month: denoted as S1-4: Read the ocean and land boundary topography data of the tth month: denoted as

[0010] S1-5: Read the global observation grid arithmetic mean data of the tth month: denoted as S1-5: Read the global observation grid arithmetic mean data of the tth month: denoted as

[0011] S2: Input data preprocessing (for details of the specific preprocessing process, see Figure 2 as shown) :

[0012] S2-1: Climate mean state observation error data preprocessing: First, update the climate mean state observation error data of the tth month. For the grid points with Nan (Nan represents no valid data) elements in the climate mean state observation error data, update the data of the grid points with Nan elements with the median of the non-Nan elements of the adjacent depths and adjacent grid points, denoted as

[0013] S2-2: Climate mean state observation error data preprocessing: Secondly, for the specified grid point (x, y), if L (x,y) ≥1 and then use the convolution kernel W1 i×j to do convolution on and update

[0014]

[0015] wherein E is a matrix with all elements being 1, x, y and z are longitude, latitude and depth grid indexes respectively, 0≤x

[0016] S2-3: CMIP mode data preprocessing: The data of each mode and each depth of CMIP is preprocessed one by one. For the simulation data of a certain depth if the corresponding topography data L (p,q) ==0, then If the corresponding terrain data L (p,q) ≥ 1 and Then use the convolution kernel W2 i×j Convolution is performed on

[0017]

[0018] Where E is a matrix with all elements being 1, p and q are grid indexes, i and j are convolution kernel sizes, i = 9, j = 11, count is the number of non-Nan elements in the convolution window, 0 ≤ r < d2, 0 ≤ p < d3, 0 ≤ q < d4, C t After calculation, it is

[0019] S2-4: For different depths, the observation data at time t is calculated according to the observation data of the specified time window in a certain way (which can be but not limited to the algebraic mean, weighted mean, geometric mean, etc. of the observation data of the specified time window), denoted as Where d2 is the total number of depths in the observation data, d3 is the number of latitude grid points, d4 is the number of longitude grid points, Lat° × d3 = 180°, Lon° × d4 = 360°.

[0020] S3: According to the depth, the data is processed in grid points (assuming the depth is depth):

[0021] S3-1: If depth ≥ d2, go to S4 and end the grid assimilation process; otherwise, go to S3-2 to continue the data assimilation calculation, the specific steps are:

[0022] S3-2: Calculate the CMIP ensemble average data C t,d of the tth observation month and the specified depth depth (abbreviated as d) participating in the grid calculation: t,d Get the ensemble average of all CMIP modes at the specified depth as the background average field, denoted as C t (abbreviated as C), which has:

[0023]

[0024] Where represents the matrix composed of elements with absolute value less than 10 in the CMIP ensemble average data C t of the tth observation month participating in the grid calculation, C t,d is reordered as

[0025] S3-3: Get the initial observation data (abbreviated as O), the initial climate mean state observation error data (abbreviated as D);

[0026] S3-4: Calculate initial background disturbance field A t,d (abbreviated as A):

[0027] A = C - repmat(mean(C T ), 1, m) T

[0028] wherein m is the number of CMIP models (i.e. m = d1); repmat(X, r1,..., rN) represents a specified scalar list r1,..., rN, which is used to describe how the copies of X are arranged in each dimension; mean(X) returns a row vector containing the mean of each column (i.e. the ensemble mean field of the background field).

[0029] S3-5: Corresponding depth gridding data preparation (see the specific process shown in Figure 4 ): Calculate the grid point coordinate matrix P t,d (abbreviated as P), the observation field Q t,d (abbreviated as Q), the analysis mean field U t,d (abbreviated as U), and the analysis disturbance field V t,d (abbreviated as V):

[0030] For grid point (x, y), if |O (x,y) | < 6 and L (x,y) ≥ 1, it is judged that the grid point needs to be calculated for data gridding, then S3-6 is entered, otherwise S3-7 is entered.

[0031] S3-6: Record the coordinate matrix of the grid points that need to be calculated:

[0032]

[0033] wherein P i is the number of rows of matrix P, record the analysis field mean data:

[0034]

[0035] wherein U i is the number of rows of U, k = x*d4+y, record the ensemble observation field data:

[0036]

[0037] wherein Q i is the number of rows of Q, at this time if the observation error data D (x,y) > 0, then

[0038] V (x,y) = D​(x,y) ,

[0039] Otherwise,

[0040]

[0041] where V i is the number of rows of the analysis perturbation field V, M i is the number of rows of the grid arithmetic mean data M read by S1-5, x, y are grid indices, 0≤x

[0042] S3-7: Calculate the observation perturbation optimization matrix Z t,d (referred to as Z):

[0043] Z0=Q T -mean(V T ) T

[0044] Z1=U T

[0045] Z (:,0) =Z0

[0046] Z (:,1) =Z1

[0047] S3-8: Define the scale vector T of the multi-scale cyclic reconstruction: determine the multi-scale vector of size size for the multi-scale cyclic reconstruction according to different typical spatial scales of the ocean, types of the processed ocean observation data, actual annual scale, and other factors, and denote it as T. Each scale in the scale vector T represents a localization radius, and the size size of the vector T can be determined according to specific conditions.

[0048] S3-9: Perform multi-scale cyclic reconstruction on a scale-by-scale basis (see FIG. 3 for the specific process): assuming that the scale index is jndex, if jndex exceeds the last scale in the scale vector, then depth = depth + 1 and go to S3-1, otherwise continue to S3-10 for data gridding. Figure 3

[0049] S3-10: Calculate the scaled set observation anomaly H t,d (referred to as H) of the tth observation month and the specified depth d (referred to as d) and the scaled gain vector N t,d (referred to as N) of the tth observation month and the specified depth d, and set the localization radius radius:

[0050]

[0051] N=R×Z0

[0052] ​H = R * (O - repmat(mean(V T ) T , 1, m))

[0053] where radius = T jndex diag(X) where X is a column vector, means to construct a diagonal matrix with the elements of X as the diagonal elements.

[0054] S3-11: assimilate each grid point, the single grid point data assimilation process is shown in Figure 5 : assuming that the grid index is index, if index exceeds the last grid point (i.e. index ≥ (d3*d4)), then jndex = jndex+1 goes to S3-9, otherwise continue S3-12 to perform data gridding.

[0055] S3-12: calculate the distance matrix F t,d,i (referred to as F):

[0056] if nx ≤ 1,

[0057] F = abs(P-index*I)

[0058] otherwise

[0059]

[0060] where nx is the maximum number of grid points in the longitude direction, index is the grid index, X = index%nx*I-P (:,0) , Y = (index / nx)*I-P (:,1) ; P (:,0) , P (:,1) are respectively the matrices composed of all the elements of the first column and the second column of matrix P;

[0061] S3-13: calculate the localization radius matrix J t,d,i (referred to as J):

[0062] J = exp(-0.5*F'*F')

[0063] where F' = F / radius.

[0064] S3-14: find the local observation matrix K t,d,i (referred to as K):

[0065] calculate the elements greater than 0 in J to obtain the local observation matrix K,

[0066] K = find(J>0)

[0067] S3-15: update the localization radius matrix J t,d,i(abbreviated as J):

[0068] J = J (K,:)

[0069] Take the value of the element in K as the row index, take out all the elements of all columns in all the specified rows in J, and update the localized radius matrix J;

[0070] S3-16: Update the background disturbance field data A t,d,i (A (index,:) is used in the following calculation formula to represent A t ): ,d,i ):

[0071] S = H (K,:) .*repmat(J,1,m)

[0072] G = (I + S T × S) -1

[0073] A (index,:) = A (index,:) × sqrt(G)

[0074] Wherein, if K is a natural number value, H (K,:) represents a row matrix composed of all column elements of the Kth row in matrix H; if K is a matrix, H (K,:) represents a matrix composed of all columns of the row specified by the row index in matrix H, and.* represents the multiplication of corresponding elements of the matrix one by one;

[0075] S3-17: Calculate the background average field data B t,d,i (B (index,:) is used in the following calculation formula to represent B t,d,i ):

[0076] B (index,:) = C (index,:) × (G × S T ) × (N (K,:) .* J)

[0077] Wherein, if K is a natural number value, N (K,:) represents a row matrix composed of all column elements of the Kth row in matrix N; if K is a matrix, N (K,:) represents a matrix composed of all columns of the row specified by the row index in matrix N, and.* represents the multiplication of corresponding elements of the matrix one by one;

[0078] S3-18: Update the background average field data B t,d,i (abbreviated as B):

[0079] B = B + mean(OT ) T

[0080] S3-19: Average data field X of the integration specified depth t,d and the perturbation data field Y t,d ; for the convenience of writing, the following two calculation formulas use to represent B t,d,i , to represent A t,d,i ;

[0081]

[0082]

[0083] wherein index is the grid point index, lat=index / d4 is the longitude grid index, lon=index%d4 is the latitude grid index, model is the CMIP model index, 0model t,d,i ; when index is cyclically increased, the complete Y t,d , X t,d ;

[0084] S3-20: After the grid point index assimilation is completed, index=index+1, and S3-11 is entered.

[0085] S4: Output grid point data

[0086] S4-1: Output the set average data field with complete spatial coverage after data grid point is completed (the final analysis field).

[0087] S4-2: Output the set perturbation data field with complete spatial coverage after data grid point is completed (the set sample is used to measure the uncertainty of analysis).

[0088] S5: Grid point is completed.

[0089] The advantages of the present application are as follows:

[0090] 1. Overcome the defects of the prior art:

[0091] 1) Avoid using the shorter period of ocean variability to represent the long-term variability (the disadvantage of the dimensionality reconstruction method).

[0092] 2) Avoid using the simple average method to represent the heat content of the unobserved area (the disadvantage of the simple average method).

[0093] 3) Avoid the simple use of linear correlation between independent information sources and heat content (the disadvantage of sea level anomaly reconstruction).

[0094] 4) Avoid the simple use of climatological values as initial guess field and background covariance (the disadvantage of traditional oceanic objective analysis method).

[0095] 2. Since the CMIP model can simulate typical ocean variability and long-term changes, it provides a reliable and independent source of information for data estimation, providing error covariance with flow-dependent characteristics. At the same time, the ensemble analysis of multiple models can reduce the influence of the model bias of a single model, and further better construct the background covariance, that is, it can better describe the correlation between different regions of the ocean.

[0096] 3. The ensemble square root filter method (EnSRF) is applied to the field of ocean observation, considering the estimation of multiple state variables and parameters, which helps to avoid the "climate drift" bias in a single model.

[0097] 4. Use the multi-scale cyclic reconstruction method to gradually reconstruct the variability of 3 typical scales. BRIEF DESCRIPTION OF DRAWINGS

[0098] Figure 1 Schematic diagram for input data for gridding method;

[0099] Figure 2 Schematic diagram of the pre-processing flow for input data for gridding method;

[0100] Figure 3 Schematic diagram of the multi-scale cyclic reconstruction flow;

[0101] Figure 4 Schematic diagram of the preparation flow for corresponding depth gridded data;

[0102] Figure 5 Schematic diagram of the single grid point data assimilation flow;

[0103] Figure 6A Visualization display diagram of the average ocean 600m depth observation data from January to March 1971;

[0104] Figure 6B Visualization display diagram of the reconstruction result of the gridding technology of the United States National Oceanic and Atmospheric Administration (NOAA) for the average ocean 600m depth observation data from January to March 1971;

[0105] Figure 6C Visualization display diagram of the reconstruction of the average ocean 600m depth from January to March 1971 by the method. DETAILED DESCRIPTION

[0106] The application will be described in further detail below with reference to the drawings, which are provided for explanation only and are not intended to limit the scope of the application.

[0107] To make the objects, contents and advantages of the present application clearer, the specific embodiments of the present application will be described in further detail below with reference to the drawings and examples, taking the ocean heat content data with a resolution of 1°x1° in January 1971 as an example. The method includes but is not limited to the following examples.

[0108] S1: reading input data:

[0109] S1-1: reading grid arithmetic mean observation data: first, for different depths, define a time window vector Window1=

[0110] [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 4, 4, 8, 8, 8, 8, 10, 10, 12, 12, 12, 12, 12, 12, 18], the specific size is 41, which respectively corresponds to the time window of 1 month for the 1st to 27th layers, and the time window of 4, 4, 8, 8, 8, 8, 10, 10, 12, 12, 12, 12, 12, 12, 18 months for the 28th to 41st layers. Then, read the observation data in the specified time window.

[0111] S1-2: reading climatic mean state observation error data: read the climatic mean state observation error data as the uncertainty data of observation, denoted as The number of depths in the climatic mean state observation error data is 41;

[0112] S1-3: reading Coupled Model Intercomparison Project (CMIP) model data: denoted as S1-3 is denoted as C 40×41×180×360 , wherein the number of CMIP models is 40, and the number of depths in each model is 41;

[0113] S1-4: reading land-sea boundary topography field: denoted as S1-4 is denoted as L 180×360 ;

[0114] S1-5: reading global observation statistical median data: denoted as S1-5 is denoted as M 1×27 ;

[0115] S2: input data preprocessing:

[0116] S2-1: Climate mean state observation error data preprocessing: First, update the observation error data of the observation month. For the grid point with Nan element, update the data of the corresponding grid point with the median of the non-Nan elements of the adjacent depth and adjacent grid points, denoted as

[0117] S2-2: Climate mean state observation error data preprocessing: Second, for a specified grid point (x, y), if L (x,y) ≥1 and then use the convolution kernel W1 i×j to convolve and update where x, y, z are the precision, dimension and depth grid index respectively, 0≤x<180, 0≤y<360, 0≤z<41, i=5, j=7 are the convolution kernel size.

[0118] S2-3: CMIP model data preprocessing: Preprocess the Historical-Run model data of each mode and each depth of CMIP one by one. For the simulation data of a certain depth if the corresponding terrain data L (p,q) ==0, then if the corresponding terrain data L (p,q) ≥1 and then use the convolution kernel W i×j to convolve :

[0119]

[0120] where E is a matrix with all elements being 1, p, q are grid indices respectively, i, j are convolution kernel sizes respectively, i=9, j=11, 0≤p<180, 0≤q<360, C t is calculated as

[0121] S2-4: For different depths, take the average value of the observation data in the specified time window as the observation data at time t, denoted as

[0122] S3: Perform grid point assimilation processing on the data layer by layer according to the depth (assuming the depth is depth, denoted as d):

[0123] S3-1: If depth≥41, then jump to S4 to end the grid point assimilation process; otherwise, go to S3-2 to continue the data assimilation calculation, the specific steps are:

[0124] S3-2: Calculate the model data C t,0: Obtain the simulation data of all CMIP modes at the specified depth as the initial guess data field, denoted as C t,0 (abbreviated as C), has:

[0125]

[0126] where C t,0 is obtained by rearranging

[0127] S3-3: Calculate the initial observation data (abbreviated as O 180×360 ), the initial uncertainty data (abbreviated as D 180×360 );

[0128] S3-4: Calculate the initial analysis field (abbreviated as A 64800×40 ): A = C-repmat(mean(C T ) T , 1, m)

[0129] where m is the number of CMIP modes, m = 40 when depth = 0; repmat(X, r1,..., rN) indicates that a scalar list r1,..., rN is specified to describe how the copies of X are arranged in each dimension; mean(X) returns a row vector containing the mean of each column.

[0130] S3-5: Calculate the grid point coordinate matrix P t,0 (abbreviated as P), the observation field Q t,0 (abbreviated as Q), the analysis field U t,0 (abbreviated as U), and the perturbation field V t,0 (abbreviated as V):

[0131] For grid point (x, y), if |O (x,y) | < 6 and L (x,y) ≥ 1, it is determined that the grid point needs to be calculated for data assimilation, and S3-6 is entered, otherwise S3-7 is entered;

[0132] S3-6: Record the coordinate matrix of the grid points that need to be calculated for assimilation:

[0133]

[0134] where P i is the number of rows of P. Record the analysis field:

[0135]

[0136] where Ui is the number of rows of U, k = x * 360 + y is the grid index, record set observation field:

[0137]

[0138] where Q i is the number of columns of Q, at this time if the disturbance field data D (x,y) > 0, then

[0139] V (x,y) = D (x,y) ,

[0140] otherwise,

[0141] V (x,y) = M 27 ,

[0142] where V i is the number of rows of V, M i is the number of rows of M, when depth = 0, after all grid calculations are completed, P 11023×3 , U 11023×40 , V 1×11023 , Q 1×11023 , M i = 27, x, y are grid indexes respectively, 0 ≤ x < 180, 0 ≤ y < 360;

[0143] S3-7: Calculate the disturbance optimization matrix Z t,0 (simplified as Z):

[0144] Z (:,0) = Q T - mean (V T ) T

[0145] Z (:,1) = U T

[0146] Z (:,0) = Z0

[0147] Z (:,1) = Z1

[0148] where, when depth = 0, Z is Z 11023×2 ;

[0149] S3-8: Define the scale vector T of multi-scale cyclic reconstruction: The grid point processing is carried out on the ocean heat content data, and the multi-scale cyclic reconstruction vector T is set to T = [3, 8, 25].

[0150] S3-9: Multiscale cyclic reconstruction by scale: Assume scale index is jndex, if jndex exceeds the last scale in scale vector, then depth = depth + 1 go to S3-1, otherwise continue S3-10 to data griding.

[0151] S3-10: Calculate H t,0 (abbreviated as H) and N t,0 (abbreviated as N), set the current scale of the range of the grid radius:

[0152]

[0153] N = R x Z0

[0154] H = R x (U - repmat(mean(V T , 1, m)) T , 1, m))

[0155] where radius = T jndex is the current grid scale, X in diag(X) is a column vector, which represents the diagonal matrix constructed with the elements of X as the diagonal elements, when depth = 0, R is R 11023×11023 , N is N 11023×1 , H is H 11023×40 ;

[0156] S3-11: Grid assimilation calculation by index: Assume grid index is index, if index exceeds the last grid (i.e. index ≥ (d3*d4)), then jndex = jndex + 1 go to S3-9, otherwise continue S3-12 to data griding.

[0157] S3-12: Calculate distance matrix F t,0,0 (abbreviated as F):

[0158] If nx ≤ 1,

[0159] F = abs(P - index*I)

[0160] Otherwise

[0161]

[0162] where nx = 360 when the resolution of input data is 1°x1°, index is the grid index, X = index%nx*I - P (:,0) , Y = (index / nx)*I - P (:,1) , P (:,0) , P (:,1) is the matrix P, respectively, the first column and the second column of all rows of the matrix elements constitute, when index = 0, F is calculated

[0163] S3-13: Compute localized radius matrix J t,0,0 (abbreviated as J):

[0164] J = exp(-0.5*F'.*F')

[0165] where F' = F / radius, when index = 0 J is computed

[0166] S3-14: Find localized observation matrix K t,0,0 (abbreviated as K):

[0167] Compute elements in J that are greater than 0 to get the localized observation matrix K:

[0168] K = find(J>0)

[0169] when index = 0 K is computed

[0170] S3-15: Update localized radius matrix J t,0,0 (abbreviated as J):

[0171] J = J (K,:)

[0172] Take the value of elements in K as row index, take out all the elements of all columns in the specified rows in J, update the localized radius matrix J;

[0173] S3-16: Update compute background disturbance field data A t,0,0 (abbreviated as A):

[0174] S = H (K,:) .*repmat(J,1,m)

[0175] G = (I+S T ×S) -1

[0176] A (0,:) = A (0,:) ×sqrt(G)

[0177] S3-17: Compute background mean field data B t,0,0 (abbreviated as B):

[0178] B (0,:) = C (0,:) ×(G×S T )×(N (K,:) .*J)

[0179] S3-18: Update background mean field data Bt,0,0 (abbreviated as B):

[0180] B = B + mean(O T ) T

[0181] S3-19: integrate the average data field X of the specified depth t,0 and the perturbation data field Y t,0 ; for the convenience of writing, the following two calculation formulas use to represent B t,0,0 , to represent A t,0,0 ;

[0182]

[0183]

[0184] wherein index is the grid point index, lat=index / 360 is the longitude grid index, lon=index%360 is the latitude grid index, and model is the CMIP model index, 0model<40. When model is sequentially increased, the complete Y t,0,0 can be obtained; when index is cyclically increased, the complete Y t,0 , X t,0 can be obtained;

[0185] S3-20: after the grid point index assimilation is completed, index=index+1, and S3-11 is entered;

[0186] S4: output the grid point data combination:

[0187] S4-1: output the complete spatial coverage set average data field after data grid point combination

[0188] S4-2: output the complete spatial coverage set perturbation data field after data grid point combination

[0189] S5: end of grid point combination.

[0190] Although the specific embodiments of the present application are disclosed for the purpose of illustrating the present application, helping to understand the content of the present application and implementing the same, it can be understood by those skilled in the art that various substitutions, changes and modifications are possible without departing from the spirit and scope of the present application and the appended claims. Therefore, the present application should not be limited to the disclosed content of the best mode, and the scope of the present application claimed is defined by the scope of the claims.

Claims

1. A method for gridding ocean observation data based on climate model, comprising the steps of: 1) reading ocean-related data of the tth month, wherein the ocean-related data comprises ocean state observation data O, climate mean state observation error data D, CMIP model data C, boundary topography data L between ocean and land, and global observation grid arithmetic mean data M; 2) Preprocessing the ocean-related data of the t-th month: first, for the grid points with Nan elements in the climate mean state observation error data D, Nan represents that there is no valid data at the grid point, and the Nan element is updated to the median of the non-Nan elements of the adjacent depth, adjacent grid points; second, if the ocean and land boundary topography data L (x,y) ≥ 1 and the climate mean state observation error data D of the grid point (x, y) is convolution update is performed on , x, y, z are the longitude, latitude and depth grid indexes respectively, 0≤x If the corresponding topography data L (p,q) = 0, If the corresponding topography data L (p,q) ≥ 1 and convolution update is performed on , p, q, r are the longitude, latitude and depth grid indexes respectively, 0≤p 3) gridding observation data layer by layer according to depth, wherein depth is depth, and the gridding is performed according to steps S3-1 to S3-20; S3-1: if depth is greater than or equal to d2, the gridding is ended; otherwise, S3-2 is entered; S3-2: calculating CMIP model ensemble average data C of the tth observation month and a specified depth d participating in gridding calculation; S3-8: defining a scale vector T of multi-scale cyclic reconstruction, wherein each scale in the scale vector T represents a local radius; S3-3: obtaining initial ocean state observation data O and initial climate mean state observation error data D according to depth depth; S3-4: calculating initial background disturbance field A = C-repmat(mean(C T ) T ,1,m); wherein m is the number of CMIP modes; S3-5: Calculate the grid point coordinate matrix P, the observation field Q, the analysis mean field U and the analysis perturbation field V which need to be data gridded; for the marine state observation data O corresponding to the grid point (x, y) (x,y) , if |O (x,y) | < 6 and L (x,y) ≥ 1, it is judged that the grid point (x, y) needs to be data gridded, then S3-6 is entered, otherwise S3-7 is entered; S3-6: Record the coordinate matrix of the calculation grid points that need to be performed: where P i is the number of rows of matrix P, recording the average data of the field U i is the number of rows of matrix U, k = x * d4 + y, A (k,:) to record the set of observation field data Q i is the number of rows of matrix Q, at this time if the observation error data D (x,y) > 0, then V (x,y) = D (x,y) , otherwise V i is the number of rows of the analysis disturbance field matrix V, M i is the number of rows of the grid arithmetic average data matrix M; S3-7: Compute the observed disturbance optimization matrix Z: Z0 = Q T - mean(V T ) T , Z1 = U T ; S3-9: performing multi-scale cyclic reconstruction step by step, wherein jndex is a scale index, if jndex exceeds the last scale in the scale vector, depth = depth + 1 is performed to return to S3-1, otherwise, S3-10 is performed; S3-10: calculating a scaled ensemble observation anomaly H of the tth observation month and the specified depth d and a scaled gain vector N of the tth observation month and the specified depth d, and setting a local radius radius; S3-11: performing assimilation calculation step by step, wherein index is a grid index, if index exceeds the last grid, jndex = jndex + 1 is performed to return to S3-9, otherwise, S3-12 is performed; S3-13: calculating a local radius matrix J of the grid of the tth observation month, the specified depth d and the index index; S3-12: Calculate the distance matrix F of the grid point of the tth observation month, specified depth d, and index index: if nx≤1, then F=abs(P-index*I); otherwise where nx is the maximum number of grid points in the longitude direction, X=index%nx*I-P (:,0) , Y=(index / nx)*I-P (:,1) ; P (:,0) , P (:,1) denote the matrices formed by the elements of the first and second columns of P, respectively; S3-14: calculating elements greater than 0 in J to obtain a local observation matrix K; S3-16: updating background disturbance field data A of the grid of the tth observation month, the specified depth d and the index index: S3-15: Update the localization radius matrix J = J (K,:) ; S3-20: after the assimilation of the grid with the index index is ended, index = index + 1 is performed to enter S3-11. where A is the element in the index row of all columns of A (index,:) = A (index,:) * sqrt(G); G = (I + S T * S) -1 , S = H (K,:) * repmat(J, 1, m); S3-17: Compute the background mean field data B for the grid point of the t-th observation month, specified depth d, and index index: where the elements of B in the index-th row of all columns of B (index,:) = C (index,:) × (G × S T ) × (N (K,:) * J); S3-18: Update background average field data B = B + mean(O T ) T ; S3-19: Integrating the average data field X for the tth observation month, specified depth d t,d and the perturbation data field Y t,d ; Reading CMIP model data as one of the bases for gridding ocean observation data, providing a covariance field with flow dependence as an independent and reliable data source.

2. The method of claim 1, wherein, The method for reading ocean state observation data of the tth month comprises the following steps: firstly, dividing ocean depth into layers, setting different time windows for different depths to obtain a time window vector; then, reading all observation data in each specified time window; and calculating observation data at time t according to the observation data in the specified time window in a specific manner, wherein the specific manner comprises algebraic mean value, weighted average value and geometric average value of the observation data in the specified time window.

3. The method of claim 1, wherein, The scaled set-observation anomaly matrix H and the scaled gain vector N are computed as: N = R x Z0, H = R x (O - repmat(mean(V T ) T ,1,m))、 represents constructing a diagonal matrix with the elements of the column vector as diagonal elements.

4. The method of claim 1, wherein, The implementation of multi-scale cyclic reconstruction is as follows: firstly, according to different marine typical spatial scales, the type of marine observation data processed, and the actual annual scale, a multi-scale vector T with a size of size is determined for multi-scale cyclic reconstruction, wherein each scale in the scale vector T represents a local radius, when the marine heat content data is processed by gridding, three scale cyclic reconstruction vectors T = [3, 8, 25] are constructed; secondly, the scale index is set as jndex, if jndex < size, the local radius radius = T jndex , J = exp (-0.5 * F ′ * F ′ ), wherein F ′ = F / radius, after the data gridding calculation of the scale T jndex is completed, jndex = jndex + 1, and the data gridding reconstruction calculation of the next scale is continued.

5. The method of claim 1, wherein, For the grid point (x, y) that meets the specified condition L (x,y) ≥ 1 and Convolution is performed on using the convolution kernel W1 i×j , and where E is a matrix with all elements being 1, x, y, z are the longitude, dimension and depth grid indexes respectively, 0 ≤ x < d3, 0 ≤ y < d4, 0 ≤ z < d2, i = 5, j = 7 are the convolution kernel sizes respectively; convolution is performed on using the convolution kernel W2 m×n , where E is a matrix with all elements being 1, p, q are grid indices, m = 9, n = 11, count is the number of non-Nan elements in the convolution window, 0≤r 6. The method of claim 1, wherein, The ocean depth is divided into 41 layers; and the time window vector Window1 = 7. The method of claim 6, wherein, ​ [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 4, 4, 8, 8, 8, 8, 10, 10, 12, 12, 12, 12, 12, 12, 18], the time window of the 1st layer to the 27th layer is 1 month, and the time windows of the 28th layer to the 41st layer are 4, 4, 8, 8, 8, 8, 10, 10, 12, 12, 12, 12, 12, 12, and 18 months respectively.

8. A climate model-based ocean observation data gridding system, characterized by, The data reading module, the data preprocessing module, and the data gridding module are included. The data reading module is configured to read marine related data of the t th month, wherein the marine related data includes grid arithmetic mean marine state observation data O, climate mean state observation error data D, CMIP model data C, boundary topography data L of the ocean and the land, and global observation grid arithmetic mean data M. Data preprocessing module. First, for the grid points in the climate mean state observation error data D with elements being Nan, where Nan indicates that there is no valid data at that grid point, update the grid point with the median of the non-Nan elements of the neighboring depth and neighboring grid points. Second, if the ocean and land boundary terrain data L of the specified grid point (x, y) (x,y) ≥1 and the climate mean state observation error data of the grid point (x, y) then perform convolution update. Here, x, y, and z are the longitude, latitude, and depth grid indices respectively, with 0 ≤ x < d3, 0 ≤ y < d4, 0 ≤ z < d2. d2 is the total number of depth layers of the ocean data, d3 is the number of grid points in the longitude direction, and d4 is the number of grid points in the latitude direction. For the simulation data of all models, with depth r and grid indices p, q in the CMIP model data C if the corresponding terrain data L (p,q) = 0, then if the corresponding terrain data L (p,q) ≥1 and then perform convolution update. Here, p, q, and r are the longitude, latitude, and depth grid indices respectively, with 0 ≤ p < d3, 0 ≤ q < d4, 0 ≤ r < d2. The data gridding module is configured to perform gridding processing on the observation data layer by layer according to depth, wherein the depth is denoted as d. S3-1: If the depth is greater than or equal to d2, the gridding is ended; otherwise, S3-2 is entered. S3-2: The CMIP model ensemble average data C of the t th observation month and the specified depth d participating in the gridding calculation is calculated. S3-3: obtaining initial ocean state observation data O and initial climate mean state observation error data D according to depth depth; S3-4: calculating initial background disturbance field A = C-repmat(mean(C T ) T ,1,m); wherein m is the number of CMIP modes; S3-5: Calculate the grid point coordinate matrix P, the observation field Q, the analysis mean field U and the analysis perturbation field V which need to be data gridded; for the marine state observation data O corresponding to the grid point (x, y) (x,y) , if |O (x,y) | < 6 and L (x,y) ≥ 1, it is judged that the grid point (x, y) needs to be data gridded, then S3-6 is entered, otherwise S3-7 is entered; S3-6: Record the coordinate matrix of the calculation grid required: where P i is the number of rows of matrix P, recording the analysis field average data U i is the number of rows of matrix U, k = x*d4 + y, A (k,:) is the element of the k-th row of all columns of matrix A; recording the set observation field data Q i is the number of rows of matrix Q, at this time if observation error data D (x,y) > 0, then V (x,y) = D (x,y) , otherwise V i is the number of rows of analysis disturbance field V, M i is the number of rows of grid arithmetic average data M; S3-7: Compute the observed disturbance optimization matrix Z: Z0 = Q T - mean(V T ) T , Z1 = U T ; S3-8: A scale vector T of the multi-scale cyclic reconstruction is defined. Each scale in the scale vector T represents a local radius. S3-9: The multi-scale cyclic reconstruction is performed scale by scale, wherein a scale index is denoted as jndex, if the jndex exceeds the last scale in the scale vector, the depth is set as depth + 1 and the process is transferred to S3-1, otherwise, S3-10 is performed. S3-10: The scaled ensemble observation anomaly H of the t th observation month and the specified depth d and the scaled gain vector N of the t th observation month and the specified depth d are calculated, and a local radius radius is set. S3-11: The assimilation calculation is performed grid point by grid point, wherein a grid point index is denoted as index, if the index exceeds the last grid point, the jndex is set as jndex + 1 and the process is transferred to S3-9, otherwise, S3-12 is performed. S3-12: Calculate the distance matrix F of the grid point of the tth observation month, specified depth d, and index index: if nx≤1, then F=abs(P-index*I); otherwise where nx is the maximum number of grid points in the longitude direction, X=index%nx*I-P (:,0) , Y=(index / nx)*I-P (:,1) ; P (:,0) , P (:,1) denote the matrices formed by the elements of the first and second columns of P, respectively; S3-13: A local radius matrix J of the grid point of the t th observation month, the specified depth d, and the index index is calculated. S3-14: The elements greater than 0 in J are calculated to obtain a local observation matrix K. S3-15: Update the localization radius matrix J = J (K,:) ; S3-16: The background disturbance field data A of the grid point of the t th observation month, the specified depth d, and the index index is updated. where A is the element in the index row of all columns of A (index,:) = A (index,:) * sqrt(G); G = (I + S T * S) -1 , S = H (K,:) * repmat(J, 1, m); S3-17: Compute the background mean field data B for the grid point at the t-th observation month, the specified depth d, and the index index: where the elements B (index,:) = C (index,:) × (G × S T ) × (N (K,:) * J) for all columns in the index-th row of B. S3-18: Update background average field data B = B + mean(O T ) T ; S3-19: Integrating the average data field X for the tth observation month, specified depth d t,d and the perturbation data field Y t,d ; S3-20: After the assimilation of the grid point with the index index is completed, the index is set as index + 1 and the process is transferred to S3-11.

9. A server, characterized by The computer program is configured to be executed by the processor, and the computer program includes instructions for performing each step in the method of any one of claims 1 to 7.

10. A computer-readable storage medium having stored thereon a computer program, characterized in that, The computer program is configured to be executed by the processor, and the computer program includes instructions for performing each step in the method of any one of claims 1 to 7.

Images