Mark L. Moran and Donald G. Albert June 1996

CRREL REPORT 96-8

Source Location and

Tracking Capability of a

Small Seismic Array

Seismic-Acoustic Source

Mine

Seismic-Acoustic Detectors
Abstract: Recordings of seismic wavefields from various sources were obtained using a small array of

vertical-component geophones under winter conditions

at Grayling, Michigan. These data were processed

using a frequency-wavenumber domain Capon minimum variance beamformer to estimate the bearing

angle and propagation velocity of the waves emitted

from the source. The cross power matrix was adaptively

estimated using a tapered block-averaging procedure.

The wave sources were sledgehammer blows on the

ground surface, .45 caliber blank pistol shots, and an

M60 tank moving at 4.5 m s–1 along a road near the

array. Reliable wavenumber spectra were obtained for

all sources. Processing results for the hammer blows

show that the dominant seismic arrival is a Rayleigh

wave traveling at roughly 220 m s–1. For the pistol

shots, two arrivals corresponding to the airwave (338

m s–1) and the air-coupled Rayleigh waves (220 m s–

How to get copies of CRREL technical publications:

Department of Defense personnel and contractors may order reports through the Defense Technical Information Center:

DTIC-BR SUITE 0944

8725 JOHN J KINGMAN RD

FT BELVOIR VA 22060-6218

Telephone 1 800 225 3842

E-mail help@dtic.mil

msorders@dtic.mil

WWW http://www.dtic.dla.mil/

All others may order reports through the National Technical Information Service:

NTIS

5285 PORT ROYAL RD

SPRINGFIELD VA 22161

Telephone 1 703 487 4650

1 703 487 4639 (TDD for the hearing-impaired)

E-mail orders@ntis.fedworld.gov

WWW http://www.fedworld.gov/ntis/ntishome.html

A complete list of all CRREL technical publications is available from:

USACRREL (CECRL-TL)

72 LYME RD

HANOVER NH 03755-1290

Telephone 1 603 646 4338

E-mail techpubs@crrel.usace.army.mil

For information on all aspects of the Cold Regions Research and Engineering Laboratory, visit our World Wide Web site:

http://www.crrel.usace.army.mil

Cover: Seismic and acoustic arrivals from a moving tank at an intelligent mine sensor array.

1) were observed. The results for these sources were

relatively insensitive to the processing parameters used.

For the moving vehicle, the dominant signals observed

were Rayleigh waves (220 m s–1). Accurate locations

were obtained for this moving source, although the

processing parameters had to be carefully selected, and

the choice of frequency parameters affected the accuracy

of the wavenumber results. Maximizing the number of

degrees of freedom and the coherence of the frequency

estimates and minimizing the variation of the coherence

across adjacent frequency bins provided the most consistently reliable strategy for obtaining accurate

wavenumber estimates for the moving vehicle. The sensitivity of the wavenumber estimates to the frequency

processing parameters seems to be related to the bias in

the phase spectra of the signals and will potentially occur

in any bearing estimation method that uses temporal

frequncy phase spectra.
Mark L. Moran and Donald G. Albert June 1996

CRREL Report 96-8

U.S. Army Corps

of Engineers

Cold Regions Research &

Engineering Laboratory

Source Location and

Tracking Capability of a

Small Seismic Array

Prepared for

OFFICE OF THE CHIEF OF ENGINEERS

Approved for public release; distribution is unlimited.
1. INTRODUCTION

Acoustic and seismic sensor array processing

has many potential applications of interest to the

Army. Recent military conflicts have clearly

demonstrated the efficiency of smart weapons

systems. A critical element in the effectiveness of

many of these systems is their ability to discern

and track a target in the presence of background

noise. The Wide Area Mine (WAM) being developed by the U.S. Army is heavily dependent upon

acoustic and seismic sensor information for bearing determination of targets. Weapon systems such

as the WAM are designed to recognize, track, and

destroy hostile military vehicles such as armor,

mobile artillery, and heavy transport vehicles by

deploying small arrays of passive acoustic and/or

seismic sensors (microphones and geophones) to

monitor local seismic and acoustic wavefields. A

cartoon of the detection–attack scenario for the

WAM device is shown in Figure 1.

Other important areas where acoustic and seismic array processing may be applied include perimeter defense and intrusion detection. In many

instances the effectiveness of line-of-sight systems

such as radar, infrared, and motion detection systems are compromised or defeated by local topography, vegetation, and obscuring meteorological

conditions. In real-world settings, acoustic and

seismic waves readily propagate along bending

ray paths and are less affected by vegetation and

airborne obscurants. This allows coverage of zones

currently outside the detection limits of systems

that rely on line-of-sight propagation paths. Acoustic and seismic array-based systems are relatively

inexpensive and may be produced and deployed

in large numbers, providing a substantial perimeter defense or detection capability.

For small array processing systems to function

reliably, it is necessary to be able to estimate accurately the wavefield’s temporal frequency phase

spectrum, even with limited spatial and temporal

wavefield samples. When using wavenumber estimation procedures based on beamforming principles, these data constraints make it difficult to

obtain reliable estimates of the spatial correlation

matrix, which is critical to the beamformer’s performance. The spectral bias and variance are the

statistical error parameters that describe the reliability of the resulting estimated frequency and

wavenumber spectra. Bias is a measure of the

energy leakage from adjacent frequency bins that

erroneously appears at an incorrect frequency;

this phenomenon is termed leakage and is caused

by the sampling distribution in the time (or space)

domain. Variance is a measure of the stability in the

estimate and is a function of the number and

degree of independence between samples available to the estimation process.

Beamforming complications may result from

source motion, acoustic-to-seismic coupling, inhomogeneous propagation environments, multiple sources, and low signal-to-noise ratios (SNRs).

Given the nature of the moving seismic–acoustic

sources being considered and the limited quantity

of data available, it is important to recognize the

nonideality of the observed signals, so that those

properties of the signal that have the highest degree of coherency across the observation array can

be exploited fully.

The combined effect of these source-dependent

characteristics and the problems arising from

sparsely sampled short-duration signal vectors

has not been previously demonstrated on the beam

response function. In this report we show that

beamforming techniques can be reliably applied

to these types of signal vectors as long as careful

attention is paid to the spectral estimation parameters used in the Fourier transform process. We

present a few rough guidelines that produce consistently reliable wavenumber estimates from Capon maximum-likelihood beamformers.

Historically, frequency wavenumber domain

(F-K) beamforming analysis has been applied to a

variety of problems. In seismology, wavenumber

beamformers using geophone arrays are routinely

used to determine source parameters for regional

and teleseismic events as well as to identify nuclear

Source Location and

Tracking Capability of a

Small Seismic Array

MARK L. MORAN AND DONALD G. ALBERT
Warhead

Seismic


Acoustic


Source

W a v e F r o n t s

W a v e F r o n t s

Target Acquisition


and Tracking Phase

Warhead Deployment Phase

Seismic


Acoustic


Source

WAM

Warhead

Computor

Sensor

Computor

Sensor

Figure 1. WAM conception of target acquisition and engagement.

explosions. The application of the Capon maximum-likelihood beamformer to teleseismic events

using the Large Aperture Seismic Array in Montana was first reported by Capon et al. (1967). In

1969 Capon derived the beamformer and noise

suppression characteristics. Beamforming is also

applied to sonar problems where passive arrays of

hydrophones locate and track submarines and

surface vehicles (Hahn 1975). Aircraft tracking

with both acoustic and radar technologies and

stellar mapping of distant radio-wave emissions

are other prominent areas of application for

beamforming technology.

The use of seismic and acoustic arrays to track

and locate targets of military interest dates from

the Vietnam War. Beginning in the late 1960s,

military efforts focused, with limited success, on

locating artillery positions. In the later 1970s,

targets such as surface vehicles were investigated (Cress 1976). Recently, bearing estimation of

slow-moving military aircraft using acoustic and

seismic array processing technologies has been

actively investigated. Of particular note are the

efforts of Nawab et al. (1985) and Lacoss et al.

(1991), who have successfully tracked both airborne and ground targets at distances of a few

kilometers.

Several researchers have investigated the impact of nonideal source attributes on bearing estimates. De Graaf and Johnson (1985) discuss the

bias of Bartlett and Capon maximum-likelihood

beamformers as a function of SNR under the assumption that the spatial correlation matrix is

known exactly and the sources are ideal stationary

plane waves. For two or more sources they show

that the bias in the wavenumber estimate is a

function of the source separation and the degree of

coherence between the sources, the SNR, the relative source power, and the number of observation

array elements. When the spatial correlation matrix must be estimated from a finite amount of

data, De Graaf and Johnson (1985) show how the

maximum-likelihood estimator becomes increasingly unstable as the data set is reduced.

Knapp and Carter (1977) and Carter (1977) discussed the impact of source motion on bias in

bearing estimation from a two- or three-element

array using a time-delay estimation procedure.

Gragg (1990) showed that there is a time dependence in the Doppler-shifted spectra observed by

a single receiver when an acoustic source has a

velocity vector that does not pass through the

receiver location. He showed that this time dependence makes the Doppler shift difficult to predict.

Unless source motion effects can be corrected,

there will be a phase bias in the Fourier estimate of

the time series. The impact of the phase bias from

source motion on wavenumber estimates remains

unquantified.

In the following sections we discuss the

beamforming equations and explicitly define the

frequency domain estimation process used to genComputer

Computer

2
erate the cross-power matrix in the beamforming

equations. A confidence interval for the variance

in the cross-power matrix is developed. This is

followed by a discussion of the experimental setting, data acquisition procedure, and seismic array characteristics. Processing results are then given

for an impulsive acoustic source, a pure seismic

source, and a moving M60 tank. We obtain accurate wavenumber estimates from the beamformer

for all three source types. Next we show that small

variations in the parameters used to estimate the

cross-power matrix contribute to wavenumber

bias effects. These bias effects are summarized for

a variety of block lengths, block overlaps, and

window types.

2. BEAMFORMING PRINCIPLES

The discussion given here follows Capon (1969)

and Johnson (1982). First, a general expression for

an F–K beamformer power function is given. The

standard Bartlett (BT) and high-resolution Capon

maximum-likelihood (ML) beamformer equations

are obtained. Following the beamformer derivation, the estimation and conditioning procedures

of the spatial correlation matrix are given.

2.1 Beam-power response function

Assume an array of M sensors have vector locations 

R

→m in the x–y plane

 

R xy

→m = + m m

















2 2

as shown in Figure 2. Each sensor is a spatial

wavefield sample point that records over time.

Such a collection of sensor data is termed the signal

vector. A wavefield containing a single plane wave

source sampled at locations 

R

→m may be represented by a signal vector with components:

 

x t st

R u

c m n t m ko () () = + m

• 



















+

→ →

 , (1)

where k

→o is the wavenumber 

 

k fc u → → o ko = 









 2π / ,

is the unit vector direction of the source, nm is the

background noise field, and c is the scalar wave

propagation velocity. In the frequency domain,

eq 1 is

W a v e

F r o n t s

x

m

y

Rm

θ k

k0

Source

Figure 2. A generic array geometry composed of M sensors located

at vector positions 

r

Rm . The pseudo-plane wave source has a

wavenumber vector 

r

ko . The observation wavenumber is 

r

k . Dots

(•) are sensor locations.

3
Xm(f) = S(f)eiko • Rm + Nm(f), (2)

where S(f) is the ideal Fourier transform of the

signal without the noise.

In matrix notation, the signal vector is

X(f) =

X0(f)

X1(f)

X2(f)

•

•

•

XM–1(f)

 . (3)

Given M spatial samples of the wave field, a

beam response function 

[ ( , )] Yf k→ is formed by

summing over each component of the signal vector and applying a phase shift that is dependent

upon an observation wavenumber k :

Y f, k = ame ∑ –ik • Rm Xm (f)

m = 0

M – 1

 . (4)

Using eq 2 and 4,

Y f, k = am(f) S(f)e ∑ i k0 – k • Rm + Nm (f)

m = 0

M – 1

(5)

where the am(f) constants are spatial filter

weights to be defined through the application of a

specific beamforming strategy, and Nm is the

noise field modified by a phase shift of eik • Rm .

Equation 5 is equivalent to an estimation of the

complex wavenumber spectrum of the sampled

wavefield. In matrix notation the beam response is

 

Yfk,

→ * 









 = A XT , (6)

where *T means complex conjugate transpose,

and

A(f) =

a0(f)e –i k • R0

a1(f)e –i k • R1

a2(f)e –i k • R2

•

•

•

aM–1(f)e –i k • RM–1

(7)

A(f) is referred to as the steering vector. The beam

power at a given frequency and observation

wavenumber is

P f, k = A*T XX*T A (8)

or

P f, k = A*T R(f)A , (9)

where R(f) is an M × M matrix composed of the

product of the signal vector and its complex conjugate transpose. The matrix R is commonly referred

to as the spatial correlation matrix. A key result of

this report will be to show that subtle changes in

the parameters used to estimate the correlation

matrix affect the bias in the beam power function.

Equation 9 is the foundation of the wavenumber estimation procedure used in this paper. Referring to eq 5, we see that the maximum of the

beam power function occurs when k k → → 0 – = 0.

When performing a narrowband F–K estimation,

eq 9 is evaluated at a single frequency over a regularly spaced grid of observation wave-numbers.

A wideband beam power function can be produced by integrating eq 9 over a band of frequencies:

 

P k df f

→ f



 



 = ∫ A*T

1

2 RA . (10)

Recently Nawab et al. (1985) formulated more

elegant methods of computing a broadband

beamformer.

2.2. Bartlett and Capon maximum-likelihood

beamformers

The BT and ML methods are special cases of eq

9 and are a consequence of defining specific values

for the am(f) weighting functions (eq 5). The BT

method is the conventional method of estimating

an F–K spectrum. The ML method attempts to

minimize energy leakage into the beam response

from regions outside the observation wavenumber.

This energy leakage problem is the spatial equivalent of frequency domain bias. Thus, am(f) is a

component of a spatial window function analogous to time series windowing.

If all am(f) = 1 then eq 9 yields the standard BT

beamformer. Under this condition, eq 8 becomes

A′ =

e –i k • R0

e –i k • R1

e –i k • R2

•

•

•

e –i k • RM–1

. (11)

The power function is

4
PBT f,k = A *TRA′ , (12)

where R is the estimated spatial correlation matrix.

In the ML beamformer, am(f) are defined by the

properties of the spatial correlation matrix so that

spatial energy leakage is minimized. This is done

by minimizing the weighted array power function

under the constraint that a fixed gain be maintained at the observation wavenumber. Johnson

(1982) sets up the problem as

∂PML f,k

∂W = ∂∂W

W*T RW = 0 ,

subject to W*T A′ = 1 . (13)

A′ is given by eq 11, and W is defined as the

weighted steering vector whose elements are

Wm = amA′m . (14)

The solution to eq 13 is (Johnson 1982, Capon 1969)

W = R–1A′

A′

*T R–1A′

 . (15)

Note that eq 15 is a function of the spatial correlation matrix. The effect of designing the weighted

steering vector using the observed signal properties is that we produce a wavenumber window

function that is adaptive to the signal-in-noise

field and the specific array geometry. The resultant beampower function is

PML f,k = W*T RW . (16)

Substituting eq 15 into 16 and using the estimation

for the cross-power matrix, we obtain

PML f,k = 1

A′

*T R

 –1

A′

, (17)

where ^ indicates an estimate based on available

data.

This is the working form of the ML beamformer

used in this paper. It should also be noted that the

ML beamformer is also termed a minimum variance method, and in the case of only one coherent

signal it is also spatially unbiased (Capon 1969), in

the sense that the bias tends toward zero as the

number of sensors becomes large. In a multiplesource wavefield the bias is substantially smaller

for ML compared with the BT method (De Graaf

and Johnson 1985).

2.3. Estimation and conditioning of the

correlation matrix

The spatial correlation matrix (R) is constructed

from the frequency domain representation of the

signal vector x. The properties of the spatial correlation matrix determine the quality of the

beampower function. To estimate R accurately,

we use a tapered overlapped block-averaging FFT

(OBAFFT). This is a widely utilized procedure for

estimating signals with intricate spectra (Welch

1967, Capon et al. 1967). The OBAFFT method

allows one to make trade-offs between the bias

and variance in the estimate of R.

It is important to guarantee the nonsingularity

of the spatial correlation matrix when applying

the ML processing procedure. This can be done by

the block averaging process or it can be forced by

applying a small amount of additive noise to the

correlation matrix (Capon 1969). A normalization

procedure that produces ones along the diagonal

of the spatial correlation matrix is also useful in

making qualitative assessments of the accuracy of

the wavenumber spectrum.

2.3.1. Tapered overlapped block-averaged FFT

The mathematical procedure used to form the

OBAFFT is as follows. The time-domain signal

vector obtained from a passive array with M sensors and components such as that of eq 1 has the

form

x =

x0(t)

x1(t)

x2(t)

•

•

•

xM–1(t)

 . (18)

Let x have N samples at intervals of ∆t. Now

define a block of this sequence with a length L,

such that L < N; further, let each block of L samples

be overlapped by a percentage p. The number of

blocks will be

nblks = 1 +

N – L

nshft

 , (19)

where nshft = INT L × 1 – p and INT is the nearest integer operator. Thus, nblks gives the total

number of overlapped segments (or blocks) that

5
are overlapped by an integer number of points.

The Fourier transform of the jth block from channel m will be

 

ˆ () ( ) – () –

X f L a x j nshft e

L

m

j

m

j if t = + •

=

∑

1 2

0

1

l

l l

l

π ∆

(20)

where the al are time-domain window taper

weights. The window tapers decrease the bias in

the estimate by ensuring that the ends of each

segment of data are gradually tapered to zero.

A variety of window taper functions may be

applied, each having a different impact on the

properties of the resulting spectral estimate. Harris (1978) presents an excellent overview on window tapering strategies for harmonic analysis

problems. Figure 3 graphically illustrates the

overlapped (OBAFFT) segmentation.

2.3.2. Formation of spatial correlation matrix

The cross-power function between all pairs of

sensors in the array is calculated for each block of

data in the OBAFFT method. The final estimate of

the spatial correlation matrix is formed by averaging over the block cross-power matrices.

Using eq 20, elements of the cross-power matrix

for the jth block are

 

R f X f Xf m,n

j

m

j

n

*j ( ) ( )* ( ) = , (21)

where m and n are sensor indices. The final estimate of the spatial correlation matrix is formed by

averaging over the block index,

 

R f nblks m,n m,n Rj

j 1

nblks

( ) =

=

∑ 1

 . (22)

The formation of the cross-power function between two vectors produces a vector that has a

phase angle that is the difference in phase between

the two constituent vectors. For waveforms produced by a stationary source and recorded by channel m, the Fourier transform at any given frequency will be a randomly oriented vector. Thus,

the OBAFFT gives, for each block of data in channel m, a series of vectors with random phase angles.

In an array of sensors recording real data, the quantity to be measured is the phase difference between

channels. The final estimate of the spatial correlation matrix given by eq 22 is therefore an average

of several phase angle difference matrices, which

becomes less variable as the number of blocks

increases. This assumes that the phase difference

between channels is not changing with time.

Taper

Data

X m

nblk 1

nblk 2

nblk 3

nblk 4

nshift

L

0

Time (s)

0.2 0.4 0.6 0.8 1.0

Figure 3. Segmentation and overlapping of a time-domain component of the

signal vector. Each tapered block is transformed into the frequency domain via

an FFT.

6
Figure 4. Spatial correlation matrix structure.

 

ˆ ( )

ˆ ( ) ˆ *

( ) ˆ ( ) ˆ *( ) ˆ ( ) ˆ *( ) ... ˆ ( ) ˆ * ( )

ˆ ( ) ˆ *

( ) ˆ ( ) ˆ *( ) ˆ ( ) ˆ *( ) ... ˆ ( ) ˆ * ( )

ˆ

–

–

R f

X fX f X fX f X fX f X fX f

X fX f X fX f X fX f X fX f

=

00 01 02 0 1

10 11 12 1 1

M

M

X fX f X fX f X fX f X fX f

X fX f X fX f X fX f

20 21 22 2 1

10 11 12

( ) ˆ *( ) ˆ ( ) ˆ *( ) ˆ ( ) ˆ *( ) ... ˆ ( ) ˆ * ( )

...

...

...

ˆ ( ) ˆ *

( ) ˆ ( ) ˆ *( ) ˆ ( ) ˆ *( ) ... ˆ

–

– ––

M

M MM

• •• •

• •• •

• •• •

X fX f M M – – ( ) ˆ * ( ) 1 1













































(23)

2.3.3. Spatial correlation matrix conditioning

The structure of the spatial correlation matrix is

given by eq 23 (Fig. 4). In the programmed implementation, each element of this matrix is normalized according to

Ri j(f) = Xi(f)Xj

*

(f)

Xi(f) X j(f)

 . (24)

If the signal vector is perfectly coherent, then the

diagonal elements of R will contain unity values,

resulting in a maximum beam power response of

1 when applying the BT method. This allows qualitative judgments to be made on the reliability of

the wavenumber estimate, which is particularly

important when dealing with nonideal signal types

or inhomogeneous propagation environments.

When using narrowband BT processing with this

normalization, a peak power response < 0.5 is

often unreliable.

The existence of the R –1 matrix is of critical

importance when applying the ML method (see eq

18). Capon (1969) notes that if the number of blocks

(N) in the estimate is less than the number of array

sensors (M), then R is of order M and rank N and

is thus singular. This is frequently the case when

processing signals of short duration. To guarantee

the nonsingularity of the estimated spatial correlation matrix we add a small amount of incoherent

noise (λ) to the elements of the spatial correlation

matrix. Capon (1969) suggests

R′ = 1 – λ R + λI . (25)

The choice of λ is determined by trial and error (the

results reported here often used values < 10–4). In

practice the accuracy of the beam response is relatively insensitive to λ, but it should be kept as

small as possible since it lowers the signal-to-noise

ratio (SNR) of the beam response.

2.3.4. Variance of the cross-power matrix

The OBAFFT method allows a degree of control

over the bias/variance trade-off in the estimate of

the cross-power matrix. The bias is mitigated by

the choice of window taper and the number of

points in each block, while the variance is a function of the number and degree of overlap between

blocks. The independence of each block is reduced

by increasing the overlap. The characteristics of

the window taper applied to each block also affect

the block‘s independence for a given degree of

overlap. This partial dependence/independence

may be exploited to decrease the variance of a

given spectral estimate.

As an illustration of the variance reduction that

can be obtained using this strategy, a signal (dt =

1/1023) was generated by superimposing two sine

waves with line spectra at 20 and 23 Hz and

unitary peak-to-peak amplitudes. The second sine

wave series is time-shifted (relative to the first

signal) by 0.0455 s (at 20 Hz in the frequency

domain this is an approximately 32° phase shift).

A unique pseudo-random Gaussian ”noise“ series

with zero mean and a maximum peak-to-peak

amplitude of 1 is added to each of the two sine

wave signals. This yields an SNR of 1 for each

spectral line in the series. The OBAFFT spectral

estimation process was applied at 20 Hz for each

sine wave plus noise signal, and the cross-power

operation was carried out. Note that the presence

of the 23-Hz spectral line in the resulting series

allows simulation of phase and amplitude bias

effects. The 20-Hz OBAFFT estimate of the cross

power between the two noisy signals was ob7
Figure 5. Histograms of cross-power amplitude and phase estimates at 20 Hz for a synthetic sine wave series

containing line spectra at 20 and 23 Hz with zero mean random noise. SNR is one. Sample populations for all

histograms are composed of 5000 estimates. All estimates utilized a Blackman window taper to reduce bias effects from

23 Hz spectral line. Cross-power histograms shown in (a) and (b) are based on estimates from blocks of 1024 data

points. Each estimate in the amplitude and phase histograms shown in (c) and (d) are based on six overlapped blocks

of 512 points each. The total series length is 1024 points, and each block overlaps the preceding block by 80%. Note

that the variance in histograms (c) and (d) is significantly smaller than the variance seen in (a) and (b).

tained for 5000 different random noise cases.

Using this sample population of cross-power values, the mean and variance of the amplitude and

phase population can be estimated.

Figure 5 shows cross-power amplitude and

phase estimate histograms using the Blackman

window function (given by W(l) = 0.42 – 0.5 cos

 

( ) /

–

2

1

π L + 0.08 cos 

( ) /

–

4

1

π L with l = 1, ... L). Figure 5a

and 5b are the amplitude and phase histograms for

the case of estimates based on a single block of 1024

points tapered with a Blackman window function.

There are two degrees of freedom in this crosspower estimate. Based on a sample population of

5000 estimates, the amplitude mean is 1.014 with a

variance of 0.0829. The phase mean is 39.36 ° with

a variance of 4.731°2. In Figure 5c and 5d the

amplitude and phase histograms are calculated

but the series is estimated using blocks overlapped

by 80%. Thus, six partially independent blocks are

used in each amplitude and phase estimate. Again

using a population of 5000 estimates, the histogram shown in Figure 5c gives an amplitude mean

8
of 1.013 with a variance of 0.05711. In Figure 5d the

phase estimate population yields a mean of 39.75°

with a variance of 3.249°2. The critical point to note

is the decrease in the variance for the case of the

overlapped six block estimates as compared with

the variance calculated for the single block estimates. Thus, it is statistically advantageous to use

OBAFFT to increase the stability of and decrease

the bias in spectral estimation.

2.4. Verification of signal-processing algorithm

The Fortran code implementing the beamforming procedures given above is discussed in detail

by Moran (1991). The algorithm was extensively

tested during development and continues to be

refined. Two tests were used to validate the procedure. The first test was numerical, using a wide

variety of synthetic signals under various signalto-noise ratios (SNRs) and with a variety of frequency bands and array geometries. The second

and most critical test came from results obtained

from field data collected under defined experimental conditions. These field data tests are based

on the impulsive sources discussed in section 5.2.

3. FIELD DATA ACQUISITION AND

EXPERIMENTAL CONDITIONS

To test the efficacy of beamforming, experimental measurements of various sources were obtained. The environmental characteristics of the

test site and the methods used to make these

recordings are discussed in this section.

3.1. Geological setting

The field tests were performed in February 1988

at Range 37 at the Michigan National Guard’s

Camp Grayling (Fig. 6). The soil in the Grayling

locale consists of glacial outwash sand and gravel

(Farrland and Bell 1982). Logs from four wells in

the vicinity show the thickness of these deposits at

116, 137, 175, and 251 m, with the bedrock below

being Mississippian in age (Lilienthal 1978). A

local driller estimated the sand thickness in the

area to be around 180 m based on his experience.

A borehole 30 m deep was drilled at the site; the

soils logged during the drilling are given in Table

1. The characteristics of eight soil samples that

were collected off the drill auger are listed in Table

2. These observations show that the site was underlain by medium sand, with some thin layers of

fine gravel and clay at depths around 4 and 24 m.

The water table was 5.8 m below the surface.

3.2. Shallow seismic velocity structure

Standard compressional (P) and shear (S) wave

refraction methods were applied at the site to

determine the seismic velocity structure. A linear

array of 24 geophones spaced 1 m apart was arranged perpendicular to the test track. Hammer

blows were then recorded at intervals off each end

of the array. The resulting distance vs. travel time

curves were analyzed using the intercept time

method to determine the seismic velocity as a

function of depth. For the P waves, vertical sledgehammer blows on a metal plate at the ground

surface were used as the source, and Mark Products Model L15-B vertical-component geophones

with resonant frequencies of 4.5 Hz were used to

record the resulting ground motion. For the S

waves, a 0.2 × 0.2 m wooden plank was aligned

perpendicular to the array and clamped to the

ground by the front wheels of a pickup truck (Fig.

7). Horizontal sledgehammer blows on each end

of the plank, which excited SH (shear horizontal)

waves of opposite polarity, were then recorded at

each source location using horizontal-component

geophones.

One difficulty encountered was caused by the

Figure 6. Location map for the seismic experiments,

showing vehicle test track, location of contractor’s sensors, and location of CRREL seismic arrays. Also shown

is the location of the snow pits, met station, frost tubes,

and borehole.

N

Vehicle


Array

Refraction


Array

Recording


Trailer

Met


Station

Vehicle Test Track

9
Table 1. Borehole soil log.*

Depth to Thickness

bottom of layer of layer

Layer (m) (ft) (m) (ft) Soil type

1 3.4 11 3.4 11 Medium sand

2 4.9 16 1.5 5 Fine gravel

3 5.8 19 0.9 3 Medium sand

4 7.6 25 1.8 6 Medium sand and water

5 7.9 26 0.3 1 Clayish fine sand

6 23.2 76 15.2 50 Medium sand and water

7 23.5 77 0.3 1 Blue clay

8 23.8 78 0.3 1 Medium and sand water

9 24.4 80 0.6 2 Blue clay

10 24.7 81 0.3 1 Brown clay

11 27.4 89 2.7 9 Medium sand and water

12 29.3 96 1.8 6 Clean sand and water

13 30.5 100 1.2 4 Medium sand and water with

clay trace in water

*From Hal Carlson, Jim’s Well Drilling, Grayling, Michigan.

presence of the frost layer. For recordings made

close to the hammer source, the first arrival was a

wave traveling directly through the frost layer, at

a speed of about 590 m s–1. Fortunately, this arrival

was highly attenuated and died out within a few

meters. We were able to pick the travel times of the

later-arriving refracted waves after this first arrival and use them in the analysis. The P-wave

velocity of the thin snow layer was measured at

about 100 m s–1 by a short line of geophones placed

above and below the snow cover and spaced 0.25

m apart. The measured first arrival travel times for

the P- and S-wave refraction experiments along

with least-squares line fits are shown in Figures 8

and 9.

Analysis of the travel times indicates a surface

layer 4.8 m thick with a P-wave velocity of 290

m s–1, underlain by a layer with a velocity of 1660

m s–1 that is identified as the water table. The shear

wave data reveal a surface layer with a velocity of

150 m s–1, 12 to 15 m thick, sloping downward at an

angle of about 9 degrees to the north, above a layer

Table 2. Physical properties of borehole soil samples.

Percentage of soil type

Moisture

Sample Depth Density content Gravel Sand Sand Sand

no. (m) (ft) Description (kg m3) (%) fine coarse medium fine Clay

1 1.5 5 Sand 1.38 3.7 7.6 1.8 34.5 63.5 0.3

2 3.0 10 Sand 1.45 3.4 4.6 2.0 33.0 63.9 1.0

3 4.0 13 Sand and gravel 1.56 2.8 18.2 10.2 35.7 35.6 1.1

4 6.1 20 Sand (below 1.45 3.4 0.4 0.3 12.9 86.7 0.4

 water table)

5 13.7 45 Sand 1.45 19.2 0.0 0.2 10.0 88.1 1.7

6 22.9 75 Blue clay — 23.2 0.0 0.1 7.1 53.4 39.4

7 24.4 80 Brown clay — 6.8 0.0 0.1 4.4 35.5 60.0

8 30.5 100 Sand with clay 1.73 19.4 0.0 0.0 15.0 75.7 9.4

The percentage of gravel, sand, and clay in each sample was determined from the grain size measurements. According to the unified

soil classification scheme, soils with grains smaller than 9.5 mm but larger than 2 mm are classified as fine gravel. Sand has grain sizes

less than 2 mm and greater than 0.076 mm and is further divided into coarse, medium, and fine at the 0.59 mm and 0.25 mm sizes.

Grains less than 0.076 mm can be either silt or clay, but are assumed here to be clay based on field observations of the samples.

10
Figure 7. Generating shear waves for the site characterization measurements. The cable leading

from the sledgehammer allows the impact time to be recorded.

100

80

60

40

20

0

100

80

60

40

20

0

0 20 40 60 80 100

Range (m)

Time (ns) Time (ns)

North

South

300

200

100

0

300

200

100

0

0 20 40 60 80

Range (m)

North

South

Time (ns) Time (ns)

Figure 9. Shear (S) wave refraction travel time vs.

distance. The straight lines are least-squares fits to the

travel time data. Shown are arrival times for waves from

hammer sources located north and south of the linear

array of geophones.

Figure 8. Compressional (P) wave refraction travel

time vs. distance. The straight lines are least-squares

fits to the travel time data. Shown are arrival times for

waves from hammer sources located north and south of

the linear array of geophones.

11
Table 3. Analysis of seismic refraction travel time.

Apparent velocity Intercept time Depth Dip

(m s–1) (ms) Velocity (m) angle

Layer South North South North (m s–1) South North (degrees)

Compressional (P) waves

1 276 301 0.0 0.0 288 4.78 4.80

2 1584 1749 32.6 32.8 1662 –0.5

Shear (S) waves

1 167 137 0.0 0.0 152 12.2 15.7

2 417 386 148.3 191.2 401

Figure 10. Compressional (P) wave downhole travel

time vs. distance. Straight-line ray paths were assumed

to calculate the slant distance from the measured sensor

depth and horizontal source offset. These lines yielded

values of 390 m s–1 down to 3 m, underlain by a velocity

of 1370 m s–1.

Slant Range (m)

Time (ms)

0

0

10

20

30

10 20 30

Slant Range (m)

Time (ms)

0

0

40

80

120

10 20 30

Figure 11. Transverse shear (S) wave downhole travel

time analysis. The least-squares lines give 180 m s–1

down to 12 m and 340 m s–1 below that depth.

with a speed of 400 m s–1. Details of the intercept

time analysis are given in Table 3. The seismic

water table depth is about 1 m less than the well log

depth, while the shear wave velocity change does

not correspond to any obvious drill log layer. The

vertical compressional and shear wave velocities

of the shallow soil were also measured by clamping a three-component geophone at various depths

in a cased, 30-m-deep borehole and recording the

travel times of P and S waves from surface hits 3 m

away from the top of the hole to the downhole

geophone. The seismograms from the two

downhole horizontal components of unknown

orientation were rotated to the radial and transverse orientations before the first arrivals were

picked. Figures 10 through 12 show the resulting

travel time curves for the P and S waves, with

least-squares line fits used to determine the velocities. These experiments indicate compressional

wave velocities of 390 m s–1 at the surface, with a

1370-m s–1 layer below, starting at a depth of

around 3 m. The measured shear wave velocity

was 180 m s–1 at the surface and 340 to 380 m s–1 at

a depth of 12 to 16 m. These results agree with the

12
Slant Range (m)

Time (ms)

0

0

40

80

120

10 20 30

Figure 12. Radial shear (S) wave downhole travel time

analysis. A velocity of 180 m s–1 down to 16 m, with 380

m s–1 below that, is indicated by the least-squares line

fit.

refraction results discussed above. The average

velocities over the upper 30 m were 1090 m s–1 for

the P waves and 250 m s–1 for the two S waves.

3.3. Environmental characterization

A number of meteorological parameters were

recorded continuously at the site of the seismic

and acoustic experiments. A Campbell Scientific

Inc. (CSI) Model 21X Micrologger was used to

collect the data (Fig. 13), which were recorded on

a cassette magnetic tape and written on a paper

tape. CSI Model 107 thermistors recorded the temperatures in the ground, in the snow, and at heights

of 0.01, 0.1, 0.3, 1, and 1.4 m above the snow. These

readings have an accuracy of ±0.2°C. A CSI Model

207 sensor provided the temperature and relative

humidity at a height of 2 m. The wind speed and

direction at a height of 3 m was recorded by Met

One Model 014A and 024A sensors with an accuracy of ±2% and ±5%, and a YSI Model 2014 barometer provided the barometric pressure to within 1

mb. Readings from the 12 sensors were taken

every minute, and the average and standard

Figure 13. The instrumentation used to record meteorological parameters. The Campbell micrologger is housed in the box at the center

of the tripod. At the top are the wind speed and direction sensors. The

radiation shield contains a thermistor and the relative humidity

sensor. Additional thermistors and a barometer are mounted on the

wooden post. A lead–acid battery below the tripod provided the

power.

13
Table 4. Meteorological data.

Temperatures (°C)

Pressure RH WS WD

Date Time Gnd Snow 0.01 m 0.1 m 0.3 m 1 m 1.4 m 2 m (mb) (%) (m s–1) (deg)

16 Feb 1630 –3.4 –3.4 –3.9 –4.2 –4.4 –4.2 –3.5 –4.5 963.6 94 3.3 201

16 Feb 1700 –3.3 –3.4 –4.0 –4.2 –4.4 –4.2 –3.5 –4.4 963.8 93 4.0 209

16 Feb 1730 –3.3 –3.4 –5.2 –5.0 –5.2 –5.0 –3.7 –4.9 963.9 92 3.9 207

16 Feb 1800 –3.3 –3.4 –5.4 –5.1 –5.3 –5.2 –4.1 –5.1 963.9 92 3.0 203

16 Feb 1830 –3.3 –3.4 –5.5 –5.1 –5.4 –5.3 –4.4 –5.2 964.0 92 2.7 222

16 Feb 1900 –3.3 –3.4 –5.5 –5.2 –5.5 –5.4 –4.6 –5.3 964.0 93 3.0 195

17 Feb 1100 –2.7 –2.7 –1.8 –1.2 –1.3 –0.8 0.0 –1.5 971.1 80 3.1 267

17 Feb 1130 –2.7 –2.7 –1.9 –1.1 –1.1 –0.8 0.3 –1.4 971.5 80 3.5 311

17 Feb 1200 –2.7 –2.7 –2.2 –1.7 –1.8 –1.5 –0.2 –1.7 971.5 81 5.1 306

17 Feb 1230 –2.6 –2.6 –1.2 –0.7 –0.7 –0.2 0.1 –1.4 971.8 76 5.3 286

17 Feb 1300 –2.6 –2.5 –2.2 –1.5 –1.7 –1.2 0.4 –1.4 971.8 76 4.5 302

17 Feb 1330 –2.5 –2.5 –2.5 –1.6 –1.6 –1.2 0.7 –1.6 971.7 74 6.4 256

17 Feb 1400 –2.5 –2.4 –1.2 –0.3 –0.1 0.3 0.9 –1.2 971.8 74 4.7 272

deviation of the values were computed and recorded every 30 minutes. The meteorological data

recorded on 16 and 17 February 1988 during the

seismic and acoustic array experiments discussed

in this report are tabulated in Table 4. Two frost

tubes were installed at the site in January: one was

located 5 m north of the met station in a topographic high, and the second was located 10 m

south of the met station in a topographic low. The

elevation difference between the two was 0.95 m.

The snow depth at the first frost tube was 0.14 m,

and frost depth exceeded the depth of the tube

(0.775 m) for all of the February readings. On 16

February, the snow depth at the second tube was

0.18 m and the frost depth was 0.60 m. Snow

depths at the site during these experiments ranged

from 0.14 to 0.20 m in most locations. The snow

had been on the ground for some period of time. A

snow pit examined on 18 February revealed the

presence of a thin wind crust at the top underlain

by rounded snow crystals with a layer of depth

hoar at the bottom (Table 5).

3.4. Array description and recording procedures

All the waveform data were digitally recorded

by a Geosource DSS-10A system on 9-track magnetic tape at 1600 bpi in the SEG B format. The

recording system has a dynamic range of 90 dB

and a bandwidth from 3 to 500 Hz at a sampling

rate of 2 kHz. Mark Products L-15B vertical- and

horizontal-component geophones with a resonant

frequency of 4.5 Hz and a sensitivity of 32 V m

s–1 were used to detect the ground motion. Globe

100C low-frequency microphones with a sensitivity of 2 V Pa–1 were also used as sensors. The

geophones and microphones were placed in an 18-

Table 5. Snow characterization, pit 2, at 1020 on 18 February 1988.

Tempera- Crystal

Depth Density ture Hardness size

Layer (mm) (kg m–3) (°C) index (mm) Crystal type

1 0 –3.5 Wind crust (irregular, broken

 crystals), planar dendrites

1 120 –3.5 60 0.5 Rounded, often branched

35 150 –5.0 5 0.5 Rounded, often branched

85 210 –6.5 10 0.5 Rounded, often elongated

2 115 –5.0 0 2.5 Depth, hoar, angular,

 stepped grains

140 160 –4.5 2.5 Depth, hoar, angular,

 stepped grains

165 –4.5 2.5 Depth hoar, angular,

 stepped grains

14
5

4

3

8 7

9

1

2

North

East

0 10 20 m

W

a

v

e

F r o n t s

6

Road

CPA

Figure 14. Grayling sensor positions and road, drawn to

scale. Relative sensor coordinates and associated channel

numbers are given in Table 6. North is referenced to the

magnetic pole. Dots (•) are sensor locations.

element crossed array with one of the legs perpendicular to the test track. Figure 14 shows the sensor

positions, and Table 6 gives the relative element

positions. The array was calibrated by vertical

hammer blows and blank pistol shots from five

locations along the test track: at the closest point of

approach of the vehicle (CPA), and at 50 and 100 m

on either side of the CPA. Shear waves from the

CPA were also recorded, as were pistol shots 50 m

from the array center at the remaining compass

points.

3.5. Characterization of the Grayling array

An important factor in qualitatively assessing

the reliability of an estimated wavenumber spectrum based on the BT or ML methods is the presence of the array factor [F k ] in the estimate. The

array factor is the zero-phase impulse response for

a given array geometry. The array impulse response provides important information on the

direction-dependent resolution capability of the

array as well as an estimate of the beam width for

an ideal plane wave. It is a relatively simple matter

to obtain such a response. The most straightforward

method is given by

F k = 1

M

e ∑ –ik • Rm

m = 0

M – 1

 , (31)

which is simply a Fourier transform of the array

element location vectors Rm. The procedure used

in this analysis is given by Kelley (1967) and differs

from eq 13 only in that it allows for the analysis of

a wideband array response based on a Gaussian

distribution of zero-phase energy.

Table 6. Seismic array element position.

Element E to W (m) N to S (m)

1 0.0 –2.0

2 0.0 0.0

3 30.0 6.0

4 0.0 12.0

5 0.0 16.0

6 1.73 –1.0

7 –5.0 3.0

8 –10.4 6.0

9 –13.8 8.0

15
Figure 15. Grayling array impulse response. Plot ranges over kx and ky

Nyquist limits. Spectral resolution at –3 dB is 0.19 m–1.

The impulse response for the 1989 Grayling

array is shown in Figure 15. The plot ranges to the

Nyquist wavenumber limits along the x and y

axes. The array factor is the x described by the –12-

dB high-amplitude ridge lines radiating from the

center of the plot. The inner contour at –3 dB

defines the spectral resolution for this array, which

is 0.19 m–1. The high-amplitude lobes seen intruding from the periphery of the plot toward the

center are grating lobes, which are a consequence

of the periodic nature of the spectral estimation

procedure.

The Grayling array geometry is not optimal

with respect to multiple source resolution or providing equal resolution capabilities from all bearing directions. Signals arriving from directions in

the vicinity of 0° or 180° and 55° or 235° will be

poorly resolved. Maximum source resolution lies

at directions perpendicular to these values. The

element spacing of the Grayling array has been

designed to allow coverage of a variety of possible

wavelengths ranging from 4 to 30 m.

4. DISCUSSION OF SOURCE

CHARACTERISTICS

To understand the direction estimates of moving seismic–acoustic sources it is instructive to

first consider the spectral properties of simpler

sources. We therefore begin by looking at impulsive acoustic and seismic wave fields generated by

a .45 caliber (blank) pistol shot and a vertically

delivered sledgehammer blow. The pistol shot

was recorded around 1500 on 16 February, and the

hammer blow around 1730 on the same day.

All estimates of the phase velocities discussed

in this section are based on a time-domain moveout method. Each waveform in the time-domain

signal vector is synchronized by specifying the

source direction (θ) and phase velocity (Vo) of the

waveforms. Since the spatial distribution of the

array elements is known, we can calculate the

relative differences in arrival times between all

elements in the array. The time shift may be expressed as

16

3.14

2.51

1.26

1.88

0.63

0.00

–0.63

–1.26

–1.88

–2.51

–3.14

–1.82 –1.45 –1.09 –0.73 –0.36 0.00 0.36 0.73 1.09 1.45 1.82

White Regions = –3 dB

Black Regions = –42 dB

Linear Grey Scale: 

East–West Component of Wavenumber (m ) –1

North–South Component of Wavenumber (m ) –1
Figure 16. Microphone subarray response to a .45 caliber blank pistol shot.

6

1

2

Elm


No.

0

– 0.8 

– 0.6 

Amp


(dB) 0

Time (s)

0.2 0.4 0.6 0.8 1.0

∆tm = xm

2 + ym2 cos αm – θ

Vo

 , (32)

where αm is the angular location of the mth sensor.

Given a digital record of an ideal plane wave

source, the accuracy of eq 32 will be limited by the

time-domain sample interval.

4.1. Impulsive acoustic source

The first field example discussed is the .45 caliber pistol firing blanks toward the array from the

closest point of approach (CPA). The time-domain

microphone and geophone response to a single

pistol shot are shown in Figures 16 and 17. The

microphone waveforms have phase velocities of

approximately 330 m s–1. The seismic waveforms

(Figure 17) have two distinct arrival phases. The

first arrival is a roughly 75-Hz waveform that has

a phase velocity of approximately 330 m s–1. It has

the distinct appearance of the acoustic impulse

observed in the microphone waveforms. Later in

the geophone time series, the waveforms shift to

25 Hz and have a phase velocity of approximately

220 m s–1.

The frequency domain spectra for these signal

vectors are given in Figures 18 and 19. Each spectrum is based on a 2000-point series sampled at

0.0005 s, segmented into 3 blocks of 1024 points,

and overlapped by 53%. A Blackman window

taper was applied to each block.

The primary features to note in the microphone

– 1.1 

0.2

– 0.6 

– 0.9 

– 0.1 

– 0.4 

– 0.6 

– 1.1 

0

Amp


(dB)

5

4

3

2

6

1

7

8

9

Elm


No.

75-Hz Air Wave 25-Hz Surface Wave

0

Time (s)

0.2 0.4 0.6 0.8 1.0

Figure 17. Time-domain vertical geophone subarray response to a .45 caliber

blank pistol shot.

17
0 50 100 150 200 0 50 100 150 200

0 50 100 150 200 0 50 100 150 200

0 50 100 150 200 0 50 100 150 200

0

0.2

0.4

0.6

0.8

1.0

0

0.2

0.4

0.6

0.8

1.0

0

0.2

0.4

0.6

0.8

1.0

– 120

– 100

– 80

– 60

– 40

– 20

0

– 120

– 100

– 80

– 60

– 40

– 20

0

– 120

– 100

– 80

– 60

– 40

– 20

0

Frequency (Hz) Frequency (Hz) Power (dB) Power (dB) Power (dB)

Power (Elm 2) Power (Elm 6) Power (Elm 1)

Figure 18. Microphone spectrum

for .45 caliber pistol shot: (left) Linear scaling and (right) logarithmic

scaling. The large spectral depression at 60 Hz is a result of notch

filtering during the recording process. The 180-Hz peak is probably a

harmonic.

0 50 100 150 200 0 50 100 150 200

0 50 100 150 200 0 50 100 150 200

0 50 100 150 200 0 50 100 150 200

0

0.2

0.4

0.6

0.8

1.0

0

0.2

0.4

0.6

0.8

1.0

0

0.2

0.4

0.6

0.8

1.0

– 120

– 100

– 80

– 60

– 40

– 20

0

– 120

– 100

– 80

– 60

– 40

– 20

0

– 120

– 100

– 80

– 60

– 40

– 20

0

Frequency (Hz) Frequency (Hz) Power (dB) Power (dB) Power (dB)

Power (Elm 2) Power (Elm 6) Power (Elm 1)

Figure 19. Spectra for three center

elements of the geophone subarray

for a .45 caliber pistol shot.

spectra are the broadband nature of the signal and

the presence of two high-amplitude bins at 45 and

70 Hz, separated by a deep spectral low around 60

Hz (produced from a notch filter). The appearance

of the geophone response is considerably more

complex and may be a result of several propagation factors such as modal effects associated with

surface wave propagation and acoustic-to-seismic

coupling. The majority of the energy in verticalgeophone spectra is contained in a band between

18
Amp


(dB)

– 1.1 

0.9

1.6

1.8

2.1

1.5

– 1.2 

0.8

0

Elm


No.

5

4

3

2

6

1

7

8

9

0 0.2 0.4

Time (s)

0.6 0.8 1.0

Figure 20. Vertical geophone subarray response to a stack of eight

hammer blows. The element number and peak-to-peak amplitude are

listed to the left of each trace. The amplitudes are arbitrarily normalized

to the trace from element 9.

10 and 30 Hz, with an additional spectral peak

near 70 Hz.

4.2. Impulsive seismic source

Figure 20 shows a series of vertical-component

geophone waveforms obtained by digitally stacking eight vertically delivered sledgehammer blows.

The source position is located on the road at a point

50 m to the southwest of the CPA at a direction of

–78°. Using eq 32 it is found that the phase velocity

is approximately 220 m s–1. The power spectra for

the three center elements from this array (Fig. 21)

indicate that the signal is broadband with peaks at

20 and 40 Hz.

The pistol and sledgehammer sources indicate

that the dominant seismic energy carrier is a surface wave propagating with a phase velocity of

approximately 220 m s–1 and having a mode frequency of between 20 and 30 Hz. Analysis of the

pistol shot records clearly shows that high-energy

acoustic sources couple strongly with the ground

to produce surface waves. This is consistent with

the previous findings of Albert (1989). In addition

to the acoustically induced seismic surface wave,

the geophone is also excited by the airborne pressure pulse. These propagation modes can interfere

and may be a source of difficulty for beamforming

processors.

4.3. Moving seismic–acoustic sources

Recordings of a tank were also made as a U.S.

Army M60 tank was driven up and down the test

track at various speeds. A malfunction with the

recording equipment prevented continuous signatures from being obtained, so a series of 1-s-long

signature snapshots was recorded. The 10 tank

signatures examined here were recorded between

1256 and 1259 on 17 February.

We now discuss the spectral properties of an

M60 tank moving at approximately 4.5 m s–1 (10

mph). The vehicle’s direction (θ) during this record

interval was approximately 260°. Figures 22 and

23 present the time-domain waveforms observed

by the vertical-geophone and microphone

subarrays, respectively. In the vertical-geophone

waveforms there is a high degree of both temporal

and spatial variation in the signals. The peak-topeak amplitude variation across the array ranges

up to 6.6 dB. Since the spatial variations between

array elements were not consistent from record to

record, it is unlikely that the variations are due to

geophone coupling to the ground or to variation in

the snow cover.

Spectra were estimated using a total signal

length of 2000 points segmented into three blocks

of 1024 points. Each block is overlapped by 53%

and tapered using a Blackman window. Figure 24

19
0 50 100 150 200 0 50 100 150 200

0 50 100 150 200 0 50 100 150 200

0 50 100 150 200 0 50 100 150 200

0

0.2

0.4

0.6

0.8

1.0

0

0.2

0.4

0.6

0.8

1.0

0

0.2

0.4

0.6

0.8

1.0

– 120

– 100

– 80

– 60

– 40

– 20

0

– 120

– 100

– 80

– 60

– 40

– 20

0

– 120

– 100

– 80

– 60

– 40

– 20

0

Frequency (Hz) Frequency (Hz) Power (dB) Power (dB) Power (dB)

Power (Elm 2) Power (Elm 6) Power (Elm 1)

Figure 21. Hammer blow spectra

from center array elements for vertical geophone subarray.

Amp


(dB)

6.1

1.6

6.8

4.7

6.4

3.0

6.6

0

Elm


No.

5

4

3

2

1

7

8

9

0 0.2 0.4

Time (s)

0.6 0.8 1.0

Figure 22. Time-domain vertical

geophone subarray response to

M60 tank moving at 4.5 m s–1.

Amp


(dB)

7.0

1.3

0

Elm


No.

2

1

6

0

Time (s)

0.2 0.4 0.6 0.8 1.0

Figure 23. Time-domain microphone subarray response to M60

tank moving at 4.5 m s–1.

20
0 50 100 150 200 0 50 100 150 200

0 50 100 150 200 0 50 100 150 200

0

0.2

0.4

0.6

0.8

1.0

0

0.2

0.4

0.6

0.8

1.0

– 120

– 100

– 80

– 60

– 40

– 20

0

– 120

– 100

– 80

– 60

– 40

– 20

0

Frequency (Hz) Frequency (Hz) Power (dB) Power (dB)

Power (Elm 2) Power (Elm 1)

Figure 24. Spectra of two

center elements of vertical

geophone subarray for an

M60 tank moving at 4.5 m

s–1. Spectra are strongly

peaked at 11 and 29 Hz.

Figure 25. Spectra of microphone subarray response to

an M60 tank moving at 4.5

m s–1. Spectra are strongly

peaked around 3 Hz and 29

Hz.

0 50 100 150 200 0 50 100 150 200

0 50 100 150 200 0 50 100 150 200

0 50 100 150 200 0 50 100 150 200

0

0.2

0.4

0.6

0.8

1.0

0

0.2

0.4

0.6

0.8

1.0

0

0.2

0.4

0.6

0.8

1.0

– 120

– 100

– 80

– 60

– 40

– 20

0

– 120

– 100

– 80

– 60

– 40

– 20

0

– 120

– 100

– 80

– 60

– 40

– 20

0

Frequency (Hz) Frequency (Hz) Power (dB) Power (dB) Power (dB)

Power (Elm 2) Power (Elm 6) Power (Elm 1)

shows the vertical-geophone power spectra of elements 1 and 2 from Figure 22. The spectra are

strongly peaked at 11 and 29 Hz. Plane-wave

modeling at these two frequencies with a phase

velocity of 220 m s–1 and the Grayling array geometry indicates that signal interference may be the

cause of the amplitude variation across the array.

Interference effects may also be present between

acoustic energy coupled near to the geophone and

surface waves excited near the source. This hypothesis is partially supported by the presence of

a strong spectral peak at 29 Hz in both the microphone and vertical-geophone spectra (Figures 24

and 25).

21
5. PROCESSING RESULTS FOR

FIELD RECORDS

This section presents ML wavenumber estimates of Grayling field data for the .45 caliber

blank pistol, the sledgehammer, and the M60 tank

moving at 4.5 m s–1. All wavenumber estimates

were formed over a regular grid space containing

2601 points. All spectral estimates contain various

degrees of bias and variance. To demonstrate the

degree to which frequency domain bias and variance can affect wavenumber estimates, we show

several results that use different OBAFFT parameters on the same M60 signal vector.

5.1. Criteria used to select processed frequencies

Several criteria were used to select the frequency

bin at which a wavenumber estimate is formed.

The most important of these processing criteria

was the average array coherence spectra, which is

formed by calculating and summing the magnitude-squared coherence spectra (Γave) between all

possible channel combinations and then normalizing by the number of combinations. This operation is given as

(33)

In general, wavenumber spectra which were

processed at frequency bins with coherences

greater than 0.65 and a coherence variance across

immediately adjacent bins of less than 0.025, produced reasonable wavenumber estimates. However, it was found that wavenumber estimates

could still vary by a considerable degree even with

average coherences of 0.9 (usually this was accompanied by high coherence variance across adjacent

bins). Such variation can be explained by recognizing that the coherence spectrum is also susceptible

to bias and variance errors.

The average array coherence spectrum is

strongly affected by the OBAFFT parameters. The

Blackman taper normally produced a smoothly

varying coherence spectrum with peaks ranging

over several frequency bins. In contrast, a boxcar

taper often produced an erratic coherence spectrum. In general, peaks in the coherence spectra

greater than 0.7, which varied smoothly through

several bins, produced acceptable wavenumber

estimates. Wavenumber estimates based on seismic signals from the moving tracked vehicles

tended to be reliable at frequencies between 9 and

15 Hz despite the occurrence of peak spectral

energies in the vicinity of 30 Hz. The 30-Hz bins

rarely produced satisfactory wavenumber estimates.

Several factors were useful in qualitatively determining the accuracy of wavenumber spectra.

The most important was the presence of a strong

array factor. This was particularly relevant when

assessing the quality of wavenumber estimates

based on small numbers of OBAFFT blocks. Another key indicator proved to be the estimate of the

phase velocity. An estimated phase velocity consistent with the surface wave velocity found by the

time-domain move-out analysis was usually a good

indicator of an accurate wavenumber estimate. In

addition, if the spatial correlation is normalized

according to eq 25 then, when using the Bartlett

method, a peak power value of 0.6 or greater often

accompanied an accurate estimate.

5.2. Impulsive sources

The first series of wavenumber estimates is

based on the vertical-geophone subarray response

to the .45 caliber blank fired from the CPA at θ =

–40°. The waveforms of this signal vector are shown

in Figure 17. Figure 26 shows the average array

amplitude and coherence spectra, based on a 2000-

point series with three blocks. Each block contains

1024 points, is overlapped by 53%, and is tapered

using a Blackman window

Recall from section 4 that there are two distinct

arrivals in this record. The wavenumber estimate

shown in Figure 27 is the Capon maximum-likelihood estimate processed at 78 Hz. This frequency

was chosen since it produces a wavelength comparable to the spacing of the center elements of the

geophone array with an acoustic propagation velocity. The time-domain interval considered was

between 0.1 and 0.25 s. This corresponds to a time

window that brackets the initial impulsive arrival

seen in Figure 17. There were 230 sample points in

this interval. The estimate was formed using a

block length of 128 points with each block overlapped 60%, giving a total of three blocks. Boxcar

windowing was applied. The coherence at 78 Hz

was 0.7. The wavenumber spectrum (Fig. 27) produces an estimated source direction of –40° and a

phase velocity of 338 m s–1. This is consistent with

the direction of the pistol shot and the acoustic

propagation velocity. Note that the convention

followed plots the true vector orientation of the

ˆ ( )• ˆ ( )

ˆ ( ) ˆ ( )

*

j i

i j

i j

∑

= +

1

1

XfX f

X fX f

M l l

l l

 

ˆ

– – –

–

Γave i=

( )( ) ∑ ∑ 2

2 1 1

1

nblks M M

nblks M

l

22
1.65

1.32

0.66

0.99

0.33

0.00

–0.33

–0.99

–1.32

–1.65

–1.65 –1.32 –0.99 –0.33 0.00 0.33 0.66 0.99 1.32 1.65

North–South Component of Wavenumber (m ) –1

White Regions = –3 dB

Black Regions = –42 dB

Linear Grey Scale: 

East–West Component of Wavenumber (m ) –1

Figure 27–29

Figure 31

Figure 33

Figure 36–41

Figure 27. ML wavenumber spectrum at 78 Hz for vertical geophone

subarray. The location of the spectral

maxima indicates a source direction

of θ = –40° and a phase velocity of

338 m s–1.

0 40 80 120 160 200

0

– 40

– 80

– 120

1.0

0.4

0.6

0.8

0.2

0

Frequency (Hz)

Power (dB) Coherence

Amplitude Spectrum

Coherence Spectrum

Surface


Wave

Air


Wave

Proc. Band

Figure 26. Average array spectra for vertical

geophone subarray to .45 caliber pistol shot.

23
Figure 28. Estimated ML wavenumber spectrum at 23 Hz for vertical geophone subarray. Location of

spectral peak gives a source direction

of θ = –37° and a phase velocity of

223 m s–1.

wavenumber. Thus, all incoming source wavenumbers are oriented 180° away from the source

direction.

Figure 28 is a 23-Hz ML wavenumber estimate

of the same pistol shot record. The time interval

considered encompasses the later portion of the

wave train between 0.2 and 0.5 s. There were 500

points in this interval. The time-domain move-out

analysis indicated a phase velocity of roughly

220 m s–1 in this interval. A block length of 256

points overlapped by 60% was used, giving four

blocks in the estimated spatial correlation matrix.

A boxcar window taper was applied. The average

array coherence in the 23-Hz bin was 0.8. The

wavenumber spectra in Figure 28 give a source

direction of θ = –37° and a phase velocity of 223 m

s–1. These values agree with the known source

direction and the time-domain estimates of the

propagation velocity.

The last beamformer response for the pistol

shot is shown in Figure 29. This is a broadband ML

estimate processed on an 800-point time interval

between 0.1 and 0.5 s. A block length of 128 points

was chosen with no block tapering. The block

overlap was 40%, giving nine blocks. The frequency band evaluated was from 15 to 78 Hz.

Selection of these parameters allowed the spectral

estimate of each block to retain a high degree of

independence with emphasis on the resolution of

specific frequency bins. In Figure 29 there are three

spectral maxima that have energy levels greater

than –3 dB. All three trend along a high-amplitude

“ridge,” indicating a source direction of roughly

–40°. The wavenumber peak with the smallest

magnitude yields a phase velocity of 223 m s–1 if

normalized by a frequency of 23 Hz. The

wavenumber peak with the intermediate magnitude yields a phase velocity of 342 m s–1 when

normalized by 78 Hz.

The above series of wavenumber spectra and

the time-domain plots clearly show the existence

of an acoustically induced seismic surface wave.

In considering the ML wavenumber spectra it

should be remembered that the wavefield was

generated by a nearly ideal stationary point source.

The OBAFFT parameters used in the ML

wavenumber spectra were chosen through trial

and error, and in many circumstances reasonable

parameters did not yield intelligible beam responses. The above set of wavenumber estimates

(Fig. 27 to 29) serve as a testament to the impact of

the frequency domain estimation methodology on

24

1.65

1.32

0.66

0.99

0.33

0.00

–0.33

–0.99

–1.32

–1.65

–1.65 –1.32 –0.99 –0.33 0.00 0.33 0.66 0.99 1.32 1.65

North–South Component of Wavenumber (m ) –1

White Regions = –3 dB

Black Regions = –42 dB

Linear Grey Scale: 

East–West Component of Wavenumber (m ) –1

Figure 27–29

Figure 31

Figure 33

Figure 36–41
Figure 29. Broadband ML wavenumber spectrum for vertical geophone subarray. The processed frequency band is 15 to 78 Hz. Highamplitude linear ridge along a

wavenumber space direction of 140°

indicates a source position of –40°.

The locations of spectral maxima

indicate waves propagating at 223

and 342 m s–1.

the resulting wavenumber spectra.

The next ML wavenumber spectra we consider

are for the vertical-geophone signal vector produced by a vertically delivered sledgehammer

blow (see Fig. 20). The source location was on the

road 50 m south of the CPA at θ = –78°. Figure 30

displays the average array power and coherence

spectra. The OBAFFT parameters used for the

power, coherence, and wavenumber spectra were

a signal length of 2000 points segmented into three

blocks of 1024 points, each overlapped by 53% and

tapered with a Blackman window. The 95% variance bounds were 

Rˆ • (0.46, 3.96) with 5.51 d.f. The

25-Hz ML wavenumber spectra based on this

record are shown in Figure 32. The average array

coherence at 25 Hz is 0.64. The location of the

spectral maximum in the wavenumber estimate

gives θ = –74° with a phase velocity of 232 m s–1.

5.3. Moving sources

The first source we consider (Figure 22) is an

M60 tank travelling along the road from southwest to northeast at 4.5 m s–1. This was recorded

when the vehicle was 300 m to the southwest of the

CPA. The OBAFFT parameters used to produce

the wavenumber, power, and coherence spectra

Figure 30. Average array spectra for vertical geophone

subarray to stack of eight hammer blows.

0 40 80 120 160 200

0

– 40

– 80

– 120

1.0

0

0.2

0.4

0.6

0.8

Frequency (Hz)

Power (dB) Coherence

Amplitude Spectrum

Coherence Spectrum

– 160

25

1.65

1.32

0.66

0.99

0.33

0.00

–0.33

–0.99

–1.32

–1.65

–1.65 –1.32 –0.99 –0.33 0.00 0.33 0.66 0.99 1.32 1.65

North–South Component of Wavenumber (m ) –1

White Regions = –3 dB

Black Regions = –42 dB

Linear Grey Scale: 

East–West Component of Wavenumber (m ) –1

Figure 27–29

Figure 31

Figure 33

Figure 36–41
Figure 31. Estimated ML wavenumber spectrum at 25 Hz for vertical geophone subarray for hammer blows. Peak amplitude gives a

source direction of θ = –74° and a

phase velocity of 232 m s–1.

Figure 32. Average array spectra for vertical geophone

subarray to M60 tank moving at 4.5 m s–1.

0 40 80 120 160 200

0

– 40

– 80

– 120

0

0.2

0.4

0.6

0.8

1.0

Frequency (Hz)

Power (dB) Coherence

Amplitude Spectrum

Coherence Spectrum

– 160

are a series length of 2000 points segmented into

three blocks of 1024 points, each overlapped by

53%, and tapered with a Blackman window function. The 95% confidence bounds are 

Rˆ • (0.46, 3.96)

with 5.51 d.f. Figure 32 shows the average array

power and coherence spectra. There is a peak in

the power spectrum at 11 Hz that has a coherence

of 0.94. The 11-Hz ML wavenumber spectrum

(Fig. 33) estimates the source direction as θ = 258°

with a phase velocity of 219 m s–1.

The result shown in Figure 33 is typical of ML

wavenumber estimates obtained from the M60 in

the 9- to 15-Hz band using the vertical-geophone

subarray. To demonstrate this, consider Figure 34,

which is a composite of direction and velocity

magnitudes estimated from ML wavenumber spectra for a consecutive series of M60 vehicle records.

The location of each direction position is sequentially labeled 1 to 16. All records in this composite

are from the same M60 drive-by series, with the

tank moving from southwest to northeast along

the road at a constant speed of 4.5 m s–1. The most

distant records in this series were at positions 1

and 16, which correspond to vehicle ranges of 500

and 300 m, respectively. At position 11 the tank is

50 m from the array’s coordinate origin.

26

1.65

1.32

0.66

0.99

0.33

0.00

–0.33

–0.99

–1.32

–1.65

–1.65 –1.32 –0.99 –0.33 0.00 0.33 0.66 0.99 1.32 1.65

North–South Component of Wavenumber (m ) –1

White Regions = –3 dB

Black Regions = –42 dB

Linear Grey Scale: 

East–West Component of Wavenumber (m ) –1

Figure 27–29

Figure 31

Figure 33

Figure 36–41
Figure 33. Example of an ML

wavenumber estimate from a vertical geophone array at 11 Hz for an

M60 tank moving at 4.5 m s–1.

Coherence was 0.94, the spectral

maximum yields a velocity of 219 m

s–1, and the source direction = 258°.

– 250 – 200 – 150 – 100 – 50 0 50 100 150 200 250

East-West Velocity (m s– 1)

North-South Phase Velocity (m s– 1)

250

200

150

100

50

0

50

100

150

200

250

225 m s– 1

North

East

CPA

16

14

13

11

1

2

3 to 8 9

10

15

12

Figure 34. Composite plot of estimated source direction and velocity

magnitudes for a series of M60 tank records observed by the vertical

geophone subarray.

27

1.65

1.32

0.66

0.99

0.33

0.00

–0.33

–0.99

–1.32

–1.65

–1.65 –1.32 –0.99 –0.33 0.00 0.33 0.66 0.99 1.32 1.65

North–South Component of Wavenumber (m ) –1

White Regions = –3 dB

Black Regions = –42 dB

Linear Grey Scale: 

East–West Component of Wavenumber (m ) –1

Figure 27–29

Figure 31

Figure 33

Figure 36–41
The composite uses a range of temporal frequencies from 9.77 to 29 Hz with the majority of

the records using the 11-Hz bin. All the estimates

used 2000-point signal lengths, segmented into

three blocks of 1024 points, each overlapping the

previous by 53%. A Blackman taper was used for

10 of the wavenumber estimates, with the remaining five estimates using a boxcar window taper.

Using the boxcar window, the 95% confidence

bounds on the crosspower matrix were

 

Rˆ • (0.42, 5.85) with 3.9 d.f. The decision to use

either a boxcar or Blackman taper was based on

the coherence spectra and the consistency of the

velocity estimates of the resulting wavenumber

spectra.

5.4. Affect of spectral estimation parameters on

beam response

In the preceding sections we demonstrated that

ML wavenumber spectra from small seismic arrays can produce accurate direction estimates with

careful choice of transformation parameters. Another result is that small differences in the frequency-domain processing parameters often have

large effects on the wavenumber estimate. Such

effects are termed a wavenumber bias. In this

section several examples of the wavenumber bias

phenomenon are demonstrated for the M60 tank.

When using reasonable OBAFFT parameters,

the wavenumber bias effect was negligible for

signal vectors observed from sledgehammer blows

or blank pistol shots. Thus, the wavenumber bias

problem is most likely associated with the characteristics of the M60. Given the difficulty in producing accurate wavenumber estimates from the vertical-geophone array in the 29-Hz peak power bin

of the M60 source, and the coincident occurrence

of large energy levels at 29 Hz in the microphone

and geophone power spectra (Fig. 25 and 26), it is

tentatively proposed that the bias effect may be

aggravated by interference between surface waves

and acoustic energy coupling near the geophone

and source movement.

5.5. Spatial frequency bias from moving sources

To show how changes in processing parameters can affect wavenumber estimation for the

moving M60, we consider a series of wavenumber

estimates using the same data set but applying

slightly different OBAFFT parameters. In each

case it will be shown how the estimation parameters lead to biased wavenumber estimates or to a

reduction in the SNR as measured by the background energy level of the wavenumber spectrum.

Consider the signal vector shown in Figure 35

for an M60 tank moving along the road at 4.5

Amp


(dB)

2.6

3.6

0.1

2.3

0.2

–

–

–

–

–

1.4

1.6

0

Elm


No.

5

4

3

2

1

7

8

9

0

Time (s)

0.2 0.4 0.6 0.8 1.0

Figure 35. Time-domain response of vertical geophone subarray to M60 tank

moving at 4.5 m s–1 along the road approximately 100 m southeast of the CPA.

These waveforms correspond to position 8 in the composite plot shown in

Figure 35.

28
m s–1 from southwest to northeast. An 11-Hz ML

wavenumber spectrum is shown in Figure 36. It

was formed using a series length of 2000 points

segmented into 20 blocks; each block had 1024

points, was tapered by a Blackman window, and

overlapped the preceding block by 95%. The 95%

confidence bounds on the variance of the

crosspower matrix was 

Rˆ • (0.38, 9.39) with 2.83

d.f. The average array coherence at 11 Hz was 0.92.

The variance estimate leads one to expect that no

significant improvement in the wavenumber spectra will result from such highly overlapped blocks.

This is not the case, as will be shown later.

The wavenumber spectrum in Figure 36 is notable in several respects. It has an extremely sharp

spectral maximum with background energy levels

as low as –45 dB. The highest side lobe in wavenumber space is 30 dB below the maximum. The

location of the spectral maximum gives a source

direction of 258° and a propagation velocity of 219

m s–1.

The ML wavenumber estimate shown in Figure

38 uses the same signal vector as that of Figure 36.

In this spectrum the OBAFFT parameters are identical except the blocks are boxcar windowed. The

variance estimate is 

Rˆ • (0.42, 5.83) with 3.89 d.f. at

the 95% confidence level. In this extremely overlapped region the variance bounds are slightly

better than those found for the Blackman window

function used in the previous example. Figure 37

is a notable wavenumber estimate in that it has an

exceptionally peaked spectral maximum, which

indicates a source direction of 258° and a phase

velocity of 219 m s–1. There are several subtle

differences between Figures 37 and 38. The first is

the occurrence of a relatively pronounced spectral

side lobe in Figure 37 that is 18 dB smaller than the

spectral maximum. In addition, the background

energy level in Figure 37 is approximately –33 dB

while that of Figure 36 was approximately –45 dB.

The next two wavenumber estimates investigate the impact of reducing the number of blocks

in the OBAFFT. In Figures 38 and 39 we process a

signal length of 2000 points, blocked into segments of 1024 points, with each overlapped by

53%, giving a total of three blocks. In the case of

Figure 38 we apply a Blackman window taper to

each block, whereas the wave number estimate

shown in Figure 39 used a boxcar window function.

In Figure 38 the spectral maximum is sharply

peaked, indicating a source direction of 258° and a

Figure 36. ML wavenumber estimate at 11 Hz for waveforms shown

in Figure 35: Blackman taper, 20

blocks (at 95% overlap), segment

length = 1024; produces a source

direction of 258° with a velocity of

219 m s–1.

29

1.65

1.32

0.66

0.99

0.33

0.00

–0.33

–0.99

–1.32

–1.65

–1.65 –1.32 –0.99 –0.33 0.00 0.33 0.66 0.99 1.32 1.65

North–South Component of Wavenumber (m ) –1

White Regions = –3 dB

Black Regions = –42 dB

Linear Grey Scale: 

East–West Component of Wavenumber (m ) –1

Figure 27–29

Figure 31

Figure 33

Figure 36–41
Figure 37. ML wavenumber estimate at 11 Hz for waveforms shown

in Figure 35: boxcar taper, 20 blocks

(at 95% overlap), segment length =

1024; produces a source direction of

258° with a velocity of 219 m s–1.

phase velocity of 219 m s–1. The highest secondary

side lobe is at –21 dB with an overall background

energy level of –36 dB. In contrast, the wavenumber

estimate shown in Figure 39 indicates a source

direction of 255° and a phase velocity of 271 m s–1.

The 3° direction difference and 52 m s–1 velocity

difference between Figures 38 and 39 is a mild

example of the wavenumber bias effect. In addition, applying the boxcar window effectively reduces the SNR of the beamformer in the sense that

the background energy level of Figure 38 increases

by 6 dB in Figure 39.

Comparing Figures 36 and 37 with 38 and 39 we

note that there seems to be significant benefit in

using a larger number of highly overlapped blocks.

Each estimate of the phase difference between two

channels for a given block of data is calculated

from a comparison of two vectors. Each vector has

a magnitude and a phase, both of which vary with

time, and each contains a small amount of random

noise. The amplitude at a given frequency will

vary slowly across blocks for a given channel,

indicating highly correlated blocks. The phase

angle variation across blocks changes (rotational

orientation) rapidly compared with the amplitude

variation, indicating a lower degree of correlation

between blocks. Thus, if one is interested in estimating the phase angle difference (i.e., calculating a spatial correlation matrix) between channels in an array and the signal duration is

fixed, it is advantageous to use a large number of

highly overlapped blocks (compare Fig. 36 and

Fig. 38).

Continuing with the same data, we contrast the

previous four satisfactory wavenumber spectra

with another pair of 11-Hz ML wavenumber estimates formed using a single block of 512 points.

This is theoretically more than twice the signal

duration needed to adequately observe an 11-Hz

waveform propagating at 220 m s–1.

The 11-Hz ML wavenumber spectrum in Figure 40 applied a Blackman window taper to the

block of data, and the 11-Hz ML wavenumber

spectrum in Figure 41 used a boxcar window

function. Figure 40 has an acceptably peaked spectral maximum indicating a source direction of 258°

and a phase velocity of 219 m s–1. However, this

estimate has a considerably higher background

energy level (–18 dB) compared with the previous

examples. The highest spectral side lobe has an

energy level of –12 dB. In a rather dramatic contrast, we note that the ML wavenumber spectrum

30

1.65

1.32

0.66

0.99

0.33

0.00

–0.33

–0.99

–1.32

–1.65

–1.65 –1.32 –0.99 –0.33 0.00 0.33 0.66 0.99 1.32 1.65

North–South Component of Wavenumber (m ) –1

White Regions = –3 dB

Black Regions = –42 dB

Linear Grey Scale: 

East–West Component of Wavenumber (m ) –1

Figure 27–29

Figure 31

Figure 33

Figure 36–41
Figure 38. ML wavenumber estimate at 11 Hz for waveforms shown

in Figure 35: Blackman taper, three

blocks (at 53% overlap), segment

length = 1024; produces a source

direction of 258° with a velocity of

219 m s–1.

Figure 39. ML wavenumber estimate at 11 Hz for waveforms shown

in Figure 35: boxcar taper, three

blocks (at 53% overlap), segment

length = 1024; produces a source

direction of 255° with a velocity of

271 m s–1.

31

1.65

1.32

0.66

0.99

0.33

0.00

–0.33

–0.99

–1.32

–1.65

–1.65 –1.32 –0.99 –0.33 0.00 0.33 0.66 0.99 1.32 1.65

North–South Component of Wavenumber (m ) –1

White Regions = –3 dB

Black Regions = –42 dB

Linear Grey Scale: 

East–West Component of Wavenumber (m ) –1

Figure 27–29

Figure 31

Figure 33

Figure 36–41

1.65

1.32

0.66

0.99

0.33

0.00

–0.33

–0.99

–1.32

–1.65

–1.65 –1.32 –0.99 –0.33 0.00 0.33 0.66 0.99 1.32 1.65

North–South Component of Wavenumber (m ) –1

White Regions = –3 dB

Black Regions = –42 dB

Linear Grey Scale: 

East–West Component of Wavenumber (m ) –1

Figure 27–29

Figure 31

Figure 33

Figure 36–41
32

Figure 41. ML wavenumber estimate at 11 Hz for waveforms

shown in Figure 35: boxcar taper,

one block, segment length = 512.

Beam is unfocused.

Figure 40. ML wavenumber estimate at 11 Hz for waveforms

shown in Figure 35: Blackman

taper, one block, segment length =

512. Spectra give a source direction of 25° and velocity of 219

m s–1.

1.65

1.32

0.66

0.99

0.33

0.00

–0.33

–0.99

–1.32

–1.65

–1.65 –1.32 –0.99 –0.33 0.00 0.33 0.66 0.99 1.32 1.65

North–South Component of Wavenumber (m ) –1

White Regions = –3 dB

Black Regions = –42 dB

Linear Grey Scale: 

East–West Component of Wavenumber (m ) –1

Figure 27–29

Figure 31

Figure 33

Figure 36–41

1.65

1.32

0.66

0.99

0.33

0.00

–0.33

–0.99

–1.32

–1.65

–1.65 –1.32 –0.99 –0.33 0.00 0.33 0.66 0.99 1.32 1.65

North–South Component of Wavenumber (m ) –1

White Regions = –3 dB

Black Regions = –42 dB

Linear Grey Scale: 

East–West Component of Wavenumber (m ) –1

Figure 27–29

Figure 31

Figure 33

Figure 36–41
origin. The vertical-geophone responses to the

M60 tank were dominated by the presence of

Rayleigh surface wave energy with spectral peaks

at 11 and 30 Hz and a propagation velocity of

220 m s–1. Most wavenumber estimates used the

highly coherent secondary maxima in the vicinity

of the 11-Hz frequency bins. It was difficult to

obtain reliable wavenumber estimates in the vicinity of the 30-Hz spectral maxima. It is hypothesized that the low coherence in the vicinity of 30

Hz is due to interference at the sensor between

acoustic and seismic waves.

Maximum-likelihood wavenumber estimates

of .45 caliber blank pistol shots exhibited two

distinct propagation modes. The first wave phase

arrives with a propagation velocity of 338 m s–1

and is interpreted as an acoustic excitation of the

ground in the vicinity of the geophone. The second

wave phase arrives with a propagation velocity of

approximately 220 m s–1 and is interpreted as a

Rayleigh surface wave that is excited by acousticto-seismic coupling near the source. Wavenumber

estimates of vertically delivered sledgehammer

blows also exhibited surface wave propagation

modes that had the same phase velocity.

In general, direction and velocity estimates obtained from wavenumber spectra showed good

agreement with known source positions and velocity estimates produced by time-domain moveout analyses. However, it was observed that

changes in the parameters used to estimate the

spatial correlation matrix could produce wavenumber bias effects. These spatial frequency bias

effects were strongly dependent on the type of

window taper applied to the time series data and

the number of blocks used in the estimate of the

spatial correlation matrix. It was observed that

low-resolution, low-bias window types produced

the highest resolution wavenumber spectra with

the least wavenumber bias and highest signal-tonoise ratios. Utilizing large numbers of blocks in

the estimate of spatial correlation matrix mitigated the impact of window taper choice. In some

cases improvements were found in the wave-number spectra even when block overlaps were as high

as 95%.

LITERATURE CITED

Albert, D.G. (1989) Preliminary analysis of acoustic-to-seismic coupling experiments at Grayling,

Michigan. Paper presented at the 17th Meeting

of NATORSG.11, 23–27 May 1988, Issy-lesin Figure 41 is essentially unfocused for the boxcar

taper.

Figure 41 is an extreme example of beamformer

failure; however, it was not uncommon to encounter similar problems in less coherent frequency

bins with more reasonable OBAFFT parameters.

The failure of the beamformer to focus is most

likely due to the highly biased nature of the frequency-domain phase angle caused by energy

leakage from neighboring frequency bins. Such

bias effects are particularly problematic with the

boxcar window. The above examples show the

extent to which the introduction of a high degree

of bias and variance in the frequency domain

estimate impacts the wave number estimate.

Considering the wavenumber results in this

section and those given in other sections, we make

several generalizations concerning the impact of

the OBAFFT parameters on ML beamformer. In

virtually all the examples discussed, the application of a strongly tapered, low-resolution, lowbias, time-domain window function (such as the

Blackman taper) produced a higher resolution

wavenumber estimate with lower background

energy levels when compared with wavenumber

estimates using high-resolution, high-bias, timedomain tapers (such as the boxcar function).

Specifying large numbers of blocks in the

OBAFFT estimation of the spatial correlation matrix mitigated the impact of the high-bias window

functions. Improvements in beamformer performance were observed even when blocks were

overlapped by as much as 95%. The unexpected

benefit of using highly overlapped blocks in the

estimate of the spatial correlation matrix may be

explained by noting that the phase angle for a

given frequency bin in each block varies rapidly

between blocks and therefore has a much smaller

degree of correlation between blocks as compared

with power spectra estimates. The variance analysis was helpful in partially explaining the

beamformer’s dependence on the OBAFFT parameters. A more complete explanation will have

to consider the phase angle bias effects as well.

6. CONCLUSIONS

Capon’s minimum variance (maximum likelihood) beamformer produces reliable direction estimates for a U.S. Army M60 tank moving at 4.5

m s–1 using an array of vertical-component geophone signals. Accurate direction estimates were

obtained at ranges from 50 to 500 m from the array

33
Moulineaux, France. USA CRREL Internal Report

1027.

Capon, J. (1969) High-resolution frequency–

wavenumber spectrum analysis. Proceedings of the

IEEE, 57: 1408–1418.

Capon, J., R.J. Greenfield and R.J. Kolker (1967)

Multidimensional maximum-likelihood processing of a large aperture seismic array. Proceedings of

the IEEE, 55(2): 192–211.

Carter, G.C. (1977) Bias in bearing estimation resulting from source motion. Journal of the Acoustical Society of America, 62(6): 1447–1448.

Cress, D.H. (1976) Seismic methods of locating

military ground targets. U.S. Army Waterways

Experiment Station, Vicksburg, Mississippi, Paper M-76-13.

De Graaf, S.R. and D.H. Johnson (1985) Capability of array processing algorithms to estimate

source bearings. IEEE Transactions, Acoustics,

Speech, and Signal Processing, ASSP-33(6): 1368–

1379.

Farrland, W.R. and D.L. Bell (1982) Quaternary

Geology of Southern Michigan, 1:500,000.

Gragg, R.F. (1990) Correlations among the modal

components of the far field of a moving source:

Low-frequency Doppler effects in shallow water.

Journal of the Acoustical Society of America, 88(5):

2326–2337.

Hahn, W.R. (1975) Optimum signal processing

for passive sonar range and bearing estimation.

Journal of the Acoustical Society of America, 58: 201–

207.

Harris, F.J. (1978) On the use of windows for

harmonic analysis with the discrete Fourier transform. Proceedings of the IEEE, January, 66: 51–83.

Johnson, D.H. (1982) The application of spectral

estimation methods to bearing estimation problems. Proceedings of the IEEE, September, 70(9):

1018–1028.

Kelley, E.J. (1967) Response of seismic arrays to

wideband signals. Massachusetts Institute of Technology, Lincoln Laboratory, Technical Note 1967-

30, June.

Knapp, C.H. and G.C. Carter (1977) Estimation of

time delay in the presence of source or receiver

motion. Journal of the Acoustical Society of America,

61(6): 1545–1549.

Lacoss, R.T., R.K. Brienzo and T.E. Jones (1991)

Characteristics of acoustic and seismic signals for

vehicle detection and direction finding. Paper presented at the 24th Meeting of the NATO AC/243,

Panel 3, Research Study Group 11, Los Angeles,

California.

Lilienthal, R.T. (1978) Stratigraphic cross-sections

of the Michigan Basin (Geological Survey Division, Lansing, Michigan). Report of Investigation

19.

Moran, M.L. (1991) Users manual for ESTK2D.FOR

and ESTK1D.FOR wavenumber estimation routines. USA Cold Regions Research and Engineering Laboratory, Special Report 91-15.

Nawab, H. S., F.U. Dowla and R.T. Lacoss (1985)

Direction determination of wideband signals. IEEE

Transactions, Acoustics, Speech, and Signal Processing, ASSP-33(4): 114–122.

Priestly, M.B. (1981) Spectral Analysis and Time

Series. New York: Academic Press.

Welch, P.D. (1967) The use of fast Fourier transform for the estimation of power spectra; a method

based on time averaging over short, modified

periodograms. IEEE Transactions, Audio and Electroacoustics, AU-15: 70–73.
1. AGENCY USE ONLY (Leave blank) 2. REPORT DATE 3. REPORT TYPE AND DATES COVERED

4. TITLE AND SUBTITLE 5. FUNDING NUMBERS

6. AUTHORS

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) 8. PERFORMING ORGANIZATION

 REPORT NUMBER

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES) 10. SPONSORING/MONITORING

 AGENCY REPORT NUMBER

11. SUPPLEMENTARY NOTES

12a. DISTRIBUTION/AVAILABILITY STATEMENT 12b. DISTRIBUTION CODE

13. ABSTRACT (Maximum 200 words)

14. SUBJECT TERMS 15. NUMBER OF PAGES

16. PRICE CODE

17. SECURITY CLASSIFICATION 18. SECURITY CLASSIFICATION 19. SECURITY CLASSIFICATION 20. LIMITATION OF ABSTRACT

 OF REPORT OF THIS PAGE OF ABSTRACT

NSN 7540-01-280-5500 Standard Form 298 (Rev. 2-89)

Prescribed by ANSI Std. Z39-18

298-102

REPORT DOCUMENTATION PAGE Form Approved

OMB No. 0704-0188

Public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, gathering and

maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this collection of information,

including suggestion for reducing this burden, to Washington Headquarters Services, Directorate for Information Operations and Reports, 1215 Jefferson Davis Highway, Suite 1204, Arlington,

VA 22202-4302, and to the Office of Management and Budget, Paperwork Reduction Project (0704-0188), Washington, DC 20503.

June 1996

Source Location and Tracking Capability of a Small Seismic Array

Mark L. Moran and Donald G. Albert

U.S. Army Cold Regions Research and Engineering Laboratory

72 Lyme Road

Hanover, N.H. 03755-1290

44

UNCLASSIFIED UNCLASSIFIED UNCLASSIFIED UL

CRREL Report 96-8

DA: 4A762784AT42

TA: SS

WU: E02

Approved for public release; distribution is unlimited.

Available from NTIS, Springfield, Virginia 22161.

Recordings of seismic wavefields from various sources were obtained using a small array of vertical-component

geophones under winter conditions at Grayling, Michigan. These data were processed using a frequencywavenumber domain Capon minimum variance beamformer to estimate the bearing angle and propagation

velocity of the waves emitted from the source. The cross power matrix was adaptively estimated using a tapered

block-averaging procedure. The wave sources were sledgehammer blows on the ground surface, .45 caliber blank

pistol shots, and an M60 tank moving at 4.5 m s–1 along a road near the array. Reliable wavenumber spectra were

obtained for all sources. Processing results for the hammer blows show that the dominant seismic arrival is a

Rayleigh wave traveling at roughly 220 m s–1. For the pistol shots, two arrivals corresponding to the airwave (338

m s–1) and the air-coupled Rayleigh waves (220 m s–1) were observed. The results for these sources were relatively

insensitive to the processing parameters used. For the moving vehicle, the dominant signals observed were

Rayleigh waves (220 m s–1). Accurate locations were obtained for this moving source, although the processing

parameters had to be carefully selected, and the choice of frequency parameters affected the accuracy of the

wavenumber results. Maximizing the number of degrees of freedom and the coherence of the frequency estimates

and minimizing the variation of the coherence across adjacent frequency bins provided the most consistently

reliable strategy for obtaining accurate wavenumber estimates for the moving vehicle. The sensitivity of the

wavenumber estimates to the frequency processing parameters seems to be related to the bias in the phase spectra

of the signals and will potentially occur in any bearing estimation method that uses temporal frequncy phase

spectra.

Beam forming Seismic array processing

Capon method Target location and tracking

Frequency-wavenumber estimation Unattended ground sensors

Maximum likelihood method Wavenumber bias effects