Calculation of cross-correlation of gray values in image pairs

The discrete correlation R of the two function values a and b over an interval n is defined as (Lemmens, 1988):

In this case a _i and b _i would be the gray values from 0 to 255 from the two images a and b, and n would be the number of samples along an image line in the chip. A correlation coefficient r can be calculated on a statistical basis by subtracting the mean value and normalizing by the standard deviation of the gray values in the chip:

The correlation estimate is done on a pair of chips, which is a subsection of the image with typical size of 32 x 32 samples. By calculating r as the chip moves across the image, one can obtain sets of r which form a correlation surface for the images (Mather, 1987). The position of a correlation peak, which represents the displacement between a detail in the two images, can be determined to a sub-pixel accuracy by means of interpolation on the correlation surface. Statistical correlation estimates were used in the previously mentioned projects on the Antarctic ice streams.

In this project, feature tracking was carried out in the Fourier domain where the correlation estimate is done on a frequency basis. The phases and amplitudes of the spatial frequencies of the gray values in the chips are calculated by the fast Fourier transform (FFT). The FFTs of the two chips are denoted F(a) and F(b) and are represented as complex images. The correlation theorem (Gonzalez and Wintz, 1987) states that the discrete correlation R between the chips can be estimated by a complex conjugate multiplication:

An inverse FFT is then performed on F(R) to obtain a cross-correlation image or correlation surface of this specific chip in the image domain (Fig. 3). The position of this peak represents the displacements between the two chips, and the determination of the position is similar to the statistical approach.

The original images were high-pass filtered before the

Fig. 3.

Examples of 128 x 128 pixel correlation surfaces in the image domain. In both cases the correlation was estimated on two chips from the same image, chosen with 15 pixels displacement in the sample direction. The left surface is estimated from chipi with high spatial frequencies, and the peak is the white spot positioned 15 pixels right of the center. The right surface is based on both high and low frequencies. The peak is smeared out by the dominant low frequencies, and the displacement cannot be measured.

correlation estimate. This removes the spatial and time differences in surface reflections and illumination on the glacier. The high-pass filtering provides an almost unchanged mean value and standard deviation for various positions of the chip on all crevassed parts of the glacier (Fig. 4). Similar correlation surfaces are therefore obtained when estimating the correlation R in the Fourier domain compared with the statistical method, where r is calculated by normalizing with the mean value and standard deviation of the chip. A higher amplitude of the low frequencies can dominate the high frequencies in the correlation surface by smoothing the surface (Fig. 3). The interpolation for subpixel accuracy will be more accurate from correlation peaks based on high frequencies, because the slope of the correlation peak will be steeper and the peak is clearly defined. The number of accepted peaks was limited by requiring that this slope be steep. By filtering away all frequencies with larger wavelength than the chip size, one also avoids some of the extra frequencies introduced as a result of the periodic assumption of the Fourier transform (Gonzalez and Wintz, 1987) when transforming the chip.

Geometrical correction of satellite images

The available satellite images SPOT HRV (23 July 1988, 153-141/0) and Landsat TM (6 August 1987, 214/4; and 31 August 1988, 215/4) had been preprocessed to level IB (SPOT Image Corp., 1985) when delivered, which means that radiometric and geometric corrections based on knowledge of the system distortion had been applied.

Each satellite image was co-registered to the coastline in a digital map, and a conformai (Helmert) transformation (Maling, 1993) between the image and the UTM system determined. The images were resampled using bilinear interpolation to 20 m pixels in the UTM projection, where a correction for the terrain displacement was also performed. The position of the satellite with respect to the scene was read from the image header file, and the off-nadir view angle was corrected for Earth curvature. The change in sample position of the nadir track for each scan line, caused by the Earth’s rotation, was described by a linear equation.

Fig. 4.

Top: Landsat TM band 4 from the accumulation area. Bottom: the same image high-pass filtered with 80% of frequencies removed.

The SPOT image was orientated using ground-control points read from a paper copy map. The Landsat TM images were co-registered to the coastline in a digital map. The rms errors of the Helmert transformation, shown in Table 1, are therefore better for the TM images.

In order to remove the remaining systematic errors, such as translation and different scale in the pair of ortho-correrted images, the images were again warped based on a Helmert transformation and bilinear resampling, with a rms error of 0.4 pixels. The images were co-registered with use of the points located on the mountains. Both Landsat TM images were warped to the SPOT image.

Construction of a DEM from digital aerial photographs

To determine how the surface elevations changed during the surge, a new terrain model was constructed from 1990 aerial photographs. A digital photogrammetric workstation (He-lava DPW) was used for image orientation and extraction of the DEM. With the use of cross-correlation, similar points in the two digital photos in a stereo model were identified and elevations calculated. The correlation process is conducted automatically, within user-specified areas and to a specified grid resolution. The resulting elevation grid can

Table. 1.

Water fluxes and chemical weathering rates at Finsterwalderbreen

Fig. 5.

Velocity contours in 1988. Points A–B are marked as in Figure 6.

be edited by inspection on the screen, providing a three-dimensional image of the stereo model and elevation contours.

Velocity

The cross-correlation estimates were carried out on three satellite images, which gave two sets of velocity vectors. The Landsat TM image from July 1987 was matched with the SPOT image from July 1988, which in turn was matched with the Landsat TM image from August 1988.

Crevasses already formed in 1987 were deformed in 1988. Where there was varying deformation within the chip size, which was the case in the accumulation area, no successful matching points were found. Near the terminus, the crevasses are covered by snow in the 1987 image. With a 1 year separation between the images, only the middle part of the glacier, along the central 6 km long mountain ridge Goldschmidtfjella (Fig. 5), gave successful matching points.

In July and August 1988, matching could be done over almost the entire glacier surface (Fig. 5). There were some difficulties near the terminus where high velocities and deformation of the crevasse pattern, together with melting of the surface, changed the pattern of the images significantly.

Typical movements during the period August 1987–July 1988 were more than 1 km, and for the 40 d period from July to August 1988 about 100 m. The average velocities are the same for both periods.

Terrain models

A DEM was constructed from ten aerial infrared photographs acquired by Norsk Polarinstitutt, Oslo, in 1990 at a scale of 1:50 000. It was developed with 10 m grid cells. There were no problems with the correlation in the surge-affected areas. The crevasses were well developed over the whole glacier in 1990, and give a good contrast in the images. Similar points are easily identified over the whole blue-ice area. In the accumulation area outside the surge-affected zone, the white snow gives no contrast in the images, and therefore similar points cannot be identified. The surface elevation of Osbornebreen constructed from the 1990 photographs has an error of about 10 m in elevation, the ground positioning about 5 m. The uncertainty is mainly a result of poor accuracy of the ground-control points (Rolstad, 1995). The method is suitable for determination of glacier surface elevations on surging glaciers, but only if accurate ground-control points are available.

The DEM from 1990 was compared with a DEM from 1966, to determine how the surge affects the glacier volume. The glacier surface is changed from 1966 to 1990 in the same way as seen on several surging glaciers (Raymond, 1987), with a reduction of the longitudinal slope. The glacier surface in the reservoir area is lowered by as much as 100 m, and the glacier front has advanced in total 2 km into the fjord, with a 100 m increase of the height at the position of the 1966 glacier front.

The terrain model constructed from the 1990 aerial photographs was inspected by calculating the surface slope in the longitudinal direction. A steep surface slope would indicate a surge front, but no such feature was detected.

Description of the velocity field

The velocity field from 1988 generally shows extending flow downstream over the entire glacier, with ice velocity increasing from stagnant to 6.0 m d^–1 at the glacier front (Fig. 5). The velocity increases steadily in the accumulation area. Along the mountain ridge, Goldschmidtfjella, the ice moves like a block with a velocity of 2.5 m d^–1. This block

Fig. 6.

Left: Aerial photograph from 1966. North is upwards. The lines show crevasse patterns which outline the same area as marked in the right photograph from 1990. The surface rivers turn to the east above point B. Right: Aerial photograph from 1990. The surface is crevassed, and an area with longitudinal and transversal crevasses is marked. The thick line shows the glacier front position in 1986 before the surge.

movement is also visible in the velocity field from 1987 to 1988. Between points B and A, the velocity decreases from 3.0 to 2.0 m d^–1 on the western part of the glacier, and increases from 3.0 to 4.0 m d^–1 on the eastern part. Before entering the glacier front area below point A, the velocity slows down from 3.5 to 3.0 m d^–1, and then increases steadily to 6.0 m d^–1 at the front. The velocity pattern on Osbornebreen is similar to what is found for the frontal zone of the continuously fast-moving nearby ice stream, Kronebreen (Rolstad, 1995). There is no distinct zone with compressive flow on either of the two glaciers. We believe that the velocity field on Osbornebreen during the active surge phase is similar to the velocity field of a tidewater glacier with continuously high ice flux.

Description of the crevasse pattern

Crevasses were present over the entire glacier in April 1987. Transverse crevasses can be seen in the satellite images from the upper area and down to point B (Fig. 5). Some longitudinal Crevasses are present on the western part of the glacier along Goldschmidtfjella, possibly as a result of shear at the lateral margin. Between B and A, longitudinal crevasses are found (Fig. 6). Further downstream of point A, a zone without crevasses is found above the transverse crevassed front zone. When not covered by snow, a very similar crevasse pattern can be seen on all the available images from the period April 1987–August 1993.

Deformation of crevasses as a result of ice flow

The crevasses are formed and later deformed by the movement of the ice (Vornberger and Whillans, 1990). An extending velocity will open transverse crevasses and close the longitudinal ones. If a surge front traveled over the entire glacier at an early stage, leaving an indistinct crevasse pattern, it could be further reduced by closing of the crevasses. The satellite images used in this project have a low spatial resolution, so narrow crevasses will not show. There is extending flow in the upper accumulation area. Along Goldschmidtfjella the ice moves like a block; the successful feature-matching over a whole year separation (August 1987–July 1988) in this area gives evidence of very moderate ice and crevasse deformation. Some longitudinal crevasses can be seen on the east part, so there is reason to believe that the crevasse pattern is unchanged here since the initial phase of the surge. Between B and A the velocity contours are no longer transverse, they have an oval shape. This will introduce both compression and transverse strain. In this area the longitudinal crevasses may be wider.

Interpretation of the crevasse pattern

Tectonic zones of crevasse patterns on various surging glaciers in Svalbard have been described and interpreted by Hodgkins and Dowdeswell (1994). Longitudinally extending tectonic regimes give rise to transverse crevasses, and longitudinally compressive regimes give longitudinal crevasses. The propagation of a surge front leaves a chaotic crevasse pattern in both the longitudinal and the transversal direction. Such a crevasse pattern is outlined between B and A in Figure 6, which has been pointed out as the position of the surge front in 1990 (Hodgkins and Dowdeswell, 1994). The upper limit of this area is near point B. As the pattern shows only transverse crevasses up-glacier from position B, there is no evident trace of a previous compression front, and point B could be the position where the surge started. This suggests that back-pressure was released below point B, and that the extending flow spread up-glacier, forming transverse crevasses in the whole accumulation area.

Terminus advance

An observation of 1.4 km terminus advance in the period August 1986–April 1987 (Dowdeswell and others, 1991) is also confirmed in the present study, which gives a velocity at the front of at least 2100 m a^–1 in the period August 1986–April 1987. This is a rather rapid advance shortly after the surge was initiated. The glacier terminus reached its maximum extension of 2 km into the fjord in July 1988, which is the same position as registered for 1990. In 1993 the front had retreated about 1 km. The crevasse pattern below B is similar in 1988 and 1990. The velocity field from 1988 shows a decrease in velocity on the eastern side from a maximum 4.0 m d^–1 down to 3.0 m d^–1 before entering the glacier front zone of continuously increasing velocity up to 6.0 m d^–1. The longitudinal slope of the DEM from 1990 was inspected to find a steeper edge which would give topographic evidence of a surge front; but no such feature was found.

Hodgkins and Dowdeswell (1994) identified two types of glacier terminus advance during surges of tidewater glaciers, depending on the position of the surge/compression front relative to a low effective-pressure zone at the glacier terminus. Steady terminus advance is when the compression front has not yet reached the low effective-pressure zone. Rapid terminus advance is when the compression front has reached the low effective-pressure zone and been eliminated because of reduced back-pressure. Osbornebreen was used as an example of the former category, described by a steady glacier front advance at an average speed of 450 m a^–1. According to the velocity field, the model where a surge front reaches the low-pressure zone and dies out could suit the situation in 1988. The decelerating area would then be the position of the compression front, now located close to the low-pressure zone. However, this deceleration is further downstream than the longitudinal crevasse pattern outlined as the surge front in the 1990 aerial photograph shown in Figure 6.

Variation in crevasse pattern and velocity as a result of bedrock topography

Bedrock topography influences the velocity and the crevasses. We inspected the area A–B on the aerial photographs from 1966 during the quiescent phase, shown in Figure 6. Transverse crevasses are found near point B, and some longitudinal crevasses outline the area where the chaotic crevasse pattern is identified on the 1990 photographs. Above point B the melt streams are turning east. According to the 1966 map, the glacier surface has an elevation difference of about 20 m in the transverse direction. This surface slope and the crevasse pattern suggest a rise in bedrock topography in this area. The 1988 velocity field shows, as mentioned, a decrease in velocity in the western part of the area A–B and an increase on the cast side (see Fig. 5). We believe the flow is rotating slightly and is moving faster as a result of thicker ice on the deep bedrock on the east side. Along the mountain ridge, both velocity fields show that the ice flows like a block (see Fig. 5). This indicates that the bedrock forms a smooth channel along the mountain ridge. The ice is not deformed as it flows along Goldschmidtfjella. When it reaches point B a bump in bedrock can form new crevasses, and further down the widening of the glacier gives a velocity/strain field which could open these longitudinal crevasses.