A polar region adaptive underwater acoustic positioning correction method based on a multi-dimensional environment coupling model

Abstract

The application discloses a polar region adaptive underwater acoustic positioning correction method based on a multi-dimensional environment coupling model, belongs to the technical field of underwater acoustic positioning, and is used for polar region underwater acoustic positioning and comprises the following steps: deploying a heterogeneous sensor network, constructing an environment perception layer, and eliminating sensor drift errors; a three-dimensional sound velocity model is established, a vortex dynamics equation, a Kriging interpolation algorithm and an adaptive Runge-Kutta method are combined to track a sound ray track; a time-frequency dual-domain joint detection algorithm is designed, a lattice IIR filter bank is combined to realize multipath suppression; a 9-dimensional robust UKF algorithm framework is constructed, and multi-source data fusion is realized in combination with federal filtering; and a geomagnetic matching and strain energy density early warning mechanism is integrated. Compared with the prior art, the application improves the positioning accuracy and stability, can more accurately predict a sound propagation path, improves positioning precision, realizes dynamic real-time calculation, can quickly and accurately process changing environmental data, and is suitable for complex polar region environments.

Description

TECHNICAL FIELD

[0001] The application discloses a polar adaptive underwater acoustic positioning correction method based on a multi-dimensional environment coupling model and belongs to the technical field of underwater acoustic positioning. BACKGROUND

[0002] In the field of ocean exploration and research, underwater acoustic positioning technology is a key means to realize underwater target positioning and tracking. With the gradual expansion of ocean development activities to the polar region, the complex marine environment in the polar region poses unprecedented challenges to underwater acoustic positioning technology.

[0003] Traditional underwater acoustic positioning methods are mostly designed based on the assumption of an ideal marine environment and have been widely used in conventional sea areas at middle and low latitudes. However, the polar region has unique hydrological conditions, such as low temperature, strong salinity gradient, and sea ice coverage. Low temperature leads to a decrease in seawater sound speed, which is significantly different from the sound speed distribution in conventional sea areas. The strong salinity gradient causes complex bending and refraction of sound ray propagation paths, which greatly deviates from the prediction of traditional models. The presence of sea ice not only reflects and scatters sound waves but also forms a complex reverberation environment, seriously interfering with the propagation and reception of underwater acoustic signals.

[0004] When existing underwater acoustic positioning technology is used in the polar region, the positioning accuracy decreases sharply and cannot meet the demand for high-precision positioning in activities such as polar ocean scientific research, resource exploration, and underwater facility maintenance. In addition, current positioning correction methods do not fully consider the coupling effect of multi-dimensional environmental factors in the polar region, and cannot adjust the positioning algorithm in real time according to the complex and variable environment in the polar region, resulting in poor reliability and stability of the positioning results. Therefore, it is urgent to develop an adaptive underwater acoustic positioning correction method that can adapt to the special environment in the polar region and comprehensively consider the coupling effect of multi-dimensional environmental factors, which is of great significance for improving the ocean exploration capability in the polar region and ensuring the smooth development of related activities. SUMMARY

[0005] The purpose of the present application is to provide a polar adaptive underwater acoustic positioning correction method based on a multi-dimensional environment coupling model to solve the problem that it is difficult to effectively deal with severe multipath interference, sudden changes in sound speed profile, and reference point drift when facing ice-covered sea areas in the prior art.

[0006] A polar adaptive underwater acoustic positioning correction method based on a multi-dimensional environment coupling model, comprising:

[0007] S1. Deploy an ice-floating buoy array, a vertical CTD chain, and an underwater acoustic beacon to form a heterogeneous sensor network, construct a spatiotemporally aligned environment perception layer through GPS / INS combined navigation, synchronous acquisition of temperature, salinity, and depth data, and acoustic signal emission, and eliminate the drift error of the piezoelectric ice stress sensor by using a dynamic calibration model;

[0008] S2. Establish a three-dimensional sound velocity model with horizontal temperature gradient compensation term, integrate vortex dynamics equations and Kriging interpolation algorithm to achieve high-precision sound velocity reconstruction, and track sound ray trajectories through adaptive Runge-Kutta method to cope with sudden changes in sound velocity in polar regions;

[0009] S3. Design a time-frequency dual-domain joint detection algorithm, using Gabor transform to extract the time delay spread and frequency domain attenuation characteristics of the ice layer reflection signal, and combining it with a lattice IIR filter bank to achieve multipath suppression;

[0010] S4. Construct a 9-dimensional robust UKF algorithm framework, integrate particle swarm optimization and quasi-Newton method, process non-Gaussian noise through Gaussian mixture model, and combine federated filtering to achieve multi-source data fusion.

[0011] S5. Establish an ice-water coupled dynamic model, compensate for buoy displacement through thermal expansion deformation formula, integrate geomagnetic matching and strain energy density early warning mechanism, and realize real-time correction of the reference coordinate system under ice layer drift.

[0012] S1 includes:

[0013] S1.1. The ice buoy array adopts a hexagonal topology layout. Each buoy integrates a dual-frequency GPS receiver, a MEMS inertial navigation module, and an ice thickness sensor. The dual-frequency GPS receiver receives L1 and L2 bands. The ice thickness sensor is based on the principle of electromagnetic induction. Data synchronization is achieved through the LoRa wireless network to achieve μs-level time synchronization. The NTPv4 protocol is used to correct clock drift.

[0014] Each vertical CTD chain is configured with multiple sensor nodes, which are connected by armored cables;

[0015] The underwater acoustic beacon is deployed at a depth of 50m below the seabed and is stabilized using an anchor chain system. The signal transmitted by the underwater acoustic beacon uses a nonlinear frequency modulation (NLFM) sequence.

[0016]

[0017] In the formula, s(t) represents the complex value of the nonlinear frequency modulation sequence at time t, A is the amplitude of the signal, j is the imaginary unit, f0 is the center frequency of the signal, k1 is the coefficient related to linear frequency modulation, k2 is the coefficient related to nonlinear frequency modulation, τ is the time delay, and T1 is the time period.

[0018] S1.2. Constructing a spatiotemporally aligned environment awareness layer includes establishing a spatiotemporal registration model:

[0019] t t =t0+Δt GPS +Δt b ;

[0020] In the formula, t tis the final alignment time after space-time registration, t0 is the original time, Δt GPS is the GPS time correction, Δt b is the sound propagation time delay compensation.

[0021] S1 includes:

[0022] S1.3. Based on the improved isolated forest algorithm, the temperature and salinity data are online cleaned to eliminate the sensor drift error, including establishing a dynamic calibration model:

[0023]

[0024] In the formula, T cal (z) is the calibrated temperature, T0 is the initial temperature, ΔT sensor is the sensor temperature difference, which is calibrated by the reference station on the ice in real time, z0 is the calibration decay depth, z is the real-time depth, e is the natural constant;

[0025] The clock synchronization of the underwater acoustic beacon adopts the two-way time transfer method to eliminate the clock deviation:

[0026]

[0027] In the formula, Δt is the time difference, T1, T2, T3 and T4 are the signal receiving and transmitting time stamps after space-time registration, δ1 and δ2 are the two measurement distances of the receiving and transmitting signals at two places, c is the sound speed, and

[0028] S1.4. Deploying a piezoelectric ice stress sensor array to measure the elastic modulus of the ice layer, and performing wavelet denoising processing on the data:

[0029]

[0030] In the formula, E ice is the elastic modulus of the ice layer, F is the external force on the ice body, A is the cross-sectional area of the ice body under force, ΔL is the change of the length of the ice body after being subjected to force, L0 is the original length of the ice body when it is not subjected to force, T e is the measured temperature after calibration, T ref is the reference temperature.

[0031] S2 includes:

[0032] S2.1. Establish a three-dimensional sound speed model, perform horizontal gradient field reconstruction, and introduce a vortex dynamics model to calculate the horizontal sound speed gradient:

[0033]

[0034] In the formula, is the horizontal gradient operator, R is an empirical constant, Γ is a regulating parameter, T is the temperature, and S is the salinity, is the temperature gradient in the horizontal direction, is the salinity gradient in the horizontal direction;

[0035] A 100m x 100m grid is generated using the Kriging interpolation algorithm, and the variogram γ(h) is calculated as:

[0036]

[0037] where h is the distance between two points;

[0038] S2.2. Adaptive step size Runge-Kutta method for solving the sound ray equation:

[0039]

[0040] where r is the depth of the sound ray, s is the arc length parameter along the sound ray, and c(z) is the sound speed at z;

[0041] The step size control condition is:

[0042]

[0043] where Δs n is the step size of the nth iteration, ε is the allowable error, c n is the high-order method solution, c n-1 is the low-order method solution, Δs max is the maximum step size.

[0044] S2 includes:

[0045] S2.3. Constructing the Mahalanobis distance detection statistic:

[0046] D 2 = (c obs -c model ) T ∑ -1 (c obs -c model );

[0047] where D is the Mahalanobis distance, c obs is the observed value of the sound speed, c model is the mean of the observed value of the sound speed, and Σ is the covariance matrix;

[0048] When D is greater than a preset value, the sound speed field reconstruction is triggered, and the dynamic update strategy uses full grid update every 5 minutes, or detects a temperature jump greater than 0.1℃ / min, and starts local grid encryption update;

[0049] S2.4. Introducing the digital elevation model DEM to correct the sound propagation equation at (x, y, z):

[0050]

[0051] where p is the sound pressure and γ(x,y) is the seabed slope function obtained by multibeam echo sounder.

[0052] S3 includes multipath suppression on the sound propagation equation obtained by S2.4;

[0053] S3.1. Multipath signal feature extraction, a double-channel Gabor transform is designed for time-frequency analysis:

[0054] G(t,f) = ∫r(τ)g(τ-t)e -j2πfτ dτ;

[0055] where G(t,f) is the result of Gabor transform, r(τ) is the original signal function, g(τ-t) is the window function, f is the signal frequency, and the multipath signal includes all the signals obtained by the sensors in S1;

[0056] Adaptive bandwidth design is adopted:

[0057]

[0058] where σ t is the time-bandwidth parameter, E is the feature quantity, and the ice layer reflection discrimination criterion is designed according to the main path signal and the ice layer reflection path;

[0059] S3.2. Multipath suppression filter design, the transfer function of the lattice IIR filter bank is constructed as:

[0060]

[0061] where H(z) is the transfer function of the lattice IIR filter bank, M is the total order of the filter, m is the m-th order, α m (k) and β m (k) are the m-th order polynomial coefficients at time k, z -1 is the delay operator;

[0062] Polynomial coefficient update rule:

[0063] α m (k+1) = α m (k) + μe(k)Re{s(k-m)};

[0064] where μ is the step factor, e(k) is the error signal, s(k-m) is the value of the input signal at time k-m, and Re{} is the real part operator.

[0065] S3 includes:

[0066] S3.3. Establishing motion platform compensation model:

[0067]

[0068] Where f is motion platform compensation, f is observation value, v is receiving platform motion velocity, and θ is signal incidence angle. comp obs

[0069] The velocity estimation of receiving platform motion velocity adopts adaptive α-β filter:

[0070]

[0071] Where v is velocity observation value at k moment, z is position measurement value at k moment, x is position estimation value at k moment, and Δt is time interval. k k

[0072] S3.4. Multipath signal angle of arrival estimation optimization, using MUSIC algorithm to improve spatial spectrum estimation:

[0073]

[0074] Where P(θ) is MUSIC spectrum, a(θ) is steering matrix corresponding to θ, E is noise subspace, and the array manifold adds ice layer reflection phase correction term. n

[0075] S4 includes:

[0076] S4.1. Extending state space model to 9 dimensions:

[0077] X = [x, y, z, v x , v y , v z , a bias , c bias , Δt1] T ;

[0078] Where x, y, and z are three-dimensional coordinate values, v x , v y , and v z are three-dimensional velocity values, a bias is accelerometer zero offset, c bias is sound speed error, obtained from multipath suppressed sound propagation equation, and Δt1 is clock bias.

[0079] S4.2. In UKF algorithm, the mean and covariance are calculated by Sigma point sampling, and the state prediction is:

[0080] ​​​​​​

[0081] where, is the prior estimate at time k, x k-1|k-1 is the state vector, u k is the control vector at time k;

[0082] The covariance prediction is:

[0083]

[0084] where, P k|k-1 is the prior covariance matrix at time k, F k-1 is the state transition matrix at time k-1, P k-1|k-1 is the covariance matrix, Q k is the process noise covariance at time k;

[0085] The particle swarm optimization (PSO) algorithm is introduced to find the optimal UKF parameters by updating the position and velocity of the particles. The particle velocity update is:

[0086]

[0087] where, is the ion position, w is the inertia weight, c1, c2 are learning factors, r1, r2 are random numbers, pbest i is the individual optimal position of the particle, gbest is the global optimal position;

[0088] The particle position update is:

[0089]

[0090] The dynamic error compensation matrix D k is introduced to correct the positioning results in real time:

[0091]

[0092] where, is the posterior estimate at time k, K k is the Kalman gain, z k is the observation vector, h() is the observation function.

[0093] S4 includes:

[0094] S4.3. The non-Gaussian noise processing is performed on to construct a mixed Gaussian model:

[0095]

[0096] where p(v1) is the probability density function of the random vector v1 for the Gaussian mixture model, ω i is the weight of the i-th Gaussian component, μ i is the mean of the i-th Gaussian component, ∑ i is the covariance matrix of the i-th Gaussian component;

[0097] is the probability density function of the i-th Gaussian component, and the weight update rule is:

[0098]

[0099] where, is the weight of the j-th Gaussian component of the kk+1 observation data point, r kk is the kk observation data point;

[0100] S4.4. A multi-source data fusion architecture is performed, and a federated filtering structure is designed:

[0101]

[0102] where, is the inverse matrix of the global estimation error covariance matrix, β ii is the weighting coefficient of the local filter of the ii sensor, the weight of the ii sensor, is the state estimation value of the ii sensor;

[0103] β ii is adaptively adjusted according to the signal-to-noise ratio:

[0104]

[0105] where SNR is the signal-to-noise ratio of the sensor signal, j1 is the total number of sensor signals, ii is the ii sensor signal, and the sensors include accelerometers and clocks.

[0106] S5 includes dynamic reference compensation, ice-water coupled dynamic modeling, and establishment of the buoy motion differential equation:

[0107]

[0108] where m is the mass of the buoy, F current is the water flow force, F wind is the force exerted by the wind on the surface of the buoy, F Coriolis is the Coriolis force, k 1 is the elastic coefficient, x 1 is the displacement of the buoy relative to a certain equilibrium position;

[0109] In S4, the On the basis of the above, when the early warning index is triggered or a geomagnetic anomaly occurs, the reference coordinate system is corrected in real time in combination with the displacement amount of the buoy;

[0110] S5.1. Ice layer crack early warning subsystem, develop an early warning index U based on strain energy density d :

[0111]

[0112] In the formula, V is the total volume of the ice layer, σ i is the stress in the small volume unit, ε i is the strain in the small volume unit, ΔV i is the volume change amount of the ice layer in the small volume unit, the stress and strain are obtained through E ice ;

[0113] S5.2. Establish a geomagnetic anomaly matching database:

[0114]

[0115] In the formula, M error is the geomagnetic anomaly matching error, B obs (k') is the geomagnetic intensity and direction information actually measured by the sample point, B model (k') is the geomagnetic vector value at the k'th sample point calculated according to the model, and n' is the total number of sample points.

[0116] Compared with the prior art, the present application has the following beneficial effects: the present application greatly improves the accuracy and stability of positioning by constructing a three-dimensional sound speed field model, designing a double-threshold time-frequency joint detection algorithm, and establishing a robust UKF algorithm fused with a particle swarm optimization; the present application can more accurately predict the sound propagation path by establishing an ice layer reflection feature library in combination with the constructed three-dimensional sound speed field model, thereby improving the positioning accuracy; the robust UKF algorithm fused with the particle swarm optimization of the present application introduces a dynamic error compensation matrix, realizes dynamic real-time solving, can quickly and accurately process the continuously changing environmental data, and is suitable for complex polar region environments. BRIEF DESCRIPTION OF DRAWINGS

[0117] Figure 1 is the technical flowchart of the present application. DETAILED DESCRIPTION

[0118] In order to make the purpose, technical scheme and advantages of the present application clearer, the technical scheme in the present application will be described clearly and completely below. Obviously, the described embodiments are part of the embodiments of the present application, rather than all the embodiments. Based on the embodiments in the present application, all other embodiments obtained by those skilled in the art without creative labor fall within the scope of protection of the present application.

[0119] The technical flowchart of the present application is shown as Figure 1 indicated, including deploying an ice floating buoy array, a vertical CTD chain and an underwater acoustic beacon composed of a heterogeneous sensor network, constructing a spatio-temporal aligned environment perception layer through GPS / INS combined navigation, simultaneous acquisition of temperature-salinity-depth data and acoustic signal emission, and eliminating sensor drift error using a dynamic calibration model; establishing a three-dimensional sound speed model containing a horizontal temperature gradient compensation term, fusing a vortex dynamics equation and a Kriging interpolation algorithm to realize high-precision sound speed reconstruction, tracking sound ray trajectory through an adaptive Runge-Kutta method to cope with polar sound speed mutation; designing a time-frequency dual-domain joint detection algorithm, extracting time delay expansion and frequency domain attenuation characteristics of ice layer reflection signals using Gabor transform, and realizing multipath suppression combining a lattice IIR filter bank; constructing a 9-dimensional robust UKF algorithm framework, fusing particle swarm optimization and quasi-Newton method, processing non-Gaussian noise through a mixture Gaussian model, and realizing multi-source data fusion combining federated filtering; establishing an ice-water coupled dynamics model, compensating buoy displacement through a thermal expansion deformation formula, integrating geomagnetic matching and strain energy density early warning mechanism, and realizing real-time correction of the reference coordinate system under ice layer drift.

[0120] An adaptive underwater acoustic positioning correction method in the polar region based on a multi-dimensional environment coupling model, comprising:

[0121] S1. Deploying an ice floating buoy array, a vertical CTD chain and an underwater acoustic beacon composed of a heterogeneous sensor network, constructing a spatio-temporal aligned environment perception layer through GPS / INS combined navigation, simultaneous acquisition of temperature-salinity-depth data and acoustic signal emission, and eliminating sensor drift error using a dynamic calibration model;

[0122] S2. Establishing a three-dimensional sound speed model containing a horizontal temperature gradient compensation term, fusing a vortex dynamics equation and a Kriging interpolation algorithm to realize high-precision sound speed reconstruction, tracking sound ray trajectory through an adaptive Runge-Kutta method to cope with polar sound speed mutation;

[0123] S3. Designing a time-frequency dual-domain joint detection algorithm, extracting time delay expansion and frequency domain attenuation characteristics of ice layer reflection signals using Gabor transform, and realizing multipath suppression combining a lattice IIR filter bank;

[0124] S4. Constructing a 9-dimensional robust UKF algorithm framework, fusing particle swarm optimization and quasi-Newton method, processing non-Gaussian noise through a mixture Gaussian model, and realizing multi-source data fusion combining federated filtering;

[0125] S5. Establishing an ice-water coupled dynamics model, compensating buoy displacement through a thermal expansion deformation formula, integrating geomagnetic matching and strain energy density early warning mechanism, and realizing real-time correction of the reference coordinate system under ice layer drift.

[0126] S1 includes:

[0127] S1.1. The ice buoy array adopts a hexagonal topology layout, each buoy is integrated with a dual-frequency GPS receiver, a MEMS inertial navigation module, and an ice thickness sensor. The dual-frequency GPS receiver receives L1 and L2 bands. The ice thickness sensor is based on the principle of electromagnetic induction. Data synchronization is achieved through LoRa wireless network to realize μs-level time synchronization. The NTPv4 protocol is used to correct clock drift.

[0128] Each vertical CTD chain is configured with multiple sensor nodes and connected through armored cables.

[0129] The underwater acoustic beacon is deployed at a depth of 50m from the seabed and is stabilized by an anchor chain system. The underwater acoustic beacon emits a signal using a nonlinear frequency modulation sequence NLFM as follows:

[0130]

[0131] In the formula, s(t) represents the complex value of the nonlinear frequency modulation sequence at time t, A is the amplitude of the signal, j is the imaginary unit, f0 is the center frequency of the signal, k1 is the coefficient related to linear frequency modulation, k2 is the coefficient related to nonlinear frequency modulation, τ is the time delay, and T1 is the time period.

[0132] S1.2. Building a spatio-temporal alignment environment perception layer includes establishing a spatio-temporal registration model:

[0133] t t =t0+Δt GPS +Δt b ;

[0134] In the formula, t t is the final aligned time after spatio-temporal registration, t0 is the original time, Δt GPS is the GPS time correction, and Δt b is the acoustic propagation time delay compensation.

[0135] S1 includes:

[0136] S1.3. Based on the improved isolated forest algorithm, online cleaning of temperature and salinity data is performed to eliminate sensor drift errors, including establishing a dynamic calibration model:

[0137]

[0138] In the formula, T cal (z) is the calibrated temperature, T0 is the initial temperature, ΔT sensor is the sensor temperature difference, which is calibrated in real time by the reference station on the ice, z0 is the calibration decay depth, z is the real-time depth, and e is the natural constant.

[0139] The clock synchronization of underwater acoustic beacon adopts the two-way time transfer method to eliminate clock bias:

[0140]

[0141] where Δt is the time difference, T1, T2, T3 and T4 are the time stamps of the signal transmission and reception after space-time registration, δ1 and δ2 are the measured distances at the two signal transmission and reception sites, c is the sound speed,

[0142] S1.4. Deploy a piezoelectric ice stress sensor array to measure the elastic modulus of the ice layer, and perform wavelet denoising on the data:

[0143]

[0144] where E ice is the elastic modulus of the ice layer, F is the external force on the ice body, A is the cross-sectional area of the ice body under force, ΔL is the change in length of the ice body after being subjected to force, L0 is the original length of the ice body when it is not under force, T e is the measured temperature after calibration, T ref is the reference temperature.

[0145] S2 includes:

[0146] S2.1. Establish a three-dimensional sound speed model, perform horizontal gradient field reconstruction, and introduce a vortex dynamics model to calculate the horizontal sound speed gradient:

[0147]

[0148] where is the horizontal gradient operator, R is an empirical constant, T is the adjustment parameter, T is the temperature, S is the salinity, is the horizontal temperature gradient, is the horizontal salinity gradient;

[0149] A 100m×100m grid is generated using the Kriging interpolation algorithm, and the variogram γ(h) is calculated as:

[0150]

[0151] where h is the distance between two points;

[0152] S2.2. Adaptive step size Runge-Kutta method to solve the sound ray equation:

[0153]

[0154] where r is the depth of the sound ray, s is the arc length parameter along the sound ray, and c(z) is the sound speed at z;

[0155] The step size control condition is:

[0156]

[0157] where Δs n is the step size of the nth iteration, ε is the tolerance error, c n is the high-order method solution, c n-1 is the low-order method solution, Δs max is the maximum step size.

[0158] S2 includes:

[0159] S2.3. Constructing the Mahalanobis distance detection statistics:

[0160] D 2 = (c obs - c model ) t ∑ -1 (c obs - c model );

[0161] where D is the Mahalanobis distance, c obs is the observed value of the sound speed, c model is the mean of the observed value of the sound speed, and Σ is the covariance matrix;

[0162] When D is greater than a preset value, the sound speed field reconstruction is triggered, and the dynamic update strategy adopts full grid update every 5 minutes, or when a temperature mutation greater than 0.1℃ / min is detected, local grid encryption update is started;

[0163] S2.4. Introducing a digital elevation model DEM to correct the sound propagation equation at (x, y, z):

[0164]

[0165] where p is the sound pressure, and γ(x, y) is the seabed slope function obtained by a multi-beam echo sounder.

[0166] S3 includes multi-path suppression of the sound propagation equation obtained in S2.4;

[0167] S3.1. Multi-path signal feature extraction, a double-channel Gabor transform is designed for time-frequency analysis:

[0168] G(t, f) = ∫r(τ)g(τ-t)e -j2πft dτ;

[0169] where G(t, f) is the result of Gabor transform, r(τ) is the original signal function, g(τ-t) is the window function, f is the signal frequency, and the multi-path signal includes all signals obtained by the sensors in S1;

[0170] Adaptive bandwidth design:

[0171]

[0172] where σ t is the time-bandwidth parameter, E is the eigenvalue, and the ice layer reflection discrimination criterion is designed according to the main path signal and the ice layer reflection path;

[0173] S3.2. Multipath suppression filter design, the transfer function of the lattice IIR filter bank is constructed as:

[0174]

[0175] where H(z) is the transfer function of the lattice IIR filter bank, M is the total order of the filter, m is the m-th order, α m (k) and β m (k) are the m-th order polynomial coefficients at time k, z -1 is the delay operator;

[0176] The polynomial coefficient update rule is:

[0177] α m (k+1) = α m (k) + μe(k)Re{s(k-m)};

[0178] where μ is the step factor, e(k) is the error signal, s(k-m) is the value of the input signal at time k-m, and Re{} is the real part operator.

[0179] S3 includes:

[0180] S3.3. Establishing a motion platform compensation model:

[0181]

[0182] where f comp is the motion platform compensation, f obs is the observation value, v is the receiving platform motion speed, and θ is the signal incidence angle;

[0183] The speed estimation of the receiving platform motion speed uses adaptive α-β filtering:

[0184]

[0185] where v k is the speed observation value at time k, z k is the position measurement value at time k, is the position estimation value at time k, Δt is the time interval, and α and β are dynamic adjustment parameters;

[0186] S3.4. Multipath signal angle of arrival estimation optimization, using MUSIC algorithm to improve spatial spectrum estimation:

[0187]

[0188] In the formula, P(θ) is the MUSIC spectrum, a(θ) is the steering matrix corresponding to θ, E n is the noise subspace, and the array manifold is added with ice layer reflection phase correction term.

[0189] S4 includes:

[0190] S4.1. Extend the state space model to 9 dimensions:

[0191] X = [x, y, z, v x , v y , v z , a bias , c bias , Δt1] T ;

[0192] In the formula, x, y and z are three-dimensional coordinate values, v x , v y and v z are three-dimensional velocity values, a bias is the accelerometer zero offset, c bias is the sound speed error, obtained from the sound propagation equation after multipath suppression, and Δt1 is the clock bias;

[0193] S4.2. In the UKF algorithm, the mean and covariance are calculated by Sigma point sampling, and the state prediction is:

[0194]

[0195] In the formula, is the prior estimation at time k, x k-1|k-1 is the state vector, and u k is the control vector at time k;

[0196] The covariance prediction is:

[0197]

[0198] In the formula, P k|k-1 is the prior covariance matrix at time k, F k-1 is the state transition matrix at time k-1, P k-1|k-1 is the covariance matrix, and Q k is the process noise covariance at time k;

[0199] The particle swarm optimization (PSO) algorithm is introduced, and the optimal UKF parameters are found by updating the position and velocity of the particles. is:

[0200]

[0201] In the formula, is the ion position, w is the inertial weight, c1, c2 is the learning factor, r1, r2 is a random number, pbest i is the particle individual optimal position, gbest is the global optimal position;

[0202] Particle position update is:

[0203]

[0204] Dynamic error compensation matrix D is introduced k , real-time correction of positioning results:

[0205]

[0206] In the formula, is the posterior estimation at time k, K k is the Kalman gain, z k is the observation vector, h() is the observation function.

[0207] S4 includes:

[0208] S4.3. Non-Gaussian noise processing is performed on , and a mixed Gaussian model is constructed:

[0209]

[0210] In the formula, p(v1) is the probability density function of the random vector v1 of the mixed Gaussian model, ω i is the weight of the i-th Gaussian component, μ i is the mean of the i-th Gaussian component, ∑ i is the covariance matrix of the i-th Gaussian component;

[0211] is the probability density function of the i-th Gaussian component, and the weight update rule is:

[0212]

[0213] In the formula, is the weight of the j-th Gaussian component of the kk+1 observation data point, r kk is the kk observation data point;

[0214] S4.4. Multi-source data fusion architecture is performed, and a federal filter structure is designed:

[0215]

[0216] wherein, is the inverse of the global estimation error covariance matrix, β ii is the weighting coefficient of the local filter of the ii-th sensor, the weight of the ii-th sensor, is the state estimation value of the ii-th sensor;

[0217] β ii is adjusted adaptively according to the signal-to-noise ratio:

[0218]

[0219] wherein, SNR is the signal-to-noise ratio of the sensor signal, j1 is the total number of sensor signals, ii is the ii-th sensor signal, and the sensors include an accelerometer and a clock.

[0220] S5 includes dynamic reference compensation, ice-water coupled dynamic modeling, and establishment of buoy motion differential equations:

[0221]

[0222] wherein, m is the mass of the buoy, F current is the water flow force, F wind is the force exerted by the wind on the surface of the buoy, F Coriolis is the Coriolis force, k 1 is the elastic coefficient, x 1 is the displacement of the buoy relative to a certain equilibrium position;

[0223] On the basis of the obtained in S4, when the early warning index is triggered or a geomagnetic anomaly occurs, the reference coordinate system is corrected in real time in combination with the displacement amount of the buoy;

[0224] S5.1. Ice layer crack early warning subsystem, develop early warning index U d based on strain energy density:

[0225]

[0226] wherein, V is the total volume of the ice layer, σ i is the stress within a small volume unit, ε i is the strain within a small volume unit, ΔVx is the volume change amount of the ice layer within a small volume unit, and the stress and strain are obtained through E ice ;

[0227] S5.2. Establish a geomagnetic anomaly matching database:

[0228]

[0229] In the formula, M error is the geomagnetic anomaly matching error, B obs (k') is the geomagnetic intensity and direction information actually measured at the sample point, B model (k') is the geomagnetic vector value at the k'th sample point calculated according to the model, and n' is the total number of sample points.

[0230] The ice buoy array adopts a hexagonal topology layout, and the length of the hexagon is 500 m. Each buoy is integrated with a dual-frequency GPS receiver, a MEMS inertial navigation module, and an ice thickness sensor. The dual-frequency GPS receiver receives L1 and L2 bands, and the carrier phase accuracy is 2 mm. The MEMS inertial navigation module has a gyro zero bias stability of 0.5° / h. The ice thickness sensor is based on the electromagnetic induction principle. Data synchronization is achieved through LoRa wireless network to realize μs-level time synchronization. The NTPv4 protocol is used to correct clock drift. Each vertical CTD chain is configured with 20 sensor nodes, and the node spacing is 5 m. The nodes are connected through armored cables.

[0231] A three-dimensional sound speed field model containing a horizontal temperature gradient compensation term is constructed, a double-threshold time-frequency joint detection algorithm is designed, a robust UKF algorithm fused with a particle swarm optimization is established, the limitations of underwater acoustic positioning technology in dealing with complex environments in ice-covered sea areas are broken through, the accurate grasp of the underwater acoustic signal propagation path in ice-covered sea areas is realized, the ice layer reflection misjudgment is reduced, and the dynamic real-time solution of positioning data is realized. An adaptive underwater acoustic positioning correction method based on a multi-dimensional environment coupling model is established, which significantly improves the accuracy, reliability and stability of underwater acoustic positioning in polar ice-covered sea areas, and provides more accurate and reliable positioning technology support for polar ocean exploration, resource exploration and other activities.

[0232] The above examples are only used to illustrate the technical solutions of the present application, and not to limit it. Although the present application has been described in detail with reference to the foregoing examples, those skilled in the art should understand that they can still modify the technical solutions recorded in the foregoing examples, or make equivalent replacements for part or all of the technical features, and these modifications or replacements do not make the essence of the corresponding technical solutions deviate from the scope of the technical solutions of the embodiments of the present application.

Claims

1. A polar region adaptive acoustic positioning correction method based on a multi-dimensional environment coupling model, characterized in that, Comprise: S1. Deploy the heterogeneous sensor network composed of ice buoy array, vertical CTD chain and underwater acoustic beacon, build the spatio-temporal aligned environmental perception layer through GPS / INS integrated navigation, synchronous acquisition of temperature-salinity-depth data and acoustic signal emission, and eliminate the drift error of piezoelectric ice stress sensor by using dynamic calibration model; S2. Establish a three-dimensional sound speed model containing a horizontal temperature gradient compensation term, fuse the vortex dynamics equation and Kriging interpolation algorithm to realize high-precision sound speed reconstruction, and track the sound ray trajectory by adaptive Runge-Kutta method to cope with the sudden change of sound speed in polar region; S3. Design a joint detection algorithm in time-frequency domain, extract the time delay expansion and frequency domain attenuation characteristics of ice layer reflection signal by Gabor transform, and realize multipath suppression by combining lattice IIR filter bank; S4. Construct a 9-dimensional robust UKF algorithm framework, fuse particle swarm optimization and quasi-Newton method, process non-Gaussian noise by Gaussian mixture model, and realize multi-source data fusion by combining federated filtering; S5. Establish an ice-water coupled dynamics model, compensate the buoy displacement by thermal expansion deformation formula, integrate the geomagnetic matching and strain energy density warning mechanism to realize real-time correction of the reference coordinate system under the drift of ice layer.

2. The polar region adaptive acoustic positioning correction method based on a multi-dimensional environment coupling model according to claim 1, characterized in that, S1 includes: S1.

1. The ice buoy array adopts a hexagonal topology layout, each buoy is integrated with a dual-frequency GPS receiver, a MEMS inertial navigation module, and an ice thickness sensor. The dual-frequency GPS receiver receives L1 and L2 bands. The ice thickness sensor is based on the electromagnetic induction principle, and data synchronization is achieved through a LoRa wireless network Level time synchronization, NTPv4 protocol is used to correct clock drift The vertical CTD chain is equipped with multiple sensor nodes connected by armored cable; The deployment depth of underwater acoustic beacon is 50m from the seabed, and the position is stabilized by anchor chain system. The underwater acoustic beacon emission signal adopts nonlinear frequency modulation sequence NLFM: ； wherein denotes the complex value of the non-linear frequency modulation sequence at time is the amplitude of the signal, is the imaginary unit, is the center frequency of the signal, is a coefficient related to the linear frequency modulation, is a coefficient related to the non-linear frequency modulation, is a time delay, is a time period;​ S1.

2. Building a spatio-temporal aligned environmental perception layer includes establishing a spatio-temporal registration model: ； wherein, is the final aligned time after spatio-temporal registration, is the original time, is the GPS time correction, is the acoustic propagation time delay compensation.

3. The polar region adaptive acoustic positioning correction method based on a multi-dimensional environment coupling model according to claim 2, characterized in that, S1 includes: S1.

3. Based on the improved isolated forest algorithm, the temperature and salinity data are cleaned online to eliminate sensor drift error, including establishing a dynamic calibration model: ； wherein, is the calibrated temperature, is the initial temperature, is the sensor temperature difference, calibrated in real time by the reference station on ice, is the calibrated decay depth, is the real-time depth, is the natural constant; The clock synchronization of underwater acoustic beacon adopts bidirectional time transfer method to eliminate clock deviation: ； wherein is the time difference, , , and is the signal transmission-reception time stamp after space-time registration, , is the two measured distances at the two signal transmission-reception sites, and c is the sound speed. S1.

4. Deploy piezoelectric ice stress sensor array to measure ice elastic modulus, and perform wavelet denoising on the data: ； wherein is the elastic modulus of the ice layer, is the external force on the ice body, is the cross-sectional area of the ice body under force, is the change in length of the ice body after being under force, is the original length of the ice body when not under force, is the temperature at the time of calibration, is the reference temperature.

4. The polar region adaptive acoustic positioning correction method based on a multi-dimensional environment coupling model according to claim 3, characterized in that S2 Comprise: S2.

1. Establish a three-dimensional sound speed model, reconstruct the horizontal gradient field, and introduce the vortex dynamics model to calculate the horizontal sound speed gradient: ； wherein is a horizontal gradient operator, is an empirical constant, is a tuning parameter, is temperature, is salinity, is a horizontal temperature gradient, is a horizontal salinity gradient; The Kriging interpolation algorithm is used to generate 100 m x 100 m grids, and the variogram is calculated is: ； wherein is the distance between two points; S2.

2. Adaptive step Runge-Kutta method is used to solve the sound ray equation: ； wherein is the depth of the sound ray, is the arc length parameter along the sound ray, is the sound speed at the point The step control condition is: ； wherein is the first iteration step size, is the tolerance, is the high order method solution, is the low order method solution, is the maximum step size.

5. The polar region adaptive acoustic positioning correction method based on a multi-dimensional environmental coupling model according to claim 4, wherein S2 Comprise: S2.

3. Construct Mahalanobis distance detection statistic: ； wherein is the Mahalanobis distance, is an observed value of the sound velocity, is a mean value of the observed values of the sound velocity, is a covariance matrix; When D is greater than the preset value, the sound speed field reconstruction is triggered, and the dynamic update strategy adopts every 5 minutes full grid update, or when the temperature mutation is greater than 0.1℃ / min, the local grid encryption update is started; S2.

4. Introduction of a digital elevation model, DEM, correction At the sound propagation equation: ； wherein is the sound pressure, is the sea floor slope function acquired by a multi-beam bathymeter.

6. The polar region adaptive acoustic positioning correction method based on a multi-dimensional environmental coupling model according to claim 5, characterized in that, S3 includes multipath suppression of the sound propagation equation obtained by S2.4; S3.

1. Multi-path signal feature extraction, design double-channel Gabor transform for time-frequency analysis: ； wherein is the result of the Gabor transform, is the original signal function, is the window function, is the signal frequency, the multipath signal comprising all signals obtained by the sensors in S1; Adaptive bandwidth design is adopted: ； wherein is a time-bandwidth parameter, is a characteristic quantity, which is designed to discriminate the ice layer reflection path from the main path signal according to an ice layer reflection discrimination criterion. S3.

2. Multi-path suppression filter design, construct lattice IIR filter bank transfer function: ； wherein is the transfer function of a lattice IIR filter bank, is the total filter order, is the order, and are the time polynomial coefficients, is the delay operator; The polynomial coefficient update rule is: ； wherein is a step size factor, is an error signal, is the value of the input signal at time is the value of the input signal at time is the take real operator.

7. The polar region adaptive acoustic positioning correction method based on a multi-dimensional environmental coupling model according to claim 6, characterized in that S3 Comprise: S3.

3. Establish a motion platform compensation model: ； wherein is the motion platform compensation, is the observation, is the received motion platform velocity, is the signal angle of incidence; The speed estimate of the platform motion velocity employs an adaptive filtering: ； wherein is a speed observation at time instant is a position measurement at time instant is a position estimate at time instant is a time interval , is a dynamic adjustment parameter S3.

4. Multi-path signal angle of arrival estimation optimization, MUSIC algorithm is used to improve spatial spectrum estimation: ； wherein is the MUSIC spectrum, is the steering matrix corresponding to is the noise subspace, is the noise subspace, the array manifold is added with a layer of ice to reflect the phase correction term.

8. The polar region adaptive acoustic positioning correction method based on a multi-dimensional environment coupling model according to claim 7, characterized in that S4 Comprise: S4.

1. Extend the state space model to 9 dimensions: ； wherein , and are three-dimensional coordinate values, , and are three-dimensional velocity values, is an accelerometer bias, is a sound speed error, obtained from the sound propagation equation after multipath mitigation, is a clock bias; S4.

2. In the UKF algorithm, the mean and covariance are calculated by Sigma point sampling, and the state prediction is: ； wherein is a priori estimate of the time instant, is the state vector, is is the control vector at the time instant The covariance prediction is: ； wherein is the a priori covariance matrix at time instant is the state transition matrix at time instant is the covariance matrix is the process noise covariance at time instant The PSO algorithm is introduced to find the optimal UKF parameters by updating the position and speed of particles is: ； wherein is an ion position, is an inertia weight, , is a learning factor, , is a random number, is a particle individual optimal position, is a global optimal position; Particle position update is: ； Introducing a dynamic error compensation matrix Real-time correction of the positioning result: ； where is the a posteriori estimate of the time instant, is the Kalman gain, is the observation vector, is the observation function.

9. The polar region adaptive acoustic positioning correction method based on a multi-dimensional environmental coupling model according to claim 8, characterized in that S4 Including: S4.

3. To perform non-Gaussian noise processing, construct a mixture Gaussian model: ； wherein is a probability density function of a random vector according to a Gaussian mixture model, is a weight of the th Gaussian component, is a mean of the th Gaussian component, is a covariance matrix of the th Gaussian component; is the probability density function of the ith Gaussian component, and the weight update rule is: ； wherein is is the weight of the th Gaussian component of the th observation data point; and th observation data point; and S4.

4. Multi-source data fusion architecture is carried out, and the federal filter structure is designed: ； ； wherein is the inverse of the global estimation error covariance matrix, is the weight of the local filter of the th sensor, is the weight of the th sensor, is the state estimate of the th sensor; Adaptation in terms of signal-to-noise ratio: ； wherein is a signal-to-noise ratio of the sensor signal, is a total number of sensor signals, is the sensor signal, the sensor comprising an accelerometer and a clock.

10. The polar region adaptive acoustic positioning correction method based on a multi-dimensional environment coupling model according to claim 9, characterized in that, S5 includes dynamic benchmark compensation, ice-water coupled dynamic modeling, and establishment of buoy motion differential equation: ； wherein is the buoyancy mass, is the water flow force, is the force exerted by the wind on the surface of the buoy, is the Coriolis force, is the elastic coefficient, is the displacement of the buoy with respect to a certain equilibrium position; On the basis of the S4 obtained When the early warning index is triggered or the geomagnetic anomaly occurs, the reference coordinate system is corrected in real time in combination with the displacement amount of the buoy. S5.

1. Ice layer crack early warning subsystem, development of early warning indicators based on strain energy density : ； wherein is the total volume of the ice layer, is the stress within a small volume element, is the strain within a small volume element, is the volume change of the ice layer within a small volume element, stress and strain are derived by wherein S5.

2. Establish a geomagnetic anomaly matching database: ； wherein is the geomagnetic anomaly matching error, is the geomagnetic vector value at the i-th sample point calculated according to the ice-water coupling dynamics model, is the geomagnetic vector value at the i-th sample point calculated according to the ice-water coupling dynamics model, is the geomagnetic vector value at the i-th sample point calculated according to the ice-water coupling dynamics model, is the total number of sample points.

Images