GNSS-acoustic seabed positioning data weighting method based on spatiotemporal covariance adaptive modeling

CN122883227APending Publication Date: 2026-10-09HARBIN INSTITUTE OF TECHNOLOGY (SHENZHEN) (INSTITUTE OF SCIENCE AND TECHNOLOGY INNOVATION HARBIN INSTITUTE OF TECHNOLOGY SHENZHEN)
View PDF 0 Cites 0 Cited by

Patent Information

Application Number
CN202611381272.X
Authority / Receiving Office
CN · China
Patent Type
Applications(China)
Current Assignee / Owner
Filing Date
2026-09-07
Publication Date
2026-10-09

AI Technical Summary

Technical Problem

[0005]现有技术的核心缺陷在于:观测残差被假设为独立同分布,即协方差矩阵为乘以单位矩阵

Benefits of technology

本发明针对传统GNSS-声学海底定位方法中观测误差时空相关性被忽视、密集采样数据过度支配解算、不同应答器误差被平均化以及经验定权主观性强等缺陷,提出一种基于时空协方差自适应建模的数据定权方法。

✦ Generated by Eureka AI based on patent content.

Smart Images

  • Figure CN122883227A_ABST
    Figure CN122883227A_ABST
Patent Text Reader

Abstract

The application discloses a GNSS-acoustic seabed positioning data weighting method based on space-time covariance adaptive modeling. In view of the problems of traditional methods, such as neglecting space-time correlation of observation errors, overfitting of dense sampling and strong subjectivity of empirical weighting, the application firstly acquires multi-source observation data and pre-processes; then, a double-parameter space-time covariance matrix containing a time correlation length and a cross-transponder correlation coefficient is constructed, and a non-diagonal element is dynamically adjusted according to a time sampling density and a transponder identification, so that adaptive weighting of dense data and sparse data is realized; then, covariance parameters and a sound speed field smoothness hyperparameter are brought into a Akaike information criterion framework for synchronous and objective optimization; finally, the generalized least square is solved by substituting the coordinates of the seabed reference point. The application can accurately depict the space-time structure of errors, significantly improve the positioning accuracy and spatial consistency, and is completely adaptive and does not need manual intervention, and is suitable for different sea areas and observation platforms.
Need to check novelty before this filing date? Find Prior Art

Description

Technical Field

[0001] This invention relates to marine geodesy and data processing technology, and in particular to a GNSS-acoustic (Global Navigation Satellite System and Acoustic Joint) seabed positioning data weighting method based on spatiotemporal covariance adaptive modeling. Background Technology

[0002] Traditional seabed positioning relies on shipborne Global Navigation Satellite Systems (GNSS) to determine the position of transducers on the sea surface, using acoustic ranging to obtain the propagation time from the transducer to the seabed transponder, and then jointly calculating the transponder coordinates. However, factors such as the spatiotemporal heterogeneity of the seawater sound velocity field, the propagation characteristics of transducer positioning errors, and the differences in the sound ray paths of different transponders severely limit further improvements in positioning accuracy.

[0003] In practical operations, high-precision GNSS receivers (sampling rates up to 10-20 Hz), fiber optic or microelectromechanical inertial measurement units (IMUs, attitude measurement accuracy better than 0.01°), multibeam echo sounders (such as the SeaBat T50-P, with up to 512 beams and a coverage angle of up to 140°), acoustic transducers (center frequency 10-30 kHz), sound velocity profilers (SVP / CTD), and remotely operated vehicles (ROVs) typically work in concert. The surface platform (survey vessel or buoy) continuously acquires GNSS data at frequencies of several hertz, while the acoustic transducers emit sound waves at intervals of several seconds and receive echoes from the transponders. This high-frequency sampling method generates a dense sequence of observations in time, with strong correlations in the errors between adjacent observations. Simultaneously, the acoustic paths of different transponders gradually separate at deeper water levels, making the observation error structure of each transponder independent.

[0004] Overview of technical principles: The observation equation for GNSS-acoustic seabed positioning can be expressed as: in For the first The two-way propagation time of the second observation For a depth-dependent function of sound speed, The spatial position of the transducer at the time of transmission (determined by GNSS+IMU). The coordinates of the seabed transponder are to be determined. To observe the residuals.

[0005] The core flaw of existing technology lies in: observation residuals. They are assumed to be independent and identically distributed, i.e., their covariance matrix is... Multiply by the identity matrix. However, in reality, transducer positioning errors propagate continuously over time, resulting in non-zero covariance between observation residuals at adjacent time points; differences in the acoustic paths of different transponders lead to varying error variances. This invention constructs a spatiotemporal covariance matrix to accurately characterize the aforementioned correlation structure, achieving objective and adaptive determination of data weights.

[0006] GNSS-acoustic seabed positioning is a key technology for constructing seabed geodetic benchmarks. Shipborne GNSS on the sea surface determines the transducer's location, the acoustic transducer transmits sound waves to a seabed transponder and records the round-trip propagation time, and the transponder's coordinates are obtained by solving the propagation time equation. In existing technologies, the weighting methods for observation data are mainly divided into two categories: equal-weighting and empirical weighting. The specific implementations of these two methods in the GNSS-acoustic field are as follows: (1) Equal weighting treatment (equal weighting method) In the GNSS-acoustic field, equal-weighting refers to assigning the same weights to all acoustic observations, i.e., the weight matrix is ​​the identity matrix. Its theoretical basis is the assumption that all observations have the same prior variance. In practice: the order of observation times is disregarded, and all propagation time observations are treated equally; transponder numbers are disregarded, and data from all transponders are mixed for calculation; the solution uses ordinary least squares, with the objective function being the minimization of the sum of squared residuals. The advantages of this method are its simplicity and low computational cost, but its disadvantages are significant: the large number of observations in densely sampled areas (such as when the ship is slow or stationary) can over-dominate the solution results, leading to an underestimation of the variance in the coordinate solution and resulting in overfitting.

[0007] (2) Judgment based on experience (judgment based on experience) Empirical weighting is based on prior geometric information or physical quantities to construct weight functions. Common formulas include: inverse distance weighting (weights are inversely proportional to the square of the distance), sine elevation weighting (weights are inversely proportional to the square of the sine of the beam incidence angle), and combined weighting. In the GNSS-acoustic field, empirical weighting is often combined with iterative weighted least squares. The specific process is as follows: first, obtain the initial coordinates using the equal weighting method; calculate the distance and elevation angle for each observation based on the initial coordinates; calculate the weights according to the empirical formula, and recalculate; iterate until convergence. Experimental studies show that empirical weighting can reduce the root mean square error of coordinates by 10% to 20% compared to the equal weighting method, but it still cannot handle time correlation and the error independence of different transponders.

[0008] like Figure 1 As shown, the traditional GNSS-acoustic data processing workflow can be broken down into the following steps: Data Acquisition and Preprocessing: The shipborne GNSS receiver acquires dual-frequency carrier phase and pseudorange data at a sampling rate of 1–10 Hz. PPK (Dynamic Post-Processing Differential) or PPP (Precise Point Positioning) is used to calculate centimeter-level transducer position timing. The acoustic transducer emits acoustic pulses at a fixed period (e.g., 2-second intervals), and the transmission and reception times are recorded to calculate the two-way propagation time. The IMU measures the ship's roll, pitch, and heading for lever arm effect correction: the transducer position equals the GNSS antenna position plus the rotation matrix multiplied by the lever arm vector. A sound velocity profiler (SVP / CTD) projects layered sound velocities for ray bending correction.

[0009] Establishment of the observation equation: For the k-th observation, the calculation of the theoretical propagation time involves the equivalent average sound speed and geometric distance. After linearization, the error equation is obtained, which expresses the difference between the observed propagation time and the initial calculated value as a linear function of the increment of the unknown parameters.

[0010] Weight matrix construction (equal weight / empirical weight): Equal weight method: The weight matrix is ​​an identity matrix. Empirical weight method: The weight matrix is ​​a diagonal matrix, and the diagonal elements are calculated based on distance or pitch angle, commonly the square of the reciprocal of the sine of the pitch angle.

[0011] The normal equation is constructed by multiplying the normal matrix by the increment of the unknown parameters, which equals the constant term. The normal matrix is ​​the transpose of the design matrix multiplied by the weight matrix and then by the design matrix. The constant term is the transpose of the design matrix multiplied by the weight matrix and then by the vector of observed values ​​minus calculated values.

[0012] Parameter solving (least squares): The increment of the unknown parameter is equal to the inverse of the normal matrix multiplied by a constant term, and the parameter is updated. It usually converges after 2 to 3 iterations.

[0013] Accuracy assessment: The unit weighted variance is the weighted sum of squared residuals divided by the degrees of freedom, and the parametric covariance matrix is ​​the unit weighted variance multiplied by the inverse of the normal matrix.

[0014] Sonic field smoothing constraint (optional): If the acoustic field is solved jointly, a smoothing constraint is introduced, and the normal equation becomes the original normal equation plus the smoothing constraint matrix multiplied by the smoothing parameter.

[0015] Output results: Output the coordinates and standard deviation of the seabed transponder, as well as the residual sequence for quality check.

[0016] The main drawbacks of existing technologies are: 1. Ignored Correlation: Transducer positioning errors propagate continuously over time, resulting in strong temporal correlation between observations at adjacent moments. Different transponders, situated along different sound ray paths, exhibit varying error structures due to the spatiotemporal heterogeneity of sound speeds. Ignoring this correlation leads to excessive dominance of the solution results by densely sampled data, resulting in overfitting and a significant underestimation of the variance of the calculated coordinates.

[0017] 2. Reduced spatial resolution: The differences in error characteristics among different transponders are averaged out, weakening the spatial consistency of the solution results. For large-scale multi-transponder arrays, the relative accuracy between reference points deteriorates, failing to meet millimeter-level relative positioning requirements.

[0018] 3. High subjectivity: The weight function of empirically determined weights relies on manually set empirical parameters and lacks objective adaptive adjustment capabilities. When the sound velocity profile changes significantly in different sea areas and seasons, the performance of the same empirical formula fluctuates greatly and its universality is poor.

[0019] It should be noted that the information disclosed in the background section above is only for understanding the background of this application, and therefore may include information that does not constitute prior art known to those skilled in the art. Summary of the Invention

[0020] The main objective of this invention is to overcome the deficiencies in the aforementioned background technology and provide a GNSS-acoustic seabed positioning data weighting method based on spatiotemporal covariance adaptive modeling.

[0021] To achieve the above objectives, the present invention adopts the following technical solution: A weighting method for GNSS-acoustic seabed positioning data based on spatiotemporal covariance adaptive modeling includes the following steps: S1. Acquire multi-source observation data and perform preprocessing. The multi-source observation data includes high-precision GNSS data, acoustic transducer data, multibeam bathymetry data, and ROV underwater robot data. S2. Construct a two-parameter spatiotemporal covariance matrix to characterize the temporal correlation of observation errors and the independence of errors of different transponders. Based on the temporal sampling density of the observation data and the transponder identifier, adaptive weighting is achieved by automatically reducing the weight of densely sampled data and automatically increasing the weight of sparse data through the off-diagonal elements of the covariance matrix. S3. Incorporate the covariance parameter and the sound velocity field smoothness hyperparameter into the Akaike Bayesian Information Criterion (ABIC) framework, and achieve synchronous objective optimization of all hyperparameters through a multi-dimensional optimization strategy. S4. Substitute the optimized covariance matrix and the optimal smoothness parameter into the generalized least squares solver to solve for the coordinates of the seabed reference point and output the positioning results.

[0022] A computer program product includes a computer program that, when executed by a processor, implements the aforementioned GNSS-acoustic seabed positioning data weighting method based on spatiotemporal covariance adaptive modeling.

[0023] The present invention has the following beneficial effects: This invention addresses the shortcomings of traditional GNSS-acoustic seabed positioning methods, such as neglecting the spatiotemporal correlation of observation errors, excessive dominance of dense sampling data in the solution, averaging of errors from different transponders, and strong subjectivity in empirical weighting. It proposes a data weighting method based on adaptive modeling of spatiotemporal covariance.

[0024] The core of this invention lies in constructing a two-parameter covariance matrix containing time correlation length and cross-responder correlation coefficient, and incorporating the covariance parameter and model smoothness parameter into the Akaike Bayesian Information Criterion (ABIC) framework for synchronous optimization, thereby achieving objective and adaptive determination of data weights. Through this two-parameter spatiotemporal covariance matrix, this invention can accurately characterize the continuous propagation characteristics of observation errors at adjacent time points and the independent error structure of different transponders due to differences in sound path, thus accurately reflecting the true distribution characteristics of observation noise.

[0025] Based on this, the off-diagonal elements of the covariance matrix can be dynamically adjusted according to the temporal sampling density and transponder identifiers of the observation data: in densely sampled regions, the time interval between adjacent observations is small, resulting in a larger covariance value. In generalized least squares solutions, the equivalent weights are automatically reduced, effectively suppressing overfitting. Sparse data, on the other hand, receive relatively higher weights, achieving data-driven adaptive weighting. Simultaneously, the covariance parameter and the sound velocity field smoothness hyperparameter are incorporated into the ABIC framework for synchronous and objective optimization, completely avoiding the subjectivity and uncertainty of manual parameter tuning. This allows the method to adapt to changes in error characteristics across different sea areas, seasons, and observation platforms.

[0026] Compared with existing technologies, this invention significantly improves the accuracy of seabed positioning. Simulation and measured data verification show that the root mean square error of the solution is reduced by 30% to 40% compared with the equal weight method and by 15% to 20% compared with the empirical weighting method. Spatial consistency is greatly improved, the baseline length repeatability standard deviation of the multi-transponder array solution is improved by about 40%, and the relative accuracy between each reference point is significantly improved. Moreover, the method is fully automated and does not require human experience intervention, which has strong universality and engineering practical value.

[0027] Other beneficial effects of the embodiments of the present invention will be further described below. Attached Figure Description

[0028] Figure 1 This is a schematic diagram of the existing GNSS-acoustic seabed positioning data processing flow.

[0029] Figure 2 This is the overall flowchart of the GNSS-acoustic seabed positioning data weighting method based on spatiotemporal covariance adaptive modeling of the present invention.

[0030] Figure 3This is a detailed flowchart of the GNSS-acoustic seabed positioning data weighting method according to an embodiment of the present invention. Detailed Implementation

[0031] The embodiments of the present invention will be described in detail below. It should be emphasized that the following description is merely exemplary and not intended to limit the scope and application of the present invention.

[0032] In traditional GNSS-acoustic seabed positioning methods, equal-weighting or empirical weighting methods often ignore the spatiotemporal correlation of observation errors, leading to excessive dominance of densely sampled data in the solution, averaging of errors from different transponders, and a decrease in overall solution stability. To address these issues, this invention aims to construct a two-parameter covariance matrix that can simultaneously characterize temporal and cross-transponder correlations to accurately reflect the true distribution characteristics of observation noise. It adaptively adjusts the off-diagonal elements of the covariance matrix based on data sampling density and ray path differences, achieving automatic weight reduction for densely sampled data and reasonable weighting for sparse data. The covariance parameter and sound velocity field smoothness parameter are incorporated into a unified Akaike Bayesian information criterion optimization framework to achieve synchronous and objective optimization of hyperparameters. Furthermore, while ensuring solution stability, spatiotemporal covariance modeling improves the independent error modeling capability of different transponders, ultimately enhancing the accuracy and reliability of seabed reference point coordinate calculation. Therefore, this invention proposes a GNSS-acoustic seabed positioning data weighting method based on adaptive spatiotemporal covariance modeling.

[0033] See Figure 2 and Figure 3 This invention provides a method for weighting GNSS-acoustic seabed positioning data based on spatiotemporal covariance adaptive modeling, comprising the following steps: Step S1: Acquire multi-source observation data and perform preprocessing. The multi-source observation data includes high-precision GNSS data, acoustic transducer data, multibeam bathymetry data, and ROV underwater robot data.

[0034] In some embodiments, the data acquisition and preprocessing in step S1 includes: all devices synchronizing with microsecond-level clocks using GNSS second pulses and UTC time; performing dynamic post-processing differential calculations on GNSS data to obtain centimeter-level transducer position timing, and adding the rod arm vector after rotation matrix transformation to the GNSS antenna position coordinates to obtain the transducer position coordinates; performing ray bending correction on acoustic data, correcting the observed round-trip propagation time to the theoretical propagation time along the actual ray path; removing outliers from multibeam data; and performing time alignment and underwater trajectory gross error detection on ROV data to generate standardized observation values ​​and initial residual values.

[0035] Step S2: Construct a two-parameter spatiotemporal covariance matrix to characterize the temporal correlation of observation errors and the independence of errors of different transponders. Based on the temporal sampling density of the observation data and the transponder identifier, adaptive weighting is achieved by automatically reducing the weight of densely sampled data and automatically increasing the weight of sparse data through the off-diagonal elements of the covariance matrix.

[0036] In some embodiments, the process of constructing the two-parameter spatiotemporal covariance matrix in step S2 includes: constructing a time difference matrix based on the sampling time of each observation, where each element of the matrix is ​​the absolute value of the difference between two corresponding observation times; simultaneously, constructing a transponder identification matrix based on the transponder number corresponding to each observation, where the identification element is 0 when two observations belong to the same transponder and 1 when they belong to different transponders; dividing each element of the time difference matrix by the time correlation length parameter and taking the negative, then performing an exponential decay operation with the natural constant as the base to obtain the time decay factor matrix; performing a power operation with the cross-transponder correlation coefficient parameter as the base and each element of the transponder identification matrix as the exponent to obtain the cross-transponder correlation factor matrix; multiplying the time decay factor matrix and the cross-transponder correlation factor matrix element by element, and then multiplying by the variance scaling parameter to generate the two-parameter spatiotemporal covariance matrix, thereby realizing the joint modeling of time decay and the independence between transponders.

[0037] In some embodiments, the adaptive weighting is achieved through the following mechanism: In densely sampled regions, the time interval between adjacent observations is small, the corresponding element in the time decay factor matrix approaches 1, and the corresponding element in the covariance matrix is ​​large. In generalized least squares calculation, the weight of the observation is determined by the inverse of the covariance matrix; when the covariance value is large, the corresponding equivalent weight automatically decreases. In sparsely sampled data, the time interval between observations is large, the corresponding element in the time decay factor matrix approaches 0, and the corresponding element in the covariance matrix is ​​small, so the corresponding equivalent weight automatically increases. Observations from different transponders are independently modeled after exponentiation of the cross-transponder correlation coefficient parameter. When this parameter is small, the covariance between observations from different transponders approaches 0, achieving decoupling of the error structure.

[0038] Step S3: Incorporate the covariance parameter and the sound velocity field smoothness hyperparameter into the Akaike Bayesian Information Criterion (ABIC) framework, and achieve synchronous objective optimization of all hyperparameters through a multi-dimensional optimization strategy; In some embodiments, the joint optimization process of the Akaike Bayesian information criterion framework in step S3 includes: constructing an objective function of the Akaike Bayesian information criterion using the covariance parameter and the sound velocity field smoothness hyperparameter as variables to be optimized. This objective function is composed of the following sums: the number of observations minus the number of unknown parameters multiplied by the logarithm of the weighted residual variance; the logarithm of the determinant of the matrix obtained by multiplying the transpose of the design matrix by the inverse of the covariance matrix, then by the design matrix, plus the smoothness hyperparameter multiplied by the smoothing constraint matrix; and the logarithm of the determinant of the covariance matrix. The algorithm calculates the number of parameters and constant terms. Within the parameter domain, it performs multidimensional optimization using grid search or gradient descent. For each candidate parameter combination, it sequentially calculates the time difference matrix, constructs the covariance matrix, calculates the weighted residual variance and the above terms, and obtains the Akaike Bayes information criterion value. It iterates through all candidate parameter combinations or along the gradient descent direction, selecting the parameter combination that minimizes the Akaike Bayes information criterion value as the optimal hyperparameter, thus achieving the synchronous and objective determination of the covariance parameter and the smoothness hyperparameter. The iteration terminates when the change in the objective function is less than the preset convergence threshold.

[0039] In some embodiments, the process of achieving synchronous objective optimization of hyperparameters through a multi-dimensional optimization strategy in step S3 further includes: The variance scaling parameter does not participate in the ABIC iterative optimization. Instead, it is directly calculated analytically given the covariance parameter and the smoothness hyperparameter. A structure matrix is ​​defined such that the covariance matrix equals the variance scaling parameter multiplied by the structure matrix. The elements of the structure matrix are determined by the product of the power operation of the cross-responder correlation coefficient and the time decay exponent. In each ABIC iteration, the structure matrix is ​​constructed given the current hyperparameter combination, and the generalized least squares method is used to obtain the estimated values ​​of the unknown parameters. Then, the residual vector is calculated based on the difference between the observed value vector and the design matrix multiplied by the estimated value of the unknown parameters. The transpose of the residual vector is multiplied by the inverse of the structure matrix, then multiplied by the residual vector again, and divided by the total number of observations to obtain the estimated value of the variance scaling parameter. Finally, the variance scaling parameter is substituted into the covariance matrix for the calculation of the next ABIC expression. After substituting the analytical expression of the variance scaling parameter into the ABIC expression, ABIC is rewritten as a function that depends only on the covariance parameter and the smoothness hyperparameter, with the variance scaling parameter being automatically optimized as a latent variable. The covariance matrix sparsity acceleration mechanism includes: utilizing time decay sparsity, where the exponential decay term is less than a threshold when the time interval between two observations exceeds a preset multiple of the time correlation length, and the correlation can be ignored. Only elements within this time interval are calculated, reducing the storage complexity from the square of the total number of observations to the order of the total number of observations multiplied by the effective bandwidth; and utilizing a block diagonal structure, where the covariance between different transponders can be ignored when the cross-responder correlation coefficient is lower than a preset threshold, and the covariance matrix degenerates into a block diagonal matrix.

[0040] Step S4: Substitute the optimized covariance matrix and the optimal smoothness parameter into the generalized least squares solver to solve for the coordinates of the seabed reference point and output the positioning result.

[0041] In some embodiments, the generalized least squares solution in step S4 includes: multiplying the inverse of the optimal covariance matrix by the transpose of the design matrix and then by the observation vector to obtain a constant term; simultaneously, multiplying the transpose of the design matrix by the inverse of the covariance matrix and then by the design matrix, adding the smoothness hyperparameter and multiplying by the smoothing constraint matrix to construct the normal matrix; obtaining the optimal estimate of the seabed reference point coordinate increment by solving the inverse of the normal matrix and multiplying by the constant term, and updating the parameters after superimposing the increment onto the initial coordinates; repeating the above process until the parameter change is less than the convergence threshold; calculating the posterior covariance matrix based on the inverse of the normal matrix, and the square root of its diagonal elements is the posterior standard deviation of the coordinates, used for accuracy assessment.

[0042] In some embodiments, the positioning results output in step S4 include: seabed reference point coordinates and their posterior standard deviation, wherein the posterior standard deviation is obtained by taking the square root of the diagonal elements of the posterior covariance matrix; optimal covariance parameter estimates, including time correlation length parameters, transponder correlation coefficient parameters and smoothness hyperparameters; digital elevation model after multibeam terrain correction; and three-dimensional position sequence after ROV trajectory smoothing.

[0043] In some embodiments, the method further includes covariance expansion of the multibeam bathymetry data, the process of which includes: treating the bathymetry data of each beam of the multibeam bathymetry as an independent observation, constructing an extended covariance function for the observation residual correlation of adjacent beams within the same transmission cycle; the function first calculates the absolute value of the difference between the incident angles of each beam to obtain an angle difference matrix, and performs a power operation with the beam angle correlation coefficient as the base and each element of the angle difference matrix as the exponent to obtain an angle correlation factor; simultaneously, for the observations of different transmission cycles, the time decay factor is calculated by dividing the time difference by the time correlation length, taking the negative, and then exponentially decaying; the angle correlation factor is multiplied by the time decay factor, and then multiplied by the multibeam variance scaling parameter to generate a multibeam covariance matrix, thereby realizing the joint modeling of the spatial correlation caused by beam angle differences and the time correlation between different transmission cycles.

[0044] In some embodiments, the method further includes covariance expansion of the multi-source sensors of the ROV underwater robot, the process of which includes: constructing a time covariance function for continuous trajectory points output by the DVL and inertial navigation system, generating a time smoothing constraint matrix by dividing the time difference by the ROV time correlation length, taking the negative, and then exponentially decaying; introducing a cross-sensor correlation coefficient for discrete positioning points obtained by USBL, scaling the covariance matrix of the DVL-INS trajectory through this coefficient, and establishing a correlation covariance sub-matrix between USBL observations and continuous trajectories; constructing a spatial covariance function to constrain the consistency of feature point clouds obtained by hyperspectral cameras and image sonar; using the covariance sub-matrices corresponding to each sensor as the main diagonal block of the block diagonal matrix, and then calculating the cross covariance matrix as the off-diagonal block according to the calibration parameters and spatiotemporal alignment relationship between different sensors, and adding the block diagonal matrix and the cross covariance matrix to obtain the unified covariance matrix of ROV multi-source observations.

[0045] The proposed GNSS-acoustic seabed positioning data weighting method based on spatiotemporal covariance adaptive modeling has the following main technical advantages: By constructing a two-parameter spatiotemporal covariance matrix composed of time correlation length and cross-transponder correlation coefficient, it accurately characterizes the continuous temporal propagation characteristics of observation errors and the independent error structure of different transponders. This enables adaptive weighting that automatically adjusts weights according to data sampling density, effectively suppressing overfitting from dense sampling. Simultaneously, it incorporates the covariance parameter and sound velocity field smoothness parameter into the Akaike Bayesian information criterion framework for synchronous and objective optimization, completely eliminating the need for human experience intervention. Compared to existing methods, it significantly improves solution accuracy, greatly enhances spatial consistency, and possesses strong generalization capabilities and engineering practicality across sea areas, seasons, and platforms.

[0046] The following further describes specific embodiments, algorithm examples, and test results of the present invention.

[0047] A weighting method for GNSS-acoustic seabed positioning data based on spatiotemporal covariance adaptive modeling, the specific process of which is as follows: Figure 3 As shown. In this method, a time-related length is constructed. Correlation coefficient across transponders The two-parameter covariance matrix is ​​calculated, and the covariance parameter and model smoothness parameter are incorporated into the Akaike Bayesian Information Criterion (ABIC) framework for simultaneous optimization, achieving objective and adaptive determination of data weights. The main steps are as follows: 1. Construct the spatiotemporal covariance function Let the observed travel time residual vector be... Each observation corresponds to a specific sampling time and transponder number. A covariance matrix C is constructed, whose elements are determined by both temporal correlation and cross-transponder correlation: in: This is the variance scaling parameter (which can be estimated in conjunction with other hyperparameters); The time correlation length controls the decay rate of the correlation between observation errors at adjacent time points. The larger the value, the stronger the time correlation. The cross-responder correlation coefficient (values ​​from 0 to 1) is the coefficient when the observation... and When they belong to the same transponder ,otherwise ; This represents the absolute value of the time interval between two observations. This function has a clear physical meaning: the correlation between observations from the same transponder decreases exponentially with increasing time interval; the correlation between observations from different transponders decreases through... Independent modeling The smaller the value, the more independent the errors of different transponders.

[0048] Estimate: Variance scaling parameter Not participating in ABIC iterative optimization, but given Under the given conditions, it can be directly calculated using analytical formulas. Let the structure matrix be... satisfy ,in ,but The maximum likelihood estimate is: in This represents the total number of observations. In each ABIC iteration, the following steps are performed in order: ① Given the current hyperparameters... Constructing a structure matrix And solve the generalized least squares problem to obtain... ① Estimate the value; ② Calculate the residual vector ③ Calculate according to the above formula The estimated value; ④ Substitute into the covariance matrix This is used for the next evaluation of the ABIC expression. After substituting the parsing expression into the ABIC expression, ABIC can be rewritten to depend only on The function, Automatic optimization is performed as a latent variable.

[0049] 2. Data-driven adaptive weighting Construct a time difference matrix based on the travel time sampling time series. Simultaneously construct the transponder identification matrix This leads to the generation of the data covariance matrix: in This represents the Hadamard product (element-wise multiplication). Densely sampled regions have short time intervals between observations and large covariance values, resulting in automatically lower equivalent weights in the generalized least squares solution; sparse data receives relatively higher weights. This achieves an adaptive weighting effect of automatically reducing weights for dense data and automatically increasing weights for sparse data.

[0050] In the application of the measured data of this invention, high-precision GNSS observation equipment continuously collects ship positions at a sampling rate of over 10 Hz, with an interval of only 0.1 seconds between adjacent GNSS epochs, resulting in a high temporal correlation in transducer positioning errors. Using the covariance function of this invention, the covariance between adjacent observations is accurately modeled, and its equivalent weight decreases exponentially as the time interval decreases, thereby avoiding the overfitting problem of densely sampled segments.

[0051] 3. Joint optimization of covariance and smoothness (ABIC framework) Covariance parameter Hyperparameters related to sound velocity field smoothness By incorporating the ABIC criterion, optimal synchronization can be achieved. The ABIC expression is: in To design the matrix, It is a smoothing constraint matrix (usually a second-order difference matrix). The number of unknown parameters. This represents the weighted residual variance. The goal is to find the parameter combination that minimizes ABIC within the parameter domain using grid search or gradient descent. This framework ensures the objectivity and statistical optimality of hyperparameter selection.

[0052] The weighted residual variance (i.e., an unbiased estimate of the unit-weighted variance) is calculated using the following formula: ,in This is the residual vector. In hyperparameter optimization, this is a constant and does not affect the selection of the optimal parameter; it can be omitted in actual calculations.

[0053] Principles for setting the search range and step size of hyperparameters: ① (Time correlation length): Search range 0.1 ~ 60 seconds, typical step size 0.5 ~ 2 seconds. Time correlation can be ignored when it is less than the GNSS sampling interval (0.1s); it is considered to be fully correlated when it is greater than the total observation duration (60s).

[0054] ② (Cross-responder correlation coefficient): Search range 0.01~0.99, typical step size 0.05~0.1. 0 indicates complete independence, 1 indicates complete correlation, with an interval of 0.05 to balance accuracy and efficiency.

[0055] ③ (Smoothness parameter): Search range 10 -6 ~10 2 A logarithmic uniform distribution (3-5 points per order of magnitude) is used. Spanning multiple orders of magnitude, logarithmic uniform sampling can more effectively cover the effective area.

[0056] Convergence Criterion: Iteration terminates when one of the following conditions is met: ① Absolute convergence condition: ② Relative convergence condition: The maximum number of iterations is set to 500 (grid search) or 200 (gradient descent).

[0057] The sparsity of the covariance matrix C accelerates computation: This invention fully utilizes the following sparsity properties of the covariance matrix for acceleration: ① Time decay sparsity: when At that time, the exponent term The correlation is negligible. Therefore, only the following calculations are performed. For elements within the range, the storage complexity is from Down to ,in For effective bandwidth (typically) ② Diagonal segmented structure: when At this time, the covariance between different transponders is negligible, and the covariance matrix degenerates into a block diagonal matrix. ③ Cholesky decomposition acceleration: Sparse Cholesky decomposition is used for calculation. For an observation scale of n=5000, the computational complexity is from Down to .

[0058] 4. Covariance expansion based on multibeam observation data Multibeam bathymetry systems (such as the SeaBat T50-P) can acquire hundreds of beams in a single ping. Seabed measurement points between adjacent beams exhibit strong spatial correlation, and temporal correlation also exists between different pings. This invention treats the bathymetry data of each beam as an independent observation. For adjacent beams within the same ping, the correlation of their observation residuals is modeled using an extended covariance function. in The beam angle correlation coefficient is... For the first The incident angle of each beam. This extension enables multibeam bathymetry data to be jointly processed with GNSS-acoustic observations, improving the overall consistency of seabed topography and the accuracy of edge beams.

[0059] 5. Combining covariance expansion with ROV underwater robots ROVs (Remotely Operated Vehicles) are equipped with high-precision DVL (Doppler velocimeter, velocity accuracy better than 0.5%), fiber optic inertial navigation (attitude accuracy better than 0.01°), USBL (ultra-short baseline positioning system, positioning accuracy of 0.5%~1% of slant range), hyperspectral cameras, 2D imaging sonar, and lidar, among other sensors. When ROVs operate underwater, there is a strong temporal correlation between adjacent points on their continuous trajectory, and the observation errors between different sensors are also correlated.

[0060] This invention, tailored to the characteristics of ROV data, designs a hierarchical covariance structure: for continuous trajectory points output by DVL+INS integrated navigation, a time covariance function is used. Smoothing constraints are applied; for discrete positioning points obtained by USBL, cross-sensor correlation coefficients are used. A correlation is established with the DVL+ inertial navigation trajectory; for feature point clouds obtained from hyperspectral cameras and image sonar, consistency constraints are applied using a spatial covariance function. Finally, a unified covariance matrix for ROV multi-source observations is constructed. in This is the cross-covariance matrix between different sensors, whose off-diagonal elements are determined by the calibration parameters and spatiotemporal alignment relationships between the sensors. This covariance structure enables the joint calculation of ROV underwater position, attitude, imagery, and sonar data with surface GNSS-acoustic observations in strict consistency.

[0061] 6. Seabed Positioning Calculation and Output Substituting the optimal covariance matrix into the generalized least squares solver, we obtain the optimal estimate of the seabed reference point coordinates: Simultaneously output the posterior covariance matrix. Used for accuracy assessment.

[0062] The final output of this method includes: seabed reference point coordinates, posterior standard deviation, and estimated optimal covariance parameter. The results include a digital elevation model with multibeam sonic logging correction and a three-dimensional position sequence with smoothed ROV trajectories. These findings can be directly applied to monitoring seafloor plate movement, deep-sea resource exploration, and high-precision navigation for autonomous underwater vehicles.

[0063] In a specific example, the main method flow is as follows: ① Data acquisition phase: High-precision GNSS observation equipment (sampling rate 10-20 Hz), IMU (200 Hz), multibeam echo sounder (512 beams per ping), acoustic transducer (2-5 second cycle), and ROV underwater robot (equipped with DVL, inertial navigation, USBL, hyperspectral camera, image sonar, and lidar) acquire data simultaneously. All equipment is synchronized with microsecond-level clocks using GNSS PPS and UTC time.

[0064] ② Preprocessing stage: GNSS data is processed by PPK calculation to obtain centimeter-level transducer position timing; multibeam data is processed by acoustic ray bending correction and outlier removal; ROV data is processed by DVL + inertial navigation combined navigation to obtain underwater continuous trajectory; USBL data is processed by gross error detection.

[0065] ③ Covariance Modeling Stage: The spatiotemporal covariance matrix of GNSS-acoustic observations, the beam angle correlation covariance matrix of multi-beam observations, and the ROV multi-sensor covariance matrix are constructed respectively. The covariance matrices are correlated through cross-covariance terms to form a unified ultra-large-scale covariance matrix.

[0066] ④ ABIC Joint Optimization Stage: The covariance parameter... With the smoothness parameter of the sound velocity field By incorporating the ABIC framework and simultaneously optimizing through 3D mesh search (or gradient descent), the computational complexity of ABIC is O(n log n). For observations of tens of thousands, sparse matrix techniques and Cholesky decomposition are used to accelerate the process.

[0067] ⑤ Joint Solution Stage: Based on the optimal covariance matrix and optimal smoothness parameters, a generalized least squares equation is constructed to solve for the seabed reference point coordinates, multibeam topography corrections, and ROV trajectory corrections. The solution process employs a parallel computing strategy, processing massive amounts of beam data in blocks.

[0068] ⑥ Output and Application Stage: Outputs centimeter-level seabed reference point coordinates, sub-decimeter-level multibeam digital elevation model of terrain, decimeter-level ROV smoothed trajectory, and various optimal covariance parameters to guide subsequent marine mapping tasks.

[0069] In summary, this invention provides a GNSS-acoustic seabed positioning data weighting method based on spatiotemporal covariance adaptive modeling. Its important innovative contributions and technical advantages include: Two-parameter spatiotemporal covariance structure: Introducing time-dependent length Correlation coefficient across transponders It simultaneously characterizes the temporal correlation of observation errors and the independence of errors from different transponders. This covariance function has a clear physical interpretation and flexible adaptive capability, and its computational efficiency is suitable for on-orbit processing on satellites / ships.

[0070] Data-driven adaptive weighting method: Based on the temporal sampling density and transponder identifiers of the observed data, the off-diagonal elements of the covariance matrix are dynamically adjusted to achieve automatic weight reduction of densely sampled data, weight enhancement of sparse data, and independent modeling of different transponder errors, which significantly suppresses overfitting and improves spatial resolution.

[0071] A joint optimization framework based on covariance and smoothness using ABIC: This framework integrates the covariance parameter... With model smoothness hyperparameter By incorporating the unified Akaike Bayesian information criterion, synchronous and objective optimization of all hyperparameters is achieved, avoiding the subjectivity and uncertainty of manual parameter tuning, while improving the model's generalization ability to different sea areas.

[0072] Compared with the prior art, the significant technical advantages of the present invention are reflected in the following aspects: Significantly improved accuracy: By accurately characterizing the spatiotemporal correlation of observation errors, the weight of densely sampled data is automatically reduced, effectively avoiding overfitting. Simulation and experimental data verification show that the root mean square error (RMSE) of the solution is reduced by 30%~40% compared to the equal weight method and by 15%~20% compared to the empirical weighting method.

[0073] Spatial consistency is significantly improved: cross-responder correlation coefficient By modeling the errors of different transponders independently, error coupling between different transponders is effectively avoided. The baseline length repeatability standard deviation of multi-transponder array solutions is improved by about 40%, and the spatial resolution is greatly improved.

[0074] Fully adaptive, requiring no human experience: The covariance parameter is adaptively determined from the data through the ABIC framework without human intervention. It can adapt to the changes in error characteristics under different sea areas, seasons, and observation platforms, and has strong universality and engineering applicability.

[0075] This invention also provides a storage medium for storing a computer program, which, when executed, performs at least the methods described above.

[0076] This invention also provides a control device, including a processor and a storage medium for storing a computer program; wherein the processor executes the computer program by performing at least the method described above.

[0077] This invention also provides a processor that executes a computer program, at least performing the methods described above.

[0078] The storage medium can be implemented by any type of non-volatile storage device, or a combination thereof. The non-volatile memory can be read-only memory (ROM), programmable read-only memory (PROM), erasable programmable read-only memory (EPROM), electrically erasable programmable read-only memory (EEPROM), magnetic random access memory (FRAM), flash memory, magnetic surface memory, optical disc or CD-ROM; magnetic surface memory can be disk storage or magnetic tape storage. The storage media described in the embodiments of this invention are intended to include, but are not limited to, these and any other suitable types of memory.

[0079] In the several embodiments provided by this invention, it should be understood that the disclosed systems and methods can be implemented in other ways. The device embodiments described above are merely illustrative. For example, the division of units is only a logical functional division, and in actual implementation, there may be other division methods, such as: multiple units or components can be combined, or integrated into another system, or some features can be ignored or not executed. In addition, the coupling or direct coupling or communication connection between the various components shown or discussed can be through some interfaces, and the indirect coupling or communication connection between devices or units can be electrical, mechanical, or other forms.

[0080] The units described above as separate components may or may not be physically separate. The components shown as units may or may not be physical units, that is, they may be located in one place or distributed across multiple network units. Some or all of the units may be selected to achieve the purpose of this embodiment according to actual needs.

[0081] In addition, in the various embodiments of the present invention, each functional unit can be integrated into one processing unit, or each unit can be a separate unit, or two or more units can be integrated into one unit; the integrated unit can be implemented in hardware or in the form of hardware plus software functional units.

[0082] Those skilled in the art will understand that all or part of the steps S of the above method embodiments can be implemented by hardware related to program instructions. The aforementioned program can be stored in a computer-readable storage medium. When the program is executed, it performs steps S of the above method embodiments. The aforementioned storage medium includes various media capable of storing program code, such as mobile storage devices, read-only memory (ROM), random access memory (RAM), magnetic disks, or optical disks.

[0083] Alternatively, if the integrated units of this invention are implemented as software functional modules and sold or used as independent products, they can also be stored in a computer-readable storage medium. Based on this understanding, the technical solutions of the embodiments of this invention, or the parts that contribute to the prior art, can be embodied in the form of a software product. This computer software product is stored in a storage medium and includes several instructions to cause a computer device (which may be a personal computer, server, or network device, etc.) to execute all or part of the methods described in the various embodiments of this invention. The aforementioned storage medium includes various media capable of storing program code, such as mobile storage devices, ROM, RAM, magnetic disks, or optical disks.

[0084] The methods disclosed in the several method embodiments provided by this invention can be arbitrarily combined without conflict to obtain new method embodiments.

[0085] The features disclosed in the several product embodiments provided by this invention can be arbitrarily combined without conflict to obtain new product embodiments.

[0086] The features disclosed in the several method or device embodiments provided by the present invention can be arbitrarily combined without conflict to obtain new method or device embodiments.

[0087] The above description, in conjunction with specific preferred embodiments, provides a further detailed explanation of the present invention. It should not be construed that the specific implementation of the present invention is limited to these descriptions. For those skilled in the art, various equivalent substitutions or obvious modifications can be made without departing from the concept of the present invention, and all such modifications, achieving the same performance or application, should be considered within the scope of protection of the present invention.

Claims

1. A weighting method for GNSS-acoustic seabed positioning data based on spatiotemporal covariance adaptive modeling, characterized in that, Includes the following steps: S1. Acquire multi-source observation data and perform preprocessing. The multi-source observation data includes high-precision GNSS data, acoustic transducer data, multibeam bathymetry data, and ROV underwater robot data. S2. Construct a two-parameter spatiotemporal covariance matrix to characterize the temporal correlation of observation errors and the independence of errors of different transponders. Based on the temporal sampling density of the observation data and the transponder identifier, adaptive weighting is achieved by automatically reducing the weight of densely sampled data and automatically increasing the weight of sparse data through the off-diagonal elements of the covariance matrix. S3. Incorporate the covariance parameter and the sound velocity field smoothness hyperparameter into the Akaike Bayesian Information Criterion (ABIC) framework, and achieve synchronous objective optimization of all hyperparameters through a multi-dimensional optimization strategy. S4. Substitute the optimized covariance matrix and the optimal smoothness parameter into the generalized least squares solver to solve for the coordinates of the seabed reference point and output the positioning results.

2. The GNSS-acoustic seabed positioning data weighting method based on spatiotemporal covariance adaptive modeling according to claim 1, characterized in that, The data acquisition and preprocessing in step S1 includes: all devices synchronizing with microsecond-level clocks using GNSS second pulses and UTC time; performing dynamic post-processing differential calculations on GNSS data to obtain centimeter-level transducer position timing, and adding the rod arm vector after rotation matrix transformation to the GNSS antenna position coordinates to obtain the transducer position coordinates; performing ray bending correction on acoustic data, correcting the observed round-trip propagation time to the theoretical propagation time along the actual ray path; removing outliers from multibeam data; and performing time alignment and underwater trajectory gross error detection on ROV data to generate standardized observation values ​​and initial residual values.

3. The GNSS-acoustic seabed positioning data weighting method based on spatiotemporal covariance adaptive modeling according to claim 1, characterized in that, The process of constructing the two-parameter spatiotemporal covariance matrix in step S2 includes: A time difference matrix is ​​constructed based on the sampling time of each observation, where each element of the matrix is ​​the absolute value of the difference between two corresponding observation times. At the same time, a transponder identification matrix is ​​constructed based on the transponder number corresponding to each observation. When two observations belong to the same transponder, the identification element takes a value of 0, and when they belong to different transponders, the value takes a value of 1. The time decay factor matrix is ​​obtained by dividing each element of the time difference matrix by the time correlation length parameter and taking the negative, and then performing an exponential decay operation with the natural constant as the base; the cross-responder correlation factor matrix is ​​obtained by performing an exponential operation with the cross-responder correlation coefficient parameter as the base and each element of the responder identification matrix as the exponent. The time decay factor matrix and the cross-responder correlation factor matrix are multiplied element-wise, and then multiplied by the variance scaling parameter to generate a two-parameter spatiotemporal covariance matrix, thereby achieving joint modeling of time decay and the independence between transponders.

4. The GNSS-acoustic seabed positioning data weighting method based on spatiotemporal covariance adaptive modeling according to claim 3, characterized in that, The adaptive weighting is achieved through the following mechanism: the time interval between adjacent observations in the dense sampling region is small, the value of the corresponding element in the time decay factor matrix is ​​close to 1, and the value of the corresponding element in the covariance matrix is ​​large. In the generalized least squares solution, the weight of the observation is determined by the inverse of the covariance matrix. When the covariance value is large, the corresponding equivalent weight is automatically reduced. In sparsely sampled data, the time interval between observations is large, the value of the corresponding element in the time decay factor matrix approaches 0, the value of the corresponding element in the covariance matrix is ​​small, and the corresponding equivalent weight automatically increases; the observations of different transponders are independently modeled after exponentiation of the cross-transponder correlation coefficient parameter. When the value of this parameter is small, the covariance between the observations of different transponders approaches 0, thus achieving decoupling of the error structure.

5. The GNSS-acoustic seabed positioning data weighting method based on spatiotemporal covariance adaptive modeling according to claim 1, characterized in that, The joint optimization process of the Akaike Bayesian information criterion framework described in step S3 includes: Using the covariance parameter and the sound velocity field smoothness hyperparameter as the variables to be optimized, the objective function of the Akaike Bayesian information criterion is constructed. The objective function is composed of the following terms: the number of observations minus the number of unknown parameters multiplied by the logarithm of the weighted residual variance, minus the design matrix transpose multiplied by the inverse of the covariance matrix and then multiplied by the design matrix, plus the logarithm of the determinant of the matrix obtained by multiplying the smoothness hyperparameter by the smoothing constraint matrix, plus the logarithm of the determinant of the covariance matrix, and a constant term; Within the parameter domain, multidimensional optimization is performed using grid search or gradient descent. For each set of candidate parameter combinations, the time difference matrix is ​​calculated sequentially, the covariance matrix is ​​constructed, the weighted residual variance is calculated, and the above items are obtained to obtain the numerical value of the Akaike Bayesian information criterion. By traversing all candidate parameter combinations or iterating along the gradient descent direction, the parameter combination that minimizes the Akaike Bayesian information criterion is selected as the optimal hyperparameter, thus achieving the synchronous and objective determination of the covariance parameter and the smoothness hyperparameter. The iteration terminates when the change in the objective function is less than the preset convergence threshold.

6. The GNSS-acoustic seabed positioning data weighting method based on spatiotemporal covariance adaptive modeling according to claim 1, characterized in that, The process of achieving synchronous objective optimization of hyperparameters through a multi-dimensional optimization strategy in step S3 also includes: The variance scaling parameter does not participate in the ABIC iterative optimization. Instead, it is directly calculated analytically given the covariance parameter and the smoothness hyperparameter. A structure matrix is ​​defined such that the covariance matrix equals the variance scaling parameter multiplied by the structure matrix. The elements of the structure matrix are determined by the product of the power operation of the cross-responder correlation coefficient and the time decay exponent. In each ABIC iteration, the structure matrix is ​​constructed given the current hyperparameter combination, and the generalized least squares method is used to obtain the estimated values ​​of the unknown parameters. Then, the residual vector is calculated based on the difference between the observed value vector and the design matrix multiplied by the estimated value of the unknown parameters. The transpose of the residual vector is multiplied by the inverse of the structure matrix, then multiplied by the residual vector again, and divided by the total number of observations to obtain the estimated value of the variance scaling parameter. Finally, the variance scaling parameter is substituted into the covariance matrix for the calculation of the next ABIC expression. After substituting the analytical expression of the variance scaling parameter into the ABIC expression, ABIC is rewritten as a function that depends only on the covariance parameter and the smoothness hyperparameter, with the variance scaling parameter being automatically optimized as a latent variable. The covariance matrix sparsity acceleration mechanism includes: utilizing time decay sparsity, where the exponential decay term is less than a threshold when the time interval between two observations exceeds a preset multiple of the time correlation length, and the correlation can be ignored. Only elements within this time interval are calculated, reducing the storage complexity from the square of the total number of observations to the order of the total number of observations multiplied by the effective bandwidth; and utilizing a block diagonal structure, where the covariance between different transponders can be ignored when the cross-responder correlation coefficient is lower than a preset threshold, and the covariance matrix degenerates into a block diagonal matrix.

7. The GNSS-acoustic seabed positioning data weighting method based on spatiotemporal covariance adaptive modeling according to claim 1, characterized in that, The generalized least squares solution in step S4 includes: multiplying the inverse of the optimal covariance matrix by the transpose of the design matrix and then by the observation vector to obtain a constant term; simultaneously, multiplying the transpose of the design matrix by the inverse of the covariance matrix and then by the design matrix, adding the smoothness hyperparameter and multiplying by the smoothing constraint matrix to construct the normal matrix; by solving the inverse of the normal matrix and multiplying by the constant term, the optimal estimate of the seabed reference point coordinate increment is obtained, and the increment is superimposed on the initial coordinates to update the parameters; repeating the above process until the parameter change is less than the convergence threshold; calculating the posterior covariance matrix based on the inverse of the normal matrix, and the square root of its diagonal elements is the posterior standard deviation of the coordinates, used for accuracy assessment.

8. The GNSS-acoustic seabed positioning data weighting method based on spatiotemporal covariance adaptive modeling according to claim 1, characterized in that, The positioning results output in step S4 include: the coordinates of the seabed reference point and its posterior standard deviation, the posterior standard deviation being obtained by taking the square root of the diagonal elements of the posterior covariance matrix; the estimated optimal covariance parameters, including the time correlation length parameter, the transponder correlation coefficient parameter, and the smoothness hyperparameter; the digital elevation model after multibeam terrain correction; and the three-dimensional position sequence after ROV trajectory smoothing.

9. The GNSS-acoustic seabed positioning data weighting method based on spatiotemporal covariance adaptive modeling according to claim 1, characterized in that, The method also includes covariance expansion of multibeam bathymetry data. The process includes: treating the bathymetry data of each beam as an independent observation, constructing an extended covariance function for the correlation of observation residuals of adjacent beams within the same transmission cycle; this function first calculates the absolute value of the difference in incident angles of each beam to obtain an angle difference matrix, and then performs a power operation with the beam angle correlation coefficient as the base and each element of the angle difference matrix as the exponent to obtain the angle correlation factor; simultaneously, for observations in different transmission cycles, the time decay factor is calculated by dividing the time difference by the time correlation length, taking the negative, and then exponentially decaying; the angle correlation factor and the time decay factor are multiplied together, and then multiplied by the multibeam variance scaling parameter to generate a multibeam covariance matrix, thereby realizing the joint modeling of spatial correlation caused by beam angle differences and time correlation between different transmission cycles.

10. The GNSS-acoustic seabed positioning data weighting method based on spatiotemporal covariance adaptive modeling according to claim 1, characterized in that, The method also includes covariance expansion of the multi-source sensors of the ROV underwater robot, the process of which includes: For continuous trajectory points output by DVL and inertial navigation integration, a time covariance function is constructed, and a time smoothing constraint matrix is ​​generated by dividing the time difference by the ROV time correlation length, taking the negative, and then exponentially decaying. For the discrete positioning points obtained by USBL, a cross-sensor correlation coefficient is introduced. This coefficient is used to scale the covariance matrix of the DVL-INS trajectory to establish the correlation covariance submatrix between USBL observations and continuous trajectories. For feature point clouds obtained by hyperspectral cameras and image sonar, a spatial covariance function is constructed to constrain the consistency of the feature point clouds. The covariance submatrices corresponding to each of the above sensors are used as the main diagonal blocks of the block diagonal matrix. Then, the cross covariance matrix is ​​calculated based on the calibration parameters and spatiotemporal alignment relationship between different sensors and used as the off-diagonal block. The block diagonal matrix and the cross covariance matrix are added together to obtain the unified covariance matrix of ROV multi-source observation.