- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Measurements of low frequency acoustic backscatter...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Measurements of low frequency acoustic backscatter from the sea surface 1991
pdf
Page Metadata
Item Metadata
Title | Measurements of low frequency acoustic backscatter from the sea surface |
Creator |
Hill, Steven |
Publisher | University of British Columbia |
Date Created | 2011-01-31 |
Date Issued | 2011-01-31 |
Date | 1991 |
Description | The overall objective of this thesis was to predict, model and measure low frequency acoustic backscatter from the sea surface zone (SSZ). In particular, the objectives were fourfold: to relate the acoustic backscatter Doppler spectrum to the directional waveheight spectrum (DWS) through a perturbation analysis; to develop instrumentation suitable for measuring the properties of acoustic backscatter from the SSZ; to design and implement signal processing hardware and software to process raw data from the instrument; and to deploy the instrument and make measurements to test the validity of the predictions of the theoretical development. A theoretical framework was developed to enable a test of the acoustic analogue of the Coastal Ocean Dynamics Applications Radar (CODAR) technique, using beamforming techniques to simulate the CODAR antennas. Expressions relating the CODAR antenna outputs to the output of an array of omnidirectional acoustic point sensors were developed, and mathematical algorithms and techniques were derived to extract information about the DWS of surface gravity waves from acoustic Doppler backscatter measurements with the array. Models were developed and implemented, showing the expected form of the power spectral density of the acoustic Doppler backscatter seen by single omnidirectional receivers, and the expected form of data products of the beamformed array. An acoustic instrument — the Upward-Looking Sonar Array System (ULSAS) — was developed for stand-alone, remotely controlled operation in both bottom-situated and deep-water, surface-tethered configurations. This device can collect and store large quantities of acoustic data from a multi-element array, under the control of a distant operator over a radio link. The bottom-situated version was deployed in the coastal waters of British Columbia, and the deep water version was deployed in the recent Surface Wave Processes (SWAPP) experiment. A preliminary test of the acoustic CODAR technique was made, yielding information consistent with the known wind and wave field. The form of the non-directional part of the extracted DWS followed approximately the expected k⁻⁴ shape for k values above saturation. Beamforming results using frequency-domain data show that the Doppler-shifted acoustic backscatter is directional in nature. These are the first results of this kind to be reported. The deep-water version of ULSAS was tested for the first time during the SWAPP cruise. In spite of a problem limiting the power output of the projector, estimates of the surface scattering strength parameter over angles of incidence less than 45° were made, showing some surprising departures from the Chapman-Harris empirical formula for S₅ , and interesting angular structure. Measurements of the ambient noise field were also made under calm conditions and during 14 kt winds. |
Subject |
Backscattering -- Measurement Acoustic Surface Waves -- Measurement |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | Eng |
Collection |
Retrospective Theses and Dissertations, 1919-2007 |
Series | UBC Retrospective Theses Digitization Project |
Date Available | 2011-01-31 |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0053270 |
Degree |
Doctor of Philosophy - PhD |
Program |
Oceanography |
Affiliation |
Science, Faculty of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
URI | http://hdl.handle.net/2429/31022 |
Aggregated Source Repository | DSpace |
Digital Resource Original Record | https://open.library.ubc.ca/collections/831/items/1.0053270/source |
Download
- Media
- UBC_1991_A1 H54.pdf [ 8.33MB ]
- UBC_1991_A1 H54.pdf
- Metadata
- JSON: 1.0053270.json
- JSON-LD: 1.0053270+ld.json
- RDF/XML (Pretty): 1.0053270.xml
- RDF/JSON: 1.0053270+rdf.json
- Turtle: 1.0053270+rdf-turtle.txt
- N-Triples: 1.0053270+rdf-ntriples.txt
- Citation
- 1.0053270.ris
Full Text
M E A S U R E M E N T S O F L O W F R E Q U E N C Y A C O U S T I C B A C K S C A T T E R F R O M T H E S E A S U R F A C E By Steven Hill B. Sc. (Physics) University of British Columbia M. Sc. (Marine Biology) University of Victoria A T H E S I S S U B M I T T E D IN P A R T I A L F U L F I L L M E N T O F T H E R E Q U I R E M E N T S F O R T H E D E G R E E O F D O C T O R OF P H I L O S O P H Y in T H E F A C U L T Y O F G R A D U A T E S T U D I E S D E P A R T M E N T O F O C E A N O G R A P H Y We accept this thesis as conforming to the required standard T H E U N I V E R S I T Y O F B R I T I S H C O L U M B I A December 1990 © Steven Hill, 1991 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department The University of British Columbia Vancouver, Canada DE-6 (2/88) Abstract The overall objective of this thesis was to predict, model and measure low frequency acoustic backscatter from the sea surface zone (SSZ). In particular, the objectives were fourfold: to relate the acoustic backscatter Doppler spectrum to the directional waveheight spectrum (DWS) through a perturbation analysis; to develop instrumentation suitable for measuring the properties of acoustic backscatter from the SSZ; to design and implement signal processing hardware and software to process raw data from the instrument; and to deploy the instrument and make measurements to test the validity of the predictions of the theoretical development. A theoretical framework was developed to enable a test of the acoustic analogue of the Coastal Ocean Dynamics Applications Radar (CODAR) technique, using beamforming techniques to simulate the CODAR antennas. Expressions relating the CODAR antenna outputs to the output of an array of omnidirectional acoustic point sensors were developed, and mathematical algorithms and techniques were derived to extract information about the DWS of surface gravity waves from acoustic Doppler backscatter measurements with the array. Models were developed and implemented, showing the expected form of the power spectral density of the acoustic Doppler backscatter seen by single omnidirectional receivers, and the expected form of data products of the beamformed array. An acoustic instrument — the Upward-Looking Sonar Array System (ULSAS) — was developed for stand-alone, remotely controlled operation in both bottom-situated and deep-water, surface-tethered configurations. This device can collect and store large quantities of acoustic data from a multi-element array, under the control of a distant operator over a radio link. The bottom-situated version was deployed in the coastal waters of British Columbia, and the deep water version was deployed in the recent Surface Wave Processes (SWAPP) experiment. A preliminary test of the acoustic CODAR technique was made, yielding information consistent with the known wind and wave field. The form of the non-directional part of the extracted DWS followed approximately the expected k~4 shape for k values above saturation. Beamforming results using frequency-domain data show that the Doppler-shifted acoustic backscatter is directional in nature. These are the first results of this kind to be reported. The deep-water version of ULSAS was tested for the first time during the SWAPP cruise. In spite n of a problem limiting the power output of the projector, estimates of the surface scattering strength parameter over angles of incidence less than 45° were made, showing some surprising departures from the Chapman-Harris empirical formula for S, , and interesting angular structure. Measurements of the ambient noise field were also made under calm conditions and during 14 kt winds. in Table of Contents Abstract ii List of Tables ix List of Figures xi Acknowledgement xv 1 Introduction 1 1.1 M o t i v a t i o n 1 1.2 Object ives 2 1.3 Previous Work Us ing Acous t ic Techniques to Measure Sea Surface Propert ies 3 1.4 R A D A R Methods 4 1.4.1 T h e C O D A R sys tem 4 1.5 The A c o u s t i c C O D A R 5 1.6 Ou t l ine of the Thesis 5 1.7 Summary of the contents of this Thesis 6 1.7.1 S tandard mathemat ica l notat ion 6 2 Acoustic CODAR Theory 7 2.1 Region of V a l i d i t y of the Theory 7 2.2 Geomet ry used i n the Ana lys i s , and the Surface Height Descr ip t ion 8 2.2.1 Geomet ry 8 2.2.2 Surface height descript ion 10 2.3 T h e Acous t i c Wave 10 2.4 T h e Per tu rba t ion Expans ion 12 2.4.1 Firs t -order terms 13 2.4.2 Second-order terms 13 2.4.3 Expressions for the first and second order acoustic pressure 15 iv 2.5 The Acous t i c Doppler Backscat ter Spec t rum 16 2.5.1 Scat ter ing direct ion and the wavevector k2 16 2.5.2 Reduc t ion to a backscatter spec t rum 17 2.6 T h e Acous t i c Backscat ter Spec t rum in terms of the Waveheight Di rec t iona l Spec t rum. . . 18 2.6.1 F i r s t order 18 2.6.2 Second order 19 2.7 N o r m a l i z a t i o n 22 2.7.1 F i r s t order 22 2.7.2 Second order 23 3 The Directional Waveheight Spectrum from Acoustic Doppler Backscatter 25 3.1 Beamforming 25 3.1.1 T h e 4 element sub-array 25 3.2 D a t a P roduc t s - the b„ Defined 29 3.2.1 W r i t i n g the bn in terms of cross-spectra 30 3.3 Re la t ing the bn(r). 6) to the Di rec t iona l Waveheight Spec t rum 32 3.3.1 T h e inverse problem 32 3.3.2 M o d e l l i n g the direct ional waveheight spect rum 33 3.3.3 T h e first-order bn 33 3.3.4 T h e second-order bn 34 3.4 So lv ing the Integral Equa t ion for Second Order 6 n 37 3.4.1 T h e constant bands approximat ion 37 3.4.2 M a t r i x form of the equat ion 39 3.5 So lv ing the M a t r i x Equa t ion 39 3.5.1 T h e straightforward approach 39 3.5.2 Linear least-squares theory 40 3.5.3 A d d i n g constraints to the model 41 3.5.4 C o m p u t i n g the covariance ma t r ix 45 4 Modelling of Acoustic Doppler Backscatter 47 4.1 Omnid i r ec t i ona l spectra using the P h i l l i p s wave spec t rum 47 4.1.1 So lv ing the integral 49 v 4.1.2 The first-order spectrum 50 4.1.3 Smearing 50 4.1.4 Results for the omnidirectional receiver (Phillips spectrum) 52 4.2 Model bn using the Phillips Spectrum 53 4.2.1 Some examples of computed bn functions 54 4.3 Omnidirectional spectra using the JONSWAP non-directional model and a variable direc- tional factor 54 4.3.1 Description of the directional spectrum used in this model 57 4.3.2 The directional factor 58 4.3.3 The total directional spectrum 60 4.3.4 Computations of typical spectra, and comparisons to the truncated Fourier series representation 60 4.3.5 Omnidirectional acoustic doppler spectra developed using the directional JON- SWAP model 65 4.4 Application of the JONSWAP/directional wave spectrum model to the computation of model bn functions 68 4.5 Tests of the inversion methods using model input 68 4.5.1 5 m/s windspeed 68 5 Instrument Development 74 5.1 ULSAS I 74 5.2 ULSAS II 74 5.2.1 Overview of system operation 75 5.3 Measurements with ULSAS I and II 77 6 Saanich Inlet Measurements. 79 6.1 Typical power spectrum 79 6.2 Computed 6n functions 82 6.3 Tests of inversion methods on real acoustic data 82 6.3.1 The straightforward method 82 6.3.2 Using iterative improvement and the "nastiness" function 85 6.4 Beamforming in the frequency domain 86 vi 6.5 Discussion of Saanich Inlet results 90 6.5.1 Subsequent deployments and instrument development 91 7 SWAPP Measurements 93 7.1 Selection of data for processing 93 7.2 Environmental sensors . . 96 7.3 Typical hydrophone data 96 7.4 Surface scattering strength parameter (5 S) variation with grazing angle 100 7.4.1 Corrections for depth variations 102 7.4.2 Mean squared pressure computations 102 7.4.3 Backscatter measurements using correlation analysis 104 7.4.4 Backscattering measurements using both methods 107 7.5 Ping by ping time series 109 7.5.1 Voltage statistics, and returned power levels 110 7.6 Frequency spectra of Ping data and ambient noise 112 7.6.1 Ambient noise spectra 113 8 Summary and Conclusions 115 Appendices 118 A Computation of F7 using "REDUCE" 118 B Expressions for the data products B„. 119 C A Detailed Description of ULSAS II 120 C . l General system description 120 C . l . l Mechanical description of the deep water version of U L S A S II 120 C. l .2 Deployment and recovery 124 C.2 Generation of the Acoustic Signal 125 C.2.1 Sine wave generation 125 C.2.2 Analog processing 128 C.2.3 Acoustic projector and power amplifier 128 C.3 Acoustic Receiving System 128 vii C.3.1 Receiving hydrophones 128 C.3.2 Analog signal processing. 129 C.4 Data Aquisition and Storage 129 C.4.1 Sample/hold and multiplexer 129 C.4.2 Digital audio ( P C M ) system 130 C.4.3 P C M Clocks 132 C.4.4 The Video recorder 132 C.4.5 Data sampling configuration • 133 C.5 Communication and Control. 134 C.5.1 Communications protocol 135 C.6 Local Microcontrollers 139 C.6.1 Main controller 139 C.6.2 Ping/Sample controllers 141 C.7 Power Supplies 141 C. S Environmental Sensors 142 C.8.1 Azimuth 142 C.S.2 Pressure 142 C.8.3 Acceleration 142 C. 8.4 Vertical tilt 143 D Pre-SWAPP calibrations of ULSAS II 144 D. l A D C System Calibration 144 D. l . l Results 144 D.2 Calculation of source level of Argotech projector 145 D. 3 Instrument problems during the S W A P P cruise 150 E Details on data processing. 152 E . l Prehminary Data Processing 152 E . l . l Hardware 152 E.1.2 Software 155 Bibliography 156 V l l l List of Tables 4.1 Expected dj(qi) values for 5 m/s windspeed and principal direction of 35° 71 4.2 Results of straighforward inversion of bn models for 5 m/s windspeed and principal direc- tion of 35° 71 4.3 Estimates of the principal direction and the factor s from results of inversion of model data 72 4.4 Results of straighforward inversion of 6n models for 5 m/s windspeed and principal direc- tion of 35°. Noise at SNL of 20 dB added 73 4.5 Estimates of the principal direction and the factor s from results of inversion of model data. Noise at SNL of 20 dB added 73 6.6 Results of straighforward inversion of real bn functions. A co-elevation angle of 45° was used 84 6.7 Estimates of the principal direction and the factor s from the results in Table 6.6 85 6.8 Results of straighforward inversion of real 6„ functions. A co-elevation angle of 31° was used 85 6.9 Estimates of the principal direction and the factor s from the results in Table 6.8 86 6.10 Results of inversion of real bn functions using iterative improvement. A co-elevation angle of 45° was used 86 6.11 Estimates of the principal direction and the factor s from the results in Table 6.10. . . . 87 6.12 Summarized results of 4-element beamforming. 6max is the mean direction from which maximum power was returned, and 0 m i n is the mean direction for minimum power. . . . 90 7.13 Statistics of ULSAS environmental sensors for deployment at 06:00 on 17/3. N = 376, covering a time of 6.2 minutes 96 7.14 Incidence angle and Bragg frequency for each range gate for a typical collection starting at cycle 30, and a carrier frequency of 393.75 Hz 110 C.15 Counter reload values and corresponding frequencies 127 ix C.16 Signal gain as a function of attenuator value 128 C.l7 MUX channel assignments 131 C.18 Sampling frequency vs 5/ setting 136 C.19 Amplifier gain versus RA setting 137 C. 20 68HC11 MUX channel assignments 140 D. 21 Calibration voltage input for different gain settings 144 D.22 Calibration measurements for different gain settings (H = high level, L = low level input). 145 D.23 Calibration offsets and factors for the 4 available instrumentation amplifier gains 145 D.24 Parameter estimates for power amp input attenuation setting of 64 147 D.25 Parameter estimates for power amp input attenuation setting of 32 148 D.26 Parameter estimates for power amp input attenuation setting of 16 148 D.27 Parameter estimates for power amp input attenuation setting of 8 149 D.28 Parameter estimates for power amp input attenuation of 08H 149 x List of Figures 2.1 Geometry of the experiment, looking down on the sea surface. Plane waves are incident from below in the x direction at co-elevation angle 6. The scattered wave (k,) is di- rected downward with co-elevation angle 9S and azimuthal angle ips. k2 is the horizontal projection of k J : and k ; is its vertical component 9 3.2 4 element hydrophone array. The reference direction is through H4 to H2 26 3.3 Dipole formed by two hydrophones 27 3.4 Wavevectors q and qi in r — s space. Curves of constant r) are plotted for the inner sidebands. The curves are symmetric about the r axis 36 3.5 Bands of constant dj(q) (dashed lines) and constant r/ (solid lines) for the outer Doppler sidebands in (r, s) space. The curves are symmetrical about both the r and s axes. . . . 38 4.6 One-sided power spectra as seen by omnidirectional receivers. The dashed line is for windspeed of 10 m/s, and the solid line is for 5 m/s 52 4.7 One-sided power spectra as seen by omnidirectional receivers. Both are calculated using a windspeed of 5 m/s: the dashed line shows the spectrum for a co-elevation angle of 60°; the solid line is the spectrum using a co-elevation angle of 30° 53 4.8 Computed bn functions using the Phillips spectrum. Windspeed is 10 m/s, co-elevation angle is 45°, and principal wave direction is 137°. The horizontal axis is the normalized Doppler frequency 55 4.9 Computed b„ functions using the Phillips spectrum. Windspeed is 10 m/s, co-elevation angle is 45°, and principal wave direction is 210°. The horizontal axis is the normalized Doppler frequency 56 4.10 Fetch (km) around instrument site in Saanich Inlet 59 4.11 Mean depth (m) around instrument site in Saanich Inlet 60 4.12 Family of JONSWAP spectra for different wind directions. Dashed line is Phillips spec- trum for the same windspeed 61 xi 1 4.13 Directional JONSWAP spectra (solid line) together with TFS representation(dashed line) for k values 0.02 through 0.05 62 4.14 Directional JONSWAP spectra (solid line) together with TFS representation (dashed line) for lc values 0.06 through 0.20 63 4.15 Directional JONSWAP spectra (solid line) together with TFS representation (dashed line) for k = 0.30 64 4.16 Power spectra as seen by omnidirectional receivers. JONSWAP directional model used, with 0.1 Hz smearing width 65 4.17 Power spectra as seen by omnidirectional receivers. JONSWAP directional model used, with 0.2 Hz smearing width 66 4.18 Power spectra as seen by omnidirectional receivers. JONSWAP directional model used, with 0.15 Hz smearing width, and noise at 10 dB SNR added 66 4.19 Power spectra as seen by omnidirectional receivers. JONSWAP directional model used, with cos4-proportional directional factor. 67 4.20 Power spectra as seen by omnidirectional receivers. JONSWAP directional model used, with coss-proportional directional factor 67 4.21 Model 6„ functions using the JONSWAP/directional spectrum. Winds to 160° (long fetch). The horizontal axis is the normalized Doppler frequency 69 4.22 Model b„ functions using the JONSWAP/directional spectrum. Winds to 225° (short fetch). The horizontal axis is the normalized Doppler frequency 70 5.23 Side view (not to scale) of ULSAS I deployment in Saanich Inlet 78 6.24 Average of 37 windowed, 512 point power spectra for hydrophone 5. Horizontal axis is Doppler frequency (Hz). Expected Bragg frequency regions are between the two sets of dotted lines 80 6.25 Computed bn functions for one four-element subset of hydrophones. The horizontal axis is the normalized Doppler frequency v 83 6.26 4-eIement subarray directional response for the band -.95 Hz to -.80 Hz 88 6.27 4-element subarray directional response for the band -.75 Hz to -.55 Hz 88 6.28 4-element subarray directional response for the band .55 Hz to .75 Hz 89 6.29 4-element subarray directional response for the band .80 Hz to .95 Hz 89 xn C.30 Backscatter from near-vertical incidence to 30° grazing angle, for ten consecutive pings. Hydrophone 6, Georgia Strait deployment 92 7.31 Weather information from SAIL system during ULSAS deployments on the SWAPP cruise. Numbers and horizontal bars at the top of the figure refer to ULSAS deployments. . . . 94 7.32 Weather information from SAIL system during deployment 7 95 7.33 ULSAS environmental sensor data (Compass, X-tilt and Y-tilt) for a 6 minute period at 06:00 17/3 (Deployment 7) 97 7.34 ULSAS environmental sensor data (Pressure and Accelerometer) for a 6 minute period at 06:00 17/3 (Deployment 7) 98 7.35 Typical transducer voltage time series during one ping 98 7.36 Transducer voltage time series of multiple pings covering a 10 s period. Note the slow variation of the mean level 99 7.37 12 overlapped transducer voltage time series from consecutive pings, for returns from near vertical to grazing angles of 30° 100 7.38 Geometry of the scattering strength parameter (S,) measurements. Z is the depth of the pressure sensor, and the projector and hydrophone arrays are 2 m apart 101 7.39 Average amplitude of mean squared pressure (solid line) and correlation magnitude (bro- ken line) vs time after beginning of ping for one hydrophone. These data are averages of 1200 pings 103 7.40 Computed values of Sj as a function of 6 using mean squared pressure for hydrophones 2,5,8 and 9 at 00:35 on 17/3/90. Windspeed was 8 m/s. In this and the next 4 figures, the dashed line is the predicted value for Ss using the Chapman-Harris model 104 7.41 Computed Ss values for data from 00:35 17/3. Results of mean squared pressure compu- tations on left, correlation analysis on right 108 7.42 Computed Ss values for data from 06:00 17/3. Results of mean squared pressure compu- tations on left, correlation analysis on right 108 7.43 Computed Ss values for data from 22:00 17/3. Results of mean squared pressure compu- tations on left, correlation analysis on right 109 7.44 Probability plots of 200 consecutive values of real and imaginary voltages for gate 2 (0 = 31.1°), referred to a normal distribution I l l xm 7.45 Mean amplitude of complex voltage for 5 range gates, plotted against 0. The dashed line is approximately the amplitude of ambient noise 112 7.46 Frequency spectra for hydrophone No. 7 for two different regions of the ping. Left plot is cycles 25-57; right plot is cycles 37-69 113 7.47 Frequency spectrum of ambient noise for hydrophone No. 5. Data are averages of 1000 unwindowed 128-point series. Left-hand plot is from File A7F22 (calm); right-hand plot is for File A8F05 (14 let winds) 114 C.48 Side view of ULSAS deployment 121 C.49 Side view of the hydrophone/projector mount 122 C.50 Side view of the electronics/battery case assembly 123 C.51 "Black-box" representation of ULSAS II electronic subsystems 125 C.52 Integer representation of sine wave (upper plot), and its power spectrum (lower plot). The horizontal axis of the lower plot is proportional to frequency, where the fundamental frequency of the sinusoid is at index 1 126 C.53 Typical transmitting voltage response of the Argotech projector 129 C. 54 Typical frequency response of the pre-amps, covering the range 1 Hz to 50 kHz. Right- hand plot is a blow-up of the region between 100 Hz and 1 kHz. Note the different vertical scales in these plots 130 D. 55 Example of data used in non-linear fit to sinusoid, together with fitted curve 146 D. 56 Residuals after non-linear fit to sinusoid 147 E. 57 Data processing hardware components and interconnections 153 xiv Acknowledgement This work was done in conjunction with the Ocean Acoustics group at the Institute of Ocean Sciences, Sidney, B.C. The Canadian Panel on Energy Research and Development supported some of the work under project number 67132, and some of the field work was supported by the U.S. Office of Naval Re- search as part of the Program on Surface Reverberation. I have been supported personally by the Science Council of B.C. through a G.R.E.A.T. scholarship, and by a National Research Council Postgraduate Scholarship. 1 would like to thank Dr. David Farmer for his unswerving support and enthusiasm through a very long project, and Dr. Stephen Pond for his timely and thoughtful advice. In addition, the personnel at Arctic Sciences Ltd., particularly David Lemon, were extremely helpful during the early part of this project. Finally, 1 would like to thank my wife Marni, and my children Annwyn and Patrick for their patience, understanding and good cheer through a longer than anticipated return to school. xv Chapter 1 Introduction This thesis describes theory and modelling of low frequency (LF) acoustic backscatter from the sea surface, the development and testing of instrumentation to measure acoustic backscatter, and the devel- opment and testing of algorithms used to infer properties of the sea surface from such measurements. 1.1 Motivation LF ocean acoustics has become more important in recent years with the development of such tools as large-aperture towed arrays, and new lightweight piezoelectric materials which make ship-borne large arrays practical and possible. Interest in LF acoustics is driven by many disparate factors, including military requirements, i.e. active sonar, seismic exploration of the sea-bottom, and the detection of nuclear tests. In addition, LF acoustics is of great interest to oceanographers, because of the its potential as a remote sensing tool. A clear understanding of the processes involved in LF acoustic scattering from the sea surface is essential to the proper design and operation of sonar systems, and to the validity and usefulness of LF reflection and propagation models. Conversely, once the mechanisms are well understood, the process can be turned around, and the characteristics of scattered acoustic signals can be used to infer various properties of the sea surface. In particular, measurements at low frequencies (i.e. a few hundred Hz) and at windspeeds greater than about 15 kt have shown that acoustic backscatter from the sea surface is not well described by existing models or physical understanding of the scattering process. It has been suggested that the sea surface, from an acoustic point of view, does not consist simply of the air-water boundary, but rather a zone, the Sea Surface Zone (SSZ), which extends some meters down into the water column from the actual air-water interface. Many physical processes which may have effects on the scattering and reflection of acoustic waves take place in this zone. These processes include the transport of entrapped air in the form of bubbles due to breaking waves, turbulence associated with vertical velocities caused by surface wave motion or breaking internal waves, and small scale circulation such as Langmuir cells. 1 Chapter 1. Introduction 2 One approach to understanding the acoustical effects of these mechanisms is to concentrate on one mechanism which we might be able to understand and model very realistically, such as acoustical scatter- ing from the air-water boundary. Then, discrepancies between measurements and predictions might be ascribed to the other mechanisms, which are poorly understood and even more poorly measured in the SSZ. The instrument which is described in this thesis is the Upward-Looking Sonar System (ULSAS), which was designed first as a tool to infer the directional waveheight spectrum(DWS) of surface gravity waves from measurements of LF acoustic backscatter. In addition, ULSAS is an excellent instrument for making measurements of Doppler backscatter spectra, the directional properties of LF backscatter from the sea surface and of ambient noise in the LF region. Although considerable work on the modelling of acoustic scattering has been done, there are very few measurements of LF scatter in the literature. ULSAS was conceived and designed to make some of these important measurements, and to test the validity of certain models of the scattering process. 1.2 Objectives In general, the objectives of this thesis are to predict, model, and measure LF acoustic backscatter from the SSZ. In particular, there are four objectives. The first objective is to derive expressions for the first and second order acoustic Doppler backscatter using a perturbation analysis, and to relate the Doppler backscatter spectrum to the DWS of surface gravity waves. This derivation follows methods developed previously in the RADAR field (described further below). Models of the expected form of the Doppler backscatter spectrum are derived using the perturbation theory and existing models of the DWS, and algorithms and processes for inferring the DWS from the acoustic Doppler backscatter spectrum are developed. The second objective is to design and develop instrumentation suitable for measurement of acoustic Doppler backscatter, as well as directional and other properties of both acoustic backscatter from the SSZ and ambient noise. The third objective is to deploy instruments and collect data suitable for testing the validity of the models and techniques developed under the first objective, and for making measurements of general properties of acoustic backscatter from the SSZ and of ambient noise. The fourth objective is to design and implement signal processing hardware and software to process raw data collected by the instrument so that tests of the theory can be carried out, and measurements Chapter 1. Introduction 3 of the properties of acoustic backscatter from the SSZ can be made. 1.3 Previous Work Using Acoustic Techniques to Measure Sea Surface Properties. Although some study of acoustic scattering from the surface went on during the Second World War, when "surface reverberation" was of great interest to sonar engineers and operators, Eckart, [16] in 1953 was the first to apply rigorous mathematical analysis to the surface scattering problem. Using Helniholtz's theorem and the Kirchhoff approximation at the air/sea boundary, he derived some of the expected properties of monochromatic sound waves reflected from the ocean surface. A fundamental conclusion of his work was that much more information about the sea surface would be revealed by long- wave scattering (i.e. acoustic wavelengths much greater than the rnis waveheight) than by short wave scattering. Following Eckart, many other workers have made contributions, some of which have been reviewed by Fortuin [17]. Although Fortuin points out that "the backscattering contains all statistical information about the surface", most measurements reported in the literature have been of forward- scattered continuous-wave(CW) sound: probably because of the increased signal to noise ratio for this configuration, and because most theoretical models were of bistatic (receiver and transmitter separated by a significant distance) configurations looking at forward scattered sound. Examples of these kinds of measurements include work by Williams [59], who also gives a concise review of related work; and Roderick and Cron [47], who report on both laboratory and ocean measurements. Scrimger [48] has written an exhaustive review of the literature (up to 1983) on the measurement of surface gravity wave properties using forward scattered CW sound. There is a large body of theoretical work on scattering (both acoustic and electromagnetic) from rough surfaces. Workers at Bell Labs have been involved in this research, notably Harper and Labianca, who have developed a general perturbation theory for long-wave scatter from a moving rough surface [27]: and applied inverse scattering methods to show how the directional waveheight spectrum might be estimated using a bistatic arrangement [26]. DeSanto and Brown [15] carried out a thorough review of analytical techniques used in the analysis of multiple scattering from rough surfaces. Recently, Kuo [25] and Mellen [40] have presented modelling results for sea-surface scattering and propagation loss based on perturbation methods. Experimental results continue to be reported in the literature up to the present time, although very few results are available at LF: some interesting measurements were reported by Lynch et al [38], Chapter 1. Introduction 4 who derived estimates of the non-directional waveheight spectrum using the phase fluctuation spectra of surface interacting rays observed during a preliminary ocean tomography experiment carried out as part of the Marginal Ice Zone Experiment (MIZEX '84). Pinkel and Smith [44] were able to estimate three-dimensional wavenumber-frequency spectra using Doppler processed returns from a surface grazing sonar. The former experiment involved forward scattered returns from a 224 Hz projector over a roughly 50 km path length, and the latter measured backscattered returns from the sea surface approximately 1 km from a 75 kHz projector. The scattering mechanisms must therefore be different in these two cases: in the former it is direct scattering of pulses from the rough sea surface which causes the observed phase fluctuations; and in the latter the Doppler-shifted backscatter is thought to be due to moving bubble clouds which are present because of wave breaking events. 1.4 R A D A R Methods. Remote sensing of sea state using HF (around 25 MHz) RADAR has a long history, also beginning in the Second World War. Crombie [13] in 1955 was the first to report the connection between Doppler shifted sidebands observed in RADAR sea returns and the ocean wave regime. The field was reviewed in 1972 by Barrick [1], and later in 1978 by Barrick and Lipa [4]. Both bistatic and backscatter methods have been used with some success to derive information about the non-directional gravity wave spectrum. Broad- beam (Barrick and Lipa [3]) and narrow-beam (Lipa ei. al. [34]) measurements have been described. Broche et. al [9] used VHF RADAR to estimate sea states and surface currents, and skywave RADAR (Maresca and Georges [39]) has been used to estimate the non-directional waveheight spectrum at a distance of 1600 km from the transmitter/receiver. 1.4.1 The CODAR system. Over the past decade, considerable work has been carried out using the CODAR (Coastal Ocean Dy- namics RADAR) system. This backscatter technique, fully described by Barrick and Lipa [3],uses a 25 MHz broad-beam crossed loop ground wave Doppler RADAR. Analysis methods for the extraction of wave spectral information from narrow-band RADAR returns were developed by Lipa and Barrick [31] and later extended by them to CODAR data [32]. Recently Wyatt [60] has given some modifications to their technique for narrow-beam RADAR which extend the region of applicability of the method to shorter wavelengths. All of these techniques are applicable, with Chapter 1. Introduction 5 some modifications, to acoustic scattering. Although CODAR has been used to measure wave spectra, there have been very few published results. The most notable success of the CODAR system has been in the field of surface current measurement (Lipa and Barrick [33]), where it has been used to map surface currents over.large areas of ocean. 1.5 The Acoustic CODAR One of the applications of ULSAS is as an "acoustic CODAR", designed to measure the directional waveheight spectrum using the Doppler shifted backscattered returns from the sea surface above the instrument. Short pulses from a low-frequency projector are backscattered from the sea surface and detected by an array of omnidirectional hydrophones. Beamforming techniques are used to simulate a CODAR crossed loop antenna. Acoustic scatter from the sea surface is mathematically simpler than RADAR scatter, since it is a scalar rather than a vector process. However, in the acoustic case there is the added complication of a finite (and usually large) angle between the incident and backscattered waves and the surface. Moreover, the presence of bubbles will influence backscatter at higher windspeeds, as mentioned previously. Thus such a device may serve as a useful probe of the departure of backscatter from predicted levels at higher windspeeds, providing due account is taken of the special angular effects mentioned above. At low windspeeds, the system offers the opportunity of remote wave spectrum measurement. 1.6 Outline of the Thesis. Some of the work in this thesis is based on the work of others, especially the work involving the acoustic CODAR. My contributions to the work are detailed in the following paragraphs. First, I recalculated and verified the form of the first and second order acoustic Doppler backscatter spectrum, following initial work by Barrick [2], incorporating some corrections and modifications. Sec- ond, 1 applied a scheme for deriving the DWS of surface gravity waves from CODAR measurements to the acoustic case, developing processing algorithms and software for processing data from the acoustic CODAR. Third, I developed numerical techniques for the solution of the acoustic CODAR equations with noisy data based on regularization and iterative improvement. Fourth, I developed an extensive array of models of the expected Doppler backscatter spectra and other data products based on the perturbation analysis for first and second order backscatter, using two different models for the DWS of Chapter 1. Introduction 6 surface gravi ty waves. F i f t h , I conceived, designed, and helped b u i l d two new configurations of U L S A S in a stand-alone, remotely operated mode: one configuration for bo t t om deployment; and the second configuration suitable for midwater deployment in the deep open ocean. F ina l l y , I have collected and analysed data f rom U L S A S deployments i n several locations. A l t h o u g h the p r imary deep ocean deploy- ment was par t ia l ly foiled by low transmit ter power, a pre l iminary test of the acoustic C O D A R theory, and measurements of parameters of interest in acoustic backscatter and L F ambient noise were carr ied out. 1.7 Summary of the contents of this Thesis In Chapter 2, a per turba t ion analysis of acoustic surface backscatter, which is the theoretical basis of the acoustic C O D A R technique, is developed. In Chap te r 3, the derivat ion of acoustic C O D A R data products is described together wi th the relat ionship between the data products and the direct ional waveheight spect rum, and methods of solving the result ing system of equations. The methods and results of mode l l ing investigations are given in Chapte r 4. T h e instrument development part of this work is summarized i n Chapte r 5. Chapter 6 gives details and interpretat ions of measurements made w i t h U L S A S in the bot tom-mounted configuration in coastal waters, and Chapte r 7 covers deployments and measurements made i n the deep open ocean. Chap te r 8 gives conclusions and recommendations for further work. Deta i led technical informat ion, especially concerning ins t rument development and operat ion, and signal processing hardware and software, can be found in the Appendices . 1.7.1 Standard mathematical notation Throughou t this paper, I w i l l adopt the fol lowing nota t ional conventions: lower case letters (e.g. a) w i l l denote scalars; bold-face lower case letters (e.g. a) w i l l represent vectors; and bold-face upper case letters (e.g. A) w i l l represent matrices. Whenever there are exceptions to these conventions they w i l l be stated expl ic i t ly . Ang le brackets around a quanti ty (e.g. < a >) i m p l y the ensemble average in theory; in practice this notat ion implies the sample average, usually meaning t ime-averaging over the sampl ing per iod. Transposi t ion of ma t r i x quantit ies is denoted by a superscript T ' , e.g. Q T is the transpose of Q. Fina l ly , the operator '* ' means the complex conjugation of a quanti ty (e.g. a* is the complex conjugate of a) . Chapter 2 Acoustic CODAR Theory In this chapter, I will develop the theory used to describe acoustic backscatter from the sea surface. This theory is based on a perturbation expansion of the acoustic pressure field at the sea-air boundary. The technique was first used by Rice [46] to describe the scattering of electromagnetic waves from a "slightly rough" surface, and the development shown here was originally carried out by Barrick [2]. It is included because all of the subsequent development of the CODAR theory rests on it, and because several errors in Barrick's original development are corrected. 2.1 Region of Validity of the Theory. Roughness (or, conversely, smoothness) is most commonly defined by the Rayleigh [45] roughness pa- rameter, which is clearly defined by deSanto and Brown [15]. In this case we can consider the system smooth enough for the perturbation expansion to first order in C to be valid if: 2fcoCcos0 < | , (2.1) where C is the rms waveheight, ko is the acoustic wavenumber, and 9 is defined in Figure 2.1. Because explicit conditions for the validity of this Rayleigh-Rice perturbation theory are unavailable except for a few surfaces with special profiles, there has always been some uncertainty as to when it can be safely used. Recently Jackson and Winebrenner [21] have compared numerical results for electromagnetic scattering from a sinusoidal boundary using both the Rayleigh-Rice perturbation theory and a new theory, said to be more rigorous and exact (Nieto-Vesperinas and Garcia [41]). Their results showed that the two theories gave identical results when carried out to fifteenth order, and were analytically identical to fifth order. Thorsos and Jackson [51] examined the validity of the perturbation approximation in rough surface scattering, comparing numerically computed results with theoretical calculations using the perturbation approximation to first and second order, for a Gaussian roughness spectrum. They found first that, 7 Chapter 2. Acoustic CODAR Theory 8 except for very small values of &oC (~ 0.10) , the perturbation approximation had to be carried out to second order for good agreement between calculated and computed scattering strengths. Secondly, they looked at the effect of the surface correlation length /, defined by: S(k)elkxdk = (C2) f{x) • and/(/) = e"1. (2.2) -cc Here S(k) is the nondirectional waveheight spectrum of the surface, and k is the wavenumber of surface waves. They found that the correlation length has important effects on the accuracy of the perturbation approximation when using a Gaussian roughness spectrum. Even for small &oC> large values of k$l (5.0 or greater) will lead to inaccuracies in the calculations, especially in the backscatter directions. Addressing the first of these problems, it will be seen below that the expansion of the acoustic pressure field is carried out to second order, so that the theory used should be applicable at least within the range of (,' defined by 2.1. Later work by Thorsos [50] using a Pierson-Moskowitz roughness spectrum (a much more realistic approximation to the sea surface than the Gaussian spectrum) shows that second-order perturbation theory is accurate to values of &(," as high as 1.79, far beyond the region of accuracy for the Gaussian spectrum. As to the correlation length restriction, Thorsos [pers. comm.] has stated that this restriction on the validity of perturbation theory is not applicable to realistic (i.e. power law) roughness spectra, but is a feature of the artificial Gaussian roughness spectrum. 2.2 Geometry used in the Analysis, and the Surface Height Description 2.2.1 Geometry The co-ordinate system used in the analysis is shown in Figure 2.1. The origin of the z axis is at the level of the mean ocean surface. Acoustic energy is incident from below and to the left of the z axis, at the co-elevation angle 9. The -fz-axis is the horizontal projection of the direction of propagation of the incident plane wave. An additional set of polar co-ordinates with the -z axis as the polar axis is defined for scattered waves, such that <ps is the azimuthal angle of the scattered wave, and 9S is its co-elevation angle. Vectors k2 and kz are the horizontal and vertical components, respectively, of the scattered wavevector ks. Chapter 2. Acoustic CODAR Theory 9 Figure 2.1: Geometry of the experiment, looking down on the sea surface. Plane waves are incident from below in the x direction at co-elevation angle 6. The scattered wave (k,) is directed downward with co-elevation angle 6S and azimuthal angle ips. k2 is the horizontal projection of ka, and k̂ is its vertical component. Chapter 2. Acoustic CODAR Theory 10 2.2.2 Surface height description Following the usual convention (e.g. Leblond and Huntley [29]), small displacements C(r,0 in the z direction about the mean surface z — 0, where r = (x.y) is the horizontal spatial vector and t is the time, are described as a superposition of an infinite number of plane waves of different lengths, frequencies, amplitudes and directions of propagation: C(r,<) = ^H(k,w)e'( k ' r - a "). (2.3) k,u' In this case the boldface uppercase letter in H(k,w) does not represent a matrix, but the Fourier coefficients of the surface height. Since £ is real, we have: H*(k,w) = H(-k , -w) . Now we hypothesize some fundamental space and time scales L, T and define the wavenumber and frequency through the following: k = (kx,ky) = (am, an) ; u = si ; where m,n, and / are integers, a — — . and Li 2n In this development, we also assume that the waves are propagating in deep water (i.e. the depth of the water is much greater than the wavelength), so that we use the deep-water dispersion relation for surface gravity waves (Leblond and Mysak [30]): w2 = g\k\. (2.4) 2.3 The Acoustic Wave The acoustic pressure wave P(r, z,i) incident on the surface can be modelled as a plane wave which we take, without loss of generality, to be propagating in the x — z plane: Chapter 2. Acoustic CODAR Theory 11 pog'(*o^ sin 8 + k0z cos 8-ui0t) / where and Ao is the acoustic wavelength. PQ is the source amplitude. The wave reflected from a flat pressure release surface at z = 0 is given by: p t(k(jX sin 8 — ICQZ COS 8~ <*>o0 In addition there is a fluctuating component reflected from the surface roughness, which we can also express as an infinite Fourier series: ] T A ( k , w ) e ' ( k - r - ^ - ( a , + w ° : ) < ) k,a. Again the boldface letter in A(k.w) does not represent a matrix, but the Fourier coefficients of the pressure field due to the surface roughness. Summing all these parts, we arrive at an expression for the acoustic pressure at any point (x,y,z) in the medium as: P(r,z,t) = 2iPosm(kozcos0)e<koXSlne-w"^ + ^A(k,u;)e'lkT-k'z-^+^ (2.5) The boundary condition for the acoustic wave at the pressure release surface z — C(r,f) is given by (Clay and Medwin [12]): P(r,C(r,t),0 = 0. (2.6) Finally, the acoustic pressure must satisfy the wave equation c1 where c = LJo/k0; whence, if u <C u>o '• Chapter 2. Acoustic CODAR Theory 12 = iJkT=k% ;fcg<fc2. (2.7) From this point on, the factor e~lWoi will be dropped from the equations, as it does not enter into the calculations further. 2.4 The Perturbation Expansion We can relate H ( k , w ) to A ( k , w ) through a perturbation expansion of A(v,z,i) at z = £(r,<). We assume that the following quantities are small: k0C ; CC Effectively what we are doing in the following is a Taylor series expansion of the boundary condition 2.6 about z = 0, keeping terms to order ( only. That Taylor series expansion can be written as: .dP P(r,Cty= P(T,0,t)+Cdz r,o,« We expand H ( k . w ) and A ( k . w ) to second order as follows: H ( k , w ) = H ( 1 ^ ( k , a ; ) + H(2\k,w) + 0(3) and A ( k , w ) = A^(k,u) + A ^ k . w ) + 0(3) Since both the acoustic and waveheight Fourier coefficients exist only because of the roughness, there is no zeroth order term in either expression. The perturbation expansion is stopped at second order, although it is necessary to go to third order to get all the intensity terms to fourth order. However, since this additional term (a product of first and third order terms) only comes in at the frequency of the first-order Bragg peaks, it is ignored here. Next we let z = C(r,<) in 2.5, using 2.3 for C(r!*) m the boundary condition 2.6. We will also expand some quantities: sin(&oC c o s 0 ) = ôC c o s ^ + 0(3) i a n c i Chapter 2. Acoustic CODAR Theory 13 e(-ikzQ = l-*fc,C + 0(2). Keeping only up to second-order quantities, 2.6 becomes: 2«P0Jfc0 cos6]T (H<X> + H<2>) e ' ( ( k + k ° ) + £ (A' 1 ' + A ^ ) X k,u> k,u> [ 1-ik, ^ H ( 1 ) ( k 1 ) W l ) C ' ( k ' - r - w ' « ) ] e i ( k r - t ) = 0 (2.8) \ ki . w i / 2.4.1 First-order terms. Now we collect terms of like order from 2.8, beginning with the first order: 2 » P o * o c o s f f ^ H f 1 ) e , « k + k ° ) - r - t t ' t > - | - 5 3 A ( 1 ) e , ( k - r - u ' ' ) = 0 (2.9) k,u> k.u' Next we multiply 2.9 by e''^2'1""2'', then integrate the equation with respect to x.y and / over the fundamental scales L and T. This procedure has the effect of removing the summations by selecting out only certain non-zero values of k and u. For example, it can be shown that: J J H ( 1 ) ( k , W ) e , ( ( k + k » - k 2 ) ' r ) d r = H ( 1 ) ( k , W ) ^ where kg = (ko sin 6, 0) is the horizontal projection of ko, and 6^~k° =6(k-(k2-ke)) is the Dirac delta-function. As a result of the integration, 2.9 becomes: A(1>(k2,w2) = -27.Pofcocos0H^(k2 - k » , w 2 ) . (2.10) 2.4.2 Second-order terms. Collecting terms of second order from 2.8 gives: Chapter 2. Acoustic CODAR Theory 14 2«P„*o cos 6 £ H ( 2)(k, w )e'(( k + k e ) r - l i ") + ^ A(2\k,io)e'(kr-^ - k,u' k.u k.w k, , u . Again multiplying by e !( k 2' r a ' 2 ^, integrating over the fundamental length and time scales, and the substituting for A^1) from 2.10, the second-order equation becomes: A ( 2 ) (k 2 ,w 2 ) = -2iP0Jfeocos0HW(k2 - k» ,w 2 ) + (2),' k,o.' k i ,wi (2.11) We could eliminate one of the double summations through the use of the delta functions, but will leave 2.11 as it is for the moment. Now H(2>(k2,w2) is the second order waveheight Fourier coefficient derived by Weber and Barrick [58]: H<2>(k2lu2) = 53 J2 T ^ k ^ s k ^ ^ O H ^ ^ C k ^ O H ' ^ C k ! ^ , ) ^ 1 ^ 1 k.oj k, ,0"! where fh is the hydrodynamic coupling coefficient of Weber and Barrick [58], given by: (2.12) 7 f c(k 1w,ki 1w 1) = - k + h + g(kki - k • ki) gko + ui2 gk2 - co2 (2.13) Consider 2.11 above. If we substitute 2.12 for H^2^ in the first term, then the double summation can be factored out after some re-arrangement. After this substitution and re-arrangement, the double summations are almost the same. In the second term, we are summing over all k, so we can replace k by k + kg, provided we make the following changes: Chapter 2. Acoustic CODAR Theory 15 H ( 1 ) ( k - k s , w ) - f HW(k,w) ka ks—kg ^ = y/kl - P - Ĵfcg - (Jb̂ + ib0 sin _ Jfc2 After this substitution, 2.11 becomes: A ( 2 ) (k 2 ,w 2 ) = -27P0fco cos 6>]T 53 7T(k,w,ki,wi) x k .a1 k i ,u> i H ^ 1 » ( k , ^ ) H ( 1 ) ( k 1 , W ] ) ^ ^ " k k 1 / a : 2 + U J l (2-H) where IT = ~fh-r tja and (2.15) Ta = A* = ^/*g cos2 0 - kt(kx + 2k0 sin 9) - lb*. (2.16) 7a is called the acoustic coupling coefficient. We now have expressions for the first and second order acoustic Fourier coefficients in terms of the first-order waveheight Fourier coefficients and the coupling coefficients. 2.4.3 Expressions for the first and second order acoustic pressure. Using 2.10 and 2.14, we can write equations for the first and second order acoustic pressure in an analogous way to 2.3 for the surface height description, i.e.: P^(v,z:t) = ^ A<1'2>(k2)w2)e,(k'-r-t**-<"'i'1) k 2 )t*-'2 where '1,2' denotes either first or second order terms. For the first order, we have: P ( 1 ) (r,z,*) = -2 JP ofc ocos0 53 H ( 1 ) ( k 2 - k « , w 2 ) e ' ( k 2 r - * ^ - a ' 2 ' ) (2.17) We can write the second order equation in an analogous way to 2.17 as follows: Chapter 2. Acoustic CODAR Theory 16 J° ( 2 )(r>M) = -2iPofcocos0 £ H^ 2 )(k 2 - k S ) w 2 ) e , ( k 3 - r - * « 2 - W a t ) (2.18) 1̂2,(̂2 where HJ/> is actually: H ( r 2 ) ( k 2 - k t f , W 2 ) = ^ £ 7 r ( k , " , k i , w O H f l ) ( k , u ^ (2.19) k,w ki 2.5 The Acoustic Doppler Backscatter Spectrum The waveheight directional spectrum S(k,w) is defined in terms of the waveheight Fourier coefficients in the usual way: (^) 2 (^) 5 (k,u,0 k = k 1 ) o ; = W l (H(k,w)n'(ki,vi))= { (2.20) 0 k ^ k i , OJ ^ ui The acoustic spectrum G(k,w) is defined in a similar way: | (3 i ) - (? i )G(k ) W ) k = k I ) W = W ] (A(k,w)A*(ki,w,)) 0 k ^ ki,w ^ wi (2.21) Note that as L, T get large, we can take: —> dkx or dky ; (jj^J ~* d2k = dkTdky ; and ——> cL;. Then (A(k,w)A*(k,w)) = G(k,w)d2kcL; gives the average acoustic power scattered from the surface at k.w in the band d2kduJ. 2.5.1 Scattering direction and the wavevector k 2 The wavevector k 2 defines the direction of the given plane wave component of the scattered field, projected onto the horizontal plane. The z-component. of the scattered wave is given by 2.7 as: Chapter 2. Acoustic CODAR Theory 17 We can write k 2 in terms of the scattering angles 0S and <ps, using polar co-ordinates with the — z axis as the polar axis (see Figure 2.1 above). Now k 2 is just the horizontal projection of the scattered plane wave, which has wavenumber ko, so |k 2 | = kos'm0s. We can write all the components of the scattered wavenumber k, as: (k2)x = fc0 sin 0, cos (2.22) (k2)y = k0 sin 0, sin ifs k, = ko cos 6f Using these relations, we can write d 2 k 2 in terms of the polar co-ordinates using the Jacobian of the transformation: F i n a l l y : rf2k2 = d'-'k. d9sdtps d6sd<f, = Hky)x d(k2)x de, dip, d£sd(ps d2k2 - k\ cos0s sin OsdO}d<ps (2.23) and we can rewrite the acoustic power spectrum 2.21) as (A(k 2,w 2)A*(k2,w 2)) = G(k2,Lj2)k'o cos9s sin 9sd0sdtpsdu>2. But sin 0sdOsd<f>s = d2Q,,, the increment of solid angle. Therefore the above expression gives the average acoustic power at a spatial frequency k 2 and temporal frequency w2 scattered into the solid angle increment in the direction specified by 0S and <ps. 2.5.2 Reduction to a backscatter spectrum We are now ready to calculate the power received at a known direction from a finite piece of surface at some distance from the receiver. For a surface area A' x Y (both large compared with XQ), which we will call the "scattering patch", at a distance R such that R^> X or Y , looking in the backscatter direction (where 9S — 0 and ips — w), we have the projected area X x Y: A T cos 0 R2 Chapter 2. Acoustic CODAR Theory 18 This gives the total solid angle subtended by the surface area, and the power backscattered from the area is G(k 2 , ^2)^0 c o s 2 0 To convert this expression to a backseattering cross-section per unit surface area, we must divide by the power which would be incident on an area A' x Y at a distance R. from a source level Po: 4TTP 2 ' Dividing by it, we can now write a temporal spectrum of the backscattering cross-section as: <ra(e,v2)=47rk2°'fffG(k2,u2) (2.24) where k 2 = (—kD sin 6,0) — —kg; and kz — k0cos9 2.6 The Acoustic Backscatter Spectrum in terms of the Waveheight Directional Spectrum. Now that we have an expression for the acoustic backscatter Doppler spectrum (Equation 2.24) and expressions for the first and second order acoustic Fourier coefficients in terms of the waveheight Fourier coefficients (Equations 2.10 and 2.14), we can write expressions for the first and second order backscat- tered Doppler spectra in terms of the waveheight directional spectra. The basic relation used in this development is a generalization of 2.24 and 2.21. cr^\0,.2) = 4J^H ( £ j (A (^ ) (k2 , W 2 )A*( 1 . 2 ) (k2 , W 2)) (2.25) where the '1,2' determine whether we are computing the first or second order backscattering cross- section. 2.6.1 First order Substituting 2.10 into 2.25 and using 2.20 we can write the first-order acoustic backscatter Doppler spectrum as a function of the waveheight directional spectrum: XY R2 du>2 ff£°(0.W2) = 16**0 cos 4 0S ( 1 ) ( -2k „ ,w 2 ) . (2.26) Chapter 2. Acoustic CODAR Theory 19 This is Bragg scattering. For backscatter, the first order acoustic return due to roughness is scatter from that part of the surface wavefield with wavenumber 2kg (or wavelength Ao/2sin#) travelling either directly toward or away from the source/receiver. Observed Doppler shift due to first order scatter will be two narrow bands symmetrically positioned about the carrier. The Barrick and Weber [6] definition for the directional waveheight spectrum in deep water is: S(k, w) = \ (S (k)«(w - v ^ * ) + S(-k)6{u + Vg~k)) . (2.27) Using this in 2.26 gives: o^\e,w2) = 8*k40 cos4 9 \s(2ke)6 (w2 + yfigkl) + S{-2ke)6 (w2 - y/2^j where S(k) is the directional waveheight spectrum in deep water. 2.6.2 Second order We begin by substituting 2.14 into 2.25, which results in the following: 1 3 4 (2.28) where the notation has been simplified by (for example): £ = £ ; 3 k3A'3 1T4)* = 7T(k3,W3,k4,Lj4) ; and Hi = # ( 1 ) (k i ,^ ) The (•) term is simplified using an identity for the quadruple product of Gaussian random variables A,B,C and D: (ABCD) = (AB)(CD) + (AC)(BD) + (AD)(BC) Chapter 2. Acoustic CODAR Theory 20 The assumption of Gaussianity is valid at least to first order (see Kinsman [23]), especially in deep water. We use this relation to find: (•) = 2 ( i/( 1)(k, w)^ 1)*(k, w))(//( 1)(k ], W l)^ 1)*(k 1, W l)) = *(T)\¥)2 S { K ' U ) S { K L ' W L ) Using this relation, we can rewrite the equation for the second order Doppler backscatter spectrum: o<?\0,u2) = 32nk4ocos40Y, E |7T (k, W ,k 1 ) W l)|M y 2TTV /2-7T T ^ ( k ^ O S t k ! , ^ ! ) ^ ^ ^ 1 First we eliminate the ^ w by taking rf- —*• du, —• J* and using 2.27, giving: 2 /2?r <r<2)(0,w2) = 16jr*«cos 40 E E E ^ ( k . m v ^ A . k i . w i ) m = ± ! k k.,^ To eliminate the second £^ we again use 2.27. Then, for fixed ki the ^2 will only pick out values satisfying the delta functions, i.e. u\ — ±y/gki, so that we can rewrite the expression once again as a{2)(e,uj2) = 8xk$co&U 53 £ £ | 7 r ( k , m V ^ , k 1 ) m i v ^ ) | 2 r ^ x m,ma=±l k kj V J 5(77?k)5(n71k1)*, ̂ ,> <5a,2v We now have to eliminate the double summation over k and ki, which we do by noting that k 2 — kg = —2kg and using the delta function involving the k values to eliminate the ^Ikj ^ v picking out only those values of ki for which ki = -2kg - k We are left with the single summation over k, which we transform to an integration by letting ( x - ) 2 —>• <f2k and Ylk ~~* Jk- To make the integration symmetric and to take into account the delta function condition on k i , we use a new integration variable s, defined through: Chapter 2. Acoustic CODAR Theory 21 k = - kg + s ki = - kg - s The final result is (replacing u>2 with OJ): o[2){e.u) = Specs' m = ±l m i = ± l • / - c o 5'(mk)S'(miki)6(w — rtiy/gk' — JVI \fgk~l) k.m^/gT•.kum]^ygl^) x (2.29) with the second order Bragg condition: k + ki = -2kg (2.30) implied. The delta function has also been written using different notation. To gain some physical insight into second order backscatter, consider 2.29 and 2.30. First we see that there is a double scatter involved. The surface wavevectors of the scattering wave pair must be such that their vector sum is equivalent to Bragg scattering as in the first order case. The Doppler shift, (u) will depend on the wavenumbers k and k\ and on the direction of travel of those waves. For each wave pair satisfying 2.30 there are 4 possible Doppler shifts, depending on the values of m and mi. We will see later that these 4 different combinations are associated with mutually exclusive regions of the Doppler spectrum. At any particular value of ui, the strength of the backscattered return is proportional to the sum of the products S(mk)S(rr?]ki) for all wave pairs in the scattering patch satisfying both 2.30 and the delta function constraint. These products are also multiplied by the known coupling coefficient The Bragg condition allows us to simplify and rewrite the coupling coefficient JT = fh + ija by noting that kx + (ki)x = — 2ky sin 6 ; and ky+(kl)y = 0. The acoustic coupling coefficient ja (Equation 2.16) simplifies to IIYI 2. Chapter 2. Acoustic CODAR Theory 22 T a = ^fc2Cos20 + k - k j . The hydrodynamic coupling coefficient. 7/, (Equation 2.13) can also be rewritten as (2.31) 7h k + ki- (kk] - k • ki) |"w2 + UJ2B ??U72i ̂ fkkx (2.32) where uiB is the Bragg frequency, given by uiB = gk2 = 2gko sin 9. 2.7 Normalization The expressions are now simplified and non-dimensionalized by normalizing some quantities. First we define normalized, dimensionless wavevectors, frequencies, Doppler spectra, and directional waveheight spectra through the following: 2&o sin 6 rj - — u>B Z(q) = 16^sin 46»5(k) &A = uBo-a (2.33) 2.7.1 First order The first order normalized Doppler backscatter spectrum (Equation 2.28) becomes m=±l (2.34) Chapter 2. Acoustic CODAR Theory 23 2.7.2 Second order First we normalize our integration variable s, and change to integrating over normalized variable q through: We transform the integration to one over q (q = |q|) and <p through d2s = 4kl sin"' 6qdqd<p. Finally, the coupling coefficients (which have the dimension of k) are normalized and made non- dimensional in the same way as the wave-vectors, i.e.: r 7T \ T — 2ko sin 9 and we arrive at the normalized, non-dimensional expression for the second-order Doppler backscatter spectrum 2.29: s kg + k x = 7T + q-2Ao sin 6 2*o sin 6 2 In terms of horizontal polar coordinates, we have: Z(mq)2(miqi )5(»j - m^fq - miyfq~[). (2.35) The second-order Bragg constraints are now written as: q + qi = —x ; and T] - rriy/q + mi y/qj. (2.36) The normalized coupling coefficients are given by: (2.37) Chapter 2. Acoustic CODAR Theory 24 Ih 2*o sin 6 (991 - q • q i ) mniiy/qqt [77 rf + 1' n2 (2.38) Chapter 3 The Directional Waveheight Spectrum from Acoustic Doppler Backscatter 3.1 Beamforming As will be seen in Chapter 4, the receiving section of the acoustic instrument is an equispaced, two dimensional and square 3 by 3 array of hydrophones, each with an omnidirectional receiving response. Given such an array, different antenna patterns can be formed by selecting the appropriate weighted sum of the individual hydrophone outputs. However, for the purposes of wave spectra measurement, we need an antenna pattern which is: as directional as possible; expressible as a truncated Fourier series; and which can be swept in azimuth without excessive variation in directional response. Antennas with beam patterns of the form cos2n((t/> — <fi)/2), where <j> is the azimuthal position of the distant source and %l< is the antenna steering angle, fulfill these requirements. Application of the CODAR technique [3] to the acoustic case requires that we have an antenna which has an azimuthal power response of the form cos\W-4>)/2). 3.1.1 The 4 element sub-array It is possible to form an antenna whose beam pattern varies as cos4((^ — ip)/2) using an omnidirectional receiver and two receivers whose responses vary in azimuth as coscj) and sin<j!>. Multiplying the latter by cos '̂ and sin ip respectively and adding the omnidirectional response gives the following sum (multiplied by 1/2 for convenience): i ( l + cos <j> cos ip + sin <j> sin ip) = (̂1 + cos(ip - <£)) = cos2 ~ ^) (3.39) This expression gives the antenna voltage output as a function of <f> and rp; squaring to find the power response gives the required cos4 pattern. We can mathematically construct these three receivers given a 4-element square array of omnidirec- tional hydrophones arranged as shown in Figure 3.2. Two of these hydrophones, for example HI and H3, form a dipole, shown in Figure 3.3 with the co-ordinate system and various angles defined. 25 Chapter 3. The Directional Waveheight Spectrum from Acoustic Doppler Backscatter 26 Figure 3.2: 4 element hydrophone array. The reference direction is through H4 to H2. We can look at the transmit pattern of the dipole. since it will be the same as the receive pattern. First define some position O in space in terms of spherical co-ordinates (r, 6.9) where 9 is the co-elevation angle. Let the two hydrophones emit pressure waves of the form: P(r,t) where ko po e'(^<-*or) T HE Ao is the acoustic wavenumber. Then if ro d rj RS ro — — sin 9 sin <j> ; and d „ . f 3 ss ro + — sin psin ^. Now if we form Pi - P3 at O, we get: 2 » / V ^ - * ° r ° ) . fk0d . . \ "1 ~ "3 ~ sin I — - sin 0 sin <j> ) . ro \ 2 J Let d = Xo/h, with h >̂ 1; the separation of the hydrophones must be much less than one wavelength at the carrier frequency (this is a fundamental requirement of the array). Then, since sin x ts X for ar <C 1, Chapter 3. The Directional Waveheight Spectrum from Acoustic Doppler Backscatter 27 we have: Figure 3.3: Dipole formed by two hydrophones. 2fP 0e , U' (~* o r o ) T . „ . Pi, = P] - P 3 ~ T sm 9 sin <? r0 h - TT = 22P— sin # sin 0 h where P = ro For sufficiently large h, the response of this hydrophone pair is very close to sinusoidal, with a loss in sensitivity (compared to a single omnidirectional hydrophone) given by 2ns\n6/h. Using a similar argument, we can show the pressure field at O formed from P 2 — P 4 to be: 2iPoe%(-wi~kor°) 7r Pc = Pt — P4 « — sin 9 cos 6 r0 h - 7T = 2zP— sin 9 cos <p. h In this way the receivers with coŝ i and sinfli azimuthal responses are constructed. The omnidirectional response is simply formed by adding the pressure outputs of all the array elements together. For example, adding P] and P 3 together gives: Pi + P 3 ss 2P cos sin 9 sin ^ and for small x, cos(i) ss 1 — x1 j2. Using this approximation we get: Chapter 3. The Directional Waveheight Spectrum from Acoustic Doppler Backscatter 28 , TTsin ( , . . , Pi + P 3 « P 2 - — r — sir/ <p • n\ 2 TT Sill P Similarly for P 2 and P 4 we get: P2 + PA * -P 2̂ - l ^ " 1 1 ^ - J cos" ^ Adding all the pressures together results in an output independent of azimuth, but with a slight loss in sensitivity: E 4 ~ / f TT sin 0 Note that the omnidirectional response is also 7r/2 radians out of phase with the dipole responses. An antenna with a cos4 azimuthal response can be formed by a linear combination of Pa, Pb, and Pc after mathematically correcting the sub-antenna responses for amplitude and phase mismatches: P„ — T-P„ = - Pi oc 1 (omnidirectional); Pt, = —Pb = —(Pi — P3) sin 9sm <j> ; and 27T 27T Pc = A p c = A ( p 2 _ p4) oc s-me cos <j>. ATT ITT Combining these as in 3.39 to form the antenna pattern and then squaring to get the power response F2 of the 4 element array gives: F2(^,<j>,9) = i ( l + sin(?cos(V;-0))2. This expression has the required response for 9 — TT/2 but not for other co-elevation angles. If we multiply the omnidirectional response by sin#, then the power response of the array to a signal originating from co-elevation a and azimuth <j> when the antenna is tuned for co-elevation 9 and azimuth ip is given by F2(i{>,<j>,9,a)= i(sin0 + sinacos(t/)-</>))2 ; (3.40) Chapter 3. The Directional Waveheight Spectrum from Acoustic Doppler Backscatter 29 and for a = 0 F2(i',<t>,0,6) = sin2 0 cos 2 a „„„4 (3-41) The voltage output of a hydrophone is directly proportional to the pressure input. Therefore we replace Pj by Vj, the voltage output from the jth hydrophone. With the hydrophones arranged as in Figure 3.2, we construct the 4-element antenna tuned for co-elevation angle 0 and azimuth ip as a linear combination of the hydrophone voltages: F(Q< 4>) = I (. ,T^-7f/h2 + ^ 4>( Vi - Vs) + £ cos iP(V2 - VA, 2. \4 - TT- sm 0jhl i t 2TT (3.42) Note that the Vj are assumed to be complex quantities, so that we have both magnitude and phase information as a function of time for each hydrophone. 3.2 Data Products - the bn Defined The power output of the 4 element antenna described in Equation 3.41, which I will denote as ̂ ( n , 6, ip), can be defined as: VA(V,Q,IP) sm20 f+* / 1 2* y_, 4>-<$> VA{ri,9,<p)d<j>. (3.43) Here we are assuming that the antenna is tuned for co-elevation 6 and azimuth ip. The antenna output is written as a function of n, the normalized Doppler frequency, so the implicit assumption is that we have transformed from the time domain to the frequency domain (see section 3.2.1 later). The definitions of the normalized narrow beam backscatter Doppler spectra a A (first and second order) are given in 2.34 and 2.35. We can rewrite 3.43 using an exact trigonometric identity: 4' - 4> £ anfn(<f>)fn(4>), (3.44) n = -2 whe cos(n<f>) if n > 0 sin(—n<p) if n < 0 Chapter 3. The Directional Waveheight Spectrum from Acoustic Doppler Backscatter 30 and 0.125 n = ± 2 On = < Substituting 3.44 into 3.43 yields: 0.5 77 = ± 1 0.375 T?. = 0 < 7 A ( 7 7 , M ' ) = S-^~ E anfnW / fnWMvJ,*)**- (3-45) »! — _*> J — TT n = - 2 Now let us approximate the antenna output ^ as a truncated Fourier series: j 2 M ' / . f f ^ ) * ^ E bn(V,0)fnW- (3-46) Z7T ^ — ' n = - 2 Then the 6„ are given by the definition of the Fourier coefficients as: b„(r,.,e)=—[ VA(ri,6,4>)U{<i>)d<p (3.47) where \ 1 if 77 = 0 fr. = < ( 1 if 77 ^ 0 By inspection of 3.45 and 3.46 we can write bn(Vl6) = an sin2 0 f <rA(ri,0s<l>)f„(<j>)d<f>. (3.48) J — TJ But we have expressions for a A (both first and second order) in terms of the directional waveheight spectrum, so that Equation 3.48 provides a link between the bn and the desired information, and 3.47 is used to derive the bn from the measured antenna output. 3.2.1 Writing the bn in terms of cross-spectra. Given time series of digitized complex hydrophone voltages V}, we transform to the frequency domain through the discrete Fourier transform(DFT), so that: Vj{t)D-^ Vjiri)- (3-49) These time-averaged complex spectra(V}(n)) can then be used to form cross-spectra of the form (77) • V£(r}) , which are the actual data products used in either forming antenna outputs or carrying out the Chapter 3. The Directional Waveheight Spectrum from Acoustic Doppler Backscatter 31 process of es t imat ing the direct ional waveheight spect rum. The co- and quad-spectra are defined in the usual way: Jv (Vjiv) • Vk'(v)) =N Cjk(n) + iNQjkin) ; (3.50) where the V indicates that these values are TV-averaged estimates of the corresponding true quantities. T h e Cjk are the co-spectra and the Qjk are the quadrature spectra. Defini t ion 3.50 immediate ly leads to: Qu = 0 Qi, = -Qu (3.51) The antenna voltage t ime series described in 3.42 can be wri t ten in the frequency domain s imply by replacing any Vj by Vj , since F is a linear combinat ion of Vj's. Then the power spec t rum can be wri t ten as F • F". W o r k i n g out this result is a tedious exercise in algebra, resul t ing i n . a series of terms each consist ing of a product of constants (which can be imaginary) and t r igonometr ic terms m u l t i p l y i n g a sum of cross-spectral terms. Since the power spec t rum is by definition real, these sums of cross-spectral terms tu rn out to be either sums of co-spectra ( i f the constant mul t ip ier is real) or quadrature spectra (if the constant mul t ip l ie r is imaginary) . T h e b„(r),6) are defined in 3.47, and the power spect rum of the antenna output is defined by: &A(V,0,1>) = (F(V>6><I>)-F*(V,O,1>))- (3-52) We know from the expressions for the F(w,9,iJ>) that the integrand in Equa t ion 3.47 for the bn(rj,9) w i l l involve products of t r igonometr ic terms like sin(2<?!>) c o s 4 ( © ) , etc. Cons tan t terms and terms involv ing only n (like the cross-spectral terms) w i l l move outside the integral , and we w i l l be left w i t h only s imple t r igonometr ic terms under the definite integral , which can then be solved analy t ica l ly (in fact, many of these terms w i l l integrate to zero). T h e final result is that the data products used i n the est imation of the di rect ional waveheight spect rum, the bn(n,9), w i l l consist of weighted sums of various combinat ions of the t ime averaged cross-spectral estimates. In order to ensure that no mistakes were made, the algebraic reduct ion of the expression for F • F* using 3.42 was carr ied out using a p rogram called ' R E D U C E ' , a 'L ISP ' -based language designed for Chapter 3. The Directional Waveheight Spectrum from Acoustic Doppler Backscatter 32 general algebraic computations. The results of this manipulation can be seen in Appendix A. To simplify the appearance of the results, the following substitutions were made: yi sin 6 . K2 = T T - C O S ^ ; and K3 - ^-sin^ The resulting expression P4 can be used directly to find the antenna output voltage as a function of azimuth by substituting for the K< and putting in the appropriate values for azimuth in the trigono- metric expressions. Using P4 in Equation 3.47 results, after solution of the definite integrals containing trigonometric functions, in an expression for the bn(rj.6) in terms of weighted sums of cross-spectra. These expressions can be seen in Appendix B. 3.3 Relating the bn(n,6) to the Directional Waveheight Spectrum 3.3.1 The inverse problem There has been a great deal of work done on the solution of integral equations of the first kind using numerical methods (see for example a review by Oldenburg [42]). Such equations are of the form: (7(a) = J K(a,b)f(b)db ; (3,53) where estimates of g(a) (the data) are available, and the so-called "kernel" K(a,b) is a known function: usually it is some model of a physical process. The object is to determine the form of f(b). Let's look at the general form of the equation we have to solve. Using 2.35 in 3.48 and also using the delta function to eliminate the integration over q, we get (ignoring the -̂dependence and constant multipliers): bn(ri) R» [+* K(r,,<p)Z(q')Z(q'1)dtp, (3.54) where K(n,ip) includes such things as the coupling coefficient, and q' and q'j are found through the delta function and the Bragg conditions and are uniquely defined for any fixed values of t) and <p: details can be seen in section 3.3.4. The main point is that both q' and q'j are functions of both 77 and <p, so that 3.54 can be written in the form of 3.53 as: g(a) = J K(a,b)f(a,b)db, Chapter 3. The Directional Waveheight Spectrum from Acoustic Doppler Backscatter 33 which is even more difficult to solve than 3.53. One approach to simplifying 3.54 to the form of 3.53 is to make some assumptions about the form of Z(q), i.e. to propose a. simple model for Z(q) and to convert 3.54 to an equation relating the b„ to the model parameters. It will also be necessary to linearize 3.54 somehow. This is the approach taken by Lipa and Barrick [32] in analysing CODAR data. In the following development, the normalized non-dimensional equations for narrow-beam acoustic Doppler backscatter (2.34 and 2.35) will be used, since they are simpler to work with and the results are easily transformed back into dimensional quantities. 3.3.2 Modelling the directional waveheight spectrum We begin by approximating the directional waveheight spectrum for any particular value of normalized wavenumber q by a truncated Fourier series. This is the same approximation used in standard analysis of pitch/roll buoy data (Longuet-Higgins and Cartwright [37]). j = -2 which means that the complete spectrum Z(q) is approximated by: A' 2 2TT 1 2 E djWM), (3.55) ^E E dii9i)fi(4>) i=lj=-2 for q in the range {qx,... ,q.\'}- The fj are defined in 3.44 (replacing n by j). 3.3.3 The first-order 6„ For the first-order b„ we use 2.34 and 3.55 in 3.48 to get: i , 0, an€n7Tsin2 flcot46> bn(m,6) = d„(l)(-»7i)". (3.J>6) Clearly the first-order bn are line spectra defined only for normalized frequencies n = m — ± 1 . It is in theory possible to get a complete directional waveheight spectrum by making a series of measurements at different carrier frequencies and then using this first order relationship to estimate the dn at different Bragg wavelengths. There are several difficulties with this approach. First, in order to measure the direc- tional waveheight spectrum at ocean wavelengths ranging from 1 m to 100 m would require a projector capable of efficient operation oyer the range 10 Hz to 1000 Hz. Such devices do not exist. Secondly, each measurement will require sufficient time to get reliable statistics— 8.5 minutes to get eight 256-point Chapter 3. The Directional Waveheight Spectrum from Acoustic Doppler Backscatter 34 time series at a pulse repetition frequency of 4 Hz. Tb get estimates of d„ for 10 different wavenumbers would require 1.4 hours, not including time spent in reprogramming the instrument between measure- ments. However, the surface statistics are only quasi-stationary, especially in a sea which is not fully developed (Kinsman [23]), so that the resulting directional waveheight spectrum would not in most cases be correct by the time the measurement series was completed. Finally, to get quantitative estimates of the d„ would require absolute calibration of the acoustic system, taking into account projector and receiver frequency transfer functions, path losses, etc. 3.3.4 The second-order bn For the second order 6 n , we use 2.35 and 3.55 in 3.48 to get: , , „ x o„ sin2 9 cot4 9 [' , , f°° , C , 1¥, , , ,, M » 7 , 0 ) = - 5 T d<p qdq d<p\TT\-fn(4>)x 6 (t] - rriy/q - mi y/qT) • (3.57) Simplifying the integral equation. Simplification of this integral equation is mainly due to work by Lipa and Barrick [31]. Equation 3.57 can immediately be simplified by noting that the dependence on <j> can be separated from the rest of the integral. Let us define a quantity Inji as: Injt('fi,<Pi,m,m1) = Next we can eliminate the integration over q by using the delta function for r) and the Bragg conditions seen in 2.36. The delta function is rewritten as a function of q through: 6(T]-my/q-miy/qT) = 6 (r? - h(y, <p)) ; Chapter 3. The Directional Waveheight Spectrum from Acoustic Doppler Backscatter 35 where y = yfq : and q\ = q2 + Iq cos <p + 1 (3.59) which gives h(y,<p) = m2/ + m1(j/4-f2j/2cosv?+ 1) 1 / 4 . (3.60) We can transform the integration over q to an integration over y through qdq = 2y3dy with no effects on the limits of integration. Finally we change the integration over y to one over h(y, ip) through: dy — dy dh dh where: 6y dh (y4 + 2y2cos<p+l)3'4 m(y4 + 2y2 cos ip + l ) 3 / 4 + miy(cosip + y2) Using 3.58, 3.60, and 3.61 in 3.57, we arrive at the following expression for the second-order bn: (3.61) where y is a solution to M>7,0) = an sin2 0 cot4 I Air 2 „ LI J — Tt y3Dh(y,<p)\rT\2 x m , m i = ± l E rfj(?) E rf'(9i)7«j'(^»vi,"j>wi) i= - 2 (3.62) »j-A(y,^) = 0. (3.63) For any fixed values of r] and cp, q and <?] can be found using 3.60 and 3.59. Symmetry considerations. Figure 3.4 shows the relationship of the two normalized wavevectors q and qi in Cartesian r — s space, where the r axis points in the "look" direction of the antenna. Integration of Equation 3.62 is actually carried out over r — s space. The integration is symmetric, in that for fixed n and any pair of wavevectors Chapter 3. The Directional Waveheight Spectrum from Acoustic Doppler Backscatter 36 (-5,0) i*S,0) r_* Figure 3.4: Wavevectors q and q i in r—s space. Curves of constant n are plotted for the inner sidebands. The curves are symmetric about the r axis. q and q i satisfying conditions 2.36 there is another pair of wavevectors q ' and q\ such that q ' = q i and q'j = q which also satisfy conditions 2.36. Therefore we can simplify 3.62 by doubling it and imposing the condition q < q\. Further examination of 3.62 and 2.36 shows that, with the above restriction on the magnitude of q versus that of q i , each of the four terms in the double summation over m and mi contribute in different and mutually exclusive regions of the Doppler spectrum. For example, in the region r? > 1, only the term in 3.62 with m = 1 and mi = 1 contributes to the integral. Linearizing the integral. Equation 3.62, because it involves the product of two sets of unknowns dj and di, is non-linear. However, under certain conditions it can be linearized. Figure 3.4 shows contours of constant T) for the inner Doppler sidebands, and Figure 3.5 shows these curves for the outer Doppler sidebands. Note that when \TJ\ ss 1, q is small, and 91 ss 1. Then di(qi) ss Let us define: r-m+A bn(m,e)= bn(r1,e)dr)=bi;m = ±l. (3.64) The integration over 6„ is necessary because the measurement system cannot return a delta-function exactly. Now using 3.64 in Equation 3.56 we can write: Chapter 3. The Directional Waveheight Spectrum from Acoustic Doppler Backscatter 37 = 46'.(C)f,4.;for m = ± L ( 3 - 6 5 ) a;e|7rsin" 0 cot 0 We can extend the region of applicability of this linearization by making the argument that for qi ^> q, i.e. for small wavelengths, the directional waveheight spectrum Z(qi) is saturated, and thus follows a gj~4 amplitude dependence. Therefore, given d|(l). we can write rf/(?i)w^;for9l«l. (3.66) Finally, using symmetry as described above, and 3.65 and 3.66, we can rewrite Equation 3.62 as: *~ m.,m1=±lJ-* ?1 Yl di(q) Yl ^jj^-Inji('P,<Purn,ml) (3.67) j = -2 / = -2 Q , e ' where the condition q < q\ is implicit. Because of the linearization method used, the solution of the integral equation is restricted to cases of long wave/short wave interactions. The short waves are considered to be wind-driven and fully saturated, i.e. have reached the peak power density possible for their wavelength (see Kinsman [23]), so that we can use the equilibrium spectrum to describe their power spectral density, at least over a small range of wavelengths near the Bragg wavelength. Lipa and Barrick [32] restricted the use of their analysis technique to regions where the normalized Doppler frequency is less than 0.2. Then the minimum value of qi will be 0.8 compared to a maximum value of 0.2 for q. Thus the worst-case ratio of long wavelength to short wavelength will be 4:1. 3.4 Solving the Integral Equation for Second Order bn. 3.4.1 The constant bands approximation To solve Equation 3.67 we must evaluate the integral over <p. The solution involves dividing r — s space up into annular regions concentric about (r,s) = (|,0). Within these bands dj(q) is assumed to be constant, so it can be brought outside the integral. Figure 3.5 shows these bands, together with contours of constant rj, for the outer Doppler sidebands. Note that the first constant-r? contour is entirely within Chapter 3. The Directional Waveheight Spectrum from Acoustic Doppler Backscatter 38 t S Figure 3.5: Bands of constant dj(q) (dashed lines) and constant 77 (solid lines) for the outer Doppler sidebands in (r, s) space. The curves are symmetrical about both the r and s axes. the first wavenumber band, whereas the next contour spans two bands, etc. Let us divide up the q region into I annuli, each of width A,(not necessarily equal). Then we can write 3.67 as : 1 2 Mi?, ") = E E KM*, WAn) (3-68) i=lj=-2 where Uijn(T},6) = ^ j - / d^p ¥jDh(y,<p) \TT\2 Y] 6'^ 7nj/(v?,y>i,m,mi) ; ~ J-n Hi j _ _ 2 °/w and g, - A, < 9 < q{ + A,-. (3.69) Notice that the sum over m,mi has disappeared, since m and m\ are uniquely denned by 77 when we require q < qi- For any value of |?7| there will be two sets of b„(rj), corresponding to positive and negative Doppler shifts. The t/,j„ will be different for these two sets, but the dj will be the same. Thus there are ten pieces of data to be used to estimate 5 parameters. Including both the inner(r7 < 1) and outer(n > 1) Doppler sidebands increases the number of data to 20. If we now extend the system to include P values of TJP, which is the way the data are organized, 3.68 becomes: / 2 bn(riP,6) = E E Ui:n{Vp,0)dj(qi) (3.70) i=lj=-2 Chapter 3. The Directional Waveheight Spectrum from Acoustic Doppler Backscatter 39 For any choice of inner or outer sidebands, we then have a system of linear equations relating 10 x P data points to 5 x 7 parameters to be estimated. 3.4.2 Matrix form of the equation Equation 3.70 can be written as a matrix equation: b = Uw (3.71) where b is a column vector of length 10 x P; w is a column vector of parameters to be solved for, 5x7 in length; and U is a known matrix with 5 x 7 columns and 10 x P rows. The elements of U are the Uijn defined in 3.69. 3.5 Solving the Matrix Equation. Now that our original equation 3.57 for the second order bn has been transformed to to an integral equation of the first kind and then to a matrix equation 3.71 via the constant-bands approximation, we can set about finding the "best" way to solve it. 3.5.1 The straightforward approach. The simplest approach is to treat the task as a simple matrix inversion problem. First we convert 3.71 to a square system: U T b = U T U w Then, if U r U has an inverse, we can solve for w directly: w = ( U T U ) _ 1 U T b (3.72) The usual problem with this approach is that the matrix U T U tends to be ill-conditioned (Golub and VanLoan [18]). Small changes in the data lead to large changes in the solution, which is not good because we expect the data to contain statistical fluctuations about some "true" values, and we want the solution to remain stable while the data are allowed to take values over this statistical range. Generally, such ill-posed matrix problems require more complex methods of solution, often termed regularizaiion. Varah ([55],[56]) presents several of these techniques, and we will use a variation of one of them in the following. Chapter 3. The Directional Waveheight Spectrum from Acoustic Doppler Backscatter 40 3.5.2 Linear least-squares theory. The next "best" method to 3.72 when dealing with data which have known statistics is linear least- squares theory. A good description and theoretical development of this technique can be found in Jenkins and Watts [22], Appendix A4.1. With this method, we model our linear system as: b = Uw + e (3.73) where e is an error vector: the data are now assumed to have some fluctuations about their "true" values. We assume that these errors are jointly Gaussian in their distribution, and that: <e> = 0. Then the covariance of e is given by: Covariance (e) = {̂ (e — (e)) (e — (e))T^ = <-T> Rewriting 3.73 as: b = (b) + e where (b) is the "true" value of the data, then: e = b - (b) and the error covariance matrix V is given by: V=(ee r )=((b-(b)) (b-(b)) T ) (3.74) Therefore the error covariance matrix is equal to the data covariance matrix. The covariance matrix is discussed in more detail in section 3.5.4. Now it can be shown (Jenkins and Watts [22]) that the generalized least squares estimates w are those which minimize the following expression: ( b - U w f v ^ b - U w ) . (3.75) Differentiating 3.75 with respect to w and setting the result equal to zero gives the least squares estimator for w: Chapter 3. The Directional Waveheight Spectrum from Acoustic Doppler Backscatter 41 w = ( U r V _ 1 U ) 1 U T V - 1 b . (3.76) This method also depends on the matrix U T V _ 1 U being non-singular and well-conditioned. A check on the validity of the solution can be made under the assumption that the estimator w is the true value of w, thus satisfying 3.71 exactly. Then we can compute the error vector e using the observed data b and the estimator w in 3.73. The variable p2 defined by: />2 = e T V _ 1 e (3.77) is chi-square distributed with n degrees of freedom (where n is the dimension of b) and standard nor- malization (p2) = n. Then we can ask the question whether p- is close enough to the expected value to be statistically acceptable. By this 1 mean that we choose a confidence limit {e.g. 80%), then compute the maximum acceptable chi-square for the given values of n and confidence level, and see if p1 is less than that value. If so, the solution is statistically acceptable; if not, we throw it out. 3.5.3 Adding constraints to the model. In fact, we know more about the physics of the problem than is specified by the model U . We can improve the estimator w by including this a priori knowledge into our estimation procedure as additional constraints on the solution. This is essentially a regularization technique, discussed in [56], [43], and [52]. The technique I will describe here is a modification of methods first used by Long and Hassleman [36] and later discussed in a more general context by Lanzano [28] and van Schooneveld (in [54]). First, we know that the normalized directional waveheight spectrum Z{q.<f) is a positive quantity for all q and so that an estimator of the "non-positiveness" of the spectrum is: f+ (Z(q,4>)-\Z(q,<p)\fd4. (3.78) J —TT Next, we might have some a priori knowledge about the form of the spectrum: perhaps from preliminary processing of the data. If we denote this a prion guess as Zrj(q,<f>), then an estimator for the departure of the estimated Z(q,<p) from Zo(q,<i>) is: [+\z(q,<f>)-Zo(q,<j>)fd<t). (3.79) Chapter 3. The Directional Waveheight Spectrum from Acoustic Doppler Backscatter 42 Long and Hassleman [36] combine the statistical requirement 3.77 and the a priori: information 3.78 and 3.79 into what they term the "nastiness function" c, given by: where a and /? are free parameters which are chosen to give varying weights to the three constraints to be satisfied by the solution. Differentiating c and setting the result equal to zero will give the best estimator for the waveheight spectrum under the conditions of the statistical requirement, the a priori knowledge, and the chosen values for the weights. Applying the regularization technique. To simplify the development, we will consider the directional waveheight spectrum at a single value of normalized wavenumber q, so that w is a vector of dimension 5 and b is a vector of dimension 10. Then from 3.55 we can write: <; = e T V - 1 e + J — TT f+ [l3(Z(qJ)-\Z(q.4>)\f + Q(Z(q,<f))-Z0(q,<p))2]d<p (3.80) (3.81) where $ = {sin 2<f>, sin <j>, 1, cos <p, cos 26}T. Using 3.81 in 3.78 we get: (Z(q,<p)-\Z(q,<i>)\fd<f>=~ (3.82) Now suppose we define the absolute value as follows: Then we can rewrite 3.82 as: J+\z{qA)-\Z{qA)\fd4 7T ( * T w ) 2 ( l - H{<f>))d4 (3.83) and 3.79 is easily written as: (3.84) where WQ is the a priori quess. Putting 3.83 and 3.84 in 3.80 gives: Chapter 3. The Directional Waveheight Spectrum from Acoustic Doppler Backscatter 43 ? = e T v - ' e +^y * [ / ? ( $ T w ) 2 ( l - J f f ( ^ ) ) + | ( $ T ( w - w 0 ) ) 2 ] # . (3.85) To find the best estimate of w we need to find the total derivative of c with respect to w, given by: 6c = dc dc de <9w de dvf Now from 3.73 we have: (3.86) Differentiating 3.85 with respect to e and w gives: de and dc_ 1 >+* dw = z\ J [P (*Tw) (] - H W ) * + § (*T(w - w0)) *] d</>. Using 3.87. 3.88 and 3.89 in 3.86 gives the desired expression for the total derivative of?: (3.87; (3.88) (3.89) 6c = 2 U T V - 3 e + - ^ y+ [/?(<J>Tw)(l - / / (^) )$+ | (c l ) T (w-wo))$] d<j> 6w. (3.90) Now we set [•] to zero to minimize 6c and rearrange to put all known quantities on the right hand side, giving : / + ir [2/?$Tw(l-^((?)) + a $ T w ] « J > ^ - | - 4 v T 2 U T V - 1 U w = 4 v T 2 U T V - 1 b + •ir Q / _ + ; ( $ T W O ) ^ . (3.9i) which gives a system of linear equations to solve for w (5 equations in 5 unknowns). If we write 3.91 in matrix notation, we have: (P + Q)w = z where z is the known right hand side vector, P is given by: (3.92) Chapter 3. The Directional Waveheight Spectrum from Acoustic Doppler Backscatter 44 Pij= l+ *i*j[20(l-H(<l>)) + a]d<l> ; (3.93) J — 7T and Q is given by: Q = 4 T T 2 U T V - 1 U (3.94) These elements can be evaluated numerically. H(<p) is evaluated using the most recent estimate for w. In practice, the technique is apphed iteratively, as follows: 1. Get a first estimate of w using 3.76. 2. Fill vector wo with a priori information. 3. Choose initial values for a and /?. 4. Compute the right, hand side of 3.91(z). 5. Solve the system 3.92 for a new estimate of w using 3.93. 6. Let: W o l d ~l~ w n t m ^new — ^ to prevent overshoot. 7. Recompute the covariance matrix to reflect the new estimate of the errors. 8. Check for convergence, i.e. P i < f for some epsilon. 9. If not convergent, then go to 4. 10. If convergent then get a new estimate of e, using w „ e u , in 3.73. Compute p2 using 3.77 and compare to desired value. 11. If p2 is too big, reduce a or j3 or both, and go to 4. 12. If p2 is small enough, we have a "good" estimate of w under the applied constraints. Chapter 3. The Directional Waveheight Spectrum from Acoustic Doppler Backscatter 45 3.5.4 Computing the covariance matrix. There are a number of ways the covariance m a t r i x can be computed, depending on how the da ta are defined. Fo l lowing the development of Long[35], let us define a da ta vector x w i th co- and quadrature spectra as the elements of the vector. For a 4-element array of hydrophones there are potent ia l ly 32 elements in x, but from 3.51 we see that there are only 16 independent and non-zero co- and quad-spectra. We can a rb i t ra r i ly define x as: X = ( C n , C12, C13, C14, C22, C23, C24, C33, C34, C44, Ql2, Ql3: Q l 4 , <523, ^24r Q 3 4 ) T • (3.95) Us ing 3.42 and 3.52 in 3.48 we can write: where b is the vector of b„ values, and B is a known m a t r i x of coefficients. Let x be one real izat ion of the data . T h e n b = Bx is one real izat ion of the vector of 6 n values. Now let x = x + 6k where x is the "true" value of the data. T h e n we can wri te b as a Taylor series expansion to first order in 6k: (3.96) b = b(x) + B6x (3.97) since Now to lowest order in 6k we can say: (b) = b(x) and therefore Chapter 3. The Directional Waveheight Spectrum from Acoustic Doppler Backscatter 46 b-(b) = B*x. (3.98) Using 3.98 in 3.74 allows us to write the covariance matrix as: V = B((5x6xT)BT. (3.99) The statistics of <5x are well known. Each of the elements of the matrix {^X^XT) is given by one of the following relations (Jenkins and Watts [22]): COl)\Cpq , Grs } — (C'prCqg + QprQqs ~¥ Cps CqT 4" Qps Qqr) V COll{Cpq , Qrs } — (CprQqS — QprCqS — CpsQqr -\- QpsCqT) V COv{Qpq,Qrs} = {CprCqS + QprQqs ~ CpsCqr - QpsQqr)"'1 (3.100) where v is the number of degrees of freedom of the cross-spectral estimates. Each of the quantities on the right-hand sides of 3.100 is an element of the vector x, the "true" data. In principle, we can find the elements of x given the "true" b, through: x = ( B T B ) _ 1 B T b (3.101) Therefore, assuming that a solution w to the system 3.71 is the correct solution, 3.71 yields the "true" value of b, which can be used to compute the covariance matrix V through 3.101, 3.100, and 3.99. In practice we cannot solve 3.101: B T B is mathematically singular. This fact is not surprising, since examination of 3.101 shows that it describes a linear system of 5 equations in 16 unknowns! This problem has come about because in following the CODAR processing scheme we have synthesised a three-element CODAR antenna out of a 4 element monopole array. Fortunately our synthesized CODAR antennas are so nearly perfect (see section 3.1.1 above) we can use derivations by Lipa and Barrick[33] of the covariance matrix elements in terms of the bn (see their equations (16), (18) and (19)). In section 3.5.2 the N-averaged measured bn values are assumed to be the "true" values and are used to compute the covariance matrix. This matrix is also used as a first guess for V in the procedure of section 3.5.3, but thereafter we make the assumptions that the vector b resulting from using the latest estimate w in 3.71 is the "true" value used to compute V. Chapter 4 Modelling of Acoustic Doppler Backscatter In order to get some idea of what measured data should be expected to look like, given the models developed in the preceding chapters, I wrote some software to produce the expected backscatter power spectra, and the expected form of the B„ functions. 4.1 Omnidirectional spectra using the Phillips wave spectrum First I will develop a model for the power spectrum of acoustic Doppler backscatter seen by an omnidi- rectional receiver. We will use the normalized equations for acoustic Doppler spectra given in 2.34 and 2.35. The full directional waveheight spectrum is given by S(k)G(8), where S(k) is the non-directional waveheight. spectrum, and G(B) is the directional factor. The usual normalizations apply, that is: /•CO / s(k)dk = < c 2 > Jo / Gifi)dO = 1 * J — 7T There is an implicit assumption here that the two parts of the spectrum are independent of each other, which is in fact not usually the case. However, the assumption is adequate for the simple modelling we are going to do here. For the non-directional part of the wave spectrum Z(q) we will use the Phillips spectrum (see LeBlond and Mysak [30]), with a simple cos4 model of the directional factor: 2(q) = 1 3 ^ \ 2 J V - V C (4.102) I 0 q < qc 2wj0*0 sin 0 Here g is gravitational acceleration, U J O is windspeed, and <f>* is the principal direction of propagation of the sea. 4.102 is valid for rn. = +1, and for m = —1, we have <j> —+ <p + VT to indicate the reverse direction 47 Chapter 4. Modelling of Acoustic Doppler Backscatter 48 of propagation. Now we can rewrite 2.35 using this definition of the wave spectrum: ^ 2 ) ( M ) mq{5 + mi q° 5 (cos <p + q)\ *x / J . I (m 1 — 1 J , * COS / . (m-l)-K ,» 4 I 9 + <P COS (4.103) In order to get this formulation, 1 have used simplifying relations given in 3.59, 3.60, 3.61, and 3.63. Now 4.103 is a narrow-beam spectrum, i.e the spectrum which would be recorded by an instrument looking at a small surface patch distant from the receiver. The line between the scattering patch and the receiver defines an azimuth angle of zero. Thus the spectrum seen by an omnidirectional receiver will be the integral of 4.103 over all azimuth angles: *i 2 ) (M) cot4 e ? _ J > / > ' ^ ( T £ ) ' ( » > - 2.5 m,mi =± 1 ,0.5 q]5 + miq°-5(cos0 + q)\ 1 x / I i J i (m-l)ir j»\ / „ / i A , ( n n-l)* 4 / W + <P + 2 V \ 4 y + VI + 2 COS 2 (4.104) The tilde indicates that this is an omnidirectional spectrum. The cos4 terms can be expanded as shown in 3.44, and we can take the integral with respect to ip inside, since the only place that ip appears is in the cos4 terms. Then the cos4 part of 4.104 can be written as: 2 2 AQI = E 5Z aJa'fj(4>* ~ ; = -2 J=-2 The fj and /; are orthogonal functions, so j = / for a nonzero result, and 4:105 becomes „ (m — 1 )JT l = -2 where 9*1 32 / = 0 ( m i — l)vT <t>l) \ / = ± 1 — / = ± 2 64 ' (4.106) Chapter 4. Modelling of Acoustic Doppler Backscatter 49 It turns out that 4.106 has an exact solution: Ao: = — ((mmi cos(^ -</>) + A)2 - 7.5) . (4.107] Finally 4.104 can be rewritten using Atn as: ' m,m] = ±l n -1 |m« 1 1 - 5 - |-m 1g 0 - 5 (cos^ + ?)| A01. (4.108) Some simplifications can be made due to symmetry considerations, which will reduce the computa- tional difficulty significantly. First, due to the Bragg constraint 2.36 on q and qj we can impose the condition that q < q\ and multiply the integral by 2. Second, we note that all functions of angle in the integrand are symmetric with respect to so we only need integrate from 0 to TT and multiply by a further factor of 2. Finally, since the integral is symmetric about the frequency origin, we only need to solve it for positive rj. 4.1.1 Solving the integral Values of q and q\ are found for fixed <p using the Newton-Raphson method. For q < q\ the value of rj forces the values given to m and m\. The starting point of the integral is determined by the value of » and the requirement for an initial guess for the Newton-Raphson technique. For the inner sidebands(n < 1), we start the integration at <f> = TT, and use: Vq= \ (\/2 - n2 + mr,) to get the (exact) initial value for q. For the outer sidebands (rj > 1), we start the integration at <p = 0, and use: to get the (exact) initial value for q. A further simplification exists for the outer sidebands. From Figure 3.5 we see that the curves of constant n are not always closed in the half-space defined by q < qi, so the integration may not proceed all the way to TT. For n > \/2 we can define the upper limit of integration by: p m a x = TT — arccos Chapter 4. Modelling of Acoustic Doppler Backscatter 50 As the integration proceeds, the previous value of q is used as the initial guess in the Newton-Raphson procedure. For each <f> in the integration, values for q and q\ are found. The value of (̂ i is found using: ql-qj-l COS <?!>!=: . 4>i must be in the third quadrant because of the condition q < q\ and the Bragg constraint 2.36. 4.1.2 The first-order spectrum. All the above computations are for the second-order spectrum. From 2.34, the first order Doppler spectrum is given by: m=±l Then substituting 4.102 for Z, integrating over all azimuth directions, and taking n = 1 (thus m — — 1) gives: $H0, l)=°-£- cot4 e £ cos4 d*. (4.109) But the value of the integral is 3VT/4, independent of <f>*, so the first order spectrum is a delta function at i] = 1, with height given by: o-̂ V, 1) = 0.00125 cot4 6. (4.110) 4.1.3 Smearing. Because the spectral resolution of the measurement system is finite, the spectrum seen by any receiver will not be the same as the model. In order to account for this in our model spectra, we need to "smear" the model results out in frequency space in some way which, reflects the properties of the measurement system and the geometry of the measurements. The first effect to take into account is the finite pulse length. Following the notation in Burdic [10], the pulse shape is given by: F{t) = A0 cos(w0i)rect (4111) which is a gated sinusoid of frequency w, where the gate length is given by r. In frequency space, this signal and its power spectrum look like: F(u) = AoT (cos(woO) ® ? ^rect E(u) = FF* Chapter 4. Modelling of Acoustic Doppler Backscatter 51 A2 j j s i n 2 ((u - u0)r) (4.112) (w - w o ) 2 The operator '<g>' is the convolution operator, and 'T' is the Fourier transform operator. By looking at the -20dB point of the ratio: we can get some idea of the amount of "smearing" of the measured spectra due to the finite pulse width. Typical pulse lengths used in my measurements were about 17 ms. The -20dB half-width of the energy ratio for this pulse width is at u> — WQ = 0.59, or about 0.1 Hz. There is at least one other source of frequency smearing, again due to the finite pulse width. The outgoing pulse interacts with the surface over an annulus of finite width, giving rise to a variation in the co-elevation angle, and thus a variation in the Bragg frequency. For example, at 45 degrees co-elevation angle, with the receivers at 50 m, at 400 Hz carrier, the variation in the Bragg frequency is 0.05 Hz. Of course most of the energy is coming from interactions near the center of the annulus, since the whole pulse interacts there, whereas only very small portions of the pulse interact at the extreme edges of the annulus, so this estimate of "smearing" is over-estimated. In any case, this source of "smearing" is at most 1/4 of that due to the spectral content of the pulse itself. To account for the effect of both sources of "smearing", we need to convolve our model spectra with a smearing function in the frequency domain, as: I chose to use a Gaussian in frequency space as a smearing function, with the choice of the parameter p given by the desired smearing halfwidth of 0.1 Hz. Therefore: and thus p = 0.2928. For spectra in normalized frequency, we simply replace w by u/ug in 4.114, which leads to: and pB — 0.2928/wj5. In applying the smearing function, the energy content of the original spectrum is conserved by applying an appropriate normalization. E(u>o) X'(OJ) - G(u>)® X(u) (4.113) (4.114) Chapter 4. Modelling of Acoustic Doppler Backscatter 52 -10 -70 1 1 1 1 0 1 2 3 Normalized Frequency Figure 4.6: One-sided power spectra as seen by omnidirectional receivers. The dashed hue is for wind- speed of 10 m/s. and the solid line is for 5 m/s. 4.1.4 Results for the omnidirectional receiver (Phillips spectrum). In Figures 4.6 and 4.7 the spectra do not disappear in the region between the first order and second order peaks: they actually continue down below the boundaries of the plot, but the plotting routines used do not draw the continuations of the lines to the boundary. Figure 4.6 shows one side of the power spectrum for two cases: the solid line is for a windspeed of 5 m/s; and the dashed line is for a windspeed of 10 m/s. Both these spectra are for co-elevation angles of 30 degrees. First note the Bragg peak at a normalized frequency of 1. These peaks are the same for both cases because the waves causing Bragg reflections are fully saturated for both windspeeds. The second order energy is in the two sidebands. Note how the two curves coincide at high and low frequencies: this behaviour is a consequence of the seas being fully saturated at the wavelengths which produce second-order interactions at these frequencies. As windspeed increases, the sideband energy increases, and energy is seen closer to the Bragg frequency, since longer waves (and thus smaller q and 9i values) are now present. Eventually as windspeed increases, the second order energy will swamp out the first order peak. Chapter 4. Modelling of Acoustic Doppler Backscatter 53 -10 -70 0 1 2 3 Normalized Frequency Figure 4.7: One-sided power spectra as seen by omnidirectional receivers. Both are calculated using a windspeed of 5 m/s: the dashed line shows the spectrum for a co-elevation angle of 60°; the solid line is the spectrum using a co-elevation angle of 30°. Figure 4.7 shows the effect of changing the co-elevation angle. Both spectra are for a windspeed of 5 m/s. The sohd line is the same spectrum as we saw in Figure 4.6; a co-elevation angle of 30 degrees. The frequency axis is different, because now the frequencies are normalized with respect to the Bragg frequency for a 60 degree co-elevation. We see that increasing the co-elevation angle increases the frequency of the first order peak, and reduces the overall spectral power. In addition, the "peakiness" of the outer sideband is reduced, and the separation between first and second order energy is somewhat increased. 4.2 Model bn using the Phillips Spectrum. We can use 4.102 in 2.35 to compute model b„ as defined in 3.48. Carrying through the derivation of this integral is very similar to the exercise for the omnidirectional spectra seen above, so I will simply write down the result: 6„(77,0) = a„ VT cos ;2 e cot 2 e Chapter 4. Modelling of Acoustic Doppler Backscatter 54 \mq{ 5 + mi9°- 5(cos<£ + q)\ *x 2 2 E E fMlJiWVnii^fa.mM) (4.115) j = -2 J=-2 where the function Injj is defined in 3.58. For the first order portion of the bn, we put 4.102 in 2.34 to compute 3.48. This results in: .02a„ cos2 e cot2 e ^ , , , (m+I)VT x r j = -2 But the /j are orthogonal, so that: , , „, 9 6» J U , ,„ + l)jr, [* r , J X , , , , , , bn(m,e) = !-— J2 a i l M + 2 ) fr,(4>)fj(<f>)d4>- (4-116) J — TT where e„ is given in 3.47. Thus we can write: , , ., .027ra2 cos2 6 cot2 9 . , (7n+l)vT. , . , , „ ^ K{m,9) = 2—-; c n / „ ( 0 * + i ^ - ) . (4.117) 4.2.1 Some examples of computed bn functions. Figures 4.8 and 4.9 show some computed bn functions for the same windspeed and carrier frequency, but different principal wave propagation directions. The y-axis of these plots is logarithmic, with the magnitudes of all plots normalized by the height of the first order peak of the bo spectrum. Thus a change of 1 in the magnitude is equal to 10 dB. Note thai, the b0 function is just the omnidirectional spectrum reflected about the frequency origin and multiplied by a constant 27raosin26'. The shape and sign of the second order parts of the bn change in a complicated way depending on the wave propagation direction. However, the sign of the first order and immediately adjacent second order parts of the spectrum change in a simple way with wave propagation direction, as can be seen from an examination of 4.117. Considering the negative frequency bn for simplicity, we see that it is proportional to /„(<?>*). Thus a cursory examination of plotted bn functions can immediately reveal a rough estimate of the principal wave propagation direction. In addition, we note that the 60,62, and 6_2 are symmetrical about the frequency origin, while the 61 and 6_i are anti-symmetrical. 4.3 Omnidirectional spectra using the JONSWAP non-directional model and a variable directional factor. In order to more closely approximate the spectra which we might see in an actual geographical location, I developed a model using a more sophisticated non-directional wave spectrum and a variable directional Chapter 4. Modelling of Acoustic Doppler Backscatter USPO: 10.0 fl/S. PHI: 137.0 SNR: 0.0 l/OB. THETfi: 45.0 FREO: 338.6 HZ. 55 N • 2 Figure 4.8: Computed bn functions using the Phillips spectrum. Windspeed is 10 m/s, co-elevation angle is 45°, and principal wave direction is 137°. The horizontal axis is the normalized Doppler frequency. Chapter 4. Modelling of Acoustic Doppler Backscatter 56 •USPO: 10.0 n/S. PHI: 210.0 SNR: 0.0 1/OB. THETR: 45.0 FREQ: 398.6 HZ. 1 Figure 4.9: Computed bn functions using the Phillips spectrum. Windspeed is 10 m/s, co-elevation angle is 45°, and principal wave direction is 210°. The horizontal axis is the normalized Doppler frequency. Chapter 4. Modelling of Acoustic Doppler Backscatter 57 factor. Much of the earlier data in this project was collected in Saanich Inlet, so the model developed reflects the geography of that region. 4.3.1 Description of the directional spectrum used in this model The JONSWAP spectrum. The J O N S W A P spectral form is a fairly recent formulation for the non-directional wave spectrum under conditions of limited fetch and limited depth of water over which the wind blows. The J O N S W A P spectrum is given by (see LeBlond and Mysak[30]): ag 5 / u> (2VT) 4,..5 exp - -4 \ui (u - um): where v = exp I —rt— 2(7- and cr = (4.118) 0~a i f W < Ulm [ crj i f w > um The spectrum is thus a function of the parameters a. um. A, cra, and o-j. a and uim depend on the fetch, through the following relations: o = 0.076(2TT)5 ( '10 -0.22 5.887 w h e r e T , = tanh 10.52 1 9 " 0.37 tanh 0.077 (4) 0.25 \ V tanh (0.52 ( £ ) 0.375 (4.119) and / is the fetch, « 1 0 is the windspeed, and z is the mean depth in the upwind direction. The fetch and depth-dependent expression for T, is one of the Bretschneider [8] equations. The J O N S W A P spectrum described above is a function of frequency. We convert i i to a wavenumber spectrum through the normalization: /•OO /"OO / J(u)dw = / J{k)kdk =< C2 > • Jo Jo Thus Chapter 4. Modelling of Acoustic Doppler Backscatter 58 We use the deep water dispersion relation (w2 = c gk) to solve this relation: some error is introduced because the mean depths in some directions are relatively shallow. However, since we do not expect to generate very long waves in this enclosed site, the errors will be very small. Then we find that: Geographical input to the model. The mean depth and fetch as a function of direction (in degrees true) around the array site were determined at angular increments of 5° starting at due North. These data are shown in Figures 4.10 (fetch) and 4.11 (depth). The directions of greatest fetch are roughly Northwest and Southwest. Strong winds are likely to blow from the NW during the summer months, and strong SW winds are sometimes associated with the passage of low pressure fronts over the area. The mean depths are greater in these directions as well, although the depth shallows when moving toward the NW from the SW. The best conditions for maximum waveheight would probably be strong, consistent winds from the SW. 4.3.2 The directional factor. A directional spectrum can be developed by letting the contribution to the spectrum for a particular direction (p be given by the JONSWAP spectrum using the mean depth for the direction (<j) — ir) and an effective fetch Fe which is given by where F{4>) is the fetch in the direction (<f> — ir) and /3 is the downwind direction. There is usually an additional directional factor applicable to the directional spectrum. This factor accounts for the concentration of spectral power along the wind direction. Typical spreading factors are described by functions like: (4.120) Chapter 4. Modelling of Acoustic Doppler Backscatter 59 -8 -6 -4 - -2 0 2 Km East of array site Figure 4.10: Fetch (km) around instrument site in Saanich Inlet. Chapter 4. Modelling of Acoustic Doppler Backscatter 60 Or -20 I " 4 0 JS -60 - CL -S -80 g -loo 2 -120 -140 50 100 150 200 250 Direction (deg. true) 300 350 Figure 4.11: Mean depth (m) around instrument site in Saanich Inlet. The higher the value of V , the more concentrated is the spectral power in directions close to the wind direction. In our case, we use similar functions for G{6), but the normalization is expressed as: G{<f>)d<p = 1 ; since we don't allow spectral power in upwind directions. 4.3.3 The total directional spectrum. Using the above results we can describe the total directional spectrum as: 5( M ) = J(k,<p)G(o) with the usual normalization: /•OO ft J J J(k,<j>)G(<j>)kd<pdk=<C2 > This type of spectrum has been used by Seymour[49] with some success to predict the nondirectional spectrum seen in fetch-restricted areas. 4.3.4 Computations of typical spectra, and comparisons to the truncated Fourier 6 e r i e s representation. Figure 4.12 shows a family of JONSWAP spectra for the geometry described above. All the spectra were computed assuming a windspeed of 10 m/s, but with differing wind directions. The differing fetch and Chapter 4. Modelling of Acoustic Doppler Backscatter 61 Windspeed 10 m/s Min. wave: 2 m Max. wave: 80 m 102 I- K \ I s I N IO1 - . 10° - u o CL « i o - i . y Q. CO IO"2 - IO' 3 - JQ-4 I ' ' ' • ' ' I I I I I I I 1 I I 1 10_1 10° Wavenumber Figure 4.12: Family of JONSWAP spectra for different wind directions. Dashed line is Phillips spectrum for the same windspeed. mean depth cause differing peak wavenumbers and differing spectral power in the small wavenumber region, but at the high wavenumber range the spectra all come together, as would be expected. For comparison, the Phillips spectrum for a fully developed sea is plotted in the same figure with a dotted line. The maximum wavelength and minimum wavelength are 80 m and 2 m respectively. Note that the wavenumber of the peak spectral power of the JONSWAP spectum is roughly 0.3, corresponding to a wavelength of about 21 m, in contrast to the cutoff wavenumber of 0.098 (corresponding to a wavelength of 64 m) predicted by the Phillips spectrum. Note also that the JONSWAP spectrum predicts more energy at higher wavenumbers than does the Phillips spectrum. This effect has been observed: the "Phillips constant" decreases with increasing fetch[19]. A "typical" directional spectrum was computed, using a 10 m/s windspeed with the downwind direc- tion being 150° true, the approximate direction of greatest fetch. At the same time, the corresponding truncated Fourier series (TFS) representation of the spectrum, given in 3.55 was calculated. The results Chapter 4. Modelling of Acoustic Doppler Backscatter 62 Vintf 10 m/t • I M ««| tewtnd. K- .03 - W i n 10 m/t 9 I M ««| 4*<ra»ind. K« M 0 M 100 ISO tOO C M 300 440 DirecUo* (d<| Imi ) Vu»4:I0 m/s • I M *«4 «m«tad. X - A t 60 100 I M tOO *W 300 S M D4r«cUoa Lmi) I lot I Figure 4.13: Directional JONSWAP spectra (solid line) together with TFS representation(dashed line) for k values 0.02 through 0.05 are plotted as a function of direction for several values of normalized k in Figures 4.13, 4.14, and 4.15. There are several things to note. First, the form of the synthetic directional spectrum is highly "peaky" for small wavenumbers, whereas it becomes quite broad for larger wavenumbers. There is only a narrow region for which the fetch is appreciable; thus all the long waves will come from that narrow band of directions. Secondly, the TFS does a very poor job of representing the directional spectrum in the long wave, "peaky" part of the spectrum, but does approximate the spectrum fairly well in the short wave, saturated region which has a broad directional distribution. The TFS is good for estimating the princi- pal direction of the sea, and numerical results indicate that the estimate of RMS waveheight using this representation is also good. It is not clear from this test how well the TFS spectral form might do when the directional spectrum is bimodal. Chapter 4. Modelling of Acoustic Doppler Backscatter 63 Figure 4.14: Directional JONSWAP spectra (solid line) together with TFS representation (dashed line) for A; values 0.06 through 0.20 Chapter 4. Modelling of Acoustic Doppler Backscatter Wind: 10 m/s @ 150 deg downwind. K=.30 j I I I I I i i 0 50 100 150 200 250 300 350 Direction (deg true) Figure 4.15: Directional JONSWAP spectra (solid line) together with TFS representation (dashed li for k = 0.30 Chapter 4. Modelling of Acoustic Doppler Backscatter 65 J0«««* 10 m/t Hi t n o M 0.1 Hz t V n n , J O N I W 10 » / • I K 0»—»« 01 « t s -w U t - J Figure 4.16: Power spectra as seen by omnidirectional receivers. JONSWAP directional model used, with 0.1 Hz smearing width. 4.3.5 O m n i d i r e c t i o n a l a c o u s t i c d o p p l e r s p e c t r a d e v e l o p e d u s i n g t h e d i r e c t i o n a l J O N - S W A P m o d e l Using the model for the normalized omnidirectional acoustic Doppler spectrum developed in Section 4.1, but using the directional JONSWAP spectrum should give us a more accurate picture of the results we should expect to find in Saanich Inlet. The derivation and computation of this spectrum is very similar to, although more complex and computationally intensive than, the model using the Phillips spectrum, so I will not go through the development. Omnidirectional spectra were calculated for several different conditions. Unless otherwise specified, the spectra described here were computed for a carrier frequency of 400 Hz at a co-elevation angle of 45 degrees, with a windspeed of 10 m/s, and have.been "smeared" as described in 4.1.3. Doppler spectra for 315° downwind (the most restricted fetch and broadest directional distribution) and for 150° downwind (greatest fetch and most "peaky" directional distribution) are shown in Figure 4.16. Note that the spectral level for the inner sidebands is actually higher for the restricted fetch case, whereas the outer sideband peak level is higher for the 150° case. The separation between the first order and second orders is decreased in the long-fetch case because of the greater energy available at small wavenumbers. Figure 4.17 shows the same spectra as described above, except that the width of the smearing function has been doubled. As expected, considerable detail in the spectra is lost, and the first and second order overlay each other to a greater extent. Chapter 4. Modelling of Acoustic Doppler Backscatter 66 JOS*** 10 « . / • >is Oo.—tia o.: ta i J O N S W A P 10 •>/> I S O D M - M 0.1 Hi io-' 1 0 " Figure 4.17: Power spectra as seen by omnidirectional receivers. JONSWAP directional model used, with 0.2 Hz smearing width. 0.15 Hi l o n o r I K K . 10 J . I Figure 4.18: Power spectra as seen by omnidirectional receivers. JONSWAP directional model used, with 0.15 Hz smearing width, and noise at 10 dB SNR added. Figure 4.18 is the same as Figure 4.16 except that random normally distributed noise with a signal to noise ratio of 10 dB has been added to the model. Note how the addition of such noise adds spurious detail to the spectrum and also creates false asymmetries, especially in the first order peaks. In Figures 4.19 and 4.20 we see the effect of narrowing the directional width of the wave spectrum by changing the directional factor from the cos2 proportionality seen in Figure 4.16 to a cos4 proportionality (Figure 4.19) and a cos8 proportionality (Figure 4.20). The main effects are: the width of the second order portion of the spectrum is reduced, particularly the outer sidebands; and the peak second order spectral level is increased somewhat for both inner and outer sidebands. First order spectral levels are unaffected. Chapter 4. Modelling of Acoustic Doppler Backscatter 67 J O N S * * * 10 m / l U S O U . H . I M 0.1 Hj V ~ C C « . " < J0MSKA» 10 « / l ISO O o — ~ J 0.1 " J C o . " ' Figure 4.19: Power spectra as seen by omnidirectional receivers. JONSWAP directional model used, with cos4-proportional directional factor. JOWSWtf 10 ->/• }IS 0u.i.,ii.a 0.1 " J %mm» C i " l « W S « A » 10 ISO O——«< 0.1 «» C " « Figure 4.20: Power spectra as seen by omnidirectional receivers. JONSWAP directional model used, with cos8-proportional directional factor. Chapter 4. Modelling of Acoustic Doppler Backscatter 68 4.4 Application of the JONSWAP/directional wave spectrum model to the computation of model fcn functions. Model bn functions were computed using methods similar to those outlined 4.2, using the JONSWAP/directional wave spectrum. Figures 4.21 and 4.22 show model bn functions for wind directions in two quadrants at the same windspeed (10 m/s). The wind directions are downwind directions: thus Figure 4.21, with winds blowing toward 160°, has the longest fetch, whereas Figure 4.22, with winds blowing toward 210°, has a short, restricted fetch. As before, all the spectral levels within each plot are normalized to the bo first order peak level, so that absolute spectral levels cannot be compared from plot to plot. The difference in the shape of the functions is mainly that the separation between the first and second order peaks is much reduced for the long-fetch cases, indicating the increased spectral power at short wavenum- bers. The resulting functions are quite different from those computed using the Phillips model for the nondirectional wave spectrum. For example, compare Figure 4.21 with Figure 4.8. These have winds at the same speed and roughly the same quadrant, but the spectral levels of the JONSWAP-based bn functions are much lower than those based on the Phillips spectrum, so that the second-order portion of the spectrum is well below the first-order part, and the separation between the first and second orders is much greater. The explanation is most likely seen in Figure 4.12, where we see that, although generally the JONSWAP spectral levels are higher, there is less energy at the low wavenumber end, and the low wavenumber cutoff is more realistic. In addition, the directional factor used in these spectra was the cos2 type, less intensely peaked around the principal direction. 4.5 Tests of the inversion methods using model input Inversion methods developed in Chapter 3 were tested using model input to evaluate their performance. 1 used bn functions computed using the Phillips spectrum as the simplest possible input to the tests. Tests of "smeared" data were carried out on data using different windspeeds, and with and without noise. The use of the direct inversion was compared to the inversion using the linear least-squares method, and finally to regularization and iterative improvement. 4.5.1 5 m/s windspeed A test data set was constructed using model bn functions for 5 m/s windspeed at a principal direction of 35°. The co-elevation angle was 45°. In addition, the data were "smeared" with a 0.1 Hz Gaussian. Chapter 4. Modelling of Acoustic Doppler Backscatter 69 Figure 4.21: Model bn functions using the JONSWAP/directional spectrum. Winds to 160° (long fetch). The horizontal axis is the normalized Doppler frequency. Chapter 4. Modelling ol" Acoustic Doppler Backscatter 70 • USPO: 10.0 n/S. PHI: ??S.O SNR: 0.0 l/DB. THE IB: <S.O FRtO: 400.0 HZ. N . 2 Figure 4.22: Model 6„ functions using the JONSWAP/directional spectrum. Winds to 225° (short fetch). The horizontal axis is the normalized Doppler frequency. Chapter 4. Modelling of Acoustic. Doppler Backscatter 71 First I computed the expected value of the dj(qj) (see equation 3.55): Ii d-i do di d2 .096 +0 000 +0.000 +0.000 +0.000 +0.000 .132 +0 000 +0.000 +0.000 +0.000 +0.000 .144 +0 000 +0.000 +0.000 +0.000 +0.000 .158 +0 000 +0.000 +0.000 +0.000 +0.000 .173 +1 761 +4.299 +5.621 +6.139 +0.641 .189 +1 227 +2.997 +3.919 +4.280 +0.447 .207 +0 853 +2.083 +2.723 +2.974 +0.310 .227 +0 590 + 1.440 + 1.883 +2.057 +0.215 .250 +0 404 +0.987 +1.290 + 1.409 +0.147 .275 +0 273 +0.666 +0.870 +0.951 +0.099 .306 +0 179 +0.437 +0.572 +0.624 +0.065 .344 +0 112 +0.273 +0.357 +0.390 +0.041 Table 4.1: Expected dj(qj) values for 5 m/s windspeed and principal direction of 35°. The results of the inversion process using the straightforward method (see Section 3.5.1) are as follows: d-2 d-i do di do 0 096 0.198 0.373 0.363 0.533 0.072 0 132 0.383 0.840 0.993 1.200 0.139 0 144 0.605 1.378 1.708 1.969 0.220 0 158 0.805 1.868 2.357 2.668 0.293 0 173 0.907 2.130 2.713 3.042 0.330 0 189 0.877 2.083 2.668 2.976 0.319 0 207 0.747 1.792 2.304 2.559 0.272 0 227 0.576 1.396 1.799 1.993 0.210 0 250 0.414 1.011 1.305 1.444 0.151 0 275 0.280 0.687 0.888 0.982 0.102 0 306 0.173 0.426 0.551 0.610 0.063 0 344 0.089 0.221 0.284 0.316 0.032 Table 4.2: Results of straighforward inversion of b„ models for 5 m/s windspeed and principal direction of 35°. Visual inspection of these results shows that they are reasonably representative of the input spectrum, especially in the region where actual wave energy exists (q < qc)- As we would expect, energy has "leaked" into the low wavenumber region because of the smearing of the theoretically calculated fc„ functions, but the peak in wave energy in the do values is close to the right place, although attenuated in value. We can compare the resulting estimates of the directional factor s and the estimates of the peak direction using the computed dj in Table 4.2 above. Lipa and Barrick [3] developed some equations Chapter 4. Modelling of Acoustic Doppler Backscatter 72 <t>* 4>l s'i +35.00 +35.00 + 17.37 +5.84 +35.00 +35.00 +5.62 +4.52 +34.99 +35.00 +4.75 +4.29 +35.00 +35.00 +4.47 +4.20 +35.00 +35.00 +4.34 +4.15 +34.99 +35.00 +4.26 +4.11 +35.00 +35.00 +4.21 +4.08 +35.01 +35.00 +4.18 +4.05 +35.00 +35.00 +4.16 +4.03 +35.00 +35.00 +4.15 +4.01 +35.00 +35.00 +4.17 +4.00 +35.00 +35.00 +4.21 +3.99 Table 4.3: Estimates of the principal direction and the factor s from results of inversion of model data. relating these quantities (see their equations (21) through (26)). I found some errors in their equations: the correct versions are given below. — — , C _ i do do d± . c i = d_zl do ' do P^yjci + Cl! ; P2 = xjc'i + C. fC-i\ 2/xi then <p\ — arctan —— ) ; s\ = 2 2 Ci J ' 2-in 1 t (C-2 — arctan ( — 2̂ 2 V O a n d s 2 = 2 + 3 , 2 + y ^ + 28/,2+_4 ^ The computed values of s and principal direction using these formulas are given in Table 4.3. Clearly the inversion process is returning direction information (which depends on relative amplitude) more faithfully than absolute amphtude information. Addition of noise to the model bn functions changes the situation considerably. 1 took the same functions used in the tests above, and added Gaussian distributed random noise so that the signal to noise level was +20 dB. The resulting functions were then re-submitted to the inversion process. The results can be seen in Table 4.4. An inspection of these results shows clearly that the addition of noise to the 6 n functions causes the straightforward method of inversion to fail. The estimates of principal direction and spreading factor s Chapter 4. Modelling of Acoustic Doppler Backscatter 73 Ii d-2 d-i do di d2 0.096 -.2936 0.7276 0.2161 0.6155 0.4034 0.132 0.5355 0.4359 1.511 0.9135 -.4449 0.144 -.3596 2.485 -.03841 3.893 -.2423 0.158 0.1514 2.454 3.367 0.7576 2.100 0.173 -4.736 9.603 -2.992 8.735 1.342 0.189 1.453 0.2603 3.626 1.846 -.5344 0.207 3.298 -2.969 8.076 -3.838 0.0405 0.227 1.730 0.5467 4.993 -4.036 4.167 0.250 -2.058 6.895 -3.598 6.432 1.077 0.275 -1.289 2.193 -.2236 2.323 -.0708 0.306 -.3985 0.8677 -.1114 1.536 0.2441 0.344 0.6131 -.5214 0.9956 -.6485 0.3010 Table 4.4: Results of straighforward inversion of 6n models for 5 m/s windspeed and principal direction of 35°. Noise at SNL of 20 dB added. 4>l <f>2 •s2 +49.77 -18.02 -3.66 -56.76 +25.51 +64.86 + 1.01 +4.88 +32.55 -61.99 -2.03 -6.14 +72.84 +2.06 + 1.23 +6.22 +47.71 -37.09 -3.71 +40.02 +8.03 +55.10 +0.69 '+4.64 -142.28 +44.65 +0.86 +4.51 + 172.29 + 11.27 +1.38 +9.30 +46.99 -31.19 -8.44 +6.41 +43.35 -46.57 -2.33 -8.86 +29.46 -29.26 -2.29 -12.02 -141.20 +31.93 + 1.44 +6.79 Table 4.5: Estimates of the principal direction and the factor s from results of inversion of model data. Noise at SNL of 20 dB added. are given in Table 4.5. Both the magnitude information and the directional information have been lost. In a way this result is not surprising, due to the "garbage in/garbage out" syndrome. Unfortunately, when the iterative improvement method is applied to this model input, it fails to converge, and the situation is not hopeful. This topic will be discussed further when we apply the inversion techniques to real data. Chapter 5 Instrument Development It should be stated at the outset that much of the work involved in this investigation concerned instrument development, modification and testing. This chapter is only a brief summary of that work. Details of instrument design and operation can be found in Appendix C. 5.1 ULSAS I The first version of the upward-looking sonar system (ULSAS I) was designed specifically as an in- strument to test the acoustic CODAR method for estimating the DWS of surface gravity waves. This instrument was designed and built by Arctic Sciences Ltd. of Sidney, B.C., and is described in a report by Lemon et. al. [14]. ULSAS I was first deployed in Saanich Inlet, B.C. in 1985, and was recovered, modified, and redeployed several times in the period 1985 through mid-1988. Figure 5.23 shows a side view of the deployment of ULSAS I; it was connected to a laboratory in the Institute of Ocean Science in Sidney, B.C. by an underwater armoured cable which carried power and data. ULSAS I sat on the bottom in about 50 m of water on a sand/silt bottom which was generally flat, and exposed to reasonably long fetches in the North- and Southwesterly directions. The instrument suffered from several major flaws. First, there was considerable analog preprocessing of the received acoustic signals before digitization and storage, so that it was difficult to determine the causes of perceived problems with the digital data. Second, the armoured cable through which ULSAS I was connected to shore made deployment difficult and risky, and in fact nearly resulted in the total loss of the instrument during one deployment. Because of the cable, the choice of locations in which ULSAS I could be deployed was extremely limited. 5.2 ULSAS II To eliminate some of the problems with ULSAS I, it was decided to make some major modifications, and it was at this point that 1 took over the instrument development. Some of the changes which were 74 Chapter 5. Instrument Development 75 implemented in the new instrument ( U L S A S II) inc luded: 1. T h e analog preprocessing of received da ta was l imi ted to amplif icat ion and broad-band filtering only. 2. D a t a was recorded using commercia l digi ta l audio equipment, a l lowing for d igi ta l rather than analog signal processing. 3. A completely redesigned mul t ichannel synchronous system for recording acoustic Doppler signals was implemented. 4. A remote communica t ion facility using packet switched radio hardware and microcontroller-based software interface was conceived and bu i l t , implement ing fail-safe defaults. 5. The power sys tem was redesigned for full bat tery power a l lowing up to 4 weeks deployment. 5.2.1 Overview o f system operation. U L S A S operates at a carrier frequency (/o), typ ica l ly in the range 400 - 1200 H z , and the information content of the signal backscattered from the ocean surface is in the modu la t ion of the carrier frequency, i.e. its Doppler shifts and ampl i tude modulat ions , which are caused by reflections f rom the mov ing sea surface. The acoustic projector produces short pulses which are t ransmi t ted omnidirect ional ly . Each pulse consists of an integral number of cycles of a s inusoidal signal at the carrier frequency (/o)- Successive pulses are phase-coherent, that is, the t ime between them is an integral mul t ip le of the carrier per iod ( l / / o ) . The sequence of a pulse of carrier frequency followed by a "dead t ime" is called a ping. T h e dead t ime between pulses is required to give the pulse time to t ravel out from the source, scatter from the surface, and be received back at the source loca t ion . T h u s the t ime between successive pings is l imi ted by the speed of sound and the depth of the source. For shallow depths there must be an addi t iona l delay between pings to allow the reverberation level to decrease well below the expected signal level . T h e m a x i m u m ping rate (or pulse repeti t ion frequency — P R F ) used in this work is approximately 4 H z . T h e backscattered, ampl i tude modula ted and Doppler shifted pulses are received at an array of 9 hydrophones, amplif ied, and digi t ized so as to preserve the phase of the returned signal at each hydrophone locat ion. T h e method used is delayed sampl ing of the signals. If we take the accepted Chapter 5. Instrument Development 76 definition of a time-varying signal x(t) as product of a carrier term and a complex envelope function (see for example Burdic [10]): x(t) = ?K [X(t)el'2*'0i] = A'/(I)COS(2IT/OO + XQ(i)sin(27r/o<) where X(t) = Xr(t) - iXQ(t) then it can easily be shown that if we take consecutive samples separated by a time equal to 1/4 of the carrier period, we get a good estimate of the complex envelope function: *(0l,=.o - * * W U 0 + ^ = A ' / ( O U 0 - 1 XQ(i)\t=to+j_ * X(t0) where io = n//o> " integer, and also provided that the envelope function does not change significantly over the time l/(4/o)- The time of the digitization relative to the time the pulse was transmitted will determine the location on the surface from which the backscattered signal originates. Because of the finite pulse width, the area on the surface which produces backscatter at a particular digitization time is an annulus of variable width. Many digitizations (one per ping) corresponding to a particular surface location produce a time series of complex voltage with a digitization rate equal to the PRF. There will be 9 time series, one for each hydrophone. These time series are as long as possible, consistent with the assumptions of weak stationarity used in spectral analysis of surface waves — in practice the time series are limited in length to 15 - 20 minutes, which gives a total number of around 4000 points in the time series. Analysis of these time series using methods described in previous chapters gives information about the surface target strength and (possibly) the directional spectrum of surface gravity waves., The bulk of the design of hardware and software for ULSAS II had to do with the generation of the coherent output pulses, and the data acquisition of signals from multiple channels at the appropriate times, and with the development of a system of control by which the user could set up data collections according to the kinds of data products desired. In addition to the acoustic system, there was also a monitoring system, measuring such variables as tilt of the hydrophone array, compass orientation, depth, and the state of the power supplies. Sensors to detect flooding were also built into the instrument. Chapter 5. Instrument Development 5.3 Measurements with ULSAS I and II Detai ls of the design and operat ion of U L S A S II can be found in A p p e n d i x C . T h e first version of U L S A S II was a b o t t o m situated design, s imi lar to U L S A S I. A descript ion of the measurements made w i t h bo t tom deployments of both U L S A S 1 and U L S A S II is given i n Chap te r 6 of this thesis. Because of problems w i t h backscatter f rom the bo t tom, i t was decided that U L S A S II would be modified once again for deployments in the deep ocean. M a n y changes had to be made to the ins t rument to allow this k i n d of deployment, and most of the details of the physical configuration of U L S A S II to be found in A p p e n d i x C refer to the deep-water, surface-tethered version of U L S A S 11. The deep water version of U L S A S 11 was taken along as an unplanned experiment on the Sur- face Waves Processes P r o g r a m ( S W A P P ) cruise of the Canad ian Scientific Ship "Par izeau" in Febru- a r y / M a r c h of 1990. T h i s was essentially a "ship of oppor tuni ty" cruise for the deep water version of U L S A S II, and extensive tests and measurements were made dur ing that cruise. Some of these measure- ments, a l though of a pre l iminary nature, are given in Chapte r 7. Chapter 5. Instrument Development 78 NOT TO SCALE Figure 5.23: Side view (not to scale) of ULSAS I deployment in Saanich Inlet. Chapter 6 Saanich Inlet Measurements. Many sets of measurements were made in Saanich Inlet over a period of roughly three years, while we sought to better the quality of the output and understand more fully the source of apparent problems. In this chapter, 1 will show some data from one particular set of measurements. The data were collected on January 8, 1986 at 09:40 PST. There were strong Southwest winds blowing at the time, and the sea surface was rough with many whitecaps. Conditions were probably just about as rough as they can get in Saanich Inlet under normal wintertime conditions. These data were collected at a carrier frequency of 415.55 Hz, with a pulse repetition rate of almost exactly 4 Hz, and a pulse length of 16.8 ms. The post-multiplication filter was set to a cutoff frequency of 2 Hz, so these data are particularly free of high frequency energy aliased into the Doppler region. 6.1 Typical power spectrum. Complex voltages from all nine hydrophones were collected for 10000 pings. These data were broken up into 37 half-overlapped 512 point subseries. Each 512 point subseries was multiplied by a window function (Blackman-Harris -74 dB window: see [20]) before transformation to the frequency domain. The resulting spectra were then used to compute a cross-spectral matrix as in 3.50. Using the notation of that equation, the estimate of the power spectrum for the ith transducer is the cross-spectral product NCH(I]). The estimated power spectrum for hydrophone 5 is plotted in Figure 6.24. Spectra for the other 8 hydrophones are essentially the same. The first thing to note is that the spectrum is not symmetric about the zero frequency axis. Theo- retical analysis and modelling results do not predict an asymmetrical spectrum for an omnidirectional receiver, so this result is surprising. It is not atypical: many of the spectra collected using this version of the instrument were asymmetrical. The second thing to note is that the shape of the spectrum does not resemble model spectra very closely, especially if the Bragg region indicates the location of first order energy. Comparing Figure 6.24 to, for example, Figure 4.16 reveals very little symmetry unless we 79 Chapter 6. Saanich Inlet Measurements. 80 Figure 6.24: Average of 37 windowed, 512 point power spectra for hydrophone 5. Horizontal axis is Doppler frequency (Hz). Expected Bragg frequency regions are between the two sets of dotted lines. Chapter 6. Saanich Inlet Measurements. 81 assume that the position of the first order peaks is in error. If we take the highest peak on the negative side of the spectrum to be the first order peak, then one could argue that there is a corresponding first order peak on the positive side. Then, with the exception of the abnormally high second order energy on the positive Doppler side, we could say that the measured spectrum resembles a model spectrum which has been "smeared" considerably more than we did in the modelling in Chapter 4. How can we account for this shift in the first order peak? The assumed first order peaks are in the 0.744 Hz to 0.595 Hz region, corresponding to co-elevation angles between 25° and 42°. But we know from the experimental conditions that the surface area sampled was bounded by co-elevation angles from 37° to 54°, which corresponds to Bragg frequencies between 0.724 Hz and 0.838 Hz for this carrier frequency. Determination of the co-elevation angles requires the geometry of the deployment (i.e. the instrument depth and the tidal height), the timing of the receiver on/off signals, and the local sound speed. Given these inputs, the co-elevation angles are uniquely determined. Variation of the sound speed is a negligible effect, and in order for the observed shift in first order frequency to be due to uncertainty in the instrument depth, the instrument would have to be about 28 m deeper than we think it is, which is clearly impossible. We are left with only four possibilities: the theory with respect to first order interactions at the surface is wrong; the instrument timing setup information does not reflect the actual measurement conditions; the analog front end of the instrument is somehow distorting the incoming signal to cause a "squeezing-in" in the frequency domain; or the peaks that we think are the first order peaks are not them. The first possibility can be thrown out, because the theory has been amply verified experimentally both in the acoustic case (see for example Roderick and Cron [47] and Koenigs [24]) and with RADAR (Broche [9]). As for the second possibility, the hardware and software was extensively tested to verify correct operation. However, the instrument is not accessible during data acquisition, so we cannot totally rule out some timing variations during actual operation. The third possibility is linked to the second, since a "squeezing" operation in the frequency domain is equivalent to a "stretching" operation in the time domain, which again implies that what we think is happening in the time domain perhaps is not so. We are left with the fourth possiblity, which we can investigate by an analysis of the data. Chapter 6. Saanich Inlet Measurements. 82 6.2 Computed 6n functions. Using the cross-spectral matrix resulting from the processing of the frequency domain data, we can compute the bn functions as outlined in 3.2.1. Figure 6.25 shows the results of this computation for one of the four element subsets of the hydrophones. bn computed for other subsets were essentially the same. These plots are similar to the model plots: the spectral levels have all been normalized to the energy of the first order bo. These plots are all in normalized frequency, and I have used the assumed Bragg frequency rather than the expected Bragg frequency in the normalization. Except for the absence of clear separation between first and second orders, and the addition of a lot of noise, these functions resemble the modelled bn fairly well: see for example Figure 4.21. As discussed in Chapter 4, a cursory examination of these functions should reveal some information about the principal wave propagation direction. Comparison of equation 4.117 with the figure shows that <f>' is in the second quadrant, and 2<f>* is in the fourth quadrant, so that cA* is somewhere between 135 and 180° relative to the reference direction for the hydrophone subset in question. Now the reference direction, which is a line from hydrophone 9 to hydrophone 1, was 107° magnetic, or roughly 130° true, so the reference direction for the hydrophone subset is 220° (see Figure 3.2). Azimuth angles increase in an anti-clockwise sense, so the principal wave direction must be between 40° and 85°, which is consistent with winds blowing from the southwest. 6.3 Tests of inversion methods on real acoustic data. The 6n data described above can be used as input to the inversion routines. I tested these routines on the data using the expected value of the Bragg frequency for normalization, and the results were reasonably good, even though there is no clear separation between the first and second order regions. I also tried the inversion methods using the assumed value of the Bragg frequency from visual interpretation of the power spectra. Results of these tests are shown in the following sections. 6.3.1 The straightforward method The first test was run on the data using a co-elevation angle of 45°, and the solution is shown in Table 6.6. The results look reasonably well-behaved. I ran a regression to the equation log(d0(g,)) = mlog(c/,) + b (6.122) Chapter 6. Saanich Inlet Measurements. File «5 Uindowed 37 by 512 pts THETR: 32.0 FREQ: 415.6 HZ. Figure 6.25: Computed 6n functions for one four-element subset of hydrophones. The horizontal the normalized Doppler frequency TJ. Chapter 6. Saanich Inlet Measurements. 84 9i d-2 d-i d0 di d-> 0.0442 0.5573 -.04567 0.2194 -0.5043 -.08749 0.0567 0.3274 -.03932 0.09606 -.02048 -.002262 0.0632 0.2037 -.01942 0.05644 -.003668 0.05154 0.0702 0.2139 -.01160 0.06169 -.009450 0.04818 0.0777 0.1564 -.01669 0.05315 -.01348 0.02464 0.0858 0.09761 -.01121 0.02723 -.006382 0.02687 0.0946 0.06474 -.01014 0.01938 -.0002516 0.02122 0.1040 0.04650 -.009276 0.02629 -.002408 0.001723 0.1141 0.04610 -.003143 0.02244 -.0005226 0.0006578 0.1250 0.02503 -.003870 0.008342 -.002232 -.005397 0.1327 0.01609 -.003389 0.003516 0.00002529 -.009953 0.1495 0.01277 -.0004725 0.001918 0.001238 -.003912 Table 6.6: Results of straighforward inversion of real bn functions. A co-elevation angle of 45° was used. using the last 8 values of the do(<7t), since the first values correspond to wavelengths which are probably not realistic for Saanich Inlet. The parameter rn, which should be around 4, was estimated to be 4.6 with an r2 of 0.85. Estimates of the principal direction and spreading factor s were also made, and those results can be seen in Table 6.7. The estimates of the spreading parameter S\ and s2 are very different: the si estimates indicate very weak directionality, whereas the estimates of s2 indicate a strongly peaked directional spectrum. The reference direction for the array used in these calculations was 220° true, then for Southwesterly winds, we would expect a principal direction somewhere in the range 155° to 200° (or —205° to —160°). The estimated values for <j>\ are mostly in the third quadrant, although with a great deal of variation. These estimates correspond more closely to Southerly rather than Southwesterly winds. However, the 0% values are tightly grouped together, with a mean of 46°. On examination, it appears that this result is an artifact of the data. The values of d-i are much too large; they should always be less than the do values. Recalling the equation for 02 from 4.121, we see that if d-2 S> d2, 4>2 will always be around 45°. I re-ran the inversion using a co-elevation of 31°, which was the co-elevation angle estimated from the apparent position of the Bragg peaks in the power spectra. These results can be seen in Table 6.8, with the estimated principal direction and s factors in Table 6.9. These data are actually worse than those using the expected co-elevation angle: note especially that the magnitudes of the d-2, d - i , d+i and d +2 values far exceed the magnitude of the do parameters, which is clearly not possible, and gives rise to the negative signs of the estimated spreading parameters. Thus, although the appearance of the power spectrum seems to indicate Bragg peaks in a location different from where we expect to find them, Chapter 6. Saanich Inlet Measurements. 85 tt tt *] Sn -137.84 +49.46 +0.37 -32.92 -117.51 +45.20 +0.60 -16.16 -100.70 +37.90 +0.42 -14.06 -129.17 +38.65 +0.28 -15.08 -128.93 +40.52 +0.51 -21.20 -119.65 +37.30 +0.62 -14.08 -91.42 +35.93 +0.71 -15.35 -104.55 +43.94 +0.45 +64.50 -99.44 +44.59 +0.15 -298.16 -119.97 +51.08 +0.73 -19.80 -89.57 +60.87 +1.86 -9.37 -20.89 +53.52 +1.06 -7.78 Table 6.7: Estimates of the principal direction and the factor s from the results in Table 6.6. "i d-i do di d2 0.0541 -147.4 -8.494 5.797 -4.201 24.97 0.0711 -105.9 -4.067 2.009 -2.305 19.34 0.0800 -71.75 -2.674 0.9037 -1.719 14.30 0.0898 -43.10 -1.715 0.5257 -1.199 10.14 0.1004 -23.73 -1.252 0.1652 -.8291 6.323 0.1120 -11.83 -.9871 0.4113 -.6922 4.252 0.1247 -7.741 -.3151 0.8633 -.4389 3.357 0.1386 -1.993 -.2101 0.4071 -.3266 2.101 0.1538 -1.120 -.1516 0.2006 -.1998 1.372 0.1705 0.1820 -.1035 0.1481 -.1439 0.5473 0.1892 0.01266 -.07147 -.003759 -.07017 0.3160 0.2100 0.1949 -.03649 -.02286 -.03922 0.2275 Table 6.8: Results of straighforward inversion of real bn functions. A co-elevation angle of 31° was used. analysis of the data shows that these cannot be the true Bragg peaks. 6.3.2 Using iterative improvement and the "nastiness" function I tested the effect of iterative improvement on the data, using a co-elevation of 45° in the process, since it has been demonstrated to be the correct value to use. The results can be seen in Table 6.10, with the estimated principal directions and s factors in Table 6.11. Values of 2 and 200 were used for a and 0 respectively, meaning that the positivity of the spectrum was highly emphasized over the a priori guess at the directional factor. Inspection of the derived solution shows that iterative improvement has cleaned up the estimated spectrum a bit: for example, note in Table 6.11 that the estimated spreading factors are now uniformly positive, indicating that the relative magnitudes of the d±2 to the do have Chapter 6. Saanich Inlet Measurements. 86 tt tt *i *2 -116.32 -40.19 +8.95 -4.90 -119.54 -39.83 -14.24 -4.44 -122.74 -39.36 -4.64 -4.29 -124.96 -38.38 -4.02 -4.28 -123.51 -37.54 -2.56 -4.16 -125.04 -35.12 -6.30 -4.76 -144.32 -33.28 +0.91 -6.51 -147.25 -21.74 + 1.82 -7.68 -142.81 -19.61 +3.34 -6.83 -144.27 +9.20 +2.98 -13.20 -134.47 + 1.15 -2.16 -4.28 -137.07 +20.29 -13.65 -5.82 Table 6.9: Estimates of the principal direction and the factor s from the results in Table 6.8. <1< d-2 d-i do di d2 0.0442 0.01310 -.1338 0.2438 -.1987 -.06388 0.0567 0.03351 -.05145 0.1278 -.1249 0.06693 0.0632 -.01546 -.03609 0.1113 -.05091 0.06663 0.0702 0.05624 -.02604 0.08732 -.06346 0.05687 0.0777 0.07195 -.03120 0.07894 -.05556 0.02905 0.0858 0.04090 -.02122 0.04123 -.04081 0.03237 0.0946 0.01071 -.008749 0.02767 -.02258 0.02093 0.1040 0.01315 -.01462 0.03685 -.01035 -.008142 0.1141 0.01536 -.01029 0.03100 -.009422 0.0008622 0.1250 -.008005 -.008066 0.02065 -.009899 -.005887 0.1367 -.004434 -.003754 0.01058 -.004669 -.006663 0.1495 -.004822 0.001183 0.006263 -.002331 0.00334 Table 6.10: Results of inversion of real 6„ functions using iterative improvement. A co-elevation angle of 45° was used. been corrected. A fit of the data to 6.122 gives a value of -3.26 for m, with r 2 = 0.86. <j>\ values are now more tightly grouped together, with a mean value of about 210° (actually —150°), which is closer to the expected direction than the unimproved estimates. The <f>*2 values are too variable to be of much use. 6.4 Beamforming in the frequency domain I wanted to test the operation of ULSAS as a superdirective array, using the beamforming equations developed in section 3.2.1. The expression P4 can be used directly to find the antenna output voltage as a function of azimuth by substituting for the K,- and putting in the appropriate values for azimuth in the trigonometric expressions. Chapter 6. Saanich Inlet Measurements. 87 tt tt *2 -146.04 +84.21 +1.93 +3.58 -157.61 +13.30 +2.24 +5:88 -144.67 -6.53 +0.78 +6.13 -157.69 +22.34 +1.29 +9.47 -150.68 +34.01 +1.35 + 10.47 -152.53 +25.82 +2.52 + 16.60 -158.82 +13.55 +1.56 +8.60 -125.30 +60.88 +0.64 +4.59 -132.48 +43.39 +0.58 +5.15 -140.83 -63.17 +0.90 +5.04 -141.20 -73.18 +0.79 +7.51 +153.09 -27.65 +0.53 +9.77 Table 6.11: Estimates of the principal direction and the factor s from the results in Table 6.10. Spectral data were first 'balanced', so that the energy contained in a narrow band around the esti- mated Bragg frequency was the same for all hydrophones. Then the antenna outputs for 4 element arrays were computed for several frequency bands both above and below the Bragg frequency, for both positive and negative Doppler shifts. There are six possible sets of 4-element subarrays using all nine hydrophones in the array, so that six separate but not necessarily independent, estimates of the directionality of surface returns can be made. Figures 6.26 through 6.29 show the resulting antenna returns. In general, the results are good. However, there is considerable variation in the patterns from sub-array to sub-array. Results from negative Doppler appear to be slightly better than those from positive Doppler, and results from the outer sidebands appear to be better than those from the inner sidebands. Thus the 'best' set of results, in terms of the similarity of the patterns, is seen for the frequency range -0.95 Hz to -0.80 Hz, and the worst results are those from the frequency range from 0.55 Hz to 0.75 Hz. For example, in this case sub-arrays 4 and 6 are considerably different from sub-arrays 1,2 and 5, and sub-array 3 has an almost circular pattern with little directionality. In addition to the plots, the direction of maximum and minimum returned power was also computed for all six sub-arrays for each frequency region. The results are summarized statistically in Table 6.12 below. These results are self-consistent, i.e. the direction of the maximum for positive Doppler bands is the direction of the minimum for negative Doppler bands, and vice-versa. In addition, the difference between the maximum and minimum directions is roughly 180°, as would be expected for waves propagating in Chapter 6. Saanich Inlet Measurements. LEGEND O- SETl O- SET Z 6 - SET 3 + » SET 4 X - SET 5 o m SET 6 Figure 6.26: 4-element subarray directional response for the band -.95 Hz to -.80 Hz I.Ok Figure 6.27: 4-element subarray directional response for the band -.75 Hz to -.55 Hz Chapter 6. Saanich Inlet Measurements. 89 Figure 6.29: 4-element subarray directional response for the band .80 Hz to .95 Hz. Chapter 6. Saanich Inlet Measurements. 90 Frequency(Hz) 0mi„±s.6. @max $min .55 to .75 320 ± 27 137 ± 2 7 183 .80 to .95 340 ± 13 157 ± 2 0 183 -.75 to -.55 165 ± 2 2 332 ± 2 1 167 -.95 to-.80 166 ± 13 333 ± 17 167 Table 6.12: Summarized results of 4-element beamforming. 6max is the mean direction from which maximum power was returned, and 0min is the mean direction for minimum power. one principal direction. The reference direction of the array for this data set was 220° true, which results in a maximum power return direction of 249° true, corresponding well to the wind direction, reported at the time to be Southwest. We would expect the maximum power return for positive Doppler to be in the upwind direction, since positive Doppler is due to advancing waves. Conversely, the maximum power return for negative Doppler would be in the downwind direction, from retreating waves. 6.5 Discussion of Saanich Inlet results. The results of the analyses shown above are tantalizing: they show that it is possible to get some information about the wave field from the acoustic backscatter measurements. Iterative improvement seems to improve the analysis somewhat as well. However, although many data sets were collected and analysed over a long period, it was not possible to definitively prove or disprove the usefulness of the techniques developed here, for two main reasons. First of all, Saanich Inlet is a protected inlet in which sea conditions are never worse than severe chop. Thus it is not really a good "test bed" for a technique which is to be used in the open ocean. Secondly, the configuration of ULSAS used for making all the Saanich Inlet measurements was one in which the "front end" of the instrument was complex and included a great deal of analog processing before digitization of the data. Consequently it was hard to understand what the actual acoustic returns seen by the hydrophone were like. A major problem with the digital data was that, for any particular range gate, there was more often than not a large DC offset to the data which was consistent ping-to-ping. Consequently it was not possible to use anywhere near the maximum available gain without saturating the signal. The useful signal was then a small fluctuation of a large DC offset. Often the variance of the signal was less than 4% of the dynamic range of the ADC system. A great deal of time was spent in analysing the data and thinking about the problem to try and understand the source of the large DC offset. It seemed that the most likely cause was bottom (or sub-bottom) reflections, but there was no way that this hypothesis could be tested with the ULSAS Chapter 6. Saanich Inlet Measurements. 91 configuration used in Saanich Inlet. 6.5.1 Subsequent deployments and instrument development. In order to address the first of the problems outlined above, it was decided to deploy ULSAS in the Strait of Juan de Fuca, near the Sheringham Point lightstation, where there would be improved chances of getting some larger waves and swell. However, problems during deployment nearly resulted in the loss of the entire system, and no measurements were made at that location. As a result of this near-disaster, it was decided to make significant changes to the instrument, so that it would be self-contained, and so that "front-end" processing would be minimized. The result was the instrument described in Chapter 5, except that it was not designed for mid-water deployment. This configuration of ULSAS was deployed in Georgia Strait, just North of Departure Bay. Measurements made there revealed the source of the large DC offsets seen in earlier data. Figure 6.30 shows the measured returns for 10 pings for one hydrophone, from ranges of 56 m to 102 m, so it includes surface returns from near vertical incidence out to a grazing angle of 30°. The returns look like a slowly varying component which is consistent ping-to-ping, with a more rapidly fluctuating component superimposed on it. The most likely explanation is that the slowly varying, consistent part of the return is due to reflections off layers in the sub-bottom — since the frequency used is so low, there will be significant penetration of the bottom — and the more rapidly fluctuating component is due to backscatter from the sea surface combined with ambient noise. Before beginning instrument design, an analysis of expected reverberation levels from both surface scatter and bottom scatter, using data from the literature, was carried out. This analysis indicated that backscatter from the bottom could be neglected; however, no thought was given to possible backscatter from layers in the sub-bottom. The geometry permits specular reflections in that case, which can give rise to large returns, even in comparison to the direct surface return. The implication of this interpretation is that the dynamic range of the ADC system again will be under-used, since the constant DC offsets will saturate if they are amplified too much. Assuming that this interpretation is correct, the only solution is to deploy in midwater with a surface-tethered system, in deep water. This is the configuration used for the SWAPP measurements. Figure 6.30: Backscatter from near-vertical incidence to 30° grazing angle, for ten consecutive pings. Hydrophone 6, Georgia Strait deployment. Chapter 7 SWAPP Measurements The SWAPP cruise was really a "ship of opportunity" for the deep-water version of ULSAS II. Because of the short time available, modifications made to ULSAS II were not tested before it went on board CSS "Parizeau". Despite a number of problems, ULSAS II performed well: a total of 9 deployment/recovery sequences, under a variety of conditions, were carried out without incident. Unfortunately, a serious prob- lem with the acoustic projector precluded making any useful measurements of the Doppler backscatter spectrum; however, useful data were collected, and estimates of the surface scattering strength parameter S, and of the ambient noise spectrum were made. Details of these and other measurements are found in this chapter. Calibration information and a "post-mortem" of the difficulties encountered during SWAPP can be found in Appendix D. Before being able to examine data from ULSAS II and compute parameters of interest, considerable work had to be put into the development of data processing hardware and software, primarily because of the amount of data stored by ULSAS II and the rate at which that data is read off the video tape. Details of that work can be found in Appendix E. 7.1 Selection of data for processing Since the purpose of SWAPP was to investigate sea-surface processes during storm conditions, the weather during the SWAPP cruise was disappointing. Windspeed and direction and other environmental information during SWAPP are shown in Figure 7.31. The highest windspeeds seen were around 18 m/s, and these were not sustained for any length of time. In fact, conditions were unusually calm for the time of year and location. The files 1 have analysed come mostly from deployment number 7, the weather for which is shown in Figure 7.32. There were consistent 8-10 m/s winds during the first half of this deployment, after which the windspeed dropped rapidly to 2-4 m/s until the end of the deployment. I also analysed an ambient sound record made early in deployment 8, when the winds were at 14 knots, with some whitecaps. 93 Chapter 7. SWAPP Measurements 94 Chapter 7. SWAPP Measurements 16- Surf. Temp(C) 12- 00:00:00 06:00:00 12:00:00 18:00:00 00:00:0 17/3 17/3 17/3 17/3 18/3 Figure 7.32: Weather information from SAIL system during deployment 7. Chapter 7. SWAPP Measurements 96 7.2 Environmental sensors Typical data from the environmental sensors in ULSAS are shown in Figures 7.33 and 7.34. These data are from a deployment at 06:00 on 17/3, so they correspond with the end of the higher windspeed period, and thus with the roughest sea conditions we might expect during this deployment. The statistics of these data are given in Table 7.13. We can see from these data and statistics that the bungy cord was doing a reasonably good job of decoupling the array from the surface motions, and that the array was very level and stable. Although the vertical motions were very small, they were still sufficient that they had to be taken into account when processing data: these corrections are discussed later in this chapter. The compass data shows that the array was remarkably stable in a rotational sense as well, considering that there is no force acting other than friction to maintain any particular rotational alignment. Sensor Min Max Mean Std. Dev. Compass (Deg) 281.9 356.0 314.7 16.76 X-tilt (Deg) -3.33 -1.62 -2.30 0.366 Y-tilt (Deg) • -2.63 -0.08 -1.52 0.435 Pressure (dBar) 54.92 56.44 55.44 0.299 Acceleration (g) 0.936 1.054 1.001 0.019 Table 7.13: Statistics of ULSAS environmental sensors for deployment at 06:00 on 17/3. N - 376, covering a time of 6.2 minutes. 7.3 Typical hydrophone data A typical time series of hydrophone data during a "ping" sequence is shown in Figure 7.35. The outgoing pulse is easily identified on the far left, with the next ping beginning on the far right. The outgoing pulses have saturated the receivers. The surface return is the next easily identifiable feature, beginning about 72 ms after the leading edge of the outgoing pulse. Tb the right of the outgoing pulse we see some reflections from other parts of ULSAS, notably the flotation on the instrument housing roughly 20 m below the hydrophone array. The strong return between the surface return and the outgoing pulse (around 60 ms) cannot be accounted for by returns from any part of ULSAS other than cabling (which cannot account for the strength of the return), but is not present in all ping records. Data from more than about 170 ms after the outgoing pulse can be assumed to be representative of ambient sound. To the right of the surface return we see energy which comes from the surface at high grazing angles, but is not the specular return. These are the data we are interested in. Note that the signal level is not Chapter 7. SWAPP Measurements 97 0 1 2 3 4 5 6 TimeCmin) 0 1 2 3 4 5 6 Tlme(mirO Figure 7.33: ULSAS environmental sensor data (Compass, X-tilt and Y-tilt) for a 6 minute period at 06:00 17/3 (Deployment 7). Chapter 7. SWAPP Measurements 0 9 0 l 1 1 1 1 1 1 0 1 2 3 4 5 6 • Time(min) Figure 7.34: U L S A S environmenta l sensor da ta (Pressure and Accelerometer) for a 6 minute period 06:00 17/3 (Deployment 7). 35000 300 Time after ping (ms) Figure 7.35: T y p i c a l transducer voltage t ime series dur ing one p ing . Chapter 7. SWAPP Measurements 99 Figure 7.36: Transducer voltage time series of multiple pings covering a 10 s period. Note the slow variation of the mean level. much higher than the ambient sound level, because we were not able to put very much acoustic energy into the water due to the projector problems discussed in Appendix D. Figure 7.36 shows a time series of pings covering a 10 s period. Note the small but significant slow variation of the mean level due to the VLF response of the hydrophone/pre-amplifier to pressure fluctuations. This variation correlates well with the pressure sensor output over the same period. Note how the strength of the surface return also seems to be fluctuating with the low frequency changes in the mean, because the phase of the surface return is changing as the distance between the mean surface and the hydrophones varies. Since we have a synchronous system, we can detect such small vertical movements by phase fluctuations in the return from the mean surface, rather like in an interferometer. Figure 7.36 points out two requirements in data processing: we must remove the mean from any data extracted on a ping-to-ping basis, or the fluctuations in the mean will overwhelm any other fluctuations in the data; and we will have to correct for pressure fluctuations in order to maintain approximately constant distance (in the sense of the phase of the measurement) between the mean surface and the hydrophones. Figure 7.37 shows overlapped multiple time series of the output voltage from one hydrophone (No.5) for 12 consecutive pings. Backscatter from near vertical to grazing angles of 30° are covered by the time Chapter 7. SWAPP Measurements 100 5000 -10000 1 : 1 1 1 1 1 140 160 180 200 220 240 1/4 Cycles since Ping start Figure 7.37: 12 overlapped transducer voltage time series from consecutive pings, for returns from near vertical to grazing angles of 30°. range chosen. This figure is to be compared with 6.30, which showed a similar set of results for a bottom deployment. We can see from a comparison of these two figures that the consistent pattern within each ping seen in the data from the bottom deployment is not present in data from the mid water deployment, as was predicted. Since there are no stationary scatterers (like sub-bottom layers) near the hydrophones in the midwater deployment, there cannot be any ping-to-ping consistency in the returns. 7.4 Surface scattering strength parameter ( S A ) variation with grazing angle The first parameter we are interested in estimating is the surface scattering strength parameter (5,), which is defined in Urick [53] as: S. = 10 log hcattered (7.123) ^incident where 7 refers to the intensity of sound either scattered from the surface or incident on the surface. Now it can be shown that, for non-directional projector and receiver in the surface scattering mode, that 7.123 can be rewritten as: S5 = RL — SL-r 30 log r - 10 log ( c r v T ) (7.124) Chapter 7. SWAPP Measurements 101 R s T 2 m J I. Z Projector Figure 7.38: Geometry of the scattering strength parameter (Ss) measurements. Z is the depth of the pressure sensor, and the projector and hydrophone arrays are 2 m apart. where r is the pulse length, c is the speed of sound in seawater, and r is the slant range to the scattering point. Because the transmitted pulses are not delta-functions, the measured intensity at any time is actually an average of the intensity of sound scattered from different regions of the surface. Figure 7.38 shows the geometry of these measurements. Suppose the measurement of the intensity of returned sound is made at time to. Then we define the following relationships, referring to Figure 7.38: The point S on the surface is centered on a pulse, the point R is the closest point on the surface from which some of the energy of that pulse can come, and the point T is the furthest point on the surface from which some of that energy can come from. We define the slant range r a s « i from the figure, and the angle 6, which is the complement of the grazing angle, as the angle between s\ and the vertical. We are interested in estimating the value of S, and its variation with the angle 6 as a function of sea state. ct0 = n + r 2 « i + *2 = r j + r2 + — ti+t2 = rx-rr2 + cT. Chapter 7. SWAPP Measurements 102 7.4.1 Corrections for depth variations. The hydrophone/projector array moved vertically in the water, although its motion was damped by the bungy cord. In order to ensure that measurements of scattering from the surface were coming from the same point on the surface, it was necessary to correct for these motions. For each ping, as data was being read in, the depth measurement was examined and compared to a "current" value. If the difference between the two was greater than c/8/o, then the "current" value was incremented or decremented (depending on the sign of the depth change) by c/8/o, and the pointer to the beginning of the next series was incremented or decremented by one, again depending on the sign of the depth change. For typical values of c and /o, a depth change of roughly 0.5 m was required to cause this correction to be implemented. 7.4.2 Mean squared pressure computations. Measurements of mean squared pressure were first used to arrive at estimates of S,. For these com- putations, the digital data from each transducer were processed as follows. First, within each ping, 256-point series were extracted for each transducer: these series started at a specified time after the beginning of the ping, and were corrected for depth variations as described above. Sampling at 4 times per cycle, these time series would be 64 cycles long, or about 2/3 of a ping when using the most common measurement setup. The mean value of these series, for each transducer, was computed and subtracted from each value in the series. The values at each point in each series were then squared, and a running sum of the points in all the series was kept. The results shown here are averages of 1200 pings, or about 5 minutes of data. At the end of the averaging period the 256 point averaged time series of the squared digitized voltage output for each of the 9 hydrophones were corrected for ambient noise. Assuming that the projector was (unfortunately) not powerful enough to produce any measurable backscatter from grazing angles of less than 30°, or 9 > 60°, returns from this region were assumed to be ambient noise. The average of these returns, for each transducer, was calculated and subtracted from the time series. An example of one such averaged, ambient-corrected time series is the solid line in Figure 7.39. Note how the amplitude very quickly drops to zero after the surface return. The region between the surface return and the ambient noise region is the region of interest. To compute Ss, the amplitude of the mean squared digitized hydrophone voltage is used as D, in Chapter 7. SWAPP Measurements 103 6000.000 4000.000 2000.000 1 1 1 1 '• X. v *» „ A '-' j 0.075 0.100 0.125 0.150 Time(s) Figure 7.39: Average amplitude of mean squared pressure (solid line) and correlation magnitude (broken line) vs time after beginning of ping for one hydrophone. These data are averages of 1200 pings. D.144. This expression for RL — SL is used in 7.124. The slant range r is easily calculated from the geometry, and is a function of the time of the measurement being processed. The expressions used for r and ^ are as follows: t = 4/o 0 = arccos (7.126) where n r is the integer index of the data point, assuming that n r = 0 is the beginning of the ping. Figure 7.40 is an example of the computed S, as a function of 0 for several of the hydrophones. These data are from 00:35 on 17/3/90 (File A7F02), very near the start of deployment 7. For comparison 1 have plotted the Chapman-Harris [11] empirical scattering function for windspeed 8 m/s and frequency 393.75 Hz as the dashed line on the plot. It can be seen that the computed values of S, are well above predicted values out to 50°, and then fall below predicted values for greater values of 0. It is difficult to compare the measured results with the predictions because SWAPP measurements were from a source very near the surface, so that measurements at any instant will include returns from a wide range of Chapter 7. SWAPP Measurements 104 5 -45 30 40 50 60 Angle(deg) Figure 7.40: Computed values of Ss as a function of 6 using mean squared pressure for hydrophones 2.5.8 and 9 at 00:35 on 17/3/90. Windspeed was 8 m/s. In this and the next 4 figures, the dashed line is the predicted value for S, using the Chapman-Harris model. grazing angles, due to the finite pulsewidth. For example, the calculated value of Ss at 35° actually includes backscatter from the range 0° < 6 < 35°. This effect is less pronounced at high values of 0, but all of our measurable returns come from the region 6 < 60°. Discussion of measurements of S, will be postponed until after the next section. 7.4.3 Backscatter measurements using correlation analysis. Because the measurement system is coherent and because the signal to noise ratio for the backscatter was so poor, it was decided to try and exploit the coherence to extend the range of useful measurements of 5, beyond where the mean square pressure method fails. The idea is developed in the following paragraphs. First, the transmitted pulse can be written as a sequence of complex numbers: when sampled at exactly 4 times the carrier frequency, assuming the initial phase is zero. Po is some constant value, Ng is the number of cycles in the transmitted ping, and n = 0 is the start of the ping. (7.127) Chapter 7. SWAPP Measurements 105 Let the (real-valued) measured data be the x(n). We can define a new complex variable g(m) as the cross-correlation of the measured data with the transmitted pulse: where N is some number greater than 4Ng. This complex number y(m) is a measure of how the received signal correlates with the transmitted pulse, and the magnitude of this number should relate to the strength of the backscattered pulse. Ambient noise, which will have a non-zero average mean squared pressure, should have an averaged correlation magnitude of zero. It can be shown that g(m) can be computed more efficiently through the use of the Fast Fourier Transform (FFT). If G(l), X(l) and P(l) are the Fourier transforms of g(m), x(n), and p(n) respectively, then we can write: and g(m) is easily recovered from G(l) through the inverse Fourier transform. P*(l) only needs to be computed once. 128 point series were used for both the pulse and the input data series. Thus the input data series were 32 cycles long, or about 80 ms long under the most typical measurement conditions. For each ping the input data series x(n) for each hydrophone were transformed as pure real series using the FFT, and then multiplied by P'(l). The resulting complex series, TV points long, were averaged. When the required number of pings had been averaged, the series were inverse transformed and normalized as described below. The magnitudes of the numbers in the resulting series can then be interpreted as the received signal due to backscatter and used in computations of Ss. Because the use of the FFT implies that all series are repeated infinitely in time, the computation in 7.129 is actually a circular correlation, meaning that contributions from one end of the series can add to contributions from the other end. The usual solution to this problem is zero-padding of the input series. In this case, p(n) is zero over most of its length, so we can just throw out some values of g{m) instead. For example, in the cases looked at here, Ng was 7, so p(n) was non-zero for only 28 points. With p(n) defined exactly as in 7.127 and TV = 128, the first 100 points of g(m) will be correct. In order to relate the values of g(m) to the received signal due to backscatter, we need to relate the magnitude of the g(m) to that of the x(k), by requiring that energy be conserved. We require that: N-l (7.128) G(i) = X(l)P*(l) (7.129) Chapter 7. SWAPP Measurements 106 N-l N-l m) k = 0 m = 0 E *2(*) = E» Now from Parsevals relation we know that A ' - l , A ' - l m=0 1=0 A ' - l N 1 = 0 If we substitute the Fourier transforms in the above, we get (simplifying the summation notation some- what): J2ff2(m) = ^ ^ ^ a - ( i ) e ' ^ ' ^ a ; * ( A - ) e ' ^ u - ^ ? > ( r ) e ' ^ ' r ^ p * ( s ) e ' ^ (7.130) m I j k r ' = i ^ L E E E rU)*~(k)p(r)p*(s) E e x p - h+r~s) j k r / v In the computation of 7.131, we are taking an average over many pings. Everything in 7.131 other than the input signal is a constant, and (x(j)x"' (k)) is zero except when j = k. Thus 7.131 becomes: <E = ^ E E Ev* 2( f c)>pW(«) E e x p ( ' | r '( r - *)) <7-131) m k r s I ^ ' But the in 7.131 is also zero unless r = s, so finally 7.131 becomes (using 7.127): £ > V ) > = *^<I> 2(*)> (7-132) m k We can use .7.132 to relate the g(l) to the output voltage of the hydrophones so we can calculate real values of Ss. To get RL - SL in 7.124, we use 2Poy/Ngg(l)g'(l)/N as Dt in D.144. An example of the value of the normalized value of g(l) can be seen in Figure 7.39. There are several things to note in this figure. First, non-zero values of correlation magnitude extend much further out than do those of the mean-squared pressure: this extension is partly due to the "smearing" of the surface return by the correlation operation, but also to the sensitivity of the correlation analysis to low-level backscatter. Essentially the correlation analysis is giving us some signal processing gain. The second notable feature Chapter 7. SWAPP Measurements 107 of this plot is the null at about 0.12 s: this null is common to most of these analyses, although its position may vary. Beyond the null, the magnitude of the correlation increases again, and then tends to zero as the backscattered returns become too weak to detect. 7.4.4 Backscattering measurements using both methods. Data from 3 collections during deployment 7 were used in this comparison: one set near the beginning of the deployment, at 00:35 17/3/90 (File A7F02), when the winds had been blowing from the South at 8-10 m/s for about 8 hours; one set at 06:00 17/3/90 (File A7F04), just before the Southerly 8-10 m/s winds began to drop off; and one set from 22:00 17/3/90 (File A7F22), when the winds had been calm for roughly 6 hours. Sea conditions during these 3 collections were quite different: the first set had choppy seas with whitecaps and no swell; the second set had choppy seas, whitecaps, and a low SW swell; and the final set had smooth, oily seas. The correlation method takes what is essentially a weighted running average of the data. The width of this average is 4Ng data points. In order to compare the estimates of Ss using mean square pressure measurements to those using the correlation method, it was necessary to take a similar running average of the mean-square pressure data. The results of these computations, using data from 6 hydrophones, are shown in Figures 7.41 through 7.43. The dashed line in these figures is the computed Ss using the Chapman-Harris [11] empirical model. There are several points to note on examination of these figures. The main point is that the computed values of Sa are significantly higher (as much as 20 dB) than those estimated with the Chapman-Harris model. This increase is valid for both methods of estimation for 0 < 45°, and for all angles using the correlation method. Secondly, the two methods yield very comparable values, although the correlation method has much more hydrophone to hydrophone variability. The slope of the variation in Ss with angle using the correlation method seems to be close to that of the Chapman-Harris model, but there is unexplained structure in the angular variability of Ss. Thirdly, the correlation method does seem to allow useful estimates of Sa to be made in regions where the mean square pressure backscatter signal is "hidden" in the ambient noise. Finally, there is no clear pattern of variablilty of Ss with sea state. For example, one would expect that S, for near vertical incidence would be higher when the sea is smooth than when it is rough, but it does not appear to be so, comparing Figure 7.41 to Figure 7.43. It may be that the sea conditions were simply not different enough to cause any noticeable effects on Ss: for example, the Chapman-Harris Chapter 7. SWAPP Measurements 108 30 40 50 60 30 40 50 60 Angle(deg) AngJetdeg) Figure 7.41: Computed S, values for data from 00:35 17/3. Results of mean squared pressure computa- tions on left, correlation analysis on right. Angte(cteg) Angte(deg) Figure 7.42: Computed S, values for data from 06:00 17/3. Results of mean squared pressure computa- tions on left, correlation analysis on right. Chapter 7. SWAPP Measurements 109 30 40 50 60 30 40 50 60 Ang&deQ) Angle(den) Figure 7.43: Computed Ss values for data from 22:00 17/3. Results of mean squared pressure computa- tions on left, correlation analysis on right. model predicts a 2.5 dB change in Ss at 45° for a change in windspeed from 8 m/s to 2 m/s. This change is small enough that it may be hidden in the hydrophone to hydrophone variablility seen in these data. 7.5 Ping by ping time series. In order to look at changes in backscatter with time, ping-by ping time series of complex voltages were extracted from the raw data. It is useful and efficient to look at several areas on the sea surface within each ping. Those areas should be separated in time by one ping length, so that the returns will be independent of one another. First the raw data had to be filtered, at least to remove the DC variability from ping to ping, by filtering in the frequency domain. For each ping and for each hydrophone, 32 point time series of complex voltages were extracted from the raw data, as described in 5.2.1. The time step for these series was then l//o, so the length of each series was 32/fo, or 81 ms for the most commonly used carrier frequency of 393.75 Hz. These series were then Fourier transformed, and multiplied by a simple frequency domain filter. The frequency resolution of the filter was obviously very gross, i.e. 12.3 Hz for f0 = 393.75 Hz. The filter chosen was a band-pass filter, where frequencies from ±12.3 Hz to ±73.8 Hz were passed, but energy at higher frequencies and at DC was removed. This filter characteristic was chosen for a number of reasons: first, to remove the ping-to-ping DC variability from the data; secondly, to remove high frequency energy which cannot be related to sea-surface changes; and finally, the filter Chapter 7. SWAPP Measurements 110 must be left open enough to allow changes occurring over the time scale of the ping length to pass, so that we can get information from more than one part of the surface with each ping. The resulting complex frequency domain series was then inverse Fourier transformed back to the time domain, and individual data points were taken and stored. Points 2,9,16,23, and 30 were saved for each hydrophone (note the separation of 7 cycles, or one ping length, between points), for each ping. This filtering and decimation process reduces data quantity by a factor of about 40, and yields data which we can use to examine the long term variability of backscattering from the sea surface at 5 "range gates", each of which corresponds to mutually exclusive annular regions on the surface. Once the software to carry out this operation is developed and tested, it is easy to vary the characteristics of the filter, and the regions on the surface to be examined. The time series of complex voltages can be analysed in the frequency domain to give information about Doppler-shifted backscattered sound in the band ± 2 Hz. 7.5.1 Voltage statistics, and returned power levels. Each range gate time series results from reflected power from a particular annular region of the sea surface, together with any volume reverberation (assumed to be very small at the low frequencies used), and ambient noise. Most data were taken starting at cycle 30 after the start of ping. The incidence angles and Bragg frequencies corresponding to the gated time series are shown in Table 7:14. 0 and / are the values at the center of the reflected ping. Note that the first gate includes mostly power from the direct reflection from the sea surface. The values of 6 and / vary widely for the first few gates, since the ping is 7 cycles long, and the instrument is quite close to the surface. Gate No. Cycle No. @max e fmin /max / 1 30 0.0 16.0 0.0 0.0 AHA 0.0 2 37 16.0 39.2 31.1 .4774 .7224 .6525 3 44 39.2 49.5 45.1 .7224 .7922 .7645 4 51 49.5 56.0 53.1 .7922 .8272 .8123 5 58 56.0 60.6 58.5 .8272 .8479 .8388 Table 7.14: Incidence angle and Bragg frequency for each range gate for a typical collection starting at cycle 30, and a carrier frequency of 393.75 Hz. Barrick and Snider [5] have shown that sea-echo voltage signals are Gaussian (normally distributed) random variables, so long as the sea height can be represented as a Gaussian random variable. Thus as a first test of the validity of the filtered, range-gated voltage time series, I looked at the statistics Chapte r 7. SWAPP Measurements 111 -10000 -5000 0 5000 10000 -10000 -5000 0 5000 10000 REAL PART IMAGINARY PART Figure 7.44: Probability plots of 200 consecutive values of real and imaginary voltages for gate 2 (6 — 31.1°), referred to a normal distribution. of the real and imaginary parts of these data. Figure 7.44 shows a typical example of these voltages, from a range gate where backscattered power was reasonably high. These plots are of the proportions of the voltage values compared to the expected proportions for a standard normally distributed random variable. Thus if the voltages are normally distributed, they will lie on a straight line. We see that the voltages do fall more or less on a straight line, but we also note that they are not zero-mean, which means that the frequency domain filter, while it will remove the mean value of the 32 point data series within any ping, cannot remove ping-to-ping variability which may show up as a small offset in the mean value of the time series. Additionally, we can see that the outliers in the data do not conform to a normal distribution, possibly because the assumption of normality cannot really extend to very large values in the return, since the assumption of Gaussianity for the sea height obviously cannot be extended to the "tails" of the distribution. I have stated before that the projector could only be run at 1/16 of its maximum possible output power: consequently, the strength of the backscattered echo was low. One indication of this lack of signal strength can be seen in the plots of mean squared pressure (Figure 7.39). Computed statistics of the range-gated data also show this problem. Figure 7.45 shows the variation in the mean amplitude of the complex voltages in the 5 range gates, plotted against 6, for a data collection made at 06:00 on 17/3/90 (File A7F04), when the sea state was as rough as it got during deployment 7. These points are Chapter 7. SWAPP Measurements 112 -o- 0 0 20 40 60 Angle (deg) Figure 7.45: Mean amplitude of complex voltage for 5 range gates, plotted against 9. The dashed line is approximately the amplitude of ambient noise. averages of 2560 pings, or over 10 minutes of data. The dashed line in this figure is the approximate level of ambient noise. Clearly, only range gate 2 has any appreciable power from surface returns in it. The consequence of this lack of reflected power is that we are unlikely to see any information about the sea surface in the time series from range gates 3 to 5. 7.6 Frequency spectra of Ping data and ambient noise. Software was developed to enable 128-point real-valued time series from any location within pings to be collected foT all available hydrophones, transformed to the frequency domain, and averaged over any number of pings, all in real-time (i.e. as the data tapes are being played back). Using these data we can analyse the frequency content of the signal at any point in the ping. Using Jo = 393.75 Hz and sampling 4X per cycle, a 128 point time series will enable frequency spectra from DC to 787.5 Hz to be collected, with a frequency resolution of 12.3 Hz. Figure 7.46 shows two of such spectra for one hydrophone (No. 7) for two different regions in the ping: cycle 25-57, which includes the direct surface return; and cycle 37-69, which starts just after the direct surface return ends. These power spectra are averages of 2000 spectra, or about 8.3 minutes of data collection. The data was from File A7F04, when the sea state was rough. We can see, as expected, very significant energy at the carrier frequency in the plot for cycles 25-57. We also see a lot of energy between the carrier frequency and the first harmonic, and between the carrier and DC. Comparing the two plots, we see some significant differences. Now the spectrum for cycles 25-57 will include energy from the surface returns that we see around the carrier in the plot for cycles Chapter 7. SWAPP Measurements 113 170 170 100 100 0 100 200 300 400 500 600 700 600 0 100 200 300 400 500 600 700 800 Frequency (Hz) Frequency (Hz) Figure 7.46: Frequency spectra for hydrophone No. 7 for two different regions of the ping. Left plot is cycles 25-57: right plot is cycles 37-69. 37-69. The differences between the two plots must be due to the energy in cycles 25-37, which is the region of the direct surface return. These data seem to show that the energy in the ping is smeared considerably by backscatter from the surface. It is not clear what mechanism exists to do this smearing: Doppler shifts from the surface are predicted by theory to be only on the order of a few Hertz. The width of the main peak around the carrier frequency is very roughly what would be expected for a 7 cycle pulse at 393.75 Hz: we would expect a width of around 60 Hz for this peak. 7.6.1 Ambient noise spectra. Data collected with the transmitter disabled can be used to look at the spectrum of ambient sound. In order to look at the broadest possible bandwidth of ambient sound, a carrier frequency of 1102.5 Hz was used, sampling at 4X per cycle. Because of the high sampling rate only 4 hydrophones could be sampled. Results presented here are typical of these four hydrophones. The data analysed were collected at 22:00 on 17/3/90 (File A7F22) and at 23:30 on 18/3/90 (File A8F05). The first was a period of very calm seas, with a low Southwesterly swell and virtually no wind blowing; the second collection was made during a period of 14 knot Easterly winds, with some whitecaps on a low Westerly swell. Figure 7.47 shows the resulting spectra: these are averages of 1000 unwindowed spectra (about 4 minutes), with a frequency resolution of 34.5 Hz. The spectra have been corrected for the frequency response of the pre-amps (see Appendix C). Chapter 7. SWAPP Measurements 114 120 1000 1500 Frequency (Hz) 1000 1500 2000 2500 Frequency (Hz) Figure 7.47: Frequency spectrum of ambient noise for hydrophone No. 5. Data are averages of 1000 unwindowed 128-point series. Left-hand plot is from File A7F22 (calm); right-hand plot is for File A8F05 (14 kt winds). The shape of these spectra follows generally accepted models of deep-water ambient noise (see for example Urick [53]). The absolute levels of the spectra are high, but these plots are not corrected for self-noise, which is estimated to be about 65 dB. Incorporating this correction brings the levels close to values reported in the literature. Comparison of the two spectra shows the expected increase in spectrum level in the higher windspeed regime. This increase is especially noticeable in the 200 - 600 Hz region, roughly where the peak of the broad spectrum of ambient noise due to wind has been reported in the literature. C h a p t e r 8 S u m m a r y a n d C o n c l u s i o n s A theoretical framework has been developed to enable a test of the acoustic CODAR technique to be carried out with an acoustical array, using beamforming techniques to simulate the CODAR antennas. Expressions relating the CODAR 3-element antenna outputs to the output of an array of omnidirectional acoustic point sensors have been developed. Mathematical algorithms and techniques have been derived to extract information about the directional waveheight spectrum (DWS) of surface gravity waves from acoustic Doppler backscatter measurements with the array. These algorithms, or data inversion methods, include techniques for extracting information in the presence of noise. Data processing routines to implement these techniques have been written and tested. To test algorithms and gain insight into the form of data products we might expect from real mea- surements, several models have been developed mathematically and in software to show the expected form of the power spectral density of the acoustic Doppler backscatter seen by single omnidirectional receivers, and the expected form of data products of the beamformed array. These models use 2 different representations of the DWS of surface gravity waves. The models were used both to test data processing routines and as a point of reference in assessing the usefulness of real data. Starting from an earlier version (ULSAS 1), an acoustic instrument (ULSAS II) has been developed for stand-alone, remotely controlled operation in both bottom-situated and deep-water, surface-tethered configurations. Using commercially available digital audio equipment for data digitization and storage, this device can collect and store large quantities of acoustic data from a multi-element array, under the control of a distant operator over a radio link. ULSAS is designed to remain deployed for long periods of time (up to 3 weeks of battery power are available), and can collect data automatically according to preset criteria, or under the direct control of the operator. Fail-safe modes of operation are built in, ensuring that data will be collected in the event of a breakdown of the communications link. The bottom-situated version has been deployed in the coastal waters of British Columbia, and the deep water version, which is the result of recent major design changes, was deployed on a "ship of opportunity" basis in the recent SWAPP experiment. 115 Chapter 8. Summary and Conclusions 116 A preliminary test ol" the acoustic CODAR technique has been made using the earliest version of ULSAS in Saanich Inlet, near Sidney, B.C. Difficulties with the instrument and the deployment site precluded a full, unambiguous test of the technique, but some encouraging results have been obtained. Application of the data inversion methods to samples of the Saanich Inlet data yield information consis- tent with the known wind and wave field. The principal wind and wave direction can be ascertained, and the form of the non-directional part of the DWS follows approximately the expected k~A shape for k val- ues above saturation. Beamforming results using frequency-domain data show that the Doppler-shifted acoustical backscatter is directional in nature. These are the first results of this kind to be reported. Results of a deployment of the bottom-situated configuration of ULSAS II showed that anomalously strong scatter from the bottom (or possibly sub-bottom) exists at the frequencies used, making bottom deployment of the instrumentation only marginally useful. Signal processing hardware and software were developed to allow real-time processing of the large amount of data available from ULSAS II. A method of estimating the surface scattering strength param- eter Ss using a correlation of the received signal with the outgoing pulse was developed, in an attempt to extend measurements of Sf into regions of low signal to noise ratio. Techniques developed and tested in this phase of the work may be applicable to later versions of ULSAS in which much of the data processing will take place in the instrument, allowing a remote user to get measurements of parameters of interest directly from ULSAS in real-time. The deep-water version of ULSAS 11 was tested for the first time during the SWAPP cruise. Except for a major problem with the acoustic projector and its driving system, the instrument performed well: valuable experience was gained in the deployment and operation of ULSAS at sea. Before this deployment, ULSAS had always been operated from a land base, and was deployed on a long-term basis. During the 3.5 weeks of the SWAPP cruise, 9 successful deployments of ULSAS were made. In spite of the problem with the projector, which limited its power output, to only 1/16 of the maximum available, we were able to make estimates of the surface scattering strength parameter over angles of incidence less than 45°. These results showed some surprising departures from the Chapman-Harris [11] empirical formula for S, , and interesting angular structure. Measurements of the ambient noise field were also made under calm conditions and during 14 kt winds. The SWAPP cruise was an excellent "test bed" for the deep water version of ULSAS II. This instru- ment has been shown to be basically sound in design and construction, and with a few relatively minor changes and "bug fixes" should be able to make some very useful measurements of low frequency acoustic Chapter 8. Summary and Conclusions 117 scattering processes in the sea surface zone, in addition to finally testing the validity and usefulness of the acoustic CODAR technique. Appendix A Computation of F2 using "REDUCE" This is the output of " R E D U C E " when computing the expression for antenna power. The Vi are the complex transducer voltages, and the Kj are defined in Chapter 3. 118 Appendix B Expressions for the data products B„. The following are the expressions for the Bn as a function of the transducer voltages V,-, where the hydrophone locations are shown in Figure 3.2. h is defined in Chapter 3. It is simplest to write these expressions in terms of the voltage cross-products (here Si indicates the real part of a complex number, and 3 the imaginary part): Let Pa and Qij Then P, ; = />•,• ; Qij — Qjii and Qi, = 0. Using these relations, the Bv can be written as: Bi = — (P22 + P44-P11- P33 - 2P2A + 2P13) 1 07T (B.133) B ^ = » / ^ ' " f ^ ( - Q l 2 + <?23 + Qu + <?34 + 2Q 2 4) 2 ( 4 - ^ ) (B.134) (B.135) 73_! = , / ' S m f 2 . ( Q 1 2 + Q 2 3 + Ql4 - Q34 + 20i3) 2 (4 - fr) (B.136) /j2 5-2 = (Pl2 + ^34 ~ P23 ~ Pu) (B . l 37' 119 Appendix C A Detailed Description of ULSAS I I The first section in this description of ULSAS II describes mechanical aspects of the deep water, surface- tethered version of the instrument, as well as giving some information about its deployment at sea, and a general overview of how the acoustic and electronic components work. The remainder of the sections discuss various aspects of the acoustic and electronic components of the system in more detail. Other than the details of physical configuration and deployment which obviously have to do with the deep water version, this description applies equally well to both the bottom mounted and deep water versions of ULSAS II. C . l General system description. C . l . l Mechanical description of the deep water version of ULSAS I I . Figure C.48 shows a scaled side view of ULSAS II, as deployed during SWAPP. The surface buoy was a "doughnut" buoy with a 2 m aluminium tripod mounted on it. The top of the tripod formed a radar reflector. The doughnut buoy was chosen because it gave considerable buoyancy at the surface, and the attached tripod was a good place to mount a large splashproof box which held electronic gear for communication between the user aboard ship and the instrument suspended below the buoy. This arrangement also improved our radio range by allowing the radio antenna to be nearly 2 m above sea level. In addition to the box, a xenon flasher, 2 RDF beacons, a Coast Guard approved buoy light, and an Argos satellite positioning transmitter were attached in various ways to this package. In order to decouple the hydrophone array from surface motion, a 15 m bungy cord with a spring constant of 120 Nt • m - 1 , with a safety line attached, was connected from the surface buoy to a 27 m length of Nilspin cable, which in turn connected to the top of the hydrophone array. The net weight in water of the complete package was designed to be 160 kg, which would stretch the bungy cord to a length of 28 m. Thus the depth of the hydrophone array was designed to be 55 m. Next in line is the hydrophone array/projector mount, shown in Figure C.49. This mount was an 120 Appendix C. A Detailed Description of ULSAS II - - S u r f o c e b u o y 50 m Q — Hydrophones and projector* 20 n U L S A S , b a t t e r i e s . V C R R e l e a s e i s 30 n b e l o w U L S A S 121 Figure C.48: Side view of ULSAS deployment. open aluminum frame, with 9 receiver hydrophones arranged in a 3 by 3 array close to the top. The inter-element spacing was approximately 23 cm. The acoustic projector was mounted at the bottom of the frame, so that the separation between the hydrophones and the projector was 2 m. This separation ensured that the hydrophones were not in the near field of the projector. The hydrophones were mounted in cylinders cut from high density open cell polyurethane foam, which in turn were held by friction in short sections of aluminium pipe. This mounting method was chosen to isolate the hydrophones from mechanical vibrations in the mounting frame or cables. There was also a small pressure case containing environmental sensors (compass, tiltmeters, pressure sensor and accelerometer) and associated electronics attached to the frame at the bottom. Connecting the hydrophone array to the electronics/battery package was a 20 m length of Nil- spin cable, to which a bundle of 13 electrical cables was attached. The Nilspin cable had 3 "hard eyes"incorporated in it, so that the hydrophone array could be picked up in a series of lifts during recovery. This method of recovery had to be used because the array could not have any flotation on it. Next was the electronics/battery package, seen in Figure C.50, which was also an open aluminium Appendix C- A Detailed Description of ULSAS II 122 Figure C.49: Side view of the hydrophone/projector mount. frame, to which were attached a large spherical pressure case and two cylindrical pressure cases. The large spherical pressure case contained most of the electronics, and several batteries. The smaller cylindrical pressure case contained a video tape recorder (VCR), and the large cylindrical pressure case contained the batteries used for the power amplifier main power supply. A lifting frame was attached to the top of this assembly. For flotation, 12 plastic floats, 30 cm in diameter, were attached in the space between the ball and the top of the lifting frame. This flotation resulted in roughly 75 kg net buoyancy to the whole assembly. A submersible xenon flasher was attached to this package at the top of the lifting frame. Below this package was approximately 30 m of Nilspin cable, attached to a release frame with an acoustic release and a safety pressure release connected in parallel. Above it were 3 plastic floats for flotation. The bottom weights were partially full barrels of concrete weighing roughly 240 kg in water, attached to the release frame. Measurements with weight attached showed that the net weight of the Appendix C. A Detailed Description of ULSAS II 123 Figure C.50: Side view of the electron)cs/battery case assembly. Appendix C. A Detailed Description of ULSAS II 124 whole assembly in water was about 173 kg, so the Bungy cord was stretched to 29 m, and thus the hydrophone array was at roughly 56 m below sea level, very close to the. design depth. C.l.2 Deployment and recovery. The instrument was deployed from the well deck on the CSS "Parizeau". Deployment proceeded in stages, starting with the surface buoy. When it was over the side and the ship had moved slightly away, the hydrophone array was picked up from the bottom, slung over the side, and released when it was totally submerged. The electronics/battery package was picked up and put over the side next. This stage had to be completed quickly, because the hydrophone array was hanging from it by the cable bundle, and we had to avoid any strain or chafing on the electrical cables. Then the hydrophone array hung between the two buoyant parts, the surface buoy and the electronics/battery package. Since the cable from the surface buoy was longer, the array still hung upside down in the water. Next the release frame was lifted up and a weight attached to it. This assembly was lifted up, and the ship backed off slowly to ensure that all components were strung in a line. At that point the release frame/weight assembly was let go, and the string of components slowly sank to end up suspended below the surface buoy. The depth of the pressure sensor on the small pressure case on the hydrophone frame was typically around 56 m. Recovery was somewhat more difficult, especially in rough weather. The acoustic release was trig- gered, and the excess buoyancy on the electronics/battery package brought it to the surface, usually fairly near the doughnut buoy. Depending on which was closest to the ship, we then picked up either the surface buoy or the electronics/battery package. Most of the time the surface buoy came aboard first. When both buoyant components had been picked up and were on deck, the hydrophone array frame then hung below the ship by the cable bundle connecting it to the electronics/battery package. It was brought on board in stages, by attaching the lift point to the nearest eye in the cable, and lifting to the next eye, which was then stopped off. The lift point was then transferred to that eye and the process repeated until the array frame could be brought aboard in one lift. This scheme worked very well — there were 11 deployment/recovery sequences carried out, and the most damage that was done was the loss of one RDF beacon due to a hard bang of the surface buoy against the side of the ship. Credit for this success is due partly to the mooring design and construction, but mainly to the skill of the ships officers and crew. Appendix C. A Detailed Description of ULSAS II 125 Remote PC f \ UHF Link \ ) 3uoy/UHF Link ~7F U/W cable C 6BHC11AB Serial I/O 1 ROM/ RAM B u s Hyjr. *5 Digital I/O rr I V Control I r = f \ , , K Stat* Audio Uachine ' Clock, etc. Serial I/O Analog in Clocks Proj- ector PCM VCR Video s ) Figure C.51: "Black-box" representation of ULSAS II electronic subsystems. C.2 Generation of the Acoustic Signal. This section begins the discussion of the electronics of ULSAS II. Figure C.51 is a block diagram of the various subsystems of ULSAS II which will be discussed in this and the following sections. C.2.1 Sine wave generation. To generate the gated sinusoidal output pulse, a 256 step lookup table holding the 8-bit 2's complement integer representation of one full cycle of a sine wave is stored in read only memory (ROM). Figure C.52 shows a plot of one cycle of the integer approximation to a sine wave, and its power spectrum. Note that the maximum sideband level, due to the inaccuracy of the integer representation, is about 60 dB below the signal level. The ROM used is a 2732, capable of holding 4096 bytes of data: it has twelve address lines, designated AO through All . The contents of the ROM are output at a rate controlled by a programmable divider, Appendix C. A Detailed Description of ULSAS II 126 Figure C.52: Integer representation of sine wave (upper plot), and its power spectrum (lower plot). The horizontal axis of the lower plot is proportional to frequency, where the fundamental frequency of the sinusoid is at index 1. Appendix C. A Detailed Description of ULSAS II 127 PRELOAD VALUE (HEX) FREQUENCY (Hz) EC 275.63 FO 344.53 F2 393.97 F4 459.38 F6 551.25 FA 918.75 FD 1837.5 FF 5512.5 Table C.15: Counter reload values and corresponding frequencies. which divides the master clock by an integer number to produce a sine wave at the frequency closest to that desired by the operator. The master clock signal, at 1.411 MHz, goes to a programmable 8-bit up-counter which is preloaded with an operator-selected value every time it overflows. The carry-out line of this counter is then used to clock another 8-bit counter which counts up from zero. The outputs of the second counter are used as the lowest 8 bits of the ROM address (AO - A7). Thus as the counter increments from 0 to 255, the contents of the ROM are output sequentially. As an example of the selection of the preload value for the first counter, suppose the operator wants a frequency of 400 Hz. The following relationships: 1.411 x 106 % Nd x 256 x 400 ; 1.411 x 106 = Nd x 256 x f0 ; and Preload value = 256 — Nd govern the choice of jVd(which must be an integer) and the preload value. The integer which most nearly satisfies the first relation above is 14. Therefore the preload value will be 256 — 14 = 242 (F2 in hexadecimal), and the actual frequency produced will be 393.97 Hz (from the second relation). Table C.15 shows some preload values (in hexadecimal) with their corresponding frequencies. To produce a ping, the sinewave output is turned on and off, by controlling a "page bit". The ROM used is a 4 Kbyte ROM, so it has many more memory locations than are needed to hold the 256-byte sine wave. Address lines A9 - All are held low, and A8 is used as the "page bit". ROM locations 100 - IFF hold the sine wave, and locations 000 - OFF are loaded with zeroes. Then the output of the ROM depends on the state of address bit A8. If this bit is a one, then the output of the ROM will be the sine wave, and if it is zero, the ROM will continuously output zeroes. Thus a ping consisting of 6 sinewave cycles followed by 94 cycles of no signal is produced by turning A8 on for 6 cycles and then off for 94 cycles. Appendix C. A Detailed Description of ULSAS II 128 VALUE(HEX) GAIN (dB) FF 0 80 -6 50 -10 IA -20 Table C.16: Signal gain as a function of attenuator value. C.2.2 Analog processing. The ROM output goes to a digital to analog converter (DAC) which produces a ± 5V full-scale analog output. To reduce sideband levels still further, this signal is passed through a linear 4 pole low. pass filter with a Bessel transfer characteristic and a cutoff frequency of 1.2 kHz. The filtered signal then passes through a programmable attenuator before being input to the power amplifier. By varying the attenuation of the input waveform, the power output of the acoustic projector can be controlled (see Section C.2.3 below). Examples of gain as a function of the digital attenuator setting are shown in Table C.16. C.2.3 Acoustic projector and power amplifier. The acoustic projector used on ULSAS II in the deep-water version was an Argo Technology Model 201 low frequency source. This device has a useful operating range of 400 to 3000 Hz, with a resonant frequency of 2 kHz. The typical transmitting voltage response is shown in Figure C.53. This projector is driven through a 30:1 step-up autotransformer by an Amcron M600 power amplifier (PA), which has a maximum power output of 1000 watts into a 4 ohm load over a bandwidth of 1 Hz to 15 kHz, and maximum total harmonic distortion rated at 0.05% of the output voltage. As mentioned in the previous section, the power output, of the PA, and thus of the acoustic projector, is controlled by varying the PA input level. Projector output can be disabled by turning off the PA power supplies. C.3 Acoustic Receiving System. C.3.1 Receiving hydrophones. The nine receiving hydrophones are ITC-1001 omnidirectional devices with a typical open circuit receiv- ing response of -187 dB re 1 V / i P a - 1 at 400 Hz. The hydrophone outputs are amplified approximately 27 dB by pre-amps potted into cables close to the hydrophones. The typical frequency response of the Appendix C. A Detailed Description of ULSAS II 129 o > 0) c Q -O 3. 03 1.0 Frequency kHz Figure C.53: Typical transmitting voltage response of the Argotech projector. pre-amps is given in Figure C.54. Outside of the passband, the frequency response of the pre-amps falls off at about 6 dB/octave. C.3.2 Analog signal processing. The pre-amp outputs pass through programmable gain amplifiers before digitization. Amplifier gain can be set to 0 dB. 12.04 dB, 24.08 dB or 36.12 dB, so the maximum overall system gain, from hydrophone to digitization, is about 63 dB. There is also a programmable switch at the output of each amplifier, which switches the output to either the amplifier or analog ground. This switch is used in system tests to measure internal electrical noise level. C.4 Data Aquisition and Storage. C.4.1 Sample/hold and multiplexer. All signals pass to a bank of sample and hold (S&H) modules which are switched in parallel. At each designated sample time, the S&H modules are put into the "hold" mode, so that all signals are sampled synchronously. Digitization of the S&H outputs then proceeds in a user-determined sequence at 44.1 kHz. Droop rate of the modules used is specified at 500 p.Y- ms, so that the maximum droop in any Appendix C. A Detailed Description of ULSAS II 130 1 10 100 1000 10000 100000 100 158 251 398 631 1000 Frequency (Hz) Frequency (Hz) Figure C.54: Typical frequency response of the pre-amps, covering the range 1 Hz to 50 kHz. Right-hand plot is a blow-up of the region between 100 Hz and 1 kHz. Note the different vertical scales in these plots. channel, assuming that all 15 MUX channels are being sampled consecutively at 44.1 kHz, is 170 pV. For comparison, a 1 bit change in the ADC (Analog to Digital Converter) output, given that the 16 bits of accuracy span a range of ±5V, is about 150 pV. Analog data from the hydrophones and other sensors passes via the S&H bank to a 16 channel multiplexer (MUX). Channel 0 is never used, so there are up to 15 channels available for data. The MUX addresses used and their corresponding data sources are shown in Table C.17. C.4.2 Digital audio (PCM) system. Analog to digital conversion of the multiplexed signals is carried out using a modified digital audio PCM system manufactured by SONY(TM). The acronym "PCM" stands for Pulse Code Modulation, a method of encoding digital data using a sequence of pulses. This device will digitize two channels of audio signals (Left and Right) with 16 bit resolution at 44.1 kHz per channel. The major modification made to the PCM is that only the Left channel is used for the digitization of analog data. Right channel storage is allocated to information about the data source and the sequencing of the data. The least significant 4 bits of the 16-bit Right channel word hold the MUX address of the data in the Left channel 16 bit word, and the upper 12 bits are used to hold a count which increments for every new sample taken. For example, if the sampling setup requests that transducers 3,5,6, and 8 be sampled, then a list of data (in hexadecimal) might look like: Appendix C. A Detailed Description of ULSAS II 131 ADDR.ESS(HEX) DATA SOURCE 01 Xducer 1 02 Xducer 2 03 Xducer 3 04 Xducer 4 05 Xducer 5 06 Xducer 6 07 Xducer 7 08 Xducer 8 09 Xducer 9 OA Compass OB X axis tilt OC Y axis tilt OD Pressure OE Accelerometer OF Ground reference Table C.17: MUX channel assignments. RIGHT L E F T C028 data 8 C033 data 3 C035 data 5 C036 data 6 C038 data 8 C043 data 3 C045 data 5 C046 data 6 C048 data 8 C0S3 data 3 Note how the first three hex digits of the Right channel work are a count that increments once for each set of samples taken. The last digit is the MUX address of the channel whose data follows in the left channel word. Other PCM modifications are to the Left channel analog input circuitry, including: • DC coupling of signals rather than AC coupling; .Appendix C A Detailed Description of ULSAS II 132 • removal of low-pass anti-aliasing filter. The PCM is designed to digitize audio signals at a 44.1 KHz digitization rate, so it has a "brick-wall" low-pass anti-aliasing filter with a cutoff frequency of around 22 KHz in front of the ADC. However, the signals we will be presenting to the ADC are coming from a multiplexer which will be switching from channel to channel at 44.1 kHz. If we left the anti-aliasing filter in, the channel to channel variability of the signal, which contains critical information, would be lost or seriously affected. • modification of gain of input amplifiers to result in a ±5V full-scale range for the ADC. C.4.3 P C M Clocks. The PCM system produces 3 clock signals which are used by the local controller (see Section C.6.2) to synchronize carrier frequency generation, ping production, and data collection. B C L K The two 16 bit wide words from the left and right channels are transmitted as a serial data stream, one after the other. The serial bit rate is determined by BCLK, and must be 44.1 kHz x 16 x 2 = 1.441 MHz. 2WCLK This clock is BCLK divided by 16 — it is the word rate of the data output by the PCM, which is 88.2 kHz. Proper interpretation of the serial data requires this clock, since it indicates where 16 bit words begin and end in the serial data stream. W C L K This clock has a frequency exactly half that of 2WCLK. It indicates the source of the data, either Left or Right channel — When WCLK is high, the serial data line contains Right channel data, and when it is low, data is from the Left channel. C.4.4 The Video recorder. The serial data stream from the left and right channels is converted by the PCM to a video signal which is recorded on a regular VHS type video recorder (VCR). The only modification of the VCR was to Appendix C. A Detailed Description of ULSAS II 133 disable the automatic dropout compensation circuitry — if we have data dropouts, we want to know about it! In addition to the video information, the VCR. can store two HI-FI audio channels and one monoaural audio channel. Only one VCR was used: with one T-160 videotape loaded a maximum of 8 hours of data can be recorded, using the extended play feature of the VCR. VCR functions were controlled through the remote control input circuitry, which is accessed directly, bypassing the infrared link. Only the POWER ON/OFF, PLAY, RECORD, and STOP functions were used. C.4.5 Data sampling configuration. The variables which must be specified for a data collection are: 1. Carrier frequency (fo)- This is determined by the choice of a counter preload value (Nd) as described in Section C.2.1. Recall that f0 is given by: / 0 = 1.411 MHz / (256 x Nd). 2. Number of transmitter 'on' cycles at fo per ping. (Ng) 3. Number of transmitter 'off: cycles at fo per ping (Nc). Then the ping duration in seconds can be calculated from: Ping length = (Ng + Nc)/f0. 4. Number of pings to produce (Np). 5. Number of times per carrier frequency cycle to sample (5'/). This number can range from 1 to 64. For example, sampling at 4 times per cycle means that samples are taken at times: t0,t0 + 1/(4 x /„),!„ + 2/(4 x /o),*0 + 3/(4 x / 0) where to — A''/fo ; and Ar is an integer. The next sample to be taken in the above sequence would be at time (A' + 1)/fo- Sampling frequency is fundamentally limited by the maximum digitization rate of 44.1 kHz, as will be discussed further below. 6. Total number of samples to collect (Cc). This number is dependent on all the others listed above. In order to collect samples for all pings transmitted, it should be set equal to: Cc = (Ng + Nc) x NP x Sf. Appendix C. A Detailed Description of ULSAS II 134 However, the ping production and sample collection processes are independent of each other — each is controlled by a separate microcontroller — so it is possible that ping production could continue with no data collection, or vice versa. 7. MUX addresses to digitize for each sample. These will range from 1 to 15, and up to 15 can be specified in any order. Let Nm be the number of MUX addresses to be digitized per sample. In determining the maximum sampling rates allowable, the following inequality must be adhered to: fo x Sj x Nm < 44.1kHz This inequality can be expressed more simply as: (Sf x Nm)/Nd < 8 The operator can choose and set both current and default values for all of the above variables, but there is no error checking carried out by the system on these settings. C.S Communication and Control. ULSAS II is self-contained, and is deployed at a shallow depth (typically 50 m). The user communicates with the instrument over a radio link between a surface buoy and a surface station. The radios used were 1COM IC-U16A UHF transceivers operating at 412.2125 MHz, and communication at 1200 baud was via a packet switching technique using AX.25 Level 2 protocol, which means simply that data are transmitted in short blocks or packets which include error checking information. In the event of errors in transmission, the packet will be retransmitted until it is error-free. Communication is initiated by the issuing of a CONNECT request to the buoy terminal node controller (TNC) by the surface station TNC. After the two TNC's are connected, their packet-switching operations are transparent to communication between the main controller in ULSAS II and the user's computer aboard ship. An AEA Model PK-88 TNC was used on the surface station, and a Heathkit HK-21 TNC was used in the buoy, where low-power operation is a necessity. Data transfer between ULSAS II and the surface buoy is via a two-wire current loop in an underwater cable attached to the wire linking the buoy to ULSAS II. Commercially available hardware designed for the SAIL (Shipboard Ascii Instrumentation Loop) system was used for the current loop supply and for the RS-232/current loop converters. Information was transmitted over the loop as an RS-232C signal at 1200 baud. Appendix C. A Detailed Description of ULSAS II 135 C.5.1 Communications protocol. Communication between the surface station (call sign MASTER) and the buoy (call sign BUOY) is initiated by the surface station issuing a connect request: C BUOY Connection is then established by the two TNC :s, each of which will issue a message, i.e.: *** CONNECTED TO BUOY At this point both TNC's enter the converse mode, in which information received is "packetized" and transmitted so that the users at either end are unaware of the process. The main controller in ULSAS II will only respond to certain commands. All of them have several things in common: • all commands begin with an ampersand ("@") character; • all alphanumerics in command strings are in upper case; • after the ampersand, all command strings begin with a 3 letter mnemonic followed by a blank space; • all numeric parameters are sent as bytes in hexadecimal code, separated by commas; and • all command strings end with a carriage return character. Most commands are used to modify parameters in the main controller's memory. These "MODIFY" commands usually have a corresponding "RETURN" command, which will return the current parameters in memory. Valid command syntax. In the following descriptions, blank spaces and commas are necessary for correct syntax, and the ">" character is used to indicate a carriage return character. Modify and return commands Ping parameters QMDP A' sL0, NgEI, NCL0, NCZI,A^LO, 7VpHI, Nd> Appendix C. A Detailed Description of ULSAS II 136 This command will modify the ping parameters to be used in default collections — these parameters will be stored in EEPROM in the main controller. Since Ng, A r c , and Np are words, they must be passed as two bytes, LO and HI. There is another command: OMCP NgLO,NgEI,NcLO,NcEI,NpLQ,NpEI,Nd> which will modify the current data collection parameters, not the default parameters. Successful execution of either of these commands will return the message: MDP ***CMD OK*** In addition, both commands have their "return parameters" version: flRDP > SMCP > which will return the parameters currently in memory. For example the response to the first of these would be: MDP NgLQ, NgEI, NcL0 ,NCEI, NPLO, NPEI, Nd> Note that the ampersand is not returned. Sample parameters. QMDS Rc,CCLO,CCVLED,CCEI,Sj , M u x l , M u x 2 , . . . ,Muxl5> There are several things to note about this command. First, the parameter Cc is three bytes long, since it can be a very large number. Second, the sample frequency (Sf) byte has a small range of valid settings: Sf VALUE SAMPLING AT 07 1 per cycle 06 2 per cycle 05 4 per cycle 04 8 per cycle 03 16 per cycle 02 32 per cycle 01 64 per cycle Table C.18: Sampling frequency vs Sf setting. Third, all 15 available MUX address values must be included in the list. However, not all the addresses need be valid. The last address must have the most significant bit set to indicate the end of the list. Appendix C. A Detailed Description of ULSAS II 137 As an example, suppose that we wanted to sample only transducers 4, 1, 9, and 3 (in that order). Then the command would look like (using XX for all bytes other than the MUX addresses: fiMDS XX,XX,XX,XX,XX,04,01,09,83,00,00,00,00,00,00,00,00,00,00,00 This command also has the "modify current" version: 0MCS Rc ,CcL0,6'eMED,CCHI,Sj ,Muxl,Mux2,.. . ,Muxl5> and the "return parameters" versions: SRDS > ARCS > which will return a string like: MDS i?c,CcL0,CcMED,CcHI,5/,Muxl,Hux2, . . . ,MuxlS> Gain values. The receiver amp gains and the power amp gain can be modified using this command, which has the form: QMDG RA,PAL0,PAHI> where RA is the receiver amp gain, which can be set as follows: RA VALUE GAIN (dB) 00 0 01 12.04 02 24.12 03 36.12 Table C.19: Amplifier gain versus RA setting. PALO and PAHI are the least significant nibble and most significant nibble of the ping signal attenuator (see Table C.16). For example, if an attenuator value of C4 is required, PALO and PAHI would be set as 04 and 0C respectively. This command also has its "modify current" version: flMCG RA,PAL0,PAHI> and its return parameters versions: ORDG > fflRCG > which will return a string like: MDG RA,PAL0,PAHI> Appendix C. A Detailed Description of ULSAS II 138 Enables The receivers and the acoustic projector can be enabled/disabled for data collection. The form of this command is: OHCE Xmtr.Rcvr > where the value of Rcvr,Xmtr will be either 00 (disable) or 01 (enable). There is no default version of this command, because the default mode of both enables is forced on. There is the corresponding return parameters version of this command: QRCE > which returns: MCE Rcvr,Xmtr> Default collect time. The time between default collections can be modified using the command: fiMDC Time> where Time is the time in hours between default collections. This parameter can also be returned using: ORDC > Real time clock commands There are two commands: SSRT Min,Hour,DayMon,Month,Year> This command sets the real-time clock, and is self-explanatory. The other command is the return real time command: QXRT > which returns a string of the form: XRT Min,MinA.Hr.HrA,Day,DayMon,Month,Year> where the appended 'A' indicates the current alarm time, and Day is day of the week (Sunday = 1). Parameter return commands ADC parameters The command: fflADC > will cause the ADC on the 68HC11 to digitize all 8 channels and return the values in a string of the form: ADC ChO,Ch1,Ch2,Ch3,Ch4,Ch5,Ch6,Ch7> Appendix C. A Detailed Description of ULSAS II 139 I/O registers The command: ORIO > will return a string of the form: RIO 01,02,03,04,11,I2> where 01 through 04 are the contents of output registers 1 through 4, and II and 12 are the contents of input registers 1 and 2. Run data collection commands. Once the ping and sample parameters and gains are set prop- erly, data collections can be run using any one of the following commands: QG0R > This command will run a data collection using the current parameters, with the VCR on. fiGOB > As above, but with the recorder off. This mode might be used to check for saturation of the receivers using the peak detector output. QGDR > Here the recorder is on, and default rather than current parameters are used. 0GDB > As above, but with the recorder off. C.6 Local Microcontrollers. C.6.1 Main controller. Local control of all functions in ULSAS II is carried out by a Motorola 68HC11A1 microcontroller. This chip has a number of useful features which make it ideally suitable for this task: 512 bytes of EE-PROM for semi-permanent storage of data; a serial interface; an 8-bit 8-channel ADC/MUX; and on-board RAM, timer, and real-time interrupt. There are many tasks carried out by this device. First, two-way communication with the surface buoy (and thus with the user at a surface station) is through the serial interface. The communication syntax is shown in Section C.5.1 above. The transmit line (data from the 68HC11 to the surface buoy) also goes to a modem for conversion to an frequency shift key (FSK) signal which can be stored on either of the HI-FI channels of the VCR. In this way, if Appendix C. A Detailed Description of ULSAS II 140 CHANNEL PARAMETER 0 Peak detector 1 X-axis tilt, 2 Y-axis tilt, 3 Compass 4 Temperature 5 Main 12V battery voltage 6 Analog +12V battery voltage 7 Pressure Table C.20: 68HC11 MUX channel assignments. the VCR is running in the RECORD mode, a record is kept of all the responses of the main controller fo user requests. Second, because the system may sometimes wait for a long time, power conservation is necessary. Most power supplies are controlled by the 68HC11, and are off when the system is in idle mode. This issue is discussed further in Section C.7. Third, the 68HC11 handles time management. Accurate time is kept by a separate real-time clock with a self-contained power supply. This clock can also produce alarm signals at preset, times. When the system is in idle mode, that is, not communicating with the user or not carrying out data collections, the default, mode of operation of the clock, as programmed by the 68HC11, is to provide an alarm at, the beginning of each hour. The controller responds to this alarm by "waking up"and listening for communication from the surface buoy. If no valid command is received during a 5 minute period after the hour the system returns to the idle mode. Additionally, every time a valid command is received, the clock is reset, to provide an alarm after 20 minutes, so that, unless a data collection is underway, the system will power down after 20 minutes if no more commands are received. Another function for which the clock is used is the default data collection. In the event that com- munication between the surface station and ULSAS II is lost due to radio failure, cable breakage, etc., the system will collect data automatically. If no communication takes place for a user-specified number of hours, the controller will initiate a data collection, using the default parameters stored in EEPROM. Fourth, the on-board 8 channel ADC is used to return information about important environmental and system parameters to the user. For example, during deployment it is necessary to know the angles of tilt of the instrument. One of the commands to the controller causes it to sample all 8 channels of the ADC and return their values. The channel assignments are given in Table C.18. The peak detector is connected to the output of hydrophone 5 , and is reset at the beginning of each Appendix C. A Detailed Description of ULSAS II 141 ping. These data can be used to provide the user with some feedback as to the choice of gain settings for the receiver amplifiers and the projector output. The rest of the parameters are self-explanatory. Finally, the mechanics of data collection are orchestrated by the main controller. The ping and sample controllers (see Section C.6.2) must be loaded with the data collection parameters; gains must be set, and devices enabled or disabled; power supplies are turned on at appropriate times; and the VCR must be controlled. During the data collection, data are output to the user indicating the status of the collection. C.6.2 Ping/Sample controllers. These are two 8749 microcontrollers which operate asynchronously to both the 68HC11A1 and the master clock of the PCM system. One controls ping production, i.e. the provision of the correct driving signal to the PA/projector. The other controls sampling of data, i.e. the generation of MUX addresses at the appropriate times, and the production of the correct digital "tag" for the Right channel PCM information. The only connection between the two is an output from the ping controller indicating the correct time to sample, according to the. sample frequency (Sj) setting. These controllers are loaded with data and set to free-run by the 68I1C11. Each will send a message to the 68HC11 when they are finished their assigned task. C.7 Power Supplies. All power in ULSAS II comes from deep cycle gel-type lead-acid batteries. There is ample power to run all systems in data collection mode for at least 8 hours, and in idle mode for about 4 weeks. The battery supplies are divided into groups. 1. Main 12V battery — supplies +5V to digital circuits through a DC-DC converter, and also supplies 12V power to the VCR. This supply is always on, although the +5V supply to the PCM and the ping/sample controller board is only on during data collection. 2. Analog power supplies — these are + 12V and -12V batteries which supply power to all the analog circuitry, including the PCM. This supply is usually off, and is switched on during data acquisition. 3. Relay control battery — this is a 12V battery which is always connected. It is used to engage battery on/off relays under the control of the 68HC11. Appendix C. A Detailed Description of ULSAS 11 142 4. Power amplifier supplies — this is a set of 5 different supplies which are only switched on if the power amplifier is in use, i.e. if the acoustic projector is producing pings. The supplies include: ± 65VDC, ± 24 VDC, and ± 6VDC. C.8 Environmental Sensors. There were 5 sensors used to monitor the location and movement of the hydrophone array frame. These sensors were mounted inside a 6 inch aluminium pressure case, which was mounted on the bottom of the hydrophone array frame (see Figure C.49), and were connected to the main instrument housing by a multiconductor underwater cable. C.8.1 Azimuth. The compass was a solid state fluxgate sensor manufactured by KVH Industries, Middletown, Rhode Island. The accuracy of this compass is specified to be ±1° , over tilts of up to 45° from the vertical. The output of this sensor is an analog voltage of from 0 to 3.59 V, where 10 mV equals one degree. The response time is resistor programmable; I chose the fastest response time available, which was 2.7 seconds. C.8.2 Pressure. The pressure sensor was a solid state sensor made by SENSYM Inc. of Sunnyvale, California. It measured gauge pressure from 0 to 100 psig (recall that 1 psig = 0.6895 dbar). The accuracy was specified to be ± 0.3%FS, or ± 0.2 dbar, and repeatability was specified at ± 0.14 dbar. Output was an analog voltage from 1 V (corresponding to 0 psig) to 6 V (corresponding to 100 psig), which was "bucked off" using an instrumentation amplifier to a 0V to +5V range. C.8.3 Acceleration. The accelerometer used was a General Oceanics Model 6020, having a range of ± 1 g, with specified accuracy of ± 0.1%, or ± 0.002 g. The frequency response of this device was D.C. to 5 Hz, and the output was an analog voltage in the range ± 2.5 V. Appendix C. A Detailed Description of ULSAS II 143 C.8.4 Vertical tilt. Tilt in the A* and Y directions about the vertical axis were measured using plumb-bob type tiltmeters; Model CP-17-0601-1, manufactured by Humphrey. Inc., San Diego, California. The resolution of these devices is specified to be 0.2° over a range of ±45°. Output was an analog voltage in the range 0 to +5V, with 2.5V representing 0° tilt. The orientation of the A" and Y axes was defined according to the right-hand rule, with the vertical axis pointing away from the earth; the V axis was aligned with the compass reference direction. Appendix D Pre-SWAPP calibrations of ULSAS II D. l A D C System Calibration. In order to calibrate the digitization system, end to end, I injected known signals into the hydrophone input for channel 5. Since the digitization is carried out by the PCM electronics for all inputs, it really doesn't matter which channel is chosen. Signals generated by a Wavetec programmable function generator were passed to the programmable amplifier in channel 5, and then through the multiplexer to the PCM hardware. The PCM video output was stored on videotapes for later analysis to provide calibration factors. The gains of the programmable amplifiers had been measured previously, but as a check, signals were input and calibration files were collected at all amplifier settings. The inputs used for the different gain settings are given in Table D.21. The input signals were 1 Hz square waves with variable amplitude. Thus we could measure the calibration on both positive and negative inputs. D . l . l Results The recorded data for square wave input were collected on computer and analysed. Data from the positive side and negative side of the square wave pulse were separated, and the maximum, minimum, standard deviation, and mean were calculated for these two sets for each test. The results are given in Table D.22. These results can be used to compute a calibration equation in the form: Vi = (Di - Og) * Gg (D.138) GAIN Input Voltage 1 4.0V 4 1.0V 16 0.25V 64 0.0625V Table D.21: Calibration voltage input for different gain settings. 144 Appendix D. Pre-SWAPP calibrations of ULSAS II 145 GAIN 11/L Min Max Mean Std Dev N 1 H 25565 25692 25566.3 5.9 450 1 L -24950 -24928 -24939.7 4.2 500 4 H 25439 25439 25439.0 0.0 450 4 L -24859 -24819 -24837.3 6.8 700 16 H 25439 25439 25439.0 0.0 400 16 L -24927 -24759 -24804.4 15.1 600 64 H 25057 25184 25061.6 23.7 500 64 L -25537 -24503 -24660.9 74.4 400 Table D.22: Calibration measurements for different gain settings (H : high level, L = low level input). GAIN Offset(0,) Factor(G,) 1 313.3 158.40E-6 4 300.8 39.78E-6 16 317.3 9.9516E-6 64 200.3 2.514E-6 Table D.23: Calibration offsets and factors for the 4 available instrumentation amplifier gains. where V; is the voltage out of the hydrophone/preamp system for the i'h hydrophone, and D{ is the corresponding digital output from the ADC. Os is the digital correction for the zero offset, which is due almost totally to the PCM input electronics, since all other components of the data acquisition system have been carefully adjusted to eliminate zero offsets. The calibration factors(Gj) for the 4 available gain settings are given in Table D.23. The subscript 'g' refers to the gain setting of the amplifiers. These calibrations can be used to relate any digital output to the input voltage, when those voltages pass through the instrumentation amplifiers. Voltages from devices which do not pass through the instrumentation amplifiers (such as the environmental sensors) can use the calibration factors for Gain of 1. D.2 Calculation of source level of Argotech projector. Because one of the main purposes of this work was to measure the surface scattering strength parameter Ss, it was necessary to measure the projector source level at the hydrophones directly. In order to do so. several data collections were made with the post-amplifier gain set to 0 dB, and with different settings of projector attenuation. The hydrophones and the projector are separated by 2 m, and according to proximity criteria for acoustic measurements given in Bobber [7], at this separation the hydrophone array is in the far field relative to the projector. Since we were measuring at a rate of 4 times per cycle, 4/o, and the separation between the receivers Appendix D. Pre-SWAPP calibrations of ULSAS II 146 5000 •3 IOOO eo Jj -1000 > -3000 3000 -5000 0 5 10 15 20 25 30 35 Figure D.55: Example of data used in non-linear fit to sinusoid, together with fitted curve. and projector is slightly different for each hydrophone, we cannot measure the ping amplitude directly. However, it is possible to fit the measured data to the expected form of the transmitted pulse, and estimate desired parameters using a nonlinear fitting technique. Each ping had 7 cycles, but 1 removed the first and last cycle of data to eliminate turn-on and turn-off transients from the analysis, so there were 20 data points to use for the fit, for each hydrophone, for each experiment. ] used the SYSTAT nonlinear fit package to fit the data to the following model: where the Xn are the data points, and the mean A, amplitude B, and phase C are the parameters to be estimated from the nonlinear fit. In order to get useful estimates for the phase, I had to modify the loss function from its usual squared error form, as follows: LOSS = (C > 0 AND C < 2TT) X (Xn - ESTIMATE)2 + (C < 0 OR C > 2ir) x 1025 (D.140) This loss function constrained the estimated phase to be in the range 0 to 27r. SYSTAT allows two different minimization algorithms to be used: the Quasi-Newton method, or the Simplex method. I tried both methods on several data sets, and found no differences between the estimates of the parameters. Four data sets were analyzed, each with a different value of the power amp input attenuation. The post-amplifier gain was 0 dB in all cases. The results are given in Tables D.24 through D.27. "ASE" (D.139) Appendix D. Pre-SWAPP calibrations of ULSAS II 147 150 100 cs 5 50 "3 s •e •S3 0 u OS -50 -100 0 5 10 15 20 25 30 35 Figure D.56: Residuals after non-linear fit to sinusoid. Hydrophone Amplitude ASE Phase ASE Mean ASE 1 3769 21 32.5 0.28 -82 15 2 9269 55 28.9 0.28 -95 42 3 9443 58 28.8 0.34 -128 41 4 8964 69 25.7 0.46 -100 49 5 9604 58 27.8 0.34 30 40 6 8931 67 23.2 0.40 -109 47 7 8462 60 25.7 0.40 -97 42 8 8887 61 29.0 0.40 -41 43 9 8590 56 20.1 0.34 -112 39 Table D.24: Parameter estimates for power amp input attenuation setting of 64. in these tables refers to the Asymptotic Standard Error as estimated by SYSTAT. Figure D.55 shows some sample data together with the fitted curve, and Figure D.56 shows the residuals after fitting. In the worst case, only 5% of the data value is "unexplained" by the fitting procedure. I then computed a linear relation between the computed amplitudes of the fitted sinusoid and the 4 values of transmitter attenuation for all nine hydrophones: Df = KfxA + Rf (D.141) where 'J4' is the attenuation setting of the transmitting system, and Df is the estimated digital amplitude for the ith hydrophone. These relationships are summarized in Table D.28. These calibration runs were all carried out at. an amplifier gain of 0 dB, so that extension of the calibrations for higher amplifications has to take that into account (see equation D.143). Note that in all nine cases the correlation coefficient Appendix D. Pre-SWAPP calibrations of ULSAS II US Hydrophone Amplitude ASE Phase ASE Mean ASE 1 1880 44 32.5 1.32 - 2 4646 47 29.0 0.57 - .3 4741 47 28.8 0.57 - 4 4488 48 25.9 0.63 - 5 4823 28 27.9 0.34 -• 6 4485 45 23.2 0.57 - 7 4233 46 25.5 0.63 -76 15 8 4439 22 29.3 0.29 - 9 4302 45 20.0 0.57 - D.25: Parameter estimates for power amp input attenuation setting Hydrophone Amplitude ASE Phase ASE Mean ASE 1 966 9.1 31.5 0.52 -115 6.5 2 2372 4.8 28.8 0.11 -89 3.4 3 2419 4.2 28.3 0.11 -101 2.9 4 2300 8.9 26.0 0.23 -81 6.3 5 2442 6.4 27.7 0.17 -70 4.5 6 2296 7.0 23.3 0.17 -82 5.0 7 2159 7.2 25.8 0.17 -79 5.1 8 2267 4.6 28.7 0.11 -41 3.3 9 2197 4.6 19.7 0.11 -78 3.0 Table D.26: Parameter estimates for power amp input attenuation setting of 16. was 0.9999 or better. Results of this calibration procedure indicate that hydrophone 1 was malfunctioning in the sense that, the gain was much reduced. Further analysis of the output from this hydrophone also revealed anomalous frequency response, so that its output could not be used in later analysis. Because much of the analysis of SWAPP data involved computation of backscattering coefficients, an expression relating the ratio of the source level (SL in dB) and the receive level (RL in dB) is required. Given the calibration results, the source level can be given by the following equation: SL = 20 log VtA - 20 log(V2) - RR, where V)A is the output from the ith hydrophone (not including the preamp gain) at, a transmitter attenuation setting 'A', and RRj is the open circuit, voltage response of the hydrophone in question. The factor 201og(\/2) is necessary because the calibration described above is for the amplitude of the sinusoid, and we must divide the estimated voltage by \/2 to reduce to rms voltage. Now the hydrophone output, voltage can be given as a function of the preamp gain P,, the amplifier gain setting '</', and the transmitter attenuation setting 'A''. Combining equations D.138 and D.141, and ignoring offsets, since we remove the means of any time series used in these calculations, we get: Appendix D. Pre-SWAPP calibrations of ULSAS II 149 Hydrophone Amplitude ASE Phase ASE Mean ASE 1 477 4.3 34.1 0.52 -52 3.0 2 1194 3.0 28.8 0.17 -33 2.2 3 1217 2.7 28.7 0.11 -25 1.9 4 1158 3.2 26.0 0.17 -48 2.3 5 1241 2.3 27.9 0.11 -9 1.6 6 1152 3.5 23.0 0.17 -28 2.5 7 1087 2.2 25.1 0.11 -43 1.6 8 1149 3.2 29.3 0.17 -25 2.3 9 1105 2.9 19.6 0.17 -18 2.0 Table D.27: Parameter estimates for power amp input attenuation setting of 8. Hydrophone Factor (Kf) Offset (Rf) 1 58.6 14.0 2 144.0 50.2 3 146.7 54.4 4 139.2 52.3 5 149.3 48.8 6 138.7 56.2 7 131.5 39.0 8 138.1 43.9 9 133.5 44.2 Table D.28: Parameter estimates for power amp input attenuation of 08H. A _ Kf xAxGi Pi Finally, we can write an equation for the source level as: SL = 20log(A7* x A) - 20log P ; - 79.02 - RR;. (D.142) The constant factor 79.02 includes both the 201og\/2 term and the 201ogG'i term. Similarly, we can write an equation for the receive level (RL): RL = 20 log Vt - RR, where V{ is the raw voltage output of the hydrophone (before pre-amp). Again using equation D.138 we can rewrite this as: RL = 20 log Di + 20 log Gg - 20 log P{ - RR, (D.143) Appendix D. Pre-SWAPP calibrations of ULSAS II 150 where D, is the digital output of the i,h hydrophone. Determinations of backscattering strength will involve the computation of the sum 'RL - SL', which is proportional to the logarithm of the ratio 'RL/SL'. Combining equations D.142 and D.143 gives: RL - SL = 20 log Di + 20 log Gg - 20 \og(KA x A) + 79.02. (D.144) Note that this relation eliminates the need to include variables such as preamp gain and open circuit voltage response in the computation. As a rough end-to-end test of the transmitter calibrations. 1 used equation D.142, measured preamp gains, published data on the OCV of the lTC-1001 hydrophones, and the data from Table D.28 to estimate the actual source level of the Argotech projector. Using values of A — 64, RR = -187 dB, Pi — 23, we get SL = 160 dB. Now the nominal specification for the transmitting voltage response of the Argotech projector at 400 Hz is 8 dB re 1 microbar/V at 1 m, or 108 dB re 1 micropascal/V. Our measured input voltage (peak-to-peak) at a gain of 64 was 460 V, or 325 Vrms. Thus the expected SL based on these numbers is 108 + 20log(325) = 158 dB. I estimate that the difference of 2 dB is well within the uncertainty due to the expected variability of the values used in this computation. D.3 Instrument problems during the SWAPP cruise In addition to the unusually calm weather, there were a number of problems with ULSAS during SWAPP. First, we noticed early on in the cruise that the hydrophones were responding to the pressure changes caused by the small vertical motions of the array. This characteristic was unexpected, as the measured response of the pre-amplifiers rolls off at 6 dB/octave below 100 Hz. In order to fix the problem, we had to insert capacitors in the cables between the hydrophones and the pre-amplifiers, since the pre-amps were potted onto the cables and were not accessible. It was not possible to totally remove the effect, but we were able to reduce the very low frequency (VLF) response of the hydrophone/pre-amps to tolerable levels. We also had many problems with the underwater connections - through some design flaw in the connectors, O-rings were being damaged and small leaks were occurring in many of the connectors. Since there were so many connectors used between the hydrophone array and the instrument housing 20 m below, this was a serious problem which was never completely eliminated. The most serious problem was that we could not operate the projector at any more than 1/4 of its maximum driving voltage (thus 1/16 of its maximum available output power) because of apparent Appendix D. Pre-SWAPP calibrations of ULSAS II 151 high voltage leakage somewhere in the path between the power amplifier and the projector. Thus the maximum SL attainable with the system as deployed on SWAPP was roughly 160 dB, or about 13 dB below the theoretically obtainable maximum with this projector at 400 Hz. The high voltage leakage problem had occurred during laboratory tests of the projector, and was traced to a faulty connector, which was subsequently replaced. However, the problem re-occurred during SWAPP, was never solved during the cruise, and proved to be a fundamental limitation in our ability to get useful data from the deployments of ULSAS during SWAPP. Finally, the environmental sensors were malfunctioning intermittently. They seemed to be turning off from time to time, and were actually non-functional for many data collections. This problem was never solved satisfactorily during SWAPP, but later in the cruise the problems seemed to occur less often. The net result of all these problems is that we were not getting useful data until almost the end of the cruise. Appendix E Details on data processing. E.l Preliminary Data Processing E.l.l Hardware The data processing hardware is shown in Figure E.57. There are two sources of data, which are discussed separately below: PCM data, which are measurements from the hydrophones and environmental sensors; and communications data, which are a record of the output of the local controller to the user during data acquisition sequences. P C M data As we have seen in Chapter 5, the PCM system has the capability of storing very large amounts of data; the storage rate is 176,400 bytes/s, so that an 8-hour VHS tape can store 5.08 Gigabytes of data. This quantity of data is difficult, to deal with, both in terms of the. rate at which it is read off the videotapes during playback, and the storage space required (for example on a hard disk) if a small computer is used to take the data in and store it for later analysis. Because of the data rate, a small computer would be kept very busy just reading the data in and storing it to disk. I have developed hardware and software to do that task on AT-class personal computers, and it requires the use of DMA (direct memory access) hardware, FIFO buffers as well as large buffers in RAM just to be able to get the data off the videotape and onto disk. Part of the problem with data rates was taken care of during design. When I calculated the data rates and storage requirements for operation at 400 Hz with 9 hydrophones and 5 environmental sensors, I realized that, even sampling at 4 times the carrier frequency, the data rate is only 96,000 bytes/s, or 54.4% of the capacity of the PCM system. The reader may recall from Appendix C that we were not allowed to use a MUX address of zero. Thus whenever the MUX address is zero, that is, whenever the least significant 4 bits of the right channel PCM data are all zero, we know that this word and the subsequent left channel word can be ignored. Thus when operating at 400 Hz with 14 channels of data, 152 Appendix E. Details on data processing. 153 VCR Video Hi/Fi MODEM circuit 286 PC PCM 601ESD (modified) Serial Clocks Data Reducer Serial + Clocks, RS-232^^"4^ Intfc. f T O A „ J ISA Bus DSP-56 Figure E.57: Data processing hardware components and interconnections. the PCM data stream consists of groups of 14 pairs of valid left/right data separated by roughly 13 pairs of invalid left/right data, with MUX addresses of zero. It is a relatively simple matter to take the serial data and clocks coming from the PCM and screen out the invalid data, thus reducing the data rate considerably. However, in spite of this reduction the data rates are still too high to do much other than store the data to disk, even using a fast personal computer and high speed hard disk. The quantity of data which has to be stored for a given sampling period has also been reduced, but at 400 Hz carrier with 14 channels sampling 4 times per carrier cycle, it still requires about 2.7 Mbyte of disk space to store one minute of data! Since in DOS systems the largest disk partition available is 32 Mbytes, the time period for data series is limited to roughly 12 minutes. Aside from this limitation, there remains the problem of what to do with the large amount of data once it is stored to disk: even the simplest processing on so large a data set will be very time consuming. The solution to the data acquisition and processing problem is to use one of the available single- chip digital signal processing (DSP) systems to carry out the data acquisition and initial processing in real-time, which accomplishes several things. Appendix E. Details on data processing. 154 • Initial error checking, for example checking that MUX addresses are in the proper sequence and that the data count is incrementing properly, must be the first stage of any processing of the real- time data. The DSP chip takes care of this task while reading in the data, immediately halving the quantity of data which must be saved or processed, since we don't need to save the data from one of the two PCM channels. • Preprocessing of the data, such as averaging, filtering, binning, etc. are all tasks for which DSP chips have been optimized. Many of these operations further reduce the quantity of data. • The DSP system and the host computer are separate entities, connected by an interface which is used to load programs and data (constants, parameters, etc.) from the host to the DSP, and to pass error-checked and preprocessed data from the DSP to the host. Thus the host is only required to load a much reduced amount of data from the DSP, and has plenty of time for storage, display and post-processing of the data it receives. 1 chose to use a fixed-point processor board based on the Motorola DSP56001, which is a second generation DSP whose basic data unit is a 24-bit integer interpreted as a signed fraction. This board (manufactured by Spectrum Signal Processing Inc. in Burnaby, B.C.) has 48 kwords of fast memory, two serial ports, one parallel port, and a host interface. One of the two serial ports is used to get data from the PCM via the data rate reducer. The DSP56001 has a clock frequency of 20.5 MHz, and most of its instructions execute in 2 clock cycles. At the maximum PCM data rate of 88.2 kHz, the DSP56001 can execute roughly 100 instructions between PCM words. Given the high degree of parallelism inherent in the chip design, there remains time to do plenty of data pre-processing. User interface data As discussed in Chapter 5, the ouput of the local controller to the user over the radio interface is also saved as FSK data on the Hi-Fi audio channels of the VCR. Therefore all the parameters of the data collection, including time, environmental sensor snapshots, and instrument setup are available for examination during processing of the data. This information also acts as a label for the PCM data stored on the video channels of the tape. To decode these data, the audio channel output was fed to a single-chip MODEM for FSK to TTL-level conversion, to level shifters for T T L to RS-232 conversion, and then to a serial interface on the personal computer. These data can be examined using any terminal emulation program, or can be read directly from the serial interface during processing of the PCM data. Appendix E. Details on data processing. 155 E.1.2 Software Software for processing of PCM data by the DSP56001 board was written in Motorola assembler on the PC. Spectrum Signal Processing provides a library of "C" functions and procedures for downloading and running programs on the DSP board, and for the transfer of data between the host PC and the DSP. Programs for storage, display and analysis of data on the PC were written primarily in "C", although some analysis was done using FORTRAN programs from Numerical Recipes [57], The commercial programs "SYSTAT" and "SYGRAPH" were also used for some statistical analysis and display of data. Bibliography D.E. Barrick. Remote sensing of sea state by RADAR. In V.E. Derr, editor. Remote Sensing of the Troposphere, chapter 12, U.S. Gov't Printing Office, 1972. D.E. Barrick. Technical report on the second-order acoustic doppler spectrum backscattered from ocean surface gravity waves. Technical Report, Ocean Surface Research, 1984. D.E. Barrick and B.J. Lipa. A compact transportable HF RADAR system for directional coastal wave field measurements. In M.D. Earle and A. Malahoff, editors, Ocean Wave Climate, pages 153- 201, Plenum Press, New York, 1977. D.E. Barrick and B.J. Lipa, Ocean surface features observed by HF coastal ground-wave RADARS: a progress review. In M.D. Earle and A. Malahoff, editors. Ocean Wave Climate, pages 129-152, Plenum Press, New York, 1977. D.E. Barrick and J.B. Snider. The statistics of HF sea-echo Doppler spectra. IEEE Transactions on Antennas and Propagation, AP-25(l):19-28, 1977. D.E. Barrick and B.L. Weber. On the nonlinear theory for gravity waves on the ocean's surface, part 2: interpretation and applications. Journal of Physical Oceanography, 7(1):11-21, 1977. R.J. Bobber. Underwater Eleciroacoustic Measurements. Peninsula Press, Los Altos, California, 1988. C. Bretschneider. Shore Protection Manual, Vol 1. Technical Report, Coastal Engineering Research Center, Fort Belvoir, Virginia, USA, 1973. P. Broche, P. Forget, J.C. deMaistre, J.L. Devenon, and M. Crochet. VHF RADAR for ocean surface current and sea state remote sensing. Radio Science, 22(l):69-75, 1987. W.S. Burdic. Underwater Acoustic System Analysis. Prentice-Hall, Englewood Cliffs, New Jersey, 1984. RP. Chapman and J.H. Harris. Surface backscattering strengths measured with explosive sound sources. Journal of the Acoustical Society of America, 34:1592, 1962. C. S. Clay and H. Medwin. Acoustical Oceanography. Wiley-Interscience, Toronto, 1977. D. D. Crombie. Doppler spectrum of sea echo at 13.56 mc/s. Nature, 175:681-682, 1955. S. Oberski D. Lemon, D. Knight and J. Marko. Development, of techniques for acoustical measure- ment of ocean wave spectra. Interim Report No. 2. Measurement System Design. Technical Report, Arctic Sciences Ltd., 1986 Mills Rd., Sidney, B.C., 1984. J. DeSanto and G.S. Brown. Analytical techniques for multiple scattering from rough surfaces. Progress in Optics, 23:, 1986. C. Eckart. The scattering of sound from the sea surface. Journal of the Acoustical Society of America, 25(3):566-570, 1953. 156 Bibliography 157 [17] L. Fortuin. Survey of literature on reflection and scattering of sound waves at the sea surface. Journal of the Acoustical Society of America, 47(5): 1209—1228, 1970. [18] G.H. Golub and C.F. VanLoan. Matrix Computations. The Johns Hopkins University Press, Bal- timore, Maryland, 1983. [19] The SWAMP Group. Ocean Wave Modelling. Plenum Press, New York; 1985. [20] F.J. Harris. On the use of windows for harmonic analysis with the discrete fourier transform. Proceedings of the IEEE, 66(l):51-83, 1978. [21] D.R. Jackson and D.P. Winebrenner. Comparison of perturbation theories for rough-surface scat- tering. Journal of the Acoustical Society of America, 83(3):961-969, 1988. [22] G.M. Jenkins and D.G. Watts. Spectral Analysis and its Applications. Holden-Day, Oakland, California, 1968. [23] B. Kinsman. Wind Waves. Their Generation and Propagation on the. Ocean Surface. Dover Publi- cations Inc., New York, 1965. [24] P.D. Koenigs. Amplitude and spectral characteristics of low grazing angle acoustic backscatter. Technical Report NUSC TD 7875. Naval Underwater Systems Center, 1987. [25] E.Y.T. Kuo. Sea surface scattering and propagation loss: review, update, and new predictions. IEEE Journal of Oceanic Engineering, 13(4):229-234, 1988. [26] F.M. Labianca, Estimation of ocean-surface directional-frequency spectra: the inverse scattering problem. Technical Report OSTP-67 FL, Bell Laboratories, 1978. [27] F.M. Labianca and E.Y. Harper. Sideband structure of sound from a harmonic point source scat- tered by a rough surface moving over an upward-refracting ocean. Journal of the Acoustical Society of America, 61(2):378-389, 1977. [28] P. Lanzano. Inversion algorithms for geophysical problems. Technical Report Memorandum Report 6138, Naval Research Laboratory, 1987. [29] P.H. LeBlond and D.A. Huntley. Ocean Wave Directionality: a review of current knowledge. Tech- nical Report, National Research Council of Canada, 1984. [30] P.H. LeBlond and L.A. Mysak. Waves in the Ocean. Elsevier, Oxford, 1978. [31] B.J. Lipa and D.E. Barrick. Analysis methods for narrow-beam high-frequency RADAR sea echo. Technical Report. ERL- 420-WPL56, NOAA, 1982. [32] B.J. Lipa and D.E. Barrick. Extraction of sea state from HF RADAR sea echo: mathematical theory and modelling. Radio Science, 21(1 ):81—100, 1986. [33] B.J. Lipa and D.E. Barrick. Least-squares methods for the extraction of surface currents from CODAR crossed-loop data: application at arsloe. IEEE Journal of Oceanic Engineering, OE- 8(4):226-253, 1983. [34] B.J. Lipa, D.E. Barrick, and J.W. Maresca. HF RADAR measurements of long ocean waves. Journal of Geophysical Research, 86(C5):4089-4102, 1981. [35] R.B. Long. The statistical evaluation of directional spectrum estimates derived from pitch/roll buoy data. Journal of Physical Oceanography, 10:944-952, 1980. Bibliography 158 [36] R.B. Long and K. Hasselmann. A variational technique for extracting directional spectra from multi-component wave data. Journal of Physical Oceanography. 9:373-381, 1979. [37] M.S. Longuet-Higgins and D.E. Cartwright. Observations of the directional spectrum of sea waves using the motions of a floating buoy. In Ocean Wave Spectra, pages 111-136, Prentice-Hall, 1963. [38] J.F. Lynch, R.C. Spindel, C. Chiu, J.H. Miller, and T.G. Birdsall. Results from the 1984 marginal ice zone experiment preliminary tomography transmissions: implications for marginal ice zone, arctic, and surface wave tomography. Journal of Geophysical Research, 92(C7):6869-6885, 1987. [39] J.W. Maresca and T .M. Georges. Measuring rms waveheight and the scalar ocean wave spectrum with HF skywave RADAR. Journal of Geophysical Research, 85(C5):2759-2771, 1980. [40] R.H. Mellen. On underwater sound scattering by surface waves. IEEE Journal of Oceanic Engi- neering, 14(3):245-247, 1989. [41] M. Nieto-Vesperinas and N. Garcia, A detailed study of the scattering of scalar waves from random rough surfaces. Opt. Acta., 28:1651-1672, 1981. [42] D.W. Oldenburg. An introduction to linear inverse theory. IEEE Transactions on Geoscience and Remote Sensing, GE-22(6):665-674, 1984. [43] D.L. Phillips. A technique for the numerical solution of certain integral equations of the first kind. SIAM J. Assoc. Comp. Mach., 9:84-97, 1962. [44] R. Pinkel and J.A. Smith. Open ocean surface wave measurement using Doppler sonar. Journal of Geophysical Research, 92(C12): 12967-12973, 1987. [45] J.W. Strutt(Lord Rayleigh). The Theory of Sound. Dover Press, New York, 1945. [46] S.O. Rice. Reflection of electromagnetic waves from slightly rough surfaces. Commun. Pure Appl. Math., 4:351-378, 1951. [47] W.I. Roderick and B.F. Cron. Frequency spectra of forward-scattered sound from the ocean surface. Journal of the Acoustical Society of America, 48(3):759-766, 1970. [48] J.A. Scrimger. Analysis of sea surface conditions by surface scattered sound. Technical Report, Jasco Research Ltd., 1983. [49] R.J. Seymour. Estimating wave generation on restricted fetches. Journal of Waterway Port Coastal and Ocean Division, 103 (WW2):251-263, 1977. [50] E.I. Thorsos. Acoustic scattering from a "pierson-moskowitz" sea surface. Journal of the Acoustical Society of America, 88(l):335-349, 1990. [51] E.I. Thorsos and D.R. Jackson. The validity of the perturbation approximation for rough surface scattering. 1987. Presented at 113th meeting of the Acoustical Society of America. [52] S. Twomey. On the numerical solution of fredholm integral equations of the first kind by the inversion of the linear system produced by quadrature. SIAM J. Assoc. Comp. Mach., 10:97-101, 1963. [53] R. J. Urick. Principles of Underwater Sound. McGraw-Hill, New York, 1983. [54] C. van Schooneveld. Inverse problems: a tutorial survey. In Y.T . Chan, editor, Underwater Acoustic Data Processing, NATO ASI Series E, Kluwer Academic Publishers, Dordrecht, 1989. Bibliography 159 '[55] J .M. Varan. Pitfalls in the numerical solution of linear ill-posed problems. SIAM J. Sci. Stat. Comput., 4(2):164-176, 1983. [56] J .M. Varah. A practical examination of some numerical methods for linear discrete ill-posed prob- lems. SIAM Review, 21 (1): 100-111, 1979. [57] S. Teukolsky W. Press, B. Flannery and W. Vetterling. Numerical Recipes: The Art of Scientific Computing. Cambridge University Press, New York, 1986. [58] B.L. Weber and D.E. Barrick. On the nonlinear theory for gravity waves on the ocean's surface, part i: derivations. Journal of Physical Oceanography, 7(1):3—10, 1977. [59] R.G. Williams. Estimating ocean wind wave spectra by means of underwater sound. Journal of the Acoustical Society of America, 53(3):910—920, 1973. [60] L. Wyatt. The measurement of the ocean wave directional spectrum from HF RADAR Doppler spectra. Radio Science, 21(3):473-485, 1986.
Cite
Citation Scheme:
Usage Statistics
Country | Views | Downloads |
---|---|---|
China | 13 | 1 |
United States | 9 | 1 |
Canada | 3 | 0 |
Iran | 3 | 0 |
France | 2 | 0 |
United Kingdom | 1 | 0 |
City | Views | Downloads |
---|---|---|
Beijing | 9 | 0 |
Unknown | 5 | 20 |
Ashburn | 4 | 0 |
Mountain View | 3 | 1 |
Brooks | 3 | 0 |
Guangzhou | 2 | 0 |
Hangzhou | 2 | 0 |
Weston-super-Mare | 1 | 0 |
Sunnyvale | 1 | 0 |
Redmond | 1 | 0 |
{[{ mDataHeader[type] }]} | {[{ month[type] }]} | {[{ tData[type] }]} |
Share
Share to: