Wave spectrum calculation method based on GNSS wave measuring buoy

This application claims priority to Chinese Patent Application No. 202310517281.7, filed on May 10, 2023, which is hereby incorporated by reference in its entirety.

The present disclosure relates to the field of ocean wave directional spectrum calculation technology technologies, and in particular, to a wave spectrum calculation method based on GNSS wave measuring buoys.

At present, basic principles of wave measurement applied on GNSS wave buoys in the market are mainly divided into two types:

1. Position information of the buoy is directly obtained through precise positioning. The mainstream positioning technology includes precision single point positioning, which requires the purchase of high-precision satellites, orbital products, and corresponding commercial software, resulting in high costs and difficulty in large-scale deployment; it can real time dynamic measure, but this technology relies on signal reference stations near the shore, which limits the ability of buoys to conduct offshore surveys.

2. Three-dimensional velocity information of the wave carrier is obtained firstly and then integrate to obtain the displacement of the waves. There are still buoy mounted accelerometers in the market for wave measurement, but this method has a low sampling frequency and uses an instantaneous velocity at that moment when integrating velocity into displacement. There is a large error in displacement calculation, resulting in inaccurate calculation of wave spectrum and wave parameters.

In summary, the application of traditional GNSS measurement methods in ocean buoys wave measurement has many limitations, such as limited observation distance, high cost, cumbersome calculation and average accuracy, and the need for correction services. Therefore, an improvement is needed.

The purpose of the present disclosure is to provide a wave spectrum calculation method based on GNSS wave buoys, which solves the many limitations of traditional GNSS measurement methods applied to ocean buoy wave measurement, such as limited observation distance, high cost, cumbersome calculation and average accuracy, and a need for correction services. The calculation accuracy is high, and no additional correction service is required, it is more suitable for wave measurement. The wave spectrum algorithm in this method can adapt to more complex observation environments, and the results of the velocity measurement algorithm can be brought into a matched wave spectrum calculation method. The algorithm is simple and stable.

In order to achieve the above objectives, the present disclosure provides a wave spectrum calculation method based on GNSS (global navigation satellite system) wave buoys, the variometric differential positioning algorithm is used to calculate a relative displacement of the buoy in a short time through high-frequency sampling results; and the wave spectrum calculation method includes the following steps:

step 1: placing GNSS drift buoys in the sea, collecting buoy position data; step 2: obtaining three-dimensional velocity data of ocean waves with the variometric differential positioning algorithm; step 3: performing high pass filtering to velocity data of a sliding window within a preset time M2 to remove noise signals and trend terms, retaining wave motion data; step 4: calculating cross spectra and wave directional spectra with Longuet-Higgins method.

In an embodiment, in step 2, GNSS drift buoys drift on the ocean at a preset time M1, data measured during a drifting process of the GNSS drift buoys is preprocessed and transmitted to ship or shore base through Beidou or Iridium communication so as to achieve real-time transmission and processing of the data; the three-dimensional velocity data is obtained through the variometric differential positioning algorithm.

In an embodiment, the variometric differential positioning algorithm is a differential observation equation obtained by subtracting standard original current carrying observation equation of adjacent epochs; various epochs dada of multiple satellites are substituted into a variometric differential positioning model, a plurality of observation equations for to-be-solved three-dimensional velocity are obtained with the variometric differential positioning algorithm; the three-dimensional velocity of the buoys is obtained by the observation equation with a least squares method;

the standard original current carrying observation equation is formula (1):

In an embodiment, in step 3, two adjacent epochs (t, t+1) are taken, the two epochs (t, t+1) are respectively substituted into the standard original current carrying observation equation to obtain observation values corresponding to the two epochs; the two equations are subtracted, an ionospheric error term is reduced to second order by applying ionospheric free combination to obtain the following ionospheric time single difference observation equation;

the ionospheric time single difference observation equation is formula (2):

In an embodiment, in the geocentric cartesian coordinate system, Δρ r s (t,t+1) is affected by geometric changes [Δρ r s (t,t+1)] OR caused by satellite orbit motion and Earth rotation, as well as affected by solid tides and ocean loads [Δρ r s (t,t+1)] EtOl , as shown in formula (3):

in the formula, Δε r s (t,t+1) is the noise in a single time difference.

In an embodiment, the GNSS drift buoys generate a 3D direction displacement Δξ r (t,t+1) in unit time; a calculation formula is formula (4) when a high-frequency calculation is used:

Z r s is a zenith angle of the satellite relative to the GNSS drift buoys, and the high frequency >1 Hz.

In an embodiment, a formula (6) can be obtained through formulas (1), (2), (3), (4) and (5)

In an embodiment, in step 4, GNSS satellite antenna receives satellite ephemeris file and transmits the ephemeris file to a data processing module; the data processing module reads the ephemeris file, inputs original data consisting of carrier phase observations, pseudorange, and clock deviation of each epoch in the ephemeris file into a simultaneous equation system, the wave direction spectrum is obtained through the Longuet-Higgins method, and the three-dimensional velocity of the buoys is calculated;

three-dimensional information calculated of the buoys is denoted as ϕη tt , η x , η y , then a cross spectral matrix can be written as formula (7):

w is a frequency of power spectral density at the expected frequency; k is wavenumber at the expected frequency.

In an embodiment, the matrix of formulas (7) and (8) is inverted to obtain the three-dimensional velocity of the buoys, and a calculation formula is formula (9):

In an embodiment, in step 1, the GNSS drift buoys include a housing, the housing is provided with a magnetic switch, the housing is provided with a lithium battery pack, a control motherboard, a GNSS satellite antenna, a satellite communication module, a data processing module, and a storage module, the magnetic switch is electrically connected to the control motherboard; the control motherboard is connected to the GNSS satellite antenna, satellite communication module, data processing module, magnetic switch, and storage module; the lithium battery pack is connected to the control motherboard, GNSS satellite antenna, satellite communication module, data processing module, and storage module, and the storage module can store data.

The present disclosure provides a wave spectrum calculation method based on GNSS wave buoys, which solves many limitations of traditional GNSS measurement methods applied to ocean buoy wave measurement, such as limited observation distance, high cost, cumbersome calculation and average accuracy, and the need for correction services. Compared with other methods for measuring buoy displacement, it has higher theoretical accuracy and does not require additional correction services, making it more suitable for wave measurement. The wave spectrum algorithm in this method can adapt to more complex observation environments and incorporate the results of the velocity measurement algorithm into a matched wave spectrum algorithm. The algorithm is simple and stable.

FIG. 1 is a flowchart of a wave spectrum calculation method based on GNSS wave buoys in the present disclosure.

In order to make the purpose, technical solution, and advantages of the present disclosure clearer and understood, the following will provide further detailed explanations of the present disclosure in combination with embodiments. It should be understood that the specific embodiments described here are only intended to explain the present disclosure and are not intended to limit the present disclosure.

Referring to FIG. 1 , the present disclosure provides a wave spectrum calculation method based on GNSS wave buoys, a relative displacement of the buoy in a short time is calculated with a variometric differential positioning algorithm through high-frequency sampling results.

Specifically, the wave spectrum calculation method includes the following steps:

The GNSS drift buoy includes a housing, on which a magnetic switch is provided. The housing is provided with a lithium battery pack, a control motherboard, a GNSS satellite antenna, a satellite communication module, a data processing module, and a storage module. The magnetic switch is electrically connected to the control motherboard. The control motherboard is connected to the GNSS satellite antenna, satellite communication module, data processing module, magnetic switch, and storage module. The lithium battery pack is connected to the control motherboard, GNSS satellite antenna, satellite communication module, data processing module, and storage module. The storage module is capable of temporarily storing data for a storage period of not less than 30 days.

This GNSS drift buoy has a miniaturized volume and a reasonable center of mass design that does not resonate with a dominant frequency of waves; which can be deployed by one single person, with a length, width, and height of 50 cm×50 cm×50 cm. The buoy has a well-controlled volume, reasonable design, strong wave following ability, and production deployment can save more manpower and material resources.

Step 2: obtaining three-dimensional velocity data of ocean waves with the variometric differential positioning algorithm.

The GNSS drift buoys drift on the ocean at a preset time of M1 (days can be set as needed), data measured during a drifting process of the GNSS drift buoys is preprocessed and transmitted to ship or shore base through Beidou or Iridium communication so as to achieve real-time transmission and processing of the data; the three-dimensional velocity data is obtained through the variometric differential positioning algorithm.

Specifically, during the drifting process of the GNSS drift buoy, GNSS satellite antenna receives raw data from the satellite (namely satellite and orbital products).

The raw data is transferred to a data processing module, which uses the variometric differential positioning algorithm to process the raw data and obtain attitude data of the buoys, which is the three-dimensional velocity data.

Attitude data of the buoys is introduced into the wave spectrum algorithm to obtain ocean environment data of the wave spectrum.

Then, through Beidou or Iridium communication, it is transmitted to ship-based or shore bases so as to achieve real-time data transmission and processing; at the same time, the raw data collected by GNSS drift buoys is stored in a storage module for subsequent other purposes.

The variometric differential positioning algorithm in the present disclosure is different from traditional positioning models of single point positioning and differential positioning. The variometric differential positioning algorithm only requires obtaining standard GPS broadcast products (orbit and ephemeris) and data collected by a separate dual frequency GPS receiver to perform data analysis. The variometric differential positioning algorithm in the present disclosure does not require complex technology or centralized data analysis, and can therefore be embedded in a GPS component.

Data collected from standard GPS broadcast products (orbit and ephemeris) and the separate dual frequency GPS receiver are inputted into the variometric differential positioning algorithm. A specific analysis method of the variometric differential positioning algorithm is as follows:

the variometric differential positioning algorithm is a differential observation equation obtained by subtracting standard original current carrying observation equation of adjacent epochs; various epochs dada of multiple satellites (at least four) are substituted into a variometric differential positioning model. By optimizing error term of the model, a plurality of observation equations for to-be-solved three-dimensional velocity are obtained; the three-dimensional velocity of the buoys is obtained by the observation equation with a least squares method when comminating a plurality of observation equations.

the standard original current carrying observation equation is formula (1):

Step 3: performing high pass filtering to velocity data of a sliding window within a preset time M2 to remove noise signals and trend terms, retaining wave motion data.

The wave spectrum of the sliding window is calculated every 20-30 minutes (i.e. a preset time M2). The sampling frequency is above 2.5 Hz.

Specifically, assuming the influence of cycle jumps is ignored, two adjacent epochs (t, t+1) are taken, and the two epochs (t, t+1) are respectively substituted into the standard original current carrying observation equation to obtain the observation values corresponding to the two epochs. The two equations are subtracted, an ionospheric error term is reduced to second order by applying ionospheric free combination to obtain the following ionospheric time single difference observation equation.

The ionospheric time single difference observation equation is formula (2):

In the geocentric cartesian coordinate system, Δρ r s (t,t+1) is affected by geometric changes [Δρ r s (t,t+1)] OR caused by satellite orbit motion and Earth rotation, as well as solid tides and ocean loads [Δρ r s (t,t+1)] EtOl , as shown in formula (3):

In this coordinate system, the GNSS drift buoys will generate a 3D displacement Δξ r (t,t+1) (up and down, east-west, north-south) in unit time, if high frequency (>1 Hz) is used, the pseudorange variation is almost the same per unit time. According to formula (3), formula (4) can be obtained:

A formula (6) can be obtained through formulas (1), (2), (3), (4) and (5),

ΔΦ r s is a difference of a carrier phase observation value of the receiver relative to the satellite; Δξ γ is a displacement scalar difference in three-dimensional direction of the receiver; Δδt r is a difference of clock error of the receiver; Δp r s is a difference of pseudorange between the satellite and the receiver;

Δδt s is a difference of clock error of the satellite;

Δp r s is difference of errors influence by other factors;

Δm r s is a difference of error caused by multipath effects;

Δε r s is the noise difference in a single time difference.

The present disclosure requires data from at least four satellites in adjacent epochs to solve (least square method) through a system of simultaneous equations (formulas 1 to 6). The more satellite data obtained, the better of the result.

Specifically, the GNSS satellite antenna receives the satellite ephemeris file and transfers it to the data processing module. The data processing module reads the ephemeris file and inputs the original data consisting of carrier phase observations, pseudorange, and clock bias of each epoch in the ephemeris file into a system of simultaneous equations. The wave direction spectrum is obtained through the Longuet-Higgins method, and the three-dimensional velocity of the buoys is calculated;

three-dimensional information calculated of the buoys is denoted as ϕη tt , η x , η y , then a cross spectral matrix can be written as formula (7):

where

All elements in the cross spectral matrix are real numbers, which indicates the frequency of the power spectral density at the desired frequency; and the wavenumber is at the expected frequency.

The matrix is inverted to obtain the wave direction spectrum, and the calculation formula is formula (9):

where

The present disclosure provides a wave spectrum calculation method based on GNSS wave measuring buoy, which solves limitations of traditional GNSS measurement methods applied to ocean buoy wave measurement, such as limited observation distance, high cost, cumbersome calculation and average accuracy, and a need for correction services. Compared with other methods for measuring buoy displacement, it has higher theoretical accuracy and does not require additional correction services, it is more suitable for wave measurement. The wave spectrum algorithm in this method can adapt to more complex observation environments and incorporate results of velocity measurement algorithm into a matched wave spectrum algorithm. The algorithm is simple and stable.

Of course, there may be a plurality of other embodiments of the present disclosure, and those skilled in the art may make various corresponding changes and deformations based on the present disclosure without departing from the spirit and essence of the present disclosure. However, these corresponding changes and deformations should all fall within the protection scope of the claims of the present disclosure.