Measurements of high-frequency acoustic scattering from glacially-eroded rock outcrops

5 MB·

by  Derek R. Olson

Measurements of high-frequency acoustic scattering from glacially-eroded
rock outcropsa)
Derek R. Olson,1 Anthony P. Lyons,2 and Torstein O. Sæbø3
1)Applied Research Laboratory, Pennsylvania State University, State College, PA,
16804
2)University of New Hampshire, Durham, NH 03824
3)Norwegian Defence Research Establishment, Kjeller N-2027, Norway
(Dated: 15 April 2016)
Measurementsofacousticbackscatteringfromglacially-erodedrockoutcropsweremadeoffthecoastofSande-
fjord, Norway using a high-frequency synthetic aperture sonar (SAS) system. A method by which scattering
6 strengthcanbeestimatedfromdatacollectedbyaSASsystemisdetailed,aswellasamethodtoestimatean
1 effective calibration parameter for the system. Scattering strength measurements from very smooth areas of
0 the rock outcrops agree with predictions from both the small-slope approximation and perturbation theory,
2 and range between -33 and -26 dB at 20◦ grazing angle. Scattering strength measurements from very rough

r areas of the rock outcrops agree with the sine-squared shape of the empirical Lambertian model and fall
p
between -30 and -20 dB at 20◦ grazing angle. Both perturbation theory and the small-slope approximation
A
are expected to be inaccurate for the very rough area, and overestimate scattering strength by 8 dB or more

4 for all measurements of very rough surfaces. Supporting characterization of the environment was performed
1 in the form of geoacoustic and roughness parameter estimates.

] PACS numbers: 43.30.Hw, 43.30.Gv, 43.30.Fn, 43.30.Pc
h
p
-
o I. INTRODUCTION McKinneyandAnderson(1964)at100kHz,andSoukup
e and Gragg (2003) at 2-3.5 kHz. A report from the
g
. Determinationoftherelationshipbetweentheacoustic Applied Physics Laboratory, University of Washington
s field scattered by the seafloor and environmental param- (APL-UW, 1994) presents model curves that were fit to
c
i eters is crucial to understanding and predicting acous- scattering strength measurements taken at a rocky site.
s
y tic interaction with the ocean environment. An impor- Previous measurements of scattering strength fall be-
h tant step in this process is to perform measurements of tween -15 to -22 dB at 20◦ grazing angle, with the ex-
p seafloor scattering in conjunction with measurements of ceptionoftheAPL-UWmeasurementfrom‘roughrock’,
[ seafloor geoacoustic and roughness properties. The scat- which is approximately -8 dB at the same angle. These

tered field is typically characterized by the differential measurementstendtodecreasemonotonicallywithgraz-
3
v scattering cross section per unit area per unit solid an- ing angle, apart from typical statistical fluctuations, and
8 gle,σ,whichwillbeabbreviatedhereas‘scatteringcross some systematic ripples in Soukup and Gragg (2003).
7 section’ or ‘cross section,’ keeping in mind that it is di- Some of these measurements likely suffer from bias. The
0 mensionless. It is a system-independent quantity that measurements by Eyring et al. (1948) likely include the
6 characterizes the angular and frequency dependence of effect of multiple interactions with the sea surface and
0
the second moment of the acoustic pressure field due to floor, and may represent an overestimate of scattering
.
1 scattering(JacksonandRichardson,2007; Pierce,1994). strength. One of the of the authors of the APL-UW re-
0 In terms of acoustic scattering measurements, rock port (Jackson, 2015) has expressed concerns regarding
6
seafloors have received little attention to date, with five the reliability of the calibration for the measurements
1
existing scattering strength measurements reported in on which the models were based. The model curves
:
v theliterature. Totheauthors’knowledge,detailedacous- from the APL-UW report exceed the maximum possible
Xi tic scattering measurements of rock seafloors, coupled Lambert’s law curve, and would violate energy conser-
withmeasuredgroundtruthhaveneverbeenmade. Scat- vation unless the scattering cross section is azimuthally
r
a tering strength measurements of rock seafloors without anisotropic or if enhanced backscattering (Ishimaru and
quantitative geophysical parameters have been made by Chen, 1990; Thorsos and Jackson, 1991) were present.
Eyring et al. (1948) at 24 kHz, Urick (1954) at 55 kHz,
The present work addresses the paucity of scatter-
ing measurements from rock seafloors by presenting es-
timates of scattering strength obtained from glacially-
eroded rock outcrops, accompanied by characterization
a)Copyright(2016)AcousticalSocietyofAmerica. Aversionofthis of geoacoustic and roughness properties. These out-
articlehasbeenpublishedinJ.Acoust. Soc. Am. 139,1833(2016)
crops contain two contrasting roughness characteristics
and may be found at http://link.aip.org/link/?JAS/139/1833.
This article may be downloaded for personal use only. Any other that allow model-data comparisons to be made under
userequirespriorpermissionoftheauthorandtheAcousticalSo- differentconditions. Acousticbackscatteringdataat100
cietyofAmerica. kHz were collected off the coast of Sandefjord, Norway
2
bytheNorwegianDefenceResearchEstablishment(FFI) i is the imaginary unit (Jackson and Richardson, 2007,
aboard the HU Sverdrup II using the HISAS 1030 syn- Chap. 9). The bulk and shear moduli, K and µ re-
thetic aperture sonar (SAS) system from a HUGIN au- spectively, are related to c˜ and c˜ through the standard
p t
tonomousvehicle(Midtgaardet al.,2011; Fossumet al., formulae: c˜ = (cid:113)(K+ 4µ)/ρ and c˜ = (cid:112)µ/ρ (Mavko
2008). This sonar has not been calibrated in terms of p 3 b t b
et al., 2003).
itsreceiversensitivitys orsourcestrengths ,whichare
r 0
required for scattering strength estimates. The product Geoacousticparametersoftheareawerenotmeasured,
of these two parameters was estimated by comparison of butboundswerecomputedbyusinganeffectivemedium
measured data to a model which used measured input approximation combined with previously measured min-
parameters. Roughness estimates of rock outcrops were eral compositions (Neumann, 1976, 1980). Crystalline
obtained using a digital stereo photogrammetry system. igneous rock is composed of randomly-oriented crystal
Geoacousticparameterswereestimatedusinganeffective grainsofindividualminerals(LeMaitre,2005). Aneffec-
mediummodelwithpreviouslymeasuredmineralcompo- tive medium model is used to replace the heterogeneous
sition of the bedrock in the area. granular material with a homogeneous material with
Characterization of the environment, including esti- properties that reflect the aggregate effect of the crystal
mates of geoacoustic and roughness parameters of the grains. Thenarrowestboundsfortheaggregatebulkand
bedrock is summarized in Section II. Section III gives an shear moduli without knowledge of grain shapes or dis-
overview of the acoustic scattering experiment, and de- tributionsareattainedbytheHashin-Shtrikman-Walpole
tailsthedataprocessingandcalibrationtechnique. Scat- (HSW) bounds for multiphase composites (Mavko et al.,
teringstrengthresultsarepresentedanddiscussedinSec- 2003). For each mineral component, its bulk modulus,
tion IV, with conclusions given in Section V. shear modulus, density, and volume fraction βi are re-
quired. The two isotropic moduli are computed from
anisotropic crystalline mineral elastic properties using a
II. ENVIRONMENTAL CHARACTERIZATION Voigt average (den Toonder et al., 1999) and the bulk
density is computed using a simple volume average.
A. Geoacoustic characterization Due to the slight porosity of crystalline igneous rock,
one of the components of the effective medium is water,
the volume fraction of which was not measured by Neu-
ThebedrocksurroundingLarvikand Sandefjord, Nor-
mann(1976,1980). Measurementsofporosityingranite,
way and their coastline is composed of monzonite, a
asimilarcrystallineigneousrock,weremadebyTullborg
crystalline intrusive igneous rock (Petersen, 1978; Neu-
and Larson (2006) in Sweden and by Norton and Knapp
mann, 1980; Le Maitre, 2005). This material supports
(1977) in the United States. These measurements range
both compressional and shear waves, and like most nat-
between 0.61% and 2.60% with a mean of 1.02% and a
ural materials, contains intrinsic dispersion and atten-
standard deviation of 0.43%. The porosity mean is used
uation (Mavko et al., 2003). Wave propagation within
to compute the mean wave speeds and the bulk density,
monzonite has been modeled as an elastic medium with
and the standard deviation is used to compute model
frequency-independentcomplexwavespeeds. Thismodel
uncertainties.
contains transverse and longitudinal waves, but is dis-
persionless with linear dependence of attenuation on fre- Attenuationat100kHzcannotbeestimatedusingthe
quency. Since linear attenuation over a limited fre- HSW bounds. Instead, an estimate for attenuation is
quency range results in logarithmic dispersion via the based on measurements of water-saturated granite mea-
Kramers-Kronig relations it does not satisfy causality sured by Coyner and Martin (1990). Using the resonant
(Futterman, 1962; Milton et al., 1997). Although the bar technique, the measured quality factor at 100 kHz,
model is acausal, it is commonly used to characterize the center frequency of the HISAS 1030, was approxi-
the seafloor especially when little information is avail- mately 30. This estimate results in δp ≈ 0.02. Shear
able (Jackson and Richardson, 2007, Chap. 9). The wave attenuation in saturated granite was found by the
frequency-dependence of interface scattering is primar- same researchers to be approximately twice that of com-
ily dependent on the roughness spectrum, and is weakly pressional wave attenuation in this frequency range, so
dependentonsoundspeeddispersion. Thusincludingap- δt =2δp.
propriatedispersionintheelasticmodelwouldlikelynot The isotropic averaged elastic moduli (K and µ), den-
be observable when comparing scattering measurements sity (ρ), and volume fraction (β) for each mineral in the
to models. Sandefjord area can be found in Table I, with the results
Parametersfortheelasticmodelarec˜ andc˜,thecom- of the Voigt average for isotropic bulk and shear moduli,
p t
pressional and shear phase speeds, δ and δ , the com- anddensity. InputstotheVoigtaverageforeachmineral
p t
pressional and shear loss parameters (imaginary compo- are from Gebrande (1982). These moduli were used to
nent of the complex sound speed divided by the phase estimate HSW bounds on compressional velocity, shear
speed), and ρ , the bulk density. The complex wave velocity and density, and are reported in Table II, along
b
speeds, c and c , are related to the phase speeds by with the attenuation parameters. Since detailed infor-
p t
c = c˜ (1+iδ ) and c = c˜(1+iδ ) respectively, where mation regarding mineral composition is rarely available
p p p t t t
3
for a given region, Table II also includes tabulated pa- TABLEI.Mineralcompositionoftheexperimentalsite,along
with their volume fractions, β, isotropic bulk modulus, K,
rameters for generic rock from Table IV.2 in APL-UW
isotropic shear modulus, µ, and density, ρ. Mineral fractions
(1994), and generic granite (a similar crystalline igneous
havebeenaveragedfromthreesitesinNeumann(1976,1980),
rock)fromTable5.2inBourbi´eetal.(1987). Thegranite
and elastic constants used to estimate isotropic moduli are
parameters are the mean of ranges reported in Bourbi´e
from Gebrande (1982).
et al. (1987), which are of 4500 - 6000 m/s for c , 2500
p
- 3300 for ct, and 2500 - 2700 kg/m3 for ρb. The pa- Mineral Name K [GPa] µ [GPa] ρ [kg /m3] β [%]
rameters for generic rock generally perform quite poorly, Water 2.20 0.00 1000 1.20
although parameters for generic granite are much closer Albite 79.17 22.14 2630 51.05
to the estimates using detailed minerology. This com- Anorthite 110.84 42.88 2760 8.97
parison indicates that tabulated values for geoacoustic Orthoclase 81.71 31.85 2560 23.91
properties may be adequate for modeling purposes if the Nepheline 80.43 32.86 2620 3.32
lithology is known. Diopside 166.31 65.77 3310 3.45
Enstatite 184.70 75.99 3200 0.08
Olivine 129.60 81.05 3224 2.88
Magnetite 161.00 91.30 5180 2.05
B. Roughness characterization Apatite 125.80 49.65 3218 1.35
An experiment to characterize the roughness of rock
TABLEII.Summaryofgeoacousticparameterestimatesand
outcrops was performed in May 2013 near Sandefjord,
their bounds. The values for c˜ , c˜ and ρ were calculated
Norway at 59◦4’26.2”N, 10◦15’42.1”E. Roughness mea- p t b
using the Hashin-Shtrikman-Walpole bounds (Mavko et al.,
surements of in-air (subaerial) glacially-eroded rock out-
2003). AlsoincludedareestimatesforgenericrockfromTable
croppings called roches moutone´es were obtained using IV.2inAPL-UW(1994),andaveragevaluesforgenericgran-
digital stereo photogrammetry. These roughness mea- ite(whichissimilartomonzonite)fromTable5.2inBourbi´e
surements provided inputs to the effective acoustic sys- et al. (1987)
tem calibration described below, and inputs to approx-
imate scattering models that were compared with mea- Lower Upper Generic Generic
Parameter Mean
sured data. These types of outcrops have two contrast- Bound Bound Rock Granite
ing roughness characteristics: a gently-undulating sur- c˜ [m/s] 5945 6842 6393 3600 5300
p
face where the ice flowed onto the outcrop (stoss), and c˜t [m/s] 3198 3353 3276 1900 2900
a stepped surface where the ice flowed off of the outcrop ρb [kg/m3] 2696 2720 2708 2500 2600
(lee). The stoss side has been shaped through the mech- δp 0.02 0.02 0.02 0.0018 0.01
δ 0.04 0.04 0.04 0.085 0.05
anism of glacial clast abrasion, whereby sediment grains t
(clasts) trapped beneath glaciers gouge and scrape the
underlying bedrock (Scholz and Engelder, 1976; Alley
et al., 1997). The resulting surface is characterized by
large-scale undulations that follow the glacial flow pat-
tern, and small-scale scratches, or striae, from individ-
ual clasts or sediment grains. The stepped leeward side
is formed when hydraulic fracturing dislodges blocks of
rock delineated by the internal joint structure, a process
termedglacialquarryingorplucking(Hallet,1996; Iver-
son, 2012; Zoet et al., 2013).
Roughness measurements are reported from two in-air
surfaces that contain glacial abrasion and plucking re-
spectively. They were chosen because they are represen-
tative of other surfaces with the same geomorphology in
the region. To illustrate the roughness characteristics of
thesesurfaces,photographsareshowninFig.1forglacial
abrasion, and Fig. 2 for glacial plucking. These surfaces
are the in-air expression of the same geomorphological
featuresstudiedbelowinSectionIII.Similaritieswiththe FIG. 1. (color online) Photograph of in-air (subaerial)
SASimageryofthesefeaturesisexploredbelow. Theice glacially abraded roughness measured in air. Ice flow in this
flowed from North to South at the roughness measure- area was from the bottom left of the image to the top right.
ment site (Mangerud, 2004). Additional measurements The large-scale undulations can be seen, and are represented
and analysis can be found in Olson (2014). intheheightfieldsinFig.3. Thedistancebetweenthesurface
shown in the bottom left and upper right is approximately
In this research, acoustic measurements were per-
2.5 m. The diagonal line in the upper part of the image is a
formed on submerged outcrops, which preserved the
shadow cast by the measurement system.
glacial features. Roughness measurements were per-
4
is the resolution perpendicular to baseline in the hori-
zontal plane, and ∆z is the precision of the depth esti-
mate. Since a 7×7 correlation window was used for the
stereo-matching algorithm in the photogrammetry pro-
cessing, a conservative estimate of the realistic image-
plane resolutions is seven times the nominal resolution,
(1.12,1.11)mmand(2.49,2.45)mmforthehigh-andlow-
resolution modes respectively. Surface features as small
asthenominalsystemresolutionareobservable,butfea-
tures at scales smaller than the realistic resolution are
subject to an effective low-pass filter due to the correla-
tion window. These resolutions correspond to Nyquist-
Shannon spatial frequencies of approximately 449 m−1
and 202 m−1 for the high- and low-resolution modes re-
spectively. The cameras were calibrated using a black
andwhitecheckerboardpatternattachedtoaglassplate.
Images were processed to obtain height fields using the
OpenCV library (Kaehler and Bradski, 2015), and the
intrinsic and extrinsic camera parameters were obtained
using the camera calibration toolbox (Bouguet, 2008).
DetailsonstereophotogrammetrycanbefoundinWong
(1980) and an example of an underwater system can be
found in Lyons and Pouliquen (2004).
FIG.2. (coloronline)Photographofglaciallypluckedrough-
ness. Ice flow in this area is out of the page. The typical
stepped characteristics of the leeward side of a roche mou- 1. Data analysis
tone´e canbeseenhere. Onelegofthestereophotogrammetry
framecanbeseenintheupperleftquadrantoftheimage,of Two-dimensional roughness power spectra were esti-
which 1.3 m is visible.
matedfrommeasuredheightfieldsusingThomson’smul-
titaperapproach(Thomson,1982; PercivalandWalden,
1993). Thisapproachusesthediscreteprolatespheroidal
formed on subaerial outcrops, which are subject to ad-
sequences (DPSS) as orthogonal window functions. One
ditional erosion through chemical weathering (Nichol-
of the advantages of the DPSS windows is that for a
son, 2009) over the approximately 12,000 years since the
given value of the equivalent noise bandwidth, N ,
glaciersretreatedfromthisarea(Mangerudet al.,2011). BW
several orthogonal windows are available. Power spec-
Chemical erosion creates small pits in the rock surface
tracomputedwithorthogonalwindowswereincoherently
as grains are dissolved due to oxidation and exposure to
averaged to produce a power spectrum estimate with re-
an acidic environment. These features are the primary
duced variance compared to single realizations (Percival
source of difference between roughness characteristics of
and Walden, 1993). Two-dimensional window functions
submerged and subaerial roches moutone´es. It is argued
were obtained by taking the inner product of two one-
below that the chemical weathering pits do not affect
dimensionalwindowfunctions. Tofurthermitigatespec-
thescatteredacousticfieldatthesefrequencies,andthat
tral leakage, a least-squares plane was subtracted from
subaerial roughness measurements are acceptable as in-
the height field before windowing and spectrum estima-
put parameters to scattering models used to predict the
tion. The NBW parameter was set at six, and seven
cross section from submerged outcrops.
windows in each direction were used, for a total of 49
Thestereophotogrammetrysystemusedinthesemea-
orthogonal windows.
surements consisted of two Nikon D7000 digital single
lens reflex cameras with Nikkor 28 mm f/2.8 D lenses
mounted to an aluminum frame. Baseline separations
were set at approximately 0.5 m and 1 m, which enabled 2. Roughness results
the system to operate at heights of both 1 m and 2 m
from rock surfaces, which will be called the high- and Two roughness measurements of a glacially abraded
low-resolution modes respectively. The nominal system rough surfaces are presented in Fig. 3. These measure-
resolutions determined by the camera resolution, cam- mentsweremadeofthesamerocksurface,butatthetwo
era separation, focal length, and mean distance from the different heights of the roughness measurement system,
rock surface are (∆x,∆y,∆z) = (160,158,77.4) µm for 2 m for Fig. 3(a), and 1 m for Fig. 3(b). A mean plane
the high-resolution mode and (356,357,176) µm for the was subtracted from each of the measurements before
low-resolution mode, where ∆x is the resolution paral- plotting. Because of this operation it is difficult to use a
lel to the baseline direction in the horizontal plane, ∆y globalcoordinatesystemforallroughnessmeasurements,
5
andcoordinatesarereferencedtotheirmeanvalues. The
region of overlap between the two measurements is indi-
cated by the dashed box in Fig. 3(a).
Roughness is displayed in a color representation in
which the color bar communicates surface height, and
surface slope information modifies the black/white bal-
ance (lightness), as if the surface were illuminated by a
lightsource. Thelightnessmodificationtothecolorbaris
biasedtolightervalues,sinceblackispartofthecolorbar
and white is not. Thus pixels may appear to have colors
that are not part of the color bar (e.g. white glints) due
toalargemodificationofitslightnessvalue. Thisvisual-
ization scheme is used because the dynamic range of the
color scale is dominated by the large amplitude features
and cannot resolve the low-amplitude, short wavelength
features captured by the photogrammetry system.
After a mean plane was subtracted, the root mean
square (rms) height of Fig. 3a is 13.1 mm and the rms
slope is 0.32. In Fig. 3b, the rms height is 2.01 mm and
the rms slope is 0.37. In both images, glacier flow is ap-
proximatelyinthenegativeydirection. Longwavelength
undulations in both the along- and across-flow direction
are evident at scales of approximately 50 cm. At wave-
lengthsontheorderofafewcm,theroughnessisprimar-
ily caused by scratches parallel to the glacier flow direc-
tion. These striae are likely caused by individual clasts
being dragged across the rock surface by the glacier. At
wavelengthslessthan5mmthereexistpitsinthesurface
that are likely the result of post-glacial chemical weath-
ering(Nicholson,2009)thatdonotseemtobeshapedin
any preferred direction.
Post-glacialweatheringisthelargestsourceofdiscrep-
ancy between submerged and subaerial bedrock rough-
ness characteristics. From Fig. 3(b), weathering pits are
common, although their prevalence is exaggerated some-
whatbytheshadingscheme. Therearealsoseveralsmall
cracks running through both height fields. The cracks
are likely pre-existing joints in the rock material that FIG.3. (coloronline)Roughinterfaceresultsfromaglacially
were widened by freeze-thaw cycles (Nicholson, 2009). abraded surface in (a) the low-resolution mode, and (b) the
high-resolution mode. The glaciers flowed in the negative y
Althoughpits10mmwideand1mmdeepexist,theyare
direction. The color bar corresponds to height reference to
relatively rare, with most pits less than 2-3 mm in hori-
the surface mean, and the brightness, or black/white infor-
zontalextent,andshallowerthan0.5mm. Itislikelythat
mation communicates the surface slope. The dashed box in
roughness parameters estimated for wavelengths larger
(a) indicates the portion of the surface that is shown in (b).
than 2-3 mm are applicable to underwater (submereged)
abraded surfaces. Since first-order perturbation theory
(Thorsos and Jackson, 1989) states that roughness at appear to be isotropic and are mostly restricted to di-
scales less than λ/2 (7.5 mm for the HISAS 1030 sonar) ameters of less than 5 mm, the high-spatial frequency
cannotaffectthescatteredfield,andperturbationtheory isotropic regime would seem to represent the effect of
is expected to be accurate for these surfaces, chemical chemical weathering. Since the weathering pits are con-
weathering pits likely do not affect the scattered field. fined to a region of spatial frequencies above the highest
Power spectra estimated for the height fields depicted Bragg spatial frequency accessible by the HISAS sonar
in Fig. 3 are presented in Fig. 4. The 2D power spectra (135 m−1), they likely do not contribute significantly to
areplottedasafunctionofhorizontalspatialfrequencies, the scattered field. So long as attention is restricted to
u and u . Both spectra are anisotropic at low spatial spatial frequencies less than 200 m−1, the spectral char-
x y
frequencies, but are isotropic at high spatial frequencies acteristicsofsubaerialroughnesscanbeusedtorepresent
with the division at approximately 200 m−1. This di- spectral characteristics of submerged roches moutone´es.
vision between large and small scales corresponds to a The low-wavenumber anisotropy takes the form of
5 mm wavelength. Since the chemical weathering pits broad peaks in the spatial frequency domain centered
6
FIG.5. (coloronline)Roughnessmeasurementresultsforthe
pluckedroughinterface. Thecolorscalecommunicatessurface
height,andthesurfaceslopehasbeenincludedaschangesto
grayscale value to accentuate low-amplitude roughness not
resolvedbythecolorscale. Onlythelowresolutionmeasure-
mentispresentedbecausethehighresolutionmodedoesnot
include enough facets to obtain a proper sample size.
superimposed at small scales. After a mean plane was
subtracted, the rms height of this surface was 45.6 mm,
andthermsslopeis5.2. Notethatsomeofthesteepfaces
appear to be quite smooth. This artifact results from a
FIG.4. (coloronline)Decibelversionofthetwodimensional
shortcoming of wide-baseline photogrammetry in which
power spectra in m4 of abraded surfaces shown in Fig. 3
thestereocorrespondencealgorithmfailsforsteepslopes,
as a function of horizontal spatial frequencies, u and u .
x y
and the missing areas are interpolated. The small-scale
Thepowerspectrumoftheabradedsurfacemeasuredinlow-
resolution mode is presented in (a), and the spectrum mea- roughness appears to be isotropic, and lacks the parallel
suredbythehigh-resolutionmodeisin(b). Anglesmentioned striae exhibited by glacially abraded surfaces. It is likely
inthetextaremeasuredcounterclockwisefromthe+u axis. that the small-scale roughness reflects the shape of the
x
preexistinginternaljointsurfacebeforeglacialquarrying,
and is not the result of glacial abrasion (Iverson, 2012).
at the origin, and at aspect angles, u = tan−1(u /u ), The two-dimensional power spectrum of this surface
φ y x
of 0◦, -21◦, and 77◦ for Fig 4(a), and 0◦ and 77◦ for the is shown in Fig. 6. The plucked spectrum is anisotropic
Fig.4(b). Bytheprojection-slicetheorem(Fergusonand overmuchofthespatialfrequencydomain. Atlowspatial
Wyber, 2005), angles in the spatial-frequency domain frequencies, it has peaks at u of 90◦, and ±23◦. These
φ
correspond to angles in the spatial domain. Since these directional peaks have most of their energy at the ori-
peaksarecenteredattheorigin,theirwidthssetthecor- gin and extend to spatial frequencies of 150 m−1. These
relation scale of the surface in that particular direction. peaks likely result from the large-scale facet structure of
Thewidthoftheanisotropicpeakat0◦likelycorresponds the plucked interface. Aside from the directional peaks,
to the correlation scale of small scratches perpendicular there is a background isotropic spectrum that decays as
to the direction of glacier flow. The anisotropic feature a power law, likely representing the isotropic small-scale
at 77◦ is present in both measurements and likely corre- roughness on each facet.
sponds to large-scale undulations approximately parallel Parameters for an isotropic two-dimensional power
to glacier flow. The peak at -21◦ in Fig. 4(a) likely cor- spectrum are required for the effective acoustic system
responds to the undulations present between the upper calibrationstepdetailedbelow,andasinputsforscatter-
left and lower right corners of Fig. 3(a). ingmodels. Averagingoverazimuthistypicallyonlyper-
Aroughsurfaceformedthroughglacialpluckingispre- formedforisotropicspectra. Inthiscasetheazimuthally
sentedinFig.5withthesamevisualizationschemeasin averaged spectrum is expected to reflect the behavior of
Fig.3. Thissurfaceiscomposedofapproximatelyplanar scatteringcrosssectionmeasurementsaveragedoversev-
polygonal facets at large scales, with low-level roughness eral azimuth angles, which may be anisotropic. Based
7
TABLE III. Roughness Spectrum model parameters and un-
certainty.
Parameter Abraded Plucked
φ2 [m4−γ2] 4.847×10−8 6.083×10−3
γ 2.73 4.36
2
η [%] 4.5 6.2
Φ
FIG.6. (coloronline)Decibelversionofthetwo-dimensional
powerspectruminm4ofthepluckedinterfaceshowninFig.5
as a function of horizontal spatial frequencies, u and u .
x y
Angles mentioned in the text are measured counterclockwise
from the +u axis.
x
on perturbation theory, the roughness spectrum com-
ponents responsible for backscattered power are at the
Bragg wavenumbers, 2k cosθ, or spatial frequencies or FIG. 7. (color online) Azimuthally averaged power spectrum
w
2u cosθ, where k = 2πu = 2πf/c is the wavenum- from abraded and plucked surfaces with power-law fit. The
w w w w
powerspectrumfromboththelow-andhigh-resolutionmodes
berinwater,c isthesoundspeedinwater,andf isthe
w of the photogrammetry system are plotted for the abraded
acoustic frequency. For the angles covered by the exper-
surfaces.
imentalgeometry, thiscorrespondstospatialfrequencies
between 105 m−1 and 135 m−1. The azimuthally av-
eraged spectrum for the small- and large-scale abraded
Uncertaintyinthemodelparameterestimatesisama-
spectrum, and the plucked spectrum are shown on a log-
jor contributor to uncertainty in the resulting scattering
log scale in Fig. 7. Low wavenumbers that are biased
strength measurements reported below. Parameters are
by the apodization functions are not shown here. The
estimated using a least-squares fit to the power-law in
small-andlarge-scaleabradedspectramatchveryclosely
log-log space, which means that parameter uncertainty
in their overlapping spatial-frequency domains, which is
is a function of the residual sum of squares, and number
expected. The plucked spectrum exceeds the abraded
of points used in the estimate. An additional source of
spectra by more than three orders of magnitude in the
uncertaintyisthevariationofΦasfunctionofazimuthal
lowspatialfrequencydomain,butismuchcloserinpower
spatial frequency, u . The total relative variance η2 in
at high spatial frequencies. φ Φ
thespectrumestimateisestimatedbyη2 =η2 +η2 ,
Apower-lawmodelisfittoazimuthally-averagedspec- Φ LS aniso
whereη2 istherelativevarianceduetotheleastsquares
tra, of the form LS
fit, and
Φ(u )=φ /uγ2 (1)
r 2 r
where Φ(ur) is the 2D power spectrum, ur is the radial η2 = 1 u(cid:90)max(cid:104)Φ2(uφ,ur)−Φ2(ur)(cid:105)uφ du
spatial frequency, φ is the spectral strength, and γ is aniso u −u Φ2(u ) r
2 2 max min r
the spectral slope. The subscripts on γ2 and φ2 indicate umin
(2)
thattheseparametersarefor2Dpowerspectra,following
theconventioninJacksonandRichardson(2007). Power-
lawfitsaredisplayedinFig.7,andestimatedparameters represents the azimuthal variability in the power spec-
can be found in Table III. All spatial frequency compo- trum. It is the variance of the power spectrum over
nentsofthespectrumthatwerenotbiasedfromapodiza- u divided by the squared mean, and averaged over the
φ
tion or the photogrammetry processing were used in the Bragg spatial frequencies accessible by the measurement
fit. Althoughasinglepower-lawisfittothepluckedspec- system. A parameter for total uncertainty of the spec-
trum, it appears to have some curvature in log-log space trum over the Bragg spatial frequencies is used rather
andcouldconvergewiththeabradedspectraathighspa- than uncertainty of individual parameter estimates. Un-
tial frequencies if higher-resolution data were available. certainties can be found in Table III.
8
III. ACOUSTIC SCATTERING EXPERIMENTS array is related to the seafloor height. Since phase is
wrappedto2πradians,discontinuitiesinthephasediffer-
Acousticbackscatteringmeasurementswereperformed encemustbedetectedinthepresenceofnoise. Typically,
off the coast of Larvik and Sandefjord, Norway in April an n×n window is used to provide an estimate seafloor
2011 (Midtgaard et al., 2011). This experiment was per- height with reduced variance at the cost ofreduced reso-
formed by the Norwegian Defence Research Establish- lution (Sæbø et al., 2013; Sæbø, 2010). Bathymetry es-
ment (FFI) aboard the HU Sverdrup II. Data were col- timatescanbeusedinscatteringstrengthmeasurements
lectedusingtheHUGIN1000HUSautonomousunderwa- by providing an estimate of the local seafloor slope and
ter vehicle (AUV) equipped with the HISAS 1030 inter- global grazing angle at each pixel. The assumption of
ferometric SAS system (Fossum et al., 2008). This sonar a constant or planar seafloor is often used in scattering
has not been calibrated in terms of its open-circuit re- strength measurements (Jackson et al., 1986a), which is
ceiversensitivity,s ,orthesourcestrength,s ,although severely violated in the case of the rock outcrops stud-
r 0
beam patterns were measured. These parameters must ied in this research. The local seafloor slope is estimated
be determined in order to estimate the absolute seafloor usingaweightedaverageversionofafinitedifferenceop-
scattering cross section. erator on the SAS bathymetry to reduce high-frequency
noise. The weights are computed using an m-point least
Theseafloorinthisareaconsistedofroches moutone´es
squaresfittoaquadraticfunction, withthederivativeof
surrounded by a sediment of cobble with a mud matrix.
thepolynomialcomputedanalytically(SavitskyandGo-
The water depth at the experimental site was approx-
lay, 1964). Second-order polynomials were used to esti-
imately 30 m, with nominal vehicle altitude of 10 m
mateallslopes,withwindowsizesof55pixelsforplucked
from the seafloor. The sound speed profile was slightly
surfaces, and 23 pixels for abraded surfaces. The larger
upward-refracting with the surface sound speed at ap-
window size was used to suppress effects of steps on lee-
proximately 1452 m/s, and approximately 1456 m/s at
wardsurfacesofrochesmoutone´es andfocusonthemean
the seafloor. The change in sound speed of the lower 10
trend.
m of the water column was a maximum of 3 m/s, and
therefore refractioneffectson the localgrazing angle can
be ignored. Further details on the experiment can be
found in Midtgaard et al. (2011). B. Estimating the scattering cross section from synthetic
aperture sonar data
A. Synthetic aperture sonar overview Estimation of scattering strength from SAS systems
is similar in principle to estimation using other real-
Syntheticaperturesonar(SAS)isadatacollectionand aperture sonar systems designed to measure scattering
beamforming technique in which a transmitter and a re- strength (Jackson et al., 1986a; Williams et al., 2002),
ceiver array move along a track, which is typically linear with the exceptions of nearfield beamforming and the
or circular. Acoustic energy is transmitted at regular in- abilitytoestimateseafloorslope. Typically,beamformed
tervals, and the resulting scattered field is sampled, also time series are incoherently averaged over independent
at regular intervals. The received field from many trans- areas of the seafloor, and then the sonar equation is in-
missions can be concatenated to form a synthetic array verted for scattering strength for each intensity sample.
thathasalengthmanytimesthatofthephysicalreceiver SinceSASisanear-fieldimagingtechnique,transmission
array. A SAS image can be formed using a variety of losschangessignificantlyalongthesyntheticarrayforall
beamformingtechniques,thesimplestandmostrobustof ranges. Consequently, sonar equation terms that vary
which is the backprojection, or delay and sum algorithm as a function of sensor location, i.e. spherical spreading
(Hawkins, 1996). Synthetic array lengths scale linearly and attenuation, are removed before the image is beam-
with pixel range and images have a constant Cartesian formed. Terms that do not vary as a function of sensor
resolution over the whole scene. The array lengths can position, such as the ensonified area and calibration pa-
be quite long and in most situations imaged objects are rameters,areremovedafterimageformation. Theresult-
in the near field. Aligning successive physical aperture ing calibrated SAS pixel intensity values correspond to
locations to form a synthetic array is challenging in the theunaveragedcrosssection,σ˜,forazimuthallyisotropic
ocean environment due to hydrodynamic forces on the scatterers.
sonar platform. Typically, an inertial navigation system A sonar equation is used that relates pixel intensity
combined with displaced-phase center acoustic naviga- to the scattering cross section. It is equivalent to Eq.
tion is used to estimate the sonar element locations to (G.11) in Jackson and Richardson (2007) but adapted
within a fraction of a wavelength (Bellettini and Pinto, to SAS quantities. This equation is valid so long as the
2002). MoreinformationonSASprocessingcanbefound scattering cross section and vertical beam pattern are
in Hansen (2011) and references therein. both slowly varying within the system resolution. Let
Interferometric SAS is possible with two vertically- v be the complex matched-filtered voltage sensed by
ij
separatedreceiverarraysonasonarplatform. Thephase the jth receiver element, and delayed such that it corre-
difference between SAS image pixels formed from each spondstotheithpixel. Letxbethealong-trackposition
9
of the pixel location, and y be the cross-track (ground that since SAS is a nearfield imaging algorithm, a given
range) position. The complex value of the ith pixel, q , patch of the seafloor is always in the nearfield of the
i
is the output of the delay-and sum beamformer, and is a synthetic array. The Fourier transform relationship be-
weighted sum over the synthetic array, with the weights tween the beam pattern in the far field and the aperture
determined by transducer patterns, and by tapering ap- weightingfunctionisusedtocomputeδx,sincethebeam
plied to reduce side-lobes. When corrected for propaga- pattern in the focal plane of a focused near-field array is
tion, vertical beam pattern effects, and coherent gain, q identical to the far-field beampattern (Mast, 2007). The
i
is defined by range resolution, δy is determined by the transmitted
bandwidth and spectral weighting of the received pulse,
 −1
Nr and is equal to 3.25 cm.
(cid:88)
qi = wjbtx(φij)brx(φij) × Therelativesignalgain,Γ,isdefinedhereasthepower
j ratio between the partially coherent array gain and the
(3)
  coherent array gain (Cox, 1973; Carey and Moseley,
(cid:88)Nr wjbtx(φij)brx(φij)|btx(rθi2ijje)2bαrrxi(jθij)|vij 1w9h9e1n;thCearreecye,iv1e9d98s)ig.nIatlicshaarflauccteturiazteisngthqeuabniatsityobinsesrtveaedd
j
of a point scatterer. In this work, partial coherence is
where w is the processing weight applied to the jth re- due to two mechanisms, 1) phase fluctuations due to un-
j
ceiver, α is the attenuation of seawater, r is the dis- certainty in the sensor positions and oceanographic con-
ij
tance (slant range) from pixel i to receiver j, b (φ ) ditions, and2)amplitudeandphasefluctuationsduethe
tx ij
andb (φ )arethetransmitandreceivehorizontalbeam scattering characteristics of random rough surfaces. The
rx ij
patterns, and b (θ ) and b (θ ) are the vertical beam firstmechanismisincludedasasourceofuncertaintyand
tx ij rx ij
patterns. The variables φ , and θ , are the horizon- discussed at the end of the next subsection.
ij ij
tal and depression angles from the sensor j and pixel The spatial coherence of rough surface scattering is a
i in the sonar’s coordinate system. The HISAS 1030 complex topic and a rigorous treatment is outside the
sonarhasasingletransmitterandmultiplereceivers. The scope of this work. However, an intuitive argument is
phase-center approximation (Bellettini and Pinto, 2002) giventhatΓisthesameforallpixels. Ingeneral,theco-
is employed to work with an equivalent monostatic con- herenceofthefieldduetoscatteringfromaroughsurface
figuration with colocated transmitters and receiver. The has a Fourier transform relationship with the covariance
unaveraged scattering cross section, σ˜, is related to qi of the pressure at the scattering surface (Mallart and
through Fink, 1991). Theoretical models of coherence typically
employ an isotropic point scatterer model (see Jackson
cos(θ−θ )
|q |2 =(s s )2Γσ˜δxδy 0 (4) and Moravan (1984) and references therein), or the van
i r 0 cosθ
Cittert-Zernike theorem (vCZT) (Born and Wolf, 1999;
where s is the receiver voltage sensitivity in V/Pa, s is Mallart and Fink, 1991). Under the experimental condi-
r 0
thesourcelevelinPa-m,andthecosinetermsresultfrom tionsinthiswork,bothmodelspredictthatthecoherence
converting from the transducer array coordinate system alongthereceiverapertureisequaltotheautocorrelation
to the seafloor coordinate system, with θ the depres- of the transmitting aperture function. This relationship,
0
sion angle of the main response axis of the transmitter whichwewillcallthevCZT,dependsonthetwoequiva-
as discussed in Appendix G of Jackson and Richardson lent assumptions that the covariance of the surface field
(2007). Thequantitiesδxandδy arethealongtrackand behaveslikeaDiracdeltafunction,orthatthescattering
range resolutions of the system respectively. These reso- cross section is isotropic. For the ground ranges studied
lutions are defined as the width of rectangular windows inthiswork,thegrazinganglesalongthesyntheticarray
that have the same integrated power as the ambiguity vary by less than 0.1◦. Therefore we can restrict at-
function along a particular direction, similar to the con- tention to azimuthal coherence, and only the aziumthal
cept of equivalent noise bandwidth. This definition of dependence of the scattering cross section is of import.
resolution differs from others, such as the 3 dB down Since azimuthally isotropic scatterers have already been
points, and is a simple way to define a sonar equation. assumed above, it is an appropriate assumption to use
The local seafloor slope at each pixel is used to define for coherence as well. Since the scattering cross section
the local grazing angle. Ensembles are formed by group- is not anisotropic for all situations, this result may not
ing local grazing angles into 1◦ bins and estimating the always be applicable. However, the intuitive argument
sample mean. canbemadethatifthebistaticscatteringcrosssectionis
The along-track resolution, δx, is nominally 3.14 cm isotropicovertheazimuthanglessubtendedbythetrans-
fortheparametersoftheHISAS1030andprocessingpa- mitter and receiver arrays, then the spatial coherence is
rameters used for these data. Due to vehicle motion, the practicallyequaltotheautocorrelationofthetransmitter
valueofδxfluctuatesbyasmuchas5%basedoncalcula- aperture.
tionsofthealongtrackresolutionforeachpixel. Tosim- UsingthevCZT,Γcanbedemonstratedtobethesame
plifydataprocessing,thenominalvalueofδxisused,and forallpixels. Thetransmitandreceivearraysarealways
its variability is included in uncertainty analysis. Note matchedinlengthduetothephasecenterapproximation
10
(Bellettini and Pinto, 2002). Additionally, the transmit face. Apartfromdiscretefeaturessuchasdropstonesand
and receive aperture weights have functional forms that edges, the image amplitude appears to slowly decrease
scale as a function of x/y = tan(φ), where φ is the hor- as the absolute value of the ground range from the sonar
izontal angle from the sensor location to a given pixel. track increases. Since all propagation and system effects
Therefore range dependence in both the partial coherent have been removed, this slow change in image amplitude
gain and coherent gain is canceled out by the definition can be attributed to the decrease in scattering strength
of Γ. Under the assumption of azimuthal isotropy, Γ is with grazing angle, which is confirmed in the scattering
identicalforallpixelsandmaybenormalizedoutduring strength plots below. The surrounding sediment is com-
the calibration procedure described below. Note that if posed of approximately 4 cm cobble in a mud matrix,
the system had been calibrated using a technique such as indicated by diver samples. The surface of the cali-
asreciprocity,thevalueforΓwouldneedtobeexplicitly bration rock scatters sound at approximately 7 dB lower
calculated based on the system parameters. than surrounding cobble, which is a further indication
that the surface is quite smooth. The calibration rock
surface, like most of the exposed bedrock in this area,
has been formed through glacial abrasion, as indicated
C. System calibration
by the smooth image texture and low amplitude of the
scattered field.
Thetwoparameterscharacterizingsystemcalibration,
The relative scattering strength is averaged over all
s and s , have not been measured, but are required for
r 0
of these aspects for the system calibration. In IIB, the
estimating the scattering cross section. An effective cali-
bration technique is applied that estimates, s=(s s )2, roughnesspowerspectraofglaciallyabradedsurfaceswas
r 0
seentobeanisotropicatscalesontheorderoftheBragg
byfittingthescatteringstrengthofanareaoftheseafloor
wavenumber for this sonar system. To minimize bias
toavalidscatteringmodelwithknowninputs. Notethat
due to anisotropy, the calibration procedure compares
this technique also normalizes out any bias that is con-
azimuthally-averaged scattering strength to the small-
stant for all pixels and measurement locations, such as
slopemodelcomputedwithazimuthally-averagedrough-
therelativesignalgain,Γ,discussedabove. Thismethod
ness power spectrum parameters. Azimuthal variability
ofcalibrationhasbeenappliedpreviouslyusingaseafloor
of the roughness spectrum at the range of Bragg spa-
scatteringmodelbyDwan(1991)fortheGloriaside-scan
tial frequencies was used to compute error bars on the
sonar, andusedforscatteringstrengthmeasurementsby
scattering model, and thus the estimate of s.
Lyonset al.(1994). Asimilartechniqueusingman-made
targets,suchasmetalspheres,hasbecomestandardtech- The first-order small-slope approximation (SSA)
nique (Foote et al., 2005) for sonar calibration. (Voronovich, 1985) was used to compute a model cross
Arockoutcropwithverylow-amplituderoughnesswas sectioncurvefromwhichsisestimated. Initsfirst-order
chosenbasedoninspectionofSASimagerytoprovidean form, the SSA can be decomposed into the product of a
effectivesystemcalibration. Theapparentlow-amplitude term depending on the geoacoustic properties, and the
roughness enabled the use of a first-order scattering Kirchhoffintegral(JacksonandRichardson,2007, Chap.
model. SAS images of the calibration rock taken at var- 13). The Kirchhoff integral was numerically evaluated
ious azimuthal angles are presented in Fig. 8. This out- using the algorithm developed by Gragg et al. (2001)
crop extends approximately 0.5 m above the surround- and Drumheller and Gragg (2001). The SSA is accurate
ing sediment. The color scale corresponds to the decibel in the respective regions of validity of the simpler small-
value of the relative scattering cross section. Note that roughnessperturbationapproximationandKirchhoffap-
the decibel quantity here and in Fig. 10 is dimensionless proximation (Jackson and Richardson, 2007, Chap. 13).
and thus appears without a reference unit. The vertical Although these bounds have not been established for
axescorrespondtothedistancealongthesyntheticaper- power-law surfaces, perturbation theory for fluid inter-
ture,andhorizontalaxescorrespondtothegroundrange faces and a von K´arm´an spectrum has been shown to be
to the sonar, with negative values representing distance accurate with kwh ≈ 1, where h is the rms roughness
fromtheportsideofthevehicle,andpositivevaluesfrom (Thorsos et al., 2000). The rms roughness of the surface
the starboard side. Boxes indicate areas that were used presented in Fig. 3(b) is 2.09 mm, which corresponds to
toestimatescatteringstrength. Inthecasethatmultiple kh=0.86,andthermsslopeis0.37. Giventheseconsid-
boxesintersect,theirunionisusedtodefineamorecom- erations, the SSA is expected to be accurate for glacially
plicated boundary that avoids double counting of pixels. abraded surfaces.
Ensonificationdirectionisnominallyfromthelinewhere The elastic SSA model requires 9 inputs, φ , γ , c ,
2 2 p
the ground range is zero, and can be from the left or c , c , δ , δ , ρ , and ρ , the density of water. It was
t w p t b w
right of the images depending on the sign of the ground argued in IIB that roughness characteristics of the sur-
range values. The calibration rock outcrop is approxi- faces from which the acoustic calibration measurements
mately10mx5minhorizontalextent, andappearedto were made are similar to the roughness characteristics of
be flat and smooth based on the image texture, as sup- glaciallyabradedsurfacesattheroughnessmeasurement
ported by interferometric SAS bathymetry. Dropstones, site. Therefore, the power spectrum estimates presented
likely deposited by glaciers, are present on the rock sur- in IIB are assumed to represent the roughness power