A marine buoy multi-source fusion positioning method and system based on factor graph optimization

Abstract

The application discloses a marine buoy multi-source fusion positioning method and system based on a factor graph optimization, which acquires GNSS observation data, IMU measurement data and marine environment auxiliary data; a factor graph containing state nodes, GNSS position factors, IMU pre-integration factors and marine dynamics constraint factors is constructed, and wave and current theories are used to constrain buoy movement; in view of multipath effects, marine surface reflection geometry and marine root mean square wave height are combined to calculate a multipath weighting factor, and a GNSS covariance matrix is adaptively adjusted; according to IMU data, a sea state level is discriminated, and an edge window length and a trigger interval of incremental smoothing solving are adaptively linked and adjusted; through adjacent buoy ranging information, collaborative constraints are constructed, and based on Mahalanobis distance and chi-square distribution threshold value detection, abnormalities are detected and local reconstruction is performed. The application effectively suppresses marine surface multipath interference, slows down the accumulation of calculation errors during signal interruption, and realizes high-availability continuous positioning under limited computing power.

Description

Technical Field

[0001] This invention relates to the field of marine navigation and positioning technology, specifically to a multi-source fusion positioning method and system for marine buoys based on factor graph optimization. Background Technology

[0002] Ocean buoys are crucial data acquisition platforms for marine hydrological environment monitoring, resource exploration, and maritime traffic safety. To achieve precise spatial matching of various marine physical parameters in target sea areas, the buoys themselves must possess high-precision real-time position coordinates. Therefore, providing stable and reliable positioning services for ocean buoys is a fundamental prerequisite for conducting various marine scientific research and engineering applications.

[0003] Currently, in routine marine observation applications, the industry commonly employs a combined positioning technology that combines Global Navigation Satellite System (GNSS) and Inertial Measurement Unit (IMU) to obtain buoy position information. Traditional combined positioning schemes often utilize basic filtering algorithms or standard map optimization models, directly fusing absolute position observation data provided by GNSS with relative motion data calculated by the IMU, aiming to achieve continuous navigation state estimation of the buoy in dynamic sea surface environments.

[0004] However, in actual marine environments, due to the strong reflective properties of the water surface, satellite signals are highly susceptible to wave reflections when reaching the buoy antenna, resulting in significant multipath errors. Existing fusion positioning schemes typically adopt conventional multipath suppression strategies used in terrestrial scenarios, without specifically building models and performing robust optimizations to address the unique geometric relationships of sea surface reflections and variations in sea surface roughness in marine environments. This leads to the system's inability to accurately assess the actual impact of multipath interference on satellite observation data. Observation data containing reflection errors fails to undergo reasonable adaptive weighting, causing these errors to be directly introduced into the positioning solution system. Ultimately, this reduces the reliability of satellite observation data acquired by the buoy, and the overall accuracy of multi-source fusion positioning significantly decreases under actual sea conditions. Summary of the Invention

[0005] To address the shortcomings of existing technologies, this invention provides a multi-source fusion positioning method and system for marine buoys based on factor graph optimization. This method solves the problems of multipath interference of satellite signals caused by strong reflection of sea waves under dynamic sea conditions, and the rapid divergence of pure inertial calculation errors when signals are interrupted due to violent swaying.

[0006] To achieve the above objectives, the present invention provides the following technical solution: The first aspect of the present invention provides a multi-source fusion positioning method for ocean buoys based on factor graph optimization, comprising:

[0007] The data acquisition module performs multi-source observation data acquisition and preprocessing, acquiring GNSS observation data, IMU measurement data, and marine environmental auxiliary data; it aligns the timestamps of all observation data with GPS time as a reference; and it classifies satellite observations into three levels—high quality, medium quality, and low quality—based on the carrier-to-noise ratio, marking low-quality observations as suspected multipath pollution on the sea surface.

[0008] The factor graph construction module constructs factor graph state nodes and multi-source factors based on the preprocessed multi-source data; it constructs navigation state nodes that include the rotation matrix from the carrier coordinate system to the navigation coordinate system, three-dimensional position vector, three-dimensional velocity vector, accelerometer deviation, and gyroscope deviation; and it constructs multi-source factors that include GNSS position factors, IMU pre-integration factors, IMU deviation constraint factors, and ocean dynamics constraint factors.

[0009] Among them, the wave motion constraint factor is constructed using the local mean sea level height, the predicted tidal height, and the amount of sea surface undulation caused by waves, which constrains the vertical displacement of the buoy to the vicinity of the sea surface; the ocean current drift constraint factor is constructed using the predicted ocean current velocity vector and the ocean current constraint covariance matrix, which constrains the horizontal motion velocity of the buoy to the vicinity of the ocean current velocity.

[0010] The factor graph construction module adaptively weights the GNSS observation data for sea surface multipath effects based on sea state parameters; it determines the grazing angle of sea surface reflection by the satellite elevation angle and the buoy antenna height above the sea surface, and calculates the sea surface reflection coefficient based on the Fresnel reflection coefficient of a smooth sea surface, the root mean square wave height of the sea surface, and the signal carrier wavelength; it constructs multipath weighting factors based on the sea surface reflection coefficient, the multipath penalty coefficient, and the ratio of the observed carrier-to-noise ratio to the reference carrier-to-noise ratio threshold; it constructs a diagonal matrix using the multipath weighting factors corresponding to the three-dimensional coordinate axes, and adjusts the GNSS covariance matrix for GNSS factor construction;

[0011] The adaptive optimization module acquires sea state parameters and performs adaptive incremental optimization based on these parameters. It calculates the intensity of buoy motion based on the squared mean of the difference between the magnitude of the triaxial accelerometer measurements and the scalar value of gravitational acceleration within the evaluation window, and classifies sea state levels accordingly. It adaptively adjusts the marginalization window duration, optimizes the trigger interval, and adjusts the variance parameters of wave constraints based on the sea state level. The iSAM2 incremental smoothing algorithm is used to solve the optimization objective function, and fixed-lag marginalization is performed on historical state nodes that exceed the marginalization window.

[0012] In the presence of multiple buoys, the collaborative correction module uses the ranging information between adjacent buoys to perform multi-buoy collaborative positioning correction; it constructs a multi-buoy collaborative distance constraint factor using the distance difference between position estimates and ranging information; when the GNSS signal of this buoy is completely interrupted for more than a preset time, and the GNSS signal of at least one neighboring buoy is available, the collaborative distance constraint factor is added to the factor graph of each buoy to participate in joint optimization.

[0013] The anomaly detection module performs anomaly drift detection and factor graph reconstruction on the positioning sequence and outputs the ocean buoy positioning results; it calculates the Mahalanobis distance between the current position estimate and the position predicted based on the previous moment's state and the ocean current model, and determines an anomaly when the Mahalanobis distance exceeds the threshold set based on the chi-square distribution; it marks the state nodes before and after the anomaly moment as nodes to be reconstructed and removes them from the current factor graph, and reconstructs the local factor graph with the nearest reliable state node as the new prior constraint and performs local optimization.

[0014] A second aspect of the present invention provides a multi-source fusion positioning system for ocean buoys based on factor graph optimization, comprising:

[0015] The data acquisition module is used to acquire GNSS observation data, IMU measurement data, and marine environmental auxiliary data, and to perform data acquisition and preprocessing operations.

[0016] The factor graph construction module is used to construct navigation status nodes and multi-source measurement factors, and to perform factor graph construction and sea surface multipath adaptive weighting operations.

[0017] The adaptive optimization module is used to perform adaptive incremental optimization solutions based on sea state levels, including sea state discrimination, marginal window adjustment, and iSAM2 optimization operations.

[0018] The collaborative correction module is used to perform collaborative positioning correction using distance measurement information between multiple buoys, and to perform collaborative constraint construction and joint optimization operations.

[0019] The anomaly detection module is used to detect anomaly localization results and perform local reconstruction of the factor graph, performing anomaly detection and reconstruction operations.

[0020] This invention provides a multi-source fusion positioning method and system for ocean buoys based on factor graph optimization. It has the following beneficial effects:

[0021] 1. This invention determines the glancing angle of sea surface reflection by combining the satellite elevation angle and the buoy antenna height above the sea surface, then calculates the sea surface reflection coefficient by combining the root mean square wave height of the sea surface, and finally constructs a multipath weighting factor by combining the observation carrier-to-noise ratio to dynamically adjust the GNSS covariance matrix for the construction of GNSS factors. This achieves effective suppression of multipath interference errors on the sea surface in the marine environment, thereby improving the reliability of satellite observation data acquired by the buoy and the overall accuracy of the final multi-source fusion positioning.

[0022] 2. This invention utilizes linear wave theory to constrain the vertical displacement of the buoy to near the sea surface and uses ocean current forecast velocity vectors to constrain the horizontal motion velocity of the buoy to near the ocean current velocity to construct an ocean dynamics constraint factor. Then, the ocean dynamics constraint factor and the IMU pre-integration factor are jointly integrated into the factor graph framework, which provides auxiliary constraints when using pre-integration for short-term position propagation during GNSS signal interruption, thereby effectively mitigating the error accumulation caused by pure inertial calculation and maintaining the continuity of positioning services.

[0023] 3. This invention calculates the intensity index of buoy motion based on the magnitude of triaxial accelerometer measurements and the scalar value of gravitational acceleration to classify sea state levels. It adaptively adjusts the marginalization window duration and optimizes the trigger interval according to the sea state level. At the same time, for abnormal positioning results where the Mahalanobis distance exceeds the threshold set based on the chi-square distribution, it reconstructs the local factor graph and performs local optimization, thereby achieving a dynamic balance between positioning accuracy and computational resource consumption. This satisfies the real-time requirements of buoy calculation and prevents abnormal errors from accumulating in global optimization. Attached Figure Description

[0024] Figure 1 This invention provides a flowchart of the technical route for a multi-source fusion positioning method for ocean buoys based on factor graph optimization, as shown in the embodiments of the present invention.

[0025] Figure 2 A comparison diagram of the positioning trajectories in a two-dimensional plane for different methods provided in the embodiments of the present invention;

[0026] Figure 3 A time series comparison chart of positioning errors of different methods provided in the embodiments of the present invention;

[0027] Figure 4 A comparison curve of service availability for different positioning methods provided in embodiments of the present invention;

[0028] Figure 5 This is a performance comparison chart of the sea state adaptive edge-mapping strategy and the fixed window strategy provided in the embodiments of the present invention. Detailed Implementation

[0029] The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0030] Please see the appendix Figure 1-5 This invention provides a marine buoy positioning device, which includes a GNSS receiver, a MEMS inertial measurement unit, an embedded processor, an inter-buoy communication module, and a power supply unit. The embedded processor runs a software program for a multi-source fusion positioning method for marine buoys based on factor graph optimization.

[0031] This embodiment also provides a multi-source fusion positioning system for ocean buoys based on factor graph optimization. The multi-source fusion positioning system for ocean buoys based on factor graph optimization includes a data acquisition module, a factor graph construction module, an adaptive optimization module, a collaborative correction module, and an anomaly detection module.

[0032] The data acquisition module is used to acquire GNSS observation data, IMU measurement data, and marine environmental auxiliary data, and to perform data acquisition and preprocessing operations.

[0033] The factor graph construction module is used to construct navigation status nodes and multi-source measurement factors, and to perform factor graph construction and sea surface multipath adaptive weighting operations.

[0034] The adaptive optimization module is used to perform adaptive incremental optimization solutions based on sea state levels, including sea state discrimination, marginal window adjustment, and iSAM2 optimization operations.

[0035] The collaborative correction module is used to perform collaborative positioning correction using distance measurement information between multiple buoys, and to perform collaborative constraint construction and joint optimization operations.

[0036] The anomaly detection module is used to detect anomaly localization results and perform local reconstruction of the factor graph, performing anomaly detection and reconstruction operations.

[0037] Combining the above-mentioned factor graph-based multi-source fusion positioning system for ocean buoys, the macroscopic overall workflow of the factor graph-based multi-source fusion positioning method for ocean buoys described in this invention includes the following steps:

[0038] Step S1: The data acquisition module performs multi-source observation data acquisition and preprocessing to obtain GNSS observation data, IMU measurement data, and marine environmental auxiliary data.

[0039] Step S2: The factor graph construction module constructs factor graph state nodes and multi-source factors based on the preprocessed multi-source data.

[0040] Step S3: The factor graph construction module performs adaptive weighting of the sea surface multipath effect on the GNSS observation data based on sea state parameters.

[0041] Step S4: The adaptive optimization module obtains sea state parameters and performs adaptive incremental optimization based on the sea state parameters.

[0042] Step S5: In the presence of multiple buoys, the collaborative correction module uses the ranging information between adjacent buoys to perform multi-buoy collaborative positioning correction.

[0043] Step S6: The anomaly detection module performs anomaly drift detection and factor graph reconstruction on the positioning sequence and outputs the ocean buoy positioning results.

[0044] Specifically, step S1 includes the following:

[0045] Step S101: Acquire GNSS observation data. The GNSS receiver mounted on the buoy acquires satellite observation data at a sampling rate of 1Hz. The satellite observation data includes pseudorange observations, Doppler observations, and the carrier-to-noise ratio of each satellite. and elevation angle The GNSS receiver obtains the position solution using the weighted least squares method. and the covariance matrix corresponding to the location solution .

[0046] Step S102: Acquire IMU measurement data. The MEMS inertial measurement unit mounted on the buoy acquires triaxial accelerometer measurements at a sampling rate of 100Hz to 200Hz. and the measurement value of the three-axis gyroscope The IMU measurement model is as follows:

[0047]

[0048]

[0049] in, This represents the rotation matrix from the navigation coordinate system to the vehicle coordinate system. Represents linear acceleration in the navigation coordinate system. Represents the gravitational acceleration vector. This represents the true angular velocity in the carrier coordinate system. Indicates accelerometer deviation. Indicates gyroscope deviation. This indicates white noise from the accelerometer. This represents the white noise of the gyroscope.

[0050] Step S103: Acquire marine environmental auxiliary data. The data acquisition module collects marine environmental auxiliary data, including: predicted tidal height. Sea state parameters and predicted ocean current speed in the area where the buoy is located. .

[0051] Step S104: Data Time Reference Unification and Signal Quality Assessment. Unify the time reference for all data sources, aligning the timestamps of all observation data with GPS time as the reference. Assess the signal quality of GNSS observations based on the carrier-to-noise ratio. Satellite observations are categorized into three levels: high quality, medium quality, and low quality, corresponding to... dB-Hz, dB-Hz and dB-Hz. Low-quality observations are flagged as suspected multipath pollution on the sea surface.

[0052] Specifically, step S2 includes the following:

[0053] Step S201: Construct the state nodes of the factor graph. Navigation state nodes at each time step k. Defined as:

[0054]

[0055] in, This represents the rotation matrix from the vehicle coordinate system to the navigation coordinate system. Represents a three-dimensional position vector. This represents a three-dimensional velocity vector.

[0056] Step S202: Construct GNSS location factors. When GNSS observations are available, the residuals of the GNSS factors... Defined as:

[0057]

[0058] in, This represents the GNSS covariance matrix after adaptive weighting for sea surface multipath effects in step S3. Indicates Mahalanobis distance .

[0059] Step S203: Construct the IMU pre-integration factor. Pre-integrate all IMU measurements between adjacent GNSS observation times into a single relative motion constraint to construct the IMU pre-integration factor. :

[0060]

[0061] in, This represents the IMU pre-integration residual. Let represent the pre-integral covariance matrix. The IMU bias is modeled as a slowly varying random walk process, and the constraint factor for the IMU bias is... for:

[0062]

[0063] in, and These represent the covariance of the random walk of the accelerometer bias and the covariance of the random walk of the gyroscope bias, respectively.

[0064] Step S204: Construct ocean dynamics constraint factors. The wave motion constraint factor, based on linear wave theory, constrains the vertical displacement of the buoy to near the sea surface. for:

[0065]

[0066] in, Represents the vertical position component of a state node. Indicates the local average sea level height. This represents the predicted tidal height at time k. This indicates the amount of sea surface undulation caused by waves. This represents the variance of the wave constraint. Wave undulation. Data is obtained through actual measurements by wave sensors mounted on the buoy, or based on significant wave height. Estimate the statistical range, take The ocean current drift constraint factor constrains the buoy's horizontal motion speed to near the ocean current speed. for:

[0067]

[0068] in, Represents the horizontal velocity component of the state node. This represents the predicted ocean current velocity vector at time k. This represents the ocean current constraint covariance matrix. The value of the ocean current constraint covariance matrix is ​​set according to the accuracy of ocean current forecasting.

[0069] Step S3: The factor map construction module adaptively weights the sea surface multipath effects on the GNSS observation data based on sea state parameters. For each visible satellite s, the weight of the sea surface multipath influence of satellite s is calculated, including the following steps:

[0070] Step S301, Grazing angle of reflection from the sea surface From satellite elevation angle and the height of the buoy antenna above the sea surface Confirmed, the formula is:

[0071]

[0072] in, This represents the horizontal distance from the reflection point of satellite s to the buoy, which is approximately... .

[0073] Step S302, Sea surface reflectance Based on the Fresnel reflection formula and the sea surface roughness factor, the formula is as follows:

[0074]

[0075] in, This represents the Fresnel reflection coefficient of a smooth sea surface. represents the sea surface root mean square wave height, This represents the carrier wavelength of the GNSS signal. The root mean square wave height of the sea surface. By Yoshihiko The estimated root-mean-square wave height at sea surface is approximately: .

[0076] Step S303: Based on the sea surface reflection coefficient and carrier-to-noise ratio, construct the multipath weighting factor for each satellite. The formula is:

[0077]

[0078] in, This is the multipath penalty coefficient, which ranges from 5 to 20. The reference carrier-to-noise ratio threshold is set to 40 dB-Hz. The noise ratio of satellite s is the observation carrier-to-noise ratio.

[0079] Step S304: Adjust the GNSS covariance matrix using the multipath weighting factor. The formula is as follows:

[0080]

[0081] in, This represents the GNSS covariance matrix before adjustment. This represents the weighted covariance matrix. Used for constructing GNSS factors in step S2. When a satellite's carrier-to-noise ratio is below 20 dB-Hz and its reflection coefficient is... When the value is greater than 0.6, the observations of the corresponding satellite will be removed from the GNSS solution.

[0082] Specifically, step S4 includes the following:

[0083] Step S401: Calculate the intensity index of buoy motion based on IMU acceleration measurement. The formula is:

[0084]

[0085] in, This indicates the number of IMU sampling points within the evaluation window. This represents the triaxial accelerometer measurement value at the i-th sampling point. The modulus representing the accelerometer measurement value This represents the scalar value of gravitational acceleration. The scalar value of gravitational acceleration is taken as... Based on the index of the intensity of buoy movement Sea state classification is based on the following rules: calm sea states correspond to... Moderate sea state corresponds to Corresponding to severe sea conditions .

[0086] Step S402: Adaptively adjust the edge-shading window according to the sea state level. Optimize trigger interval Wave constraint variance The adjustment rule is as follows: when the sea state is determined to be calm, the window is edged out. Set to 120 seconds to optimize the trigger interval. Set to 1.0s, wave constraint variance Set to 0.25 When the sea state is determined to be moderate, the marginal window... Set to 60 seconds to optimize the trigger interval. Set to 1.0s, wave constraint variance Set to 1.0 When the sea state is determined to be severe, the marginal window is adjusted. Set to 30 seconds, optimize trigger interval. Set to 2.0s, wave constraint variance Set to 4.0 As shown in Table 1.

[0087] Table 1: Adaptive Adjustment Rules for Sea State Levels

[0088]

[0089] Step S403: Solve the factor graph optimization problem using the iSAM2 incremental smoothing algorithm. The complete optimization objective function is as follows: for:

[0090]

[0091] in, This represents the prior factors of the initial state. These are the initial state estimates. Let be the initial state covariance. This represents the set of all state nodes in the factor graph. Whenever a new GNSS observation arrives or the trigger interval is optimized... Upon expiration, the new factor is added to the factor graph, and the iSAM2 incremental smoothing algorithm is invoked to perform an incremental update. For factors outside the marginal window... Historical state nodes will be subject to fixed-latency edgedification, which will exceed the edgedification window. Historical state nodes are removed from the optimization graph and replaced with prior factors to control the size of the factor graph and computational complexity.

[0092] Step S404: During GNSS signal interruption, short-term position propagation is performed using IMU pre-integration, with auxiliary constraints provided by ocean dynamics constraint factors. The maximum duration of IMU autonomous estimation is set according to sea state level: 15s for calm sea states, 10s for moderate sea states, and 5s for severe sea states. If the maximum estimation time is exceeded, the current positioning result is marked as low confidence.

[0093] Specifically, step S5 includes the following:

[0094] Step S501: When there is ranging information between adjacent buoys in the buoy network, acquire the ranging information between adjacent buoys. The ranging value between adjacent buoy A and buoy B is denoted as... Ranging information is obtained through acoustic ranging equipment between buoys or by estimating radio signal strength. The acoustic ranging accuracy is 0.1 to 1.0 m, and the radio ranging accuracy is 1 to 5 m. Communication between buoys uses time-division multiple access, with a communication cycle of 1 to 5 seconds.

[0095] Step S502: Construct multi-buoy cooperative distance constraint factors The formula is:

[0096]

[0097] in, and These represent the estimated positions of buoy A and buoy B at time k, respectively. This represents the variance of the distance measurement error.

[0098] Step S503: The cooperative distance constraint factor is added to the factor graph of each buoy to participate in joint optimization. When the GNSS signal of a buoy is completely interrupted, the positioning information of neighboring buoys is transmitted through the cooperative distance constraint factor to provide external position constraints for the interrupted buoy. The triggering condition for cooperative correction is: the GNSS interruption of this buoy exceeds 3 seconds, and the GNSS positioning of at least one neighboring buoy is available.

[0099] Step S6: Perform abnormal drift detection and factor graph reconstruction on the buoy positioning result sequence.

[0100] Step S601: Calculate the Mahalanobis distance between the current estimated position and the ocean dynamics predicted position. The formula is:

[0101]

[0102] in, This indicates the current position estimate. This indicates the position based on the previous time-instance state and the ocean current model prediction. The covariance matrix represents the predicted location.

[0103] Step S602, when the Mahalanobis distance Exceeding the threshold If this occurs, the current location result is determined to be abnormal. Threshold Based on the chi-square distribution, the 95th quantile with 3 degrees of freedom is taken, i.e. .

[0104] Step S603: After detecting an anomaly, perform local reconstruction of the factor graph. First, mark the time before and after the anomaly. The state node within the last second is the node to be reconstructed. The value is between 5 and 10 seconds. The node to be reconstructed is removed from the current factor graph. The local factor graph is reconstructed using the nearest reliable state node as the new prior constraint, and local optimization is performed. After reconstruction, the local optimization results are reintegrated into the global factor graph.

[0105] Step S604: Record the anomaly detection results and factor graph reconstruction history as system health status indicators. If anomaly detection is triggered more than 5 times consecutively, the system will issue a sensor check alarm.

[0106] Step S7: Conduct experimental verification and evaluate beneficial effects, specifically including:

[0107] Step S701: Configure the experimental verification environment. Data collection was conducted at a buoy observation station in a selected sea area, covering both calm and moderate sea states. The buoy was equipped with a u-blox F9P GNSS receiver with a sampling rate of 1Hz and an ICM-20948 MEMS inertial measurement unit with a sampling rate of 200Hz. The buoy antenna was positioned 2.5m above the sea surface. The reference ground truth was provided by a high-precision GNSS / INS integrated navigation system mounted on the buoy, with a horizontal positioning accuracy better than 0.1m. The total experimental data duration was 1200s.

[0108] Step S702: Set the comparison method and evaluation indicators. Four comparison schemes are set: Method A is a GNSS single-point positioning scheme, using only the position solution output by the GNSS receiver; Method B is a standard FGO scheme, employing loosely coupled fusion of GNSS position factors and IMU pre-integration factors, with a fixed marginalization window set to 60s, and no ocean dynamic constraints; Method C is a real-time FGO scheme, enabling IMU autonomous calculation and iSAM2 incremental optimization based on Method B, with a fixed marginalization window set to 60s; Method D is the multi-source fusion positioning method provided by this invention. The evaluation indicators include horizontal root mean square error, three-dimensional root mean square error, and positioning service availability. And the computation time for a single optimization.

[0109] Step S703: Comparison of experimental results. Method A has a horizontal root mean square error (RMSE) of 8.67m, a three-dimensional root mean square error (RMSE) of 11.36m, and a positioning service availability of 72.7%. Method B has a horizontal RMSE of 3.85m, a three-dimensional RMSE of 5.31m, a positioning service availability of 72.7%, and a single optimization calculation time of 12.3ms. Method C has a horizontal RMSE of 6.83m, a three-dimensional RMSE of 7.93m, a positioning service availability of 100.0%, and a single optimization calculation time of 5.8ms. Method D has a horizontal RMSE of 4.71m, a three-dimensional RMSE of 5.24m, a positioning service availability of 100.0%, and a single optimization calculation time of 4.5ms. Method D's horizontal RMSE is 45.7% lower than Method A and 31.0% lower than Method C. Method D achieves positioning accuracy close to that of batch processing method B while maintaining 100.0% service availability, as shown in Table 2.

[0110] Table 2: Comparison of Positioning Accuracy and Service Availability of Different Methods

[0111]

[0112] Step S704: Verify the beneficial effects of the sea state adaptive edge-mapping strategy. Under a 60s edge-mapping window, the horizontal root mean square error of sea state adaptive edge-mapping is 4.3m, a 36.8% reduction compared to the fixed edge-mapping's 6.8m; the single-cycle optimization time of sea state adaptive edge-mapping is 4.5ms, a 25.0% reduction compared to the fixed edge-mapping's 6.0ms. Under severe sea conditions, the single-cycle optimization calculation time of sea state adaptive edge-mapping is reduced by more than 50% compared to the fixed window method, meeting the real-time requirements of limited buoy computing resources; under calm sea conditions, the positioning accuracy of sea state adaptive edge-mapping is improved by more than 15% compared to the short-window scheme.

[0113] Step S705: Verify the beneficial effects of the ocean dynamics constraint factor and multi-buoy cooperative positioning correction. During GNSS signal interruption, the ocean dynamics constraint factor provides physical constraints on wave motion and ocean current drift for the buoy, reducing the position drift rate during IMU autonomous estimation by more than 40%. In a 10s GNSS interruption scenario, the horizontal positioning error is reduced from more than 8m in existing methods to less than 5m. At an error threshold of 10m, the positioning service availability of method D reaches 92%, method C reaches 78%, and method A reaches 55%. When single-buoy GNSS is completely interrupted, through multi-buoy cooperative positioning correction, the positioning service availability is improved from less than 60% in existing methods to more than 85%.

[0114] Step S706: Verify the beneficial effects of the sea surface multipath adaptive weighting mechanism and anomaly detection. The sea surface multipath adaptive weighting mechanism dynamically adjusts the weights of each satellite observation using a sea surface reflection geometry model. Under moderate sea state conditions, the horizontal positioning accuracy is improved by more than 25% compared to the fixed weighting method. Based on the anomaly drift detection and factor map reconstruction mechanism using ocean dynamics priors, the recognition rate of anomaly positioning results reaches over 95%. The overall technical solution of this invention improves the comprehensive positioning accuracy of buoys under moderate sea state conditions from over 10m in existing methods to within 5m, and the availability of positioning services from below 50% to over 85%.

[0115] For the GNSS receiver weighted least squares solution method, time division multiple access communication control logic, Mahalanobis distance and chi-square distribution numerical calculation model, root mean square error statistical calculation formula, etc. involved in the embodiments of the present invention, those skilled in the art can refer to surveying and navigation standards and basic mathematics manuals to implement them. The relevant operation process and underlying calculation principle are well known technologies in this field and will not be described in detail here.

[0116] Although embodiments of the invention have been shown and described, it will be understood by those skilled in the art that various changes, modifications, substitutions and alterations can be made to these embodiments without departing from the principles and spirit of the invention, the scope of which is defined by the appended claims and their equivalents.

Claims

1. A multi-source fusion positioning method for ocean buoys based on factor graph optimization, characterized in that, Includes the following steps: Step S1: Perform multi-source observation data acquisition and preprocessing to obtain GNSS observation data, IMU measurement data, and marine environmental auxiliary data; The marine environmental auxiliary data includes tidal forecast height, sea state level parameters, and ocean current forecast speed in the sea area where the buoy is located; Step S2: Based on the preprocessed multi-source data, construct the factor graph state nodes and multi-source factors; Construct state nodes of the factor graph. The navigation state node at each time step includes the rotation matrix from the vehicle coordinate system to the navigation coordinate system, the three-dimensional position vector, the three-dimensional velocity vector, the accelerometer deviation, and the gyroscope deviation. Construct GNSS location factors. When GNSS observations are available, use the Mahalanobis distance as the residual of the GNSS factor by using the GNSS covariance matrix after adaptive weighting for sea surface multipath effects. An IMU pre-integration factor is constructed, which pre-integrates all IMU measurements between adjacent GNSS observation times into a single relative motion constraint. The IMU bias is modeled as a slowly varying random walk process. The constraint factor of the IMU bias is constructed using the random walk covariance of accelerometer bias and the random walk covariance of gyroscope bias. Ocean dynamics constraint factors are constructed, among which the wave motion constraint factor is based on linear wave theory and is constructed using the local mean sea level height, the predicted tidal height, and the sea surface undulation caused by waves, constraining the vertical displacement of the buoy to be near the sea surface; the ocean current drift constraint factor is constructed using the predicted ocean current velocity vector and the ocean current constraint covariance matrix, constraining the horizontal motion velocity of the buoy to be near the ocean current velocity. Step S3: Perform adaptive weighting of the GNSS observation data for sea surface multipath effects based on sea state parameters; Step S4: Obtain sea state parameters and perform adaptive incremental optimization solution based on the sea state parameters; The intensity index of buoy motion is calculated based on the square mean of the difference between the magnitude of the triaxial accelerometer measurement and the scalar value of gravitational acceleration within the evaluation window. The sea state is classified according to the buoy motion intensity index, and the marginalization window duration, trigger interval and wave constraint variance parameters are adaptively and dynamically adjusted according to the sea state classification results. Step S5: In the presence of multiple buoys, use the ranging information between adjacent buoys to perform multi-buoy cooperative positioning correction; Step S6: The anomaly detection module performs anomaly drift detection and factor graph reconstruction on the positioning sequence and outputs the ocean buoy positioning results.

2. The multi-source fusion positioning method for ocean buoys based on factor graph optimization according to claim 1, characterized in that, Step S1 includes: The GNSS receiver on the buoy acquires satellite observation data, including pseudorange observations, Doppler observations, carrier-to-noise ratios and elevation angles of each satellite, and obtains the position solution and the corresponding covariance matrix by weighted least squares method. The buoy acquires IMU measurement data. The MEMS inertial measurement unit on board the buoy acquires triaxial accelerometer and triaxial gyroscope measurements and establishes an IMU measurement model that includes linear acceleration, gravitational acceleration vector, true angular velocity, accelerometer bias, gyroscope bias, and white noise in the navigation coordinate system. Data time reference unification and signal quality assessment: all observation data timestamps are aligned with GPS time as the reference; satellite observations are divided into three levels—high quality, medium quality, and low quality—based on the carrier-to-noise ratio, and low-quality observations are marked as suspected multipath pollution on the sea surface.

3. The multi-source fusion positioning method for ocean buoys based on factor graph optimization according to claim 1, characterized in that, Step S3 involves calculating the sea surface multipath impact weights for each visible satellite, including: The glancing angle of the reflection from the sea surface is determined by the satellite elevation angle and the height of the buoy antenna above the sea surface; The sea surface reflection coefficient is calculated based on the Fresnel reflection coefficient of a smooth sea surface, the root mean square wave height of the sea surface, and the signal carrier wavelength. Based on the sea surface reflection coefficient, multipath penalty coefficient, and the ratio of the observed carrier-to-noise ratio to the reference carrier-to-noise ratio threshold, a multipath weighting factor for each satellite is constructed.

4. The multi-source fusion positioning method for ocean buoys based on factor graph optimization according to claim 3, characterized in that, Step S3 further includes: Construct a diagonal matrix using the multipath weighting factors corresponding to the three-dimensional coordinate axes; Multiply the unadjusted covariance matrix with the diagonal matrix to obtain the weighted covariance matrix, which is used to construct the GNSS factor in step S2. When the carrier-to-noise ratio of a satellite is lower than the preset carrier-to-noise ratio value and the reflection coefficient is greater than the preset reflection coefficient value, the observations of the corresponding satellite will be removed from the solution.

5. The multi-source fusion positioning method for ocean buoys based on factor graph optimization according to claim 1, characterized in that, Step S4 includes: The iSAM2 incremental smoothing algorithm is used to solve the optimization objective function, which includes the set of initial state prior factors, IMU pre-integration factors, bias constraint factors, GNSS position factors, wave motion constraint factors, and ocean current drift constraint factors. Whenever a new observation arrives or the optimization trigger interval expires, an incremental update is performed. Historical state nodes that exceed the marginalization window are marginalized with fixed hysteresis, removed, and replaced with a priori factors. During signal interruption, short-term position propagation is performed using pre-integration, and auxiliary constraints are provided in conjunction with the ocean dynamics constraint factor. The maximum duration is set autonomously based on the sea state level. If the maximum duration is exceeded, the positioning result is marked as low confidence.

6. The multi-source fusion positioning method for ocean buoys based on factor graph optimization according to claim 1, characterized in that, Step S5 includes: To obtain ranging information between adjacent buoys communicating using time division multiple access in a buoy network; By using the distance difference between buoy position estimates and the ranging information, and combining the ranging error variance, a multi-buoy cooperative distance constraint factor is constructed. If the signal of this buoy is completely interrupted for more than a preset time and the signal of at least one neighboring buoy is available, the cooperative correction trigger condition is met, and the cooperative distance constraint factor is added to the factor graph of each buoy to participate in joint optimization.

7. The multi-source fusion positioning method for ocean buoys based on factor graph optimization according to claim 1, characterized in that, Step S6 includes: Calculate the Mahalanobis distance between the current position estimate and the position predicted based on the previous time-instance state and the ocean current model; When the Mahalanobis distance exceeds the threshold set based on the chi-square distribution, the current positioning result is determined to be abnormal; Upon detecting an anomaly, a local reconstruction of the factor graph is performed. The state nodes before and after the anomaly are marked as nodes to be reconstructed and removed from the current factor graph. The local factor graph is reconstructed using the nearest reliable state node as the new prior constraint, and local optimization is performed. After reconstruction, the local factor graph is reconnected to the global factor graph. Record the anomaly detection results and factor graph reconstruction history as output indicators of system health status.

8. A multi-source fusion positioning system for ocean buoys based on factor graph optimization, employing the multi-source fusion positioning method for ocean buoys based on factor graph optimization as described in any one of claims 1-7, characterized in that, include: The data acquisition module is used to acquire GNSS observation data, IMU measurement data, and marine environmental auxiliary data, and to perform data acquisition and preprocessing operations. The factor graph construction module is used to construct navigation status nodes and multi-source measurement factors, and to perform factor graph construction and sea surface multipath adaptive weighting operations. The adaptive optimization module is used to perform adaptive incremental optimization solutions based on sea state levels, including sea state discrimination, marginal window adjustment, and iSAM2 optimization operations. The collaborative correction module is used to perform collaborative positioning correction using distance measurement information between multiple buoys, and to perform collaborative constraint construction and joint optimization operations. The anomaly detection module is used to detect anomaly localization results and perform local reconstruction of the factor graph, performing anomaly detection and reconstruction operations.

Images