Continent-wide 2-d Co-seismic Deformation of the 2015 Mw 8.3 Illapel, Chile Earthquake Derived from Sentinel-1a Data: Correction of Azimuth Co-registration Error

In this study, we mapped the co-seismic deformation of the 2015 Mw 8.3 Illapel, Chile earthquake with multiple Sentinel-1A TOPS data frames from both ascending and descending geometries. To meet the requirement of very high co-registration precision, an improved spectral diversity method was proposed to correct the co-registration slope error in the azimuth direction induced by multiple Sentinel-1A TOPS data frames. All phase jumps that appear in the conventional processing method have been corrected after applying the proposed method. The 2D deformation fields in the east-west and vertical directions are also resolved by combing D-InSAR and Offset Tracking measurements. The results reveal that the east-west component dominated the 2D displacement, where up to 2 m displacement towards the west was measured in the coastal area. Vertical deformations ranging between´0.25 and 0.25 m were found. The 2D displacements imply the collision of the Nazca plate squeezed the coast, which shows good accordance with the geological background of the region.


Introduction
On 16 September 2015, an earthquake of Mw 8.3 [1] occurred 48 km from the west of Illapel, Chile, causing extensive damages to the area and leading to a large tsunami [2].In this earthquake, at least 15 people died, 34 people injured, and tens of thousands of people became homeless.A total of 13 aftershocks ě Mw 6.0 soon afterwards occurred within 200 kilometers around the epicenter [1].Chile is in one of the most earthquake-prone regions in the world and has been attacked by more than a dozen of quakes ě Mw 8.0 (see Figure 1) since 1900.The Nazca plate is moving at an average velocity of 74 mm/year east-northeast, and subducting beneath the continent at the Peru-Chile Trench of the South American plate [3].The 2015 Mw 8.3 Illapel earthquake occurred as the result of thrust faulting on the interface between the Nazca and South American plates in central Chile [1,4].The USGS report [4] shows that the epicenter of this event locates at 31.57˝S, 71.654 ˝W with a focal depth of about 25 km, and its dimension is typically about 230 ˆ100 km (length ˆwidth).
Currently, few geodesy measurements (such as GPS, spirit leveling) are available for scientists to further study this event.Thanks to a new generation of SAR systems, such as the Sentinel-1A (S1A) [5,6] equipped with the Terrain Observation by Progressive Scans (TOPS) [7], the continent-wide Interferometric Synthetic Aperture Radar (InSAR) measurements becomes possible.Previous studies [8][9][10][11][12] have demonstrated the capability and great potential of S1A TOPS in ground surface deformation monitoring.In contrast to the conventional strip-map or ScanSAR data, multiple consecutive S1A TOPS frames (defined by ESA, where one zip file of S1A TOPS data distributed in Scihub [13] denotes a single frame) can generate long seamless strip data which is 250 km wide.Thus, it extremely benefits the study of large continent-wide event, such as the 2015 Mw 8.3 Illapel earthquake.
The S1A TOPS, as a new mode of InSAR, also has some key technical issues that need further study, especially for co-registering multiple S1A TOPS data frames.The Doppler Centroid of TOPS SAR data varies strongly along tracks, so a co-registration accuracy of at least 0.001 SLC pixel is required to eliminate the effect [7].Currently, the Spectral Diversity (SD) [14] method and the enhanced SD (ESD) [15], implemented in GAMMA [16] or GMTSAR [17], could perform well for a single frame data case, but cannot resolve the problem of co-registering multiple S1A TOPS data frames well.
In this study, we firstly proposed an improved Spectral Diversity (ISD) method to co-register the multiple S1A TOPS Interferometric Wide Swath (IWS) data frames.We then map the co-seismic LOS deformation associated with the 2015 Mw 8.3 Illapel earthquake using both ascending and descending S1A data and the methods of D-InSAR and Offset Tracking, respectively.Subsequently, the D-InSAR and Offset Tracking measurements are assessed in terms of precision and used together to retrieve the east-west and vertical co-seismic deformations of this earthquake.Finally, we discussed implications of ISD co-registration method and co-seismic measurements, as well as the potential contribution of the S1A data to further continent-wide event studies.[8][9][10][11][12] have demonstrated the capability and great potential of S1A TOPS in ground surface deformation monitoring.In contrast to the conventional strip-map or ScanSAR data, multiple consecutive S1A TOPS frames (defined by ESA, where one zip file of S1A TOPS data distributed in Scihub [13] denotes a single frame) can generate long seamless strip data which is 250 km wide.Thus, it extremely benefits the study of large continent-wide event, such as the 2015 Mw 8.3 Illapel earthquake.
The S1A TOPS, as a new mode of InSAR, also has some key technical issues that need further study, especially for co-registering multiple S1A TOPS data frames.The Doppler Centroid of TOPS SAR data varies strongly along tracks, so a co-registration accuracy of at least 0.001 SLC pixel is required to eliminate the effect [7].Currently, the Spectral Diversity (SD) [14] method and the enhanced SD (ESD) [15], implemented in GAMMA [16] or GMTSAR [17], could perform well for a single frame data case, but cannot resolve the problem of co-registering multiple S1A TOPS data frames well.
In this study, we firstly proposed an improved Spectral Diversity (ISD) method to co-register the multiple S1A TOPS Interferometric Wide Swath (IWS) data frames.We then map the co-seismic LOS deformation associated with the 2015 Mw 8.3 Illapel earthquake using both ascending and descending S1A data and the methods of D-InSAR and Offset Tracking, respectively.Subsequently, the D-InSAR and Offset Tracking measurements are assessed in terms of precision and used together to retrieve the east-west and vertical co-seismic deformations of this earthquake.Finally, we discussed implications of ISD co-registration method and co-seismic measurements, as well as the potential contribution of the S1A data to further continent-wide event studies.

ISD Method
Compared with the conventional strip-map InSAR, the TOPS mode InSAR has difficulty in azimuth co-registration due to a strong azimuth-variant Doppler centroid [7].The interferometric phase bias caused by the Doppler centroid induced azimuth mis-registration (or azimuth offset) can be expressed as [7,18]: where f DC is the Doppler centroid of the target point, d azi is the azimuth co-registration error, and PRF is the Pulse Repetition Frequency of focused SLC data.We should note that f DC varies linearly within a TOPS burst in the along-track dimension.Thus, a constant d azi will produce a series of azimuth phase ramps imposed onto an interferogram, resulting in phase discontinuity (or phase jump) on the burst interface [19].
To suppress the phase jump at burst overlap to ˘1.5 ˝, a co-registration accuracy of at least 0.001 SLC pixel is required [7,14].Currently, an overall mis-registration error is usually corrected with constant azimuth offset d azi (see Equation ( 1)) based on the method of SD [18] or ESD [14].For a single S1A IWS frame [13], which contains three sub-swaths and in each sub-swath there are typically nine burst, the number of burst is small (about 3 ˆ(9 ´1) = 24 of burst overlaps), and these conventional methods can work well.However, conventional methods with a constant correction would be ineffective in the case of multiple consecutive S1A TOPS frames [20], because large variations of the Doppler centroid will result in considerable variations of azimuth mis-registration (or azimuth offset) for long imaging duration.More importantly, obvious phase jump would also occur on overlap area between burst, especially at the ends of the long data taken.In order to address this problem, an improved Spectral Diversity method is proposed in this study.Rather than treating d azi as a constant, the new method compensates the variation of the azimuth mis-registration by a linear model: where d 0 azi is the intercept, k is the change rate of azimuth mis-registration error, and t is the acquisition time difference between the current and the first SLC line.By substituting Equation (2) into Equation (1), one gets: where s " PRF 2π¨f DC is a scale factor converting phase to SLC azimuth offset.Given N observations (N is the quantity of burst overlaps) we get a linear system in matrix form: where x " rd 0 azi , ks T is the unknown vector, and B is the design matrix: The observation vector L " rL 1 , L 2 , ¨¨¨, L N s T , where L i " φ i ¨si pi " 1, ¨¨¨, Nq is the azimuth offset of the i-th burst overlap, with φ i being the average of the phase differences of the current and the next burst interferogram of the i-th burst overlap and s i " PRF 2π¨f i DC being the scale factor of i-th burst overlap.
When N > 2, Equation ( 4) becomes over-determined problem.To resolve the parameters x, we use an iterative weighted least square estimator [21][22][23], to account for possible gross error in observation vector L. When the estimation of x was obtained, the azimuth co-registration error for each azimuth line could be calculated by Equation ( 2) and used to refine the SLC resample look-up table in S1A TOPS SLC co-registration, which will be discussed in next section.

Procedure of S1A TOPS SLC Co-Registration with ISD
To meet the requirement and to counteract the topographic effect, the co-registration of wide-swath SLC images of S1A TOPS SLC [16] generally involves a sequentially DEM-assisted method (Step 1), conventional Cross-Correlation (CC) (Step 2) method, and a SD method (Step 3).Our ISD method can be easily embedded into the S1A SLC co-registration procedure for multiple S1A TOPS frames.Figure 2 shows the flow chart of the SLC co-registration, and the three key steps are summarized as follows.

‚
Step 1: initial co-registration.Estimate initial SLC images offset and set up SLC resampling look-up table [24] with orbital state vectors and an external DEM, e.g., SRTM DEM, and resample the slave SLC into master imaging geometry.

‚
Step 2: co-registration by the CC method.Determine the image offsets between master and the resampled slave SLC from Step 1 by the CC method, and use the image offsets to refine the look-up table.The original slave SLC is, again, resampled into master imaging geometry.

‚
Step 3: co-registration by ISD method.Calculate the mean phase differences of all bursts overlaps between the master and the resampled slave SLC from Step 2. Then the constant azimuth offset d 0 azi and the change rate k are resolved by ISD method.The look-up table is refined with an azimuth offset calculated from Equation ( 2) and the estimated parameters and, finally, output the resampled slave SLC using the ultimately-refined look-up table.
Remote Sens. 2016, 8, 376 4 of 12 When N > 2, Equation ( 4) becomes over-determined problem.To resolve the parameters x , we use an iterative weighted least square estimator [21][22][23], to account for possible gross error in observation vector L .When the estimation of x was obtained, the azimuth co-registration error for each azimuth line could be calculated by Equation ( 2) and used to refine the SLC resample lookup table in S1A TOPS SLC co-registration, which will be discussed in next section.

Procedure of S1A TOPS SLC Co-Registration with ISD
To meet the requirement and to counteract the topographic effect, the co-registration of wideswath SLC images of S1A TOPS SLC [16] generally involves a sequentially DEM-assisted method (Step 1), conventional Cross-Correlation (CC) (Step 2) method, and a SD method (Step 3).Our ISD method can be easily embedded into the S1A SLC co-registration procedure for multiple S1A TOPS frames.Figure 2 shows the flow chart of the SLC co-registration, and the three key steps are summarized as follows.

•
Step 1: initial co-registration.Estimate initial SLC images offset and set up SLC resampling lookup table [24] with orbital state vectors and an external DEM, e.g., SRTM DEM, and resample the slave SLC into master imaging geometry.

•
Step 2: co-registration by the CC method.Determine the image offsets between master and the resampled slave SLC from Step 1 by the CC method, and use the image offsets to refine the lookup table.The original slave SLC is, again, resampled into master imaging geometry.
• Step 3: co-registration by ISD method.Calculate the mean phase differences of all bursts overlaps between the master and the resampled slave SLC from Step 2. Then the constant azimuth offset

Data Processing
Ascending and descending tracks S1A data are used in this study (see Figure 1b).The details of the SAR data are shown in Table 1.The methods of Differential InSAR (D-InSAR) and Offset tracking will be exploited to estimate the Light of Sight (LOS) deformations.Finally, 2-D co-seismic displacement fields of this event will be derived.

Data Processing
Ascending and descending tracks S1A data are used in this study (see Figure 1b).The details of the SAR data are shown in Table 1.The methods of Differential InSAR (D-InSAR) and Offset tracking will be exploited to estimate the Light of Sight (LOS) deformations.Finally, 2-D co-seismic displacement fields of this event will be derived.

Co-Registering Multiple Frames of S1A TOPS Data
The SAR co-registration follows the steps in Section 2.2.The 90 m SRTM DEM (v4.1) is used to estimate the initial look-up table [24] and remove the topographic phase.Thanks to the precise orbital state vector [6], an accuracy of a few pixels could be achieved [16] in the preliminary co-registration using orbital information.Then, the CC method is used to estimate the image offsets from evenly-distributed 64 ˆ192 (range ˆazimuth) tie-points with matching windows size of 128 ˆ256 (range ˆazimuth) pixels.All offsets are used to fit an offset polynomial that will subsequently be used to refine the look-up table.Due to limited accuracy of the CC method (~1/32 pixel) [25,26], the required co-registration accuracy of 0.001 SLC pixel could still not be achieved (see Table 2).Therefore, azimuth offset estimation with the proposed ISD is required for multiple S1A TOPS frames co-registration.We firstly calculated the observation vector L, the azimuth co-registration error, from each burst overlap.For descending data, three sub-swathes contain a number of (3 ˆ(27 ´1) = 84) burst overlaps.Since some burst overlaps situate in low coherent area (average coherence < 0.75), such as sea or vegetation areas, 69 valid burst overlaps are used to fit the azimuth co-registration error.For ascending data, sub-swath 1 is located in the sea area, which is discarded directly.Sub-swathes 2 and 3 provide a total of 34 burst overlaps for fitting the azimuth co-registration error.The final fitted linear models are, respectively, d des azi " 0.01320 ´2.1698 ˆ10 ´4 ¨t and d asc azi " 0.01477 `1.772 ˆ10 ´4 ¨t for descending and ascending tracks.The fitting results are also shown in Figure 3e.The azimuth co-registration accuracy with ISD is listed in Table 2, indicating that an accuracy of 0.001 SLC pixel has been achieved.Then the SLC resampling look-up table is refined by the fitted linear model.Finally, the slave SLC is precisely co-registered to the master SLC.To demonstrate the effectiveness of our ISD method, the datasets are also co-registered with SD method implemented in GAMMA [16].From comparisons in Figure 3a-d, we can see that obvious phase jumps still exist in the interferogram by the conventional SD method.c) mark the locations of phase jump; the solid green triangle, solid and dashed green line in (e) denote the samples of burst overlap azimuth offset, the fitted azimuth offset adopted by, the constant offset adopted by SD for the descending track, containing three S1A TOPS IWS frames with a number of 84 burst overlaps; the solid blue circles, solid and dashed blue lines in (e) represent the counterparts for the ascending track, containing two S1A TOPS IWS frames with 51 burst overlaps.

D-InSAR and Offset Tracking Processing
The TOPS D-InSAR processing, somewhat like traditional D-InSAR procedures [27], can be conducted after the SLC co-registration.As the Doppler centroid variations are different within bursts, we do not apply the azimuth common band filter in the interferogram [16].The topographic phase is removed with 90 m SRTM, thus generating the differential interferogram, and multi-looking factors of 20 × 4 (range × azimuth) is applied to set the output result with a resolution of about 50 × 50 m.The improved adaptive phase noise filter [28] is utilized to smooth the interferogram.Subsequently, the smoothed interferogram is unwrapped by network programming-based phase unwrapping method [29].Low coherent area (coherence ≤ 0.75) is masked out in the phase unwrapping process.Eventually, the LOS deformations from both the ascending and descending tracks are obtained (see Figure 4a,b).the solid green triangle, solid and dashed green line in (e) denote the samples of burst overlap azimuth offset, the fitted azimuth offset adopted by, the constant offset adopted by SD for the descending track, containing three S1A TOPS IWS frames with a number of 84 burst overlaps; the solid blue circles, solid and dashed blue lines in (e) represent the counterparts for the ascending track, containing two S1A TOPS IWS frames with 51 burst overlaps.

D-InSAR and Offset Tracking Processing
The TOPS D-InSAR processing, somewhat like traditional D-InSAR procedures [27], can be conducted after the SLC co-registration.As the Doppler centroid variations are different within bursts, we do not apply the azimuth common band filter in the interferogram [16].The topographic phase is removed with 90 m SRTM, thus generating the differential interferogram, and multi-looking factors of 20 ˆ4 (range ˆazimuth) is applied to set the output result with a resolution of about 50 ˆ50 m.The improved adaptive phase noise filter [28] is utilized to smooth the interferogram.Subsequently, the smoothed interferogram is unwrapped by network programming-based phase unwrapping method [29].Low coherent area (coherence ď 0.75) is masked out in the phase unwrapping process.Eventually, the LOS deformations from both the ascending and descending tracks are obtained (see Figure 4a,b).Apart from the D-InSAR measurements, the Offset Tracking [19] measurements in the slant range and azimuth direction are also achieved by conventional CC method [25] using the coregistered SLC images.The sliding step size of offset estimate window is set to be 20 × 4 (range × azimuth) pixels to generate the same resolution result with D-InSAR.To reduce the noise, center weighted median filter [30] with a window size of 11 × 11 (range × azimuth) is employed to smooth the offset measurement, in which the correlation coefficient of the matching windows is used as the weight factor.From our offset tracking measurement implementation, we find that the azimuth deformation is difficult to be detected.This is also confirmed by Multiple-Aperture Interferometry (MAI) technology, whose expected accuracy for S1A TOPS IWS data is about 27 cm [25].Thus, we do not use the azimuth offset measurements in our further study.The range offset measurements for ascending and descending orbits are shown in Figure 4d,e.

Precision Assessment and 2-D Displacement Retrieval
The precisions of both D-InSAR and Offset Tracking measurements are assessed by evaluating the local standard deviations (STDs) on sliding window with a size of about 1 × 1 km pixel-by-pixel [31].The linear trend is fitted and subtracted from the sliding window to get stationary signals for Apart from the D-InSAR measurements, the Offset Tracking [19] measurements in the slant range and azimuth direction are also achieved by conventional CC method [25] using the co-registered SLC images.The sliding step size of offset estimate window is set to be 20 ˆ4 (range ˆazimuth) pixels to generate the same resolution result with D-InSAR.To reduce the noise, center weighted median filter [30] with a window size of 11 ˆ11 (range ˆazimuth) is employed to smooth the offset measurement, in which the correlation coefficient of the matching windows is used as the weight factor.From our offset tracking measurement implementation, we find that the azimuth deformation is difficult to be detected.This is also confirmed by Multiple-Aperture Interferometry (MAI) technology, whose expected accuracy for S1A TOPS IWS data is about 27 cm [25].Thus, we do not use the azimuth offset measurements in our further study.The range offset measurements for ascending and descending orbits are shown in Figure 4d,e.

Precision Assessment and 2-D Displacement Retrieval
The precisions of both D-InSAR and Offset Tracking measurements are assessed by evaluating the local standard deviations (STDs) on sliding window with a size of about 1 ˆ1 km pixel-by-pixel [31].The linear trend is fitted and subtracted from the sliding window to get stationary signals for calculating the STDs.The local STDs maps of the measurements from D-InSAR and Offset Tracking in the slant range are shown in Figure S1, which will be used as the weight factor to derive the 2D displacement of this earthquake.
Theoretically, three dimensional (3D) deformations [32] could be resolved with LOS and along-track measurements from two or more viewing geometries [33,34].However, the north-south direction displacement is insensitive in D-InSAR measurements in general, and beyond the detectability of S1A TOPS IWS data in this earthquake in particular [25].Thus, the north-south deformation is compromisingly neglected in this study.The unit projection vector of LOS deformation for ascending and descending orbits are S asc " rs e , s n , s u s T " r´0.607, ´0.170, 0.755s T and S des " rs e , s n , s u s T " r0.608, ´0.168, 0.776s T , respectively.By assuming the north-south deformation is null and exploiting the projection vectors S asc and S des to establish observation equation, two dimensional (2D) deformations, i.e., east-west and up=down (or vertical) are derived by combing the D-InSAR and Offset Tracking measurements from ascending and descending orbits [35].We acknowledge that neglecting the north-south deformation will induce errors in the east-west and vertical deformation estimation, but for the Illapel earthquake it might be very slim, because the azimuth deformation itself is very small.The final 2D deformations are shown in Figure 4c,f.

Results
We got the co-seismic LOS deformations and the 2-D, i.e., east-west and vertical displacements with the elaborated data processing strategy presented in previous section.Figure 4a,d show the co-seismic LOS deformations extracted by D-InSAR and Range Offset Tracking with descending track InSAR data, and Figure 4b,e are the co-seismic LOS deformations with ascending track InSAR data.Due to the effect of decorrelation, there are still some areas along the Andes Mountains that has been masked out.
Comparing the D-InSAR and the Offset Tracking measurements, we can find that they have the same LOS deformation pattern and the same deformation magnitude, showing remarkably good agreement between each other.To make more detailed comparisons, the profiles along A-B and C-D are displayed as insets in Figure 4d,e.These results illustrate quite good consistency in the LOS D-InSAR and Offset Tracking measurements for S1A TOPS mode data.Comparing the D-InSAR measurements from both descending and ascending track, we can see that they have quite a similar deformation pattern, which implies that the horizontal displacement dominates the ground deformation; we can also see that most of the deformation with large magnitude locates on the upper part along profile A-B near the coast in descending track, while it is opposite in ascending track, implying that obvious vertical displacement has occurred and the upper part sank and the lower part uplifted.
Figure 4c shows the east-west displacement, where up to 2 m and towards the west occurred in the coastal area.Vertical displacement (Figure 4f) associated with the ground uplift and subsidence ranging between [´0.25, 0.25] m is also detected.It should be noted that east-west and vertical displacements are only calculated over the common area of all measurements (Figure 4a,b,d,e), which ensures that there are redundant observations for alleviating the effect of atmospheric [36] in the 2-D deformation inversion.

Discussion
To improve the accuracy of co-registering multiple S1A TOPS data frames, an ISD method is proposed and embedded into the S1A SLC co-registration procedure.Comparing with the conventional SD or ESD methods, the ISD method is more suitable for multiple S1A TOPS data frames co-registration by compensating the variation of azimuth co-registration error using a linear model.However, some issues shall be carefully settled to ensure the accuracy of ISD.Firstly, similar to the SD or ESD methods, the ISD method is subject to small azimuth offset measuring.The prerequisite for using ISD method is that the azimuth offset should be within ˘0.75 SLC pixels [19] Therefore, the azimuth offset shall firstly be truncated to this limit by other methods before applying the ISD method, e.g., the conventional CC method.Secondly, as shown in Equations ( 1) and (2), the observation, i.e., azimuth co-registration error d azi on each burst overlap, relates linearly to the phase offset φ.To obtain a high coherent phase offset, the burst overlap's phase difference of subsequent bursts shall be multi-looked with a large factor, e.g., 100 ˆ4 (range ˆazimuth) pixels and, subsequently, be filtered to reduce phase noise [28,31].The multi-looking factor depends on the phase quality, and we found that a factor of 100 ˆ4 is empirically enough.Additionally, it should be noted that the north-south deformation would cause an azimuth offset [15] too, which will superimpose on the normal azimuth co-registration error.If only a part of burst overlaps have north-south deformation, our ISD method would automatically exclude those observations by an iterative weighted least square estimator [21], as those observations possibly manifest as gross errors with respect to the normal azimuth offset observations.Furthermore, given the prior information we could also manually remove the observations contaminated by north-south deformation, and the final ISD accuracy would be improved in such a situation.
The 2015 Mw 8.3 Illapel earthquake occurred on the interface between the Nazca and South American plates (see Figure 1).The co-seismic deformation of this earthquake is mapped with S1A TOPS data in this study.Comparing our results with those obtained by the European Space Agency's Project on Sentinel-1 INSAR Performance Study with TOPS Data [37], we found that ours are continuous and clear, while their D-InSAR interferograms still have obvious phase jumps, similar to the results in Figure 3c, for both ascending and descending tracks.Affecting by the discontinuities of the D-InSAR measurements, their decomposed vertical displacement has considerable leap.All of these errors would provide false information for the interpretation of earthquake mechanism.On the contrary, there are no phase jumps in our results (see Figure 4a,b).In addition to the D-InSAR measurements, we also extract the Offset Tracking measurements in the LOS direction.Both D-InSAR and Offset Tracking LOS measurements show good consistency (see Figure 4a,b,d,e).This confirms that the SLC images have been precisely co-registered by our ISD method.In addition, we have derived the 2-D deformations from the LOS D-InSAR and Offset Tracking measurements.The results show that the coastal area moved westwards by up to 2 m, and the vertical displacement (see Figure 4f) associated with ground uplift and subsidence ranging between [´0.25, 0.25] m has also been detected.Additionally, by measuring the deformation area, we also found the dimension of the event is consistent with the result from USGS [4].Additionally, the 2-D displacements, a recent study [38] has tried to recover the 3-D displacements of the event, in which they use the along-track interferometry [38] to measure the azimuth deformation.By comparing them, we found that our results are in line with theirs.However, the along-track measuring [38] can only be applied to the burst overlap region, which occupies only about 15% of the whole frame.While, for the non-burst overlap region, occupying about 85% of the frame, interpolation processes must be implemented, which will, of course, induce errors.
Chile is prone to earthquakes that have caused severe destruction in the past century.The Nazca plate moves to the east-northeast at a varying rate (approximately 80-65 mm/year from north to south) relative to a fixed South American plate, resulting in complex geologic processes along the Nazca subduction zone [1].At the latitude of the 2015 Mw 8.3 Illapel earthquake, many massive earthquakes have occurred, including the 2010 Mw 8.8 Maule earthquake in central Chile [3,39].The 2D results show that the coastal area has moved westwards by up to 2 m accompanying with a little vertical deformation.This implies the collision of the Nazca plate squeezed the coast and the east-west displacement dominating the co-seismic deformation.The result is in good accordance with the geological background of the region.All those measurements would benefit for the geophysical parameters inversion of the 2015 Mw 8.3 Illapel earthquake.S1A satellite is now in normal on-orbit operation.As the first SAR satellite of Sentinel-1 constellation, S1A systematically acquires SAR image over land area with a revisit of 12 days in Europe and 24 days for other regions [40].S1A is controlled to follow an orbital tube of 100 m radius [41], which guarantees short spatial perpendicular baseline.These two merits reduce the decorrelation, which has been demonstrated in previous studies [27].As a C-band SAR system, S1A also has a limited detectable deformation gradient capability with D-InSAR technology [42].However, the increasing slant range resolution (up to 2 m) enhances the detectability along the slant range with the SAR Offset Tracking method, and makes it a good supplement for D-InSAR measurement.Additionally, the S1A acquires data with a wide swath (250 kilometer wide for IWS mode), which would be more favorable to study continent-wide event than the conventional SAR satellite (<100 km), like the ENVISAT and the ALOS.The advantages of S1A for large earthquake study have been presented in this work.The advent of S1A will boost geo-hazard monitoring and crisis management [43,44].

Conclusions
In this study, we mapped the co-seismic deformation of the 2015 Mw 8.3 Illapel earthquake using multiple S1A TOPS frames data.In order to achieve a co-registration accuracy of 0.001 pixel, we proposed an ISD method by compensating the variation of azimuth co-registration error in the azimuth direction with a linear model, and successfully embedded it into the S1A TOPS SLC co-registration procedure.The co-registration results of the Illapel earthquake demonstrate that our ISD method can be used for co-registering multiple frames of S1A TOPS data, which will be beneficial for continent-wide event study with S1A TOPS data.
The co-seismic LOS deformation of the Illapel earthquake is measured by both the D-InSAR and the Offset Tracking method.Comparing the D-InSAR and the Range Offset measurements, we found that they are consistent with each other, which reveals the enhanced range deformation detectability of S1A TOPS data.Thus, the range offset tracking can act as a good supplement for D-InSAR in the situation where interferometric phases are turned to be decorrelated because of large displacement.The 2D deformations (east-west and vertical directions) are also resolved based on the D-InSAR and Offset Tracking LOS measurements.Results show that the coastal area has moved westward by up to 2 m, and the vertical deformation varies between [´0.25, 0.25] m, implying the collision of the Nazca plate squeezed the coast.
Acknowledgments: This work is supported by the National Natural Science Foundation of China (41574005, 41474007 and 41404013) and the Fundamental Research Funds of Central South University (2014zzts049).The S1A data is from the ESA project and is downloaded from the Sentinel-1 Scientific Data Hub.Bathymetry data used in Figure 1 is provided by TOPEX, UCSD.Figures are produced by General Mapping Tools (v5.1.2).The GAMMA software (v20150702) is partially used in InSAR processing of S1A data.The author would like to thank the four anonymous reviewers for their very constructive remarks and suggestions.
Author Contributions: Bing Xu and Zhiwei Li conceive the research work and Bing Xu wrote the first draft of the paper.Guangcai Feng, Zeyu Zhang, Qijie Wang, Jun Hu, Xingguo Chen contributed to experiment implementation and result interpretation.All authors contributed to paper writing and revision.

Conflicts of Interest:
The authors declare no conflict of interest.

Figure 1 .
Figure 1.Tectonic settings of the 2015 Mw 8.3 Illapel earthquake and SAR data coverage.The main shock is presented by the yellow star.(a) shows earthquakes occurred along the interface between the Nazca and the South American plates (USGS, 2015a) since 1900.Red stars and dark purple dots represent earthquakes ≥Mw 8.0 and <Mw 8.0, respectively; (b) shows the topography of the area marked by black rectangle in (a) and the aftershocks occurred between 17 September 2015 and 1 October 2015.Red dots and dark purple dots represent earthquakes ≥Mw 6.0 and <Mw 6.0, respectively.The coverage of the S1A images is shown by blue rectangles.One rectangle denotes one S1A frame.The red line denotes the plate boundaries.

Figure 1 .
Figure 1.Tectonic settings of the 2015 Mw 8.3 Illapel earthquake and SAR data coverage.The main shock is presented by the yellow star.(a) shows earthquakes occurred along the interface between the Nazca and the South American plates (USGS, 2015a) since 1900.Red stars and dark purple dots represent earthquakes ěMw 8.0 and <Mw 8.0, respectively; (b) shows the topography of the area marked by black rectangle in (a) and the aftershocks occurred between 17 September 2015 and 1 October 2015.Red dots and dark purple dots represent earthquakes ěMw 6.0 and <Mw 6.0, respectively.The coverage of the S1A images is shown by blue rectangles.One rectangle denotes one S1A frame.The red line denotes the plate boundaries.
the change rate k are resolved by ISD method.The look-up table is refined with an azimuth offset calculated from Equation (2) and the estimated parameters and, finally, output the resampled slave SLC using the ultimately-refined look-up table.

Figure 2 .
Figure 2. Flowchart of S1A TOPS SLC co-registration.The procedure of S1A co-registration is somewhat like a 3× iteration, which gradually improves the accuracy of SLC co-registration.

Figure 2 .
Figure 2. Flowchart of S1A TOPS SLC co-registration.The procedure of S1A co-registration is somewhat like a 3ˆiteration, which gradually improves the accuracy of SLC co-registration.

Figure 3 .
Figure 3.Comparison of the co-registration results by SD and ISD method.(a) Differential interferogram in the descending track co-registered by the SD method; (b) differential interferogram in the descending track co-registered by the ISD method; (c) enlarged view of the rectangle area marked by c in (a); (d) enlarged view of the rectangle area marked by d in (b); (e) variation of the azimuth co-registration error and its linear fitting.Black arrows in (a) and (c) mark the locations of phase jump; the solid green triangle, solid and dashed green line in (e) denote the samples of burst overlap azimuth offset, the fitted azimuth offset adopted by, the constant offset adopted by SD for the descending track, containing three S1A TOPS IWS frames with a number of 84 burst overlaps; the solid blue circles, solid and dashed blue lines in (e) represent the counterparts for the ascending track, containing two S1A TOPS IWS frames with 51 burst overlaps.

Figure 3 .
Figure 3.Comparison of the co-registration results by SD and ISD method.(a) Differential interferogram in the descending track co-registered by the SD method; (b) differential interferogram in the descending track co-registered by the ISD method; (c) enlarged view of the rectangle area marked by c in (a); (d) enlarged view of the rectangle area marked by d in (b); (e) variation of the azimuth co-registration error and its linear fitting.Black arrows in (a) and (c) mark the locations of phase jump;the solid green triangle, solid and dashed green line in (e) denote the samples of burst overlap azimuth offset, the fitted azimuth offset adopted by, the constant offset adopted by SD for the descending track, containing three S1A TOPS IWS frames with a number of 84 burst overlaps; the solid blue circles, solid and dashed blue lines in (e) represent the counterparts for the ascending track, containing two S1A TOPS IWS frames with 51 burst overlaps.

Figure 4 .
Figure 4. Co-seismic deformation of the earthquake.(a) and (b) are the D-InSAR LOS deformations from ascending and descending geometries, respectively.A positive sign means the ground moves towards the sensor; (c) east-west deformation.(+: towards east); (d) and (e) are the LOS deformations determined by Offset Tracking; (f) vertical deformation (+: Up).The yellow star indicates the location of the main shock.The red line denotes the boundary of the Nazca and South American plates; the inner figure in (d) shows the deformation profiles by D-InSAR and Offset Tracking along section A-B for comparisons, and the inner figure in (e) shows that of the profiles along section C-D.

Figure 4 .
Figure 4. Co-seismic deformation of the earthquake.(a) and (b) are the D-InSAR LOS deformations from ascending and descending geometries, respectively.A positive sign means the ground moves towards the sensor; (c) east-west deformation.(+: towards east); (d) and (e) are the LOS deformations determined by Offset Tracking; (f) vertical deformation (+: Up).The yellow star indicates the location of the main shock.The red line denotes the boundary of the Nazca and South American plates; the inner figure in (d) shows the deformation profiles by D-InSAR and Offset Tracking along section A-B for comparisons, and the inner figure in (e) shows that of the profiles along section C-D.

Table 1 .
Information of SAR data.
1No. of Total Burst = No. of Frame ˆNo. of Sub-Swath per Frame ˆNo. of Burst per Sub-Swath.