Acoustic Clutter Removal

With the widespread availability of commercial and recreational underwater vehicles, maintaining underwater security in ship channels and harbors is more challenging than ever before. Acoustic clutter from moored and moving vessels, bottom objects, and docks make threat detection technically difficult and costly. The objective of this article is to demonstrate how basis vector properties of maximum length sequences may be used in a bistatic geometry to coherently subtract clutter from the received signal, thereby transforming a difficult signal-in-clutter detection problem into a more tractable signal-in-noise detection problem. Signal processing steps are described in Part 1. Steps that are critical, or have not been well documented previously, are in italics. In Part 2, the performance of a notional traffic monitoring system with underwater coverage is simulated for a segment of the Hong Kong Ship Channel. Although the focus of this work is on underwater applications, the algorithm may apply in other situations where small targets must be detected in the presence of much larger objects. Active sonar, medical imaging, and ground penetrating radar are among the potential applications.


I. INTRODUCTION
M AXIMUM length sequences, or m-sequences, have been used to probe the dynamics of oceans since Steinberg and Birdsall's 1964 pioneering experiment in the Florida Straits [1]. In the following 50 years, oceanographic studies using m-sequences and near-stationary acoustic sources and receivers were conducted at longer range, eventually leading to the Heard Island experiment, tomography, ATOC, and others [2]- [9]. More recent experiments were conducted in geometries and at frequencies related to coastal and harbor defense [10]- [15].
Collectively, past research has shown that, when motion from the source and receiver is small relative to an acoustic wavelength, refracted and bottom reflected acoustic paths are remarkably stable, subject only to changes in oceanographic conditions. Energy interacting with the ocean surface, on the other hand, may be affected by waves and surface currents and exhibit significant time and frequency spreading.
The objective of this article is to demonstrate how the basis vector properties of m-sequences allow coherent subtraction of clutter from the received signal, thereby transforming a difficult signal-in-clutter detection problem into a more tractable signal-in-noise detection problem. The result is a relatively new approach to harbor and coastal defense using fixed multistatic sources and receivers with low-power m-sequence transmissions.
Part 1 of this article is organized according to the signal processing steps leading to clutter removal and detection: fundamental equations, signal generation, demodulation, correlation, Doppler processing, and clutter removal. Signal processing steps that are critical, or have been poorly documented previously, are in italics. The theoretical development from Part 1 is applied to a notional harbor defense system for the Hong Kong Ship Channel in Part 2.

A. Fundamental Equations
Maximum length sequences, m-sequences, are strings of +1s and 0s derived from coefficients of prime polynomials. The sequence length, L, is related to the degree of the polynomial, L = 2 N − 1, where N is the degree of the polynomial. The theory has been well documented by Golomb [16] and will not be repeated here. However, if 0s are replaced with -1s, and the m-sequence is multiplied by arctan(sqrt(L)) [17], then the following very important relationship is obtained for circular correlation: where n = mod(n, L), " * " indicates conjugation, and "·" indicates multiplication.
From (1), m-sequences may be viewed as orthogonal basis vectors that span a sample space of dimension L. Within this article, the term "m-sequences" will refer to modified sequences with the properties of (1).
The circular correlation, shown in the following equation, is the projection of the input sample vector s onto the m-sequence m(n): If a segment of the correlation Γ(p) is identified as clutter, the clutter may be removed by setting the nonclutter portion of the correlation to zero, transforming the modified correlation Γ (p) back to the time domain, correcting for Doppler frequency shift, correcting time compression/expansion and, finally, subtracting the clutter segment from the original signal vector. The division to obtain the signal vector with clutter removed is facilitated by the property in (1) If the clutter in (3) is from a moving vessel, it may be identified visually or by radar. If the clutter is from bottom objects, distant shores, or moored vessels, it can be identified by its consistency in range, bearing, and Doppler. An example is provided in Part 2.
It should be noted that subtraction of the time-domain clutter signal removes all the energy associated with clutter, including the cross-correlation energy in the ambiguity function floor (addressed later in this article).
Equations (1)-(3) are the basis for coherent subtraction, keeping in mind that relationships apply only to individual m-sequences. This requirement leads to time multiplexing msequences, addressed in Section II-B.

B. Signal Generation
Unlike previous work, the transmitted signal is viewed as a carrier frequency, phase modulated by time-multiplexed msequences. Using parameters from Part 2 of this article, the length of a single time-multiplexed m-sequence for the following example is L = 2 11 − 1 = 2047 samples. Forty eight m-sequences are time multiplexed to form one time-multiplexed sequence. Specifically, the first sample from each of the 48 sequences is ordered sequentially, 1 to 48, to form the first "digit." Similarly, the second sample from each of the 48 sequences is ordered to form the second digit; and so on until 2047 digits are formed, each with 48 samples. A modulated carrier is illustrated in Fig. 1. The sample rate is 40 kHz, resulting in eight samples per carrier cycle at 5 kHz and six carrier cycles per digit. With this rather unconventional view of m-sequence modulation, each time-multiplexed sequence can be demultiplexed and processed individually, preserving the relationships of (1)- (3).
Phase changes in Fig. 1 occur where the carrier frequency is near zero amplitude. Changing phase near zero carrier amplitude minimizes the transient energy introduced by discontinuities in the m-sequence phase modulation. It also minimizes strain on the acoustic projector that must reproduce the time-domain waveform with minimal distortion. The signal bandwidth is nominally the reciprocal of the digit duration, about 833 Hz. As a final modification to the transmitted (or received) signal, the bandpass signal is filtered to 50% of its initial bandwidth (approximately 415 Hz) to smooth the discontinuity at each m-sequence phase change. Without filtering, interpolation of the m-sequence modulation, which is required for processing Doppler-shifted signals, introduces unacceptable noise. A segment of the transmitted signal, after filtering, is shown in Fig. 2.

C. Modeling the Received Sampled Signal
Each of the received demultiplexed signals is modeled as the sum of delayed and Doppler-shifted replicas of the transmitted m-sequence. Ignoring bandpass filtering and noise and considering a single demultiplexed signal where K is the number of sound-reflecting objects (typically in the hundreds when reverberation is modeled), ω 0 is the carrier radian frequency, and τ k (t) is the time-varying sound propagation delay along the kth source-reflector-receiver path. Doppler is modeled by a Maclaurin series expansion of the time-varying path delay where τ 0k is the initial source-reflector-receiver propagation delay for path k, c 0 is the speed of sound, and ν k is the relative velocity of the moving reflector. When the reflector velocity is constant over the signal duration, the approximation in (5) is exact. Substituting the Maclaurin series approximation into (4) and considering only a single source-reflector-receiver path, the received signal is Switching to discrete notation, the time variable t may be replaced by the product of a sample index n and the time between samples 1/(f 0 · S) seconds, where S is the number of samples per carrier cycle. The time delay from the start of sequence transmission to the initial signal reception at the receiving array τ 0k is determined by the sound speed and the propagation path length. A phase modulated demultiplexed signal with Doppler, received on path k, is denoted as follows: The optimum detector for a signal in white Gaussian noise and with unknown time delay and Doppler is an array of filters each matched to a Doppler-shifted replica of the transmitted signal [19]. It will be convenient to write the Doppler shift v k /c 0 as the product of a matched filter bin index and the bandwidth of each matched filter bin. The bin index M is a positive or negative integer. The matched filter bandwidth for Doppler bin M is the reciprocal of the m-sequence duration. Writing the m-sequence duration in terms of the carrier frequency and the number of carrier cycles Q · L, the Doppler term in (7) is Substituting into (7), the Doppler-shifted demultiplexed and sampled received signal for path k, before filtering and at Doppler frequency index M , is With delay and Doppler set to zero, (8) is the basis for generating the transmitted signal. With appropriate Doppler, delay, and amplitude, (8) is also the basis for simulating the direct arrival, surface and bottom clutter, and moving targets.

D. Demodulation
The carrier frequency is removed from (8) by writing the equation in complex form, i.e., cos x = (e jx + e −jx )/2, multiplying by the complex exponential exp(−j(2πn/S)), and then low-pass filtering to remove the 2ω 0 component. With these operations, the received signal with carrier frequency removed is A simple averaging filter, with averaging over one-half carrier cycle, was used to remove the 2ω 0 component. Our choice of a low-pass filter is somewhat intuitive since one-half carrier cycle is exactly one cycle at 2ω 0 , so the 2ω 0 component is exactly removed by simple averaging. Phase is maintained by forward and reverse filtering (MATLAB (sin(x)/x) 2 "filtfilt" function). The Doppler remaining after demodulation [line 1 of (9)] is removed in subsequent processing.

E. Autocorrelation for Doppler-Shifted Sequences
Earlier we stated that the optimum detector for a signal in white Gaussian noise and with unknown delay and Doppler was a bank of filters each matched to the Doppler-shifted transmitted signal. Although this is correct, matched filters centered at different frequencies have a different number of samples due to Doppler time compression and expansion. Since the autocorrelation properties of m-sequences require exactly L = 2 N − 1 samples, where N is the order of the sequence generating polynomial and L is the sequence length, the usual bank of matched filters is not practical. Matched filtering is therefore implemented by frequency-shifting replicas of the received signal to zero Doppler, resampling to correct each replica for the appropriate Doppler time compression or expansion, and finally correlating each of the time-multiplexed sequences with a zero Doppler reference.
Following (3), correlation is performed in the frequency domain using Parseval's theorem as follows: While the Doppler search computation over S · Q sequences may seem formidable, the autocorrelation properties of the m-sequence are preserved. The autocorrelation for a signal at zero Doppler is shown in Fig. 3(a). The autocorrelation floor is more than 250 dB below the peak value, indicating that the oneto-one correspondence between a point on the autocorrelation function and corresponding time-domain m-sequences has been preserved. The S · Q points on the autocorrelation peak are from S · Q time-multiplexed m-sequences. The autocorrelation of a Doppler-shifted signal, received in frequency bin 12, frequency shifted to bin 0, resampled, then correlated, is shown in Fig. 3(b). Because of computational noise in the resampling interpolation algorithm, the autocorrelation floor for the signal with Doppler is 80 dB below the peak. (While this floor could be reduced by further filtering and a higher sample rate, the 80 dB range is sufficient for the example presented in Part 2.)

F. Doppler Search
An "ambiguity function" is obtained by frequency-shifting and time-correcting multiple replicas of the received signal, one replica for each Doppler. In Fig. 4, the viewing angle has been altered to show the autocorrelation function floor, approximately 80 dB below the autocorrelation peak. The cross-correlation floor is, on average, 10logL below the autocorrelation peak for all Doppler bins. The cross-correlation floor limits detection when high amplitude clutter is present.

G. Removing Clutter
In his dissertation under Prof. T. Birdsall, Chang [18] demonstrated that the one-to-one relationship between a value of the msequence complex autocorrelation function and a time-domain m-sequence could be used to remove clutter at zero Doppler. Because of frequency shifting and time compression/expansion, a different approach is required to remove Doppler-shifted clutter sources. As summarized in (1)-(3), the approach is to select segments of the ambiguity function containing clutter, then  reverse the order of the steps outlined previously. Expanding on our previous explanation, the reverse steps are as follows: 1) considering the correlation function magnitude for each Doppler, identify clutter that exceeds the background by a threshold (typically the background will be a floor that is approximately 10 log L below the peak clutter arrival); 2) in the complex correlation function, zero all nonclutter areas below the threshold, saving only the clutter; 3) inverse transform the clutter correlation to the time domain; 4) transfer the time-domain function to zero Doppler by frequency shifting and resampling; 5) subtract the time-domain clutter from the zero Doppler received signal vector; 6) recompute the ambiguity function using the "cleaned" received signal; 7) repeat the clutter removal process until a target is detected or until the cleaned signal vector approaches ambient noise. While this clutter removal process is computationally intensive, in practice, only a few Doppler bins are expected to contain Doppler-shifted clutter, and the entire range-Doppler plane need not be processed. To illustrate the Doppler removal process, we consider a notional fixed bistatic system in Part 2.

H. Other Considerations
When a segment of the ambiguity function containing clutter is selected, the segment also contains the cross-correlation floors of the reference signal correlated with other clutter and with the target. Typically, these cross-correlation floors are each 10 log L below their respective peaks. Our simulations have shown that sections of these cross-correlations are sufficiently small that they do not detract from coherent subtraction. This is an area for future study.

A. Notional System for Ship Channel Monitoring
The Hong Kong Ship Channel is one of many possible locations for our notional system. As shown in Fig. 5, a promontory provides a physical acoustic barrier between the acoustic projector and the horizontal line array receiver. Separating the source and receiver eliminates the direct acoustic path, which, because of its high energy, can limit the performance of bistatic systems. (Unlike historic "ping-then-listen" sonars, m-sequence transmissions are of low power and long duration. "Listening" continues throughout the signal transmission.) The geometry also places moving clutter sources near perpendicular to the source and receiver paths, minimizing the Doppler channels and array beams that need to be processed.
In Fig. 5, the red lines near the far shore and on the shipping channel boundary indicate regions where acoustic clutter from the channel bottom is expected. Ten clutter points with random amplitudes are generated at both the far shore and the dredged shipping channel boundary within the 5°of the receiving array beam. In addition to bottom clutter from the channel boundary and far shore, reverberation is modeled by 33 clutter points distributed between the intersection of the source and receiver beam patterns and the far shore. Reverberation energy is determined by bottom and surface backscattering coefficients associated with a rocky bottom and by the scattering area. Considering the channel boundary, far shore clutter, and reverberation, 53 clutter sources are included in the received beam.
Two moving vessels are modeled. The outbound ship is traveling at 20 kn with a +15-dB target strength. Doppler-shifted direct and surface-reflected acoustic paths are included, to and from the ship, resulting in four closely spaced acoustic arrivals at the receiver. Similar modeling represents the unknown inbound submerged vessel traveling just beyond the ship channel. The speed of the unknown vessel is 6 kn, target strength -15 dB, and depth 6 m. The objective is to detect the submerged vessel in the presence of reverberation, bottom clutter, and the much larger +15 dB moving ship.
"Simulation" is a word that has been used to describe a variety of systems with widely varying complexity. As shown in Fig. 6, the simulation in this work separates the signal processing and environmental algorithms in such a way that the environmental portion of code can be replaced by the "real" ocean, thereby simplifying the transition from simulation to experiment.
The m-sequence is designed, usually with some trial and error, in the first block. Sequence duration is chosen to exceed the out-and-back bistatic distance. Using the time-multiplexing described in Part 1, there are 2047 digits per sequence with 48 time-multiplexed sequences per digit. Bursts of five timemultiplexed sequences are transmitted. The resulting duration of a single time-multiplexed sequence is over 2 s, long enough for acoustic propagation across the ship channel and back. A carrier frequency of 5 kHz was chosen based on the availability of reasonably priced directional linear acoustic projectors.
In the environmental simulation, sound speed is isovelocity. Doppler is included for refracted and surface-reflected paths for shipping and for the unknown vessel. Ambient noise spectrum level is 70 dB, representative of harbors [20]. The line array beam pattern is modeled as a pie-shaped segment with coherent gain for clutter and incoherent gain for noise.
Simulations with noise, clutter, and shipping are shown in Fig. 7. Axes are source-object-receiver Doppler (+4 to -4 kn), source-object-receiver distance (0-3600 m), and received energy after beamforming and m-sequence processing (10log|A| 2 ). Noise alone is shown in Fig. 7(a). Reverberation (beginning at the intersection of source and receiver area coverage), clutter from the dredged shipping channel, and clutter from the far  shore have been added to noise in Fig. 7(b). In Fig. 7(c), both an outgoing ship, traveling at 15 kn with a + 15 dB target strength, and a submerged inbound underwater vehicle, traveling at 6 kn with a −15 dB target strength, have been added.
The ambiguity function before coherent subtraction is shown in Fig. 8(a). Although reverberation, bottom clutter, and the submerged −15-dB inbound vessel are present, they are masked by the ambiguity function floor generated by the outbound ship. In Fig. 8(b), multiple acoustic arrivals from the +15-dB outbound vessel have been coherently subtracted from the received signal, and the ambiguity function reprocessed. Correlations from the unidentified vessel, bottom clutter, and reverberation remain. After the second iteration [see Fig. 8(c)], boundary clutter residuals have been coherently subtracted, and only the −15-dB 6-kn target remains. The ambiguity function floor is near the "noise only" value in Fig. 7(a), and the −15-dB underwater vessel is easily detectable.

B. Summary
With the widespread availability of recreational and military underwater vehicles, maintaining underwater security in ship channels and harbors is more challenging than ever before. Although experiments are needed, this work demonstrates that coherent subtraction may provide a practical method for reducing clutter in a complex ship channel or harbor environment. Careful placement of the acoustic source and line array receiver can reduce signal processing complexity, minimize equipment cost, and minimize the personnel required for system operation and maintenance. While the focus of this work is on underwater applications, the algorithm may apply in other situations where small targets must be detected in the presence of much larger objects. Active sonar, medical imaging, and ground penetrating radar are among the potential applications.

APPENDIX SONAR EQUATION ANALYSIS FOR A NOTIONAL SYSTEM
A preliminary sonar equation analysis served as the basis for choosing projector power, receiving array beamwidth, signal duration, and m-sequence gain. The analysis is presented in the following material.
Considering the placement of our notional source and receiver, the distance to the farther shore is approximately 1750 m. An m-sequence that avoids data folding would have a minimum duration that allows propagation over twice this distance (out and back). With 1500-m/s sound speed, the sound travel time is approximately 2.4 s. An m-sequence transmitted on a 5000-Hz carrier with 6 cycles/digit and 2047 digits has duration of 2.45 s and allows sonar returns from the opposite shore in a few center beams.
From the sonar equation, the energy returned from the target is the projector source level, less the source-target-receiver propagation loss, plus the target strength, beamformer gain, and m-sequence processing gain. We choose a source level of 180 dB and a triangular geometry with source, target, and receiving array at the vertices of an isosceles triangle with 1500-m sides, again consistent with our notional harbor. One-way propagation loss is assumed 17 log R and target strength −15 dB. An array of 32 elements is assumed with a coherent signal gain of 20 log 32 = Background interference is the energy sum of ambient noise, energy scattered from the surface and bottom, and energy from the projector received in beamformer sidelobes. We estimate each interference component separately.
Considering ambient noise, the noise spectrum level in harbors at 5 kHz is approximately 70 dB. The bandwidth of our signal, after bandpass filtering, is approximately 400 Hz. We assume an incoherent noise gain in the receiving array and an incoherent signal gain. Substituting To determine the scattered energy from the surface and bottom, we determine the scattering area of a single complex sample at the target range and in the target beam. The range extent of a single demultiplexed ambiguity function sample is c 0 /(8f 0 ) and is less than a meter. From [20], the beamwidth of a 32 element (16 wavelength) array at 30 o off broadside is approximately 4 o . This angle must be multiplied by 2 to account for the "back angle" associated with the line array cone of revolution beam pattern. Considering the array beamwidth and the range extent of a single complex sample, the scattering area is less than 1 m 2 .
We assume coherent addition of the scattered energy samples at theN − element hydrophone array, and incoherent addition of L scattered energy samples in the matched filter. We choose −30 dB as our bottom scattering strength estimate, recognizing that directional bistatic scattering strength at low grazing angles is still somewhat of a research issue.
Substituting, our estimate for bottom-scattered energy is Since the backscattering coefficient for the ocean surface is approximately 20 dB less than the bottom coefficient, the contribution to total clutter energy will be neglected. Scattered energy from the surface and bottom, not considering discrete scatterers, should not be a limitation.
Finally, we consider energy from ambiguity function sidelobes. This interference can arise from discrete clutter objects within the receiver beam pattern. Vessels moored near the opposing shore, bottom topography, navigation aids, and moving vessels are examples of objects providing significant sonar target strength. Because of surface currents, even stationary objects can provide Doppler spreading capable of masking a target. As a potentially limiting clutter source, we consider a moving vessel within the navigation channel. Considering the ship channel geometry, a moving vessel with a +15-dB target strength is assumed at a range of 500 m from both the acoustic source and receiving array. The clutter arrival is assumed to be coherent with a signal gain of 20 log L and ambiguity function sidelobes Clearly, our −15-dB target would not be detectable in the presence of surface traffic unless the interference can be reduced by coherent subtraction.