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 <\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 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)-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 (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 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 (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, | [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