A highly robust adaptive direction of arrival estimation method

Abstract

This invention relates to direction-of-arrival (DOA) estimation methods, specifically a highly robust adaptive ODA method that addresses the technical problems of low accuracy, poor robustness, and instability in time-varying environments caused by impulse noise interference in traditional MUSIC algorithms. This invention employs a logarithmic compression operator to compress the array signal, significantly improving the overall signal-to-noise ratio, suppressing impulse noise interference, and enhancing the accuracy of ODA estimation. Furthermore, it introduces an adaptive adjustment mechanism for compression parameters, avoiding the mismatch between compression intensity and noise environment that may occur with fixed-parameter compression. This mechanism can handle impulse noise interference from different sources, improving the robustness and stability of ODA estimation in complex time-varying environments. Simultaneously, the dual-branch ODA estimation framework based on the PAST and EWMA algorithms ensures both the real-time adaptive adjustment of compression parameters and improved ODA accuracy.

Description

Technical Field

[0001] This invention relates to direction-of-arrival (DOA) estimation methods, specifically to a highly robust adaptive DOA estimation method. Background Technology

[0002] Direction-of-arrival (DOA) estimation is a technique in array signal processing that determines the target's location by analyzing the phase difference and intensity changes of the signal. It is widely used in radar target detection, acoustic target detection, localization, and tracking. Compared with traditional algorithms such as beamforming, minimum variance distortionless response, and rotation-invariant factor estimation of signal parameters, the Multiple Signal Classification (MUSIC) algorithm, through orthogonality analysis of the signal and noise subspaces, exhibits superior performance in multi-target resolution and complex environments with low signal-to-noise ratios, and has therefore attracted widespread attention.

[0003] The high-resolution direction-of-arrival (DOA) estimation performance of the traditional MUSIC algorithm relies on two core assumptions: 1) the Gaussianity assumption, which states that the noise in the input signal follows a zero-mean stationary normal distribution; and 2) the orthogonality assumption, which states that the signal subspace and noise subspace in the input signal are orthogonal. However, in an airborne acoustic environment, noise often fails to meet these assumptions. On the one hand, the transmission loss of signals of different frequencies in the air varies, causing the energy distribution of the input signal to change dynamically with the propagation distance. On the other hand, random interference such as natural wind noise, vehicle horns, and other short-duration impulse noise exhibits significant non-stationarity and non-Gaussianity, leading to an asymmetric long-tailed distribution of noise amplitude. This directly violates the Gaussianity assumption of the MUSIC algorithm regarding the input signal, resulting in increased estimation bias in the noise subspace, which in turn reduces the accuracy of DOA estimation and weakens its robustness and stability. Summary of the Invention

[0004] The purpose of this invention is to solve the technical problems of low accuracy, poor robustness and stability of direction-of-arrival estimation in time-varying environments caused by impulse noise interference in traditional MUSIC algorithms, and to provide a highly robust adaptive direction-of-arrival estimation method.

[0005] To achieve the above objectives, the technical solution adopted by the present invention is as follows:

[0006] A robust adaptive direction-of-arrival estimation method, characterized by the following steps:

[0007] Step 1: Construct a logarithmic compression operator for time-varying non-stationary and non-Gaussian noise; then, given the compression candidate parameter set for strong noise state, moderate noise state and stationary noise state, as well as the compression parameter update conditions, the range of noise subspace drift, the initial compression parameters of the logarithmic compression operator, the initial signal subspace estimation, the initial covariance matrix and the initial noise subspace projection matrix.

[0008] Step 2: Acquire the array signal. Based on the initial signal subspace estimation, use the PAST algorithm to process the array signal to obtain the signal subspace estimation. Calculate the noise subspace projection matrix based on the signal subspace estimation. Then, calculate the noise subspace drift based on the noise subspace projection matrix and the initial noise subspace projection matrix. If the noise subspace drift is within the range of noise subspace drift variation, proceed to Step 3; otherwise, calculate the instantaneous covariance matrix based on the array signal, and then proceed to Step 6.

[0009] Step 3: Compress the array signal according to the logarithmic compression operator corresponding to the initial compression parameters to obtain the compressed array signal, and calculate the instantaneous covariance matrix based on the compressed array signal; then, project the compressed array signal onto the noise subspace according to the noise subspace projection matrix to obtain the projected signal, and calculate the projection residual; perform JB test on the projection residual to obtain the JB test probability value of the current noise distribution, and calculate the current composite objective function based on it;

[0010] Step 4: Based on the JB test probability value of the current noise distribution, determine the current noise state, and select the corresponding compression candidate parameter set from the compression candidate parameter set given in Step 1 for strong noise state, medium noise state and stationary noise state according to the current noise state. Then, according to the logarithmic compression operator corresponding to the candidate parameter in the compression candidate parameter set, obtain the noise distribution JB test probability value and composite objective function of each candidate parameter in the compression candidate parameter set according to the method in Step 3.

[0011] Step 5: Based on the JB test probability value of the noise distribution and the composite objective function of each candidate parameter in the compressed candidate parameter set, as well as the current JB test probability value of the noise distribution and the current composite objective function, determine whether the candidate parameters in the compressed candidate parameter set meet the compression parameter update conditions. If the compression parameter update conditions are met, take the candidate parameter corresponding to the minimum value of the composite objective function in the compressed candidate parameter set as the initial compression parameter, and then execute Step 6; if the compression parameter update conditions are not met, then execute Step 6.

[0012] Step 6: Based on the instantaneous covariance matrix and the initial covariance matrix, the EWMA algorithm is used to recursively obtain the covariance matrix. Then, eigenvalue decomposition and MUSIC spatial spectrum peak search are performed on the covariance matrix to obtain and output the direction of arrival.

[0013] Step 7: Use the signal subspace estimate as the initial signal subspace estimate, the noise subspace projection matrix as the initial noise subspace projection matrix, and the covariance matrix as the initial covariance matrix. Then return to step 2 to obtain the next frame of array signal until the direction of arrival estimation is completed.

[0014] Further, in step 1, the logarithmic compression operator is:

[0015]

[0016] Where f(x(t)) is the logarithmic compression operator, x(t) is the array signal, t represents the current frame; sgn() is the sign function, α(t) and β(t) are the initial compression strength parameter and the initial compression region parameter, respectively, and ε is an infinitesimal quantity of the logarithmic compression operator;

[0017] The conditions for updating the compression parameters are as follows:

[0018] If the noise distribution JB test probability value of the candidate parameters in the compressed candidate parameter set and the composite objective function satisfy the following conditions, then the compressed parameters should be updated:

[0019] ,and

[0020] or,

[0021]

[0022] in, , These are the current composite objective function and the composite objective function of the candidate parameters in the compressed candidate parameter set, respectively. , These are candidate parameters from the compression candidate parameter set, namely the compressibility strength parameter and the compression region parameter; , These are the JB test probability values ​​of the current noise distribution and the JB test probability values ​​of the noise distribution of the candidate parameters in the compressed candidate parameter set, respectively, ε. J ε p These are the threshold values ​​for the composite objective function and the probability value of the JB test, respectively.

[0023] The compression candidate parameter sets for the strong noise state, the medium noise state, and the stationary noise state are given as follows:

[0024] The candidate set of compression parameters for the strong noise state is: the candidate set of compression parameters for the compression intensity parameter remains unchanged, and the compression region parameter is... ;

[0025] The candidate set of compression parameters under moderate noise conditions is: the candidate set of compression parameters with the compression intensity parameter remaining constant and the compression region parameter being... ;

[0026] The candidate sets of compression parameters for stationary noise states are as follows: the candidate sets of compression intensity parameters and compression region parameters are respectively... , ;

[0027] in, Adjust the step size of the compressibility strength parameter for the current steady noise state. , , These are the parameter adjustment step sizes for the current strong noise state, medium noise state, and steady noise state, respectively. .

[0028] Furthermore, in step 4, the specific method for determining the current noise state is as follows:

[0029] If the current noise distribution has a JB test probability value p JB <0.01, or the JB test probability value of the current noise distribution p JB If the value is less than 0.01 and the noise subspace drift δ(t) is greater than 1.1, then it is a strong noise state;

[0030] If the JB test probability value of the current noise distribution is 0.01 ≤ p JB <0.05, or the JB test probability value of the current noise distribution is 0.01 ≤ p JB If the noise level is less than 0.05 and the noise subspace drift is 1.1 ≥ δ(t) > 0.9, then it is considered a medium noise state.

[0031] If the current noise distribution has a JB test probability value p JB ≥0.05, or the JB test probability value p of the current noise distribution. JB If the noise level is ≥0.05 and the noise subspace drift δ(t) ≤0.9, then it is a stationary noise state.

[0032] Furthermore, the current composite objective function is calculated using the following formula:

[0033]

[0034] Where J0 is the core objective function for noise compression. Penalty constraints for the objective function;

[0035] J JB This is the JB test normalization term for the noise distribution. ; , These represent the current skewness and current kurtosis, respectively, and ω1, ω2, and ω3 represent J... JB , , The weights;

[0036] These are candidate parameters from the set of candidate compression parameters for the compressed region parameters. β serves as a reference value for adjusting the parameters of the compression region. max Let λ1 and λ2 be the maximum values ​​in the set of candidate compression parameters for the compression region, respectively. , The constraint strength, An indicator function for adjusting the step size of the current compression region parameters. , These are the compression region parameters for the next frame of the array signal. Adjust the step size for the parameters of the current compression region.

[0037] Furthermore, step 2 specifically involves:

[0038] Step 2.1: Obtain the array signal. Based on the initial signal subspace estimation, use the PAST algorithm to process the array signal. Through step-by-step iterative processing, the signal subspace estimation is obtained using the following formula:

[0039]

[0040]

[0041]

[0042] Among them, U s (t) represents the signal subspace estimate, U s (t-1) represents the initial signal subspace estimate; , is the conjugate transpose matrix of the initial signal subspace estimation. The product of the array signal x(t); k(t) is the gain vector. Let be the conjugate transpose of k(t);

[0043] λ is the forgetting factor of the PAST algorithm. Let be the conjugate transpose of y(t), and P(t) and P(t-1) be the inverses of the autocorrelation matrices of y(t) and y(t-1), respectively.

[0044] Step 2.2: Based on the signal subspace estimation, and utilizing the orthogonality between the signal and noise subspaces, calculate the noise subspace projection matrix using the following formula:

[0045]

[0046] Among them, P n (t) is the projection matrix of the noise subspace, I N Let N be an N-order identity matrix, where N is the number of array elements used to acquire the array signal. For U s The conjugate transpose of (t);

[0047] Step 2.3: Based on the noise subspace projection matrix and the initial noise subspace projection matrix, calculate the noise subspace drift using the following formula:

[0048]

[0049] Where δ(t) is the noise subspace drift, P n (t-1) is the initial noise subspace projection matrix. It is the Frobenius norm;

[0050] Step 2.4: If the noise subspace drift is within the range of noise subspace drift, proceed to step 3; otherwise, calculate the instantaneous covariance matrix based on the array signal, and then proceed to step 6.

[0051] Further, in step 6, the covariance matrix is ​​calculated using the following formula:

[0052]

[0053] Where R(t) is the covariance matrix, R(t-1) is the initial covariance matrix, and λ R R is the forgetting factor of the covariance matrix. x (t) is the instantaneous covariance matrix.

[0054] Further, in step 3, the JB test probability value of the current noise distribution is obtained by the following formula:

[0055]

[0056] Where JB(t) is the JB statistic of the JB test, and N tot This represents the total number of signals in the array. Let represent the cumulative distribution function of a chi-square distribution with 2 degrees of freedom. The array signal acquired for the nth array element, where n is an integer and 0 ≤ n ≤ N. σ is the sample mean of the array signal. x denoted as the sample standard deviation of the array signal.

[0057] Furthermore, step 5 also includes:

[0058] If the compression parameter update condition is met, then the adjustment step size of the current noise state's candidate compression parameter set is updated as follows:

[0059]

[0060] If the compression parameter update condition is not met, then the adjustment step size of the current noise state's compression candidate parameter set is updated as follows:

[0061]

[0062] in, , These are the adjustment step sizes for the compression intensity parameter and compression region parameter of the next frame, respectively. Adjust the step size for the current compressive strength parameter. , These represent the maximum adjustment step sizes for the candidate parameter sets of the compressibility strength parameter and the compression region parameter, respectively. , These are the minimum values ​​of the adjustment step size for the candidate parameter sets of the compressibility strength parameter and the compression region parameter, respectively.

[0063] Furthermore, step 1 also includes setting the number of frozen frames for updating compression parameters;

[0064] Step 3 specifically involves:

[0065] The array signal is compressed using the logarithmic compression operator corresponding to the initial compression parameters to obtain a compressed array signal. The instantaneous covariance matrix is ​​then calculated based on the compressed array signal. It is determined whether the number of frames of the array signal has reached the number of frozen frames for updating the compression parameters. If so, the compressed array signal is projected onto the noise subspace according to the noise subspace projection matrix to obtain a projected signal, and the projection residual is calculated. The JB test is performed on the projection residual to obtain the JB test probability value of the current noise distribution, and the current composite objective function is calculated based on it. Then, step 4 is executed. If not, step 6 is executed.

[0066] Further, in step 1, the range of the noise subspace drift is [0.3, 1.2], the initial signal subspace estimate is a random matrix, the initial covariance matrix R(t-1) = 0.01 × I, and the initial noise subspace projection matrix P n (t-1)=0.01×I, where I is the identity matrix;

[0067] In step 5, the The value is 0.08-0.2. It is 0.01-0.1.

[0068] Compared with the prior art, the present invention has the following beneficial effects:

[0069] 1. The present invention provides a highly robust adaptive direction-of-arrival estimation method, which uses a logarithmic compression operator to compress the array signal. By approximately linearly amplifying the low-amplitude signal and significantly compressing the high-amplitude signal, the overall signal-to-noise ratio of the array signal is significantly improved, the quality of array signal acquisition in low signal-to-noise ratio and complex sound field environments is improved, thereby suppressing impulse noise interference and improving the accuracy of direction-of-arrival estimation.

[0070] 2. The present invention provides a highly robust adaptive direction-of-arrival estimation method, which introduces an adaptive parameter adjustment mechanism of a logarithmic compression operator to avoid the mismatch between compression intensity and noise environment that may exist when compression is performed with fixed parameters. It can cope with impulse noise interference from different sources and improve the robustness and stability of direction-of-arrival estimation under complex time-varying environments.

[0071] 3. This invention provides a highly robust adaptive direction-of-arrival (DOA) estimation method based on a dual-branch DOA estimation framework using the PAST (Projection Approximation Subspace Tracking) and EWMA (Exponentially Weighted Moving Average) algorithms. The PAST algorithm is used for adaptive adjustment of compression parameters, avoiding complex covariance matrix estimation and eigenvalue decomposition, thus ensuring the real-time performance of the adaptive adjustment. The EWMA algorithm is used for DOA estimation, assigning different weights to array signals at different times. The exponential decay of these weights rapidly reduces the impact of noise statistical characteristic drift and strong impulse noise, thereby providing a high-quality covariance matrix for DOA estimation and further improving its accuracy.

[0072] 4. The present invention provides a highly robust adaptive direction-of-arrival estimation method. The composite objective function adopts multi-index fusion and penalty constraints to evaluate the non-Gaussianity and non-stationarity of the noise residual after adaptive adjustment of the compression parameters. This provides a comprehensive criterion for whether to update the compression parameters, which can improve the noise statistical characteristics and make them meet the noise stationarity and Gaussianity requirements of the MUSIC algorithm. Attached Figure Description

[0073] Figure 1 This is a flowchart of a method according to an embodiment of the present invention;

[0074] Figure 2 The above are comparison diagrams of the logarithmic compression operator under different values ​​of compression parameters in step 1 of the present invention, where (a) and (b) are comparison diagrams of the logarithmic compression operator under different values ​​of compression intensity parameter and compression region parameter, respectively.

[0075] Figure 3This is a schematic diagram of the array signal obtained in step 2.1 of an embodiment of the present invention;

[0076] Figure 4 Comparison images of directions of arrival obtained by using the traditional MUSIC algorithm, the MUSIC algorithm based on M estimation, the MUSIC algorithm based on FLOM, and the method of this embodiment;

[0077] Figure 5 for Figure 4 Comparison of results from cumulative distribution function analysis;

[0078] Figure 6 This chart compares the success rates of the traditional MUSIC algorithm, the MUSIC algorithm based on M estimation, the MUSIC algorithm based on FLOM, and the method described in this embodiment. Detailed Implementation

[0079] The highly robust adaptive direction-of-arrival estimation method proposed in this invention will be further described in detail below with reference to the accompanying drawings and specific embodiments. Those skilled in the art should understand that these embodiments are merely used to explain the technical principles of this invention and are not intended to limit the scope of protection of this invention.

[0080] A robust adaptive direction-of-arrival estimation method, such as Figure 1 As shown, it includes the following steps:

[0081] Step 1: Construct a logarithmic compression operator for time-varying, non-stationary, and non-Gaussian noise. Then, given the compression candidate parameter set for strong noise, moderate noise, and stationary noise states, as well as the compression parameter update conditions, the range of noise subspace drift, the initial compression parameters of the logarithmic compression operator, the initial signal subspace estimation, the initial covariance matrix, the initial noise subspace projection matrix, and the number of frozen frames for updating the compression parameters, and set the frame counter to zero.

[0082] The noise subspace drift ranges from [0.3, 1.2], the initial signal subspace estimate is a random matrix, the number of frozen frames L for compression parameter updates is 10, the initial covariance matrix R(t-1) = 0.01 × I, and the initial noise subspace projection matrix P n (t-1)=0.01×I, where I is the identity matrix.

[0083] The logarithmic compression operator is:

[0084]

[0085] Where f(x(t)) is the logarithmic compression operator, x(t) is the array signal, t represents the current frame; sgn() is the sign function, α(t) and β(t) are the initial compression strength parameter and the initial compression region parameter, respectively, and ε is an infinitesimal of the logarithmic compression operator.

[0086] like Figure 2 As shown, the logarithmic compression operator exhibits strong sensitivity to small signals. For small signals with amplitudes close to 0, it will be amplified approximately linearly; while for signals with larger amplitudes, the curve slope tends to be gentler, meaning it can effectively compress signals with larger amplitudes. Simultaneously, the compression intensity and region of noise are adjusted by the compression intensity parameter and compression region parameter, respectively. The larger the compression intensity parameter and the smaller the compression region parameter, the greater the amplification intensity for signals with smaller amplitudes, and the more concentrated the effective region. The logarithmic compression operator, through nonlinear stretching of signals with different amplitudes, can effectively suppress strong impulse noise interference under low signal-to-noise ratio conditions, causing the noise in the compressed signal to re-follow a normal distribution. Different values ​​of the compression intensity parameter and compression region parameter represent different compression ranges and compression amplitudes. If fixed compression parameters are used, there may be situations where the compression intensity is too large, suppressing the effective signal, or insufficient, failing to suppress impulse noise interference. To improve robustness under different time-varying environments, this embodiment introduces an adaptive adjustment mechanism for the compression parameters in steps 2-5.

[0087] The conditions for updating compression parameters are:

[0088] If the noise distribution JB test probability value of the candidate parameters in the compressed candidate parameter set and the composite objective function satisfy the following conditions, then the compressed parameters should be updated:

[0089] ,and

[0090] or,

[0091]

[0092] in, , These are the current composite objective function and the composite objective function of the candidate parameters in the compressed candidate parameter set, respectively. , These are candidate parameters from the compression candidate parameter set, namely the compressibility strength parameter and the compression region parameter; , These are the JB test probability values ​​of the current noise distribution and the JB test probability values ​​of the noise distribution of the candidate parameters in the compressed candidate parameter set, respectively, ε. J ε p These are the threshold values ​​for the composite objective function and the JB test probability, respectively.

[0093] The JB test probability value of the noise distribution reflects the Gaussianity of the noise residual. The Gaussianity of the noise residual is a core premise of the orthogonality assumption between the signal and noise subspaces in the MUSIC algorithm. If the Gaussianity of the noise can be significantly improved through the logarithmic compression operator, it should also be considered an effective parameter update. Therefore, in this embodiment, if the candidate parameters in the compression candidate parameter set satisfy the reduction of the composite objective function and the Gaussianity of the noise residual does not significantly degrade, or the Gaussianity of the noise residual is significantly improved, it indicates that its compression effect is better than the initial compression parameters, and the compression parameters are updated.

[0094] In other embodiments, to enable the array signal under strong noise conditions to quickly escape the non-Gaussian environment, when the compression region parameter increases, it satisfies... ,and This allows for compression parameter updates, enabling rapid increases in compression region parameters to enhance the stretching effect on signals with smaller amplitudes. Simultaneously, to ensure effective compression of the array signal, a more stringent judgment is applied to the decrease in compression region parameters under various noise conditions; only when certain conditions are met... ,and Only then will compression parameter updates be allowed.

[0095] The compression candidate parameter sets for strong noise, medium noise, and stationary noise states are given as follows:

[0096] The candidate set of compression parameters for the strong noise state is: the candidate set of compression parameters for the compression intensity parameter remains unchanged, and the compression region parameter is... ;

[0097] The candidate set of compression parameters under moderate noise conditions is: the candidate set of compression parameters with the compression intensity parameter remaining constant and the compression region parameter being... ;

[0098] The candidate sets of compression parameters for stationary noise states are as follows: the candidate sets of compression intensity parameters and compression region parameters are respectively... , ;

[0099] in, Adjust the step size of the compressibility strength parameter for the current steady noise state. , , These are the parameter adjustment step sizes for the current strong noise state, medium noise state, and steady noise state, respectively. .

[0100] In other words, under strong noise conditions, the compression strength parameter remains unchanged, while the compression region parameter remains the same or increases significantly; under moderate noise conditions, the compression strength parameter remains unchanged, while the compression region parameter remains the same or increases slightly; under steady noise conditions, the compression strength parameter and the compression region parameter search for the optimal parameters in both directions.

[0101] Step 2: Acquire the array signal. Based on the initial signal subspace estimation, use the PAST algorithm to process the array signal to obtain the signal subspace estimate. Then, calculate the noise subspace projection matrix based on the signal subspace estimate. Next, calculate the noise subspace drift based on the noise subspace projection matrix and the initial noise subspace projection matrix. If the noise subspace drift is within the range of its variation, proceed to Step 3; otherwise, calculate the instantaneous covariance matrix based on the array signal, and then proceed to Step 6. Specifically:

[0102] Step 2.1, obtain as follows Figure 3 The array signal shown is processed using the PAST algorithm based on the initial signal subspace estimation. Through step-by-step iterative processing, the signal subspace estimation is obtained using the following formula:

[0103]

[0104]

[0105]

[0106] Among them, U s (t) represents the signal subspace estimate, U s (t-1) represents the initial signal subspace estimate; , is the conjugate transpose matrix of the initial signal subspace estimation. The product of the array signal x(t); k(t) is the gain vector. Let be the conjugate transpose of k(t);

[0107] λ is the forgetting factor of the PAST algorithm. Let be the conjugate transpose of y(t), and P(t) and P(t-1) be the inverses of the autocorrelation matrices of y(t) and y(t-1), respectively. When this step is executed for the first time, P(t-1) = 0.01 × I.

[0108] In this embodiment, a 6-element linear array is used to acquire array signals, wherein the effective acoustic signal frequency is 1kHz, the pulse noise is Gaussian pulse noise, the pulse width is 100ms, and the interval is 1.02s.

[0109] The PAST algorithm can transform the problem of tracking and updating the signal subspace into a problem of minimizing the objective function. It uses historical signal subspace estimation to approximate the current projection coefficients, transforming the complex eigenvalue decomposition problem into a linear least squares regression problem. By using the least squares method for recursion, it significantly reduces computational complexity and improves computational efficiency.

[0110] Step 2.2: Based on the signal subspace estimation, and utilizing the orthogonality between the signal and noise subspaces, calculate the noise subspace projection matrix using the following formula:

[0111]

[0112] in, The projection matrix of the noise subspace. It is an N-order identity matrix, where N is the number of array elements used to acquire array signals; For U s The conjugate transpose of (t).

[0113] satisfy and ,in, Represents the noise subspace, i.e., any matrix is... After projection, the signal components will be eliminated, leaving only the noise components.

[0114] Step 2.3: Based on the noise subspace projection matrix and the initial noise subspace projection matrix, calculate the noise subspace drift using the following formula:

[0115]

[0116] in, This represents the drift amount in the noise subspace. Let be the initial noise subspace projection matrix. This represents the Frobenius norm (or F-norm for short).

[0117] Step 2.4: If the noise subspace drift is within the range of noise subspace drift, proceed to step 3; otherwise, calculate the instantaneous covariance matrix based on the array signal, and then proceed to step 6.

[0118] When the statistical characteristics of the noise are stationary, P n The change of P(t) is slow, and the value of δ(t) is small; when the noise changes abruptly or there is strong impulse noise interference, P nThe noise will change significantly, and the value of δ(t) will jump rapidly. The larger the value of δ(t), the more drastic the change in noise characteristics. Therefore, in this embodiment, δ(t) is used as a sensitive noise indicator to achieve real-time tracking of noise dynamics. In this embodiment, the range of noise subspace drift is set to [0.3, 1.2]. When δ(t) < 0.3, the noise is considered to be stable, and no logarithmic compression is performed. When δ(t) > 1.2, the environmental noise is considered to be extremely severe or the evaluation result of this frame has a large error, and no logarithmic compression is performed. When δ(t) ∈ [0.3, 1.2], the statistical characteristics of the noise environment are considered to have changed significantly, the stability of the array signal is disrupted, and there may be impulse interference or noise power jumps. It is necessary to further examine the Gaussianity of the noise fluctuation and update the compression parameters.

[0119] Step 3: Compress the array signal according to the logarithmic compression operator corresponding to the initial compression parameters to obtain the compressed array signal, and calculate the instantaneous covariance matrix based on the compressed array signal; determine whether the number of frames of the array signal has reached the number of frozen frames L for updating the compression parameters. If so, project the compressed array signal to the noise subspace according to the noise subspace projection matrix to obtain the projected signal, and calculate the projection residual; perform JB test on the projection residual to obtain the JB test probability value of the current noise distribution, and calculate the current composite objective function based on it, and then execute step 4; if not, execute step 6.

[0120] The JB test probability value for the current noise distribution is obtained using the following formula:

[0121]

[0122] in, Let N be the JB statistic for the JB test. tot This represents the total number of signals in the array. Given the current skewness, For the current peak value, The cumulative distribution function represents a chi-square distribution with 2 degrees of freedom; The array signal acquired for the nth array element, where n is an integer and 1 ≤ n ≤ N. σ is the sample mean of the array signal. x denoted as the sample standard deviation of the array signal.

[0123] The current composite objective function is calculated using the following formula:

[0124]

[0125] in, Let J0 be the current composite objective function, and J0 be the core objective function for noise compression. Penalty constraints for the objective function;

[0126] J JB This is the JB test normalization term for the noise distribution. ω1, ω2, and ω3 are respectively J JB , , The weights;

[0127] These are candidate parameters from the candidate parameter set for the compressed region parameters. β serves as a reference value for adjusting the parameters of the compression region. max Let λ1 and λ2 be the maximum values ​​in the set of candidate compression parameters for the compression region, respectively. , The constraint strength, An indicator function for adjusting the step size of the current compression region parameters. , These are the compression region parameters for the next frame of the array signal. Adjust the step size for the parameters of the current compression region.

[0128] First, for the time-varying logarithmic Gaussian noise contained in the compressed array signal obtained by logarithmic compression, if only the JB test probability value is used as the objective function, the JB test probability value is prone to jumps due to instantaneous impulse interference from the noise, thus leading to misjudgments in parameter adjustment. Second, for the MUSIC algorithm, the time-varying logarithmic Gaussian noise has circular symmetry, which directly determines the orthogonality of the signal subspace and the noise subspace and whether the noise satisfies the algorithm assumptions. Moreover, the noise is non-stationary, which can easily lead to an imbalance in the covariance of the real and imaginary parts, violating the orthogonality assumption between the signal subspace and the noise subspace. Finally, in order to balance the optimization real-time performance and convergence stability of adaptive compression parameter adjustment, it is necessary to avoid the compression parameters from oscillating due to the instantaneous high-intensity fluctuations of impulse noise, or getting trapped in local optima, which would lead to a degradation in the compression effect. Therefore, in this embodiment, the composite objective function includes multi-index fusion and penalty constraints, which can comprehensively measure the non-Gaussianity and non-stationarity of the noise residual in the compressed array signal. During the adjustment of compression parameters, the greater the deviation of the compression region parameters from their conditional reference values, the greater the penalty, thereby rapidly increasing the composite objective function value when the core objective function of noise compression is small. At the same time, asymmetric penalties are applied to the rise and fall of compression region parameters. When the compression region parameters fall, a penalty of 1.5 times is applied, and when they rise, a penalty of 1 times is applied, to avoid the reduction in compression effect caused by the blind fall of compression region parameters.

[0129] Step 4: Determine the current noise state based on the JB test probability value of the current noise distribution, and select the corresponding compression candidate parameter set from the compression candidate parameter sets given in Step 1 for strong noise state, medium noise state and stationary noise state based on the current noise state. Then, according to the logarithmic compression operator corresponding to the candidate parameter in the compression candidate parameter set, obtain the noise distribution JB test probability value and composite objective function of each candidate parameter in the compression candidate parameter set according to the method in Step 3.

[0130] The specific method for determining the current noise state is as follows:

[0131] If the current noise distribution has a JB test probability value p JB <0.01, or the JB test probability value of the current noise distribution p JB If the value is less than 0.01 and the noise subspace drift δ(t) is greater than 1.1, then it is a strong noise state;

[0132] If the JB test probability value of the current noise distribution is 0.01 ≤ p JB <0.05, or the JB test probability value of the current noise distribution is 0.01 ≤ p JB If the noise level is less than 0.05 and the noise subspace drift is 1.1 ≥ δ(t) > 0.9, then it is considered a medium noise state.

[0133] If the current noise distribution has a JB test probability value p JB ≥0.05, or the JB test probability value p of the current noise distribution. JB If the noise level is ≥0.05 and the noise subspace drift δ(t) ≤0.9, then it is a stationary noise state.

[0134] Step 5: Based on the JB test probability value of the noise distribution and the composite objective function of each candidate parameter in the compressed candidate parameter set, as well as the current JB test probability value of the noise distribution and the current composite objective function, determine whether the candidate parameters in the compressed candidate parameter set meet the compression parameter update conditions. If the compression parameter update conditions are met, take the candidate parameter corresponding to the minimum value of the composite objective function in the compressed candidate parameter set as the initial compression parameter, and update the adjustment step size of the compressed candidate parameter set in the current noise state as follows: Set the frame counter to zero, then execute step 6; if the compression parameter update condition is not met, update the adjustment step size of the current noise state's compression candidate parameter set to: Then, increment the frame counter and proceed to step 6. , These are the adjustment step sizes for the compression intensity parameter and compression region parameter of the next frame, respectively. Adjust the step size for the current compressive strength parameter. , These represent the maximum adjustment step sizes for the candidate parameter sets of the compressibility strength parameter and the compression region parameter, respectively. , These are the minimum adjustment step sizes for the candidate parameter sets of the compressibility strength parameter and the compression region parameter, respectively. In this embodiment, It is 0.08. It is 0.015.

[0135] Step 6: Based on the instantaneous covariance matrix and the initial covariance matrix, the covariance matrix is ​​obtained recursively using the EWMA algorithm. Then, eigenvalue decomposition and MUSIC spatial spectrum peak search are performed on the covariance matrix to obtain and output the direction of arrival.

[0136] The covariance matrix is ​​calculated using the following formula:

[0137]

[0138] in, Let covariance matrix be the variance matrix. Let λ be the initial covariance matrix. R The forgetting factor of the covariance matrix. Let be the instantaneous covariance matrix.

[0139] While the PAST algorithm avoids complex covariance estimation and eigenvalue decomposition, directly updating the noise subspace projection matrix in real time through gradient descent iterations to meet the real-time requirements of adaptive algorithms, a high-quality covariance matrix is ​​still needed for the high-precision requirements of direction-of-arrival (DOA) estimation. Therefore, this embodiment, based on the PAST algorithm, employs the EWMA algorithm to provide a high-quality covariance matrix for the MUSIC algorithm, adapting to the time-varying log-Gaussian noise contained in the compressed array signal and meeting the engineering application requirements of adaptive DOA estimation. The EWMA algorithm assigns different weights to the array signal at different times, rapidly reducing the impact of noise statistical characteristic drift and strong impulse noise through exponential decay of the weights.

[0140] The formula for calculating the covariance matrix is ​​rewritten to obtain:

[0141]

[0142] Where R(1) is the covariance matrix of the first frame array signal, R x (ta) is the instantaneous covariance matrix of the array signal in the ta-th frame, where a is an integer and 0≤a≤t-2.

[0143] It is evident that, within the EWMA algorithm framework, the weight of historical data decays exponentially with the increase of array signal frames. Even when the current array signal is subjected to strong impulse noise or abrupt changes in noise statistical characteristics, causing a drastic jump in the covariance matrix, its impact can be rapidly attenuated. Furthermore, due to the accumulation of long-term statistical characteristics of noise in historical data, through λ... R The attenuation weights retain this trend and can also effectively prevent the covariance matrix from deviating from the true statistical characteristics due to the current array signal anomalies, thereby improving the reliability of covariance matrix estimation under time-varying non-Gaussian noise and providing high-quality input for eigenvalue decomposition and spatial spectrum peak search of the MUSIC algorithm.

[0144] Step 7: Use the signal subspace estimate as the initial signal subspace estimate, the noise subspace projection matrix as the initial noise subspace projection matrix, and the covariance matrix as the initial covariance matrix. Then return to step 2 to obtain the next frame of array signal until the direction of arrival estimation is completed.

[0145] The algorithm employs traditional MUSIC algorithms, MUSIC algorithms based on M estimation, MUSIC algorithms based on FLOM (Fractional Order Moments), and this embodiment. Figure 3 The direction of arrival (DOA) of the array signal is estimated, and the obtained DOA is as follows: Figure 4 , Figure 5 As shown, the direction of arrival obtained using the method in this embodiment has the best effect, with a 90th percentile < 3.63°. Figure 6 As shown, for Figure 3 Success rate analysis was performed on the directions of arrival obtained by various algorithms. The success rate of the method in this embodiment was 0.7405, which is the best. The success rate is defined as the absolute error between the estimated direction of arrival and the true direction of arrival < 2°.

Claims

1. A highly robust adaptive direction-of-arrival estimation method, characterized in that, Includes the following steps: Step 1: Construct a logarithmic compression operator for time-varying, non-stationary, and non-Gaussian noise; Then, given the compression candidate parameter sets for strong noise state, medium noise state and stationary noise state, as well as the compression parameter update conditions, the range of noise subspace drift, the initial compression parameters of the logarithmic compression operator, the initial signal subspace estimate, the initial covariance matrix and the initial noise subspace projection matrix; Step 2: Obtain the array signal. Based on the initial signal subspace estimation, use the PAST algorithm to process the array signal to obtain the signal subspace estimation. Then, calculate the noise subspace projection matrix based on the signal subspace estimation. Then, the noise subspace drift is calculated based on the noise subspace projection matrix and the initial noise subspace projection matrix. If the noise subspace drift is within the range of noise subspace drift, then proceed to step 3; Otherwise, calculate the instantaneous covariance matrix based on the array signal, and then proceed to step 6; Step 3: Compress the array signal according to the logarithmic compression operator corresponding to the initial compression parameters to obtain the compressed array signal, and calculate the instantaneous covariance matrix based on the compressed array signal; then, project the compressed array signal onto the noise subspace according to the noise subspace projection matrix to obtain the projected signal, and calculate the projection residual. Perform a JB test on the projected residuals to obtain the JB test probability value of the current noise distribution, and calculate the current composite objective function based on it; Step 4: Based on the JB test probability value of the current noise distribution, determine the current noise state, and select the corresponding compression candidate parameter set from the compression candidate parameter set given in Step 1 for strong noise state, medium noise state and stationary noise state according to the current noise state. Then, according to the logarithmic compression operator corresponding to the candidate parameter in the compression candidate parameter set, obtain the noise distribution JB test probability value and composite objective function of each candidate parameter in the compression candidate parameter set according to the method in Step 3. Step 5: Based on the JB test probability value of the noise distribution and the composite objective function of each candidate parameter in the compressed candidate parameter set, as well as the current JB test probability value of the noise distribution and the current composite objective function, determine whether the candidate parameters in the compressed candidate parameter set meet the compression parameter update conditions. If the compression parameter update conditions are met, take the candidate parameter corresponding to the minimum value of the composite objective function in the compressed candidate parameter set as the initial compression parameter, and then execute Step 6; if the compression parameter update conditions are not met, then execute Step 6. Step 6: Based on the instantaneous covariance matrix and the initial covariance matrix, the covariance matrix is ​​obtained recursively using the EWMA algorithm. Then, the covariance matrix is ​​subjected to eigenvalue decomposition and MUSIC spatial spectrum peak search to obtain and output the direction of arrival. Step 7: Use the signal subspace estimate as the initial signal subspace estimate, the noise subspace projection matrix as the initial noise subspace projection matrix, and the covariance matrix as the initial covariance matrix. Then return to step 2 to obtain the next frame of array signal until the direction of arrival estimation is completed.

2. The highly robust adaptive direction-of-arrival estimation method according to claim 1, characterized in that, In step 1, the logarithmic compression operator is: ； Where f(x(t)) is the logarithmic compression operator, x(t) is the array signal, t represents the current frame; sgn() is the sign function, α(t) and β(t) are the initial compression strength parameter and the initial compression region parameter, respectively, and ε is an infinitesimal quantity of the logarithmic compression operator; The conditions for updating the compression parameters are as follows: If the noise distribution JB test probability value of the candidate parameters in the compressed candidate parameter set and the composite objective function satisfy the following conditions, then the compressed parameters should be updated: ,and ; or, ； in, , These are the current composite objective function and the composite objective function of the candidate parameters in the compressed candidate parameter set, respectively. , These are candidate parameters from the compression candidate parameter set, namely the compressibility strength parameter and the compression region parameter; , These are the JB test probability values ​​of the current noise distribution and the JB test probability values ​​of the noise distribution of the candidate parameters in the compressed candidate parameter set, respectively, ε. J ε p These are the threshold values ​​for the composite objective function and the probability value of the JB test, respectively. The compression candidate parameter sets for the strong noise state, the medium noise state, and the stationary noise state are given as follows: The candidate set of compression parameters under strong noise conditions is: the candidate set of compression parameters for which the compression intensity parameter remains constant and the compression region parameter is... ; The candidate set of compression parameters under moderate noise conditions is: the candidate set of compression parameters with the compression intensity parameter remaining constant and the compression region parameter being: ; The candidate sets of compression parameters for stationary noise states are as follows: the candidate sets of compression intensity parameters and compression region parameters are respectively... , ; in, Adjust the step size of the compressibility strength parameter for the current steady noise state. , , These are the parameter adjustment step sizes for the current strong noise state, medium noise state, and steady noise state, respectively. .

3. The highly robust adaptive direction-of-arrival estimation method according to claim 2, characterized in that, In step 4, the specific method for determining the current noise state is as follows: If the current noise distribution has a JB test probability value p JB <0.01, or the JB test probability value of the current noise distribution p JB If the value is less than 0.01 and the noise subspace drift δ(t) is greater than 1.1, then it is a strong noise state; If the JB test probability value of the current noise distribution is 0.01 ≤ p JB <0.05, or the JB test probability value of the current noise distribution is 0.01 ≤ p JB If the noise level is less than 0.05 and the noise subspace drift is 1.1 ≥ δ(t) > 0.9, then it is considered a medium noise state. If the current noise distribution has a JB test probability value p JB ≥0.05, or the JB test probability value p of the current noise distribution. JB If the noise level is ≥0.05 and the noise subspace drift δ(t) ≤0.9, then it is a stationary noise state.

4. The highly robust adaptive direction-of-arrival estimation method according to claim 3, characterized in that, In step 3, the current composite objective function is calculated using the following formula: ； Where J0 is the core objective function for noise compression. Penalty constraints for the objective function; J JB This is the JB test normalization term for the noise distribution. ; , These represent the current skewness and current kurtosis, respectively, and ω1, ω2, and ω3 represent J... JB , , The weights; These are candidate parameters from the set of candidate compression parameters for the compressed region parameters. β serves as a reference value for adjusting the parameters of the compression region. max Let λ1 and λ2 be the maximum values ​​in the set of candidate compression parameters for the compression region, respectively. , The constraint strength, An indicator function for adjusting the step size of the current compression region parameters. , These are the compression region parameters for the next frame of the array signal. Adjust the step size for the parameters of the current compression region.

5. The highly robust adaptive direction-of-arrival estimation method according to claim 4, characterized in that, Step 2 is as follows: Step 2.1: Obtain the array signal. Based on the initial signal subspace estimation, use the PAST algorithm to process the array signal. Through step-by-step iterative processing, the signal subspace estimation is obtained using the following formula: ； ； ； Among them, U s (t) represents the signal subspace estimate, U s (t-1) represents the initial signal subspace estimate; , is the conjugate transpose matrix of the initial signal subspace estimation. The product of the array signal x(t); k(t) is the gain vector. Let be the conjugate transpose of k(t); λ is the forgetting factor of the PAST algorithm. Let be the conjugate transpose of y(t), and P(t) and P(t-1) be the inverses of the autocorrelation matrices of y(t) and y(t-1), respectively. Step 2.2: Based on the signal subspace estimation, and utilizing the orthogonality between the signal and noise subspaces, calculate the noise subspace projection matrix using the following formula: ； Among them, P n (t) is the projection matrix of the noise subspace, I N Let N be an N-order identity matrix, where N is the number of array elements used to acquire the array signal. For U s The conjugate transpose of (t); Step 2.3: Based on the noise subspace projection matrix and the initial noise subspace projection matrix, calculate the noise subspace drift using the following formula: ； Where δ(t) is the noise subspace drift, P n (t-1) is the initial noise subspace projection matrix. It is the Frobenius norm; Step 2.4: If the noise subspace drift is within the range of noise subspace drift, proceed to step 3; otherwise, calculate the instantaneous covariance matrix based on the array signal, and then proceed to step 6.

6. The highly robust adaptive direction-of-arrival estimation method according to claim 5, characterized in that, In step 6, the covariance matrix is ​​calculated using the following formula: ； Where R(t) is the covariance matrix, R(t-1) is the initial covariance matrix, and λ R R is the forgetting factor of the covariance matrix. x (t) is the instantaneous covariance matrix.

7. The highly robust adaptive direction-of-arrival estimation method according to claim 6, characterized in that, In step 3, the JB test probability value of the current noise distribution is obtained by the following formula: ； Where JB(t) is the JB statistic of the JB test, and N tot This represents the total number of signals in the array. Let represent the cumulative distribution function of a chi-square distribution with 2 degrees of freedom. The array signal acquired for the nth array element, where n is an integer and 0 ≤ n ≤ N. σ is the sample mean of the array signal. x denoted as the sample standard deviation of the array signal.

8. The highly robust adaptive direction-of-arrival estimation method according to claim 7, characterized in that, Step 5 also includes: If the compression parameter update condition is met, then the adjustment step size of the current noise state's candidate compression parameter set is updated as follows: ； If the compression parameter update condition is not met, then the adjustment step size of the current noise state's compression candidate parameter set is updated as follows: ； in, , These are the adjustment step sizes for the compression intensity parameter and compression region parameter of the next frame, respectively. Adjust the step size for the current compressive strength parameter. , These represent the maximum adjustment step sizes for the candidate parameter sets of the compressibility strength parameter and the compression region parameter, respectively. , These are the minimum values ​​of the adjustment step size for the candidate parameter sets of the compressibility strength parameter and the compression region parameter, respectively.

9. The highly robust adaptive direction-of-arrival estimation method according to claim 8, characterized in that, Step 1 also includes setting the number of frozen frames for updating compression parameters; Step 3 specifically involves: The array signal is compressed using the logarithmic compression operator corresponding to the initial compression parameters to obtain a compressed array signal. The instantaneous covariance matrix is ​​then calculated based on the compressed array signal. It is determined whether the number of frames of the array signal has reached the number of frozen frames for updating the compression parameters. If so, the compressed array signal is projected onto the noise subspace according to the noise subspace projection matrix to obtain the projected signal, and the projection residual is calculated. The JB test is performed on the projection residual to obtain the JB test probability value of the current noise distribution, and the current composite objective function is calculated based on it. Then, step 4 is executed. If not, proceed to step 6.

10. A highly robust adaptive direction-of-arrival estimation method according to claim 9, characterized in that, In step 1, the range of the noise subspace drift is [0.3, 1.2], the initial signal subspace estimate is a random matrix, the initial covariance matrix R(t-1) = 0.01 × I, and the initial noise subspace projection matrix P n (t-1)=0.01×I, where I is the identity matrix; In step 5, the The value is 0.08-0.

2. It is 0.01-0.1.

Images