Noise Suppression for Direction of Arrival Estimation in Co-located MIMO Sonar

Noise suppression capacity in multiple-input multiple-output (MIMO) sonar signal processing is derived under the assumption of white Gaussian noise. However, underwater noise mainly includes white Gaussian noise and colored noise. There exists a certain correlation between the noise signals received by each MIMO sonar array element. The performance of traditional direction-of-arrival (DOA) estimation methods decreases obviously in complex marine noise. In this paper, we propose a marine environment noise suppression method for MIMO applied to multiple targets’ DOA estimation. The noise field can be decomposed into a symmetric noise component and an asymmetric noise component. We use the covariance matrix imaginary component to pre-estimate the signal sources, then use the dimension reduction transformation to reconstruct the real component of the covariance matrix. The Toeplitz technique is utilized to reduce the correlation of the reconstructed covariance matrix. Thus, the subspace decomposition-based techniques such as multiple signal classification (MUSIC) can be used for multiple targets’ DOA estimation. To reduce the computational complexity of the methods, search-free direction-finding techniques such as the estimation of signal parameters via rotational invariance techniques (ESPRIT) can be utilized. As a result, the proposed methods can achieve better direction-finding performance in the condition of limited snapshots with lower computational cost. The corresponding Cramer-Rao bound (CRB) is deduced and the signal-to-noise ratio (SNR) gain obtained by dimension reduction processing is discussed. Simulation results also show the superiority of the proposed method over the existing methods.


Introduction
Multiple-input multiple-output (MIMO) sonar is a system that consists of the transmit array and the receive array. Both arrays employ multiple sensors. The transmit array emits orthogonal waveforms and the receive array completes the echo signal acquisition [1][2][3]. The MIMO array can be mainly classified into two types. One type is equipped with separated sensors [4][5][6][7][8]. This kind of MIMO sonar system utilizes widely separated sensors to gain spatial diversity. The large aperture arrays at both the transmitter and receiver enable the MIMO sonar to view different aspects of a target [9]. The other type is equipped with co-located sensors [10][11][12]. This kind of MIMO sonar employs arrays of closely spaced sensors to view the same side of one target from the same angle. In this scenario, it is assumed that the source is the point target in the far field of MIMO sonar. Additionally, the transmitted signals are normally assumed to be narrowband. Hassanien et al. [13] assumed the point target signal model. Under this circumstance, since the receiver gathers the multiple independent waveforms, the diversity of the waveform is used to increase the receiver virtual aperture. sensitive to snapshots. It is worth mentioning that the probability of the target resolution of this method increases faster than traditional methods. Simultaneously, this method avoids the freedom loss caused by decoherent processing. We can simulate complex marine environments by adding colored noise to white Gaussian noise. In the actual noise field, the proportion of the symmetric noise component is much larger than that of the asymmetric noise component [38,39]. Considering that the imaginary part of the covariance matrix has no relationship with the symmetric noise component, we can remove the real part of the covariance matrix to avoid the effect of the symmetric noise component on DOA estimation. The real part of the covariance matrix is reconstructed by using the method of dimension reduction transformation and the replacement principle of the matrix imaginary component, and the interference of the bilateral spectrum interference is avoided. Thus, we obtain the reconstructed covariance matrix. Moreover, Toeplitz is utilized to reduce the correlation of the covariance matrix. This operation can reduce the influence of limited snapshots and complex marine noise on the covariance matrix. The new covariance matrix we obtained is closer to the ideal covariance matrix. Therefore, the phenomenon of fuzzy division between the signal subspace and the noise subspace can be avoided. To further reduce the computational complexity, dimension reduction processing is once again used to obtain a lower dimension covariance matrix. Above all, we get a novel noise suppression method for MIMO sonar DOA estimation. This method is based on dimension reduction transformation and the Toeplitz decoherence technique. For the sake of simplicity, we abbreviate this method as RC-STIM. Then, DOA estimation methods, such as MUSIC can be used to pre-estimate the target angle and carry out the final DOA estimation. Due to the established rotational invariance property at the MIMO sonar array, ESPRIT can be used to replace MUSIC. This search-free method can obviously reduce the computational complexity. In view of the advantageous property of MUSIC, we try to combine the advantages of the two methods to reach a balance between performance and complexity. As a result, the proposed algorithms can achieve better DOA estimation performance at lower computational cost, in limited snapshots, simultaneously.
The paper is organized as follows. In Section 2, the signal model of MIMO sonar is briefly introduced. Section 3 presents the influence of the noise component and the noise suppression MIMO sonar model. Several performance parameters, such as computational complexity, Cramer-Rao Bound and SNR gain are presented in Section 4. The simulation results which show the superiority of the proposed MUSIC-based and ESPRIT-based RC-STIM MIMO sonar DOA estimation techniques are presented in Section 5, followed by conclusions drawn in Section 6.

MIMO Sonar Signal Model
Consider a uniform linear array (ULA) of co-located MIMO sonar, equipped with a transmit array of M sensors and a receive array of N sensors, with half-wavelength spacing employed for both the transmit and receive arrays. The two arrays are collinear and the genetic centers are overlapped. The sensors are assumed to be in close proximity so that all the elements can detect a far-field target at the same spatial angle. Assume that there are K targets, satisfying the condition 2K < MN. A T = [a t (θ 1 ), a t (θ 2 ), . . . , a t (θ K )] is a M × K dimensional matrix composed of the K transmit array steering vector. A R (θ) = [a r (θ 1 ), a r (θ 2 ), . . . , a r (θ K )] is a N × K dimensional matrix composed of the K receive array steering vector. The effectiveness of the proposed method is not influenced by whether the transmit and receive array element numbers are equal or not. For simplicity, but without loss of generality, take the co-located MIMO array with equal numbers of transmit and receive array elements as an example. The convolution operation is performed on the positions of the transmit and receive array elements [40], we can obtain the positions of the MIMO sonar virtual array elements, which are shown in Figure 1. Under the condition of far-field, the coordinates of each virtual array element are equal to the convolution of the corresponding transmit and receive array elements. The coordinates of the virtual array elements can be expressed as VP i = (m + n − 1)d t , m = 1, . . . , M, n = 1, . . . , N. In this scenario, several overlapped positions of virtual elements will appear due to the co-located location of the transmit and receive array elements. Each column in Figure 1 represents a set of virtual elements in the same location. the transmit and receive array elements [40], we can obtain the positions of the MIMO sonar virtual array elements, which are shown in Figure 1. Under the condition of far-field, the coordinates of each virtual array element are equal to the convolution of the corresponding transmit and receive array elements. The coordinates of the virtual array elements can be expressed as In this scenario, several overlapped positions of virtual elements will appear due to the co-located location of the transmit and receive array elements. Each column in Figure 1 represents a set of virtual elements in the same location. The M transmit elements are used to transmit M orthogonal waveforms. Therefore, by discretization of P snapshots, the representation of the transmitted signals is a unit energy baseband pulse matrix. We assume that there is no Doppler shift between the targets and each element. Simultaneously, the effects of channel fluctuation, dielectric absorption, echo distortion and propagation loss on the echo signals are ignored. Then, we operate under the assumption that there is no multipath among the different sources' emissions, i.e., the propagation is nondispersive. The narrowband echo signals of K far-field targets obtained at the receiver can be modeled as where k β is the reflection coefficient of the target located at an unknown spatial angle k θ . Spatial noise N , consisting of white Gaussian noise and colored noise, is a N P × dimensional matrix.
T ( ) stands for transpose. The distribution model of array elements and target signals is plotted in Figure 2a. As can be seen from this figure, the vertical direction of the line array is the reference direction 0°. The angle increases with rotating to the right until 90°.
(a) The M transmit elements are used to transmit M orthogonal waveforms. Therefore, by discretization of P snapshots, the representation of the transmitted signals is a unit energy baseband pulse matrix. We assume that there is no Doppler shift between the targets and each element. Simultaneously, the effects of channel fluctuation, dielectric absorption, echo distortion and propagation loss on the echo signals are ignored. Then, we operate under the assumption that there is no multipath among the different sources' emissions, i.e., the propagation is nondispersive. The narrowband echo signals of K far-field targets obtained at the receiver can be modeled as where β k is the reflection coefficient of the target located at an unknown spatial angle θ k . Spatial noise N, consisting of white Gaussian noise and colored noise, is a N × P dimensional matrix. ( ) T stands for transpose. The distribution model of array elements and target signals is plotted in Figure 2a. As can be seen from this figure, the vertical direction of the line array is the reference direction 0 • . The angle increases with rotating to the right until 90 • .
Sensors 2019, 19, x 4 of 21 the transmit and receive array elements [40], we can obtain the positions of the MIMO sonar virtual array elements, which are shown in Figure 1. Under the condition of far-field, the coordinates of each virtual array element are equal to the convolution of the corresponding transmit and receive array elements. The coordinates of the virtual array elements can be expressed as In this scenario, several overlapped positions of virtual elements will appear due to the co-located location of the transmit and receive array elements. Each column in Figure 1 represents a set of virtual elements in the same location. The M transmit elements are used to transmit M orthogonal waveforms. Therefore, by discretization of P snapshots, the representation of the transmitted signals is a unit energy baseband pulse matrix. We assume that there is no Doppler shift between the targets and each element. Simultaneously, the effects of channel fluctuation, dielectric absorption, echo distortion and propagation loss on the echo signals are ignored. Then, we operate under the assumption that there is no multipath among the different sources' emissions, i.e., the propagation is nondispersive. The narrowband echo signals of K far-field targets obtained at the receiver can be modeled as where k β is the reflection coefficient of the target located at an unknown spatial angle k θ . Spatial noise N , consisting of white Gaussian noise and colored noise, is a N P × dimensional matrix. The components of spatial noise include white Gaussian noise and colored noise. It can be regarded as the superposition of discrete planar waves generated by several noise sources. Compared to actual noise, the error decreases as the number of noise sources increases. According to statistical characteristics, if we add as many noise sources as possible, we can gain an environmental noise field that is closer to the true value. As shown in Figure 2b, spatial noise can be decomposed into several narrowband sources. φ denotes the azimuth of the noise source. Assume that the spatial noise consists of Q narrowband noise sources with a frequency of f . The power and azimuth of the th q noise source is 2 q σ and q φ , respectively. Obviously, we can get The noise waveform is represented by ( ) q n t . Therefore, we can obtain the noise waveform received by the th m array element: . Therefore, we can define 1 ′ N and 2 N as symmetric noise and define Δ N as asymmetric noise. Further, we can draw the conclusion that random noise can be decomposed into a symmetric noise component and an asymmetric noise component. Figure 3 shows the power distribution and decomposition of the noise component. Therefore, subfigure (a) depicts the power of the noise component, it contains colored noise and white Gaussian noise term. Subfigure (b) shows the symmetric noise component. It can be seen from subfigure (b) that the azimuth and energy of noise are statistically symmetric around the receive array center. Obviously, we can draw that symmetrical noise component account for a much larger proportion than the asymmetric noise component; this may lead to a direct way of suppressing noise. The components of spatial noise include white Gaussian noise and colored noise. It can be regarded as the superposition of discrete planar waves generated by several noise sources. Compared to actual noise, the error decreases as the number of noise sources increases. According to statistical characteristics, if we add as many noise sources as possible, we can gain an environmental noise field that is closer to the true value. As shown in Figure 2b, spatial noise can be decomposed into several narrowband sources. φ denotes the azimuth of the noise source. Assume that the spatial noise consists of Q narrowband noise sources with a frequency of f . The power and azimuth of the qth noise source is σ 2 q and φ q , respectively. Obviously, we can get φ q = φ Q−q+1 . The noise waveform is represented by n q (t). Therefore, we can obtain the noise waveform received by the mth array element: Q ∑ q=1 n q (t)e −jk T n (φ q )P m . Where k n (φ q ) = 2π f v n (φ q )/c, v n (φ q ) denotes the direction vector of the noise sources, c is the sound velocity, P m is the position vector of each element. For simplicity, we take two symmetric noise sources, N 1 and N 2 as examples. The azimuths are φ and −φ. The power is σ 2 N 1 and σ 2 N 2 , σ 2 N 1 > σ 2 N 2 . We can get σ 2 , i.e., N 1 can be divided into N 1 with the power of σ 2 N 2 and N ∆ with the power of ∆ σ 2 N 1 . Therefore, we can define N 1 and N 2 as symmetric noise and define N ∆ as asymmetric noise. Further, we can draw the conclusion that random noise can be decomposed into a symmetric noise component and an asymmetric noise component. Figure 3 shows the power distribution and decomposition of the noise component. Therefore, subfigure (a) depicts the power of the noise component, it contains colored noise and white Gaussian noise term. Subfigure (b) shows the symmetric noise component. It can be seen from subfigure (b) that the azimuth and energy of noise are statistically symmetric around the receive array center. Obviously, we can draw that symmetrical noise component account for a much larger proportion than the asymmetric noise component; this may lead to a direct way of suppressing noise. By using matched filtering on echo signals generated by the orthogonal signals S , we can obtain the autocorrelation function at the output end.
Because of the orthogonality between the echo signals, the covariance matrix of the transmitting signal can be simplified into the identity matrix M M ×

I
. By substituting the matrix into Equation (2), the receive signal matrix can be reformulated as By using matched filtering on echo signals generated by the orthogonal signals S, we can obtain the autocorrelation function at the output end.
Because of the orthogonality between the echo signals, the covariance matrix of the transmitting signal can be simplified into the identity matrix I M×M . By substituting the matrix into Equation (2), the receive signal matrix can be reformulated as By vectorizing Equation (3), we can obtain the output end signal sample vector after matched filtering.
obeys the spatial colored noise distribution. n 2 = vec(N 2 S H ) obeys the complex Gaussian distribution with zero mean and covariance matrix σ 2 (4) can be expressed as The receiving signal matrix is composed of L snapshots.
By calculating the echo signal covariance matrix of Equation (7), we can obtain where Q is the covariance matrix of the spatial noise component, which contains spatial colored noise and white Gaussian noise. The uncorrelation between the target reflection coefficients enables simplifying the matrix into a diagonal matrix.

Problem Formulation
The main goal is to suppress the influence of noise in the direction of arrival estimation, especially the symmetric component of the spatial noise. In an actual acoustic environment, the symmetric noise component in the noise component has a large influence on the real part of the covariance matrix and has little effect on the imaginary part. By using the imaginary part DOA estimation method, the real part of the echo signal covariance matrix is removed, and the influence of symmetric noise on the direction-finding performance of the algorithm is reduced. To further analyze the covariance matrix, we consider the idea of dimension reduction transformation and obtain Sensors 2019, 19, 1325 It can be rewritten as a matrix The signal covariance matrix R s can be modeled as The signal covariance matrix can be expressed as (R s ) a,b = G( K ∑ k=1 β 2 k e jπ(b−a) sin θ k )G H ; the imaginary part of the signal covariance matrix can be obtained as Considering Equation (15), the two parts of the equation can be regarded as the covariance matrix of the signal matrix with the incident angle of θ k and −θ k . It is worth noting that the imaginary part of the signal covariance matrix contains the correct target direction-finding information and the false target direction-finding information. The two sets of azimuths are symmetric about the normal direction of the array.
The noise component can be considered as the superposition of mutually independent plane waves generated by several noise sources. The larger the number of noise sources, the closer the noise is to the real noise. In this paper, we decompose the spatial noise component into several narrowband noises with the quantity of P and the center frequency of f . We can rewrite the noise as n = P ∑ p=1 a n (ϕ p )n p (t), (16) To simplify the calculation, in Figure 2, we assume that there are only two receive array elements and two symmetric noise sources, each of which is independent with each other and has the same power δ 2 . The azimuth angle of the two noise sources relative to the array normal is θ and −θ. The output signals with matched filtering of the two receive array elements can be respectively modeled as The spatial correlation function of X 1 (t) and X 2 (t) can be expressed as where (19) can be simplified as It is worth noting that the imaginary part of γ 12 is zero. Therefore, the covariance matrix of the symmetric noise component is a real symmetric matrix. Moreover, the imaginary part of the covariance matrix is independent of the symmetric noise components in the noise component. From Equation (20), we observe that the interference of symmetric noise components in the noise component can be suppressed by eliminating the real part of the covariance matrix.
Let the real part of the covariance matrix R YY is R R , and the imaginary part is R I . Suppose that Λ = [A, A c ]. We can define a new matrix where ( ) c denotes conjugate. As [A, A c ] and diag(R BB , −R BB ) are all full column rank matrices and the rank is 2K. Therefore, the rank of ∆R is 2K. The power spectrum obtained from the imaginary part contains 2K peaks, where K pseudo-peaks are symmetric with the real peaks. The angles of the peaks are expressed as Ψ = [ψ 1 , ψ 2 , . . . , ψ K , −ψ 1 , −ψ 2 , . . . , −ψ K ], the real angle estimation is Among them, the first K elements are the real target locations. Suppose that the matrix I is the permutation matrix generated by noise, the locations of each element in Θ are changed by one transformation and we can get the following relationship Ψ = ΘI. By constructing the array manifold matrix of the imaginary part by A c By substituting Equation (21) into Equation (22), we can get when A I = Λ, we can get I = E. The diagonal elements in R are sorted correctly. In this scenario, the first K diagonal elements are positive and the corresponding K target reflected signal power. When A I = Λ, the diagonal sort is not correct in R. The matrix is still a diagonal matrix after two transformations, but the positions of the elements on the diagonal change. The positive signal power estimated value corresponds to the correct signal angle. By confirming the correct signal angle estimation and signal power estimation according to the element positions of the diagonal in Equation (23), we extract the positive value of the diagonal in R as the signal power σ 2 k (k = 1, 2, · · · , K). The angle information in Ψ corresponding to the position is the target direction of the estimation. Construct the estimated power matrix containing the signal as ∆ s = diag σ 2 1 , σ 2 2 , . . . , σ 2 K . By considering the idea of dimension reduction transformation, we can obtain T , which can be rewritten as a matrix We reconstruct the covariance matrix aŝ By taking the real part of D∆ s D H , we can obtain . Therefore, the covariance matrix of the echo signal with noise suppression can be reconstructed as Because of the numerical difference and complex noise in the array elements, the covariance matrix is an approximate Toeplitz matrix. In this case, when we do eigen-decomposition on the covariance matrix, the division of the signal subspace and noise subspace is not clear, which will lead to a performance decline in the DOA estimation of MUSIC algorithm. This phenomenon is particularly pronounced in the case of low SNR, the direction-finding result shows worse performance compared with the traditional MUSIC method. Therefore, we can modify the real part R R and imaginary part R I of the reconstructed covariance matrix R by Toeplitz respectively.
By averaging the diagonal elements of the matrix R ∆ , we can obtain the covariance matrix of the ideal situation. Symbol m and n denote the rows and columns of M + N − 1 dimension matrix, r mn denotes the element of the mth row and the nth column. The modified matrix elements can be expressed asr We can obtain the modified R R by Equation (25). In the same way, by averaging the diagonal elements of matrix R I , we can obtain the modified R I . However, it is worth noting that, since the MIMO sonar system normally has a large virtual array number, the computational complexity of the conventional subspace-based method grows up greatly. We use the reduced dimension transformation method to improve computational efficiency. Based on Equation (12), we define W = G H G, and get W = diag(1, 2, · · · , min(M, N), · · · , min(M, N) , · · · , 2, 1). We denote P = W − 1 2 G H by the reduced dimension transformation method and then use the transformation P for the echo signal after matched filtering. We can obtain the dimension reduction matrix z(t) = Py(t). Moreover, the covariance matrix R z can be expressed as Through the above transformation, we obtain the covariance matrix R z , which efficiently reduces the computational complexity. It is worth pointing that since the application of the spatial spectrum estimation algorithm is in the DOA estimation, it is necessary to know the number of signal sources. The signal sources can be estimated according to the distribution of the covariance matrix eigenvalues [41]. On the other hand, considering the MIMO sonar configuration and the dimension reduction transformation, the DOF of the proposed method is DOF = M + N − 2, i.e., the number of detectable signal sources.
Up to now, we have achieved the noise suppression method for DOA estimation in co-located MIMO sonar. We show the major steps of the proposed method as follows: i.
Remove the real part of the covariance matrix R YY . The imaginary part R I is used for DOA pre-estimation. ii.
Construct G matrix, then use the sorting refactoring method and the dimension reduction transformation to reconstruct the real part of R YY . We can obtain the new covariance matrix R. iii.
Reduce the coherence of R by the Toeplitz method. iv.
Construct W and P from G, and then use P left multiply R and use P H right multiply R.
Through the dimension reduction transformation, we can get the covariance matrix R z . v.
Compute the noise subspace and the signal subspace, finally estimating DOAs.
The choice of the twice DOA estimation methods in the steps will be discussed in the following section.

Performance Analysis
In this section, we analyze the performance of the noise suppression technique for MIMO sonar DOA estimation. For simplicity, RC-STIM is used to represent the noise suppression method we proposed.
In the following subsection, we assume a uniform linear transmit array consists of M = 8 omnidirectional sensors. The spacing of the elements is half a wavelength. The receive array is the same as the transmit array; the centers of the two arrays are overlapped. We compare the aforementioned methods in terms of the computational complexity, Cramer-Rao bound (CRB) and SNR gain methods.

Computational Complexity Analysis
The RC-STIM method we proposed needs to be combined with the appropriate DOA estimation algorithm to complete direction-finding. Considering the characteristics of the MIMO sonar array, we adopt the MUSIC algorithm and the ESPRIT algorithm. For achieving the effect of suppressing noise, we need to perform two DOA estimation operations. In order to obtain the computational complexity of the proposed algorithm, we can analyze the complexity of the traditional MUSIC algorithm (hereinafter referred to as MUSIC). The computational complexity of the MUSIC algorithm is mainly affected by the covariance matrix calculation, eigenvalue decomposition (EVD) and spectral peak search calculation. Compared with the MUSIC algorithm, the proposed method requires higher computational complexity. On the other hand, the dimension reduction transformation can help to reduce the computational complexity. For convenience, we abbreviate the MUSIC-based dimension reduction method as RC-MUSIC. While O{(M + N − 1) 2 } presents the MUSIC method which needs higher computational cost than RC-MUSIC. n is the number of searching steps. Hence, the total computational complexity of the proposed method is Similarly, we can get the computational complexity of the ESPRIT-based RC-STIM algorithm. The abbreviated definition of RC-ESPRIT is the same as above. Different from the MUSIC algorithm, the ESPRIT algorithm avoids the computational load brought by spectral peak search and the calculating speed is more prominent than the MUSIC-based methods. Hence, the total computational complexity of the proposed method is Through all these equations, we can see that the ESPRIT-based methods outperform the MUSIC-based methods in computational complexity. In order to achieve noise suppression, we perform DOA pre-estimation, which leads to an increase in computational complexity. In this case, we use the dimension reduction transformation and search-free direction-finding techniques to further reduce the computational complexity.

Cramer-Rao Bound
In our framework, a useful statistical bound for evaluating the limiting DOA estimation performance is the CRB. In this section, we analyze the CRB of MIMO sonar DOA estimation accuracy. It can be derived from the considered reconstructed covariance matrix, i.e., the model expressed by Equation (27). It is worth mentioning that the proposed MIMO sonar model with the reconstructed covariance matrix differs from the traditional MIMO sonar model given in Equation (4) due to the fact that the reduced dimension transformation method is used in Equation (27). Because of the irrelevance of noise and signals, the covariance matrix of noise can be written as According to Refs. [42,43], we derive the CRB estimation as where s(t l )s H (t l ) and stands for the Schur-Hadamard matrix product.

SNR Gain
Assume that the noise source signals that generate colored noise are irrelevant, therefore, the covariance matrix of colored noise can be modeled as γ and we assume that the imaginary part of the main diagonal is diag(γ 1 , γ 2 , . . . , γ MN ). We can obtain the SNR of the received signals from Equation (8) The SNR after dimension reduction can be obtained from Equation (27). It is worth noting that the processing of the covariance matrix reconstruction helps to reduce the symmetric noise component. Hence, we can obtain We can achieve the extra SNR gain from the dimension reduction processing.

Simulation Results
Throughout our simulations, we adopt the co-located MIMO array. We assume that the transmit and receive arrays are both uniform linear arrays with half wavelength spacing of each element and M = N = 8. Suppose that the direction of the two targets are θ 1 = −15 • , θ 2 = 20 • . The additive noise is spatial colored noise and the production process has been shown in Figure 3. We use several simulations to compare the performances of the proposed DOA methods, especially MUSIC and ESPRIT.
In this paper, we use the proposed covariance matrix reconstruction and dimension reduction method to improve the noise suppression ability. At this time, the subspace decomposition-based techniques such as MUSIC and the search-free direction-finding techniques such as ESPRIT can be used for multiple targets' DOA estimation. In all simulations, unless otherwise stated, all methods are computed based on 500 independent runs.
We evaluate the DOA estimation performance of our algorithms and present the root-mean-square error (RMSE) as where MONT denotes the number of Monte Carlo experiments,θ (n) k denotes the estimated angle in the nth experiment. Figure 4 shows the RMSE for the MUSIC-based and ESPRIT-based DOA estimators versus SNR for all the method test results. The noise environment is spatial colored noise and white Gaussian noise. A total number of snapshots L = 50 is used. It can be seen from this figure that the MUSIC-based MIMO sonar with the proposed RC-STIM method outperforms the traditional MIMO sonar MUSIC algorithm at the low SNR region while the opposite occurs at the high SNR region, the same situation to the other two methods based on ESPRIT. This means that the traditional algorithm is more sensitive to colored noise at the low SNR region, while the covariance matrix reconstruction operation helps to reduce the influence of the symmetric colored noise part on the DOA estimation. It can also be observed that the two methods based on ESPRIT show a similar DOA estimation accuracy. This is because under limited snapshots, the estimation performance is limited by ESPRIT. In practice, if we need to achieve a balance between the computational complexity and the DOA estimation accuracy, we can combine the MUSIC and ESPRIT algorithms, this will be considered in the following section. It can be seen from Figure 5 that in the white Gaussian noise environment the noise suppression performance of ESPRIT-based RC-STIM algorithm is prominent at the low SNR region, while the superiority of MUSIC-based algorithms is feasible at the high SNR region. Therefore, the ESPRIT-based algorithms can be selected when the white Gaussian noise is the main noise component. It can also be observed from Figure 5 that the covariance matrix reconstruction and dimension reduction processing can efficiently reduce the sensitivity of the algorithm to white Gaussian noise. Furthermore, we consider that the DOA estimation is considered to be solved when the following is satisfied It can be seen from Figure 5 that in the white Gaussian noise environment the noise suppression performance of ESPRIT-based RC-STIM algorithm is prominent at the low SNR region, while the superiority of MUSIC-based algorithms is feasible at the high SNR region. Therefore, the ESPRIT-based algorithms can be selected when the white Gaussian noise is the main noise component. It can also be observed from Figure 5 that the covariance matrix reconstruction and dimension reduction processing can efficiently reduce the sensitivity of the algorithm to white Gaussian noise. It can be seen from Figure 5 that in the white Gaussian noise environment the noise suppression performance of ESPRIT-based RC-STIM algorithm is prominent at the low SNR region, while the superiority of MUSIC-based algorithms is feasible at the high SNR region. Therefore, the ESPRIT-based algorithms can be selected when the white Gaussian noise is the main noise component. It can also be observed from Figure 5 that the covariance matrix reconstruction and dimension reduction processing can efficiently reduce the sensitivity of the algorithm to white Gaussian noise. Furthermore, we consider that the DOA estimation is considered to be solved when the following is satisfied Furthermore, we consider that the DOA estimation is considered to be solved when the following is satisfied where ∆θ = |θ 2 − θ 1 |, andθ k denotes the estimation of θ k . The source resolution probability versus SNR for all the methods we proposed is shown in Figure 6. It can be seen from the figure that the probability of the source resolution for each method starts to grow at a certain point as the SNR increases.
All methods, except the ESPRIT-based RC-STIM method, exhibit a 100% correct source resolution probability at a high SNR value. It can be seen from Figure 6 that the ESPRIT-based RC-STIM method has the highest SNR threshold. This is because the direction-finding precision of the ESPRIT-based algorithms is obviously affected by the preliminary DOA estimation processing and the limitation of algorithm DOA estimation accuracy will influence the reconstruction of the covariance matrix. The SNR thresholds of the ESPRIT method, the MUSIC method and the MUSIC-based RC-STIM method are lower than the aforementioned method, at about 0 dB. Moreover, except for the two RC-STIM algorithms, the target resolution probability of the remaining two algorithms is close to zero as SNR decreases to −10 dB and −20 dB. The ESPRIT-based RC-STIM algorithm is better than these two methods, the target resolution probability will not decrease to zero. Finally, the MUSIC-based RC-STIM method remains at a target resolution probability above 50%, i.e., the best performance of the source resolution probability.
ˆ, 1, 2,..., 2 where 2 1 θ θ θ Δ = − , and ˆk θ denotes the estimation of k θ . The source resolution probability versus SNR for all the methods we proposed is shown in Figure 6. It can be seen from the figure that the probability of the source resolution for each method starts to grow at a certain point as the SNR increases. All methods, except the ESPRIT-based RC-STIM method, exhibit a 100% correct source resolution probability at a high SNR value. It can be seen from Figure 6 that the ESPRIT-based RC-STIM method has the highest SNR threshold. This is because the direction-finding precision of the ESPRIT-based algorithms is obviously affected by the preliminary DOA estimation processing and the limitation of algorithm DOA estimation accuracy will influence the reconstruction of the covariance matrix. The SNR thresholds of the ESPRIT method, the MUSIC method and the MUSIC-based RC-STIM method are lower than the aforementioned method, at about 0dB . Moreover, except for the two RC-STIM algorithms, the target resolution probability of the remaining two algorithms is close to zero as SNR decreases to 10dB − and 20dB − . The ESPRIT-based RC-STIM algorithm is better than these two methods, the target resolution probability will not decrease to zero. Finally, the MUSIC-based RC-STIM method remains at a target resolution probability above 50%, i.e., the best performance of the source resolution probability. It is worth noticing that, in the case of low SNR, the MUSIC-based RC-STIM algorithm has a strong suppression of colored noise, while the ESPRIT-based RC-STIM algorithm can effectively reduce the impact of the white Gaussian noise. Considering the computational complexity and noise suppression, we can try to combine the proposed two methods. The computational complexity of the proposed methods has been discussed above. To access a higher DOA estimation accuracy and lower computational complexity to improve the algorithm performance, we utilize ESPRIT to carry out the preliminary DOA estimation before covariance matrix construction processing, and take the advantage of MUSIC in final DOA estimation to ensure the target direction-finding precision. Another option is to reverse the two methods. From the above, the computational complexity of this method is It is worth noticing that, in the case of low SNR, the MUSIC-based RC-STIM algorithm has a strong suppression of colored noise, while the ESPRIT-based RC-STIM algorithm can effectively reduce the impact of the white Gaussian noise. Considering the computational complexity and noise suppression, we can try to combine the proposed two methods. The computational complexity of the proposed methods has been discussed above. To access a higher DOA estimation accuracy and lower computational complexity to improve the algorithm performance, we utilize ESPRIT to carry out the preliminary DOA estimation before covariance matrix construction processing, and take the advantage of MUSIC in final DOA estimation to ensure the target direction-finding precision. Another option is to reverse the two methods. From the above, the computational complexity of this method is Figure 7 shows the computational complexity of the MUSIC-based and ESPRIT-based DOA estimations. It can be seen from this figure that the MUSIC-ESPRIT-based method of Equation (37) with a covariance matrix construction operation and dimension reduction processing has an advantage in computational complexity, the computational load is similar to the ESPRIT-based methods. Moreover, with the increasing of snapshots, this advantage will gradually decrease. Therefore, under limited snapshot conditions, the performance loss due to a limited snapshot number is proportional to the covariance matrix dimension. The lower the covariance matrix dimension, the less performance loss. As a result, the method we proposed gains better performance than traditional algorithms versus the limited snapshot number. On the other hand, the method of Equation (38) shows a slight advantage in improving computational complexity. However, this method better preserves the performance of MUSIC algorithm in the preliminary DOA estimation. In the processing of preliminary DOA estimation, the precision of direction-finding directly affects the accuracy of covariance matrix reconstruction. Figure 7 shows the computational complexity of the MUSIC-based and ESPRIT-based DOA estimations. It can be seen from this figure that the MUSIC-ESPRIT-based method of Equation (37) with a covariance matrix construction operation and dimension reduction processing has an advantage in computational complexity, the computational load is similar to the ESPRIT-based methods. Moreover, with the increasing of snapshots, this advantage will gradually decrease. Therefore, under limited snapshot conditions, the performance loss due to a limited snapshot number is proportional to the covariance matrix dimension. The lower the covariance matrix dimension, the less performance loss. As a result, the method we proposed gains better performance than traditional algorithms versus the limited snapshot number. On the other hand, the method of Equation (38) shows a slight advantage in improving computational complexity. However, this method better preserves the performance of MUSIC algorithm in the preliminary DOA estimation. In the processing of preliminary DOA estimation, the precision of direction-finding directly affects the accuracy of covariance matrix reconstruction. Since the computational complexity of the algorithms is related to the snapshots, we reduce the snapshots to 10 L = . In this way, we can observe the effect of the snapshots on the proposed algorithms. Figure 8 shows the RMSE versus SNR respectively. Compared with Figure 4, the RMSE of the four corresponding algorithms in Figure 8 is not significantly different. Therefore, the snapshots have little effect on the DOA estimation accuracy of the algorithms. Obviously, Method (37) shows better performance than Method (38). In the low signal area, Method (37) outperforms all the other methods. With the increase of SNR, Method (37) gradually lost this advantage. The target resolution probability of the four algorithms in Figure 6 in the condition of 10 L = are shown in Figure 9. In this figure, the target resolution probabilities of the shown methods are all reduced, compared with Figure 6. The ESPRIT-based RC-STIM method failed to reach a 100% correct source resolution at SNR 20dB = . Meanwhile, the SNR threshold of the remaining three algorithms is about 15dB , which is much higher than the situation in Figure 6. The target resolution probability of the two traditional algorithms decreases to 0 at SNR 10dB = − and SNR 5dB = − . The ESPRIT-based RC-STIM method shows better performance, which will not go down to 0. The MUSIC-based RC-STIM method outperforms all the other methods. In brief, the target resolution probability of all the methods are affected by snapshots, meanwhile, the RC-STIM methods show a lower sensitivity to snapshots. Since the computational complexity of the algorithms is related to the snapshots, we reduce the snapshots to L = 10. In this way, we can observe the effect of the snapshots on the proposed algorithms. Figure 8 shows the RMSE versus SNR respectively. Compared with Figure 4, the RMSE of the four corresponding algorithms in Figure 8 is not significantly different. Therefore, the snapshots have little effect on the DOA estimation accuracy of the algorithms. Obviously, Method (37) shows better performance than Method (38). In the low signal area, Method (37) outperforms all the other methods. With the increase of SNR, Method (37) gradually lost this advantage. The target resolution probability of the four algorithms in Figure 6 in the condition of L = 10 are shown in Figure 9. In this figure, the target resolution probabilities of the shown methods are all reduced, compared with Figure 6. The ESPRIT-based RC-STIM method failed to reach a 100% correct source resolution at SNR = 20 dB. Meanwhile, the SNR threshold of the remaining three algorithms is about 15 dB, which is much higher than the situation in Figure 6. The target resolution probability of the two traditional algorithms decreases to 0 at SNR = −10 dB and SNR = −5 dB. The ESPRIT-based RC-STIM method shows better performance, which will not go down to 0. The MUSIC-based RC-STIM method outperforms all the other methods. In brief, the target resolution probability of all the methods are affected by snapshots, meanwhile, the RC-STIM methods show a lower sensitivity to snapshots.    (37) and (38). For comparison, we consider different snapshots, i.e., 10 L = and 50 L = . Method (37) has better performance compared with Method (38), regardless of the value of snapshot. Method (37) successfully reaches a 100% correct source resolution at SNR 5dB = , both at 10 L = and 50 L = . This situation is obviously better than the algorithms in Figure 9, however, it is slightly worse than the algorithms in Figure 6. This means that the performance of Method (37) is not sensitive to changes of snapshots, which helps to provide better performance in limited snapshot conditions. Meanwhile, the target resolution probability of Method (38) is not sensitive to snapshots in a low SNR situation. As SNR increases, the SNR threshold of 50 L = is obviously better than 10 L = . Compared to the aforementioned four algorithms, the SNR threshold of Method (38) is higher at snapshots 50 L = , and basically consistent with them in limited snapshots. As a result, Method (37) outperforms all the mentioned methods in limited snapshots.    (37) and (38). For comparison, we consider different snapshots, i.e., 10 L = and 50 L = . Method (37) has better performance compared with Method (38), regardless of the value of snapshot. Method (37) successfully reaches a 100% correct source resolution at SNR 5dB = , both at 10 L = and 50 L = . This situation is obviously better than the algorithms in Figure 9, however, it is slightly worse than the algorithms in Figure 6. This means that the performance of Method (37) is not sensitive to changes of snapshots, which helps to provide better performance in limited snapshot conditions. Meanwhile, the target resolution probability of Method (38) is not sensitive to snapshots in a low SNR situation. As SNR increases, the SNR threshold of 50 L = is obviously better than 10 L = . Compared to the aforementioned four algorithms, the SNR threshold of Method (38) is higher at snapshots 50 L = , and basically consistent with them in limited snapshots. As a result, Method (37) outperforms all the mentioned methods in limited snapshots.  Figure 10 shows the target resolution probability of Methods (37) and (38). For comparison, we consider different snapshots, i.e., L = 10 and L = 50. Method (37) has better performance compared with Method (38), regardless of the value of snapshot. Method (37) successfully reaches a 100% correct source resolution at SNR = 5 dB, both at L = 10 and L = 50. This situation is obviously better than the algorithms in Figure 9, however, it is slightly worse than the algorithms in Figure 6. This means that the performance of Method (37) is not sensitive to changes of snapshots, which helps to provide better performance in limited snapshot conditions. Meanwhile, the target resolution probability of Method (38) is not sensitive to snapshots in a low SNR situation. As SNR increases, the SNR threshold of L = 50 is obviously better than L = 10. Compared to the aforementioned four algorithms, the SNR threshold of Method (38) is higher at snapshots L = 50, and basically consistent with them in limited snapshots. As a result, Method (37) outperforms all the mentioned methods in limited snapshots.  (37) and (38) in the case of different snapshot number. Figure 11 presents the effect of array element number on the performance of the algorithms. Obviously, the performance of the algorithm is proportional to the number of array elements, i.e., for the same algorithm, the larger the number of array elements, the higher the direction-finding accuracy of this algorithm. This situation can be attributed to the diversity gain. In the scenario of 6 M N = = and low SNR, the MUSIC-based RC-STIM method outperforms the other two methods. As the SNR increases, when SNR reaches about 5dB − , we can gain better DOA estimation through the traditional MUSIC method. However, under the condition of 10 M N = = , Method (37) has the best DOA estimation accuracy, until the SNR is greater than 5dB − .  Figure 11. The RMSE versus SNR for MUSIC-based and ESPRIT-based DOA estimation in the case of different elements.

Conclusions
The influence of symmetric noise component on algorithm performance for multiple targets DOA estimation is proposed. A group of methods for noise suppression is introduced. The essence Figure 10. The probability of the target resolution versus SNR for Methods (37) and (38) in the case of different snapshot number. Figure 11 presents the effect of array element number on the performance of the algorithms. Obviously, the performance of the algorithm is proportional to the number of array elements, i.e., for the same algorithm, the larger the number of array elements, the higher the direction-finding accuracy of this algorithm. This situation can be attributed to the diversity gain. In the scenario of M = N = 6 and low SNR, the MUSIC-based RC-STIM method outperforms the other two methods. As the SNR increases, when SNR reaches about −5 dB, we can gain better DOA estimation through the traditional MUSIC method. However, under the condition of M = N = 10, Method (37) has the best DOA estimation accuracy, until the SNR is greater than −5 dB.  (37) and (38) in the case of different snapshot number. Figure 11 presents the effect of array element number on the performance of the algorithms. Obviously, the performance of the algorithm is proportional to the number of array elements, i.e., for the same algorithm, the larger the number of array elements, the higher the direction-finding accuracy of this algorithm. This situation can be attributed to the diversity gain. In the scenario of 6 M N = = and low SNR, the MUSIC-based RC-STIM method outperforms the other two methods. As the SNR increases, when SNR reaches about 5dB − , we can gain better DOA estimation through the traditional MUSIC method. However, under the condition of 10 M N = = , Method (37) has the best DOA estimation accuracy, until the SNR is greater than 5dB − .  Figure 11. The RMSE versus SNR for MUSIC-based and ESPRIT-based DOA estimation in the case of different elements.

Conclusions
The influence of symmetric noise component on algorithm performance for multiple targets DOA estimation is proposed. A group of methods for noise suppression is introduced. The essence

Conclusions
The influence of symmetric noise component on algorithm performance for multiple targets DOA estimation is proposed. A group of methods for noise suppression is introduced. The essence of the aforementioned methods is covariance matrix reconstruction and dimension reduction. Unlike the previous methods that directly utilize the covariance matrix to do DOA estimation, a priori processing, i.e., pre-estimation, is used to decrease the effect of the symmetric noise component.
Dimension reduction processing can help to reduce the computational complexity. MUSIC and ESPRIT can then be used for multiple targets' DOA estimation. The ESPRIT-based methods are more advantageous in terms of reducing computational complexity. Simultaneously, Method (37) allows us to achieve a better DOA estimation performance in the lower SNR region. Performance analysis of the proposed noise suppression methods and comparison to the existing MIMO sonar techniques are given. The computational complexity of the proposed methods can be controlled by the selection of the snapshots and the prior DOA estimation algorithm. The performance of the MUSIC-based RC-STIM method and Method (37) outperform the other methods. Simulation results show the superiority of the two proposed methods over the existing methods.
In future work, we intend to further optimize the noise suppression algorithm model, reduce the impact of the pre-estimation result on the DOA estimation accuracy. On the other hand, the pre-estimation result can be utilized to design the transmit beam space to improve SNR. Thus, the DOA estimation accuracy can be improved.