A physical surface
noise model
Billy D. Jones
11-Apr-2008
The spatial cross-correlation function allows one to determine the important measurement scales for the problem of interest. In this section we present a simple physically motivated ocean surface noise model for various multipole sources of interest and integrate to obtain the spatial cross-correlation function for the independent vertical and horizontal orientations.
Keeping with the tenet of the principle of superposition nicely put in the opening line of [1]:
“Ambient noise in deep water can be modeled more or less satisfactorily as a linear superposition of independent plane waves propagating in all directions,”
and noting that in shallow water the principle of superposition still holds just now one needs to carefully add up all of the discrete (normal) and continuous (leaky) modes (see [2]) in the sound channel of interest, the vertical and horizontal two-point spatial cross-correlation functions at frequency f = c/l and hydrophone separation d, for assumed plane wave source Green’s functions from the ocean surface, as nicely developed in chapter 10 of Burdic [3], are given by
and
respectively, where |N|2 is the surface noise angular spectral density whose physical model is described just below and the vertical (horizontal) array element spacing is dz (dr). The angles in these equations can be simply described by
These integration limits are over the upward hemisphere (from
the viewpoint of the array looking up at the ocean surface) most easily thought
about in terms of the elevation angle interval being 
. Note the following relations, useful when integrating the
above formulas. We are using identical notation for the angles as shown in Fig.
10-6 of p. 286 in [3].
The physical surface noise model in terms of a multipole source (m = 0, horizontal impulse; 1, isotropic noise; 2, dipole field; 3, longitudinal quadrupole; 4, longitudinal octupole; etc.) for the assumed azimuthally symmetric noise angular spectral density, |N(f)|2, in the R12 integrands above is given by [3]
where dIH is the received differential intensity of hydrophone H at depth z and horizontal range r, dWH is the differential solid angle subtended at H by a differential annular ring on the ocean surface at elevation angle f. This annular ring is projected onto a surface normal to the direction of f—hence the factor of sin(f). I0 is the source intensity per unit area of surface, and |g(f)|2 is the angular radiation pattern density of the annular ring at elevation angle f.
Quite interestingly, note the geometrical effect of how the r and z dependence exactly cancel in this formula for |N(f)|2. The ocean surface area increase of integrating over all surface elements is exactly canceled by the decrease due to the geometrical spreading of the wave as it propagates towards hydrophone H. This exact cancelation is broken if attenuation is included because more energy is lost as the wave propagates along, however at frequencies below 10 kHz this effect is very small.
The angular radiation pattern density |g(f)|2 of a multipole source at elevation angle f on the ocean surface is physically modeled by [3,4]
where
m = 0: impulsive field at f = 0 (distant shipping, low frequency)
m = 1: monopole field (isotropic noise, low-med freq. transition region)
m = 2: dipole field (surface noise, med frequency)
m = 3: longitudinal quadrupole (fine structure)
m = 4: longitudinal octupole (hyperfine structure)
etc.
Given this angular radiation pattern |g(f)|2, the noise angular spectral density |N(f)|2, defined above, is shown in the next figure for the first five multipoles. Note the m = 0 low frequency mode being peaked in the horizontal direction relevant especially for long-range low frequency propagation and also continental slope propagation towards deeper depths [Ref. 4, p. 515] where disappearance of the horizontal noise notch has been noted.
In order to incorporate noise notch effects (continuous modes dropping out for horizontal propagation; see p. 500 of [4]) it is useful to introduce an impulsive noise angular spectral density given by a Dirac delta function at an arbitrary elevation f = f0:
In Appendix A we show that this impulsive angular density at f0 = 0 (i.e. horizontal) is the same as the above multipole mode at m = 0, an original result that shows all of Burdic’s factors of 4p etc. are right on.
Tidying up the notation, it is convenient to normalize the spatial cross-correlation functions by the average ambient noise intensity received at the center of the array. Also, to include beamforming steering effects, a possible delay t between channels is allowed. Thus we define
where 
is the real part of
the corresponding expression. Substituting the formula for |g(f)|2 into the formula for |N(f)|2 and then integrating the R12
formulas above, leaves the analytic normalized single frequency vertical and
horizontal spatial cross-correlation functions r12.
These analytic formulas for r12
are placed in Appendix A. Just below we plot r12
using these formulas for various hydrophone separations and beam-steering delays
t. The fact that we can obtain the
formulas for r12 in closed analytic
form is one of the primary advantages of this simple but physical surface noise
model. Analytic work aids understanding.
Given the surface noise angular spectra in the figure above, the vertical and horizontal cross-correlation functions are as follows in the next four figures just below.
Note the ordering between the multipoles is opposite for the vertical and horizontal array cases. Given a multipole source at the surface, in the vertical array case, for higher multipoles, the first zero of the cross-correlation function is smaller (less correlated) but the oscillations above dz ~ l (higher frequency effect) are larger (more correlated). The exact opposite is seen for the horizontal array case. In general at zero delay (t = 0), for a source at the surface, the overall picture is that the horizontal array has larger correlation lengths (first zero crossing) with smaller high frequency oscillations in the cross-correlation function, whereas the vertical array has smaller correlation lengths with larger high frequency oscillations.
Beam steering (nonzero t)
can reverse these conclusions as is well known. The third figure just below for
the vertical array with nonzero t shows
a dipole changing into a lower order multipole (larger correlation length) by
steering the beam to look up towards the source, and vice versa when it looks
down away from the source. This makes intuitive sense in that if the source is
in the broadside beam then the correlations will be larger, and if it is not in
the broadside beam then if you steer the beam towards the source, this should
increase the correlation length. Beam steering for the horizontal array is not
as dramatic, but if the beam is steered to endfire then the correlations should
decrease (since this direction is perpendicular to the direction to the
source), and they do as measured by the correlation length in the fourth figure
just below.
Beamforming of the surface noise model
Given these noise model results it is of interest to apply them to a realistic vertical array such as is being developed for this project. The goal is to demonstrate the feasibility. The array gain for an eight element vertical line array is developed in this section given the surface noise dipole field vertical cross-correlation function derived in the previous section.
In general, as developed in chapters 11 and 14 of [3], the average
output noise power of a vertical line array with M discrete elements 
given the noise
angular spectral density of the last section |N(f,g)|2
with elevation angle f and azimuthal
angle g is given by
where GM is the array pattern function which is the Fourier transform of the normalized array aperture function gM (written here for M discrete even numbered elements):
Azimuthal symmetry has been assumed and f0 is the steering elevation angle
for an assumed incoming plane wave with frequency f = c/l and elevation angle f.
The average output noise power for a single
omnidirectional sensor 
is the same formula
with GM set to unity:
The array gain in the signal-to-noise ratio (SNR) of the system—in this case for the extra noise that can be nulled out by beam steering to elevation angle f0 given an incident unit amplitude plane wave at elevation angle f—is therefore
Before integrating we enlarge our dipole surface noise field model to include the bottom image as well as direct arrival signal (constants that cancel in the array gain ratio are not written in this formula):
where Rb is the plane wave pressure reflection coefficient, given by the well-known Rayleigh formula (although Frisk [2] on p. 40 notes that it was originally derived by George Green of Green’s function fame):
where 
is the seabed/seawater
density ratio and 
is the index of
refraction of the bottom with 
being complex to
include seabed attenuation (
, with 
the seabed attenuation
in nepers per unit distance). Except near the critical or intromission angles, the
Rb amplitude is nearly a constant as shown by the next two figures
for typical [5] fast and slow seabeds respectively. The blue curves being Rb
amplitude, gold curves bottom loss (inverse amplitude in dB), and pink
curves Rb phase for completeness. The figures are drawn at the same
scale for clarity. Note how in this typical example the slow bottom is lossier
than the fast bottom over all elevation angles f.
Since the blue curves are nearly constant except near critical and intromission
angles we will assume |Rb| is constant in the f integrations for array gain in the simple
model below.
Note that a general reflection coefficient with attenuation in
the bottom is approximately invariant under frequency scalings—change frequency
and the reflection coefficient (bottom loss) does not change. This scaling is
exact if the attenuation coefficient ab
scales exactly linearly with
frequency, i.e. iff
. In practice, for sediment attenuation, this scaling is
found to hold fairly accurately over an amazingly wide range of frequencies
2.5-400 kHz. See Fig. 5.27 on p. 167 of [5] which is a log-log plot that shows . A final comment is that the attenuation k shown in these two figures below has units
of dB/m/kHz, and if this scaling being discussed here is exact, this translates
to a constant k (but the constant changes
depending on the bottom type). The relationship between ab in nepers/m and k
in dB/m/kHz is simple and well known:
where the assumed unit of the respective quantity in this formula is shown in parentheses.
Now we perform the vertical noise power integration for a vertical line array with M even-numbered discrete elements. As stated above, we assume azimuthal symmetry and use a physical dipole-field surface noise model that includes direct and bottom bounce paths. Substituting in the relations above, integrating over z, and simplifying gives
where the sin(f) factor on the far right comes from the
dipole-field noise model and the cos(f)
comes from the solid angle measure dW. Note
that the GM sum was done in closed form. Also note that this final form is the same whether M is even
or odd. This is useful to know in comparing results with Burdic [3] p. 324,
because Burdic uses an odd number of elements (M = 9) whereas our array has an
even number of elements (M = 8). This final form above agrees with Burdic’s
result[i]
(Burdic’s 
is our; our ‘
’ is for ‘vertical’, not frequency—frequency dependence is
implicit in our notation).
The result for the noise power at a single omni phone is
simple. Using the formula for 
above we have
Taking the ratio of these previous two results and subsequent
logarithm as in the array gain formula above, setting M = 8 as for our array,
and choosing the array design frequency (d = l/2,
green curve) and frequencies an octave below (d = l/4, red curve) and above (d = l,
blue curve) gives the following three colored curves in the left figure below. For
comparison the 100-element result is shown in the right figure at the same
scale. These are polar plots of array gain in dB over elevation steering angle f0. The dashed curve is the
surface noise (with bottom bounce) dipole field angular spectrum 
plotted in relative dB
for understanding.
We see the largest gain is near horizontal steering (f0 = 0) for all three frequencies,
where the noise is quietest. Note however (1) the high frequency (blue) curve
has more of a noise notch than the other two, which in this case is a spatial
processing artifact, and (2) if one picks the curve with largest array gain
across all steering angles f0 
[-p/2, p/2] one
sees high frequencies (blue curves) better at positive (upward looking)
elevations, design frequencies (green curves) better near horizontal, and low
frequencies (red curves) better at aft endfire near f0 = -p/2. This shows how adaptive beamforming can be
made to work—the analogy is precise since the phase in the steering vector is
directly proportional to frequency which is what is being adapted as the colors
change in these curves.
A final note is that the high frequency blue curves scale roughly like the well-known 10 log(M) formula when comparing the 8 and 100 element results. This gives an approximately 10 log(100/8) ~ 11 dB gain in AG across all steering angles. However, for the design frequency green curves, near horizontal steering, this gain in AG is almost twice that at 20 dB.
8-element VLA (same scale) 100-element VLA
Polar plot array gain (dB) vs. elevation steering angle for three frequencies
Theory
Given the principles of superposition and energy conservation we can go a long ways. Thus we emphasize the plane wave reflection coefficient (whose information is a nice summary of energy conservation) and also discuss its inversion from pressure measurements of a point source in a homogeneous fluid layer bounded by arbitrarily horizontally stratified media [2] for understanding and also a possible alternative/comparative method.
Energy conservation
First, energy conservation: for an arbitrary acoustic plane wave traveling in a homogeneous fluid (water) at grazing angle f (f positive is down here; it is the negative of the elevation angle, however in this section we use the same symbol f as used for elevation in the other sections of this paper) with constant density r and sound speed c impinging on a flat boundary of a second fluid (modeled bottom) with density rb, sound speed cb, and attenuation ab, energy is conserved in the time-averaged sense at the instant when the wave strikes the boundary and splits into a reflected and refracted (transmitted) wave. After that, given attenuation in the transmitted wave, extra energy is lost to the bottom. However, the reflection coefficient has all that information in it and knows how much energy will be lost to the bottom even before it has been lost!
We say “energy conservation” but here in all its jargon by this is really meant:
In the direction normal to the flat surface (lets say “z”), the magnitude of the time average of the acoustic energy flux in the incident wave is equal to the sum of the magnitudes of the acoustic energy flux in the reflected and transmitted waves.
In an equation (since energy flux is energy per unit time per unit area, i.e. intensity, we use the symbol I) energy conservation here is
Consult Frisk [2] for example or any ocean acoustics (or quantum mechanics or electromagnetism) textbook for the setup. The incident, reflected, and transmitted plane waves are written out with Rb being the pressure reflection coefficient and Tb = 1 + Rb (continuity of pressure gave this relation) being the pressure transmission coefficient. Continuity of pressure and normal acoustic velocity are demanded on physical grounds, and then the above equation for Iz conservation can now be written as follows.
We have generalized Frisk’s discussion by including
attenuation in the bottom, ab.
Here, z > 0 is down into the
bottom, so we see the exponential factor is decaying as long as 
, which is how one can see that the relation is ill-defined if attenuation in the water
is included and 
. In practice this is not satisfied so things are
well-defined because a/ab ~ 0.001. In the plots below we
have a nonzero a (seawater
attenuation). This is added to the above formulas by taking 
.
Now we ask whether Eq. (1) above holds for nonzero ab or not. Given the decaying
exponential factor we see that it does not for nonzero ab unless z = 0. This makes sense because heat is
lost to the bottom as it propagates with nonzero ab.
Energy is always conserved, it just goes to different places. Also, with no
attenuation in the bottom, i.e. ab
= 0, energy flux is conserved for all z (normal to the direction of the boundary)
except for angles such that total
internal reflection is obeyed, because then the evanescent lateral wave leads
to an imaginary 
term and energy is lost to the bottom in a similar way as
discussed for the nonzero ab case.
All is well; the equation has all the information in it.
So, with nonzero ab, we now show that Eq. (1) does indeed exactly hold at z = 0 (the instant the wave impacts with the boundary and splits in two) over all ab and grazing angles f. We show this with a couple plots. The fact that the pink lines are at unity is a statement of conservation of energy in this notation. The blue curve is actually the plane wave intensity reflection coefficient. The dense‑slow and dense‑fast bottom parameters (except for a and ab) are the same as for the reflection coefficients plotted above. Interestingly enough notice how the dense‑slow and dense‑fast results look rather different except for the case with a very large attenuation in the bottom; then they are not identical, but they are similar. Attenuation smooths out the critical and intromission angles.
The point: bottom loss BL = -20 log |Rb| is very general and tells us how much energy gets taken out of the incident wave and transmitted into the bottom. Also, given bottom loss measurements, one can invert these equations and, as a function of grazing angle, gain information about the amount of energy flux transmitted into the bottom (if any) and also the magnitude of the bottom parameters rb, cb, and ab. Hence the emphasis on measuring bottom loss as a function of grazing angle in this project.
Inversion of a
point source in a homogeneous fluid with arbitrarily bounded flat boundaries
As alluded to, given the principle of superposition, we start with a point source in a homogeneous medium and work from there. This discussion follows the work of Frisk in §6.3 of [2] and Frisk et al. in [6] where they introduce the technique and derive the following important result for inhomogeneous media (my bolding)
“The key point is that the effect of the boundaries can be incorporated into the theory in an exact manner using plane-wave reflection coefficients even for the case of an inhomogeneous ocean.”
With this as motivation, we now write the general Green’s function for the case of a point source in a homogeneous ocean (constant density r and sound speed c) bounded by flat but other than that arbitrary media. Let the plane-wave surface and bottom pressure reflection coefficient (of Rayleigh form) be given by RS and RB respectively. Given this, Frisk derives the following general relation. Keep in mind his conventions have z > 0 being down from the ocean surface, water depth is given by h, the source is at (0, z0) and receiver at (r, z). We keep his exact notation:
In this equation the horizontal wavenumber is kr
which is the integration variable, and the vertical wavenumber is kz—actually
kz in this equation can be thought of as a shorthand for 
which of course has to
be integrated over. The beauty of this chapter is Frisk then proceeds to (1)
expand the denominator in a geometrical series and show how it exactly leads to
the infinite ray theory solution in this environment, and then (2) in
the full form of equation 2 sets RS ® -1, c1 >
c (fast bottom), and r1 >
r, the famous Pekeris model, and
proceeds to do a contour integration in the complex kr plane showing
how the discrete normal modes come out exactly, and how the branch line
integral leads to the continuous improper modes which in an asymptotic
expansion valid at high frequencies become an infinite sum of leaky modes,
with a formula looking very similar to the normal mode solution, except that
the leaky modes are multiplied by exponentials that decay with range—hence their
name. They are thus rightly so de‑emphasized in long range propagation
problems. Quite a chapter!
So it is all there (in equation 2 above), but our point here is to invert the equation as Frisk et al. do in [6] except we will not assume that the incident and reflected wave are separable, because this will not hold for a distributed source such as ocean surface noise. The algebra is simple, Hankel transforming equation 2 above leaves
Assuming surface noise is a sheet of point sources just below a pressure release sea surface (a good approximation except at very high sea states and/or extremely high frequencies) allows us to set RS ® -1 leaving the final form
Note that the grazing angle and frequency dependence in this formula for the pressure reflection coefficient, RB, is hidden in the relations
Given the exactness behind the fundamentals of this
equation, the alternative method idea is to measure pressure across a near
horizontal glider path (constant z) and obtain many ensembles of p(r; z, z0)
data for a source at z0 very near the surface (a parameter that can
be fit for). Then perform the required Hankel transform with well-known efficient
numerical projection techniques [7] for the given ensemble. Repeat and finally
perform an average over all ensembles leaving a final estimate 
. Given this, bottom parameters can be inverted for as
discussed above, or also the average modal attenuation coefficient follows
directly form BL according to Kornhauser and Raney’s classic formula [8, 9]:
In this formula, SL is surface loss, which vanishes if RS = -1 as assumed above. However, if RS is known or if it is fit for in the same way as for source depth z0 discussed above, then it might be useful in this formula. Dc is the cycle distance which includes ray displacement if known and/or measured. This is where the phase of the reflection coefficient can be measured:
Transmission loss also has a bottom loss term that can be used to compare with separate TL measurements as in Arvelo [10]:
Where Nc is the number of complete cycles which does not necessarily have to include a bounce.
Ambient noise inversion:
a simple method with a littoral glider
A field-tested simple robust method for ambient noise
inversion using conservation of energy principles is discussed in Arvelo’s
recent paper “Robustness and constraints of ambient noise inversion” [10]. His
work is based on Harrison and Simon’s paper “Geoacoustic inversion of ambient
noise: A simple method” [11] from
The method is based on an estimate of the surface noise spatial cross-correlation function discussed in the first section of the current paper. Given our array is an eight element VLA we will specialize to the vertical cross-correlation function which sets Arvelo’s array tilt angle to 90°. From Arvelo’s paper (using the notation in the first section of the current paper for continuity) in Harrison and Simon’s ray-based derivation we have
where the cross-correlation
function has been evaluated at the origin to estimate the average ouput ambient
noise power of a single omnidirectional sensor (beamforming is discussed next).
As discussed in Arvelo’s paper, in the above equation, the path length of the
ray from the source at the surface directly to the receiver (direct path) is 
, and the path length of the ray from the source and one
bounce off the bottom and then to the receiver (bottom bounce) is 
. Each factor of 
(and the direct path
incident wave) gets an attenuation factor 
and each factor of 
gets an attenuation
factor 
. These path lengths are approximately (exact for
straight-line propagation)
with the source at 
(usually set to zero
for surface noise applications, but in 100 m of water it is possible to have a
1% effect hidden here), the receiver at z, and water depth given by h as in the
previous section. The angle the bottom bounce ray makes with the seabed is 
, and the direct path elevation angle from receiver to source
is 
.
Now we include beamforming as above in the current paper, with a modification according to the above Arvelo formula because it uses point-source Green’s functions instead of simple plane waves which we used above. Also we integrate over all elevation angles to account for bottom as well as surface arrival paths as we did in the beamforming examples above. The noise power array output becomes
where GM is the pattern function of a discrete M-element VLA as discussed in previous sections. Then the upward and downward steered beam powers are subtracted according to [10, 11]
Recall [2] 
~ 0.0001 << 1,
but at 40 kHz (phase I array design frequency), and after 200 m of a worst case
scenario in 100 m of water gives 
, and now this term cannot be neglected. That is why Arvelo
emphasizes making measurements where 
vanishes, or at least
is very small, i.e. put the receiver (glider) near the bottom.
This formula in the blue box is simply a statement of conservation of energy [11] in that no matter how complex the environment or waveguide features, the ratio of the upward and downward moving energy flux locally measured is the power reflection coefficient by definition. Acoustic energy flux principles are discussed at length in the theory section of this current paper.
This final equation boxed in blue is the main result of our measurement method discussed at length in [10] with simulated examples for a robustness/sensitivity analysis, and also real-world East China Sea ASIAEX 2001 data with promising results of BL and TL matching well with independent high-fidelity geoacoustic inversion results [14]. We list here the three main lessons learned in [10] since Arvelo is on the team of this project, and his methods/lessons-learned are being used with the glider platform-mounted 8‑element PVDF vertical array discussed in earlier status reports. The three main lessons:
1. Array element location must be known within one-fifteenth of an acoustic wavelength.
2. Electronic and flow noise must be more than 10 dB below the omnidirectional ambient noise level across the array.
3. Nearby loud ship interference is particularly damaging to the accuracy of the bottom loss estimation.
The first point especially will
be challenging with a glider where the underwater horizontal position is not
accurately known (it is known within a meter or better), however if one does
observe Arvelo’s plots upon which this conclusion was made (less than 1 dB of
error in bottom loss), we see the 3 cm result scales according to a “better
than 
” rule, with the BL error being at most 3 dB. This is
workable especially with an improved glider control board which is under development.
Appendix A
The results of the analytic integrations of the formulas for r12 defined in the first section of the body of the paper are listed in this appendix. They were obtained using Mathematica v. 6.0.1.0. Jn(x) is indeed the nth-order Bessel function of the first kind. From these results we see the 0‑pole field and impulsive noise field at f0 = 0 are equivalent as mentioned in the body of the paper.
In integrating the horizontal array formulas note from the
above useful relations that sin(f) =
cos(y) cos(q). Also, the m = 0 formulas are actually singular, but a
well-defined limiting procedure has been found (an original result of this
paper) which sets m ® e,
a small positive infinitesimal number, performs the analytic calculations, and
then at the end of the day takes the limit e
®
0. For this special m = 0 case, singular e-dependent
multiplicative factors (actually just a simple pole, i.e. 1/e) cancel in the numerator and denominator of
the formula for 
, leaving a finite well-defined limit as e ® 0:
, seen to be the same as 
in the impulsive noise model above. Recap: horizontal impulsive noise and the 0-pole
field at the surface are equivalent.
