Imageodesy on MPI & GRID for Co-seismic Shift Study Using Satellite Optical Imagery

1. Knowledge discovery from massive data processing for earthquake study Imageodesy on MPI & GRID for Co-seismic Shift Study Using Satellite Optical Imagery Jian Guo Liu and Jinming Ma Department of Earth Science and Engineering Imperial College LondonSouth Kensington campus, London SW7 2AZ, UK j.g.liu@imperial.ac.uk

2. Outline of the presentation • Introduction: data mining software development for geohazard study • Imageodesy: the principle and implementation on MPI and GRID • Case study: Co-seismic shift of the Ms 8.1 Kunlun earthquake • Conclusions

3. Acknowledgement This research is part of DiscoveryNet project (GR/R67750/01) supported by EPSRC e-science pilot project grant. Computing Centre of Imperial College London provided parallel processing facilities and technical support. MIT Phase Correlation website provided free access and technical support for software development. ER Mapper image processing software has been used for data visualisation and analysis. Xinjiang Bureau of Seismology is acknowledged for providing field photos and some reference materials of Kunlun earthquake.

4. 1. Introduction Geohazard study in DiscoveryNet project Geohazard monitoring and assessment is one of the major application areas of DiscoveryNet project aiming to scientific knowledge discovery. Remote sensing A typical characteristic of geohazards is that movement and displacement will be produced. For instance, an earthquake produces co-seismic displacement while a landslide is characterised by slope failure and mess movement. Remotely sensed imagery data enable detection and measurement of these changes and thus provide vital information for hazard assessment and prevention. Algorithm and software development Fast imageodesy algorithms and software for high accuracy change detection have been developed and integrated into DiscoveryNet workbench. The massive data processing is possible only with advanced algorithms using powerful MPI/GRID computing.

5. 2. Imageodesy: the principle The ‘Imageodesy’ technique (Crippen 1992; Crippen & Blom 1996) is capable of measuring the horizontal shift of image features, at a sub-pixel level accuracy, through normalized cross-correlation (NCC) between pre- and post-event images. This technique is a complementary to the well-established radar interferometry technique that is sensitive to vertical deformation. The major technical challenge of high accuracy imageodesy is the huge demand on computing to handle with massive data processing. We have developed software to implement the imageodesy on MPI parallel processor based on the conventional template normalised cross-correlation (NCC) algorithm. Furthermore we developed new algorithm and software based on phase correlation algorithm.

6. 2. Imageodesy: FNCC algorithm The NCC is defined as: The Fast NCC algorithm (FNCC) reduces large quantity of repeated operations and further improvement including conditional jumps and smart sampling avoids unnecessary operations over homogeneous image features. Thus the processing can be speed up by 5-10 times depending on the size of processing windows and the image characteristics.

7. Image “before” Image “after” Read dataset Read dataset Set searching window Set calculation window Move calculation window Maximum correlation N Y Shift-Y Shift-X Correlation coefficient 2. Imageodesy: FNCC MPI implementation With an improved FNCC algorithm, operating on a MPI UNIX parallel computer with 24 processors, it takes 10 hours to complete imageodesy processing for one pair of cross-event Landsat-7 ETM+ Pan imagery data. The image size is 3.75GB after interpolating to 3m pixel size.

8. 2. Imageodesy: FNCC GRID implementation • FNCC is a neighbourhood processing, each image line (the minimal data unit) must carry its neighbour image lines with it in order to conduct the processing. This wipe off 50% efficiency of FNCC and increases data communication by m times, where m is the search window size. • The experiments on GRID using a small image (512512) completed much slower than local processing using a single PC (2GHz processor). Submission larger images of a few thousands lines and columns to the GRID simply blocked the processing pipe line and failed to complete the task. • The current status of GRID is not sufficient for the massive neighbourhood processing of FNCC imageodesy. The future of GRID for dealing with the type of processing of imageodesy lies on very fast high throughput network.

9. Image “before” Image “after” Read Dataset Read Dataset Hamming Windowing Hamming Windowing FFTW FFTW Phase Correlation Inverse FFTW Delta Y Delta X Correlation coefficient 2. Imageodesy: Phase Correlation algorithm We have implemented phase correlation algorithm for imageodesy operating on single UNIX and PC. By transforming the image data within a matching window into frequency domain via FFT, the phase correlation can pinpoint the best matching position directly as the peak of the overlap between the frequency distribution of the two images, without the time consuming searching. Phase Correlation Imageodesy algorithm scheme

10. 3. Case study: Co-seismic shift of Ms 8.1 Kunlun earthquake Kunlun Earthquake On 14 November 2001, at 09:26:18 UTC, an Ms 8.1 earthquake occurred in the East Kunlun Mountains, along the Kusai Lake segment of Kunlun fault. A 400 km surface rupture zone of E-W to WNW-ESE orientation was produced and, according to the field observations of Chinese scientists, the fault displacement was as large as 16.3 m (Lin et al. 2002, 2003). Remote sensing study Landsat TM and SPOT images have already been used to locate and map the visible surface rupture features (Fu & Lin 2003). SAR interferometry (InSAR) would be an ideal technique to reveal the stress field of the earthquake and to provide two-dimensional quantitative measurements of the fault movement. The lack of high quality cross-event ERS SAR fringe pairs for this region and in this particular time unfortunately hindered the use of interferometry.

11. 3. Case study data: ETM+ imagery Landsat ETM+ imagery data was chosen because of its large coverage, improved 15 m resolution panchromatic band, and the availability of suitable pre- and post-earthquake, almost cloud-free scenes. The output images of imageodesy are: X-shift, Y-shift and R (the NCC coefficient).

12. 3. Case study results:the left-lateral movement of Kunlun fault

13. 3. Case study results:Kusai Lake scene Histograms of smoothing filtered X-shift image of KL scene with 0.7 correlation threshold. Left: The whole scene: the high peak on the left is at 0.7 and the shoulder on the right is centred at 2.5. Middle: The south side of the Kunlun fault zone, the peak is at 2.5. Right: The north side of the Kunlun fault zone, the peak is at 0.7. The range of the left-lateral shift is 1.5~8.1 m, while shift is most commonly 5.4 m. The maximum net left-lateral displacement can be as great as 13 m.

14. 3. Case study results:Kusai Lake scene The co-seismic shift vectors of the Kusai Lake area overlaid on the post earthquake ETM+ Pan image. The vectors were derived from X and Y-shift images, with 371371 window averaging, 20% cut-off for elimination of extreme values, and a 0.8 NCC coefficient criterion.

15. 3. Case study results:Kusai Lake scene The co-seismic shift vectors of the Kusai Lake area overlaid on the post earthquake ETM+ Pan image. The vectors were derived with X-shift compensated by -3.6 m.

16. The earthquake fault zone from space

17. 3. Case study results:Buka Daban scene Histograms of smoothing filtered X-shift image of BD scene with 0.7 correlation threshold. Left: The whole scene: the high peak on the left is at 0.35 and the low peak on the right is centred at 1.7. Middle: The south side of the southern branch, the peak is at 1.7. Right: The north side of the northern branch, the peak is at 0.3. The range of the left-lateral shift is 1.0 to 8.2 m, while the most representative figure is 4.2 m.

18. 3. Case study summary Kusai Lake scene: Fault is consistently in W-E to WNW-ESEdirection. Left-lateral shift range 1.5 m~8.1 m, average 4.8 m, the maximum 13 m. Buka Daban scene:With the splayed nature of the fault in this section, the displacement patterns become complicated and stepwise. The average left-lateral shift over the broad fault zone is 4.6 m and ranges from 1.0 m to 8.2 m. Both sides of the faults moved toward the east. the south side of the fault has been displaced significantly to the right (east) relative to the largely stable, or slightly right-shifting, northern block. The relative movement of the fault is left-lateral and the south side of the fault is the active block.

19. 4. Conclusions • As an essential function of remote sensing data mining for geohazard study in DiscoveryNet project, FNCC imageodesy technique has been implemented on the workbench operating on MPI. • The FNCC implementation onGRID yields a disappointing performance because the demand for data communication increases dramatically when the neighbourhood processing of imageodesy is distributed to many nodes. • The phase correlation based algorithm is currently not operational on MPI or GRID processing mode. It is only more efficient when the forward and inverse FFT operations in phase correlation take less time than searching in FNCC. • Our imageodesy results present the first regional 2-D picture of the co-seismic displacement of Kunlun fault as the result of the Ms 8.1 Kunlun earthquake, in a vast area, of about 320 km W-E and 180 km N-S. It is an important scientific knowledge discovery.

20. Thank You! Picture provided by Xinjiang Bureau of Seismology