A wave physical constraint-based buoy wind speed dynamic compensation method and system

Abstract

The present application relates to the technical field of offshore buoy wind speed measurement, and specifically relates to a buoy wind speed dynamic compensation method and system based on wave physical constraints. The method comprises: obtaining inertial measurement unit data, global navigation satellite system data, anemometer data and buoy design parameters of a buoy; calculating the natural frequency and damping ratio of the vertical motion of the buoy according to the buoy design parameters; calculating the wave main frequency and viscous damping coefficient according to the inertial measurement unit data and the global navigation satellite system data; fusing the inertial measurement unit data and the global navigation satellite system data by using a Kalman filtering algorithm according to the natural frequency, the damping ratio, the wave main frequency and the viscous damping coefficient to obtain optimal estimation of the buoy motion state; and calculating the real wind speed based on the optimal estimation of the buoy motion state and the anemometer data. The present application solves the problem that the dynamic compensation of the buoy wind speed is not accurate enough under extreme sea conditions in the prior art.

Description

Technical Field

[0001] This invention relates to the field of marine buoy wind speed measurement technology, specifically a dynamic compensation method and system for buoy wind speed based on wave physics constraints. Background Technology

[0002] Buoys are the core carriers for observing near-surface wind speed at sea. The anemometers on board are inevitably affected by the six-degree-of-freedom motion coupling of the buoy, including the additional velocity error caused by the translation of the buoy and the deviation of the wind speed vector direction caused by the attitude tilt. In order to obtain reliable and true wind speed data, the motion state of the buoy must be accurately calculated and motion compensation must be completed.

[0003] Current mainstream compensation schemes mostly adopt a navigation architecture that combines inertial measurement units (IMUs) and global navigation satellite systems (GNSS) to obtain the buoy's motion state. Under normal sea conditions, this can meet basic compensation requirements. However, it has two shortcomings under extreme sea conditions such as typhoons: First, large buoy swaying can easily cause the GNSS carrier phase to lose lock. After the lock is lost, the position and velocity errors of pure inertial integration accumulate and diverge rapidly over time, and the accuracy of state estimation decays rapidly. Second, under severe impact, small-range high-precision IMUs are prone to range saturation, resulting in distorted output data, which further aggravates motion calculation errors and ultimately leads to a serious decrease in wind speed compensation accuracy.

[0004] Therefore, there is an urgent need for a dynamic compensation method for buoy wind speed that can stably maintain the accuracy of state estimation under extreme sea conditions. Summary of the Invention

[0005] To address the shortcomings of existing technologies, this invention provides a method and system for dynamic compensation of buoy wind speed based on wave physics constraints, solving the problem that dynamic compensation of buoy wind speed is not accurate enough under extreme sea conditions in existing technologies.

[0006] To achieve the above objectives, one aspect of the present invention provides a dynamic compensation method for buoy wind speed based on wave physics constraints. The method includes: acquiring inertial measurement unit (IMU) data, global navigation satellite system (GNSS) data, anemometer data, and buoy design parameters of the buoy; calculating the natural frequency and damping ratio of the buoy's vertical motion based on the buoy design parameters; calculating the wave dominant frequency and viscous damping coefficient based on the IMU data and the GNSS data; fusing the IMU data and the GNSS data using a Kalman filter algorithm based on the natural frequency, the damping ratio, the wave dominant frequency, and the viscous damping coefficient to obtain an optimal estimate of the buoy's motion state; and calculating the actual wind speed based on the optimal estimate of the buoy's motion state and the anemometer data.

[0007] This invention calculates the natural frequency and damping ratio of vertical motion by acquiring buoy design parameters, extracts real-time wave dominant frequency and viscous damping coefficient based on inertial measurement unit (IMU) data and global navigation satellite system (GNSS) data, embeds wave physical constraints into a Kalman filter framework, and utilizes the buoy's own dynamic characteristics to constrain the state recursion process, thus suppressing the divergence of inertial integral error during GNSS lockout. Simultaneously, by updating the dual-weight measurement, it fully leverages the observation information from the IMU and GNSS, completely eliminating the interference of buoy translation and rotation on wind speed measurement, significantly improving the compensation accuracy and data continuity of offshore buoy wind speed observations under severe sea conditions such as typhoons.

[0008] Optionally, acquiring the buoy's inertial measurement unit data, global navigation satellite system data, anemometer data, and buoy design parameters includes: acquiring the buoy's main inertial measurement unit data, auxiliary inertial measurement unit data, global navigation satellite system data, anemometer data, and buoy design parameters; and obtaining inertial measurement unit data based on the main inertial measurement unit data and the auxiliary inertial measurement unit data.

[0009] This invention achieves adaptive weighted fusion of the two sets of inertial data by synchronously acquiring data from primary and secondary inertial measurement units with complementary ranges and accuracies, and then calculating confidence weights based on window statistical variance and signal amplitude. Under normal sea conditions, the high-precision measurement values ​​of the primary IMU are prioritized, while under extreme sea conditions, the primary sensor automatically switches to the secondary IMU with wide-range effective data when the primary sensor becomes saturated and distorted. The entire process is smooth with no data jumps caused by hard switching, avoiding data truncation failure of a single inertial unit under the impact of large waves, thus improving the integrity, continuity, and reliability of inertial measurement data under all sea conditions.

[0010] Optionally, obtaining inertial measurement unit (IMU) data based on the primary IMU data and the secondary IMU data includes: extracting primary acceleration data, primary velocity data, secondary acceleration data, and secondary angular velocity data from the primary IMU data and the secondary IMU data, respectively; performing sliding sampling on the primary acceleration data and the primary velocity data within a predefined time window to obtain multiple time-series window acceleration data and time-series window angular velocity data; calculating the window acceleration variance and the window endpoint acceleration comprehensive amplitude based on the time-series window acceleration data, and calculating the window angular velocity variance and... The window endpoint angular velocity comprehensive amplitude is calculated; the main acceleration confidence weight is calculated based on the window acceleration variance and the window endpoint angular velocity comprehensive amplitude, and the main velocity confidence weight is calculated based on the window angular velocity variance and the window endpoint angular velocity comprehensive amplitude; the main acceleration data and the auxiliary acceleration data are weighted and fused according to the main acceleration confidence weight to obtain fused temporal acceleration data; the main velocity data and the auxiliary angular velocity data are weighted and fused according to the main velocity confidence weight to obtain fused temporal angular velocity data; inertial measurement unit data is formed based on the fused temporal acceleration data and the fused temporal angular velocity data.

[0011] This invention separates the acceleration and angular velocity data from the primary and secondary IMUs, calculates the dynamic confidence weight of the primary sensor using the variance and instantaneous amplitude of statistical data through a sliding window, and synthesizes unified inertial data by adaptively weighting and fusing the acceleration and angular velocity separately. The weights of the primary and secondary sensors are smoothly allocated according to signal fluctuations and saturation levels. Under normal sea conditions, the high-precision advantage of the primary IMU is fully preserved. When large waves cause saturation distortion in the primary sensor, its weight is automatically reduced, and the effective data from the secondary IMU with a wide range is adopted. There are no data jumps caused by hard switching, and inertial measurement noise and truncation distortion are suppressed, which improves the smoothness, reliability and measurement accuracy of the raw inertial data under all sea conditions.

[0012] Optionally, the calculation of the wave dominant frequency and viscous damping coefficient based on the inertial measurement unit data and the global navigation satellite system data includes: extracting vertical acceleration component data from the inertial measurement unit data; calculating the wave dominant frequency based on the vertical acceleration component data; extracting horizontal velocity time-series data and signal quality indication parameters from the global navigation satellite system data; and calculating the viscous damping coefficient using the horizontal velocity time-series data when the signal quality indication parameters meet preset availability conditions.

[0013] This invention solves for the real-time wave dominant frequency by extracting the vertical acceleration after coordinate transformation from fused inertial data, while simultaneously extracting GNSS horizontal velocity and signal quality parameters. It identifies the horizontal viscous damping coefficient online only when the positioning signal meets the standard, and obtains sea state dynamic parameters in real time based on the actual motion response of the buoy. This avoids the defect of mismatch between fixed model parameters and actual sea waves. When the signal deteriorates, it uses historical effective damping values ​​to ensure continuous parameter output, providing dynamic physical constraints for the filtering model and improving the fit of Kalman filter state prediction and the stability of buoy motion estimation under extreme sea states.

[0014] Optionally, the step of fusing the inertial measurement unit data and the global navigation satellite system data using a Kalman filter algorithm based on the natural frequency, the damping ratio, the wave dominant frequency, and the viscous damping coefficient to obtain the optimal estimate of the buoy's motion state includes: setting a state vector of a Kalman filter containing three-dimensional position, three-dimensional velocity, and three-dimensional attitude angles; constructing an attitude transfer matrix using the angular velocity in the inertial measurement unit data as the control input, and constructing a vertical motion discrete-time state transfer matrix and a horizontal motion discrete-time state transfer matrix based on the natural frequency, the damping ratio, and the viscous damping coefficient; and using the vertical motion discrete-time state transfer matrix and the water... The state prior estimate is obtained by performing time prediction on the state vector using the translational motion discrete state transition matrix and the attitude transition matrix, combined with the wave dominant frequency; a first observation matrix is ​​constructed based on the position and velocity in the global navigation satellite system data; the state prior estimate is updated by a first measurement based on the first observation matrix to obtain the state estimate after the first measurement update; a second observation matrix is ​​constructed based on the acceleration in the inertial measurement unit data; the state estimate after the first measurement update is updated by a second measurement based on the second observation matrix to obtain the state estimate after the second measurement update; and the optimal estimate of the buoy's motion state is obtained based on the state estimate after the second measurement update.

[0015] This invention constructs a filtered state vector containing three-dimensional position, velocity, and attitude. It builds attitude, vertical, and horizontal motion transfer matrices based on IMU angular velocity and buoy dynamic parameters, respectively. Combined with wave dominant frequency, it completes state-time prediction with physical constraints. Then, it sequentially constructs GNSS and IMU acceleration observation matrices to perform two-level step-by-step measurement updates. The filtering is constrained by wave dynamics laws. When GNSS is lost, it suppresses the divergence of inertial integral error. It fully utilizes the observation information of the two types of sensors in stages to correct state deviations, thereby improving the accuracy and robustness of buoy three-dimensional motion state estimation under extreme sea conditions.

[0016] Optionally, the step of constructing an attitude transfer matrix using the angular velocity in the inertial measurement unit data as control input, and constructing a vertical motion discrete-time state transfer matrix and a horizontal motion discrete-time state transfer matrix based on the natural frequency, the damping ratio, and the viscous damping coefficient includes: constructing a vertical motion continuous-time system matrix based on the natural frequency and the damping ratio; discretizing the vertical motion continuous-time system matrix using the matrix exponentiation method to obtain the vertical motion discrete-time state transfer matrix; constructing a horizontal motion discrete-time state transfer matrix using the viscous damping coefficient; and constructing an attitude transfer matrix based on the quaternion kinematic equations and using the angular velocity in the inertial measurement unit data as control input.

[0017] This invention constructs a vertically continuous system matrix based on the buoy's natural frequency and damping ratio. It then uses the matrix exponential method to accurately discretize the vertical transfer matrix, ensuring numerical stability. Based on online identification of the viscous damping coefficient, a horizontal transfer matrix is ​​built. Finally, a singular attitude transfer matrix is ​​constructed by combining the IMU angular velocity and quaternion kinematic equations. This fully conforms to the hydrodynamic motion law of the buoy, avoids the recursive errors caused by Euler angle gimbal lock and approximate discretization, and makes the filtered state evolution closely match the real wave response. This improves the accuracy, continuity, and numerical stability of the recursive calculation of buoy attitude, horizontal and vertical motion under extreme sea conditions.

[0018] Optionally, the step of obtaining a prior state estimate by performing time prediction on the state vector using the vertical motion discrete-time state transition matrix, the horizontal motion discrete state transition matrix, and the attitude transition matrix, combined with the wave dominant frequency, includes: recursively calculating the state vector using the vertical motion discrete-time state transition matrix, the horizontal motion discrete state transition matrix, and the attitude transition matrix to obtain a predicted state value; adjusting the Kalman filter process noise covariance matrix according to the wave dominant frequency to obtain an adaptive process noise covariance matrix; recursively calculating the state covariance matrix using the adaptive process noise covariance matrix, the vertical motion discrete-time state transition matrix, the horizontal motion discrete state transition matrix, and the attitude transition matrix to obtain a prior covariance matrix; and forming a prior state estimate based on the predicted state value and the prior covariance matrix.

[0019] This invention recursively obtains state prediction values ​​using three types of transition matrices: vertical, horizontal, and attitude. An adaptive matrix is ​​generated based on the noise covariance matrix of the real-time wave main frequency scaling adjustment process. A prior covariance matrix is ​​obtained by recursively combining the various transition matrices. The state prediction values ​​and the prior covariance matrix together constitute a complete prior estimate. The confidence weight of the filter on the dynamic model is dynamically adjusted according to the real-time wave intensity. In large waves, the process noise is moderately amplified to weaken the model constraints, while in calm sea conditions, the noise is reduced to strengthen the model constraints. This adapts to time-varying wave disturbances and improves the adaptive capability and estimation accuracy of the prior prediction of buoy motion state under different sea conditions.

[0020] Optionally, the step of recursively deriving the state vector using the vertical motion discrete-time state transition matrix, the horizontal motion discrete state transition matrix, and the attitude transition matrix to obtain the state prediction value includes: recursively deriving the vertical position and vertical velocity of the buoy using the vertical motion discrete-time state transition matrix to obtain vertically recursed data; recursively deriving the horizontal position and horizontal velocity of the buoy using the horizontal motion discrete state transition matrix to obtain horizontally recursed data; recursively deriving the attitude angles in the state vector according to the attitude transition matrix to obtain recursed attitude angles; and forming the state prediction value based on the vertically recursed data, the horizontally recursed data, and the recursed attitude angles.

[0021] This invention independently recursively derives the vertical motion data, horizontal motion data, and attitude angles of the buoy using vertical, horizontal, and attitude transfer matrices, and then integrates the three types of recursive results to synthesize a complete state prediction value. It relies on the independent hydrodynamic laws of different degrees of freedom of the buoy to solve the problem in modules, avoiding model interference and computational redundancy caused by overall coupled recursion. The dynamic constraints of each channel do not interfere with each other, reducing the computational burden of filtering numerical calculations and ensuring that the motion recursion of each dimension closely matches the real response of the ocean waves. This improves the computational efficiency and the accuracy of multi-dimensional estimation of the buoy's multi-degree-of-freedom motion state prediction.

[0022] Optionally, adjusting the Kalman filter process noise covariance matrix according to the wave dominant frequency to obtain an adaptive process noise covariance matrix includes: obtaining a preset reference frequency and the Kalman filter process noise covariance matrix; calculating a noise scaling factor based on the relative deviation between the wave dominant frequency and the reference frequency; and adjusting the diagonal elements corresponding to the velocity components in the process noise covariance matrix using the noise scaling factor to obtain the adaptive process noise covariance matrix.

[0023] This invention retrieves the covariance matrix of the reference frequency and the nominal process noise, calculates the noise scaling factor based on the deviation between the wave dominant frequency and the reference frequency, and generates an adaptive covariance matrix by scaling only the diagonal elements of the velocity components within the matrix. It dynamically adjusts the degree of trust of the filter in the velocity prediction model according to the intensity of sea wave vibration. During high-frequency and large waves, it amplifies velocity noise and reduces model weight, while during calm sea surfaces, it reduces noise and strengthens model constraints. It accurately matches the characteristics of real-time sea state disturbances, thereby improving the Kalman filter's adaptability to time-varying sea wave environments and the accuracy of buoy velocity state prediction.

[0024] Another aspect of the present invention provides a buoy wind speed dynamic compensation system based on wave physics constraints, comprising: a processor, an input device, an output device, and a memory, wherein the processor, the input device, the output device, and the memory are interconnected, wherein the memory is used to store a computer program, the computer program including program instructions, and the processor is configured to call the program instructions to execute the buoy wind speed dynamic compensation method based on wave physics constraints as described in any of the preceding aspects of the present invention.

[0025] The present invention provides a buoy wind speed dynamic compensation system based on wave physics constraints, which has a compact structure, stable performance, high integration and simple composition. It can stably execute the buoy wind speed dynamic compensation method based on wave physics constraints provided in the preceding part of the present invention, further improving the overall applicability and practical application capability of the present invention. Attached Figure Description

[0026] Figure 1 This is a flowchart of a buoy wind speed dynamic compensation method based on wave physics constraints according to an embodiment of the present invention. Figure 2 This is a schematic diagram of a buoy wind speed dynamic compensation system based on wave physics constraints, according to an embodiment of the present invention. Detailed Implementation

[0027] Specific embodiments of the present invention will now be described in detail. It should be noted that the embodiments described herein are for illustrative purposes only and are not intended to limit the invention. In the following description, numerous specific details are set forth in order to provide a thorough understanding of the invention. However, it will be apparent to those skilled in the art that these specific details are not necessary to practice the invention. In other instances, well-known circuits, software, or methods have not been specifically described to avoid obscuring the invention.

[0028] Throughout this specification, references to "an embodiment," "an embodiment," "an example," or "an example" mean that a particular feature, structure, or characteristic described in connection with that embodiment or example is included in at least one embodiment of the invention. Therefore, the phrases "in an embodiment," "in an embodiment," "an example," or "an example" appearing in various places throughout the specification do not necessarily refer to the same embodiment or example. Furthermore, specific features, structures, or characteristics can be combined in one or more embodiments or examples in any suitable combination and / or sub-combination. Moreover, those skilled in the art will understand that the illustrations provided herein are for illustrative purposes and are not necessarily drawn to scale.

[0029] Please see Figure 1 In order to solve the problems in the prior art, in an alternative embodiment, such as Figure 1The buoy wind speed dynamic compensation method based on wave physics constraints, as shown, includes the following steps: Step S1: Obtain the buoy's inertial measurement unit data, global navigation satellite system data, anemometer data, and buoy design parameters.

[0030] The acquisition of inertial measurement unit data, global navigation satellite system data, anemometer data, and buoy design parameters specifically includes the following sub-steps: Step S101: Obtain the main inertial measurement unit data, auxiliary inertial measurement unit data, global navigation satellite system data, anemometer data, and buoy design parameters of the buoy.

[0031] In this embodiment, multi-source sensor data are synchronously acquired under the same time reference: the main inertial measurement unit (IMU) outputs data from a small-range, high-precision MEMS inertial measurement element, including triaxial angular velocity and triaxial acceleration (e.g., accelerometer range ±6g, gyroscope range ±500° / s); the auxiliary IMU outputs data from a large-range, low-precision MEMS inertial measurement element, also including triaxial angular velocity and triaxial acceleration (e.g., accelerometer range ±30g, gyroscope range ±2000° / s), with a range greater than that of the main IMU; the entire... The data from the global navigation satellite system is output by a multi-frequency, multi-system high-precision receiver, including three-dimensional position, three-dimensional velocity, and accompanying signal quality indication parameters, including the position accuracy factor (PDOP) and carrier-to-noise ratio (C / N0). The anemometer data is output by a three-dimensional ultrasonic anemometer, including the relative wind speed vector in the sensor coordinate system, the direction of which is defined as the direction of the wind relative to the instrument. The buoy design parameters are provided by the buoy design manual or the factory test report, including the displacement volume, waterline area, total buoy mass, added mass coefficient, and vertical damping coefficient.

[0032] The reason for setting up two inertial measurement units (IMUs) with different ranges is that, under normal sea conditions, the buoy's motion amplitude is relatively small, and the smaller-range primary IMU can provide a more accurate measurement signal. However, under extreme sea conditions such as typhoons, the buoy is subjected to severe wave impacts, and its acceleration and angular velocity may far exceed the upper limit of the normal range. This causes the output signal of the smaller-range primary IMU to be saturated and truncated, resulting in completely distorted measurement data. At this time, although the larger-range secondary IMU has relatively lower accuracy, its wider range allows it to remain within the linear measurement range under the same extreme sea conditions, outputting usable data that has not been truncated. Through the redundant configuration of the two IMUs and subsequent adaptive weighted fusion, the high-precision data of the primary IMU can be prioritized under normal sea conditions, while the weights can be automatically switched to the secondary IMU under extreme sea conditions, thus maintaining the continuous availability of inertial measurement data under all sea conditions. During the data preprocessing stage, the zero-bias correction values ​​and coordinate rotation matrices obtained from the factory calibration are applied to the main inertial measurement unit and the auxiliary inertial measurement unit respectively. The output data of the two are unified into the same coordinate system and the zero bias is subtracted. The zero-bias correction values ​​and coordinate rotation matrices are obtained by joint calibration of the six-degree-of-freedom vibration table before leaving the factory. The zero bias of the accelerometer and the zero bias of the gyroscope of the inertial measurement unit are not included in the subsequent online state estimation in order to reduce the filter dimension and avoid the risk of numerical drift caused by unobservable states.

[0033] Step S102: Obtain inertial measurement unit data based on the main inertial measurement unit data and the auxiliary inertial measurement unit data.

[0034] The process of obtaining inertial measurement unit data based on the main inertial measurement unit data and the auxiliary inertial measurement unit data specifically includes the following sub-steps: Step S10201: Extract the main acceleration data, main velocity data, auxiliary acceleration data, and auxiliary angular velocity data from the main inertial measurement unit data and the auxiliary inertial measurement unit data, respectively.

[0035] In this embodiment, three-axis acceleration vectors and three-axis angular velocity vectors are separated from the data stream output by the primary inertial measurement unit (PMU), and are denoted as primary acceleration data and primary velocity data, respectively. Simultaneously, three-axis acceleration vectors and three-axis angular velocity vectors are separated from the data stream output by the secondary inertial measurement unit (TIMU), and are denoted as secondary acceleration data and secondary angular velocity data, respectively. All extracted acceleration and angular velocity data are defined in the buoy body coordinate system unified by the coordinate rotation matrix. During the extraction process, the timestamp information in the synchronous data stream is used to ensure that the primary acceleration data and secondary acceleration data, and the primary velocity data and secondary angular velocity data at the same moment are strictly aligned in time, providing a data foundation for time registration for subsequent sliding window processing and weighted fusion.

[0036] Step S10202: Perform sliding sampling on the main acceleration data and the main velocity data in a predefined time window to obtain multiple time window acceleration data and time window angular velocity data.

[0037] In this embodiment, a fixed-length sliding time window is first set, with the window length determined by the number of sampling points (e.g., 50 sampling points, corresponding to a 0.5-second time span). Then, with the current moment as the end of the window, sliding sampling is performed on the main acceleration data and the main velocity data respectively: For the main acceleration data, the magnitude of the triaxial acceleration is taken at each sampling moment, and arranged in chronological order to form a time-series window acceleration data sequence; similarly, the magnitude of the triaxial angular velocity is taken at each sampling moment, and arranged in chronological order to form a time-series window angular velocity data sequence. The window slides forward by one sampling step (e.g., one sampling interval) each time. That is, at each new sampling moment, the data of the current moment is added to the end of the window, and the data of the earliest moment at the beginning of the window is removed from the window, always keeping the window containing the latest fixed number of continuous sampling points. This results in multiple time-series window acceleration data and time-series window angular velocity data that are continuously updated over time. These window data retain the local change information of the main acceleration and the main velocity over time, which is used for subsequent calculation of the statistical characteristics within the window. The sliding sampling process is performed only on the data from the master inertial measurement unit, and its window length and sliding step remain constant throughout the fusion process to ensure the consistency of statistical feature calculation.

[0038] Step S10203: Calculate the window acceleration variance and the combined amplitude of the window endpoint acceleration based on the time-series window acceleration data, and calculate the window angular velocity variance and the combined amplitude of the window endpoint angular velocity based on the time-series window angular velocity data.

[0039] In this embodiment, for time-series acceleration data, the arithmetic mean of the acceleration magnitudes of all sampling points within the window is first calculated as the window acceleration mean. Then, the deviations between the acceleration magnitude of each sampling point within the window and the window acceleration mean are squared and averaged to obtain the window acceleration variance. This variance reflects the dispersion of the main acceleration data within the window time range; a larger variance indicates more severe fluctuations in the main acceleration and more significant signal noise. Simultaneously, the instantaneous magnitudes of the three-axis components of the main acceleration at the endpoint (i.e., the current moment) in the time-series acceleration data are taken as the window endpoint acceleration composite amplitude, which reflects the total intensity level of the main acceleration signal at the current moment. Similarly, for time-series window angular velocity data, the same calculation method is used. First, the arithmetic mean of the angular velocity magnitudes of all sampling points within the window is calculated as the mean of the window angular velocity. Then, the variance of the window angular velocity is calculated to characterize the dispersion and noise level of the protagonist velocity data. At the same time, the instantaneous magnitude of the three-axis components of the protagonist velocity at the end of the time-series window angular velocity data is taken as the comprehensive amplitude of the window's end angular velocity, which reflects the overall intensity level of the protagonist velocity signal at the current moment. The above-mentioned window acceleration variance and the comprehensive amplitude of the window's end acceleration are calculated independently, and the window angular velocity variance and the comprehensive amplitude of the window's end angular velocity are also calculated independently, without affecting each other.

[0040] Step S10204: Calculate the confidence weight of the principal acceleration based on the variance of the window acceleration and the comprehensive amplitude of the acceleration at the end of the window, and calculate the confidence weight of the principal velocity based on the variance of the window angular velocity and the comprehensive amplitude of the angular velocity at the end of the window.

[0041] In this embodiment, the confidence weight of the main acceleration is calculated jointly by the variance term and the amplitude term. The variance term reflects the degree of fluctuation of the main acceleration data. The greater the fluctuation, the lower the data confidence. The amplitude term reflects the position of the current main acceleration amplitude relative to the upper limit of the range. The higher the amplitude, the closer the sensor is to the saturation region and the lower the confidence.

[0042] The confidence weights of the principal accelerations satisfy the following formula:

[0043] in, The main acceleration confidence weight of the main inertial measurement unit. Let Variance be the window acceleration. The preset acceleration variance threshold can be determined through Allan variance analysis of the static noise of the primary inertial measurement unit. In this embodiment, the variance value corresponding to 60% of the upper limit of the acceleration range of the primary inertial measurement unit is taken. The total acceleration amplitude at the end of the window. This is an auxiliary function for acceleration amplitude. (Acceleration amplitude auxiliary function) The definition of is:

[0044] in, The upper limit of the acceleration range of the main inertial measurement unit (determined by the factory parameters of the main inertial measurement unit, which is 6g in this embodiment). The physical meaning of this amplitude auxiliary function is as follows: when the main acceleration amplitude does not exceed 70% of the upper limit of the range, the main acceleration data is in a good linear working range, and the amplitude term does not produce a weight reduction; when the amplitude is between 70% and 100% of the upper limit of the range, the sensor gradually enters the nonlinear transition region, and the amplitude term weight gradually decreases from 1 to 0.5 according to a linear law; when the amplitude reaches or exceeds the upper limit of the range, the sensor is completely saturated, the amplitude term weight drops to 0, and the main acceleration confidence weight returns to zero. At this time, the main acceleration data is completely unreliable, and the acceleration data of the auxiliary inertial measurement unit completely takes over.

[0045] The angular velocity confidence weights are calculated independently using the same calculation structure, and the angular velocity confidence weights satisfy the following formula:

[0046] in, The confidence weight of the main velocity of the primary inertial measurement unit. The variance of the window angular velocity. The preset angular velocity variance threshold can be determined through Allan variance analysis of the static noise of the main inertial measurement unit. In this embodiment, the variance value corresponding to 60% of the upper limit of the angular velocity range of the main inertial measurement unit is taken. The combined amplitude of the angular velocity at the end of the window. This is an auxiliary function for the magnitude of angular velocity. The definition of is:

[0047] in, The upper limit of the angular velocity range of the main inertial measurement unit (determined by the factory parameters of the main inertial measurement unit; in this embodiment, it is 500° / s). The confidence weight of the main acceleration and the confidence weight of the main velocity are calculated independently; the saturation state of acceleration only affects... The protagonist's speed saturation state only affects The two do not interfere with each other.

[0048] Step S10205: The main acceleration data and the auxiliary acceleration data are weighted and fused according to the main acceleration confidence weight to obtain fused time-series acceleration data.

[0049] In this embodiment, the main acceleration confidence weight of the main inertial measurement unit is calculated. The confidence weight of the auxiliary acceleration of the auxiliary inertial measurement unit Values That is, the sum of the confidence weights of the two accelerations is always 1. During fusion, the main acceleration data of the main inertial measurement unit is combined with the main acceleration confidence weights. Multiply by combining the auxiliary acceleration data from the auxiliary inertial measurement unit with the auxiliary acceleration confidence weights. The two data points are multiplied and then added together axis by axis to obtain the fused temporal acceleration data. This fusion process is performed independently along the X, Y, and Z axes, but all three axes share the same set of weighting coefficients. and This means that the acceleration data of the three axes at the same time are fused using the same fusion ratio. The fused three-axis acceleration data is the fused temporal acceleration data at the current time. This data retains the high precision characteristics of the primary inertial measurement unit under normal sea conditions and the wide range characteristics of the secondary inertial measurement unit under extreme sea conditions. Moreover, the weights between the primary and secondary units transition continuously as the signal quality of the primary inertial measurement unit changes, without any hard switching.

[0050] Step S10206: The protagonist velocity data and the auxiliary angular velocity data are weighted and fused according to the protagonist velocity confidence weight to obtain fused temporal angular velocity data.

[0051] In this embodiment, the confidence weight of the master inertial measurement unit's main velocity is calculated. The confidence weight of the auxiliary angular velocity of the auxiliary inertial measurement unit Values That is, the sum of the confidence weights of the two angular velocities is always 1. During fusion, the main velocity data of the master inertial measurement unit is combined with the main velocity confidence weights. Multiply by combining the auxiliary angular velocity data of the auxiliary inertial measurement unit with the confidence weight of the auxiliary angular velocity. The two data points are multiplied and then added together axis by axis to obtain the fused temporal angular velocity data. This fusion process is performed independently along the X, Y, and Z axes, but all three axes share the same set of weighting coefficients. and This means that the angular velocity data of the three axes at the same time are fused using the same fusion ratio. The fused three-axis angular velocity data is the fused time-series angular velocity data at the current time. This data retains the high precision characteristics of the primary inertial measurement unit under normal sea conditions and the wide range characteristics of the secondary inertial measurement unit under extreme sea conditions. Moreover, the weights between the primary and secondary inertial measurement units transition continuously with the changes in the signal quality of the primary inertial measurement unit, without any hard switching.

[0052] Step S10207: Inertial measurement unit data is formed based on the fused temporal acceleration data and the fused temporal angular velocity data.

[0053] In this embodiment, the obtained fused temporal acceleration data and fused temporal angular velocity data are aligned and combined according to the same timestamp. That is, the triaxial acceleration data and triaxial angular velocity data at the same moment are merged into a set of inertial measurement unit (IMU) data frames, forming IMU data. This IMU data stream contains triaxial fused acceleration and triaxial fused angular velocity, and is output at the same sampling frequency as the original sensor (100Hz in this embodiment), serving as the sole source of IMU data in all subsequent steps. Through the above processing, the data from the two physical IMUs have been fused into one IMU. This IMU automatically inherits the high-precision characteristics of the primary IMU under normal sea conditions and the wide-range characteristics of the secondary IMU under extreme sea conditions. The transition between the two is continuous and smooth throughout, without any hard switching logic, thus ensuring the continuity and stability of subsequent filtering and data fusion under all sea conditions.

[0054] Step S2: Calculate the natural frequency and damping ratio of the buoy's vertical motion based on the buoy design parameters.

[0055] In this embodiment, the buoy design parameters include drainage volume, waterline area, total buoy mass, added mass coefficient, and vertical damping coefficient: drainage volume characterizes the buoy's static buoyancy; waterline area determines the hydrostatic restoring force per unit vertical displacement, which is the core parameter of vertical motion stiffness; the total buoy mass and added mass together constitute the equivalent inertial mass of vertical motion; the added mass coefficient reflects the added inertial effect generated by the surrounding water when the buoy moves in the water; the vertical damping coefficient reflects the degree of viscous energy dissipation during the buoy's vertical motion, and its value directly affects the calculation accuracy of the damping ratio.

[0056] The formula for calculating the natural frequency is:

[0057] in, The natural frequency of the buoy's vertical motion. The density of seawater (taken in this embodiment) ), Let gravitational acceleration be (take) ), The surface area of ​​the buoy at the waterline. The total mass of the buoy. Add mass to the vertical motion of the buoy. The added mass satisfies the calculation relationship:

[0058] in, The additional mass factor, provided in the buoy design manual, is set to 0.3 in this embodiment, meaning the additional mass is 30% of the total buoy mass. The still water restoring stiffness coefficient satisfies the calculation relationship:

[0059] in, This is the hydrostatic restoring stiffness coefficient, representing the hydrostatic restoring force generated per unit vertical displacement. A higher damping ratio indicates greater vertical stability of the buoy. The formula for calculating the damping ratio is:

[0060] in, The damping ratio for the vertical motion of the buoy. This is the vertical damping coefficient. The static water restoring stiffness coefficient is... The total mass of the buoy. Add mass to the vertical motion. Vertical damping coefficient. This coefficient can be obtained directly from the buoy design manual. If not provided in the manual, it can be obtained through calibration via a buoy free decay test. This coefficient reflects the damping force generated by the viscosity of the surrounding water and wave radiation effects during buoy movement, and is a key parameter connecting the buoy's hydrodynamic characteristics and motion response model. Taking a certain type of 10-meter diameter buoy as an example: the displacement volume is approximately... The waterline surface area is approximately Total mass approximately The added mass coefficient is taken as 0.3, and the vertical damping coefficient is approximately... Substituting into the above formula, the natural frequency is calculated. (Corresponding to a natural period of approximately 5.7 seconds), damping ratio The inherent frequency Damping ratio Once calculated, it remains unchanged throughout the filtering process and serves as a fixed parameter for the vertically damped simple harmonic oscillator model in subsequent step S4.

[0061] Step S3: Calculate the wave dominant frequency and viscous damping coefficient based on the inertial measurement unit data and the global navigation satellite system data.

[0062] The calculation of the wave dominant frequency and viscous damping coefficient based on the inertial measurement unit data and the global navigation satellite system data specifically includes the following sub-steps: Step S301: Extract the vertical acceleration component data from the inertial measurement unit data.

[0063] In this embodiment, three-axis acceleration data is extracted from the inertial measurement unit data, and the roll, pitch, and bow angles from the optimal estimate of the buoy's motion state output in the previous filtering cycle are used to transform the three-axis acceleration from the body coordinate system to the navigation coordinate system. Since the buoy experiences roll and pitch motions on the sea surface, the Z-axis of the body coordinate system does not always coincide with the geographic vertical direction; therefore, the Z-axis acceleration in the body coordinate system cannot be directly used as the vertical acceleration. During the transformation, a rotation matrix from the body coordinate system to the navigation coordinate system is first constructed based on the roll, pitch, and bow angles. This rotation matrix is ​​constructed using the ZYX Euler angle sequence. The bow angle primarily affects the azimuth direction of the horizontal acceleration component and has no significant impact on the amplitude of the vertical acceleration component. Then, the three-axis acceleration data in the body coordinate system is multiplied left by this rotation matrix to obtain the three-axis acceleration data in the navigation coordinate system. The third axis component is then directly extracted from this data as the vertical acceleration component. It should be noted that the optimal attitude estimate from the previous filtering cycle is used for coordinate projection here, rather than the prior prediction from the current cycle, to avoid introducing causal circular dependencies during the filtering recursion process. Through this coordinate projection processing, even if the buoy tilts significantly under extreme sea conditions such as typhoons, the extracted vertical acceleration component data can still accurately reflect the buoy's vertical acceleration, providing reliable input data for subsequent wave frequency extraction.

[0064] Step S302: Calculate the wave main frequency based on the vertical acceleration component data.

[0065] In this embodiment, vertical acceleration component data is used as input, and a recursive autoregressive model with a forgetting factor is employed for real-time spectral estimation to extract the dominant wave frequency of the current sea state. The autoregressive model represents the current vertical acceleration as a linear combination of the vertical accelerations at several past moments plus the forecast error. The model order is determined offline to be 4 using the Akaike Information Criterion to achieve a balance between spectral resolution and computational efficiency.

[0066] A forgetting factor is introduced during the recursive process, assigning weights to historical data in an exponentially decaying manner. This reduces the impact of data further removed from the current time on parameter estimation, enabling the algorithm to track the time-varying characteristics of sea state. The forgetting factor ranges from 0.90 to 0.98; in this embodiment, it is set to 0.95. The coefficients of the autoregressive model are updated in real time using a recursive least squares method with a forgetting factor. A coefficient update is performed every time a new vertical acceleration sampling point is received.

[0067] The power spectral density is calculated based on the updated autoregressive coefficients, and the peak value of the power spectral density curve is searched within the typical effective frequency range of ocean waves (e.g., 0.1~3 rad / s). The frequency corresponding to the peak value is determined as the dominant wave frequency at the current moment. The dominant wave frequency is output once every 1 second for adaptive adjustment of the process noise covariance matrix in subsequent steps. The calculation of power spectral density and dominant frequency values ​​is performed each time based on the updated autoregressive coefficients at the current moment, ensuring that the dominant frequency output can keep up with changes in sea state in a timely manner. Under resonance conditions, the frequency corresponding to the peak value of the response spectrum still represents the dominant wave frequency and will not shift due to the amplification of the frequency response function amplitude. Therefore, no special processing is required for resonance conditions.

[0068] Step S303: Extract horizontal velocity timing data and signal quality indication parameters from the Global Navigation Satellite System data.

[0069] In this embodiment, eastward and northward velocities are extracted from Global Navigation Satellite System (GNSS) data as horizontal velocity time-series data. Simultaneously, signal quality indication parameters are extracted from the accompanying output. These parameters include the Position Precision Factor (PDOP) and the Carrier-to-Noise Ratio (C / N0). The PDOP reflects the quality of the spatial geometry of the satellites currently participating in the positioning process; a smaller value indicates a more uniform satellite distribution and higher positioning accuracy. The C / N0 reflects the signal reception strength of each satellite; a higher value indicates better signal quality. The extracted horizontal velocity time-series data is updated frame-by-frame according to the GNSS receiver's original output frequency, with each data point accompanied by a corresponding timestamp to ensure consistency with the time reference in subsequent steps.

[0070] Step S304: When the signal quality indication parameter meets the preset availability conditions, the viscous damping coefficient is calculated using the horizontal velocity timing data.

[0071] In this embodiment, the system first determines whether the GNSS data at the current moment meets preset availability conditions based on signal quality indication parameters. The preset availability conditions are a position accuracy factor (PDOP) less than 5 and a carrier-to-noise ratio (C / N0) greater than 35 dB-Hz. When both conditions are met simultaneously, the GNSS signal is deemed available, and online identification of the viscous damping coefficient is initiated. During identification, horizontal velocity time-series data is used as input. This data includes eastward and northward velocities, corresponding to the values ​​in the formula... and The identification is based on a pure dynamic model of the buoy's horizontal motion. This model contains only velocity damping and process noise terms, excluding IMU acceleration terms, to avoid mathematical contradictions caused by combining IMU measurements (which already include the actual damping effect) with the damping term to be identified. Specifically, it is expressed as follows:

[0072]

[0073] in, , These are the horizontal accelerations in the east and north directions, respectively. , These are the viscous damping coefficients to be identified in the east and north directions, respectively. , The horizontal velocities are respectively eastward and northward. , Process noise represents the driving effect of unmodeled external forces such as wind, waves, and currents on the horizontal motion of the buoy, as well as residual errors introduced by model simplification. In recursive least squares identification, process noise does not need to be directly measured or pre-assigned, and its statistical characteristics (zero mean, finite variance) are automatically adapted by the algorithm during the recursive process.

[0074] The above continuous-time dynamic model is adapted according to the GNSS sampling period. (In this embodiment, the time is 0.1 seconds) Precise discretization is performed to obtain the discrete recurrence relation:

[0075]

[0076] in, , For discretized recursive coefficients, This represents the discrete time step number. , The eastward and northward velocities are the GNSS measurements at the current moment. , These are the eastward and northward velocities measured by GNSS at the previous moment. , This refers to the process noise after discretization.

[0077] When identifying For input, For the output, a recursive least squares method with a forgetting factor is used for real-time estimation. ; Identification Afterwards, through The viscous damping coefficient is obtained by reverse calculation. ; The identification process and The two processes are completely identical and performed independently. The forgetting factor ranges from 0.95 to 0.99, and in this embodiment, it is set to 0.98. This causes the weight of historical data to decay exponentially with the power of the forgetting factor. The further away the observation data is from the current time, the smaller its contribution to the parameter estimation. This allows the identification results to track the slow time-varying damping coefficient caused by changes in sea state.

[0078] Discretized recursive coefficients during the identification process and Physical boundary constraints are applied, forcibly clamping the value to a preset minimum and maximum value. In this embodiment, the clamping interval is taken as... ,in The preset maximum allowable damping coefficient (0.5s in this embodiment) is used. -1 The lower limit is approximately 0.951; clamping to this interval physically means ensuring the damping coefficient obtained through reverse calculation. , The values ​​are always non-negative and do not exceed a preset upper limit to ensure that the horizontal motion state transition matrix is ​​always a stable matrix, avoiding divergence in identification results caused by GNSS velocity measurement noise. The identified viscous damping coefficients are stored in a cache according to timestamps. When the GNSS signal quality deteriorates and no longer meets the usability conditions, the latest effective damping coefficient in the cache is used as the current damping coefficient value for the construction of subsequent horizontal motion discrete state transition matrices.

[0079] Step S4: Based on the natural frequency, the damping ratio, the wave dominant frequency, and the viscous damping coefficient, the Kalman filter algorithm is used to fuse the inertial measurement unit data and the global navigation satellite system data to obtain the optimal estimate of the buoy's motion state.

[0080] The process of fusing the inertial measurement unit data and the global navigation satellite system data using a Kalman filter algorithm to obtain the optimal estimate of the buoy's motion state, based on the natural frequency, the damping ratio, the wave dominant frequency, and the viscous damping coefficient, specifically includes the following sub-steps: Step S401: Set the state vector of the Kalman filter, which includes three-dimensional position, three-dimensional velocity, and three-dimensional attitude angle.

[0081] In this embodiment, the state vector of the Kalman filter is set to 9 dimensions, specifically including three-dimensional position, three-dimensional velocity, and three-dimensional attitude angle, as shown below:

[0082] in, For state vectors, , , These represent the positions in the east, north, and sky directions within the navigation coordinate system. , , These represent the speeds in the east, north, and sky directions, respectively, within the navigation coordinate system. , , These are the roll angle, pitch angle, and bow angle, respectively, indicated by superscript. This represents the vector transpose. The position and velocity in the state vector are both defined in the Northeast-Eastern-Upper-Southern (ENU) navigation coordinate system, i.e. The axis points due east. The axis points due north. The axis points to the zenith (positive); attitude angles are adopted in the ZYX Euler angle sequence (first rotating around the X-axis). Then rotate around the Y-axis Finally, rotate around the Z-axis The state vector is defined to describe the rotational relationship between the buoy's body coordinate system and the northeast-to-northeast navigation coordinate system. The accelerometer and gyroscope zero biases of the IMU have been deducted during the data preprocessing stage through factory calibration parameters, and therefore are not included in the online estimation of the state vector. This reduces the filter dimension, computational burden, and avoids the risk of numerical drift caused by the unobservability between zero bias and attitude.

[0083] Step S402: Construct an attitude transfer matrix using the angular velocity in the inertial measurement unit data as the control input, and construct a vertical motion discrete-time state transfer matrix and a horizontal motion discrete-time state transfer matrix based on the natural frequency, the damping ratio, and the viscous damping coefficient.

[0084] The process of constructing an attitude transfer matrix using the angular velocity from the inertial measurement unit data as control input, and constructing a vertical motion discrete-time state transfer matrix and a horizontal motion discrete-time state transfer matrix based on the natural frequency, the damping ratio, and the viscous damping coefficient, specifically includes the following sub-steps: Step S40201: Construct a vertical motion continuous-time system matrix based on the natural frequency and the damping ratio.

[0085] In this embodiment, based on the inherent frequency Damping ratio A continuous-time system matrix for vertical motion is constructed. The vertical motion of the buoy can be modeled as damped simple harmonic motion, and its equations of motion describe the dynamic behavior of the buoy under the combined action of hydrostatic restoring force (proportional to the vertical displacement and pointing towards the static equilibrium position) and damping force (proportional to the vertical velocity and opposite in direction to the velocity). The vertical position is taken as... and vertical velocity The state variable is the vertical motion state vector. Then the continuous-time state equation for vertical motion Can be written as ,in Vertical process noise represents external disturbances such as unmodeled wave forcing.

[0086] Vertical motion continuous-time system matrix Satisfy the following formula:

[0087] In this matrix, the first row represents the vertical velocity as the rate of change of the vertical position, i.e. ; in the second row of the matrix The acceleration term corresponding to the restoring force in still water (i.e.) ), The acceleration term corresponding to the damping force (i.e.) This matrix fully describes the inherent dynamic characteristics of the buoy's vertical motion.

[0088] Step S40202: The vertical motion continuous-time system matrix is ​​discretized using the matrix exponentiation method to obtain the vertical motion discrete-time state transition matrix.

[0089] For vertical motion continuous-time system matrix The vertical motion discrete-time state transition matrix is ​​obtained by using the matrix exponentiation method for precise discretization:

[0090] in, Let be the sampling period of the Kalman filter (0.01 seconds in this embodiment), and let the damped natural frequency satisfy . , For a damped natural frequency, The natural frequency of the buoy's vertical motion. Let be the damping ratio of the buoy's vertical motion. This calculation precisely maps the continuous-time system matrix to the discrete-time state transition matrix. This ensures that the state evolution of the system between sampling times conforms to the analytical solution of damped simple harmonic motion.

[0091] Compared to the Euler approximation discretization, the matrix exponential method offers the advantage of unconditional numerical stability, preventing numerical divergence even with large sampling steps. This discrete-time state transition matrix contains no IMU-measured acceleration terms and consists solely of the buoy's natural frequency. Damping ratio and sampling period The decision was made to fully preserve the physical characteristics of the damped simple harmonic oscillator model. Due to the use of precise discretization, the recursive relationship remains numerically stable under any parameter conditions. Even in the extreme case where the GNSS is completely unlocked and the recursion relies solely on the vertical physical model, the vertical position and vertical velocity can still be constrained within a physically reasonable range, avoiding the divergence problem caused by the secondary accumulation of errors over time due to pure IMU integration.

[0092] Step S40203: Construct the discrete state transition matrix for horizontal motion using the viscous damping coefficient.

[0093] Based on the viscous damping coefficient and Construct discrete state transition matrices for horizontal motion in the east and north directions, respectively, where The viscous damping coefficient is the eastward viscous damping coefficient. This is the northward viscous damping coefficient, which is the latest valid value in the buffer when the GNSS signal quality deteriorates and no longer meets the usability conditions.

[0094] Taking the x-direction as an example, the continuous-time dynamic model of horizontal motion with viscous damping is as follows: The relationship between its position and velocity satisfies The continuous-time model is precisely discretized using the matrix exponentiation method consistent with the vertical direction, yielding the discrete state transition matrix for eastward horizontal motion. The velocity components in this matrix decay exponentially, and the decay rate is determined by the viscous damping coefficient. Control; position component is velocity during the sampling period The integral accumulation within the period is determined by the damping decay rate and the sampling period. As it approaches 0 (i.e., the undamped limit), the matrix degenerates into a purely kinematic recursive form.

[0095] Similarly, construct the discrete state transition matrix for northward horizontal motion. Its form is similar to Completely identical, only Replace with The discrete state transition matrix for horizontal motion does not contain any IMU horizontal acceleration terms and is driven entirely by a dynamic model with viscous damping. The damping force is a natural recursion of the horizontal velocity, and there is no redundant calculation.

[0096] Step S40204: Based on the quaternion kinematic equations, construct an attitude transfer matrix using the angular velocity in the inertial measurement unit data as the control input.

[0097] In this embodiment, the attitude transfer matrix is ​​constructed based on quaternion kinematic equations, using the angular velocity components in the inertial measurement unit data. As a control input, , , These represent the angular velocity components along the X, Y, and Z axes of the inertial measurement unit (IMU) data. Quaternions are used for attitude recursion because they eliminate the singularity problem of Euler angles at large angles (e.g., pitch angles approaching ±90°), maintaining numerical stability in attitude recursion even under extreme sea conditions with significant buoy swaying. Attitude quaternions are defined as follows:

[0098] in, The real part of the quaternion. , , For the imaginary part of the quaternion, The quaternion represents the vector transpose and satisfies the normalization constraint. .

[0099] The kinematic equations for quaternions are:

[0100] in, The time derivative of the attitude quaternion. This is the quaternion multiplication operator. This is the quaternion representation of the angular velocity of the inertial measurement unit, i.e., a quaternion with a real part of 0 and an imaginary part consisting of the three angular velocity components. The continuous-time quaternion kinematic equations are then precisely discretized using the matrix exponentiation method to obtain the discrete-time state transition matrix of the attitude quaternion. The expression is:

[0101] in, The attitude quaternion discrete-time state transition matrix is... The sampling period of the Kalman filter. angular velocity The corresponding 4×4 antisymmetric matrix:

[0102] because As an antisymmetric matrix, its matrix exponent can be accurately calculated using a closed-form analytical method, avoiding the truncation error inherent in numerical matrix exponent calculations. After discretization, the state recurrence relation of the attitude quaternions is as follows: ,in Let be the attitude quaternion at the current moment. Let be the attitude quaternion from the previous moment. This represents the discrete-time step number. It should be noted that the above quaternion discrete state transition matrix... Numerical recursion is used only for attitude angles to avoid the singularity problem of Euler angles at large angles. In the covariance recursion, the attitude components are linearized and approximated using the first-order Jacobian matrix corresponding to the kinematic equation of the Euler angles to match the dimensional definition of the 9-dimensional state vector. This processing method conforms to the standard framework of extended Kalman filtering.

[0103] After each step of the recursion is completed, the quaternion needs to be... Convert to ZYX Euler angles (roll angles) Pitch angle Bow roll angle The transformation relationship is as follows:

[0104]

[0105]

[0106] in, It is the arctangent function in the fourth quadrant. It is an arcsine function.

[0107] The obtained Euler angles are used to update the attitude components in the Kalman filter state vector, and to construct the coordinate system rotation matrix required for acceleration coordinate transformation.

[0108] This invention employs an open-loop attitude recursive architecture driven solely by the three-axis angular velocity of the inertial measurement unit, adapting to marine conditions characterized by high-frequency, large-amplitude buoy oscillations. Buoy attitude exhibits high-frequency oscillation characteristics, and GNSS and accelerometer measurements cannot effectively constrain high-frequency attitude. Furthermore, wind speed compensation calculation relies solely on instantaneous attitude angles and does not require long-term attitude accumulation. This open-loop architecture eliminates the need for real-time correction from external observations, enabling stable output of instantaneous high-precision attitude and improving the real-time performance and engineering adaptability of wind speed compensation calculations under extreme sea conditions.

[0109] It should be noted that the above quaternion discrete state transition matrix Numerical recursion is used solely for attitude angles to avoid the singularity problem of Euler angles at large angles. In the covariance recursion, the attitude components are linearized using the first-order Jacobian matrix corresponding to the Euler angle kinematic equations to match the dimensional definition of the 9-dimensional state vector. In this invention, quaternions are only used as numerical integration tools for the attitude kinematic equations. The attitude components in the state vector are always represented by Euler angles. The covariance recursion is linearized based on the state quantities (Euler angles). In engineering practice, both are equivalent to converting the quaternion recursion results to Euler angle space and then performing linearization analysis, which conforms to the standard framework of extended Kalman filtering. The accuracy of this open-loop attitude recursion scheme depends on the zero-bias stability and noise level of the inertial measurement unit's angular velocity. For short- to medium-term marine wind speed measurement tasks (within several hours), the angular velocity integral attitude error of a high-precision inertial measurement unit is within an acceptable range, sufficient to support the accuracy requirements of wind speed compensation. For applications requiring long-term attitude maintenance, magnetometer or differential global navigation satellite system heading observations can be further introduced as supplementary measurements. Such extensions do not constitute a substantial change to this invention. The core innovation of this invention lies in embedding wave physical constraints into linear motion state prediction. Since the wind speed compensation accuracy is less sensitive to attitude error than to translational velocity error, this open-loop attitude recursion architecture is reasonable and feasible in engineering.

[0110] Step S403: Using the vertical motion discrete-time state transition matrix, the horizontal motion discrete state transition matrix, and the attitude transition matrix, and in conjunction with the wave dominant frequency, time prediction is performed on the state vector to obtain a priori state estimate.

[0111] Specifically, the process of obtaining a priori state estimation by performing time prediction on the state vector using the vertical motion discrete-time state transition matrix, the horizontal motion discrete state transition matrix, and the attitude transition matrix, combined with the wave dominant frequency, includes the following sub-steps: Step S40301: The state vector is recursively derived using the vertical motion discrete-time state transition matrix, the horizontal motion discrete state transition matrix, and the attitude transition matrix to obtain the state prediction value.

[0112] The process of recursively deriving the state vector from the state vector using the vertical motion discrete-time state transition matrix, the horizontal motion discrete state transition matrix, and the attitude transition matrix to obtain the state prediction value specifically includes the following sub-steps: Step S4030101: The vertical position and vertical velocity of the buoy are recursively derived using the vertical motion discrete-time state transition matrix to obtain the vertical recursive data.

[0113] Extract the vertical position from the optimal estimate (posterior estimate) of the buoy's motion state at the previous moment. and vertical velocity As the recursive input, where This represents the discrete-time step number. The discrete-time state transition matrix for vertical motion is used. By performing time prediction recursion, we obtain the prior predicted values ​​of the vertical position and vertical velocity at the current moment, which are the data after vertical recursion.

[0114] This recursive process does not include any IMU acceleration measurements; it is entirely driven by a damped simple harmonic oscillator physical model, with vertical position and vertical velocity determined by the buoy's natural frequency. Damping ratio The determined dynamic laws evolve naturally. The vertical recursive data obtained will serve as the prior estimate of the vertical component at the current moment, participating in the assembly of the subsequent complete state prediction vector.

[0115] Step S4030102: Use the horizontal motion discrete state transition matrix to recursively deduce the horizontal position and horizontal velocity of the buoy to obtain the horizontal recursive data.

[0116] Extract the eastward position from the optimal estimate (posterior estimate) of the buoy's motion state at the previous moment. Eastward speed Northward position and northbound speed As the recursive input, where This represents the discrete time step number. The discrete state transition matrix for eastward horizontal motion is used. Discrete state transition matrix for northward horizontal motion The horizontal position and horizontal velocity in the east and north directions are predicted and recursively calculated over time to obtain the prior predicted values ​​of the position and velocity in the east and north directions at the current time, which are the data after horizontal recursion.

[0117] This recursive process does not include any IMU horizontal acceleration measurements; it is entirely driven by a viscous damped dynamic model. The eastward and northward directions are recursively calculated independently, with horizontal position and velocity determined by the viscous damping coefficient. and The determined exponential decay law evolves naturally. The horizontal recursive data obtained will serve as the prior estimate of the horizontal component at the current moment, participating in the assembly of the subsequent complete state prediction vector.

[0118] Step S4030103: The attitude angles in the state vector are recursively derived based on the attitude transition matrix to obtain the recursively derived attitude angles.

[0119] Extract the roll angle from the optimal estimate (posterior estimate) of the buoy's motion state at the previous moment. Pitch angle and bow angle As the recursive input, where, This represents the discrete-time step number. Since the attitude recursion is performed internally in quaternion form, the Euler angles from the previous time step are first converted into attitude quaternions. Using attitude quaternions discrete-time state transition matrix Perform time prediction recursion on attitude quaternions, and after the recursion is completed, output the predicted quaternions. Convert to ZYX Euler angles to obtain the roll angle at the current moment. Pitch angle and bow angle The prior predicted value, i.e. the attitude angle after recursion.

[0120] Step S4030104: Based on the vertical recursive data, the horizontal recursive data, and the recursive attitude angle, a state prediction value is formed.

[0121] In this embodiment, the recursively obtained prior predictions of vertical position and vertical velocity, eastward and northward position and velocity, and roll, pitch, and yaw angles are assembled into a complete 9-dimensional state prediction vector according to the definition order of the Kalman filter state vector. This state prediction vector is the prior prediction value of the current moment. The linear motion part (position and velocity) of this state prediction vector is entirely driven by the damped simple harmonic oscillator model and the viscous damped dynamics model, while the attitude angle part is driven by the angular velocity of the inertial measurement unit (IMU), without introducing IMU acceleration as a control input. After assembly, this state prediction vector will serve as the prior estimate of the current moment, forming a complete prior estimate pair with the prior covariance matrix obtained in subsequent steps, for use in subsequent measurement update stages.

[0122] Step S40302: Adjust the process noise covariance matrix of the Kalman filter according to the wave dominant frequency to obtain an adaptive process noise covariance matrix.

[0123] The process noise covariance matrix of the Kalman filter is adjusted according to the dominant wave frequency to obtain the adaptive process noise covariance matrix, which specifically includes the following sub-steps: Step S4030201: Obtain the preset reference frequency and the process noise covariance matrix of the Kalman filter.

[0124] In this embodiment, the preset reference frequency is an empirically calibrated value, representing the typical wave frequency under moderate sea states, taken as 1.0 rad / s (corresponding to a wave period of approximately 6.28 seconds). This value is determined based on statistical analysis of a large amount of measured sea state data, representing the wave frequency range with the highest frequency throughout the year. The process noise covariance matrix of the Kalman filter is a 9×9 diagonal matrix, where the diagonal elements are the nominal process noise variance of each state component under calm sea states, and the off-diagonal elements are all 0, indicating that the process noise between each state component is uncorrelated. The nominal process noise variance is jointly calibrated through offline system identification and Allan variance analysis: the variance of the position component comes from the static statistics of GNSS positioning noise, the variance of the velocity component comes from the statistical value of the sum of GNSS velocity measurement noise and IMU velocity random walk, and the variance of the attitude angle component comes from the error statistics after integrating the IMU angular velocity random walk. The above calibration work is completed in advance during the filter design stage and stored in the system as fixed parameters. This nominal matrix will be adaptively scaled according to the wave dominant frequency in subsequent steps to match the excitation intensity of the actual sea state.

[0125] Step S4030202: Calculate the noise scaling factor based on the relative deviation between the wave dominant frequency and the reference frequency.

[0126] The noise scaling factor satisfies the following formula:

[0127] in, This is the process noise scaling factor. This is a preset adjustment coefficient used to control the scaling intensity. It can be determined through offline parameter debugging. In this embodiment, it is taken as... The scaling factor is designed based on the following: the higher the dominant wave frequency, the faster the oscillation frequency of the buoy's motion, the lower the tracking ability of the damped simple harmonic oscillator model for high-frequency disturbances, and the greater the instantaneous error of the model prediction. Therefore, it is necessary to appropriately amplify the process noise to reduce the confidence in the dynamic model. Conversely, the lower the dominant wave frequency, the more quasi-static the buoy's motion tends to be, and the more accurate the model prediction. Therefore, it is necessary to reduce the process noise to strengthen the model constraint. The above scaling is achieved through... It directly affects the diagonal elements of the velocity component in the process noise covariance matrix, enabling the filter's confidence level in velocity prediction to adaptively adjust with wave frequency characteristics. To prevent numerical anomalies caused by excessively large or small scaling factors, a lower bound protection is applied to the scaling factor. and protection of the upper realm ,in, The minimum allowable scaling factor is preset (0.1 in this embodiment). To preset the maximum allowable scaling factor (10 in this embodiment), we ensure that the scaled process noise covariance matrix always maintains positive semidefiniteness.

[0128] Step S4030203: Adjust the diagonal elements corresponding to the velocity components in the process noise covariance matrix using the noise scaling factor to obtain an adaptive process noise covariance matrix.

[0129] In this embodiment, the diagonal elements corresponding to the velocity components in the nominal process noise covariance matrix are scaled using the calculated noise scaling factor; that is, the eastward velocity variance, northward velocity variance, and vertical velocity variance are all multiplied by the scaling factor. The process noise variances corresponding to the position and attitude angle components remain unchanged because changes in the wave dominant frequency directly affect the prediction accuracy of the dynamic model for velocity state, while the process noise of position and attitude angles is not directly related to the wave dominant frequency. The scaled velocity variance, together with the unscaled position and attitude angle variances, constitutes the adaptive process noise covariance matrix at the current moment. This matrix will be used for the time prediction recursion of the prior covariance matrix in subsequent steps, enabling the filter to automatically adjust its confidence level in the physical model based on real-time sea conditions.

[0130] Step S40303: Obtain the prior covariance matrix by recursively deriving the state covariance matrix based on the adaptive process noise covariance matrix, the vertical motion discrete-time state transition matrix, the horizontal motion discrete state transition matrix, and the attitude transition matrix.

[0131] In this embodiment, the vertical motion discrete-time state transition matrix is ​​first... Discrete state transition matrix for eastward horizontal motion Discrete state transition matrix for northward horizontal motion And the Euler angle kinematic Jacobian matrix corresponding to the attitude components, according to the Kalman filter state vector The components are arranged in order to form a complete unit. State transition matrix Since the state variables of the vertical, horizontal (east / north) and attitude subsystems are physically independent, the off-diagonal elements across subsystems in the assembled complete state transition matrix are all 0, meaning that the state predictions between different subsystems do not affect each other. Within the same subsystem, there is a non-zero mapping relationship between the position and velocity components, determined by the dynamic model, which conforms to the state transition law of linear systems. For the attitude components, because their numerical recursion uses quaternions while their covariance recursion uses the Euler angle kinematics Jacobian matrix, the linearization mapping method of the attitude subsystem differs from other subsystems, but the overall assembly logic is consistent, conforming to the standard processing method of extended Kalman filtering.

[0132] The prior covariance matrix is ​​recursively derived using the time prediction formula of the standard Kalman filter:

[0133] in, For the current moment (the th) The prior covariance matrix of (step) The complete state transition matrix after assembly ( matrix), For the previous moment (the first moment) The posterior covariance matrix of (steps), This represents the matrix transpose operation. Let be the adaptive process noise covariance matrix at the current moment. This represents the discrete-time step number. This recursion transforms the uncertainty of the state estimate from the previous time step into the state transition matrix. Propagate to the current time step and superimpose the process noise covariance matrix at the current time step. The uncertainty of the current state prediction (i.e., the prior covariance matrix) is obtained. The recursively obtained prior covariance matrix, together with the state prediction value, forms a complete prior estimate pair, which is used in the subsequent measurement update process.

[0134] Step S40304: Form a prior state estimate based on the predicted state value and the prior covariance matrix.

[0135] In this embodiment, the predicted state values ​​and the prior covariance matrix are combined into a complete prior estimate pair. This prior estimate pair fully describes the statistical characteristics of the buoy's motion state at the current moment before the measurement update. The predicted state values ​​provide the optimal predictions for each state component, while the prior covariance matrix quantifies the degree of uncertainty of these predictions. After combination, this prior estimate pair will serve as input for subsequent GNSS measurement updates and IMU accelerated measurement updates. The measurement update process will calculate the posterior state estimate based on the difference between this prior information and the sensor observations, thereby correcting the predicted state values.

[0136] Step S404: Construct a first observation matrix based on the position and velocity data from the Global Navigation Satellite System.

[0137] In this embodiment, a first observation matrix is ​​constructed based on the three-dimensional position and three-dimensional velocity in the GNSS data. The GNSS observation vector is a six-dimensional vector, containing the position and velocity in three directions in the navigation coordinate system. The first observation matrix is ​​a 6×9 matrix used to map the GNSS observation vector to the corresponding position and velocity components in the state vector; the positions in the matrix corresponding to the position and velocity of the state vector are 1, and the remaining positions are 0, meaning that this matrix directly extracts the three-dimensional position and three-dimensional velocity from the state vector as the predicted GNSS observation values. After construction, this observation matrix will be used in step S405 to calculate the GNSS observation residual and Kalman gain, and simultaneously continuously adjust the GNSS observation noise covariance matrix according to the GNSS signal quality indication parameters to achieve adaptive weighting of GNSS measurements.

[0138] Step S405: Based on the first observation matrix, perform a first measurement update on the state prior estimate to obtain the state estimate after the first measurement update.

[0139] In this embodiment, the prior state estimate formed in step S40304 is used as input, and Kalman filtering measurement updates are performed using the first observation matrix and GNSS observation vector constructed in step S404. First, the GNSS observation residual is calculated, which is the difference between the actual GNSS observation value and the predicted state observation value obtained by projecting through the first observation matrix. Then, the Kalman gain is calculated based on the prior covariance matrix, the first observation matrix, and the GNSS observation noise covariance matrix. The GNSS observation noise covariance matrix is ​​continuously adjusted according to the signal quality indication parameters (position accuracy factor PDOP and carrier-to-noise ratio C / N0)—its diagonal elements approach the nominal value when the signal quality is good, and automatically increase when the signal quality deteriorates to weaken the correction effect of GNSS measurements. Subsequently, the calculated Kalman gain is used to correct the predicted state value, resulting in the updated state estimate after GNSS measurement, while simultaneously updating the covariance matrix to reflect the reduced uncertainty. When the GNSS signal is completely unlocked, the Kalman gain approaches zero, and the GNSS measurement does not produce any correction effect on the state estimate; the state estimate relies entirely on the physical model recursion.

[0140] Step S406: Construct a second observation matrix based on the acceleration data from the inertial measurement unit.

[0141] In this embodiment, a second observation matrix is ​​constructed based on the triaxial accelerations in the inertial measurement unit data. First, a rotation matrix is ​​constructed using the prior attitude angles output by the time prediction. The specific force acceleration in the body coordinate system is projected onto the navigation coordinate system, and the gravity vector is subtracted to obtain the actual motion acceleration observation value in the navigation coordinate system. The second observation matrix is ​​a 3×9 matrix used to map the acceleration predicted by the physical model to the acceleration observation vector. Its non-zero elements are determined by the physical models of vertical and horizontal motion. In the vertical channel, the partial derivative of the predicted acceleration with respect to the vertical position is the square of the negative natural frequency, and the partial derivative with respect to the vertical velocity is twice the damping ratio multiplied by the natural frequency. In the horizontal channel, the partial derivative of the predicted acceleration with respect to the horizontal velocity is the negative viscous damping coefficient in the corresponding direction. These partial derivative terms are derived from the natural frequency and damping ratio calculated in step S2, and the viscous damping coefficient identified in step S304, respectively. Once constructed, the observation matrix will be used to calculate the observation residuals and Kalman gain in subsequent IMU acceleration measurement updates; the IMU acceleration observation noise covariance matrix will also be continuously adjusted according to the acceleration confidence weights, and its diagonal elements will automatically increase when the confidence of the main IMU decreases in order to weaken the measurement correction effect.

[0142] Step S407: Based on the second observation matrix, perform a second measurement update on the state estimate updated by the first measurement to obtain the state estimate updated by the second measurement.

[0143] In this embodiment, the state estimate updated by the first measurement obtained in step S405 is used as input, and a second Kalman filter measurement update is performed using the second observation matrix and IMU acceleration observation vector constructed in step S406. First, the IMU acceleration observation residual is calculated, which is the difference between the IMU measured acceleration (after deducting gravity and converting to the navigation frame) and the acceleration predicted by the physical model obtained by projecting through the second observation matrix. This residual reflects the deviation between the physical model prediction and the sensor measurement. Then, the Kalman gain is calculated based on the covariance matrix updated by the first measurement, the second observation matrix, and the IMU acceleration observation noise covariance matrix, where the IMU acceleration observation noise covariance matrix is ​​continuously adjusted according to the acceleration confidence weight. Subsequently, the calculated Kalman gain is used to further correct the state estimate updated by the first measurement, resulting in the state estimate updated by the second measurement. Simultaneously, the covariance matrix is ​​updated to reflect the further reduction in uncertainty. This step allows the difference between the physical model prediction and the IMU measured acceleration to be used to correct the position and velocity components in the state estimate, effectively compensating for the recursive bias of the physical model.

[0144] Step S408: Obtain the optimal estimate of the buoy's motion state based on the updated state estimate from the second measurement.

[0145] In this embodiment, the updated state estimate from the second measurement is output as the optimal estimate of the buoy's motion state at the current moment. Simultaneously, the updated covariance matrix from the second measurement is used as the optimal estimated covariance matrix for the current moment, which is then used for time prediction recursion in the next filtering cycle. Thus, the Kalman filter completes a full recursion cycle, sequentially going through three stages: time prediction (state prediction and covariance prediction), GNSS measurement update, and IMU accelerometer measurement update. The output optimal estimate of the buoy's motion state includes the buoy's three-dimensional position, three-dimensional velocity, and three-dimensional attitude angles in the navigation coordinate system at the current moment.

[0146] Step S5: Calculate the actual wind speed based on the optimal estimate of the buoy's motion state and the anemometer data.

[0147] In this embodiment, the true wind speed is calculated based on the optimal estimate of the buoy's motion state output in step S408 and the anemometer data obtained in step S101. First, the three-dimensional velocity and three-dimensional attitude angle of the buoy's center of mass in the navigation coordinate system are extracted from the optimal estimate of the buoy's motion state. Simultaneously, the installation position vector of the anemometer on the buoy (defined in the body coordinate system) is obtained. A rotation matrix from the body coordinate system to the navigation coordinate system is constructed using the attitude angle. The additional velocity generated by the buoy's rotation at the anemometer's installation position is calculated by the cross product of the angular velocity and the installation position vector (in the body coordinate system), and then transformed to the navigation coordinate system using the rotation matrix. The transformed rotational additional velocity is added to the translational velocity of the buoy's center of mass to obtain the total additional motion velocity of the anemometer's measurement center in the navigation coordinate system. The relative wind speed vector output by the anemometer (in the body coordinate system) is transformed to the navigation coordinate system using the same rotation matrix. According to the principle of velocity vector synthesis, the true wind speed is equal to the vector sum of the relative wind speed and the total additional motion velocity of the anemometer, i.e.:

[0148] in, The true wind speed vector in the navigation coordinate system. The relative wind speed vector in the navigation coordinate system. This represents the total additional velocity of the anemometer in the navigation coordinate system. Physically, the three-dimensional ultrasonic anemometer measures the velocity vector of the wind relative to the instrument; the true wind speed equals the vector sum of the relative wind speed and the sensor's velocity. The final output true wind speed includes three-dimensional components in the navigation coordinate system. Its magnitude can be obtained by taking the square root of the sum of the squares of the three components, and its direction is determined by all three components. This true wind speed completely eliminates the influence of the buoy's own six degrees of freedom motion (translation and rotation) on wind speed measurement, maintaining high accuracy even under conditions of intense buoy movement.

[0149] In the entire technical solution of this invention, the buoy is regarded as a damped oscillating system driven by wave forces. Its vertical motion naturally satisfies the damped simple harmonic motion equation determined by its own hydrodynamic parameters, namely its natural frequency and damping ratio. Its horizontal motion satisfies the attenuation law controlled by the viscous damping coefficient. These physical laws can independently describe the buoy's motion evolution trend without relying on any external sensor observations. Based on this objective fact, this invention completely abandons the traditional architecture of using inertial measurement unit acceleration integral to drive state recursion, and instead uses a pure physical model to drive state prediction. The inertial measurement unit acceleration is only used as a measurement value to correct the model prediction deviation, fundamentally avoiding the error of repeatedly subtracting the restoring force from the inertial measurement unit measurement value that already includes the damping effect. At the same time, the wave dominant frequency, as a characterization index of the external excitation intensity, is only used to adjust the filter's confidence in the physical model, i.e., the process noise covariance, rather than changing the model itself, thus ensuring the purity of the physical model. This filtering architecture, which combines physical model prediction with sensor measurement correction, ensures that even in extreme sea conditions where the global navigation satellite system is completely unlocked and the inertial measurement unit acceleration is distorted due to range saturation, the buoy's vertical motion can still be constrained within a physically reasonable bounded range by its natural frequency and damping ratio, and its horizontal motion can still be constrained by the viscous damping coefficient to be a stable exponential decay. The state estimation will not exhibit the divergence phenomenon of pure inertial integral error accumulating twice over time. Furthermore, the adaptive fusion of dual inertial measurement units ensures the continuity of inertial data, and the wave dominant frequency adaptively adjusts the filtering parameters to match the real-time sea conditions. Ultimately, it achieves continuous, stable, and high-precision estimation of the buoy's six degrees of freedom motion state under all sea conditions, and completes accurate motion compensation for wind speed measurement data.

[0150] like Figure 2 As shown, in another aspect, the present invention also provides a buoy wind speed dynamic compensation system based on wave physics constraints, comprising: a processor, an input device, an output device, and a memory, wherein the processor, the input device, the output device, and the memory are interconnected, wherein the memory is used to store a computer program, the computer program includes program instructions, and the processor is configured to call the program instructions to execute the relevant steps of the relevant embodiments of the buoy wind speed dynamic compensation method based on wave physics constraints of the present invention.

[0151] This invention provides a buoy wind speed dynamic compensation system based on wave physics constraints. The functional components can be integrated into a single processing unit, exist as individual physical entities, or be integrated into a single unit. The integrated components can be implemented in hardware or as software functions.

[0152] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, and not to limit them. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some or all of the technical features therein. Such modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the scope of the technical solutions of the embodiments of the present invention, and they should all be covered within the scope of the claims and specification of the present invention.

Claims

1. A dynamic compensation method for buoy wind speed based on wave physics constraints, characterized in that, The method includes: Acquire data from the buoy's inertial measurement unit, global navigation satellite system, anemometer, and buoy design parameters; Calculate the natural frequency and damping ratio of the buoy's vertical motion based on the buoy design parameters; The wave dominant frequency and viscous damping coefficient are calculated based on the inertial measurement unit data and the global navigation satellite system data. Based on the natural frequency, the damping ratio, the wave dominant frequency, and the viscous damping coefficient, the Kalman filter algorithm is used to fuse the inertial measurement unit data and the global navigation satellite system data to obtain the optimal estimate of the buoy's motion state. The actual wind speed is calculated based on the optimal estimate of the buoy's motion state and the anemometer data.

2. The buoy wind speed dynamic compensation method based on wave physics constraints according to claim 1, characterized in that, The acquisition of inertial measurement unit data, global navigation satellite system data, anemometer data, and buoy design parameters of the buoy includes: Acquire data from the buoy's main inertial measurement unit, auxiliary inertial measurement unit, global navigation satellite system, anemometer, and buoy design parameters; Inertial measurement unit data is obtained based on the data from the primary inertial measurement unit and the data from the secondary inertial measurement unit.

3. The buoy wind speed dynamic compensation method based on wave physics constraints according to claim 2, characterized in that, The process of obtaining inertial measurement unit data based on the primary inertial measurement unit data and the secondary inertial measurement unit data includes: Extract the main acceleration data, main velocity data, auxiliary acceleration data, and auxiliary angular velocity data from the main inertial measurement unit data and the auxiliary inertial measurement unit data, respectively. Sliding sampling is performed on the main acceleration data and the main velocity data in a predefined time window to obtain multiple time-series window acceleration data and time-series window angular velocity data; The window acceleration variance and the combined amplitude of the window endpoint acceleration are calculated based on the time-series window acceleration data, and the window angular velocity variance and the combined amplitude of the window endpoint angular velocity are calculated based on the time-series window angular velocity data. The confidence weight of the principal acceleration is calculated based on the variance of the window acceleration and the comprehensive amplitude of the acceleration at the end of the window, and the confidence weight of the principal velocity is calculated based on the variance of the window angular velocity and the comprehensive amplitude of the angular velocity at the end of the window. The primary acceleration data and the secondary acceleration data are weighted and fused according to the primary acceleration confidence weight to obtain fused time-series acceleration data; Based on the confidence weight of the protagonist's velocity, the protagonist's velocity data and the auxiliary angular velocity data are weighted and fused to obtain fused temporal angular velocity data; Inertial measurement unit data is formed based on the fused temporal acceleration data and the fused temporal angular velocity data.

4. The buoy wind speed dynamic compensation method based on wave physics constraints according to claim 1, characterized in that, The calculation of wave dominant frequency and viscous damping coefficient based on the inertial measurement unit data and the global navigation satellite system data includes: Extract the vertical acceleration component data from the inertial measurement unit data; The wave dominant frequency is calculated based on the vertical acceleration component data; Extract horizontal velocity timing data and signal quality indication parameters from the Global Navigation Satellite System data; When the signal quality indication parameter meets the preset availability conditions, the viscous damping coefficient is calculated using the horizontal velocity timing data.

5. The buoy wind speed dynamic compensation method based on wave physics constraints according to claim 1, characterized in that, The step of fusing the inertial measurement unit data and the global navigation satellite system data using a Kalman filter algorithm based on the natural frequency, the damping ratio, the wave dominant frequency, and the viscous damping coefficient to obtain the optimal estimate of the buoy's motion state includes: Set the state vector of the Kalman filter to include three-dimensional position, three-dimensional velocity, and three-dimensional attitude angle; An attitude transfer matrix is ​​constructed using the angular velocity in the inertial measurement unit data as the control input, and a vertical motion discrete-time state transfer matrix and a horizontal motion discrete-time state transfer matrix are constructed based on the natural frequency, the damping ratio, and the viscous damping coefficient. By using the vertical motion discrete-time state transition matrix, the horizontal motion discrete state transition matrix, and the attitude transition matrix, and in conjunction with the wave dominant frequency, a state prior estimate is obtained by performing time prediction on the state vector. A first observation matrix is ​​constructed based on the position and velocity data from the global navigation satellite system. Based on the first observation matrix, the state prior estimate is updated by a first measurement to obtain the state estimate after the first measurement update; A second observation matrix is ​​constructed based on the acceleration data from the inertial measurement unit. Based on the second observation matrix, the state estimate updated by the first measurement is updated by the second measurement to obtain the state estimate updated by the second measurement. The optimal estimate of the buoy's motion state is obtained based on the updated state estimate from the second measurement.

6. The buoy wind speed dynamic compensation method based on wave physics constraints according to claim 5, characterized in that, The process of constructing an attitude transfer matrix using the angular velocity from the inertial measurement unit data as control input, and constructing a vertical motion discrete-time state transfer matrix and a horizontal motion discrete-time state transfer matrix based on the natural frequency, the damping ratio, and the viscous damping coefficient, includes: Construct a vertical motion continuous-time system matrix based on the natural frequency and the damping ratio; The vertical motion discrete-time system matrix is ​​discretized using the matrix exponentiation method to obtain the vertical motion discrete-time state transition matrix; Construct a discrete state transition matrix for horizontal motion using the aforementioned viscous damping coefficient; Based on the quaternion kinematic equations, an attitude transfer matrix is ​​constructed using the angular velocity from the inertial measurement unit data as the control input.

7. The buoy wind speed dynamic compensation method based on wave physics constraints according to claim 5, characterized in that, The step of obtaining a priori state estimation by performing time prediction on the state vector using the vertical motion discrete-time state transition matrix, the horizontal motion discrete state transition matrix, and the attitude transition matrix, combined with the wave dominant frequency, includes: The state vector is recursively derived using the vertical motion discrete-time state transition matrix, the horizontal motion discrete state transition matrix, and the attitude transition matrix to obtain the state prediction value. The adaptive process noise covariance matrix is ​​obtained by adjusting the process noise covariance matrix of the Kalman filter based on the wave dominant frequency. The prior covariance matrix is ​​obtained by recursively deriving the state covariance matrix based on the adaptive process noise covariance matrix, the vertical motion discrete-time state transition matrix, the horizontal motion discrete state transition matrix, and the attitude transition matrix. A prior state estimate is formed based on the predicted state value and the prior covariance matrix.

8. The buoy wind speed dynamic compensation method based on wave physics constraints according to claim 7, characterized in that, The process of recursively deriving the state vector from the state vector using the vertical motion discrete-time state transition matrix, the horizontal motion discrete state transition matrix, and the attitude transition matrix to obtain the state prediction value includes: The vertical position and vertical velocity of the buoy are recursively derived using the vertical motion discrete-time state transition matrix to obtain the vertical recursive data. The horizontal position and horizontal velocity of the buoy are recursively derived using the horizontal motion discrete state transition matrix to obtain the horizontal recursive data; The attitude angles in the state vector are recursively derived from the attitude transition matrix to obtain the recursively derived attitude angles. Based on the vertically recursed data, the horizontally recursed data, and the recursed attitude angles, a state prediction value is formed.

9. The buoy wind speed dynamic compensation method based on wave physics constraints according to claim 7, characterized in that, The process noise covariance matrix obtained by adjusting the Kalman filter process noise covariance matrix according to the wave dominant frequency includes: Obtain the preset reference frequency and the process noise covariance matrix of the Kalman filter; The noise scaling factor is calculated based on the relative deviation between the wave's dominant frequency and the reference frequency; The diagonal elements corresponding to the velocity components in the process noise covariance matrix are adjusted using the noise scaling factor to obtain an adaptive process noise covariance matrix.

10. A buoy wind speed dynamic compensation system based on wave physics constraints, characterized in that, include: The system includes a processor, an input device, an output device, and a memory, all interconnected. The memory stores a computer program, which includes program instructions. The processor is configured to invoke the program instructions to execute a buoy wind speed dynamic compensation method based on wave physics constraints as described in any one of claims 1 to 9.

Images