A Tidal Prediction Method Based on Kriging Interpolation

Abstract

The present invention belongs to the field of spatial interpolation technology, and relates to a tidal prediction method based on a Kriging interpolation algorithm. The historical tidal height data of known tide gauge stations after preprocessing are subjected to harmonic analysis to obtain the tidal characteristic information of each tide gauge station. The distance between each tide gauge station is calculated, and the semivariance function of the harmonic constant of the tide gauge station and the distance relationship are fitted, and then the semivariance function of the harmonic constant is determined. The semivariance function is used as a weight to perform interpolation calculation on the position to be predicted, and its harmonic constant is determined. Then, the tidal height information at any time is predicted according to the tidal harmonic forecasting method. Finally, it is proved through experiments that the present invention proposes a feasible spatial interpolation tidal prediction method.

Description

Technical Field

[0001] The invention belongs to the technical field of spatial interpolation and relates to a tidal prediction method based on a Kriging interpolation algorithm. Background Art

[0002] The regular rising and falling movement of seawater in the ocean due to the gravitational effects of the moon and the sun is called tide. Tide, like other ocean hydrological conditions, is one of the necessary factors for people to engage in ocean-related activities. Mastering the laws of ocean tides can ensure the safe navigation of ships. For example, when ships enter and leave some ports with large tidal changes, in order to avoid them running aground or breaking anchor, they must wait until the tide level meets the requirements. In addition, marine military activities, marine fisheries, marine environmental protection and tidal energy utilization are all closely related to tidal laws. Therefore, the study of tidal prediction is of great significance.

[0003] As early as 1687, Isaac Newton first proposed the theory of equilibrium tides, which explained many of the most basic phenomena of tides. However, the complex actual ocean conditions made many phenomena inconsistent with the theory of equilibrium tides. In 1775, Laplace established the fluid dynamics equations for describing tidal motion by assuming that the earth was completely surrounded by an ocean of equal depth. Subsequently, Thomson designed a harmonic analysis method in 1868. In 1886, Darwin further harmonically expanded the tidal force and obtained the frequencies of the main tidal components. In 1921, Dudson further expanded the tidal force into a purely harmonic expansion, thereby obtaining quite accurate calculation results. In 1928, he obtained 60 quite accurate harmonic constants of tidal components by using as little calculation as possible. So far, the harmonic analysis method has been widely used in tidal prediction applications. At present, most of the research on tidal prediction is based on the historical data of a given tide gauge station, and the future tidal height of the observation point is predicted by harmonic analysis and other methods, which belongs to time series analysis. However, considering that there are no tide gauges in many places in the ocean and near many ports, historical data on tide heights in most sea areas is missing, which makes it impossible to make tidal predictions.

[0004] In the 1990s, satellite launch technology was well developed. Satellite altimeter data represented by TOPEX / Poseidon (T / P) and Jason were widely used in ocean tide research. Based on the hydrodynamic equations, satellite altimeter data and tide station data were assimilated using data fusion and assimilation methods to obtain some ocean tide models. These ocean tide models can provide tidal characteristic information for any area in the ocean in the form of two-dimensional grid data, and tide prediction can be performed based on tidal characteristic information. The TPXO series model established by Oregon State University in the United States is based on the Laplace tidal equation and uses the least squares method to assimilate T / P and Jason satellite altimeter data. The Space Research Center of the University of Texas in the United States established the CSR series model based on empirical algorithms and T / P satellite data. The National Astronomical Observatory of Japan established the NAO series model based on the two-dimensional nonlinear shallow water equation and the blending method, and assimilated 5 years of T / P satellite data. The FES series model is a global finite element ocean tide model developed by the French Tidal Group. The University of Hamburg in Germany established the HAMTIDE series model based on the generalized inversion method and the least squares method. Based on the response method, the Technical University of Denmark assimilated T / P and Jason satellite data to establish the DTU series of models. NASA Goddard Space Flight Center established the GOT series of models using T / P and ERS satellite data. The above tidal models have better tidal prediction effects in deep water areas in the center of the ocean. However, the prediction effects in shallow water areas are different and the prediction errors are large. On the one hand, this is because in shallow water areas, tides are affected by some nonlinear factors, such as seabed friction, ocean current convergence, etc. On the other hand, different models use different assimilation methods, assimilated satellite altimeter data, and tide station data are also different. For places with tide station data or abundant tide stations, the prediction is good and the assimilation effect is good. In response to the above problems, there is an urgent need for a method to predict tides for stations without historical tide station data based on the historical data of surrounding tide stations. Summary of the invention

[0005] The present invention proposes a method for tidal prediction for locations without historical data of tide gauge stations based on the Kriging interpolation algorithm and the use of the tidal harmonic analysis method. The historical tidal height data of known tide gauge stations after preprocessing are harmonized and analyzed to obtain the tidal characteristic information of each tide gauge station. The distance between each tide gauge station is calculated, and the semivariance function of the harmonic constant of the tide gauge station and the distance relationship are fitted, and then the semivariance function of the harmonic constant is determined. The semivariance function is used as a weight to perform interpolation calculation on the position to be predicted, and its harmonic constant is determined, and then the tidal height information at any time is predicted according to the tidal harmonic forecasting method. Finally, it is proved through experiments that the present invention proposes a feasible spatial interpolation tidal prediction method.

[0006] The technical solution of the present invention:

[0007] A tidal prediction method based on Kriging interpolation uses the historical tidal height data and longitude and latitude of multiple tide gauge stations in the past at least 24 hours as the input of the Kriging interpolation model, and predicts the tidal height value at any time at the desired location as the output. The specific steps are as follows:

[0008] (1) For the tidal height data provided by the tide gauge station, data preprocessing is performed. For missing or abnormal data, that is, when the value is Nan, the data mean is used for interpolation to fill the data and extract the longitude and latitude information of the tide gauge station;

[0009] (2) Perform a harmonic analysis on the preprocessed tidal height data, where the tidal level expression is:

[0010]

[0011] Where h(t) is the tidal height at time t; S0 is the average sea level during the observation period; J is the tidal number; σ j V is the tidal angle rate; j is the tidal phase angle; j is the intersection factor, u j is the intersection correction angle of the astronomical phase angle, and the two represent the correction values ​​of the average amplitude and phase angle caused by the 18.61-year change of the lunar orbit; H j , g j The tidal harmonic constants are collectively called the tidal harmonic constants, which reflect the response of the ocean to this frequency external force. The tidal harmonic constants H are calculated using the harmonic analysis toolkit T_TIDE by R. Pawlowicz. j and g j ; They reflect the ocean's response to this frequency external force. This response is determined by the dynamic properties of the ocean itself. Since the ocean environment changes very slowly, the harmonic constants have great stability for a certain sea area and can be sufficiently approximated to be constants in a relatively short period of time.

[0012] (3) Perform Kriging interpolation on the calculated harmonic constants of some tide gauge stations to solve the harmonic constants of the points to be solved. The Kriging interpolation algorithm is also called the spatial autocovariance optimal interpolation method. It is an optimal interpolation method based on the semivariogram function theory and structural analysis, and is suitable for regional variables with spatial correlation. Considering that the tides within a certain sea area are correlated with each other in terms of geographic spatial distribution, it is assumed that they all have spatial correlation and all random errors have second-order stationarity. The expression of the Kriging interpolation algorithm is:

[0013]

[0014] In the formula, is the estimated value at the point (x0, y0), that is, z0 = z(x0, y0); λi is the weight coefficient; z i is the attribute value where there is observation data;

[0015] Use the weighted sum of the data of all known points in space to estimate the value of the unknown point, and its weight coefficient is to satisfy the estimated value at the point (x0, y0) The coefficient with the smallest difference from the true value z0, that is At the same time, it satisfies the conditions for unbiased estimation Assume that the spatial attribute z is uniform, and for any point (x, y) in space, there is the same expectation c and variance σ 2 ; That is, for any point (x, y) Var[z(x,y)]=σ 2 , select the ordinary Kriging interpolation method; Kriging interpolation is to determine the law of the change of the research object with the spatial position based on the sample points, so as to infer the attribute value of the unknown point; this law is the semivariance function, which is defined as:

[0016]

[0017] Among them, the covariance of any two points is C ij = Cov(z i ,z j ), so the expression for solving the weight coefficient of the ordinary Kriging interpolation algorithm is:

[0018]

[0019] Among them, φ is the Lagrange multiplier, and the weight coefficient λ is obtained in the above calculation i The system of equations; written in matrix form is:

[0020]

[0021] The weight coefficient can be solved by inverting the matrix. The only unknown is the semivariance function. According to the first law of geography, variables that are close in space have similar attributes. The semivariance function expresses the similarity of attributes. The spatial similarity is expressed by distance. Define z i (x i ,y i ) and z j (x j ,y j ) Kriging interpolation assumptions ij With d ijThere is a functional relationship, which can be linear, quadratic, exponential, or logarithmic. In order to confirm this relationship, we need to first analyze the observed data set {z(x1,y1),z(x2,y2),…,z(x n-1 ,y n-1 ),z(x n ,y n )} Calculate the distance d between any two points ij and semivariogram function r ij , and get n 2 (d ij ,r ij ) data pairs, plot all d and r into a scatter plot, find an optimal fitting curve to fit the relationship between d and r, and get the functional relationship r = r (d). Then for any two points (x i ,y i ),(x j ,y j ), first calculate the distance d ij , and then according to the obtained functional relationship, we can get the semivariance r of these two points ij .

[0022] After the semivariogram function is determined, the weight coefficient can be further determined according to formula (6), and then the variable value of the point to be determined can be determined by formula (4).

[0023] (4) After calculating the harmonic constant information of the desired point using the Kriging interpolation algorithm in step (3), the tidal harmonic prediction method can be used to predict the tidal height value at any time at the desired location.

[0024] The basic principle of tidal harmonic forecasting method is the inversion of tidal harmonic analysis, that is, determining the partial tide by the tidal harmonic constant, and then approximating the actual tide by the superposition of these partial tides. In this study, the astronomical partial tide is mainly considered, and the main 8 partial tides are selected for tidal analysis and forecasting. The angular frequencies of these 8 partial tides are as follows:

[0025] Table 1 Main tidal angle frequencies

[0026]

[0027] Where cph (cycles per hour) is the number of cycles per hour.

[0028] Beneficial effects of the present invention: The present invention provides a method for predicting tides using Kriging interpolation, which (1) can perform Kriging interpolation on the tidal height in the entire sea area based on the historical tidal height data of some tide gauge stations in a certain sea area, with high prediction accuracy; (2) can perform harmonic analysis on the historical tidal height data of some tide gauge stations in a certain sea area to obtain the harmonic constant, perform Kriging interpolation on the harmonic constant of the entire sea area, and then predict the tidal height at any period in the sea area, with high prediction accuracy. BRIEF DESCRIPTION OF THE DRAWINGS

[0029] Figure 1 A framework structure diagram of a tidal prediction method based on Kriging interpolation according to the present invention;

[0030] Figure 2 Check the histogram for data distribution;

[0031] Figure 3 Normal QQPlot distribution diagram for data distribution test;

[0032] Figure 4 is the semivariogram;

[0033] Figure 5 is the semivariogram fitting graph;

[0034] Figure 6 This is a comparison chart of the hourly tide heights interpolated by Kriging and the observed values ​​at Shipu throughout the day on January 1, 2021;

[0035] Figure 7 This is a comparison chart of the hourly tidal height predicted by the harmonic constant of Kriging interpolation and the observed value at Shipu throughout the day on July 23, 2022. DETAILED DESCRIPTION

[0036] The specific implementation of the present invention is further described below in conjunction with the accompanying drawings and technical solutions.

[0037] The overall process of the present invention is as follows Figure 1 As shown. The historical tidal height data and station information of 39 tide gauge stations in the East China Sea in 2021 were obtained from the National Marine Information Center. The historical tidal height data of each station were preprocessed, and the data were interpolated and filled in by using the data mean method. The time and tidal height values ​​in the tidal height data were converted into the format required for subsequent calculations. The longitude and latitude information of the tide gauge stations were extracted.

[0038] To conduct the Kriging interpolation experiment, we first determine the interpolation area based on the longitude and latitude information of the tide gauge station, and select the sea area [120.6°E~122.8°E, 27.7°N~31.5°N] for surface interpolation. Here, the tide gauge station Shipu (121.95°E, 29.2°N) is used as an unknown point, and the other tide gauge stations are used as sample points for the experiment. The spatial correlation between each sample point is calculated, and the calculation form of the semivariogram function formula (5) is:

[0039]

[0040] in is the distance between sample points, n is the number of pairs of sample points separated by d, and z is the attribute value of the point. First, find the distance between all sample points. different distances, and then sort all the distances from small to large into n groups, with a group spacing of Finally, the average value of the semivariance function corresponding to each group of distances and the average value of the distance corresponding to the semivariance function are calculated, and the relationship between the semivariance function and the distance is fitted.

[0041] The ordinary kriging interpolation method requires that the data obey the normal distribution, and the data distribution needs to be tested. Figure 2 Data distribution test histogram shows that the distribution in the histogram is approximately bell-shaped, the mean 4.0596 and the median 3.9933 are very close, the skewness value 0.14068 is around 0, and the kurtosis value 2.8065 is around 3. Figure 3 Data distribution test Normal QQPlot distribution diagram shows that the data distribution is close to a straight line with the standard normal distribution. So the experimental data is subject to normal distribution.

[0042] Explore spatial structure, using pairs of sample locations, examine spatial autocorrelation in the data and plot it as a variogram cloud. Figure 4 The semivariogram shows that close pairs of points have smaller differences. As the distance increases, the square of the difference increases and the value on the ordinate is higher. It can be seen that the experimental data has spatial autocorrelation, and the ordinary Kriging interpolation method is applicable.

[0043] Fit the semivariogram model and use the relationship between distance and semivariogram to get the best fitting theoretical semivariogram model. Here we use the exponential model. The best fitting effect is Figure 5 As shown. After obtaining the semivariogram function, the selected area is interpolated. By performing ordinary Kriging interpolation on the tidal height data of the sample tide gauge station on January 1, 2021, the tidal height value of Shipu tide gauge station is obtained. Figure 6 This is a comparison chart between the ordinary Kriging interpolation results and the actual observed values ​​at Shipu.

[0044] The data of each tide station are minute-by-minute tidal height values, and there are 525,600 data throughout the year. Given the latitude information of the tide station, the start time of the harmonic analysis, and the tidal height data, specify the specific 8 astronomical tides mentioned above, and use the T_Tide toolkit to perform tidal harmonic analysis. The harmonic constants of the 8 astronomical tides of each tide station can be obtained, and each tide station has 16 harmonic constants. Repeat the above ordinary Kriging interpolation method, and perform surface interpolation on these 16 harmonic constants in the sea area of ​​[120.6°E~122.8°E,27.7°N~31.5°N]. The harmonic constants at Shipu can be obtained as shown in the following table:

[0045] Table 2 Main tidal harmonic constants

[0046]

[0047] The main tidal harmonic constants at Shipu are obtained by ordinary kriging interpolation. The tidal heights of any future period can be predicted based on the harmonic constants and the longitude and latitude information of Shipu. Given the main tidal angular frequencies in Table 1, the latitude information of the tide gauge station, the start time of the harmonic forecast and the length of the prediction time series, as well as the main tidal harmonic constant values ​​in Table 2, the T_Tide toolkit is used for harmonic forecasting. A comparison chart of the hourly tidal heights predicted for the entire day on July 23, 2022 based on the interpolated harmonic constants and the observed values ​​can be obtained, as shown in Figure 7 . Method 1 calculates the tidal height at Shipu by ordinary Kriging interpolation of the tidal heights at other tide gauges; Method 2 calculates the harmonic constant by interpolation, and then predicts the tidal height at Shipu by the harmonic forecasting method. The root mean square errors (RMSE) of the two methods are 0.07961 and 0.13844, respectively, both of which are relatively small, which shows that the calculation results of the two methods are relatively accurate. Although Method 1 is more effective, Method 1 requires interpolation based on historical tidal height data for each calculation, which is complicated to operate and cannot predict tidal heights. Method 2 first harmonically analyzes a large amount of historical tidal height data to obtain the harmonic constant, and then performs an interpolation calculation on the harmonic constant. It can predict the tidal height at any period and is more practical.

Claims

1. A tidal prediction method based on Kriging interpolation, which uses the historical tidal height data and longitude and latitude of multiple tide gauge stations in the past at least 24 hours as the input of the Kriging interpolation model, and predicts the tidal height value at any time at the desired location as the output, characterized in that: The specific steps are as follows: (1) For the tidal height data provided by the tide gauge station, data preprocessing is performed. For missing or abnormal data, that is, when the value is Nan, the data mean is used for interpolation to fill the data and extract the longitude and latitude information of the tide gauge station; (2) Perform a harmonic analysis on the preprocessed tidal height data, where the tidal level expression is: Where h(t) is the tidal height at time t; S0 is the average sea level during the observation period; J is the tidal number; σ j V is the tidal angle rate; j is the tidal phase angle; j is the intersection factor, u j is the intersection correction angle of the astronomical phase angle, and the two represent the correction values ​​of the average amplitude and phase angle caused by the 18.61-year change of the lunar orbit; H j , g j The tidal harmonic constants are collectively called the tidal harmonic constants, which reflect the response of the ocean to this frequency external force. The tidal harmonic constants H are calculated using the harmonic analysis toolkit T_TIDE by R. Pawlowicz. j and g j ; (3) Perform Kriging interpolation on the calculated tidal harmonic constants of the tide gauge station to obtain the harmonic constants of the points to be determined. The expression of the Kriging interpolation algorithm is: In the formula, is the estimated value at the point (x0, y0), that is, z0 = z(x0, y0); λ i is the weight coefficient; z i is the attribute value where there is observation data; Use the weighted sum of the data of all known points in space to estimate the value of the unknown point, and its weight coefficient is to satisfy the estimated value at the point (x0, y0) The coefficient with the smallest difference from the true value z0, that is At the same time, it satisfies the conditions for unbiased estimation Assume that the spatial attribute z is uniform, and for any point (x, y) in space, there is the same expectation c and variance σ 2 ; That is, for any point (x, y), E[z(x, y)] = E[z] = c, Var[z(x, y)] = σ 2 , select the ordinary Kriging interpolation method; Kriging interpolation is to determine the law of the change of the research object with the spatial position based on the sample points, so as to infer the attribute value of the unknown point; this law is the semivariance function, which is defined as: Among them, the covariance of any two points is C ij = Cov(z i ,z j ), so the expression for solving the weight coefficient of the ordinary Kriging interpolation algorithm is: The semivariance function expresses the similarity of attributes, and the spatial similarity is expressed by distance. Define z i (x i ,y i ) and z j (x j ,y j ) Kriging interpolation assumptions ij With d ij There is a functional relationship; in order to confirm this relationship, we need to first analyze the observed data set {z(x1,y1),z(x2,y2),…,z(x n-1 ,y n-1 ),z(x n ,y n )} Calculate the distance d between any two points ij and semivariogram function r ij , and get n 2 (d ij ,r ij ) data pairs, fitting the semivariogram function; (4) The tidal harmonic constant of the target point is determined by the Kriging interpolation algorithm in step (3), and then the tidal harmonic prediction method is used to predict the tidal height value at any time at the target location.

2. The tidal prediction method according to claim 1, characterized in that: In the process of calculating the harmonic constant, the actual tide level is composed of many partial tides superimposed, so formula (1) is equivalent to: From the cosine theorem, we know that In the formula, a j =R j cosθ j , b j =R j sinθ j , the amplitude without introducing the intersection factor Delay angle without introducing intersection correction angle Combining formula (1), we have

Images