Underwater single-beacon Gaussian Newton positioning method and system under linear acoustic profile condition

Abstract

The invention belongs to the technical field of underwater positioning, and discloses an underwater single-beacon Gaussian Newton positioning method and system under a linear sound profile condition. According to the method, underwater acoustic signal propagation time observation noise is modeled into Gaussian distribution, a maximum likelihood estimation problem is constructed, and then the maximum likelihood estimation problem is converted into a nonlinear least square estimation problem; and solving the gradient of the observed quantity to the position coordinates of the AUV based on a chain rule, iteratively solving a nonlinear least square estimation problem based on a Gauss-Newton method, and obtaining the position coordinates of the AUV after iterative convergence. Compared with other single beacon positioning methods which do not consider sound propagation refraction, the method of the invention can obtain a more accurate position calculation result under a linear sound profile condition, and enhances single beacon positioning, especially the practical application capability of the method in deep sea areas.

Description

Technical Field

[0001] The invention belongs to the field of underwater positioning technology, and in particular relates to an underwater single-beacon Gauss-Newton positioning method and system under linear acoustic profile conditions. Background Art

[0002] Underwater vehicles (AUVs) are widely used in underwater search and rescue, underwater operations, and marine environment data collection. AUVs need to obtain their own position coordinates in real time to ensure the accuracy of motion control and the validity of the collected data. Since electromagnetic wave signals attenuate rapidly underwater, hydroacoustic positioning has become the most widely used and effective positioning method in underwater environments. Compared with the common underwater long baseline positioning and ultra-short baseline positioning technologies, the positioning method that combines a single hydroacoustic beacon with dead reckoning information only requires a single hydroacoustic beacon with a known position to provide positioning assistance, which can greatly reduce the time and economic cost of long baseline positioning. The beacon can be fixed to the seabed or carried on a surface ship, and the system has high flexibility. In view of the above advantages of single beacon positioning, it has received widespread attention from researchers.

[0003] The current underwater single beacon positioning method usually assumes that the water acoustic speed is constant, and directly uses the distance between the AUV and the beacon as the observed value. In actual underwater positioning applications, affected by factors such as seawater temperature, salinity and depth, the water acoustic speed varies with space, especially in the depth direction, the sound speed value changes with a large gradient, which will cause obvious sound propagation refraction. Ignoring the influence of sound propagation refraction will cause a large distance calculation error, which will in turn lead to a large single beacon positioning error, affecting the actual application capability of the single beacon system. Studies have shown that for deep-sea environments, the water acoustic speed profile can be approximated as linear, that is, the sound speed value can be approximated as a linear function of depth. In other words, studying the single beacon positioning method under the condition of linear sound speed profile has broad application prospects. By more accurately describing the trajectory of the sound propagation channel under the condition of linear sound profile and the signal propagation time observation model, the error of single beacon positioning can be reduced to a certain extent. Summary of the invention

[0004] In order to overcome the problems existing in the related art, the embodiments disclosed in the present invention provide an underwater single-beacon Gauss-Newton positioning method and system under linear acoustic profile conditions.

[0005] The technical solution is as follows: A Gauss-Newton positioning method for a single underwater beacon under linear acoustic profile conditions, comprising:

[0006] S1, take any point in the positioning area as the origin, set the east, north and sky directions as X, Y and Z axes respectively, and establish the underwater single beacon positioning geodetic coordinate system; take the AUV centroid as the origin, and set the front, right and bottom directions as x, y and z axes respectively, and establish the AUV body coordinate system;

[0007] S2, constructing a dead reckoning model based on the velocity of the AUV in the body coordinate system detected by the DVL and the attitude of the AUV detected by the attitude sensor;

[0008] S3, the hydroacoustic beacon periodically emits hydroacoustic signals, the beacon position coordinates are known, the AUV is equipped with a hydrophone, the propagation time of the hydroacoustic signal is measured, the sound velocity profile is modeled as a linear sound profile, and the known hydroacoustic profile parameter information and the hydroacoustic beacon position coordinates are combined to construct a hydroacoustic observation model under the linear sound profile condition, and the relationship between the signal propagation time and the AUV position is obtained;

[0009] S4, integrating AUV dead reckoning information and using all acoustic signal propagation time observations within the time window width N, constructs a virtual long baseline positioning model under linear acoustic profile conditions;

[0010] S5, calculates the gradient matrix of the propagation time of the underwater acoustic signal relative to the position variable, uses the Gauss-Newton method to iteratively solve the positioning problem, and obtains the position estimate of the AUV after the iteration converges.

[0011] In step S2, a dead reckoning model is constructed, including:

[0012] The relationship between the AUV position variables at time k and at time ki,i≥1 is:

[0013]

[0014] In the formula, x k ,y k ,z k is the real position coordinate of AUV at time k, x k-i is the x-direction position coordinate of the AUV at time ki, y k-i is the y-direction position coordinate of the AUV at time ki, z k-i is the z-direction position coordinate of the AUV at time ki, Δt is the discrete time interval, i is the time index variable, and v x,k-t is the x-direction velocity of AUV at time kt, v y,k-t is the y-direction velocity of AUV at time kt, v z,k-t is the z-direction velocity of the AUV at time kt.

[0015] In step S3, a hydroacoustic observation model under the condition of linear acoustic profile is constructed, and the expression is:

[0016]

[0017] Where P b is the coordinate vector of the hydroacoustic beacon position, x bis the east coordinate of the hydroacoustic beacon, y b is the north coordinate of the hydroacoustic beacon, z b The celestial coordinates of the hydroacoustic beacon;

[0018] The sound velocity profile is modeled as linear, and the sound velocity value is a linear function of the depth:

[0019] c=az+b

[0020] Where c is the depth at depth z, a is the sound velocity gradient, z is the depth, and b is the sound velocity at sea level;

[0021] When the sound velocity profile is known, the sound velocity at the beacon is known and is denoted by c b , the speed of sound at the AUV position at time k is recorded as c k , the propagation time of the underwater acoustic signal at time k is recorded as t k , the relationship between horizontal distance and signal propagation time is:

[0022]

[0023] Among them, r k is the horizontal distance between the AUV and the beacon, and e is the natural base;

[0024] Squaring both sides of the above equation yields:

[0025]

[0026] Simplifying gives:

[0027]

[0028] Substituting the sound velocity profile parameters into the equation:

[0029]

[0030] definition:

[0031]

[0032] In the formula, m k is a virtual observation, in which all parameters are known and the value of each parameter depends on the signal propagation time observation value and the sound profile parameter at the current moment;

[0033] The relationship between the AUV position and the propagation time of the underwater acoustic signal is obtained as follows:

[0034] 0=(x k -x b ) 2 +(y k -y b ) 2+(z k -z b ) 2 -m k (az k +b).

[0035] In step S4, a virtual long baseline positioning model under the linear acoustic profile condition is constructed, including:

[0036] Define the time window width as N, use the propagation time of the underwater acoustic signal received from kN to k, solve the position coordinates of the AUV at k, and for the ki moment, apply the established observation equation of the underwater acoustic signal propagation time as follows:

[0037] 0=(x k-i -x b ) 2 +(y k-i -y b ) 2 +(z k-i -z b ) 2 -m k-i (az k-i +b),

[0038] for i=1,2…N

[0039] In the formula, m k-i is the virtual observation at time ki.

[0040] Substitute the AUV’s dead reckoning model, that is, the position relationship between time k and time ki, and we get:

[0041] 0=(x k -x b,i ) 2 +(y k -y b,i ) 2 +(z k -z b,i ) 2 -m k-i (az k +b i ), for i=1,2…N

[0042] in,

[0043]

[0044] In the formula, x b,i is the x-axis position coordinate of the ith virtual beacon, y b,i is the y-axis position coordinate of the i-th virtual beacon, z b,i is the z-axis position coordinate of the ith virtual beacon, bi is the surface sound velocity corresponding to the i-th virtual beacon; further, define:

[0045] x b,0 =x b

[0046] y b,0 =y b

[0047] z b,0 =z b

[0048] b0=b

[0049] In the formula, x b,0 ,y b,0 ,z b,0 is the virtual beacon coordinate at time k-0, b0 is the virtual surface sound speed at time k-0;

[0050] Combined with the observation equation at time k, all N+1 observation equations are expressed as follows:

[0051] 0=(x k -x b,i ) 2 +(y k -y b,i ) 2 +(z k -z b,i ) 2 -m k-i (az k +b i ),

[0052] for i=0,1,2…N.

[0053] In step S5, the positioning problem is iteratively solved using the Gauss-Newton method, including:

[0054] The propagation time of the underwater acoustic signal contains observation noise, and the observation value of the propagation time of the underwater acoustic signal is recorded as:

[0055]

[0056] In the formula, is the observed value of the propagation time of the underwater acoustic signal, t k is the actual underwater acoustic signal propagation time, v k is the observation noise of the underwater acoustic signal propagation time, and its mean is modeled as 0 and the variance is Gaussian random variable; all N+1 underwater acoustic signal propagation time observations are combined into a vector All N+1 real underwater acoustic signal propagation times are combined into a vector tk , all N+1 underwater acoustic propagation time observation noises are integrated into a vector v k ,as follows:

[0057]

[0058] The observation noise v k Modeled as Gaussian noise with a mean of 0 and a variance matrix denoted as Q k , then the maximum likelihood estimation problem for single beacon positioning is:

[0059]

[0060] In the formula, is the maximum likelihood estimate of the AUV position at time k, arg max is the maximization problem, exp{} is the exponential function, and p k is the position of the AUV at time k, the superscript T represents the vector matrix transpose, t k (p k ) is the real propagation time of the underwater acoustic signal, and its value is a function of the position coordinates; the variance matrix Q k It is obtained according to the observed variance of the propagation time of the underwater acoustic signal at each moment, as follows:

[0061]

[0062] In the formula, σ k-N is the standard deviation of the observed variance at time kN, σ k-N+1 is the standard deviation of the observed variance at time k-N+1, σ k is the standard deviation of the observed variance at time k, diag(a) is a diagonal matrix with the elements of a as diagonal elements. According to the monotonicity of the logarithmic function, the above optimization problem is transformed into:

[0063]

[0064] Where arg min( ) is the minimization function;

[0065] definition:

[0066]

[0067] Among them, S k for Choleskey decomposition, the optimization problem is transformed into:

[0068]

[0069] This problem is a nonlinear least squares optimization problem, and the Gauss-Newton method is used to solve it as follows:

[0070]

[0071] in, is the position estimate obtained after j+1 iterations, is the position estimate obtained at the jth iteration, is the Gauss-Newton position increment for the jth iteration, defined as:

[0072]

[0073] In the formula, is the gradient matrix of the signal propagation time to the position coordinates obtained at the jth iteration, is the acoustic signal propagation time calculated based on the AUV position estimated in the current iteration cycle, Defined as:

[0074]

[0075] In the formula, is the partial derivative of the propagation time of the underwater acoustic signal with respect to the position coordinate in the x direction, is the partial derivative of the underwater acoustic signal propagation time with respect to the position coordinate in the y direction, is the partial derivative of the underwater acoustic signal propagation time with respect to the position coordinate in the z direction, is the x, y, and z position estimate obtained at the jth iteration.

[0076] Furthermore, the gradient matrix of the position coordinates Each element is calculated as follows:

[0077] According to the hydroacoustic beacon observation equation:

[0078] 0=(x k -x b,i ) 2 +(y k -y b,i ) 2 +(z k -z b,i ) 2 -m k-i (az k +b i ),

[0079] for i=0,1,2…N

[0080] Two sides x k ,y k ,z k Find the partial derivatives respectively, and we get:

[0081]

[0082] In the formula, is the partial derivative of the virtual observation with respect to the propagation time of the underwater acoustic signal, are the partial derivatives of the underwater acoustic signal propagation time with respect to the AUV position;

[0083] According to m k-i The definition of , we get:

[0084]

[0085] Then we have:

[0086]

[0087] get:

[0088]

[0089] for i=0,1…N.

[0090] Furthermore, the acoustic signal propagation time is calculated based on the AUV position estimated in the current iteration cycle. The method is:

[0091] According to the equation:

[0092]

[0093] get:

[0094] n i α i =(α i -1) 2

[0095] Where n i With α i are auxiliary parameters, and their definitions are:

[0096]

[0097] Solving the above quadratic equation yields:

[0098]

[0099] Right now:

[0100]

[0101] According to the above equation, we can get two different t k-i , since the real underwater acoustic signal propagation time is a positive number, select the one that is greater than 0; substitute the position coordinates obtained in the jth iteration into the above equation to obtain the underwater acoustic signal propagation time vector

[0102] Furthermore, the iterative convergence criterion is set as follows: When the second norm of is less than a certain threshold ε, the iteration is considered to converge and the single beacon positioning position estimate is output

[0103] Another object of the present invention is to provide an underwater single-beacon Gauss-Newton positioning system under linear sound profile conditions, which implements the underwater single-beacon Gauss-Newton positioning method under linear sound profile conditions. The system addresses the problem of refraction of the sound propagation channel caused by the change of water sound velocity with depth under linear sound profile conditions, and obtains the AUV position coordinates through a positioning solution method based on the Gauss-Newton algorithm.

[0104] Furthermore, the system is provided with an AUV, wherein the AUV is equipped with a hydrophone, a DVL, a posture sensor and a CTD; the AUV detects the sound velocity profile parameters of the experimental area through the CTD; the hydroacoustic beacon periodically broadcasts a hydroacoustic signal, and the AUV periodically observes its velocity in the body coordinate system through the DVL carried by the AUV, calculates the velocity of the AUV in the geodetic coordinate system in combination with the AUV attitude angle observed by the posture sensor, performs dead reckoning, and records and stores the dead reckoning data; after receiving the hydroacoustic signal emitted by the hydroacoustic beacon, the AUV obtains the hydroacoustic signal transmission time in combination with the known hydroacoustic signal emission time. propagation time; when the AUV receives a certain number of hydroacoustic signals, a virtual long baseline array is constructed in combination with the stored dead reckoning data; the virtual long baseline positioning problem is constructed in combination with the relationship between the hydroacoustic propagation time and the AUV position under the linear acoustic profile; the hydroacoustic signal propagation time observation noise is modeled as a Gaussian distribution, a maximum likelihood estimation problem is constructed, and then it is converted into a nonlinear least squares estimation problem; the gradient of the observed quantity with respect to the AUV position coordinates is solved based on the chain rule, and the nonlinear least squares estimation problem is iteratively solved based on the Gauss-Newton method, and the AUV position coordinates are obtained after the iteration converges.

[0105] In combination with all the above-mentioned technical solutions, the beneficial effects of the present invention are as follows: the present invention focuses on the problem of sound propagation refraction caused by the spatial variation of sound speed in underwater single beacon positioning, and designs a single beacon positioning position solution method under linear sound profile conditions based on the Gauss-Newton method. The underwater vehicle (AUV) is equipped with a hydrophone, a Doppler velocimeter, an attitude sensor and a temperature-salinity-depth meter; the AUV detects the sound velocity profile parameters of the experimental area through the temperature-salinity-depth meter; the hydroacoustic beacon broadcasts hydroacoustic signals periodically, and the AUV uses the Doppler velocimeter carried by the AUV in combination with the attitude sensor to perform dead reckoning and record and store the dead reckoning data; after receiving the hydroacoustic signal emitted by the hydroacoustic beacon, the AUV calculates the propagation time of the hydroacoustic signal; when the AUV receives a certain number of hydroacoustic signals, it constructs a virtual long baseline array in combination with the stored dead reckoning data; combining the relationship between the hydroacoustic propagation time and the AUV position under the linear sound profile, a virtual long baseline positioning problem is constructed; the hydroacoustic signal propagation time observation noise is modeled as a Gaussian distribution, a maximum likelihood estimation problem is constructed, and then it is converted into a nonlinear least squares estimation problem; the gradient of the observed quantity with respect to the AUV position coordinates is solved based on the chain rule, and the nonlinear least squares estimation problem is iteratively solved based on the Gauss-Newton method, and the AUV position coordinates are obtained after the iteration converges. Compared with other single-beacon positioning methods that do not consider the refraction of sound propagation, the method of the present invention can obtain more accurate position solution results under linear sound profile conditions, enhance single-beacon positioning, especially its practical application capability in deep sea areas.

[0106] The present invention takes into account the refraction of the underwater acoustic propagation channel that is widely present in real underwater acoustic positioning applications, and adopts a more reasonable and accurate underwater acoustic observation model, which can obtain better ranging accuracy and positioning accuracy, improve the effectiveness and value of the data collected by the AUV, and enhance the operation and application value of the AUV itself. The present invention takes into account the refraction problem of the underwater acoustic propagation channel that is widely present in single beacon positioning.

[0107] Affected by the temperature, salinity and depth changes of seawater, the gradient of the water acoustic velocity in the depth direction is large, which will cause obvious refraction of the sound propagation channel. This problem is widely present in all underwater positioning methods based on water acoustic ranging, including single beacon positioning. The current common processing method is to ignore the refraction of the channel and directly use a constant sound velocity value, or introduce the concept of effective sound velocity and directly regard the effective sound velocity as unknown. Ignoring the refraction of the channel will cause ranging errors and positioning errors, while directly regarding the effective sound velocity as unknown will ignore the valuable physical information of the water acoustic velocity profile, which will also cause performance loss and positioning errors. It is urgent to study the single beacon positioning method considering the refraction of sound propagation. The present invention directly and accurately describes the water acoustic observation model under the linear sound profile, directly processes the observation model considering the propagation refraction, and uses the Gauss-Newton method to iteratively solve the AUV position, which can obtain better position estimation performance and solve the technical problems that people are eager to solve. BRIEF DESCRIPTION OF THE DRAWINGS

[0108] The accompanying drawings herein are incorporated in and constitute a part of the specification, illustrate embodiments consistent with the present disclosure, and together with the description, serve to explain the principles of the present disclosure;

[0109] Figure 1 It is a flow chart of an underwater single beacon Gauss-Newton positioning method under linear acoustic profile conditions provided by an embodiment of the present invention;

[0110] Figure 2 It is a schematic diagram of the underwater single beacon positioning geodetic coordinate system and the AUV body coordinate system provided by an embodiment of the present invention;

[0111] Figure 3 It is the simulation trajectory diagram;

[0112] Figure 4 The sound velocity profile used for the simulation;

[0113] Figure 5 This is a comparison chart of positioning errors of different positioning methods;

[0114] Figure 6 This is a graph of the number of iterations at each moment of the Gauss-Newton method. DETAILED DESCRIPTION

[0115] In order to make the above-mentioned objects, features and advantages of the present invention more obvious and easy to understand, the specific embodiments of the present invention are described in detail below in conjunction with the accompanying drawings. In the following description, many specific details are set forth to facilitate a full understanding of the present invention. However, the present invention can be implemented in many other ways different from those described herein, and those skilled in the art can make similar improvements without violating the connotation of the present invention, so the present invention is not limited by the specific implementation disclosed below.

[0116] The innovation of the present invention is that the present invention innovatively introduces the relationship between horizontal distance and signal propagation time under the condition of linear sound velocity profile, and innovatively derives the hydroacoustic observation model under the condition of linear sound profile. The present invention integrates the dead reckoning information of AUV, and constructs a virtual long baseline positioning method under the condition of linear sound profile under the framework of virtual long baseline positioning. The present invention designs an AUV position solution method based on the Gauss-Newton method: first initialize the iteration, calculate the gradient of the hydroacoustic signal propagation time of the current iteration step for the AUV position, construct the gradient matrix, and then calculate the Gauss-Newton descent direction, correct the AUV position coordinates of the current iteration step, and obtain the position coordinates of the AUV after the iteration converges. The method of the present invention takes into account the refraction of the hydroacoustic propagation channel in single beacon positioning, and can obtain more accurate position estimation results and enhance the practical application ability of AUV. The present invention designs a positioning solution method based on the Gauss-Newton algorithm for the problem of refraction of the sound propagation channel caused by the change of hydroacoustic speed with depth under the condition of linear sound profile, so as to obtain better positioning performance than the method that does not consider the influence of sound propagation refraction.

[0117] Embodiment 1, as Figure 1 As shown, the underwater single beacon Gauss-Newton positioning method under linear acoustic profile conditions provided by an embodiment of the present invention includes:

[0118] S1, take any point in the positioning area as the origin, set the east, north and sky directions as X, Y and Z axes respectively, and establish the underwater single beacon positioning geodetic coordinate system; take the AUV centroid as the origin, and set the front, right and bottom directions as x, y and z axes respectively, and establish the AUV body coordinate system;

[0119] S2, constructing a dead reckoning model based on the velocity of the AUV in the body coordinate system detected by the DVL and the attitude of the AUV detected by the attitude sensor;

[0120] S3, the hydroacoustic beacon periodically emits hydroacoustic signals, the beacon position coordinates are known, the AUV is equipped with a hydrophone, the propagation time of the hydroacoustic signal is measured, the sound velocity profile is modeled as a linear sound profile, and the known hydroacoustic profile parameter information and the hydroacoustic beacon position coordinates are combined to construct a hydroacoustic observation model under the linear sound profile condition, and the relationship between the signal propagation time and the AUV position is obtained;

[0121] S4, integrating AUV dead reckoning information and using all acoustic signal propagation time observations within the time window width N, constructs a virtual long baseline positioning model under linear acoustic profile conditions;

[0122] S5, calculates the gradient matrix of the propagation time of the underwater acoustic signal relative to the position variable, uses the Gauss-Newton method to iteratively solve the positioning problem, and obtains the position estimate of the AUV after the iteration converges.

[0123] Exemplarily, in step S1, the established geodetic coordinate system and body coordinate system are as follows: Figure 2 shown.

[0124] In step S2, a dead reckoning model is constructed, and the horizontal coordinates of the AUV at time k are represented as:

[0125]

[0126] Among them, x k represents the x-coordinate at time k, y k represents the y-direction coordinate at time k, z k represents the z-direction coordinate at time k, and the velocity of the AUV detected by DVL at time k in the body coordinate system is recorded as:

[0127]

[0128] The AUV attitude angle detected by the attitude sensor is recorded as:

[0129]

[0130] In the formula, are the forward, rightward and downward speeds of the AUV respectively, θ k ,ψ k They are the roll angle, pitch angle and heading angle of the AUV respectively.

[0131] Combining the velocity of the AUV in the body coordinate system and the AUV attitude angle, the velocity of the AUV in the earth coordinate system at time k can be obtained as:

[0132]

[0133] In the formula, R k is the rotation matrix, defined as:

[0134]

[0135] According to the position and speed of AUV at time k-1, the position coordinates of AUV at time k are obtained as follows:

[0136] P k =P k-1 +v k-1 Δt

[0137] Where Δt represents the discrete time interval of dead reckoning;

[0138] Expressing the above in component form is:

[0139] x k =x k-1 +v x,k-1Δt

[0140] y k =y k-1 +v y,k-1 Δt

[0141] z k =z k-1 +v z,k-1 Δt

[0142] Then the relationship between the AUV position variables at time k and time ki (i ≥ 1) is obtained as follows:

[0143]

[0144]

[0145] In the formula, x k-i is the x-direction position coordinate of the AUV at time ki, y k-i is the y-direction position coordinate of the AUV at time ki, z k-i is the z-direction position coordinate of the AUV at time ki, Δt is the discrete time interval, i is the time index variable, and v x,k-t is the x-direction velocity of AUV at time kt, v y,k-t is the y-direction velocity of AUV at time kt, v z,k-t is the z-direction velocity of the AUV at time kt.

[0146] In step S3, a hydroacoustic observation model under the condition of linear acoustic profile is constructed, and the expression is:

[0147]

[0148] Where P b is the coordinate vector of the hydroacoustic beacon position, x b is the east coordinate of the hydroacoustic beacon, y b is the north coordinate of the hydroacoustic beacon, z b The celestial coordinates of the hydroacoustic beacon;

[0149] The sound velocity profile is modeled as linear, and the sound velocity value is a linear function of the depth:

[0150] c=az+b

[0151] Where c is the depth at depth z, a is the sound velocity gradient, z is the depth, and b is the sound velocity at sea level;

[0152] When the sound velocity profile is known, the sound velocity at the beacon is known and is denoted by c b , the speed of sound at the AUV position at time k is recorded as c k , the propagation time of the underwater acoustic signal at time k is recorded as tk , the relationship between horizontal distance and signal propagation time is:

[0153]

[0154] Among them, r k is the horizontal distance between the AUV and the beacon, and e is the natural base;

[0155] The present invention innovatively squares both sides of the above equation to obtain:

[0156]

[0157] Simplifying gives:

[0158]

[0159] Substituting the sound velocity profile parameters into the equation:

[0160]

[0161] definition:

[0162]

[0163] In the formula, m k is a virtual observation, in which all parameters are known and the value of each parameter depends on the signal propagation time observation value and the sound profile parameter at the current moment;

[0164] The relationship between the AUV position and the propagation time of the underwater acoustic signal is obtained as follows:

[0165] 0=(x k -x b ) 2 +(y k -y b ) 2 +(z k -z b ) 2 -m k (az k +b).

[0166] In step S4, a virtual long baseline positioning model under the linear acoustic profile condition is constructed, including:

[0167] Define the time window width as N, use the propagation time of the underwater acoustic signal received from kN to k, solve the position coordinates of the AUV at k, and for the ki moment, apply the established observation equation of the underwater acoustic signal propagation time as follows:

[0168] 0=(x k-i -x b ) 2+(y k-i -y b ) 2 +(z k-i -z b ) 2 -m k-i (az k-i +b), for i=1,2…N

[0169] In the formula, m k-i is the virtual observation at time ki.

[0170] Substitute the AUV’s dead reckoning model, that is, the position relationship between time k and time ki, and we get:

[0171] 0=(x k -x b,i ) 2 +(y k -y b,i ) 2 +(z k -z b,i ) 2 -m k-i (az k +b i ), for i=1,2…N

[0172] in,

[0173]

[0174] In the formula, x b,i is the x-axis position coordinate of the ith virtual beacon, y b,i is the y-axis position coordinate of the i-th virtual beacon, z b,i is the z-axis position coordinate of the ith virtual beacon, b i is the surface sound velocity corresponding to the i-th virtual beacon; further, define:

[0175] x b,0 =x b

[0176] y b,0 =y b

[0177] z b,0 =z b

[0178] b0=b

[0179] In the formula, x b,0 ,y b,0 ,z b,0 is the virtual beacon coordinate at time k-0, b0 is the virtual surface sound speed at time k-0;

[0180] Combined with the observation equation at time k, all N+1 observation equations are expressed as follows:

[0181] 0=(x k -x b,i ) 2 +(y k -y b,i ) 2 +(z k -z b,i ) 2 -m k-i (az k +b i ),

[0182] for i=0,1,2…N.

[0183] In step S5, the positioning problem is iteratively solved using the Gauss-Newton method, including:

[0184] In real applications, the propagation time of the underwater acoustic signal contains observation noise. The observed value of the propagation time of the underwater acoustic signal is recorded as:

[0185]

[0186] In the formula, is the observed value of the propagation time of the underwater acoustic signal, t k is the actual underwater acoustic signal propagation time, v k is the observation noise of the underwater acoustic signal propagation time, and its mean is modeled as 0 and the variance is Gaussian random variable; all N+1 underwater acoustic signal propagation time observations are combined into a vector All N+1 real underwater acoustic signal propagation times are combined into a vector t k , all N+1 underwater acoustic propagation time observation noises are integrated into a vector v k ,as follows:

[0187]

[0188]

[0189] The observation noise v k Modeled as Gaussian noise with a mean of 0 and a variance matrix denoted as Q k , then the maximum likelihood estimation problem for single beacon positioning is:

[0190]

[0191] In the formula, is the maximum likelihood estimate of the AUV position at time k, arg max is the maximization problem, exp{} is the exponential function, p k is the position of the AUV at time k, the superscript T represents the vector matrix transpose, t k (p k ) is the real propagation time of the underwater acoustic signal, and its value is a function of the position coordinates; the variance matrix Q k It is obtained according to the observed variance of the propagation time of the underwater acoustic signal at each moment, as follows:

[0192]

[0193] In the formula, σ k-N is the standard deviation of the observed variance at time kN, σ k-N+1 is the standard deviation of the observed variance at time k-N+1, σ k is the standard deviation of the observed variance at time k, diag(a) is a diagonal matrix with the elements of a as diagonal elements. According to the monotonicity of the logarithmic function, the above optimization problem is transformed into:

[0194]

[0195] Where arg min( ) is the minimization function;

[0196] definition:

[0197]

[0198] Among them, S k for Choleskey decomposition, the optimization problem is transformed into:

[0199]

[0200] This problem is a nonlinear least squares optimization problem, and the Gauss-Newton method is used to solve it as follows:

[0201]

[0202] in, is the position estimate obtained after j+1 iterations, is the position estimate obtained at the jth iteration, is the Gauss-Newton position increment for the jth iteration, defined as:

[0203]

[0204] In the formula, is the gradient matrix of the signal propagation time to the position coordinates obtained at the jth iteration, is the acoustic signal propagation time calculated based on the AUV position estimated in the current iteration cycle, Defined as:

[0205]

[0206] In the formula, is the partial derivative of the propagation time of the underwater acoustic signal with respect to the position coordinate in the x direction, is the partial derivative of the underwater acoustic signal propagation time with respect to the position coordinate in the y direction, is the partial derivative of the underwater acoustic signal propagation time with respect to the position coordinate in the z direction, is the x, y, and z position estimate obtained at the jth iteration.

[0207] Prioritize the gradient matrix of the position coordinates Each element is calculated as follows:

[0208] According to the hydroacoustic beacon observation equation:

[0209] 0=(x k -x b,i ) 2 +(y k -y b,i ) 2 +(z k -z b,i ) 2 -m k-i (az k +b i ),

[0210] for i=0,1,2…N

[0211] Two sides x k ,y k ,z k Find the partial derivatives respectively, and we get:

[0212]

[0213] In the formula, is the partial derivative of the virtual observation with respect to the propagation time of the underwater acoustic signal, are the partial derivatives of the underwater acoustic signal propagation time with respect to the AUV position;

[0214] According to m k-i The definition of , we get:

[0215]

[0216] Then we have:

[0217]

[0218] get:

[0219]

[0220] for i=0,1…N.

[0221] Prioritize the acoustic signal propagation time calculated based on the AUV position estimated in the current iteration cycle The method is:

[0222] According to the equation:

[0223]

[0224] for i=0,1,2…N

[0225] get:

[0226] n i α i =(α i -1) 2

[0227] Where n i With α i are auxiliary parameters, and their definitions are:

[0228]

[0229] Solving the above quadratic equation yields:

[0230]

[0231] Right now:

[0232]

[0233] According to the above equation, we can get two different t k-i , since the real underwater acoustic signal propagation time is a positive number, select the one that is greater than 0; substitute the position coordinates obtained in the jth iteration into the above equation to obtain the underwater acoustic signal propagation time vector

[0234] Set the iterative convergence criterion to when the position update amount When the second norm of is less than a certain threshold ε, the iteration is considered to converge and the single beacon positioning position estimate is output

[0235] The pseudo code of the underwater single-beacon Gauss-Newton positioning method under the linear acoustic profile condition provided by the embodiment of the present invention is shown in Table 1.

[0236] Table 1 Pseudo code of Gauss-Newton positioning method for underwater single beacon under linear acoustic profile conditions

[0237]

[0238]

[0239] Example 2: The underwater single-beacon Gauss-Newton positioning system under linear acoustic profile conditions provided in the embodiment of the present invention focuses on the refraction problem of the acoustic propagation channel under linear sound velocity profile conditions, and designs a single-beacon positioning position solution method based on the Gauss-Newton method. The AUV of the system is equipped with a hydrophone, a DVL, a posture sensor, and a CTD; the AUV detects the sound velocity profile parameters of the experimental area through the CTD; the hydroacoustic beacon periodically broadcasts hydroacoustic signals, and the AUV periodically observes its velocity in the body coordinate system through the DVL it carries, and calculates the AUV velocity in the geodetic coordinate system in combination with the AUV attitude angle observed by the posture sensor, performs dead reckoning, and records and stores the dead reckoning data; after receiving the hydroacoustic signal emitted by the hydroacoustic beacon, the AUV obtains the propagation time of the hydroacoustic signal in combination with the known hydroacoustic signal emission time; When V receives a certain number of hydroacoustic signals, a virtual long baseline array is constructed in combination with the stored dead reckoning data; a virtual long baseline positioning problem is constructed in combination with the relationship between the hydroacoustic propagation time and the AUV position under the linear acoustic profile; the observation noise of the hydroacoustic signal propagation time is modeled as a Gaussian distribution, a maximum likelihood estimation problem is constructed, and then it is converted into a nonlinear least squares estimation problem; the gradient of the observed quantity with respect to the AUV position coordinates is solved based on the chain rule, and the nonlinear least squares estimation problem is iteratively solved based on the Gauss-Newton method, and the AUV position coordinates are obtained after the iteration converges.

[0240] To further illustrate the relevant effects of the embodiments of the present invention, the following simulation experiments are conducted: This embodiment simultaneously demonstrates the positioning performance of three single beacon positioning methods that do not consider the refraction of sound propagation. The sound velocity value is calculated using the average value of the entire linear sound velocity profile. Combined with the observed propagation time of the hydroacoustic signal, the distance observation value between the AUV and the beacon can be calculated, and the distance is used to estimate the position coordinates of the AUV. Method 1 is based on weighted least squares to solve the virtual long baseline positioning problem (denoted as R-WLS), method 2 is based on the semidefinite programming method in the literature Yanshen Du, Ping Wei, and Huaguo Zhang. Semidefinite programming approaches for source localization problems. In 2014 IEEE International Conference on Communiction Problem-solving, pages 339–342, 2014. to solve the virtual long baseline positioning problem (denoted as R-SDR1), and method 3 is based on the semidefinite programming method in the literature KW Cheung, WK Ma, and HCSo. Accurate approximation algorithm for to a-based maximum likelihood mobile location using semidefinite programming. In 2004 IEEE International Conference on Acoustics, Speech, and Signal Processing, volume 2, pages ii–145, 2004. to solve the virtual long baseline positioning problem (denoted as R-SDR2). The method of the present invention directly uses the underwater sound propagation time as the observed quantity, constructs a nonlinear least square estimation problem, and solves it based on the Gauss-Newton method, which is abbreviated as T-GN.

[0241] The simulation parameters are set as follows: AUV dead reckoning interval is 0.1 seconds, acoustic signal transmission interval is 15 seconds, and acoustic beacon is fixed at [0 0 100] T m, and the total simulation time is 500 seconds. The AUV dynamics model is used to generate the spiral diving motion trajectory of the AUV. The control input of the AUV is horizontal rudder 8°, vertical rudder -8°, and propeller speed 1000 rpm. The initial position of the AUV is [100100 10] T m, the forward initial velocity is 1.2 m / s, and the initial velocity in other directions is 0. The AUV motion trajectory obtained by simulation is as follows Figure 3The standard deviation of the simulated DVL velocity observation noise in three directions is [0.01 0.01 0.01] T m / s, the standard deviation of the attitude sensor roll angle and pitch angle observation is 0.03°, the standard deviation of the heading angle observation is 0.3°, and the standard deviation of the underwater acoustic signal propagation time is 100us. The simulated linear sound velocity profile function is c = -0.22z + 1542m / s, and its profile shape is as follows Figure 4 shown.

[0242] The Bellhop hydroacoustic simulation toolbox is used to generate hydroacoustic sound lines and real hydroacoustic signal propagation time during the AUV movement. The Bellhop model parameters used are shown in Table 2.

[0243] Table 2 Bellhop environmental parameters

[0244] Signal transmission frequency 10KHz depth 200m Seafloor modeling Acoustic elastic half space Minimum and maximum sound angle -90°，90° Maximum sound line distance 1.6 km Maximum sound depth 200 meters Number of beams 9001

[0245] For all methods, the width of the sliding window is set to 15, that is, the current moment and the 15 acoustic ranging observations before the current moment are used to solve the position coordinates of the AUV at the current moment. For the proposed Gauss-Newton method, the iterative convergence threshold is set to ε = 10 -5, the maximum number of iterations is set to 100. The initial position of the Gauss-Newton algorithm at each moment is set to the AUV coordinates estimated at the previous moment, and the initial position of the AUV at the initial moment is the position of the beacon plus a Gaussian random number with a standard deviation of 1. The R-SDR1 and R-SDR2 methods are solved by using the CVX toolbox of Matlab software, using the sedumi solver, and the solution accuracy is set to the highest, see the reference Michael Grant and Stephen Boyd.CVX: Matlab software for disciplined convex programming, version 2.0beta.https: / / cvxr.com / cvx, September 2013. and the reference Michael Grant and Stephen Boyd.Graphimplementations for nonsmooth convex programs, Recent Advances in Learning and Control (a tribute to M.Vidyasagar), V.Blondel, S.Boyd, and H.Kimura, editors, pages 95-110, Lecture Notes in Control and Information Sciences, Springer, 2008.http: / / stanford.edu / ~boyd / graph_dcp.html.; the positioning error and average positioning error of the four methods are used as performance evaluation criteria, and the calculation formulas of positioning error and average positioning error are:

[0246]

[0247]

[0248] Among them, K represents the number of single beacon positioning settlements; the positioning errors of the four methods are as follows: Figure 5 The average positioning error is shown in Table 3. The number of iterations of the proposed Gauss-Newton method at each moment is shown in Figure 6 shown.

[0249] according to Figure 5 As shown in Table 3, the positioning accuracy of the Gauss-Newton method proposed in the present invention is significantly better than that of the traditional single beacon positioning method that does not consider the refraction of sound propagation. Its positioning error is lower than that of other comparison methods at most times. This shows that the refraction of the sound propagation channel has a certain degree of influence on the performance of single beacon positioning, and the method proposed in the present invention can eliminate this influence. Figure 6 It can be seen that the method of the present invention can achieve a faster convergence speed, and the maximum number of iterations does not exceed 10 times. The method of the present invention has good real-time performance and can meet the real-time application requirements of single beacon positioning.

[0250] Table 3 Comparison of average positioning errors of different methods

[0251] Positioning method Average error (meters) T-GN 0.908090 R-WLS 2.396516 R-SDR1 1.518397 R-SDR2 2.374167

[0252] The above description is only a preferred specific implementation manner of the present invention, but the protection scope of the present invention is not limited thereto. Any modifications, equivalent substitutions and improvements made by any technician familiar with the technical field within the technical scope disclosed by the present invention and within the spirit and principles of the present invention should be covered within the protection scope of the present invention.

Claims

1. A single underwater beacon Gauss-Newton positioning method under linear acoustic profile conditions, characterized in that: The method includes: S1, take any point in the positioning area as the origin, set the east, north and sky directions as X, Y and Z axes respectively, and establish the underwater single beacon positioning geodetic coordinate system; take the AUV centroid as the origin, and set the front, right and bottom directions as x, y and z axes respectively, and establish the AUV body coordinate system; S2, constructing a dead reckoning model based on the velocity of the AUV in the body coordinate system detected by the DVL and the attitude of the AUV detected by the attitude sensor; S3, the hydroacoustic beacon periodically emits hydroacoustic signals, the beacon position coordinates are known, the AUV is equipped with a hydrophone, the propagation time of the hydroacoustic signal is measured, the sound velocity profile is modeled as a linear sound profile, and the known hydroacoustic profile parameter information and the hydroacoustic beacon position coordinates are combined to construct a hydroacoustic observation model under the linear sound profile condition, and the relationship between the signal propagation time and the AUV position is obtained; S4, integrating AUV dead reckoning information and using all acoustic signal propagation time observations within the time window width N, constructs a virtual long baseline positioning model under linear acoustic profile conditions; S5, calculates the gradient matrix of the propagation time of the underwater acoustic signal relative to the position variable, uses the Gauss-Newton method to iteratively solve the positioning problem, and obtains the position estimate of the AUV after the iteration converges.

2. The underwater single beacon Gauss-Newton positioning method under linear acoustic profile conditions according to claim 1, characterized in that: In step S2, a dead reckoning model is constructed, including: The relationship between the AUV position variables at time k and at time ki,i≥1 is: In the formula, x k ,y k ,z k is the real position coordinate of AUV at time k, x k-i is the x-direction position coordinate of the AUV at time ki, y k-i is the y-direction position coordinate of the AUV at time ki, z k-i is the z-direction position coordinate of the AUV at time ki, Δt is the discrete time interval, i is the time index variable, and v x,k-t is the x-direction velocity of AUV at time kt, v y,k-t is the y-direction velocity of AUV at time kt, v z,k-t is the z-direction velocity of the AUV at time kt.

3. The underwater single beacon Gauss-Newton positioning method under linear acoustic profile conditions according to claim 1, characterized in that: In step S3, a hydroacoustic observation model under the condition of linear acoustic profile is constructed, and the expression is: Where P b is the coordinate vector of the hydroacoustic beacon position, x b is the east coordinate of the hydroacoustic beacon, y b is the north coordinate of the hydroacoustic beacon, z b The celestial coordinates of the hydroacoustic beacon; The sound velocity profile is modeled as linear, and the sound velocity value is a linear function of the depth: v=az+b Where c is the depth at depth z, a is the sound velocity gradient, z is the depth, and b is the sound velocity at sea level; When the sound velocity profile is known, the sound velocity at the beacon is known and is denoted by c b , the speed of sound at the AUV position at time k is recorded as c k , the propagation time of the underwater acoustic signal at time k is recorded as t k , the relationship between horizontal distance and signal propagation time is: Among them, r k is the horizontal distance between the AUV and the beacon, and e is the natural base; Squaring both sides of the above equation yields: Simplifying gives: Substituting the sound velocity profile parameters into the equation: definition: In the formula, m k is a virtual observation, in which all parameters are known and the value of each parameter depends on the signal propagation time observation value and the sound profile parameter at the current moment; The relationship between the AUV position and the propagation time of the underwater acoustic signal is obtained as follows: 0=(x k -x b ) 2 +(y k -y b ) 2 +(z k -z b ) 2 -m k (a) k +b)。 4. The underwater single beacon Gauss-Newton positioning method under linear acoustic profile conditions according to claim 2, characterized in that: In step S4, a virtual long baseline positioning model under the linear acoustic profile condition is constructed, including: Define the time window width as N, use the propagation time of the underwater acoustic signal received from kN to k, solve the position coordinates of the AUV at k, and for the ki moment, apply the established observation equation of the underwater acoustic signal propagation time as follows: 0=(x k-i -x b ) 2 +(y k-i -y b ) 2 +(z k-i -z b ) 2 -m k-i (az k-i +b),for i=1,2…N In the formula, m k-i is the virtual observation at time ki; Substitute the AUV’s dead reckoning model, that is, the position relationship between time k and time ki, and we get: 0=(x k -x b,i ) 2 +(y k -y b,i ) 2 +(z k -z b,i ) 2 -m k-i (az k +b i ),for i=1,2…N in, In the formula, x b,i is the x-axis position coordinate of the ith virtual beacon, y b,i is the y-axis position coordinate of the i-th virtual beacon, z b,i is the z-axis position coordinate of the ith virtual beacon, b i is the surface sound velocity corresponding to the i-th virtual beacon; further, define: x b,0 =x b and b,0 =and b With b,0 =from b b0=b In the formula, x b,0 ,y b,0 ,z b,0 is the virtual beacon coordinate at time k-0, b0 is the virtual surface sound speed at time k-0; Combined with the observation equation at time k, all N+1 observation equations are expressed as follows: 0=(x k -x b,i ) 2 +(y k -y b,i ) 2 +(z k -z b,i ) 2 -m k-i (a) k +b i ), for i=0,1,2…N.

5. The underwater single beacon Gauss-Newton positioning method under linear acoustic profile conditions according to claim 2, characterized in that: In step S5, the positioning problem is iteratively solved using the Gauss-Newton method, including: The propagation time of the underwater acoustic signal contains observation noise, and the observation value of the propagation time of the underwater acoustic signal is recorded as: In the formula, is the observed value of the propagation time of the underwater acoustic signal, t k is the actual underwater acoustic signal propagation time, v k is the observation noise of the underwater acoustic signal propagation time, and its mean is modeled as 0 and the variance is Gaussian random variable; all N+1 underwater acoustic signal propagation time observations are combined into a vector All N+1 real underwater acoustic signal propagation times are combined into a vector t k , all N+1 underwater acoustic propagation time observation noises are integrated into a vector v k ,as follows: The observation noise v k Modeled as Gaussian noise with a mean of 0 and a variance matrix denoted as Q k , then the maximum likelihood estimation problem for single beacon positioning is: In the formula, is the maximum likelihood estimate of the AUV position at time k, arg max is the maximization problem, exp{} is the exponential function, and p k is the position of the AUV at time k, the superscript T represents the vector matrix transpose, t k (p k ) is the real propagation time of the underwater acoustic signal, and its value is a function of the position coordinates; the variance matrix Q k It is obtained according to the observed variance of the propagation time of the underwater acoustic signal at each moment, as follows: In the formula, σ k-N is the standard deviation of the observed variance at time kN, σ k-N+1 is the standard deviation of the observed variance at time k-N+1, σ k is the standard deviation of the observed variance at time k, diag(a) is a diagonal matrix with the elements of a as diagonal elements. According to the monotonicity of the logarithmic function, the above optimization problem is transformed into: Where arg min( ) is the minimization function; definition: Among them, S k for Choleskey decomposition, the optimization problem is transformed into: This problem is a nonlinear least squares optimization problem, and the Gauss-Newton method is used to solve it as follows: in, is the position estimate obtained after j+1 iterations, is the position estimate obtained at the jth iteration, is the Gauss-Newton position increment for the jth iteration, defined as: In the formula, is the gradient matrix of the signal propagation time to the position coordinates obtained at the jth iteration, is the acoustic signal propagation time calculated based on the AUV position estimated in the current iteration cycle, Defined as: In the formula, is the partial derivative of the propagation time of the underwater acoustic signal with respect to the position coordinate in the x direction, is the partial derivative of the underwater acoustic signal propagation time with respect to the position coordinate in the y direction, is the partial derivative of the underwater acoustic signal propagation time with respect to the position coordinate in the z direction, is the x, y, and z position estimate obtained at the jth iteration.

6. The underwater single beacon Gauss-Newton positioning method under linear acoustic profile conditions according to claim 4, characterized in that: Gradient matrix of position coordinates Each element is calculated as follows: According to the hydroacoustic beacon observation equation: 0=(x k -x b,i ) 2 +(y k -y b,i ) 2 +(z k -z b,i ) 2 -m k-i (a) k +b i ), for i=0,1,2…N Two sides x k ,y k ,z k Find the partial derivatives respectively, and we get: In the formula, is the partial derivative of the virtual observation with respect to the propagation time of the underwater acoustic signal, are the partial derivatives of the underwater acoustic signal propagation time with respect to the AUV position; According to m k-i The definition of , we get: Then we have: get:

7. The underwater single beacon Gauss-Newton positioning method under linear acoustic profile conditions according to claim 5, characterized in that: The propagation time of the underwater acoustic signal is calculated based on the AUV position estimated in the current iteration cycle. The method is: According to the equation: get: n i α i ＝(α i -1) 2 Where n i With α i are auxiliary parameters, and their definitions are: Solving the above quadratic equation yields: Right now: According to the above equation, we can get two different t k-i , since the real underwater acoustic signal propagation time is a positive number, select the one that is greater than 0; substitute the position coordinates obtained in the jth iteration into the above equation to obtain the underwater acoustic signal propagation time vector 8. The underwater single beacon Gauss-Newton positioning method under linear acoustic profile conditions according to claim 1, characterized in that: Set the iterative convergence criterion to when the position update amount When the second norm of is less than a certain threshold ε, the iteration is considered to converge and the single beacon positioning position estimate is output 9. An underwater single-beacon Gauss-Newton positioning system under linear acoustic profile conditions, characterized in that: The system implements the underwater single-beacon Gauss-Newton positioning method under linear sound profile conditions as described in any one of claims 1 to 8. The system obtains the AUV position coordinates through a positioning solution method based on the Gauss-Newton algorithm to address the problem of sound propagation channel refraction caused by the change of underwater sound velocity with depth under linear sound profile conditions.

10. The underwater single beacon Gauss-Newton positioning system under linear acoustic profile conditions according to claim 9, characterized in that: The system is provided with an AUV, which is equipped with a hydrophone, a DVL, a posture sensor and a CTD; the AUV detects the sound velocity profile parameters of the experimental area through the CTD; the hydroacoustic beacon periodically broadcasts hydroacoustic signals, and the AUV periodically observes its velocity in the body coordinate system through the DVL carried by the AUV, calculates the velocity of the AUV in the geodetic coordinate system in combination with the AUV attitude angle observed by the posture sensor, performs dead reckoning, and records and stores the dead reckoning data; After receiving the acoustic signal emitted by the acoustic beacon, the AUV obtains the acoustic signal propagation time based on the known acoustic signal emission time; When the AUV receives a certain number of hydroacoustic signals, it constructs a virtual long baseline array in combination with the stored dead reckoning data; Combined with the relationship between the underwater acoustic propagation time and the AUV position under the linear acoustic profile, a virtual long baseline positioning problem is constructed; the observation noise of the underwater acoustic signal propagation time is modeled as a Gaussian distribution, and a maximum likelihood estimation problem is constructed, which is then converted into a nonlinear least squares estimation problem; the gradient of the observation quantity with respect to the AUV position coordinates is solved based on the chain rule, and the nonlinear least squares estimation problem is iteratively solved based on the Gauss-Newton method, and the AUV position coordinates are obtained after the iteration converges.

Images