Real-time azimuth tracking model-based algorithm based on horizontal circular array

Abstract

The invention relates to a real-time azimuth tracking model-based algorithm based on a horizontal circular array. The real-time azimuth tracking model-based algorithm comprises the following steps: S1, determining equipment parameters and water area environment parameters for measuring a moving target; s2, constructing an azimuth angle state space equation and a calculation formula of a sensor model; s3, measuring the background noise level of the water area, and giving a noise coefficient of a measurement equation; s4, setting an azimuth angle initial value by observing the approximate direction of the moving target by a measurer, or estimating the azimuth angle initial value through a conventional beam forming method; and S5, performing calculation through an unscented Kalman filter azimuth tracking algorithm framework, and solving the azimuth angle of the moving target in real time. According to the unscented Kalman filtering framework, the input is sound pressure real-time data, the output is a moving target azimuth angle, the online real-time performance of operation is ensured, meanwhile, a measurement equation based on a plane wave propagation model is constructed according to the spatial position of a horizontal circular array, and by means of the isotropic property of the circular array, the real-time performance of the horizontal circular array is improved. And the reliability of estimation results of the algorithm in all directions is ensured.

Description

Technical Field

[0001] The invention relates to the technical field of underwater moving target azimuth tracking, in particular to a real-time azimuth tracking model-based algorithm based on a horizontal circular array. Background Art

[0002] The passive sonar target tracking problem has always been the key point of sonar signal processing. Because the premise of sonar target tracking is to complete the target azimuth estimation, and then design the tracking algorithm. Or use batch processing to take snapshot data one by one, perform offline segmented target azimuth estimation, and realize moving target tracking. And most algorithms are based on linear arrays. The directivity index of linear arrays is very low in the end-fire direction, while circular arrays are isotropic and have consistent directivity index in all directions, which is more suitable for azimuth tracking. Therefore, the online real-time algorithm technology for moving target tracking that combines the two traditional methods and uses circular arrays as the formation came into being.

[0003] Traditional moving target tracking technology estimates the initial arrival azimuth through beamforming, then establishes the target's velocity-related motion equation and the state space model of the physical characteristics, and realizes moving target tracking through the Kalman filter method. Alternatively, the target azimuth is solved by snapshots one by one through batch processing. Both methods face some important problems. The first method is that in practical applications, the target's motion model cannot be effectively modeled; the second method is that the algorithm lacks real-time performance, and the length of the snapshot directly affects the performance of the estimation. In addition, due to the complexity of the actual ocean environment, some snapshots will cause large estimation errors and "wild value" points.

[0004] Both of the above methods face some important problems. The first method is that in practical applications, the motion model of the target cannot be effectively modeled; therefore, the motion equation is inaccurate, making the Kalman filter algorithm inaccurate. The second method is that the algorithm lacks real-time performance, resulting in slow calculation results, and the length of the snapshot directly affects the performance of the estimation. Due to the complexity of the actual marine environment, some snapshots will cause large estimation errors and "wild value" points. Summary of the invention

[0005] Due to its isotropy, the horizontal circular array has an advantage over the linear array in horizontal azimuth estimation. Therefore, the purpose of the present invention is to develop a model-based azimuth tracking algorithm based on the horizontal circular array. According to the assumption that the azimuth angle of the target does not change much in a short snapshot, a state space model is introduced, and combined with the measurement equation, according to the nonlinear characteristics of the measurement equation and the azimuth angle, an unscented Kalman filter framework is introduced to form azimuth tracking of the moving target. In this way, the advantages of the two methods in the prior art can be effectively integrated. On the one hand, the establishment of the motion equation can be avoided, and the plane wave propagation model is used. This is because the sound pressure signal can be directly obtained. On the other hand, the state space model is introduced as a priori information to avoid the appearance of "wild value" points. The iterative algorithm operates under the framework of the unscented Kalman filter, and ensures the online real-time nature of the estimation. The unscented transformation simplifies the nonlinear relationship operation between the azimuth angle and the sound pressure, ensuring the low complexity of the estimation.

[0006] The technical solution of the present invention is:

[0007] A real-time azimuth tracking model-based algorithm based on a horizontal circular array comprises the following steps:

[0008] S1. Determine the equipment parameters and water environment parameters for measuring moving targets;

[0009] S2, constructing the calculation formula of the azimuth state space equation and the sensor measurement equation;

[0010] S3. Measure the background noise level of the water area and provide the noise coefficient of the measurement equation;

[0011] S4. Setting the initial azimuth angle by observing the approximate direction of the moving target by the surveyor, or estimating the initial azimuth angle by conventional beamforming method;

[0012] S5. Calculate the azimuth angle of the moving target in real time by using the unscented Kalman filter azimuth tracking algorithm framework.

[0013] In step S1, the equipment parameters include the radius of the horizontal circular array, the number of horizontal circular array elements, and the sampling frequency; the water environment parameters include the underwater sound speed and the frequency domain range of the underwater target;

[0014] At the same time, the moving target is divided into a single-frequency target, a narrowband target or a broadband target, and the frequency point selection and calculation method is as follows:

[0015] (1) For single-frequency targets, the frequency corresponding to the single frequency point is directly used for algorithm calculation to achieve target tracking;

[0016] (2) For narrowband targets, 5 to 10 frequency points in the narrowband frequency band are uniformly selected, and the algorithm operation is performed at each frequency point. The average of the iterative estimates of each frequency point is taken as the tracking result;

[0017] (3) For broadband targets, select 30 to 50 frequency points in a manner that the frequency points are densely selected in the low frequency band and sparsely selected in the high frequency band. Perform algorithm calculations on the frequency points and take the average of the iterative estimates of each frequency point as the tracking result.

[0018] In step S2, the azimuth state space equation is constructed based on the short-term azimuth invariance characteristic of the moving target, as shown in the following equation:

[0019]

[0020] In the formula, θ represents the azimuth angle, represents the rate of change of azimuth angle, t is time, θ t represents the azimuth at time t, represents the rate of change of the azimuth at time t, ΔT represents the short-term snapshot time, ω represents the state-space model noise, and ω(t) represents the model noise that varies with time;

[0021] The sensor measurement equation is a measurement equation based on the geometric characteristics of the horizontal circular array, as shown in the following equation:

[0022]

[0023] Where M is the number of elements in the horizontal circular array, π is the pi constant, f is the frequency value corresponding to the frequency point, r is the radius of the horizontal circular array, c is the water acoustic velocity, Indicates the azimuth angle and the rate of change of the azimuth angle, with a range of (0,360).

[0024] In step S5, the specific steps include:

[0025] S51, using random perturbation to obtain the state space model noise ω, the azimuth state space equation established in step S2 is input into various parameters (including the azimuth initial value θ 0 , time parameter t, short-time snapshot time ΔT) and then calculate the output azimuth angle θ, input the azimuth angle θ and other sensor parameters into the sensor measurement equation f(x) and then add the sensor noise v to calculate the frequency domain phase difference y of the adjacent hydrophone receiving signals;

[0026] S52, initial azimuth angle θ 0 Perform UT transformation, enter each parameter in the azimuth state space equation, calculate and weight it to get the weighted azimuth The weighted azimuth Perform UT transformation and transform the transformed The frequency domain phase difference of the hydrophone receiving signal is calculated and weighted by the sensor measurement equation f(x) in step S2 and other sensor parameters.

[0027] S53, parameter y and parameter Output parameter e after difference calculation Where e is defined as the residual, which is the difference between the actual y output by the measurement equation and the parameter θ calculated by the state space model. The parameter Ke is then output after the parameter e is multiplied by the Kalman gain value K(t). Ke represents the angle difference that should be compensated for the angle given by the state space equation calculated by the Kalman gain according to the phase residual information e given by the measurement equation and the state equation.

[0028] S54, weighted azimuth The sum of the parameters Ke is calculated and the real-time calculated azimuth angle is output.

[0029] The beneficial effects of the present invention are:

[0030] 1. The unscented Kalman filter framework of the present invention has real-time sound pressure data as input and azimuth of moving targets as output, thus ensuring the online real-time performance of the calculation;

[0031] 2. The present invention constructs a measurement equation based on a plane wave propagation model according to the spatial position of the horizontal circular array, incorporates it into an unscented Kalman filter framework, and utilizes the isotropic property of the circular array to ensure the reliability of the algorithm's estimation results in all directions. BRIEF DESCRIPTION OF THE DRAWINGS

[0032] Figure 1 The present invention is a workflow diagram of a real-time azimuth tracking model-based algorithm based on a horizontal circular array;

[0033] Figure 2 Schematic diagram of the process of obtaining the noise coefficient of the measurement equation in step S2 of the present invention;

[0034] Figure 3 This is a measurement example of the background noise level in the water area in step S3 of the present invention (a time-frequency diagram of the background noise in a water area of ​​Qiandao Lake in a quiet state in January 2024);

[0035] Figure 4 This is a calculation example of estimating the initial value of the azimuth angle by the conventional beamforming method in step S4 of the present invention (the azimuth angle estimation result of a water area of ​​Qiandao Lake by the conventional beamforming method in January 2024);

[0036] Figure 5 This is a MATLAB data processing platform interface diagram of a real-time azimuth tracking model-based algorithm based on a horizontal circular array of the present invention;

[0037] Figure 6This is an interface diagram showing the processing results of the algorithm of the present invention (a comparison diagram between circular array beamforming and Kalman filtering when the sound source moves toward the buoy). DETAILED DESCRIPTION

[0038] The present invention will be further described below in conjunction with the accompanying drawings and embodiments:

[0039] To solve the problem that the azimuth tracking of underwater moving targets requires real-time online processing while most algorithms use offline batch processing, the present invention designs an unscented Kalman filter azimuth tracking algorithm for moving targets based on a horizontal circular array. By using an unscented Kalman filter and introducing a state-space model and a plane wave measurement equation, real-time estimation of the azimuth of the moving target can be achieved based on the sound pressure data received by the circular array.

[0040] In the framework of unscented Kalman filtering, the real-time azimuth tracking algorithm based on horizontal circular array of the present invention mainly includes three modules: a state space equation module with state space model as the core, a measurement equation module with sensor model as the core, and an unscented Kalman filtering module with Kalman gain weighting as the core. Figure 1 As shown, the present invention first adds untraceable transformation to the first two modules to deal with the nonlinear mapping problem, and finally calculates the new azimuth in real time by calculating the Kalman gain weighting, thereby forming a real-time azimuth tracking algorithm.

[0041] Specifically, the real-time azimuth tracking model-based algorithm based on the horizontal circular array provided by the present invention adopts the unscented Kalman filter azimuth tracking algorithm framework to measure the azimuth of the moving target, and includes the following steps:

[0042] S1. Determine the equipment parameters and water environment parameters for measuring moving targets;

[0043] S2, constructing the calculation formula of the azimuth state space equation and the sensor measurement equation;

[0044] S3. Measure the background noise level of the water area and provide the noise coefficient of the measurement equation;

[0045] S4. Setting the initial azimuth angle by observing the approximate direction of the moving target by the measurement personnel (the direction can be roughly given by the human sense), or estimating the initial azimuth angle by conventional beamforming method;

[0046] S5. After the preparations from step S1 to step S4 are completed, the azimuth angle of the moving target can be calculated in real time by using the unscented Kalman filter azimuth tracking algorithm framework.

[0047] Furthermore, in step S1, the equipment parameters include the radius of the horizontal circular array, the number of elements of the horizontal circular array, and the sampling frequency, and the water environment parameters include the underwater sound speed and the frequency domain range of the underwater target;

[0048] At the same time, the moving target is divided into a single-frequency target, a narrowband target or a broadband target, and the frequency point selection and calculation method is as follows:

[0049] (1) For single-frequency targets, the frequency corresponding to the single frequency point is directly used for algorithm calculation to achieve target tracking;

[0050] (2) For narrowband targets, 5 to 10 frequency points in the narrowband frequency band are uniformly selected, and the algorithm operation is performed at each frequency point. The average of the iterative estimates of each frequency point is taken as the tracking result;

[0051] (3) For broadband targets, select 30 to 50 frequency points in a manner that the frequency points are densely selected in the low frequency band and sparsely selected in the high frequency band. Perform algorithm calculations on the frequency points and take the average of the iterative estimates of each frequency point as the tracking result.

[0052] According to the motion characteristics of the azimuth of the moving target at sea, a state space model is set. In the present invention, it is assumed that the azimuth does not change much in a short period of time, so the state space model is a unit matrix.

[0053] Furthermore, in step S2, the azimuth state space equation is constructed based on the short-term azimuth invariance characteristic of the moving target, as shown in the following equation:

[0054]

[0055] In the formula, θ represents the azimuth angle, represents the rate of change of azimuth angle, t is time, θ t represents the azimuth at time t, represents the rate of change of the azimuth at time t, ΔT represents the short-term snapshot time, ω represents the state-space model noise, and ω(t) represents the model noise that varies with time;

[0056] Furthermore, in step S2, the sensor measurement equation is a measurement equation based on the geometric characteristics of the horizontal circular array, as shown in the following equation:

[0057]

[0058] Where M is the number of elements in the horizontal circular array, π is the pi constant, f is the frequency value corresponding to the frequency point, r is the radius of the horizontal circular array, c is the water acoustic velocity, Indicates the azimuth angle and the rate of change of the azimuth angle, with a range of (0,360).

[0059] Furthermore, in step S3, the specific steps of measuring the noise coefficient of the equation and measuring the background noise level of the water area are:

[0060] S31, measuring the noise coefficient of the equation, comprising the following steps:

[0061] (1) Determine the processing length of the sampling window, where the number of sampling points is equal to the sampling frequency fs;

[0062] (2) Determine the overlap ratio of the sampling window length;

[0063] (3) Perform Fourier transform on the data within each window length to obtain the spectrum, and take the absolute value of the amplitude value corresponding to each frequency point and square it to obtain the power spectrum of the data within each window length;

[0064] (4) taking the average value of the data in each window length to obtain the average power spectrum;

[0065] (5) Averaging all the frequency values ​​in the average power spectrum and taking the square root to obtain the noise coefficient of the measurement equation;

[0066] S32, measuring the background noise level of the water area, comprising the following steps:

[0067] (1) First calculate the average power spectrum value of each frequency point to obtain the decibel level. The average power spectrum can be taken in the form of 10*log10;

[0068] (2) Calculate the background noise level of the measured water area. The specific calculation method is: the average power spectrum value of each frequency point (10*log10) + the sensitivity decibel level of the hydrophone + the decibel level of the preamplifier = the background noise level of the measured water area.

[0069] Furthermore, in step S4, an initial azimuth angle θ is estimated by an artificial or conventional beamforming method. 0 , or the way to estimate the initial azimuth angle by conventional beamforming method is as follows:

[0070] (1) Selecting a processing time T, performing Fourier transform on the data within the processing time T, and obtaining a frequency spectrum within the processing time T;

[0071] (2) Take the amplitude value corresponding to the target frequency point. Assuming there are M array elements, we get an M×1 data vector, defined as X data ;

[0072] (3) Define the driving vector. The driving vector of the horizontal circle array is:

[0073] a=[e -j×2×π×f×r×cos(θ) / c ,e -j×2×π×f×r×cos(θ-360 / M) / c ,...,e -j×2×π×f×r×cos(θ-360(M-1) / M) / c ] T ;

[0074] (4) Transpose and multiply the data vector obtained in the second step by the driving vectors in each direction (0 to 360 degrees) defined in the third step, that is, In the formula, H represents conjugate transpose; take the absolute value of the result corresponding to each direction (0 degrees to 360 degrees) and square it, that is, Get the power decibel value in each direction;

[0075] (5) Take the angle corresponding to the maximum power decibel value as the initial azimuth angle θ estimated by the conventional beamforming method 0 .

[0076] Furthermore, in step S5, the specific steps include:

[0077] S51, using random perturbation to obtain the state space model noise ω, the azimuth state space equation established in step S2 is input into various parameters (including the azimuth initial value θ 0 , time parameter t, short-time snapshot time ΔT) and then calculate the output azimuth angle θ,

[0078] The azimuth angle θ and other sensor parameters are input into the sensor measurement equation f(x) and then added with the sensor noise v to calculate the frequency domain phase difference y of the adjacent hydrophone receiving signals;

[0079] S52, initial azimuth angle θ 0 Perform UT transformation, enter each parameter in the azimuth state space equation, calculate and weight it to get the weighted azimuth Then the weighted azimuth Perform UT transformation and transform the transformed The frequency domain phase difference of the hydrophone receiving signal is calculated and weighted by the sensor measurement equation f(x) in step S2 and other sensor parameters.

[0080] S53, parameter y and parameter Output parameter e after difference calculation Where e is defined as the residual, which is the difference between the actual y output by the measurement equation and the parameter θ calculated by the state space model. The parameter Ke is then output after the parameter e is multiplied by the Kalman gain value K(t). Ke represents the angle difference that should be compensated for the angle given by the state space equation calculated by the Kalman gain according to the phase residual information e given by the measurement equation and the state equation.

[0081] S54, weighted azimuth The sum of the parameters Ke is calculated and the real-time calculated azimuth angle is output.

[0082] Furthermore, in step S51, the frequency domain phase difference y of the adjacent hydrophone received signals is the sum of the output value of the azimuth angle θ substituted into the sensor measurement equation f(x) and the sensor measurement noise v, as shown in the following formula:

[0083]

[0084] In the formula, y represents the frequency domain phase difference of the signals received by adjacent hydrophones, M is the number of elements in the horizontal circular array, and the circular array is numbered from 0 to M-1, then: y 1 represents the frequency domain phase difference between the signals received by hydrophone No. 1 and hydrophone No. 0, y M-1 represents the frequency domain phase difference between the signals received by the M-1th hydrophone and the M-2th hydrophone, π is the pi constant, f is the frequency value corresponding to the frequency point, r is the radius of the horizontal circular array, c is the water acoustic velocity, θ is the azimuth, v is the sensor measurement noise, t is time, and v(t) represents the measurement noise that changes with time.

[0085] Furthermore, in step S52, based on the initial azimuth angle θ 0 Will (At high sampling rates, the instantaneous rate of change is very small, then x=θ 0 ) is converted into 5 sigma points after UT transformation, and the corresponding weights of each sigma point are [w 0 ,w 1 ,w 2 ,w 3 ,w 4 ], calculate each sigma point through the azimuth state space equation to get the angle value sigma of the corresponding sigma point i (i=0,…,4), and then through the weight [w 0 ,w 1 ,w 2 ,w 3 ,w 4 ] is weighted to obtain the weighted (priori estimate), which is specifically expressed as:

[0086]

[0087] Among them, for the initial azimuth angle θ 0 The formula for UT transformation and weight calculation of 5 sigma points is as follows:

[0088] Sigma 0 =μ x

[0089]

[0090] In the formula, Indicates the azimuth angle and the rate of change of azimuth angle, μ x represents the mean of the azimuth angle, S x Represents the covariance matrix R = x H The Cholesky decomposition result of x, (S x ) 1 For S x The first column of the matrix, (S x ) 2 For S x The second column of the matrix, n+κ is the empirical value, in this paper, κ is 1 and n+κ is 3;

[0091] Among them, the angle value sigma of each sigma point i When (i=0, ..., 4) is substituted into the azimuth state-space equation, the state-space model noise ω is not considered, as shown in the following formula:

[0092]

[0093] In the formula, θ represents the azimuth angle, represents the rate of change of azimuth angle, t is time, θ t represents the azimuth at time t, represents the rate of change of the azimuth at time t, and ΔT represents the short-term snapshot time.

[0094] Furthermore, in step S52, based on the weighted azimuth Will (At high sampling rates, the instantaneous rate of change Very small, at this time ) is transformed into 5 sigma points (with the initial azimuth angle θ 0 The UT transformation method and weights are the same), and each sigma point is calculated through the measurement equation f(x) to obtain the frequency domain phase difference y of the adjacent hydrophone receiving signals at the corresponding sigma point. y is calculated through the same weight [w 0 ,w 1 ,w 2 ,w 3 ,w 4 ] is weighted to obtain the weighted The specific expression is:

[0095]

[0096] The frequency domain phase difference y of the signals received by adjacent hydrophones is based on the weighted azimuth Convert UT into 5 sigma points and convert the angle value sigma of each sigma point i Substitute (i=0, ..., 4) into the sensor measurement equation f(x) to calculate, ignoring the sensor measurement noise v, and obtain the frequency domain phase difference y, as shown in the following formula:

[0097]

[0098] In the formula, y represents the frequency domain phase difference of the signals received by adjacent hydrophones, M is the number of elements in the horizontal circular array, and the circular array is numbered from 0 to M-1, then: y 1 represents the frequency domain phase difference between the signals received by hydrophone No. 1 and hydrophone No. 0, y M-1 It represents the frequency domain phase difference between the received signals of the M-1th hydrophone and the M-2th hydrophone, π is the pi constant, f is the frequency value corresponding to the frequency point, r is the horizontal circular array radius, c is the water acoustic velocity, and θ is the azimuth angle.

[0099] Furthermore, in step S53, the Kalman gain value is calculated as follows:

[0100]

[0101] In the formula, K(t) represents the value of Kalman gain, which changes with time t. represents the azimuth and the rate of change of azimuth, represents the phase residual information given by the measurement equation and the state equation, R represents the covariance matrix, R xε represents the covariance matrix of the azimuth angle and the phase residual given by the measurement equation and the state equation, R εε Represents the phase residual autocorrelation matrix given by the measurement equation and the state equation;

[0102] At the same time, ξ is used to represent the residual information of the sigma point, and the expression is as follows:

[0103]

[0104] In the formula, ξ represents the residual information of the sigma point, U represents the value of the sigma point after being mapped to the nonlinear measurement equation, w represents the weighted value, i represents the index number of the sigma point, t represents the time, ξ i Represents the residual information of the i-th sigma point, w i Represents the weighted value of the i-th sigma point, U i Represents the i-th sigma point by sigma i The value after mapping to the nonlinear measurement equation, namely U i =f(sigma i), (t|t-1) means predicting the value at the next time t at time t-1, U i (t|t-1) means that the i-th sigma point predicts the U at the next time t at time t-1 i value;

[0105] On this basis, ξ is used to represent the residual information of the sigma point instead of ε, and the content in the Kalman gain is replaced by:

[0106]

[0107] In the formula, (t|t-1) means the value predicted at the next time t at time t-1, R xε represents the covariance matrix of the azimuth angle and the phase residual given by the measurement equation and the state equation, R εε Represents the phase residual autocorrelation matrix given by the measurement equation and the state equation; R xξ Represents the covariance matrix of the residuals of the azimuth and sigma points, R ξξ Represents the autocorrelation matrix of the residuals at the sigma point, w i Represents the weighted value of the i-th sigma point, sigma i represents the i-th sigma point, ξ' represents the transpose of ξ, i Represents the residual information of the i-th sigma point.

[0108] In the present invention, in step S3, the process of obtaining the noise coefficient of the measurement equation is as follows: Figure 2 As shown: The first step is to determine the processing window length num_window, for example, 1 second, that is, one second is a sampling window length, and the number of sampling points is equal to the sampling frequency fs. The second step is to confirm the window length overlap rate overlap, for example, the overlap rate is 0.5, which represents that the first 50 percent of the next sampling window is consistent with the last 50 percent of the previous sampling window. The third step is to perform Fourier transform on the data in each window length to obtain the spectrum, and take the absolute value and square the amplitude value corresponding to each frequency point to obtain the data power spectrum in each window length. The fourth step is to take the average value of the data in each window length to obtain the average power spectrum. The fifth step is to average all the frequency point values ​​in the average power spectrum and take the square root to obtain the measurement equation noise coefficient;

[0109] The process of measuring the background noise level of the water area is as follows: first, the average power spectrum is taken in the form of 10*log10 to obtain the decibel level, and then the background noise level of the measured water area is calculated based on the sensitivity decibel level of the hydrophone and the decibel level of the preamplifier. The specific calculation method is: the average power spectrum value of each frequency point (10*log10) + the sensitivity decibel level of the hydrophone + the decibel level of the preamplifier = the background noise level of the measured water area. Figure 3The time-frequency diagram of the background noise level in a water area of ​​Qiandao Lake in January 2024 is given (where the sensitivity decibel level of the hydrophone is -192dB and the decibel level of the preamplifier is 30dB).

[0110] In step S4, the approximate direction of the moving target can be estimated manually by a surveyor, or the azimuth angle can be estimated by a conventional beamforming method. Figure 4 The azimuth estimation results of a water area in Qiandao Lake in January 2024 using the conventional beamforming method are given in the figure. It can be seen from the figure that the initial azimuth angle θ given at the highest peak is 0 It is 63.4 degrees.

[0111] The present invention uses the data processing platform of MATLAB software to realize the algorithm simulation task, and completes the calculation formula construction of the azimuth state space equation and the sensor measurement equation in step 2 through the platform, and the establishment of the unscented Kalman filter azimuth tracking algorithm in step S5. The implementation function of the unscented Kalman filter algorithm can be based on the following: Figure 5 The MATLAB language shown in the figure is used to write the program. On this basis, the equipment parameters and water environment parameters based on the horizontal circular array in step S1, the noise coefficient of the measurement equation in step S3, and the initial value of the azimuth angle θ obtained in step S4 are input. 0 , the azimuth of the moving target can be estimated in real time.

[0112] In the specific implementation, under the scenario of sound source frequency of 1270Hz, underwater sound speed of 1500m / s, circular array radius of 1.5m, number of circular array elements of 16, sampling frequency of 51200Hz, the tracking data of the sound source azimuth moving from 47 degrees to 60 degrees is processed in real time, and the following is obtained: Figure 6 The Kalman filter algorithm results are compared with the conventional beamforming method. Figure 6 It can be seen that the Kalman filter algorithm processing results obtained by introducing the state space model in the present invention are compared with the conventional beamforming method, which effectively suppresses the appearance of "wild value" points, improves data reliability, and has good real-time performance.

[0113] In summary, the present invention has developed a model-based azimuth tracking algorithm based on the horizontal circular array. According to the assumption that the azimuth angle of the target does not change much in a short snapshot, a state space model is introduced, and combined with the measurement equation, according to the nonlinear characteristics of the measurement equation and the azimuth angle, an unscented Kalman filter framework is introduced to form azimuth tracking of the moving target.

[0114] The present invention is characterized in that:

[0115] (1) The horizontal circular array is isotropic and has greater application potential in azimuth estimation than the linear array. In addition, the algorithm of the present invention performs Gaussian approximation (traceless transformation) on the relationship between the azimuth angle and the received sound pressure, which is an online iterative calculation, thereby ensuring the real-time performance of the result.

[0116] (2) Compared with the traditional azimuth tracking method which only uses the plane wave measurement equation for matched filtering, or only introduces the Kalman filter method after combining the motion equation with the state equation, the difference of the present invention is that the plane wave measurement equation is combined with the state equation, the nonlinear unscented Kalman filter is introduced, the azimuth is directly iteratively estimated, thereby avoiding secondary calculation, and the state space model is introduced as prior information, thereby avoiding the appearance of "wild value" points, and the calculation efficiency and reliability are better;

[0117] (3) The algorithm of the present invention is a new model-based algorithm for real-time azimuth tracking based on a horizontal circular array, and has good application prospects in real-time estimation of the azimuth angle of a moving target.

[0118] It should be noted that the above embodiments are merely preferred embodiments of the present invention and are not intended to limit the protection scope of the present invention. All equivalent modifications made on the basis of the above embodiments belong to the protection scope of the present invention.

Claims

1. A real-time azimuth tracking model-based algorithm based on a horizontal circular array, characterized in that: The steps include: S1. Determine the equipment parameters and water environment parameters for measuring moving targets; S2, constructing the calculation formula of the azimuth state space equation and the sensor measurement equation; S3. Measure the background noise level of the water area and provide the noise coefficient of the measurement equation; S4. Setting the initial azimuth angle by observing the approximate direction of the moving target by the surveyor, or estimating the initial azimuth angle by conventional beamforming method; S5. Calculate the azimuth angle of the moving target in real time by using the unscented Kalman filter azimuth tracking algorithm framework.

2. The real-time azimuth tracking model-based algorithm based on a horizontal circular array according to claim 1, characterized in that: In step S1, the equipment parameters include the radius of the horizontal circular array, the number of horizontal circular array elements, and the sampling frequency; the water environment parameters include the underwater sound speed and the frequency domain range of the underwater target; The moving targets are classified into single-frequency targets, narrow-band targets or broadband targets, and the frequency selection and calculation method are as follows: (1) For single-frequency targets, the frequency corresponding to the single frequency point is directly used for algorithm calculation to achieve target tracking; (2) For narrowband targets, 5 to 10 frequency points are uniformly selected from the narrowband frequency band, and the algorithm is calculated at each frequency point. The average of the iterative estimates of each frequency point is taken as the tracking result. (3) For broadband targets, select 30 to 50 frequency points in a manner that the frequency points are densely selected in the low frequency band and sparsely selected in the high frequency band. Perform algorithm calculations on the frequency points and take the average of the iterative estimates of each frequency point as the tracking result.

3. The real-time azimuth tracking model-based algorithm based on a horizontal circular array according to claim 1, characterized in that: In step S2, the azimuth state space equation is constructed based on the short-term azimuth invariance characteristic of the moving target, as shown in the following equation: In the formula, θ represents the azimuth angle, represents the rate of change of azimuth angle, t is time, θ t represents the azimuth at time t, represents the rate of change of the azimuth at time t, ΔT represents the short-term snapshot time, ω represents the state-space model noise, and ω(t) represents the model noise that varies with time; The sensor measurement equation is a measurement equation based on the geometric characteristics of the horizontal circular array, as shown in the following equation: Where M is the number of elements in the horizontal circular array, π is the pi constant, f is the frequency value corresponding to the frequency point, r is the radius of the horizontal circular array, c is the water acoustic velocity, Indicates the azimuth angle and the rate of change of the azimuth angle, with a range of (0,360).

4. The real-time azimuth tracking model-based algorithm based on a horizontal circular array according to claim 1, characterized in that: In step S3, the specific steps of measuring the noise coefficient of the equation and measuring the background noise level of the water area are: S31, measuring the noise coefficient of the equation, comprising the following steps: (1) Determine the processing length of the sampling window, where the number of sampling points is equal to the sampling frequency fs; (2) Determine the overlap ratio of the sampling window length; (3) Perform Fourier transform on the data within each window length to obtain the spectrum, and take the absolute value of the amplitude value corresponding to each frequency point and square it to obtain the power spectrum of the data within each window length; (4) taking the average value of the data in each window length to obtain the average power spectrum; (5) Averaging all the frequency values ​​in the average power spectrum and taking the square root to obtain the noise coefficient of the measurement equation; S32, measuring the background noise level of the water area, comprising the following steps: (1) First calculate the average power spectrum value of each frequency point to obtain the decibel level; (2) Calculate the background noise level of the measured water area. The specific calculation method is: average power spectrum value of each frequency point + sensitivity decibel level of the hydrophone + decibel level of the preamplifier = background noise level of the measured water area.

5. The real-time azimuth tracking model-based algorithm based on a horizontal circular array according to claim 1, characterized in that: In step S4, the method of estimating the initial value of the azimuth angle by using a conventional beamforming method is as follows: (1) Selecting a processing time T, performing Fourier transform on the data within the processing time T, and obtaining a frequency spectrum within the processing time T; (2) Take the amplitude value corresponding to the target frequency point. Assuming there are M array elements, we get an M×1 data vector, defined as X data ; (3) Define the driving vector. The driving vector of the horizontal circle array is: a=[and -j×2×π×f×r×cos(θ) / c ,And -j×2×π×f×r×cos(θ-360 / M) / c ,...,And -j×2×π×f×r×cos(θ-360(M-1) / M) / c ] T ; (4) Transpose and multiply the data vector obtained in the second step by the driving vectors in each direction (0 to 360 degrees) defined in the third step, that is, In the formula, H represents conjugate transpose; take the absolute value of the result corresponding to each direction (0 degrees to 360 degrees) and square it, that is, Get the power decibel value in each direction; (5) Take the angle corresponding to the maximum power decibel value as the initial value of the azimuth angle estimated by the conventional beamforming method.

6. The real-time azimuth tracking model-based algorithm based on a horizontal circular array according to claim 3, characterized in that: In step S5, the specific steps include: S51, using random perturbation to obtain the state space model noise ω, after inputting various parameters into the azimuth state space equation of step S2, the output azimuth θ is calculated, the azimuth θ and other sensor parameters are input into the sensor measurement equation f(x) of step S2, and then the sensor noise v is added to calculate the frequency domain phase difference y of the adjacent hydrophone receiving signals; S52, perform UT transformation on the initial azimuth angle θ0, calculate and weight the weighted azimuth angle after inputting each parameter into the azimuth angle state space equation The weighted azimuth Perform UT transformation and transform the transformed The frequency domain phase difference of the hydrophone receiving signal is calculated and weighted by the sensor measurement equation f(x) in step S2 and other sensor parameters. S53, parameter y and parameter The parameter e is output after the difference calculation, and then the parameter e is multiplied by the Kalman gain value to output the parameter Ke; S54, weighted azimuth The sum of the parameters Ke is calculated and the real-time calculated azimuth angle is output.

7. The real-time azimuth tracking model-based algorithm based on a horizontal circular array according to claim 6, characterized in that: In step S51, the frequency domain phase difference y of the signals received by adjacent hydrophones is the sum of the output value of the sensor measurement equation f(x) substituted into the azimuth angle θ and the sensor measurement noise v, as shown in the following formula: Where y represents the frequency domain phase difference between the signals received by adjacent hydrophones, M is the number of elements in the horizontal circular array, and the circular array is numbered from 0 to M-1. Then, y1 represents the frequency domain phase difference between the signals received by hydrophone No. 1 and hydrophone No. 0, and y M-1 represents the frequency domain phase difference between the signals received by the M-1th hydrophone and the M-2th hydrophone, π is the pi constant, f is the frequency value corresponding to the frequency point, r is the radius of the horizontal circular array, c is the water acoustic velocity, θ is the azimuth, v is the sensor measurement noise, t is time, and v(t) represents the measurement noise that changes with time.

8. The real-time azimuth tracking model-based algorithm based on a horizontal circular array according to claim 6, characterized in that: In step S52, the initial azimuth angle θ0 is converted into 5 sigma points after UT transformation, and the corresponding weights corresponding to each sigma point are [w0, w1, w2, w3, w4], and each sigma point is calculated by the azimuth state space equation to obtain the angle value sigma of the corresponding sigma point i (i=0,…,4), and then weighted by weights [w0,w1,w2,w3,w4] to get the weighted The specific expression is: Among them, the UT transformation of the initial azimuth angle θ0 and the weight calculation formula of the five sigma points are as follows: sigma0=μ x In the formula Indicates the azimuth angle and the rate of change of azimuth angle, μ x represents the mean of the azimuth angle, S x Represents the covariance matrix R = x H The Cholesky decomposition result of x, (S x )1 is S x The first column of the matrix, (S x )2 is S x The second column of the matrix, n+κ, is the empirical value; Among them, the angle value sigma of each sigma point i When (i=0, ..., 4) is substituted into the azimuth state-space equation, the state-space model noise ω is not considered, as shown in the following formula: In the formula, θ represents the azimuth angle, represents the rate of change of azimuth angle, t is time, θ t represents the azimuth at time t, represents the rate of change of the azimuth at time t, and ΔT represents the short-term snapshot time.

9. The real-time azimuth tracking model-based algorithm based on a horizontal circular array according to claim 7, characterized in that: In step S52, based on the weighted azimuth Will After the same UT transformation, it is converted into 5 sigma points. After each sigma point is calculated through the measurement equation f(x), the frequency domain phase difference y of the adjacent hydrophone receiving signals at the corresponding sigma point is obtained. y is weighted by the same weight [w0, w1, w2, w3, w4] to obtain the weighted The specific expression is: Among them, the frequency domain phase difference y of the adjacent hydrophone receiving signals is based on the weighted azimuth Convert UT into 5 sigma points and convert the angle value sigma of each sigma point i Substitute (i=0, ..., 4) into the sensor measurement equation f(x) to calculate, ignoring the sensor measurement noise v, and obtain the frequency domain phase difference y, as shown in the following formula: Where y represents the frequency domain phase difference between the signals received by adjacent hydrophones, M is the number of elements in the horizontal circular array, and the circular array is numbered from 0 to M-1. Then, y1 represents the frequency domain phase difference between the signals received by hydrophone No. 1 and hydrophone No. 0, and y M-1 It represents the frequency domain phase difference between the received signals of the M-1th hydrophone and the M-2th hydrophone, π is the pi constant, f is the frequency value corresponding to the frequency point, r is the horizontal circular array radius, c is the water acoustic velocity, and θ is the azimuth angle.

10. The real-time azimuth tracking model-based algorithm based on a horizontal circular array according to claim 6, characterized in that: In step S53, the Kalman gain value is calculated as follows: In the formula, K(t) represents the value of Kalman gain, which varies with time t, represents the azimuth and the rate of change of azimuth, Represents the phase residual information given by the measurement equation and the state equation, R represents the covariance matrix, R xε represents the covariance matrix of the azimuth angle and the phase residual given by the measurement equation and the state equation, R εε Represents the phase residual autocorrelation matrix given by the measurement equation and the state equation; ξ is used to represent the residual information of the sigma point, and the expression is as follows: In the formula, ξ represents the residual information of the sigma point, U represents the value of the sigma point after being mapped to the nonlinear measurement equation, w represents the weighted value, i represents the index number of the sigma point, t represents the time, ξ i Represents the residual information of the i-th sigma point, w i Represents the weighted value of the i-th sigma point, U i Represents the i-th sigma point by sigma i The value after mapping to the nonlinear measurement equation, namely U i =f(sigma i ), (t|t-1) means predicting the value at the next time t at time t-1, U i (t|t-1) means that the i-th sigma point predicts the U at the next time t at time t-1 i value; Use ξ to represent the residual information of the sigma point instead of ε, and the content in the Kalman gain is replaced by: In the formula, (t|t-1) means the value predicted at the next time t at time t-1, R xε represents the covariance matrix of the azimuth angle and the phase residual given by the measurement equation and the state equation, R εε Represents the phase residual autocorrelation matrix given by the measurement equation and the state equation; R xξ Represents the covariance matrix of the residuals of the azimuth and sigma points, R ξξ Represents the autocorrelation matrix of the residuals at the sigma point, w i Represents the weighted value of the i-th sigma point, sigma i represents the i-th sigma point, ξ' represents the transpose of ξ, i Represents the residual information of the i-th sigma point.

Images