DOA Estimation of Underwater Targets via Improved Monopulse Methodwith Sonar Array
Jiani Wu, Changchun Bao
College of Meteorology and Oceanography, National University of Defense Technology, ChinaAbstract
In this paper, a DOA estimation method using improved monopulse to resolve two targets in a single beam
is proposed. We first establish an echo model of two targets with sonar array. Then we derive the maximum
likelihood angle estimation based on a numerical optimization solution (Levenberg-Marquardt method). Afterthat an improved monopulse method is proposed to estimate the DOA of the targets according to the maximum
likelihood estimation principle. Finally, the simulation results prove that, the improved monopulse method
performs very well in many aspects, including small estimation error and less computational complexity.
Keywords: DOA estimation; two targets; sonar array; Monopulse

1INTRODUCTION

In the field of sonar detection, the most commonly used method for DOA estimation of underwater targets is theconventional beamforming (CBF) algorithm. It has the advantages in small computation, easy implementation,high robustness and wide application range. However, due to the limitation of Rayleigh limit, the method
cannot resolve multiple targets located within one beam. Based on array signal processing, the spatial spectrum
analysis can effectively achieve the multi-target super-resolution, and the MUSIC algorithm[1] and the ESPRIT
algorithm[2] are the most commonly used methods. Based on these two methods, scholars have proposed
many improved methods, such as root-MUSC[3], LS-ESPRIT[4], etc.. These algorithms have high accuracy andgood robustness, but require that the signal sources are narrow-band and uncorrelated. When the conditions
are not met, the performance of these algorithms will drop rapidly and sometimes even fail. In the 1980s, a
new class of DOA estimation algorithms called subspace fitting algorithms, including the maximum likelihood(ML) algorithm[5] and the weighted subspace fitting (WSF) algorithm[6], was proposed. The subspace fitting
algorithm has excellent estimation performance. But it has large computational complexity and is inconvenient
for real-time engineering processing.
In order to reduce the amount of calculation, this paper studies DOA estimation of two targets in a single
beam with sonar array. We first establish an echo model of two targets received by the array, and develop
the maximum likelihood angle estimation based on Levenberg-Marquardt(LM) method. Since it has a slow
convergence rate, then we propose a dual-monopulse method to estimate the DOA of the targets according to
the maximum likelihood estimation principle. The method is expected to gain advantages in DOA estimation
accuracy and moderate requirement on the system complexity as it is based on the dual-monopulse system. The
core idea the proposed method is that, two monopulse systems are formed and each of them points at oneof the targets, and iterations are then carried out to reduce the DOA estimation errors. In the Section 5, theperformance of this method is demonstrated via Monte Carlo simulations.

2TWO-TARGET ECHO MODEL OF SONAR ARRAY

Consider the uniform linear array consisting of N elements, the single snapshot received by the sonar array inthe case of two targets isx = A1 s (u1 ) + A2 s (u2 ) + n(1)3185

wheres (ui ) = [e− j(N−1)dπuiλ,e− j(N−3)dπuiλ,...,ej(N−3)dπuiλ,ej(N−1)dπui Tλ]is the array steering vector,ui = sin θi , i = 1, 2, θ1 and θ2 are the angles of the two targets with respect to thenormal of the array. λ is the wavelength, d is the distance between adjacent array elements. A1 , A2 ∈ C1×1 are
the complex amplitudes of the two target echoes, respectively. n ∈ CN×1 represents the thermal noise received
by array, which is independently and identically distributed white Gaussian noise with zero mean and a varianceof σ 2 .

3MAXIMUM LIKELIHOOD DOA ESTIMATION

According to the distribution of the noise in (1), one can conclude in the following Gaussian density functionof the array measurement with respect to A1 , A2 , u1 , u2 ,p(xx; u, a ) = [1/(πσ 2 )]N exp{−[1/σ 2 ](xx − S a )H (xx − S a )}(2)

where a = [ A1 A2 ], S = [ s (u1 ) s (u2 ) ]. The maximum likelihood estimates can then be obtained by maximizing the above formula, which equals to minimizing (xx − S a )H (xx − S a ). For notationalh convenience, we nor- i∆malize x as √1 x = √1 (A1 s (u1 ) + A2 s (u2 ) + n ) = A1 s 1 +A2 s 2 +nn0 , in which s 1 s 2 = √1N s (u1 ) √1N s (u2 )NN10√= S and n = N n . Denote Q = (xx − S a )H (xx − S a ) , then A1 , A2 , u1 , u2 can be estimated by setting the partial∆differentiations of Q with respect to them to 0, i.e.,H∂ Q/∂ A1 = −ssH
1 x + A1 + A2 s 1 s 2 = 0(3)H∂ Q/∂ A2 = −ssH
2 x + A1 s 2 s 1 + A2 = 0

(4)

0H
H0H
∂ Q/∂ u1 = 2Re[−AH1 s 1 x + A1 A2 s 1 s 2 ] = 0

(5)

0H
H0H
∂ Q/∂ u2 = 2Re[−AH2 s 2 x + A2 A1 s 2 s 1 ] = 0(6)which conclude in the following estimates according to (3) and (4),Â1 =H H
sH
1 x − s2 xs1 s2
21 − ||ssH
1 s 2 ||(7)Â2 =H HsH
2 x − s1 xs2 s1
21 − ||ssH
2 s 1 ||(8)However, it is difficult to obtain the expressions of û1 , û2 straightforwardly based on these equations[7]. Letω = (u1 , u2 ), and consider Q = kxx − S̄S a k2 . minQ(ω) becomes a least square problem. So we make use of anumerical optimization solution to obtain the value of û1 , û2 .
ω) =The Levenberg-Marquardt method[8] is used to solve this problem. For the least square problem minQ(ω0 ω
1 F ω 2T F (ωωωωωωFωkF(ω)k,theLM
methoddefines
that:g(ω)=∇Q(ω)=J(ω)).
WhereJ(ω)=F(ω)=[∇F(ω),...,
12TJTF
T−1FωddJIJJJωF∇F m (ω )] . And the search direction k is defined as k = −(J k k + µk )kk , where k = (ω k ),k =F (ωω k ),and µk > 0. According to the iteration rule ω k+1 = ω k + αk d k , then we derive the iteration sequenceω k }. The last value is the minimum point of optimization problem. The step factor αk is determined by the{ωArmijo search method[9].

3186We propose LM-based maximum likelihood angle estimation(LM-ML)[10].The process of the LM-ML is showedin the following:
Step one: choose the initialization value of ω as ω 0 , set the permissible error 0 ≤ ε ≤ 1, the maximum numberof iterations K and µ0 > 0, let k = 0;
ω k ) = J k T F k , when kgg(ωω k )k < ε or k ≥ K, stop and output the ω k as theStep two: calculate J k and get g (ωoptimal solution;
Step three: calculate d k , and get d k = −(JJ k T J k + µk I )−1 J k T F k ;Step four: obtain the step factor αk by the Armijo search method;
Step five: set ω k+1 = ω k + αk d k , k = k + 1. Update µk and substitute ω k+1 into Eq.(7) and (8) to update A1 , A2 ,then go to step two.
After several iterations, we obtain the minimum point ω k = (u1k , u2k ), which is the DOA estimation of theLM-ML.

4DUAL-MONOPULSE BASED DOA ESTIMATION METHOD

Since the LM-ML has a slow convergence rate and a long computation time, we propose a dual-monopulse
method for angle estimation following the guidelines of the above derivations. In the method, two monopulse
systems are constructed and the angles of the two unresolved targets are estimated via iterations[11]. It can beconcluded from Eq. (7) and (8) that,
Â1 =

H H
Ps2 ⊥ s 1 )H x
(PsH
1 x − s2 xs1 s2== wΣ1 x21 − ||ssH
Ps2 ⊥ s 1 k2kP
1 s 2 ||(9)Â2 =H H
Ps1 ⊥ s 2 )H x
(PsH
2 x − s1 xs2 s1== w Σ2 x
21 − ||ssH
Ps1 ⊥ s 2 k2kP
2 s 1 ||

(10)

where P s1 ⊥ = I − s 1 s 1 H and P s2 ⊥ = I − s 2 s 2 H are the orthogonal projection matrices with respect to s 1 and s 2 ,respectively, which satisfy
Ps1 ⊥ s 2 k2 = kP
Ps2 ⊥ s 1 k2 = 1 − |ss2 H s 1 |2kP

(11)

Moreover,w Σ1 =

Ps2 ⊥ s1 )H(P,Ps2 ⊥ s 1 k2kPw Σ2 =Ps1 ⊥ s2 )H(PP s1 ⊥ s 2 k 2kP(12)It can be easily obtained that w Σ1 H s 1 = 1, w Σ1 H s 2 = 0 and w Σ2 H s 1 = 0, w Σ2 H s 2 = 1 .
Thus Â1 , Â2 can be deemed as the sum beams of the monopulse system I and II, respectively, which are denotedas Σ1 , Σ2 in the rest of the paper. wΣ1 , wΣ2 are the weight vectors of Σ1 , Σ2 . The formulation of the sum beams
helps to enhance the target signal of interest and suppress the signal of the other target at the same time.
According Eq.(5) and (6), one can obtain that,H

H

H∆1 = s 01 x − A2 s 01 s 2 = s 01 (xx − A2 s 2 )

(13)

HHH
∆2 = s02 x − A1 s02 s1 = s02 (xx − A1 s1 )(14)

Substituting Eq.(9) and (10) into Eq.(13) and (14) respectively yields

3187

H
∆1 = s 01

H∆2 = s02x−

Ps1 ⊥ s 2 )H x
(Px−

Ps2 ⊥ s 1 )H x
(PPs1 ⊥ s 2 k2kP

Ps2 ⊥ s 1 k2kP!
s2= w ∆1 x

(15)

= w ∆2 x(16)!
s1whereHw ∆1 = s 01

Hw ∆2 = s 02I−

Ps1 ⊥ s 2 )Hs 2 (PI−Ps2 ⊥ s 1 )Hs 1 (P!
H= s 01 −2

Ps1 ⊥ s2 kkP!
H= s 02 −Ps2 ⊥ s 1 k2kP

s01Hs2
Ps1 ⊥ s 2 )H(P2PskP s1 ⊥ 2 k(17)s02Hs1Ps2 ⊥ s 1 )H(PPs2 ⊥ s 1 k2kP

(18)

It can be easily concluded that w ∆1 H s 1 = 0 , w ∆1 H s 2 = 0 , w ∆2 H s 2 = 0 , w ∆2 H s 1 = 0Thus ∆1 , ∆2 can be deemed as the difference beams of monopulse system I and II respectively, and w ∆1 , w ∆2 are
the weight vectors of ∆1 , ∆2 . The difference beams of the two systems help to suppress the jamming signal.
The iterative formulas of monopulse system I and II are  −1 ∆1−1 ∆2
u1t = u10 + k, u2t = u20 + k(19)Σ1Σ2where k is the slope of the monopulse response curve (MRC) and its value is derived according to literature[12] when the angle estimation error is very small. The formula for setting the value of MRC isN−12d 2k = [ss (u) s (u) + s (u) s (u)]/ss (u)ss(u) = 8 ∑πn( ) /NλN−100

HH 00

H(20)

n=− 2According to the above derivations and Eq.(19), the dual-monopulse system for angle estimation can be constructed as that shown in Fig.1.
According to Eq.(19), this system estimates the directions of the two unresolved targets based on the current
beam pointing orientations. Then the beam orientation is updated according to the refreshed DOA estimates.
The iteration is terminated when a predefined convergence criterion is satisfied. In this paper, the convergencecriterion is set when the difference between the new DOA estimates and the beam pointing orientations is lessthan a certain iterative threshold in a single cycle. In Fig.1, u10 , u20 are the sine values of the current beam
pointing angles and u1t , u2t are the sine values of the DOA estimates of the two targets, γ is the iterativethreshold.

5SIMULATION

Consider a vertical uniform linear array whose inter-element distance equals half-wavelength. The number of
λarray elements is N=100, then the beamwidth is θ3dB = 0.886 Nd≈ 1.016◦ . The iterative threshold γ = 10−8 .
SNR is defined as the ratio of the direct signal power to the noiseq variance of each array element, i.e., SNRi =|Ai |2. The root mean square error (RMSE) defined as RMSEi =σ23188E[(θ̂i − θi )2 ] is introduced for DOA estimation

Figure 1. The process of the dual-monopulse based DOA estimation method

Figure 2. Iterations of dual-monopulse method and LM-ML for different SNR
precision evaluation. In order to enhance the generality of the conclusions, the target azimuth, the included
angle and the RMSE are normalized with respect to the beam width, the unitary azimuth of the target isθ̄i = θi /θ3dB , the unitary inter-target angle is ∆θ = |θ1 − θ2 |/θ3dB , the unitary root mean square error (URMSE)
is URMSE =RMSE/θ3dB .
The influences of the SNR and the azimuth on the performance of the proposed method are studied in thefollowing simulations. When the influence of a certain factor on the performance is studied, the other parameters
are fixed on typical values. In each situation, 1000 Monte Carlo simulations are carried out to obtain thestatistical performances.
5.1 Relationship between the Number of Iterations and the SNR
We compare the number of iterations of the dual-monopulse method with that of LM-ML algorithm. In thisexperiment, we assume the radian of the two target azimuths to be θ1 = −0.007, θ2 = 0.008 , respectively, andthe corresponding unitary azimuths are θ̄1 = −0.41, θ̄2 = 0.47. Fig.2 shows the relationship curve between theiterations of two methods and the SNR with SNR = SNR1 = SNR2 .
Fig.2 shows that the number of iterations of the dual-monopulse method is less than that of the LM-ML method
and it decreases with increasing SNR.

3189

Figure 3. Iterations of dual-monopulse method and LM-ML for different unitary inter-target angle

Figure 4. URMSE1 derived from two methods for different SNR
5.2 Relationship between the Number of Iterations and the Inter-Target Angle
In this experiment, we assume the SNR of the two signals as SNR1 = SNR2 = 30dB, Fig.3 shows the relationshipcurve between the iterations of two methods and the unitary inter-target angle ∆θ .
Fig.3 shows that the number of iterations of the dual-monopulse method decreases significantly when the intertarget angle increases. And when the inter-target angle is smaller than 0.7 times the beamwidth, the LM-ML
method cannot reach the convergence condition after 3000 iterations. From the above results one can conclude
that, compared with the LM-ML method, the dual-monopulse method converges fast and has low computationalcomplexity. So it can estimate DOA of the targets more quickly and efficiently.
5.3 Relationship between the DOA Estimation Accuracy and the SNR
We compare the DOA estimation accuracy of the dual-monopulse method with that of the LM-ML method. In
this experiment, we assume the radian of the two target azimuths to be θ1 = −0.007, θ2 = 0.008 , respectively.
Fig.4 shows the relationship curve between the URMSE1 of the first target of two methods and the SNR with
SNR = SNR1 = SNR2 .
Fig.4 shows that the RMSE decreases and the angle estimation accuracy improves with increasing SNR. The3190Figure 5. URMSE1 derived from two methods for different unitary inter-target angle
angle estimation error of the dual-monopulse method is close to that of the LM-ML method.
5.4 Relationship between the DOA Estimation Accuracy and Inter-Target Angle
In this experiment, we assume the SNR of the two signals as SNR1 = SNR2 = 30dB, Fig.5 shows the relationshipcurve between the URMSE1 of two methods and the unitary inter-target angle ∆θ .
Fig.5 shows that the URMSE increases when the inter-target angle decreases. The phenomenon accords with theexperience of most sonar users that, when the inter-targer angle decreases, the two targets get closer and become
more difficult to be resolved, and the precision of DOA estimation decreases accordingly. Moreover, the angle
estimation error of the dual-monopulse method is close to that of the LM-ML method, which approximates
CRLB[10].

6CONCLUSIONSThis paper studied the problem of resolving two unresolved targets with sonar array, and a dual-monopulse
based DOA estimation method is proposed. The simulation results indicate that, this method works effectively,
and it gains advantages in small estimation error. Moreover, compared with the maximum likelihood angle
estimation algorithm, the proposed method is faster and less computationally intensive.

REFERENCES[1] R. O. Schmidt. Multiple emitter location and signal parameter estimation[J]. IEEE Transactions on Antennasand Propagation, 1986, 34(3):276-280.
[2] R. Roy, T. Kailath. ESPRIT: A subspace rotation approach to estimation of parameters of cissoids innoise[J]. IEEE Transactions on Acoustics Speech, and Signal Processing, 1986,34(10): 1340-1342.
[3] B. D. Rao, K. V. Hari. Performance analysis of Root-MUSIC[J].IEEE Transactions on Signal Processing,
1989, 37(12):1939-1949.
[4] R. Roy, T. Kailath. ESPRIT estimation of signal parameter via rotational invariance techniques[J]. IEEE
Transactions on Signal Processing,1989, 37(7):984-995.
[5] B. Ottersten, M. Viberg, P. Stoica. Exact and large sample ML techniques for parameter estimation and3191

detection in array processing[J]. In Haykin,Litva,and Shepherd, Editor, Radar array processing ,SpringerVerlag,1993,99-151.
[6] J. A. Cadzow. A high resolution direction-of-arrival algorithm for narrow-band coherent an incoherent
sources[J]. IEEE Transaction on Acoustics Speech, and Signal Processing, 1988, 36(7):965-979.
[7] W. D. White, Low-angle radar tracking in the presence of multipath. IEEE Transactions on Aerospace andElectronic Systems, 10(6), 1974, 835-852[8] N. Yamashita, M. Fukushima, On the Rate of Convergence of the Levenberg-Marquardt Method[J]. Computing, 2001, 15: 239-246.
[9] Ma Changfeng, Optimization method and its Matlab programming. Beijing: Science Publisher, 2010[10] J. N. WU, Y. G. Chen. Maximum Likelihood Angle Estimation of Two Targets within a Single Array
Radar Beam[J]. Journal of national university of defense technology, 2016, 6.
[11] J. N. WU, Z. Y. Xiong. Resolution of multiple unresolved targets via dual monopulse with arrayradar[C].2014 11th European Radar Conference,Rome, Italy, Oct. 2014.
[12] Nickel U., ”Overview of Generalized Monopulse Estimation”, IEEE A&E Systems Magazine, June 2006,
21(6):27-36

3192