Target azimuth is an important parameter for the identification, detection, positioning, and tracking of underwater targets (Luo and Shen, 2021; Chen et al., 2023; Xie et al., 2023; Zhao et al., 2023). Array signal processing has been shown to be effective for the direction of arrival (DOA) estimation. The methods can be classified as traditional beamforming, subspace-based, deconvolution (DCV), and transform domain.

The most commonly employed traditional beamforming method is the conventional beamforming (CBF). However, CBF has the disadvantages of wide beamwidths and poor spatial resolution due to the Rayleigh limit. Several high-resolution methods have been proposed to overcome these such as the maximum entropy algorithm (Burg, 1975) and the minimum variance distortionless response (MVDR) algorithm (Capon, 1969). These methods improve the spatial resolution but suffer from poor performance when used to detect broadband signals in actual ocean environments. This is because they are sensitive to signal mismatch and cannot estimate the DOA of coherent sound sources (Somasundaram, 2012). The steered minimum variance (STMV) algorithm was proposed to overcome these problems (Krolik and Swingler, 1989). STMV has better spatial resolution for coherent acoustic sources and fast convergence, but poor robustness (Somasundaram et al., 2015). The Rayleigh limit was overcome with the subspace-based algorithm multiple signal classification (MUSIC) (Schmidt, 1986). Subsequent subspace-based algorithms such as estimation of signal parameters via rotational invariance techniques (ESPRIT) (Roy and Kailath, 1989), root-multiple signal classification (RMUSIC) (Rao and Hari, 1989), maximum likelihood (ML) (Stoica and Nehorai, 1989), and weighted subspace fitting (WSF) (Bengtsson and Ottersten, 2001) provide improved performance but have sensitivity and snapshot deficiency problems when used in practical applications (Baggeroer and Cox, 1999). Another issue is that the number of acoustic sources is usually unknown and this makes it difficult to estimate the signal and noise subspaces. Further, existing algorithms can only be used to estimate DOA for incoherent or weak-coherent acoustic sources, making detection of coherent signals difficult in the actual ocean environments.

The above algorithms are either sensitive to array element errors or limited to array aperture. Various studies on the formation structure have been introduced to improve the performance of the algorithms (Zhang et al., 2021; Zhou et al., 2022; Yang, 2023; Ye et al., 2023). Additionally, in recent years, there has been a growing focus on researching robust high-resolution beamforming algorithms. Deconvolution (DCV) algorithms have attracted widespread attention for underwater acoustic applications. DCV was initially considered with both uniform linear arrays and circular arrays (Yang, 2017; Yang, 2018). It was shown that the performance is better than CBF. The super-directivity performance of DCV with a small-sized array was verified using the SwellEx96 horizontal array (Yang, 2019). However, these DCV methods are only suitable for arrays with a shift-invariant point spread function (PSF) beam pattern, such as a horizontal line array or circular array. Therefore, new DCV methods were developed for shift-variant PSF beam patterns. A DCV method based on non-negative least squares (NNLS) and an improved NNLS method called extended Richardson-Lucy (Ex-RL) were presented which provide high resolution, robustness, and excellent array gain (Sun et al., 2019; Sun et al., 2020).

Transform domain methods were originally developed to estimate seismic wave velocity and azimuth (Cheng et al., 2018). More specifically, the frequency-wavenumber (f-k) power spectrum can be obtained using the space-time two-dimensional Fourier transform of the seismic signal. Then, the f-k power spectrum can be converted into the transform domain to extract the velocity and azimuth of the waves. While transform domain methods have been widely used in seismic exploration (Zywicki and Rix, 1999), there have been few DOA estimation applications. The least squares line fitting (LSLF) algorithm was employed to obtain the slope of the local peak-energy line in the f-k domain and then the sum of the points on this line was used as an estimate of the energy output of the azimuth spectrum (Li et al., 2019). However, this method is sensitive to outliers in the image since it minimizes the sum of the squares of the distances from the points to the line. Thus, the performance can be degraded significantly, particularly in low signal-to-noise ratio (SNR) environments or when there are multiple adjacent targets. In this paper, a broadband high-resolution DOA estimation method based on the generalized weighted Radon transform (GWRT) is proposed. The array signal is converted into the f-k domain by solving the conditional wavenumber spectrum function (CWSF) and then the mathematical relationship between the spatial distribution of broadband signal energy in the f-k domain and target azimuth is obtained. To improve performance, image gradient information is utilized as weights for the GWRT, and a linear integral mathematical model is derived by the GWRT processing in the f-k domain. The resulting model contains the complete image information in the f-k domain. This is then converted into the azimuth-energy domain to realize high-resolution DOA estimation. The proposed method does not require prior knowledge of the number of sources or signal pre-estimation. In addition, it is not sensitive to outliers in the image and the results in the f-k domain provide higher transform gain and better robustness. Both simulation and sea-trial experiments are conducted to validate the proposed method. The results obtained indicate that the proposed method has better performance and offers several advantages compared with existing approaches as follows.

The remainder of this paper is organized as follows. Section 2 introduces the broadband signal model and DOA estimation using the CWSF is presented. In Section 3, we derive the expression of the mathematical model for DOA estimation using GWRT. The performance of the proposed method is evaluated via simulation and compared with other DOA methods in Section 4. The results of the sea-trial experiments are given in Section 5. Finally, Section 6 provides a summary of the paper.

We consider a linear sensor array with M receivers uniformly spaced at a distance d. The signal is assumed to be from a broadband source located at the far field of the array. This signal has a look direction θ and arrives at the array as a plane wave. Then, θ is the azimuth angle of the target, which is defined as the anticlockwise angle between the horizontal array and the target as shown in Figure 1 . Therefore, the signal received by element m at time t can be expressed as

where   m = 1 , 2 , ⋯ , M , s t is the incident signal, n m t is the noise received by element m, which is uncorrelated with s t , a m   is the array manifold factor of element m which is equal to exp ( j 2 π f ( m − 1 ) d cos θ / c ) , and   f   and c   are frequency and the speed of sound, respectively.

As the actual data processing is based on a discrete-time model, the received signal should be sampled at   f s . Therefore, the signal received by element m can be expressed as

where L is the number of samples which is an integer. The discrete Fourier transform (DFT) of x m l gives the corresponding frequency spectrum

Then, the frequency-wavenumber spectrum of the received signal can be obtained by applying the DFT to X m f with Q points in the spatial domain as

where k is the wavenumber which is an integer in the range − Q / 2 ≪ k < Q / 2 , Q is an integer with Q ≫ M , S f   is the frequency spectrum of s t , N m f is the frequency spectrum of n m t .

From Equation (4), the frequency-wavenumber power spectral density can be expressed as

where E { · }   denotes expectation, superscript H denotes conjugate transpose, and σ 2   is the noise power. The frequency power spectra Ψ f   can be expressed as

To mitigate the impact of high-frequency attenuation and enhance the outcomes of high-frequency components within Ψ k , f , the CWSF (Beall et al., 1982) is employed to derive the conditional wavenumber spectral density Ψ ( k | f ) from Equations (5) and Equations (6). Ψ ( k | f ) can be expressed as

where ω v = | sin ( π M v ) Q sin ( π v ) | is a periodic function with period 1. It is symmetric about v = 0 , and has its maximum value when   v = 0 . The first zero points of ω v   are v 0 = ± 1 / M , hence the mainlobe width is 2 / M .

In this paper, we only consider ω v for a single cycle. From Equation (7), Ψ ( k | f ) has its maximum value when v = f d cos θ c − k Q = 0 , and mainlobe width 2 Q / M . Figure 2 gives Ψ ( k | f ) for a broadband signal. This shows that the target energy is concentrated in the mainlobe, and the peak-energy points in the mainlobe are on a straight line r passing through the origin. The slope   ε   of the line   r can be expressed as ε = k f = Q d cos θ c . Therefore, ε is a linear function of   cos θ .

As mentioned above, θ   can be accurately estimated using the slope of the line r in Ψ ( k | f ) , thus realizing high-resolution DOA estimation. However, there will be sidelobels and perhaps outliers in Ψ ( k | f ) due to the windowing effect of the DFT and the random fluctuation noise which makes determining r difficult. To solve these issues, morphological grayscale reconstruction is used to extract regional maxima in Ψ ( k | f ) and obtain the reconstructed matrix Ψ new ( k | f ) (Vincent, 1992; Vincent, 1993). This method utilizes erosion and dilation operations based on a structuring element to reconstruct or eliminate specific regions in an image, which removes most of the outliers and significantly reduces the sidelobes.

The generalized Radon transform (GRT) is commonly used to extract information from images (Radon, 1986; Hansen and Toft, 1996; Ramm and Katsevich, 2020; Sun et al., 2021). However, it only considers amplitude information and ignores gradient information, which is not ideal. Image gradients provide the rates and directions of change for the pixels, which is useful information. Therefore, we propose a method based on the GWRT to achieve the integration of a multivariate function over a given path. Compared with the GRT, the GWRT makes full use of the image gradient information as its weights, thus providing better results. For a two-dimensional Euclidean space, the GWRT can be defined for a continuous image as (Alpatov et al., 2015)

where p is a vector containing the parameters of the line, and   c p ,   x , y , and   ϕ x , y are, respectively, a known line, space coordinates, and the intensity of points on a line in the two-dimensional image   ϕ . ∇ ,   〈 〉, and n → in Equation (8) are, respectively, the gradient operator, scalar product operation, and the unit normal vector which is perpendicular to the line c p . As the transform will be obtained from discrete-time data, the GWRT in discrete form is used rather than the integral form. The GWRT of Ψ new ( k | f ) in discrete form can be expressed as

where   f min   and   f max  are the lower and upper limits of frequency employed, respectively, with   f max  ≤ c / 2 d . Equation (9) converts Ψ new ( k | f ) into a one-dimensional matrix U weighted θ which reflects the energy distribution related to the parameter   θ . The position of the maximum value of U weighted θ is the DOA estimation.

As mentioned above, the slope of the peak-energy line in Ψ ( k | f ) is a function of the signal direction θ . Therefore, the lines in Ψ ( k | f ) must be discernible so targets from different directions can be distinguished. In other words, the difference in coordinates on the k-axis for   f max  must be equal to or greater than the mainlobe width of Ψ ( k | f ) which means

where   cos θ 1   and  cos θ 2   represent the directions of two targets. Then, the difference in directions should satisfy

Equation (11) is a function of   f max , so this frequency should be large to obtain high-resolution performance. Note that the mainlobe of Ψ new ( k | f ) will be much narrower than that of Ψ ( k | f ) due to the morphological grayscale reconstruction operation. Therefore, the resolution of the GWRT will be less than 2 c M f max d , which confirms that the proposed method has high-resolution performance.

The steps of the proposed method are as follows.

The performance of the proposed method is evaluated for a scalar towed array measurement system. Consider a line array of 32 receivers uniformly spaced at a distance d=0.25 m. The proposed method is compared with three commonly used DOA estimation methods, namely CBF, STMV, and DCV.

To evaluate the performance of the proposed method in practical applications, DOA estimation results were obtained for a towed line array with 32 elements uniformly spaced at 2 m. The experiments were conducted in Huanghai Sea, China in the summer of 2014. The water depth was approximately 40 m and the towed array was about 20 m above the sea floor. The recorded data suggests that the received signals include signals from passing vessels, experimental vessels, and the broadband pulses transmitted by the experimental vessels. The sampling frequency was equal to 8 kHz, with a total of 120 data frames, each comprising 4096 samples. The data within the frequency range of 1500 to 3000 Hz was processed using the CBF, STMV, DCV, and GWRT methods and the bearing time records (BTRs) are given in Figure 8 . Additionally, the low-frequency analysis and recording (LOFAR) for the data from one element is shown in Figure 9 .

The BTRs for the CBF method in Figure 8A show that due to the Rayleigh limit and the ambient ocean noise, the mainlobe width for each target is relatively wide so it is not possible to distinguish the targets located in the directions around  cos θ = 0.78   and cosθ=0.88. Figure 8B indicates that the STMV method has a narrower mainwidth but still fails to distinguish the two targets. The BTRs for the DCV and GWRT methods in Figures 8C, D , respectively, have much clearer backgrounds than with the CBF and STMV methods. The target trajectories are clearly distinguishable with a much narrower mainlobe width for each target. Figure 8 also shows a set of broadband pulse signals in the direction of around cos θ = 0.30 . There are two clear focused points in the red circle in Figure 8D which are not as well distinguished by other methods. This indicates that the GWRT method has a lower background noise level and thus better weak target detection and anti-interference capability.

Figure 10 gives the BTRs for the four methods at 60 s. These results indicate that the target located in the direction of around cos θ = 0.78 cannot be distinguished by the CBF and STMV methods due to the strong inference from the target located in the direction of around cos θ = 0.88 . Conversely, both the DCV and GWRT methods clearly distinguish these targets. The GWRT method has a mainlobe width similar to that of the DCV method but the sidelobe levels are lower. Thus, it is better able to suppress the interference due to strong targets and noise which makes it easier to detect weak targets. Therefore, the proposed method provides better high-resolution performance in multiple target environments.

A generalized weighted Radon transform to estimate the DOA for broadband targets was proposed. The GWRT was used on the conditional wavenumber spectrum density to convert image information in the f-k domain to the azimuth-energy domain for high-resolution DOA estimation. Simulation and sea-trial results were presented which show that the proposed method is simple and effective and does not require a priori information. It is not sensitive to the outliers and thus provides good robustness even with position errors. Furthermore, it produces a narrow mainlobe with low sidelobe levels which results in good performance when there are multiple targets and the target SNR is low. However, the proposed method is only applicable for broadband signals and it is not suitable for real-time applications. Therefore, a short-time model for DOA estimation with narrow-band signals will be considered as future work.

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

ML: Conceptualization, Validation, Writing – original draft, Methodology. DS: Supervision, Writing – review & editing, Funding acquisition, Project administration. TG: Supervision, Writing – review & editing, Validation. YL: Supervision, Writing – review & editing, Data curation, Validation. JM: Data curation, Supervision, Validation, Writing – review & editing.