A deep-sea multi-source fusion global navigation method based on delay backtracking compensation

Abstract

This invention proposes a deep-sea multi-source fusion global navigation method based on time delay backtracking compensation, comprising the following steps: decomposing the actual state of the submersible into a nominal state and an error state; using the output of the inertial measurement unit to perform high-frequency integration extrapolation on the nominal state to obtain the predicted value of the nominal state; and recursively deriving the prediction error covariance matrix; acquiring the velocity measurement of the Doppler log, calculating the velocity measurement residual according to its working mode and introducing ocean current compensation; performing Kalman update to correct the nominal state to obtain the locally compensated nominal state; establishing a historical state buffer; when receiving an acoustic position measurement with a timestamp, calculating its time delay and retrieving the corresponding historical state; performing Kalman update to obtain the corrected historical state; and re-integrating and recursively deriving to the current time to obtain the current nominal state, which serves as the final global pose estimate of the submersible. This method can enable the submersible to maintain high-precision and high-stability navigation in complex deep-sea environments.

Description

Technical Field

[0001] This invention relates to the field of deep-sea submersible navigation and positioning technology, and in particular to a deep-sea multi-source fusion global navigation method based on time delay backtracking compensation. Background Technology

[0002] Deep-sea hydrothermal sampling operations typically take place at complex seabeds thousands of meters deep. In the deep-sea environment, electromagnetic waves attenuate rapidly, resulting in the complete loss of GNSS signals. Maintaining high-precision global positioning is a key challenge during the process of a bottom-diving submersible being released from its mother ship and descending to within several kilometers of the target seabed.

[0003] In the deep-sea environment, there are no GNSS signals, and submersible navigation relies on the fusion of IMU, DVL, and acoustic positioning systems. However, pure IMU dead reckoning suffers from rapidly increasing position drift over time, DVL is prone to abnormal velocity jumps when tracking hydrothermal plumes or lost bottoms, and USBL / LBL acoustic positioning suffers from second-level transmission delays and low update rates. Furthermore, the traditional EKF loose coupling method does not accurately compensate for acoustic delays, treating delayed observations as instantaneous quantities, which easily leads to trajectory trailing and positioning jumps. It also lacks robust fault tolerance for strong interferences such as underwater acoustic multipath and hydrothermal plumes, and cannot provide continuous, stable, and high-precision global positioning in the complex deep-sea environment. Summary of the Invention

[0004] In view of this, the present invention proposes a deep-sea multi-source fusion global navigation method based on time delay backtracking compensation, which can provide submersibles with continuous, smooth, and trajectory-free global pose output, maintain high-precision and high-stability navigation in complex deep-sea environments, and meet the positioning requirements of deep-sea operations.

[0005] The technical solution of this invention is implemented as follows: In a first aspect, this invention provides a deep-sea multi-source fusion global navigation method based on time delay backtracking compensation, comprising the following steps: S1, decompose the actual state of the submersible into nominal state and error state, establish the kinematic equation of nominal state, and use the output of the inertial measurement unit to perform high-frequency integration extrapolation on nominal state to obtain the predicted value of nominal state. S2, based on the kinematic equation of the nominal state, performs first-order linearization to obtain the error state transition matrix and noise driving matrix. After discretization, the error covariance matrix of the current time is recursively predicted using the error covariance of the previous time step. S3. During the acoustic position measurement signal transmission gap, the velocity measurement of the Doppler log is acquired. The velocity measurement residual is calculated according to its working mode and ocean current compensation is introduced. Kalman update is performed to correct the nominal state and obtain the nominal state after local compensation. S4. Establish a historical state buffer to store the nominal state, error covariance matrix, and raw data of the inertial measurement unit at each time. When an acoustic position measurement with a timestamp is received, calculate its time delay and retrieve the corresponding historical state. Perform Kalman update to obtain the corrected historical state. Then, use the raw data of the inertial measurement unit to reintegrate from the corrected historical state to the current time to obtain the current nominal state, which serves as the final global pose estimate of the submersible.

[0006] Based on the above technical solutions, preferably, step S1 includes the following sub-steps: S11 defines the true state vector of the submersible, which includes three-dimensional position, three-dimensional velocity, attitude quaternion, three-dimensional zero bias of accelerometer and three-dimensional zero bias of gyroscope. S12 decomposes the actual state of the submersible into a nominal state and an error state, and the nominal state and the error state are combined to form the actual state through manifold addition. S13. Based on the real state vector, establish the kinematic integral equation corresponding to the nominal state. With the sampling period of the inertial measurement unit as the step size and the angular velocity and specific force output by the inertial measurement unit as the driving variables, perform high-frequency integral extrapolation on the nominal state according to the kinematic integral equation to obtain the predicted value of the nominal state at each moment.

[0007] Based on the above technical solutions, preferably, step S2 includes the following sub-steps: S21, Define the error state vector, which includes three-dimensional position error, three-dimensional velocity error, three-dimensional attitude small angle error, three-dimensional zero bias error of accelerometer and three-dimensional zero bias error of gyroscope; S22, obtain the Jacobian matrix for the kinematic equation of the nominal state to obtain the error state transition matrix in continuous time, and construct the noise driving matrix to map the accelerometer white noise, gyroscope white noise, accelerometer zero-bias random walk noise and gyroscope zero-bias random walk noise to the corresponding error state space. S23, using the sampling period of the inertial measurement unit as the discrete time step, the error state transition matrix is ​​converted into a discrete time state transition matrix using a second-order approximation. Based on the pre-calibrated continuous-time process noise power spectral density matrix, the discrete process noise covariance matrix is ​​obtained by approximation using the noise driving matrix and the sampling interval. S24. Multiply the discrete-time state transition matrix with the error covariance matrix of the previous time step and the transpose of the discrete-time state transition matrix, and add the discrete process noise covariance matrix to obtain the prediction error covariance matrix at the current time step.

[0008] Based on the above technical solutions, preferably, step S3 includes the following sub-steps: S31, acquire the three-dimensional velocity measurement value in the carrier coordinate system output by the Doppler log and the returned lock status word, identify the current working mode according to the lock status word, the working mode includes bottom tracking mode and water tracking mode; S32, if it is the bottom tracking mode, calculate the rotation matrix from the carrier coordinate system to the navigation coordinate system based on the attitude quaternion in the current nominal state, transform the carrier speed measured by the Doppler log to the navigation coordinate system through the rotation matrix to obtain the navigation speed, and subtract the navigation speed in the nominal state from the navigation speed to obtain the speed measurement residual. S33, If it is water tracking mode, the velocity measurement residual calculated by bottom tracking mode is added to the current estimated ocean current velocity estimate to obtain the compensated velocity measurement residual, and the value of the observation noise covariance matrix is ​​replaced with a preset limit value, which is greater than the value of the observation noise covariance matrix in bottom tracking mode. S34. Determine the corresponding velocity measurement matrix based on the velocity measurement residual under the corresponding working mode. Perform Kalman update based on the currently predicted error covariance matrix, velocity measurement matrix and corresponding observation noise covariance matrix to obtain the optimal estimate of the error state, which is used to correct the nominal state and obtain the nominal state after local compensation.

[0009] Based on the above technical solutions, preferably, step S4 includes the following sub-steps: S41, Establish a historical state buffer, which stores the nominal state, error covariance matrix and raw data of the inertial measurement unit at each time. S42, when a global position measurement is received from the acoustic positioning system, the total delay is calculated based on the timestamp of the global position measurement and the current system time; S43, based on the timestamp of the global position measurement, search for the matching index in the historical state buffer, and extract the nominal state and error covariance matrix of the corresponding historical moment; S44. Based on the difference between the global position measurement and the position in the historical nominal state, the position measurement residual is calculated. The position measurement matrix is ​​determined according to the definition of the position measurement residual. Using the error covariance matrix, the position measurement matrix and the noise covariance matrix of the acoustic measurement at the historical moment, Kalman update is performed to obtain the optimal estimate of the error state, which is used to correct the historical nominal state and obtain the corrected historical state. S45, from the next moment of the historical moment to the current moment, traverse the original data of each frame of the inertial measurement unit in the historical state buffer in turn. Starting from the corrected historical state, recalculate the nominal state at each moment by integrating frame by frame according to the kinematic integral equation of the nominal state, and update the state value of the corresponding moment in the historical state buffer. S46, the nominal state obtained at the current moment is used as the final global pose estimate of the submersible.

[0010] Based on the above technical solutions, preferably, the Kalman update includes the following sub-steps: Calculate the innovation covariance matrix based on the prediction error covariance matrix, measurement matrix, and observation noise covariance matrix at the corresponding time. The Kalman gain is obtained by multiplying the prediction error covariance matrix by the transpose of the measurement matrix and then by the inverse of the innovation covariance matrix. The Kalman gain is multiplied by the measurement residual at the corresponding time to obtain the optimal estimate of the error state. The optimal estimate of the error state is then injected into the nominal state at the corresponding time through manifold addition to correct the position, velocity, attitude and inertial zero bias of the nominal state. The error state vector is reset to zero. The updated posterior error covariance matrix is ​​obtained by subtracting the Kalman gain multiplied by the measurement matrix from the identity matrix and then multiplying it by the prediction error covariance matrix.

[0011] Based on the above technical solutions, preferably, step S5 is also included for robust fault tolerance before measurement update, including the following sub-steps: S51, Before each measurement update, calculate the Mahalanobis distance based on the new information covariance matrix and the corresponding measurement residuals; S52, compare the Mahalanobis distance with the preset first threshold. If the Mahalanobis distance is less than or equal to the first threshold, determine that the current measurement is a normal measurement, accept the current measurement and perform Kalman update, and update the accelerometer zero bias and gyroscope zero bias of the inertial measurement unit online. S53, if the Mahalanobis distance is greater than the first threshold and less than or equal to the preset second threshold, it is determined that there is strong interference in the current measurement, the observation noise covariance matrix is ​​enlarged proportionally, and Kalman update is performed. S54. If the Mahalanobis distance is greater than the second threshold, the current measurement sensor is determined to be faulty. Kalman update is not performed, and the system enters the inertial traverse mode. Using the pre-calibrated accelerometer zero bias and gyroscope zero bias, pure inertial integration is performed based on the angular velocity and specific force output of the inertial measurement unit to maintain the recursive estimation of the submersible's position, velocity, and attitude until the sensor signal returns to normal.

[0012] Secondly, the present invention provides a deep-sea multi-source fusion global navigation system based on time delay backtracking compensation, implemented using a deep-sea multi-source fusion global navigation method based on time delay backtracking compensation, comprising: The nominal state prediction module is used to decompose the actual state of the submersible into a nominal state and an error state, establish the kinematic equation of the nominal state, and use the output of the inertial measurement unit to perform high-frequency integral extrapolation on the nominal state to obtain the predicted value of the nominal state. The error covariance prediction module is used to perform first-order linearization of the kinematic equations based on the nominal state to obtain the error state transition matrix and the noise driving matrix. After discretization, the error covariance matrix of the current time is recursively predicted using the error covariance of the previous time step. The local velocity compensation module is used to acquire the velocity measurement of the Doppler log during the gap in the acoustic position measurement signal transmission, calculate the velocity measurement residual according to its working mode and introduce ocean current compensation, perform Kalman update to correct the nominal state, and obtain the nominal state after local compensation. The global backtracking compensation module is used to establish a historical state buffer, storing the nominal state, error covariance matrix, and raw data of the inertial measurement unit at each moment. When an acoustic position measurement with a timestamp is received, its time delay is calculated and the corresponding historical state is retrieved. Kalman update is performed to obtain the corrected historical state, and the raw data of the inertial measurement unit is used to re-integrate from the corrected historical state to the current moment to obtain the current nominal state, which serves as the final global pose estimate of the submersible.

[0013] Thirdly, the present invention provides a terminal device, the terminal device comprising: a memory, a processor, and a program for a deep-sea multi-source fusion global navigation method based on time delay backtracking compensation stored in the memory and executable on the processor, the program for the deep-sea multi-source fusion global navigation method based on time delay backtracking compensation being configured to implement the steps of the deep-sea multi-source fusion global navigation method based on time delay backtracking compensation.

[0014] Fourthly, the present invention provides a computer-readable storage medium storing a deep-sea multi-source fusion global navigation method and system program based on time delay backtracking compensation. When the program for the deep-sea multi-source fusion global navigation method based on time delay backtracking compensation is executed, the deep-sea multi-source fusion global navigation method based on time delay backtracking compensation is implemented.

[0015] The deep-sea multi-source fusion global navigation method based on time delay backtracking compensation of the present invention has the following advantages over the prior art: (1) By adopting time delay backtracking compensation based on error state Kalman filtering and multi-source information level fusion, the IMU high-frequency motion information, DVL velocity constraints and acoustic global positioning are effectively combined to accurately compensate for the second-level time delay caused by underwater acoustic signal transmission. This eliminates the trajectory trailing and positioning jump effects caused by improper delay processing in traditional filtering. In the complex environment of the deep sea, continuous, smooth and stable global pose output is achieved, which improves the accuracy and stability of global navigation. (2) Through adaptive switching and speed correction of DVL bottom tracking and water tracking, it can automatically adapt to the entire process of the submersible from sea surface diving, mid-level navigation to near-bottom operation, complete the smooth switching of ground speed and water speed and ocean current compensation, effectively suppress the short-term drift of pure inertial navigation, and improve the environmental adaptability and continuous working capability of the navigation system. (3) By constructing a historical state buffer and executing a backtracking correction and repropagation mechanism, the delayed acoustic position observations are accurately projected to the corresponding historical time for correction, and then re-propelled to the current time along the buffered IMU data, ensuring the physical continuity of state updates, avoiding trajectory mutations caused by delayed correction, and improving the smoothness and consistency of global positioning. (4) Through the detection of anomalies in Mahalanobis distance and the fault tolerance mechanism of inertial crossing, strong interferences such as hydrothermal plumes and underwater acoustic multipaths can be quickly identified and graded. When the sensor fails for a short time, the low drift pure inertial navigation is maintained by relying on the converged high-precision IMU zero bias, which improves the anti-interference ability and operational safety of the system in extreme deep-sea environments. Attached Figure Description

[0016] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0017] Figure 1 This is a flowchart of the deep-sea multi-source fusion global navigation method based on time delay backtracking compensation according to the present invention; Figure 2 This is a schematic diagram of the repropagation timing of acoustic time delay backtracking compensation in the deep-sea multi-source fusion global navigation method based on time delay backtracking compensation of the present invention. Detailed Implementation

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

[0019] like Figure 1-2 As shown, the present invention provides a deep-sea multi-source fusion global navigation method based on time delay backtracking compensation, comprising the following steps: S1 decomposes the actual state of the submersible into a nominal state and an error state, establishes the kinematic equation of the nominal state, and uses the output of the inertial measurement unit to perform high-frequency integral extrapolation on the nominal state to obtain the predicted value of the nominal state.

[0020] Step S1 includes the following sub-steps: S11 defines the true state vector of the submersible, which includes three-dimensional position, three-dimensional velocity, attitude quaternion, three-dimensional zero bias of accelerometer and three-dimensional zero bias of gyroscope. The true state is a 15-dimensional vector, including three-dimensional position. p ∈ R 3 ,speed v ∈ R 3 Posture Quaternions q ∈ S 3 accelerometer zero bias b a ∈ R 3 Zero bias of the gyroscope b g ∈ R 3 .

[0021] S12, decompose the actual state of the submersible into a nominal state and an error state, and the nominal state and error state are combined to form the actual state through manifold addition; the expression is: ; In the formula, X t This represents the actual state of the submersible. X In nominal condition, δX This is an error state.

[0022] S13. Based on the real state vector, establish the kinematic integral equation corresponding to the nominal state. With the sampling period of the inertial measurement unit as the step size and the angular velocity and specific force output by the inertial measurement unit as the driving variables, perform high-frequency integral extrapolation on the nominal state according to the kinematic integral equation to obtain the predicted value of the nominal state at each moment. Among these, the kinematic equations of the nominal state are established, and the angular velocity measured by the IMU is used. oh m And comparison a m As the driving variable, high-frequency pure integral extrapolation is performed, and the derivatives of the nominal state position, velocity, attitude quaternions, accelerometer zero bias, and gyroscope zero bias are calculated respectively. The expressions are as follows:

[0023] In the formula,R ( q ) is a rotation matrix obtained by quaternion transformation. g This is the local gravity vector. This is quaternion multiplication. To expand the angular velocity difference into a quaternion with a pure imaginary part.

[0024] By integrating the above steps, we obtain the predicted value X of the nominal state at each time step. k|k-1 .

[0025] In this embodiment, high-frequency IMU data is used to provide a continuous and smooth basic motion trajectory, laying the time reference for subsequent fusion. The decomposition of nominal state and error state separates nonlinear integral and linear filtering, avoiding the error caused by directly linearizing the nonlinear system using EKF.

[0026] S2, based on the kinematic equations of the nominal state, performs first-order linearization to obtain the error state transition matrix and noise driving matrix. After discretization, the error covariance matrix of the current time step is recursively predicted using the error covariance of the previous time step.

[0027] Step S2 includes the following sub-steps: S21, Define the error state vector, which includes three-dimensional position error, three-dimensional velocity error, three-dimensional attitude small angle error, three-dimensional zero bias error of accelerometer and three-dimensional zero bias error of gyroscope; The error state vector is a 15-dimensional vector. δX ∈R 15 :

[0028] In the formula, δp For positional error, dv For speed error, dth For small angle attitude error, To achieve three-dimensional zero bias error of the accelerometer, This represents the three-dimensional zero-bias error of the gyroscope.

[0029] S22, obtain the Jacobian matrix for the kinematic equation of the nominal state to obtain the error state transition matrix in continuous time, and construct the noise driving matrix to map the accelerometer white noise, gyroscope white noise, accelerometer zero-bias random walk noise and gyroscope zero-bias random walk noise to the corresponding error state space. By linearizing the nominal state kinematic equations to the first order, we obtain the continuous-time error state transition matrix, which is expressed as follows: ; In the formula, I n×n It is an n-dimensional identity matrix, 0 n×nFor an n-dimensional zero matrix, [·] × Let the antisymmetric matrix of a three-dimensional vector be denoted as , for any a =[ a 1, a 2, a 3] T ,but: ; Noise driving matrix G The expression is:

[0030] In the formula, n a For accelerometer white noise, n g For gyroscope white noise, n ba For accelerometer zero-bias random walk noise, n bg This is the gyroscope's zero-bias random walk noise.

[0031] S23, using the sampling period of the inertial measurement unit as the discrete time step, the error state transition matrix is ​​converted into a discrete time state transition matrix using a second-order approximation. Based on the pre-calibrated continuous-time process noise power spectral density matrix, the discrete process noise covariance matrix is ​​obtained by approximation using the noise driving matrix and the sampling interval. The expression for the error covariance matrix is: ; In the formula, Φ k This is the discrete-time state transition matrix, which describes how the system error state evolves within a time interval. F The continuous-time error state transition matrix reflects the derivative coupling relationship between position, velocity, and attitude errors, Δ t The IMU sampling interval I It is a 15-dimensional identity matrix, representing the self-maintained part of the error state.

[0032] The expression for the noise covariance of a discrete process is: ; In the formula, Q d The discrete-time process noise covariance matrix is ​​quantized at the IMU sampling interval Δ. t During this period, the cumulative uncertainty caused by random noise injected into the IMU system state; G The noise-driven matrix determines how the IMU's raw physical noise is distributed across the 15-dimensional error states. Q cThis is the noise power spectral density matrix for a continuous-time process; the parameters are obtained from the IMU manual calibration. ; In the formula, The power spectral density of the accelerometer white noise reflects the influence of high-frequency vibration and electronic noise on velocity estimation. The power spectral density of the gyroscope's white noise reflects the instantaneous drift during attitude estimation. The power spectral density of the zero-bias random walk noise of the accelerometer is typically simulated as a random walk. The power spectral density is the zero-bias random walk noise of the gyroscope.

[0033] S24. Multiply the discrete-time state transition matrix with the error covariance matrix of the previous time step and the transpose of the discrete-time state transition matrix, and add the discrete process noise covariance matrix to obtain the prediction error covariance matrix at the current time step. The expression for the prediction error covariance matrix at the current time is: ; In the formula, P k-1|k-1 Let be the error covariance matrix from the previous time step, representing the system's assessment of the uncertainty in estimating the current state. P k|k-1 The error covariance matrix for the predicted time reflects the increase in uncertainty due to noise accumulation after inertial integration.

[0034] This embodiment quantifies the propagation law of uncertainty in the IMU integration process, providing a statistical basis for the subsequent weight allocation of Kalman filtering.

[0035] S3: During the gap in the acoustic position measurement signal transmission, acquire the velocity measurement of the Doppler log, calculate the velocity measurement residual according to its working mode and introduce ocean current compensation, perform Kalman update to correct the nominal state, and obtain the nominal state after local compensation.

[0036] During the acoustic position measurement signal transmission gap, i.e. between two acoustic updates, DVL velocity measurement updates are performed at a higher frequency.

[0037] Step S3 includes the following sub-steps: S31, acquire the three-dimensional velocity measurement value in the carrier coordinate system output by the Doppler log and the returned lock status word, identify the current working mode according to the lock status word, the working mode includes bottom tracking mode and water tracking mode; S32, if it is the bottom tracking mode, calculate the rotation matrix from the carrier coordinate system to the navigation coordinate system based on the attitude quaternion in the current nominal state, transform the carrier speed measured by the Doppler log to the navigation coordinate system through the rotation matrix to obtain the navigation speed, and subtract the navigation speed in the nominal state from the navigation speed to obtain the speed measurement residual. The expression for calculating the velocity measurement residual in bottom tracking mode is: ; In the formula, z dvl The velocity measurement residual represents the difference between the velocity observed by DVL and the nominal velocity calculated by the inertial navigation system. The original velocity measurement of DVL in the carrier coordinate system usually includes forward, lateral and vertical velocities; The rotation matrix from the carrier system to the navigation system is calculated based on the quaternions in the current nominal state and is used to project the DVL velocity into the global coordinate system. The nominal velocity is the current global velocity of the submersible calculated from the IMU integration in step S1.

[0038] Velocity measurement residuals and error state δ X A linear relationship exists: ; In the formula, H dvl The DVL measurement matrix maps the 15-dimensional error state vector to a 3-dimensional velocity observation space. In ESKF, this typically takes the form of... This reflects the impact of velocity and rotation errors on the observation, ν bt The measurement noise in the bottom tracking mode has a covariance matrix R. bt The smaller value reflects the high accuracy of multiplying the Kalman gain by the measurement residual DVL at the corresponding time. S33, If it is water tracking mode, the velocity measurement residual calculated by bottom tracking mode is added to the current estimated ocean current velocity estimate to obtain the compensated velocity measurement residual, and the value of the observation noise covariance matrix is ​​replaced with a preset limit value, which is greater than the value of the observation noise covariance matrix in bottom tracking mode. When the submersible is in the middle depth water, exceeding the bottom tracking range of the DVL (Digital Velocity Scale), the DVL measures the velocity relative to the water layer. In this case, an ocean current velocity estimate needs to be introduced for compensation. ; In the formula, To estimate ocean current velocity, in water-tracking mode, DVL measures the velocity relative to the water layer, so the ocean current velocity must be added to obtain the global velocity relative to the ground. The ocean current velocity estimate can be updated in real time through filter state estimation, or long-term calibration can be performed using low-frequency information from the acoustic positioning system.

[0039] Velocity measurement residuals and error state δ X A linear relationship exists: ; In the formula, ν wt The measurement noise in water tracking mode has a covariance matrix R. wt Much larger than R bt This is to reduce the negative contribution of this mode to global accuracy.

[0040] S34. Determine the corresponding velocity measurement matrix based on the velocity measurement residual under the corresponding working mode. Perform Kalman update based on the currently predicted error covariance matrix, velocity measurement matrix and corresponding observation noise covariance matrix to obtain the optimal estimate of the error state, which is used to correct the nominal state, reset the error state and obtain the nominal state after local compensation.

[0041] In this embodiment, by adaptively identifying the DVL working mode and adopting different noise covariance and ocean current compensation strategies, the DVL velocity information can be effectively utilized in both mid-water tracking and near-bottom tracking. The bottom tracking mode provides high-precision velocity constraints, while the water tracking mode, although less accurate, is sufficient to suppress pure inertial navigation divergence. The two modes can be seamlessly switched to ensure the continuity and robustness of the navigation system in complex deep-sea environments.

[0042] S4. Establish a historical state buffer to store the nominal state, error covariance matrix, and raw data of the inertial measurement unit at each time. When an acoustic position measurement with a timestamp is received, calculate its time delay and retrieve the corresponding historical state. Perform Kalman update to obtain the corrected historical state. Then, use the raw data of the inertial measurement unit to reintegrate from the corrected historical state to the current time to obtain the current nominal state, which serves as the final global pose estimate of the submersible.

[0043] Step S4 includes the following sub-steps: S41, Establish a historical state buffer, which stores the nominal state, error covariance matrix and raw data of the inertial measurement unit at each time. S42, when a global position measurement is received from the acoustic positioning system, the total delay is calculated based on the timestamp of the global position measurement and the current system time; It should be noted that, N The buffer capacity is usually taken as... N = tmax / Δ t ,in, t max To maximize the possible acoustic delay; When received with timestamp t obs Acoustic position measurement p acous At that time, calculate the delay: ; In the formula, t For time delay, t now This refers to the current system time. S43, based on the timestamp of the global position measurement, search for the matching index in the historical state buffer, and extract the nominal state and error covariance matrix of the corresponding historical moment; S44. Based on the difference between the global position measurement and the position in the historical nominal state, the position measurement residual is calculated. The position measurement matrix is ​​determined according to the definition of the position measurement residual. Using the error covariance matrix, the position measurement matrix and the noise covariance matrix of the acoustic measurement at the historical moment, Kalman update is performed to obtain the optimal estimate of the error state, which is used to correct the historical nominal state and obtain the corrected historical state. The expression for the position measurement residual is: ; In the formula, p pred ( k τ (This refers to the historical nominal state) X kτ The positional component in the middle.

[0044] The expression for the position measurement matrix is: ; Acoustic measurement noise covariance matrix R a The Kalman gain is multiplied by the measurement residual at the corresponding time point and calibrated by the USBL system.

[0045] Calculate the Kalman gain: ; In the formula, P kτ Here is the error covariance matrix; Historical moment error state estimation: X kτ = K · z acous ; Correcting historical nominal status: ; S45, from the next moment of the historical moment to the current moment, sequentially traverse the raw data of each frame of the inertial measurement unit in the historical state buffer, taking the corrected historical state as the starting point, and recalculate the nominal state at each moment by recursively integrating frame by frame according to the kinematic integral equation of the nominal state, and update the state value of the corresponding moment in the historical state buffer.

[0046] S46, the nominal state obtained at the current moment is used as the final global pose estimate of the submersible.

[0047] This embodiment solves the problem of large time delay in underwater acoustic communication. By using the mechanism of retrospective correction and repropagation, the delayed global position observation is naturally and smoothly integrated into the current state, avoiding trajectory jumps caused by direct correction, and ensuring the physical continuity and global accuracy of state updates. Compared with traditional interpolation or direct delay compensation, this method can make full use of the cached high-frequency IMU data to recover the complete motion trajectory from the historical time to the current time.

[0048] The Kalman update includes the following sub-steps: Calculate the innovation covariance matrix based on the prediction error covariance matrix, measurement matrix, and observation noise covariance matrix at the corresponding time points; the expression is: S=HPH T +R; In the formula, H is the measurement matrix, P is the prediction error covariance matrix at the corresponding time, and R is the observation noise covariance matrix.

[0049] Multiplying the prediction error covariance matrix by the transpose of the measurement matrix, and then by the inverse of the innovation covariance matrix, yields the Kalman gain; the expression is: K = PH T S -1 ; The Kalman gain is multiplied by the measurement residual at the corresponding time to obtain the optimal estimate of the error state. The optimal estimate of the error state is then injected into the nominal state at the corresponding time through manifold addition to correct the position, velocity, attitude and inertial zero bias of the nominal state. The optimal estimate of the error state is expressed as follows: ; In the formula, z This represents the corresponding measurement residual.

[0050] Reset all components of the error state vector to zero. Then, subtract the Kalman gain multiplied by the measurement matrix from the identity matrix, and multiply by the prediction error covariance matrix to obtain the updated posterior error covariance matrix. The expression is: ; This embodiment provides the standard measurement update framework of ESKF, which optimally integrates sensor observations and predicted states, while ensuring that the corrected state satisfies manifold constraints. Error state reset avoids repeated accumulation, and posterior covariance update provides accurate uncertainty input for prediction at the next time step.

[0051] It also includes step S5 for robust fault tolerance before measurement updates, which includes the following sub-steps: S51, Before each measurement update, calculate the Mahalanobis distance based on the new information covariance matrix and the corresponding measurement residuals; the expression is: ; In the formula, d 2 The standardized Mahalanobis distance squared is a scalar that measures the degree to which the current measurement residual deviates from its predicted distribution. Unlike the Euclidean distance, it is normalized by the covariance matrix, eliminating the influence of correlation and dimensions in each dimension.

[0052] S52, compare the Mahalanobis distance with the preset first threshold. If the Mahalanobis distance is less than or equal to the first threshold, determine that the current measurement is a normal measurement, accept the current measurement and perform Kalman update, and update the accelerometer zero bias and gyroscope zero bias of the inertial measurement unit online. S53, if the Mahalanobis distance is greater than the first threshold and less than or equal to the preset second threshold, it is determined that there is strong interference in the current measurement, the observation noise covariance matrix is ​​enlarged proportionally, and Kalman update is performed. The observation noise covariance matrix is ​​scaled up proportionally, and its expression is: ; Performing a Kalman update based on the amplified observation noise covariance matrix reduces the confidence in the measurement and prevents the filter from over-converging to erroneous values.

[0053] S54, if the Mahalanobis distance is greater than the second threshold, the current measurement sensor is determined to be faulty, Kalman update is not performed, and the inertial traverse mode is entered. Using the pre-calibrated accelerometer zero bias and gyroscope zero bias, pure inertial integration is performed based on the angular velocity and specific force output of the inertial measurement unit to maintain the recursive estimation of the submersible's position, velocity and attitude until the sensor signal returns to normal; where the recursive formula is the kinematic equation of the nominal state.

[0054] In this embodiment, during the normal cruise phase before entering the plume region, the ESKF has converged the IMU drift parameters online through high-precision constraints of DVL and acoustics. When the sensor is detected to have failed due to the plume, the system automatically triggers the inertial crossing mode, using the calibrated zero bias to maintain pure inertial integration. Since the zero bias has reached extremely high precision, even within several minutes without any external compensation, the position divergence rate can be maintained at an extremely low level until the submersible passes through the plume region and regains sensor lock.

[0055] This embodiment provides a multi-level anomaly detection and processing mechanism. The Mahalanobis distance takes into account system uncertainty and sensor noise, and can effectively identify strong interferences such as plumes and multipath reflections. The inertial traverse mode utilizes the previously converged IMU zero bias, which can maintain low drift pure inertial navigation for several minutes without any external sensors, greatly enhancing the robustness and safety of the system in extreme environments.

[0056] This embodiment takes a certain type of deep-sea submersible performing a sampling task in a hydrothermal vent area at a depth of 3000 meters as an example; The system first initializes, with the submersible receiving initial GNSS latitude and longitude coordinates from the mother ship at the sea surface. Initial bias values ​​for the IMU's accelerometers and gyroscopes are then calibrated using a static base. The nominal state and error covariance matrix are initialized, and the historical state buffer is cleared.

[0057] During the mid-water descent phase, the submersible begins its descent. Once it exceeds the DVL bottom tracking range, the DVL automatically enters water tracking mode. The system receives IMU data at 200Hz, executes steps S1 and S2, continuously integrates the nominal state, and recursively calculates the covariance. Simultaneously, it receives DVL velocity measurements at 5Hz, executing step S3: introducing real-time estimated ocean current velocity to compensate for the velocity residual, and setting the observation noise covariance R. wt =10 R bt During this stage, the positioning error remains within tens of meters, mainly relying on the IMU; DVL is only used to suppress divergence.

[0058] In the near-bottom approach phase of the bottom tracking mode, when the submersible descends to approximately 100 meters from the bottom, the DVL locks onto the seabed echo and automatically switches to bottom tracking mode. At this point, the observation noise covariance in step S3 is switched to a smaller value. R bt The velocity measurement residual directly reflects the ground velocity. The Kalman update quickly converges to zero bias in the IMU, and the position error is reduced to within the meter level. The submersible begins to cruise in search of hydrothermal vents.

[0059] The surface support vessel provides the absolute position of the submersible every 2 seconds via USBL, but underwater acoustic propagation causes a measurement delay. Each time a USBL measurement is received, the system executes step S4 to calculate the delay, retrieves the state and covariance before the delay from the historical buffer, performs a Kalman update at that historical moment, calculates the Kalman gain, obtains the error estimate, corrects the historical state, and starts from that historical moment to reintegrate and recursively advance to the current moment using the cached IMU data to obtain a smooth and jump-free final pose. When the submersible enters the hydrothermal vent plume region, the DVL beam becomes abnormal due to scattering by suspended particles, and the USBL beam is also affected by multipath interference. In step S5, the Mahalanobis distance detector calculates the Mahalanobis distance in real time. If the Mahalanobis distance is less than or equal to the first threshold, the system updates normally. When the plume causes the residual to increase and the Mahalanobis distance is between the first and second thresholds, the system automatically amplifies the observation noise covariance and still performs the update, but with a reduced weight. If the Mahalanobis distance is greater than the second threshold, the current measurement sensor is deemed to be faulty, the measurement is rejected, and the system enters the inertial traverse mode. Pure inertial integration is performed using only the precisely calibrated IMU zero bias, with a position drift of less than 2 meters per minute. Once the submersible exits the plume region, the sensor signal recovers, the Mahalanobis distance returns to the normal range, and the system automatically exits the inertial traverse mode and continues normal fusion.

[0060] The current nominal status output in step S4 is sent to the navigation control module of the submersible in real time for path tracking, obstacle avoidance, and landing alignment.

[0061] Secondly, the present invention also provides a deep-sea multi-source fusion global navigation system based on time delay backtracking compensation, implemented using a deep-sea multi-source fusion global navigation method based on time delay backtracking compensation, comprising: The nominal state prediction module is used to decompose the actual state of the submersible into a nominal state and an error state, establish the kinematic equation of the nominal state, and use the output of the inertial measurement unit to perform high-frequency integral extrapolation on the nominal state to obtain the predicted value of the nominal state. The error covariance prediction module is used to perform first-order linearization of the kinematic equations based on the nominal state to obtain the error state transition matrix and the noise driving matrix. After discretization, the error covariance matrix of the current time is recursively predicted using the error covariance of the previous time step. The local velocity compensation module is used to acquire the velocity measurement of the Doppler log during the gap in the acoustic position measurement signal transmission, calculate the velocity measurement residual according to its working mode and introduce ocean current compensation, perform Kalman update to correct the nominal state, and obtain the nominal state after local compensation. The global backtracking compensation module is used to establish a historical state buffer, storing the nominal state, error covariance matrix, and raw data of the inertial measurement unit at each moment. When an acoustic position measurement with a timestamp is received, its time delay is calculated and the corresponding historical state is retrieved. Kalman update is performed to obtain the corrected historical state, and the raw data of the inertial measurement unit is used to re-integrate from the corrected historical state to the current moment to obtain the current nominal state, which serves as the final global pose estimate of the submersible.

[0062] It should be noted that this system corresponds to the aforementioned deep-sea multi-source fusion global navigation method based on time delay backtracking compensation. All implementation methods in the above method embodiments are applicable to the embodiments of this system and can achieve the same technical effect.

[0063] Those skilled in the art will recognize that the units and algorithm steps of the various examples described in conjunction with the embodiments disclosed herein can be implemented in electronic hardware, or a combination of computer software and electronic hardware. Whether these functions are implemented in hardware or software depends on the specific application and design constraints of the technical solution. Those skilled in the art can use different methods to implement the described functions for each specific application, but such implementations should not be considered beyond the scope of this invention.

[0064] Those skilled in the art will understand that, for the sake of convenience and brevity, the specific working process of the system and modules described above can be referred to the corresponding process in the foregoing method embodiments, and will not be repeated here.

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

[0066] The units described 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 can be selected to achieve the purpose of this embodiment according to actual needs.

[0067] In addition, the functional units in the various embodiments of the present invention can be integrated into one processing unit, or each unit can exist physically separately, or two or more units can be integrated into one unit.

[0068] If the aforementioned functions are implemented as software functional units and sold or used as independent products, they can be stored in a computer-readable storage medium. Based on this understanding, the technical solution of this invention, essentially, or the part that contributes to the prior art, or a portion of the technical solution, 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 steps 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 USB flash drives, portable hard drives, ROM, RAM, magnetic disks, or optical disks.

[0069] Furthermore, it should be noted that in the system and method of the present invention, it is obvious that the components or steps can be decomposed and / or recombined. These decompositions and / or recombinations should be considered equivalent solutions of the present invention. Moreover, the steps performing the above series of processes can naturally be executed in the order described, but are not necessarily required to be executed in chronological order; some steps can be executed in parallel or independently of each other. Those skilled in the art will understand that all or any step or component of the method and apparatus of the present invention can be implemented in any computing device (including processors, storage media, etc.) or network of computing devices, in hardware, firmware, software, or a combination thereof. This is something that those skilled in the art can achieve by using their basic programming skills after reading the description of the present invention.

[0070] Therefore, the object of the present invention can also be achieved by running a program or a set of programs on any computing system. The computing system can be a known general-purpose system. Therefore, the object of the present invention can also be achieved simply by providing a program product containing program code implementing the method or apparatus. That is, such a program product also constitutes the present invention, and the storage medium storing such a program product also constitutes the present invention. Obviously, the storage medium can be any known storage medium or any storage medium developed in the future. It should also be noted that in the apparatus and method of the present invention, it is obvious that the components or steps can be decomposed and / or recombined. These decompositions and / or recombinations should be considered equivalent to the present invention. Furthermore, the steps performing the above series of processes can naturally be performed in the order described, but are not necessarily required to be performed in chronological order. Some steps can be performed in parallel or independently of each other.

[0071] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A deep-sea multi-source fusion global navigation method based on time delay backtracking compensation, characterized in that, Includes the following steps: S1, decompose the actual state of the submersible into nominal state and error state, establish the kinematic equation of nominal state, and use the output of the inertial measurement unit to perform high-frequency integration extrapolation on nominal state to obtain the predicted value of nominal state. S2, based on the kinematic equation of the nominal state, performs first-order linearization to obtain the error state transition matrix and noise driving matrix. After discretization, the error covariance matrix of the current time is recursively predicted using the error covariance of the previous time step. S3. During the acoustic position measurement signal transmission gap, the velocity measurement of the Doppler log is acquired. The velocity measurement residual is calculated according to its working mode and ocean current compensation is introduced. Kalman update is performed to correct the nominal state and obtain the nominal state after local compensation. S4. Establish a historical state buffer to store the nominal state, error covariance matrix, and raw data of the inertial measurement unit at each time. When an acoustic position measurement with a timestamp is received, calculate its time delay and retrieve the corresponding historical state. Perform Kalman update to obtain the corrected historical state. Then, use the raw data of the inertial measurement unit to reintegrate from the corrected historical state to the current time to obtain the current nominal state, which serves as the final global pose estimate of the submersible.

2. The deep-sea multi-source fusion global navigation method based on time delay backtracking compensation as described in claim 1, characterized in that, Step S1 includes the following sub-steps: S11 defines the true state vector of the submersible, which includes three-dimensional position, three-dimensional velocity, attitude quaternion, three-dimensional zero bias of accelerometer and three-dimensional zero bias of gyroscope. S12 decomposes the actual state of the submersible into a nominal state and an error state, and the nominal state and the error state are combined to form the actual state through manifold addition. S13. Based on the real state vector, establish the kinematic integral equation corresponding to the nominal state. With the sampling period of the inertial measurement unit as the step size and the angular velocity and specific force output by the inertial measurement unit as the driving variables, perform high-frequency integral extrapolation on the nominal state according to the kinematic integral equation to obtain the predicted value of the nominal state at each moment.

3. The deep-sea multi-source fusion global navigation method based on time delay backtracking compensation as described in claim 2, characterized in that, Step S2 includes the following sub-steps: S21, Define the error state vector, which includes three-dimensional position error, three-dimensional velocity error, three-dimensional attitude small angle error, three-dimensional zero bias error of accelerometer and three-dimensional zero bias error of gyroscope; S22, obtain the Jacobian matrix for the kinematic equation of the nominal state to obtain the error state transition matrix in continuous time, and construct the noise driving matrix to map the accelerometer white noise, gyroscope white noise, accelerometer zero-bias random walk noise and gyroscope zero-bias random walk noise to the corresponding error state space. S23, using the sampling period of the inertial measurement unit as the discrete time step, the error state transition matrix is ​​converted into a discrete time state transition matrix using a second-order approximation. Based on the pre-calibrated continuous-time process noise power spectral density matrix, the discrete process noise covariance matrix is ​​obtained by approximation using the noise driving matrix and the sampling interval. S24. Multiply the discrete-time state transition matrix with the error covariance matrix of the previous time step and the transpose of the discrete-time state transition matrix, and add the discrete process noise covariance matrix to obtain the prediction error covariance matrix at the current time step.

4. The deep-sea multi-source fusion global navigation method based on time delay backtracking compensation as described in claim 3, characterized in that, Step S3 includes the following sub-steps: S31, acquire the three-dimensional velocity measurement value in the carrier coordinate system output by the Doppler log and the returned lock status word, identify the current working mode according to the lock status word, the working mode includes bottom tracking mode and water tracking mode; S32, if it is the bottom tracking mode, calculate the rotation matrix from the carrier coordinate system to the navigation coordinate system based on the attitude quaternion in the current nominal state, transform the carrier speed measured by the Doppler log to the navigation coordinate system through the rotation matrix to obtain the navigation speed, and subtract the navigation speed in the nominal state from the navigation speed to obtain the speed measurement residual. S33, If it is water tracking mode, the velocity measurement residual calculated by bottom tracking mode is added to the current estimated ocean current velocity estimate to obtain the compensated velocity measurement residual, and the value of the observation noise covariance matrix is ​​replaced with a preset limit value, which is greater than the value of the observation noise covariance matrix in bottom tracking mode. S34. Determine the corresponding velocity measurement matrix based on the velocity measurement residual under the corresponding working mode. Perform Kalman update based on the currently predicted error covariance matrix, velocity measurement matrix and corresponding observation noise covariance matrix to obtain the optimal estimate of the error state, which is used to correct the nominal state and obtain the nominal state after local compensation.

5. The deep-sea multi-source fusion global navigation method based on time delay backtracking compensation as described in claim 4, characterized in that, Step S4 includes the following sub-steps: S41, Establish a historical state buffer, which stores the nominal state, error covariance matrix and raw data of the inertial measurement unit at each time. S42, when a global position measurement is received from the acoustic positioning system, the total delay is calculated based on the timestamp of the global position measurement and the current system time; S43, based on the timestamp of the global position measurement, search for the matching index in the historical state buffer, and extract the nominal state and error covariance matrix of the corresponding historical moment; S44. Based on the difference between the global position measurement and the position in the historical nominal state, the position measurement residual is calculated. The position measurement matrix is ​​determined according to the definition of the position measurement residual. Using the error covariance matrix, the position measurement matrix and the noise covariance matrix of the acoustic measurement at the historical moment, Kalman update is performed to obtain the optimal estimate of the error state, which is used to correct the historical nominal state and obtain the corrected historical state. S45, from the next moment of the historical moment to the current moment, traverse the original data of each frame of the inertial measurement unit in the historical state buffer in turn. Starting from the corrected historical state, recalculate the nominal state at each moment by integrating frame by frame according to the kinematic integral equation of the nominal state, and update the state value of the corresponding moment in the historical state buffer. S46, the nominal state obtained at the current moment is used as the final global pose estimate of the submersible.

6. The deep-sea multi-source fusion global navigation method based on time delay backtracking compensation as described in claim 5, characterized in that, The Kalman update includes the following sub-steps: Calculate the innovation covariance matrix based on the prediction error covariance matrix, measurement matrix, and observation noise covariance matrix at the corresponding time. The Kalman gain is obtained by multiplying the prediction error covariance matrix by the transpose of the measurement matrix and then by the inverse of the innovation covariance matrix. The Kalman gain is multiplied by the measurement residual at the corresponding time to obtain the optimal estimate of the error state. The optimal estimate of the error state is then injected into the nominal state at the corresponding time through manifold addition to correct the position, velocity, attitude and inertial zero bias of the nominal state. The error state vector is reset to zero. The updated posterior error covariance matrix is ​​obtained by subtracting the Kalman gain multiplied by the measurement matrix from the identity matrix and then multiplying it by the prediction error covariance matrix.

7. The deep-sea multi-source fusion global navigation method based on time delay backtracking compensation as described in claim 6, characterized in that, It also includes step S5 for robust fault tolerance before measurement updates, which includes the following sub-steps: S51, Before each measurement update, calculate the Mahalanobis distance based on the new information covariance matrix and the corresponding measurement residuals; S52, compare the Mahalanobis distance with the preset first threshold. If the Mahalanobis distance is less than or equal to the first threshold, determine that the current measurement is a normal measurement, accept the current measurement and perform Kalman update, and update the accelerometer zero bias and gyroscope zero bias of the inertial measurement unit online. S53, if the Mahalanobis distance is greater than the first threshold and less than or equal to the preset second threshold, it is determined that there is strong interference in the current measurement, the observation noise covariance matrix is ​​enlarged proportionally, and Kalman update is performed. S54. If the Mahalanobis distance is greater than the second threshold, the current measurement sensor is determined to be faulty. Kalman update is not performed, and the system enters the inertial traverse mode. Using the pre-calibrated accelerometer zero bias and gyroscope zero bias, pure inertial integration is performed based on the angular velocity and specific force output of the inertial measurement unit to maintain the recursive estimation of the submersible's position, velocity, and attitude until the sensor signal returns to normal.

8. A deep-sea multi-source fusion global navigation system based on time delay backtracking compensation, implemented using the deep-sea multi-source fusion global navigation method based on time delay backtracking compensation as described in any one of claims 1-7, characterized in that, include: The nominal state prediction module is used to decompose the actual state of the submersible into a nominal state and an error state, establish the kinematic equation of the nominal state, and use the output of the inertial measurement unit to perform high-frequency integral extrapolation on the nominal state to obtain the predicted value of the nominal state. The error covariance prediction module is used to perform first-order linearization of the kinematic equations based on the nominal state to obtain the error state transition matrix and the noise driving matrix. After discretization, the error covariance matrix of the current time is recursively predicted using the error covariance of the previous time step. The local velocity compensation module is used to acquire the velocity measurement of the Doppler log during the gap in the acoustic position measurement signal transmission, calculate the velocity measurement residual according to its working mode and introduce ocean current compensation, perform Kalman update to correct the nominal state, and obtain the nominal state after local compensation. The global backtracking compensation module is used to establish a historical state buffer, storing the nominal state, error covariance matrix, and raw data of the inertial measurement unit at each moment. When an acoustic position measurement with a timestamp is received, its time delay is calculated and the corresponding historical state is retrieved. Kalman update is performed to obtain the corrected historical state, and the raw data of the inertial measurement unit is used to re-integrate from the corrected historical state to the current moment to obtain the current nominal state, which serves as the final global pose estimate of the submersible.

9. A terminal device, characterized in that, The terminal device includes: a memory, a processor, and a program for a deep-sea multi-source fusion global navigation method based on time delay backtracking compensation stored in the memory and executable on the processor. The program for the deep-sea multi-source fusion global navigation method based on time delay backtracking compensation is configured to implement the steps of the deep-sea multi-source fusion global navigation method based on time delay backtracking compensation as described in any one of claims 1 to 7.

10. A computer-readable storage medium, characterized in that, The storage medium stores a deep-sea multi-source fusion global navigation method and system program based on time delay backtracking compensation. When the program is executed, it implements the deep-sea multi-source fusion global navigation method based on time delay backtracking compensation as described in any one of claims 1-7.

Images