Interactive comment on “ Estimation of three-dimensional crustal movements in the 2011 Tohoku-oki , Japan earthquake from TerraSAR-X intensity images ”

General comments: Liu et al. first use TSX images to derive 2-D offset fields with a previous proposed method, and then construct 3-D surface displacements in the 2011 Tohoku-oki, Japan earthquake basing on the three sets of detected 2-D movements. The results compared with the GPS observations indicate that the application is acceptable. Although the methods used in this study are not new, the work is very well. I consider the manuscript is good and appropriate for the NHESSD.


Introduction
The M w 9.0 Tohoku-Oki earthquake occurred on 11 March 2011 off (oki in Japanese) the Pacific coast of northeastern (Tohoku) Japan and caused large tsunamis and widespread devastation.According to the GPS observation network system (GEONET) of the Geospatial Information Authority (GSI) in Japan, significant crustal movements were observed over a wide area (Ozawa et al., 2011).Because information about crustal movements is important for estimating the fault model after an earthquake and for evaluating seismic hazards, crustal movements are monitored daily at GPS ground stations around the world.Hashimoto (2013) summarized the crustal deformation before, during, and after the 2011 Tohoku-Oki earthquake using data sets from GEONET and sea-floor GPS-acoustic sites that are operated by the Japan Coast Guard.However, GPS observation systems have not been established in all parts of the world, especially in developing countries.Thus, estimating crustal movements from earth observation satellites is considered to be efficient and important.
Interferometric synthetic-aperture radar (InSAR) is a useful method for detecting spatially coherent ground information (Zebker, 2000).Many researchers have detected daily and coseismic displacements based on InSAR (Stramondo et al., 2002;Prats et al., 2008;Wicks et al., 2013).InSAR can detect displacements at the centimetre level despite its metrelevel image (pixel) resolution, but displacements can only be obtained in the slant-range direction.Another drawback of InSAR is that the method is applicable only for relatively small coherent displacements.The cross-correlation (pixeloffset) method is an alternative for measuring the relative two-dimensional (2-D) displacements between two images at a pixel-to-sub-pixel scale, such as in detecting coseismic deformation (Michel et al., 1999;Tobita et al., 2006;Leprince et al., 2007).Actual 3-D surface displacements can be obtained by these methods using SAR image pairs that are captured from ascending and descending orbits (Fujiwara et al., 2000;Fialko et al., 2001;Michele et al., 2010).In contrast, a combination of two InSAR analyses results in estimates of two-and-a-half-dimensional (2.5-D) displacements.In in- versions to estimate fault-slip distributions, GPS records are required in addition to InSAR analyses (Zhang et al., 2007;Beavan et al., 2011;Jiang et al., 2013) to estimate 3-D displacements.
The Tohoku-Oki earthquake produced significant crustal movements over a large area and dramatic surface changes due to mainly repeated tsunamis.The ground displacements have been calculated using both InSAR and pixel-offset techniques using TerraSAR-X (TSX) data (Martinez et al., 2012).We have also employed an improved pixel-offset method for a descending pair of high-accuracy georeferenced TSX data to estimate the absolute ground displacements (Liu and Yamazaki, 2013).We calculated the cross-correlation between the pre-and post-event TSX images for individual non-damaged buildings because the coastal lowlands in the Sendai area were severely affected by the tsunamis (Liu et al., 2013); thus, cross-correlations were difficult to obtain for the ground surface.The maximum error between the GPS records and the detected 2-D displacements from the TSX imagery was less than 0.5 m in spite of the pixel size (1.25 m), which demonstrates the high accuracy of the proposed technique.
In this paper, three pairs of pre-and post-event TSX intensity images for the 2011 Tohoku-Oki earthquake taken in ascending and descending paths were used to detect 2-D crustal movements using the improved pixel-offset method based on the displacements of non-damaged buildings.The actual 3-D movements were then estimated from the three sets of detected results from the different acquisition conditions.The accuracy of the 3-D detection results was demonstrated and compared with the results from three GEONET station records.

Study area and image data used
The study area focuses on Miyagi Prefecture, which is located in the coastal zone of Tohoku, Japan, and is shown in Fig. 1a.This area was severely affected by the 2011 Tohoku-Oki earthquake tsunami.Six temporal TSX images, which were taken before and after the earthquake along three differ- ent paths, were used to estimate the 3-D crustal movements.The acquisition conditions for pairs A, B, and C are shown in Table 1.The descending Pair C was used in the previous study (Liu and Yamazaki, 2013).
All of the images were captured with HH polarization in the StripMap mode.Due to the different incidence angles, the resolutions of the ascending pair in both the azimuth and range directions were approximately 3.0 m, whereas those of the descending pair were 3.0 m in the azimuth direction and 3.5 m in the range direction.We used the standard Enhanced Ellipsoid Corrected (EEC) products at processing level 1B, which were projected to a WGS 84 reference ellipsoid with a resampled square pixel size of 1.25 m.The pixel-offset method is usually applied to slant-range images.However, geocoded products were used in this study and thus the method was applied to both the azimuth and ground range directions.
Because the six images cover somewhat different areas, the target area was defined as the region where they overlapped, which is marked by the black rectangle in Fig. 1a.Two pre-processing steps were applied to the images before extracting the crustal movements.First, the six TSX images were transformed to a sigma naught (σ 0 ) value, which represents the radar reflectivity per unit area in the ground range.Then, an enhanced Lee filter was applied to the original SAR images to reduce speckle noise.The application of an adaptive filter improves the correlation coefficient of two SAR images.To minimize any loss of information from the intensity images, the window size of the filter was set as 3 × 3 pixels.The pixel localization corrected by the GPS orbit determination was used directly; thus, the accuracy of the proposed method depends on the geolocation accuracy of the products (Breit et al., 2010).Based on our previous study, the accuracy of the 2-D detection for StripMap images is approximately 0.3 m.Several recent studies demonstrated that the geolocation accuracy of TSX images can be improved to a few centimetres after the correction of atmospheric refractivity and Earth surface dynamics (Eineder et al., 2011;Cong et al., 2012).
After pre-processing, three colour-composite images of the ascending and descending pairs were generated as shown in Fig. 1b-d, where the pre-event image was assigned to the green and blue bands, and the post-event image was assigned to the red band.The common frame shows the target area, which includes three GEONET stations, Rifu, Natori, and Watari, that are located from north to south.

GEONET and field survey
The GSI has established approximately 1300 GPS ground control stations, called GEONET (GNSS Earth Observation Network System), throughout Japan.Real-time observations of crustal movements are made possible by continuous observations at these stations, where radio waves from the GPS satellites are constantly received.Thus, ground movements in Japan have been continuously monitored by GEONET since 1993 (Sagiya, 2004).There are three types of GEONET GPS receiving towers; examples of two types are shown in the photographs in Fig. 2, which were taken during the authors' field survey in January 2012.These ground control stations were located in the study area of this paper.
More than 5.3 m of horizontal movement and 1.2 m of vertical movement were observed by the GEONET system after March, due to many aftershocks, was more than 20 cm.Although GEONET is one of the densest GPS networks in the world, the distance between neighbouring stations is more than 20 km.Hence, it is difficult to estimate the distribution of crustal movement in detail using only GPS data.In addition, several GEONET stations stopped operating due to power outages, and some were structurally damaged by shaking and tsunami surges.
Our field surveys showed that the Natori GEONET station, which is located at Yuriage Junior High School, was hit by the tsunami.Signs of the tsunami surge were observed at the antenna tower and nearby buildings.This GEONET station stopped operating just after the mainshock and was restored more than 1 month afterwards on 18 April 2011.Inundation of the Natori station can be confirmed from the aerial photos shown in Fig. 4a-b, which were taken by GSI.Significant surface changes are also observed in the colour composite of the pre-and post-event TSX intensity images shown in Fig. 4c. Figure 4d shows a panoramic photo taken during our field survey on 4 August 2011, in which a ship carried by the tsunami is still on the ground 5 months after the disaster.The Watari GEONET station, which was not affected by the tsunami, also stopped operating from 12 to 15 March 2011, most likely because of a power outage.The normal crosscorrelation was applied to this area as an experiment.The size of the search window was set as 10 pixels larger than the moving window in each direction.Although a correlation coefficient matrix between −1 and 1 was obtained, the maximum value due to tsunami damage was −0.04.Thus, it is difficult to evaluate whether the results obtained from the normal method are reliable for detecting displacements.

Two-dimensional displacement detection
The 2-D crustal movements within the target area were extracted by the previously proposed pixel-offset method (Liu and Yamazaki, 2013).The entire target area was first divided into a 2.5 × 2.5 km 2 mesh.In each sub-area (2.5 × 2.5 km 2 ), solid (reinforced concrete or steel) buildings with strong radar-reflection areas greater than 150 m 2 were extracted by a simple segmentation approach.As shown in Fig. 5a, the wall of a building in the sensor direction shows a strong reflection (called backscatter) due to the side-looking nature of SAR.Thus, we extract this type of high-backscatter area as a building object.The threshold value of the backscatter coefficient used in the segmentation was set at −2.5 dB, so pixels with coefficients greater than −2.5 dB were grouped as a building object.Several natural elements or other structures were incorrectly extracted as "building objects" due to their high reflections.However, these objects also have unique reflection patterns that are similar to buildings.Thus, they were regarded as "building objects" in the displacement detection.Then, potentially non-damaged buildings were extracted by comparing the building locations in the two images.If a solid building was present in the post-event image within 5 pixels of the location of the target building in the pre-event image, the target building was regarded as potentially non-damaged.The displacement of a potentially non-damaged building between the two TSX intensity images was calculated by an area correlation method.Due to the large extent of inundation, the extraction of potentially non-damaged buildings could avoid calculations of obviously changed areas and save processing time.To improve the location accuracy, the TSX images surrounding unchanged structures were resampled to 0.25 m pixel −1 by cubic convolution to one-fifth of the original pixel size.Thus, changes of building shapes could be detected at a sub-pixel level.
A small target area was selected from the pre-event intensity image at the location of a potentially unchanged building object.A larger search area surrounding the target area was then selected in the post-event intensity image.The target area was overlain by the search area and then shifted within it.The cross-correlation coefficient was calculated for each shift.The location of the maximum value for different shifts was defined as the final offset of the target building in the search area as shown in Fig. 5b.The detected displacement was considered to be valid only if the correlation coefficient of a non-damaged building between the two TSX images was greater than 0.8; the building is then considered to be nondamaged.The crustal movements of the entire area were obtained in the same way.Finally, the average movement of the non-damaged buildings in that sub-area was determined and assumed to represent the crustal movement.To ensure the reliability of the results, only sub-areas that contained more than five building displacements were used.
This building-based method was applied to the three pairs of TSX intensity images.Because the displacements were calculated from the building locations, the proposed method could not be applied to mountainous areas that did not contain buildings.The detected 2-D movements from the three image pairs are shown with vectors in Fig. 6a.Displacements were detected by the proposed method for most of the subareas.Because of the different incidence angles for paths B and C, the detected displacements from path B are smaller than those from path C. The different acquisition conditions also influence the layover shapes of the same building in the different pairs.However, the detections were performed on the image pair with the same conditions.Thus, the geometrical characteristics of the structures do not affect the accuracy of the detection.The displacements of the three subareas that contain the GEONET stations were all detected; their values are shown in Table 2.The detected displacements show that the largest movement was observed around the Rifu GEONET station, which was the closest station to the focal region.

Three-dimensional (3-D) displacement detection
The surface displacement is a vector in 3-D space with three components, D E , D N , and D Z , in the east, north, and vertical directions respectively.A SAR sensor observes these displacements in the slant-range and azimuth directions.Because the EEC products used in this study have been geocoded by the cylindrical map projection, the proposed pixel-offset method can detect the 2-D movements in the east and north directions, which were transformed from the slantrange and azimuth directions.The vertical displacement is transformed to one in the slant-range direction and then decomposed into the east and north directions in the ground range.Figure 7 shows the schematic layout of an actual 3-D crustal movement and its 2-D shift on a ground range SAR image.The relationship between a crustal movement (D) in the east, north, and vertical directions and its shift (M) in the ground-range (R) and azimuth (A) directions in an EEC product is represented by Eq. ( 1): where α is the heading angle clockwise from north, and θ is the SAR incidence angle.By rotating the relationship in the east and north directions, we obtain the relationship between the 2-and 3-D displacements as (2) The detected 2-D displacement from one pair of TSX images is not sufficient to estimate the 3-D movement because the two equations contain three unknowns (D E , D N , and D Z ).However, two pairs of different observation sets can form four individual equations; thus, a 3-D estimation is possible from them.Although there were some time differences between the three pre-event images, no significant geodetic movement occurred during these periods.Thus, the locations in the pre-event images were considered to be unchanged.There was a difference of 6 days between the three post-event images; however, GPS records indicate that the crustal movements in this period were small enough to be ignored.Hence, the displacements between these three image pairs can also be considered to be negligible; thus, it is possible to estimate the actual 3-D crustal movements from them.An additional factor would be introduced if significant displacements occurred due to the aftershocks.Because the number of equations is greater than the number of unknown displacements in the three directions, a regression analysis was applied to calculate the most suitable displacement values.The estimations using two and three pairs results was carried out to clarify the relationship between the combination condition and the accuracy.

Estimation using two date pairs
A 3-D displacement can be estimated from two sets of 2-D detection results by regression using Eq. ( 3): where D i is an independent variable (the actual geodetic movements in the east, north, and vertical directions) obtained at a point, and M ij is a dependent variable (the shifts in the SAR image in the east and north directions) obtained from the 2-D analysis for the j th pre-and post-event acquisition pair.Three combinations exist when two pairs of 2-D estimation results are used: one with ascending path A and descending path B, one with path A and descending path C, and one with paths B and C. The regression analysis was applied to Table 3.Comparison of the estimated 3-D displacements from the regression analysis with those from GPS records.Three two-set combinations were evaluated using a 2.5 km mesh, and two threeset combinations were assessed both by the mesh and individual buildings (units: m).Bold numbers represent the minimum value of RMS errors.The horizontal displacements from these combinations were stable, which means the displacements of neighbouring grids show similar amplitudes and directions.The vertical displacements from the set with ascending path A and descending path B were more stable than those from the other two sets.The vertical displacements estimated from descending paths B and C were much greater than those from the combinations of the ascending and descending paths.Thus, the vertical arrows in Fig. 6d are plotted at twice the scale of those in Fig. 6b-c.The set that used ascending path A and descending path C resulted in the smallest vertical displacements.Many of the vertical displacements from this set were directed upward, although the GPS data indicated subsidence.These inconsistent results were caused by errors in the 2-D detections.For example, the horizontal displacement detected around the Rifu station from path A is smaller than the observed value, whereas the horizontal displacement from path C is greater than the observed value.Thus, the estimated vertical displacement using these two pairs shows an upward movement despite the subsidence that was observed by the GPS station.
The estimated results for the sub-areas that include the GPS ground stations from the two different path pairs are shown in Table 3 along with the GPS records for comparison.The GPS records obtained on 1 April 2011, which is the midpoint of the post-event acquisitions, were used as a reference.The root mean square (RMS) errors were calculated for the differences between the analysis results and the GPS records for the three directions.The RMS errors for the first set were smallest in the sub-areas of the Rifu and Watari stations.The RMS error for the second set was smallest in the sub-area that contains the Natori station.The RMS errors for the third set were large in all of the sub-areas.These results confirmed that the combinations of ascending and descending paths yield better results than only the descending paths.In addition, the combination of pairs A and C with different incidence angles (35 and 21 • ) gave better results than using pairs A and C, which had similar incidence angles (35 and 37 • ).

Estimation using three date pairs
We also performed an estimation using three sets of detection results by a regression analysis using Eq. ( 4), which represents six equations with three unknowns.Because not all of the sub-areas contain results from all three pairs, the subareas that have results from at least two pairs were calculated.Based on the comparisons shown in Table 3, the combinations of ascending and descending paths had higher accuracies than those with two descending paths (B and C).Thus, only estimation results from the ascending and descending combinations were used hereafter in cases with two pairs.
The estimated results using the three sets are shown in Fig. 8a.The specific pairs used for the estimations in each sub-area are shown in Fig. 8b.The displacement amplitudes in the horizontal and vertical directions over the entire target area were determined by Kriging interpolation as shown in Fig. 8c and d.The 3-D movements were determined for 73 sub-areas (approximately 456 km 2 ).The displacements for 52 subareas were estimated from the three sets of estimation results, and the other 21 sub-areas had the same values as those shown in Fig. 6.The horizontal movements ranged from 2.62 to 3.55 m, while the vertical movements ranged from −0.65 to 0.19 m. Figure 8 confirms that the largest movements occurred near the northeastern coast, and the values decreased to the south.The detected movements point increas- ingly eastward from north to south.These trends matched the observed GPS data from GEONET (Ozawa et al., 2011;Hashimoto, 2013).The subsidence was also greatest along the coast and decreased westward.This trend is also consistent with the observed GPS data.However, the results were inconsistent in several sub-areas.A comparison of the estimated results with the number of pairs used confirmed that the results for the sub-areas that were estimated from three image pairs were more stable than those estimated from only two pairs.
The comparison of the 3-D results near the three GPS ground stations using three sets of estimation results are also shown in Table 3.The maximum difference between our results and the GPS measurements was approximately 0.25 m in the north direction at Natori station.The RMS errors for the three sub-areas in 3-D space were all less than 0.15 m.Compared with the displacements that were estimated using two sets of results, the RMS errors at Rifu and Watari stations are slightly larger, but the estimation using three sets resulted in more stable displacements because more equations were employed in the regression.However, the combination pairs A and B is enough to estimate the 3-D displacements with high accuracy.

Conclusions
In this study, the 3-D crustal movements caused by the 2011 Tohoku-Oki earthquake were estimated from 2-D displacements that were detected from TerraSAR-X intensity image pairs from three different (one ascending and two descending) paths.First, 2-D crustal movements in Miyagi Prefecture were obtained from the three sets of pre-and postevent SAR image pairs using a building-based 2-D crosscorrelation method.The 3-D movements in a 2.5 km grid were then estimated from two or three sets of 2-D results based on the SAR observation geometry.Because three unknown displacements (independent variables) exist in the three directions, which is fewer than the number of equations for dependent variables, a regression analysis was used for the estimation.
The different combinations of satellite paths were discussed, and the results were evaluated using the GPS data from three GEONET stations.The combination of the ascending path and the descending path with a 21 • incidence angle had the smallest RMS error.Overall, the combination of all three paths provided the most stable and reliable estimates for the entire target area.Most of the estimated results in the horizontal direction were stable and showed similar trends to the GPS data.However, the results in the vertical direction were unstable in several sub-areas.Compared with the observation data from the three GEONET stations, the errors in the horizontal and vertical directions were less than 0.3 m.
The proposed method can estimate 3-D coseismic displacements at a sub-pixel level using only SAR images.The results are useful in seismic source modelling and the method does not require GPS data.Since the proposed method focuses on building objects, it works only in urban areas and requires high-resolution SAR images.The standard crosscorrelation method can complement the proposed method for the areas without buildings.Although the results obtained in this study are promising, we will test the proposed method for other earthquake events with smaller surface displacements to further validate its applicability.

Figure 1 .
Figure 1.(a) Study area with the six TerraSAR-X images in Miyagi Prefecture.(b-d) Colour composites of the three pairs of pre-and post-event TSX images in paths A, B, and C.

Figure 2 .
Figure 2. Photos of different types of GNSS-based ground control stations taken during the field surveys of the study area.

Figure 3 .
Figure 3. GPS-observed displacement vectors at 18 ground stations in Miyagi Prefecture during the period from 1 March to 1 April 2011 in the horizontal and vertical directions; the black rectangle shows the common target area used in this study.

Figure 4 .
Figure 4. (a-b) Pre-and post-event aerial photographs taken by the Geospatial Information Authority of Japan (GSI) around the Natori GEONET station.(c) A colour composite of the pre-and post-event TerraSAR-X (TSX) intensity images taken in path C at the same location.(d) A handheld panoramic ground photo of the orange area in (b) that was taken during the field survey on 4 August 2011.

Figure 5 .
Figure 5. (a) Schematic image of the backscattering intensity of a building captured by a side-looking SAR sensor.(b) An example of displacement detection from a non-damaged building by the crosscorrelation method.

Figure 6 .
Figure 6.(a) Detected 2-D displacement vectors in each sub-area plotted on the pre-event TSX image from path A (the top, middle, and bottom vectors represent the results from paths A, B, and C respectively).(b-d) Estimated displacement vectors in each sub-area from the pairs of paths A and B, A and C, and B and C, plotted on the pre-event TSX image.

Figure 7 .
Figure 7. 3-D view and views of the relationship between the 3-D movement and SAR observation from the top and range directions.
all three combinations.The 3-D displacements in each subarea were estimated and are shown in Fig.6b-d.The red and blue vectors indicate movements in the horizontal and vertical directions respectively.Because the vertical displacements were much smaller than the horizontal displacements, they are shown at a different scale.

Figure 8 .
Figure 8.(a) Estimated 3-D displacement vectors obtained from the three pre-and post-event TSX pairs.(b) Image pairs used for the estimation.(c-d) Displacement amplitudes in the horizontal (c) and vertical (d) directions obtained by Kriging interpolation.

Table 1 .
Acquisition conditions of six pre-and post-event TerraSAR-X data used in this study.

Table 2 .
Detected 2-D displacements from three image pairs in the sub-areas that contain the GEONET ground stations (units: m).