Bam, Iran – Earthquake detection with the ASAR instrument of ENVISAT

 The purpose of this article

The devastating major earthquake which occured on December 23, 2003 with a moment magnitude of 6.6, that caused a huge number of human casualties, as well as major destructions in Bam and the surrounding Kerman province in Iran, had been the topic of several geophysical studies in the past 10 years. Furtunately, and due to the scientific effort of the seismological societey, the basic processes that induced the earthquake are well known by now. This short article should provide some useful information for processing spaceborne Radar data for the detection of earthquakes and faults. For this purpose I have chosen ASARAdvanced Synthetic Aperture Radar data from the ENVISAT mission before and after the event to determine the structural differences and rupture zones induced by the earthquake.

The data sets provided by ESA are:

subset_9_of_ENVISAT-ASA_IMS_1PXPDE20040211_061300_000000142024_00120_10194_0013_ML.dim

subset_10_of_ENVISAT-ASA_IMS_1PNUPA200331203_061259_000000162022_00120_09192_0099_ML.dim

Which processing steps were applied to the product?

As post-processing is very important for Radar data, several corrections and image processing algorithms have to be applied on the data. According to the product name the Single-Look Complex (SLC) image data is a high-resolution, narrow swath products based on data acquired at one of seven subswaths. It is intended for SAR image quality assessment, calibration and interferometric or wind/wave applications. A small number of corrections and interpolations are performed like the preprocessing (for all data the same; i.e. I/Q bias removal, I/Q gain imbalance correction, I/Q non-orthogonality correction), a Doppler Centroid Estimation and a Range Doppler Image Formation. To obtain the Doppler Centroid Estimation for the high-resolution Image Mode of the product, the Multi-Look Cross Correlation Method is used.

Figure 1: Multilook processed image of Bam, Iran from 3rd December 2003, with the airport in the north east of Bam

NEST Multilook Operator

The original SAR images with inherent speckle noise were processed with the Multilook Operator to reduce the speckle noise by incoherently combining several images as if they corresponded to different looks of the same scene and to produce an image with nominal image pixel size (azimuth to range ratio of 5:1, reprocessed to 1:1). As a result, the final images approximately have square pixel spacing.

Range and Azimuth Direction obtained from the metadata

The original intensity image has num_output_lines=5000, num_samples_per_line=1200 with an azimuth_looks=1.0 and range_looks=1.0. The processed image has num_output_lines=1000 with azimuth_looks=5.0. Therefore the azimuth direction is from north to south whereas the range direction is from west to east.

Why is the automatic coregistration necessary?

The process accurately co-register the master image (11th February 2004) with the slave image (3rd December 2003). For this purpose the master image is used as a reference image and the slave image is overlayed in a way that each pixel of the slave image corresponds to the same geoposition pixel in the master image. Due to the fact that both images are complex images and that the scene described by the images is in a flat area the build-in algorithm ensures a matching accuracy of better than 0.05 pixels. For complex images the registration is achieved by maximizing the coherence between master and slave images at a series of imagettes defined across the image. After the registration, a warp function is created which maps the master-slave Ground Control Points (GCP) of the slave image into pixels in the master image.

What information can be obtained out of the residual file?

The Create Stack Operator allows collocating two spatially overlapping products. Collocating two products implies that the pixel values of one product (the slave) are resampled into the geographical raster of the other (the master). When two products are collocated, the band data of the slave product is resampled into the geographical raster of the master product. In order to establish a mapping between the samples in the master and the slave rasters, the geographical position of a master sample is used to find the corresponding sample in the slave raster. The calculation of the new pixel value is performed by weighting the 36 (for the chosen settings) surrounding pixels.

Figure 2: Cubic convolution resampling method for 16 (4×4) sourrounding pixels

For the parse 0, every valid GCP is used to resample the projected output product (by Cubic Convolution) of the slave image. The respective RMS errors in x- (row) and y- (column) direction are listed in the residual file. All GCP-pairs of the master and slave image that exceed the RMS threshold defined previously are excluded from the 2nd parse of the same routine. As a result, the RMS mean gets better by each use of the Cubic Convolution. As the RMS values below show, the accuracy increases significantly from Parse 0 to Parse 1, but only very poor from Parse 1 to Parse 2. As the desired accuracy of 0.05 was already reached in Parse 1, the process could have been stopped after the first parse.

Figure 3: Automatic co-registered image with calculation of the coherence and Multilook processed with post earthquake event as master image and prior earthquake as slave image. Interessiting in this image are the areas indicated in red because it might represent a crack/surface rupture induced by the earthquake.

Figure 4: Multilook processed image with overlayed orthorectified optical (RGB) satellite imagery of Bam and the surrounding region Kerman. The green areas are vegetation and the greyish to brownish structures are builidings and roads. The satellite image was obtained 2009, but the extent of the city Bam remained nearly constant over the time. The red arrows indicate the satellite azimuth and range direction at the time the radar image was acquired.

 

 

The definition of coherence is the similarity of the co-registered GCP-pairs (by Coarse Registration and Fine Registration by performing a cross-correlation of master and slave imagettes) and a computation of the best sub-pixel shift of GCP such that the slave imagette at the new GCP position gives the maximum coherence with the master imagette according to Powell’s method. Assuming I1 and I2 be RxC imagettes and I*2 is the complex conjugate of I2, coherence is computed by:

Method:

The coherence in this scene is computed with a 10×2 (azimuth x range due to the 5:1 looks ratio) sliding window in two steps:

1. First for each pixel in the imagette, a 10×2 window centered at the pixel is determined for both master and slave imagettes, and coherence is computed for the two windows using equation above.
2. Average coherences computed for all pixels in the imagette to get the final coherence for the imagette.
   The low coherence areas (dark spots) correspond mainly to building structures like houses and roads as well as vegetation, whereas the high coherence areas (bright spots) correspond to the surrounding desert in the scene. Due to the earthquake, a phase-shift occured on the building structures (destroyed buildings) and the vegetation differs naturally with each image. Therefore the coherence value is very low for the city area whereas the surrounding desert isn’t effected that much, because the phase-shift is quite low (except for the new crack) due to a lack of horizontal or vertical displacement and therefore low differentiation.

The low coherence areas (dark spots) in Figure 3 correspond mainly to building structures like houses and roads as well as vegetation, whereas the high coherence areas (bright spots) correspond to the surrounding desert in the scene. Due to the earthquake, a phase-shift occured on the building structures (destroyed buildings) and the vegetation differs naturally with each image. Therefore the coherence value is very low for the city area whereas the surrounding desert isn’t effected that much, because the phase-shift is quite low (except for the new crack) due to a lack of horizontal or vertical displacement and therefore low differentiation. See Figure 4 for a more natural looking image of the scene.

Why do we have to remove the topographic phase?

InSAR is normally applied for estimating topography. Here we are interesting in the displacement only, meaning that we need to subtract the phase difference caused by topography. That can be done if a DEM is available or by means of a third SAR image. The topography around Bam is rather flat, and therefore the interferogram doesn’t change much.

The interferogram and estimate of the displacement

This interferogram (Figure 5) represents contours of ground displacements towards and away from the Envisat satellite that occurred during the earthquake (each fridge corresponds to half the wavelength of the radar = 2.8cm). There is a two-quadrant pattern consistent (in this descending orbit image) with a near-vertical strike-slip fault oriented north-south (Figure 6). The ground moved approximately 30cm towards the satellite in the south-east quadrant, and approximately 15cm away from the satellite in the north-east quadrant.

Figure 5: Phase-Interferogram with topographic correction of the scene and applied colour palette spectrum.spd

Figure 6: Interferogram with overlayed optical image data and indicated faults. The new fault (A’-B’) occurred approximately 4km west of the known fault (A-B)

Nico Trebbin is a geophysicist and future satellite application engineer who studied at the LMU/TUM in Munich and is currently working at DLR