3-D surface deformation of the 2000 Usu eruption measured by matching of SAR images

A subpixel level offset estimation technique was applied to SAR (Synthetic Aperture Radar) images of Usu volcano, Japan. Using RADARSAT SAR amplitude images acquired from ascending and descending orbits, we generated a map of three-dimensional displacement vectors showing ground deformation associated with criptdome formation in April, 2000. Moreover, displacement velocities and volume changes were estimated. From April 3-29, inflation was observed with vertical and horizontal displacements each exceeding 20m. During the period of major activity from April 3-5, uplift and N-S spreading velocities reached 3.3m/day and 2.7m/day, respectively, while volume change was +3.4×10⁶m³/day. Although this technique is less precise than other methods and is not very sensitive to local deformations, it provides several key advantages including (1) personnel are not required to enter hazardous areas, (2) it provides wide spatial coverage, and (3) it does not depend on coherence, making it more widely applicable than InSAR.

Usu volcano, located on the southern rim of Lake Toya in the southwestern part of Hokkaido Island, Japan (Figure 4), became active in March 2000. Eruptions continued for several months on the northwest flank of Mt. Usu and formed domes, craters, pyroclastic cones and faults. In this paper, we estimate crustal deformations and volume changes using RADARSAT SAR amplitude images collected during the eruption. SAR interferometry (InSAR) can be a powerful tool for studying volcanic activity [e.g., Massonnet et al., 1995; Rosen et al., 1996; Amelung et al., 2000; Fujiwara et al., 2000]. At Usu, however, the magnitude of uplift (~70m at maximum) and snow coverage prevented effective use of the InSAR technique.

Michel et al. [1999] produced a map of major fault ruptures from the Landers earthquake using SAR amplitude images. Because ground deformation in range was below their data noise level, they measured deformation in azimuth only.

We employed a similar technique at Usu volcano but we measured both the range and azimuth components. Moreover, offset measurements from two different line-of-sight (LOS) directions (i.e., in descending and ascending orbit) enabled us to map three-dimensional components of deformation.

1. Use two SAR images
2. Use amplitude (not phase, not interferometry)
3. Window size of cross-correlation is 128x128 pixel (about 900x630 m).
4. 1-sigma uncertainty is 0.15 pixel.
5. This corresponds to 0.7 to 1.1 m displacement error

Two SAR amplitude images are coregistered using tie-points. In this step, tie-points in deformation areas were excluded in order to avoid the influence of deformation on the coregistration. We selected image pairs with short baselines to avoid stereoscopic effects. Residual offsets in range and azimuth correspond to ground deformation in the LOS and azimuth directions, respectively.

The pixel spacing of the RADARSAT fine mode SAR is ~4.6m in range and ~5.6m in azimuth. Subpixel level estimation of offsets [Tobita et al., 1999] between tie-points enabled us to measure the offsets at a resolution of 1/32 pixel and an accuracy of better than 0.2pixels (i.e., ~1m in this case). Unlike the usual coregistration process for InSAR, in which fewer than 100 tie-points are used, our goal was to measure offsets at many thousands of locations in order to produce a deformation image. A large offset signal (~5 pixels) due to deformation was expected, making it difficult to distinguish good data from bad. To overcome these difficulties, we used a large (128´128 pixel) correlation window [Tobita et al., 1999]. This reduced the occurrence of bad offset data to lower than 5%.

Table 1 shows the SAR data used in this study. Offsets between amplitude images from the same observation mode were measured. LOS azimuth direction (to satellite at surface) is 99° and 261° in F4 and F3 mode, respectively. The standard error of offset measurements is 0.15 pixels, which corresponds to about 0.75m deformation.

Table 1. RADARSAT SAR Data

Data Aqc. Time(JST)	Obs.Mode	Orbit	Incident Angle (deg)
05:37 04/03/2000	F4	descending	44.5
17:37 04/05/2000	F3	ascending	42.5
05:37 04/27/2000	F4	descending	44.5
17:37 04/29/2000	F3	ascending	42.5

Figure 6 shows residual offsets of two images between (a) April 3 and April 27, 2000 and (b) April 5 and April 29, 2000. Blue dots denote tie-points and red lines from the dots indicate the estimated offsets. Large offsets due to surface displacement are observed on Usu's northwest flank. The summit region shows only small offsets. Maximum offset magnitudes are (a) 14.4m (azimuth) and 18.6m (range) for F4 mode; and (b) 10.9m (az.) and 16.3m (range) for F3 mode. Approaching movement is dominant in both images of (a) and (b), implying that most of the ground movement is vertical (upward) rather than horizontal. However, some spreading displacement can be seen in the N-S direction in both images as well.

Offsets measured from just one direction provide range and azimuth components of displacements only. Offsets measured from two directions, however, provide four independent estimates of 3-D displacement projections. From these four components a 3-D vector deformation field (up, north, east) can be synthesized.

Measured offsets in azimuth and range, dA and dR, are given by inner products between the displacement vector M and unit vectors as follows:

dA1 = M1 o T1
dR1 = M1 o S1
dA2 = M2 o T2
dR2 = M2 o S2 ------ (1)

where the subscripts, 1 and 2, denote variables from two directions and/or from two periods. T and S are unit vectors representing measurement directions, azimuth (the satellite track angle) and range (from target to satellite). There are 4 observations (dA1, dR1, dA2, and dR2), and 6 unknowns for each image point. To further constrain the solution, we introduce the following assumptions:
M2 = k M1 ------ (2)
where k is a single proportionality coefficient over the image. This amounts to the assumption that the directions of displacement do not change in time and the ratio of displacement magnitudes of period 1 to those of period 2 is constant (= k) in space. Moreover, the displacement from April 27-29 is assumed to be zero (see Figure 7).

We define M1 and M2 to be displacement vectors from April 5-29 (F3) and from April 3-27 (F4), respectively. The velocity vector of displacement from April 3-5 (2.5 days) is then

V=(M2-M1)/2.5=(k-1)M1/2.5 (m/day) ------ (3)

At each point in the image, therefore, the number of unknowns decreases to 4, but one of the unknowns, k, is constant over an image. Thus a whole image with n measured points has 4n observations and 3n+1 unknowns. Unknowns are estimated using a least square adjustment.

In this process, corresponding points between images 1 and 2 must be identified. To accomplish this, amplitude images were simulated using GSI's 50-m spacing digital elevation model (DEM) and then coregistered to the real images. The latitude and longitude of pixels in images 1 and 2 are given through the DEM. Thus, points that have the same coordinates can be extracted from the two images.

The coefficient k was estimated to be 1.6. This means that the deformation magnitude from April 3-5 (i.e., 2.5 days) is equal to 60% of that from April 5-29 (i.e., 24 days). This is in good agreement with the ratio obtained from GPS measurements: 57% for north component and 63% for east component. GPS displacements are shown in Figure 7.
The vector deformation fields presented in this section are estimates of M2 (i.e., from April 3-27).

Figure 8 shows the 2-dimensional distribution of the horizontal components of ground deformation. The main area of deformation is centered at about 140°48.4'E, 42°33.1'N. A smaller area of deformation lies northeast of the main area, centered near 140°48.8'E, 42°33.3'N. This bimodal distribution is consistent with the distribution of craters identified by photo-interpretation [Sekiguchi and Takazawa, 2000]. In the main area, a N-S spreading displacement is dominant but there is some E-W spreading as well. Relative displacements of the N-S spreading exceed 20m at the center of the main area. An area with horizontal components greater than 1m in magnitude is also observed to the southwest of the main area, centered near 140°47.8'E, 42°32.6'N.

Figure 9 shows the two-dimensional distribution of vertical displacements. There is no subsidence, only uplift. This figure suggests a bimodal distribution as well. The main area of uplift is at the same position as the main area of horizontal motion. Maximum uplift is greater than 22m.

From equation (3) and k=1.6, the deformation velocity vector V for April 3-5 are estimated to be 0.15 times the displacement vector M2 for April 3-27. Uplift and N-S spreading velocities therefore reach 3.3m/day and 2.7m/day, respectively, for this period in the main area.

Having the 3-D deformation vector field allows us to calculate volume changes due to uplift. The volume change from April 3-27 is +22.9×10⁶m³, while the volume change from April 5-29 is +14.3×10⁶m³. From these, the average rate of volume change from April 3-5 is +3.4×10⁶m³/day.

We generated a map of three-dimensional displacement vectors and estimated displacement velocities and volume changes. From April 3-29, 2000, inflation was observed with vertical and horizontal displacements each exceeding 20m. The rates of displacements and volume change were 3.3m/day (uplift), 2.7m/day (N-S spreading) and +3.4×10⁶m³/day during the period of major activity from April 3-5. The surface deformation measured by this study is due to inflation of subsurface magma.

It is important to note that our technique cannot detect local deformation at scales less than the correlation window, in this case 900X630m. Heliborne laser altimeter data collected by the Public Works Research Institute of Japan recorded a peak uplift of 69m from March 31 to April 26, 2000. To check quantitative consistency between the altimeter and our SAR technique, we compared the two datasets after correcting for resolution and period differences.

We averaged the 1-m spacing digital altimeter data (i.e., height change) in a rectangle window. Window sizes of 900X630m, 450X315m, and 225X158m were tried, corresponding to correlation windows of 128X128, 64X64, and 32X32 pixels in a RADARSAT SAR image, respectively. The 69m uplift peak (narrow with steep slopes) is centered in these simulations. The altimeter data were also reduced by a factor of 0.63, derived from the north components of the GPS measurements in Figure 7, because the altimeter data cover a longer period of deformation.

Table 2. Comparison between SAR and Altimeter Uplift after Correcting for Resolution (window size) and Period Differences in a Local Uplift Area.

SAR image offset		Laser altimeter		Difference
correlation win-size (pix)	uplift(m)	averaging win-size (m)	uplift(m)	(m)
128X128	21.7	900X630	20.7	1.0
64X64	28.9	450X315	30.2	-1.3
32X32	37.4	225X158	35.8	1.6

Comparison of heliborne laser altimeter and SAR offsets shows good agreement to within 2m (Table 2). This implies that a measured displacement in our technique is an average in a correlation window (900X630m in size). Smaller correlation windows are more sensitive to local deformations but include more bad data, thereby requiring greater effort to distinguish good data from bad. This is a trade-off.

The advantages of our technique are: (1) personnel are not required to enter hazardous areas, (2) it provides wide spatial coverage, (3) it is possible to obtain three-dimensional components of displacement, (4) it is possible to estimate displacement velocities, and (5) it does not depend on coherence, making it more widely applicable than InSAR. The disadvantages are: (1) it is less precise than InSAR and GPS, (2) it provides less temporal resolution than GPS, and (3) it is difficult to perform real-time monitoring.

Acknowledgments. We are grateful to S. Matsuzaka and D. A. Phillips for their comments and advice. We thank M. Shimada for providing information of RADARSAT data and advice. SAR data were provided by National Space Development Agency of Japan. RADARSAT R-DATA (c) Canadian Space Agency (2000). Public Works Research Institute and Hokkaido Government kindly provided us digital data of heliborne laser altimeter of Mt. Usu.

• Tobita, M., Mak. Murakami, H. Nakagawa, H. Yarai, S. Fujiwara, P. A. Rosen (2001), 3-D surface deformation of the 2000 Usu eruption measured by matching of SAR images, Geophys. Res. Lett. Vol. 28 , No. 22 , 4291-4294.

• Amelung, F., S. Jonsson, H. Zebker, and P. Segall, Widespread uplift and 'trapdoor' faulting on Galapagos volcanoes observed with radar interferometry, Nature,407, 993-996, 2000.
• Fujiwara, S., T. Nishimura, M. Murakami, H. Nakagawa, M. Tobita, and P. A. Rosen, 2.5-D surface deformation of M6.1 earthquake near Mt Iwate detected by SAR interferometry, Geophys. Res. Lett., 27, 2049-2052, 2000.
• Geological Survey of Hokkaido, Usu eruption in 2000 (in Japanese), GSH News, 2000, 5, vol. 16, (ISSN1345-1138), 4 pp., 2000.
• Massonnet, D., P. Briole, and A. Arnaud, Deflation of Mount Etna monitored by spaceborne radar interferometry, Nature, 375, 567-570, 1995.
• Michel, R., J-P Avouac, and J. Taboury, Measuring ground displacements from SAR amplitude images : application to the Landers earthquake, Geophys. Res. Lett., 26, 875-878, 1999.
• Rosen, P. A., S. Hensley, H.A. Zebker, F.H. Webb, and E. J. Fielding, Surface deformation and coherence measurements of Kilauea Volcano, Hawaii, from SIR-C radar interferometry, J. Geophys. Res., 101, 23109-23125, 1996.
• Sekiguchi, T., and S. Takazawa, Repeat photointerpretation of evolving eruptive activity around Usu volcano (in Japanese), Abstract of 2000 Meeting of Japan Cartographers Assoc., 32-33, 2000.
• Tobita, M., S. Fujiwara, M. Murakami, H. Nakagawa, and P. A. Rosen, Accurate offset estimation between two SLC images for SAR interferometry (in Japanese), J. Geodetic Soc. of Japan, 45, 297-314, 1999.