A Velocity-independent DOA Estimator of Underwater Acoustic Signals via an Arbitrary Cross-linear Nested Array

In this paper, an estimator for underwater DOA estimation is proposed by using a cross-linear nested array with an arbitrary cross-angle. The estimator removes the acoustic velocity factor by deriving the geometric relation of the cross-linear array. Compared with conventional DOA estimation algorithms via a linear array, this estimator eliminates systematic errors caused by the uncertainty factor of the acoustic velocity in the underwater environment. In comparison with the conventional acoustic velocity-independent algorithm, this estimator uses the nested array and does not have to estimate the value of the acoustic velocity, which improves the performance of DOA estimation. Moreover, the proposed method is only slightly inferior to some of the conventional algorithms if the acoustic velocity is estimated accurately. Numerical simulations are provided to validate the analytical derivations and corroborate the improved performance in underwater environments where the actual acoustic velocity is not accurate.

(1) The number of identifiable sources with ULAs is limited to the number of array elements minus one; (2) The accuracy of DOA estimation is also limited by the inherent size of the array aperture; (3) Unable to deal with the DOA estimation problem in some special environments effectively, such as underwater communication with inaccurate acoustic velocity; (4) The estimated angular range is limited to the specific side that the sensor array is pointing at.
The 1D-DOA estimation under the condition of unknown velocity is essentially a multi-parameter estimation problem, which is similar to the 2D-DOA estimation under the condition that the speed of sound is known. To solve this type of problem, numerous various array geometries have been proposed, such as uniform planar arrays [1,12], uniform circular arrays [10,13], L-shaped arrays [21,24], double L-shaped arrays [8], and so on. Among them, L-shaped arrays are more widely used for their relative simplicity in the data processing. In addition, sparse linear arrays (SLAs), such as co-prime [11,26], and nested [4,5,16], are significantly better than conventional ULAs in terms of estimation accuracy.
At present, to solve the 1D-DOA estimation problem under the condition of unknown acoustic velocity, we had given some feasible velocity-independent (VI) methods, VI-ESPRIT [15] and FVI-ESPRIT [14]. However, both of them are via ULAs, which are inferior to SLAs in terms of accuracy. Further, these methods can theoretically only accomplish DOA estimation on a specific side of the coordinate axis the array is located on. Finally, these methods do not provide an explicit matching scheme for the multi-source case.
Inspired by the respective properties of L-shaped arrays and SLAs, we arranged a cross-linear nested array in this paper. On this basis, a velocity-independent ESPRIT using Toeplitz Matrix Reconstruction (TVI-ESPRIT) algorithm is proposed. The Toeplitz covariance matrices of the two linear arrays are reconstructed, and then the geometric relationship of the crossed line arrays is exploited to eliminate the impact of the acoustic velocity factor. Finally, the effective DOA estimation is accomplished by the ESPRIT algorithm.
Compared to the previous work, the contributions of this paper are summarized as: (1) We derive the principle of eliminating the acoustic velocity factor by the crosslinear array with an arbitrary angle and give an approach to estimate the underwater acoustic velocity and a matching scheme for the rotation factor in the case of multiple sources; (2) We analyze in detail eight different spatial regions of the signal arrival direction and give a unified DOA estimation method; (3) We combine the advantages of L-shaped arrays and sparse line arrays to ensure the validity and accuracy of DOA estimation in environments with inaccurate acoustic velocities.

Data Model and ESPRIT Algorithm
Assume K far-field, uncorrelated narrowband acoustic signals from distinct direction {θ 1 , θ 2 , . . . , θ K } impinge on a nested linear array (NLA) with M sensors. Considering that each sensor collects N snapshots, the signal observed matrix X ∈ C M×N of NLA can be expressed as The second order ULA where A(θ ) = [a(θ 1 ), . . . , a(θ K )] ∈ C M×K is the manifold matrix made up of the array steering vectors a(θ k ), k = 1, . . . , K . S ∈ C K ×N denotes the signal sampling matrix and N ∈ C M×N is the Additive White Gaussian Noise matrix with zero mean and σ 2 n variance of the entries. The NLA usually refers to the combined array composed of two ULAs connected in series, as shown in Fig. 1. The number of elements in the first-order ULA is M 1 , and the interval between sensors is d 1 (Usually half of the signal wavelength). And the number of elements in the second-order ULA is M 2 , and the interval between sensors is d 2 = (M 1 + 1)d 1 . Hence, the sensor layout set D (in ascending order) of the NLA can be expressed as Then, the covariance matrix R x of the NLA can be expressed as where R s ∈ C K ×K is the correlation matrix of the signal. The observed covariance matrix R x is then converted to a vector by the property vec { AY B} = B T ⊗ A vec{Y }: In the above equation, represents the Khatri-Rao product and ⊗ denotes the Kronecker product. p = σ 2 1 , σ 2 2 , . . . , σ 2 K T ∈ C K ×1 represents the power vector of the signals and 1 n is the vectorized identity matrix. We can regard the z as the observed signal of a single snapshot from the coherent source p. The matrix A * A analogues to an extended manifold matrix with the equivalent sensors located at S = {D(n 1 ) − D(n 2 ) : n 1 , n 2 = 1, 2, . . . , where the location set D of the real sensors is defined in Eq. (2) and note that M = M 1 + M 2 .
In the set S, the elements are not arranged in ascending order, and some of them are duplicated. We should sort the elements in the vector z and remove redundancy to get a vector that matches the order of ULAs [16]. Then a new equivalent signal vector is formed,z ∈ C M sr ×1 . The number of equivalent sensors is M sr = max S−min S d 1 . A r ∈ C M sr ×K is its manifold matrix and 1 n denotes an all-zero vector except 1 at the M sr +1 2 th entry. In general, to obtain the maximum extended aperture (or degree of freedom, DOF) of the array, the optimal secondary nested array is laid out as the following arrangement: If the value of M is set to an even number (here, M sr = M 2 −2 2 + M), then the manifold matrix A r in Eq. (6) can be expressed as Meanwhile, the location set S r of the sensors is Since the rank of the signal matrix p in Eq. (6) is 1, which no longer satisfies the uncorrelated condition. A common approach is to reconstruct a Toeplitz matrix R t p to recover the rank of the matrix so that it is equal to the number of sources K . The size of the Toeplitz matrix is M t p = ( M 2 4 + M 2 ). Note thatz(n) denotes the nth element of the vectorz.
The ESPRIT algorithm is then used to obtain the rotational phase information of the array. First, by performing the eigenvalue decomposition (EVD) of the matrix R t p , the signal subspace U s ∈ C M tp ×K spanned by K larger eigenvalues corresponding eigenvectors can be obtained. Then, the signal subspace is divided into the signal subspace (U s1 ∈ C (M tp −1)×K and U s2 ∈ C (M tp −1)×K ) of two subarrays, respectively.
Their corresponding steering matrices can be expressed as A rotation matrix Φ relationship between the steering submatrices A r 1 and A r 2 exists because of the unit space ( = d 1 ) between the two subarrays.
In addition, due to the correspondence between the signal subspace U s and the steering matrix A r , a non-singular unitary matrix T ∈ C K ×K should be defined, which satisfies the following relation: Further, we can get the following result: By using the property that the diagonal elements of matrix Φ are equal to the eigenvalues of matrix Ψ , the diagonal elements of matrix Φ can be estimated. And the arrival angle of the signal can be obtained, too. However, the real-time velocity c real is in the range of 1450−−1550 m/s. In general, we set the expected acoustic speed to c 0 = 1500 m/s, which is an error between this and the actual acoustic velocity. And the bias of DOA estimation increases with the increase of acoustic velocity error (see Eq. (26)).

Proposed Method
In this section, a DOA estimation method is proposed to eliminate the underwater acoustic velocity error. We use a cross-linear nested array with the crossed angle δ, as shown in Fig. 2.
The sensor at the intersection of two linear arrays is regarded as the reference array element O. We divide the one-dimensional space into eight parts by taking the axes where the two linear arrays are located and their normal lines as the boundary. Suppose that the coordinate axes of the two linear arrays are, respectively, X-axis and Y-axis, and the angle between the incoming direction of the kth signal (k = 1, 2, . . . , K ) and the positive direction of the X-axis normal is x nk , and the angle between the incoming direction and the positive direction of the Y-axis normal is y nk . Note that x nk , y nk ∈ [0, π].
When the signal is incoming from area 1, then y nk − x nk = δ; when the signal is incoming from area 2, then y nk − x nk = δ; when the signal is incoming from area 3, then y nk + x nk = δ; when the signal is incoming from area 4, then x nk − y nk = δ; when the signal is incoming from area 5, then x nk − y nk = δ; when the signal is incoming from area 6, then x nk − y nk = δ; when the signal is incoming from area 7, then x nk + y nk = 2π − δ; when the signal is incoming from area 8, then y nk − x nk = δ; Summarizing the above four cases, we get the following conclusions: DOA in the area 4,5,6 x nk + y nk = 2π − δ DOA in the area 7 (16) When the incident signal reaches the reference element O first, the incoming direction of the X-axis or Y-axis (θ xk , θ yk ) is negative. Otherwise, when the incident signal reaches the reference element O last, the direction of incoming wave (θ xk , θ yk ) is positive.
When the signal is incoming from area 1, then θ xk = x nk , θ yk = y nk ; when the signal is incoming from area 2, then θ xk = x nk , θ yk = y nk ; when the signal is incoming from area 3, then θ xk = −x nk , θ yk = y nk ; when the signal is incoming from area 4, then θ xk = −x nk , θ yk = −y nk ; when the signal is incoming from area 5, then θ xk = −x nk , θ yk = −y nk ; when the signal is incoming from area 6, then θ xk = −x nk , θ yk = −y nk ; when the signal is incoming from area 7, then θ xk = x nk , θ yk = −y nk ; when the signal is incoming from area 8, then θ xk = x nk , θ yk = y nk .
Combining the above analysis and Eq. (16), the relationship between the wave directions coming up on the X and Y axes in the four areas can be obtained: We can generalize the conclusion in Eq. (17) that the relation between the arrival angles of the same signal to two linear arrays is sin θ yk = sin (θ xk + δ) k = 1, 2, . . . , K By performing the eigenvalue decomposition of the matrices Ψ X and Ψ Y in Eq. (15), respectively, their rotation matrices Φ X and Φ Y can be obtained. Then, we list the rotational phase of the two axes separately, for the kth signal (k = 1, 2, . . . , K ).
Ideally, let φ x = arg (Φ xk ) and φ y = arg Φ yk . Figure 2 and Eq. (19) demonstrate that the incident region and the phase state exhibit the following properties in Table  1. From it, the objective information about the direction of signal incidence can be obtained (e.g., φ x ≥ 0 means that the incident region is in the set {①, ②, ⑦, ⑧}).
By Eq. (19), we have arg Φ yk From such a process, it is clear that the acoustic velocity is eliminated by the reduction of the fraction. Thus, it eliminates the error caused by acoustic velocity in DOA estimation at the same time. Next, substitute Eq. (18) into Eq. (20), and then the kth DOA value is estimated (X -axis as the main axis).
where the arccotangent function arccot(·) has the value domain of [0, π] that corresponds to the incident region ①, ②, ⑦ and ⑧ in Fig. 2. For a single signal, the algorithm has completed DOA estimation. However, in the case of multiple signals in the channel, some constraints need to be added to Eq. (20) to complete their pairing. For the receiver sensors, the acoustic velocity received by each of them is similar (from 1450 m/s to 1550 m/s). After deducting Eq. (19), we can get the calculation expression of the kth acoustic signal velocity aŝ Then, a total of K ! paired combinations of arg Φ yk and arg Φ xk are available. For each combination, Eq. (22) can be used to find the acoustic velocities ĉ k , k = 1, 2, . . . , K of K signals. Finally, the group with the smallest variance of the acoustic velocity estimated corresponds to the correct matching combination.
Since the proposed velocity-independent algorithm works on Toeplitz matrix reconstruction technology, it is named TVI-ESPRIT method, which is summarized in Algorithm 1.
Moreover, the Cramer-Rao Bound (CRB) of DOA estimation for the arrays is derived from the literature [7].
where γ denotes the ratio form of the signal-to-noise ratio, and x m ,y m denote the positions of the array elements along the x and y coordinate axes. And the total CRB

Input:
A cross-linear nested array with M sensors at a cross-angle of δ; The receiving matrix of array signals with N snapshots, X ∈ C M×N ; Output: 1: Calculate and vectorize the signal covariance matrices (R x and R y ) of two nested line arrays by Equations (3) and (4), respectively. z x and z y can be obtained; 2: Reorder them and get rid of the redundancy and obtainz x andz y as Equation (6); 3: Construct the Toeplitz matrices applicable to the two axes, respectively, by Equation (10); 4: Utilize ESPRIT algorithm to calculate the matrices x and y of two linear arrays by Equations (11) and (15) of 1D DOAs can be expressed by the following equation.
According to Eqs. (24) and (25), it can be concluded that the CRB of 1D DOAs decreases as the number of snapshots or SNR increases. And the array structure also has an impact on the CRB, including the number and location of sensors.

Numerical Results
In this section, we evaluate the performance of the proposed TVI-ESPRIT algorithm and compare it with the other ESPRIT algorithms. These algorithms include GLS-ESPRIT [22] algorithm via a ULA, SS-ESPRIT [16] algorithm via a NLA and VI-ESPRIT [15] algorithm via L-shaped ULA.
To be fair, the array used by all the algorithms includes M = 11 sensors, and the cross-angle of the cross-linear nested array is set δ = 90 • , which corresponds to the L-shaped array used by the VI-ESPRIT algorithm.
Unless otherwise specified, all numerical results are obtained for the case of sources with f = 10 kHz center frequency by means of W = 1000 Monte Carlo trials with f s = 25 kHz sampling frequency and N = 200 snapshots. Nextly, the SNR of the signal is set to 0 dB in the Additive White Gaussian Noise situation, and the incoming directions are θ = {30 • , 60 • }.
The expected acoustic velocity is set to 1500 m/s. The unit spacing between the sensors is half of the expected wavelength. Note that c = c real − c 0 , where c 0 = 1500 m/s. After deducting, the theoretical expression of the error caused by acoustic velocity error is where c real represents the real-time acoustic velocity. Error c is used as the reference of simulation experiment to reflect the influence of velocity error on estimation accuracy.

Effect of array structure on estimation performance
In the first trial, the impact of the structure on the estimation accuracy is theoretically analyzed in terms of CRB. Figure 3a and b shows that the nested array structure is superior to the ULA in terms of theoretically maximum estimation accuracy, no matter the cross-array and the linear array. In addition, the accuracy of the estimation of the crossed array structure is inferior to that of the linear array for the same number of array elements. Therefore, the NLA can theoretically obtain the optimal estimation accuracy without considering the effect of conditions such as acoustic velocity errors.
However, the effect of acoustic velocity errors on the estimation can be eliminated with the help of the cross-array geometric relationship (as demonstrated in Sect. 3). Thus, the TVI-ESPRIT algorithm proposed in this paper chooses a cross-nested array structure to combine their respective advantages.

Comparison of Algorithms with Different SNR
In the second trial, the other methods are compared in terms of RMSE (Root Mean Square Error) with the TVI-ESPRIT algorithm under different SNR conditions. RMSE is defined as follows: The acoustic velocity error is set as 0 m/s and 20 m/s under each SNR condition. As shown in Fig. 4a, TVI-ESPRIT algorithm is inferior to SS-ESPRIT algorithm in terms of accuracy but outperforms GLS-ESPRIT and VI-ESPRIT algorithm. According to Fig. 4b, all methods via the ULA inevitably produce systematic errors caused by inaccurate acoustic velocity, no matter how large the SNR value is. Due to the inaccuracy of the acoustic velocity, their error curves overlap highly with those of the theory. In addition, the TVI-ESPRIT algorithm proposed in this paper greatly improves the performance compared with the VI-ESPRIT algorithm.

Comparison of Algorithms with Different Number of Snapshots
In the third trial, the above algorithms are compared with respect to the number of sampling snapshots. Figure 5a demonstrates that the TVI-ESPRIT algorithm performs better than the other, although its performance is not as good as that of the SS-ESPRIT without velocity error. And in Fig. 5b, the RMSEs of the algorithms via the ULA are overlapped with the velocity error curves, no matter how large their numbers of snapshots are. Obviously, their major error in estimation originates from the uncertainty in the velocity. In addition, the accuracy of the proposed algorithm is generally better than that of VI-ESPRIT while they are both via an L-shaped array (δ = 90 • ). The conclusion of this experiment is similar to that of the previous trial.

Comparison of Algorithms with Inaccurate Acoustic Velocity
In the last trial, the above algorithms are compared in terms of RMSE with respect to specific acoustic velocity error. The other simulation experiment parameters are the same as those of the overall setting. From the two results ( Fig. 6a and b), the estimation error of the ESPRIT algorithms via the ULA increases as the absolute value of the acoustic velocity error | c| increases. Meanwhile, the estimation accuracy of the two kinds of VI algorithms is almost not affected by the acoustic velocity error. Moreover, the TVI-ESPRIT outperforms the VI-ESPRIT algorithm in terms of estimation accuracy, especially in low SNR conditions. Then, in the case of low SNR, the estimation errors of the conventional algorithms via the ULA mainly originate from noise and inaccurate velocity. While at high SNR, their estimation errors mainly come from acoustic velocity. It is reflected by the fact that with the increment of SNR, their RMSE curves keep approaching the Error c curve caused by the theoretical error due to the acoustic velocity error. Finally, our TVI-ESPRIT is slightly inferior to SS-ESPRIT with a small velocity error, but it still outperforms other algorithms in the other cases.

Conclusion
In this paper, a velocity-independent algorithm via an arbitrary cross-linear nested array is proposed. The proposed method employs the rotational phase matrix of two crossed linear arrays to eliminate the acoustic velocity factor. Also, it extends the aperture of the array by arranging nested arrays to obtain a better estimation accuracy. In addition, it reaches the estimable angular range of (−π, π] compared to (−π/2, π/2] of ULAs. Simulation results verify that the proposed method is not susceptible to acoustic velocity in contrast to the other ESPRIT algorithms (e.g., GLS-ESPRIT). Finally, the proposed method outperforms the conventional VI-ESPRIT method in terms of estimation accuracy for the same cases.