Slip analysis of the Kokoxili earthquake using terrain-change detection and regional earthquake data

On 14 November 2001, the Kunlun fault experienced a M_w 7.8 earthquake, which became known as the Kokoxili earthquake. Since the earthquake, several studies have been conducted along the faulted region in order to describe the slip direction, amount of slip, deformational style, and principal stress orientations (Lin et al., 2002, 2003, 2004; Andronicos et al., 2004a, 2004b; Ozacar and Beck, 2004). In this study, we used terrain-change detection using satellite imagery to model slip along the rupture zone of the 2001 Kokoxili earthquake. The resulting vector displacement fields were then analyzed using a strain inversion algorithm in order to determine the principal extensional and shortening axes. We compared both the resulting displacement fields and the infinitesimal strain axes to the results of the previous studies in order to understand the pattern of deformation along the rupture zone. To first order, the deformation field is influenced by the fault's segmented geometry.

The Kokoxili earthquake created a 400-km-long surface rupture and as much as 16.3 m of left-lateral strike slip (Van der Woerd et al., 2002; Lin et al., 2002, 2003) (Fig. 1). This slip is most evident as fault scarps developed on alluvial fans and terraces ⁰¹. Lin et al. (2003) divided the portion of the fault that slipped in 2001 into four segments based on morphological characteristics, geological structures, and displacement distributions (Fig. 1). The four segments are (from west to east): Buka Daban Peak, Hongshui River, Kusai, and Kunlun Pass. All segments exhibit sinistral displacement, but each segment also displays small amounts of either transtensional or transpressional deformation. For instance, the northern part of the Hongshui River segment displays a small transpressional component, which is manifested in mole tracks where the fault bends to the north (Lin et al., 2004). The Buka Daban Peak and Kusai segments both display transtension as indicated by the occurrence of many extensional cracks (Lin et al., 2003).

The field observations by Lin et al. (2003) are in agreement with a regional stress analysis by Andronicos et al. (2004a, 2004b). Their work used focal mechanisms for M_w > 5.5 earthquakes from the Harvard Central Moment Tensor catalogue that occurred between 1977 and 2004. For each event, the trend and plunge of the axes of greatest (P) and least compression (T) were extracted and sorted into five categories: thrust, transpressive, strike slip, transtensive, and normal. When Andronicos et al.'s (2004a, 2004b) results are combined with the results from Lin et al. (2003), the deformation types observed along the fault correlate with the fault's segmented geometry (Fig. 1). For instance, Lin et al. (2003) observed tensional cracks and sinistral displacement ⁰¹ in the Kusai segment north of Kusai Lake, and Andronicos et al. (2004a, 2004b) showed transtensive stresses north of the fault and strike slip south of the fault for that same segment (Fig. 1). This correlation is also depicted in an earthquake rupture model for the Kokoxili earthquake designed by Ozacar and Beck (2004), who used inverted teleseismic P waveforms, differing types of focal mechanisms (ranging from pure strike slip to strike slip with compressional and tensional components), and source models.

The coseismic slip resulting from the Kokoxili earthquake can be measured with optical remote-sensing data. To do this, image processing algorithms are needed to measure apparent offsets in the geographic locations of corresponding pixels in two (or more) images of the same portion of Earth's surface taken at different times. These interimage pixel offsets define vectors with orientations that indicate the direction of terrain displacement and lengths that denote the magnitude of that displacement. Previous work with Satellite Probatoire d'Observation de la Terre (SPOT) images has shown the feasibility of using optical imagery for lateral displacement change detection using Fourier methods (Van Puymbroeck et al., 2000; Dominguez et al., 2003; Michel and Avouac, 2002; and Feigl et al., 2002). For example, Dominguez et al. (2003) were able to resolve a coseismic displacement of 5–11 m (11 m near the epicenter) along a major thrust fault associated with the 1999 Chi Chi earthquake in the Central Range of Taiwan with this approach.

Two types of pixel offsets exist: true and apparent. Apparent pixel offsets can originate from a variety of sources and are not related to terrain changes. These apparent pixel offsets can be caused by topography, the geometry of the imaging system, the attitude of the spacecraft, and image decorrelation due to seasonal or atmospheric effects. The sources of apparent pixel offset need to be quantified and their effects minimized in order to obtain a reliable terrain-change measurement. The most important apparent pixel offsets to minimize are those that arise from incomplete knowledge of the geometry of the imaging system. The effect of this type of apparent pixel offset is strongest in the distal areas of the examined region, which we discuss in the results section, but can be mitigated by using image pairs that have a maximum of 3° difference in look angle (Schiek, 2004; Van Puymbroeck et al., 2000).

Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) is a push-broom instrument, which contains three separate imaging subsystems that cover the visible and near infrared spectrum (VNIR) with 15 m resolution, the shortwave infrared spectrum (SWIR) with 30 m resolution, and the thermal infrared spectrum (TIR) with 90 m resolution. The VNIR subsystem consists of a solid-state focal plane that records imagery one line at a time in a 60-km-wide swath as the satellite passes over Earth. It includes two Schmidt reflecting-refracting telescopes, both with a focal length of 329 mm. One telescope is nadir (straight down) looking, and the other points 27.6° off-nadir, allowing for stereo viewing and the extraction of digital elevation models (DEMs) (Abrams et al., 2002). All images used in this study are VNIR band 3n (0.78–0.86 µm wavelengths) images from ASTER level 1B scenes. This band was chosen because it has the highest resolution and because a sensor model (Yu et al., 2003), which describes the interior and exterior orientations of the imaging system (e.g., Mikhail et al., 2001), is available. ASTER level 1B scenes were chosen due to their availability in this region, and because they have already been preprocessed to remove some artifacts (i.e., stripe removal, but not parallax) and have been georef-erenced. For more information about ASTER, refer to Abrams et al. (2002) or http://asterweb.jpl.nasa.gov/.

For three locations along the Kunlun fault, pairs of time-separated ASTER images were acquired (Fig. 1; ⁰²). In the first, area A, we examined an image pair that spans the 2 yr time interval from March 2000 to July 2002. This area lies at the eastern end of the Kusai segment (Fig. 1; ⁰²). The images we analyzed for area A show no well-defined fault scarps along the mapped trace of the Kunlun fault. Area B lies in the western part of the Hongshui River segment (Fig. 1; ⁰²), and here we use images that cover a 2 mo interval from 4 November 2001 to 22 December 2001. Area C is located in the western part of the Kusai segment (Fig. 1; ⁰²). The images acquired for area C span more than a year, from October 2000 to December 2001. The images used for both areas B and C include well-defined fault scarps that formed during the 2001 earthquake.

Our analysis required a series of image preprocessing steps. The goal of the preprocessing is to minimize apparent pixel offsets, thus enhancing the true pixel offset due to terrain change. The principal preprocessing step comprised ortho-rectification of the input imagery using the Yu et al. (2003) sensor model, image parameters found in the metadata of each input image, an shuttle radar topography mission (SRTM) DEM of the investigated area, and the Leica Photogrammetry Suite in ERDAS IMAGINE 8.7. See ⁰² for a summary of the image parameters and results of the orthorectification. The ground-control points (GCP) listed in ⁰² were chosen in areas outside of the faulted area, where no change is thought to have occurred. The Leica Photogrammetry Suite performs orthorectification through a series of triangulation calculations between the DEM, image, and ground-control points. After the triangulation procedure, the image is given corrected coordinates, which are corrected for the curvature of the Earth, and which reflect the true geographic location of pixels. For details about the protocol used to complete our orthorectification, we refer the reader to Schiek (2004).

Our method requires that the input image pair is initially co-registered to within ±1 pixel. By linking the images using ENVI 4.0 software, and comparing the same geographic location in corresponding images, we were able to evaluate the quality of the orthorectification. Area A was the only set of images not to pass the ±1 pixel criterion. In this case, we then attempted to co-register them again using ENVI. The result of this co-registration also did not meet the ±1 pixel criterion. This poor co-registration is due to the large pointing-angle difference between the two images ⁰². However, as seen in the results (Fig. 2), this image pair was still sufficient for producing a well-constrained displacement vector direction, but gave an unacceptable vector magnitude.

The factor of 10² was determined empirically and maximized the performance of the water-level deconvolution as an adaptive filter.

A computer program, Fast Fourier Transform Change Detection (FFT_CD), was written to operate on the input images using a moving window scheme, computing a shift for each image segment. Each segment is analyzed separately and an aggregate pixel shift vector, indicating the direction and magnitude of terrain change for that region, is computed. Experimentation shows that the optimum window size is ∼27 × 27 pixels (Schiek, 2004). The results of the FFT_CD analysis are depicted as pixel shift vectors overlain on the “before” image of the input image pair (Fig. 2). Figure 3 gives a direct comparison of the measurements made from this analysis to the measurements from the Lin et al. (2003) field measurements. The displacements and uncertainties we obtained are the average and standard deviation, respectively, of the magnitudes of the vectors in the vector map. The slip azimuth was obtained using rose diagrams (Fig. 4) and depicts the direction of slip along the fault. The azimuthal direction displaying the greatest population of vectors is considered to be the slip azimuth (Fig. 4). This slip direction is not to be confused with the direction of structural trends, such as mole tracks and extensional cracks. The structural trends occur ∼45° from the measured slip azimuth (Fig. 2).

The resulting displacement fields were analyzed with a strain inversion program, StrainSim Pro v. 2.5.7 (Allmendinger, 2005), to determine the strain field as well as the average strain axes responsible for the observed displacement field. We then compared each displacement field and corresponding strain field to the observations made by Lin et al. (2002, 2003) ⁰¹ and the stress analysis by Andronicos et al. (2004a, 2004b) (Fig. 1). Figure 5 displays maps of strain axes calculated from StrainSim Pro 2.5.7 using the FFT_CD vector coordinates. These strain axes were used to map the structural trends displayed in Figure 2.

StrainSim Pro computes finite strain by importing a displacement field comprising the x-y coordinates of the vector starting points and their corresponding magnitudes (Fig. 2). The program creates two grids with this information. One uses the x-y coordinates as the initial position of pixels before deformation (red points in Fig. 5). The other grid is generated by using the magnitude values to compute vector end points, e.g., the final position of pixels after deformation (blue points in Fig. 5). A strain ellipse is computed by building a triangular irregular network (TIN) for both grids (Fig. 5) and comparing the two. After calculating the strain field, the average strain axes are determined by averaging the entire strain field.

All image pairs were processed using the FFT_CD algorithm to yield terrain-change vector maps (Fig. 2). Analysis of the area A image pair resulted in 2.91 ± 0.61 m of average left-lateral slip at an azimuth of 92° (Fig. 4A). The sense of slip and azimuth of slip are consistent with Lin et al.'s (2002, 2003) results (Fig. 2A; ⁰¹). Assessing the accuracy of the amount of slip we obtained with FFT_CD in area A, however, is less straightforward. The closest of Lin et al.'s (2002, 2003) field sites is site 7 (Fig. 1; ⁰¹). They measured several offset features there that had displacements ranging from 3.3 m on a road to 6.8 m on a stream channel ⁰¹. The first measurement of 3.3 m was, nevertheless, similar to our result. The underestimated amount of displacement in area A is likely due to the relatively poor quality of the orthorectification of the input imagery ⁰²; however, the azimuth of the displacement seems well constrained.

The average strain field in area A is compressional, and the strain ellipse has an extensional direction of 140° and a shortening direction of 50° (Figs. 2A and 5A). The compressional strain field measured for this area is compatible with the thrust-type greatest compression axes computed from earthquake data by Andronicos et al. (2004a, 2004b) (Fig. 1). In detail, we find the strain field for area A to be a complex pattern of extension and shortening (Fig. 5A). This complex pattern is compatible with the field observations of Lin et al. (2002, 2003) within the eastern Kusai Lake segment (e.g., site 7 in ⁰¹ and Fig. 1). For example, at site 7, located to the east of area A, Lin et al. (2002, 2003) observed mole tracks, which are indicative of compression. Our analysis predicts compression, and, therefore, features such as mole tracks, in this same area (Figs. 2A and 5A).

The investigation of area B obtained an average left-lateral displacement of 11.39 ± 0.56 m at an azimuth of 93° (Fig. 4B). The strain inversion yielded a compressional strain field with an extensional axis orientation of 125° and a shortening direction of 35° (Figs. 2B and 5B). Our measured displacement was much greater than the 4.6–4.8 m of displacement measured in a gully by Lin et al. (2002, 2003) (site 2 in ⁰¹; Fig. 1). Our overestimation of the amount of slip is mostly likely because the FFT_CD analysis considers a large swath (∼30 km wide; Fig. 2B) surrounding the fault zone and not just the proximal area of the fault where Lin et al. (2002, 2003) conducted their field survey. Thus, our estimate includes both displacement along the fault, as well as displacement at some distance away from the fault (Fig. 2B).

Overall, area B displays transpression along the rupture zone, which is consistent with the transpressional state of stress obtained from the analysis of earthquake data (Fig. 1) (Andronicos et al., 2004a, 2004b). In detail, a complex distribution of compression and extension is observed in our results for area B (Fig. 5B). This complexity is borne out by field observations in the eastern Buka Daban River and western Hongshui River segments (sites 1 and 2 in ⁰¹ and Fig. 1) (Lin et al., 2002, 2003). In these areas, a variety of features are reported, including mole tracks, extensional cracks, and discontinuous shear zones (Lin et al., 2003).

The terrain-change analysis of the area C imagery detects sinistral slip along an azimuth of 92° (Fig. 4C). The resulting strain field is extensional, and the average strain ellipse has an extensional direction of 128° and a shortening direction of 37° (Figs. 2C and 5C). Extension is expected given the extensional cracks observed in the field at site 3 ⁰¹ (Lin et al., 2002, 2003) and the earthquake-derived extensional stress orientations (Fig. 1) (Andronicos et al., 2004a, 2004b). The average amount of slip obtained from our analysis, 6.18 ± 0.6 m, correlates well with the 5.7 m of displacement on an extensional crack measured by Lin et al. (2002, 2003) (site 3 in ⁰²; Fig. 1).

The Fourier terrain-change detection and strain analysis successfully detected the slip from the 2001 Kokoxili earthquake. Here, we interpret our results in terms of large-scale and small-scale processes along the Kunlun fault. The large-scale observations explain the relationship between the relative plate motion of the India-Asia collision zone and the Kunlun fault. The small-scale observations may give insight into the surface rupture characteristics of the Kunlun fault.

The motion of India with respect to Asia in the Tibetan-Himalayan collision zone is at an azimuth of N15°E, which causes SW-NE compression in the Tibetan Plateau (Fig. 1). According to Tapponnier et al. (2001), the material in the plateau is extruded to the southeast along left-lateral faults, such as the Kunlun fault, in response to the collision. The SW-NE compressional direction and the SE extensional direction match the average strain axes we computed in our strain analysis (Fig. 2). We obtained an azimuth range of 35°–50° for the compressional direction and an azimuth range of 125°–140° for the extensional direction (Fig. 2). The computed compressional and extensional axes are similar to the overall stress field of the Tibetan-Himalayan collision zone, but the range of the computed numbers for compression and extension gives insight into stress changes along the Kunlun fault.

Although we detected terrain change consistent with left-lateral strike slip in all portions of the Kunlun fault that we studied, our results from areas A and C suggest that the Kunlun fault becomes more compressional toward the west. This is likely due to the northward bend the Kunlun fault takes, which strikes between 91.5° and 92.5° (Fig. 1) (Lin et al., 2003; Ozacar and Beck, 2004). This change in strike, in the face of the large-scale plate motion, results in variations in the state of stress at different sections along the fault. As the fault bends northward, the strain axes are rotated northward causing the compressional axes to be almost perpendicular to the fault. For instance, in area A, at the easternmost end of the study area where the fault strikes at 91.5°, the extensional direction is 140°, and the compressional direction is 50° (Figs. 1 and 2A). By contrast, area B, located in the westernmost part of the study area, has an extensional direction of 125° and a compressional direction of 35° (Figs. 1 and 2B). However, the fault strikes at 92.5° through this area, and the more northward strike and the counterclockwise rotation of the stress axes yield a more transpressive regime in area B than in area A.

This transition to progressively more compressional stress along the Kunlun fault as a result of the change in strike and the rotation of the stress field is also manifested in the topography. In the eastern part of the study region, where areas A and C are located, the fault lies on the northern edge of a wide basin (Fig. 1). This basin becomes narrower to the west, terminating ∼250 km to the west of area B, where the terrain becomes significantly more mountainous. The change in strike of the Kunlun fault has developed a very large-scale, northward restraining bend, which is evident in the transition from transpressive strike slip at the longitude of area B to the fold-and-thrust belt defined by the north Kunlun thrust system in the Kunlun Shan further west (Figs. 1 and 2; ⁰¹).

In detail, however, our results show complexity in the deformation patterns along the Kunlun fault not attributable to the northward bend alone. For instance, in area C, in between the transpressive zones at A and B, we observe strike slip with a significant transtensional component (Figs. 1, 2, and 5). In addition, along-strike variations in the magnitudes and orientations of the stress axes within areas A, B, and C show similar complexity (Figs. 1, 2, and 5). As with many large strike-slip faults, the surface expression of the Kunlun fault zone likely is composed of a series of en echelon strands (Lin et al., 2004), and compressional and tensional features can develop between and along these en echelon faults. Stepwise en echelon faulting is consistent with the field observations of Lin et al. (2003), who observed features such as mole tracks and extensional cracks ⁰¹. Mole tracks occur between, and perpendicular to, right-stepping en echelon faults and strike N50°W (Figs. 2, 5, and 6) (Lin et al., 2003, 2004). This orientation is consistent with our inferred slip azimuth, and is also what is expected from Andronicos et al.'s (2004a, 2004b) stress orientations and our derived strain orientations. Similarly, extensional cracks occur along, and perpendicular to, the left-stepping en echelon faults and strike N40°W (Figs. 2,5, and 6). From the calculated stress field in Figure 5, the predicted structural trends of mole tracks and extensional tracks is shown in Figures 2 and 5. Mole tracks will exist in areas where the compressional axis (blue) is greater than the extensional axis (red) (Fig. 5), and extensional cracks exist where the compressional axis (blue) is less than the extensional axis (red) (Fig. 5). If the two axes are equal, then strike-slip movement is expected. Figure 6 demonstrates how these features can exist between, and along, the right- and left-stepping en echelon fault strands.

The results from this study employing FFT_ CD are consistent with the results of other studies conducted along the Kunlun fault (i.e., Lin et al., 2002, 2003; Andronicos et al., 2004a, 2004b). The observed pattern of increasing transpression to the west, with small pockets of transtension along the northward bend, is explained by the segmented geometry of the Kunlun fault. Our technique offers a way to observe deformation along a fault in remote areas, such as the Kunlun fault, where other geodetic data may not be available. To achieve higher accuracy and precision in future studies, higher-resolution imagery, with more accurate orbital parameters, is necessary. The greater resolution will enable the FFT_CD procedure to measure smaller amounts of terrain change over smaller time intervals. More precise orbital parameters will greatly diminish uncertainties in the computed terrain-change vector maps, and increase the signal-to-noise ratio between true and apparent terrain change in a scene.

We would like to thank Aaron Velasco, Vladik Kreinovich, Chris Andronicos, Jean-Philippe Avouac, and Sébastien Leprince for their assistance in this project.