Wave direction spectrum wavelet transformation method of single-point single-frequency GNSS buoy system

Abstract

The invention belongs to the technical field of GNSS (Global Navigation Satellite System) buoy systems, and particularly relates to a wave direction spectrum wavelet transform method of a single-point single-frequency GNSS buoy system, wavelet transform can provide multi-scale analysis of high time resolution of a high-frequency part and high frequency resolution of a low-frequency part, so that better signal resolution capability can be kept under different frequencies, and the signal resolution ratio of the GNSS buoy system is improved. The method is very suitable for wave analysis under complex sea conditions.

Description

Technical Field

[0001] The present application belongs to the technical field of GNSS buoy systems, and in particular relates to a wavelet transform method for a wave directional spectrum of a single-point single-frequency GNSS buoy system. Background Art

[0002] At present, most commonly used methods for inverting the directional spectrum of waves rely on multi-point measurement arrays, such as multi-sensor (accelerometer, magnetometer, compass combination) array buoys. However, these methods are not only complex and costly, but also require a high sea area during deployment, making it difficult to achieve efficient deployment over a large area. The present invention proposes a method for directly using a single-point GNSS buoy to observe the directional spectrum of waves, and provides detailed settings for the deployment of on-site equipment.

[0003] GNSS buoys are ocean observation equipment that have been used in wave measurement in recent years because of their high-precision positioning and real-time monitoring capabilities. However, most of the existing wave directional spectrum inversion methods based on GNSS buoys use traditional Fourier transform technology. Although these methods can provide wave directional spectrum information with a certain degree of accuracy, they are not good enough in complex sea conditions. Especially in the case of coexistence of multiple wave systems or strong nonlinear wave interactions, the Fourier transform method cannot effectively extract instantaneous wave direction information, resulting in reduced measurement accuracy. Summary of the invention

[0004] Based on the above problems, this application proposes a single-point GNSS buoy wave directional spectrum inversion method based on wavelet transform to improve the real-time inversion capability of wave directional spectrum under complex sea conditions, enhance measurement accuracy and reduce equipment and deployment costs. Its technical solution is: A wavelet transform method for wave directional spectrum of a single-point single-frequency GNSS buoy system comprises the following steps: S1. Establish a carrier phase measurement model; S2. Establish a difference-based observation equation to solve the receiver position change and clock bias, and construct the error equation of carrier phase time difference TDCP; S3. Carrier phase time difference TDCP cycle slip detection and repair; S4. Convert the data coordinate format to obtain the wavefront displacement; S5. Data optimization processing; S6. Wavelet transform of directional spectrum; S7. Adaptive wavelet scaling Calculation of S8. Use kernel density to estimate the directional distribution spectrum; S9. Calculation of wave direction spectrum.

[0005] Preferably, in step S1, the measurement model of the carrier phase is expressed as: ; : The wavelength of the frequency signal, Indicates frequency coding; : Receiver With satellite between Carrier phase observations on frequency; : The geometric distance between the receiver and the satellite; : speed of light; and : The clock difference between the receiver and the satellite; : Integer ambiguity of carrier phase; : ionospheric delay error; : Tropospheric delay error; : Random measurement noise, including multipath effects and other residual errors.

[0006] Preferably, in step S2, in the absence of a cycle slip, the carrier phase observation values ​​of consecutive epochs are differentiated to obtain: ; : The carrier phase change between two epochs; : The geometric distance difference between the receiver and the satellite, indicating the change in the receiver position between two epochs; : The clock error difference between the receiver and the satellite, representing the impact of the clock error change between two epochs on the measurement; and : The changes of ionospheric and tropospheric delay errors, which are the changes of the signal passing through the ionosphere and troposphere between two epochs; : Random measurement error variation between two epochs, including multipath effects and other residual errors.

[0007] Preferably, in step S2, the differential satellite distance equation is as follows: ; is the unit direction vector between the satellite and the receiver, and is calculated as: ; in, and Respectively represent the epoch The position coordinate vectors of the satellite and the receiver at time ; Insert the differential satellite distance equation into the carrier phase differential observation equation to construct the TDCP error equation: ; in, represents the residual error; is the observation vector, calculated as follows: ; The unknowns in this equation include two: : the change in position of the receiver between two epochs; : Change of receiver clock bias; Calculation of receiver speed: By the change in position between two epochs , calculate the instantaneous speed of the receiver : ; in, and is a time point between two epochs.

[0008] Preferably, in step S4, the coordinate conversion step is: S41. Determine the coordinates of the station center: First, according to the latitude and longitude position of the GNSS buoy, determine its origin in the NEU coordinate system; S42. Calculate the transformation matrix: To transform from the geocentric coordinate system to the NEU coordinate system, we first need to construct a rotation matrix from the XYZ coordinate system to the NEU coordinate system; assuming that the geographical location of the GNSS buoy is ,in is the latitude, is the longitude, then the rotation matrix is: ; This matrix can convert XYZ data in the geocentric coordinate system to the NEU coordinate system; S43. Apply the rotation matrix: Use the above rotation matrix to transform the displacement data in the XYZ coordinate system Convert to displacement data in NEU coordinate system : ; The displacement of the GNSS buoy in the east, north and vertical directions is obtained.

[0009] Preferably, in step S5, the data optimization process includes data segmentation, data filtering and data cleaning; Data segmentation: set the data collection interval and the minimum threshold of valid data; Data Filtering: Filter Transfer Function It is expressed as: ; in: is the cutoff frequency, is the signal frequency, is the order of the filter; Data cleaning includes the following steps: S5-31 Calculate the mean and standard deviation: For wavefront displacement data , first calculate its mean and standard deviation : ; ; in, represents the number of samples; S5-32 Eliminate abnormal data beyond five times the standard deviation: The data points were screened by setting a threshold of five times the standard deviation; for each data point , if the following conditions are met: ; Then the data point is considered an outlier and needs to be removed; S5-33 Data update: After eliminating abnormal data, the remaining valid data is retained for subsequent wavefront displacement calculation and analysis.

[0010] Preferably, in step S6, the displacement time series collected by a given GNSS buoy is , wavelet transform is expressed as: ; in: is the time shift parameter, indicating that during the transformation, the sequence The temporal offset of the local features of is the scale parameter, which is related to the frequency of the signal Associated, through the frequency scale conversion formula Make adjustments; It is called the Fourier wavelength factor and its calculation formula is: ; is the dimensionless frequency parameter; Indicates that The complex conjugate of is the mother wavelet, which is used to localize the time and frequency information of the signal. Morlet wavelet is selected for analysis and is defined as: ; The cross wavelet transform is defined as: ; in, and are the wavelet transform results of the east-west and north-west displacement data, yes The complex conjugate of ; The time point after wavelet transformation is recorded as , Use frequency To express it, the local frequency is obtained and time The direction of wave propagation under: ; in, and Real part of the cross wavelet transform of the vertical displacement and the easting and northing displacements.

[0011] Preferably, in step S7, the adaptive wavelet scale The calculation formula is as follows: ; in, It represents the local instantaneous frequency of a signal, describing the main oscillation frequency of the signal at a specific point in time; Signal The analytical signal is: ; in, , represents the imaginary unit, is the Hilbert transform, the instantaneous frequency can be obtained by analyzing the phase of the signal calculate: ; ; Indicates plural The angle between the real axis and the complex plane is the phase angle of the signal; represents the Hilbert transform With the original signal The inverse tangent function of the ratio of .

[0012] Preferably, in step S8, the formula for kernel density estimation is: ; in: Indicates in the direction The estimated density at ; Indicates time The number of local wave direction data points on ; is the bandwidth of the kernel function, which determines the smoothness of the estimate; For in time Local wave direction data observed at is the kernel function: ; ; in, is the standard deviation of the local wave direction data; For each frequency , in the direction Compute the kernel density estimate on: ; All sample points Sum the kernel functions to get the direction The density estimate at , the estimated density function is standardized: ; Thus, the directional distribution function The following conditions are met: .

[0013] Preferably, in step S9, the directional spectrum Describe the waves at different frequencies and direction The energy distribution on the vertical wave spectrum and directional distribution function The product is: ; in Determined by the mean of the wavelet energy over time.

[0014] Compared with the prior art, the present invention has the following beneficial effects: 1. Non-stationary signal processing: Wavelet transform is particularly suitable for processing non-stationary signals, such as transient fluctuations and multi-scale changes commonly seen in ocean waves, and can effectively capture the local time-frequency characteristics of waves.

[0015] 2. Multi-resolution analysis: Wavelet transform can provide multi-scale analysis with high time resolution in the high-frequency part and high frequency resolution in the low-frequency part, which enables it to maintain good signal resolution at different frequencies and is very suitable for wave analysis under complex sea conditions.

[0016] 3. Real-time performance: Wavelet transform shows high flexibility and accuracy in processing real-time signals, and can better reflect the instantaneous state of waves, especially in complex wave conditions such as the coexistence of wind waves and swell waves. BRIEF DESCRIPTION OF THE DRAWINGS

[0017] Figure 1 Diagram of the anchoring scheme for a single-point single-frequency GNSS buoy to observe waves.

[0018] 1-GNSS buoy, 2-floating ball; 3-weight block; 4-rubber rope; 5-PP rope.

[0019] Figure 2 The NEU displacement time series observed by a single-point single-frequency GNSS buoy, (a) represents the north direction, (b) represents the east direction, and (c) represents the vertical direction.

[0020] Figure 3 is the wavelet energy distribution.

[0021] Figure 4 The wavelet transform results of the single-frequency GNSS buoy NEU displacement calculation. (a) represents the three-dimensional wave directional spectrum; (b) represents its directional distribution function. DETAILED DESCRIPTION

[0022] The following will be combined with the embodiments of the present invention and the drawings of the specification to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0023] The present invention proposes a wavelet transform method for inverting the wave directional spectrum based on a single-point single-frequency GNSS buoy, aiming to obtain the wave directional spectrum data through the GNSS buoy, and improve the accuracy and real-time performance of data processing through the wavelet transform technology. The technical solution of the present invention mainly includes the structural design of the buoy system, the data collection method, the specific application method of the wavelet transform and the implementation steps of the data processing.

[0024] 1. Single-point single-frequency GNSS buoy: The core device of the present invention is a single-point single-frequency GNSS buoy, which is used to collect dynamic data of sea surface waves and transmit and process them in real time. The buoy system has multiple integrated functions to ensure its efficient and stable operation in complex sea conditions. Its main components and improved designs are as follows: GNSS module: A high-precision single-frequency GNSS module is used to measure the speed of the buoy in real time, and the wave surface displacement is further derived by analyzing the motion trajectory. The GNSS module supports high-frequency sampling to ensure accurate capture of wave data, and provides a low-power working mode to extend the endurance of the buoy.

[0025] Buoy body: Made of corrosion-resistant, lightweight, high-strength materials, it has good wind and wave resistance and impact resistance, ensuring long-term stable operation in extreme marine environments. The buoy body is designed with a compact internal space structure, integrating GNSS antennas and high-capacity data storage devices, and adopts a multi-layer waterproof sealing structure to ensure that electronic equipment is not corroded by seawater. The buoy is also designed with a buoyancy system with automatic adjustment function, which can automatically adjust the buoy's buoyancy according to sea conditions to ensure the stability of data collection.

[0026] Energy system: A multi-source power supply system is used, including high-efficiency solar panels and large-capacity lithium battery packs, to ensure the buoy's ability to provide power autonomously under long-term unattended conditions. The energy system's management module can dynamically allocate energy use according to environmental conditions and realize intelligent switching between solar energy and lithium batteries. In addition, the system is equipped with a high-performance backup battery to ensure continuous power supply in extreme weather or long periods of no sunlight. To extend the life of the buoy, the system is also designed with a low-power mode to reduce power consumption when no measurements or communications are being performed.

[0027] Computational control system: The buoy's built-in efficient computational control module integrates data acquisition, processing, storage and system monitoring functions. The module can process the data collected by the GNSS module in real time according to the preset sampling frequency and strategy, and make intelligent adjustments based on environmental conditions and system status. The computational control system supports edge computing, which can preliminarily analyze wave data and generate wave surface displacement and velocity information, thereby reducing the computational burden on shore-based data centers.

[0028] Communication module: The buoy integrates a variety of wireless communication methods, which can transmit data to the shore-based data center in real time through satellite communication or 4G / 5G network. The module has an automatic switching function, which can intelligently select the optimal communication method according to the network coverage of the sea area where the buoy is located to ensure real-time transmission and continuity of data. In order to improve communication efficiency, the system is also designed with data compression and package transmission functions to effectively reduce communication costs.

[0029] Environmental adaptability design: The buoy design takes into account the variability of the complex marine environment. It has the ability to withstand strong winds and strong waves, and can work stably under harsh conditions such as high salinity and humidity. The external material of the buoy has self-cleaning properties to prevent the attachment of salt and microorganisms in seawater, reduce maintenance requirements and extend the service life of the buoy.

[0030] 2. Anchoring scheme for single-point single-frequency GNSS buoy wave observation: The present invention proposes an anchoring design suitable for a single-point single-frequency GNSS buoy. Through a reasonable connection method and structural arrangement, the relative motion of the buoy under the action of waves can be effectively measured, while avoiding knotting, twisting or drifting of the buoy system due to currents and waves, thereby improving the accuracy and stability of the buoy observation data.

[0031] (a) Connection between GNSS buoy and buoyancy ball: The GNSS buoy is connected to the buoyant ball by a 20-meter-long rubber elastic rope. The rubber elastic rope has good elasticity and can absorb the up and down movement of the buoy under the action of waves and currents, avoiding excessive or rapid movement of the buoy due to wave impact, thereby maintaining the relative stability of the buoy.

[0032] (b) Design of buoyancy ball: The buoyancy ball has a buoyancy of 30kg and is used to provide additional buoyancy support to ensure that the connection between the buoy and the rubber rope remains stable on the water surface. The buoyancy ball can help reduce the drift amplitude of the buoy system when it is affected by waves and stabilize the position of the buoy on the water surface.

[0033] (c) Connection between the buoyancy ball and the seabed gravity anchor: The buoyancy ball is connected to the gravity anchor on the seabed by a polypropylene rope. The polypropylene rope has good corrosion resistance and strength, and has a certain buoyancy to prevent it from sinking or entangled in the water. To ensure the firmness of the connection, the length of the polypropylene rope is designed to be twice the water depth, which can ensure the free movement range of the buoy when the water depth changes, and at the same time provide sufficient anchoring depth to enhance the stability of the buoy.

[0034] (d) Design of gravity anchor: The gravity anchor has a grip of 100kg and can be effectively fixed to the seabed to prevent the buoy system from drifting or moving due to currents or waves. The design of the gravity anchor ensures that it has good grip performance in various complex seabed environments, thereby providing a stable anchor point for buoy observation.

[0035] (e) Design to prevent rope tangling and twisting: In order to prevent the polypropylene rope from getting knotted or twisted under the action of currents or waves, which would affect the measurement accuracy of the buoy, the present invention designs multiple auxiliary floats and weights to ensure the stability of the rope: A buoy with 5kg buoyancy is tied 5 meters from the seabed. The function of the buoy is to reduce the friction between the rope and the seabed by increasing the buoyancy of a certain section of the polypropylene rope, thereby preventing the rope from getting knotted or entangled in strong currents or waves.

[0036] A 1kg weight is added 5 meters from the water surface. The function of the weight is to provide a certain vertical tension for the buoy system, reduce the floating tendency of the buoy and the polypropylene rope caused by the waves on the water surface, maintain the relative position of the buoy on the water surface, and thus improve the stability of the measurement.

[0037] 3. Data collection: The present invention adopts a single-frequency L1 band GNSS receiving system. The L1 band is the main navigation signal frequency band of the GNSS satellite system, has the ability of global coverage, and can stably obtain the three-dimensional motion data of the buoy under complex sea conditions.

[0038] Sampling frequency: The sampling frequency of the satellite epoch is set to 2 Hz, which means data is collected twice per second. This frequency can capture the dynamic changes of sea surface waves and ensure the continuity and accuracy of the buoy's motion trajectory. The sampling frequency of 2 Hz is reasonable in most ocean wave observation applications and can fully reflect the periodicity and amplitude changes of waves.

[0039] The main contents of the collection are: Carrier Phase Data: GNSS receivers accurately estimate the relative position between the receiver and the satellite by analyzing the phase changes of the carrier signal.

[0040] Broadcast Ephemeris: Broadcast Ephemeris contains data such as satellite orbit, time calibration, signal status, etc. GNSS receivers use this information to correct their own time and position measurements.

[0041] 4. Data processing, single-frequency TDCP (Time-Differenced Carrier Phase): The measurement model of carrier phase can be expressed as the following formula: ; : The wavelength of the frequency signal, Indicates frequency coding; different frequencies of GNSS have different wavelengths, such as and The present invention adopts , : Receiver With satellite between Carrier phase observations on frequency; : The geometric distance between the receiver and the satellite; : speed of light; and : The clock difference between the receiver and the satellite; : The integer ambiguity of the carrier phase; it is an integer representing the uncertainty of the carrier phase during signal propagation; : Ionospheric delay error; an error that affects the propagation of GNSS signals due to refraction of the signal as it passes through the Earth’s ionosphere; : Tropospheric delay error, the delay caused by the signal passing through the troposphere of the Earth's atmosphere; : Random measurement noise, including multipath effects and other residual errors.

[0042] In the absence of cycle slips, the carrier phase observations of consecutive epochs are differentiated to obtain: ; : The carrier phase change between two epochs; : The geometric distance difference between the receiver and the satellite, indicating the change in the receiver position between two epochs; : The clock error difference between the receiver and the satellite, representing the impact of the clock error change between two epochs on the measurement; and : The changes of ionospheric and tropospheric delay errors, which are the changes of the signal passing through the ionosphere and troposphere between two epochs; : Random measurement error variation between two epochs, including multipath effects and other residual errors.

[0043] Further get the satellite distance difference: ; is the unit direction vector between the satellite and the receiver, and is calculated as: ; in, and Respectively represent the epoch The position coordinate vectors of the satellite and receiver at time .

[0044] The position information of the receiver in consecutive epochs can be expressed as ,in is the change in position of the receiver between two epochs. Inserting the differential satellite range equation into the carrier phase differential observation equation yields the error equation for the time differential carrier phase method TDCP: ; in, represents the residual error; is the change in receiver clock bias.

[0045] is the observation vector, calculated as follows: ; The unknowns in this equation include two: : the change in position of the receiver between two epochs; : Change of receiver clock bias; When the number of observed satellites is greater than 4, the least squares method can be used to solve the position change of the receiver and receiver clock offset changes .

[0046] Calculation of receiver speed: By the change in position between two epochs , calculate the instantaneous speed of the receiver : ; in, and is a time point between two epochs.

[0047] When solving the TDCP error equation, the satellite position is calculated using the broadcast ephemeris, while the receiver position is determined using single-point positioning based on pseudorange measurements. These quantities are primarily determined by the direction vector Affects the accuracy of speed measurement. As long as their error is within 10 meters, the impact on speed measurement can be controlled at the order of millimeters per second.

[0048] The satellite's clock drift rate is also calculated from the broadcast ephemeris. Typically, the clock drift rate is less than 0.001 nanoseconds / second, which affects the velocity measurement by less than 0.3 millimeters / second.

[0049] Ionospheric and tropospheric delays generally change slowly. When the sampling frequency reaches 1 Hz or higher, the ionospheric and tropospheric delays between adjacent epochs show strong correlation. Most atmospheric delays can be eliminated by the difference method, and the residual error is treated as noise.

[0050] When constructing the TDCP observation equation, it is assumed that no cycle slips occur between epochs. Therefore, in practical applications, cycle slip detection and repair should be performed. The present invention proposes to use a Doppler-assisted method to detect and repair cycle slips in GNSS data. This method uses Doppler shift measurements, which are less sensitive to cycle slips than phase measurements, thereby more effectively detecting and removing cycle slips.

[0051] 5. Convert the data coordinate format to obtain the wave surface displacement: In the data processing process, the velocity in the geocentric coordinate system (XYZ) is first obtained through the single-point GNSS buoy. By integrating the velocity, the displacement in the XYZ coordinate system can be obtained. However, the calculation of waves requires the displacement in the NEU (North, East, Up) coordinate system, so the XYZ coordinate system needs to be converted to the NEU coordinate system in order to obtain the displacement of the buoy in the east (E), north (N) and vertical (U) directions.

[0052] Coordinate transformation steps: (a) Determine the coordinates of the station center: First, determine the origin of the GNSS buoy in the NEU coordinate system based on its latitude and longitude positions; (b) Calculate the transformation matrix: To transform from the geocentric coordinate system to the NEU coordinate system, we first need to construct a rotation matrix from the XYZ coordinate system to the NEU coordinate system; assuming that the geographic location of the GNSS buoy is ,in is the latitude, is the longitude, then the rotation matrix is: ; This matrix can convert XYZ data in the geocentric coordinate system to the NEU coordinate system; (c) Apply the rotation matrix: Use the above rotation matrix to transform the displacement data in the XYZ coordinate system Convert to displacement data in NEU coordinate system : ; The displacement of the GNSS buoy in the east, north and vertical directions is obtained.

[0053] 6. Data segmentation: In the data processing process of this application, in order to calculate a set of wave data every half an hour, the data segmentation technology is used, and the minimum threshold of valid data is set. If there are less than 2000 data points in a data set, the data set will be regarded as invalid data and excluded.

[0054] Time interval: Split the original time series data into half-hour intervals. For each half-hour of data, form an independent data block, and the timestamp is the starting time of the data segment, which is used for subsequent wave calculation and analysis.

[0055] Check the number of data points: In each half-hour data block, check the number of data points. The sampling frequency is 2Hz, so there should be 3600 data points in half an hour in theory. If the number of data points in a data block is less than 2000, the data block is marked as invalid data and discarded.

[0056] The valid data blocks will be used for subsequent wave calculations and analysis.

[0057] 7. Data filtering: In order to process the long-period error caused by systematic drift in single-frequency GNSS buoy data, the present invention uses a Butterworth high-pass filter to filter the data to ensure the accuracy of the calculation of the wave directional spectrum. The Butterworth filter has a smooth frequency response characteristic, which can maximize the retention of high-frequency components while ensuring smooth transition of the signal, thereby effectively removing systematic drift and low-frequency noise.

[0058] (a) Filter design Filter order: Select a second-order Butterworth filter. The higher the order, the more obvious the filtering effect and the narrower the transition band. For the GNSS buoy data used in the present invention, it is recommended to use a second-order Butterworth high-pass filter to ensure accuracy.

[0059] Cut-off frequency: According to the sampling frequency of the GNSS buoy, the cut-off frequency of the high-pass filter is set to 0.03 Hz, which can effectively filter out the systematic drift without affecting the high-frequency characteristics of the waves.

[0060] Filtering formula Butterworth filter transfer function It is expressed as: ; in: is the cutoff frequency (corresponding to 0.03 Hz), is the signal frequency, is the order of the filter, this application uses .

[0061] 8. Data cleaning: In the data cleaning and quality control process of this application, a gross error elimination method based on five times the standard deviation is used. Considering the large fluctuations of waves in the marine environment, the use of five times the standard deviation can effectively eliminate abnormal data without excessively filtering normal fluctuations. The specific steps are as follows: S5-31 Calculate the mean and standard deviation: For wavefront displacement data , first calculate its mean and standard deviation : ; ; in, represents the number of samples; S5-32 Eliminate abnormal data beyond five times the standard deviation: The data points were screened by setting a threshold of five times the standard deviation; for each data point , if the following conditions are met: ; Then the data point is considered an outlier and needs to be removed; S5-33 Data Update: After removing the abnormal data, the remaining valid data is retained for subsequent wave surface displacement calculation and analysis. This method can filter out obvious abnormal data while retaining the natural wave changes.

[0062] 9. Wavelet transform of directional spectrum: Continuous Wavelet Transform (CWT) is a time-frequency analysis tool that is particularly suitable for non-stationary signals, such as waves in the ocean. Traditional Fourier transform assumes that the signal is stationary when processing it, which is often difficult to meet in complex ocean environments. In contrast, wavelet transform can provide localized information of the signal in both the time domain and the frequency domain, thereby capturing the frequency characteristics of wave signals at different time scales. Methods based on wavelet transform are particularly suitable for complex wave directional spectrum inversion, and can describe the nonlinear characteristics and directional distribution of waves in more detail.

[0063] Single-point GNSS buoys can continuously collect motion data of sea surface waves, but due to the singleness of their observation points, traditional methods are difficult to accurately calculate the directional spectrum. The present invention uses wavelet transform to extract more detailed wave frequency and direction information from single-point data and realize the inversion of directional spectrum.

[0064] Displacement time series collected for a given GNSS buoy (Time series such as Figure 2As shown), the wavelet transform is expressed as: ; in: It's a signal In scale and time shift Wavelet transform of .

[0065] is the time shift parameter, which is the center position of the wavelet function on the time axis. It can be understood as the center point of the "observation window". It controls the position of the wavelet function on the time axis relative to the original signal and represents the local characteristics of different time points in the transformation. The wavelet transform is achieved by translating the wavelet function to different Position, so as to perform local analysis on different parts of the time series, through which the signal details at different time points can be observed.

[0066] is the scale parameter, which is related to the frequency of the signal Associated, through the frequency scale conversion formula Adjust it to give different levels of attention to the overall (low frequency) or local (high frequency) structure; It is called the Fourier wavelength factor and its calculation formula is: ; is a dimensionless frequency parameter, this application ; .

[0067] Indicates that The complex conjugate of It is the mother wavelet, which is used to localize the time and frequency information of the signal. Morlet wavelet is usually selected for analysis and is defined as: ; when When the scale parameter is large, The value of changes slowly, indicating that the wave packet of the wavelet function becomes wider, corresponding to the low-frequency component; When the scale parameter is small, The value of changes rapidly, indicating that the wave packet of the wavelet function becomes narrower, corresponding to the high-frequency component; when When changes, the entire wavelet function moves along the time axis, corresponding to the signal analysis at different time points. It is the integral variable in the convolution operation, which represents the actual time point of the signal. It represents the wavelet function and the signal on the entire time axis. The matching process.

[0068] Wavelet transform is actually a transformation of the wavelet function Move along the time axis (i.e. change ), and then for each The signals near the location are analyzed to extract the local features of the signal at each time point.

[0069] The cross wavelet transform is defined as: ; in, and are the wavelet transform results of the east-west and north-west displacement data, yes The conjugate complex number of ; the cross wavelet transform reveals the common energy and relative phase of two signals in time-frequency space.

[0070] The time point after wavelet transformation is recorded as , Use frequency To express it, the local frequency is obtained and time The direction of wave propagation under: ; in, and Real part of the cross wavelet transform of the vertical displacement and the easting and northing displacements.

[0071] 10. Scale and frequency relation:

[0072] The present invention proposes a signal time-frequency analysis method based on adaptive wavelet transform (AWT), which can dynamically adjust the wavelet scale according to the local characteristics of the signal, thereby realizing accurate time-frequency analysis of complex and non-stationary signals. Compared with the traditional fixed-scale wavelet transform, the method of the present invention can more effectively capture the instantaneous frequency changes of the signal and improve the time-frequency resolution.

[0073] Adaptive wavelet transform introduces a dynamic scale adjustment mechanism to adaptively adjust the wavelet scale according to the local characteristics of the signal. Specifically, at each time point, the local instantaneous frequency of the signal is Compute wavelet scaling. Adaptive wavelet scaling The calculation formula is as follows: ; in, It represents the local instantaneous frequency of a signal, describing the main oscillation frequency of the signal at a specific point in time; Signal The analytical signal is: ; in, , represents the imaginary unit, is the Hilbert transform, the instantaneous frequency can be obtained by analyzing the phase of the signal calculate: ; In the complex plane, Can be used To express: ; Indicates plural The angle between the real axis and the complex plane is the phase angle of the signal; represents the Hilbert transform With the original signal The inverse tangent function of the ratio of .

[0074] 11. Use kernel density to estimate the directional distribution spectrum: Kernel density estimation is a nonparametric method for estimating probability density functions. It estimates the overall density by smoothly distributing each data point to its neighboring area. For the directional distribution function The present invention uses kernel density estimation to estimate the local direction data to be processed.

[0075] The formula for kernel density estimation is: ; in: Indicates in the direction The estimated density at ; Indicates time The number of local wave direction data points on ; is the bandwidth of the kernel function, which determines the smoothness of the estimate; For in time Local wave direction data observed at is the kernel function: ; For each observation point The distribution has a bell-shaped density function. Bandwidth Controls the smoothing of the kernel density estimate. If If it is too small, it will lead to overfitting and the estimation result will be too sharp. If it is too large, it will lead to underfitting and the estimation result will be too smooth. It can be determined by cross-validation or rule law (such as Silverman rule), and the present invention uses the following rule: ; in, is the standard deviation of the local wave direction data; For each frequency , in the direction Compute the kernel density estimate on: ; All sample points Sum the kernel functions to get the direction The density estimate at , the estimated density function is standardized: ; Thus, the directional distribution function The following conditions are met: .

[0076] 12. Wave direction spectrum calculation:

[0077] Direction spectrum Describe the waves at different frequencies and direction The energy distribution on the vertical wave spectrum and directional distribution function The product is: ; in It is determined by the mean value of the wavelet energy over time. The wavelet energy distribution is as follows: Figure 3 As shown in Figure 4 As shown in (a), the direction distribution function is Figure 4 As shown in (b).

[0078] 13. Result transmission: In the present invention, the wave data collected and processed by the single-point single-frequency GNSS buoy is used to calculate the wave directional spectrum through wavelet transform technology. In order to improve the efficiency of data transmission, the buoy only sends the processed key result parameters to the server instead of transmitting the original GNSS data, thereby significantly reducing the transmission bandwidth requirements and costs. The key data parameters transmitted include: Time: Each transmitted data carries an accurate timestamp to mark the directional spectrum result calculated at that moment.

[0079] Latitude and longitude: The real-time geographic location (latitude and longitude) of the buoy, which is used to determine the specific location of the measurement point and the corresponding sea conditions.

[0080] Directional spectrum matrix: The wave directional spectrum matrix calculated using wavelet transform technology contains the wave energy distribution and propagation direction at different frequencies. This is the core information of the directional spectrum and can accurately describe the state of the sea surface waves at that moment.

[0081] Main direction: The main propagation direction of the waves, which is the main energy propagation direction calculated based on the directional spectrum matrix.

[0082] The patented technology was verified on site in the coastal waters of Shandong. Figure 2 The figure shows the NEU displacement time series on July 20, 2024 observed by a single-point single-frequency GNSS buoy, where (a) represents the northward displacement time series, (b) represents the eastward displacement time series, and (c) represents the vertical displacement time series.

[0083] Significantly reduced costs: This invention uses a single-point single-frequency GNSS buoy to observe the three-dimensional directional spectrum of ocean waves, which greatly reduces the complexity and overall cost of the system compared to the traditional buoy matrix system. Accurate directional spectrum information can be obtained by using only a single-point low-cost buoy, avoiding the high deployment and maintenance costs of a multi-buoy system, making large-scale observations more economical and practical.

[0084] The amount of calculation is significantly reduced, and the real-time performance is improved: the present invention uses single-frequency GNSS data for directional spectrum inversion. Compared with dual-frequency or multi-frequency GNSS data, the amount of calculation is greatly reduced, which significantly shortens the data processing time and the overall calculation cycle, thereby improving the real-time performance of the system. The buoy can calculate the directional spectrum of the waves in a short time and quickly transmit the results to the shore-based server to meet the needs of real-time ocean monitoring and early warning.

[0085] Advantages of wavelet transform: Compared with traditional Fourier transform, wavelet transform can analyze in both time and frequency domains, and is particularly suitable for processing non-steady-state processes and short-period wave characteristics in ocean waves. Through wavelet transform, the dynamic changes of waves at different scales can be extracted more accurately, improving the accuracy of directional spectrum analysis in complex ocean conditions.

[0086] Effectively capture multi-scale wave information: Through wavelet transform technology, the present invention can extract energy distribution of different frequency bands from ocean wave data. Especially when multiple wave systems coexist, the propagation direction and energy characteristics of each wave system can be analyzed more accurately, thereby improving the accuracy of directional spectrum calculation.

[0087] GNSS systems do not require regular calibration compared to IMUs: An important advantage of GNSS buoys over systems based on inertial measurement units (IMUs) is that they do not require regular calibration of direction and position. IMU systems will experience cumulative errors and drifts during long-term observations and must be calibrated regularly, while GNSS systems rely on satellite signals for positioning and do not require direction calibration, thereby reducing equipment maintenance requirements and work complexity. In addition, GNSS systems maintain precise positioning and direction tracking for a long time, making them suitable for long-term ocean monitoring and remote deployment, greatly improving system reliability.

[0088] In summary, the present invention has significant technical advantages in reducing costs, reducing calculation complexity, improving real-time performance, enhancing system stability and calibration-free operation, and is particularly suitable for application scenarios such as large-scale ocean observation and real-time early warning.

[0089] The above shows and describes the basic principles, main features and advantages of the present invention. Those skilled in the art should understand that the present invention is not limited by the above embodiments. The above embodiments and descriptions are only preferred examples of the present invention and are not intended to limit the present invention. Without departing from the spirit and scope of the present invention, the present invention may have various changes and improvements, which fall within the scope of the present invention. The scope of protection of the present invention is defined by the attached claims and their equivalents.

Claims

1. A wavelet transform method for the wave directional spectrum of a single-point single-frequency GNSS buoy system, characterized in that: The following steps are involved: S1. Establish a carrier phase measurement model; S2. Establish a difference-based observation equation to solve the receiver position change and clock bias, and construct the error equation of carrier phase time difference TDCP; S3. Carrier phase time difference TDCP cycle slip detection and repair; S4. Convert the data coordinate format to obtain the wavefront displacement; S5. Data optimization processing; S6. Wavelet transform of directional spectrum; S7. Adaptive wavelet scaling Calculation of S8. Use kernel density to estimate the directional distribution spectrum; S9. Calculation of wave direction spectrum.

2. The wavelet transform method for the wave directional spectrum of a single-point single-frequency GNSS buoy system according to claim 1 is characterized in that: In step S1, the measurement model of the carrier phase is expressed as: ； : The wavelength of the frequency signal, Indicates frequency coding; : Receiver With satellite between Carrier phase observations on frequency; : The geometric distance between the receiver and the satellite; : speed of light; and : The clock difference between the receiver and the satellite; : Integer ambiguity of carrier phase; : ionospheric delay error; : Tropospheric delay error; : Random measurement noise.

3. The wavelet transform method for the wave direction spectrum of the single-point single-frequency GNSS buoy system according to claim 2 is characterized in that: In step S2, in the absence of a cycle slip, the carrier phase observations of consecutive epochs are differentiated to obtain: ； : The carrier phase change between two epochs; : The geometric distance difference between the receiver and the satellite, indicating the change in the receiver position between two epochs; : The clock error difference between the receiver and the satellite, representing the impact of the clock error change between two epochs on the measurement; and : The changes of ionospheric and tropospheric delay errors, which are the changes of the signal passing through the ionosphere and troposphere between two epochs; : Random measurement error variation between two epochs.

4. The wavelet transform method for the wave directional spectrum of a single-point single-frequency GNSS buoy system according to claim 2 is characterized in that: In step S2, the differential satellite distance equation is as follows: ； is the unit direction vector between the satellite and the receiver, and is calculated as: ； in, and Respectively represent the epoch The position coordinate vectors of the satellite and the receiver at time ; Insert the differential satellite distance equation into the carrier phase differential observation equation to construct the TDCP error equation: ； in, represents the residual error; is the observation vector, calculated as follows: ； The unknowns in this equation include two: : the change in position of the receiver between two epochs; : Change of receiver clock bias; Calculation of receiver speed: By the change in position between two epochs , calculate the instantaneous speed of the receiver : ； in, and is a time point between two epochs.

5. The wavelet transform method for the wave directional spectrum of a single-point single-frequency GNSS buoy system according to claim 1, characterized in that: In step S4, coordinate conversion step: S41. Determine the coordinates of the station center: First, according to the latitude and longitude position of the GNSS buoy, determine its origin in the NEU coordinate system; S42. Calculate the transformation matrix: To transform from the geocentric coordinate system to the NEU coordinate system, you first need to construct a rotation matrix from the XYZ coordinate system to the NEU coordinate system; Assume that the geographical location of the GNSS buoy is ,in is the latitude, is the longitude, then the rotation matrix is: ； This matrix can convert XYZ data in the geocentric coordinate system to the NEU coordinate system; S43. Apply the rotation matrix: Use the above rotation matrix to transform the displacement data in the XYZ coordinate system Convert to displacement data in NEU coordinate system : ； The displacement of the GNSS buoy in the east, north and vertical directions is obtained.

6. The wavelet transform method for the wave directional spectrum of a single-point single-frequency GNSS buoy system according to claim 1, characterized in that: In step S5, data optimization processing includes data segmentation, data filtering and data cleaning; Data segmentation: set the data collection interval and the minimum threshold of valid data; Data Filtering: Filter Transfer Function It is expressed as: ； in: is the cutoff frequency, is the signal frequency, is the order of the filter; Data cleaning includes the following steps: S5-31 Calculate the mean and standard deviation: For wavefront displacement data , first calculate its mean and standard deviation : ； ； in, represents the number of samples; S5-32 Eliminate abnormal data beyond five times the standard deviation: The data points were screened by setting a threshold of five times the standard deviation; for each data point , if the following conditions are met: ； Then the data point is considered an outlier and needs to be removed; S5-33 Data update: After eliminating abnormal data, the remaining valid data is retained for subsequent wavefront displacement calculation and analysis.

7. The wavelet transform method for the wave directional spectrum of a single-point single-frequency GNSS buoy system according to claim 1, characterized in that: In step S6, the displacement time series collected by a given GNSS buoy is , wavelet transform is expressed as: ； in: is the time shift parameter, indicating that during the transformation, the sequence The temporal offset of the local features of is the scale parameter, which is related to the frequency of the signal Associated, through the frequency scale conversion formula Make adjustments; It is called the Fourier wavelength factor and its calculation formula is: ； is the dimensionless frequency parameter; Indicates that The complex conjugate of is the mother wavelet, which is used to localize the time and frequency information of the signal. Morlet wavelet is selected for analysis and is defined as: ； The cross wavelet transform is defined as: ； in, and are the wavelet transform results of the east-west and north-west displacement data, yes The complex conjugate of ; The time point after wavelet transformation is recorded as , Use frequency To express it, the local frequency is obtained and time The direction of wave propagation under: ； in, and Real part of the cross wavelet transform of the vertical displacement and the easting and northing displacements.

8. The wavelet transform method for the wave directional spectrum of a single-point single-frequency GNSS buoy system according to claim 1, characterized in that: In step S7, the adaptive wavelet scale The calculation formula is as follows: ； in, It represents the local instantaneous frequency of the signal, describing the main oscillation frequency of the signal at a certain point in time; Signal The analytical signal is: ； in, , represents the imaginary unit, is the Hilbert transform, where the instantaneous frequency is obtained by analyzing the phase of the signal calculate: ； ； Indicates plural The angle between the real axis and the complex plane is the phase angle of the signal; represents the Hilbert transform With the original signal The inverse tangent function of the ratio of .

9. The wavelet transform method for the wave direction spectrum of the single-point single-frequency GNSS buoy system according to claim 1, characterized in that: In step S8, the formula for kernel density estimation is: ； in: Indicates in the direction The estimated density at ; Indicates time The number of local wave direction data points on ; is the bandwidth of the kernel function, which determines the smoothness of the estimate; For in time Local wave direction data observed at is the kernel function: ； ； in, is the standard deviation of the local wave direction data; For each frequency , in the direction Compute the kernel density estimate on: ； All sample points Sum the kernel functions to get the direction The density estimate at , the estimated density function is standardized: ； Thus, the directional distribution function The following conditions are met: 。 10. The wavelet transform method for the wave directional spectrum of a single-point single-frequency GNSS buoy system according to claim 1, characterized in that: In step S9, the directional spectrum Describe the waves at different frequencies and direction The energy distribution on the vertical wave spectrum and directional distribution function The product is: ； in Determined by the time mean of the wavelet energy.

Images