I. Introduction
The technique of power spectrum analysis has for a number of years been a useful tool for communications engineers as well as other scientists (Reference Blackman and TukeyMarks and Pierson, 1952; Reference PanofskyPanofsky, 1955) in extracting periodic variations from space or time series. In this paper the power spectrum technique is used to analyze remote sensing data of the Arctic ice pack, in the form of undersea and surface profiles, collected by the U.S. Navy. Since a visual examination of both under-ice and surface profiles suggest that the profiles are wave-like in nature, the analysis of a random sample using power spectrum techniques should reveal certain spatial periodicities. Intuitively such preferred spacings might be expected, inasmuch as if the ice behaves in an essentially elastic manner then the spacings between fractures and ridges will to some degree be controlled by the thickness and physical properties of the ice.
In addition to indicating the dominant periodicities for several samples of pack-ice profiles, power spectra of intersecting profiles suggest that the ridges do not exhibit directional isotropy. Moreover, the lineation or directionality of the pressure ridges is approximately determined by the analysis of two intersecting profiles. Also, using surface profiles, we find that the spectral plots may be quite useful in determining ice types by the amplitude of spectral components in different frequency regions.
II. Data reduction and limitations
The data used for this analysis were obtained by the U.S. Navy and consist of (1) upward looking sonar profiles taken in the Beaufort and East Chukchi Seas in February 1960, and (2) airborne laser profiles taken in April 1968 approximately 300 km north of Point Barrow.
A. Surface profile data
The U.S. Naval Oceanographic Office on 10 April 1968 sponsored several test flights at different altitudes over arctic pack ice with a “Spectra Physics Geodolite 3A” laser profiler. The profiler is a continuous wave source of coherent light of 632.8 nm which transmits an amplitude-modulated laser beam. The phase delay between the transmitted and reflected light at the modulated frequency is measured, providing a distance measurement as accurate as 3 cm (depending on the modulation time) at 300 m, the altitude at which the profiler was flown. The beam has an angular divergence of only 10−4 rad so that the cone size is effectively negligible. The profile was recorded in the form of analog strip charts.
The most significant source of error in the profile data is the variation of the aircraft altitude. This variation however is of a relatively low frequency and can be approximately removed by high-pass filtering. The profile height variation due to yaw, roll and pitch of the airplane are also typically of a relatively low frequency and can thus to a large extent also be removed by high-pass filtering. Errors due to variation in the laser beam’s speed over the ground will be present, but since the profiles examined consist of lengths of only about 500 m each we expect these variations to be small compared to the mean aircraft speed (about 93 ms−1).
B. Sonar profile data
The sonar data of interest in this report were obtained by the U.S. Naval Undersea Research and Development Center at San Diego, California. Sonar transducers were mounted on submarines that cruised beneath the arctic pack ice and a sound pulse was beamed vertically upward. The returning echoes were recorded on “Edo” sounders, similar to those used by ships for depth recording. These sounders record the profile data on a specially treated paper strip chart driven at constant speed beneath an electrostatic stylus. The stylus moves in such a manner that the outgoing sonar pulse is synchronized with the start of the stylus across the face of the paper. When a returning pulse is received, the stylus is activated and sparks, burning a mark on the strip-chart. In this fashion, an analog profile of the under-ice surface is drawn, with much vertical exaggeration, as the submarine cruises beneath the ice pack. These strip-chart profiles were the major source of data for this study.
The limitations of the sonar data stem primarily from the large cone diameter of the sonar beam. In particular at depths of 45 m and 105 m, from which data for this study were taken, the sonar beam “illuminated” areas of 70 m2 and 390 m2 respectively. Although this indicates poor resolution, the data suggest that in practice, the records consist primarily of vertical downward reflections. This is true because although some energy travelling upward at an angle can be reflected from a surface at a distance away from the zenith, and, with the right combination of angle of reflectance and refraction, can be returned to the transducer, this does not occur in a continuous fashion and examination of the records shows that such reflections have an attenuated and anomalous appearance. However, even with a careful examination of the records to eliminate anomalous reflections, there may still be considerable deformation of the actual profile due to the cone diameter. Consequently, one should not rely on the profiles to determine absolute thickness nor to produce realistic spectra at high frequencies. However, in the calculation of the spectrum, one is not concerned with absolute thickness but only with the periodicities. Also the spectral peaks of interest occur at lower frequencies which will be representative of the ice even if the cone width deforms the profile somewhat.
There are errors, however, that will affect the spectrum analyses. The first and most important factor is the existence of small changes in the direction of movement of the submarine (in order to avoid deep-draft ridges) which are quite common when the submarine is cruising in shallow water. We have considered in this paper only tracks where either the submarine was at a relatively great depth and did not undergo a series of small course corrections, or tracks where there were no deep ridges to avoid. The other factor is the variation of the ship speed which is typically of the order of ±0.5 m s−1. However, since this is small compared to the overall velocity (approximately 5 to 5.5 m s−1) it is effectively negligible.
C. Digitizing procedure
Since all the submarine sonar and airborne laser data were recorded in analog form, digitization was necessary before numerical analysis by digital computer could be initiated. The process of digitizing was accomplished by means of a semi-automatic line-follower digitizer that operated in the following manner: A time-base generator programmed a mechanical arm to move across a strip-chart segment at a constant speed along its x-direction. A cross-hair, mounted on the arm and moveable in the y-direction, was controlled manually using a potentiometer. The voltage across the potentiometer reproduced, in analog form, the original data. This analog wave-form was then processed by an analog/digital converter and the output punched on paper tape. There are several errors inherent in this system, but they are believed to be of the same order of magnitude or less than those inherent in the data recording discussed above.
III. Power spectrum analysis technique
A. Theory
In this section we briefly describe the theory of power spectra. For simplicity we will consider an infinite, continuous time series. In practice, of course, one has a finite digitized sample and a finite data window whose transform is convoluted with the desired “true” power spectrum. The techniques of using appropriate filters as is done in our computer programs to extract the “true” power spectra is beyond the scope of this paper, and is discussed at length in Reference Blackman and TukeyBlackman and Tukey (1958).
Given a time series X(t), the time average,
where
is called the auto-covariance function, with τ being called the lag time.
The quantity
is called the power spectrum. Power-spectrum estimates for given time series are usually plotted in terms of relative amplitude per unit frequency versus incremental frequency band- widths. The frequency units are denoted by wave-number with the lighest wave-number corresponding to the Nyquist frequency 1/(2ΔT) where ΔT is the sampling interval. In this work, the basic units are distance, not time, therefore the frequencies discussed are spatial frequencies; the reciprocal of spatial frequency is the wavelength, which is an easier unit to visualize and therefore is often given as an additional scale on power-spectra plots.
From the above definition of the power spectrum it is easily shown that C(f) is also given by
This expression is most useful for interpretation because the right-hand side of Equation (5) is simply the square of the Fourier transform of η(t). This procedure is, in fact, simply one way of carrying out a Fourier analysis. Since the auto-covariance function consists of an average over many data points, its spectrum more nearly represents the “true” spectrum (Reference Blackman and TukeyBlackman and Tukey, 1958) of the time series. Calculating the auto-covariance function is also a way of smoothing the data. From the general theory of power spectra one can predict, from the number of data points and the number of lags, the confidence limits (Reference Blackman and TukeyBlackman and Tukey, 1958; Reference Jenkins and WattsJenkins and Watts, 1968) of a given spectrum about the “ideal” spectrum assuming a stationary Gaussian process. These confidence limits are quite useful in determining the significance of a given spectral peak. The confidence limits may be determined from the equivalent degrees of freedom per spectral point. In the calculations presented in this paper a Hamming spectral window is used so that the equivalent number of degrees of freedom per spectral point (Reference Jenkins and WattsJenkins and Watts, 1968) is 2.67 T/M where T is the length of sample analyzed and M is the maximum value of lag time.
It should also be pointed out that in the computational program used in this paper, the auto-covariance function is computed by the formula
where ΔT is the sampling interval so that, for a given stationary time series, ρ(τ) will be effectively independent of the number of data points used in the sample. This is significant because the total area under the power spectrum is proportional to the total variance. (This follows by inverting Equation (3) and taking τ = 0. See for example Parseval’s equality in Reference Jenkins and WattsCarrier and others (1966).) Consequently, for an ideal stationary time series the power- spectrum magnitudes (as a function of frequency) will be independent of the number of sample points and, therefore, areas under different power spectra may be compared regardless of the different sample lengths with a greater area corresponding to a greater variance and hence a greater ice roughness.
B. Model calculations
Power-spectrum computations are often applied to time series that have a Gaussian distribution of points about a mean value. It is clear from plotting distributions that the sonar and laser data are not always Gaussian and, although non-Gaussian series are not precluded from analysis, it was felt to be instructive to examine certain simple mathematical functions by power-spectrum analysis for comparison.
For example, power-spectrum computations were made for a sine wave series of 300 data points, a truncated sine wave of 300 points and a series of 1 080 random numbers (from a table of random numbers). These spectra are illustrated in Figure 1. The sine-wave spectrum plot displays a strong peak where expected, as well as the mathematical artifacts inherent in the analysis of a finite record. The truncated sine wave, made by replacing all values falling beneath the mean axis with a d.c. level, has a waveform that comes closer to being a model for the real sonar data than either the sine wave or the random numbers. The spectral plot of the truncated sine wave has a peak at the same point as the original sine wave from which it was generated, but the plot also shows distinct harmonics generated by the truncating process. The spectral plot of the random number series, a perfectly Gaussian function, is, within the confidence limits expected, a straight line.
Fig. 1. Power-spectra computations for three mathematical models. The frequency in cycles per data-point spacing is given by (0.005) L where L is the wave number.
The spectrum of the truncated sine wave, while demonstrating that a simple model approximating the real sonar profile waveform generates spectral peaks where expected, also shows that the model generates harmonics which are only one order of magnitude below the desired peak. The artifacts due to the finite nature of the record and the filters in the computer program are there, but, as in the case of the sine wave, are two orders of magnitude below the signal peak and so are not likely to be confusing.
IV. Examination of sonar and laser profiles
A. Selection of sonar profiles
Two continuous sonar records, from the East Chukchi Sea and the Beaufort Sea respectively, were considered for statistical analysis. The East Chukchi Sea data were generally unsatisfactory for analysis. This was primarily because the water is shallow in the portion of sea from which the sample was taken, and consequently the direction of the submarine was constantly changed to avoid deep-draft ridges. This not only biased the sample against deeper-draft ridges but also caused the track to be non-linear in both space and time, which of course greatly decreases the confidence limits of any statistical analysis. We will consider in this paper only three spectral plots of sonar data (illustrated in Figure 2): two spectra from the Beaufort Sea which best illustrate the significant spectral peaks, and one spectrum of essentially open water from the Chukchi Sea for comparison. The open water section would be a linear path because there were no deep ridges to avoid. A more complete set of spectra is contained in LeSchack and others (unpublished).
Fig. 2. Power spectra of under-ice prof iles in the Beaufort and Chukchi Seas. The frequency scale is related to the wave number scale by f 0.0424 L (m−1) for a and b, and f 0.116 L (m−1) for c, where L is the wave number.
In selecting sections of data for analysis care was taken to select apparently homogeneous samples since this is a condition for the computation of a power spectrum. Visual examination of both records indicated that distinct ice “provinces” exist, with obvious changes in overall statistical variation or homogeneity. The province boundaries were picked visually and served to delineate the length of sample to be used for each computation.
B. Power spectrum analysis—sonar data
The power spectra in Figures 2-a and 2-b exemplify the spectra which may be obtained from homogeneous samples of under-ice profiles. The two profiles (denoted by 227–228 and 228–229) also have the advantage of being linear profiles intersecting at approximately right angles, as the submarine made a right-angle turn at 228.
The spectrum in Figure 2-a has a marked peak at 82.5 m and at 40 and 35 m. The spectrum in Figure 2-b is somewhat rougher as fewer data points were used but nevertheless has a marked peak at 56 m. We will discuss in Section 3 the hypothesis that the peak at 82.5 m in Figure 2-a and the 56 m peak in Figure 2-b correspond to the same physical ridge structure crossed in different directions. The important point here is that if the ridge structure were isotropic, then the auto-correlation taken along a profile in any direction over the same region would be approximately independent of direction (Reference McDonald and KatzMcDonald and Katz, 1969). If this were the case the spectral peaks would remain approximately the same and would certainly not shift in any uniform way. The fact that the spectral peaks appear to shift reinforces the theory that they do in fact represent a given periodic ridge structure.
It is also instructive to compare Figure 2-C with Figures 2-A and 2-B. Figure 2-C represents a very smooth profile and it is clear that the magnitude of the area under the power spectrum in Figure 2-C is less than the areas under the spectra in Figures 2-A and 2-B which represent rougher profiles. This is certainly to be expected since a rough profile has a greater variance than a smooth profile.
C. Determination of ridge orientation
Consider the following situation, using data obtained from record segments 227–228 and 228–229 as shown in Figures 2-A and 2-B. At 228, the submarine made a right-angle turn, thus changing its relative direction with respect to any ridge orientation. It was observed that 13.3 km before the turn and 6.4 km after the turn a substantial section of open water was traversed. Since there was very little open water in this portion of the traverse, it is quite possible that this water is a continuation of the same open lead. If this is true, the lead would form the hypotenuse of a right triangle whose sides are 227–228 and 228–229. It is also anticipated from knowledge of pack ice dynamics that ridge lineation might occur roughly parallel to the lead.
The geometry of this traverse, as seen in Figure 3, suggests that a determination of ridge orientation can be made by comparison of power spectra computed from straight line courses differing by a known angle. From power spectra computed for straight line paths, e.g. 227–228, 228–229, certain dominant periods appear in the spectral plot and are believed to correspond to actual ridge separations. If there were no preferred orientation then these spectra should not vary appreciably from one another. If there is a preferred orientation, however, the spectral peaks, if real, should shift when the profile direction is changed. For example, when in Figure 3 the angle θ is less than 45°, the spatial period corresponding to a particular physical ridge separation will be larger along leg 227–228 than 228–229. (By physical ridge separation we mean the average distance between ridges measured perpendicular to lineation.) Specifically, if P 1 is the wavelength of a given ridge formation ensemble obtained from the power spectrum of the leg 227–228 and P 2 is the wavelength of the same ridge formation ensemble obtained from the power spectrum along leg 228–229, then the magnitude of θ is given by
Fig. 3. Beaufort Sea Traverse segments 227–228 and 228–229 intersect at a right angle. Open water was encountered prior to 227 and subsequent to 229. The ice-ridge geometry shown above is suggested.
From power-spectra computations of segments 227–228 and 228–229, it has been determined that θ = 34°. As can be seen from the spectra, the peak occurring in the 227–228 record at wave number 29 (period = 82.5 m) has a corresponding peak approximately at wave number 41 (period = 56 m) in the record 228–229 (see Fig. 2). The ratio of P 2 to P 1 is therefore 56/82.5 (= tan 34°). In Figure 4 the spatial period scales of the spectral plots have been adjusted by a two-block averaging of the raw space series data for profile 227–228 and recalculation of the spectra for profile 228–229 with a maximum wave number of 150. This procedure modifies the frequency scales so that the significant peak pairs mentioned above line up.
Fig. 4. The power spectra for segments 227–228 and 228–229 (shown in Figure 2) have been recomputed by calculating the 228–229 spectrum with a maximum wave number of 150 and performing a two-block average on the 227–228 raw data series and then computing the spectrum with a maximum wave number of 100. The frequency at wave number L is (0.000 28 L) m−1 for the 228–229 spectra and (0.000 21 L) m−1 for the 227–228 spectra. There is now overall agreement between the spectral estimates suggesting a technique for determination of ice-ridge orientations while navigating beneath the ice. 227–228 = thick line, 228–229 = thin line.
The block-averaging process may be regarded as smoothing followed by decimation. The smoothing multiplies the spectra by some power transfer function (Reference Blackman and TukeyBlackman and Tukey, 1958) and the decimation folds high-frequency spectral components into lower ones. In Figure 5 we illustrate the power transfer functions for two- and three-block averaging in the spectral region obtained after decimation. To the extent that the aliasing caused by decimation may be neglected, the block-averaged spectrum may be converted to be a “true” spectrum by dividing it by the value of the transfer function at a particular wave number. Block averaging is preferable to straight decimation because the smoothing decreases the high-frequency spectral components prior to their being folded into lower frequencies.
Fig. 5. Power transfer functions for two-block and three-block averaging procedure. To convert a block-averaged spectrum to the true spectrum, the block-averaged spectrum at wave number L is divided by the appropriate transfer function at wave number L. (Because of the particular computational program used in this paper a factor of 1/k, that is not present in Reference Blackman and TukeyBlackman and Tukey (1958), is included in the transfer function for a k-block averaging procedure.)
With regard to profile 228–229, in our computational program, the effect of using a higher maximum wave number does not modify the true spectrum other than to give a better frequency resolution and to change the confidence limits.
It will be noticed in Figure 4 that there is now an overall agreement between the spectral peaks (and absence of peaks as at wave number 82) of the two plots. This appears to verify that these peaks are, in fact, real, and that the ridges have a definite lineation with an orientation of approximately 34°. This comparison also indicates that lineations of different frequencies are essentially parallel. That these segments are in fact of similar ridge structure is further borne out by their very similar depth-distribution plots (Reference LeSchack, LeSchack, Hibler and MorseLeSchack and others, unpublished). The computed ridge angle is of the same order of magnitude as would be suggested by the geometry shown in Figure 3.
It should be noted that there is, however, an angular ambiguity when only two legs are used for this analysis, i.e. the ridge orientation could also have an angle of —θ and still produce the same spatial period ratios. This ambiguity could be eliminated if a third leg, intersecting the other two, were traversed. This direction finding of ridge orientation has potential application for arctic submariners since the chances of encountering refrozen leads or open water might be enhanced by travelling perpendicular to known ridge orientations, which could otherwise not be easily determined.
We mention in passing that one should remember when comparing spectral peaks, as in Figure 4, that the relative amplitude is plotted on a logarithmic scale so that the peaks can be much more dominant than they appear by casual observation. Also regarding the “confidence” limits of a given point—that is the possible deviation of a spectrum for a finite sample from the ideal spectrum that could be obtained if one had an infinite homogeneous sample—the 80% confidence limits are constant factors for a given spectrum in proportion to the magnitude of the ideal peak. This is important because even though the magnitudes of peaks with higher lag number are down by orders of magnitude from lower lag number peaks, the smaller peaks can still be quite significant because the confidence limits are also smaller.
D. Analysis of laser profiles
For surface profiles, two adjacent but different ice profile types were selected. A section of laser profile data whose profile path was also covered by visual photography was analyzed. By correspondence with the visual photograph it could be ascertained that segment 3–4 of the profile represented mostly multi-year ice, and segment 4–5 represented mostly young ice. The two time series were filtered with a high-pass filter to remove the aircraft altitude variation, and then analyzed. The high-pass filtering was carried out by Fourier transforming all the points and then removing all wavelengths larger than 180 m. This type of filtering produces some error especially near the end-points and consequently the exact transfer function deviates from a sharp cut-off at 180 m. However, both the young and multi-year ice profiles were of the same length and were filtered in the same manner, so that for comparison purposes the transfer functions due to the high-pass filtering can be approximated by a high-pass filter passing all frequencies above about 1/200 m−1 and dropping to a zero pass value rapidly for lower frequencies. The high-frequency spectral values should be valid in any case, with the low frequencies being modified, but still useful for comparison between spectra. Plots of the superimposed power spectra are shown in Figure 6. The plot marked long period spacing was obtained by three-block averaging the raw data. As discussed previously the true spectra may be obtained from the block-averaged spectra by dividing by the appropriate transfer function given in Figure 5. The block-averaging procedure gives a better view of the low- frequency components with a minimum of aliasing.
Fig. 6. Power spectra of surface ice profiles from both young and multi-year ice as observed during the winter in the central Polar Basin. The frequency scale is related to the wave number scale by f = 1.667 L (m−1) for the short-period spacing plot and f = 0.556 L (m−1) for the long-period spacing plot where L is the wave number. The long-period spectrum has been obtained by a three-block average. See Figure 5 for the power transfer function.
The most noticeable difference between these spectra is the greater amplitude of all spatial periods shorter than 20 m in the young ice profile. On the other hand, for wavelengths longer than 20 m, the multi-year ice has distinctly greater amplitude. Physically, this indicates that the older, more weathered ice is more undulating than the young ice, while the younger ice, although flatter, has greater high frequency roughness. This is, in fact, what is observed in the field.
Although it is difficult to draw definite conclusions from the examination of only two profiles, the analysis technique described appears promising. The laser profile data are logistically much simpler to obtain than sonar data and, because they are obtained by aircraft, are closer to being synoptic than are data recorded from a submarine. It is believed that improvements in ice type determination and ice morphology studies could be made from airborne laser profile data processed along the lines described.
Conclusion
In conclusion, from the power spectra of sonar data it does appear that there are certain significant real spectral peaks which correspond to a periodic rather than random ridge structure. Moreover, as the direction is changed, these peaks shift, which indicates that there is a possibility of a preferred ridge orientation. A more careful study of the area under spectral peaks at different frequencies might yield information about the structural nature of ice ridges as well as the physical mechanisms involved in hummocking.
From the laser profile spectra at high frequencies a distinction between older multi-year ice and first-year ice may be made. This suggests that appropriate filtering of laser profiles could be used to determine ice types.
Acknowledgements
The authors are grateful for the assistance of the following persons at the U.S. Naval Oceanographic Office, Washington, D.C.: Mr Walter I. Wittmann, Head, Sea Ice Branch, and his colleagues, Mr Richard Ketchum and Mr P. Martin, who provided all the raw data and helpful suggestions during the conduct of the work.
Thanks also go to Dr W. A. Weeks at U.S.A.C.R.R.E.L. for a number of helpful discussions, and to Mr S. F. Ackley and Mr S. J. Mock for comments on the manuscript.
The submarine sonar data were collected under the guidance of Dr Waldo Lyon of the U.S. Naval Undersea Research and Development Center, San Diego, California.
This research was reported under Office of Naval Research contract N0014–70–C–0110 (NR–307–322) to the Development and Resources Transportation Corporation and under ARPA order #1615 to U.S.A.C.R.R.E.L.
