Introduction

In a sonar system, one is often searching for targets, e.g. submarines, mines, or fish. The sonar system gathers sounds from its acoustic sensors searching for either echoes or sounds emitted by the target. These received sounds are the observations, and for active sonar, are collected into sets of observations related to the time of transmission of the sonar waveform. The sounds related to the broadcast of a single transmission is called a ping history. These pings occur sequentially in time, so one naturally has a sequence of observations (sound recordings) indexed by the time that a ping was transmitted.

The sonar system decision space includes hypotheses about the target’s presence, location, velocity and classification. The observations YkYk size 12{Y rSub { size 8{k} } } {} at ping k contain information about the sonar decision space, but are also influenced, and often dominated by other sounds, such as noise and reverberation. Echoes, noise and reverberation are significantly influenced by the propagation properties of the ocean. These environmental effects are important when making useful inferences about the target echoes that may or may not be present in the sonar ping history.

In most sonar systems today, environmental interference effects, such as noise and reverberation, are treated as random variables. The sonar processing designer develops algorithms that make detections and estimates target states by assuming a statistical model of the echo and interference, choosing environmental interference model parameters (amplitude, covariance, autocorrelation, etc.) and then computing a detection decision or state estimate.

The environmental effects are usually estimated as part of the target state decision process, or the processing algorithm is constructed to be invariant to the environmental effects.

The primary detection processing method for current active sonars is to process the ping history with a bank of matched filters. The filters are constructed so that each filter is constructing the cross-correlation between the transmitted waveform and the pre-whitened ping history.

A monostatic sonar has the source and receiver in the same location, and hence the receiver cannot realistically capture the ping history until the waveform transmission is complete. To simplify the problem, we will look for targets that are far away from the sonar, so that the echo reception occurs when the reverberation level has fallen below the background noise level. In this way, we are dealing with target echoes that are essentially embedded in background noise only.

Decision Space

We use the symbol φφ size 12{φ} {} to designate the target absent hypothesis. The other hypotheses concern the location of a single target. We then use as a decision space the composite space D=φ∪HτD=φ∪Hτ size 12{D=φ union H"" lSup { size 8{τ} } } {}where HτHτ size 12{H"" lSup { size 8{τ} } } {}is the space constructed from all target present hypotheses. The target hypothesis space consists of those locations around the sonar that generate echoes embedded in noise only. These location hypothesis hh size 12{h} {} are part of the decision space h=(x,y)∈Hτh=(x,y)∈Hτ size 12{h= \( x,y \) in H rSup { size 8{τ} } } {}specified by the set of ordered pairs h=(x,y)h=(x,y) size 12{h= \( x,y \) } {} such that

h=(x,y)∈Hτh=(x,y)∈Hτ size 12{h= \( x,y \) in H rSup { size 8{τ} } } {}if Rmin<(x−xr)2+(y−yr)2<RmaxRmin<(x−xr)2+(y−yr)2<Rmax size 12{R rSub { size 8{"min"} } < sqrt { \( x - x rSub { size 8{r} } \) rSup { size 8{2} } + \( y - y rSub { size 8{r} } \) rSup { size 8{2} } } <R rSub { size 8{"max"} } } {}

Where we are assuming that the depth of the target is small when compared to its (x,y)(x,y) size 12{ \( x,y \) } {} coordinates, the receiver is located at (xr,yr)(xr,yr) size 12{ \( x rSub { size 8{r} } ,y rSub { size 8{r} } \) } {}. RminRmin size 12{R rSub { size 8{"min"} } } {}is the range at which the echo is noise, not reverberation limited, and RmaxRmax size 12{R rSub { size 8{"max"} } } {}is the farthest range of interest. For this problem, hh size 12{h} {} is an index into the target range from the sonar.

Sonar Receiver Model

The sonar transmits the waveform m(t)m(t) size 12{m \( t \) } {}for each ping. In most sonar transmitters, the transmitted waveform is narrow-band, that is, the waveform bandwidth is much smaller than its center frequency, ff size 12{f} {}. This is true because efficient sonar transmitters use resonant mechanical and electrical components to provide maximum electrical to sound power transfer. An approximation therefore is to model the transmitted waveform as an amplitude modulated carrier:

m(t)=sin(2πft)w(t)m(t)=sin(2πft)w(t) size 12{m \( t \) ="sin" \( 2π ital "ft" \) w \( t \) } {}, t=(0,T)t=(0,T) size 12{t= \( 0,T \) } {}

We will assume that the target is motionless, so that Doppler effects can be ignored. We will assume that the sonar receiver is a single sensor, with no directionality characteristics. For each target location hypothesis h=(x,y)h=(x,y) size 12{h= \( x,y \) } {} we know approximately the received echo time series:

g ( t ∣ h ) = Bm ( t − 2R / c ) g ( t ∣ h ) = Bm ( t − 2R / c ) size 12{g \( t \lline h \) = ital "Bm" \( t - 2R/c \) } {}

The amplitude BB size 12{B} {} is related to the propagation loss out to the target hypothesis location, and the reflection characteristics of the target. The time delay 2R/c2R/c size 12{2R/c} {} corresponds to the time it takes for the transmission waveform to reach the target and return to the sonar. RR size 12{R} {} is the range to the target and c is the effective speed of sound, when including refraction and boundary reflections.

The received echo is band-limited to approximately the same frequency band as the transmission. The receiver bandwidth may be greater than the transmitted bandwidth due to Doppler frequency shifts, but for the present, we are assuming that the target is not moving. Sonar receivers use heterodyne techniques to reduce the data storage of the ping history. The sonar receiver multiplies the ping history by a carrier signal e−j2πfte−j2πft size 12{e rSup { size 8{ - j2π ital "ft"} } } {}to shift the positive frequency part of the received echo closer to DC. The resulting signal is then low pass filtered to eliminate the shifted negative frequency part of the ping history. Since the original ping history was real, the negative frequency part of the signal spectra carries no additional information. The result is a complex signal with a lower bandwidth, but retains all of the echo related information of the original ping history. This heterodyne process can be done in the analog or digital domain.

A target echo passing through the heterodyne part of the sonar receiver becomes:

r ( t ∣ h ) = Ae jθ w ( t − 2R / c ) r ( t ∣ h ) = Ae jθ w ( t − 2R / c ) size 12{r \( t \lline h \) = ital "Ae" rSup { size 8{jθ} } w \( t - 2R/c \) } {}

The phase shift θθ size 12{θ} {} corresponds to the phase shift due to heterodyne operation; the uncertainty in propagation conditions; and the summation of multi-path arrivals with almost the same time delay, etc.

We will assume that the target echo amplitude, AejθAejθ size 12{ ital "Ae" rSup { size 8{jθ} } } {},is a complex Gaussian random variable with zero mean and with standard deviation σ2(h).σ2(h). size 12{σ rSup { size 8{2} } \( h \) "." } {}We are modeling the echo as having the same waveform as the transmission, but with an uncertain phase and amplitude. This is assuming that the target echo amplitude obeys Swerling target type I statistics with unknown phase.

While receiving an echo, we will also receive ambient noise, q(t)q(t) size 12{q \( t \) } {}, which we will assume to be complex Gaussian noise, with constant power spectral density over the receiver’s bandwidth. The noise power spectral density over the receiver’s bandwidth BW is assumed to be N0N0 size 12{N rSub { size 8{0} } } {}Pascals^2/Hertz.

We receive L complex valued samples, yy size 12{y} {} as the ping history after heterodyning. For hypothesis h, the observation in discrete time is:

y(kΔt)=q(kΔt)+r(kΔt−2R/c∣h)y(kΔt)=q(kΔt)+r(kΔt−2R/c∣h) size 12{y \( kΔt \) =q \( kΔt \) +r \( kΔt - 2R/c \lline h \) } {} for , k=1,…Lk=1,…L size 12{k=1, dotslow L} {}

where ΔtΔt size 12{Δt} {} is the digital sample rate after heterodyning, and q is a sample of the noise and reverberation. Note that r(kΔt−2R(h)/c∣h)=0r(kΔt−2R(h)/c∣h)=0 size 12{r \( kΔt - 2R \( h \) /c \lline h \) =0} {} when 2R(h)/c>kΔt,2R(h)/c>kΔt, size 12{2R \( h \) /c>kΔt,} {} because the echo is delayed. The delay for hypothesis hh size 12{h} {} in samples, is given by

D ( h ) = 2R cΔt D ( h ) = 2R cΔt size 12{D \( h \) = left [ { {2R} over {cΔt} } right ]} {}

where [x] is the nearest integer to x. We choose the sample rate ΔtΔt size 12{Δt} {} to be small enough to satisfy the Nyquist sampling criteria for the received echo. We will assume that the non-zero part of the echo is N samples long.

Active Sonar Likelihood Function

We will find that the natural logarithm of the measurement likelihood ratio simplifies the detection expression considerably:

log L ( y ∣ h ) L ( y ∣ φ ) = log ( det ( C ) ) − log ( det ( C + C r ) ) + y H C − 1 − ( C r + C ) − 1 y log L ( y ∣ h ) L ( y ∣ φ ) = log ( det ( C ) ) − log ( det ( C + C r ) ) + y H C − 1 − ( C r + C ) − 1 y size 12{"log" { {L \( y \lline h \) } over {L \( y \lline φ \) } } ="log" \( "det" \( C \) \) - "log" \( "det" \( C+C rSub { size 8{r} } \) \) +y rSup { size 8{H} } left (C rSup { size 8{ - 1} } - \( C rSub { size 8{r} } +C \) rSup { size 8{ - 1} } right )y} {}

The determinant of CC size 12{C} {}is det(C)=det(N0IL−Dmin)det(R)=(N0)L−Dmindet(R)det(C)=det(N0IL−Dmin)det(R)=(N0)L−Dmindet(R) size 12{"det" \( C \) ="det" \( N rSub { size 8{0} } I rSub { size 8{L - D rSub { size 6{"min"} } } } \) "det" \( R \) = \( N rSub {0} size 12{ \) rSup {L - D rSub { size 6{"min"} } } } size 12{"det" \( R \) }} {}

It is convenient to partition C+CrC+Cr size 12{C+C rSub { size 8{r} } } {}into sub-matrices compatible with r.

C + C r = R 0 0 0 0 N 0 I D h − D min 0 0 0 0 N 0 I N + σ A ( h ) 2 ww H 0 0 0 0 N 0 I L − D ( h ) − N C + C r = R 0 0 0 0 N 0 I D h − D min 0 0 0 0 N 0 I N + σ A ( h ) 2 ww H 0 0 0 0 N 0 I L − D ( h ) − N size 12{C+C rSub { size 8{r} } = left [ matrix { R {} # 0 {} # 0 {} # 0 {} ## 0 {} # N rSub { size 8{0} } I rSub { size 8{D left (h right ) - D rSub { size 6{"min"} } } } {} # 0 {} # 0 {} ## 0 {} # 0 {} # N rSub {0} size 12{I rSub {N} } size 12{+σ rSub {A \( h \) } rSup {2} } size 12{ bold "ww" rSup {H} } {} # size 12{0} {} ## size 12{0} {} # size 12{0} {} # size 12{0} {} # size 12{N rSub {0} size 12{I rSub {L - D \( h \) - N} }} {} } right ]} {}

The determinant of C+CrC+Cr size 12{C+C rSub { size 8{r} } } {}is

det ( C + C r ) = det ( R ) det ( N 0 I D ( h ) − D min ) det ( N 0 I N + σ A ( h ) 2 rr H ) det ( N 0 I L − D ( h ) − N ) det ( C + C r ) = det ( R ) det ( N 0 I D ( h ) − D min ) det ( N 0 I N + σ A ( h ) 2 rr H ) det ( N 0 I L − D ( h ) − N ) size 12{"det" \( C+C rSub { size 8{r} } \) ="det" \( R \) "det" \( N rSub { size 8{0} } I rSub { size 8{D \( h \) - D rSub { size 6{"min"} } } } \) "det" \( N rSub {0} size 12{I rSub {N} } size 12{+σ rSub {A \( h \) } rSup {2} } size 12{ bold "rr" rSup {H} } size 12{ \) "det" \( N rSub {0} } size 12{I rSub {L - D \( h \) - N} } size 12{ \) }} {}

Which becomes

det ( C + C r ) = ( N 0 ) L − D min − N det ( N 0 I N + σ A ( h ) 2 rr H ) det ( R ) det ( C + C r ) = ( N 0 ) L − D min − N det ( N 0 I N + σ A ( h ) 2 rr H ) det ( R ) size 12{"det" \( C+C rSub { size 8{r} } \) = \( N rSub { size 8{0} } \) rSup { size 8{L - D rSub { size 6{"min"} } - N} } "det" \( N rSub {0} size 12{I rSub {N} } size 12{+σ rSub {A \( h \) } rSup {2} } size 12{ bold "rr" rSup {H} } size 12{ \) "det" \( R \) }} {}

Using Sylvester’s determinant theorem:

det ( I + AB ) = det ( I + BA ) det ( I + AB ) = det ( I + BA ) size 12{"det" \( I + bold "AB" \) =" det" \( I + bold "BA" \) } {}

We obtain:

det ( C + C r ) = ( N 0 ) L − D min ( 1 + σ A ( h ) 2 r H r N 0 ) det ( R ) det ( C + C r ) = ( N 0 ) L − D min ( 1 + σ A ( h ) 2 r H r N 0 ) det ( R ) size 12{"det" \( C+C rSub { size 8{r} } \) = \( N rSub { size 8{0} } \) rSup { size 8{L - D rSub { size 6{"min"} } } } \( 1+ { {σ rSub {A \( h \) } rSup {2} size 12{r rSup {H} } size 12{r}} over {N rSub {0} } } size 12{ \) "det" \( R \) }} {}

Which equals

det ( C + C r ) = ( N 0 ) L − D min ( 1 + σ A ( h ) 2 N 0 ) det ( R ) det ( C + C r ) = ( N 0 ) L − D min ( 1 + σ A ( h ) 2 N 0 ) det ( R ) size 12{"det" \( C+C rSub { size 8{r} } \) = \( N rSub { size 8{0} } \) rSup { size 8{L - D rSub { size 6{"min"} } } } \( 1+ { {σ rSub {A \( h \) } rSup {2} } over { size 12{N rSub {0} } } } size 12{ \) "det" \( R \) }} {}

Now we partition the observed ping history y into

y = y R y N1 y h y N2 y = y R y N1 y h y N2 size 12{y= left [ matrix { y rSub { size 8{R} } {} ## y rSub { size 8{N1} } {} ## y rSub { size 8{h} } {} ## y rSub { size 8{N2} } } right ]} {}

so that:

C − 1 − ( C r + C ) − 1 = 0 0 0 0 0 0 0 0 0 0 1 N 0 I N − ( N 0 I N + σ A ( h ) 2 ww H ) − 1 0 0 0 0 C − 1 − ( C r + C ) − 1 = 0 0 0 0 0 0 0 0 0 0 1 N 0 I N − ( N 0 I N + σ A ( h ) 2 ww H ) − 1 0 0 0 0 size 12{C rSup { size 8{ - 1} } - \( C rSub { size 8{r} } +C \) rSup { size 8{ - 1} } = left [ matrix { 0 {} # 0 {} # 0 {} # 0 {} ## 0 {} # 0 {} # 0 {} # 0 {} ## 0 {} # 0 {} # { {1} over {N rSub { size 8{0} } } } I rSub { size 8{N} } - \( N rSub { size 8{0} } I rSub { size 8{N} } +σ rSub { size 8{A \( h \) } } rSup { size 8{2} } bold "ww" rSup { size 8{H} } \) rSup { size 8{ - 1} } {} # {} ## 0 {} # 0 {} # 0 {} # 0{} } right ]} {}

The Woodbury matrix identity states:

( A + UCV ) − 1 = A − 1 − A − 1 U ( C − 1 + VA − 1 U ) − 1 VA − 1 ( A + UCV ) − 1 = A − 1 − A − 1 U ( C − 1 + VA − 1 U ) − 1 VA − 1 size 12{ \( A+ bold "UCV" \) rSup { size 8{ - 1} } =A rSup { size 8{ - 1} } - A rSup { size 8{ - 1} } U \( C rSup { size 8{ - 1} } + bold "VA" rSup { size 8{ - 1} } U \) rSup { size 8{ - 1} } bold "VA" rSup { size 8{ - 1} } } {}

So that

( N 0 I N + σ A ( h ) 2 ww H ) − 1 = 1 N 0 I N − 1 N 0 2 ww H 1 N 0 + 1 σ A ( h ) 2 ( N 0 I N + σ A ( h ) 2 ww H ) − 1 = 1 N 0 I N − 1 N 0 2 ww H 1 N 0 + 1 σ A ( h ) 2 size 12{ \( N rSub { size 8{0} } I rSub { size 8{N} } +σ rSub { size 8{A \( h \) } } rSup { size 8{2} } bold "ww" rSup { size 8{H} } \) rSup { size 8{ - 1} } = { {1} over {N rSub { size 8{0} } } } I rSub { size 8{N} } - { {1} over {N rSub { size 8{ {} rSub { size 6{0} } } } rSup {2} } } { { size 12{ bold "ww" rSup {H} } } over { size 12{ { {1} over {N rSub {0} } } size 12{+ { {1} over {σ rSub {A \( h \) } rSup {2} } } }} } } } {}

And hence

y H C − 1 − ( C r + C ) − 1 y = 1 N 0 σ A ( h ) 2 / N 0 ( 1 + σ A ( h ) 2 N 0 ) ∣ w H y h ∣ 2 y H C − 1 − ( C r + C ) − 1 y = 1 N 0 σ A ( h ) 2 / N 0 ( 1 + σ A ( h ) 2 N 0 ) ∣ w H y h ∣ 2 size 12{y rSup { size 8{H} } left (C rSup { size 8{ - 1} } - \( C rSub { size 8{r} } +C \) rSup { size 8{ - 1} } right )y= { {1} over {N rSub { size 8{0} } } } { {σ rSub { size 8{A \( h \) } } rSup { size 8{2} } /N rSub { size 8{0} } } over { \( 1+ { {σ rSub { size 8{A \( h \) } } rSup { size 8{2} } } over {N rSub { size 8{0} } } } \) } } lline w rSup { size 8{H} } y rSub { size 8{h} } rline rSup { size 8{2} } } {}

So that the log likelihood function becomes

log L ( y ∣ h ) L ( y ∣ φ ) = − log ( 1 + σ A ( h ) 2 N 0 ) + 1 N 0 σ A ( h ) 2 / N 0 ( 1 + σ A ( h ) 2 N 0 ) ∣ w H y h ∣ 2 log L ( y ∣ h ) L ( y ∣ φ ) = − log ( 1 + σ A ( h ) 2 N 0 ) + 1 N 0 σ A ( h ) 2 / N 0 ( 1 + σ A ( h ) 2 N 0 ) ∣ w H y h ∣ 2 size 12{"log" { {L \( y \lline h \) } over {L \( y \lline φ \) } } = - "log" \( 1+ { {σ rSub { size 8{A \( h \) } } rSup { size 8{2} } } over {N rSub { size 8{0} } } } \) + { {1} over {N rSub { size 8{0} } } } { {σ rSub { size 8{A \( h \) } } rSup { size 8{2} } /N rSub { size 8{0} } } over { \( 1+ { {σ rSub { size 8{A \( h \) } } rSup { size 8{2} } } over {N rSub { size 8{0} } } } \) } } lline w rSup { size 8{H} } y rSub { size 8{h} } rline rSup { size 8{2} } } {}

Putting the noise variance inside the observation vector yields for the log likelihood function

log L ( y ∣ h ) L ( y ∣ φ ) = − log ( 1 + σ A ( h ) 2 N 0 ) + σ A ( h ) 2 / N 0 ( 1 + σ A ( h ) 2 N 0 ) ∣ w H y h N 0 ∣ 2 log L ( y ∣ h ) L ( y ∣ φ ) = − log ( 1 + σ A ( h ) 2 N 0 ) + σ A ( h ) 2 / N 0 ( 1 + σ A ( h ) 2 N 0 ) ∣ w H y h N 0 ∣ 2 size 12{"log" { {L \( y \lline h \) } over {L \( y \lline φ \) } } = - "log" \( 1+ { {σ rSub { size 8{A \( h \) } } rSup { size 8{2} } } over {N rSub { size 8{0} } } } \) + { {σ rSub { size 8{A \( h \) } } rSup { size 8{2} } /N rSub { size 8{0} } } over { \( 1+ { {σ rSub { size 8{A \( h \) } } rSup { size 8{2} } } over {N rSub { size 8{0} } } } \) } } lline w rSup { size 8{H} } { {y rSub { size 8{h} } } over { sqrt {N rSub { size 8{0} } } } } rline rSup { size 8{2} } } {}

The log-likelihood ratio shows that the assumed echo shape, wHwH size 12{w rSup { size 8{H} } } {}, the assumed energy of the signal σA(h)2σA(h)2 size 12{σ rSub { size 8{A \( h \) } } rSup { size 8{2} } } {}, the ambient noise density N0N0 size 12{N rSub { size 8{0} } } {}, and the expected location in time of the echo needed to select yhyh size 12{y rSub { size 8{h} } } {}are the parameters required to evaluate the log-likelihood function. The assumed echo shape and energy can vary by hypothesis h, but the noise properties N0N0 size 12{N rSub { size 8{0} } } {}have been assumed to be the same for each decision hypothesis h.

Signal to Noise Ratio of the Detector

The magnitude squared term, ∣wHyhN0∣2∣wHyhN0∣2 size 12{ lline w rSup { size 8{H} } { {y rSub { size 8{h} } } over { sqrt {N rSub { size 8{0} } } } } rline rSup { size 8{2} } } {}, is a matched filter, where the observations are cross-correlated with the signal template. The observations are normalized by the noise variance before cross correlation, which is a form of pre-whitening. The term σA(h)2N0σA(h)2N0 size 12{ { {σ rSub { size 8{A \( h \) } } rSup { size 8{2} } } over {N rSub { size 8{0} } } } } {} is the Energy to Noise Density Ratio (ENR) of the detection problem. Note that the signal to noise ratio of the matched filter output can be written as:

E ∣ w H Aw N 0 ∣ 2 E ∣ w H q N 0 ∣ 2 = σ A ( h ) 2 N 0 E w H qq H w N 0 = σ A ( h ) 2 N 0 w H Iw = σ A ( h ) 2 N 0 = ENR E ∣ w H Aw N 0 ∣ 2 E ∣ w H q N 0 ∣ 2 = σ A ( h ) 2 N 0 E w H qq H w N 0 = σ A ( h ) 2 N 0 w H Iw = σ A ( h ) 2 N 0 = ENR size 12{ { {E left lbrace lline w rSup { size 8{H} } { {Aw} over { sqrt {N rSub { size 8{0} } } } } rline rSup { size 8{2} } right rbrace } over {E left lbrace lline w rSup { size 8{H} } { {q} over { sqrt {N rSub { size 8{0} } } } } rline rSup { size 8{2} } right rbrace } } = { { { {σ rSub { size 8{A \( h \) } } rSup { size 8{2} } } over {N rSub { size 8{0} } } } } over {E left lbrace { {w rSup { size 8{H} } ital "qq" rSup { size 8{H} } w} over {N rSub { size 8{0} } } } right rbrace } } = { { { {σ rSub { size 8{A \( h \) } } rSup { size 8{2} } } over {N rSub { size 8{0} } } } } over {w rSup { size 8{H} } ital "Iw"} } = { {σ rSub { size 8{A \( h \) } } rSup { size 8{2} } } over {N rSub { size 8{0} } } } = ital "ENR"} {}

This result is a general result for matched filters. A matched filter’s SNR is the Energy to Noise Density Ratio of the problem. The energy of a signal being detected is related to the average amplitude and the duration of the signal. The SNR output of the matched filter is independent of the details of the waveform being detected, only the signal energy and the noise spectral density determine the matched filter response.

The likelihood function is dependent only the ENR as well.