CMel An
a .W. wml 0 “A
) Ay Ai ‘4
Mississippi 39522-5001 April 1985 BAe
So WHO! i TR 279
DOCUMENT
COLLECTION
Description, Analysis, and Prediction of
Sea-Floor Roughness Using Spectral Models
CHRISTOPHER G. FOX
Advanced Technology Staff
Approved for public release; distribution unlimited.
| EC Prepared under the authority of
/ Commander,
YS Naval Oceanography Command
FOREWORD
Knowledge of sea-floor characteristics is required to interpret
deep-sea processes correctly. This technical report discusses the
use of spectral models in describing, analyzing and predicting sea-floor
roughness.
ip Commanding Officer
om
0 0303- 00652075
WES
WHOL
sss ssn ss sss
SECURITY CLASSIFICATION OF THIS PAGE (When Data Entered)
REPORT DOCUMENTATION PAGE pa a RS a ot ala
T. REPORT NUMBER 2. GOVT ACCESSION NO, 3. RECIPIENT'S CATALOG NUMBER
TR-279 TR 279
4. TITLE (and Subtitle) S. TYPE OF REPORT & PERIOO COVERED
Description, Analysis, and Prediction of
Sea Floor Roughness Using Spectral Models
6. PERFORMING ORG. REPORT NUMBER
- AUTHOR(ae) 8. CONTRACT OR GRANT NUMBER(2)
Christopher G. Fox
. PERFORMING ORGANIZATION NAME AND ADORESS 10. PROGRAM ELEMENT, PROJECT, TASK
AREA & WORK UNIT NUMBERS
Advanced Technology Staff
Naval Oceanographic Office
NSTL, MS 39522-5001
+ CONTROLLING OFFICE NAME ANDO ADORESS 12. REPORT DATE
Naval ceanographic Office
th ee
4. MONITORING AGENCY NAME & AOORESS(If different from Controlling Office) 1S. SECURITY CLASS. (of thie report)
Unclassified
1Sa. DECLASSIFICATION/ OOWNGRADING
SCHEDULE
16. DISTRIBUTION STATEMENT (of this Report)
Approved for public release; distribution unlimited
17. OISTRIBUTION STATEMENT (of the abetract entered In Block 20, If different from Report)
18. SUPPLEMENTARY NOTES
19. KEY WORDS (Continue on reverse aide if neceeeary and identify by block number)
Roughness, frequency spectrum, marine geology and geophysics, acoustic
bottom interaction, multi-beam sonar, bathymetry, numerical modelling.
20. ABSTRACT (Continue on reverse side If neceseary and identity by block number)
A method has been developed which allows a valid statistical model of
the variability of oceanic depths to be derived from existing digital
bathymetric soundings. The bathymetry of the world ocean has been
mapped using a variety of acoustic sounding instruments and traditional
| contouring methods. The bathymetric contours represent a low-frequency,
_ deterministic model of the seafloor. To describe the higher frequency
_ (continue on reverse)
DD , "on™, 1473 Gam) oe 1 NOV 68 1S raceverc
SANROMOZ CEOS 3660) SECURITY CLASSIFICATION OF THIS PAGE (When Data Entered)
ee D
SECURITY CLASSIFICATION OF THIS PAGE (When Deta Entered)
variability, or roughness, of the seafloor requires the development of
an appropriate statistical method for generating a valid stochastic
model. The smooth contoured surface (often preserved as a geographic
grid of depths), when supplemented by such a roughness model, provides a
complete description of the relief. The statistical model of seafloor
roughness is also a valuable tool for predicting acoustic scattering and
bottom loss, and in addition contains a wealth of geological information
for interpreting deep-sea processes.
To allow the variability of depths to be described as a function of
scale (spatial frequency), the amplitude spectrum is employed as the
fundamental statistic underlying the model. Since the validity of the
amplitude spectrum depends upon the assumption of a statistically sta-
tionary sample space, a computer algorithm operating in the spatial
domain was developed which delineates geographic provinces of limited
statistical heterogeneity. Within these provinces, the spectral model
is derived by fitting the amplitude estimates within the province with
one or several two-parameter power law functions, using standard regression
techniques.
The distribution of the model parameters is examined for a test area
adjacent to the coast of Oregon (42°N-45°N, 130°W-124°W), which includes
a variety of geologic environments. The distribution of roughness
corresponds generally with the various physiographic provinces observed
in the region. Within some provinces, additional complexities are
apparent in the roughness model which cannot be inferred by simply
studying the bathymetry. These patterns are related to a variety of
geological processes operating in the region, such as the convergence of
the continental margin and the presence of a propagating rift on the
northern Gorda Rise.
In many cases, the roughness statistics are not constant when observed
in different directions, due to the anisotropic nature of the seafloor
relief. A simple model is developed which describes the roughness
statistics as a function of azimuth. The parameters of this model
quantify the anisotropy of the seafloor, allowing insight into the
directionality of the corresponding relief-forming processes. Finally,
the model is used to successfully predict the roughness of a surface at
scales smaller than those resolvable by surface sonar systems. The
model regression line (derived from a hull-mounted sonar) is compared
to data from deep-towed sonars and bottom photographs. The amplitude
of roughness is predicted to within half an order of magnitude over
five decades of spatial frequency.
S/M 0102- LF- 014-6601
SECURITY CLASSIFICATION OF THIS PAGE(When Date Entered)
DESCRIPTION, ANALYSIS AND PREDICTIONS OF SEA-FLOOR ROUGHNESS
USING SPECTRAL MODELS
Christopher Gene Fox
Submitted in partial fulfillment of the
requirements for the degree of
Doctor of Philosophy
in the Graduate School of Arts and Sciences
COLUMBIA UNIVERSITY
1985
ABSTRACT
Description, Analysis and Prediction of Sea-—Floor Roughness
Using Spectral Models
CHRISTOPHER GENE FOX
A method has been developed which allows a valid statistical model
of the variability of oceanic depths to be derived from digital bathy-
metric soundings. Existing bathymetric contour charts represent low-
frequency, deterministic models of the sea floor. To describe the
higher frequency variability, or roughness, of the sea floor requires
the development of a valid stochastic model. The statistical model of
sea-floor roughness is also a valuable tool for predicting acoustic
scattering and in addition contains a wealth of geological information
for interpreting deep-sea processes.
To allow the variability of depths to be described as a function of
scale (spatial frequency), the amplitude spectrum is employed as the
fundamental statistic underlying the model. Since the validity of the
amplitude spectrum depends upon the assumption of a statistically sta-
tionary sample space, a computer algorithm operating in the spatial
domain was developed which delineates geographic provinces of limited
statistical heterogeneity. Within these provinces, the spectral model
is derived by fitting the amplitude estimates with one or two power law
functions.
The distribution of the model parameters is examined for a test area
adjacent to the coast of Oregon. The distribution of roughness corre-
sponds generally with the various physiographic provinces observed in
the region. Additional complexities are apparent in the roughness model
which can not be inferred by simply studying the bathymetry. These pat-
terns are related to geological processes operating in the region.
In many cases, the roughness statistics are not constant when
observed in different directions, due to the anisotropic nature of the
sea-floor relief. A simple model is developed which describes the
roughness statistics as a function of azimuth. The parameters of this
model quantify the anisotropy of the sea floor, allowing insight into
the directionality of the corresponding relief-forming processes.
Finally, the model is used to successfully predict the roughness of a
surface at scales smaller than those resolvable by surface sonar
systems. The model regression line (derived from a hull-mounted sonar)
is compared to data from deep-towed sonars and bottom photographs. The
amplitude of roughness is predicted to within half an order of magnitude
over five decades of spatial frequency.
ye fal wins rm ial
imo capes
| sea ih el ges eld
‘ Ratan ‘ut hserlee te We
re que its iy’ tag th yay ak
ssi ayes een ipa ris et e007
adr So an Wgendenhin la on ai iseaebanacith jen thie, i versal
! sa a eal eds" amy) bere
ee We malware btsbiud' “i arate weap |
chipt asuvedFirnsaess “Bie wh a Raniovier he Lobe
Frinalis pntwin Esti’ “ia otitis, aes 95" cities ane
PE satis ae Gam iy an ihataaher as si ay iahew old ea
thet | ettie® wa: nape oe ie ——, erty a basa:
Rate rind Sih o ait Heal: anit pana wonea eee % vad pag Ad
ef prgeay ott Atta “ann la adh ar bone ak a7
nie ee pak 4 Salad bene tonre ui arene ‘, shoes bi -
Ups Ste, 2 aio Dn rags etal eth oo
- eo Wcull hes a. o mee tt 08 ye ee i ad ale ”
ee 7 q re " a“ i ‘
] sa ails iw ioe
2° - = ay ie
Table of Contents
IACKNOWLEDGEMEN TS icretslorere:cueleieveiel sioleroleletercioveloleteteloleloieiceiolctetcterelevererereveretlxX
LIST OF SYMBOL Sicrereporevetelel elelichelevotonslelcvelelelevelolcloreloicveleneiorerelchenelelereielenerererXs
LIST OF TEL USTRATIONS ivevciere cc lercvevcrote clever everevevalelelevevexereveie ave evovsie ew:
Wo ISpy COHN CN Io go g 60000000 0000000 DDD00D000000000000000000000!
Oo NADH NNO Bo ogo 06600000000 d OSC 000 00000000000 00000000000"%
JH EREVIOUSMWORKGielololeterclevalelelolotorclelelelalolerereverelclelotetercrenctatevovelavercycvcien?/
/, STATISTICAL CONSIDERATIONS.....ccccccccccccccvccccccsccell
A. Statistical Measures of Roughness....ccccccccccccccoll
B. Validity of Measurement Over Large AreaS....ccccee+.1l6
C. Functional Representation of Spectra....ccccccccccceld
D. Extension of Model into High Frequencies.....cccccce28
Ew Effect of Linear Featuresc.ccecccceecccecccccecececea]
F. SUMMA TY oie o1eie.clelolelorelese cele) sseierulerevelelcroreree) everareravevaretnlevaroveicvers he:
vi
Si
DELINEATION OF STATISTICALLY HOMOGENEOUS PROVINCES......45
General SPhilOSOphy,cjicrercielciclelelelelelevelerelclelelotelelalolelelelereleverevevee>
Generation of Amplitude Spectra....ccccccccccccccccce48
Physical Interpretation of the Amplitude Spectrum...50
TheePhase| (Specie rumss.c)c/cclelelclelelelale/elcleleleleiololeleleleleloleleleicteter).O
Creation of the Model and Interpretation............63
ANISOTROPY OF SURFACES <\<jcicie/ovcleielelelel elec creleieloieteleleleveteretoisiersieienio
Theoretical Mode ilicicjcverelevelelelolelelejlelelcleleycvoloverctolevereretevonevoreleten io
Comparison of Theoretical Functional Forms with
Multibeam Sonar Datiacs os ccc cc cclelcicicicicicicicicicicicie celslelelele G7,
Estimation of Two-Dimensional Spectra from
Randomly Oriented Ship Track.....cccccccccccccccesselO3
PREDICTION OF HIGH FREQUENCY ROUGHNESS......ccccccccce lJ
Source of Error in Spectral Estimates.....ccccceeeelO7
Propagation of Error to High Frequency Estimates...11l
Comparison of Surface-Ship Sonar Results to
Deep-Towed Sonar Results and Results from
Bottom PHOEOLTAPHY,ccierelelolcielclelc elelelelelelclelcvcreheiclerolclerelerereier lite
vii
8. SUMMARY AND CONCLUSIONS... ..cccccccccccccccccccccccccoll/]
9. REEERENGE Siercloielererctolereleloleleielelereleleloieiclelcielelelelelelelctelelereleieteleveverenli2O
APPENDIX
APPENDIX
APPENDIX
APPENDIX
APPENDIX
Algorithms for Performing Power Law
Regression AnalysiS..cccccccccccvccccccccecl 24
FORTRAN Software for Weighted log-log
BL Biss fae olieieyeveveue ciaveteioleue sere) syevexerexoieveyeretavene euaseterevererol
FORTRAN Software for Iterative linear-
linear FeLi a vex avvaieusvere areleueteve ieiaucheverateloreieneveleneiererati ele:
Algorithm and Performance Tests for
Amplitude Spectrum Province Picker.........135
FORTRAN Software for Province Picker.......149
Linear Combination of Sinusoids of
Identical Wavelength....ccccccccccccccccccel/8&
FORTRAN Software for Iterative Sinusoid
FTG ara tarav'elleiri'ollovs latleveqeverevovekere ove talopeveloyetoiotote stereietolerever lO
FORTRAN Software for Computing Pre-
Whitened SPCCELacielejelciclelelcleleloicicleloloielcloistelotereterer Oo,
Azimuthally Dependent Variations of Spectra
Derived from Artificially Generated
Anisotropic Surfaces.....ccccccscccccccccee 202
DESCRIPTION OF SLE RMSareteleseicvolcielerevcicicle re) ciclorersleteleleiciciereretetoleielercieiere Le
viii
Acknowledgments
This study was conducted while the author was simultaneously a grad-
uate student at Lamont-Doherty Geological Observatory (L-DGO) and an
employee of the U.S. Naval Oceanographic Office (NAVOCEANO). Original
funding was provided through the long-term training program of NAVOCEANO
and I thank all who are involved with the program for providing such a
fine opportunity. Additional funding for the effort at L-DGO was pro-
vided through the Office of Naval Research, and special thanks are due
to Murray MacDonald and Gerald Morris. The remainder of the project was
performed with the support of NAVOCEANO and, in particular, Tom Tavis
and my colleagues on the Advanced Technology Staff.
The original concept of the study, many of the approaches used, and
the recognition of the value of the study to underwater acoustics are
all from the fertile mind of Tom Davis, NAVOCEANO. Special thanks to my
major professor, Dennis Hayes, L-DGO, for recognizing and pursuing the
geological aspects of the study, and guiding me through the graduate
school maze under such unusual circumstances. Also thank you to the
staff of L-DGO (especially Mia Leo and Mary Ann Garland), my friends and
fellow graduate students at L-DGO, and my friends and colleagues at
NAVOCEANO, in particular Terry Blanchard for computer assistance.
Data were provided from several sources and many thanks are due to
all of the following:
John Farre Columbia University SEAMARC=1
Jeff Fox Univ. of Rhode Island SEABEAM
Geology Branch NAVOCEANO Ba thyme try
Geomagnetics Division NAVOCEANO Magnetics
Richard Gregory NAVOCEANO SASS
Mary Linzer Scripps Inst. of Ocean. Deep-Tow
Arthur Nowell Univ. of Washington Bottom Photo
Mark Wimbush Univ. of Rhode Island Bottom Photo
In addition, to Drs. Davis and Hayes, many ideas were provided by
Jim Cloutier of NAVOCEANO, David Berman of the Naval Research Labora-
tory, and Benoit Mandelbrot of IBM Watson Research Center.
Word processing chores were most graciously performed at L-DGO by
Carol Elevitch and Erika Free.
Illustrations were laboriously constructed by Donna Waters
(NAVOCEANO) and David Johnson (L-DGO).
List of Symbols
a coefficient of spatial frequency
Agrb, initial estimates of a,b for iterative regression
a,b,x estimated values based on regression estimate
Aa,Ab correction increments for a,b
a(6),b(8) azimuthally dependent function of a,b
A amplitude
b exponent of spatial frequency (spectral slope)
Cy cross product at lag k
D Fractal dimension
E expected value
£ (8) ,£,(8) azimuthally dependent function of a,b
f (x) general function
F(s) Fourier transform of f(x)
FFT Fast Fourier Transform algorithm
k lag
n number of observations
P(|)
Var
Y(k)
number of values in ensemble average
conditional probability distribution function
power
final residuals in least squares technique
root mean square
spatial frequency
mean level of a(6)
amplitude of cosine component of a(0)
variance of random distribution
general independent variable
random variable
observation i
statistical mean of random variable X
general dependent variable
scaling factor
residual error
autocovariance as a function of lag
spatial wavelength
xi
d! apparent wavelength
v minimized residuals
p(k) autocorrelation as a function of lag
8 azimuth
8, azimuth perpendicular to trend (strike) of
surface
2
has the Fourier transform
xii
5-5
5-6
List of Dllustrations
Scale dependence of stationarity of the
MEAN c\aicoinisielelalelaleiohetelclevevovelolelcisicteleleleloielorersielelelcleveleleletctercverevererl (oO
Method of Davis (1974) for defining homo-
geneous provinces in a designated wave band ..........20
Typical amplitude spectrum of sea-floor topog-
F-47118 7500000000000000000000000000000000000000000000000043)
Spectral estimates derived by T.H. Bell
(CIEE) G6 06050000 OO0DD00DOD0OO000OO000 000000 000000000044
White-noise level of an amplitude spectrum..........-.3l
Effect of directional sampling of corru-
gated SULLAC EC cisicvercicie wis leveloielcleleicieleleieieieleleleleletercveleielelelercvereterets ss
Bathymetric chart of the Mendocino Frac-
ture ZONE iavevera a wie Tinie Ge 6 eters Ocerehornlalies ev olavereveleveleiorevelevereyolelerererevets0
North-south profile of Mendocino Fracture
Zone and amplitude spectrum....cccccccccccccccccccccce rs
East-west profile of Mendocino Fracture Zone
and amplitude spectrum...cccccccccccccccccccccccccccee3d&
Comparison of Figures 4-8 and Uo goq00000D00000000000e")
Importance of prewhitening of amplitude
BPSCELalcrcrcvclevcicie ciclelelelelolelelelolololcloiolelelorelelelelerelelelololelololoteloleleretel ot
Relationship of spectral slope parameter
to the aspect ratio of sinusoids at dif-
ferent frequenciesS....ccccccccccccccccccccccccccccccscI4
Microrelief map of sea floor collected
August 16, VOB Vie scatarei Wrovevensialle vevenetaccpove aleve eleleVe otwvelerorcreieieieictereOile
Microrelief map of sea floor collected
July 20, UDG Wig scerareveverstoversits ersietacveiovevevere erave ei elevetevereleletelevereteters ale,
Composite spectra from sea-floor microrelief.......20+59
Typical phase spectrum and statistical distri-
bution of phase angles from sea-floor relief..........61
Distribution of roughness in the vicinity
of the Gorda Rise, Northwest Pacific Ocean...........63a
xili
6-2
6-3
6-4
6-8
6-9
6-10
6-11
Schematic illustration (from Hey, 1977) of
the magnetic anomaly pattern resulting from
a PLOpagatingert He ajeiercieisloleieleleleiclorcleleiclelelelolelejcicicleleielsiclerere n OS
Magnetic anomaly chart of the Gorda Rise.....cccccceee/0
Amplitude spectrum of a long bathymetric
profile spanning non-homogeneous reliefs....ccccccccce/3
Typical amplitude spectrum of relief from
the Tufts Abyssal Plain, Northwest Pacific Ocean......74
Typical amplitude spectrum of relief from
the Continental Rise, east coast U.S...cccccccccccccccclO
Location of bathymetric profiles showing
similar amplitude spectra....ccccccccccccccccccccccccc/&
Graphic representation of an elementary
theoretical model for anisotropic surfaces......cccccc82
Radial sampling pattern used in anisot-
ropy SEUAL OSS sia eye lateratav cle lovelotelcvavelertetalotoeVetelveleVeleleeve vee evéveleiere Oo
Distribution of spectral parameters versus
azimuth for theoretical anisotropic surface.......0..e8)
Location of two-dimensional SASS bathymetry
from the Gorda RUG Cicieoveiercicuerevetoveravalevevevercvesstotederabeve evalerelewicieed.
Graphic representation of bathymetry from
the Gorda IRS Cie aiererorelsieiej eles ececcvevela¥oralctetercteretskelloveletel everereuarereceiel
Distribution of spectral parameters versus
azimuth for Gorda Rise bathymetry..ccccecscccccccceece92
Comparison of model spectrum to measured
spectrum for “worst case” profile....ccccccccccccccc co 094
Location of two-dimensional SASS bathym-
etry from the upper continental rise.....ccccscceeceee9D
Graphic representation of bathymetry from
the upper continental rise....cccccccccccccccccccccec cdl
Distribution of spectral parameters versus
azimuth for upper continental rise bathymetry.........98
Location of SEAMARC - 1 derived bathymetry
from the Carteret CanyOniclersivicieleleiclelolelcleteveloleleloicicleleiereverere lOO
xiv
B-l
B-2
B-5
B-6
Distribution of spectral parameters ver-
sus azimuth for Carteret Canyon bathymetry..........-l0l
Two-dimensional amplitude spectrum from
the Gorda Rise with model spectrum.....c.ccccccscceece oe lO4
Results of ensembling spectral estimates
from sixteen parallel profiles from a multi-
beam sonar SYSCEM. cccccccciocvccccccleciciececcccccoosssl lO
Composite amplitude spectrum of bathymetry
profiles of various scales.....ccccccccccccccccvcvccclld
Frequency response of band-pass filters
used to estimate the amplitude spectrum
discretely in the spatial domain.......cccccccccccccol3/
Uniformly distributed random noise with
corresponding amplitude spectrum... .cccccccccccccccs ol 40
Normally distributed random noise with
corresponding amplitude spectrum...cccccccccccccccce oll
Inverse FFT generated white noise with
corresponding amplitude spectrum......cccccccccccccccl42
Province picker output with white-noise
AN DU Ciererevcle! cleleVelelelelololevolelclolelololoieterclereleloverelelole cloleheleveleleleleleloronl 4
Province picker output with mixed input
of known Siionial’sicielclclelclelelelctetelevelelevelcie/eleVelelele) efolololeloletolorerere)l 40
Province output with input of digital
balthymetry’s/cls «c\c cic. o1s)0\clole'e)e\e elelelele) elslclolelc eleleleiciele elejelelej lao
Graphic representation and spectral param-
eters versus azimuth for artificial surface
composed of two identical orthogonal trends..........204
Graphic representation and spectral param-
eters versus azimuth for artificial surface
composed of two orthogonal trends of differ-
ent spectral levels....cccccccccccccccccccccccccsces so 200
Graphic representation and spectral param-
eters versus azimuth for artificial surface
xv
composed of two orthogonal trends of differ-
entmspecitralmisilopelaioisietereieleleislolelciclelelcleleleicleleieteietelsicicitererezOS
Graphic representation and spectral param-
eters versus azimuth for artificial surface
composed of two orthogonal trends of differ-
ent spectral slope and intercept....cccccccccccccsccere2l™
Location of two-dimensional SASS bathym-
etry from the Mendocino Fracture Zone.....cccccccccee2ll
Graphic representation and spectral param-
eters versus azimuth for Mendocino Fracture
ZONE teleloleveleleleleleleieloleterereverelelelelererciclereleleioierereievavelerersielevereieniemima les
ier
Denon
hae hs a
A i \ bs
HAS scien
a
Hawi peer |
oa
Hebe oaay visinanas We er sy TOE
area: a pesinapanseenl csPoeiban dit)
A Poa ae. 08 ea ig eee
Ta Lestnooy hea: seldeciaamaet
| coer
ee et ee ame
rs nea ae Le HG Khe TaN A ie ae eetamico!
aa we " is ows Past yapiibe sui bie
ipl) bil Dtad i* qe le ‘eka ‘ib sie
‘ei untae gate
Beet oe hows eee ws q apy rev
‘i ia qinaantae’ Wms ne wedi RY big eee
prea sa 2 a Linaumiucieady WA ae "4 Se Soke tlh *
Ne ay
iat Bate
Hacc is COR NE aipuk: wih | iteron ian?
va, ¢ SRN ¢ Balt MS a! Sees ‘iM, yy We un,
Savi ay sah bie or nN hee +e of
colt iter elie
Aen Men ata ie MAREE ER IN tye Wie it or
' eC EW os 4'e a Piet en es
WF ahs Ch et ers bik hid 4 epee ene
Junius at sorta E hid yshinsngpro paint
cavapollentel eo iach Deans
ine ane iad ake? etsy ee oh
a ; a8 tigi Zi:
‘Aa heewke
1. Introduction
The contour map of oceanic depths, or bathymetric chart, has pro-
vided a fundamental tool for inferring deep-sea processes, both sedi-
mentological and tectonic. The identification of the major features of
the sea floor led to more elaborate geophysical studies and such unify-
ing theories as sea-floor spreading. The efforts of numerous institu-
tion. world-wide have produced comprehensive contour charts of global
bathymetry, which have formed the basis of further geophysical survey
efforts.
The method of contouring has the mathematical equivalent of fitting
a continuous surface to discrete three-dimensional data. In the case of
bathymetric contours, the discrete data are usually in the form of vari-
ably spaced and randomly oriented ship tracks along which are discrete
soundings at some interval. Modern soundings are normally estimated
from surface ships using acoustic sounders in which the two-way travel
time of a pulse of sound is interpreted as a depth at a discrete point.
Geometric spreading of the sound with depth causes a large area of the
sea floor to be insonified (often called the “footprint”), and the meas-
ure of the first significant return of sound to the ship is interpreted
as a sample of the shallowest depth from this large insonified area.
The uncertainties in the precise location of the measured depth, the
sound velocity structure of the water column and the positioning of the
ship at sea, combine to introduce a noise component into the output sig-
nal. The random noise of measurement, in addition to the uncertainty of
interpolation between often widely spaced samples and ship tracks,
imposes a physical limit on both the vertical (depth) accuracy and the
spatial resolution possible in utilizing bathymetric contouring methods.
Despite its inherent uncertainties, a bathymetric contour chart
does provide an absolute estimate of depth at all points in space. In
the terminology of numerical modelling, a contour chart can be consid-
ered a “deterministic model” of sea-floor topography. Because such a
model is analogous to fitting a least-squares surface through randomly
spaced, noisy data, it is by nature a smoothed representation of the
actual “submarine topography. The degree of smoothing required (or
equivalently the cutoff frequency for low-pass filtering of spatial
frequency) depends upon the accuracy, resolution, and density of data
used in the contouring. Van Wyckhouse (1973) demonstrated the deter-
ministic aspect of bathymetric charts with the creation of SYNBAPS
(Synthetic Bathymetric Profiling System), a data base containing depths
on an evenly spaced grid. The grid cell spacing of 5° of latitude and
longitude used for SYNBAPS, seems to represent a workable estimate of a
suitable interpolation interval for deterministic modelling. Recent
work at the U.S. Naval Oceanographic Office has extended this 5' grid
world-wide.
There are also some practical considerations, in addition to the
physical limitations, in determining submarine topography to high spa-
tial frequencies. Survey instruments such as deep-towed sleds incorpo-
rating stereo photography, narrow-beam sounders and side-scan sonar, are
able to map small areas of the sea floor with high resolution. However,
it is practically impossible to extend these surveys globally, due to
the operational difficulties and the great expense of these methods.
Another practical difficulty of extending deterministic, gridded models
to a smaller grid cell spacing is the storage requirements of the data.
For example, the 5° gridded data set extended from 70°S to 70°N being
developed by the Naval Oceanographic Office will require approximately
five million storage locations. To extend this grid to the order of 100
meters spacing (assuming this was determinable globally) would require
approximately 5 x 1010 storage locations.
In light of the physical and practical limits of describing high
spatial frequency sea-floor topography deterministically, the apparent
alternative is some probabilistic (stochastic) model describing the dis-
tribution of features in a statistical sense rather than determining the
exact locations of depths. To be useful, this stochastic model must
describe, with a reasonable number of parameters, most of the true vari-
ability of a region of sea floor. Using a stochastic representation of
high frequency feature distribution in combination with the lower fre-
quency deterministic models being developed, an essentially complete
description of the sea floor is possible. The development of such a
stochastic model of sea-floor topography is the intent of this
dissertation.
ie Aarts a renin f s tha) a see = senna k
34) ,,
eres ay. tS
oa ' wean in
tan Pote seit cle
ae "bs
=a wee - , sienna
bg i J
Teh
2. Applications
In addition to the intellectual satisfaction of completing the des-
eription of sea-floor depths, a stochastic model of high spatial fre-
quency submarine topography, or surface roughness, has many practical
scientific and engineering applications. Many of the most direct appli-
cations are in the field of underwater acoustics. Clay and Medwin
(1977) provide one of the best physical descriptions of the interaction
of an acoustic signal with a rough surface. Matthews (1980) reviews the
importance of relative scale in acoustic bottom interaction. Briefly,
the relative spatial frequencies of the incident acoustic signal and the
bottom roughness determine whether the surface acts as a reflector or as
a scatterer of energy. This relationship illustrates the necessity of
describing bottom roughness in terms of spatial frequency.
The importance of acoustics to marine geophysical surveying systems
cannot be overstated. Sea-floor bottom loss of sonar and seismic sys-
tems, ranging of side-scan sonars, and signal strength of outer beams on
multibeam sounders all depend heavily on bottom roughness. In naviga-
tion applications, the roughness of the sea floor has an impact on the
backscattering of Doppler sonar, as well as determining the background
noise for navigation by bottom features. A major application for the
U.S. Navy of bottom roughness information is as environmental input into
long-range acoustic propagation models, used in submarine and surface
ship tracking. For engineering applications, the a priori knowledge of
the spatial frequency composition of the bottom could be used in
computer-aided design of future sensing systems, as well as aid in main-
taining coherent signals for underwater communication systems.
Another field with significant applications for sea-floor roughness
information is physical oceanography. This requirement led to the
extensive work of T. H. Bell. Planetary Rossby waves are strongly
affected by bottom roughness in long wavelengths (Rhines and Bretherton,
1973). The propagation of long surface waves such as tides and tsunamis
Bee also affected (Rhines, 1977). Bell (1973, 1975a) showed that the
interaction of deep ocean currents and bottom topography may lead to the
generation of internal gravity waves in the oceans, a major influence in
ocean dynamics as well as submarine operations.
Another geophysical application is in the general field of survey
design. Davis (1974) has formulated a method which, with a knowledge of
the spectral content of the field being measured, allows a predetermined
survey accuracy to be attained. Briefly, the two-dimensional (or in
some applications, three-dimensional) spectral content estimates are
used in algorithms which prescribe preferred track spacing, sampling in-
tervals and track orientation. The method has been successfully applied
to gravity, magnetic and physical oceanographic surveys. The availabil-
ity of an adequate spectral content model for submarine topography would
make this technology available for bathymetric survey design.
In the detection of anomalous features in a field, the “normal”
background variability must be removed by filtering to aid detection.
This method has been applied successfully in magnetic anomaly detection
and theoretically could be applied to bathymetric anomaly detection.
McDonald, Katz and Faas (1966) applied this concept to submarine detec-
tion, specifically in the search for the nuclear submarine Thresher.
There are many applications of sea-floor topography modelling to
the geologic interpretation of seafloor processes. Deterministic models
illustrate that sea-floor relief is often decidedly anisotropic. This
investigation will provide insight into the systematic nature of such
spatial patterns and at several scales. Hayes and Conolly (1972) demon-
strated the value of spectral techniques in studying the complex topog-
raphy of the Australian-Antarctic Discordance. Since many geological
processes result in linear features of characteristic wavelengths, an
accurate, frequency-dependent, stochastic model could provide the capa-
bility of decomposing such features and delineating their geographic
distributions. It is anticipated that many important but totally
unforeseen relationships will be discovered in exploiting this higher
frequency model in the same manner that numerous fundamental relation-
ships were discovered through the generation of global bathymetric
charts.
anes Kec aretsen is ves, dase
Pai AGE
Whe
3. Previous Work
Although a great deal of quantitative geomorphology has been done
on terrestrial landscape, as well as on lunar and planetary landscapes,
relatively little quantitative study has been done on sea-floor morphol-
ogy. Only since the late 1950's, when acoustic sonar systems became
commercially available, has it been possible to attempt such statistical
studies of bathymetry. The continuing refinement of sonar and naviga-
tional systems has given recent investigators even greater opportunities
for success. Equally important has been the development of large dig-
ital computers and efficient statistical algorithms such as the fast
Fourier transform, which allow sophisticated statistical analyses on
large volumes of data, which were impractical until recently.
The work of previous investigators is somewhat sparse and inconsis-
tent, and is reviewed only briefly here. More information is given in
later sections as it becomes pertinent to various aspects of the study.
Some of the earliest studies were done by Agapova (1965), who generated
mean, variance, skewness, and kurtosis statistics from measurements of
slopes of a transect of the mid-Atlantic ridge. Heezen and Holcombe
(1965) were able to describe physiographic provinces over a large por-
tion of the North Atlantic Ocean. After rejecting spectral and filter-
ing techniques, these authors developed a method of comparing the spac-
ing of adjacent peaks and troughs. In effect, the method calculates the
average distribution of slopes without regard to spatial frequency.
Krause and Menard (1965) studied depth distributions from 15 pro-
files in the east Pacific Ocean and found them to be normally distrib-
uted. In studying normalized autocorrelation functions derived from the
same profiles, the authors found no regularity with respect to lag. The
profiles were used for delineating provinces on the basis of height ver-
sus width ratio of the abyssal hills. Smith et al. (1965) continued the
work of Krause and Menard (1965) and concluded that the distribution of
depths is Gaussian when observed in distinct wavebands of 2 to 16 nauti-
cal miles. Larson and Spiess (1970) later used a deep-towed instrument
package to study the distribution of slopes in a small area of the east-
ern North Pacific. Krause, Grim and Menard (1973) generated simple
cumulative frequency plots of slopes in two areas of the East Pacific
Rise. They found a very consistent power law form to describe these
Ai Pronto and concluded that marine geomorphology could be described
beers a few parameters.
: Neidell (1966) generated spectral estimates of bathymetric, mag-
netic and gravity profiles from the Atlantic and Indian Oceans. All
spectra showed a "red-noise” character, that is, a decrease of power
with increasing spatial frequency. The comparison of the various
spectra was shown to be a valuable tool in investigation of complex
geophysical problems. McDonald and Katz (1969) in their study of the
directional dependence of roughness developed a polar autocorrelation
function. Hayes and Conolly (1972) used spectral analysis as an inter-
pretive tool in an area south of Australia. Distinct linear trends were
interpreted by projecting randomly oriented tracks into north-south and
east-west profiles and investigating the consistency of the resulting
spectra.
Clay and Leong (1974) rigorously described the relationship between
surface roughness and the coherence of acoustic reflections. The dura-
tions of returned pulses from hull-mounted sonars were used to estimate
RMS roughness of microtopography (<0.2 km). Histograms of spectral
estimates (periodograms) were generated by hand and shown to map con-
sistently in an area southwest of Spain.
Bell (1975b), also analyzing data from the eastern North Pacific
Ocean, used power spectral techniques to generate composite estimates
from several sources. He discovered a functional relationship for power
versus spatial frequency of the form ax (power law form) to be consis-
tent over a large range of enon scales. In a later paper (Bell,
1978), this relationship was shown to hold for an enlarged data base,
which is also confined to the same geographic area. The importance of
anisotropy was recognized, and an initial look at the aspect ratios of
submarine features was presented.
Berkson (1975) generated spectra from a variety of bottom types and
attempted to fit these with a power law form. Although a wide variety
of coefficients were calculated, the power law form seemed to be consis-
tent over many types of topography. Akal and Hovem (1978) generated
two-dimensional spectra of sea-floor roughness from two sets of stereo-
pair bottom photographs and a contoured bathymetric chart. They noted a
remarkable consistency of form in all three spectra. Matthews (1980)
developed a deterministic approach to describe bottom roughness. The
North Atlantic and North Pacific Oceans were divided into 30 x 30 nau-
tical mile squares and the maximum relief calculated. These cells were
then grouped by range of relief (0-1100 m, 1100-1900 m, >1900 m) and the
results mapped. Recently, Naudin and Prud'homme (1980) quantitatively
described bottom morphologies from several areas based on multibeam
sonar data collected by the SEABEAM system. Most recently, Berkson and
Matthews (1983) have extended the work of Berkson (1975) and included
estimates of sediment-basaltic interface roughness.
10
4. Statistical Considerations
In order to develop a useful model of sea-floor roughness, one must
select, from a seemingly infinite variety of available methods, an
approach which is both suitable and tractable. In addition to consider-
ing the nature of the sea floor itself, one must also consider such
aspects as data resolution, computer storage capacity, statistical
validity, and compatibility with various applications. Often, the
choice of a particular method involves trade-offs between several of
these considerations. In the following sections, many of these funda-
mental considerations are addressed, and an initial approach to devel-
oping this particular model is presented.
Statistical Measures of Roughness
If one considers surface roughness to be the variability of heights
(or depths), the realm of statistics offers a multitude of measures to
describe the roughness of a surface. Perhaps the most fundamental sta-
tistic available to describe roughness would be one of the standard
measures of data dispersion, such as the root mean square, standard
deviation, or variance. These measures have the advantage of producing
a simple parameter to describe the variability of depths in a given sam-
ple. The major disadvantage of such simple measures is that they do not
| provide information for roughness in terms of wavelength, and the sta-
tistic decivedudecends upon the sample spacing and length of sampling as
well as the actual roughness of the surface.
11
To illustrate the importance of having control over frequency
dependence, consider two examples. If a generally flat surface which
contains a high frequency roughness component of wavelength 4, were sam-
pled at spacing A or any integer multiple thereof, each sample would
fall at precisely the same depth and yield a variance of zero. Mathe-
matically, all values of X would be identical, therefore the mean
= n
».4 =+ ie 2a would be equal to all X‘s and therefore
roa (Xq - x)?
Var (Xx) = eal = 0
Although this example presents an extreme case, it is obvious that to
assure an accurate measure of the variance at a given frequency, the
sample spacing must be less than that corresponding wavelength, and the
sample length long enough to sample all portions of the cycle.
A more important shortcoming of these standard dispersion measures
occurs at the longer wavelengths. Consider a generally smooth but
broadly sloping surface. Examples from the deep sea might be a conti-
mental rise or an abyssal fan. Since these dispersion measures record
the average variation of individual samples from the sample mean, it is
obvious that a relatively long sample would span a greater range of
depths and produce a larger dispersion statistic than a smaller sample
located in the same data. In this case, the measure of roughness would
depend largely on the sample length.
Many acoustic models of bottom interaction use the more general
dispersion measure of the root mean square (RMS) roughness. In these
acoustic models, this value represents the RMS variation of depth about
some predicted value, normally the smoothed bathymetry.
RMS roughness = [
where 2S represents a predicted value of depth at point i. By measuring
the roughness relative to a predicted depth, rather than a simple mean
as in the case of the standard deviation, the effect of long wavelength
slopes on sampling is reduced. This improved measure does not, however,
provide control over the distribution of roughness with frequency.
Since the reflection or scattering of an incident acoustic signal on a
surface is dependent upon spatial frequency, this often used parameter
appears inadequate.
Another possible measure of roughness is termed the “roughness
coefficient” by Bloomfield (1976), and has the form
2
nl iy - x, _,)
it oe tol
ao ee i
121 6% > ©)
Because this measure (also referred to as the von Neumann ratio and the
Durbin-Watson statistic) is normalized by the total sum of squares of
the residuals, the dependence on sample length is minimized. However,
because this measure also depends on the squared differences of adjacent
points, it measures only the variability of the signal at a wavelength
corresponding to the sampling interval. This roughness coefficient then
is essentially a ratio of the high frequency variability of a signal to
its long-term trend, an interesting statistic, but not adequate for a
complete stochastic model.
Several statistical functions exist which describe the sample vari-
ability as a function of discrete data spacings or lags (see Chatfield,
13
1980). Perhaps the simplest function of this type is the mean cross-—-
product of terms at given lags k
CK = =k ee Xq - Xg—ee = EL(X_)(X4-K)]
This function (called simply the autocorrelation function in electrical
engineering literature; see Bracewell, 1978) is both unnormalized and
uncentered (the mean is not removed). Because it depends on absolute
magnitude values (in this case, the total water depth), this simple
measure can be improved for the purposes of roughness modelling by
removing the sample mean, which yields the autocovariance function
v(k) = E[(Xq - X) (X4-~ - X)]
It is obvious that at a lag of k=0, (y(0)) is simply the sample vari-
ance. By normalizing the autocovariance by the sample variance y(0),
one derives the normalized autocorrelation function
p(k) = ¥(K)/\ (9)
Notice the equivalence between the normalized autocorrelation function
at lag k=l, (p(1)) and the roughness coefficient described previously.
In comparing the roughness of two different samples, it is undesirable
to normalize by the sample variance, therefore, the autocovariance func-
tion provides the best measure of variability with frequency (centered,
but unnormalized).
14
Fully describing the variability of a process with its autocovari-
ance function requires values at all n lags. By taking the Fourier
transform of either the autocovariance or the unnormalized autocorrela-
tion function, the process can be expressed more concisely as a function
of frequency. This measure is the well known power spectral density
function and can be estimated directly from the data with Fourler trans-
forms. Besides providing a concise expression for the roughness as a
function of frequency, the power spectrum also has many useful proper-
ties which are described in Chapter 6 of Bracewell (1978). One theorem
of particular interest is the derivative theorem which states that if a
function f(x) has the transform F(s), then its derivative f'(x) has the
transform i27sF(s). In this application, given the power spectrum of
depths as a measure of roughness, the distribution of slopes (first
derivative of depth) can be directly calculated.
An additional measure of bottom roughness, which is particularly
favored by those interested in acoustic modelling of bottom interaction,
is the two-point conditional probability distribution function
P(hy yhg|ry,F2)- This function defines the probability of measuring two
heights (h, and hj) given two distance vectors (ry and FQ). Two consid-
erations make this approach intractable. First, the description of the
function requires a large number of parameters to be retained. Sec-
ondly, complete two-dimensional data are required to adequately generate
the function. This type of survey data is rarely available and only in
areas which have been surveyed using multibeam sonar.
The power spectral density function appears to be the best choice
of statistical measure for bottom roughness, and it will underlie the
stochastic model generated in this study. For convenience, the ampli-
15
tude spectrum (square root of power spectrum) will be used. When
properly normalized, the amplitude spectrum allows the amplitude of
component sinusoids to be expressed in simple length units, which can be
interpreted more easily than units of length-squared. The method of
calculation will follow Davis (1974) with proper windowing, prewhiten-
ing, etc. As will be discussed in later sections, the simple and
consistent functional form of spectra of topography allows relatively
easy description and manipulation of the model. Also, recent work by
Brown (1982, 1983) has shown the value of working in the frequency
domain in modelling scattering from rough surfaces.
Validity of Measurement Over Large Areas
In generating the variance, autocovariance, or power spectrum from
a discrete sample, only one realization of an infinite number of possi-
ble realizations within the population is observed. The degree to which
this realization is valid over the entire population depends upon the
degree of homogeneity of the process. The condition of spatial homoge-
neity is generally known as “stationarity” and the population being
described referred to as a “stationary process”. The term “process” is
used in the statistical sense of the variation of data with either time
or space. In the case of sea-floor topography, the depths vary as a
function of space.
Stationarity is normally defined in two ways (see Chatfield, 1980;
Popoulis, 1962). A spatial series is said to be “strictly” (or "“first-
order”) stationary if the joint distribution of the process does not
depend upon position. This implies that the mean and variance do not
16
depend on position. The less restrictive definition of “weakly” (or
“second-order") stationary processes requires the mean to be constant
and the autocovariance function to depend only on lag, not on position.
This second definition has somewhat greater application, however, it is
still too restrictive for general use with a spatial process as varied
as submarine topography.
Consider a broadly sloping surface, such as an abyssal fan. The
mean depth in this province would by definition vary with position, and
therefore would not satisfy the first requirement for stationarity. Yet
if the process is homogeneous in higher spatial frequencies, one might
prefer to treat this province as a homogeneous area for modelling. In
practice the existing deterministic model could be used to describe the
low-frequency trend. The sample data could be high-pass filtered to
remove the non-stationary trend before generation of a spectrum.
The presence of a low-frequency trend in samples of geophysical
data is quite typical. In almost any length sample of a natural proc-
ess, there are frequency components present with wavelengths greater
than the sample length. In most natural systems, there is a finite
limit on the rate of change of the process, causing the longer wave-
length components to be also of greater amplitude. This typical spec-
trum of natural processes was termed a “red-noise” spectrum by Shapiro
and Ward (1960), an analogy to the red color of low frequency visible
light.
Although the red=-noise spectrum is the usual form in natural sys-
tems, Figure 4-1 illustrates schematically that the presence of non-
stationary components (in this case the statistical mean) can occur at
any frequency and is dependent upon the horizontal extent of the sample.
17
C
Figure 4-1
Schematic illustration of the scale dependence of station-
arity. All traces are derived from the same function (B),
however the inferred stationarity of the mean would be dif-
ferent depending upon the scale of observation. Observed
at the scale of frame A, the mean is obviously non-station-
ary. However, the mean appears stationary at the larger
scale of B or the smaller scale of C. In most natural sys-
tems, the non-stationary components occur in the partially
resolved low spatial frequencies. This difficulty can be
handled analytically by defining homogeneous provinces on
the basis of stationarity and prewhitening of the spectra.
18
In order to generate a stochastic model of sea-floor topography, we must
ensure that the process is stationary within some limit and over some
defined area, with respect to the statistical parameters being calcu-
lated. Ideally, this is accomplished by identifying homogeneous prov-
inces based upon these criteria. Davis (1974) developed such a “prov-
ince picker” for use in geophysical survey design, and his method is
illustrated in Figure 4-2.
By actively defining provinces which are weakly stationary in the
frequency band of interest, one improves the validity of the statistical
measures generated within each province. In addition, by delineating
province boundaries, one can also alleviate the problems associated with
the least-squares or averaging nature of Fourier transforms. In gener-
ating a power spectrum from a data set, the operator must select the
length and location of the sample series to be transformed. The result-
ing spectrum will reflect the average frequency composition over the
sample. If the sample spanned two distinct statistical processes, the
result would provide the average composition of the two provinces, and
would accurately represent neither. By confining one's samples within
the boundaries of a homogeneous province, one insures a valid and repre-
sentative statistic. These concepts will be discussed in detail in
Chapter 5.
Functional Representation of Spectra
One property which makes the Fourier transform, both continuous and
discrete, such a powerful analytical tool is its ability to express a
spatial process in the spatial frequency domain both exactly and com-
Band-Pass Filter (A) =(B)
Low - Pass Filter (C) =(D)
Contoured (D) = (E)
Figure 4-2 Schematic illustration of the method developed by Davis
(1974) for delineating homogeneous roughness provinces
within a given spatial waveband. The original signal (A)
is band-pass filtered at a pre-selected wave band of inter-
est to produce B. The now centered data of B is full-wave
rectified by taking the absolute value of all terms to pro-
duce C. This output is low-pass filtered (smoothed) to
yield the continuous function of D. Finally, this function
is contoured at some selected interval (not necessarily
linear) to delineate the boundaries of provinces X, Y, and
Z in frame E.
20
pletely. This complete correspondence between domains allows detailed
analyses to be made on the spectra with the assurance that there is an
exact analog in the spatial domain. In order to maintain this corres-
pondence, it is necessary to retain the complete transform, (both ampli-
tude and phase components) in the spatial frequency domain. In the case
of discrete data and transforms, it would be necessary to retain all of
the degrees of freedom present in the original data set.
In the creation of a probabilistic model, it makes little sense to
retain as much information in the model as was present in the original
data. Presumably, one would prefer to use the original data as a deter-
ministic model. Also, the purpose of the model is to describe the gen-
eral high frequency structure of the sea floor, rather than to analyze
in detail the individual frequency components. Finally, the restric
tions of computer storage space require a limited number of parameters
in the model.
All of the above considerations argue strongly for a severe paring
of information in the frequency domain model. The phase spectrum, which
requires fully half of the information in the spatial frequency domain,
defines the origin in space of all component sinusoids of the amplitude
spectrum and adds very little insight into the general structure of the
sea-floor. An analysis of bathymetric phase spectra presented in Chap-
ter 5 will show that the phase is in fact randomly distributed. For the
purpose of this model, no phase spectra will be retained, as none of the
previously stated applications require phase information.
All measured data contain a component of random noise. Remotely
sensed data are especially susceptible to measurement noise, the partic-
ular noise problems in measuring oceanic depths having been mentioned in
21
the introduction. The presence of noise in the spatial data manifests
itself in two ways in spatial frequency spectra. Inaccuracies of verti-
cal measurements in the spatial domain result in the presence of a hori-
zontal “white-noise level” in all amplitude or power spectra. This
problem will be treated in detail in the following section. Uncertainty
in the location of features in the spatial domain results in the scat-
tering of amplitude estimates about the true frequency spectrum. These
distinctions in the sources of error are somewhat artificial since the
vertical and horizontal uncertainties are interdependent.
Figure 4-3 reproduces a typical amplitude spectrum of depths. The
red-noise character of the distribution as well as the scattering of
amplitude values is apparent. The spectrum was derived from data col-
lected by the U.S. Naval Oceanographic Office using SASS (Sonar Array
Subsystem), and represents the highest resolution bathymetric informa-
tion currently available from a surface ship (see Glenn, 1970). Control
over relative horizontal location (navigational accuracy) of soundings
is especially good due to the use of large, stable surveying platforms
(in this case, USNS Dutton) and SINS (Ship's Inertial Navigation
System). The degree of scattering of amplitude estimates would pre-
sumably be greater in less sophisticated systems.
Several methods come to mind to smooth this somewhat noisy spec-
trum. A simple moving average taken over the amplitude estimates in the
frequency spectrum would smooth the data. However, information would be
lost from the high and low frequency extremes of the spectrum, while the
density of data in the intermediate frequencies would not be reduced.
The use of spectral windows or lag windows could be used both to smooth
the spectrum and decrease the data density. The use of data windows in
22
DEPTH (FATHOMS)
104
103
102
10!
AMPLITUDE (FATHOMS/CYCLE/DATA INTERVAL)
Figure 4-3
i120 169225) 28 S793 4490506 mt 562 O18 NNN674
CUMULATIVE DATA POINTS DOWNTRACK
10-2 10°! 10° 10!
NORMALIZED FREQUENCY (CYCLES/DATA INTERVAL)
Illustration of a typical amplitude spectrum of sea floor
topography. A profile of the data, collected by SASS on
the Gorda Rise, is presented above. The spectrum shows
clearly the power law form as a linear fit on log-log plot.
Possible extrapolations of the trend are shown as dashed
lines beyond the fundamental and Nyquist frequencies.
23
the smoothing of spectra is discussed in many texts (see for example
Bloomfield, 1976) and will not be discussed in detail here. The spec-
trum presented in Figure 4-3 was in fact generated using the method of
Davis (1974) which utilizes prewhitening as well as low-pass filtering
of the spectral estimates.
Another method of smoothing this rough spectrum is through the use
of regression techniques. Calculating a continuous mathematical func-
tion to describe the distribution of amplitudes would produce a smooth
representation of the spectrum while, depending upon the complexity of
the function used to fit the data, greatly reducing the number of param-
eters retained in the model. This least-squares representation will be
used in this study.
By describing the spectrum with a simple, continuous mathematical
function, one can more easily take advantage of the many symmetry prop-
erties of the Fourier transform described by Bracewell (1978). For
example, Rayleigh's Theorem (or Parseval's Formula for discrete series)
states that the integral of the power spectrum equals the integral of
the squared modulus of the function, or
fP \£(x)|? ax = f° |F(s)|? as
=O —o
This is equivalent to stating that the total energy in one domain is
exactly equal to the total energy in the other. If one were interested
in the total energy in a particular band of frequency (in order to study
the bottom interaction of sound of a particular wavelength, for exam-
ple), the high and low frequencies of the band pass would become the
limits of integration of the power spectrum. With the spectrum repre-
24
sented as a simple function, this definite integral could be evaluated
analytically to derive the RMS variability for a discrete waveband.
Having decided to use functional representations as the basis of a
stochastic model of sea-floor roughness, the selection of a suitable
functional form for the spectra becomes crucial. In order to minimize
the size of the model, the simplest functional form which is justified
by the data should be the best. An examination of Figure 4-3 as well as
many other spectra presented later, would suggest the use of a simple
straight line fit to the data. In light of the scatter in this already
smoothed data, no higher order functional form is justified.
The normal form for fitting a straight line to data is the estima-
tion of the coefficients a and b in the equation
a a
y=bxta
Notice in Figure 4-3, however that the data are plotted versus logarith-
mic scales. The regression equation would therefore be written,
log A =b log s + log a
where A = amplitude and s = frequency. This equation can be rewritten
in terms of A as,
Aza. sd
This inferred relationship between amplitude and frequency is often
termed a “power law” or “power function” relationship. Appropriate
25
regression techniques must be selected in order to assure a proper
regression fit to the power law form. Of prime concern is whether the
residuals are minimized in log-log space or linear-linear space. The
methods used in this study with accompanying theory and computer soft-
ware are presented as Appendix A.
The power law form seems to represent the best model for describing
sea-floor topography in the spatial frequency domain. Its simplicity
permits the frequency structure of a sample profile to be described
using only two parameters, a considerable reduction of the original,
deterministic data. Extensive work by Benoit Mandelbrot (1982) has
produced a theoretical basis for this consistent relationship, formu-
lated in terms of fractal dimension, a parameter functionally related to
coefficient b above (see Berry and Lewis, 1980). Bell (1975b) discov-
ered the same power law relationship in data from the Pacific Ocean that
including lower spatial frequencies than those of interest in this study
(see Fig. 4-4). Notice that Migure 4-4 plots power spectral density,
rather than amplitude as in Figure 4-3 and other example spectra.
Because the vertical axis in both plots is logarithmic, this exponentia-
tion appears graphically as a linear transformation. Mathematically,
the squaring of amplitude represents a simple doubling of the slope,
1.e., multiplication of the b term by a factor of two.
Bell (1975b) finds a fairly consistent relationship at many size
scales, with the slope of the log transformed power spectrum of 6 =
=2. Although this value probably represents a good mean estimate, the
examination of many spectra in this study will show an accountable var-
fation in this value. Berkson (1975) also discovered a significant
variation of regression coefficients for spectra generated from differ-
26
F(K), m*/(cyc/km)
Figure 4-4
BALMINO ET AL. (1973)
WARREN (1973)
BRETHERTON (1969)
COX & SANDSTROM (1962 )
ABYSSAL HILLS
10°§ 1073 107% to 861072 10° §=610° ~~ §610! 102
K, cyc/km
Spectral estimates collected by T. H. Bell (1975b) from
various sources. The linear distribution of the data
illustrates the power law form of submarine topography
spectra for a wide range of scales. Notice that this spec-
trum is plotted versus power (amplitude squared) which
results in a doubling of the ordinate and the slope of an
amplitude spectrum. Refer to Bell (1975b) for data
sources.
27
ent bottom types. Analyses presented in Chapter 5 will link this
“universal” relationship
a a
Aza.s
to the absence of stationary provincing techniques.
Extension of Spectra into High Frequencies
The component of random noise present in all empirical data,
includes uncertainties from many possible sources in the total meas-
urement system. There are uncertainties associated with the measuring
device itself, for example, the errors in the timing of the arrival of a
sonar pulse. The interference of external sources, such as radio waves,
may also affect the measuring instruments. The nature of the environ-
ment between the detection device and the process being sampled may
introduce error, such as the variability of sea water sound velocity in
bathymetric surveying. The truncation of significant digits in the
storage of digital data introduces a finite level of noise, often called
“round-off error.” All of these sources combine to form a total noise
level for a measured data set.
With the exception of round-off error, which is calculable, the
source of these random errors can not be decomposed. However, using
spectral techniques the level of the total noise can be estimated. It is
a well-known property (see Bracewell, 1978, Chapter 16, for a complete
discussion) that the spectrum of a randomly generated signal varies
28
about a constant value. This makes intuitive sense, since one would not
expect any particular frequencies to dominate a random series. Again in
analogy to the spectrum of visible light, this random component of a
signal is often called “white noise”, reflecting the equal contribution
of all frequencies.
In the same way that Rayleigh’s or Parseval's theorems can be used
to relate energy in one domain to energy in the other, the level of
noise in a spatial signal can be estimated from the amplitude of the
noise level of an amplitude spectrum. The noise level in the spatial
domain is normally expressed as a simple dispersion measure of the vari-
ability of the data, in this case, the oceanic depths. The following
formula relates the white-noise level of the power spectrum to the root
mean square.
RMS = ¥ P/na
where n = number of data points in series; and P = the mean power level
of the white noise. In the case of sonar systems, knowledge of this RMS
level defines the resolving capability of the system for a given signal
level. Using these simple spectral techniques, the resolving power of
the various sounding systems in use today can be calculated and
compared.
The red-noise structure of natural systems has been mentioned and
illustrated previously (see Fig. 4-3). The interaction of natural sig-
nals with instrument noise takes the form of a decreasing red-noise
spectrum of the signal “intersecting” the horizontal white-noise level.
In the spatial domain perspective, lower frequency features tend to have
29
higher amplitude and therefore the signal in these frequencies is “vis-
ible” above the RMS noise. The frequency at which the spectrum of the
signal intersects the noise level represents the highest frequency that
is being resolved in the system.
Figure 4-5 illustrates these concepts on a spectrum of SASS data.
Notice that Figure 4-5 shows an obvious noise level while the spectrum
shown in Mgure 4-3 does not. Both sets of data were collected using
the same sounding system within days of each other, and it is expected
that the instrument noise in each is approximately equal. The differ-
ence is therefore that the data in Figure 4-3 represent a rougher area
in which the signal energy in the highest frequencies is sufficient to
maintain the resolution of the signal above the noise. The noise level
is present in both sample spectra, however, it was never intersected in
the sample from higher energy sea floor. The ability of a sonar to
resolve horizontal features depends not only on the horizontal resolving
power limitations of the instrument (normally limited by the size of the
“footprint”), but also by the amplitude of features present in the sea
floor.
In light of the limitations of noise in all measurement systems, a
fundamental question in the development of a stochastic model of sea-
floor roughness from widely available surface sonar, is how far into the
high frequencies the model can be extended. One could argue that due to
the persistent exponential form of the spectra noted in this study, as
well as the work of Bell (1975b), that this functional form can simply
be mathematically extrapolated into higher frequencies. Further justi-
fication might be provided by the generation of spectra from deep-towed
instrument packages, bottom photographs, and direct observations, to
30
8
8
eee
SIRI IN PL ELL LL
DEPTH (FATHOMS)
Sas
35 71 106 141 177° 212 247 283 «#4318 353 ©6388 «9424
CUMULATIVE DATA POINTS DOWNTRACK
AMPLITUDE (FATHOMS/CYCLE/DATA INTERVAL)
10°73 10-2 10°! 10° 10!
NORMALIZED FREQUENCY (CYCLES /DATA INTERVAL)
Figure 4-5 Illustration of an amplitude spectrum of sea floor topog-
raphy which has encountered the "white-noise" level in the
data. The spectrum of random noise has the form of a hori-
zontal line. This noise level, which was derived from data
collected by SASS, corresponds to a random component in the
spatial data with root mean square dispersion of 1.9
fathoms.
31
provide spot checks of the high frequency structure. This approach will
be pursued in Chapter 7.
Effect of Linear Features
In generating a Fourier transform of topography within a stationary
province, a one-dimensional statistic is generated from a two-
dimensional surface. Provided the surface is isotropic, that is that
there is no significant directional dependence, statistics derived from
a one-dimensional sample would be equally valid for any orientation.
Even a cursory examination of a deep-sea bathymetric chart reveals
clearly that the sea floor, at least in the lower spatial frequencies
presented in a contour chart, is quite anisotropic. There is extensive
evidence that bathymetric lineations aieotertat to some extent in the
higher spatial frequency roughness of interest to this study. To com-
pletely model sea-floor roughness, it is essential to account for ‘any
major directional dependence of the topographic features.
Figure 4-6 illustrates the effect of sampling a simple two-
dimensional periodic function in differing directions. Notice that the
true wavelength (A) of the feature is sampled only when the sampling is
exactly perpendicular to strike (8 = 0°). Any oblique angle produces an
apparent wavelength (A') which is greater than A. At ® = 90° (a sample
taken parallel to strike), the feature is not expressed at all (\' =©@).
The apparent wavelength is related to the true wavelength by
At = [costo] . A
32
Maes ve
i i h i fi
vii fou los Te A by
eine
Figure 4-6 Schematic illu of the effect of directional sam-
pus °¢ on spati a1 tre enue ncy estima ae The figure illus-
ae inc of snoaen eC aeyaien th of a
perio Sie ae ae perce ampled at any angle other an
ctly pe reper ndicular to at ike.
Notice that the amplitude of the feature is not affected, provided a
full wave form is sampled.
The effect of directional sampling in the spatial frequency domain
can be calculated using the previous relationship in combination with
the similarity theorem of Fourier transforms (see Bracewell, 1978,
p.- 101). The similarity theorem states that given the transform pair
f(x) > F(s), then
f(ax) > |a|~! F(s/a)
Applying the geometry for directional sampling of linear features, i.e.,
a= |cos 6|
£(|cos 8| . x) D |cos ellen - F (s/cos 9)
Because |cos 6| must always be less than one, the effect of oblique
sampling in the frequency domain is to shift the amplitude peak to lower
frequencies, narrow its width, and increase its amplitude. This is a
result of the fact that a signal of equal height but lower frequency has
greater power than its higher frequency counterpart. It is important to
note that these theorems assume an infinite continuous signal. In ana-
lyzing finite length signals, it is necessary to normalize the spectrum
by the sample length, that is, divide each amplitude estimate by the
number of values in the time series. This normalization preserves the
true amplitude of the individual waveforms, and allows comparisons
between samples of different length. The effect of anisotropy on sam-
pling a surface with continuous spectra is developed in a later section.
34
Several authors have recognized the importance of anisotropy of the
sea floor to power spectral studfes. Hayes and Conolly (1972) recog-
nized distinct peaks in their spectral analysis of the bathymetry of the
Australian-Antarctic Discordance. These groups were normalized to two
linear trends by projecting the data to simulate north-south and east-
west samples. Distinct trends were successfully identified, however, in
longer wavelengths than the high spatial frequency band of interest in
this study. Bell (1978) also recognized the importance of anisotropic
features in his study of abyssal hills in the north Pacific Ocean. One
result reported in his study of the aspect ratios of features in this
province is that the degree of anisotropy tends to decrease in the
higher spatial frequencies. Whether or not this is a true relationship
or the result of resolving limitations will be discussed further in
Chapter 6.
Figure 4-7 shows the sample locations for two spectra from the lin-
ear Mendocino Fracture Zone in the northeastern Pacific Ocean. Figure
4-8 presents the profile and amplitude spectrum for line A-A', which was
sampled perpendicular to strike, reflecting the long wavelength fracture
zone. Figure 4-9 is the corresponding plots for line B-B', which paral-
lels strike and is located in the zone of disturbance. The exponential
form of both spectra is evident, however Figure 4-10, which shows a com-
posite of the trend of both ‘spectra, illustrates the differences in
slope (exponent) of the two spectra. Segment A-A’ contains more power
in the low spatial frequencies, while segment B=-B' contains more power
in the high spatial frequencies.
Both spectra are valid representations of the spatial frequency
distribution in this physiographic province, but neither spectrum indi-
35
126°16'W-
40°
US IO ee 30’
126° 38'W
40°
126° 38'W 126°16'W
Figure 4-7 Bathymetric chart from the Mendocino Fracture Zone of the
eastern Pacific Ocean showing the location of the data used
in Figures 4-8 and 4-9. Notice that both samples lie
within the same bottom environment, but that A-A' jis sam-
pled across the long axis of the fracture while B-B' jis
sampled along the axis.
36
DEPTH (FATHOMS)
8
A 35 70 105 140 175 210 246 281 316 351 386 421 A!
CUMULATIVE DATA POINTS DOWNTRACK
105
z
> 104
[a4
ud
—
z
<I
= 10
Ta)
~
uw
O10?
>
O
“N
(Zp)
S 10!
Se
=
=
= 10°
lu
ra)
=)
=
: 10°!
10-2
10-3 10-2 10-1 10° 10!
NORMALIZED FREQUENCY (CYCLES/DATA INTERVAL)
Figure 4-8 Profile and amplitude spectrum of SASS data collected along
line A-A' in Figure 4-7.
37
DEPTH (FATHOMS)
%
B 63 127 190 254 317 381 444 508 571 634 698 761 B!
CUMULATIVE DATA POINTS DOWNTRACK
AMPLITUDE (FATHOMS/CYCLE/DATA INTERVAL)
10-3 10-2 10°! 10° 10!
NORMALIZED FREQUENCY (CYCLES/DATA INTERVAL)
Figure 4-9 Profile and amplitude spectrum of SASS data collected along
line B-B' in Figure 4-7.
38
10°
104
103
10°
10°!
AMPLITUDE (FATHOMS/CYCLE/DATA INTERVAL)
10-2
10-3 10-2 10°! 10° 10!
NORMALIZED FREQUENCY (CYCLES/DATA INTERVAL)
Figure 4-10 Composite of Figures 4-8 and 4-9 illustrating the anisot-
ropy of the sea floor in the Mendocino Fracture zone.
Although both spectra retain their power law form, there is
an obvious difference in slope and intercept. Profile B-B'
appears to contain more high spatial frequency energy,
while profile A-A' contains more energy in lower spatial
frequencies.
39
vidually describes the roughness in all directions. The need for a
directionally dependent function is obvious. The two-dimensional
Fourier transform might appear to be appropriate, since it describes the
two-dimensional roughness of a surface. However, its calculation
requires a complete two-dimensional grid of data values which is gener-
ally unavailable. The double Fourier transform, well described in Davis
(1973), is calculated from two orthogonal one-dimensional spectra. This
method, especially without the retention of the phase spectra, can not
unambiguously identify the orientation of trend. Both the two-dimen-
sional and double Fourier transforms require a large two-dimensional
matrix to be retained in the model.
McDonald and Katz (1969) describe a method for estimating autocor-
relation functions as a function of azimuth 9. A similar approach will
be attempted here for use with amplitude spectra. By studying the azi-
muthally dependent distribution of the coefficients derived from the
power law regression, a and b, it is anticipated that some functional
form or forms will be revealed to allow modelling of the entire process
as a function of both spatial frequency (s) and true azimuth (8). If
a=f,(®) and b=f,() then
F(8,8) = £,(6) . s*b(°)
With this functional form, input into the model of simple compass
direction for a given location would return the unique spatial-
frequency-dependent function coefficients for the amplitude spectrum in
that particular direction. Evaluation of these 9 dependent functions
should provide a simple measure of the degree and direction of bottom
anisotropy.
40
The derivation of the functional form of this 8 dependence would
best be investigated and verified in areas where complete two-
dimensional measures of depth are available. Areas surveyed with multi-
beam sonar systems such as SASS and SEABEAM are ideal for this purpose.
If a functional form is determined to be consistent, discrete samples
from randomly oriented tracklines could be fitted with this functional
model and F(®,s) estimated. Such studies will be presented in Chapter 6
and Appendix E.
One complication that could arise due to anisotropic roughness is
in the use of the province picker for the delineation of stationary
provinces. In a highly lineated area, the power level of particular
marrow frequency bands would be directionally dependent. The province
picker, however, measures total energy (integrated power) rather than
discrete power, and as mentioned previously, the peak shift due to
oblique sampling of a linear feature does not affect the total energy
measured. Unless a major spectral peak is shifted beyond the low cut-
off frequency of the band pass used in the province picker, the results
should not be adversely affected by directional sampling problems. Even
this problem seems unlikely to arise as thus far in all spectra sampled,
none have shown unusual low spatial frequency peaks such as are seen by
Hayes and Conolly (1972), and which are probably unique to a few tec-
tonic provinces such as the Antarctic-Australian Discordance.
41
ipkhew: ae ihn ¢ vt ‘¢ wha ei, very ie 1
‘eal Aleit, Me id sa Ne both f fl Aa Areouael 7 hah hai ‘“ ae ; r
ety el at: i ‘ail wae cae 7 ne mee ay
A eine We avy iw pea Nh orca he aah et re a | fu 7 ;
ait Dh mee | PE waa aantetiad ry
Wry dna! Hari Ae wah le sd boak oll ane
eo mmoaune ike omy hme atin’ |
tide ea Lain § CaN) ARR AN:
i emt qewe A dime i ieee aie pin
mest Gt we = ae year igi iin 1 bee
tones? mi
- Bibbstoauy 0 ley barat tbe ah a Jeena |
Rierienenie: sy gee S200 SM
nen bisa On butt pha he thot Lon Jae ov hain
ee obgen elk, ot ded ony to, nt, odd a Na
YN). eet aerial Say
Oh Ssciuge cae A wih btw (3
alte thins soca hat iat: i. “a
dep eaberins 4 ee are nen bi
op aT wll be tig May
a Re, Ae inn aval. ayer, wh wai
Diahiess MR mrss Camerata ad ston cote |
‘i “tials tii (apwniy ar tao CL a heen s see sit
hua iO ini citi ee
mi ik | i hone
y ‘Muee teh nat Stee sg ie gana cov ‘
> iy a oy Wa there Bua lial te
Dic vee rey ie re a ae susie,
ae ‘ait h
at, terest tc Re ata le
Pat | if 4 ;
Riau ith omfokg oan, gl, wy pia ai ait
a
eee « oun wy By! (yma Tonoldapyth hoasatia vi
‘ive + oe rhe ot Cras ag pce a deaideey ; Ke ate otal I: 46%)
phonon aso Ra the au ely) ke wotsa oH ‘otettig: pean enti
ay ey + ae ether ed vy, ite sey
Ww OR ee yp sali asieupeatt AaSiage, wat ‘i
“ods wah m3 a i Ra Le yi iy ot Girt she poetry) wht haa.
h iy a
oubba range aalfatamphmaga i ated, cad on tte wonton
aide
pit ine 1A : Hay ee wie a ¢ nem mt sh acihalae i he al
rade Whee hin.) Bae bina) Wi
wienegldl peasery ae MaaRe Ret edhe
ea ia
Summary
In summarizing the preceding sections, a tentative strategy can be
formulated for approaching the stochastic modelling of sea-floor
roughness.
(1) Delineations of homogeneous provinces, using a province
picker, such as that developed by Davis (1974), would insure
some degree of stationarity, and therefore validity, for the
statistics describing the area. As will be discussed in the
following chapter, this province-picking algorithm must be
based on the same statistical measure that underlies the
model, that is, the frequency spectrum.
(2) Generation of amplitude spectra from available data within the
delineated provinces would follow. Sample lengths for spectra
generation would be confined to within the province boundaries
defined in stage 1. One-dimensional spectra would describe
the roughness in several orientations to provide input for
later modelling of topographic anisotropy.
(3) Regression analyses would be performed on these amplitude
spectra to generate the coefficients of the functional form
chosen to represent the spectra. Preliminary indications are
that this form will be one or several power law relationships,
each of which require only two coefficients to describe.
(4) Anisotropy of the bottom would be modelled by studying the
variation of the coefficients a and 6 of the power law model
as a function of direction 8. In areas where complete infor-
42
(5)
mation is available in two dimensions, such as a SASS or
SEABEAM survey area, it might be possible to verify through
regression analysis a simple ® dependent function to describe
a and b. These functional forms, £,(9)=a and £,,(8)=b, could
then define the best model for estimating model coefficients
in areas where only randomly oriented track data is available.
Extension of the functional representations beyond the spatial
frequency at which the surface sonar encountered its white-
noise level, would be attempted by studying the spectra of
deep-towed sonar and bottom photographs. It is anticipated
that a general functional form, perhaps a simple extrapolation
of the power law form, will be discovered. In the many areas
where high resolution bathymetry is not available, this
extrapolated function should provide the best available esti-
mate of spatial frequency structure at very short wavelengths.
(6) Incorporation of the model into existing data sets would be
the final development stage. This roughness model would be
calculated for 5° grid cells and merged with the existing
gridded data sets being developed at the Navai Oceanographic
Office. All grid cells within a homogeneous province would be
represented by coefficients generated from data located any-
where within that province. By integrating this stochastic
model with the existing deterministic models of oceanic
depths, an essentially complete description of the sea floor
will be contained in a single data base. This data base, when
combined with similar gridded models of sea water sound veloc-
ity, sediment column sound velocity, sediment thickness, and
43
other geophysical data, would provide a comprehensive environ-
mental data base for further acoustical, oceanographic, geo-
logical, and geophysical modelling and interpretation.
The following chapters generally follow the approach presented
above, beginning with a more rigorous look at the delineation of sta-
tionary provinces. The method used for generating valid amplitude
spectra are fully described in Davis (1974), and are reviewed only
briefly. The regression techniques used in the study are described in
detail in Appendix A. The problems of anisotropy and extension of the
model into high spatial frequencies are discussed in detail in later
chapters. Throughout, interpretation of the results in terms of geo-
logical processes are presented.
44
oe ge ih ae Sila) i,
ie ZRiihae manent SS eeadegnet, v6 ti
| Scanoniil nares ie HME meee
“ae at My atic cree Thien “hy
bene anes ma, ‘Pane
Pens se 2 ny Sees ide tag
Mpc, te AN 2a ean
a A tia eo ‘eat, a at
ial pets nt cola mt ; oPaae sep
a 1 oF apne som we i + rato Ae PgR
| i. nao ibe : wi ay oe neg Re HR MR
H Povaiit ihn TD "4 ‘pevetad: oe bales vs tli
ea al il nih " we Th oui rams oan ¥ ut ie cig thee! wa vol
REY cep is. Sh wid Lota spticien: se ies! one be i
Aye
by
a Ys he gah Brine
ae
rr MR oh a a
Phe ue Rita Chad 1 ai Aw: uF i inkption
AP08), RRL iy lin MM RaRR
5. Delineation of Statistically Homogeneous Provinces
General Philosophy
The importance of defining statistically homogeneous provinces was
described in the previous chapter. The validity of frequency spectra,
and indeed nearly all statistical measures, requires a stationary sample
space. Unfortunately, truly stationary phenomena are usually only
available to theoreticians and experimentalists. The statistics of most
natural phenomena vary in either time, space, or both. It is therefore
essential in attempting to describe statistically non-stationary phenom-
ena, to delineate areas in which the statistic being generated varies
minimally and only within defined limits. In order to accomplish this
preliminary procedure of “province picking,” it is necessary to design a
detection algorithm which takes account of the phenomenon being des-
cribed and the statistic being used.
To illustrate this principle of matching the province detector to
the statistic being generated, we begin with an elementary statistic.
As an example, assume that it is necessary to describe the areal distri-
bution of height of the people of Africa. Assume for this discussion
that the desired significance of the mean requires samples of at least
10,000 individuals. One approach might be simply to divide the conti-
nent into regular square areas and randomly select 10,000 heights from
each area to generate a mean. The means so generated should indeed rep-
resent the populations of these arbitrary squares.
45
To illustrate the effect of non-stationarity on the validity of the
measured statistic, assume that in a particular sample square, the
southern region is inhabited exclusively by Pygmies, averaging only 4
feet tall, while the northern region is inhabited by Watusis, averaging
7 feet tall. The me thod just described would predict a single popula-
tion, averaging 5'6” height, inhabiting the area. In fact, very few of
the individuals in the population are near this height and our statistic
has failed to perform its intended function; to describe the areal dis-
tribution of heights of the population.
In order to confine our sampling to relatively homogeneous sample
spaces, it is necessary to detect the boundary between independent pop-
ulations prior to final sampling. The best method for accomplishing
this: “province picking” is to measure the mean crudely with much smaller
samples, say ten individuals or even one individual over smaller areas.
Even these crude measures could be sufficient to define the large gra-
dient in population mean across the boundary. Notice that the same sta-
tistical measure (the mean) is used to describe the population and to
define the province boundary. After the provinces are delineated, sam-
pling within homogeneous provinces ensures statistical validity, at
least to the degree that stationarity was confined in the province pick-
ing procedure. An additional operational advantage of the province
picking procedure is that very large areas of stationary means might be
detected which would require only one random sample of 10,000 individ-
uals to describe a large area, rather than conducting several repetitive
samplings using the arbitrary grid technique.
When one uses more advanced statistical measures to describe the
earth, it becomes necessary to design more complex procedures for homo-
46
geneous province detection. The method used by Davis (1974) was
described in Figure 4-2. In this application, the design of optimum
survey spacing for marine gravity data collection, the statistic used
was the total RMS energy in a particular spatial frequency band. The
location of this band in frequency was dictated by later applications of
the marine gravity data. Construction of a digital model of oceano-
graphic sound speed requires delineating provinces in space and time
based on statistics describing the shape of oceanographic profiles (T.M.
Davis, personal communication, 1983).
In order to delineate stationary provinces for the description of
sea-floor roughness using frequency spectra, it becomes necessary to
make a crude estimate of the amplitude spectrum discretely in the spa-
tial domain. Recall that in transforming to the frequency domain, sta-
tionarity has already been assumed, and therefore Fourler transform
techniques are not appropriate. The method used in this study takes
advantage of the relationship between band-limited energy in the spatial
and frequency domains (Parseval's Formula), and the inferred power law
form of amplitude spectra of topographic surfaces. Just as amplitude
spectra represent the amplitude of component sinusoids at discrete fre-
quencies, an equivalent estimate can be made in the spatial domain by
band-pass filtering the frequency of interest and evaluating its ampli-
tude. While the frequency domain estimate represents the least squares
average amplitude over the entire leneen of sample, the amplitude: can be
estimated discretely in the spatial domain using the Hilbert Transform.
This is very similar to the method developed by Davis (1974) for a sin-
gle frequency band.
47
In order to estimate the full spectrum, it is necessary to evaluate
the amplitude at several frequencies, spanning the range of the desired
spectrum. This is accomplished by convolving a bank of band pass fil-
ters, centered at different frequencies, with the data and evaluating
the amplitude of the band-limited signals discretely. Knowing the
“power law” functional form of the spectrum in advance, one can fit the
several amplitude versus frequency estimates at discrete points in
space, using the iterative regression technique described in Appendix A.
The regression coefficients a and b, mow available at every point along
the profile, are often highly variable and must be smoothed. Also,
because the two parameters are statistically correlated, it is prefer-
able to use the exponent of frequency (b) and total band-limited RMS as
detection parameters. Just as the presence of white noise at high fre-
quencies must not be included in the regression analysis of the ampli-
tude spectrum, amplitude estimates at the noise level in the spatial
domain are also ignored. The method is described in detail in Appendix
B, along with the results of various performance tests on signals of
known properties used to calibrate the sensitivity of the detector.
Generation of Amplitude Spectra
Having delineated statistically homogeneous segments of data on the
basis of their estimated frequency spectra, the next step is to generate
amplitude spectra from these segments. Were the spatial domain esti-
mates adequate, it would not be necessary to generate the spectra at
all. However, due largely to instabilities in the slope (b) parameter,
those estimates are not adequate and true FFI's must be run to estimate
48
the model parameters. Also, the Fourier transform method has the addi-
tional advantage of estimating the amplitude at many more frequencies
than the ten bands arbitrarily selected for the spatial domain
algorithn.
Since the proper generation of spectra is the basis for the entire
model, great care has been taken to ensure that the best estimate of the
true amplitude spectrum are obtained. The techniques used are described
in detail by Davis (1974) and will only be reviewed here. The computer
software used in this study was modified from programs provided by T.M.
Davis and is presented as Appendix D.
In using a finite length sample to represent an infinite series,
the observer has in effect multiplied the infinite series by another
infinite series consisting of zeros beyond the sample and ones at all
sample locations. The multiplication of this so-called “boxcar” func-
tion in the spatial domain, causes the true transform of the signal to
be convolved with the boxcar's transform, a sinc function, in the fre-
quency domain (see Bracewell, 1965). The presence in the frequency
domain of side lobes on the sinc function, causes energy to be “leaked”
into adjacent frequencies during convolution. Because of the red-noise
character of spectra of sea-floor topography, this “leakage” tends to
transfer energy artificially from lower to higher frequency.
Although the use of tapered windows rather than boxcar sampling
tends to reduce leakage, the preferred technique uses the method of pre-
whitening. Tapered windows have a sinc function transform with eaduee’
sidelobes and a broadened mainlobe, which reduce spectral leakage at the
expense of spectral resolution. In prewhitening, a specially designed
high-pass filter is convolved with the signal, modifying it such that
49
its spectrum appears flat, or white, rather than red. When this pre-
whitened signal is then passed through the Fourler transform, there is
no preferential transfer of energy in either frequency direction. To
obtain the “corrected” spectrum which approximates that of the true
(Anfinite) signal, the prewhitened spectrum is divided by the impulse
response of the prewhitening filter. This operation is equivalent to
deconvolving the filter in the spatial domain.
The importance of proper prewhitening can not be overstated. Leak-
age of energy into high frequencies would cause a consistent underesti-
mation of the magnitude of b (spectral slope), and degrade the ability
of the model to estimate high frequency roughness (i.e., overestimation
of amplitude at high frequencies). Figure 5-1 illustrates prewhitening
by showing a raw spectrum, prewhitened spectrum, and corrected spectrum
on one plot. Further examples can be found in Davis (1974).
Physical Interpretation of Spectral Model Parameters
Before examining the distribution of the spectral roughness model in
selected study areas, it is worthwhile to discuss the physical meaning
of the model parameters a and 6. The proportionality constant a in the
expression
a
Azae° 3?
where A = amplitude
s = spatial frequency
50
g
oh
N
Q
DEPTH (METERS)
Bg
(o7)
NORMALIZED AMPLITUDE (KILOMETERS)
2 4 6 8 10 12 14
DISTANCE (KILOMETERS)
7
ae
CORRECTED \ Wee
MH
AML
PREWHITENED
10°! 10° 10! 102
FREQUENCY (CYCLES/KILOMETER)
Figure 5-1 Illustration of the importance of prewhitening of amplitude
spectra. The raw spectrum is the result obtained by simply
performing a spectral analysis on raw data. The prewhitened
spectrum is the result of. performing the same analysis on
data which has been convolved with a special high-pass fil-
ter. The corrected spectrum results from the quotient of
the prewhitened spectrum and the frequency spectrum of the
high-pass filter, and represents an estimate of the spectrum
of the corresponding infinite signal from which the sample
data was derived.
51
represents a simple scaling factor for roughness. For a given frequency
(s) and exponent (d), the amplitude (A) is proportioned to a. Due to
the method of calculation of this expression, the actual value of a cor-
responds to the amplitude of the component sinusoid with a wavelength of
one kilometer and is usually expressed in meters or kilometers. This
particular normalization was selected because the one kilometer wave-
length falls within the sampling of most surface sonar data. For exam-
ple, the required 0.5 km sample rate would be obtained with a 1 minute
sonar ping rate on a ship traveling 30 km/hr (or 16 knots). The value
of a does not necessarily correspond to any particular features in the
signal, but only to the amplitude of the component sinusoid.
The interpretation of the exponential parameter (b) is somewhat
less intuitive. For the case of b= 0, the amplitude of all component
frequencies is constant and equal to a. This ig the well-knowm “white
noise” associated with random series such as instrument noise. Such a
value for b would customarily be interpreted as instrument noise in any
spectra from sea-floor topography. Values of b > O imply that ampli-
tudes increase at shorter wavelengths, a condition that has never been
observed in bathymetric data. What is consistently observed is the case
where b < 0, the previously mentioned “red-noise” spectrum, in which
amplitudes of component sinusoids increase with decreasing spatial fre-
quency (longer wavelength). This indicates simply that broader features
have greater height.
An interesting special case occurs when b=-1. The expression
52
becomes
or
or
A/A =a
Simply stated, for b = -1, the ratio of amplitude to wavelength (or
height to width) is constant and equal to a at all scales. This
condition was termed “self-similarity” by Mandlebrot (1982) and corre-
sponds to a fractal dimension of D= 1.5.
One might visualize the special case b = -l as a signal which
appears identical at all scales of observation. In another sense, the
signal appears equally “rough” at all scales, the magnitude of roughness
being prescribed by the magnitude of a. In cases where -l < b < 0, the
ratio of height to width tends to decrease at longer wavelengths,
although the absolute amplitude does increase. Such signals appear
rougher at high frequencies. The converse case of b < - 1, implies that
the ratio of height to width increases at longer wavelengths, and there-
fore the signal appears smoother at high frequencies. Figure 5-2 sum-
marizes these relationships.
53
LOWER RELATIVE HIGHER
SPECTRAL FREQUENCY ASPECT FREQUENCY
SLOPE COMPONENT RATIOS COMPONENT
(D) (A/d)
(a)<G)
Low High
Figure 5-2 Relationship of spectral slope parameter (b) to aspect ratio
(A/A) of sinusoid at different frequencies. For b= -1, the
aspect ratio remains constant in all frequencies. For b<-l,
the aspect ratio increases at lower frequencies. For b>-1,
the aspect ratio increases at higher frequencies.
54
Although the physical interpretation of the model parameters a and
b is clear, an important question remains as to the geological signifi-
cance of these terms. Why do certain areas of the sea floor have a par-
ticular representative spectrum, and why do all spectra seem to show
such a consistent power law form over large ranges of spatial frequency?
An obvious hypothesis is that the spectral form reflects the unique
interaction of the relief-forming processes and the materials being
affected. For example, the formation of new sea-floor crust at oceanic
ridge crests affects the relief of the new sea floor at all spatial fre-
quencies. If the relief-forming process is uniform over some geographic
region and interval of geologic time, there is no reason to suppose a
change in the statistics of the surface being constructed, although its
deterministic shape might change. Conversely, if there is a change in
the reldief-forming process (such as the spreading rate) or material
(perhaps a change in the properties of the magma source), it is likely
that the resulting relief would also be affected.
Many geological environments represent a composite of several
relief-forming processes (tectonic, sedimentary, erosional) and several
types of material. Such composite reliefs should result in an amplitude
spectrum reflecting the composite spectra of these several processes and
materials. If each style of relief is dominant over a different spatial
frequency band, and each component spectrum conforms to the power law
functional form observed in one-component cases, the composite spectrum
should appear as a set of straight line segments on a plot of log ampli-
tude versus log frequency. Examples of such composite spectra will be
shown in a later section.
55
It is difficult to prove a direct relationship between statistical
relief and combined process and material, but an interesting insight can
be gained from a study of a highly variable environment, sedimentary
microtopography. In 1981, Mark Wimbush of the University of Rhode
Island deployed stereo camera on a structure on the upper continental
rise northeast of Cape Hatteras. Stereo-pair bottom photographs were
taken at an interval of twenty-seven days. Time-lapse photography
showed that between the dates of these stereo-pair photographs, the fine
scale sea-floor relief beneath the structure was altered both by biolog-
ical activity and episodic bottom current events.
Two microrelief maps were generated from the stereo-pair images and
these are illustrated in Figures 5-3 and 5-4. Transects of heights were
taken at 0.5 cm intervals (labelled A, B, C, D) across the surface.
Amplitude spectra all showed the power law form found in spectra at
lower frequencies, and in addition the spectral parameters showed no
significant differences in spite of the gross change in the surface (see
Figure 5-5). This implies that the two surfaces merely represent two
realizations of the same statistical process. In terms of frequency
domain analysis, only the phase spectrum is altered by the redistribu-
tion of features, not the amplitude spectrum. This simple experiment
does not unequivocally prove a causal relationship between statistical
relief and process, but it does provide an encouraging result.
The Phase Spectrum
To reconstruct a profile or surface from its frequency domain repre-
sentation, it is not sufficient to model only the amplitude of each com-
56
CONTOUR INTERVAL: 2.5 cm
0 2 4 6 8 10cm
Reet |
LONGITUDE: 72°14.2°W
” Tepograghic Mep Prepered For The
UNIVERSITY OF RHODE ISLAND GRADUATE SCHOOL OF OCEANOGRAPHY
Figure 5-3 Contour representation of a bottom stereo-pair photograph
collected on August 16, 1981. Dotted lines represent tran-
sects used for spectral model generation.
57
S meri es i at
AG D
DATE OF PHOTOGRAPHY: 7- mn Se peed LATITUDE: 36°33.3'N
DEPTH OF PHOTOGRAPHY: es a cea m LONGITUDE: 72°14.2'W
CONTOUR INTERVAL: 2.5 cm 0 2 4 6 8 10cm
Fee tS |
Topographic Map Prepsred For The
UNIVERSITY OF RHODE (BLAND GRADUATE SCHOOL OF OCEANOGRAPHY
As Part Of The
GH ENERGY GENTHIC BOUNDARY LAVER EXPERIMENT
Figure 5-4 Contour representation of a bottom stereo-pair photograph
collected at the same location as that shown in Figure 5-3,
but 27 days earlier, on July 20, 1981. The spectral models
generated along indicated transects were not significantly
different from those generated from transects of Figure 5-3.
58
NORMALIZED AMPLITUDE (KILOMETERS)
10? 103 104 105 10°
SPATIAL FREQUENCY (CYCLES/KILOMETER)
Figure 5-5 Amplitude spectrum derived from profile A illustrated in
Figure 5-3. Regression lines represent model spectra from
profiles A,B,C, and D. Differences between these spectra
are within estimation error, indicating no significant
variation in time for the microtopography at this location
nor significant anisotropy within each sample.
59
ponent frequency. One must, in addition, define the position of each
component sinusoid relative to some geographic origin. The location of
each sinusoid in space is expressed by its phase relative to this geo-
graphic origin, and the composite of all component frequencies with
their corresponding phases represents the phase spectrum. Since sines
and cosines are trigonometric functions, phase is normally expressed as
an angle between - 180° and 180°.
Results from this study show that within statistically homogeneous
provinces, the amplitude spectrum can be consistently modelled with a
single or several power law functions. Although there is some varia-
bility of the measured amplitude around the simplified model, the calcu-
lated parameters remain consistent over often large geographic areas.
However, any two sample profiles are not necessarily identical or even
statistically correlated. The differences in the spatial domain man-
ifestations of identical amplitude spectra can only be due to differ-
ences in the phase spectra.
Figure 5-6 illustrates a typical phase spectrum and the statistical
distribution of its phase angles, derived from a single bathymetric pro-
file. Several profiles were examined, which represented a variety of
geographic locations and geological environments. The variability of
phase angle with increasing frequency appeared to be random. A simple
one-sample runs test was performed on several phase spectra, and all
proved to be randomly ordered to within 95% confidence limits. The runs
test is a non-parametric method (Siegle, 1956), meaning that no proba-
bility distribution of the population is assumed. Examination of the
distribution of phase angles indicates a uniform statistical distribu-
tion; that is, a random distribution in which all phase angles are
60
PHASE ANGLE (DEGREES)
0.0000 0.3689 0.7378 1.1067 1.4756 1.8445 2.2134
FREQUENCY (CYCLES/KILOMETER)
NUMBER OF OCCURRENCES
cS)
VUGCASCACASACaA aaa wACaPawcea’
GUN GRAAUACRUBAEA AAA RERABS
NNNANANAN AN AANA
CABARBUBRWSABAsSARaASaaasayl
KN NAN ANNAN ANANSI
KAN SAA NANI
BDEMEAEABARAEADSOBAMAVSAY
LAUER CAESASCAACaES
& —KSSSSSSSS5)
PAN NVNAN ANNAN SANA SNA SANS
KAR ANANANAANANA NASA SAAN
RAE MASE EMMADE RA
INSANE
TMVBABaATADBABASBBSSS
w RANNNAVLSVLVVANANNAANNSANANSSNASN
©) 7] SMAI
WALSALL AGRA SESAS
4+ Fess cn
WABAVAAAACaA CAF VCas
~ Rrvawawewvae nevus
IINSNENS NESS NENTS NN NNN ENS)
RENE NENGNENS NES NENGNENE NENG
NAN AN AASAAANSNAN ASN SAAS
ICN ANN ANAAAN ANI
ENENENTNT NS NEN NE NEN NENG]
MBABRAMAABAW
INSNENSNSNSNENENS NONE
ICNANASN ANNAN AS SANS
RANANRAN SANNA ANANS
INN ANNA AN AN ANNAN
— ([RRANRANANANSANNANARSY
8 Trae uw SNS
KRNANAN ANY
fo)
nm
3s eae
S RUT TT
8
8
8
&
3
8
8
0 60
HASE ANGLE (DEGREES)
mo)
Figure 5-6 Typical phase spectrum from a bathymetric profile collected
in the Cascadia Basin. The upper diagram plots phase angles
versus corresponding spatial frequencies. The histogram
(below) shows the distribution of phase angles by 10-degree
class intervals. Statistical testing indicates the series
to be uniformly distributed random noise.
61
equally likely to occur for any frequency component. A Kolmogorov-
Smirnov goodness-of-fit test was performed on several phase spectra, and
mone were significantly (95%) different from the uniform distribution
(Siegle, 1956).
The spatial consistency of the amplitude spectrum and the uniformly
distributed random nature of the phase spectrum of topography indicate
that the differences in bathymetric surfaces within statistically homo-
geneous provinces simply represent multiple realizations of the same
statistical process. The observed changes in the microtopography
recorded in Figures 5-3 and 5-4 can be modelled by combining two differ-
ent random phase spectra of the same distribution with the known ampli-
tude spectrum. With the functional representation of the amplitude
spectrum for an area, a “typical” profile or surface can be produced by
generating a uniformly distributed set of random numbers to represent
the phase spectrum, and performing an inverse Fourier transformation to
the space domain.
The concept of representing the sea floor as a deterministic surface
combined with stochastic variability was introduced in an earlier sec-
tion. In that discussion, the deterministic components appeared as a
smoothed, long wavelength surface which was combined with a higher fre-
quency, stochastic roughness component. By describing the higher fre-
quency components with a spectral representation, we can now visualize
the amplitude spectrum as being determined (by modelling) and the phase
spectrum as being a purely random (stochastic) process.
62
Creation of the Model and Interpretation
Having developed the algorithms for defining quasi-stationary prov-
inces and generating valid amplitude spectra, these methods were applied
to an area off the Oregon coast. The area (42°N-45°N, 130°W-124°W) was
selected due to data availability and the variety of geologic environ-
ments represented within this relatively small area. The area includes
the continental margin (shelf and slope), Astoria deep-sea fan, Tufts
Abyssal Plain, Gorda Rise spreading center, Blanco Fracture Zone, the
Cascadia Channel and numerous seamounts. All spectral estimates were
generated from data collected on the SASS multibeam sonar system by the
U.S. Naval Oceanographic Office. Only center beam depths were used in
this portion of the analysis, in order to simplify processing.
Figure 5-7 compiles the results of the combined province-picking
and spectra-generating procedures. The areas delineated by the various
shading patterns represent stationary provinces with similar ranges of
the a statistic; that is, the amplitude of the component sinusoid at a
wavelength of one kilometer. In cases of coincident values at crossing
lines, a simple average was taken, ignoring for this analysis the
effects of anisotropy. In many cases, provinces separated in the prov-
ince-picking procedure became recombined in the final spectral model,
indicating that the provincing algorithm used is more stringent than the
levels selected for presentation.
Within each of the larger provinces, the spectral slope parameter
(b) estimates were averaged and these values shown within the provinces.
The standard deviation of all estimates within any province was less
than 0.1 in all but one case shown. There is one province in which the
63
ROUGHNESS PROVINCES
Topographic Amplitude (meters) at A= km
(a) 0.0-0.5 F 15-25
1.0-1.5 [-1.50] Spectral Slope
DATA DISTRIBUTION
Figure 5-7 Distribution of roughness in the vicinity of the Gorda Rise, Northwest Pacific Ocean.
63a
slope parameter is different in two areas within a single shading pat-
tern (bottom center of chart). Many roughness province boundaries coin-
cide with obvious physiographic province boundaries. In these cases,
the bathymetry was used to trace the province boundaries between sample
spectra. In many other cases, no boundary was obvious in the bathymetry
and such tracing was not possible; that is, what appears as an abrupt
boundary may represent simply an arbitrary contour of a continuous gra-
dient. Notice also that the illustrated bathymetry from Chase et al.
(1981) was based on a totally different data set than that used for
spectral model generation, and therefore some provinces appear in the
model which are apparently unsupported independently in the bathymetry.
The gross distribution of the roughness statistic a corresponds
fairly well with what one would expect intuitively. The roughest areas
(a > 2.5 m) are located at the Gorda Rise crest, Blanco Fracture Zone
and over most seamounts, in particular the President Jackson Seamounts.
The least rough areas correspond to the sedimentary provinces of the
Tufts Abyssal Plain, Cascadia Basin, the Astoria deep-sea fan, and the
continental shelf. The Cascadia Channel appears as an intermediate
roughness province which can be traced very easily through the Blanco
Fracture Zone and onto the Tufts Abyssal Plain. The Astoria Channel,
located to the east of the Cascadia Channel, is too narrow (<8 km) for
this analysis and thus does not appear as a separate province.
Within this gross distribution of roughness, some more subtle pat-
terns can be identified. The continental margin between the shelf and
abyssal plain appears to be banded with the topography becoming gener-
ally rougher down slope. This particular continental slope represents a
slowly converging margin between the North American Plate and the Juan
65
de Fuca-Gorda Plate. Carson (1977) and Barnard (1978) have both exam-
ined this convergent margin in areas to the north off the coast of
Washington, and both present seismic cross-sections of this area. Pre-
sumably, similar processes are at work along the Oregon margin. Barnard
(1978) infers a change of compression rate from 2.3 em/year before 0.5
mybp to a present rate of 0.7 cm/year. Deformation of Cascadia Basin
sediments is progressing westward, making the deepest areas of the mar-
gin also the most recently deformed. This westward progression of
deformation process is expressed in the bathymetry as a downslope
increase in surface roughness. Barnard (1978) classifies the slope ter-
rain into an upper slope extending to a depth of about 1500 m, and an
accretionary “borderlands” complex of en-enchelon, anticlinal ridges.
Between these ridges are sediment-filled basins. These physiographic
divisions, derived by qualitative observation of the geological struc-
ture of the region, correspond closely to the roughness model generated
by quantitative methods.
The Gorda Rise is one of the more active areas of the world sea
floor, and a corresponding complexity is evident in the derived pattern
of bottom roughness provinces. Atwater and Mudie (1973) reviewed the
tectonic history of the area. Additional history of spreading rate and
spreading direction changes can be found in Elvers et al. (1973). Much
of this tectonic history may be peripheral to this study because the
sea floor affected is now buried beneath the Tufts Abyssal Plain and
Gorda Deep-Sea Fan sediments.
One interesting feature of the bottom roughness chart presented in
Figure 5-7 is the very rough ridge crest which terminates abruptly on
either flank. The full width of the feature is about 25 km. Were
ridge- forming processes constant through time and ridge crest relief
“frozen” into the topography, the same roughness would be expected to
persist, at least in long wavelengths, on older sea floor. The other
major process affecting the ridge flanks, sedimentation, would be
expected to affect the short wavelengths initially (due to their lower
amplitude), and form a smooth transition zone, not the abrupt boundary
observed in the roughness pattern. Magnetic data from the area indicate
that this zone falls within the Bruhnes-Matayama boundary and must
therefore represent crust younger than 0.7 my. Recalling that Barnard
(1978) found evidence of a slowing of compression rate on the continen-
tal margin at a time younger than 0.5 mybp in areas to the north, it is
possible that this roughness feature reflects the same change in proc-
ess, most likely a slowing of spreading rate. Future examination of
other ridge axes should reveal whether this pattern is unique to the
Gorda Rise or present under other tectonic conditions.
Perhaps equally interesting is the abrupt termination of this ridge
crest at latitude 42°20'°N. The roughness values drop (as supported by
three tracklines) two and three roughness levels at the feature's termi-
nus. This disruption of the ridge crest falls along a trend which
encompasses President Jackson Seamounts and other seamounts to the
northwest, and a major (900 m) bathymetric deep to the southeast. The
break in the ridge crest trend also appears in the bathymetric chart.
Hey (1977) developed a “propagating rift” model to describe the
plate geometrics and magnetic anomaly pattern found on the Juan de Fuca
Ridge by the Pioneer survey. According to this model, the growing
spreading center propagates along strike as the dying spreading center
becomes inactive and is added to one of the rigid plates. As this proc-
67
ess continues through geologic time, a V-pattern of fossil spreading
centers forms a “propagator wake”, as is illustrated in Figure 5-8.
This same pattern is found in the isosynchronous magnetic anomaly pat-
tern. Since this model was proposed, similar processes have been
observed on other spreading centers, in particular the Cocos-Nazca
spreading center (Searle and Hey, 1983).
A propagating ridge crest model offers one explanation for the
abrupt termination of the Gorda Rise crest shown in Figure 5-7. In
order to test this hypothesis, a magnetic anomaly map of the Gorda Rise
was constructed and is illustrated in Figure 5-9. The chart is based on
original magnetic anomaly data collected by the U.S. Naval Oceanographic
Office. The V-pattern associated with the “propagator wake” is evident,
extending to the east and nae tmaselin the ridge crest at 42°N. The
direction of the V-pattern indicates that the ridge crest to the north
is propagating toward the south at the expense of the southern portion
of the Gorda Rise crest. The geometry of the schematic model by Hey
(1977) for the. Juan de Fuca Ridge (Mgure 5-8), represents a nearly per-
fect analog to the magnetic anomaly pattern of Gorda Rise (Figure 5-9).
It would appear that the abrupt termination of the high roughness zone
shown in Mgure 5-7, is due to its association with the tip of a prop-
agating rift system.
Another feature of interest in the distribution of the a parameters
in Figure 5-7 is the existence of east-west trends of selected roughness
provinces on the ridge flank. One might have expected roughness prov-
inces to align themselves sub-parallel to the ridge strike, perhaps
reflecting changes in processes through time being felt along the length
of the ridge axis. This is the case with the very rough central valley
ways
='<
Pil
08
RS
Figure 5-8 Schematic illustration (from Hey, 1977) of the pattern of
isochronous seafloor resulting from a southward propagating
rift. Double lines represent active spreading centers,
dashed lines are fossil spreading centers, heavy line is
active transform fault, and dotted lines are associated
fracture zones. Diagonal trend lines indicate the so-called
“propagator wake”.
- .
Pp so?
ry
s*
.
42°
aie :
20° 128° W 40° 20° 1Z77° 40‘ 20° 126° 40’ 20’ 125°
Figure 5-9 Magnetic anomaly chart of the Gorda Rise area. Positive
anomalies are shown in black; negative anomalies are shown
in white. The northeast-southwest trending Gorda Rise is
obvious as is the northwest-southeast trend of the Blanco
Fracture Zone near the northern limit of the chart. Compare
the inferred “propagator wake" (indicated by dashed lines)
to the theoretical model of a propagating rift illustrated
in Figure 5-8.
70
portion. The east-west trend of roughness indicates instead relief-
forming processes which act relatively constantly through time, but
quite variably along the ridge strike. Francheteau and Ballard (1982)
describe the along-strike variability of processes along the East
Pacific Rise, and relate the changes to distance from the ridge crest/
fracture zone intersection. The change in process is expressed by the
relative importance of fluid lava flows and pillow lava flows. These
petrologic changes are in turn associated with the elevation of the rift
valley along the ridge crest; topographic highs are associated with high
ratios of fluid lavas and topographic lows associated with pillow lavas.
Indeed the shallowest portion of the Gorda Ridge segment is located near
latitude 42°45'N, which shows relatively low roughness values as one
would expect from the sheet-like flow of fluid lavas. Ridge flank areas
become rougher adjacent to deeper axial valley segments, which might
reflect the rougher surface of pillow lavas. Although other explana-
tions for the trend of roughness on the ridge flanks could be put forth,
the distributions found in this study are consistent, although not nec-
essarily typical. Extension of the model to more thoroughly investi-
gated ridge crests should shed light on this particular hypothesis.
The patterns apparent in the distribution of a are not as evident
when one examines the distribution of the spectral slope (b). It is
clear that the “universal” value of -1 for the spectral slope inferred
by Bell (1975b) and others is not supported by this modelling effort.
The reason for this discrepancy is not readily apparent, however the
attention given to defining stationary sample space in this study does
represent one major difference in method. In order to test this hypoth-
71
esis, amplitude spectra were generated for several profiles in the Gorda
Rise area which spanned multiple stationary provinces. The effect of
these extensions was to reduce the degree of statistical homogeneity of
each profile. The computed spectral slope (b) parameters for the
resulting amplitude spectra converge consistently to the value b = -1.0
as longer profiles are tested. Figure 5-10 illustrates one such long
profile extending nearly 500 km from the Oregon coast. The profiles
analyzed by Bell (1975b) were often much longer. No formal statistical
argument for this convergence to b = -1.0 will be attempted here, how-
ever, it would appear that the concatenation of multiple profiles of
differing spectral characteristics results in a profile which resembles
a random walk. The need to define statistically homogeneous sample
spaces before generating statistics is clearly demonstrated.
Most slope values shown in Figure 5-7 cluster about -1.5 with the
exception of the smooth ridge axis segment south of 42°20'N. The two
large sedimentary provinces are represented by two values for spectral
slope. Figure 5-11 illustrates a typical amplitude spectrum from these
sedimentary provinces. The spectrum clearly separates into two distinct
straight-line segments of different slope. The average values of these
distinct slopes for all profiles is given in the corresponding boxes
(see Figure 5-7). If the hypothesis is accepted that the characteristic
spectral slope of amplitude spectra of sea-floor topography represents a
dominant relief-forming process, these spectra should represent areas
where two processes are at work, affecting the relief in different spa-
tial frequency bands. It is likely that the higher frequency band rang-
ing from A = 2.5 km and with 6 = -.6, represents the sedimentary regime
in these areas. The lower frequency process with b < - 1.4, 1s less
72
DEPTH (METERS)
a :
8
te) 100 200 300 400
DISTANCE (KILOMETERS)
SLOPE =-1.813
INTERCEPT >.418-883
10°!
10-2
104
NORMALIZED AMPLITUDE (KILOMETERS)
10-5
10°
10-3 10-2 10°! 10° 10!
FREQUENCY (CYCLES/KILOMETER)
Figure 5-10 Amplitude spectrum of a long bathymetric profile, which
encompasses several statistically homogeneous provinces.
The inclusion of non-stationary segments into the analyses
results in a spectral slope parameter which approaches b=
73
DEPTH (METERS)
= —_ _ —_ —_
(-} (=) [~) (>) (>)
b & > ou, °
NORMALIZED AMPLITUDE (KILOMETERS)
r=)
bs
Figure 5-11
8 16 24 32 4 48 56
DISTANCE (KILOMETERS)
10°! 10° 10! 10 2
FREQUENCY (CYCLES/KILOMETER)
Typical amplitude spectrum of profiles collected on the
Tufts Abyssal Plain. Spectrum shows two distinct
straight-line segments which intersect at a spatial
frequency of ~0.4 cycles/km. The longer wavelength
portion is described by = -1.847, a = 0.0294 m. The
shorter wavelength segment is described by 6 = -0.541, 4 =
0.113 m. Profiles from the Cascadia Basin province have
similar spectra, however the model parameters are slightly
different.
74
obvious, but may well represent an underlying tectonic effect due per-
haps' to the compressional nature of the area overprinting a long wave-
length relief on the smooth sediment.
This general correspondence of b < -1 in tectonic provinces, and -1l
< b < 0 for sedimentary provinces is found in many areas outside this
study area. The reason for this general relationship can only be specu-
lated upon. Recall from the previous section that b < -1 requires that
the ratio of height to width of component features, increase at longer
wavelengths. If one envisions tectonic relief forming processes, the
entire morphology is constructed and then erosional and sedimentary
processes begin smoothing small features first, progressing to con-
stantly larger scales. In sedimentary processes, one can envision a
smooth layer of sediment being affected by bottom current interaction
for example. This constructive process begins at the highest spatial
frequencies and progresses to larger scales. This is in agreement with
the relationship of -1 < b< 0, in which the ratio of height to width of
component features increases at shorter wavelengths.
The study area illustrated in Figure 5-7 was selected because of
its rich geological diversity. As such, the distribution of roughness
provinces is correspondingly complex. To allow a comparison with a more
tectonically stable area of the world ocean floor, a large number of
tracklines off the United States east coast were analyzed to compare a
passive continental margin. Although the trackline density was not ade-
quate for a complete chart to be drawn, one interesting relationship was
discovered. An extremely large area of the continental margin, compris-
ing most of the continental rise, was found to have in common a very
distinct amplitude spectrum. Figure 5-12 illustrates a typical ampli-
75
DEPTH (METERS)
) 16 32 48 64 80 96 112, 128
DISTANCE (KILOMETERS)
10-1
10-2
10-3
NORMALIZED AMPLITUDE (KILOMETERS)
10-2 10°! 10° 10! 102
FREQUENCY (CYCLES/KILOMETER)
Figure 5-12 Typical amplitude spectrum of profiles collected on the
Continental Rise, east coast U.S. All fourteen profiles
examined in this province (shown in Figure 5-13) show
nearly identical patterns. The two linear segments, in
this example intersect at ~0.3 cycles/km., with 6 =
-1.813, 4 = 0.0142 m. for the longer wavelength model, and
= -.524, a = 0.0645 m. for the shorter wavelength
model.
tude. spectrum from this large province, and shows the same two-process
mature of the spectrum in Figure 5-ll. Im fact, this spectrum is almost
identical to the measured spectra from the Tufts Abyssal Plain and the
Cascadia Basin.
Figure 5-13 shows the location of profiles with this characteristic
spectrum. Some of these profiles are over 200 nom long, and the profiles
are distributed over an area more than 1000 nm in extent. The total
standard deviation for the spectral parameters in all fourteen profiles
was only 0.15 for the slope (b) parameters and 0.02 meters for the
intercept (a) parameters. These values apply to both the lower fre-
quency and higher frequency line segments. The area may well extend
even further northeast or southwest. Profiles are only broken by sea-
mounts or deep=sea channels associated with submarine canyons. As such,
one can see that the areal extent of a roughness province with a given
level of allowed non-stationarity, can wary, econ hundreds of thousands
of square miles to the slopes of a single seamount.
77
"p-G pue E-G saunbyy uy paqeugsni i} sydeshoyoyd w07z30q B43 yO uO,QeD0] aYyQ syseW MH ceaUe
- 3FIGGSH FY JO UOLQEIO, JYyQZ Sysew H °394450 BZYydeubouR|aZN [eAeN |YyQ Aq padnpoud (gqgga) eseg
beg Bsqzowmyqeg 12346;q pappys6 aqnujw-9Ayy ayQ wos, (S4dq9W Uy) PaiNOqUOD A{ 12d; 3ewWO Ne
S} Asqemyyeg
*RuZIads apnzy (due [ed;quapy A_ueau Buymoys sajyyoud I, ujzaw{yjeq yO u0y3Qe907 ET-G aunby J
78
6. Anisotropy of Surfaces
The necessity of directionally treating anisotropic surfaces was
introduced in Chapter 4. Any statistic generated from a one-dimensional
profile of a two-dimensional surface is only valid in all directions if
that surface is isotropic. This chapter examines the importance of
anisotropy in detail. A simplified theoretical model of the effect of
anisotropy on the frequency spectra of directionally sampled profiles is
formulated and tested. Such spectra are generated for two areas of the
sea floor where complete, two-dimensional bathymetric data are available
from multibeam sonar. An identical study is performed on data from
side-sean sonar. These results are then compared to the theoretical
model. Next, the possibility of estimating such two-dimensional func-
tions from randomly oriented bathymetric profile data is discussed.
Theoretical Model
Before examining the effect of anisotropy on measured frequency
spectra derived from actual bathymetric profiles, it is instructive to
examine a very simple theoretical model of this effect. Several of the
concepts introduced in Chapter 4 will be utilized. There it was shown
that the effect of one-dimensionally sampling a sinusoid which has been
extended to two dimensions (see Figure 4-6), is to stretch the true
wavelength as,
AY = |cos™!6| Oi
79
where A" = apparent wavelength
4 = true wavelength
@ = angle of sampling (0° = perpendicular
to linear trend)
This relationship can then be combined with the similarity theorem
of Fourier Transforms to yield the transform pair
£ (|cos 6] * x) D> |cose|~! + F (s/cose)
Recall that these relationships were formulated for a single component
wave form extended to two dimensions. Actual sea-floor topography has a
spectrum which {fs continuous and conforms to a power law functional form
Amzae gb
To extend the model to the continuous case, one can envision gener-
ating a topographic profile (or other signal with continuous power law
form, such as a random walk model), and extending all points on the pro-
file to the second dimension. This surface is then sampled in various
directions and the spectrum of each profile evaluated for a and b above.
Combining the above relationships yields
£(|cos 0| * x) D |eos e|7! + a = (s/cos 0)”
or equivalently
£ (|cos 0|~! = x) D |cos 8| + a+ (s * cos @)>
or
F (s, 0) = |cos 0| + a * (s * cos 6)?
where F(s, 8) is the Fourier transform of £(|cose|~1 ° x) expressed
in terms of 6 and s. The coefficient a can now be expressed as a func-
tion of 6 as
a(6) = |cos0| * a
which can be peneralized to
a(@) = |cos(6-8,)| * a
where 6, {s the azimuth perpendicular to the linear trend. Notice
that for © = 6, that is a profile generated perpendicular to trend,
|cos(0-6,)| = 1 and F(s,0) = a * 8, the original function. For 6=0, =
£90°, |cos( - 94) | = 0° and F(s,9) = 0, that is, for profiles sampled
parallel to strike, the series is a constant (normalized to zero) and
therefore contains no energy. The parameter b is independent of 0.
Figure 6-1 illustrates the hypothetical “plane wave” surface that is
used in this simple model. As a test of the above theory, radial tran-
sects of the surface were generated and the resulting profiles input to
the standard Fourier transform routines used throughout this study.
Figure 6-2 depicts the sampling patterns used to generate the profiles.
81
Depth
Figure 6-1
Randomly Generated Plane Wave
ma HU
I)
Mf; ma ‘ah
Miff f
My)
rn SNe
vin SS Mi Mi
at He ‘yey
WT HU
64.0
Longitude
Graphic representation of an elementary theoretical model
for anisotropic surfaces. The surface is created by gener-
es a ue oi walk (with theoretical amplitude spectrum of
= ) in the x (90° azimuth) direction. These values
are met extended to the second dimension, creating a
lineated surface with strike = 0°. Viewpoint is from the
southwest.
82
37 TRANSECTS
EVERY 5°
Figure 6-2 Radial sampling pattern used to generate one-dimensional
amplitude spectra from various azimuths on a surface.
83
Figure 6-3 plots the value of a (“intercept”) and b ("slope of spec-
trum”) as a function of azimuth. The resulting spectral estimates are
fitted with the above model. Approaching azimuths of 0° and 180°, the
profiles nearly parallel the linear trend and the results are unreli-
able, due to the small number of depths available for analysis resulting
in a very narrow frequency band width used in the model regression.
The coefficients b(6) do not show any relationship to the trend, as
predicted by theory. This does not imply that directional dependence in
b(6) is never found in spectra from sea-floor profiles. If a surface
had two or more distinct signals (perhaps due to differing relief-
forming processes) superimposed, cyclical behavior in b(@) would be pos-
sible. For example, envision a hypothetical surface composed of an iso-
tropic two-dimensional signal with spectra slope b,(@), overlain by a
simple linear trend (like the one described above) with spectral slope
bg(®). Perpendicular to trend, b(8) would be some combination of b,(6)
and bg(®) depending on their relative amplitudes. Parallel to trend,
b(@) would equal just b, (8), in this example, since the linear trend
with slope bg(9) is constant in the direction parallel to strike. The
example spectra from the Mendocino Fracture Zone presented in Chapter 4
illustrate this effect. Appendix E examines a series of artificially
generated surfaces and their spectral characteristics.
In the results shown in Figure 6-3, the parameter a(@) shows the
expected relationship to cos(6-6) 0, = 90°, as predicted by theory.
This model cosine function can be used to parameterize the surface
anisotropy. The model sinusoid is generated by an iterative regression
technique (see Appendix C.2) which determines the following equation
84
= lm o
°o ° °
INTERCEPT (< )
(o}
ce)
-1.0
Figure 6-3
RANDOMLY GENERATED PLANE WAVE
& o
SLOPE OF SPECTRUM (b )
‘
i)
20 40 60 80 100 120 140 160 180
AZIMUTH (degrees)
Distribution of spectral parameters versus azimuth of samp-
ling for theoretical surface shown in Figure 6-1. The
upper series represents the slope of the spectrum in log-
log space and varies randomly around the theoretical value
of -1 as predicted by theory. Near 0° and 180°, the esti-
mates become unstable due to poor sampling. The lower ser-
jes represents the intercept of the spectrum in log-log
space (or the coefficient of frequency) and agrees well
with the sinusoid model predicted by theory. Notice the
maximum intercept, and therefore total spectral energy,
corresponds to a sample taken perpendicular to strike
(azimuth = 90°).
85
a(9) =u +v° cos(2 ° (9-8,))
the three term regression technique yields estimates for u, v, and Qo.
These terms have definite physical meaning. u represents the simple
mean roughness level of the surface, that is the mean a(@) of the signal
sampled in all directions. This could be visualized as the “isotropic”
component of the surface. v determines the amplitude of the sinusoidal
component of the regression model and represents a measure of the degree
of anisotropy of the surface. 80 estimates the normal to the true azi-
muth of the linear trend. Frequency is not estimated since the perio-
dicity of 1 cycle/180° is known.
Unfortunately, it is not possible to decompose more than one linear
trend in a surface using this method. Envision a surface consisting of
two linear trends of differing orientation (8 a and 6g), “anisotropy”
levels (va and vg), and “isotropy” levels (ug and ug)- The surface,
being a simple linear combination of the two component trends can be
expressed as
a(®@) = (ug + ug) + va © cos(2 * (6-0,)) + vp ° cos(2 * 8&-6g))
In this example, u, and ug are both presumably equal to zero for “per-
fect” linear trends. However, even in non-perfect cases in which some
energies are available parallel to strike, the u components are linearly
combined and can not be differentiated. The “anisotropy” components
also combine linearly to yield another sinusoid whose amplitude and
phase are dependent on the relative amplitudes and phases of the orig-
inal sinusoidal components. Appendix C presents a geometric proof of
this relationship.
In the case where two or more linear trends are present in a sur-
face, the form of a(@) will show a simple sinusoid with phase and ampli-
tude which can not be uniquely decomposed into component sinusoids. If
an estimate of the azimuth and amplitude of one of the trends could be
produced independently, this component could be removed and the remain-
ing component sinusoids analyzed. It should be emphasized that the
inability to decompose component trends in no way invalidates the model,
it simply complicates the interpretation of the model statistics in
terms of formation processes.
Multiple linear trends can only be decomposed in cases where the
trends are sufficiently band-limited to appear as distinct peaks in the
frequency spectrum. This approach allows decomposed trends to be
uniquely identified by spectra generated in two orthogonal directions,
as was shown by Hayes and Conolly (1972). Their work clearly shows such
trends in the large scale topography of the Antarctic-Australian
Discordance. In examining a great many spectra of small-scale (i <8 km)
topography during the present study, no significant spectral peaks have
been observed, even in topography which is highly lineated at the larger
scales. This result may be due in part to the presentation of these
spectra in log-log form.
Comparison of Theoretical Functional Forms with Multibeam Sonar Data
The previous section developed a simple theoretical model of the
effect of linear trends on frequency spectra from profiles sampled at
87
varying azimuths on an anisotropic surface. Figures 4-7-10 illustrated
the effect of linear features on two bathymetric profiles collected at
near right angles. To test the validity of the theoretical model to
actual bathymetry, ie is necessary to sample regularly around the con-
pass at a single location on the sea floor. This is only possible for
areas which have complete areal bathymetric coverage of high spatial
resolution.
The most practical instruments available for obtaining such data
are the multibeam sonar systems. These systems, which provide a “swath”
of discrete soundings on a line perpendicular to the ship's track, allow
complete coverage of an area. By conducting surveys in which the paral-
lel survey tracks are spaced so that the outer beams of adjacent tracks
overlap or are nearly juxtaposed, a complete two-dimensional survey can
be performed. Very few of such data sets are currently available. How-
ever, two data sets from contrasting geologic environments were made
available for this study (see Chapter 5, Section B). Both surveys were
conducted by the U.S. Naval Oceanographic Office using the SASS multi-
beam sonar system (Glenn, 1970). The recent acquisition of the academic
SEABEAM systems should provide more full-coverage surveys in the future.
New techniques for processing side-scan sonar data from the SEAMARC-1
system allow similar two-dimensional analyses to be performed at smaller
scales.
Figure 6-4 presents the contoured bathymetry from the Gorda Rise
area of the northeast Pacific Ocean. Contours were created automati-
cally and only appear where supported by multibeam soundings. The
coverage is generally complete with the exception of a small area on the
western margin of the chart and small gaps between swaths. This chart
¥ == .
ag Pe, SO
DW BN
=
== O
= R=
= SS ~ a
~ ‘~=
WA\NAY ° =
ANY
NEN
y™\QG58)
a SS ES
- « Ie)
2 > Gs)
rs
& SSe7 i)
zi SEN a7 y
RSS WN
Ge .
) iw
° J i
\
) f J
i ;
‘ TVW xa,
A | rf
WY e- ||
1
1
Laer
ZG
Z
Ta
GY
°o
=
Ew
Tet
(ES;
LY
4 i vA ‘i 4, fi
S) Ai e A
Figure 6-4 Index maps showing the location of two-dimensional SASS
bathymetry data used in study of azimuthal dependence of
topographic spectra. The data are in the vicinity of the
crest of the Gorda Rise in highly lineated topography.
Central coordinates are 42°53.5'N, 126°40'W. Contour
interval is 10 fathoms.
89
was created using a grid spacing of 0.25 minutes of longitude and lati-
tude (~460m x ~300m), well above the resolving capability of the SASS
system.
The area shown in Figure 6-4 was selected to test the effect of
anisotropy on one-dimensional amplitude spectra. The data set is
located in the vicinity of the ridge crest, which was identified as a
quasi-stationary province by the methods described in Chapter 5. The
lineation of the topography is obvious on the index chart and trends
approximately N25°E.
Figure 6-5 represents graphically the test area. Although the grid
spacing used for the illustration is 0.1 minutes of latitude and longi-
tude, the profiles were generated from a grid with spacing of 0.05 min-
utes (~100m), which approaches the resolving limit of the SASS system
for these water.depths and noise level conditions. Again the sampling
pattern shown in Figure 6-2 was used to generate radial profiles of 256
points. No ensemble averaging was used.
The spectral parameters a(®) and b(6) are plotted versus azimuth in
Figure 6-6. The results agree closely with the theoretical model devel-
oped in the previous section. Notice first that the parameter b(6),
plotted above in these diagrams and labelled “Slope of Spectrum", shows
no systematic fluctuation with azimuth, as predicted by theory. Note
also that the mean slope is -1.24, well below the -1.0 slope of the ran-
dom walk (Markov process) model. The large variability of this paran-
eter could be reduced by ensemble averaging of several spectral esti-
mates, created by offsetting the center of the sampling pattern. The RMS
variability would be reduced by 1/ Y N , where N is the number of esti-
mates (see Chapter 7).
SASS Bathymetry — Gorda Rise
Figure 6-5 Graphic representation of SASS bathymetry data projected
onto an evenly-spaced grid. Illustration uses 128 x 128
points spaced at 0.1 minutes of latitude and longitude
without cartographic projection. Fourier analysis was per-
formed on 256 points from a 0.05 minute grid. Measured
strike of the lineations is approximately 25° azimuth.
Viewpoint is from the northeast.
91
INTERCEPT (¢)
S
°
&
a)
SASS BATHYMETRY—GORDA RISE
SLOPE OF SPECTRUM (> )
-2
40 60 80 100 120 140 160 180
AZIMUTH (degrees)
Figure 6-6 Distribution of spectral parameters versus azimuth for
Gorda Rise spreading center bathymetry shown in Figure 6-4
and 6-5. Spectral slope parameter (above) shows no appar-
ent functional relationship to azimith as predicted by
theory (notice that the mean slope is -1.24, well below
-1.0, which would correspond to a Markov process). The
intercept parameter (below) clearly shows the effect of
seafloor anisotropy and generally conforms to the sinu-
soidal model. Model parameters are as follows: mean ampli-
tude = 1.68 m, amplitude of sinusoid = 0.51 m, azimuth of
maximum energy = 115° (perpendicular to observed strike).
92
The results for a(8) also conform reasonably well to the model pre-
diction. The values, plotted below and labelled “Intercept”, represent
the amplitude (in meters) of the Fourier component at a wavelength of 1
kilometer. The mean intercept is 1.68 meters, indicating a relatively
rough topography, and corresponding to the “isotropic” term u in the
model. The anisotropy of the surface is obvious from the large ampli-
tude (0.51 meters) of the model sinusoid (the v term). The maximum
value of the model occurs at an azimuth of a, = 115°. This corresponds
to the normal to the linear trend (which was measured as 25° in the
full-coverage chart) as expected. These values fully parameterize the
effect of surface anisotropy on the frequency domain description.
Although the functional models for b(®) and a(®) appear to be of
the proper form, there are some obvious variations in the measured par-
ameters from the model. Im order to give an intuitive impression for
the degree that the generalized model departs from the actual measured
spectrum, Figure 6-7 illustrates a “worst case” example in which the
modelled a(@) differs from the observed a(@) by ~0.5 m. The selected
profile is from 9 = 115° (see Figure 6-6). Plotted with the measured
spectrum is the regression line derived from the functional model,
rather than the least-square fit usually shown. The model-derived spec-
trum provides an excellent representation of the spectrum and appears to
fall well within the estimation noise of the spectrum, even in this
“worst case” example.
The second study area is shown in Figure 6-8. This data set repre-
sents a totally different style of sea-floor topography due to the geo-
logic environment, which is controlled by sedimentological processes
rather than the tectonic setting of the Gorda Rise. The broad trend of
93
-250
DEPTH (METERS)
oe
ae 40 60 80 100 120 140 160 180
DISTANCE (KILOMETERS)
10!
SLOPE = -1.240
INTERCEPT = 2.190 METERS
10° AZIMUTH = 115.0°
10°!
10-2
10°3
10-4
NORMALIZED AMPLITUDE (KILOMETERS)
10-5
10°
10-2 10°! 10° 10! 102
FREQUENCY (CYCLES/KILOMETER)
Figure 6-7 Example amplitude spectrum from the Gorda Rise crest sampled
at azimuth 115 degrees. The straight line segment through
the spectrum represents the model prediction, rather than a
least squares fitted line. As seen in Figure 6-6, this
particular azimuth represents a large deviation of the
fitted parameters from the model sinusoid. The model line
still appears to fall well within estimation noise, even for
this “worst case” example.
94
Y iN SS : J
y ~>
z = Dy: 7 <4 Mis ) W/Z,
LOX SN ‘ S Ne :
\ S ih NY ap a *
DS Ver = PH es (SD \ Pp
iS DD NOTA aC Apex’
HK OS Re TS ¢ <\>}
Bo U = x we AS BAO
C7. \Wo@ mS aay (\ \
y.. XFAN, {ee , ARON
~. < eZ Gi 7 tf
ee hs a; 4 i a)
yk < SNRs ) ~ 5
y K ‘.. XQ D U =
Y Ea SF 7 eG
N
y,
JNNWAS
We es
WEA
wy was
ane
De
Pi ad
Zp
+ E
Wl Ae
LEY
71°50
Ws
ENP hs
:
=
(Ge
—
‘G
SARS
PT ESSeA
INSSE WS
yp
DSW)
¥
=z
D
AS
»
=
SJ
Figure 6-8
Me 1 ON
nort Ww
< Vp
Nolte
Org &
er unr
—-o el wv
BmpOgk e
SCevgoEe
Cow wS°9
Tak ¢ UL
Law »
VUE_5-— 9
2 og ot
seas
oo Ppa
ise oo
ay es
PEces -
fo) ow
yo a@>
rs) post
(rs cj 6 WY
o 8 oe
(<j zZzrte
[e) C= = fe a
H- ALP — NV
wyTaqw .
eo 7000s
bd ie)
o
oa ee
=
cF°os
cS
2 ee
e =
SoD °
n= Pez
oseVo-
= o eg
HEgke O
2s ae
sy3O8eN
£U8, oF
“v eee .
Son, © =z
-_ =
ot 6 oe
eee) von
© >o
& “FT -
x a>aP - LPM
ers) 3
UPHPU—-— ww
Eemw~asoat
HOANPL
95
slope toward the south-southeast is very long wavelength and not
included in the analysis. The small channels traversing the slope are
due presumably to downslope transport of sediment and represent an
apparent linear trend with wavelength of approximately 10-15 kilometers.
This is at the low frequency limit of the present analysis, but might
indicate such a trend in shorter wavelengths. Higure 6-9 graphically
illustrates the data set (in this case gridded at 0.1 minutes of lat-
{tude and longitude), and shows this linear trend due to down-slope
processes. Coincidentally, the survey track was run quasi-parallel to
this trend which could further complicate interpretation.
The distribution of spectral parameters with azimuth are plotted in
Figure 6-10, in the same format as the previous plots. Notice again
that the parameter b(8) plotted above shows no functional relationship
to azimuth 8, as predicted by theory. In this case, the mean value of
the b(®8)'s is -1.05. The intercept parameter a(@) reflects the rel-
atively smooth, isotropic nature of this sample of the sea floor. In
this case, the mean amplitude of 0.097 meters is less than 62 of the
same parameter in the Gorda Rise area. The “anisotropy” term, vin this
ease is estimated at 0.0063 meters, only 1% of the value for the Gorda
Rise. With these almost isotropic conditions, the estimate of the azi-
muth of maximum energy cannot be made with any fidelity.
The final test area represents an intermediate level of both gen-
eral roughness and degree of anisotropy. The data were collected by the
SEAMARC-1 side-scan sonar system, which is a deep-towed instrument
developed by W.B.F. Ryan of Lamont=Doherty Geological Observatory. The
vehicle is towed at approximately 500 m above the sea floor and collects
data from side-scan sonar and from down-looking sonar. The depth of the
ou SASS Bathymetry — Continental Rise
1500 “ee
Depth in Fathoms
i)
q U
wecee le
'
t}
I
Hs
el
si
a
4
aH
a ¢
SS
Figure 6-9 Graphic representation of SASS bathymetry data projected
onto a 128 x 128 point grid. Grid spacing is 0.1 minutes
of latitude and longitude and is presented without cartog-
raphic projection. Visible in the surface are spikes assoc-
jated with sonar processing noise and smooth areas where
surface was interpolated between tracks. Viewpoint is from
the southwest.
97
NN.
fo)
°o
INTERCEPT (¢)
2,
>)
°.
Po)
SASS BATHYMETRY—CONTINENTAL RISE
SLOPE OF SPECTRUM (b )
40 60 80 100 120 140 160 180
AZIMUTH ( degrees)
Figure 6-10 Distribution of spectral parameters versus azimuth for con-
tinental rise bathymetry shown in Figures 6-8 and 6-9.
Spectral slope (above) shows no apparent functional rela-
tionship to azimuth and has a mean slope of -1.05. The
intercept parameter (below) reflects the relatively smooth,
isotropic nature of the surface. Model parameters are as
follows: mean amplitude = 0.097 m, amplitude of sinuscid =
0.0063 m, azimuth of maximum energy = 120°. Due to the low
ive) of anisotropy, the azimuth direction is not signif-
cant.
vehicle is continuously measured. Farre and Ryan (1984) have developed
a method of combining these sources of information into a detailed con-
tour chart of bathymetry. These contours, when evaluated for depth on
an evenly spaced grid, comprise the data set used in this study.
Figure 6-11 illustrates a contour chart of the study area, which
encompasses the Carteret Canyon, continental slope, and upper continen-
tal rise. The smaller area outlined was used to produce the spectral
parameters illustrated in Figure 6-12. The data grid uses a spacing of
only 5 meters and represents much higher resolution than the surface-
derived sonar data in the other two examples. Unfortunately, due to
round-off errors producing a white-noise level of 0.3 meters in the data
(only whole meter depths were retained), most spectra were only above
noise to the 100 m wavelength band, similar to the resolution of SASS.
Figure 6-12 illustrates the results of the azimuthal dependence of
spectra study. The slope parameter (b) shows a mean of -1.88, lower
than the other two study areas. The data also indicate what might be an
example of azimuthal dependence of b, such as that observed on the
Mendocino Fracture Zone and discussed in Chapter 4 and Appendix E. The
intercept (a) parameters are u = .87 Mm, v = .23 m and relative azimuth
(8,) = 85°, which represent a level of roughness and anisotropy inter-
mediate to the previous two examples. Because the original survey was
collected at an azimuth of 140°, the true azimuth of anisotropy is 45°.
In conclusion, the simple theoretical model of the effect of linear
trends on frequency domain statistics, can adequately describe anisot-
ropy in three test areas. The parameter b(6) appears to be independent
of azimuth in at least two areas, although quite noisy. The sinusoidal
construction of the intercept parameter as
99
“I LILYA
Jeuos ayy yo yaed ayy aqedjpuy saup, Aaway °(HeET) UeAY
pue assejy worg St e3eG “paus{ano sy Apnas yynwyze voy pasn
UO}3IaS 6 Huyssad0ud weuos UeIS-aPLS T-IUVWIWSS WOU PaAjsap
vaue uokue) yavaquej/adoys equaujquoo aux yo Auqawkyyeg TI-9 eunby4
} (
\
.
a\\
wane
\ \\y
YE y
i WB
of — Ye y
i) Z A
Ip,
Ni)
yl
Y
a
100
INTERCEPT (0)
° = rX)
o o o
5
Figure 6-12
SEAMARC-1—CONTINENTAL SLOPE
- To)
wo
SLOPE OF SPECTRUM (> )
AZIMUTH (degrees)
Distribution of spectral parameters versus azimuth for con-
tinental slope/submarine canyon bathymetry shown in Figure
6-11. Spectral slope (above) has a mean value of -1.88,
although it may also contain a cyclical component not noted
in other examples. The intercept parameters (below) repre-
sent an intermediate level of both mean roughness and
anisotropy between those of the Gorda Rise and Continental
Rise examples. Model parameters are as follows: mean
amplitude = 0.87 m, amplitude of sinusoid = 0.23 m, azimuth
of maximum energy = 85°. Due to the direction of sampling,
true azimuth on the earth corresponds to 45°.
101
a(@) =u +v¥° cos(2 * (0-89))
allows the background, or isotropic, roughness, as well as the degree of
anisotropy and its trend to be quantified, and therefore compared in
different areas.
If one assumes this simple model of anisotropy to be true, at least
in some cases, an interesting insight into the effect of scale on
anisotropy can be seen in the mathematics. The assumption of a constant
value of b(@) at all azimuths can be envisioned as a family of lines in
log-log space of constant slope whose levels vary with azimuth. The
orthogonal azimuths which represent the extremes of anisotropy, would
maintain a constant spacing in amplitude at all frequencies in log-log
Space. That is, if at a given frequency the amplitude in one direction
were twice that of the normal azimuth, that relationship would remain
constant at all frequencies.
In examining bathymetry and other geological data, it often appears
that anisotropy decreases at shorter wavelengths. Bell (1975) reached
that conclusion in studying the aspect ratio of shapes of seamounts of
different sizes. Much of the validity of this statement depends upon
how anisotropy is defined. As stated above, if a relative doubling of
amplitude in the orthogonal direction occurs at one scale, this same
doubling should occur at all scales. However, if one considers the
absolute difference in amplitudes in orthogonal directions, at long
wavelengths the difference between perhaps one and two meters of ampli-
tude represents one meter of difference, while at very short wavelengths
this same relationship might appear as the difference between one and
two centimeters. Although the proportional relationship of amplitude
102
(doubling) is constant, the absolute difference decreases exponentially,
depending upon the value of b. Since the ability to resolve amplitude
is limited, the ability to resolve anisotropy at small scales is also
limited and could lead to an erroneous conclusion concerning anisotropy
at high frequencies. These relationships do not apply in cases where b
varies regularly with azimuth, such as those discussed in Appendix E.
Estimation of Two-Dimensional Spectra from Randomly
Oriented Ship Track
An alternative method for describing the topographic roughness of a
surface over all azimuths is with the two-dimensional amplitude spec-
trum. The method involved is quite similar to that used in the one-
dimensional case, but prewhitening requires a circularly symmetric high-
pass filter, and a two-dimensional Fast Fourier Transform algorithm is
used. Perhaps most important to the practical use of this method is the
requirement for a complete two-dimensional array of data (depths) as
input.
As mentioned previously, complete areal bathymetric surveys are
available in very few areas of the world ocean. To be practical, such
surveys must use a multibeam sonar array such as SASS or SEABEAM, and
tracks must be spaced so that adjacent swaths are juxtaposed. Of the
areas presented in the previous section, the Gorda Rise survey shown in
Figure 6-4 indicated the highest degree of anisotropy and was therefore
selected for two-dimensional FFT analysis.
Figure 6-13 illustrates the results of generating a two-dimensional
amplitude spectrum from the gridded bathymetric data shown in Figure
103
TWO-DIMENSIONAL AMPLITUDE SPECTRUM
2.57122
1.71415
0.95707
0.95707
FREQUENCY (CYCLES/KILOMETER)
fo)
-1.71415
-3.5157 -2.3438 -1.1719 0.0000 1.1719 2.3438 3.5157
FREQUENCY (CYCLES/KILOMETER)
Figure 6-13 Two-dimensional amplitude spectrum from the Gorda Rise area
illustrated in Figure 6-5. Log-transformed amplitude esti-
mates, computed via a two-dimensional Fourier transform,
are represented by light contours drawn every 0.5 order of
magnitude. The heavy lines represent amplitude estimates
predicted by the four-parameter, azimuthally dependent
model derived for the area.
104
6-5. Amplitude estimates appear as irregular contours. The amplitude
values were log transformed before plotting, and these values are plot-
ted at integer increments. Contours are drawn each .5 order of magni-
tude of amplitude. The data set, and therefore its amplitude spectrum,
is oriented with the columns parallel to longitude and rows parallel to
latitude.
Plotted with the spectrum in heavy lines is the two-dimensional
spectrum as modelled from one-dimensional profiles by the method
described in the previous section. Because the gridded data base used
in the analysis was spaced evenly in latitude and longitude, the spec-
trum as a function of spatial frequency is necessarily distorted. High
frequency noise associated with the east-west oriented track lines
appears as a smearing of the contours in the vertical and horizontal.
The simple model spectrum expisines most of the variance in ampli-
tude. The lineation of the topography with a strike of 6 = 25° can be
easily seen as an elongation of the contours in the cross-strike (@ =
115°) direction. The degree of anisotropy (v) term in the model deter-
mines the elongation of the contours. The model appears to overestimate
the amplitudes in the high frequencies slightly, which would indicate a
slightly lower slope (b) value than that derived by the described
method. An improved fit results if the value b= -1.5, the value
derived for the ridge crest in Chapter 5, is used. Overall, however,
the match of the true two-dimensional spectrum with the model is quite
good. It is questionable whether the additional detail present in the
true spectrum represents a valuable signal or simply additional noise.
There are several advantages to the model proposed in this study
(derived with respect to azimuth) over the two-dimensional FFT method
105
(constructed in cartesian coordinates). The proposed model requires
only four parameters (b, u, v and ee) to describe the surface rough-
ness. The two-dimensional spectrum in this case requires a 128 x 128
array, or 16,384 parameters. In addition, the four parameters used in
the model have physical meaning attached to them which may prove useful
in comparisons of different areas. Also, the computer algorithms used
to generate the model require far fewer calculations than the direct
transform method.
Perhaps the most relevant advantage in the azimuthal model con-
struction is that the two dimensional nature of the surface can be esti-
mated from randomly oriented ship tracks. Each profile yields a one-
dimensional estimate of the amplitude spectrum at the azimuth of the
ship's heading. Such estimates can be thought of as cross sections
through the surface contoured in Figure 6-13. For example, the profile
and amplitude spectrum shown in Figure 6-7 represent a cross section of
the two-dimensional surface collected at azimuth NIL15°E. Given a suf-
ficient number of such randomly oriented tracks over an adequate range
of headings, the model can be constructed as described in this chapter.
Until a great many more multibeam surveys have been collected, this
method of estimating the anisotropy of bottom roughness will remain the
only available method over most of the world oceans.
106
' sein eek sada ete oth ‘Salas 3!
| sil led ipl rac itt hw >yoaet Wwenny eahlyies "2 ahan
operon il oweah ied). from cage tance ph ine piasah mt
de een tis pee nel it ‘ahh 9
SE ial neh can ee eb ei ‘stella did day Ne ana
ee uvepeebtndba ca tiads
ee mainte” tp yd aorta ge asia mie No Hehehe Fe ; Re
| enn nore NER bepmeanerat oa Nabil “gabe 284, por dammeyng
pleat oa sit quate! some ee i ee cd aontnete a
Wines GS AS abkadier ven Poesia! Ad
es Nb were. voted iarpablerbe es ‘NARA PO es
bien nit ie sent aq Devo eked Dates ee oq saglet
hearggh> dats ait ‘wou pdonie iweeieatle. sabe 4
ayis? site th sabe euals ee ae att
ilar wERiy: iventader: anced ia tei
a wees ninadiia & sve een ott Voit ot
ee i pidge aaue ay napa 2 a hip semoiele “eee ——
team af Whe UKE. Retest menae Mel gue ionua ott the) otek, te gute ;
qorel... “ha.” Ca sighs lenis ee ae toned peenenn an: bi
-
tae aii Lea Oyiitay ieee ge Be Lge ee aking oA zie ‘iia ©
4 Ge uvterss nO Uh aad tw ei nodal ule aad | es) ‘ee: whee
(derives <tth elipers & eet] Me en ee
. a L
i
7. Prediction of High Frequency Roughness
Having developed a spectral model of sea-floor topography based on
measurements from surface ship sonar systems, the question remains
whether this model can be extrapolated into spatial scales smaller than
those resolved by the sounding system. This question is particularly
important to underwater acoustic applications, where the acoustic fre-
quencies of interest in a scattering problem do not necessarily corre-
spond to the spatial frequencies sampled to generate the model. The
concepts of measurement noise levels were introduced in Chapter 4. In
this chapter, the effect of estimation errors on prediction will be
examined, sources of measured high frequency bathymetry introduced, and
a simple prediction test presented.
Source of Error in Spectral Estimates
The model parameters used to describe the amplitude spectrum of sea-
floor topography are derived from regression estimates of spectra from
profiles of noisy data collected in a generally non-stationary environ-
ment. As such, there is estimation error in the model parameters from
several sources, which will necessarily result in prediction errors as
the model is extrapolated to high frequencies. Although many techniques
are used to reduce these errors, some level of error will always remain.
Unfortunately, due to the variability of data quality, variable track-
line spacing, and the presence of some level of non-stationarity in a
data province, it is not possible to quantify completely the estimation
107
error associated with the model. Only by comparing estimated values
with values measured in high frequencies over many data sets and geo-
logic environments, can a reasonable statistical base be assembled to
assess the prediction capabilities of the model quantitatively.
The effect of instrument noise (illustrated in Mgure 4-5) when
encountered by the signal spectrum, has the effect of reducing the
spatial frequency range of amplitude estimates available to the regres-
sion analysis. Certain spectra examined in the course of this study
varied from complete white-noise spectra to spectra showing only two or
three amplitude estimates above the noise. Such spectra are of no use
in model generation. The length of the spectrum available for regres-
sion analysis depends as well on the length of data available to the
FFT. Long profiles from large quasi-stationary provinces produce cor-
respondingly long amplitude spectra and therefore more reliable model
estimates. As mentioned in Appendix B, a minimum of one hundred points
is allowed in a profile for analysis. The use of a higher minimum pro-
file length, while improving regression estimates, would allow more non-
stationarity in the profile and reduce the ability to resolve small
provinces.
No attempt has been made to quantify the level of non-stationarity
(as defined for this study) by examining the mean variability of the
spatial domain estimates used in province picking (see Appendix B). Due
to the nature of the sea floor, certain provinces appear over thousands
of square miles, while others fall smaller than the hundred point mini-
mum required for analysis and must be combined. Such variability prob-
ably precludes estimating the degree of non-stationarity with any accu-
racy; we can only attempt to constrain the effect with the province
108
picking procedure. In general, however, one can treat estimates from
large provinces of persistent geological processes as more reliable than
those generated in a relatively small area.
One source of error in the model can be quantified, and that is the
residual error from the regression. Using standard statistical tech-
niques (see, for example, Draper and Smith, 1981), the difference in
individual amplitude estimates from the regression model estimates can
be expressed as a root mean square. Errors in this study averaged €, =
t .03 and E1og ae = .015m for a single spectrum. These errors reflect
many of the errors associated with the model, although they cannot be
decomposed into component sources.
Another important factor determining the level of error in the
model is the number of estimates used in generating a composite spec-
trum. In the case of an anisotropic area, a variety of azimuths dis-
tributed about the compass allows a better estimate to be made. In gen-
eral, the ensembling of N time series composed of signal with noise,
results in a decrease of noise (as RMS) of 1/ Vv N. In the case of our
spectral model, the signals are the derived amplitude spectra along an
azimuth. A simple method of improving the prediction capability and
accuracy of this model is to ensemble-average the amplitude versus fre-
quency estimates from several proximal and near-parallel tracks. The
multibeam sonar provides exactly this capability and future developments
should take advantage of it.
Figure 7-1 illustrates this reduction of estimation error through
ensemble-averaging of multibeam sonar derived spectra. For this exam-
ple, sixteen parallel beams (profiles) from the SASS multibeam system
were analyzed for the statistically homogeneous province of the Gorda
109
©e0®@ eo ©
STANDARD ERROR (METERS)
@ee ee
-1.38
030
“1.34 os
[a4
a
A 020
b <
1.30 S
B
010
1.26
000
NUMBER OF POINTS USED IN ENSEMBLE AVERAGES
Figure 7-1 Results of ensembling the spectral estimates (a and b) from
sixteen parallel profiles collected by a multibeam sonar
system. Plotted points indicate the estimated parameters
for each of the sixteen profiles and various ensembles of
adjacent profile parameters. The calculated standard error
(solid line) decreases as VN (dashed line), as predicted by
standard statistical theory.
110
Rise crest. Each derived E-W profile was 17 km in length and spaced
approximately 100 m apart. Amplitude spectra were generated for each
profile, and rather than averaging each amplitude estimate, the derived
regression parameters a and b were assembled for adjacent beams in
groups of 2, 4, 8, and 16. Figure 7-1 plots for both a and b, the
estimated parameters for each of the sixteen profiles and the various
results of ensembling. The standard error of each set is plotted as a
solid line, with the theoretical :!/VN relationship shown as a dashed
line. The derived standard error for the final average of sixteen
points is necessarily zero. Similar techniques could be applied in the
province picking algorithm to improve reliability.
Propagation of Error to High Frequency Estimates
Although the errors associated with the model can not be determined
with any accuracy, it is still instructive to examine how the errors
affect the ability to predict amplitude in frequencies beyond the range
of analysis. Due to the power law form of the model,
A= a 0 ab where A = amplitude
s = spatial frequency
a a
a, b = regression parameters
the errors of estimate (error bars) are not linear. It is therefore
somewhat easier to visualize the function in log A-log s space where the
error bars are linear.
111
As previously stated, the parameter a represents the intercept (in
meters) of the function with frequency log s = 0, or wavelength ’4 = 1
km. The error associated with a (e,) appears as a constant vertical
shift of the regression line in log-log space. It is in fact a multi-
plicative factor in linear space and has the effect of multiplying or
dividing the A value by ([antilog e,! at any frequency. Since €, is
{ndependent of frequency, it is stable over large extrapolations. Since
all spectra of topography are red noise, the absolute level of estima-
tion error effectively decreases at higher frequencies (lower ampli-
tudes).
The spectral slope parameter b te not independent of frequency. In
log-log space, the dimensionless bd appears as the slope of the linear
regression line. Error associated with b (ey) at s = Q, causes an
increasing prediction error at higher or lower frequencies. The rela-
tionship in linear space is also multiplicative and depends on frequency
as multiplying or dividing the A value by [antilog (|log 3| * €,)]. The
total error of estimate requires linearly combining the two sources €.
and €,, which translates into multiplying or dividing the value A(s) by
[antilog (c€, + |log s| ° €,)]. An example of these calculations is
included in the following section.
Comparison of Surface Ship Sonar Results to Deep-Towed
Sonar Results and Results from Bottom Photography
To quantify accurately the ability of the spectral model derive
from surface ship sonar systems to predict amplitudes at high spatial
frequencies requires a large data base of small-scale bathymetry pro-
112
files. At present, such a large data base does not exist, at least for
the meter-millimeter scales of bottom photography. One area of the
world ocean (near 40°27'N and 62°20'W) has been extensively surveyed at
small scales and this is the location of the High Energy Benthic
Boundary Layer Experiment (HEBBLE).
The HEBBLE area falls into the large roughness province of the East
Coast Continental Rise which was described in Chapter 5 and illustrated
in Figures 5-12 and 5-13. A spectral model was generated for this large
area based on averaging spectral estimates from all profiles illustrated
on Figure 5-13. As mentioned previously, the amplitude spectrum for
this area consists of two segments, the lower frequency model (with a =
0.0142 a, b = -1.813) extending to wavelengths of approximately A = 3
km, and the higher frequency model (with a = 0.0602 m, b= -0.603)
extending from A = 3 km to A = 200 a.
Data collected by the Deep-Tow sonar system in the HEBBLE area were
provided by Scripps Institute of Oceanography. The Deep-Tow, which col-
lects profiles from a height of only 25-50 meters above the sea floor,
is able to sample bathymetry at a horizontal sample spacing of 5 meters.
The vehicle is positioned via a transponder navigation system which pro-
vides relative location to the beacons every five minutes (or approxi-
mately 270 meters) of track. All depths recorded by the system in this
area appear to be above the instrument noise level.
Digital height data from one stereo-pair bottom photograph col-
lected in the HEBBLE area were provided by the University of Washington.
The data set consists of six horizontal and six vertical transects sam-
pled at 1 mm spacing. All derived spectra were virtually identical in
their model characteristics.
113
Figure 7-2 presents a composite amplitude spectrum showing data
from all three sources in the HEBBLE area. The SASS derived spectrum
(from the nearest available trackline) spans wavelengths of 100 km to
200 m. The model regression lines, derived from spectra generated
throughout this large province, are shown as straight line segments.
The Deep-Tow derived spectrum, spanning wavelengths of 1 km to 10 on, is
very well predicted by the model regression line. The lower frequency
estimates begin to diverge at A= 300 m., which may be due to positioning
distortions from the transponder navigation system, which is interro-
gated at 5 minute (or ~300 meter) intervals. Finally, the bottom photo-
graph derived spectrum, spanning wavelengths of 25.6 cm to 2 mm, falls
approximately .5 orders of magnitude below the regression line.
The prediction residual of .5 orders of magnitude over 5 decades of
frequency indicates some combination of errors in parameter estimates a
and d. If all errors were in Ay it would be in error by .5 orders of
magnitude. Similarly, the error in b would be 2.1, using the relation-
ships described in the last section. The evaluation of several more
bottom photograph-derived profiles from the same region would allow a
more quantitative treatment of the model prediction error. It is pos-
sible that other high frequency spectra may scatter around the predic-
tion line, due perhaps to actual variability of the microrelief in the
area rather than error in the prediction model. Visual inspection of
the several photograph pairs taken in HEBBLE indicate that other areas
are in fact rougher than the relatively featureless photo analyzed here
(Arthur Nowell, personal communication, 1983).
In any case, the prediction of millimeter scale topographic ampli-
tude from a surface hull-mounted sonar system to within half an order of
114
10°!
NORMALIZED AMPLITUDE (KILOMETERS)
r=)
10°
Figure 7-2
SPECTRAL ESTIMATES FROM SEVERAL SOURCES
DEEP-TOW
10-2 10°! 10 0 10! 102 102 104 105 10°
SPATIAL FREQUENCY (CYCLES/KILOMETER)
Composite amplitude spectrum of bathymetric profiles derived
at various scales by different survey instruments. All data
were collected on the East Coast continental rise near
40°27'N, 62°20'W, proximal to the HEBBLE area. Regression
lines represent the mean bottom roughness model derived from
fourteen SASS profiles illustrated in Figure 5-13. Model
predicts the deep-tow derived spectrum almost exactly over
~1.5 decades of spatial frequency. The spectrum derived
from a stereo-pair bottom photograph is predicted to within
0.5 orders of magnitude over ~ decades of spatial frequency
(200 meters to 2 millimeters wavelength).
115
magnitude is a surprisingly good result. The ability to predict high
frequency roughness from a low frequency model appears to be quite pos-
sible, at least for some areas of the world ocean. Ome serious caution
must be taken into account. The successful prediction from the SASS-
derived spectrum was only possible by identifying the break in slope at
3 km. Were data only available to that 3 km wavelength, the lower fre-
quency estimates would have been used, and prediction error of nearly
six orders of magnitude incurred in estimating the millimeter scale amp-
litudes. If it is assumed that such slope breaks are due to changes in
relief-forming processes at various frequencies, then any estimates of
high frequency roughness exclusively from a lower frequency model pre-
sume a continuity of process over all intervening frequencies, and the
absence of a significant break in the power law form of the spectrum.
116
8. Summary and Conclusions
A method has been developed to allow a valid stochastic description
of sea-floor relief to be generated in a. relatively simple statistical
model. The fundamental statistic used is the amplitude spectrum of spa-
tial frequency, which is both elementary enough to be generated opera-
tionally from existing digital bathymetric data bases, and general
enough to be applied to a variety of engineering applications and scien-
tific problems. The model allows relatively large areas of the world
sea floor to be described by as few as two model parameters for simple
isotropic surfaces. The difficulties of producing a stationary statis-
tie in a non-stationary environment are minimized with the use of a spa-
tial domain provineing technique. The model also accounts for the
directional dependence of anisotropie surfaces. The results of one sim-
ple experiment indicate that the model may be extrapolated with high
fidelity to frequencies beyond the resolving capability of surface ship
sonar systems. This stochastic model when combined with lower frequency
deterministic models, such as the gridded bathymetric models developed
by NAVOCEANO, allows a complete description of the sea-floor relief.
There are numerous avenues available to improve the model. The use
of ensembles of either data or derived statistics allows the estimation
errors of the model to be reduced substantially. The availability of
multibeam sonar systems makes this improvement possible immediately.
The ensembling of the statistics from the sixteen beams of the SEABEAM
system would reduce the estimation error to one quarter the results from
a single beam. Data derived from very narrow beam sonars allow a higher
rate of sampling, and therefore a wider frequency band, to be available
117
for model generation. The collection of more complete areal bathymetric
surveys would allow further refinement of the models of sea-floor aniso-
tropy. Future improvements in deep-towed side-scan sonars and the col-
lection of additional bottom photographs would allow a quantitative
estimate to be made of the predictive ability of the model.
The utility of the model falls into two broad categories; engineer-
ing application and scientific investigation. The application of the
model to underwater acoustics is obvious. The model spectrum of the
surface relief is an important environmental factor in the scattering of
sound from the sea floor. Efforts are currently underway to develop
scattering models which utilize such stochastic environmental informa-
tion. More traditional models require the description of the bottom as
a faceted surface, which can easily be derived from the derivative spec-
trum and a knowledge of the probability distribution of depths. Even
models requiring a fully determined surface can obtain a valid realiza-
tion by combining the model amplitude spectrum with a randomly generated
phase spectrum.
In terms of scientific investigation, the method provides a new
tool for studying the earth's geological and tectonic processes. The
resulting patterns of roughness discovered on the Gorda Rise and Oregon
continental margin illustrate the ability of the method to detect inter-
esting relationships not obvious by simply studying bathymetric charts.
Comparing the distribution of roughness on a variety of spreading cen-
ters, continental margins and other submarine environments, should yield
valuable insight into the distribution of geologic processes on the sea
floor. In the case of the patterns on the Gorda Rise, the distribution
of roughness may well relate to the lava formation processes, which in
118
turn may be related to the distribution of polymetallic sulphide miner-
als formed in association with hydrothermal vents.
The ability to quantify the anisotropy of the sea floor also intro-
duces interesting possibilities for geological investigation. The dis-
tribution and trend of these measurements provides a new tool for the
quantitative investigation of deep-sea processes in various environ-
ments. Also the relationship between sea-floor spreading rate and
roughness, long assumed qualitatively, could be quantified and analyzed
in terms of the operative geological processes.
Finally, the author hopes that this study represents more than sim-
ply a study of sea-floor roughness. Much effort has been made to pre-
sent an approach to modelling natural phenomena which could be applied
to a wide variety of problems in the natural sciences. As most disci-
plines within the earth sciences represent a marriage of one of the
“pure” sciences to the study of the earth, this study represents a
crossing of natural science with engineering statistical methods. Very
little such work has been done by earth scientists in the past, perhaps
because the earth is rarely as well ordered as the controlled labora-
tories of the chemist or physicist. Due to the natural variability of
the phenomena under study, the geologist’s attempt to describe the earth
quantitatively is particularly difficult. It is hoped that the approach
and philosophy presented in this attempt will be of use to future
investigators.
119
hi sortie shit ele cs baci bowl
. worn serait ns olay
| penitent! rua i a EH
‘ge sagt) ast ons wan nnuesianpictnatg hte % oe
Saisie ‘bmi sume! \aannapsaeitceda nesta ; ;
a Waetiew ‘pel ban Ska lh wpe ‘ese nee nao setadde a eb +5 f
ve ‘eat. bat) so lee, aiyreattnge, abut cart Sag ag siaasenhagiey, eae
i" vewtipcet: iis stb hala ia-sartibe eu ite sna sian baat a
ies Se geal © Seti. “eainat wecobnaiys ieee piverabab’ bona"
lane acne ahaa a seen wl ome sau doar, se yforae
7 et Be! wan Rhee ‘ayo the wa a Pian P baie.
aps aad uo ht made have an? ; :
har: Sy gta oa? Mtoe ated sntebbiel ge wate matt rer
ee:
a ieee bail gk: njalsnntos ta ae . ih, ones heed ‘wed Citas ‘down
Peas tesusnsasact st 48 “talstats SRBC n> i aR
means oot ste is: somroyaah ety 06 es te 8
Weenie oom oh per aye eh, Rs AEM
| cusipaalinninar ennai ce pn mote iad
tgs: 8 hil te i ai ena as ms aif
> MBopsigpa tig the Mate mrchotk CF asin wee ‘a : wa
CK g enn eS. Wy tom weak sie
pay Lay Lorn s
: #3 : - i. ie 4 4h ren oe ya ve i ey J
ot cngtone ay Weld | culene
9. References
Agapova, G.P., 1965, Quantitative characteristics of bottom slope angles
in seas and oceans: Oceanology, vol. 5, no. 4, p. 135-138.
Akal, T., and J. Hovem, 1978, Two-dimensional space series analysis for
sea-floor roughness: Marine Geotechnology, vol. 3, no. 2,
p. 171-182.
Atwater, T., and J.D. Mudie, 1973, Detailed near-bottom geophysical
study of the Gorda Rise, J. Geophys. Res., v. 78, no. 35, p. 8665-
8686.
Barker, F.S., DR. Barraclough, and S.R.C. Malin, 1981, World Magnetic
Charts for 1980 - spherical harmonic models of the geomagnetic field
and its secular variation: Geophysical Journal of the Royal Astro-
nomical Society, vol. 65, p. 525-533.
Barnard, W.D., 1978, The Washington continental slope: quaternary tech-
tonics and sedimentation: Marine Geology, v. 27, p. 79-114.
Bell, T.H., 1975a, Topographically generated internal waves in the open
ocean: Jour. Geophys. Res., vol. 80, p. 320-327.
Bell, T.H., 1975b, Statistical features of sea-floor topography: Deep-
Sea Research, vol. 22, p. 883-892.
Bell, T.H., 1978, Mesoscale sea-floor roughness: Deep-Sea Research,
vol. 26A, p. 65-76.
Berkson, J.M,, 1975, Statistical properties of ocean bottom roughness
(abs.), power spectra of ocean bottom roughness: Jour. Acous. Soc.
Amer., vol. 58, no. Sl, p. 587.
Berkson, J.M., and J.E. Matthews, 1983, Statistical properties of sea-
floor roughness. In: Acoustics in the Sea-Bed, N.G. Pace (ed.),
Bath University Press, p. 215=223.
Berry, M.V., and Z.V. Lewis, 1980, On the Weierstrass-Mandelbrot fractal
function: Proceeding of the Royal Society of London, vol. A370,
po 459-484.
Blackman, R.B., and J.W. Tukey, 1958, The Measurement of Power Spectra,
Dover Publications, Inc., New York, 190 pp.
Bracewell, R., 1978, The Fourier Transform and its Applications: McGraw-
Hill Book Col, New York, NY, 391 pp.
120
Brown, Gary S., 1982, A stochastic Fourier transform approach to scat-
tering from perfectly conducting randomly rough surfaces: IEEE
Transactions on Antennas and Propogation, Vol. AP-30, No. 6, p.
1135-44.
Brown, Gary S., 1983, New results on coherent scattering from randomly
gation, Vol. AP-31, No. 1, p. 5-11.
Carson, B., 1977, Tectonically induced deformation of deep-sea sediments
off Washington and northern Oregon: mechanical consolidation, Marine
Geology, v. 24, p. 289-307.
Chase, T.E., P. Wilde, W.R. Normark, C.P. Miller, D.A. Seckins, and J.D.
Young, 1981, Offshore Topography of the Western United States
between 32° and 49° North Latitudes, Plate 2, U.S. Geological Survey
Open File Map 81-443.
Chatfield, C., 1980, The Analysis of Time Series: An Introduction:
Chapman and Hill, New York, 268 pp.
Clay, C.S., and W.K. Leong, 1974, Acoustic estimates of the topography
and roughness spectrum of the sea floor southeast of the Iberian
Peninsula. In: Physics of Sound in Marine Sediments, L. Hampton
Clay, C.S., and H. Medwin, 1977, Acoustical Oceanography: Principles and
Applications: John Wiley & Sons, Inc., 544 pp., New York.
Davis, J.C., 1973, Statistics and Data Analysis in Geology: John Wiley
& Sons, Inc., New York, 550 pp.
Davis, T.M., 1974, Theory and practice of geophysical survey design:
NAVOCEANO TR-13, Bay St. Louis, MS, 137 pp.
Draper, Norman and Harry Smith, Jr., 1981, Applied Regression Analysis,
Second Edition: John Wiley & Sons, New York, 709pp.
Elvers, D., S.P. Srivastava, K. Potter, J. Morley, and D. Sdidel, 1973,
Asymmetric spreading across the Juan de Fuca and Gorda Rises as
obtained from a detailed magnetic survey: Earth and Planet. Sci.
Lett., v. 20, p. 211-219.
Engel, A.E.J., C.G. Engel, and R.C. Havens, 1965, Chemical character-
istics of oceanic basalts and the upper mantle: Bull. Geol. Soc.
Am.; Ve 76, pe 719-734.
Farre, J.A., and W.B.F. Ryan, 1984, A 3=D view of erosional scars on the
U.S. Mid-Atlantic continental margin: AAPG submitted.
Glenn, M.F., 1970, Introducing an operational multibeam array sonar:
Int. Hydrographic Review, 47, p. 35-39.
121
Godfrey, Michael D., 1967, Prediction for non-stationary stochastic
processes, In: Bernard Harris (Ed.) Advanced Seminars on Spectral
Analysis of Time Series: John Wiley & Sons, Inc., New York, N.Y.
Hayes, D.E., and J.R. Conolly, 1972, Morphology of the southeast Indian
Ocean. In: D. E. Hayes (ed.) Antarctic Oceanology II: The
Australian-New Zealand Sector: American Geophysical Union,
Washington, D.C.
Heezen, B.C., and T.L. Holcombe, 1965, Geographic distribution of bottom
roughness in the North Atlantic: Bell Telephone Labs, Inc., and
Lamont=Doherty Geol. Obs. (Columbia Univ., Palisades, NY), 41 pp and
6 charts.
Hey, Richard, 1977, A new class of “pseudofaults” and their bearing on
plate tectonics: a propagating rift model: Earth and Planetary
Science Letters, vol. 37, p. 321-325.
Kanasewich, E.R., 1981, Time Sequence Analysis in Geophysics, University
of Alberta Press, Edmonton, 480 pp.
Krause, D.C., and H.W. Menard, 1965, Depth distribution and bathymetric
classification of some sea floor profiles: Marine Geology, vol. 3,
p- 169-193.
Krause, D.C., P.J. Grim, and H.W. Menard, 1973, Quantitative marine geo-
morphology of the East Pacific Rise: NOAA Tech. Rept. ERL 275-
AOML 10, pe 1-73.
Larson, R.L., and F.N. Spiess, 1970, Slope distributions of the East
Pacific Rise crest: SIO Ref. 70-8 (Scripps Inst. Ocean., Univ.
California, San Diego, CA), 4 pp.
Mandelbrot, Benoit, 1982, The Fractal Geometry of Nature, W.H. Freeman &
Co., San Francisco, 460 pp.
Martin, Marcel A., 1957, Frequency domain applications in data process-
ing: General Electric Co., Missile and Ordinance Systems Dept.
(Philadelphia) Document No. 57SD340, 128 pp.
Matthews, J.E., 1980, An approach to the quantitative study of sea-floor
topography: NORDA TN-47, Bay St. Louis, MS, 55 pp.
McClain, C.R., and H. Walden, 1979, On the performance of the Martin
digital filter for high- and low-pass applications: National
Aeronautics and Space Administration Technical Memorandum 80593,
Goddard Space Flight Center, Greenbelt, Md., 20771.
McDonald, M.F., and E.J. Katz, 1969, Quantitative method for describing
the regional topography of the ocean floor: Jour. Geophys. Res.,
vol. 74, no. 10, p. 2597-2607.
122
McDonald, M.F., E.J. Katz, and R.W. Faas, 1966, Depth study report on
the detectibility of submarines on the ocean bottom: Vol. III,
Electric Boat Div. (Gen. Dynamics Corp., Groton, CT), Rept. U413-66-
116.
Menard, H.W., 1964, Marine Geology of the Pacific: McGraw-Hill Book
Co., New York, NY.
Naudin, J.J., et R. Prud'homme, 1980, L'Analyse cartographique: Etude
numerique des caractéristiques morphologiques des surfaces: Sciences
de la Terre, Paris, Série Informatique Geologique, no. 4, p. 47-71.
Neidell, N.S., 1966, Spectral studies of marine geophysial profiles:
Geophysics, vol. 31, p. 122-134.
Papoulis, A., 1962, The Fourier Integral and Its Applications: McGraw-
Hill Book Co., New York, NY, 319 pp.
Rhines, P.B., 1977, The dynamics of unsteady currents. In: The Sea,
vol. 6, E. D. Goldberg, I. N. MecCave, J. J. O'Brien and J. H.
Steele, editors, John Wiley, p. 189-318.
Rhines, P., and F. Bretherton, 1973, Topographic Rossby waves in a
rough-bottomed ocean: Journal of Fluid Mechanics, vol. 6l,
p. 583-607.
Ross, Sheldon M., 1980, Introduction to Probability Models: Academic
Press, New York, 376 pp.
Scarborough, J.B., 1930, Numerical Mathematical Analysis: The Johns
Hopkins Press, Baltimore, Md., 416 pp.
Searle, R.C., and R.N. Hey, 1983, GLORIA observations of the propagating
rift at 95.5°W on the Cocos-Nazca spreading center: Journal of
Geophysical Research, vol. 88, no. B8, p. 6433-6447.
Shapiro, R., and F. Ward, 1960, The time-space spectrum of the geo-
strophic meridional kinetic energy: J. Meteor., no. 17, p. 621-626.
Smith, S.M., T.E. Chase, W.E. Farrell, I.L. Taylor, and M.E. Weed, 1965,
Distribution and statistical analysis of sea-floor relief in the
northeast Pacific: Bell Telephone Labs and Univ. of California, San
Diego, CA, Tech. Rept. 6-601347, 59 pp. and 8 charts.
Van Wyckhouse, R.J., 1973, Synthetic bathymetric profiling system
(SYNBAPS): NAVOCEANO TR=-233, Washington, DC, 138 pp.
123
Appendix A.l
Having concluded that the frequency spectrum of submarine topography
conforms to a power law functional form, it becomes critical in con-
structing a model based on this statistic to ensure a valid regression
fit to the data. Such a regression analysis is necessary both in fit-
ting the amplitude spectrum itself (amplitude as a function of fre-
quency) and in fitting the energy envelopes for various spatial fre-
quency bands for use in spatial-domain provincing. We represent our
power law function as follows:
y = f(a,x,b) = ax> where y = dependent variable (amplitude)
x = independent variable (frequency)
a,b = regression coefficients
We can use the property of power law functions that they plot linearly
on log-log axes and recast the function in log-log space as follows:
log(y) = log(a) + b log(x)
which can be easily solved by a simple linear regression. Upon deriving
the coefficients 6 and log (a), these terms can be reconverted to
linear-linear space as,
a
b=b
a = antilog (log(a))
124
to reconstruct our power law form,
y = ax
Recalling that the method of least-squares minimizes the total sum
of squares of residual distance from the observed data to the resultant
regression curve, the method just described minimizes these distances in
log-log space. The solution derived in this manner does not minimize
the total residual distance in linear-linear space, although the esti-
mates might be quite close. Methods exist for performing the least-
square fit in either log-log or linear-linear spaces, and these are dis-
cussed here. The choice of method and the use of weighting schemes
depend on the distribution of the data being fit as well as the distri-
bution of the estimation error. In all cases, the error is assumed to
reside in the amplitude estimates, rather than in the independent
variable.
In fitting the energy envelope estimates produced in the delinea-
tion of stationary provinces (see Chapter 5), the regression must be
performed in linear-linear space without weighting. The errors asso-
ciated with the envelope estimates (dependent variable), do not depend
on frequency band (independent variable) and should not be log-trans-
formed. In the case of amplitude spectra however, it was show by
Blackman and Tukey (1958) that estimation error is related to the chi-
squared distribution and that error bars remain constant in log-log
space. Under these conditions, the regression analysis is optimally
performed on log-transformed data.
125
As explained in Scarborough (1930), Article 114., the log-transfor-
mation of the dependent variable (y) causes the residuals in the least-
squares residual equations to be of unequal weight. In the case of a
power law function (which is given as an example in Scarborough (1930)
and will not be reproduced here), the weighting function is the squared
dependent variable (y2). Appendix A.2 presents an ASCII FORTRAN 77 pro-
gram for doing such a weighted regression analysis in log-log space.
Because of the chi-square distribution of errors associated with the
amplitude spectrum, however, the residuals in this case are equally
weighted and no weighting function is required.
In fitting the envelope estimates of discrete band passes for use
in “province picking,” there is no requirement for log transformation of
the data or weighting of residuals. This is the case because the error
associated with all envelope estimates is theoretically constant. In
order to derive properly the regression coefficients a and b in linear-
linear space, it is necessary to use an iterative method. The following
development is modified from Scarborough (1930), Article 115., to apply
to the power law functional form.
We can express the regression coefficients as the sum of initial
estimates and differences as follows:
where a,,bo = initial estimates
4a, db = correction factors
126
It is convenient to use the coefficients derived from performing the
linear regression on log-transformed data as the initial estimates (a>
b,) for the iteration process. If we define a new function in terms of
estimated coefficients as,
) am by
y' = £(x,a,,b,) aox
the discrete values of this approximating function will be,
y"'y = £(x1,895b 9)
y'2 = £(x9,895bo)
¥'n = £(Xq,a9sbo)
where n is the total number of data pairs used in the regression
equations.
Realizing that the coefficients a and b produce the “best fit” solu-
tions, the minimized residuals are represented as,
Vo = £(xo,a,b) = Y2
Ya | £(x, a,b) on Ya
where Y], Y2»:--ceYq are the observed dependent variables. Substituting
our approximation yields
127
A = £(x, > ay + Aa, Ls + Ab) Tiyan Dy areveieyy Th
or
vs +yi = £(x, » a, + Aa, Le + Ab),i = 1.2,...,n
Expanding the right side by Taylor's Theorem for the two variables,
a and b, ytelds
ae, af,
v, +9, = £(x,, a, bo) + dalse- » + Ab(sp mp.) > coop $b WpAyooodt
substituting aay = £(x,5 as? b)) we have
of,
af,
vi tongiriy L Tere + Ab(sp ) tise sitt Lumaladss «ast
Rearranging terms and dropping higher order derivatives yields
of, af,
Vine sags) + Ab) + y' tn ae sb Oo Ap ooonm
We define
™,* y's er mw 12 sicleeigD
where the r's are the residuals for the approximation curve
y' = £(x 229d 5)- Substituting, we can write our residual equations
128
of, of,
WA da(s) + Ab(i=—) + 4, 3; 121,2,...,n
°
{ ob 1
Since this system of equations is linear in the correction terms
4a and Ab, these terms can be derived by the method of least-squares.
At this point, one variation from the description given by
Scarborough (1930) is introduced. Because the linearization of the sys-
tem of equation is only valid for small Aa and Ab, it is possible to
derive correction terms which yield solutions that extend beyond the
local neighborhood of linearization. Under these circumstances, it is
possible that the method will not converge to a valid solution. During
any particular iteration, this non-convergence would appear as an
increase in the total sum of squares of the residuals over the previous
iterations.
To ensure a decrease in the total residuals (and therefore a con-
verging solution) during each iteration, the correction terms are scaled
by a term a to yield,
= + ad
a a
o> @>
= bb + aAb
The scaling term a is first set equal to one, and the calculated
total residuals compared to those calculated in the previous iteration.
If the residuals do not decrease, the scale factor a is halved until a
decrease in the total residuals is observed. These new A and b become
the new “initial estimates” a, and Bo» which are then used in the next
iteration. The process continues until the values of both Aa and Ab
reach some suitable minimum, and a final a and b are derived. It is
129
theoretically possible for the solution to converge to a subsidiary min-
{mum and therefore yield a poor model. However, the selection of the
initial a, and b, by a regression in log-log space makes convergence on
subsidiary maxima unlikely. An ASCII FORTRAN 77 subroutine for perform-
ing this algorithm is presented in Appendix A.3.
In practice, the regression models for the amplitude spectra (and
final roughness model) are formulated interactively with a graphic dis-
play terminal. The operator can interactively edit the bathymetric pro-
file under examination as well as control the frequency limits included
in the regression analysis. This allows the software to delete white-
noise levels, interpolation effects, and other contaminating factors
from the analysis. The interactive control also allows the operator to
detect visually spectra composed of two power law segments, and to fit
each segment individually. This interactive software is included in the
amplitude spectrum generating software listed in Appendix D.
130
ene
bie patient
oe hadnt yack
neh yen ai iy
Wy
PD Bh a hin 8
er a ah ee,
gh 4
f
! ie
ANNRNM AAANANNANNANAN|ANN|ANAN
10
50
Appendix A.2
SUBROUTINE POWWGT(A,B,FIRSTX,DELX,Y,N,ILIST,CTOFF1,CTOFF2)
THIS ROUTINE PERFORMS A BEST FIT TO DATA WITH A POWER LAW
FUNCTION OF THE FORM Y=A*X**B USING A WEIGHTING METHOD
AS. DESCRIBED IN SCARBOROUGH(1930), ART. 114.
INPUTS ARE
A=COEFFICIENT OF X
B= EXPONENT OF X
X= ARRAY OF INDEPENDENT VARIABLE VALUES
Y= ARRAY OF DEPENDENT VARIABLE VALUES
N= NUMBER OF DATA PAIRS OF X AND Y
AMIN= MINIMUM VALUE FOR A CORRECTION TO STOP ITERATION
BMIN= MINIMUM VALUE FOR B CORRECTION TO STOP ITERATION
ILIST= 1 FOR SUMMARY OF ITERATION PROCESS, = 0 , NO LISTING
PROGRAMMED BY C.G.FOX-ADVANCED TECHNOLOGY STAFF ,NAVOCEANO,4/15/83
DIMENSION X(1024),¥(1024)
COMPUTE INITIAL ESTIMATE OF A AND B BY PERFORMING A SIMPLE
LINEAR FIT ON LOG TRANSFORMED DATA
YSQR=1.
YSQSUM=0.0
XPROD=0.0
XSUM=0.0
YSUM=0.0
XSQR=0.0
X(1)=FIRSTX
DO 10 I=2,N
X(1)=X(I-1)+DELX
WRITE(6,'(80X,2F10.4)') (X(I),¥(1),1=1,20)
DO 50 J=1,N
IF(Y¥(J).LE.0.0) GO TO 50
YTEMP=AL0G10(Y(J))
XTEMP=AL0G10(X(J))
IF (XTEMP.GT .CTOFF2.0R.XTEMP.LT.CTOFF1) GO TO 50
YSQR=YTEMP*YTEMP
XPROD=XPROD+( YTEMP*XTEMP*YSQR )
XSUM=XSUM+XTEMP*YSQR
YSUM=YSUM+YTEMP*YSQR
YSQSUM=YSQSUM+YSQR
XSQR=XSQR+(XTEMP*XTEMP*YSQR )
CONTINUE
B=(( YSQSUM*XPROD ) - (XSUM*YSUM) ) /( ( YSQSUM*XSQR ) - (XSUM*XSUM) )
A=(YSUM/YSQSUM)-(B*(XSUM/YSQSUM) )
A=10.**A
RETURN
END
131
ANAND AAANAANANADANANNANN|AINN|A
C
C
C
Appendix A.3
SUBROUTINE POWFIT(A,8,X,Y,N,AMIN,BMIN,ILIST)
THIS ROUTINE PERFORMS A BEST FIT TO DATA WITH A POWER LAW
FUNCTION OF THE FORM Y=A*X**B USING AN ITERATIVE SUID
AS DESCRIBED IN SCARBOROUGH(1930), ART. 115.
INPUTS ARE
A=COEFFICIENT OF X
B= EXPONENT OF X
X= ARRAY OF INDEPENDENT VARIABLE VALUES
Y= ARRAY OF DEPENDENT VARIABLE VALUES
N= NUMBER OF DATA PAIRS OF X AND Y
AMIN= MINIMUM VALUE FOR A CORRECTION TO STOP ITERATION
BMIN= MINIMUM VALUE FOR B CORRECTION TO STOP ITERATION
ILIST= 1 FOR SUMMARY OF ITERATION PROCESS, = 0 , NO LISTING
PROGRAMMED BY C.G.FOX-ADVANCED TECHNOLOGY STAFF ,NAVOCEANO,4/15/83
DIMENSION X(10),¥(10)
COMPUTE INITIAL ESTIMATE OF A AND B BY PERFORMING A SIMPLE
LINEAR FIT ON LOG TRANSFORMED DATA
XN=FLOAT(N)
XPROD=0.0
XSUM=0.0
YSUM=0.
XSQR=0.
RIEL, '(2F10.4)') (X(I),¥(1),121,N)
DO 50 J=1,N
IF(Y(J). EQ. 0.0) GO TO 50
YTEMP=AL0G10(Y(J))
XTEMP=AL0G10(X(J))
XPROD=XPROD+( YTEMP*XTEMP )
XSUM=XSUM+XTEMP
YSUM=YSUM+YTEMP
50 XSQR=XSQR+(XTEMP*XTEMP )
BO=( (XN*XPROD) -(XSUM*YSUM) ) /( (XN*XSQR ) - (XSUM*XSUM) )
AO=(YSUM/XN ) = (BO*(XSUM/XN ) )
A0=10.**A0
COMPUTE SUM OF SQUARES OF THE RESIDUALS
POLD=0.0
DO 100 I=1,N
100 POLD= POLD+((Y(1)-F3(X(I), AO, 80) )**2)
ITERAT=0
NBIS=0
IF (ILIST.EQ.1)WRITE(6,110)
110 FORMAT(' ITERATION # OF BISECTIONS A B
*RES IDUALS**2' )
IF (ILIST.EQ.1)WRITE(6,120)ITERAT ,NBIS ,AO,80,POLD
132
NNO
(ol wi a)
NAD ANQO AND
AANANMD
aAaNM
120 FORMAT(' ',16,10X,14,6X,3(4X,F10.4))
ZERO OUT MATRIX TERMS
150 Al=0.0
B1=0.0
D1=0.0
E1=0.0
G1=0.0
COMPUTE TERMS FOR LEAST SQUARES MATRIX CONSTRUCTION
DO 200 I=1,N
PARTA=F1(X(I) ,A0,B0)
PARTB=F2(X(I) ,A0,B0)
POWF =F 3(X(I) ,A0,B0)
Al=A1+(PARTA**2)
B1=B1+(PARTA*PARTB)
D1=D1+(PARTB**2)
E1=E1+(PARTA*(Y¥(1I)-POWF ))
200 Se a eA
C1=B
COMPUTE CORRECTION TERMS FOR A AND B
DIVSOR=(A1*01-B1*C1)
ACORR=(D1*E 1-B1*G1)/DIVSOR
BCORR=(A1*GI-C1*E1)/DIVSOR
CREATE NEW A AND B
230 A=A0+ACORR
B=B0+BCORR/AO
COMPUTE NEW SUM OF SQUARES OF RESIDUALS WITH NEW ESTIMATES
PNEW=0.0
00 250 I=1,N
250 PNEW=PNEW+((Y(1)-F3(X(1I),A,B) )**2)
TEST FOR CONVERGENT SOLUTION(PNEW < POLD)
IF NOT, BISECT CORRECTIONS AND RECOMPUTE
‘EF (PNEW.LT.POLD) GO TO 300
ACORR=. 5*ACORR
BCORR=.5*BCORR
NBIS=NBIS+1
IF (NBIS.GT.10) GO TO 300
GO TO 230
TEST FOR MINIMUM CHANGE OF A AND B
300 IF(ABS(A-AO).GT.AMIN) GO TO 500
133
IF (ABS (B-B0).GT.BMIN) GO TO 500
GO TO 900
C CORRECTION TERM NOT FINE ENOUGH, START NEW ITERATION
500 ITERAT=ITERAT +1
A0=A
BO=B
POLD=PNEW
IF (ILIST.EQ.1)WRITE (6,520) ITERAT ,NBIS ,A0,B80,POLD
520 FORMAT(' ',16,10X,14,6X,3(4X,F10.4))
NBIS=0
GO TO 150
900 ITERAT=ITERAT +1
IF (ILIST.EQ.1)WRITE (6,920) ITERAT ,NBIS,A0,B0,POLD
920 FORMAT(' *,16,10X,14,6X,3(4X,F10.4))
RETURN
END
FUNCTIONS TO CALCULATE POWER LAW FUNCTION AND PARTIAL
DERIVATIVES WITH RESPECT TO A AND B
FUNCTION F1(X3,A3,B83)
CALCULATE PARTIAL OF A*X**B WITH RESPECT TO A
F 1=X 3**B3
A3=A3
RETURN
END
(=) NAOMONOOM
FUNCTION F2(X4,A4,B4)
C CALCULATE PARTIAL OF A*X**B WITH RESPECT TO B
F2=(X4**B4)*ALOG(X4)
A4=A4
RETURN
END
FUNCTION F3(X5,A5,B5)
C CALCULATE POWER LAW A*X**B
F 3=A5*X5**B5
RETURN
END
134
Appendix B.1l
Before generating amplitude spectra from a statistically non-sta-
tionary sample space, such as the sea floor, one must initially define
provinces which are relatively homogeneous with respect to the frequency
spectrum. Since generating an amplitude spectrum directly by Fourier
transform assumes stationarity over the length of the input series, this
“province picking” procedure must be performed by estimating the spec-
trum discretely in the spatial domain. Although the present application
is new, the concept of estimating spectra in the time/space domain is
not. Blackman and Tukey (1958) referred to such estimates as “pilot
spectra” and describe two methods for their calculation. Godfrey (1967)
also describes a method, very similar to that detailed here, which is
used for predicting non-stationary time series. This section gives
details of the procedure used in this study, presents performance tests
of the algorithm, and include full FORTRAN-77 software for performing
the analysis.
The initial step in processing is identical to that required for
running an FFT. The data must be projected onto a straight-line segment
and interpolated to even increments in distance. The straight line is
generated by a simple least-squares fit of the navigation track; the
best results are obtained when relatively straight line navigation is
input. Large deviations in the track greatly degrade results. Next,
points are mapped onto the nearest location on the least-squares line,
without alteration if they are within a designated “pivot distance,” or
modified according to the local gradient if beyond this distance.
Finally, the data are interpolated at a specified interval using a one-
135
dimensional cubic spline, which tends to preserve frequency content.
The full algorithm is included in SUBROUTINE MAPCTIN.
The subsequent stage of processing requires band-pass filtering of
the interpolated data at ten nearly equi-spaced frequency bands. Fil-
tering is done by sequential application of low-pass and high-pass
filters, using a non-recursive, symmetric, least-squares filter devel-
oped by Martin (1957). The frequency response of this bank of filters
is illustrated in Figure B-l. A review of the performance of the Martin
filter can be found in McClain and Walden (1979). The filter cutoffs
are designed to juxtapose at the 100% energy pass level. The total fre-
quency bank was selected to span .02 - .25 cycles/data interval, or
wavelengths of approximately 2 - .16 nautical miles for data recorded at
12 second intervals, or 10 - .8 nautical miles for one minute data. The
filtering is performed by SUBROUTINE FILTER.
The next step of the algorithm requires estimating the instanta-
neous amplitude of all ten band-passed signals. Davis (1974) used a
simple full-wave rectification, followed by a low-pass smoother to esti-
mate the energy envelope. For this study, a true Hilbert transform is
performed and manipulated to generate the energy envelope. The reader
is referred to Kanasewich (1981) for a full development. Notice that
the calculation of the envelope via Hilbert transform does require oper-
ating in the frequency domain. However, because the complete Fourier
Transform is retained throughout, the requirement of stationarity is not
applicable. The enveloping algorithm is contained in SUBROUTINE ENVEL.
The envelope generated in the previous step represents a continuous
estimate of amplitude through space, for each selected frequency band.
To estimate the full amplitude spectrum at each point on the profile,
136
% RESPONSE
FREQUENCY (cycles/data interval)
Figure B-1 Simplified frequency response of the bank of band-pass fil-
ters used to estimate amplitude spectra discretely in the
spatial domain. Actual responses include side lobe leakages
of less than 5%.
137
all ten amplitude versus frequency estimates are considered. Since the
power law form of the spectrum is known, this functional form is used in
an iterative regression procedure described in Appendix A and performed
in SUBROUTINE POWFIT. The regression coefficients vary widely along the
profile, and therefore it is desirable, in order to aid in the detection
procedure, to smooth these estimates. To save calculation time, this is
done by averaging 91 envelope estimates at each frequency before enter-
ing the regression routine. This process tends to smear. the boundaries
somewhat, but makes detection of provinces more reliable.
With smoothed estimates of the amplitude spectrum available contin-
uously along a profile, the final stage of processing is to detect sig-
nificant chances in the estimated spectra and on this basis impose prov-
inee boundaries. The regression coefficients a and b represent the
antilog of the y intercept, and the slope, respectively, of the spectrum
projected in log-log space. The spectral slope, b (which is related to
the so-called Fractal dimension) is a worthwhile parameter for province
detection. A simple algorithm is run across the slope estimates to
detect significant, rapid changes (boundaries).
The regression coefficient a is correlated to b and therefore does
mot represent an independent parameter for detection. An alternative is
to look at the total, band-limited RMS energy of the estimated spectra
for significant shifts in total energy. These RMS parameters are easily
calculated using Parseval's Formula (integrating the power spectrum)
applied to the estimated spectra. In addition to detecting rapid
changes as with the slope, the RMS is contoured to form segments of some
minimum size. Although these detectors have proven fairly reliable, it
is often necessary to interpret certain boundaries where the various
138
detectors disagree. All of the pertinent software is listed in Appendix
B.2. The program runs with a core size of approximately 540 K bytes,
and was written in ANSI Standard FORTRAN (1977) for a UNIVAC 1180/2 com-
puter. The maximum profile length is 4100 interpolated data points and
145 points are lost from the end of the signal due to filtering.
Just as an electrical engineer investigates the performance of an
instrument by operating on signals of known properties, the same tech-
nique can be used here to test the performance of the province picker
algorithm. While an engineer might use a step, ramp, or impulse func-
tion as input, random signals with known spectral forms are used in this
analysis. Many of the seemingly arbitrary choices of filters, averaging
procedures, and other design decisions incorporated into the present
province picker, were selected through feedback from performance tests
with known signals.
An obvious choice of a signal with a known amplitude spectrum is
random “white noise,” with a continuous spectrum of zero slope. There
are several means of producing such a signal. The simplest is to gener-
ate pseudo-random noise series of either a normal or uniform distribu-
tion. Figures B-2 and B-3 illustrate such signals with their spectra,
which are indeed relatively flat. Notice the large amount of scatter in
the amplitude estimates. An alternative method of generating “white
noise” is actually to use a constant amplitude spectrum and uniformly
distributed random phase spectrum, separate their real and imaginary
parts, and inverse-transform the signals using the FFI. In this manner
a perfectly flat, non-varying spectrum is assured, as is illustrated in
Figure B-4. This is the signal source used in testing for this study,
and the algorithm is presented in SUBROUTINE WHINOL.
139
UNIFORMLY DISTRIBUTED NOISE
-200
DEPTH (METERS)
NOs 212 4 nemtol 18) S202
DISTANCE (KILOMETERS)
10°!
10-2
NORMALIZED AMPLITUDE (KILOMETERS)
10-2 10°! 10° 10! 102
FREQUENCY (CYCLES/KILOMETER)
Figure B-2 An example of uniformly distributed random noise is shown in
upper profile. Its corresponding amplitude spectrum (in
lower profile) shows the flat (zero slope) form indicative
of white noise. Notice the high degree of scatter in ampli-
tude estimates. The same computer software was used as that
used for bathymetric data.
140
NORMALLY DISTRIBUTED NOISE
-100
DEPTH (METERS)
fo)
Op EZ Oa NnOink ShpaslOwe A2wil4 ja Tog) SIS) aad AY22 E24
DISTANCE (KILOMETERS)
107!
1072
r=)
rs
ra)
Fs
NORMALIZED AMPLITUDE (KILOMETERS)
°
w
107? 10-1 10° 10! 10?
FREQUENCY (CYCLES/KILOMETER)
Figure B-3 An example of normally distributed random noise in the same
format as Figure B-2. Notice the tendency of values to
cluster about the mean (normal distribution) and the large
scatter of amplitude estimates in the amplitude spectrum.
141
INVERSE TRANSFORM NOISE
1500
-825
-150 74
§25
DEPTH (METERS)
OP) QPF" GET NSS OT Gey NO Pel 2raralamenor 18 -20F 92279924
DISTANCE (KILOMETERS)
NORMALIZED AMPLITUDE (KILOMETERS)
10-2 10°! 10° 10! 102
FREQUENCY (CYCLES/KILOMETER)
Figure B-4 An example of random “white noise" generated using an
inverse Fast Fourier Transform. The method of calculation
forces a perfectly constant amplitude spectrum. This “white
noise" generator was used for performance testing of the
province picker algorithm.
142
Figure B-5 illustrates the performance of the province picker when
white noise is input. Examination of this output results in some inter-
esting insights into the nature of stochastic processes. Despite the
known property that the input signal has a perfectly flat, non-variable
amplitude spectrum, all ten band-pass filtered signals show large fluc-
tuations in amplitude. The amplitude recorded in the frequency spectrum
represents simply an average amplitude computed over the length of the
input time series. At any discrete point, the amplitude at.a frequency
could be widely removed from the average value. It was discovered
through experimentation that these fluctuations were damped when wider
band-pass filters were used; thus, the overlapping filter bank illus-
trated in Figure B-l was designed. The slope values (plotted above the
ten band-passed signals) do fluctuate about zero as expected. The
standard deviation about zero averages about 0.2 over several runs which
implies (assuming a no.mal distribution) that a change of wsiope of t0.4
can be detected with 95% confidence. Many of the decisions made in gen-
erating the regression analysis and smoothing, were designed to minimize
this fluctuation.
In order to design a province detector properly, it is necessary to
combine known signals of differing spectral slopes and RMS energies.
For the purpose of this study, the resulting signal must also have an
amplitude spectrum with power law form. One method of generating such a
signal is through a Markov chain with a probability transition matrix
which allows subsequent events either to raise or lower a constant
increment with probability of 0.5 (see Ross, 1980). This signal, which
is a special case of a random walk, fluctuates about its initial value
and has a spectral form of
143
° BAND PASS 10
BAND PASS 9
BAND PASS 8
BAND PASS 7
BAND PASS 6
BAND PASS 5
BAND PASS 4
BAND PASS 3
BAND PASS 2
BAND PASS |
INPUT SIGNAL
DISTANCE (NAA)
Figure B-5 Example of province picker output with "white noise" input
illustrated in Figure B-4. The input signal is shown at the
bottom, above which are the output of ten band-pass filters
convolved with the data. The lowest signal (Band Pass 1) jis
the lowest frequency pass, while the highest (Band Pass 10)
is the highest frequency pass. Above each band-passed sig-
nal is the energy envelope calculated by the Hilbert Trans-
form method. Notice the large variability in energy for each
band despite the constant amplitude input. At the top, the
estimated spectral exponent (slope of log-transformed spec-
tra) and the band-limited RMS energy calculated along the
profile are plotted. Standard deviation of the slope par-
ameter is .2. The algorithm was designed to minimize this
"natural" variability.
144
It is now possible with two known signals of different spectral
characteristics to study the response of the province picker as it
crosses the boundary Secsesn two such signals (provinces). Since it is
also necessary to detect changes in total RMS energy without a change of
slope, each signal can be multiplied by a constant using a corollary of
the addition theorem of Fourier transforms (see Bracewell, 1965).
Figure B-6 illustrates an artificially generated random signal
which consists of four distinct provinces. The first half of the signal
is composed of “white noise" and the second half of 1/s noise generated
via the random walk technique. Notice the rapid change in the slope
parameter from 0 to -l at the boundary, which was easily detected.
Within each half of the signal, the provinces are again divided
with the second halves (second and fourth quarters) representing a dou-
bling of the first halves (first and third quarters). Notice the
obvious change of RMS energy, although the slope parameter is unaf-
fected, illustrating the importance of detecting on two uncorrelated
parameters. Several false alarms are observed on the derivative detec-
tor, but the province boundaries derived from contouring (represented by
straight lines above the RMS energy profile) correspond to the known
boundaries in the signal. Notice that provinces one and four show the
same RMS energy level, but are delineated by their differences in spec-
tral slope.
As stated earlier, the automatic detection is an aid to province
boundary recognition, but in some cases, human intervention is needed to
resolve inconsistencies. Also, the setting of contour interval or slope
145
PROVINCE 1 PROVINCE 2
PROVINCE 3 PROVINCE 4
BAND PASS 10
BAND PASS 9
BAND PASS 8
BAND PASS 7
BAND PASS 6
BAND PASS 5
BAND PASS 4
BAND PASS 3
BAND PASS 2
BAND PASS 1
INPUT SIGNAL
Figure B-6 Example of province picker output with mixed input of known
signals. The first half of the signal is white noise (A =
as9), joined to the second half of random walk generated
noise (A = as~') the amplitude of each half is doubled at
the mid-point. The break in spectral slope parameter is
easily detected. The change in level of RMS energy is con-
toured, as shown by straight line segments above the RMS
energy profile.
146
threshold can =. modified by the investigator depending upon the amount
of resolution desired. Figure B-7 is a sample profile from real bathy-
metric data collected by the SASS system in the vicinity of the Gorda
Rise. Boundaries are less distinct than in the artificially generated
signal, as one would expect, however distinct changes of RMS energy (and
in two cases changes of slope) do define quasi-stationary provinces. In
practice, the slope parameter has offered little independent information
for province picking, and therefore, construction of a simplified prov-
ince detector based on RMS energy alone may prove to be adequate.
147
*Suaqoweued
Lapow aonpoud 03 pazAleue-uatunoy swaze, aue uawbas aout
-Oud yora 03 Bbulpuodsauu0d saipjoud e3zeg “Laray ABuaua swy
ug sabueyd Aq pauljyap aue saluepunog aduLaoud ysow °yndut
Asqawkyzeq SSyS YatiM yndyno waydid adujaoud yo aydwexy /-g aunbLy
148
AA AAA AAA AA ARAAAARAAAAAAAAAAARAAAAANRAAAAARAAARAARAAAARAANNDNANM|AN|ND
Appendix B.2
PROGRAM TO SELECT PROVINCES FOR SPECTRUM GENERATION
PROGRAMMED BY C.G. FOX,ADVANCED TECHNOLOGY STAFF ,NAVOCEANO,5/15/84
PROVINCE SELECTION IS BASED ON A SPATIAL DOMAIN ESTIMATE OF THE
FUNCTIONAL DESCRIPTION OF THE SPECTRUM BASED ON A POWER LAW MODEL
AMPLITUDE= A*FREQUENCY**B
THE SPECTRUM PARAMETERS A & B ARE ESTIMATED BY PERFORMING A
REGRESSION ANALYSIS ON THE ENERGY ENVELOPES GENERATED ON TEN
BAND PASSES OF THE DATA SPACED EVENLY IN FREQUENCY. PROVINCE
BOUNDARY SELECTION IS BASED ON CHANGES IN THE BAND LIMITED
RMS LEVEL OF EACH LOCATION ALONG THE TRACKLINE AS DERIVED
BY INTEGRATING THE RESULTS OF THE REGRESSION MODEL.
PROGRAM CONTROLS ARE ENTERED IN FREE FORMAT FOLLOWING @XQT
JTOTL= NUMBER OF TRACKLINES OF DATA TO BE ANALYZED IN THIS ANALYSIS
MAXIMUM IS CURRENTLY SET TO 50
IPLOT= 1, PRODUCE PLOT TAPE ON UNIT 10 OF INTERPOLATED DATA,
OPTIONAL BAND PASS PLOTS, PROVINCES SPECTRAL PARAMETER
ESTIMATES, AND PROVINCE BOUNDARIES.
OUTPUT CAN GO TO ZETA OR CALCOMP PLOTTER
0, NO PLOT TAPE IS PRODUCED
IBAND= 1, PLOT ALL BAND PASSED SIGNALS WITH ENVELOPES
0, NO BAND PASSES
IGNORED IF IPLOT=0
ISAVE1= 1, SAVE DATA WITHIN EACH PROVINCE ON UNIT 13 FOR LATER
ANALYSIS IN FFT PROGRAM.
0, NO DATA SAVED
ILIST= 1, PRODUCE COMPLETE LISTING OF ANALYSIS
O, ONLY RUDIMENTARY LISTING IS PRODUCED
GRID= DESIRED SPACING OF INTERPOLATED DATA IN NAUTICAL MILES
THIS SPACING IS NORMALLY SLIGHTLY LESS THAN THE ORIGINAL
SPACING OF THE RAW DATA.
KPTS= NUMBER OF RAW DATA VALUES TO BE ANALYSED IN EACH PROFILE.
READING OF THE DATA WILL CEASE WHEN THIS NUMBER IS REACHED
DATA INPUT IS CURRENTLY DESIGNED SUCH THAT FOLLOWING THIS CONTROL
CARD, THE PROGRAM EXPECTS A CARD IMAGE CONTAINING THE NAME(LESS
149
ANANANANANDNANMNND
(o> en op es Se SP SP)
THAN 20 CHARACTERS) OF A FILE IN WHICH IS CONTAINED A LIST OF
ALL FILE NAMES TO BE USED IN THE ANALYSIS. AFTER READING THIS
FILE NAME, JTOTL FILE NAMES ARE READ FROM THE FILE AND THE
PROGRAM AUTOMATICALLY OPENS THESE DATA FILES AS NEEDED. USERS
MUST INSURE THAT THE FILE NAME LIST FILE AND ALL DATA FILES
ARE AVAILABLE TO THE RUN AT EXECUTION. ONE CONVENIENT WAY
TO PROVIDE DATA TO THE RUN IS BY CREATING DATA AND LIST
ELEMENTS IN A PROGRAM FILE AND THEN CREATE AN @ADD ELEMENT TO
COPY ALL ELEMENTS INTO TEMPORARY FILES BEFORE EXECUTION
PARAMETER (ISIZ=7000)
DIMENSION X(ISIZ),Y(ISIZ),Z(ISIZ) ,AX(ISIZ) ,AY(ISIZ),AZ(ISIZ),
1WGT (240) ,K1(1000) ,K2(1000) ,K3(1000)
2,FRLOW(10) ,FRHIGH(10) ,ENVEST(10) ,TZ(1SIZ,12)
3,FRMEAN( 10)
4, ENERGY (ISIZ),TZ11(1S1Z),1Z12(1S1Z)
EQUIVALENCE (TZ(1,11),1Z11(1)),(7Z(1,12) ,1Z12(1))
CHARACTER*20 FILE(50)
CHARACTER*20 INFILE
FRHIGH=HIGH FREQUENCY FOR DATA BANDPASS
FRLOW= LOW FREQUENCY FOR DATA BANDPASS
DATA FRLOW/.002, .023,.044, .065, .086,.107,.128, .149,.170,.191/
DATA FRHIGH/.041, .062, .083,.104,.125, .146, .167, .188, .209, .230/
SET FILTER PARAMETERS FOR BANDPASSES
FILSLP=.008
NUMF IL=50
NFIT IS LENGTH OF AVERAGE OVER ENERGY ENVELOPES
NFIT=91
C2MAX=0.0
SET PARAMETERS TO CONTROL DERIVATIVE PICKS ON SPECTRAL PARAMETERS
SLWIN=.4
SLWID1=1.
SLWID2=14.
ENWIN=400.
ENWID1=20.
ENWID2=100.
PDEL DETERMINES THE INTERVAL OF RMS ENERGY USED FOR PROVINCE
DETECTION
PDEL=.03
Oat otc NUMBER OF POINTS ALLOWED WITHIN A PROVINCE
MIN=9
OSCALE=SCALING FACTOR(INCHES/UNIT) OF ORIGINAL DATA FOR PLOTTING
OSCALE=-1.
ESCALE= SCALING FACTOR FOR ENVELOPE UNITS
ESCALE=.0001
READ(5,*) JTOTL,IPLOT, IBAND,ISAVE1,ILIST,GRID,KPTS
IKONT =0
LTPS=1000
GRIDKM=GRID*1.852
IF(IPLOT.NE.1) GO TO 24
INITIALIZE PLOTTER
150
aaNnD
a-R
USE FOLLOWING CALL FOR CALCOMP PLOTTER OUTPUT
CALL PLOTS(0,0,10) |
USE FOLLOWING CALL FOR ZETA PLOTTER OUTPUT
CALL PLOTS(53,0,-10)
CALL FACTOR(.5)
CALL PLOT(1.0,1.0,-3)
READ NAME OF FILE CONTAINING LIST OF DATA FILE NAMES
24 READ(5,'(A20)') INFILE
25 OPEN(UNIT=18,FILE=INFILE,STATUS='OLD' )
READ LIST OF DATA FILE NAMES FROM FILE INFILE
00 5 IK=1,JTOTL
READ(18,'(A20)') FILE(IK)
5 CONTINUE
CLOSE(UNIT=18)
READ X,Y,Z, VALUES FOR THIS TRACK.FIRST POINT IS
ASSUMED TO BE ORIGIN,PUT END OF FILE AT THE END OF EACH TRACK.
440 IKONT=IKONT+1
ZERO OUT TEMPORARY ARRAYS
DO 441 J=1,10
DO 441 beens
441 TZ(1, ge
442 X(I)= 0. 0
IAVPLT=0
INPUT X(I)=LONGITUDE(DEC. DEG.),Y(1)=LAT,Z(1)=DEPTH
N=# OF POINTS - 1
IKNT=IKONT
JPTS=KPTS
IF(ILIST.EQ.1) WRITE(6,45) FILE(IKONT)
45 FORMAT(' OPENING FILE ',A20)
ROUTINE TO READ LAT,LON,DEPTH FROM FILE # IKNT
CALL PROVRD( JPTS,Y,X,Z,IKNT,FILE)
4 N=JPTS-1
IF(ILIST.EQ.1) WRITE(6,'(24H NUMBER OF INPUT POINTS ,I5)")N
INTERPOLATE DATA ON STRAIGHT LINE
CALL MAPCTN(GRID,X,Y,Z,N,AX,AY,AZ, JCT)
eS TE Dee en NUMBER OF INTERPOLATED POINTS ,16)')
Ci
839 C=JCT
AZAVE=0.0
PLOT INTERPOLATED DATA
DO 261 I=1,JCT
261 AZAVE=AZAVE+(AZ(1)/C)
IF(IPLOT.NE.1) GO TO 26
XL=(C*.02)+0.5
DV=50.*GRID
CALL NEWPEN(1)
CALL AXES(0.,0.,12HDISTANCE(NM),-12,XL,0.,1.,0.,DV,-1)
EV=AZAVE-(3.0/0SCALE )
DV=1.0/0SCALE
CALL AXES(0.,-3.,13HORIGINAL DATA,!#,¢.,().,!.,EV,DV,1)
CALL SYMBOL(XL+.75,-.25,.5,12HINPUT SIGNAL,0.0,12)
IF(C2MAX.LT.XL) C2MAX=XL
CALL PLOT(0.0,(AZ(1)-AZAVE)*OSCALE, 3)
00 26 [=1,JCT
AJ=I-1
C=AJ*.02
D=(AZ(1I)-AZAVE )*OSCALE
CALL PLOT(C,D,2)
26 CONTINUE
00 27 I=1,JCT
27 Z(1)=AZ(1)
IF(IPLOT.EQ.1)CALL PLOT(0.0,1.0,-3)
DO 237 IRUN=1,10
XRUN=IRUN
ENVMAX=0.0
IF(IPLOT.EQ.1.AND.IBAND.EQ.1)CALL PLOT(0.0,4.0,-3)
C BAND PASS FILTER DATA AT SPECIFIED INTERVAL
CALL FILTER(AZ,JCT,FRLOW( IRUN) ,FILSLP,NUMFIL,1,WGT)
CALL FILTER(AZ,JCT ,FRHIGH( IRUN) ,FILSLP ,NUMFIL,O,WGT )
IF(IPLOT.EQ.1) CALL PLOT(0.0,0.0,3)
WGTMN=0.0
XNORM=( (FRHIGH( IRUN)+.2)-FRLOW( IRUN) )*(1./GRIDKM)
C XNORM=1.
NTWICE=2*NUMF IL
DO 31 I=NTWICE,JCT-NTWICE
Geeel
IF(IPLOT.NE.1.0R.IBAND.NE.1) GO TO 31
Gasol
AJ=1-1
C=AJ*.02
B=(AZ(1)/XNORM)*ESCALE .
IF(I.EQ.NTWICE) CALL PLOT(C,B,3)
CALL PLOT(C,B,2)
31 AZ(1)=(AZ(1))/KNORM
(eS
IF (IPLOT.NE.1.0R.IBAND.NE.1) GO TO 32
kasi
EV=-1.5/ESCALE
DV=1./ESCALE
CALL PLOT(0.0,0.0,3)
GAUIPAXES(OMMONM LHI SL XUN OM SIERO ORR)
CALL AXES(XL,-1.5,13HFILTERED DATA, -13,3.,90.,1.5,£V,DV,0)
CALL SYMBOL(XL+.75,-.25,.5,9HBAND PASS,0.0,9)
CALL NUMBER(XL+5.75,-.25,.5,XRUN, 0.0, -1)
CALL PLOT(0.0,0.0,3)
CALL NEWPEN(2)
C COMPUTE APPROXIMATION TO ENVELOPE OF BAND PASSED DATA
32 CALL ENVEL(AZ,JCT)
152
Gece
Gast
321
322
44
836
837
840
841
845
846
850
IF(IPLOT.EQ.1) CALL PLOT(0.0,0.0,3)
ENVAVG=0.0
IF(IPLOT.NE.1.0R.IBAND.NE.1) GO TO 321
DO 44 I=NTWICE,JCT-NTWICE
AJ=I-1
C=AJ*.02
B=AZ(1)*ESCALE
IF(I.EQ.NTWICE) CALL PLOT(C,B,3)
CALL PLOT(C,B,2)
GO TO 322
DO 44 I=NTWICE,JCT-NTWICE
AZ(1)=ABS(AZ(I))
IF(AZ(1).GT.ENVMAX) ENVMAX=AZ(T)
ENVAVG=ENVAVG+(AZ(1)/( JCT-2*NUMF IL) )
ENVELOPE ENERGY IS STORED IN THE APPROPRIATE COLUMNS
OF 2-DIMENSIONAL ARRAY TZ
TZ(1, IRUN)=AZ(1)
FRDIFF =FRHIGH( IRUN) -FRLOW( IRUN)
IF(IPLOT.EQ.1) CALL PLOT(0.0,0.0,3)
IF(IPLOT.EQ.1) CALL NEWPEN(1)
DO 836 [=1,JCT
AZ(1)=Z(T)
CONTINUE
LINEAR REGRESSION OF TEN ENVELOPE ESTIMATES FOR EACH POSITION
FIRST COMPUTE AVERAGE ENVELOPE
DO 841 I=1,10
FRMEAN( I) =(( (FRLOW(I)+.017)+FRHIGH(I))/2.)
CONVERT CYCLES/DI TO CYCLES/KM
FRMEAN(I)=FRMEAN(I)*(1./GRIOKM)
FR1=(FRLOW(1)+.017)*(1./GRIDKM)
FR2=FRHIGH(10)*(1./GRIDKM)
DO 850 I=NTWICE , JCT-(NTWICE+NF IT)
DO 846 J=1,10
ENVEST(J)=0.0
FIT=FLOAT(NFIT)
DO 845 K=1,NFIT
ENVEST(J)=ENVEST(J)+TZ(I+K, J)
CONT INUE
IF(FIT.EQ.0.0) GO TO 846
ENVEST(J)=ENVEST(J)/FIT
CONTINUE
PERFORM ITERATIVE REGRESSION,. STORE B & A IN
COLUMNS 11 & 12 OF ARRAY TZ
CALL POWFIT(TZ12(I+((NFIT/2)+1)) ,TZ11(1+((NFIT/2)+1)) SFRMEAN( 1),
*ENVEST(1),10,0.000001, .01,0)
PLOT SLOPE PARAMETER B
IF(IPLOT.NE.1) GO TO 859
CALL PLOT(0.0,4.0,-3)
CALL ‘AXES(0.0,0.0,1H ,1,XL,0.,1.,0.,0.,-2)
CALL NEWPEN(2)
153
CALL AXES(XL,-2.0,5HSLOPE,-5,2.0,90.,1.0,-4.0,2.0,-1)
CALL PLOT(0.0,0.0,3)
859 AVSLP=0.0
KKNT=0
SLWID=SLWID2-SLWID1
ENWID=ENWID2-ENWID1
INWID1=IF IX (ENWID1)
LSWIDL=IFIX(SLWID1)
INWID2=IF IX (ENWID2)
LSWID2=IF IX (SLWID2)
DO 877 K=(NTWICE+1)+(NFIT/2), JCT-((NTWICE+1)+(NFIT/2))
KKNT=KKNT+1
C=(K-1)*.02
B=(TZ(K,11))/2.
AVSLP=AVSLP+(TZ(K,11)/(JCT-(4*NUMF IL+(NFIT+1))))
IF(IPLOT.NE.1) GO TO 876
IF (KKNT.EQ.1) CALL PLOT (C,B,3)
CALL PLOT(C,B,2)
876 IF(KKNT.LE.LSWID2) GO TO 877
IF (K.GE. (JCT-(NTWICE+NFIT/2+1+LSWID2))) GO TO 877
PREAV=0.0
POSTAV=0.0
C DO 8761 KNB=LSWID1,LSWID2
C ——- PREAV=PREAV+(TZ(K-KNB,11)/SLWID)
C8761 POSTAV=POSTAV+(TZ(K+KNB,11)/SLWID)
IF (ABS(POSTAV-PREAV).LT.SLWIN) GO TO 877
CALL SYMBOL(C,-2.,.2,16,0.0,-1)
CALL PLOT(C,B,3)
877 CONTINUE
C PLOT ARRAY OF RMS ENERGY
IF(IPLOT.EQ.1) CALL NEWPEN(1)
KKNT=0
AVENER=0.0
DO 878 L=(NTWICE+1)+(NFIT/2) , JCT-((NTWICE+1)+(NFIT/2))
CALCULATE BAND LIMITED RMS-FIRST SQUARE FUNCTION AND INTEGRATE
INTEGRAL ( (A*F **B )**2)=(A**2/(2*B+1)*F **(2*B+1)
XINTLO=(TZ(L,12)**2/(2.*TZ(L,11)+1.))*FR1**(2.*TZ(L,11)+1. )
XINTHI=(TZ(L,12)**2/(2.*TZ(L,11)+1.))*FR2**(2.*TZ(L, 11) +1. )
SINCE SPECTRUM IS SYMMETRIC ABOUT ZERO, CAN EVALUATE INTEGRAL
BY MULTIPLYING EVALUATED RESULT BY TWO
ENERGY(L)=2.*(XINTHI-XINTLO)
TO CALCULATE MEAN SQUARE, DIVIDE BY WIDTH
ENERGY (L)=ENERGY(L)/(2.*(FR2-FR1) )
NOW TAKE SQUARE ROOT TO DETERMINE RMS
ENERGY (L)=SORT (ENERGY (L))
878 AVENER=AVENER+(ENERGY (L)/( JCT-(4*NUMFIL+(NFIT+1))))
ENSCAL=4.*AVENER ;
IF(IPLOT.NE.1) GO TO 8782
eee AXES(0.0,0.0,10HRMS ENERGY,10,2.,90.,0.5,0.0,AVENER/2.,
ee
CALL PLOT(0.0,0.0,3)
8782 DO 8781 L=(NTWICE+1)+(NFIT/2) , JCT-((NTWICE+1)+(NFIT/2))
an
QO a anaQNn
154
KKNT=KKNT#1
C=(L-1)*.02
B=ENERGY(L)*(1./AVENER)
IF(L.EQ.(NTWICE+1)+(NFIT/2).AND.IPLOT.EQ.1) CALL PLOT(C,B,3)
IF(IPLOT.EQ.1) CALL PLOT(C,B,2)
IF(KKNT.LE.INWID2) GO TO 8781
IF(L.GE.( JCT-( (NTWICE+NFIT/2+1)+INWID2))) GO TO 8781
PREAV=0.0
POSTAV=0.0
DO 8771 KNB=INWID1,INWID2
PREAV=PREAV+(ENERGY (L-KNB) /ENWID)
8771 POSTAV=POSTAV+( ENERGY (L+KNB ) /ENWID)
IF (ABS(POSTAV-PREAV).LT.ENWIN) GO TO 8781
CALL SYMBOL(C,2.,.2,16,180.0,-1)
CALL PLOT(C,B,3)
8781 CONTINUE
DO 879 I=(NTWICE+1)+(NFIT/2) , JCT-( (NTWICE+1)+(NFIT/2) )
RMSSLP=RMSSLP+ABS ((TZ(1I,11)-AVSLP)**2) /( JCT-(4*NUMFIL+(NFIT-1)))
879 RMSENG=RMSENG+ABS ( (ENERGY (I )-AVENER)**2)/( JCT-(4*NUMFIL+(NFIT-1) ))
RMSSLP=SQRT(RMSSLP )
RMSENG=SQRT (RMSENG )
C WRITE(6,'(4(F8.3,2X))') AVSLP,RMSSLP,AVENER,RMSENG
IST=(NTWICE+1 )+(NFIT/2+1)
IEND=JCT-(IST-1)
333 CONTINUE
C PICK PROVINCE BOUNDARIES ,STORE PROVINCE NUMBER FOR EACH POINT IN X
DO 34 I=IST,IEND
X(I)=IFIX(SQRT( ENERGY (1) )/PDEL )+1
34 CONTINUE
C OUTPUT POSITIONS Or PROVINCE BOUNDARIES
C WRITE(06,60) CUT
C 60 FORMAT(' PROV.NO TYPE 1ST LAT LONG LAST LAT LONG
C 1 RMS-FILTERED CUT=',F7.5)
IPNUM=X (IST)
SUM=0.0
JST=IST
KST=1
DO 35 I=IST,IEND
JPNUM=X (I)
SUM=SUM+X (I)
IC=I-JST
IF (JPNUM.EQ.IPNUM.AND.I.LT.IEND) GO TO 35
C IF((IC+1).LT.MIN) GO TO 35
JJ=I
C COMPUTE AVERAGE PROVINCE NO.
AC=IC+1
AVE=SUM/AC
IAVE=AVE
BAVE=IAVE
AVE=AVE-BAVE
IF(AVE.GE.0.5) IAVE=IAVE+1
IPNUM=IAVE
155
K1(KST)=JST
K2(KST)=JJ
K3(KST)=IPNUM
JST=JJ
SUM=0.0
IPNUM=X ( JST )
IF(KST.EQ.1) GO TO 116
MST=KST-1
IF(IAVE.NE.K3(MST)) GO TO 116
K2(MST )=K2(KST)
GO TO 35
116 KST=KST+1
IF(KST.GT.LTPS) GO TO 117
35 CONTINUE
GO TO 118
117 WRITE(06,119)
119 FORMAT(' YOU HAVE TOO MANY PROVINCES-INCREASE DIMENSION OR PDEL')
GO TO 500
C NOW OUTPUT FINAL PROVINCE BOUNDARIES OF AT LEAST MIN. SIZE
118 K3(KST)=99
KST=KST-1
C SMOOTHE SMALL PROVINCE SEGMENTS
eae PRVFIX(K1,K2,K3,MIN,KST)
J2=
DO 36 [=1,KST
II=I+1
IF(K3(II).EQ.K3(I)) GO TO 36
JST=K1(J2)
JJ=K2(1)
IPNUM=K3(J2)
ISDLAT=AY( JST)
A=ISDLAT —
SELAT=ABS (AY ( JST) -A)*60.
ISOLNG=AX (JST)
A=ISDLNG
SELNG=ABS (AX ( JST)-A)*60.
IDLAT=AY (JJ)
A=IDLAT
EELAT=ABS(AY(JJ)-A)*60.
IDLNG=AX (JJ)
A=IDLNG
EMLNG=ABS (AX( JJ)-A)*60.
IDIFF=K2(1)-K1(1)
C COMPUTE RMS LEVEL OF HIGH PASSED DATA FOR THIS PROVINCE
C ALSO AVERAGE SPECTRAL SLOPE AND Y-INTERCEPT
C3=ABS( JJ-JST+1)
RMS=0.0
SLOPEX=0.0
XINTER=0.0
DO 37 IE=JST,JJ
SLOPEX=SLOPEX+(TZ(IE,11)/C3)
XINTER=XINTER+(TZ(IE,12)/C3)
156
37 RMS=RMS+(ENERGY(IE)/C3)
C ADJUST ENERGY LEVEL IN LOG-LOG SPACE AFTER NORMALIZATION
XINTER=10**(ALOG10(XINTER)-.5)
IPNUM=IF IX(SQRT(RMS )/PDEL )+1
IF(ILIST.EQ.1)WRITE(6,51)1, 1PNUM, ISDLAT,SELAT, ISDLNG,SELNG, IDLAT,
1EELAT, IDLNG,EMLNG,RMS , IDIFF , SLOPEX, XINTER
C LOOP TO INSERT DATA INTO FILE FOR LATER FFT RUNS
IF(ISAVE1.NE.1) GO TO 49
ENOMRK=9.999999
NPTS=JJ-JST
WRITE(13,46) FILE(IKONT),I ,SLOPEX,XINTER,GRID,NPTS
WRITE (13,50)1, IPNUM, ISDLAT, SELAT, ISDLNG, SELNG, IDLAT ,EELAT,
LIDLNG, EMLNG , RMS
WRITE(13, 47)(AZ(1Q) ,1Q=9ST, Jd)
XRMS= SRMS (AZ( JST) , NPTS )
46 FORMAT(A20,12,1X, F8. 351X693 a XeF 74), 15)
47 FORMAT ( 10F8. 6)
WRITE(13,47) ENDMRK
49 IF (IPLOT.NE.1) GO TO 248
APNUM=.5+IPNUM*.5
AJ=JST-1
BJ=JJ-1
C1=AJ*.02
C2=BJ*.02
CALL PLOT(C1,APNUM, 3)
CALL PLOT(C2,APNUM, 2)
50 FORMAT(I5,18,4(1X,14,F6.2),F10.4)
51 FORMAT(15,18,4(1X,14,F6.2),F10.4,16,2F10.6)
248 J2sII
36 CONTINUE
880 IF(IKONT.EQ.JTOTL) GO TO 500
C REPOSITION PLOTTER FOR NEXT PLOT
IF (MOD(IKONT,4).NE.0) GO TO 882
C2MAX=C2MAX+15 .0
IF(IPLOT.EQ.1) CALL PLOT(C2MAX,-44.0, -3)
C2MAX=0.0
GO TO 440
882 IF(IPLOT.EQ.1) CALL PLOT(0.0,8.0,-3)
GO TO 440
500 IF(IPLOT.EQ.1) CALL PLOT(0.0,0.0,999)
IF(ISAVE1.GT.0) CLOSE(UNIT= 13)
STOP
END
157
READER SUBROUTINE FOR PROVPICK
SUBROUTINE PROVRD(NP,XLAT ,XLONG, DEPTH, IKONT FILE)
DIMENSION XLAT (3000) ,XLONG( 3000) , DEPTH ( 3000)
CHARACTER*20 FILEN
CHARACTER*20 FILE(20)
KNTER=0
400 FILEN=FILE(IKONT )
WRITE (6,433 )FILEN
433 FORMAT(' OPENING FILE ',A20)
OPEN(UNIT=12,FILE=FILEN,STATUS='OLD' ,FORM='FORMATTED' )
DO 420 I=1,NP
438 READ(12,434,END=470) XLAT(1I),XLONG(I),IDEPTH
434 FORMAT (9X ,F12.8,F13.8,17)
KNTER=I
IF(IDEPTH.LE.0) PRINT *, "DATA POINT SKIPPED'
IF (IDEPTH.LE.0) GO TO 438
DEPTH(I )=IDEPTH*. 0018288
420 IF(I.LT.30)WRITE(6,*)XLAT(I) ,XLONG(I) ,DEPTH(T)
GO TO 490
470 NP=KNTER
490 CLOSE (UNIT=12)
WRITE (6,492 )NP
492 FORMAT(' NUMBER OF POINTS READ =',15)
RETURN
END
158
SUBROUTINE MAPCTN(GRID,X,Y,Z,N,AX,AY,AZ, ICT)
ROUTINE TO MAP DATA VALUES(Z) WITH ASSOCIATED POSITIONS
(X=LONG DECIMAL DEG.,Y=LAT DECIMAL DEG. )FROM A RELATIVELY STRAIGHT
SEGMENT OF SURVEY TRACK ONTO A STRAIGHT LINE,ADJUST THE (Z) VALUE FOR
THE AMOUNT OF SHIFT REQUIRED BY THE MAPPING AND INTERPOLATE NEW (Z)
VALUES (AZ) AT AN EQUALLY SPACED DISTANCE(GRID IN DECIMAL NAUTICAL MI.
AND ASSOCIATED POSITIONS AX,AY ALONG THIS STRAIGHT LINE.
C** N=NO.OF ORIGINAL X,Y,Z INPUT PTS., ICT= NO.OF OUTPUT PTS.
C** NOTE-SINCE THIS IS A CARTESIAN MAP WHICH DOES NOT ACCOUNT FOR LONG
C CONVERGENCE ,IT SHOULD NOT BE USED FOR SEGMENTS COVERING A LARGE
C LATITUDE RANGE.
C** NOTE-ORIGINAL XYZ DATA IS DESTROYED IN THIS ROUTINE
C*** THIS ROUTINE CALLS GINT, SPLINE ,SPLICO,SORTY
DIMENSION X(7000),Y (7000) ,Z( 7000) ,AX (7000) ,AY( 7000) ,AZ( 7000),
*DIST (7000)
00 55 I=1,7000
55 DIST(I)=0.0
JOIR=1
IF (ABS (X(N)-X(1)).GE.ABS(Y(N)-¥(1))) JDIR=0
ATER=9999.99
RAD=0. 00029089
BLONG=X (1)*60.0
BLAT=Y(1)*60.0
DO 1 I=1,N
IF (JDIR.EQ.1)GO0 TO 3
X(1)=BLONG-X (1)*60.0
Y(1)=¥(1)*60.0-BLAT
GO TO 1
3 AY(1I)=BLONG-X (1)*60.0
X(1)=¥(1)*60.0-BLAT
Y(1I)=AY(T)
1 CONTINUE
AN=N
C MAP INPUT POSITIONS ONTO STRAIGHT LINE TO PREPARE DATA FOR
C INTERPOLATION ON AN EQUAL DISTANCE BASIS
ANANNQARM
aAnxy> 8 oOOWn>
ow ow
5 D=0+¥(1)*xX(1)
Al=(C*B-D*A) / (AN*B-A**2)
A2=(C-A1*AN)/A —
IF (ABS (A2).LT.0. 000001 )A2=0. 000001
LEAST SQUARES LINE IS Y=A1+A2X
NOW SEARCH FOR A POINT LESS THAN PIVOT DISTANCE FROM TRACK LINE
TO USE FOR 1ST PIVOT AND MAP PTS.ONTO TRACK. WITH CORRECT Z VALUE
C POINTS-LT PIVOT DIST FROM TRACK WILL HAVE THEIR POSITS MAPPED ONTO
C We aees BUT THEIR Z VALUE WILL NOT BE CHANGED
UM=0.0
(owe w)
159
DO 6 I=1,N
DIST(1)=ABS((¥(1)-A2*X(1)-A1)/ (SQRT(A2**2+1.0)))
6 SUM=SUM+DIST(1)**2
PIVOT =SQRT (SUM/AN )
C REMOVE POINTS THAT MAY HAVE A BAD POSITION
XP IVOT=3. O*PIVOT
J=0
DO 42 I=1,N
IF (DIST(1).GT.XPIVOT) GO TO 42
J=J+1
X(J)=X(T)
Y(J)=¥(1)
Z(J)=Z(1)
DIST(J)=DIST(I)
42 CONTINUE
N=J
AN=N
00 7 I=1,N
IF (DIST(I)- PIVOT) 8,8,7
7 CONTINUE
8 A3= -1.0/A2
AX (1 )=(A1+A3*X (1 )-¥(1))/(A3-A2)
AY(I)= A2*AX(1) +Al
AZ(1)=Z(1)
IA=I+1
C NOW WORK BACKWARDS ON TRACK TO PICK UP POINTS THAT FAILED
C PIVOT TEST
11 J=I-1
LEO) 1254259
9 DELZ=(Z(J)-Z(1))/SQRT((X(1)-X (J) )**2 +(¥(1)-¥ (0) )**2)
AX (J )=(A1+A3*X (J )-¥(d)) /(A3-A2)
AY (J )=A2*AX (J)+A1
AZ (J )=(DELZ*SQRT ( (AX (1 )-AX (J) )**2+ (AY (1)-AY(J))**2))+AZ(T)
I=J
GO TO ll
C NOW WORK FORWARD ON TRACK TO PICK UP REMAINING PTS.
12 DO 13 I=IA,N
IF (DIST(I)- PIVOT) 14,14,15
14 AX(1)=(A1+A3*X (1 )-¥(1))/(A3-A2)
AY(T)= A2*AX(I) +Al
AZ(1)=Z(T)
GO TO 13
15 J= I-1
DELZ=(Z(1)-Z(9))/SQRT((X(1)-X(J))**2 +(¥(1)-¥(0))**2)
AX (1 )=(A1+A3*X (1)-¥(1))/(A3-A2)
AY(1)=A2*AX(I) +Al
AZ (1 )=(DELZ*SQRT ( (AX (1 )-AX (J) )**2+ (AY (1 )-AY (J) )**2) )+AZ (0)
13 CONTINUE
SORT INPUT SO THAT THE INDEPENDENT VARIABLE IS MONOTONIC AND REMOVE
CLOSELY SPACED POINTS TO CONTROL SPLINE INTERPOLATION
DAVE =0.0
DO 44 [=2,N
J=I-1
an
160
an (ok Sr) AND
am
44 DAVE=DAVE+X (1)-X (J)
DAVE =0. 3* (ABS (DAVE /AN ) )
CALL SORTY(AY,AX,AZ,Y,X,Z,N,1,1,DAVE )
IF(JDIR.LT.1) GO TO 92
DO 93 I=1,N
Z(1)=AZ(T)
X(1)=AY(1)
93 Y(1)=AX(1)
0
92 DO 95
94 CONTINUE
RESET ORIGIN TO FIRST MAPPED POINT
BLAT =BLAT +Y¥ (1)
BLONG=BLONG-X (1)
DX=X (1)
DY=Y(1)
DO 71 I=1,N
X(1)=X (1 )-DX
71 Y(1)=¥(1)-DY
NOW INTERPOLATE ALONG TRACK AT DESIRED GRID SPACING
COMPUTE APPROX.LENGTH OF LONGITUDE FOR THIS TRACK
USE APPROX. MIDLATITUDE AS BASIS AND TABLE 6(BOWDITCH)
NN=N/2
AL =(Y(NN )+BLAT )*RAD
A=(111415.13*COS (AL )-94.55*COS (3. *AL )+.012*COS(5.*AL ))/60.0
UNITS OF A ARE METERS/MINUTE OF LONGITUDE
COMPUTE APPROX.NAUTICAL MILES/MINUTE OF LONGITUDE
19 B=A*(1.0/1852.0)
CONVERT X COORDINATE OF MAPPED POSITION TO MILES AND COMPUTE DISTANCE
DOWN TRACK.
00 25 I=1,N
25 DIST(1)=SQRT ((X(1)*B)**2+Y (1 )**2)
IF (ABS (X(N)).LT.0.00001) X(N)=0.00001
THETA=ATAN2(Y(N), (X(N)*B) )
CALL GINT(GRID,N,0.0,DIST(N), ICT, DIST,Z,AX,AZ)
STORE NEW POSIT OF INTERPOLATED PT.IN AX(LONG,AY(LAT) DECIMAL DEG.
NEW INTERPOLATED VALUE OF Z IS STORED IN AZ
DO 26 I=1,ICT
AY (1T)=(AX (T)*SIN (THETA )+BLAT )/60.0
26 AX (I )=(BLONG-(AX (I )*COS(THETA)/B) )/60.0
RETURN
END
161
C FUNCTION TO CONVERT LATITUDE INTO MERIDIONAL PARTS:
FUNCTION YMP(Z)
DATA AP/0.7853981634/
Y=ABS (Z)*0. 290888209E -03
T=TAN (AP+Y*0. 5)
YM=7915. 7045*ALOG10(T )-23. 268932*SIN(Y)
YMP=YM*SIGN(1.0,Z)
RETURN
END
162
SUBROUTINE SORTY(X,Y,Z,AX,AY,AZ,K,KODE , JCODE ,GRID)
C Y=INPUT VARIABLE TO BE SORTED, X,Z=VALUES ASSOCIATED WITH Y
C K=LENGTH OF Y,IF KODE=1,VALUES OF Y WHICH ARE WITHIN 1 GRID INT OF
C PREVIOUS VALUE ARE REMOVED,IF KODE=0 ALL VALUES OF Y ARE
C RETAINED, IF JCODE=+1,Y IS SORTED IN INCREASING ORDER,IF JCODE=-1,
C Y IS SORTED IN DECREASING ORDER,OUTPUT IS SORTED VALUES OF
C Y WITH ASSOCIATED X AND Z
DIMENSION Y¥(20),X(20),Z(20),AX (20) ,AY(20) ,AZ(20)
KB = K
CODE =JCODE
Jzl
1295 alsa el
JCT=0
AY(J) = Y(T)
132 TEMP= CODE*(AY(J)-Y(I+1))
IF ((ABS (TEMP) )-GRID+0.0001) 122,120,120
TZO SIR GUEME) lal 36; 123
2 eles a te
Wao 6 dh)) Ss) obeys ite y2e ac)
123) eh =) 2
AY(J) =
AX(J) =
AZ(J) =
RU a ee
122 IF(KODE) 120,120,136
136 KD=1+2
IF (KD-KB) 124,124,139
139 K=K-1
KB=KB-1
GO TO 125
124 DO 126 JD=KD,KB
JF = J0-1
Y(JD)
X(JD)
Z(J0)
nn oO
Ws) (ule)
125 IF(JCT) 127,127,128
Weil Nid) YG)
AX(J) = X(1)
aa) 2 OU)
KT = 2
128) dander
IF(J - K) 131,133,133
131 DO 134 KA = KT,KB
JT = KA -
Y(JT)
X(JT)
IY. Aaa)
KB = KB - 1
GO To 129
“won ow
~<
=
ras
>
~
163
133 IF(JCT) 137,137,138
137 KB=KB+1
138 AY(K)
135
AX (K)
AZ(K)
DO 135
Y(T)
X(1)
Z(1)
RETUR
END
auunu
Y(KB- 1)
X(KB- 1)
Z(KB- 1)
I = 1,K
AY(1)
AX(1)
AZ(1)
164
SUBROUTINE GINT(DELX,M,XBGN,XEND,ICT,X,Y,AX,AY)
C MODIFICATION OF ORIGIONAL GINT SO THAT INPUT DATA IS NOT DESTROYED
C GENERAL 1-D SPLINE INTERPOLATION FOR MIN.STORAGE
C INTERPOLATES OVER 50 INPUT PTS WITH OVERLAP
C DELX=DESIRED INTERPOLATION INTERVAL,X AND Y ARE INPUT WITH X=INDEP.
C VARIABLE,AX AND AY ARE OUTPUT,M=LENGTH OF X, XBGN AND XEND=DESIRED
C BEGINNING AND ENDING VALUES OF AX,ICT=LENGTH OF AX
C***NOTE IF M MOD 50 IS LESS THAN 2 INTERPOLATED OUTPUT MAY STOP SHORT
C OF XEND
C THIS ROUTINE REQUIRES SPLINE AND SPLICO SUBROUTINES
DIMENSION X(1),Y(1),AX (1) ,AY(1)
ICT=(ABS (XEND-XBGN ) /DELX )+1.0
ISECT=M/50
IA=ISECT*50
IC=0
JJ=0
MM=49
IF(IA.EQ.0) MM=M
IF((M-IA).GE.2) JJ=1
C OVERLAP INTERPOLATION INTERVALS BY 2 INPUT PTS SO SPLINE ROUTINE
C IS DIMENSIONED TO MAX OF 53 PTS
JCONT=1
KA=1
KC=MM+1
IF (IA.EQ.0) KC=KC-1
K=1
53 IKT= (ABS (X(MM)-XBGN)/DELX )+1.0
IF (MM.EQ.M) IKT=ICT
ATER=9999.999
DO 26 I=KA,IKT
AJ=1-1
X INT =Ad*DELX+X BGN
IF (XINT.GT.XEND) GO TO 58
CALL SPLINE (X(K),Y(K) ,KC,XINT, YINT ATER)
AX (1) =X INT
26 AY(I)=YINT
JCONT =JCONT +1
IF(IA.EQ.0) GO TO 56
IF(JCONT.GT.ISECT) GO TO 55
JC=JCONT*50
57 IMM=MM-1
K=IMM
MM=MM+50
IF(IC.EQ.1) MM=M
KA=IKT+1
KC=JC-IMM+1
GO TO 53
55 IF(JJ.LT.1) GO TO 56
IC=1
JJ=0
JC=M
GO TO 57
58 ICT=I-1
165
RETURN
56 ICT=IKT
RETURN
END
166
SUBROUTINE SPLINE (X,Y,M,XINT,YINT,ATER)
SEE PENNINGTON REF. FOR DESCRIPTION OF THIS SUBROUTINE
DIMENSION X(1),Y(1),C(4,53)
K=0
IF (X(1)+Y¥ (M)+Y¥ (M-1)+X (M-1)+Y (M-2)-ATER) 10,3,10
10 CALL SPLICO (X,Y,M,C)
ATER= X(1)+Y¥(M)+Y(M-1)+X (M-1)+Y (M-2)
K=1
3 IF(ABS(XINT-X(1)).LT.0.00001) GO TO 1
IF (XINT-X(1)) 70,1,2
70 K=1
GO TO 7
1 YINT=¥(1)
RETURN
2 IF (ABS(XINT-X(K+1)).LT.0.00001) GO TO 4
IF (XINT-X(K+1))6,4,5
4 YINT=Y(K+1)
RETURN
5 K=K+1
IF(M-K) 71,71,3
71 K=M-1
GO TO 7
6 IF (ABS(XINT-X(K)).LT.0.00001) GO TO 12
IF (XINT-X(K))13,12,11
12 YINT=Y(K)
RETURN
13 K=K-1
GO TO 6
11 YINT=(X(K+1)-XINT )*(C(1,K)*(X(K+1)-X INT )**2+C(3,K) )
YINT=YINT+(XINT-X (K) )*(C(2,K)*(XINT-X (K) )**2+C(4,K))
RETURN
7 PRINT 101,XINT
7 CONTINUE
101 FORMAT(' CAUTION VALUE AT POSITION',F10.2,'WAS EXTRAPOLATED' )
GO TO ll
RETURN
END
167
Ww P
=~
toy)
™N
SUBROUTINE SPLICO (X,Y,M,C)
DIMENSION X(1),¥(1),C(4,53) ,D(53) ,P(53) ,E{53),A(53,3),B(53) ,Z(53)
MM=M-1
DO 2 K=1,MM
D(K) =X (K+1)-X (K)
P(K)=D(K)/6.
E(K)=(Y(K+1)-¥(K))/D(K)
00 3 K=2,MM
B(K)=E(K)-E(K-1)
A(1,2)=-1.-D(1)/D(2)
A(1,3)=0(1)/D(2)
A(2,3)=P(2)-P(1)*A(1,3)
A(2,2)=2.*(P(1)+P(2))-P(1)*A(1,2)
A(2,3)=A(2,3)/A(2,2)
B(2)=B(2)/A(2, 2)
DO 4 K=3,MM
A(K,2)=2.* (P(K-1)+P(K) )-P(K-1)*A(K-1, 3)
B(K)=B(K)-P(K-1)*B(K-1)
A(K,3)=P(K) /A(K, 2)
B(K)=B(K) /A(K, 2)
Q=D(M-2)/D(M-1)
A(M,1)=1.+0+A(M-2, 3)
A(M,2)=-Q-A(M, 1)*A(M-1, 3)
B(M)=B(M-2)-A(M, 1)*B(M-1)
Z(M)=B(M)/A(M, 2)
MN =M-2
DO 6 I=1,MN
K=M-I
Z(K)=B(K)-A(K, 3)*Z(K+1)
Z(1)=-A(1,2)*Z(2)-A(1,3)*Z(3)
DO 7 K=1,MM
Q=1./(6.*D(K))
C(1,K)=Z(K)*Q
C(2,K)=Z(K+1)*Q
C(3,K)=¥(K) /D(K)-Z(K)*P(K)
C(4,K)=¥(K+1)/D(K)-Z(K+1)*P(K)
RETURN
END
168
SUBROUTINE AXES(X,Y,IBCD,NC,AXLEN,ANG,DELTIC ,FIRSTV,DELVAL ,NDEC)
C MODIFIED CALCOMP AXIS SUBROUTINE ---RANKIN,NOV.1971
C X,Y COORDINATES OF STARTING POINT OF AXIS IN INCHES
C IBCD AXIS TITLE
C NC NUMBER OF CHARACTERS IN TITLE
C IF NC IS(-),HEADING AND TICKS ARE BELOW THE AXIS
C AXLEN FLOATING POINT AXIS LENGTH IN INCHES
C ANG ANGLE OF AXIS FROM HORIZONTAL IN DEGREES
C DELTIC DISTANCE BETWEEN TIC MARKS IN INCHES
C FIRSTV SCALE VALUE AT FIRST TIC MARK
C DEL VAL SCALE INCREMENT
C NDEC NUMBER OF DECIMAL PLACES OF TIC ANNOTATION PLOTTED(PUNCH
C -1 IF ONLY INTEGER(NO DECIMAL POINT)IS DESIRED)
C -2 IF NO ANNOTATION DESIRED IE. ONLY TICKS
DIMENSION IBCD(10)
A=1.0
KN=NC
IF(NC)1,2,2
1 A=-A
KN=-NC
2 XVAL=FIRSTV
STH=ANG*0. 0174533
CTH=COS (STH )
STH=SIN(STH)
OXB=-0.1
DYB=0.15*A-0.05
XN=X+DXB*CTH-DYB*STH
YN=Y+DYB*CTH+DXB*STH
NT IC=AXLEN/DELTIC+1.0
NT=NTIC/2
00 10 I=1,NTIC
IF(NDEC.EQ.-2) GO TO 12
C CHANGED NUMBER HEIGHT FROM .105 TO .15
CALL NUMBER(XN,YN,0.15,XVAL ,ANG,NDEC )
12 XVAL=XVAL+DELVAL
XN=XN+CTH*DELTIC
YN=YN+STH*DELTIC
IF (NT )10,11,10
11 Z=KN
OXB=-0.07*Z+AXLEN*0.5
DYB=0. 325*A-0.075
XT =X+DXB*CTH-DYB*STH
YT =Y+DYB*CTH+DXB*STH
C CHANGED HEIGHT FROM .14 TO .18
CALL SYMBOL (XT,YT,0.18,1BCD(1) ,ANG,KN)
10 NT=NT-1
CALL PLOT (X+AXLEN*CTH, Y+AXLEN*STH, 3)
IF(NDEC.EQ.-2) GO TO 14
OXB=-0.07*A*STH
DYB=0.07*A*CTH
GO TO 13
14 DXB=-0.05*A*STH
DYB=0.05*A*CTH
169
13
20
10
20
A=NTIC-1
XN=X +A*CTH*DELTIC
YN=Y+A*STH*DELTIC
00 20 [=1,NTIC
CALL PLOT(XN, YN, 2)
CALL PLOT(XN+DXB, YN+DYB, 2)
CALL PLOT(XN,YN, 2)
XN=XN-CTH*DELTIC
YN=YN-STH*DELTIC
CONT INUE
RETURN
END
FUNCTION SRMS(X,N)
DIMENSION X(1)
AVE=0.0
SRMS=0.0
DO 10 I=1,N
AVE =AVE +(X(1)/FLOAT(N))
00 20 I=1,N
SRMS =SRMS + ((X(1)-AVE )**2)
SRMS=SRMS /FLOAT (N )
SRMS =SQRT (SRMS )
RETURN
END
170
SUBROUTINE FILTER(X,NP,CUT,H,N,K,WGT )
C GENERAL ROUTINE TO HIGH OR LOW PASS A SET OF EQUALLY
C SPACED DATA USING MARTIN FILTERS.
C X=INPUT DATA AND OUTPUT DATA,NP=NO.OF POINTS IN X,CUT=NORMALIZED
C CUTOFF OF FILTER IN CYCLES/DATA INTERVAL ,H=SLOPE PARAMETER,
CN=HALF LENGTH OF FILTER,TOTAL LENGTH=2N+1,
C K=1=HIGH PASS,=0 FOR LOW PASS.
C WEIGHTS STORED IN WGT
DIMENSION X(1),WGT(1)
NA=N+1
WGT (1)=2.0* (CUT+H)
C CENTER WEIGHT STORED IN LOCATION 1
SUM=0.0
PI=3. 1415926
DO 61 I=2,NA
P=I-1
Q=1.0-16. O*H**2*P**2
IF (ABS(Q).GT.0.0001) GO TO 62
WGT (I )=SIN(2.*PI*P*(CUT+H))/(4. O*P)
GO TO 61
62 WGT(1)=((COS(2.*PI*P*H))/Q)*((SIN(2.*PI*P*(CUT+H)) )/(PI*P))
61 SUM=SUM+WGT (I)
DELTA=1.-(WGT(1)+2.*SUM)
AX=2*N4+1
IF(K.LT.1) GO TO 78
DO 65 I=2,NA
65 WGT(1)=(WGT (1)+DELTA/AX )*(-1.0)
WGT (1)=1.0-(WGT (1)+DELTA/AX )
GO TO 79
78 DO 80 I=1,NA
80 WGT(I)=WGT(I)+DELTA/AX
79 NB=NP-N
C CONVOLVE WEIGHTS WITH DATA.
DO 63 I=NA,NB
II =I+1-NA
SUM=0.0
DO 64 J=1,NA
Jl=I+J-1
J2=I-J+1
64 SUM=SUM+WGT (J)*(X(J1)+X (J2))
63 X(II)=SUM-WGT (1)*X (I)
C SHIFT FILTERED DATA TO CORRECT LOCATION AND ZERO ENDS.
II =NB+1-NA
DO 67 I=1,I1
J=II+1-1
67 X(L)=X(
I
68 X(I)=0.
DO 69 I=NC,NP
69 X(I)=0.0
RETURN
END
171
AANANQNND
AANDARMANAADNNNMA Ce)
aamM oO
10
100
110
120
200
‘250
260
SUBROUTINE ENVEL (DATA,NP)
THIS ROUTINE RECEIVES A TIME SERIES OF LENGTH NP -
IN ARRAY DATA, AND RETURNS ENVELOPE ESTIMATES BASED
ON A HILBERT TRANSFORM METHOD AS DESCRIBED IN KANASEWICH
(1981) ,P. 362-368.
C.G.FOX ,ADVANCED TECHNOLOGY STAFF ,NAVOCEANO-5/12/83
DIMENSION DATA(1)
COMPLEX XDATA(2048)
NSTART=1
NPSEG=NP
IF (NPSEG.LE.0) RETURN
PERFORM ENVELOPE CALCULATIONS IN SEGMENTS OF 2048 POINTS
IF (NPSEG.GT.2048) THEN
NPSEG=NPSEG-2048
LENGTH=2048
GO TO 100
ELSE
LENGTH=NPSEG
NPSEG=0
END IF
TRANSFER INPUT DATA TO REAL PORTION OF COMPLEX ARRAY XDATA
00 110 I=1,LENGTH
XDATA(I )=(1.0,0.0)*DATA(NSTART+(I-1))
ZERO FILL ARRAY IF LESS THAN 2048 POINTS
IF (LENGTH.EQ. 2048) GO TO 200
DO 120 I=LENGTH+1,2048
XDATA(T )=(0.0,0.0)
FOURIER TRANSFORM TO FREQUENCY DOMAIN USING FFT
CALL NLOGN(11,XDATA,-1.0)
COMPUTE HILBERT TRANSFORM IN THE FREQUENCY DOMAIN. THE
TRANSFORM IS COMPUTED AS I*SGN(OMEGA)*F (OMEGA), WHERE
OME GA=FREQUENCY F(OMEGA)=FOURIER TRANSFORM
SGN (X )=1, (X>0), =0, (X=0) , =-1(X<0)
THIS OPERATION INDUCES A NINETY DEGREE PHASE SHIFT ON
COMPLEX PLANE. THE CALCULATION IS DONE BY ZEROING X(1),
(OMEGA=0) ,TRANSFERRING REAL TO IMAGINARY,MINUS IMAG TO REAL
FOR X(2)-X(N/2),(OMEGA>O), AND TRANSFERRING MINUS REAL TO
IMAGINARY,AND IMAGINARY TO REAL FOR X(N/2+1) TO X(N),
(OMEGA<O )
XDATA(1)=(0.0,0.0)
DO 250 [=2,1024
TEMPR=REAL (XDATA(I))
TEMP I =AIMAG(XDATA(1))
XDATA(I )=CMPLX (-TEMPI, TEMPR )
DO 260 [=1025,2048
TEMPR=REAL (XDATA(T ))
TEMPI=AIMAG(XDATA(I))
XDATA(I )=CMPLX (TEMPI , -TEMPR )
INVERSE TRANSFORM SHIFTED DATA
CALL NLOGN(11,XDATA,+1.0)
CALCULATE ENVELOPE AFTER EQUATION 21.3-1 OF KANASEWICH
DUE TO SYMMETRIES IN THE TRANSFORM, IMAGINARY PORTION
OF XDATA IS NEAR ZERO AND CAN BE IGNORED
172
DO 300 I=1,2048
N=NSTART+(I-1)
300 _DATA(N)=SQRT ( (DATA(N)*DATA(N ) )+(REAL (XDATA(I))**2))
NSTART=NSTART +2048
GO TO 10
END
173
SUBROUTINE NLOGN(N,X,SIGN)
C*****NLOGN COMPUTES THE DISCRETE FOURIER TRANSFORM BY THE FAST FOURIER
C*****TRANSFORM METHOD.REFERENCE ROBINSON,PAGE 63.
C*****xTHIS SUBROUTINE CONTAINS 51 STATEMENTS.
C*****N=POSITIVE INTEGER FOR THE POWER OF 2 DESIRED.
C*****X=INPUT AND OUTPUT DATA ARRAY (COMPLEX ).
C*****STGN=-1.0 FOR FOURIER TRANSFORM,+1.0 FOR INVERSE FOURIER TRANSFORM.
DIMENSION M(15),X(2)
COMPLEX X,WK,HOLD,Q
EX =2e=N
DO 11 =1,N
TMC) = 2** (NT)
00 4L = 1,N
NBLOCK = 2**(L -1)
LBLOCK = LX / NBLOCK
LBHALF = LBLOCK / 2
K =0
DO 4 IBLOCK = 1,NBLOCK
FK =K
FLX = LX
V = SIGN * 6.2831853071796 * FK / FLX
WK = CMPLX(COS(V),SIN(V))
ISTART = LBLOCK * (IBLOCK - 1)
00 2 I = 1,LBHALF
J = ISTART + I
JH = J + LBHALF
Q = X(JH) * WK
X(JH) = X(J) - Q
X(J) = X(J) + Q
2 CONTINUE
00 3 I = 2,N
II =I
IF(K.LT.M(I)) GO TO 4
3K =K - MI)
4K =K + MII)
K = 0
DON Jes a leniex
IF (K.LT.J) GO TO 5
HOLD = X(J)
X(J) = X(K + 1)
X(K + 1) = HOLD
5 00 6 I =1,N
II =I
IF(K.LT.M(I)) GO TO 7
6 K =K - M(I)
7X =K + M(II)
IF (SIGN.LT.0.0) RETURN
DO 8 I = 1,LX
8 X(I) = X(I) / FLX
RETURN
END
ANNAAMNAND
SUBROUTINE PRVFIX(K1,K2,K3,MIN,KST )
DIMENSION K1(1000) ,K2(1000) ,K3(1000) ,11(200) ,T2(200) ,13(200)
SUBROUTINE FOR USE WITH THE PROVINCE PICKER PROGRAM
SEARCHES PROVINCE BOUNDARIES FOR SEGMENTS LESS THAN MIN
SEGMENTS ARE INCORPORATED INTO ADJACENT PROVINCES IF THEY
ARE EQUAL-SETS BOUNDARY TO INCLUDE SEGMENT IN ROUGHEST
PROVINCE IF ADJACENT PROVINCES ARE UNEQUAL-DELETES
SEGMENTS AT ENDS OF TRACKLINE
C.G.FOX-LDGO-JULY 6,1982
MINSAV=MIN
MIN=0
10 ISKIP=0
ICYCLE=0
MIN=MIN+5
IF (MIN.GT .MINSAV )MIN=MINSAV
DO 100 I=1,KST
IE GLEYCRESGE=)G0) 10/195
IX=I-ISKIP
Il =I+l
TEST FOR ADJACENT PROVINCES OF SAME NUMBER
IF(K3(1).NE.K3(II)) GO TO 30
T1(IX )=K1(1)
T2(IX )=K2(1T)
T3(IX )=K3(1)
ICYCLE=1
ISKIP=ISKIP+1
GO TO 100
TEST FOR MINIMUM PROVINCE SIZE
30 IDIFF=K2(I )-K1(1)
IF(IDIFF.LT.MIN)GO TO 40
NORMAL PROVINCE - STORE IN T ARRAY AND CONTINUE
T1(IX )=K1(1)
T2(IX )=K2(1)
T3(IX )=K3(1)
GO TO 100
LESS THAN MINIMUM SIZE-TEST NUMBER OF ADJACENT PROVINCES
IF FIRST OR LAST PROVINCE OF LINE, DELETE
40 IF(I.NE.KST.OR.IX.NE.1) GO TO 50
ISKIP=ISKIP+1
GO TO 100
TEST IF SURROUNDING PROVINCES ARE EQUAL
50 IF(T3(1X-1).NE.K3(II)) GO TO 70
IF EQUAL, EXTEND SECOND BOUNDARY OF PREVIOUS PROVINCE
TO END OF FOLLOWING PROVINCE
T2(1X-1)=K2(11)
ICYCLE=1
ISKIP=ISKIP+2
GO TO 100
ADJACENT PROVINCES ARE DIFFERENT-SET BOUNDARY CLOSER TO
LOWER VALUED PROVINCE
70 IF(T3(1X-1).LT.K3(11)) GO TO 75
T1(1X )=K1(T)
175
AANAARAARAANRANRANM
—S
on
|
oO
80
95
100
150
500
T2(1X )=K2(IT)
T3(1X )=K3(1)
IF ((K2(11)-K1(11)).GT.IDIFF )T3(1X )=K3(11)
GO TO 80
T2(1X-1)=K2(1)
T1(IX )=K1 (IT)
T2(IX )=K2(11)
T3(1X )=K3(1)
IF ((K2(1I1)-K1(11)).GT.IDIFF )T3(1X )=K3(11)
BISECT SEGMENT AND PUT HALVES INTO ADJACENT PROVINCES
K=(K2(1)+K1(1))/2
T2(1X-1)=K
T1(IX )=K
T2(1X)=K2(11)
T3(IX )=K3(1T)
ISKIP=ISKIP+1
ICYCLE=1
GO TO 100
ICYCLE=ICYCLE -1
CONTINUE
IF (ISKIP.EQ.0.AND.MIN.EQ.MINSAV) GO TO 500
KST=KST-ISKIP
DO 150 I=1,KST
K1(I)=T1(1)
K2(T)=T2(1)
K3(I )=T3(1)
K3(KST +1 )=99
G0 TO 10
RETURN
END
176
FUNCTION DECDEG(I,X )
C CONVERT DEGREES+MINUTES TO DECIMAL DEGREES
DECDEG=ABS(I)+X/60.
DECDEG=SIGN(DECDEG,1)
RETURN
END
SUBROUTINE DEGMIN(A,J,B)
C CONVERT DECIMAL DEGREES(A) TO DEGREES(J)+MINUTES(B)
J=INT(A)
B=ABS ( (A-FLOAT (J))*60. )
RETURN
END
177
oy a 4 digi’ -
ah i fh }
age
a.t ee sei pase
“s ENV CE mi
ath) TAM ne (AYE
Lae TEL fon nan) ns
Appendix C.1
The addition of two sinusoids of different amplitude and phase, but
identical wavelength results in another sinusoid of the same wavelength,
but amplitude and phase which is a function of the amplitudes and phases
of the component sinusoids. Assuming two sinusoids of differing ampli-
tudes (vy and Vp) and phases (84 and 9,), their sum can be expressed in
terms of a new sinusoid of amplitude v, and phase 96 as
vi cos(6-8,) + vps cos(8-8,) = ie ° cos(9-8,) (Cc-1)
Using the addition formula of sinusoids,
VN cos9 ° cos8 , Tvia. sind ° sind, + VR’ cosd ° coso, +
VR’ sinO « sino, = Vone cos9 ° cos8, TaVon sind ° sin, (C=-2)
Combining terms, and equating sine and cosine terms
(Jr, cos8 , +v,° cos® 5) * cos = (vo . cos6 ,) * cosd (C=3)
and
(vy ° sind, + vg ° sin?) ¢ sind = (ve ° sin®¢) ¢ sind (c=-4)
Since C-3 and C-4 are true for any value of 9, the unknown v, and 9,
¢c
are determined from:
178
v, ° cos6, + vp ° cos®g = vg * cos8> (c-5)
Va ° sind, + vg ° sind, = Vo° sin8> (C-6)
Squaring both sides and summing equations C-5 and C-6 removes
the 9, terms (cos26 + sin*0 = 1);
23 vac ° cos”6 , + vR- ° cos “0, +2e Via 2 Vp? cos6 , e cos6, +
2
MC
2 e 2 e 2 e e e e -
sin*9, + vp sin“9, + 2° v, * vg * sind, * sind, (C-7)
VN A
Again using the relationship cos26 + sin6 = 1, the addition
formula for cosines, and taking the root of the result gives the
formula for the Vc term as,
1
Vo = (v4? + vB aa ONAN SANGO (cos(8,-0,))/2 (C=8)
The unknown resulting phase angle 9¢ can be derived by dividing
(c-6) by (C-5),
Vince sind, + VR ° sin,
tano, = (C-9)
Vino? sin® , + VR ° cos,
or
2 Va ° sin® , + VR ° sin?,
Q, = tan ( ) (C-10)
Va ° cos8 , + vp ° cos8,
179
These general results show the dependence of a sinusoid of given
wavelength on a linear combination of two other sinusoids of the same
wavelength. The same combination process could be extended to any number
of component sinusoids.
Appendix C.2 contains FORTRAN -77 software to perform an iterative
fit of a sinusoid to data.
180
a rts ‘
wn. aif te aero aap out
iy ‘" + WA re 4 ain tty, ae) ary MMi t pane? 5. ’ sy tony
i aii ded Bia Re eae Siatncetp eee + “tn “ Ay bait di tia
De pilin wad tag ith enue pf ey Henne, pee ‘any
a tot hy
‘ti, | <n mrih bu pins em aia Sy te had je yt awed ; é4 vidtay!
Thon wi |
AAANDAAAAARANRAAANRANRNNANAINNA|ANH
ANAANMNNM
aNND
Appendix C.2
SUBROUTINE COSFIT(A,B,PHI,X,Y,N,AMIN, BMIN, PHIMIN, ILIST)
THIS ROUTINE PERFORMS A BEST FIT TO DATA WITH A TRIGONOMETRIC
FUNCTION OF THE FORM Y=A+B*(COS(4*PI*(X-PHI))) USING AN ITERATIVE
METHOD AS DESCRIBED IN SCARBOROUGH(1930), ART. 115.
INPUTS ARE
A=CONSTANT
B=AMPLITUDE OF SINUSOIDAL COMPONENT
PHI= ANGLE OF MAXIMUM AMPLITUDE
X= ARRAY OF DEPENDENT VARIABLE VALUES
Y= ARRAY OF DEPENDENT VARIABLE VALUES
N= NUMBER OF DATA PAIRS OF X AND Y
AMIN= MINIMUM VALUE FOR A CORRECTION TO STOP ITERATION
BMIN= MINIMUM VALUE FOR B CORRECTION TO STOP ITERATION
PHIMIN= MINIMUM VALUE FOR PHI CORRECTION TO STOP ITERATION
ILIST= 1 FOR SUMMARY OF ITERATION PROCESS, = 0 , NO LISTING
PROGRAMMED BY C.G.FOX-ADVANCED TECHNOLOGY STAFF ,NAVOCEANO, 8/30/83
DIMENSION X(1024),Y(1024)
REAL 19,J9,K9,L9
COMPUTE INITIAL ESTIMATE OF A AND B BY PERFORMING A SIMPLE
LINEAR FIT ON LOG TRANSFORMED DATA
XN=FLOAT(N)
XPRCD=0.0
XSUM=0.0
YSUM=0.0
XSQR=0.0
WRITE(6,'(2F10.4)') (X(I),¥(1),1=1,N)
FIND MEAN LEVEL AND CONVERT X TO RADIANS
DO 50 I=1,N
X(T)=X(1)/57.2957795
50 AO=A0+(Y(1)/XN)
LOCATE LARGEST POSITIVE DIFFERENCE FROM THE MEAN AND ITS PHI
00 60 I=1,N
IF((Y(I)-AO).LT.BO) GO TO 60
BO=Y(1I)-AO
PHIO=X(1)
60 CONTINUE
COMPUTE SUM OF SQUARES OF THE RESIDUALS
POLD=0.0
DO 100 I=1,N
100 POLD=POLD+((Y(1)-F3(X(1),A0,BO,PHIO) )**2)
ITERAT=0
NBIS=0
IF (ILIST.EQ.1)WRITE (6,110) :
110 FORMAT(' ITERATION # OF BISECTIONS A B
181
aNND
AND
aan (ok i a)
anNMND
*PHIO RES IDUALS**2)
IF(ILIST.EQ.1)WRITE (6,120) ITERAT,NBIS,A0,BO,PHIO,POLD
120 FORMAT(’ ',16,10X,14,6X,4(4X,F10.4))
150 A9=0.
200
ZERO OUT MATRIX TERMS
B9=0.
c9=0.
E9=0.
F9=0.
19=0.
J9=0.
K9=0.
L9=0.
(efookokoko koe)
COMPUTE TERMS FOR LEAST SQUARES MATRIX CONSTRUCTION
DO 200 I=1,N
PARTA=F1(X(I),A0,B0)
PARTB=F2(X(1I),A0,PHIO)
PARTP=F4(X(1),A0,80,PHIO)
POWF =F 3(X(I),A0,B0,PHIO)
A9=A9+(PARTA**2 )
B9=B9+(PARTA*PARTB )
C9=C9+(PARTA*PARTP )
E9=E9+(PARTB**2)
F9=F9+(PARTB*PARTP )
19=19+(PARTP**2 )
J9=J9+ (PARTA* (Y(1)-POWF ) )
K9=K9+ (PARTB*(Y(1I)-POWF ) )
L9=L9+(PARTP*(Y(1)-POWF ) )
D9=B9
G9=C9
H9=F9
COMPUTE CORRECTION TERMS FOR A AND B
D=A9*E9*19+B9*F 9*G9+C9*D9*H9 -C9*E 9*G9-B9*D9*19-A9*F9*HO
ACORR=(( (E9*19-F9*H9 )*J9 )+ ( (C9*H9-BO*I9)*KI)+( (BO*F9-E9*C9)*LO))/D
BCORR=( ( (G9*F9-D9*I9 )*J9)+( (A9*19-G9*C9)*K9)+ ( (D9*C9-AI*F9)*LI9))/D
PCORR=(( (D9*H9-G9*E9 )*J9 )+ ( (G9*B9-A9*H9 )*K9)+( (AI*E9-B9*DI)*LO) )/D
CREATE NEW A AND B
230
A=A0+ACORR
B=B0+BCORR
PHI=PHIO+PCORR
COMPUTE NEW SUM OF SQUARES OF RESIDUALS WITH NEW ESTIMATES
PNEW=0.0
DO 250 I=1,N
182
ANNND
QanNM
AND
(op)
250 PNEW=PNEW+((Y(1)-F3(X(I),A,B,PHI) )**2)
TEST FOR CONVERGENT SOLUTION(PNEW < POLD)
IF NOT, BISECT CORRECTIONS AND RECOMPUTE
IF(PNEW.LT.POLD) GO TO 300
ACORR=.5*ACORR
BCORR=.5*BCORR
PCORR=.5*PCORR
NBIS=NBIS+1
IF(NBIS.GT.10) GO TO 300
GO TO 230
TEST FOR MINIMUM CHANGE OF A AND B
300 IF(ABS(A-AO).GT.AMIN) GO TO 500
IF (ABS(B-BO).GT.BMIN) GO TO 500
IF (ABS(PHI-PHIO).GT.PHIMIN) GO TO 500
GO TO 900
CORRECTION TERM NOT FINE ENOUGH, START NEW ITERATION
500 ITERAT=ITERAT+1
AO=A
BO=B
POLD=PNEW
IF (ILIST.EQ.1)WRITE (6,520) ITERAT ,NBIS,A0,B80,PHIO,POLD
520 FORMAT(' ',16,10X,14,6X,4(4X,F10.4),
NBIS=0
GO TO 150
900 ITERAT=ITERAT+1
IF (ILIST.EQ.1)WRITE (6,920) ITERAT,NBIS,A0,B0,PHIO,POLD
920 FORMAT(' ',16,10X,14,6X,4(4X,F10.4))
RECONVERT RADIANS TO DEGREES
DO 930 I=1,N
930 X(1)=X(1)*57.2957795
PHI=PHI*57.2957795
PRINT *,A,B,PHI
RETURN
END
FUNCTIONS TO CALCULATE POWER LAW FUNCTION AND PARTIAL
DERIVATIVES WITH RESPECT TO A AND B
FUNCTION F1(X3,A3,B3)
cae PARTIAL OF FUNCTION WITH RESPECT TO A
Fl=1.
A3=A3
X3=X3
B3=B3
RETURN
183
END
FUNCTION F2(X4,A4,PHI4)
CALCULATE PARTIAL OF FUNCTION WITH RESPECT TO B
F2=COS(2.*(X4-PHI4) )
A4=A4
RETURN
END
FUNCTION F3(X5,A5,B5,PHI5)
CALCULATE FUNCTION
F 3=A5+B5*COS (2.* (X5-PHI5) )
RETURN
END
FUNCTION F4(X6,A6,B6,PHI6)
CALCULATE PARTIAL WITH RESPECT TO PHI
F4=2.*SIN(2.*(X6-PHI6) )
A6=A6
B6=B6
RETURN
END
184
Appendix D
NOTE: Many of the subroutines required by the following
routines are common to the province selection
software listed in Appendix B.2.
GENERATE AMPLITUDE SPECTRUM OF DEPTHS OUTPUT FROM PROVINCE PICKER.
PROGRAM IS A SLIGHT MODIFICATION OF THE ORIGINAL CREATED BY
T.M.DAVIS, IN WHICH THE INPUTS HAVE BEEN SET WITHIN THE CODE
AND ONLY THOSE PARAMETERS NECESSARY FOR DAILY USE HAVE BEEN LEFT
FOR INPUT BY THE USER. INPUT IS EXPECTED FROM UNIT 13, WHICH IS
THE UNIT NUMBER USED BY PROVPICK. THE PROGRAM ADDS AN INTERACTIVE
GRAPHICS SECTION WHICH ALLOWS MODIFICATION OF THE INPUT SERIES,
DELETION OF SECTIONS, AND CONTROL OF THE INDEPENDENT VARIABLE FOR
REGRESSION ANALYSIS. ROUGHNESS MODEL PARAMETERS A & B ALONG WITH
LATITUDE AND LONGITUDE INFORMATION FOR EACH PROFILE ARE OUTPUT
TO UNIT 14 FOR LATER USE IN PROVCHART PROGRAM.
PROGRAMMMED BY C.G.FOX,ADVANCED TECHNOLOGY STAFF ,NAVOCEANO 5/17/84
THE FOLLOWING COMMENTS ARE FROM THE ORIGINAL PROGRAM BY T.M.DAVIS
PROGRAM TO COMPUTE PREWHITENED,CORRECTED AND SMOOTHED AMPLITUDE
SPECTRUM CUT,H,N=PARAMETERS TO CONTROL PREWHITENING AND SMOOTHING
FILTERS ,ANORM=SPECTRUM NORMALIZATION IN TERMS OF DATA INTERVALS
XORG,YORG=LOG-LOG PLOT ORIGIN IN POWERS OF 10, ITG=HEADING
IUNIT=5 IF CONTROL DATA ON CARDS = TAPE UNIT IF ON TAPE
SAME FOR JUNIT FOR INPUT DATA
TAPE UNIT 2=QUTPUT
CODE IFEOF=1 IF EOF FOLLOWS EACH SET
FIRST SET OF FILTER PARAMETERS ARE FOR HIGH PASS PREWHITENING
SECOND SET ARE FOR LOW PASS SMOOTHING
IF N(1)=-1 COSINE TAPER IS APPLIED,=0 NO TAPER OR PREWHITENING
LEAVE HIGH OR LOW PASS PARAMETERS BLANK IF NO FILTER DESIRED
INPUT DATA LIMIT IS 2048 PTS,NP=NO.OF INPUT PTS
IF ANORM IS BLANK SET TO 1 IN PROGRAM
USES SUBROTINE NLOGN FROM ROBINSON
IN FIRST CONTROL CARD NSETS=NO.OF SETS OF DATA. THIS RUN
SET KPHA=1 IF PHASE SPECTRUM IS DESIRED
IORGN=NO.OF INPUT DATA PT.TO USE AS ORIGIN FOR PHASE SPECTRUM
IF PHASE IS DESIRED CODE IPHA=0 IF PHASE IN DEGREES OR CODE IPHA
=NO.OF DATA INTERVALS/INCH FOR PLOT IF YOU WANT PHASE IN DATA INT.
SET JCODE=1 IF ONLY PLOT OF FINAL SMOOTHED SPECTRUM DESIRED
SET JCODE=-1 FOR NO PLOT AT ALL
CODE ILIST=1 IF YOU DESIRE ONLY LISTING OF SMOOTHED SPECTRUM
**TMPORTANT-IF YOU HAVE ALREADY ADDED ZEROES TO THE BEGINNING OR END
OF YOUR INPUT DATA AND YOU REQUEST PREWHITENING OR COSINE TAPER
YOU MUST CODE IADJ=1 FOR PROPER EXECUTION
IMEAN=1 IF YOU WANT THIS MEAN REMOVED
DIMENSION DATA(4500) ,STO(2500) ,CUT(2) ,H(2) ,N(2) ,ADATA( 4100)
185
1,SDATA( 4500) , BDATA( 2050) ,CDATA( 2050) , PHASE (2050) ,SMPH( 2050)
COMMON /BOUND/I, IPNUM, ISDLAT ,SELAT, ISDLNG,SELNG, IDLAT,EELAT,
LIDLNG,EMLNG
CHARACTER*4 ITH(8) ,ITG(8)
CHARACTER*1 RESPON
C READ(5,* )NSETS , JCODE , IUNIT, JUNIT, IFEOF ,GRIDNM
NSETS=1000
JCODE =1
IUNIT=5
JUNIT=2
IFEOF=1
C GRID SPACING IS READ FROM DATA FILE
IF(JCODE.GE.0) CALL PLOTS(STO,2500,2)
CALL FACTOR(.65)
JSET=1
IF(JCODE.GE.0) CALL PLOT(1.0,0.0,-3)
READ(IUNIT ,*)NP ,KXORG,KYORG,ANORM1, (CUT(I),H(I),N(1),1=1,2),IMEAN,
1 KPHA,IPHA, IORGN, ILIST, IADJ
NUMBER OF POINTS IS READ FROM DATA FILE, MAXIMUM IS HELD TO 2000
NPTS=2000
KXORG=-2
KYORG=-7
ANORM1=2.0
CUT(1)=-1.0
H(1)=.2
N(1)=3
CUT (2)=.08
H(2)=.2
N(2)=3
IMEAN=1
KPHA=0
IPHA=0
IORGN=1
ILIST=0
IADJ=1
245 CONTINUE
APHA=IPHA
XORG=10.0**KXORG
YORG=10. 0**KYORG
IOUT =7
IF (ABS (ANORM1) .LT.1.0E-20) ANORM=1.0
CALL PRVOUT(NP,DATA, ITG,XINTER,SLOPEX ,GRIDNM)
GR IOKM=GRIDNM*1. 852
C XINTER=10** (ALOG1O(XINTER )-(ALOG10(1. /GRIDKM)+.5))
2451 CTOFF2=.5*(1./GRIDKM)
CTOFF 1=.01*(1. /GRIDKM)
CTOFF 1=AL0G10(CTOFF 1)
CTOFF 2=AL0G10(CTOFF 2)
CT ISV=CTOFF 1
CT 2SV=CTOFF 2
C SAVE DATA ARRAY IN SDATA ARRAY
DO 2452 I=1,NP
2452 SDATA(I)=DATA(I)
(on ap ap)
186
25 IF (ANORM1.GT.1. )ANORM=1./FLOAT (NP)
38 IF(CUT(1).NE.-1.0) GO TO 29
FUNDAMENTAL FREQUENCY ,H(1)=.2,N(1)=3
CUT (1)=.8/ (FLOAT (NP ) )
ACL) =<12
N(1)=3
IF (ANORM1.GT.1.) ANORM=1./(FLOAT(NP)-(2*(N(1)-1)))
29 IF(IADJ.NE.1) GO TO 383
COMPUTE INDEX OF FIRST NON-ZERO POINT
28 00 385 J=1,NP
Kl=J
IF (DATA(J).NE.0.0) GO TO 386
385 CONTINUE
WRITE (IOUT, 82) XORG, YORG,ANORM, ITG
82 FORMAT(15H PLOT ORIGIN X=,E11.4,4H Y= ,£11.4,15HNORMALIZATION =,
1£11.4,30A1)
COMPUTE INDEX OF LAST NON-ZERO POINT
386 DO 391 J=NP,K1,-1
K2=J
IF (DATA(J).NE.0.0) GO TO 392
391 CONTINUE
SHIFT DATA TO LEFT AND RECOMPUTE NO.OF INPUT POINTS
392 NP=K2-K1+1
DO 393 J=1,NP
K=K1+J-1
393 DATA(J)=DATA(K)
DEMEAN DATA AND APPLY COSINE TAPER
383 SUM=0.0
XNP=NP
DO 388 J=1,NP
388 SUM=SUM+DATA( J)
AVE =SUM/XNP
IF(N(1).GE.0.AND.IMEAN.EQ.0) GO TO 381
DO 389 J=1,NP
389 DATA(J)=DATA(J)-AVE
CALL BATHXS(DATA,NP,GRIDKM)
IF(N(1).GE.0) GO TO 381
DO 382 J=1,NP
AM=J-1
TNP=NP -1
TNP=TNP/2.0
AM=3.14159* (AM-TNP ) /TNP
382 DATA(J)=DATA(J)*0.5*(1.0+COS (AM) )
RAISE NO.OF PTS TO A POWER OF 2(M)AND STORE IN INP
381 KOUNT=1
DO 5 M=1,12
[=2**™
IF(I-NP)5,6,7
5 CONTINUE
6 INP=I
GO TO 9
7 INP=I
NPP=NP+1
187
DO 8 J=NPP, INP
8 DATA(J)=0.0
9 AINP=INP
DELR=1.0/AINP
DELK=(1.0/AINP) /GRIDKM
C INP= TOTAL NO.OF DATA PTS.TO A POWER OF 2
c DELR=NORMALIZED FREQUENCY INTERVAL IN CYCLES PER DATA INTERVAL
WRITE (IOUT ,50)DELR, INP,NP
50 FORMAT(' FREQ. INT=',£8.3,' POWER OF 2 DATA PTS=',16,' NO.OF INPUT
1 DATA PTS=',16)
WRITE ( IOUT, 560) IORGN
560 FORMAT(' ORIGIN FOR PHASE SPECTRUM IS INPUT PT.NO.',14)
C NOW STORE DATA IN COMPLEX FORM IN ADATA
53 00 2. J=1,INP
1=(J*2)-1
ADATA(I)=DATA(J)
IA=I+1
2 ADATA(IA)=0.0
c COMPUTE AND PLOT NORMALIZED AMPLITUDE SPECTRUM
CALL NLOGN(M,ADATA, -1.0)
INX=INP+1
IF (KOUNT-2) 301,301, 302
C CORRECT FFT OF PREWHITENING FILTER FOR PHASE SHIFT
302 ANN=N(1)
CALL TSHIFT(M,ADATA, ANN, DELR)
DO 303 I=1,INX,2
JX=(1+1)/2
J=I+1
BDATA( JX )=(SQRT (ADATA(I)**2+ADATA(J)**2) ) *ANORM
IF(ADATA(I).LT.0.0) BDATA( JX )=-BDATA( JX )
303 CONTINUE
GO TO 204
301 DO 201 I=1,INX,2
JX=(1+1)/2
J=I+1
201 BDATA( JX )=( SQRT(ADATA(I)**2 +ADATA(J)**2) )*ANORM
c COMPUTE ROUGH PHASE SPECTRUM FROM PREWHITENED DATA
G OR COSINE TAPERED DATA
IF (KPHA.EQ.1.AND.N(1).LE.0) GO TO 567
IF (KPHA.EQ.1.AND.KOUNT.£Q.2) GO TO 567
GO TO 568
567 ANN=IORGN-1
C CORRECT PHASE SPECTRUM FOR DESIRED ORIGIN
CALL TSHIFT(M,ADATA,ANN, DELR)
DO 569 I=1,INX,2
JX=(I+1)/2
B=JX-1
B=B*DELR*360.0
J=I+1
IF (ABS (ADATA(I)).LT.1.£-20) ADATA(I)=ADATA(I)+1.€-20
PHASE ( JX )=(ATAN2(ADATA(J) ,ADATA(I)) )*57.295779
569 IF(IPHA.GT.0.AND.B.GT.0)PHASE ( JX )=PHASE ( JX )/B
IF(IPHA.GT.0) GO TO 541
92
93
94
95
97
83
99
206
WRITE (IOUT, 570)
FORMAT (' ROUGH PHASE SPECTRUM(DEG)' )
GO TO 542
WRITE (IOUT , 543)
FORMAT(' ROUGH PHASE SPECTRUM(DATA INTERVALS )')
IF(ILIST.NE.1) WRITE (IOUT, 81) (PHASE (J),J=1,JX)
IF (KOUNT -2) 202,203, 204
WRITE (IOUT, 75)
FORMAT (28H AMP.SPECT.OF ORIGINAL DATA )
IF(ILIST.NE.1) WRITE (IOUT,81)(BDATA(J),J=1,Jx )
FORMAT (10E11.4)
GO TO 206
WRITE (IOUT, 76)
FORMAT (30H AMP.SPECT.OF PREWHITENED DATA)
GO TO 77
WRITE (IOUT,91) CUT(1),H(1),N(1)
FORMAT(' SPECT.OF PREWHITENING FILTER WITH PHASE REVERSALS CUTOF
PSY Pe G Se Bs SP Alek ats2 =',13)
DO 92 J=1,JX
BDATA(J) =BDATA(J)/ANORM
IF(ILIST.NE.1) WRITE (IOUT, 81)(BDATA(J),J=1, JX)
CORRECT SPECTRUM OF FILTERED DATA FOR PREWHITENING FILTER
D0 93 J=1,JX
IF (ABS (BDATA(J)).LT.1.£-20)BDATA(J)=BDATA(J)+1.E-20
BDATA(J)=ABS (CDATA(J)/BDATA(J))
WRITE (IOUT , 94)
FORMAT(' FINAL ROUGH CORRECTED SPECTRUM WITHOUT PHASE REVERSALS')
IF(ILIST.NE.1) WRITE (IOUT,81)(BDATA(J),J=1, Ux )
GO TO 206
CONT INUE
SMOOTH FINAL ROUGH CORRECTED SPECTRUM
COMPUTE LOW PASS WEIGHTS
K=2
IF(N(2).EQ.0) GO TO 330
STORE ROUGH SPECTRUM IN DATA
00 97 J=1,JX
DATA(J)=BDATA(J)
GO TO 98
DO 99 I=1,NA
BDATA(1)=BDATA(I )+DELTA/X
CONVOLVE LOW PASS SMOOTHING WEIGHTS WITH ROUGH SPECTRUM
NB=JX-N(K)
GO TO 235
IF (JCODE.EQ.1.AND.N(2).GT.0) GO TO 311
IF(JCODE.LT.0) GO TO 311
X=ALOG1O ( DELK/XORG)*3.0125
TPP=BDATA(2)/YORG
IF(TPP.LT.1.E-20) TPP=YORG
Y=-0.5+AL0G10(TPP )*1. 35714
CALL PLOT (XS Y, 3)
D0 42 J=3,JX-1
XJ=J-1
X=AL0G10((XJ*DELK ) /XORG)*3.0125
189
TPP=BDATA(J)/YORG
IF(TPP.LT.1.E-20) TPP=YORG
Y=-0.5+AL0G10(TPP)*1. 35714
42 CALL PLOT(X,Y,2)
311 IF (KOUNT-2)78, 79,95
78 KOUNT =2
COMPUTE WEIGHTS FOR PREWHITENING FILTER AND STORE IN BDATA
IF(N(1).LE.0) GO TO 95
K=1
98 NA=N(K)+1
BDATA(1)=2.0* (CUT (K)+H(K) )
CENTER WEIGHT STORED IN LOCATION 1
SUM=0.0
PI=3. 1415926535898
COMPUTE REMAINING WEIGHTS
DO 40 I=2,NA
P=I-1
Q=1.0-16.*H(K)**2*P**2
IF( ABS(Q).GT.0.0001) GO TO 11
BDATA(I)= SIN(2.*PI*P* (CUT (K)+H(K)))/(4.0*P)
GO TO 40
11 BDATA(I)=(( COS(2.*PI*P*H(K)))/Q)*(( SIN(2.*PI*P*(CUT(K)+H(K))))
1 /(PI*P))
40 SUM=SUM+BDATA (I)
CORRECT WEIGHTS FOR UNITY GAIN
DELTA=1.-(BDATA(1)+2.*SUM)
X=2*N(K)+1
COMPUTE FINAL CORRECTED WEIGHTS ,HIGH PASS OR LOW PASS
IF(K.EQ.2) GO TO 83
DO 41 I=2,NA
41 BDATA(I)= (BDATA(I)+ DELTA/X )*(-1.0)
BDATA(1)=1.0-(BDATA(1)+DELTA/X )
NB=NP-N(K)
CONVOLVE HIGH PASS WEIGHTS WITH ORIGINAL DATA,STORE IN ADATA
235 DO 44 I=NA,NB
IA=I-NA+1
SUM=0.0
DO 43 J=1,NA
Jl=I+J-1
J2=1-J+1
43 SUM=SUM+BDATA(J)*(DATA(J1) + DATA(J2))
44 ADATA(IA)=SUM-BDATA(1)*DATA(I)
IF(K.£Q.2) GO TO 236
STORE PREWHITENED DATA IN DATA AND FILL IN ZEROES
-NJ=NB-N(K )
DO 51 J=1,NJ
51 DATA(J)=ADATA( J)
NJ=NJ+1
DO 54 J=NJ, INP
54 DATA(J)=0.0 .
PUT WEIGHTS IN PROPER ORDER -N.O.N AND STORE IN CDATA
JJ=N(K)+1
DO 52 J=1,JJ
190
237
491
4911
J1=J+N(K )
J2=J0= (0-1)
CDATA(J1)=BDATA(J)
CDATA(J2)=CDATA(J1)
COMPUTE AMPLITUDE SPECTRUM OF PREWITENED DATA AND PLOT
GO TO 53
CONTINUE |
STORE HIGH PASS WEIGHTS IN DATA AND AMP.SPECT.OF PREWHITENED
DATA IN CDATA
JJ=2*N(K )+1
DO 102 J=1,JJ
DATA(J)=CDATA(J)
JJ=JJ+1
DO 104 J=JJ,INP
DATA(J)=0.0
DO 103 J=1,JxX
CDATA(J)=BDATA(J)
COMPUTE AMPLITUDE SPECTRUM OF PREWHITENING FILTER
KOUNT=3
GO TO 53
CONT INUE
OUTPUT SMOOTHED SPECTRUM AND PLOT
FIRST=(N(2)+1)*DELR
FIRSTK=(N(2)+1)*DELK
WRITE (IOUT, 237) FIRSTK
FORMAT (38H FINAL SMOOTHED AND CORRECTED SPECT. ,1l6HFIRST FREQUENC
1Y=,E8. 3)
CALL PARSVL (ADATA,NP,INP,AVE,ANOIS ,N(2) ,ANORM )
WRITE (IOUT, 491) AVE ,ANOIS
FORMAT(' VALUE OF AMP SPECT OF NOISE=',E10.4,' RMS WHITE NOISE LE
1VEL=',£10.4)
WRITE (IOUT, 81)(ADATA(J), J=2,1A)
CALL POWWGT(B0,B81,FIRSTK,DELK,ADATA(2) ,IA-1,0,CTOFF1,CTOFF2)
BOSAVE =B0
BISAVE=B1
CALCULATE RMS FOR WHITE NOISE WITH INTERCEPT (BO)
AVE =B0
AVE =AVE /ANORM
AVE =AVE**2
ANOISE=(AVE*INP)/(NP*INP )
ANOISE =SQRT (ANOISE )
WRITE (IOUT, 492)AVE,ANOISE
FORMAT(' AVE= ',F8.3,' RMS WHITE NOISE = ',F8.3)
IF(JCODE.LT.0) GO TO 486
ENCODE(13,601,ITH) Bl
FORMAT('SLOPE =',F6.3)
CALL SYMBOL (7.5,8.5,.20,ITH,0.,13)
ENCODE (19,602, 1TH)BO
FORMAT('INTERCEPT =',£8.3)
CALL SYMBOL (7.5,8.0,.20,ITH,0.,19)
PLOT REGRESSION LINE FROM POWWGT
CALL NEWPEN( 3)
IPLPSS=0
191
603 FIRST1=FIRSTK
C WRITE (IOUT, '(8HFIRST1= ,£11.4)' )FIRST1
X=ALOG10(FIRST1/XORG)*3.0125
AMP 1=BO*F IRST1**B1
C WRITE (IOUT, '(6HAMP1= ,£11.4)')AMP1
TPP=AMP1/YORG
IF(TPP.LT.1.E-20) TPP=ABS(TPP)
Y=-0.5+AL0G10(TPP )*1. 35714
CALL PLOT(X,Y,+3)
C CALL SYMBOL { X,Y,.1,32,0.0,-1)
FIRST 1=((1A)*DELK )
C WRITE (IOUT, '(8HFIRST1= ,£11.4)')FIRST1
X=ALOG10(FIRST1/XORG)*3.0125
AMP 1=BO*F IRST1**B1
C WRITE (IOUT, '(6HAMP1= ,£11.4)' )AMP1
TPP=AMP 1/YORG
IF(TPP.LT.1.E-20) TPP=ABS (TPP)
Y=-0.5+AL0G10(TPP )*1. 35714
CALL PLOT(X,Y,+2)
CALL SYMBOL ( X,Y,.1,32,0.0,-2)
CALL NEWPEN(1)
PLOT SPECTRAL ESTIMATES
X=ALOG10 (FIRSTK/XORG)*3.0125
TPP=ADATA( 2) /YORG
IF (TPP.EQ.0.0)TPP=1.E-15
IF(TPP.LT.1.E-20) TPP=ABS (TPP)
Y=-0.5+AL0G10(TPP)*1. 35714
CALL PLOT(X,Y,+3)
C ALL SYMBOL ( X,Y,.05,32,0.0,-1)
DO 238 J=3,IA-1
XJ=J-1+N(2)
X=ALOG10( (XJ*DELK ) /XORG)*3.0125
TPP=ADATA(J)/YORG
IF (TPP.EQ.0.0)TPP=1.E-15
IF(TPP.LT.1.E-20) TPP=ABS (TPP)
Y=-0.5+ALOG10(TPP)*1. 35714
238 CALL PLOT(X,Y,+2)
C 238 CALL SYMBOL ( X,Y,.05,32,0.0,-2)
CALL PLOT(0.0,0.0,+3)
IF (IPLPSS.EQ.1) GO TO 330
IF(LPASS.EQ.1) GO TO 330
IPLPSS=1
BO=XINTER
B1=SLOPEX
INCLUDE NEXT LINE IF A PLOT OF SPATIAL DOMAIN ESTIMATE IS DESIRED
GO TO 603
PLOT AXIS
330 IF(ANORM.EQ.1.0) GO TO 333
CALL AXES(0.,-.5, ‘NORMALIZED AMPLITUDE (IN KM)',27,9.5,90.0
wes 4) O05 Os 0s— 1)
GO TO 332
333 CALL AXES(0.,-.5, AMPLITUDE (IN KILOMETERS)',25,9.5,90.0
*,1.3514,10.0,0.0,-1)
AND
(ok ek a)
192
C
C
C
332 CALL AXES(0.,-0.5, 'FREQUENCY (CYCLES /KILOMETER)',27,12.2,0.0
POOL 25 LOO NOnOee— 9)
CALL SYMBOL (4.0,-1.4,.32,1TG,0.0, 30)
AYORG=KYORG
AXORG=KXORG
00 362 J=1,8
AA=J-1
Y=(1.35714*AA)-0. 37
FPN=AYORG+AA
362 CALL NUMBER(-.25,Y,.16,FPN,90.0,-1)
00 363 J=1,5
AA=J-1
X=(3.0125*AA )+0. 20
FPN=AXORG+AA
363 CALL NUMBER(X,-.25,.16,FPN,0.0,-1)
DO 364 J=1,4
IB=KXORG+J-1
00 364 JC=2,9
AC=JC
AA=AC*10.0**IB
X=AL0G10(AA/XORG)*3.0125
CALL PLOT(X,-.5,3)
364 CALL PLOT(X, -.46,2)
00 365 J=1,7
IB=KYORG+J-1
DO 365 JC=2,9
AC=JC
AA=AC*10.0**IB
Y=-0.5+AL0G10(AA/YORG)*1. 35714
CALL PLOT(0.0,Y,3)
365 CALL PLOT(-.04,Y,2)
CALL PLOT(0.0,-1.0,3)
CALL PLOT(0.0,-1.0,3)
+
INTERACTIVE ROUTINE TO MODIFY INPUT STRING
WRITE (6, 3641)
3641 FORMAT(T70,' ARE YOU SATISFIED WITH THE DATA SET BOUNDARIES?'/
*T70,' ENTER E TO END RUN')
READ(5,3652) RESPON
IF(RESPON.EQ.'E') GO TO 487
IF (RESPON.EQ.'Y') GO TO 3650
WRITE (6, 3642)
3642 FORMAT(T70,' WOULD YOU LIKE TO ELIMINATE THE ENTIRE SEGMENT?')
READ(5, 3652) RESPON
IF (RESPON.EQ.'N') GO TO 3644
CALL ERASE
GO TO 3658
3644 WRITE (6, 3643)
3643 FORMAT(T70,' ENTER FIRST,SECOND LOCATIONS(IN KM) FOR TRUNCATION' )
READ(5,*)FIR,SEC
NF IR=FIR/GRIDKM
IF(NFIR.LT.O) NFIR=0
193
NSEC=SEC/GRIDKM
IF (NSEC.GT.NP) NSEC=NP
C RECREATE DATA ARRAY
DO 3645 ICHG=1,NP
LOCCHG=ICHG+NFIR
DATA( ICHG )=SDATA(LOCCHG )
IF (LOCCHG.EQ.NSEC) GO TO 3646
3645 CONTINUE
C INTERPOLATE NEW LAT & LON
3646 SDIFF=FLOAT (NSEC )/FLOAT (NP)
FDIFF =FLOAT (NFIR)/FLOAT (NP )
FLAT =DECDEG( ISDLAT ,SELAT )
FLON=DECDEG( ISDLNG,SELNG)
SLAT =DECDEG(IDLAT ,EELAT )
SLON=DECDEG(IDLNG,EMLNG)
FNLAT=((SLAT-FLAT )*FDIFF )+FLAT
SNLAT=((SLAT -FLAT )*SDIFF )+FLAT
FNLON=((SLON-FLON )*FDIFF )+FLON
SNLON=((SLON-FLON )*SDIFF )+FLON
C RECONSTRUCT INTO DEGREES AND MINUTES
CALL DEGMIN(FNLAT, ISDLAT ,SELAT )
CALL DEGMIN(FNLON, ISDLNG,SELNG)
CALL DEGMIN(SNLAT, IDLAT ,EELAT )
CALL DEGMIN(SNLON, IOLNG,EMLNG)
NP=NSEC-NFIR
CALL ERASE
GO TO 2451
INTERACTIVE ROUTINE fu MODIFY REGRESSION FIT TO SPECTRUM
AND
3650 WRITE (6, 3651)
3651 FORMAT(T70,' ARE YOU SATISFIED WITH THE REGRESSION FIT?')
READ(5, 3652 )RESPON
3652 FORMAT (A1)
IF(RESPON.EQ.'Y') GO TO 3656
WRITE (6, 3653)
3653 FORMAT(T70,' ENTER LOG(LOWER FREQUENCY) ,LOG(UPPER FREQUENCY )' )
READ(5,*) CTOFF1,CTOFF2
CALL ERASE
LPASS=1
GO TO 4911
3656 CTOFF 1=CT1SV
CTOFF 2=CT 2SV
CALL ERASE
LPASS=0
WRITE (14, 399)NP, IPNUM, ISDLAT ,SELAT , ISOLNG,SELNG, IDLAT,EELAT,
LIDLNG,EMLNG, BISAVE , BOSAVE
399 FORMAT(15,18,4(1X,14,F6.2),F6.3,£9.3)
C PLOT PHASE SPECTRUM
3658 IF(KPHA.NE.1) GO TO 581
CALL PLOT(15.0,0.0,-3)
X=ALOG10 (DELK /XORG)*3.0125
Y=ABS (PHASE (2) /30.0)
194
IF(IPHA.GT.0) Y=Y*30.0/APHA
IBB=4
IF(Y.LT.8.0.AND.PHASE(2).GE.0.0) GO TO 552
IF(Y.LT.8.0.AND.PHASE(2).LT.0.0) IBB=2
IF(Y.GE.8.0) Y=8.5
CALL SYMBOL (X,Y,.05,188,0.,-1)
GO TO 553
552 CALL PLOT(X,Y, 3)
553 DO 571 J=3,JxX
XJ=J-1
X=AL0G10((XJ*DELK ) /XORG)*3.0125
Y=ABS (PHASE (J )/30.0)
IF(IPHA.GT.0) Y=¥*30.0/APHA
IBB=4
IF(Y.LT.8.0.AND.PHASE(J).GE.0.0) GO
IF(Y.LT.8.0.AND.PHASE(J).LT.0.0) IBB
IF(Y.GE.8.0) Y=8.5
CALL SYMBOL (X,Y,.05,1BB,0.,-2)
GO TO 571
554 CALL PLOT(X,Y,2)
571 CONTINUE
486 CONTINUE
COMPUTE AND PLOT SMOOTHED PHASE SPECTRUM
IF(N(2).EQ.0) GO TO 579
DO 572 I=NA,NB
IA=I -NA+1
SUM=0.0
DO 573 J=1,NA
Jl=I+J-1
J2=I-J+1
573 SUM=SUM+BDATA( J )* (PHASE (J1)+PHASE (J2) )
572 SMPH(IA)=SUM-BDATA(1)*PHASE (1)
IF(IPHA.GT.0) GO TO 544
WRITE (IOUT ,582) FIRSTK
‘582 FORMAT(' SMOOTHED PHASE SPECTRUM(DEG)-FIRST FREQ.=',E8.3)
GO TO 545
544 WRITE(IOUT,546) FIRSTK
546 FORMAT(' SMOOTHED PHASE SPECTRUM(DATA INT. )-FIRST FREQ.=',F8.5)
545 WRITE( IOUT, 81) (SMPH(J),J=2,1A)
IF(JCODE.LT.0) GO TO 581
CALL NEWPEN(1)
X=ALOG10(FIRSTK /XORG)*3.0125
Y=ABS (SMPH(2)/30.0)
IF(IPHA.GT.0) Y=Y¥*30.0/APHA
IBB=4
IF(Y.LT.8.0.AND.SMPH(2).GE.0.0) GO TO 548
IF(Y.LT.8.0.AND.SMPH(2).LT.0.0) IBB=2
IF(Y.GE.8.0) Y=8.5
CALL SYMBOL (X,Y,.05,1BB,0.,-1)
GO TO 549
548 CALL PLOT(X,Y,3)
549 DO 574 J=3,IA
XJ=J-1+N(2)
TO 554
=2
195
X=AL0G10((XJ*DELK ) /XORG)*3. 0125
Y=ABS (SMPH(J)/30.0)
IF (IPHA.GT.0)Y=Y*30.0/APHA
IBB=4
IF (Y.LT.8.0.AND.SMPH(J).GE.0.0
IF(Y.LT.8.0.AND.SMPH(J).LT.0.0
IF(Y.GE.8.0) Y=8.5—
CALL SYMBOL (X,Y,.05,IBB,0. ,-2)
GO TO 574
550 CALL PLOT(X,Y,2)
574 CONTINUE
PLOT PHASE AXIS
579 IF(JCODE.LT.0) GO TO 581
CALL NEWPEN(1)
IF (IPHA.GT.0) GO TO 547
CALL AXES(0.,0.,14HABS PHASE (DEG),14,6.,90.,1.,0.,30.,-1)
GO TO 558
547 CALL AXES(0.,0.,20HABS PHASE(DATA INT.),20,8.,90.,1.,0.,APHA, -1)
558 CALL AXES(0.,0.,1TG,+30,12.2,0.,3.0125,10.,0.,-1)
CALL SYMBOL (2.0,-.2,.14, 33HTRIANGLE INDICATES NEGATIVE PHASE,0.0,
1 33)
581 CONTINUE
IF(N(2).EQ.0) GO TO 331
COMPUTE AMPLITUDE SPECTRUM OF SMOOTHING FILTER
wo
GO TO 550
IBB=2
DELX= 0.01
D0 239 J=1,51
AJ=J-1
XJ=AJd*DELX
SUM=0.0
DO 240 I=2,NA
Al=I-1
240 SUM=SUM+2*BDATA(I)* COS(2.*PI*AI*XJ)
DATA(J)= SUM+BDATA(1)
239 CDATA(J)=XJ
WRITE(IOUT,241) CUT(2),H(2),N(2)
241 FORMAT(' AMP.SPECT.OF SMOOTHING FILTER CUTOFF=',F5.4,
We Ta Bein Gs Ue i avs 3053)
IF (ILIST.NE.1) WRITE (IOUT, 242) (CDATA(I) ,DATA(I),1=1,51)
242 FORMAT (8(F5.2,F8.3))
331 WRITE (IOUT, 243) JSET
243 FORMAT(' END OF DATA SET NO.',I3)
END FILE 12
END FILE 12
IF (JSET.EQ.NSETS) GO TO 244
JSET=JSET +1
IF(JCODE.LT.0) GO TO 245
CALL PLOT(15.0,0.0,-3)
GO TO 245
244 IF(JCODE.LT.0) GO TO 487
CALL PLOT(0.0,0.0,999)
487 END FILE 14
END FILE 14
STOP
END
196
(ol i a)
(ot ww)
THIS SUBROUTINE IS A SUPPLEMENT TO FFT1D PROGRAM
AND DRAWS A PROFILE OF BATHYMETRIC DATA
COMPILER (DIAG=3)
SUBROUTINE BATHXS (DATA,NPTS ,GRIDKM )
DIMENSION DATA(1)
CONVERT DEPTHS IN KILOMETERS TO METERS
WRITE (IOUT, 3)(DATA(I),1=1,10)
3 FORMAT(' ',10(F10.5,1X))
DO 5 I=1,NPTS
5 DATA(I)=DATA(I)*1000.
CALL PLOT(0.,10.,-3)
FIND MAXIMUM AND MINIMUM DEPTH
DE PMIN=0.0
DE PMAX=0.0
DO 10 I=1,NPTS
IF(DATA(I).LT.DEPMIN) DEPMIN=DATA(I)
10 IF(DATA(I).GT.DEPMAX ) DEPMAX=DATA(TI)
ROUND TO NEAREST 100
DE PMAX=(AINT (DEPMAX/100. )*100. )+100.
DEPMIN=(AINT (DEPMIN/100. )*100. )-100.
IF (GRIDKM.GT.0.0001) GO TO 11
DE PMAX=-1.0
DEPMIN=1.0
GO TO 12
11 IF(GRIDKM.GE.0.01) GO TO 12
DE PMAX=-20.0
DEPMIN=20.0
12 YSCALE=-(DEPMAX-DEPMIN)/4.
SETUP AXES
ae AXES(0.,0.,' DEPTH (IN METERS)',18,4.0,90.,1.00,DEPMAX, YSCALE
ae)
XDIST=FLOAT (NPTS )*GRIDKM
XSCALE =12.805/XDIST
XINC=1.0
55 XTICK=XSCALE
56 IF(XTICK.GT.1.) GO TO 60
XINC=2.*X INC
XTICK=2.*XTICK
GO T0 56
60 IF(GRIDKM.GT.0.0001) GO TO 70
CALL AXES(0.0,0.0,'DISTANCE (IN 0.1 METERS)',24,12.805,0.0,XTICK
HOW oWai)
G0 TO 90
70 IF(GRIDKM.GT.0.01) GO TO 80
CALL AXES(0.0,0.0, ‘DISTANCE (IN 100 METERS)',24,12.805,0.0,XTICK
3560 3a 3)
G0 TO 90
80 CALL AXES(0.0,0.0,'DISTANCE (IN KILOMETERS )',24,12.805,0.0,XTICK
*,0.0,XINC,-1)
90 DO 100 I=1,NPTS
YVAL =- ((DEPMAX-DATA(I))/YSCALE )
XVAL =(FLOAT (1 )*GRIDKM )*XSCALE
IF(I.EQ.1) CALL PLOT(0.0,YVAL, 3)
197
C
100 CALL PLOT(XVAL,YVAL, 2)
CALL PLOT(0.,-10.,-3)
CONVERT BACK TO KILOMETERS
00 105 I=1,NPTS
105 DATA(1I)=DATA(I)/1000.
RETURN
END
198
SUBROUTINE PRVOUT(NP,DEPTH,ITG,XINTER,SLOPE ,GRID)
DIMENSION DEPTH(5000)
COMMON /BOUND/1T, IPNUM, ISDLAT ,SELAT, ISDLNG,SELNG, IDLAT ,EELAT,
1IDLNG,EMLNG
CHARACTER*4 ITG(8)
READER ROUTINE FOR DEPTH OUTPUT FROM PROVINCE PICKER
READ(13,5,END=500) (ITG(I),1=1,6),SLOPE,XINTER,GRID,NP
5 FORMAT (6A4,F7.3,1X,£9.3,1X,F7.4,15)
NP=NP+1
ITG(7)=' :
ITG(8)=' :
WRITE(6,55,END=500) (ITG(I),1=1,6) SLOPE ,XINTER,GRID,NP
55 FORMAT (1T70,6A4,F7.3,1X,E10.5,F7.4,15)
READ(13,50)1, IPNUM, ISDLAT ,SELAT, ISDLNG,SELNG, IDLAT ,EELAT,
LIDLNG,EMLNG,RMS
50 FORMAT (15,18,4(1X,14,F6.2),F 10.4)
READ(13,6) (DEPTH(I),I=1,NP)
6 FORMAT (10F 8.6)
IF(NP.GT.2048) NP=2048
READ(13,7) CHECK
7 FORMAT (F8.2)
IF (CHECK .NE.9.999999) PRINT 8
8 FORMAT(' NINES RECORD DOES NOT CHECK‘)
GO TO 900
500 DEPTH(1)=9.999999
900 RETURN
END
199
C ROUTINE TO COMPUTE THE AVE VALUE OF THE LAST 10 PERCENT OF THE
C HARMONICS(AVE) OF AN AMPLITUDE SPECTRUM(DATA) AND USE THIS VALUE WITH
C PARSEVALS FORMULA TO COMPUTE THE RMS WHITE NOISE LEVEL(ANOIS).
NP=THE NO.OF ORIGINAL PTS INPUT TO THE FFT USED TO COMPUTE DATA ARRAY
C**
C**
C**
C*x*
Cx*
SUBROUTINE PARSVL(DATA,NP,INP,AVE ,ANOIS ,NWGT ,ANORM )
INP=THE NO.OF POINTS TO A POWER OF 2 USED FOR THE FFT
DATA= UNNORMALIZED FFT AMPLITUDE SPECTRUM WITH UNITS OF
10
11
INPUT UNITS/CYCLE/DATA INTERVAL AND IS OF LENGTH(INP/2)+1-2*NWGT
NWGT=LENGTH OF SMOOTHING FILTER YOU USED ON YOUR FFT -NWGT..0..+NWGT
DIMENSION DATA(1)
IF (ANORM.EQ.1.) GO TO 11
DO 10 I=1,NP
DATA(1 )=DATA(I) /ANORM
JX=INP/2+1-2*NWGT
JT =JX- (JX /10)
X=0. 2* INP
JX=X-2.0
T=0.13*INP
JT=T-2.0
DIF =JX-JT+1
ANP=NP
AVE=0.0
AVEPS=0.0
00 1 I=JT,JX
AVE=AVE+DATA(TI )
AVEPS=AVEPS+DATA(I )**2
AVE =AVE /DIF
AVEPS=AVEPS /DIF
C PARSEVALS FORMULA IS MEAN SQ NOISE=(INP*AVEPS )/(NP*INP)
ANOTS=SQRT (AVEPS /ANP )
IF (ANORM.EQ.1.) RETURN
00 12 [=1,NP
12 DATA(I )=DATA(1)*ANORM
RETURN
END
200
SUBROUTINE TSHIFT(M,DATA,A, DELF )
C*****THIS SUBROUTINE CONTAINS 26 STATEMENTS.
C*****TSHIFT APPLIES EFFECTS OF TIME SHIFT ON REAL AND IMAGINARY PARTS OF THE
C***x**FOURIER TRANSFORM.
Cx****M -= POWER OF 2 WHICH IS EQUAL TO THE NUMBER OF REAL OR IMAGINARY PARTS
(Chaat OF THE FOURIER TRANSFORM.
C*****DATA -= SERIES ... REAL PARTS ARE ODD INDEX, IMAGINARY PARTS ARE EVEN.
C*****xA -= YAXIS SHIFT, IE. TO CORRECT THE FOURIER TRANSFORM OF A FUNCTION
WHICH
Cae HAS BEEN SHIFTED A DATA INTERVALS IN THE +X DIRECTION,SIGN OF A IS
+,
C****x*DELF -= NORMALIZED FREQUENCY INCREMENT.
DIMENSION DATA(1)
M2=2%2%2M
PT =3. 1415926536
SFT =2. O*PI*A*DELF
00 1 I=1,M2,2
J=I+1
TRE=DATA(I)
TIM=DATA(J)
K=1/2
ARG=SFT *K
CN=COS (ARG)
SN=SIN (ARG)
DATA(T )=TRE*CN-TIM*SN
1 DATA(J)=TIM*CN+TRE*SN
RETURN
END
201
Uk Ae ote, wer
Appendix E
In Chapter 6, a simple model of an anisotropic surface was created,
and the variation of its spectral characteristics as a function of azi-
muth was derived analytically and confirmed by empirical results. The
surface was constructed by generating a random signal of known spectral
properties and extending each value of the series into the second dimen-
sion. In the resulting relationship,
Ne a(8)°sP(9)
The proportionality factor a(@) was shown to vary as a sinusoid and b(9)
remained constant over all 0. Results from multibeam-sonar-derived
bathymetry indicated that this simple model is adequate to describe some
actual surfaces on ‘the earth, such as the Gorda Rise spreading center.
Beyond this very simple model, a more elaborate surface can be gen-
erated by summing several signals of differing characteristics at a
variety of orientations. If one assumes that the simple corrugated
surface presented in the elementary model of Chapter 6 is the result of
a dominant, unidirectional process, then these composite surfaces would
represent areas of the sea floor where several geological processes have
been active, perhaps acting in different directions. The Gorda Rise,
particularly at the ridge crest, is dominated by the processes associ-
ated with crustal formation, and its relief does appear to conform to
the one-process surface model.
The relief of the Mendocino Fracture Zone shown in Mgures 4-7-10,
although only represented by spectra from two azimuths, indicates a
202
clear difference in spectral slope (b(®)), which would not be expected
in a simple surface dominated by one process. In fact, the main east-
west trend of the fracture zone is produced by tectonic processes which
are reflected in profiles collected in a north-south orientation. The
various mass-wasting processes acting down this slope produce a differ-
ent style of relief which is evident in profiles collected east-west.
The bathymetry shown in Migure 6-10 from the continental slope might
also be the results of two processes at work. This is reflected in the
possible variation of b(®) as shown in Mgure 6-11.
No attempt will be made to derive a general mathematical (geometri-
cal) form for the azimuthally-dependent spectra of these multiple com-
ponent surfaces. Such a treatment would be beyond the scope of this
eoide. However, an insight into the nature of these surfaces can be
gained by examining a few examples. All of the examples following are
generated using the same algorithms. First, two random signals of spec-
ified spectral characteristics (a and b) are generated using an inverse
FFT method. One signal becomes the initial row, and the other the ini-
tial column, of a 128 x 128 element matrix. Each element of the matrix
is then computed by summing the appropriate row and column element.
This method results in a surface which when sampled at O° azimuth, pro-
duces one of the input signals displaced by a random constant. For the
90° azimuth, the other input signal is sampled, similarly displaced. To
study the spectra of the combined signals in other azimuths, the matri-
ces were sampled and analyzed in an identical fashion to the bathymetry
grids of Chapter 6.
Four examples of two spectral component surfaces are shown in Fig-
ures E-1 through E-4. A variety of combinations were used. Figure E-l
203
-1000
Depth
Figure
S
INTERCEPT (¢ )
6
©
°o
N
o
°
°
5
o
E-1
Nan
ASSN
ARAM ANKI
TON GN Wi \
Y) NO PAWN WM) "
BORO NYA
EO ONO NL
i's ou
NATE
IAN “4 Ny A i
4 a a Die
ie a in .
a R X SI) ‘
SN ony
WI ti, t6
VOR
NAAN
(\ BAYS
big ‘ CNN
ate
Ge
i ey Z
lie
Ly
iy
My
M)
Ay TD Uf
TAM ITMTNS
i
"i,
% y
teh Gai
My igs
SN
SLOPE OF SPECTRUM (b )
80 100 120 140 160 180
AZIMUTH (degrees)
20 40. 60
Artificially generated surface composed of two identical
orthogonal trends, both with spectral slope b= -1.5, and
spectral intercept a= 1.0. Viewpoint is, from, the
southwest. Below, the variation of parameters a and b are
shown versus azimuth.
204
illustrates the characteristics of a surface in which two signals of
identical spectral slope (b = -1.5) and intercept (a = 1.0) were com-
bined in orthogonal directions. Arbitrarily assuming a viewpoint from
the southwest for all examples, the input signals represent pure signals -
in the east-west and north-south orientations. The resulting surface
shows a clear northwest-southeast trend. Another realization might
yield a northeast-southwest trend. With only a knowledge of the result-
ing surface, an investigator might infer a single process acting at 45°
or 135° azimuth. The trend actually results from the vector sum of two
orthogonal processes.
The distribution of spectral parameters with azimuth clearly
reflects the departure of this surface from a one process model. The
spectral slope parameter (b) does show the designated value of -1.5 at
azimuths of 0°, 90° and 180°, as it must. The corresponding tatwecene™
(a) parameters aiso correspond to the input value of 1.0, indicating
that these profiles represent uncontaminated samples of the input sig-
nals. At intermediate azimuths, however, both parameters are consis-
tently higher than those of the input signal. Whereas in the simple
model of Chapter 6 the slope (b) parameter remains constant with azi-
muth, the same parameter shows two clear maxima in the range 0° to 180°.
The intercept (a) Parameter also shows two maxima, rather than the
single maximum of the corrugated surface model. Figure E-2 presents a
similar surface in which the intercept (a) parameter in the east-west
direction has been increased to 1.2. The variation of spectral param-
eters with azimuth is also very similar to that shown in Figure E-l,
however the variation of both a and b is amplified.
205
-1000
Depth
o
°o
INTERCEPT (c )
5 8
°
to)
°o
Figure E-2
NS
MY ER Wns wi
RO)
4 x Hi KR
Va Pe {} AYA
MN RY) s, MONT ae
TLS 3 " ROOK XX)
RN NS ‘ uy AG a.
Nie
hg
OOS
TAIT
Y, hy Hy SNAG
oS RRC
Mi,
X *, We fing
SH A WN
NY
Mh ae a Y
WR
oe
LAAN
GAN ‘t iN NWN RY) Nis
NS NNN i
MY, hi nai
oy h
ee
oO
SLOPE OF SPECTRUM (b)
20 40 60 80 100 120 140 160 180
AZIMUTH (degrees)
Artificially generated surface composed of two orthogonal
trends, both with spectral slopes of b= -1.5, and a=1.0 in
the N-S direction, a=l1.2 in the E-W direction. Viewpoint
is from the southwest. The variation of 4 and 6 parame-
ters with azimuth are plotted below. :
206
In Figure E-3, the intercept terms (a) are again set equal, but the
apectrall slopes (b) of the component signals differ. The north-south
component has spectral slope b = -1.5, while the east-west component has
a spectral slope b = -1.0. The artificial surface was designed to mimic
the bathymetry of the Mendocino Fracture Zone. The north-south signal
dominates the surface, as reflected in an east-west trending contour
chart. The wavelength associated with the intercept term a, was arbi-
trarily selected to correspond to a wavelength of two data points.
Because of the higher angle of negative slope in the spectrum of the
north-south profile, this trend contains higher amplitudes in all fre-
quencies lower than one-cycle-per-two-data intervals. Only at the very
southern limit of the surface, where the north-south signal is rela-
tively constant, can the orthogonal trend be detected.
At higher frequencies, the spectra of the two baupoueat signals
cross over, and the east-west signal contains higher amplitudes and dou-
inates the surface topography. This somewhat complicated set of circum-
stances is not outside geological experience. For a terrestrial exam-
ple, envision a long ridge of several kilometers width and perhaps one
thousand meters height, oriented east-west. On the side of this ridge
are a series of north-south trending streams and gullies with relief of
tens-of-meters which shed the runoff from the ridge. If one were to
travel in the north-south direction, the effect of the streams would be
minimal, and the long wavelength shape of the ridge would present the
only obstacle. If one were to travel in the east-west direction on the
face of the ridge, the obstacles in the terrain would be dominated by
the lower amplitude, but higher frequency streambeds. This difficulty
of travel represents in a very direct sense the concept of surface
207
-1000
SR
RN AN
uannnned
RN
NAN
ETI
TES
0 TY ip
1
a
Ha
. (ri a ant) ven ROO nt
a \ ANNUAL \\
Depth
SLOPE OF SPECTRUM (b )
INTERCEPT (2)
0 20 40 60 80 Nati Waa ga LE
AZIMUTH (degrees)
Figure E-3 Artificially generated surface composed of two orthogonal
trends with identical intercepts of a= 1.0, but spectral
slopes of b= -1.5 in the N-S direction and b= -1.0 in the
E-W direction. Azimuthally dependent spectral parameters
are plotted below.
208
roughness. The artificial surface shown in Figure E-3, as well as the
Mendocino Fracture Zone shown in Chapter 4, are analogous to this hypo-
thetical example.
The computed spectral parameters shown below in Figure E-3 also
reflect the complexities of this surface. The a(0) parameters appear to
vary regularly with azimuth, although not following the cosine curve
derived by simple regression, the variation occurs in spite of the fact
that the input signals have identical a's. Only at exactly 90° azimuth,
does the a parameter jump to the input value. The b(8) values also fol-
low a complex pattern which is not described by the illustrated cosine
curve. Extensive analytical geometry would be necessary to reach an
understanding of these variations.
A final artificial surface is presented as Migure E-4. In this
case, both the spectral slopes (b) and intercepts (a) of the input sig-
nals are d‘fferent The construction is identical to that shown in Fig-
ure E-3, with the exception that the intercept (a) parameter in the
east-west direction has been increased to 2.0. Between the input values
at 0°, 90°, and 180°, the variation of the spectral parameters appears
even more complicated than the results from Figure E-3. The combination
of signals at oblique azimuths, or the combination of more than two sig-
nals, would result in an even more complex pattern.
With the insight gained by examining these artificially generated
surfaces, a more complete analysis of the anisotropy of the Mendocino
Fracture Zone can now be conducted. Figure E-5 presents a contour chart
of the surface used in this analysis, which represents a subarea of the
base chart shown in Figure 4-7. The digital data, collected by the SASS
multibeam sonar system, is gridded at a spacing of 0.05 minutes of lati-
209
tude and longitude. Figure E-6 shows a graphic representation of this
surface and its spectral parameters as a function of azimuth, in the
same format as the artificially generated surfaces shown in Figures E-1
through E-4.
Compare the bathymetric surface shown in Figure E-6 to the artific-
tially generated surface shown in Figures E-3 and E-4. The overall mor-
phologies are quite similar, with a longer wavelength, higher amplitude
component in the north-south direction relative to the orthogonal trend.
This similarity in morphologies is also expressed as a similarity in the
azimuthally dependent roughness models, although the parameters gener-
ated from the bathymetric surface are somewhat noisier than the artifi-
cially generated examples. In all cases, the slope parameter b is not
constant with azimuth, as was the case for the surfaces studied in Chap-
ter 6. Like the artificial surfaces of Figures E-3 and E-4, the spec-
tral slope is approximately b = -1.5 in the 0° or 180° azimuth and
approaches b = -1.0 for the 90° azimuth. The intercept parameter in
Figure E-6 reaches its maximum at ® =90°, similar to the example in
Figure E-4. In the case of the Mendocino Fracture Zone, this parameter
is approximately doubled in the east-west direction over the north-south
direction. This indicates that for wavelengths near 1 km, the Fracture
Zone surface is twice as rough for ® = 90° as for 9 = 0°. In longer
wavelengths, this relationship is reversed as evidenced by the much
higher total relief in the north-south direction. Such a reversal is
only possible for near orthogonal trends with different spectral slopes.
Sinusoidal regression lines, which comprise the basic model of Chapter
6, are included in Figure E-6 to emphasize how poorly this simple model
describes the two-trend case.
210
-1000-
po)
0
3 aI
{
A
il
1
t}
jt
1
1000 -|
}
\]
l)
1
"
i
"
it
i]
]
00 KF
oa
=
>
[a4
=
O
-1 os
(7)
wi
5
2
a
8)
So
iz 1.0
[vey
)
& 0.0
S
Zz
) 20 40 60 80 100 120 140 160 180
AZIMUTH (degrees)
Figure E-4 Artificially generated surface composed of orthogonal
trends with spectral parameters b= -1.5, a= 1.0 in the N-S
direction, and b= -1.0, a= 2.0 in the E-W direction. Azi-
muthally dependent spectral parameters are plotted below.
211
126°35'W
126°25'W
40°
126°35'W 126° 25'W
Figure E-5 Contour chart of a segment of the Mendocino Fracture Zone.
Automatically generated contours are based on SASS multibeam
sonar data gridded at 0.05 minutes of latitude and
longitude.
212
tude and longitude. Figure E-6 shows a graphic representation of this
surface and its spectral parameters as a function of azimuth, in the
same format as the artificially generated surfaces shown in Figures E-1
through E-4.
Compare the bathymetric surface shown in Figure E-6 to the artific-
ially generated surface shown in Figures E-3 and E-4. The overall mor-
phologies are quite similar, with a longer wavelength, higher amplitude
component in the north-south direction relative to the orthogonal trend.
This similarity in morphologies is also expressed as a similarity in the
azimuthally dependent roughness models, although the parameters gener-
ated from the bathymetric surface are somewhat noisier than the artific-
jally generated examples. In all cases, the slope parameter b is not
constant with azimuth, as was the case for the surfaces studied in Chap-
ter 6. Like the artificial. surfaces of Figures E-3 and E-4, the spec-
tral slope is approximately b = -1.5 in the 0° or 180° azimuth and
approaches b = -1.0 for the 90° azimuth. The intercept parameter in
Figure E-6 reaches its maximum at @ =90°, similar to the example in
Figure E-4. In the case of the Mendocino Fracture Zone, this parameter
is approximately doubled in the east-west direction over the north-south
direction. This indicates that for wavelengths near 1 km, the Fracture
Zone surface is twice as rough for 6 = 90° as for 8 = 0°. In longer
wavelengths, this relationship is reversed as evidenced by the much
higher total relief in the north-south direction. Such a reversal is
only possible for near orthogonal trends with different spectral slopes.
Sinusoidal regression lines, which comprise the basic model of Chapter
6, are included in Figure E-6 to emphasize how poorly this simple model
describes the two trend case.
213
1000
Ww
wwe
SN NWN
Depth in Fathoms
o
SLOPE OF SPECTRUM( »)
INTERCEPT (2)
0 20 40 60 80 100 120 140 160 180
AZIMUTH (degrees)
Figure E-6 Graphic representation of the Mendocino Fracture Zone
bathymetry shown in Figure E-5 is shown above (viewpoint
from the southwest). Shown below are the spectral estimates
versus azimuth for this surface. Compare this example to
the artificially generated surfaces shown in Figures E-3 and
E-4.
214
ee 4 ee *
cantly MOND
PAA EW
gadenr teh
I Yh
DESCRIPTION OF TERMS
This informal glossary is included to aid readers who represent
diverse backgrounds in geology, geophysics, acoustics, statistics, and
other fields. The descriptions of terms, phrases, and acronyms included
are intended to reflect their use within this report, rather than a
general or rigorous definition.
Amplitude
The departure of a periodic function from its mean.
Anisotropy
Condition of having different properties in different directions.
Band-1imited
Condition in which a signal contains a finite range of frequency.
With reference to Fourier analysis, this band ranges from the fun-
damental frequency to the Nyquist frequency.
Boxcar function
A function which equals unity over some finite length and zero
everywhere else. The multiplication of this function with an infi-
nite series represents mathematically the sampling of a finite
series from an infinite process. Sometimes called a rectangle
function.
Chi-square distribution
The probability distribution for the estimation orror of the
power spectrum. Due to the form of the distribution, the errors
associated with the amplitude spectrum appear constant when plotted
on log-transformed axes.
Coherent signal
A signal in which the phase relationship of the various frequency
components is retained.
Convolution —
A mathematical operation which is equivalent to multiplication in
the opposite transform domain.
Deterministic model
A numerical model in which the mathematical parameters represent
specific measurable quantities of the process under study.
Ensemble average
The mean value of a group of repetitive samplings for a given sta-
tistical measure.
FFT
Fast Fourier Transform. A numerical algorithm for estimating the
Fourier transform with a greatly reduced number of operations.
215
Fractal dimension
A topological term which refers to a dimension which which may be
either integer or fractional. The fractal dimension (D) is related
to the spectral slope (b) by b=-(5/2-D). Refer to Mandelbrot
(1982) for a complete discussion.
Functional form
Refers to the type of function or functions used as a model for the
distribution of data in a regression analysis. The term does not
apply to the calculated parameters used to describe a specific data
set.
Fundamental frequency
The lowest frequency treated in a Fourier analysis, corresponding
to the inverse of the length of data.
HEBBLE
High Energy Benthic Boundary Layer Experiment
Isotropy
Condition of having the same properties in all directions.
Leakage
In spectral analysis, the transfer of energy from one frequency
into other frequency bands.
Linear-linear space
A two-dimensional coordinate system in which both orthogonal axes
represent simple evenly-spaced scales. Usually called a rec-
tangular Cartesian coordinate system.
Log-log space
A two-dimensional coordinate system in which both orthogonal axes
are scaled by evenly spacing the logarithm of the linear scale. In
this study, all such transformations use base-ten logarithms.
Magnetic anomaly ;
The départure of the measured magnetic field from some low fre-
quency model.
Markov process
A stochastic process in which the conditional probability state is
unaffected by the historical state of the system.
Multibeam sonar
A bathymetric sounding system in which several discrete soundings
of the sea floor can be derived by a single discharge of acoustic
energy.
NAVOCEANO
United States Naval Oceanographic Office
Non-parametric statistics
Statistical theory in which the probability distribution of the
underlying data is not assumed.
216
Nyquist frequency i
The highest frequency treated in a Fourier analysis, corresponding
to the inverse of twice the data spacing.
Orthogonal
An orientation in which all axes intersect perpendicularly.
Phase
The location of a periodic function along an axis relative to
some arbitrary origin.
Planetary Rossby waves
A large-scale, stable wave motion in the global ocean.
Power
The squared amplitude of a periodic function.
Power law
A mathematical relationship of the form y=ax), This function
plots linearly in log-log space.
Prewhitening
A technique which reduces the effect of spectral leakage in ana-
lyzing non-white-noise signals.
Process (geological)
A natural continuing activity or function.
Process (statistical)
Any quantity which is defined in terms of its relationship to some
independent variable, usually time or space.
Provincing :
Techniques which divide a sample space into discrete regions based
on some predefined statistical property.
Random walk
A stochastic process of independent increments.
Regression model
A functional description of a relationship between a dependent
variable and independent variable(s), usually derived by the method
of least-squares.
Round-off error
The level of uncertainty in a data set due to the finite number of
digits retained.
SASS
Sonar Array Sub-System. A multibeam sonar system operated by the
U.S. Navy.
SEABEAM
A commercially available multibeam sonar system.
217
SEAMARC-1
A commercially available deep-towed geophysical instrument package.
Sedimentary process
Relief-forming process which involves the transportation of
suspended material.
Sinc function
The function y = sin(x)/x.
Sinusoid
A function having the form of a sine or cosine function, but with a
variable phase value.
Spatial frequency
The inverse of wavelength.
Spectral intercept
The intersection of an amplitude spectrum with the amplitude axis.
This term is specific to this report.
Spectral slope
The slope of an amplitude spectrum which has been plotted on log-
log axes. This value becomes the exponent of spatial frequency
following transformation to linear-linear space. This term is spe-
cific to this report.
Spreading center
A region of the sea floor where recent crustal material is being
formed.
Stationarity
Although rigorously defined in statistical theory, used in this
report to imply relative constancy of a particular statistical
parameter as a function of position; statistically homogeneous.
Stereo-pair bottom photography
A photogrammetric technique which allows the measurement of micro-
relief on the sea floor by analyzing photographic images of a sur-
face from offset viewpoints.
Stochastic model
A numerical model in which the mathematical parameters describe the
process under study in terms of its random variability.
Strike
The bearing of the long axis of a linear trend.
Tectonic process
Relief-forming process which involves the deformation of the
earth's crust.
White noise
A random series whose amplitude spectrum is constant with fre-
quency.
218
DISTRIBUTION LIST
CNO (OP-952, OP-212, PME-124)
ASSTSECNAV RES
ONR
COMNAVOCEANCOM
NRL (B. Adams, D. Berman, R. Feden)
USNA
NAVPGSCOL
NAVSWC
NORDA (110, 200, 220, 240, 250, 250B, 260,
270, 300, 340, 350, 360, 361, 362)
PrePNNMWN ON.W
—
NOSC
NUSC (Newport, RI)
NUSC (New London, CN-T. Bell, J. Dobler,
J. Schumacher )
NOAA/PMEL (E. Bernard, R. Burns, S. Hammond,
R. Embley)
NOAA/NGDC (T. Holcombe)
UT/ARL (H. Boehm, P. Widmar)
JHU /APL
PSU/ARL
UW/APL (D. Jackson, D. Winebrenner
UW/DGS (A. Nowell)
CU/LDGO (D. Hayes, A. Watts, J. Weissel,
S. Lewis, W, Ryan, D. Martinson
G. Mountain)
Cornell (D. Turcotte)
UC/SIO
DMA/HTC
WHOT
UH/HIG
UH/HIG (P. Taylor)
- Wr
NP NMP eNE
Pee PP pe
TR 279
Description, Analysis, and Prediction of
Sea-Floor Roughness Using Spectral Models
April 1985
Live Music Archive
Librivox Free Audio
Metropolitan Museum
Cleveland Museum of Art
Internet Arcade
Console Living Room
Books to Borrow
Open Library
TV News
Understanding 9/11