A sea surface height measurement correction method, device, computer equipment and program product

Abstract

The present disclosure provides a sea surface height measurement correction method, device, computer equipment and program product, wherein the method comprises: obtaining each time window; for any time window, performing adaptive spectral decomposition and wave system signal extraction on the vertical displacement in each three-way displacement to obtain the time window in the band-pass frequency band corresponding to the wind wave component and the swell component and the significant wave height; determining the wind wave component direction and the swell component direction according to the mutual covariance between the three-way displacements, and determining the main frequency of the wind wave component and the main frequency of the swell component according to the band-pass frequency band; determining the wave number of the wind wave component and the wave number of the swell component through the wave number iteration mode according to the projection component of the background flow field in the wind wave component direction and the swell component direction, and the main frequency of the wind wave component and the main frequency of the swell component; determining the vertical offset according to the wave number of the wind wave component, the wave number of the swell component and the significant wave height, and determining the height sequence of the sea surface by using the vertical offset and each vertical displacement.

Description

TECHNICAL FIELD

[0001] The present disclosure relates to the field of ocean technology, and in particular, to a sea surface height measurement correction method and device, computer equipment and program product. BACKGROUND

[0002] As an important observation direction in modern ocean observation, sea surface height measurement can realize the analysis of ocean circulation patterns, tidal laws and ocean eddy phenomena through the observed sea surface height changes, thereby providing a basis for understanding ocean dynamic processes. Current sea surface height measurement methods mainly include two types: one is based on the Euler system of fixed spatial reference system for measurement, typical measurement means such as satellite altimeter and tide gauge; the other is based on Lagrange system for measurement, typical measurement means such as Global Navigation Satellite System (GNSS) buoy measurement.

[0003] The sea surface height obtained based on GNSS buoy measurement will be offset from the sea surface height obtained based on satellite altimeter and tide gauge due to the change of sea surface waves, so in order to improve the accuracy of GNSS buoy measurement, the offset caused by the change of sea surface waves needs to be eliminated. However, the current offset elimination methods often have the problem of distortion and even amplification of offset estimation, resulting in low accuracy and poor effect of offset elimination, which is obviously disadvantageous. SUMMARY

[0004] The embodiments of the present disclosure at least provide a sea surface height measurement correction method, device, computer equipment and program product to eliminate the offset in GNSS buoy observation and improve the accuracy of sea surface height measurement.

[0005] In a first aspect, the embodiments of the present disclosure provide a sea surface height measurement correction method, comprising:

[0006] Obtaining each time window under sea surface buoy data; each time window comprises a plurality of continuously collected sea surface buoy data, and the sea surface buoy data comprises three-direction displacement;

[0007] For any time window, performing adaptive spectral decomposition and wave system signal extraction on the vertical displacement in each three-direction displacement to obtain the band-pass frequency band and significant wave height corresponding to the wind wave component and the swell component of the time window, respectively;

[0008] Determining the wind wave component direction and the swell component direction of the time window according to the cross-covariance between the three-direction displacements, and determining the wind wave component main frequency and the swell component main frequency of the time window according to the band-pass frequency band;

[0009] According to the projection components of the background flow field of the sea surface in the wind wave component direction and the swell component direction, and the main frequency of the wind wave component and the main frequency of the swell component, the wave number iteration mode is used to determine the wind wave component wave number and the swell component wave number corresponding to the time window.

[0010] According to the wind wave component wave number, the swell component wave number, and the significant wave height, the vertical offset of the time window is determined, and the height sequence of the sea surface is determined by using the vertical offset and each vertical displacement.

[0011] In a second aspect, the embodiments of the present disclosure further provide a sea surface height measurement correction device, comprising:

[0012] The acquisition module is configured to acquire each time window under sea surface buoy data. Each time window includes a plurality of continuously collected sea surface buoy data, and the sea surface buoy data includes three-way displacement.

[0013] The decomposition module is configured to perform adaptive spectral decomposition and wave system signal extraction on the vertical displacement in each three-way displacement for any time window, to obtain the band-pass frequency band and the significant wave height corresponding to the wind wave component and the swell component of the time window.

[0014] The determination module is configured to determine the wind wave component direction and the swell component direction of the time window according to the mutual covariance between each three-way displacement, and determine the main frequency of the wind wave component and the main frequency of the swell component of the time window according to the band-pass frequency band.

[0015] The iteration module is configured to determine the wind wave component wave number and the swell component wave number corresponding to the time window by using the wave number iteration mode according to the projection components of the background flow field of the sea surface in the wind wave component direction and the swell component direction, and the main frequency of the wind wave component and the main frequency of the swell component.

[0016] The correction module is configured to determine the vertical offset of the time window according to the wind wave component wave number, the swell component wave number, and the significant wave height, and determine the height sequence of the sea surface by using the vertical offset and each vertical displacement.

[0017] In a third aspect, the optional implementation manner of the present disclosure further provides a computer device, a processor and a memory. The memory stores machine readable instructions executable by the processor. The processor is configured to execute the machine readable instructions stored in the memory. When the machine readable instructions are executed by the processor, the steps of the first aspect are executed.

[0018] In a fourth aspect, the optional implementation manner of the present disclosure further provides a computer program product, comprising a computer program. When the computer program is executed, the steps of the first aspect are implemented.

[0019] The sea surface height measurement correction method, device, computer device and program product provided by the embodiments of the present disclosure can automatically identify the main frequency, band-pass frequency band, significant wave height and direction corresponding to the wind wave and the swell wave respectively by performing adaptive spectral decomposition and wave series signal extraction on the vertical displacement in the sea surface buoy data, so as to flexibly adapt to the real wave situation under different sea conditions and improve the accuracy of the extracted wave signal. By introducing the projection correction of the background flow field, the error caused by the frequency shift of the sea wave signal can be eliminated, and the accuracy of the deduced wave number can be improved. The use of more real and accurate wave numbers and significant wave heights for offset determination can improve the matching degree of the determined vertical displacement and the real sea condition, and obtain an accurate vertical displacement. The vertical displacement in the sea surface buoy data is corrected by using the vertical displacement, so that the centimeter-level vertical displacement correction can be realized, the systematic deviation caused by the prolonged wave peak residence time of the buoy can be effectively eliminated, and thus the precision and reliability of the offset elimination can be improved, and the effect of the buoy data height measurement can be improved.

[0020] For the effect description of the sea surface height measurement correction device, computer device and computer program product, refer to the description of the sea surface height measurement correction method, which will not be repeated here.

[0021] In order to make the above-mentioned purposes, features and advantages of the present disclosure more obvious and easy to understand, the following preferred embodiments are described in detail below, and the accompanying drawings are described as follows. BRIEF DESCRIPTION OF DRAWINGS

[0022] In order to more clearly illustrate the technical solutions of the embodiments of the present disclosure, the following will briefly introduce the drawings needed to be used in the embodiments. The drawings are incorporated into the specification and form a part of the specification. These drawings show the embodiments consistent with the present disclosure, and are used to illustrate the technical solutions of the present disclosure together with the specification. It should be understood that the following drawings only show some embodiments of the present disclosure, and therefore should not be regarded as a limitation on the scope, and other related drawings can also be obtained by those skilled in the art without creative labor.

[0023] Figure 1 A flow chart of a sea surface height measurement correction method provided by the embodiments of the present disclosure is shown;

[0024] Figure 2 A time series comparison chart of the original observation value of the GNSS buoy and the reference value of the same position sea bottom tide gauge provided by the embodiments of the present disclosure is shown;

[0025] Figure 3 A difference comparison chart provided by the embodiments of the present disclosure is shown;

[0026] Figure 4A schematic diagram of a sea surface height measurement correction device is shown.

[0027] Figure 5 A structural schematic diagram of a computer device is shown. DETAILED DESCRIPTION

[0028] To make the objectives, technical solutions and advantages of the embodiments of the present disclosure clearer, the technical solutions in the embodiments of the present disclosure will be described clearly and completely below with reference to the drawings in the embodiments of the present disclosure. Obviously, the described embodiments are only part of the embodiments of the present disclosure, rather than all the embodiments. The components of the embodiments of the present disclosure described and shown herein can be arranged and designed in various different configurations. Therefore, the following detailed description of the embodiments of the present disclosure is not intended to limit the scope of the claimed present disclosure, but only represents selected embodiments of the present disclosure. Based on the embodiments of the present disclosure, all other embodiments obtained by those skilled in the art without creative work belong to the scope of protection of the present disclosure.

[0029] In addition, the terms "first", "second", and the like in the specification and claims of the embodiments of the present disclosure and the above-described drawings are used to distinguish similar objects, and do not necessarily have to describe a specific order or sequence. It should be understood that the data used in this way can be interchanged under appropriate circumstances, so that the embodiments described herein can be implemented in an order other than that illustrated or described herein.

[0030] "Multiple or several" referred to herein means two or more. "And / or" describes the association relationship between the associated objects, which means that there can be three relationships, for example, A and / or B, which can represent the three cases of A existing alone, A and B existing together, and B existing alone. The character " / " generally represents that the front and rear associated objects are in an "or" relationship.

[0031] It is found that satellite altimeter and tide gauge are two core technical means of modern ocean observation and satellite calibration. They generally use Euler system under fixed spatial reference system to measure the change of sea surface height at a fixed position with time. With the development of global navigation satellite system technology, GNSS buoy observation technology has gradually become an important supplement to satellite altimetry calibration and nearshore tide monitoring due to its high precision, low cost and flexible deployment capability. It is widely used in satellite altimetry calibration, nearshore water level monitoring and ocean dynamic process observation. However, the sea surface height observation of GNSS buoy belongs to the observation under Lagrangian system, which records the height change of the buoy at a fixed position with wave orbit motion, which is naturally different from the observation under Euler system. For example, when there is wave on the sea surface, the sea surface buoy will rise with the wave crest and descend with the wave trough. Due to the nonlinear wave orbit effect, the residence time of the sea surface buoy at the wave crest is slightly longer than that at the wave trough, so that the average sea surface height observed by the buoy is higher than that measured by the Euler system, forming Stokes drift. This drift is often proportional to the square of wave height and wave number, and can reach several centimeters under medium or large wave conditions, which is at the same order of magnitude as the precision requirements of satellite altimetry and tide observation, so it has become a factor that cannot be ignored in sea surface height observation. With the satellite altimetry and calibration business entering the era of centimeter-level precision, how to accurately correct the above drift has become one of the key technical problems in high-precision sea surface height observation.

[0032] Some current GNSS buoy observation techniques do not display and correct the Stokes drift generated in the observation data processing, but only weaken the wave signal collected by the GNSS buoy through simple high-pass or band-pass filtering, or weaken the buoy motion effect through attitude compensation. Although these methods can partially reduce high-frequency noise, they cannot physically correct systematic bias errors caused by waves, ignore the bias problem, and result in centimeter-level deviations between the sea surface height (SSH) measured by the GNSS buoy and the measurements of satellite altimeters, tide gauges, and other measurements in the Euler system, which is difficult to meet the high-precision business requirements. Although some offset correction methods have appeared, these methods often have the following problems: 1. Assuming that the wave signal is smooth and single-peak, the observation data is often analyzed by Fourier spectrum, but due to the insufficient spectrum decomposition accuracy, it is difficult to adapt to the non-stationary and multi-wave characteristics of the superimposed wind waves and swell waves in actual sea conditions (such as wind waves and swell waves of different properties), and the frequency spectrum aliasing phenomenon is obvious, causing instability in the direction and amplitude of the offset estimation. 2. Ignoring the influence of the background flow field on the measurement, when the background flow field is strong or the wave nonlinearity is enhanced, the wave number inversion is prone to poor convergence or systematic deviation, resulting in an enlarged offset correction error. 3. Based on the assumption that the wave field has a single propagation direction, the offset correction is performed by a scalar method, which is inconsistent with the wind wave direction and swell direction in the real sea conditions, making it impossible to handle the case where wind waves and swell waves come from different directions, and the offset estimation is prone to distortion in complex sea conditions. 4. The offset is corrected by relying on manual pre-setting of the frequency band range of the wave signal and empirical calculation parameters, making it difficult to truly adapt to different sea conditions, lacking self-adaptability and robustness, and difficult to achieve long-term, stable, and automatic centimeter-level correction in actual business.

[0033] Based on the above research, the present disclosure provides a sea surface height measurement correction method, device, computer equipment and program product. By adaptively decomposing the frequency spectrum of the vertical displacement in the sea surface buoy data and extracting the wave signal, the main frequency, band-pass frequency band, significant wave height and direction corresponding to the wind wave and swell wave can be automatically identified for real-time sea surface buoy data, so as to flexibly adapt to the real wave conditions in different sea conditions and improve the accuracy of the extracted wave signal. By introducing the projection correction of the background flow field, the error caused by the frequency shift of the sea wave signal can be eliminated, and the accuracy of the deduced wave number can be improved. Using a more accurate wave number and significant wave height to determine the offset can improve the matching degree of the determined vertical offset and the real sea conditions, and obtain an accurate vertical offset. Using the vertical offset to correct the vertical displacement in the sea surface buoy data can achieve centimeter-level correction of the vertical displacement, effectively eliminate the systematic deviation caused by the prolonged wave peak residence time of the buoy, thereby improving the accuracy and reliability of the offset elimination and enhancing the effect of the buoy data height measurement.

[0034] The shortcomings of the above solutions are the result of the inventor's practical experience and careful research. Therefore, the discovery process of the above problems and the solutions proposed in this disclosure below should be considered as the inventor's contribution to this disclosure.

[0035] It should be noted that similar labels and letters in the following figures indicate similar items. Therefore, once an item is defined in one figure, it does not need to be further defined and explained in subsequent figures.

[0036] It is understood that before using the technical solutions disclosed in the various embodiments of this disclosure, users should be informed of the types, scope of use, and usage scenarios of the personal information involved in this disclosure in an appropriate manner in accordance with relevant laws and regulations, and user authorization should be obtained.

[0037] To facilitate understanding of this embodiment, a detailed description of the sea surface height measurement correction method disclosed in this disclosure will be provided first. The execution subject of the sea surface height measurement correction method provided in this disclosure is generally a terminal device or other processing device with certain computing capabilities. The terminal device can be a user equipment (UE), mobile device, user terminal, terminal, personal digital assistant device (PDA), handheld device, computer device, etc. In some possible implementations, the sea surface height measurement correction method can be implemented by the processor calling computer-readable instructions stored in the memory.

[0038] The sea surface height measurement correction method provided in this disclosure is illustrated below, taking a computer device as the execution subject as an example.

[0039] like Figure 1 The flowchart shown is a method for correcting sea surface height measurement according to an embodiment of this disclosure, which may include the following steps:

[0040] S101: Acquire data from various time windows of the sea surface buoy; each time window includes multiple continuously collected sea surface buoy data, which includes three-dimensional displacement.

[0041] The sea surface height measurement correction method provided in this application is a measurement method applied to GNSS buoy observation technology. It is used to eliminate the offset caused by wave changes during buoy observation, thereby improving the accuracy of GNSS buoy observation. The sea surface buoy data can be displacement data collected from the sea surface buoy, which is set at a position on the sea surface, and the number of sea surface buoys can be one.

[0042] The buoy data can include the buoy's three-dimensional displacement, which refers to displacement in three orthogonal directions. Specifically, it can include horizontal displacement in the horizontal direction, longitudinal displacement in the vertical direction, and vertical displacement in the direction perpendicular to the sea surface. For example, the buoy data uses... It indicates. Among them, This represents the horizontal displacement at time t. This represents the longitudinal displacement at time t. This represents the vertical displacement at time t.

[0043] In GNSS buoy observations, surface buoy data can be continuously acquired in the time domain. A time window is obtained by dividing the continuously acquired surface buoy data into windows, with each window corresponding to the same time length and being continuous in time. Each time window includes multiple continuously acquired surface buoy data points, and the surface buoy data included in different time windows may vary. The surface buoy data in each time window can be used... It means that, among them, , and , for the initial collection , and The result is obtained after preprocessing, and the preprocessing steps will be described later.

[0044] In practical implementation, data from various sea surface buoys collected continuously over the time domain can be acquired. The continuously collected sea surface buoy data can then be divided into windows and preprocessed according to a preset window length to obtain various time windows. The preset window length can be set empirically, and this embodiment does not impose a specific limitation. For example, the preset window length can be 30 minutes (min), 60 minutes, etc. The i-th time window can be represented as: That is, i is the identifier of the time window.

[0045] in, This represents the start time of the i-th time window. Indicates the preset window length. This represents the end time of the i-th time window. Each time within the time window corresponds to one data point of a sea buoy. Therefore, the in-window data sequence composed of sea buoy data within the i-th time window can be represented as: .in, Indicates the i-th time window Horizontal displacement over time; Indicates the i-th time window Longitudinal displacement over time; Indicates the i-th time window Vertical displacement over time.

[0046] In one embodiment, S101 described above can be implemented according to the following steps:

[0047] S101-1: Acquire the continuous output of sea surface buoy data in the time domain.

[0048] For example, it is possible to acquire raw buoy data continuously output in the time domain from buoys deployed on the surface of any sea area. .

[0049] S101-2: Divide the data of each sea surface buoy into windows according to the preset time length to obtain each initial time window.

[0050] Here, the preset time length is the preset window length, for example, 30 minutes. The initial time window is the time window obtained by dividing the original sea surface buoy data into windows. The sea surface buoy data in this time window has not yet been preprocessed.

[0051] For example, the continuously output sea surface buoy data can be divided into windows according to a preset time length to obtain various initial time windows. For instance, the i-th initial time window can be represented as: The in-window data sequence consisting of sea surface buoy data within the i-th initial time window can be represented as: . Indicates the i-th time window The original horizontal displacement over time, Indicates the i-th time window The original longitudinal displacement over time, Indicates the i-th time window The original vertical displacement over time.

[0052] S101-3: For any initial time window, determine the proportion of missing sea surface buoy data within the initial time window.

[0053] Here, the missing proportion is used to characterize the proportion of missing buoy data in the data sequence within the initial time window to the total number of buoy data in the window.

[0054] For example, if each initial time window is 30 minutes long and the data collection interval for each buoy is 10 seconds, then the total number of buoy data points within an initial time window should be 180. For any given initial time window, if the number of buoy data points within that window is 150, then the missing rate for that window is approximately 17%.

[0055] S101-4: In a case where the missing proportion is less than or equal to the preset proportion, the missing buoy data in the initial time window is completed to obtain a time window.

[0056] Here, the preset proportion can be set according to experience, and the embodiments of the present application do not make specific limitations. For example, the preset proportion can be 5%.

[0057] For example, for any initial time window, if the missing proportion of the buoy data of the window is greater than the preset proportion, it can be determined that the window is unqualified, and the window and each sea surface buoy data in the window are directly deleted. If the missing proportion of the buoy data of the window is less than or equal to the preset proportion, it can be determined that the window is qualified and, in a case where there is missing buoy data in the window, the missing buoy data is generated using an interpolation method to obtain filling data corresponding to the missing buoy data and filled in the window data sequence of the window to obtain a completed time window. The specific interpolation method is not limited in the present application, for example, it can be a cubic spline interpolation method, a nearest neighbor interpolation method, a linear interpolation method, etc. To ensure the rationality of the interpolation data, the same interpolation method can be used for each initial time window when interpolating. For example, according to the missing time corresponding to the missing buoy data, the effective buoy data within a time range before and after the missing time is interpolated to obtain the filling data.

[0058] It can be understood that if the missing proportion of the buoy data of the initial time window is 0, the initial time window can be directly taken as the completed time window.

[0059] After obtaining the completed time window, the completed time window can be directly taken as the preprocessed time window, or can be further preprocessed according to the manner described below to obtain a time window.

[0060] In an embodiment, for S101-4 described above, the following steps can be implemented:

[0061] S101-4-1: The missing buoy data in the initial time window is completed to obtain a completed time window.

[0062] In specific implementation, for any initial time window with a missing proportion less than or equal to the preset proportion, whether the window has missing buoy data is determined according to the collection time of each sea surface buoy data in the window data sequence of the window; if not, the window is directly taken as a completed time window, and the completed time window at this time includes each original sea surface buoy data. If yes, the missing data is filled using the interpolation method described above to obtain a completed time window, and the completed time window at this time includes the original sea surface buoy data and the filled sea surface buoy data.

[0063] S101-4-2: According to the mean value and linear slope of each three-direction displacement in the completion time window in each direction, the de-meaning and de-linear slope processing are respectively performed on each three-direction displacement in the completion time window, to obtain each processed three-direction displacement.

[0064] In a specific implementation, for any direction, the mean displacement corresponding to the completion time window in the direction is determined according to the displacement corresponding to each three-direction displacement in the completion time window in the direction , and the linear slope corresponding to the completion time window in the direction is obtained by fitting each displacement corresponding to the direction . Wherein, The mean displacement in the horizontal direction , the mean displacement in the longitudinal direction , and the mean displacement in the vertical direction may be included. The linear slope in the horizontal direction , the linear slope in the longitudinal direction , and the linear slope in the vertical direction may be included.

[0065] Then, for each direction, the mean displacement and the linear slope in the direction can be used to perform de-meaning and de-linear slope processing on each displacement in the direction, to obtain each processed displacement; the processed displacements in each direction constitute each processed three-direction displacement; and each processed three-direction displacement constitutes the zero-mean disturbance sequence corresponding to the completion time window.

[0066] Exemplarily, the zero-mean disturbance sequence corresponding to the completion time window can be represented by the following formula:

[0067] ;

[0068] Wherein, represents the zero-mean disturbance horizontal displacement at time t in the i th completion time window, represents the zero-mean disturbance longitudinal displacement at time t in the i th completion time window, represents the zero-mean disturbance vertical displacement at time t in the i th completion time window.

[0069] S101-4-3: According to the mean error of the completion time window in each direction, the three-direction displacement to be processed whose mean error does not meet the preset error standard is selected from each processed three-direction displacement.

[0070] ​​​Here, the preset error criterion can be set according to experience, and embodiments of the present application do not make specific limitations. For example, the preset error criterion can be that the mean error of each processed displacement is less than or equal to 5 times the mean error in the corresponding direction. The mean error is the standard deviation, and one complete time window corresponds to one mean error in one direction. For example, the mean error of the i th complete time window in the horizontal direction can be represented as ; the mean error in the vertical direction can be represented as ; and the mean error in the vertical direction can be represented as .

[0071] In specific implementation, for any complete time window, the mean error of the window in each direction can be determined according to the following formula:

[0072]

[0073] wherein, represents the total number of processed three-direction displacements in the i th complete time window, represents the mean displacement of each zero-mean perturbation horizontal displacement in the i th complete time window; represents the mean displacement of each zero-mean perturbation vertical displacement in the i th complete time window; represents the mean displacement of each zero-mean perturbation vertical displacement in the i th complete time window.

[0074] Further, for each processed three-direction displacement, if the processed three-direction displacement satisfies the condition shown in the following formula, it can be determined that the processed three-direction displacement is a gross error data point and needs to be processed, and therefore it is regarded as a three-direction displacement to be processed:

[0075]

[0076] It can be understood that if each processed three-direction displacement in the complete time window meets the preset error criterion, the complete time window can be directly regarded as a time window whose preprocessing is completed.

[0077] S101-4-4: Using an interpolation method, determining an optimized displacement corresponding to the three-direction displacement to be processed, and replacing the three-direction displacement to be processed with the optimized displacement to obtain a time window.

[0078] For example, for each three-direction displacement to be processed in the complete time window, an interpolation method can be used to generate an optimized displacement corresponding to the three-direction displacement to be processed. The optimized displacement includes optimized displacements in three directions. Then, the optimized displacement is used to replace the three-direction displacement to be processed to obtain a time window whose preprocessing is completed. The interpolation method used when generating the optimized displacement can be the same as or different from the interpolation method used when generating the filling data.

[0079] S102: For any time window, performing adaptive frequency spectrum decomposition and wave system signal extraction on the vertical displacement in each three-direction displacement to obtain a band-pass frequency band and a significant wave height corresponding to the wind wave component and the swell component respectively in the time window.

[0080] Here, the band-pass frequency band of the time window is used to provide a band-pass boundary for determining the propagation direction and wave number inversion of the wave signal. The significant wave height is used to determine the offset of the wave signal. The wave system signal can include a wind wave system signal and a swell wave system signal, and the wind wave component and the swell component are signal components under the wind wave system and the swell wave system. There are corresponding band-pass frequency bands and significant wave heights in each time window on the wind wave component and the swell component.

[0081] In specific implementation, for any time window, the Hilbert-Huang transform (HHT) can be used to perform adaptive frequency spectrum decomposition and multi-wave system signal extraction on each vertical displacement in the time window, and based on the decomposition results of each vertical displacement, wind wave component information and swell component information, the wind wave band-pass frequency band and the wind wave significant wave height of the time window are determined, and the swell band-pass frequency band and the swell significant wave height of the time window are determined.

[0082] In an embodiment, the above S102 can be implemented according to the following steps:

[0083] S102-1: For any vertical displacement, performing empirical mode decomposition on the vertical displacement according to a preset decomposition rule to obtain each intrinsic mode component, and performing frequency spectrum analysis on each intrinsic mode component to obtain the instantaneous amplitude and the instantaneous angular frequency corresponding to each intrinsic mode component.

[0084] Here, the preset decomposition rule can include three parts of mode number range limit, end point strategy and decomposition stop criterion. Among them, the mode number range limit can be that the number of intrinsic mode components (IMF) decomposed is between 3-8. The end point strategy can be to determine the end point mode by mirror extension or spline end point constraint, so as to avoid the end point effect. The decomposition stop criterion can be that the number of zero crossings is not more than 1.

[0085] In specific implementation, for any vertical displacement in any time window , empirical mode decomposition (EMD) can be performed on according to the preset decomposition rule to obtain The corresponding intrinsic modal components and residual terms. Understandably, the number of intrinsic modal components decomposed for different vertical displacements can be the same or different. For example, the intrinsic modal components and residual terms corresponding to the vertical displacement can be decomposed using the following formula:

[0086] ;

[0087] in, Indicates the number of intrinsic modal components. express The j-th intrinsic mode component that is decomposed, express The residual terms extracted from the decomposition. The residual terms extracted from each vertical displacement will be discarded and not used for subsequent processing.

[0088] Furthermore, for each intrinsic mode component (IMF) The instantaneous amplitude and instantaneous angular frequency of the IMF can be obtained by performing a Hilbert transform on it using the following formula:

[0089] ;

[0090] in, express The corresponding instantaneous amplitude, express The corresponding instantaneous angular frequency, express The corresponding instantaneous phase.

[0091] S102-2: Determine the modal energy and target period of each intrinsic modal component based on the instantaneous amplitude and instantaneous angular frequency.

[0092] Here, modal energy is used to represent the energy magnitude of the IMF; the target period is used to characterize the representative period of the IMF. In practice, for each IMF, the representative period (i.e., the target period) of that IMF within a time window can be defined using the following formula:

[0093] ;

[0094] in, Represents the j-th IMF in the i-th time window. The corresponding target period, Indicates the time window of the i-th time window of the median.

[0095] Furthermore, for each IMF, the modal energy of that IMF can be determined based on its corresponding instantaneous amplitude. .

[0096] S102-3: According to the modal energy and the target period, the inherent modal components of each vertical displacement are screened and grouped to obtain a modal group corresponding to the wind wave component and the swell component respectively, and based on the modal group and the instantaneous angular frequency, a band-pass frequency band and a component sequence corresponding to the wind wave component and the swell component respectively in the time window are determined.

[0097] Here, a time window, a modal group corresponding to the wind wave component , and a modal group corresponding to the swell component . The modal group includes each dominant modal of the time window on the corresponding component. The component sequence can include a wind wave component sequence corresponding to the wind wave component in the time window, and a wind wave component sequence corresponding to the swell component.

[0098] In specific implementation, the modal energy proportion of each IMF can be determined according to the modal energy of each IMF of the vertical displacement in the time window. For example, the modal energy proportion of each IMF can be determined according to the following formula:

[0099]

[0100] wherein, represents the modal energy proportion of the jth IMF in the ith time window, represents the modal energy square of the kth IMF corresponding to the vertical displacement at the time in the ith time window.

[0101] Then, according to the IMF corresponding to each vertical displacement, the IMF with a modal energy proportion greater than or equal to a preset energy proportion is determined as the dominant modal index set corresponding to the ith time window. Wherein, =0.10. At the same time, the IMF with < is regarded as noise / secondary modal and is removed.

[0102] Further, the target period and the preset period can be used to group the dominant modal index set to obtain the modal group and the modal group of the time window. The preset period can be 10 seconds (s). Specifically, and can be represented by the following formula:

[0103] ​​ .

[0104] Further, the wind wave component and the swell wave component corresponding to the time window can be determined according to the following formulae: and

[0105] .

[0106] It can be understood that the two component sequences here are the "two components corresponding to the current time window ", that is, the two components corresponding to the entire time window. After decomposition, screening and grouping are performed on the two components, two time sequences (i.e., the wind wave component sequence and the swell wave component sequence) are reconstructed in the entire window, and they both cover the entire time. Further, the significant wave height of the wind wave component corresponding to the time window can be determined according to the wind wave component sequence , and the significant wave height of the swell wave component corresponding to the time window can be determined according to the swell wave component sequence . Further, the band-pass frequency band corresponding to the wind wave component and the swell wave component respectively corresponding to the time window can be determined according to the of each IMF according to the following formulae:

[0107]

[0108] ;

[0109] wherein the band-pass frequency band corresponding to the wind wave component and the band-pass frequency band corresponding to the swell wave component ; and represent the minimum frequency and the maximum frequency in the band-pass frequency band respectively. P05 represents the 5th percentile, and P95 represents the 95th percentile.

[0110] It can be understood that, on the basis of assuming that the horizontal displacement, the vertical displacement and the vertical displacement have consistent energy response characteristics in the same frequency band, the band-pass frequency band obtained by processing each vertical displacement in the time window can also be used as the band-pass frequency band corresponding to the horizontal displacement and the vertical displacement for subsequent band-pass filtering and propagation direction calculation.

[0111] S103: determining the wind wave component direction and the swell wave component direction of the time window according to the cross-covariance between each three-direction displacement, and determining the wind wave component main frequency and the swell wave component main frequency of the time window according to the band-pass frequency band.

[0112] ​​​​Here, the wind wave component direction and the swell component direction are respectively used to indicate the average propagation direction of the wind wave and the average propagation direction of the swell in the time window. The cross-covariance can include the cross-covariance between the horizontal displacement and the vertical displacement corresponding to the wind wave component and the swell component respectively, and the cross-covariance between the longitudinal displacement and the vertical displacement corresponding to the wind wave component and the swell component respectively.

[0113] In a specific implementation, for any time window, the wind wave component and the swell component corresponding to each three-direction displacement can be obtained by filtering each three-direction displacement. Then, the cross-covariance between each three-direction displacement corresponding to the wind wave component and the swell component is determined according to the wind wave component of each three-direction displacement, and the cross-covariance between each three-direction displacement corresponding to the wind wave component and the swell component is determined according to the swell component of each three-direction displacement. The wind wave component direction of the time window is determined based on the cross-covariance corresponding to the wind wave component, and the swell component direction of the time window is determined based on the cross-covariance corresponding to the swell component. Meanwhile, the main frequency of the wind wave component and the main frequency of the swell component corresponding to the time window can be determined according to the cross-covariance corresponding to the wind wave component and the swell component. For example, the main frequency of the wind wave component and the main frequency of the swell component of the time window can be determined through frequency analysis and spectrum searching.

[0114] In an embodiment, the S103 can be implemented according to the following steps:

[0115] S103-1: Band-pass filtering each three-direction displacement to determine the wind wave component and the swell component of each three-direction displacement.

[0116] For example, for any time window, each three-direction displacement in the time window can be band-pass filtered to obtain the wind wave component and the swell component of each three-direction displacement.

[0117]

[0118] wherein wind represents the wind wave, swell represents the swell, when m is wind, represents the wind wave component of , and when m is swell, represents the swell component of . For the meanings of and when m is wind or swell, refer to the description of , which will not be repeated here.

[0119] S103-2: determining the cross-covariance between each three-direction displacement according to the covariance between the wind wave component and the swell component of each three-direction displacement, and determining the wind wave component direction and the swell component direction of the time window based on the cross-covariance. ​​

[0120] In implementation, for each of the three displacements, the standard deviation of each of the wind wave component and the swell component can be determined according to the following formula

[0121]

[0122] wherein, represents the average value of each vertical displacement in the i-th time window.

[0123] Further, the cross-covariance can be calculated according to the standard deviation. Then, the Fourier coefficients a and b corresponding to the wind wave component and the swell component respectively can be calculated according to the cross-covariance and the following formula

[0124]

[0125] wherein, represents the Fourier coefficient a of the i-th time window in the wind wave component or the swell component. represents the Fourier coefficient b of the i-th time window in the wind wave component or the swell component. represents each wind wave component or each swell component of each vertical displacement in the i-th time window; and the meanings of and can be obtained in the same way. represents the covariance between the horizontal displacement and the vertical displacement in the i-th time window. represents the covariance between the longitudinal displacement and the vertical displacement in the i-th time window. represents the variance of each horizontal displacement in the i-th time window. represents the variance of each longitudinal displacement in the i-th time window. represents the variance of each vertical displacement in the i-th time window. represents the cross-covariance between the horizontal displacement and the vertical displacement of the i-th time window, represents the cross-covariance between the longitudinal displacement and the vertical displacement of the i-th time window. Further, the directions of the wind wave component and the swell component can be determined according to the following formula

[0126]

[0127] wherein, the value of m can be wind or swell, is the direction of the wind wave component of the i-th time window,

[0128] is the direction of the swell component of the i-th time window. ​​​​​​​​​

[0129] S103-3: In the band-pass frequency range, the fast Fourier transform is performed on each vertical displacement wind wave component and swell component to determine the corresponding transform frequency of each vertical displacement wind wave component and swell component.

[0130] In specific implementation, for any time window , the frequency analysis is performed on the corresponding of each vertical displacement in the window , and the dominant frequency is extracted. wherein f represents the dominant frequency, i represents the i-th time window, and m takes the value of wind or swell.

[0131] Specifically, for any time window , the wind wave component and swell component corresponding to each vertical displacement in the window have been separated . The band-pass frequency corresponding to the wind wave component and swell component in the time window is . In the band-pass frequency , the fast Fourier transform is performed on each :

[0132] ;

[0133] wherein f represents the transform frequency corresponding to the vertical displacement wind wave component or swell component; represents the fast Fourier transform.

[0134] S103-4: Based on the frequency peak value, the wind wave component dominant frequency and the swell component dominant frequency of the time window are determined.

[0135] Specifically, in the frequency band range, the frequency peak value corresponding to the wind wave component and the swell component of the time window can be determined according to each change frequency. For example, the frequency peak value can be determined as follows:

[0136]

[0137] wherein when m takes wind or swell, and represent the frequency peak value corresponding to the wind wave component and the frequency peak value corresponding to the swell component of the time window, respectively.

[0138] S103-5: Based on the frequency peak value, the wind wave component dominant frequency and the swell component dominant frequency of the time window are determined.

[0139] In specific implementation, the wind wave component dominant frequency and the swell component dominant frequency of the time window can be determined according to =​ The dominant frequency of the wind and wave components within the time window was calculated. and the dominant frequency of the surge component .

[0140] S104: Based on the projection components of the background flow field on the wind wave component and the swell component, as well as the dominant frequencies of the wind wave component and the swell component, the wave number of the wind wave component and the wave number of the swell component corresponding to the time window are determined by wave number iteration.

[0141] Here, the background current field U refers to the background current field of the sea surface observed by GNSS buoys. This current field can be measured by the current meter at the buoy location or it can be a predicted current field provided by a public platform. The background current field U can include the velocity component of the background current in the horizontal direction. and the corresponding velocity component in the longitudinal direction , that is .

[0142] In practice, the projected component of the background flow field U in the direction of the wind and wave components can be determined. Based on this projected component and the dominant frequency of the wind and wave components, the wave number of the wind and wave components corresponding to the time window can be retrieved through wave number iteration. Simultaneously, the projected component of the background flow field U in the direction of the swell components can be determined. Based on this projected component and the dominant frequency of the swell components, the wave number of the swell components corresponding to the time window can be retrieved through wave number iteration.

[0143] In one embodiment, the above-described S104 can be implemented according to the following steps:

[0144] S104-1: Determine the initial wave number corresponding to the wind wave component and the swell component based on the dominant frequency of the wind wave component and the dominant frequency of the swell component.

[0145] In practical implementation, the calculations above can be used. and the direction of dissemination By combining the background flow field U, the wave frequency is corrected by background flow projection to obtain the natural frequencies corresponding to the wind wave component and the swell component within the time window. .

[0146] During wavenumber inversion, the input data is Direction of dissemination Background flow field Then, the projected components of the background flow field along the wind wave component and the swell component are calculated using the following formulas:

[0147] ;

[0148] in, represents the projection component of the background flow field in the direction of the wind wave component; represents the projection component of the background flow field in the direction of the swell component, which has a bias effect on the observation frequency.

[0149] Further, since the buoy observation obtains the Doppler effect of the relative flow field of the wave The Doppler effect of the wave relative flow field is contained. And the actual wave satisfies the dispersion relation , And The relationship between them is:

[0150] .

[0151] Wherein, represents the wave number corresponding to the wind wave component of the i-th time window; represents the wave number corresponding to the swell component of the i-th time window.

[0152] In order to facilitate the calculation of wave number, the initial wave number corresponding to the wind wave component and the swell component respectively can be calculated by using the deep water approximate wave number calculation method It can be understood that That is, the initial wave number corresponding to the wind wave component or the swell component in the i-th time window .

[0153] S104-2: According to the initial wave number and the projection component of the background flow field in the direction of the wind wave component and the direction of the swell component, the initial frequency corresponding to the wind wave component and the swell component respectively is determined.

[0154] In specific implementation, the initial frequency can be preliminarily corrected according to And To obtain the first round of natural frequency , which is the initial frequency. Specifically, the initial frequency can be calculated according to the following formula:

[0155] .

[0156] S104-3: The initial wave number is taken as the current iteration wave number, and the initial frequency is taken as the current iteration frequency.

[0157] In specific implementation, since the wave number is solved by the wave number iteration method in the present application, after obtaining the initial wave number and the initial frequency, the two can be taken as the current iteration wave number and the current iteration frequency respectively, so as to iterate the new wave number and the frequency.

[0158] S104-4: According to the water depth of the sea surface, the current iteration wave number and the current iteration frequency, the new wave number is determined by using the finite water depth iteration method, and according to the new wave number and the projection component, the new frequency corresponding to the wind wave component and the swell component respectively is determined.

[0159] In implementation, the linear dispersion relation between wave number and natural frequency can be represented by the following formula:

[0160] ;

[0161] where g represents the acceleration of gravity, and h represents the water depth at the buoy.

[0162] Since is a nonlinear term, the inversion of needs to be performed by a numerical method. Therefore, the finite water depth iteration method can be used for the inversion. The finite water depth iteration method can be represented by the following formula:

[0163] ;

[0164] where represents the wave number of the wind wave component or the swell component in the i-th time window obtained in the n-th iteration, and n is taken from 0. When n is taken as 0, , which is the initial wave number, , which is the initial frequency.

[0165] Solve = 0, and perform wave number iteration by the following iteration formula:

[0166] ;

[0167] where represents the wave number of the wind wave component or the swell component in the i-th time window obtained in the n+1-th iteration;

[0168] where .

[0169] In implementation, the current iteration wave number and the current iteration frequency can be used to iterate a new wave number according to the above formula, and a new frequency can be determined according to the new wave number and the formula .

[0170] S104-5: The new wave number and the new frequency are taken as the new current iteration wave number and the new current iteration frequency, and the step of determining the new wave number according to the water depth of the sea surface, the current iteration wave number and the current iteration frequency by the finite water depth iteration method is returned until the wave numbers obtained by the previous and subsequent iterations converge to a preset precision, and the wave number of the wind wave component and the wave number of the swell component corresponding to the time window are obtained.

[0171] In implementation, the new wave number obtained in a certain iteration can be taken as the initial wave number of the next iteration. and the new frequency , as the new current iteration wave number and the new current iteration frequency. Then return to the step of S104-4 until the wave numbers obtained in two adjacent iterations satisfy the following formula, determine the end of iteration, and take the latest obtained wave number as the wave number of the wind wave component and the wave number of the swell component of the time window:

[0172] ;

[0173] wherein, represents the preset precision, which may be 10-6, for example, in the case of being less than the preset precision, it can be considered that the iteration converges, and the latest obtained wave number is taken as the wave number of the time window . And, according to the formula and the solved , the wave number of the time window is determined.

[0174] For example, the initial wave number may be used to calculate the initial frequency ; the and are used to iterate , and the is used to calculate ; then, the and are used to iterate , and the is used to calculate ; the and are used to iterate , and the is used to calculate ; through continuous iteration, the is iterated, and and converge to the preset precision, the is taken as the wave number of the time window , and the natural frequency of the time window is calculated according to . In this way, the wave number is inverted by using the natural frequency and the water depth information through the linear dispersion relation, which can create a basic condition for the second-order stokes shift calculation.

[0175] S105: According to the wave number of the wind wave component, the wave number of the swell component, and the significant wave height, the vertical displacement of the time window is determined, and the height sequence of the sea surface is determined by using the vertical displacement and each vertical displacement.

[0176] Here, the vertical offset of the time window is used to indicate the overall offset of all vertical displacements in the time window in the vertical direction.

[0177] In implementation, according to the wave number and significant wave height corresponding to the wind wave component and the swell component, the offset calculation is performed according to the second-order Stokes theory to obtain two offsets of the SSH of the GNSS buoy on the wind wave component and the swell component. The SSH of the GNSS buoy can be specifically each vertical displacement in the time window in the application. The vertical offset is used to correct each vertical displacement in the time window to obtain the height sequence of the sea surface.

[0178] In an embodiment, the S105 can be implemented according to the following steps:

[0179] S105-1: determining a first offset on the wind wave component according to the wave number of the wind wave component and the significant wave height corresponding to the wind wave component, and determining a second offset on the swell component according to the wave number of the swell component and the significant wave height corresponding to the swell component.

[0180] Here, the first offset is used to indicate the offset of the three-dimensional displacement in the wind wave component direction within the time window; and the second offset is used to indicate the offset of the three-dimensional displacement in the swell component direction within the time window.

[0181] For example, according to the second-order Stokes wave theory, the asymmetry of the wave particle orbit causes the buoy to stay in the wave crest for a slightly longer time than in the wave trough, thereby causing a systematic offset of the vertical mean, and the expression of the offset is as follows:

[0182] ;

[0183] wherein, represents the first offset or the second offset of the i th time window. According to different values of m, it can be or .

[0184] Further, according to the above formula, the offsets corresponding to the wind wave component and the swell component can be calculated according to the following formula:

[0185] ;

[0186] wherein, represents the first offset of the i th time window on the wind wave component; represents the second offset of the i th time window on the swell component; represents the wave number of the wind wave component of the i th time window; represents the wave number of the swell component of the i th time window.

[0187] S105-2: Determine the vertical offset of the time window based on the first offset and the second offset.

[0188] In practice, since wind waves and swells are statistically independent, their deviation from the vertical mean plane is an additive scalar quantity. Therefore, for any given time window... The second-order Stokes offsets (i.e., the first and second offsets) corresponding to the wind and wave components and the swell components respectively within this time window can be scalarly added together to obtain... Total corrected offset (That is, vertical offset). For example, Vertical offset It can be calculated using the following formula:

[0189] .

[0190] S105-3: Using the vertical offset of the time window, the vertical displacements within the time window are offset and corrected to obtain the corrected displacements corresponding to each vertical displacement.

[0191] In practice, the vertical offset corresponding to the time window can be used. Converted to sampling corrections consistent with GNSS buoy observations. The vertical offset of the time window is used as the correction for the vertical displacement at each time point within the time window.

[0192] In practice, the SSH sequence obtained from GNSS buoy observations can be used as follows: , This includes various time windows. Utilizing each time window and its time interval ,right Each of the Perform offset modifications to obtain each The corresponding corrected displacement.

[0193] S105-4: Determine the sea surface height sequence within the time window based on each corrected displacement.

[0194] For example, point-by-point correction is performed on each vertical displacement:

[0195] ;

[0196] in, express The vertical displacement is then corrected in the Euler reference frame. Furthermore, based on the vertical displacement under each time window... According to the height sequence under each time window, a sea surface height time sequence in Euler system is determined.

[0197] In this way, the application uses the wave main frequency observed by the buoy, the background flow horizontal velocity at the buoy point, and the water depth to obtain accurate wave natural frequency and wave number by iteratively solving the equation coupling the Doppler effect and the nonlinear dispersion relationship. By adaptively decomposing the collected sea surface buoy data into multiple independent wave systems such as wind waves and swell waves, and then performing wave number inversion based on each wave system component, the accuracy of the determined wave number and significant wave height can be improved. Using the wave number and significant wave height, the second-order Stokes wave theory formula considering the effect of limited water depth is used to calculate the offset correction of each wave system; and the correction of each wave system is added up to obtain an accurate sea surface height correction value.

[0198] In summary, the application effectively improves the accuracy and stability of GNSS buoy sea surface height measurement by introducing HHT adaptive spectral decomposition, background flow correction, and second-order Stokes offset correction. Compared with traditional methods, the application does not rely on fixed frequency bands or empirical formulas, but automatically identifies the main frequency band of wind waves and swell waves using the HHT method, thereby effectively adapting to changes in different sea conditions and extracting real wave signals. On this basis, by introducing the projection correction of the background flow, the error caused by the Doppler frequency shift is eliminated, and the accuracy of the wave number inversion is improved, especially suitable for nearshore and complex flow field areas. Then, the vertical offset caused by the wave is accurately calculated using the second-order Stokes theory, realizing centimeter-level SSH correction, effectively eliminating the systematic deviation caused by the prolonged residence time of the wave peak, and obtaining Euler system SSH data matching the satellite altimeter and tide gauge, to improve the accuracy and reliability of the calibration and verification.

[0199] To further illustrate the sea surface height measurement correction method provided by the embodiments of the application, the Stokes offset correction effect is verified when the GNSS buoy sea surface height measurement is performed. In a predetermined sea area, the existing GNSS buoy, sea bottom tide gauge, and the sea surface height measurement correction method of the application are used for comparison and verification. The average water depth of the HI in the predetermined sea area is about 30 meters, the observation time lasts for 24 hours, and the sampling rate of the GNSS buoy is 1 Hz. During the observation period, the fixed-point observation data of the sea bottom tide gauge (tide gauge) at the same position as the buoy are synchronously acquired as a comparison reference. The sea surface height measured by the tide gauge is the Euler average height, which is not affected by the Stokes offset, so the true value of the correction effect of the application can be verified. As shown in Figure 2 the time sequence comparison diagram of the original observation value (Lagrangian height, not corrected by the method provided by the application) of a GNSS buoy provided by the embodiments of the application and the reference value (Euler height, unified to the GNSS height system) of the sea bottom tide gauge at the same position. Based onFigure 2 It can be seen that the two groups of sequences are basically consistent in the trend of tidal variation. The GNSS buoy sequence superimposes a large number of high-frequency disturbances, effectively retaining the instantaneous displacement signal of the sea wave, while the tide gauge, due to being placed at a fixed position on the seabed, is weakly affected by the wave, and the sequence is smoother after internal smoothing processing of the instrument. Figure 3 As shown in a difference comparison chart provided by the embodiment of the present application, the chart includes the difference between the original observation value of the GNSS buoy and the reference value of the seabed tide gauge, and the difference between the observation value after offset correction and the reference value of the seabed tide gauge. Figure 3 It can be seen that before the Stokes offset correction is added, the mean value of the original difference sequence is about 2.4 cm (the position of the highest point), indicating that there is a significant systematic deviation, and the standard deviation is 1.8 cm. After the correction amount calculated by the present application is applied, the mean value of the difference is greatly reduced to 0.6 cm, and the systematic deviation is eliminated by 75%; at the same time, the standard deviation of the difference is also reduced to 0.7 cm, and the residual error is reduced by about 61%. Therefore, it can be known that the present application can accurately eliminate the Stokes offset in the Lagrange observation of the GNSS buoy, and significantly improve the precision and stability of the sea surface height measurement.

[0200] Those skilled in the art can understand that in the above method of the specific implementation, the writing order of each step does not mean a strict execution order and does not constitute any limitation on the implementation process, and the specific execution order of each step should be determined by its function and possible internal logic.

[0201] Based on the same inventive concept, the present embodiment also provides a sea surface height measurement correction device corresponding to the sea surface height measurement correction method. Since the principle of solving problems of the device in the present embodiment is similar to the above-mentioned sea surface height measurement correction method of the present embodiment, the implementation of the device can be referred to the implementation of the method, and the repeated parts will not be described.

[0202] As shown in a difference comparison chart provided by the embodiment of the present application, the chart includes the difference between the original observation value of the GNSS buoy and the reference value of the seabed tide gauge, and the difference between the observation value after offset correction and the reference value of the seabed tide gauge. Figure 4 As shown in a schematic diagram of a sea surface height measurement correction device provided by the embodiment of the present application, it comprises:

[0203] The acquisition module 401 is configured to acquire each time window under the sea surface buoy data. Each time window comprises a plurality of continuously collected sea surface buoy data, and the sea surface buoy data comprises three-direction displacement.

[0204] The decomposition module 402 is configured to, for any time window, perform adaptive spectral decomposition and wave system signal extraction on the vertical displacement in each three-direction displacement, to obtain the significant wave height and the band-pass frequency band corresponding to the wind wave component and the swell component of the time window.

[0205] The determining module 403 is configured to determine the wind wave component direction and the swell component direction of the time window according to the cross-covariance between the respective three-way displacements, and determine the wind wave component main frequency and the swell component main frequency of the time window according to the band-pass frequency band.

[0206] The iteration module 404 is configured to determine the wind wave component wave number and the swell component wave number corresponding to the time window by means of wave number iteration according to the projection components of the background flow field of the sea surface in the wind wave component direction and the swell component direction, and the wind wave component main frequency and the swell component main frequency.

[0207] The correction module 405 is configured to determine the vertical offset of the time window according to the wind wave component wave number, the swell component wave number and the significant wave height, and determine the height sequence of the sea surface by using the vertical offset and the respective vertical displacements.

[0208] The description of the processing flow of each module in the device and the interaction flow between the modules can refer to the related description in the above method embodiments, and will not be described in detail here.

[0209] Based on the same technical concept, the embodiments of the present application also provide a computer device. Referring to FIG. 13, which is a structural schematic diagram of a computer device provided by the embodiments of the present application, the computer device comprises: Figure 5

[0210] a processor 501, a memory 502 and a bus 503. The memory 502 stores machine readable instructions executable by the processor 501, and the processor 501 is configured to execute the machine readable instructions stored in the memory 502, and when the machine readable instructions are executed by the processor 501, the processor 501 executes the above steps S101-S105.

[0211] The above memory 502 comprises an internal memory 5021 and an external memory 5022; the internal memory 5021 is also called an internal storage, and is used to temporarily store operation data in the processor 501 and exchange data with the external memory 5022 such as a hard disk, and the processor 501 exchanges data with the external memory 5022 through the internal memory 5021, and when the computer device is running, the processor 501 and the memory 502 communicate through the bus 503, so that the processor 501 executes the instructions mentioned in the above method embodiments.

[0212] The embodiments of the present disclosure also provide a computer readable storage medium, which stores a computer program, and when the computer program is run by a processor, the steps of the sea surface height measurement correction method described in the above method embodiments are executed. The storage medium can be a volatile or non-volatile computer readable storage medium.

[0213] ​The embodiment of the present disclosure further provides a computer program product, which carries program codes. The program codes include instructions for performing the steps of the sea surface height measurement correction method described in the above method embodiments. Details can be referred to the above method embodiments, and will not be repeated here.

[0214] The computer program product can be specifically implemented by hardware, software or a combination thereof. In an optional embodiment, the computer program product is specifically embodied as a computer storage medium. In another optional embodiment, the computer program product is specifically embodied as a software product, such as a software development kit (SDK) and the like.

[0215] Those skilled in the art can clearly understand that, for the convenience and brevity of description, the specific working process of the above-described device can refer to the corresponding process in the above method embodiments, and will not be repeated here. In several embodiments provided by the present disclosure, it should be understood that the disclosed devices and methods can be implemented by other means. The device embodiments described above are only schematic. For example, the division of the units is only a logical function division, and actual implementation can have another division manner. For another example, a plurality of units or components can be combined, or some features can be ignored or not executed. In addition, the coupling or direct coupling or communication connection between the units or components shown or discussed can be indirect coupling or communication connection through some communication interface, device or unit, which can be electrical, mechanical or other forms.

[0216] The units described as separate components can or can not be physically separate, and the components shown as units can or can not be physical units, that is, they can be located in one place or distributed on multiple network units. According to actual needs, some or all of the units can be selected to achieve the purpose of the present embodiment.

[0217] In addition, each functional unit in each embodiment of the present disclosure can be integrated in one processing unit, or each unit can exist physically, or two or more units can be integrated in one unit.

[0218] Finally, it should be noted that the above-described embodiments are merely specific embodiments of the present disclosure, used to illustrate the technical solutions of the present disclosure, and are not intended to limit the present disclosure. The protection scope of the present disclosure is not limited thereto. Although the present disclosure has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that any person skilled in the art can make modifications or easy changes to the technical solutions described in the foregoing embodiments, or easily think of changes or equivalent replacements for some of the technical features; and these modifications, changes or replacements do not cause the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of the present disclosure, and should be covered within the protection scope of the present disclosure. Therefore, the protection scope of the present disclosure should be subject to the protection scope of the claims.

Claims

1. A method for correcting sea surface height measurement, characterized in that, include: Acquire various time windows of sea surface buoy data; each time window includes multiple continuously collected sea surface buoy data, and the sea surface buoy data includes three-dimensional displacement; For any of the aforementioned time windows, adaptive spectral decomposition and wave system signal extraction are performed on the vertical displacement in each of the three-dimensional displacements to obtain the bandpass and significant wave height corresponding to the wind wave component and the swell component, respectively, for the time window. Based on the cross-covariance between the three-dimensional displacements, the wind and wave component directions and the surge component directions of the time window are determined, and based on the bandpass, the dominant frequencies of the wind and wave components and the surge component of the time window are determined. Based on the projection components of the background flow field on the sea surface in the directions of the wind wave component and the swell component, as well as the dominant frequencies of the wind wave component and the swell component, the wave number of the wind wave component and the wave number of the swell component corresponding to the time window are determined by wave number iteration. Based on the wave number of the wind wave component, the wave number of the swell component, and the significant wave height, the vertical offset of the time window is determined, and the height sequence of the sea surface is determined using the vertical offset and each of the vertical displacements.

2. The method according to claim 1, characterized in that, The adaptive spectral decomposition and wave system signal extraction of the vertical displacement in each of the three-dimensional displacements, to obtain the bandpass and significant wave height corresponding to the wind wave component and swell component of the time window, respectively, includes: For any given vertical displacement, empirical mode decomposition is performed on the vertical displacement according to a preset decomposition rule to obtain each intrinsic mode component, and spectral analysis is performed on each intrinsic mode component to obtain the instantaneous amplitude and instantaneous angular frequency corresponding to each intrinsic mode component; Based on the instantaneous amplitude and the instantaneous angular frequency, determine the modal energy and target period of each intrinsic modal component; Based on the modal energy and the target period, the inherent modal components of each vertical displacement are filtered and grouped to obtain the modal groups corresponding to the wind wave component and the swell component of the time window, respectively. Based on the modal groups and the instantaneous angular frequency, the bandpass and component sequence corresponding to the wind wave component and the swell component of the time window are determined, respectively. Based on the component sequences corresponding to the wind and wave components and the swell components, the significant wave heights corresponding to the wind and wave components and the swell components in the time window are determined.

3. The method according to claim 1, characterized in that, The step of determining the wind and wave component directions and the surge component directions of the time window based on the cross-covariance between the three-dimensional displacements, and determining the dominant frequencies of the wind and wave components and the surge component of the time window based on the bandpass, includes: Bandpass filtering is applied to each of the three-dimensional displacements to determine the wind wave component and swell component of each of the three-dimensional displacements. Based on the covariance between the wind and wave components and the surge components of each of the three-dimensional displacements, the cross-covariance between each of the three-dimensional displacements is determined, and based on the cross-covariance, the direction of the wind and wave components and the direction of the surge components of the time window are determined. Within the specified passband, the wind wave component and surge component of each vertical displacement are subjected to Fast Fourier Transform to determine the transform frequency corresponding to each wind wave component and surge component of each vertical displacement. Based on the transformation frequencies corresponding to the wind wave component and swell component of each vertical displacement, the peak frequency corresponding to the time window in the wind wave component and swell component is determined respectively. Based on the frequency peak, the dominant frequencies of the wind wave component and the swell component within the time window are determined.

4. The method according to claim 1, characterized in that, The step of determining the wave number of the wind and wave components and the wave number of the swell component corresponding to the time window through wave number iteration, based on the projection components of the background current field of the sea surface in the directions of the wind and wave components and the dominant frequencies of the wind and wave components and the swell components, includes: Based on the dominant frequencies of the wind wave component and the swell component, determine the initial wave number corresponding to the wind wave component and the swell component respectively; Based on the initial wave number and the projected components of the background flow field in the directions of the wind wave component and the swell component, determine the initial frequencies corresponding to the wind wave component and the swell component, respectively. The initial wavenumber is used as the current iteration wavenumber, and the initial frequency is used as the current iteration frequency; Using the finite depth iterative method, a new wave number is determined based on the sea surface depth, the current iteration wave number, and the current iteration frequency. Then, based on the new wave number and the projection component, new frequencies corresponding to the wind wave component and the swell component are determined respectively. The new wavenumber and the new frequency are used as the new wavenumber and the new frequency of the current iteration, and the process of determining the new wavenumber based on the water depth, the current wavenumber, and the current frequency using the finite depth iteration method is repeated until the wavenumbers obtained from the two iterations converge to a preset accuracy, thus obtaining the wind wave component wavenumber and swell component wavenumber corresponding to the time window.

5. The method according to claim 1, characterized in that, The step of determining the vertical offset of the time window based on the wave number of the wind wave component, the wave number of the swell component, and the significant wave height, and determining the sea surface height sequence using the vertical offset and each of the vertical displacements, includes: Based on the wave number of the wind wave component and the significant wave height corresponding to the wind wave component, a first offset on the wind wave component is determined, and based on the wave number of the swell component and the significant wave height corresponding to the swell component, a second offset on the swell component is determined. The vertical offset of the time window is determined based on the first offset and the second offset. Using the vertical offset of the time window, each vertical displacement within the time window is offset and corrected to obtain the corrected displacement corresponding to each vertical displacement. Based on each of the corrected displacements, the sequence of sea surface heights within the time window is determined.

6. The method according to claim 1, characterized in that, The various time windows for acquiring sea surface buoy data include: Acquire continuous sea surface buoy data output by the buoy in the time domain; According to a preset time length, the data of each sea surface buoy is divided into windows to obtain each initial time window; For any given initial time window, determine the proportion of missing sea surface buoy data within that initial time window; If the missing proportion is less than or equal to a preset proportion, the missing float data within the initial time window is filled in to obtain the time window.

7. The method according to claim 6, characterized in that, The step of completing the missing buoy data within the initial time window to obtain the time window includes: The missing buoy data within the initial time window is filled in to obtain the filled time window; Based on the mean and linear slope of each triaxial displacement in each direction within the completion time window, the mean and linear slope of each triaxial displacement within the completion time window are respectively processed to obtain the processed triaxial displacements. Based on the mean square error of the completion time window in each direction, the three-dimensional displacements to be processed that do not meet the preset error standard are selected from each processed three-dimensional displacement. Using an interpolation method, the optimized displacement corresponding to the three-dimensional displacement to be processed is determined, and the optimized displacement is used to replace the three-dimensional displacement to be processed to obtain the time window.

8. A sea surface height measurement correction device, characterized in that, include: The acquisition module is used to acquire various time windows of sea surface buoy data; each time window includes multiple continuously collected sea surface buoy data, and the sea surface buoy data includes three-dimensional displacement. The decomposition module is used to perform adaptive spectral decomposition and wave system signal extraction on the vertical displacement of each triaxial displacement for any given time window, so as to obtain the bandpass and significant wave height corresponding to the wind wave component and the swell component of the time window, respectively. The determination module is used to determine the wind and wave component direction and the surge component direction of the time window based on the cross-covariance between the three-dimensional displacements, and to determine the dominant frequency of the wind and wave component and the dominant frequency of the surge component of the time window based on the bandpass. The iteration module is used to determine the wave number of the wind wave component and the wave number of the swell component corresponding to the time window by means of wave number iteration, based on the projection components of the background flow field of the sea surface in the directions of the wind wave component and the swell component, as well as the main frequencies of the wind wave component and the swell component. The correction module is used to determine the vertical offset of the time window based on the wave number of the wind wave component, the wave number of the swell component, and the significant wave height, and to determine the sea surface height sequence using the vertical offset and each of the vertical displacements.

9. A computer device, characterized in that, include: The computer device includes a processor, a memory, and a bus. The memory stores machine-readable instructions executable by the processor. When the computer device is running, the processor communicates with the memory via the bus. When the machine-readable instructions are executed by the processor, they perform the steps of the sea surface height measurement correction method as described in any one of claims 1 to 7.

10. A computer program product, comprising a computer program, characterized in that, When the computer program is run by a computer device, the computer device performs the steps of the sea surface height measurement correction method as described in any one of claims 1 to 7.

Images