Using repeat-pass interferometric synthetic aperture radar (SAR) data acquired over the Belridge and Lost Hills oil fields from several satellites, we have measured and mapped ground subsidence of up to 6 cm over a 70 day period, corresponding to an annualized subsidence rate of over 30 cm per year, and have mapped differential subsidence rates over the producing fields. Our spaceborne measurements are accurate at the millimeter level.

Timely, cost-effective remote monitoring of surface deformation at the millimeter level supports assessment of subsidence induced by hydrocarbon production, property damage, and the resulting surface drainage problems. As part of a new research effort to evaluate this technology, we have tasked the Radarsat satellite to acquire additional SAR imagery over Belridge, Lost Hills, and the Rocky Mountain Arsenal for monitoring ground subsidence, observing the progress of large-scale grading operations, and generating accurate digital elevation models. Without any additional information, it is now possible to create digital terrain models from spaceborne imagery with one meter absolute accuracy (vertical) anywhere in the world. Introduction of ground control points or survey markers into the data set can greatly improve the overall accuracy of the resulting model.

This project is an example of highly leveraged research we are conducting jointly with Atlantis Scientific and is partially funded through a peer-reviewed program sponsored by the Canadian Space Directorate.

SAR interferometry offers geophysicists their first opportunity to measure surface deformation with continuous spatial coverage and millimeter precision. Since the launch of the ERS-1 satellite by the European Space Agency ESA in 1991, the topic of interferometric processing of the signals from Synthetic Aperture Radar (SAR) has gained a lot of attention in the radar Remote Sensing community. With this very sensitive measurement technique it becomes possible to measure deformation effects of the Earth's surface over large areas with a sub-centimeter accuracy. Massonnet’s demonstration of the power of the technique published in Nature in 1993 showed that this method was capable of measuring the horizontal and vertical ground displacements due to relaxation after an earthquake. Comparison of the satellite measurements with numerical models showed remarkable resemblance. Besides deformation measurements, the technique also allows the creation of digital elevation maps, even in areas dominated by cloud cover (equatorial regions) inaccessible by optical sensors.

A Synthetic Aperture Radar (SAR) is an imaging radar device which images the radar backscatter of the surface of the earth over large areas (100 km x 100 km is typical for a satellite system) with a moderate spatial resolution (10-20 m). In the satellite not only the power of the backscattered electromagnetic radiation is recorded, but also its phase. In radar interferometry the change in phase is measured between the backscatter in two different SAR images of the same target terrain, taken from slightly different positions or taken at different times. In repeat-pass SAR interferometry, two SAR images of the same area from almost identical perspective are coregistered and combined into an interferogram. The phase difference for each pixel of the interferogram is a measure of the change in distance between the scatterer and the SAR antenna between images. From the phase information, a digital elevation model (DEM) can be derived if no large scale deformation occurs between the recordings.

Very small changes of the path length can thus be measured and related to either topographic height or temporal surface deformation. Because the radar is an imaging device the information on deformation or height can in principle be obtained for every image pixel within the scene. This is a unique feature of this technique.

Comparison with other Techniques

There are now several techniques for deriving topographic maps from space imagery. The use of digital versions of conventional photogrammetry with SPOT Panchromatic imagery yields topographic maps with spatial (XY) accuracies of about 10 m and absolute vertical accuracies of 20 m. In practice, we have found that the consistency, or precision, of vertical measurements over large areas within an image can be as great as 1 meter.

With relatively little modification, digital photogrammetric techniques can be applied to side-looking radar data. Two satellite passes are still needed to form a stereo pair, and a complication is introduced by the fact that present radar satellites "look" in only one direction, so that in order to get a sufficient difference in the parallax angle required to perform radargrammetry, one must use one image from an ascending orbit and the other image from a descending orbit, with about 20 degrees difference in look direction between the two tracks. However, the great advantage of radar data is that clouds and darkness do not interfere with the acquisition of data. The resolution and accuracy of current spaceborne radar systems is about the same as that of the optical SPOT satellite. Airborne systems can obtain better spatial resolution by flying lower, but the absolute vertical accuracy remains about one meter.

However, one can also use interferometric techniques to derive a topographic map from radar data, because the instrument can be oriented to measure coherent scattering from the target surface. This still requires two satellite passes, but they must both be in the same direction to preserve target phase coherency. Because phase information is used in this technique, the vertical precision is much greater than that of radargrammetry. Because this technique can detect changes in ground elevation of less than a millimeter, it is useful for mapping subsidence, tectonic deformation before and during earthquakes and relaxation after them, as well as inflation of volcanoes due to magma movements. Although it cannot measure changes due to mudflows or human activities, because it relies on the correlation of the small-scale roughness to provide a coherent scattered signal, it can detect them. The method also breaks down when the signal is decorrelated due to a long time interval between passes, or when vegetation growth is rapid compared with the frequency of data acquisition.

Although digital photogrammetry is highly compatible with GIS systems, Land Surveyors, who are used to precision on the order of millimeters rather than tens of meters, rarely use it. This points to a fundamental difference between the three types of data:

• Conventional Survey Data and GPS: a sparse network of very precisely measured points, with all else being interpolation of varying degrees of reliability.
• Aerial Photogrammetry: contours and a sparse network of control points are known very precisely, and all else is interpolation. Based on parallax.
• Interferometry: the position of every pixel is known somewhat accurately, but no point is identified precisely. Based on phase correlation.

Photogrammetry. The most important difference between photogrammetry and radar interferometry is that photogrammetry is based on the concept of parallax, whereas interferometry actually depends on coherency and interferometric phase differences. This is translated into the need for a relatively large baseline for photogrammetry and a very small baseline for SAR interferometry. For photogrammetric applications, a major restriction is cloud cover. Since SAR is an active sensor, it creates its own illumination, at a wavelength that penetrates through clouds. A significant difference in the imagery is the layover effect, which is away from the sensor in photogrammetry and towards the sensor for INSAR.

GPS. Although GPS is an ideal system for point positioning, especially in its differential mode, one would need some 37.5 million GPS receivers to reach the same coverage as a SAR quarter scene. Furthermore, due to the incidence angle, a SAR is especially sensitive for height information, whereas the height for GPS is usually the weakest component. However, GPS permits long time series of surveys, whereas SAR can be seen as an instantaneous measurement. Another important difference between the techniques is that GPS yields absolute coordinates, whereas SAR interferometry only performs relative difference measurements. GPS and INSAR can be regarded as complementary techniques: spatial versus temporal continuity.

Optical Remote Sensing. For optical remote sensing applications, similar comments can be made as for photogrammetry. A major restriction

is being formed by cloud cover: in Europe alone only 10% of all acquired data is useful for DEM purposes, and equatorial regions have even more cloud cover. SAR is an active sensor, it creates its own illumination and is unaffected by clouds. Another problem is formed by those areas that do not exhibit good optical recognition characteristics, like the white snow coverage in polar regions. In these areas, stereoscopic parallax is hard to find. With SAR interferometry, similar problems exist when the scattering characteristics of the topography change, resulting in a loss of coherence.

Leveling. At the present time, INSAR can reach the leveling accuracies of terrestrial height measurement techniques only in ideal circumstances. However, the feature that makes INSAR so promising is the possibility to get accurate height change information at every pixel from only two observations. This nearly continuous coverage of an area located anywhere on Earth makes it a very attractive technique when compared to point-positioning geodetic techniques such as leveling or differential GPS. Furthermore, since the height information obtained by INSAR is a relative measurement, it can play a complementary role if combined with absolute GPS measurements. The practical choice for a certain height determination technique is dependent on a number of considerations: the size and type of the area of interest, the desired coverage and accuracy, the costs, and the availability of data.

Estimating Heights and Deformations

A SAR instrument emits microwave pulses to the Earth at a specific look-angle (23 degrees in the case of ERS-1). These pulses are reflected at the Earth's surface, with a small part of this backscattered energy being detected at the SAR antenna. The two-way travel time is measured, discriminating between targets in the range direction, perpendicular to the flight direction. If only the travel time of the transmitted pulse were to be measured, the achieved height accuracy would be in the order of 10-20 meters, although the ambiguity between two points at the same range distance can not be resolved.

The signal received by the SAR antenna contains not only an amplitude, which depends on the target reflectivity, but also a phase. The phase is proportional to the travel time, and therefore the distance traveled, and can be measured to a fraction of the emitted SAR wavelength.

The observed phase is actually a summation of the phases of all the scattering elements located within a resolution element. The absolute phases of pixels in a SAR image are actually quite random, and a phase plot of only one SAR image would not contain any useful information. If, however, the summation of all these scattering elements is the same for a second image - the signals are coherent - the differences between the phases in the first and second image will vary only due to path length differences.

The distance from a satellite sensor to the target can be measured very accurately, as can the baseline distance between SAR antennae. This information is used to solve geometrically for the height of the target.

Height estimates are averaged over a surface resolution element of the radar image (often termed a picture element or "pixel"), typically tens of meters in diameter. In the case of ERS-1, each pixel is 25 m on a side. RADARSAT has variable spatial resolution, approaching 6.25 m.

To obtain reasonably accurate height estimates, how well do we need to know the viewing geometry? The uncertainty in height is related to

uncertainty in sensor orientation angle.

For a spacecraft system where the range to the target is large (at least several hundred kilometers), the attitude error must be very small to achieve acceptable height accuracy. For example, to achieve a height accuracy of 2 m from a 500-km-altitude spacecraft, the attitude error must be no more than about 1 arc sec. For an L-band radar system (wavelength 26 cm) with a baseline of 1 km, this implies a relative position measurement accuracy (i.e., knowledge of the position of one antenna with respect to the other) of about 5 mm. This accuracy is difficult but feasible with modern spacecraft tracking techniques, such as GPS. Errors in the baseline orientation cause errors in absolute elevation, e.g., with respect to the center of mass of the Earth, but the relative heights within the radar swath are largely unaffected. This means that a few ground control points can be used to calibrate absolute height estimates even in the absence of arc-second attitude information.

When the relative positions of the radar scatterers within a pixel change by amounts greater than the radar wavelength, we are out of luck - phase coherence will not be maintained in the second image relative to the first, and we cannot perform a pixel-by-pixel phase comparison between the first and second image. We cannot recover topography with a single-antenna system, nor can we say anything quantitative about the nature of the change, beyond stating that it has been larger than some minimum value. However, if we had the foresight and resources to acquire our first set of data with a dual-antenna system (thereby obtaining topography on the first pass), we could at least perform a differencing operation on two independent height estimates, obtaining displacement information precise to the several-meter level.

Now consider a second case, where the only surface change is a large-scale coherent change common to several adjacent pixels. By this we mean that, within a given pixel, the position of radar scatterers has not changed to any significant degree, but the ensemble of scatterers (i.e., the entire ground surface within that pixel, as well as some adjacent ones) has moved up, down, or sideways in some correlated fashion. Now we are able to perform a phase comparison of the two images. The differential phase contains information on the range change to the radar antenna, thus providing information on the surface change, precise to a fraction of a radar wavelength, or to a few millimeters to centimeters for typical radars. To measure surface displacement to this precision requires a relatively accurate a priori estimate of topography in the region and information on the position and orientation of the antennas during each image acquisition.

One measurement of phase change based on two successive images constitutes only one component of the three-dimensional surface displacement vector, specifically the component projected onto the spacecraft-target vector. In this respect, it differs from the three-dimensional vector measured by GPS geodesy. The great advantage of SAR is that the observations are acquired with virtually complete spatial coverage, rather than the sparse observations inherent in conventional (including GPS) geodesy. In principle, all three components of surface displacement could be obtained from space; in practice, we are only likely to obtain at most two, if data from both the ascending and descending passes of the satellite are available. Additional information (e.g., from ground observations) would normally be required to determine the full three-dimensional displacement field. In practice, one or two components will be adequate for most applications.

We have said that measurement of surface displacement with SAR interferometry depends on the nature of the surface change. In fact, the conditions under which SAR interferometry gives useful information on surface change is a major topic of current research. We can restate two important necessary conditions for detecting and measuring surface change with SAR:

• Changes between successive images must not be too large; specifically the displacement gradient across a pixel must fall within some value.
• The radar-scattering characteristics within each pixel must remain similar in the time between the two image acquisitions, specifically the root-mean-square (rms) position of the surface scatterers within a pixel must remain constant within a fraction, say 10 to 20%, of the radar wavelength.

The first condition is unlikely to pose significant problems, at least in theory. When it is violated because natural phenomena generate large displacements and displacement gradients (e.g., volcanic eruptions, landslides), simple differencing of before and after topographic data should be adequate to define the changes to the degree of accuracy required for their study. Of course, obtaining such topographic data may be quite challenging in a practical sense.

The second condition is more problematic. When this condition is not met, it is termed temporal decorrelation, and it constitutes one of the major problems for SAR interferometry. Temporal decorrelation precludes or makes difficult the phase comparison of the two SAR images. Temporal decorrelation has been observed on time scales as short as a few hours in vegetated areas experiencing windy conditions. On the other hand, the initial result for the Landers earthquake demonstrated phase correlation over several months in a desert area with no storms.

Consider again two radar antennas observing the same ground swath but at different times, as in the repeat-pass mode. The measured phase at each point in each of the two radar images is equal to the sum of a propagation part, proportional to the round-trip distance, and a scattering part due to the interaction of the wave with the ground. If each pixel on the ground behaves the same for each observation, calculating the difference in the phases removes dependence on the scattering mechanism and gives a quantity dependent only on imaging geometry.

It should be noted that optimum baseline length is a function of radar wavelength and the desired result. For surface change estimates, similar imaging geometry in the before-and-after images is desirable, implying the shortest possible baseline. For topographic estimates, the baseline should be neither too short (reducing sensitivity to surface height) nor too long (causing geometric characteristics of the image to differ by too much and preventing the phases of the two images to be compared interferometrically because they are uncorrelated). Thus, an interferometer designed to recover high accuracy height estimates represents a compromise between these competing effects.

For C-band (~6-cm wavelength) ERS-1 SAR data, the optimum baseline is ~200 m, and antenna separations of 100 to ~400 m can be used to obtain good height estimates. One meter of topography gives a phase signature of 4.5° (actually less than the real noise limit of ~20° , implying that ERS-1 is not sensitive to topography at this level). However, for the same pass pair, a 1-m surface displacement yields a phase signature of 12,800° , or nearly 3000 times greater sensitivity. Since we seek to measure 1-cm surface changes, this implies that we require topographic data accurate to about 3000 x 1 cm, or +/-30 m.

However, in the differential case, a change of 1 cm in the range direction yields a phase difference of 127° , which is easily detectable.

Naturally, from a practical point of view a zero-baseline is an ideal and therefore unusual situation. Therefore, the technique of differential interferometry is applied for measuring deformations. The problem in repeat pass interferometry with a non-zero baseline is that phase differences can be caused either by spatial or by temporal path length differences.

The former is determined by the length of the effective baseline, the latter by changes in the topography. The corresponding accuracies are also dependent on the baseline for relative height differences (yielding an accuracy of some meters), whereas the topographic changes in between the two image acquisitions can be measured as a fraction of the wavelength (sub-centimeter accuracy).

In differential interferometry, a priori information on the height distributions used to predict phase differences, and then these phase differences are stripped from the interferogram, leaving only phase differences caused by surface deformations.

Two possibilities for a priori height information are an existing DEM obtained with other techniques, or a DEM obtained from a previous interferometric SAR pair, of which it is known that no deformation took place. Both techniques have been extensively demonstrated on seismic deformation, volcano deflation, glacier movement, and ground swell.

Error Sources

Since SAR interferometry only works under coherent conditions, in which the received reflections of the pulse are correlated in the two SAR images, this is evidently the most important limitation and error source. The loss of coherence, called decorrelation, can be caused by a number of sources. Decorrelation due to the registration of the two images can be distinguished into two parts. The spatial coregistration and the necessary resampling introduce noise in one of the images, resulting in a reduced coherence in the interferogram. The spectral registration error (also referred to as baseline decorrelation) is caused by the different mapping of harmonic features of the scene to the two different look directions of the satellites. For repeat pass interferometry, temporal decorrelation is a key problem in large parts of the world - major areas within a SAR scene change due to weather or antropogeneous activity. Orbit errors introduce phase errors, which translate to height errors. Atmospheric inhomogeneities can cause a delay of the electromagnetic signal, resulting in a phase bias in the interferogram, whereas changes in the dielectric constant of the terrain have a similar effect. A decorrelation effect can be directly coupled to the variance of the interferometric phase, and hence to the parameters (height or height change) which are inferred from this phase.

Other error sources can occur in the interferometric processing sequence. Phase unwrapping, the procedure through which phase ambiguities are resolved, can be quite cumbersome in rugged terrain, leading to local or global errors in the inferred height.

In differential interferometry, finally, an error in the applied a priori elevation model yields artifacts in the deformation maps, which are produced using these DEMs.

It can be noted that possibilities exist to overcome or suppress some error sources. The multi-look procedure is a trade-off in which spatial averaging (hence, a loss in image resolution) is performed to obtain better phase statistics.

Belridge Case Study

At Shell's South Belridge Field, removal of fluids from the shallow, thick, low strength Diatomite reservoir has led to substantial reservoir compaction, causing surface subsidence and well failure. In some areas over 10 feet of subsidence have occurred, and well failure rates have averaged over 3% per year, costing several million dollars per year.

A program of monitoring and research has been underway to address this issue: an updated elastoplastic rock law based on laboratory data has been generated and incorporated into a finite element geomechanics simulator, and since the late 1980s a network of level survey monuments has been used to monitor subsidence. It has been shown that water injection reduces subsidence from 7-8 inches per year to 1-2 inches per year, but does not abate well failure. The goal of this program is to generate the forecasting tools to assess the impact of primary and waterflood operating policies on subsidence and well failure, and to devise economic methods to mitigate any adverse consequences.

Diatomite reservoir thickness ranges from 800-1000 feet and is overlain by 400-600 feet of Tulare formation sands and shales and 200-400 feet of "air sands". The density of the Diatomite formation, which is composed of fragments of diatomite tests, is only 2.2-2.5 g/cc (sic), its porosity is 55-60%, and its Young's modulus only 25,000-80,000 psi (cf. >1,000,000 psi for typical sandstones).

An area in the NW quadrant of section 33, which had experienced a large number of well failures, was selected by CalResources for field calibration of their geomechanical model. This area had subsidence of 7 inches per year at the center of the subsidence bowl during primary production, decreasing to less than 1 in/year after waterflood. It had experienced 9 years of primary production from 1978, followed by 11 years of 1.25 acre waterflood. The model predicts that a proposed accelerated waterflood on a 5/8 acre spacing, with a 3/1 ratio of producing to injecting wells rather than the current 1/1 ratio will increase subsidence rates and shear strains in the area.

At present there is no demonstrated quantitative relationship between surface subsidence and well failure. Qualitatively, however, CalResources expects greater well failure rates in or adjacent to areas of high subsidence. Dogleg failures in wells suggest a shear mechanism of failure, and the model predicts highest shear beneath areas of high surface gradient, which is being monitored by means of tiltmeter arrays.

We compared the SAR Interferometric measurements of surface subsidence with the surveyed monument records. INSARresults compare favorably with the monument survey data. We found that it is indeed possible to measure the subsidence from ERS pairs accurately and cost-effectively. Profiles were extracted throught the monument survey points and the deformation images. Five ERS-1 and ERS-2 images were analyzed with variable temporal ranges. To simplify comparison, we normalized all data to a minimal 34-day deformation rate.

We observed some discrepancies in location between the SAR Interferometric deformation maps and the monument measurements, probably due to the fact that INSAR pairs were located from ephemeris data only and not georeferenced from ground control points. These spatial differences were adjusted by calculating shifts and applying them to the SAR Inteferomtric data to make it correspond more closely with the monument measurements. We observed shifts up to 750m

Applications of Digital Elevation Models (DEM) have proven to be of significant value for the Oil and Gas Exploration industry. These technologically advanced maps save industries large sums of money when planning future exploration efforts. By replacing the present day low resolution DEMs with SAR Interferometric DEMs the accuracy is increased while the cost decreases. There are multitudes of advantages for using INSAR technology to produce high grade DEMs. The fact that DEMs are generated over a wide area for a fraction of the cost makes this appealing to consumers. The Interferometric DEM below is compared with a standard USGS DEM, of higher resolution and accuracy.

The DEM below was created from 2 SAR images of the same area taken from an almost identical perspective. The images were coregistered and then combined into an interferogram in which the phase difference for each pixel is used to derive a DEM.

Interferometry has provided the same answers that surface surveying and instrumentation does, but at a fraction of the cost. The interferometric data has shown more spatial detail in the subsidence pattern than can the repeated surveying of fixed monuments. The spatial and temporal data show shifts of up to 750m in the ERS images. The estimated shifts can be estimated by using ground control subsidence measurements that will appropriately correspond to one another.

Our results show a great deal of spatial as well as temporal detail. We know already that we can obtain usable passes much more frequently than the once a year that the leveling has been done. We also are getting a very good feel for the time scale of SAR decorrelation in this type of environment. Much of the noise in our data, due both to radar speckle and to decorrelation, can be cleaned up through suitable image processing techniques, such as low-pass filtering.

In the future we plan to resolve these image calibration anomalies permanently. We also hope to be able to derive changes in surface slope and compare it with the tiltmeter measurements. We will also determine the cost of doing interferometric studies and the likely frequency of obtaining usable overpasses, as well as the actual limitations of the technique in practice.

This case study at Belridge and Lost Hills, CA is an example of highly leveraged research Shell Oil is conducting jointly with Atlantis Scientific and is partially funded through a peer-reviewed program sponsored by the Canadian Space Agency. For the fiscal year of 1997 Shell Oil has applied $90,000 to this research effort. There has been regular contact between the researchers at Shell and Atlantis Scientific. Contacts with Tim Perrott of Atlantis Scientific have addressed issues of data exchange. Atlantis Scientific has provided georeferenced and calibrated deformation maps, DEMs of Belridge and Keeks Corner test sites for both ERS and RADARSAT sensors, a USGS 7.5 minute DEM and detailed analysis and comparison between RADARSAT and the USGS DEM for Keeks Corner.

In May E.K. Biegert visited Atlantis’ offices in Ottawa for technical discussions and training fin the use of the INSAR software previously installed at Shell’s Bellaire Technology Center, including the use of different data sources, and data selection and ordering. During the visit the Belridge deformation maps were analyzed and compared to the monument measurements. There were also discussions on the future joint research project involving the analysis of an Interferomatric SAR DEM of Shell’s Rocky Mountain Arsenal test site in CO.

Within Shell we continue to publicize the excellent results we have discovered at the test sites. Presentations have been made at Shell on numerous occasions to the Geology and Geophysics research staff and at the Geoscience and Reservoir Engineering Conference at Shell which includes 450 participants drawn from all of Shell’s international locations.

We have also presented our results in a number of public forums including the 1996 RADARSAT Applications Workshop in Houston, the ERS symposium in Florence (March 1997), the American Association of Petroleum Geologists’ annual international meeting (Dallas, April 1997), and GER’97 (Ottawa, May 1997), as well as the Applied Geologic Remote Sensing ERIM Conference (Denver, November 1997). The inclusion of this technology on our webpage also provides for increased public awareness. We plan to continue this research with Atlantis Scientific and work towards a publication in the open literature.

In response to the internal desire of Shell to use the technology, we are developing training material and courses to spread awareness of this exciting capability to our domestic and international staff and are working closely with Atlantis to improve our in-house capability to perform data processing and generation of subsidence maps.