Measurement of small co-seismic deformation field from multi-temporal SAR interferometry: application to the 19 September 2004 Huntoon Valley earthquake

ABSTRACT Interferometric synthetic aperture (InSAR) has been widely applied to natural disaster monitoring. However, it has limitations due to the influence of noise sources such as atmospheric and topographic artefacts, data processing errors, etc. In particular, atmospheric effect is one of the most prominent noise sources in InSAR for the monitoring of small magnitude deformations. In this paper, we proposed an efficient multi-temporal InSAR (MTInSAR) approach to measure small co-seismic deformations by minimizing atmospheric anomalies. This approach was applied to investigate the 18 September 2004 earthquake over Huntoon Valley, California, using 13 ascending and 22 descending ENVISAT synthetic aperture radar (SAR) images. The results showed that the co-seismic deformation was ±1.5 and ±1.0 cm in the horizontal and vertical directions, respectively. The earthquake source parameters were estimated using an elastic dislocation source from the ascending and descending acquisitions. The root mean square errors between the observed and modelled deformations were improved by the proposed MTInSAR approach to about 3.8 and 1.8 mm from about 4.0 and 5.2 mm in the ascending and descending orbits, respectively. It means that the MTInSAR approach presented herein remarkably improved the measurement performance of a small co-seismic deformation.


Introduction
Synthetic aperture radar (SAR) interferometry (InSAR) utilizes SAR images acquired at different times to measure surface deformation with an unprecedented spatial resolution (e.g. Massonnet & Feigl 1998;Lu 2007). The InSAR technique can be used to derive the source parameters associated with seismic and volcanic processes, which has allowed for better investigation of the mechanisms of geologic activities (e.g. Jung et al. 2011;Zhao et al. 2012;Deffontaines et al. 2016). Compared to seismology and other approaches, InSAR provides independent and remotely sensed information of earthquakes (Weston et al. 2011(Weston et al. , 2012. However, the InSAR technique based on a repeat-pass observation is faced with the limitation that the interferometric phase is delayed by atmospheric effects when the radar signals travel through the atmosphere (mainly troposphere and ionosphere). Especially, localized phase changes due to small deformation can suffer from severe atmospheric artefacts. For example, several studies CONTACT Hyung-Sup Jung hsjung@uos.ac.kr have indicated that atmospheric artefacts lead to errors of above 10 cm in the ground deformation measurements (Zebker et al. 1997;Lu & Dzurisin 2014). If the atmospheric anomalies in the interferometric measurement are not corrected, it would be impossible for small deformation to be detected well through the InSAR method.
In the atmospheric layers, ionospheric effect can be negligible when using short-wavelength radar signals because the effect is inversely proportional to the radar frequency (Jung et al. 2013;Liu et al. 2014;Jung & Lee 2015). In contrast, the tropospheric effect is independent of wavelength, and is mainly caused by refractive index changes due to tropospheric conditions, such as the atmospheric water vapour, pressure and temperature. Generally, short-wavelength SAR data is suitable for the measurement of a small deformation field because it is more sensitive than long-wavelength SAR data, once interferometric coherences can be maintained. Consequently, while the precise measurement of small deformations can be accomplished by the short-wavelength InSAR method, tropospheric artefacts remain problematic.
Several methods have been proposed to mitigate the atmospheric artefacts using external data such as global positioning system (GPS), meteorological data, MODIS, etc. (Li et al. 2005;Ding et al. 2008). However, it is not easy for the methods to mitigate the atmospheric artefacts due to the differences in spatial resolution and acquisition time of the external data. Therefore, a multi-temporal InSAR (MTInSAR) method can be another alternative method which could mitigate atmospheric artefacts. For example, the InSAR stacking method provides a good solution by averaging N independent interferograms and then reducing the errors to 1/ ffiffiffiffi N p (Biggs et al. 2007). However, if one pair among the stacked interferograms is contaminated by severe atmospheric artefacts, errors will remain in the resulting stacked interferogram. Moreover, this method causes degradation of the temporal resolution, and only works with linear ground deformations (Ding et al. 2008).
On the other hand, MTInSAR techniques such as Permanent Scatterer SAR Interferometry (Ferretti et al. 2000;Sara et al. 2015), Small Baseline Subset (SBAS) (Berardino et al. 2002;Jung et al. 2008;Lee et al. 2012), and Temporally Coherent Point InSAR (Zhang et al. , 2014 have been widely applied to measure ground deformations with high accuracy. These approaches achieve more efficient reduction of the atmospheric artefacts using low-pass filtering in the spatial domain and high-pass filtering in the temporal domain. A feasibility test for multi-temporal method is verified for atmosphere error, topography errors mitigation . Nevertheless, these methods face difficulties in the precise measurement of abrupt deformations such as co-seismic deformations, because they often assume linear surface deformation. Moreover, since they can misestimate small co-seismic deformations as being atmospheric artefacts, the MTInSAR methods are not readily applicable to the measurement of small co-seismic deformations. In this paper, we investigated an efficient approach to improve the measurement performance of small co-seismic deformations using an MTInSAR method. For this purpose, an efficient MTInSAR approach that is appropriate for the measurement of small co-seismic deformations associated with earthquakes was proposed, utilizing ENVISAT Advanced Synthetic Aperture Radar (ASAR) data. The approach can mitigate both atmospheric artefacts and topographic errors. Moreover, the approach is able to minimize the effect of abrupt deformations, such as small co-seismic deformations, by removing and adding back the abrupt deformations using elastic dislocation modelling (Okada 1985). Finally, we discuss how earthquake source parameters were improved using the MTInSAR approach and provide a comparison with various seismic catalogues.

Methodology
In this paper, we describe the main processing flow of the proposed method and briefly summarize the multi-temporal algorithm employed. We explain the way how we can improve the measurement of small surface deformation by reducing error sources from atmospheric artefacts, topographic error, etc. The interferometric phase (Df) was defined as follows (Berardino et al. 2002): where x and r are the azimuth and slant-range pixel coordinates, λ is the radar wavelength, r is the range from master sensor to pixel, i and j are the time variables for the master acquisition and slave acquisition, respectively. The first term accounts for the phase differences between times i and j from the surface displacement (Dd) in the line-of-sight (LOS) direction. The second term accounts for topographic error (Dz), which depends on the perpendicular baseline (B ? Þ and the SAR look angle (uÞ: The third term accounts for the atmospheric phase component (Df atmo ) between the acquisitions times at master i and slave j. If shorter wavelength SAR data are utilized, such as X-or C-band images, the atmospheric phase is mostly caused by tropospheric effect. The last term accounts for other phase noise sources (Df n ), such as the temporal and spatial decorrelation and thermal noise. It is well known that MTInSAR methods can effectively mitigate atmospheric artefacts and topographic errors (Berardino et al. 2002;Jung et al. 2008;Jo et al. 2010;Lee et al. 2010). The InSAR measurement from a co-seismic interferometric pair is not good enough to measure centimetre-level deformations which occur due to small earthquakes. Globally to date, a few researches have studied by small earthquake, such as 4.8 < Mw < 5.4 earthquakes in the Zagros Mountains, Iran (Lohman & Simons 2005), and the moment magnitude (Mw) 5.6 Little Skull Mountain earthquake (depth 9.4 km) (Lohman et al. 2002), and Mw 4.7 Katanning earthquake (depth < 1 km) and Mw 4.4 Kalannie earthquake (depth < 1.5 km) in Western Australia (Dawson et al. 2008). However, compared to earthquake occurrence frequency, the researches related to InSAR method are still less (Weston et al. 2011(Weston et al. , 2012. The main contributions to this deficiency come from the fact that the InSAR measurement is perturbed by atmospheric artefacts and topographic errors, as shown in Equation (1). Thus, it is necessary to reduce the atmospheric and topographic effects from the co-seismic deformations estimated by a single InSAR measurement. As mentioned earlier, the MTInSAR method is one of the best ways to reduce the atmospheric and topographic effects. However, abrupt deformations such as small co-seismic deformations can be misinterpreted by MTInSAR methods as being atmospheric artefacts.
In order to overcome these problems and allow for the precise measurement of small co-seismic deformations (DdÞ, we need to consider the following: (1) mitigate the topographic error (Dz) and atmospheric artifacts (Df atmo ); (2) minimize the effects of abrupt deformations; and (3) effectively remove the deformations before and after the earthquake from a co-seismic deformation pair. An efficient MTInSAR approach has been proposed to satisfy the considerations. Figure 1 summarizes the detailed workflow of the proposed MTInSAR approach, which is composed of five main steps: (1) creation of a model deformation map from a co-seismic deformation pair using elastic dislocation modelling (Okada 1985); (2) removal of the model deformation from co-seismic interferometric pairs; (3) generation of a digital elevation model (DEM) error map, atmospheric artefact maps and time-series deformation maps using a refined SBAS method (Jung et al. 2008;Lee et al. 2010); (4) addition of the model deformation map into deformation maps after the earthquake; and (5) creation of an improved co-seismic deformation map from the deformation maps right before and right after the earthquake.
The proposed MTInSAR approach estimates the model deformation from a co-seismic interferometric pair using the elastic dislocation model, after which the model deformation is removed from the co-seismic interferometric pairs. All interferograms used as the input data of the proposed approach do not contain any abrupt deformations. Thus, the problem of abrupt deformations being recognized as atmospheric components should be solved through this approach. The mitigation of topographic errors and atmospheric artefacts was performed by a refined SBAS method, which is characterized by time-series deformation refinement through an iteration process, reference phase correction, the finite-difference noise compression approach, etc. (Jung et al. 2008). Especially, using iterative approach to minimize interferograms errors, it leads to consistent deformation results. Since time-series deformations retrieved by the refined SBAS method do not include the model deformations, the model deformations need to be added back into the time-series deformations. The model deformations were added into the deformations after the earthquake, through which the final time-series deformations with an earthquake event can be acquired. Finally, an improved coseismic deformation map can be created by subtracting the deformation map right before an earthquake from the deformation map right after the earthquake. Consequently, the improved co-seismic deformations can be used for the elastic dislocation model again, to obtain improved earthquake source parameters and model deformations.

Study area
The study area includes the Huntoon Valley area, as shown in Figure 2. The Huntoon Valley is a fault-bounded basin within the Excelsior Mountains, CA, U.S.A. The land surface is covered by an arid semi-desert soil with little vegetation. The valley trends north-east between the Adobe Hills to the south-east, with the Excelsior Mountains to the north-west (Wesnousky 2005). There is little sign of recent fault activity other than local occurrences of over steepening at the base of the range front and a number of well-developed triangular facets on the north-west side of the valley (Wesnousky 2005). Generally, the reported slip rate was less than 0.2 mm/year (Adams & Sawyer 1998).
On 18 September 2004, a sequence of three Mw 5.2-5.6 events occurred just east of Mono Lake, beneath the Adobe Hills in the Huntoon Valley area ( Figure 2). To distinguish the earthquake, which caused surface deformation, the NCAeqDD data-set was compiled (Waldhauser & Scaff 2008). The NCAeqDD catalogue reported that the depths and magnitudes of the three earthquakes were 3.26, 7.15, and 8.76, with Mw 5.6, Mw 5.2, and Mw 5.4, respectively. The co-seismic deformations are highly dependent on the magnitude and depth while InSAR measurements are insensitive to magnitude 5.5 occurring deeper than 6 km, according to simulations (Dawson & Tregoning 2007). Thus, the co-seismic deformation with Mw 5.6 and 3.26 km reported by the NCAeqDD catalogue can be retrieved by the InSAR measurement. Consequently, in this study, focus was given only to measurement of the co-seismic deformation field for the first and largest earthquake, Mw 5.6, among the three events which occurred on 18 September 2004.

Data processing
To measure the co-seismic deformation field which occurred due to the Huntoon Valley earthquake, C-band ENVISAT ASAR raw data were collected from ascending and descending orbit acquisitions. Thirteen and 22 SAR raw images, which were obtained from the ascending and descending acquisitions from 2003 to 2010, were used in this study. In order to generate qualified interferograms, we first selected interferometric pairs under the conditions: (1) perpendicular baseline of less than 350 m; and (2) temporal baseline of less than five years, because the study area was a moderately mountainous desert. Thirty and 40 interferograms were generated in the ascending and descending acquisitions, respectively, and they were used for a dislocation inversion. Table 1 shows a list of the qualified interferograms selected for this study. For each interferometric pair, the azimuth common band filter was applied to the master and slave raw data to improve the geometric similarity between an interferometric pair and the interferometric coherence. The range-Doppler processing was applied to the master and slave raw data for the creation of the master and slave single look complex (SLC) images. For the sub-pixel-level co-registration processing, the slave image was resampled into the coordinate of the master image, and then an interferogram was generated by a complex conjugate operation of the master and resampled slave image, after which a range common band filtering was applied to the interferogram. A differential interferogram was created by a complex conjugate operation of the interferogram and a synthetic interferogram that was simulated from the 1-arc-second the national elevation dataset. After that, the differential interferogram was multilooked by 2 £ 10 looks (»40 £ 40 m) in the range and azimuth directions. The multilook operation was applied by using a two-step strategy, as follows: (1) 1 £ 5 looks before the flat-Earth and topographic phase corrections; and (2) 2 £ 2 looks after the corrections to reduce phase noise. The differential interferogram was smoothed using an adaptive filter (Goldstein & Werner, 1998) with a window size of 32 for reduction of the phase variance, and then unwrapped with a minimum cost unwrapping algorithm (Costantini 1998). Figure 3 shows the relation between the perpendicular and temporal baselines of the multiple interferograms used to retrieve the time-series deformations. Two different SBAS from both the ascending and descending acquisitions were used for the MTInSAR approach. For this approach, the interferogram pair with the smallest temporal baseline for the ascending and descending orbits was jointly used for estimation of the model deformations. The ascending model deformation was removed from the five co-seismic deformation pairs, and the descending deformation was removed from the seven co-seismic deformation pairs (see Table 1). The refined SBAS method was respectively applied to the ascending and descending interferometric pairs, and the time-series deformation maps for the ascending and descending orbits were produced. In addition, the DEM error map and atmospheric artefact

Co-seismic deformation measurement from an InSAR interferogram
Examples of the interferograms spanning the time windows before, during, and after the earthquake are shown in Figure 4. Signals with lobe patterns persisting in the co-seismic interferograms (Figure 4(b,e)) were unlikely to be atmospheric artefacts. However, the interferogram from the ascending track (Figure 4(b)) relatively showed a clear co-seismic deformation while the descending interferogram (Figure 4(e)) was quite noisy due to atmospheric influences. Atmospheric effects were also clearly visible in Figure 4(d,e). In addition, the topography-induced error could also be deduced from the interferograms in Figure 4 (see box A). The topographical error is dependent on the perpendicular baseline. Figure 4(b,f) may be not be contaminated by topographical errors because the baselines of these interferograms were short. The short baselines can make the interferograms insensitive to any plausible errors in the DEM. In contrast, the long perpendicular baseline interferograms, such as those in Figure 4(a,c), clearly displayed topographical errors. Such noise sources can decrease the accuracy of earthquake modelling for small deformation fields.
To remove the earthquake deformation, the co-seismic interferogram pair (20040829-20041003 and 20031029-20050720) with shortest temporal baseline from ascending and descending tracks were jointly modelled using an elastic dislocation source (Okada 1985), because the joint inversion of the ascending and descending acquisitions is closer to true geometry than the separate inversion (Lohman et al. 2002). The rectangular elastic dislocation source consists of 10 model parameters: two location coordinates for the centre of the source (x, y), depth (z), length, width, three components of slip (strike, dip, tensile), and strike and dip of the dislocation plane. The downhill simplex method and Monte Carlo simulations (Press et al. 1992) were applied to estimate the optimal parameters and their uncertainties, including the root mean square errors (RMSE) between the observed and modelled interferograms as the prediction-fit criterion. The best-fit parameters are listed in Table 2. The length of the Huntoon fault with a strike of about 2418 is 7.0 km, and the depth of the fault is about 6.5 km. From Okada modelling, the estimated depth represents the centre of fault. The fault dips about 83.28 to the north-east. The earthquake examined was dominated by a strike-slip component (about 21.0 cm) and accompanied with a slight dip-slip component (¡0.9 cm). For validation purposes, the co-seismic deformations from the ascending and descending orbits were compared with those simulated using the initial fault geometric parameters, listed in Table 2. Figure 5 shows the observed ( Figure 5(a,b)), modelled ( Figure 5(c,d)), and residual ( Figures 5(e,f)) deformations from the ascending and descending tracks, respectively. The RMSE values between the observed and modelled deformations were 4.0 and 5.2 mm for the ascending and descending acquisitions, respectively. The RMSE indicates the discrepancy between the observed and modelled deformations. It cannot be assumed that the RMSEs of 4.0 and 5.2 mm are small, since the co-seismic deformation in question was very small. The magnitude of the residual deformation was not much smaller than that of the modelled deformation for the descending case ( Figure 5(c,e)). In particular, as seen in Figure 5(d,f), the magnitude of the residual deformation was very similar to that of the modelled deformation. The dominant error sources included atmospheric effects and topographic error, as mentioned earlier. The measurement performance of the co-seismic deformation needs to be improved through minimization of such noise sources.

Measurement of co-seismic deformation from MTInSAR interferograms
To reduce the high RMS misfit of 4.0 and 5.2 mm, the measurement performance of the co-seismic deformation was improved using the MTInSAR approach described in Section 3.2. Figure 6 shows the topographical error maps derived from the ascending and descending MTInSAR approaches. It can be assumed that the topographic error estimation was carried out correctly because of the good match between the two topographic error maps. Although small differences between the two error maps were observed due to the small estimation errors of the atmospheric artefacts, time-series deformation signals, reference phases, such influences were negligible. As seen in box A of Figure 6, topographic errors leading up to about 20 m were retrieved. The topographic error of about 20 m  corresponded to the deformation of about 1.9 cm for the slant-range distance of about 840 km, the incidence angle of 238, and the perpendicular baseline of 300 m. Thus, topographic error correction was essentially required to measure the centimetre-level deformation. The deformation maps right before and right after the earthquake, from the ascending and descending acquisitions, are shown in Figure 7.  Figure 7(a) was not apparent, while that caused by the earthquake was evident, as shown in Figure 7(b). The same observations could be made for Figure 7(c,d). LOS deformation of about 2 cm was visible on the north-west side of the fault in the ascending track (Figure 7(b)), while LOS deformation of less than 2 cm was also detectable on the south-east side of the faults in the descending track (Figure 7(d)). It can be recognized that the LOS deformations from the ascending and descending acquisitions apparently presented a small deformation caused by earthquake.
To compare the time-series deformations before and after the earthquake, two points (AP1, DP1) were selected, which displayed the maximum deformation for both tracks (see Figure 7(b,d)). Figure 8(a) shows the time-series deformations of the ascending pairs at point AP1, while Figure 8(b) presents those of the descending pairs at point DP1. In the ascending case, the co-seismic deformation was obtained from observation of the long time span of 301 days from 8 September 2004 to 20 July 2005. The utilization of a long time span heightens the possibility that co-seismic deformation could be distorted by the deformations before and after the earthquake. Since the Huntoon earthquake event occurred on 18 September 2004, there is a high possibility that co-seismic deformation could be distorted by the deformation which occurred after the earthquake. However, no deformations were apparent after the earthquake, as seen in Figure 8(a). Therefore, it can be assumed that the distortion was negligible. The time span used in the descending case was 35 days. This was a much shorter span than the ascending case. Because no apparent deformations were observed after the earthquake (see Figure 8(b)), it was assumed that the distortion of the co-seismic deformation signal was negligible.
We next carried out elastic dislocation modelling of the improved co-seismic deformation fields using the proposed multi-temporal approach. The best-fit parameters are listed in Table 2. The length of the Huntoon fault with a strike of about 2368 was 8.0 km, and the depth was about 4.4 km. The fault dipped about 74.98 to the north-east. The earthquake was dominated by a strike-slip component (about 27.7 cm), and was accompanied with a slight dip-slip component (¡0.3 cm). The main difference between the source parameters obtained from the InSAR and MTInSAR measurements was that the fault depth of MTInSAR was much shallower than that of InSAR, as shown in Table 2. The differences are caused by the atmospheric and topographic errors that the InSAR measurements possess. Thus, it is desirable that the error sources of InSAR should be corrected to perform elastic dislocation modelling, especially for the precise measurement of small deformations. For validation purposes, the co-seismic deformations from the ascending and descending orbits were compared with those simulated using the final fault geometric parameters, listed in Table 2. Figure 9 shows the observed (Figure 9(a,b)), modelled (Figure 9(c,d)), and residual (Figure 9(e,f)) results from the ascending and descending tracks, respectively. Through the use of this modelling method, the RMSE between the observed and modelled deformations was reduced from 4.0 and 5.2 mm to 3.8 and 1.8 mm, based on the ascending and descending co-seismic interferograms after the multi-temporal approach, respectively. The RMSE from the descending track displayed higher reduction than that from the ascending track because more SAR images were used and the time span of the co-seismic observation was much shorter in the descending track. The use of more than 20 SAR images allowed mitigation of the atmospheric artefacts (Kim et al. 2007). The magnitude of the residual deformation was relatively smaller than that of the modelled deformation for the descending case. The dominant error sources, including atmospheric effect and topographic error, were successfully minimized by the proposed MTInSAR approach. Therefore, improvement of the measurement performance of the co-seismic deformation was achieved.

Two-dimensional deformation fields
InSAR images from the ascending and descending SAR tracks can be used to calculate the displacement components, which can allow retrieval of the east-west and vertical displacement components (Wright et al. 2004). Because both the ascending and descending tracks of ENVISAT are near-polar orbits, the north-south component cannot be resolved based on two LOS InSAR observations (Wright et al. 2004;Jung et al. 2011). Let d = (dx,dy) T be a two-dimensional (2D) deformation vector in a local (east, up) reference frame. If u is the unit LOS deformation vector in the same local reference frame, then it is defined by (sinucos', cosu) T , where u is the radar incidence angle from the vertical and ' is the satellite track angle from the north. Thus, the unit vector u (u asc , u dsc ) can be calculated based on u and ' from the ascending and descending tracks, respectively. The deformation r (r asc , r dsc ) represents LOS observations from both ascending and descending tracks. The deformation vector d = ¡(u T u) ¡1 (u T r) can be obtained, where u and r are given by u = (u asc , u dsc ) T and r = (r asc , r dsc ) T . The interferograms were precisely georeferenced to a common coordinate system before 2D decomposition.
The east-west and vertical displacement components based on the improved co-seismic deformation are shown in Figure 10(a,b), respectively. The deformation signals were very clear across the north-east-trending fault that produced the earthquake. Figure 10(c,d) shows the horizontal and vertical displacements along the profile A-A' across the fault (Figure 10(a,b)). The displacement in the east-west component was between ¡1.5 and 1 cm across the fault. Meanwhile, the deformation of the vertical component was between ¡0.5 and 1 cm. From the analysis of the deformation profiles, it can be concluded that the horizontal component dominated the deformation pattern, and that the deformation showed left-lateral strike-slip movement (Figure 10(e)), which is consistent with the modelling results.

Discussion
Herein, we compared the results of the developed InSAR-based focal mechanism with those of seismologic catalogues from different organizations using seismic data (Table 3). Especially, in the case of small deformation fields caused by earthquake events, few studies have included a comparison of InSAR and the earthquake catalogues due to atmospheric noise, data incoherence, or unfavourable earthquake depths (Weston et al. 2011). In the present study, three catalogues were compiled for  (Waldhauser & Schaff 2008). CMT is a global earthquake catalogue, while CISN and NCAeqDD are local earthquake catalogues. NCAeqDD has improved resolution only for hypocentre locations in central and Northern California, achieved through the use of waveform cross-correlation and double-difference methods. The focal mechanism solutions of the catalogues were slightly different due to possible errors in the velocity model, poor distribution of seismic stations, and the use of different algorithms for parameter determination (Mellors et al. 2004).
Using best-fit parameters for the dislocation source, the earthquake Mw was calculated from the equation: Mw = 2/3 log 10 (M 0 ) ¡ 10.7 (Hanks & Kanamori 1979), where M 0 is the seismic moment. The seismic moment is equal to mAS (Table 2), where m is the shear strength of the country rock, about 3 £ 10 11 dyne/cm 2 for typical continental crust, A is the area of the fault, and S is the average displacement along the fault plane. The Mw calculated from the InSAR-derived initial and final source parameters were both Mw 5.5. This estimation is the same as that from the CISN catalogue,  but slightly larger than that made by the CMT and slightly smaller than the NCAeqDD. According to the study by Weston et al. (2011), comparison of the seismic moment values calculated from CMT and InSAR source models revealed that the results predicted through InSAR are smaller than the results from seismic moments. Even though the magnitude of NCAeqDD is M L, however, the average difference between M L and Mw is negligible in Southern California area for M L > 4.0 events (Clinton et al. 2006). Moreover, the magnitude is slightly larger than that of InSAR, the difference of magnitude between the seismic catalogues, and the SAR data-set was quite small. Among the earthquake focal parameters, the depth also plays an important role in determining the surface displacement and seismic hazard. It is generally difficult to determine depth because the arrival time data are generally observed at stations on the surface of the Earth (Bai et al. 2006). In three seismological depth solutions, global centroid moment tensor (GCMT) depth means centroid depth which is same meaning with InSAR-derived depth, while NCAeqDD and CISN depth represents hypocentre. Therefore, we cannot compare the depth information directly between InSARderived depth and hypocentre except for GCMT. However, in case of moderate earthquake, when ground deformation is identified in InSAR, hypocentre location can be reliably estimated within 1-2 km. It means that seismological depth solution is consistent with the depths determined by modelling our procedure except for GCMT depth. The depth of 4.38 km calculated using the MTInSAR approach agrees well with that of 3.2 and 7.6 km calculated from NCAeqDD and CISN, respectively, within the error of the fault dimension, while 12 km determined by CMT ( Figure 11). Especially, the CMT catalogues have very little sensitivity with depth; therefore, typically the depth is fixed (Weston et al. 2011). In this sense, the earthquake depth from CMT is difficult to resolve for shallow event.
Since remarkable improvement of the RMSE was achieved through the MTInSAR approach developed herein, the focal depth calculated from InSAR would be about 4.4 km rather than 6.5 km, which was the modelling result from one SAR interferogram. In addition, because the NCAeqDD catalogue provides a theoretically more accurate location, the focal depth from the earthquake catalogues of about 3.2 km would be more accurate, rather than the 12.0 km calculated from CMT or 7.6 km from CISN. The results showed the following: (1) that the CMT and CISN catalogues may not be accurate for source parameter estimation of small earthquakes, (2) that the NCAeqDD catalogue was the best choice among the three catalogues, and (3) that the MTInSAR approach developed herein was much better than one interferogram approach for precise estimation of the source parameters.
It should be noted that the following points need to be discussed about the proposed approach and its validation. (1) The performance test of the proposed method according to the initial earthquake correction might be required for more accurate validation. It can be done by using the Figure 11. Comparison of the catalogue-and InSAR-derived focal depths, where CMT is the Global Centroid Moment Tensor, CISN is the California Integrated Seismic Network, and NCAeqDD is the Double-Difference Earthquake Catalog for Northern California.
synthetic interferograms. (2) Although the initial co-seismic deformation modelled by the dislocation inversion was used for this study, it can be replaced into one of other simulated deformations. For example, Mogi or Yang models can be used for volcanic deformation modelling as well as the in situ data such as GPS network can be used to generate the initial deformation. (3) An iteration approach of the proposed method, which replaces the initial elastic model into the final elastic model, might further improve the measurement performance of small deformation fields, even though there is no improvement in this study. (4) If an interferometric pair with shorter temporal baseline such as Sentinel-1 interferometric wide swath mode can be available, the proposed method would have a better performance (Barra et al. 2016).

Conclusions
In this paper, a MTInSAR approach was applied for the measurement of small co-seismic deformation fields with reduction of error sources such as atmospheric error, topographic artefacts, etc. The MTInSAR approach was composed of five main steps: (1) creation of a model deformation map from a co-seismic deformation pair using elastic dislocation modelling; (2) removal of the model deformation from the co-seismic interferometric pairs; (3) generation of a DEM error map, atmospheric artefact maps and time-series deformation maps using a refined SBAS method; (4) addition of the model deformation map to the deformation maps after earthquake; and (5) creation of an improved co-seismic deformation map from the deformation maps right before and right after the earthquake. The developed approach was applied for investigation of the 18 September 2004 earthquake over Huntoon Valley, CA, U.S.A. using 13 ascending and 22 descending ENVISAT SAR images.
The results showed that the displacements of the earthquake were §1.5 and §1.0 cm in the horizontal and vertical components, respectively. Moreover, the earthquake source parameters were estimated using a dislocation source by joint modelling of the improved deformation fields from the ascending and descending geometries. The best-fit fault model struck approximately N-E with a length of about 8 km and a width of 3 km, centred at a depth of 4.4 km. The RMSEs of the final modelling based on the co-seismic deformation after multi-temporal processing were 3.8 and 1.8 mm, while those of the initial modelling based on one pair of the best-quality co-seismic interferograms were 4.0 and 5.2 mm from the ascending and descending orbits, respectively. This demonstrates that the MTInSAR approach developed herein can reduce the atmospheric artefacts and topographic error to allow enhanced retrieval of earthquake source parameters.
The source parameters derived from InSAR were compared with those obtained from three earthquake catalogues: CMT, CISN, and NCAeqDD. In three seismological depth solutions, only GCMT depth means centroid depth while, NCAeqDD and CISN depth represents hypocentre. However, in case of moderate earthquake, we could compare depth solution roughly with modelling our procedure. The focal depth of 4.38 km calculated by the MTInSAR approach showed good agreement with that of 3.2 km derived by NCAeqDD, rather than those of 7.6 and 12 km from CISN and CMT. The depth of 12 km estimated by the CMT catalogue was too deep to be believable, and varied largely from those obtained through other catalogues and InSAR-derived models. The comparison revealed the following: (1) that the CMT and CISN catalogues may not be accurate for source parameter estimation of small earthquakes; (2) that the NCAeqDD catalogue was the best choice among the three catalogues; and (3) that the MTInSAR approach proposed herein was much better than the one interferogram approach for the precise estimation of source parameters.

Disclosure statement
No potential conflict of interest was reported by the authors.