A high-precision positioning method, positioning system and device for a seabed reference station

Abstract

The application belongs to the technical field of underwater acoustic data positioning and identification, and discloses a high-precision positioning method, a positioning system and equipment for a seabed reference station, the positioning method comprising: introducing a variable cubic B-spline function to represent a disturbance sound speed time domain based on a variable sound speed field time domain and a variable sound speed structure vertical direction, constructing a step-by-step iteration-gradual correction position information and sound speed information iteration inversion model based on sound ray tracing, and gradually correcting seabed reference station coordinates and disturbance sound speed to obtain high-precision position information of the seabed reference station. The application gradually corrects seabed reference station coordinates and disturbance sound speed, effectively eliminates sound speed related systematic errors caused by sound speed structure changes, better depicts disturbance sound speed details by using the good local support of the cubic B-spline function, and effectively improves the positioning accuracy of the seabed reference station.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of underwater acoustic data positioning and identification, and particularly relates to a high-precision positioning method, a positioning system and equipment for a seabed reference station. BACKGROUND

[0002] At present, marine spatial position information is an important basis for all marine activities. As an important marine spatial information infrastructure of a country, the construction of a seabed geodetic control network puts forward higher requirements for the positioning of seabed reference stations. At present, GNSS / A (Global Navigation Satellite System / Acoustics) combined technology is widely used for positioning seabed reference stations.

[0003] Due to the wide range, complexity and variability of the marine natural environment, the sound velocity distribution in the ocean is complex and variable, and the sound velocity is one of the most significant factors restricting the positioning accuracy of underwater acoustics. Therefore, in high-precision underwater acoustic positioning, the sound velocity error needs to be properly modeled to improve the positioning accuracy. At present, sound velocity error modeling generally considers fitting the time-domain variation of the disturbed sound velocity structure by a polynomial, such as a quadratic polynomial. However, the physical background of a simple polynomial is weak, and the coherence before and after the change of the sound velocity structure is poor, which cannot well reflect the change of the sound velocity. In view of the above problems, in order to better model the sound velocity error and improve the positioning accuracy of underwater acoustics, it is necessary to study a modeling method that is more in line with the actual sound velocity variation.

[0004] Through the above analysis, the problems and defects of the prior art are that the sound velocity error modeling cannot well fit the actual situation, resulting in low positioning accuracy of seabed reference points. SUMMARY

[0005] In view of the problems existing in the prior art, the application provides a high-precision positioning method, a positioning system and equipment for a seabed reference station, and particularly relates to a high-precision positioning method for a seabed reference station considering the time-domain variation of the sound velocity structure.

[0006] The application is implemented as follows: a high-precision positioning method for a seabed reference station considering the time-domain variation of the sound velocity structure comprises the following steps:

[0007] Based on the time-domain variation of the sound velocity field and the vertical direction of the sound velocity structure, a variable cubic B-spline function is introduced to represent the time-domain variation of the disturbed sound velocity, and based on ray tracing, a step-by-step iteration-gradual correction position information and sound velocity information iteration inversion model is constructed, and the seabed reference station coordinates and the disturbed sound velocity are gradually corrected to obtain high-precision position information of the seabed reference station.

[0008] Further, the high-precision positioning method for a seabed reference station considering the time-domain variation of the sound velocity structure further comprises the following steps:

[0009] Step one, the survey ship implements continuous observation on the seafloor reference station, and obtains the survey ship track, the seafloor reference station position and the approximate position of the sound velocity profile measurement;

[0010] Step two, the collected sound velocity profile is preprocessed, and the disturbed sound velocity structure is inverted by using the cubic B-spline function method;

[0011] Step three, the seafloor reference station position and the disturbed sound velocity structure are iteratively inverted, and the gradual correction of the seafloor reference station coordinates and the sound velocity structure is realized.

[0012] Further, the step one of measuring the survey ship on the seafloor reference station to implement continuous observation, and obtaining the survey ship track, the seafloor reference station position and the approximate position of the sound velocity profile measurement comprises the following steps:

[0013] (1) The survey ship sails on the sea surface, continuously observes the seafloor reference station, collects sonar measurement data, sound velocity profile data and GNSS data, etc., and obtains the survey ship track and the sound velocity profile measurement data;

[0014] (2) Layered isogradient ray tracing is implemented, and the approximate coordinates of the seafloor point are inversely calculated.

[0015] Further, the step two of pre-processing the collected sound velocity profile, and inverting the disturbed sound velocity structure by using the cubic B-spline function method comprises the following steps:

[0016] 1) The collected sound velocity profile is preprocessed, including sound velocity measurement outlier rejection, partial sound velocity profile data extension, adjacent isogradient sound velocity layer merging, and calculation of the average sound velocity profile and the difference between each sound velocity profile and the average sound velocity profile;

[0017] 2) The cubic B-spline function is introduced to fit the sound velocity structure change.

[0018] Further, the step three of iteratively inverting the seafloor reference station position and the disturbed sound velocity structure to realize the gradual correction of the seafloor reference station coordinates and the sound velocity structure comprises the following steps:

[0019] (i) The seafloor station position is inverted by using the fixed sound velocity-position solving method, the ray tracing nonlinear system is linearized at the approximate coordinates of the seafloor reference station according to the Taylor series linearization principle, and the seafloor reference station position correction is calculated based on the robust Bayesian least square estimation principle;

[0020] (ii) The sound speed structure inversion is performed by using the fixed position-solve sound speed method. The actual sound speed structure is defined as the coupling of the reference sound speed structure and the perturbed sound speed structure. The cubic B-spline function is introduced to fit the sound speed structure variation. The ray tracing nonlinear system is linearized based on the reference sound speed profile. The sound speed structure time variation is inverted based on the robust Bayesian least square estimation principle.

[0021] (iii) The sea bottom reference station coordinates and the perturbed sound speed are gradually corrected.

[0022] Further, the sea bottom station position inversion in step (i) by using the fixed sound speed-solve position method comprises:

[0023] (i-1) The ray tracing nonlinear system is linearized at the sea bottom reference station approximate coordinates based on the Taylor series linearization principle. Specifically, the linearization is performed as follows:

[0024] t r =f 0 +A1·dX tp

[0025] wherein f 0 is the propagation time obtained by ray tracing based on the initial coordinates of the sea bottom reference station; dX tp is the sea bottom reference station position correction, A1 is the coefficient matrix of the propagation time with respect to the sea bottom reference station position, wherein the derivation with respect to the plane position can be obtained based on the ray acoustics principle according to the following formula, and the derivation with respect to the z direction can be obtained by numerical analysis method.

[0026]

[0027] wherein p is the Snell constant determined by the ray tracing process, i.e. the sound speed horizontal slowness; r is the plane distance between the sound source point and the sea bottom reference station,

[0028] (i-2) The sea bottom reference station position correction is calculated based on the robust Bayesian least square estimation principle. Specifically, the calculation is performed as follows:

[0029]

[0030] wherein P is the sonar observation value weight matrix; P x is the model parameter prior weight matrix; if the model parameter prior information is not considered, i.e. P x = 0, the method degenerates into the classical least square method; L = t r -f 0 is the observation vector. is the prior value of the model parameters. It is noted that if there are multiple stations for synchronous solution, the coordinates of the seafloor reference station can be solved by looping through the multiple stations.

[0031] Further, the fixing position-solving sound velocity in step (ii) for the sound velocity structure inversion includes:

[0032] (ii-1) The actual sound velocity structure is defined as the coupling of the reference sound velocity structure and the perturbed sound velocity structure. Considering the time-domain variation of the sound velocity field, the sound velocity structure is expressed as follows:

[0033] V(z, t e ) = V r (z) + V p (z, t e )

[0034] In the formula, represents the observation time sequence, V represents the sound velocity structure ignoring horizontal spatial heterogeneity, which is related to depth and observation time; V r represents the reference sound velocity structure, which can be a measured sound velocity profile or an average profile of multiple measured sound velocity profiles; V p represents the perturbed sound velocity structure;

[0035] (i1-2) A cubic B-spline function is introduced to fit the sound velocity structure variation, while considering the difference in sound velocity structure variation between the upper and lower layers of seawater, i.e. considering the smooth and weak characteristics of the sound velocity variation in the deep water layer within the sonar observation time window, the perturbed sound velocity is set as a piecewise function form according to the empirical depth. The entire observation time sequence is divided into m segments according to a certain time interval Δt. If the observation time t e belongs to the i-th segment, its perturbed sound velocity structure is as follows:

[0036]

[0037] In the formula, a represents the B-spline function coefficient; B represents the B-spline basis function; d represents the B-spline degree, which is taken as d = 3 in this paper; z0 represents the reference depth, which is generally taken as 1000-1600m for deep water sea areas; u i is the unitized observation time, and the calculation formula is u i = (t e -t i ) / Δt, Δt = t i+1 -t i . t i and t i+1 are the start and end observation times of the i-th segment, respectively;

[0038] (ii-3) The sound ray tracing nonlinear system is linearized based on the reference sound velocity profile, and the following formula is obtained:

[0039] tr = f 0 + A2 da

[0040] wherein, The derivative formula of the sound speed disturbance with respect to the spline coefficient is derived from the derivative principle of the cubic B-spline function; the derivative of the sound signal propagation time with respect to the disturbance sound speed can refer to the numerical calculation formula of the derivative of the discrete function, as shown in the following formula:

[0041]

[0042] wherein, k1 and k2 respectively represent the indexes of the sound speed profile layers where the sound source point and the seabed reference station are located; c k represents the incident sound speed value in the kth layer; g k represents the sound speed gradient value in the kth layer;

[0043] (ii-4) Based on the robust Bayesian least square estimation principle, the time domain variation of the sound speed structure is inverted. It should be noted that when multiple time windows are calculated, the sound speed structure inversion of the multiple windows can be completed through a loop.

[0044] Further, in step (iii), the seabed reference station coordinates and the disturbance sound speed are gradually corrected, including:

[0045] (iii-1) The sound speed structure variation obtained from the sound speed structure inversion of the fixed position-sound speed is fed back to the seabed reference station coordinate calculation, and the seabed reference station position inversion is performed again;

[0046] (iiii-2) When the coordinate increment and the sound speed variation are within the preset threshold range, the calculation is ended, and the gradual correction of the seabed reference station coordinates and the disturbance sound speed is completed.

[0047] Another object of the present application is to provide a seabed reference station high-precision positioning system considering the time domain variation of the sound speed structure, including:

[0048] The measurement approximate position acquisition module is used for measuring the continuous observation of the ship on the seabed reference station, and obtaining the measurement ship track, the seabed reference station position and the sound speed profile measurement approximate position.

[0049] The disturbance sound speed structure inversion module is used for data preprocessing of the collected sound speed profile, and the disturbance sound speed structure inversion is performed by using the cubic B-spline function method.

[0050] The reference station coordinate and sound speed structure gradual correction module is used for iterative inversion of the seabed reference station position and the disturbance sound speed structure, so that the gradual correction of the seabed reference station coordinates and the sound speed structure is realized.

[0051] Another object of the present application is to provide an apparatus characterized in that the apparatus comprises a memory and a processor, the memory stores a program, and the program is executed by the processor to enable the processor to perform the high-precision positioning method of a seabed reference station considering time-domain variation of sound velocity structure.

[0052] In combination with the above technical solutions and the technical problems solved, the advantages and positive effects of the technical solutions to be protected by the present application are analyzed from the following aspects:

[0053] Firstly, in view of the technical problems existing in the prior art and the difficulty in solving the problems, the technical solutions to be protected by the present application and the results and data in the research and development process are combined to analyze in detail and profoundly how the technical solutions solve the technical problems and bring some creative technical effects after solving the problems.

[0054] The difficulty of the defects existing in the prior art is as follows: one is how to represent the time-domain variation of complex disturbance sound velocity by using appropriate mathematical methods; and the other is how to correctly process the time-domain variation of sound velocity structure, establish a suitable calculation model, and calculate the coordinates of the seabed reference point and the variation of sound velocity.

[0055] The high-precision positioning method of a seabed reference station considering time-domain variation of sound velocity structure provided by the present application solves the problems existing in the background art, and brings the following significance: by constructing an iterative estimation model of position information and sound velocity information, gradually correcting the coordinates of the seabed reference station and the disturbance sound velocity, reducing the influence of sound velocity error on positioning accuracy, and obtaining high-precision position information of the seabed reference station.

[0056] Secondly, from the perspective of the product as a whole, the technical solutions to be protected by the present application have the following technical effects and advantages:

[0057] The present application considers the time-domain variation of the sound velocity field and the variation characteristics of the sound velocity structure in the vertical direction, introduces a cubic B-spline function to represent the time-domain variation of the disturbance sound velocity, constructs an iterative inversion model of position information and sound velocity information based on the ray tracing theory, and gradually corrects the coordinates of the seabed reference station and the disturbance sound velocity, so as to obtain high-precision position information of the seabed reference station.

[0058] Compared with the prior art, the present application has the following advantages:

[0059] The application combines the actual underwater positioning, constructs an iterative estimation model of position information and sound velocity information, and provides a high-precision positioning method of seabed reference point considering the time-domain variation of sound velocity structure.

[0060] Thirdly, the creativity of the application as the claim is also embodied in the following important aspects:

[0061] Whether the technical solution of the application solves the technical problems that people have been eager to solve but have failed to obtain success:

[0062] Sound velocity error is the core factor restricting the high-precision positioning of seabed reference station, especially influenced by complex and dynamic marine environment, the sound velocity structure has obvious time-domain variation characteristics, how to characterize and process the sound velocity disturbance is the key problem to be solved for realizing the high-precision positioning of seabed reference station. The technical solution of the application considers the time-domain variation of sound velocity structure, characterizes the sound velocity disturbance details based on cubic B-spline function, and processes the sound velocity structure error through iterative inversion, so as to weaken the influence of marine sound velocity disturbance error and realize the high-precision positioning of seabed reference station. BRIEF DESCRIPTION OF DRAWINGS

[0063] Figure 1 It is a high-precision positioning method of seabed reference station considering the time-domain variation of sound velocity structure provided by the embodiment of the application;

[0064] Figure 2 It is a high-precision positioning method of seabed reference station considering the time-domain variation of sound velocity structure provided by the embodiment of the application;

[0065] Figure 3 It is a schematic diagram of underwater acoustic positioning based on a survey ship provided by the embodiment of the application;

[0066] Figure 4 It is a high-precision positioning method of seabed reference station considering the time-domain variation of sound velocity structure provided by the embodiment of the application;

[0067] Fig. 5 is an average sound velocity profile and sound velocity variation value schematic diagram provided by the embodiment of the application;

[0068] Fig. 6 is a sound velocity correction schematic diagram based on cubic B-spline function method provided by the embodiment of the application;

[0069] Figure 7 This is a schematic diagram illustrating the results of calculating the location of the seabed station using the traditional and new methods provided in this embodiment of the invention.

[0070] Figure 8 This is a schematic diagram of the statistical results of the observation residuals of the traditional method and the new method provided in the embodiments of the present invention;

[0071] Figure 9 This is a schematic diagram of a high-precision positioning system for a seabed reference station that takes into account the time-domain variation of sound velocity structure provided in this embodiment of the invention;

[0072] Another objective of this invention is to provide a high-precision positioning system for a seabed reference station that takes into account the time-domain variations of sound velocity structure, comprising:

[0073] The approximate location acquisition module is used by the survey vessel to conduct continuous observations of the seabed reference station and acquire the survey vessel's track, the location of the seabed reference station, and the approximate location of the sound velocity profile.

[0074] The disturbed sound velocity structure inversion module is used to preprocess the collected sound velocity profiles and perform disturbed sound velocity structure inversion using the cubic B-spline function method.

[0075] The module for progressively correcting the coordinates and sound velocity structure of the reference station is used to iteratively invert the location of the seabed reference station and the disturbed sound velocity structure, thereby realizing the progressive correction of the coordinates and sound velocity structure of the seabed reference station. Detailed Implementation

[0076] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the invention.

[0077] I. Explanatory and Illustrative Embodiments. To enable those skilled in the art to fully understand how the present invention is specifically implemented, this section provides an explanatory and illustrative description of the embodiments described in the claims.

[0078] like Figure 1 As shown, the high-precision positioning method for seabed reference stations that takes into account the time-domain variation of sound velocity structure provided by the present invention includes the following steps:

[0079] S101: Starting from the theory of ray acoustics, based on the intrinsic direct sound ray model of the sound source point and the seabed reference station, using the sound ray propagation time as the observable, and following Snell's law to implement layered equal gradient sound ray tracking, and inverting and solving the coordinates of the seabed point.

[0080] S102: Construct an iterative inversion model for position and sound velocity information with "step-by-step iteration and gradual correction";

[0081] S103: gradually correct the seafloor reference station coordinates and the perturbed sound speed until the solution meets the threshold requirement, so as to obtain the seafloor reference station position and the time-domain variation of the sound speed structure.

[0082] Embodiment 1

[0083] As Figure 2 shown, the principle of the high-precision positioning method of the seafloor reference station considering the time-domain variation of the sound speed structure provided by the embodiment of the application includes:

[0084] First, based on the initial sound speed profile, the seafloor reference station coordinate tracking and position inversion are performed in combination with the basic principle of ray tracing, the sound speed perturbation model is constructed according to the residual information, the sound speed structure inversion is implemented, and it is judged whether the coordinate inversion result and the sound speed structure inversion information meet the preset threshold value. If not, the sound speed information is updated, the seafloor reference station coordinate inversion is performed again, and the sound speed structure inversion is continued. The threshold value is judged again until the threshold value is met, the iterative inversion is ended, and the seafloor station position information is output.

[0085] Among them, Figure 3 is a schematic diagram of the underwater acoustic positioning based on a survey ship provided by the embodiment of the application;

[0086] Embodiment 2

[0087] As Figure 4 shown, the high-precision positioning method of the seafloor reference station considering the time-domain variation of the sound speed structure provided by the application includes the following steps:

[0088] S1, starting from the ray acoustics theory, based on the intrinsic direct sound ray model of the sound source point and the seafloor reference station, taking the sound ray propagation time as the observation quantity, implementing layered equal-gradient sound ray tracing in accordance with Snell's law, and inversely solving the seafloor point coordinates, the specific process includes the following steps:

[0089] S1.1, the survey ship sails on the sea surface, continuously observes the seafloor reference station, collects sonar measurement time data, sound speed profile data, GNSS data and the like, and obtains the survey ship track and sound speed profile measurement (the last three measurements) and the like.

[0090] S1.2, data preprocessing is performed on the collected sound speed profile, including sound speed measurement outlier rejection, partial sound speed profile data extrapolation, and adjacent same-gradient sound speed layer merging, and the average sound speed profile and the difference between each sound speed profile and the average sound speed profile are calculated; as shown in FIG. 5.

[0091] S1.3, layered equal-gradient sound ray tracing is implemented, and the seafloor point approximate coordinates are inversely solved.

[0092] S2, a position information and sound speed information iterative inversion model of "step-by-step iteration-gradual correction" is constructed, and the specific process includes the following steps:

[0093] S2.1, First, the "fixed sound speed-solve position" method is used to perform seabed station position inversion:

[0094] First, according to the Taylor series linearization principle, the sound ray tracking nonlinear system is linearized at the approximate coordinates of the seabed reference station, specifically:

[0095] t r =f 0 +A1·dX tp

[0096] In the formula, f 0 is the propagation time obtained by sound ray tracking based on the initial coordinates of the seabed reference station ; dX tp is the seabed reference station position correction, A1 is the coefficient matrix of the propagation time with respect to the seabed reference station position, wherein the derivation with respect to the plane position can be obtained based on the basic principles of ray acoustics according to the following formula, and the derivation with respect to the z direction can be obtained by numerical analysis method.

[0097]

[0098] In the formula, p is the Snell constant determined by the sound ray tracking process, that is, the sound speed horizontal slowness; r is the plane distance between the sound source point and the seabed reference station,

[0099] Second, the seabed reference station position correction is calculated based on the robust Bayesian least squares estimation principle, specifically:

[0100]

[0101] In the formula, P is the sonar observation weight matrix; P x is the model parameter prior weight matrix; if the model parameter prior information is not considered, that is, P x = 0, the method degenerates into the classical least squares method; L = t r -f 0 is the observation vector; is the model parameter prior value. It should be pointed out that if there are multiple stations for synchronous solution, the seabed reference station coordinate solution of multiple stations can be completed through a loop.

[0102] S2.2, secondly, the "fixed position-solve sound speed" method is used to perform sound speed structure inversion:

[0103] The first step is to define the actual sound speed structure as the coupling of the reference (background) sound speed structure and the perturbation sound speed structure. Considering the time-domain variation of the sound speed field, the sound speed structure can be expressed as follows:

[0104] V(z, t e ) = V r (z) + V p (z, t e )

[0105] In the formula, represents the observation time sequence, V represents the sound speed structure ignoring horizontal spatial heterogeneity, which is related to depth and observation time; V r represents the reference sound speed structure, which can be a measured sound speed profile or an average profile of multiple measured sound speed profiles; V p represents the perturbation sound speed structure.

[0106] The second step is to introduce a cubic B-spline function to fit the sound speed structure variation, as shown in FIG. 6, while considering the difference in sound speed structure variation between the upper and lower layers of seawater, i.e., considering the smooth and weak characteristics of the sound speed variation in the deep water layer within the sonar observation time window. The perturbation sound speed is set as a segmented function according to the empirical depth. The entire observation time sequence is divided into m segments according to a certain time interval Δt. If the observation time t e belongs to the i-th segment, the perturbation sound speed structure is as follows:

[0107]

[0108] In the formula, a represents the B-spline function coefficient; B represents the B-spline basis function; d represents the B-spline degree, which is taken as d = 3 in this paper; z0 represents the reference depth, which is generally empirically taken as 1000-1600 m for deep water sea areas; u i is the unitized observation time, and the calculation formula is u i = (t e -t i ) / Δt, Δt = t i+1 -t i . t i and t i+1 are the start and end observation times of the i-th segment, respectively.

[0109] The third step is to linearize the ray tracing nonlinear system based on the reference sound speed profile, which can be obtained as follows:

[0110] t r = f 0 +A2·da

[0111] In the formula, The derivative formula of the sound speed perturbation with respect to the spline coefficient can be obtained from the derivative principle of the cubic B-spline function; the derivative of the sound signal propagation time with respect to the perturbed sound speed can refer to the numerical calculation formula of the derivative of a discrete function, as shown in the following formula.

[0112]

[0113] In the formula, k1 and k2 respectively represent the indexes of the sound speed profile layers where the sound source point and the seabed reference station are located; c k represents the incident sound speed value in the kth layer; g k represents the sound speed gradient value in the kth layer;

[0114] In the fourth step, the time-domain variation of the sound speed structure is inverted based on the robust Bayesian least square estimation principle. It should be noted that when multiple time windows are used for calculation, the sound speed structure inversion of the multiple windows can be completed through a loop.

[0115] In S3, the seabed reference station coordinates and the perturbed sound speed are gradually corrected until the calculation result meets the threshold requirement, so as to obtain the seabed reference station position and the time-domain variation of the sound speed structure. The specific process includes the following steps:

[0116] In the first step, the sound speed structure variation obtained from S2.2 is fed back to the seabed reference station coordinate calculation, and the seabed reference station position inversion is performed again;

[0117] In the second step, when the coordinate increment and the sound speed variation are within the preset threshold range, the calculation is ended, and the gradual correction of the seabed reference station coordinates and the perturbed sound speed is completed.

[0118] Embodiment 3

[0119] As Figure 9 shown, the embodiment of the present application provides a high-precision seabed reference station positioning system considering the time-domain variation of the sound speed structure, which comprises:

[0120] A measurement approximate position acquisition module 1 is used for measuring the continuous observation of a ship on a seabed reference station, and acquiring a measurement ship track, a seabed reference station position and a sound speed profile measurement approximate position.

[0121] A perturbed sound speed structure inversion module 2 is used for data preprocessing of the collected sound speed profile, and performing perturbed sound speed structure inversion by using a cubic B-spline function method.

[0122] A reference station coordinate and sound speed structure gradual correction module 3 is used for iterative inversion of the seabed reference station position and the perturbed sound speed structure, so as to realize the gradual correction of the seabed reference station coordinates and the sound speed structure.

[0123] II. Application Examples. In order to prove the creativity and technical value of the technical solutions of the present application, this part is the application examples of the technical solutions of the claims on specific products or related technologies.

[0124] Application Examples

[0125] The embodiment of the present application provides a high-precision positioning method of a seabed reference station considering time-domain changes of sound velocity structure, comprising the following steps:

[0126] Starting from the ray acoustics theory, based on the intrinsic direct sound ray model of a sound source point and a seabed reference station, taking the sound ray propagation time as an observation, implementing layered equal-gradient sound ray tracking in accordance with the Snell law, and inversely solving the coordinates of the seabed point, the essence of which follows the following rules:

[0127] |t obs -t r |=min

[0128] In the formula, t obs represents the actual measured sound signal propagation time sequence, t r represents the sound signal propagation time sequence obtained by sound ray tracking according to the ray acoustics principle, and the relationship of the sound source point, the seabed reference station coordinates and the sound velocity structure is as follows:

[0129] t r =f(X td ,X tp ,V)

[0130] In the formula, f represents a nonlinear mapping relationship of sound ray tracking, X td =(x td ,y td ,z td ) represents three-dimensional coordinates of the sound source point, X tp =(x tp ,y tp ,z tp ) represents three-dimensional coordinates of the seabed reference station, and V represents the sound velocity structure.

[0131] In a preferred embodiment, a “step-by-step iteration-gradual correction” inversion model based on robust Bayesian least squares is constructed, that is, “fixed sound velocity-solve position” and “fixed position-estimate sound velocity” are alternately iterated, the seabed station coordinates and the sound velocity structure are gradually corrected, until the solution result meets the threshold requirement, so that the position of the seabed reference station and the time-domain changes of the sound velocity structure are obtained.

[0132] In a preferred embodiment, the “step-by-step iteration-gradual correction” inversion model is established and solved as follows:

[0133] (1) Seabed station position inversion

[0134] The position inversion of the seabed reference station mainly employs a "fixed sound velocity - solution for position" method. Based on the Taylor series linearization principle, the nonlinear sound ray tracking system is linearized at the approximate coordinates of the seabed reference station, yielding the following equation:

[0135] t r =f 0 +A1·dX tp

[0136] In the formula, f 0 Based on the initial coordinates of the seabed reference station The propagation time obtained by sound ray tracking; dX tp This is the correction amount for the position of the seabed reference station. A1 is the coefficient matrix of the propagation time with respect to the position of the seabed reference station. The derivative with respect to the planar position can be obtained based on the fundamental principles of ray acoustics, and the derivative with respect to the z-direction can be obtained by numerical analysis.

[0137]

[0138] In the formula, p is the Snell constant determined by the ray tracking process, i.e., the horizontal slowness of the sound speed; r is the planar distance between the sound source point and the seabed reference station.

[0139] Based on the robust Bayesian least squares estimation principle, the solution for the seabed reference station position correction is obtained as follows:

[0140]

[0141] In the formula, P is the sonar observation weight matrix; P x Let P be the prior weight matrix of the model parameters; if we do not consider the prior information of the model parameters, i.e., P... x =0, the method degenerates into the classical least squares method; L=t r -f 0 For observation vectors; These are the prior values ​​for the model parameters. It should be noted that if multiple stations are performing simultaneous calculations, the coordinates of the seabed reference stations for multiple stations can be obtained by iterating through the data.

[0142] (2) Sound velocity structure inversion

[0143] The sound velocity structure inversion mainly employs a "fixed position - solve for sound velocity" approach. The actual sound velocity structure is a coupling of the reference (background) sound velocity structure and the perturbed sound velocity structure. Considering the time-domain variation of the sound velocity field, the sound velocity structure can be expressed as follows:

[0144] V(z,t e ) = V r (z)+Vp (z, t e )

[0145] wherein, denotes the sequence of observation time, V denotes the sound speed structure ignoring horizontal spatial heterogeneity, which is related to depth and observation time; V r denotes the reference sound speed structure, which can be a certain actual measured sound speed profile or an average profile of multiple measured sound speed profiles; V p denotes the perturbed sound speed structure. It is generally considered that the sound speed structure changes the same from the sea surface layer to the seabed layer, so V p can be expressed as a quadratic polynomial or a cubic polynomial with respect to the observation time within a certain time window, as shown in the following formula.

[0146] V p (z, t e ) = a0 + a1(t e -t0) + a2(t e -t0) 2

[0147] wherein, t e is the observation time, a i (i = 0, 1, 2) is the model parameter, and t0 is the reference time. In this method, there is no correlation between the sound speed structures corresponding to the discrete observation times, so the smoothness of the sound speed structure change is reduced. It should be noted that considering the time and spatial range of the underwater acoustic positioning operation, the sound speed profile change is mainly affected by internal waves and submesoscale ocean dynamic processes, so the window setting can be performed in steps, that is, a larger window is set first, and then a smaller window is set, and the long and short period perturbed sound speed is corrected in turn.

[0148] In a preferred embodiment, considering the characteristics of the time-domain continuous change of the sound speed structure in a broad sense, a cubic B-spline function is introduced to fit the sound speed structure change. In addition, the present application takes into account the difference in sound speed structure change between the upper and lower layers of seawater, that is, considering the characteristics of the smooth and weak change of the sound speed in the deep water layer within the sonar observation time window, the perturbed sound speed is set as a segmented function according to the empirical depth. The entire sequence of observation times is divided into m segments according to a certain time interval Δt. If the observation time t e belongs to the i-th segment, then its perturbed sound speed structure is as follows:

[0149]

[0150] wherein, a denotes the B-spline function coefficient; B denotes the B-spline basis function; d denotes the B-spline degree, which is taken as d = 3 in this paper; z0 denotes the reference depth, which is generally empirically taken as 1000-1600m for deep water sea areas; u i is the unitized observation time, and the calculation formula is u i =(te -t i ) / Δt, Δt=t i+1 -t i .t i , t i+1 are the initial and final observation time of the i-th segment, respectively.

[0151] Similarly, the ray tracing nonlinear system is linearized based on the reference sound speed profile, and the following formula is obtained:

[0152] t r =f 0 +A2·da

[0153] In the formula, The derivative formula of the sound speed disturbance with respect to the spline coefficient can be obtained from the derivative principle of the cubic B-spline function; the derivative of the sound signal propagation time with respect to the disturbed sound speed can refer to the numerical calculation formula of the derivative of the discrete function, as shown in the following formula.

[0154]

[0155] In the formula, k1 and k2 represent the indexes of the sound speed profile layers where the sound source point and the seabed reference station are located, respectively; c k represents the incident sound speed value in the k-th layer; g k represents the sound speed gradient value in the k-th layer;

[0156] Similarly, based on the robust Bayesian least square estimation principle, the time domain variation of the sound speed structure can be realized. It should be pointed out that when multiple time windows are used for calculation, the sound speed structure inversion of multiple windows can be completed through a loop.

[0157] In a preferred embodiment, the sound speed structure variation obtained by the above inversion is fed back to the calculation of the seabed reference station coordinates, and the seabed reference station position inversion is performed again until the coordinate increment and the sound speed variation are within the preset threshold range, and the gradual correction of the seabed reference station coordinates and the disturbed sound speed is completed.

[0158] III. Evidence of the effects of the embodiments. The embodiments of the present application have achieved some positive effects during research and development or use, and indeed have great advantages compared with the prior art. The following content is described in combination with data, graphs and other test processes.

[0159] The test adopts the measured data of the "deep sea comprehensive test of marine geodetic datum and new navigation technology in the South China Sea" in July 2019, the average water depth of the test area is about 3000 m, and the measuring ship is equipped with a GNSS receiver, a high-precision attitude sensor, a sound velocity profiler and a sea surface long baseline underwater acoustic positioning system and the like. The test adopts a multi-fold line sailing strategy, continuously observes the seabed datum station, collects sonar time measurement data, sound velocity profile data, GNSS data, attitude data and other sensor data, carries out experiments according to the positioning method provided by the application, uses a statistical traditional method (without applying sound velocity correction), a quadratic polynomial sound velocity correction method (QP method) and a cubic B-spline function sound velocity correction method (CBS method) provided by the application to perform calculation, and then compares the positioning accuracy of the seabed datum station of the traditional method and the application.

[0160] As shown in Figure 7 , the test shows that whether the sound velocity correction is applied, the vertical direction coordinate of the seabed datum station changes obviously; and only the accurate and appropriate sound velocity correction can not cause the deterioration of the plane direction positioning accuracy. Because the sound velocity correction is unreasonable, the component or projection of the radial direction time measurement error in the plane direction cannot be balanced or offset. The method of not correcting and only inverting the seabed datum station coordinates has an internal coincidence accuracy of the seabed datum station three-dimensional coordinates better than 10 cm, the unit weight mean error corresponding to the sonar time measurement observation value is 1.42 ms, the sound velocity correction method based on the QP method and the CBS method inverts the sound velocity structure change, and the corresponding unit weight mean error is 0.38 ms and 0.18 ms respectively, and the internal coincidence accuracy of the seabed station three-dimensional coordinates is better than 5 cm.

[0161] As shown in Figure 8 , the statistical characteristics of the sonar time measurement observation value residual are calculated, including the maximum value (Max), the mean value (Mean), the standard deviation (Std) and the root mean square error (RMS). From the mean value of the residual absolute value, the traditional method is 1.22 ms, the QP method is 0.36 ms, and the CBS method is 0.17 ms; from the root mean square error of the residual, the traditional method, the QP method and the CBS method are 1.43 ms, 0.44 ms and 0.21 ms respectively. It is not difficult to find that the root mean square error value of the residual under the inversion model with the sound velocity correction is obviously smaller; compared with the QP sound velocity correction method, the spline function used in the CBS sound velocity correction method has better smoothness and continuity, and is more in line with the actual characteristics of the smooth and continuous sound velocity structure change, so that better positioning effect is obtained.

[0162] It should be noted that embodiments of the present application can be realized by hardware, software, or a combination of software and hardware. The hardware portion can be realized by a special logic; the software portion can be stored in a memory and executed by a proper instruction execution system, such as a microprocessor or a specially designed hardware. Those skilled in the art can understand that the above-mentioned apparatus and method can be realized by executable instructions and / or included in processor control codes, such as a carrier medium, such as a magnetic disk, CD or DVD-ROM, a programmable memory, such as a read-only memory (firmware), or a data carrier, such as an optical or electronic signal carrier. The apparatus of the present application and its modules can be realized by a hardware circuit, such as a very large scale integrated circuit or a gate array, a semiconductor, such as a logic chip, a transistor, or a programmable hardware device, such as a field programmable gate array, a programmable logic device, or the like, by software executed by various types of processors, or by a combination of the above-mentioned hardware circuit and software, such as firmware.

[0163] The above description is merely a specific implementation of the present application, but the protection scope of the present application is not limited thereto. Any modification, equivalent replacement, and improvement within the technical range disclosed by the present application, and within the spirit and principle of the present application, should be included in the protection scope of the present application.

Claims

1. A high-precision positioning method for a seabed reference station considering time-domain changes in sound velocity structure, characterized in that, The high-precision positioning method of the seabed reference station considering the time-domain change of the sound velocity structure comprises the following steps: Based on the time-domain change of the sound velocity field and the vertical direction of the change of the sound velocity structure, a changeable cubic B-spline function is introduced to represent the time-domain change of the disturbed sound velocity, and based on the ray tracing, an iterative inversion model of the position information and the sound velocity information is constructed by step-by-step iteration and gradual correction, and the gradual correction of the coordinates of the seabed reference station and the disturbed sound velocity is carried out to obtain the high-precision position information of the seabed reference station. The high-precision positioning method of the seabed reference station considering the time-domain change of the sound velocity structure further comprises the following steps: Step one, the measurement ship continuously observes the seabed reference station to obtain the measurement ship track, the seabed reference station position and the sound velocity profile measurement rough position; Step two, the collected sound velocity profile is preprocessed, and the disturbed sound velocity structure is inverted by using the cubic B-spline function method; Step three, the iterative inversion of the seabed reference station position and the disturbed sound velocity structure is carried out to realize the gradual correction of the seabed reference station coordinates and the sound velocity structure; The step one of continuously observing the seabed reference station by the measurement ship to obtain the measurement ship track, the seabed reference station position and the sound velocity profile measurement rough position comprises the following steps: (1) The measurement ship sails on the sea surface to continuously observe the seabed reference station, collects the sonar measurement time data, the sound velocity profile data and the GNSS data, and obtains the measurement ship track and the sound velocity profile measurement data; (2) Layered equal-gradient ray tracing is carried out to inversely solve the rough coordinates of the seabed point; The step two of pre-processing the collected sound velocity profile and inversely calculating the disturbed sound velocity structure by using the cubic B-spline function method comprises the following steps: 1) The collected sound velocity profile is preprocessed, including sound velocity measurement outlier elimination, partial sound velocity profile data extension, adjacent same-gradient sound velocity layer merging, and calculation of the average sound velocity profile and the difference between each sound velocity profile and the average sound velocity profile; 2) The cubic B-spline function is introduced to fit the sound velocity structure change; The step three of iteratively inverting the seabed reference station position and the disturbed sound velocity structure to realize the gradual correction of the seabed reference station coordinates and the sound velocity structure comprises the following steps: (i) The fixed sound velocity-position solving method is used to invert the seabed station position, the nonlinear system of ray tracing is linearized at the rough coordinates of the seabed reference station according to the Taylor series linearization principle, and the seabed reference station position correction is calculated based on the robust Bayesian least squares estimation principle; (ii) The fixed position-sound velocity solving method is used to invert the sound velocity structure, the actual sound velocity structure is defined as the coupling of the reference sound velocity structure and the disturbed sound velocity structure, the cubic B-spline function is introduced to fit the sound velocity structure change, the nonlinear system of ray tracing is linearized based on the reference sound velocity profile, and the time-domain change of the sound velocity structure is inversely calculated based on the robust Bayesian least squares estimation principle; (iii) The seabed reference station coordinates and the disturbed sound velocity are gradually corrected; The step (i) of using the fixed sound velocity-position solving method to invert the seabed station position comprises: (i-1) The nonlinear system of ray tracing is linearized at the rough coordinates of the seabed reference station according to the Taylor series linearization principle, specifically as follows: t r = f 0 +A1·dX tp where f 0 is the initial coordinate of the seafloor reference station the propagation time obtained by ray tracing; dX tp is the position correction of the seafloor reference station, A1is the coefficient matrix of the propagation time with respect to the position of the seafloor reference station, where the derivation with respect to the plane position is obtained based on the basic principle of ray acoustics according to the following formula, and the derivation with respect to the z direction is obtained by a numerical analysis method; where p is the Snell constant, i.e., the sound speed level, determined from the ray tracing process; r is the horizontal distance between the source point and the seafloor reference station, (i-2) calculating the position correction of the seafloor reference station based on the robust Bayesian least square estimation principle, specifically: In the formula, P is the weight matrix of sonar observation value; P x is the prior weight matrix of model parameters; if the prior information of model parameters is not considered, that is, P x =0, the method for calculating the position correction of the seabed reference station based on the robust Bayesian least square estimation principle degenerates into the classical least square method; L=t r -f 0 is the observation vector; is the prior value of model parameters; it should be pointed out that if multiple stations are used for synchronous solution, the coordinates of the seabed reference station of multiple stations are solved through a loop. The fixed position-solved sound velocity in step (ii) includes: (ii-1) the actual sound velocity structure is the coupling of the reference sound velocity structure and the perturbed sound velocity structure; considering the time-domain variation of the sound velocity field, the sound velocity structure is expressed as follows: V(z, t e ) = V r (z) + V p (z, t e ) wherein, denotes the sequence of observation time, V denotes the sound speed structure ignoring horizontal spatial heterogeneity, which is related to depth and observation time; V r denotes the reference sound speed structure, which takes a certain actual measured sound speed profile or an average profile of multiple measured sound speed profiles; V p denotes the perturbed sound speed structure; (ii-2) Introducing cubic B-spline function for sound speed structure variation fitting, while considering the difference of sound speed structure variation between upper and lower layers of seawater, i.e. considering the characteristics of smooth and weak variation of sound speed in deep water layer within the time window of sonar observation, the perturbed sound speed is set as a piecewise function according to the empirical depth; the entire sequence of observation time is divided into m sections according to a certain time interval Δt; if the observation time belongs to the i-th section, the perturbed sound speed structure is as follows: ​ where a represents the B-spline function coefficient; B represents the B-spline base function; d represents the B-spline order, d=3 in this paper; z0 represents the reference depth, which is generally empirically valued at 1000-1600 m in deep water; u i is the unitized observation time, and the calculation formula is Δ=t i+1 -t i ; t i , t i+1 are the start and end observation times of the i-th section, respectively; (ii-3) linearizing the ray tracing nonlinear system based on the reference sound velocity profile to obtain the following formula: t r = f 0 +A2·da In the formula, The derivative formula of the sound speed perturbation with respect to the spline coefficient is derived from the derivative principle of the cubic B-spline function. The derivative numerical calculation formula of the sound signal propagation time with respect to the perturbed sound speed is referenced from the derivative of the discrete function, as shown in the following formula: In the formula, k1 and k2 respectively represent the index of the sound velocity profile layer where the sound source point and the seabed reference station are located; c k represents the incident sound velocity value in the kth layer; g k represents the sound velocity gradient value in the kth layer; (ii-4) inverting the time-domain variation of the sound velocity structure based on the robust Bayesian least square estimation principle; it should be noted that when there are multiple time windows for calculation, the sound velocity structure inversion of multiple windows is completed through a loop; The gradual correction of the seafloor reference station coordinates and the perturbed sound velocity in step (iii) includes: (iii-1) the sound velocity structure variation obtained by the sound velocity structure inversion of the fixed position-solved sound velocity is applied to the calculation of the seafloor reference station coordinates, and the seafloor reference station position inversion is performed again; (iii-2) when the coordinate increment and the sound velocity variation are within the preset threshold range, the calculation is ended, and the gradual correction of the seafloor reference station coordinates and the perturbed sound velocity is completed.

2. A high-precision positioning system of a seafloor reference station considering the time-domain variation of sound velocity structure, which applies the high-precision positioning method of a seafloor reference station considering the time-domain variation of sound velocity structure according to claim 1, characterized in that, The high-precision positioning system of the seafloor reference station considering the time-domain variation of the sound velocity structure includes: A measurement approximate position acquisition module is configured to measure continuous observation of a ship on a seafloor reference station, and to acquire a measurement ship track, a seafloor reference station position, and a sound velocity profile measurement approximate position. A perturbed sound velocity structure inversion module is configured to perform data preprocessing on the collected sound velocity profile, and to perform perturbed sound velocity structure inversion using a cubic B-spline function method. A reference station coordinate and sound velocity structure gradual correction module is configured to perform iterative inversion on the seafloor reference station position and the perturbed sound velocity structure, and to achieve gradual correction of the seafloor reference station coordinates and the sound velocity structure.

3. An apparatus, comprising: The device includes a memory and a processor, and the memory stores a program which, when executed by the processor, causes the processor to perform the high-precision positioning method of the seafloor reference station considering the time-domain variation of the sound velocity structure.

Images