Complete Three‐Dimensional Coseismic Deformation Field of the 2016 Central Tottori Earthquake by Integrating Left‐ and Right‐Looking InSAR Observations With the Improved SM‐VCE Method

The three‐dimensional (3‐D) deformation field associated with the 2016 Central Tottori earthquake is retrieved from advanced land observing satellite 2 interferometric synthetic aperture radar (InSAR) observations with four different viewing geometries, that is, ascending/descending tracks and left‐/right‐looking modes. The strain model and variance component estimation (SM, VCE, SM‐VCE) method is exploited and improved to integrate the InSAR observations with different viewing geometries so that the 3‐D deformation field is not affected by the inconsistent coverage of SAR footprints or the gross errors in InSAR observations. The obtained results are consistent with GNSS observations, indicating that the improved SM‐VCE method, known as SM‐RVCE in this paper, is capable of retrieving an accurate and spatially complete 3‐D deformation field for this earthquake. In addition, the precision of the InSAR observations and the estimated 3‐D deformation are quantitatively assessed by the SM‐RVCE method. Finally, on the basis of the estimated 3‐D coseismic deformation, the source parameters of this event are inverted, revealing an asperity with a maximum strike‐parallel slip of ~1.1 m concentrated at depths between 2 and 10 km. The estimated seismic moment is 2.4 × 1018 Nm, which corresponds to a Mw 6.2 event.


Introduction
At 05:07:22 (UTC) on 21 October 2016, an earthquake of magnitude 6.6 determined by the Japan Meteorological Agency (M JMA 6.6) struck Central Tottori Prefecture (Zhao et al., 2018), Western Honshu, Japan. The epicenter was located~6 km south of Kurayoshi, a city of Tottori Prefecture, at a depth of~5.6 km (USGS, 2016). Many infrastructures and residential buildings in the area with a radius of 20 km centered on the epicenter suffered serious damage. Fortunately, no people died in the event (REUTERS, 2016). In addition, this earthquake was slightly felt in the surrounding prefectures, such as Okayama and Shimane (USGS, 2016).
Highly accurate and spatially continuous three-dimensional (3-D) surface deformation can provide invaluable data for the interpretation and assessment of geo-hazards, such as earthquakes and volcanic eruptions (Morishita et al., 2016;Rocca, 2003;Wang et al., 2014). Interferometric synthetic aperture radar (InSAR), characterized by its high spatial density of measurements, large spatial coverage, and high measurement accuracy, has been widely used in the monitoring of ground surface deformation (Bawden et al., 2001;Fialko et al., 2005;Funning, Barke, et al., 2005;Gabriel et al., 1989;Goldstein et al., 1993;Hilley et al., 2004;Massonet et al., 1995;Massonnet et al., 1993). However, this technique can capture only onedimensional deformation along the line-of-sight (LOS) direction. Integrating multiple InSAR measurements with greatly varied imaging geometries can generate 3-D surface deformation in high-latitude regions (Gray, 2011), but this approach is not successful in mid-or low-latitude regions because present-day InSAR observations are basically obtained with the right-looking imaging mode of SAR satellites. Furthermore, the measurements from the pixel offset tracking (POT) (Michel et al., 1999) and multi-aperture InSAR (MAI) (Bechor & Zebker, 2006;Jung et al., 2009) techniques, typified by the complementary azimuth direction but an ungraded accuracy, can be integrated with InSAR LOS measurements to obtain significant 3-D surface deformation generated by geophysical events, such as earthquakes (Fialko et al., 2001;Xu et al., 2018), volcanic eruptions (Jo et al., 2017;Xu et al., 2017), and glacier movements Milillo et al., 2017;Zhou et al., 2019), on the basis of both ascending and descending SAR images.
Although these methods can generate surface deformation in three dimensions, they are subject to either the ill-posed configuration of multiple InSAR right-looking observations or the poor accuracy of ungraded observations (i.e., POT/MAI-derived azimuth deformation). If the left-looking mode, which has imaging geometries notably different from the right-looking mode, of SAR satellites can be exploited, highly accurate 3-D surface deformation maps can be extracted from the integration of ascending and descending InSAR LOS measurements (Wright et al., 2004). For instance, Morishita et al. (2016) combined ascending/descending and right-/left-looking advanced land observing satellite 2 (ALOS-2) InSAR measurements to retrieve the 3-D deformation generated by the volcanic activity of Sakurajima on 15 August 2015.
The ALOS-2 satellite imaged the 2016 Central Tottori earthquake along both ascending and descending tracks with left-and right-looking modes, providing invaluable data for retrieving the 3-D ground deformation. However, the footprints of several ALOS-2 images do not completely cover the seismic event, leading to difficulties in extracting the associated 3-D deformation. In this paper, we intend to estimate the spatially complete 3-D deformation field associated with the 2016 Central Tottori earthquake by integrating available left-/right-looking InSAR measurements acquired on ascending and descending tracks. Through this integration, we exploit and improve a recently proposed method (Liu, Hu, Li, Zhu, et al., 2018), known as strain model and variance component estimation (SM, VCE, and SM-VCE), to eliminate the boundary effects and gross errors caused by the inconsistent coverage of the different InSAR measurements. Moreover, the introduction of the VCE algorithm provides an opportunity to quantitatively assess the precision of different InSAR observations. To the best of our knowledge, this is the first attempt to assess the precision of InSAR observations in a posterior way without using any external data (e.g., leveling and GNSS data). Some preliminary results derived by the SM-VCE method can be found in . The 3-D deformation field derived by the SM-RVCE method proposed in this paper is superior to that presented in . In addition, a fault slip inversion with respect to this earthquake is conducted, and a detailed discussion is provided.

Data
The event that occurred in Central Tottori, Japan, on 21 October 2016 was observed by four interferometric pairs acquired by the ALOS-2 satellite from ascending and descending tracks with both left-and rightlooking modes (Table 1). The footprints of these four interferograms are shown in Figure 1a. The study area (i.e., the blue rectangle in Figure 1a) is split into three zones by these footprints, labeled as A, B and C. Note that different zones are imaged by different interferogram assemblies; for example, Zone A is imaged by ascending right-looking interferograms and by descending right-and left-looking interferograms.
Gamma software is employed to generate the LOS deformation measurements on the basis of the two-pass differential InSAR technique (Gabriel et al., 1989). In the interferometry process, a multilooking operation with 24 looks in the range direction and 30 looks in the azimuth direction is applied to suppress decorrelation noises, resulting in a final spatial resolution of approximately 60 m × 60 m. The Shuttle Radar Topographic Mission 1-arcsecond (~30 m spacing) digital elevation model is utilized to remove the topographic contribution from the interferograms. Then, decorrelation noise is minimized by an improved Goldstein filter . The minimum cost flow method (Chen & Zebker, 2002) is exploited to unwrap the interferometric phase by masking the areas with coherence smaller than 0.4. Then, reference points are selected from the extremely far field, which is relatively stable. In addition, a linear polynomial model is employed to remove the possible phase ramp associated with orbital and ionospheric errors (Hu, Wang, Li, et al., 2014). Furthermore, as the study area is located in a mountainous region and is situated along the seaside, data from the Generic Atmospheric Correction Online Service for InSAR (GACOS) are applied to correct the possible tropospheric delay (Yu et al., 2018).
The final geocoded LOS deformation maps with different viewing angles are shown in Figure 2. The two descending measurements image the whole area of interest, whereas the coverage of the two ascending measurements is incomplete but complementary. The deformation patterns of the interferograms extracted from the ascending left-looking (AsL) and descending right-looking (DesR) measurements are very similar. However, the ascending right-looking (AsR) and descending left-looking (DesL) images show distinct differences to the immediate southwest of the epicenter, which can be ascribed to the north deformation component. On the other hand, the DesL interferogram, though with a time span of approximately 2 years, is highly coherent in the study area due to the robust capability of L-band SAR to resist temporal decorrelation (Strozzi et al., 2008). In contrast, the AsL interferogram, whose time interval is shorter than that of the DesL interferogram, shows obvious decorrelation errors in the lower right part of the image, which might be caused by snow cover in the master image acquired in January as well as by residuals produced by the topography-related atmospheric delay.
The Japan mainland has been experiencing long-term tectonic deformation according to the GNSS network (http://geodesy.unr.edu/). However, the interferograms in this paper can be considered to be free of interseismic deformation since such deformation is basically mitigated by the unwrapping operation due to the very similar patterns throughout the study area. In addition, by investigating the deformation time series from the GNSS stations in the study area, little deformation is observed following the second day after the earthquake, that is, 22 October 2016. Compared with the centimeter-level coseismic deformation, this slight postseismic deformation can generally be neglected (Jung et al., 2011;Wang et al., 2015). Therefore, the errors induced by the discrepancies among the slave images' dates of the four interferograms can also be ignored in the estimation of the 3-D deformation. A detailed analysis of the GNSS observations can be found in the supporting information (i.e., Text S1 and Figures S1 and S2).

Methodology
Since the interferograms from the ascending track (i.e., AsL and AsR) cover only part of the study area, the spatially complete 3-D surface deformation field cannot be recovered by classic pixel-by-pixel methods, such as the weighted least squares (WLS) algorithm (Jung et al., 2011).
The SM-VCE method has the potential to retrieve complete and accurate 3-D surface deformation field, and it has two advantages over the WLS algorithm. First, the SM-VCE method exploits the spatial correlation of adjacent points' deformation based on the geo-kinematic model (i.e., strain model and SM) to establish the observation functions; this approach is obviously superior to the purely mathematical filtering procedure and the pixel-by-pixel-based WLS method. Second, accurate weights can be determined by the VCE algorithm with the SM-based established observation functions. However, based on the pixel-by-pixel solutions, the VCE algorithm is inapplicable since only a few redundant observations are available. Therefore, the weights can be determined only on the hypothesis that the observations are ergodic (Jung et al., 2011), which is hardly satisfied in reality.
Some gross errors induced by phase unwrapping and extreme tropospheric delays can seriously distort InSAR measurements, but they are hardly distinguishable from deformation signals. In this study, the SM-VCE method (Liu, Hu, Li, Zhu, et al., 2018) is improved by introducing the robust VCE algorithm (Yang et al., 1999) to eliminate gross errors; the improved method is known as SM-RVCE in this paper and is briefly described as follows. Detailed definitions of some variables can be found in Liu, Hu, Li, Zhu, et al. (2018).

Establishing Observation Functions Based on the Strain Model
The strain model reflects the relationship between the strain field and the deformation field on the basis of elastic theory (Vaníček et al., 2008). According to the concepts of the SM, a portion of the Earth's surface that is deformed during a geodynamic process (e.g., an earthquake) can be regarded as a homogeneous strain field Guglielmino, Bonforte, et al., 2011). We suppose that a point of interest P 0 with the ) in a window with a size of 1 × 1 km, whose 3-D position and deformation are The superscripts e, n, and u denote the east, north, and up components, respectively.
Ã T denotes the coordinate-increment vector from point P 0 to point P i , the relationship between d i and d 0 can be expressed as: where H=S + R is the displacement gradient matrix and S and R, the symmetric and antisymmetric parts, respectively, can be written as: where ξ and ω represent the parameters of the strain tensor and rigid body rotation tensor, respectively. According to equation (1), the observation function between the unknown vector l and d i can be modeled by: ð Þ ξ ee ξ en ξ eu ξ nn ξ nu ξ uu ω en ω eu ω nu Â Ã T and B i sm is the transform matrix between d i and l (Liu, Hu, Li, Zhu, et al., 2018).

Assuming that there is an observation vector
ments of point P i with the subscripts 1,2,3,and 4 denoting the AsL, AsR, DesL, and DesR measurements in this paper, respectively. The relationship between the 3-D deformation components d i and the observation vector L i can be expressed by: where B i geo is the transform matrix written as: Þ is a unit vector in the LOS direction from the ground to the satellite for the jth measurement associated with point P i . where Here, α i and θ i are the satellite heading angle (clockwise from north) and the radar incidence angle of point P i , respectively.
Combining equations (4) and (5), we obtain the following relationship between the unknown vector l and the InSAR measurements L i : and the numbers 1-12 in the top row indicate the column indexes of this matrix.
For K points surrounding P 0 , the overall observation system can be constructed on the basis of equation (6): where V with a size of 4K × 1 is the residual vector. The observation vector L with the same size as V and the coefficient matrix B with a size of 4K × 12 can be respectively expressed as:

Determining Accurate Weighting Factors From the Robust VCE Algorithm
The InSAR measurements used in this study can be divided into four groups corresponding to AsL, AsR, DesL, and DesR measurements. For simplicity, we assume that the initial weight matrix of each group of observations is equal to the unit matrix, that is, W 1 = W 2 = W 3 = W 4 = I, where I is the unit matrix with a size of K × K. L j and B j (j = 1,2,3,4) correspond to the observation and design matrixes, respectively, of the jth group, which can be extracted from L and B (i.e., equations (8) and (9), respectively).
Subsequently, the relationship between the observation residual δ and the variance component b σ 2 can be constructed as (Hu & Xiao, 2016 and Γ is the transform matrix between b σ 2 and δ (Liu, Hu, Li, Zhu, et al., 2018).
After determining the variance component b σ 2 , we can obtain the standard deviation (SD) of each measurement with respect to P 0 . Note that this is not considered in the SM-VCE method.
The variance component b σ 2 is also exploited to update the weight matrixes of the four groups of InSAR measurements: Since the performance of the VCE algorithm greatly degrades in the presence of gross errors (such as those in the AsL measurement), we propose to reupdate the determined weight factors in equation (13) on the basis of the robust VCE algorithm (Yang et al., 1999). Assuming that c W i j is the ith diagonal element of c W j , which corresponds to the weight factor of L i j , the reupdated weight factor c c W i j can be retrieved from the following equivalent weight function (Yang et al., 1999): where follows the standard normal distribution, and only~0.3% of the random error is larger than 3×SD 0 j . Therefore, the values of u i j >k 1 ¼ 3 will be taken as mistakes. k 0 = 1.5 is an empirical bound to downweight observations with relatively high corrections v i j , when the observation errors follow the normal distribution (Yang et al., 1999).
The reupdated weight matrixes c c W j in equation ((14)) are then used in equations ((10)), ((11)), ((13)), and ( (14)) to iteratively update the weights until the following equation is satisfied (Hu & Xiao, 2016 where b σ 2 max and b σ 2 min are the maximum and minimum, respectively, of the variance component b σ 2 , and ε 2 is the convergence criterion. Since b σ 2 1 is finally derived and determined to be the SD of L 1 , ε 2 can be determined based on the general accuracy of the observations. In this study, we use 1.0 × 10 −6 since the accuracy of InSAR deformation measurements is generally at the millimeter level. In addition, the accuracy of the 3-D deformation field can be calculated from the SD of each measurement according to the law of error propagation. Since independent observations should be guaranteed when using the law of error propagation, the relationship among the observations of adjacent points in equation (6) is not exploited to estimate the accuracy of the 3-D deformation. Hence, we simplify the observation function (equation (6)) as: The variance-covariance matrix of the 3-D deformation components D 0 can thus be formed as:

Fault Slip Inversion
The source geometry and fault slip distribution are determined from the generated 3-D coseismic deformation field. We subsample the data (i.e., Figure S3) using the quadtree method (Jónsson et al., 2002), which divides the data into squares until the values within each box do not exceed a certain variance threshold. The average value of the corresponding pixels is assigned to their local point. Guided by the focal mechanism solution and the distribution of aftershocks, we use a single rectangular dislocation model to model the data while assuming a homogeneous and isotropic elastic half-space with a Poisson's ratio of 0.25 and a shear modulus of 30 GPa. We set relatively tight bounds for the strike and location of the causative fault. Considering that InSAR is least sensitive to the north-south displacement, we set the lowest weight to this component and give equal weights to the other two displacement components in the inversion. First, the optimal fault model parameters with uniform slip are estimated simultaneously using a Monte Carlo-type simulated annealing algorithm, followed by a gradient-based iterative method using the output from the simulated annealing as the initial estimate (Cervelli et al., 2001). To retrieve the fault slip distribution model, we enlarge the dimensions of the fault to 30 × 30 km and discretize the fault into small patches (2 × 2 km). A Laplacian smoothing condition (i.e., Figure S4) is applied to avoid abrupt variations in the slip estimation. Under the fault geometry estimated from the first step, the slips on the patches are solved directly using a nonnegative least squares method.

Results
The spatially complete 3-D deformation field (see Figures 3a-3c) associated with the 2016 Central Tottori earthquake is derived from only InSAR LOS measurements with four distinct viewing angles (i.e., AsL, AsR, DesL, and DesR) using the SM-RVCE method. For the SM-RVCE method, the number of adjacent points in the limited range is vital for the accuracy and efficiency of the algorithm. A series of experiments with window sizes ranging from 5 × 5 to 65 × 65 pixels are carried out to assess the performance of the SM-RVCE method. The root mean square errors (RMSEs) of all three components first increase and then decrease as the window size increases (see Figure 4). This is expected since a larger window could involve more inhomogeneous pixels to estimate the 3-D deformation, which is not favored by the SM-RVCE method.
On the other hand, the computational burden on the SM-RVCE method generally increases with increasing window size. However, when the window size is small (i.e., 5 × 5 pixels), both the computational burden and the RMSE are quite large. This can be attributed to the small number of points involved in the computation, resulting in more iterations of the RVCE method and a lower accuracy of the estimated 3-D deformation. Therefore, to obtain a tradeoff between the 3-D deformation accuracy and the computational burden, a fixed window of 21 × 21 pixels (i.e., approximately 1 × 1 km) is preferred in this paper. As shown in Figures 3a-3c, the maximum deformation values on the east, north, and up components are approximately 8.7, 11.7, and 4.6 cm, respectively. The vertical coseismic deformation zone (i.e., Figure 3c) is divided into four sections. In the southeast and northwest sections, a subsidence pattern with generally epicenter-ward horizontal deformation is observed, while the opposite deformation pattern appears in the other two sections.
The 3-D deformation predicted by the model (Figures 3d-3f) generally agrees with the InSAR-derived deformation (i.e., Figures 3a-3c), and their differences are shown in Figures 3g-3i. Figure S5 exhibits the histograms of these differences, which normally obey a Gaussian distribution and have RMSEs of 1.0, 2.7, and 0.7 cm for the east, north, and vertical components, respectively. The residual on the north component is much higher than those on the east and up components. This is expected since InSAR measurements are insensitive to the deformation on the north component. In addition, we calculate the residuals for the original InSAR observations after subtracting those calculated from the estimated 3-D deformation. As shown in Figure S6, the residuals for the AsL, AsR, DesL, and DesR observations generally obey a Gaussian distribution with RMSEs of 1.8, 0.4, 0.4, and 0.4 cm, respectively. A relatively large value can be observed in the southern part of the AsL residual map due to the contamination of decorrelation noise. Two profiles across the fault trace are selected for detailed analysis. As shown in Figure 3j, some discrepancies are observed on the east component along profile D-D' and on the north component along profile E-E'. The tropospheric delay caused by the steep topographic slope and the presence of serious water vapor along the seaside can hardly be corrected by the GACOS data. The residuals due to the tropospheric delay could be responsible for these discrepancies. Since GACOS employs 0.125°(~14 km) gridded ECMWF data and continuous GNSS tropospheric delay estimates to produce global and near-real-time tropospheric delay products (Yu et al., 2018), these data can correct only the long-wavelength tropospheric delays whose length scales are longer than 14 km. Another possible reason for these discrepancies is the strong constraint in the fault slip inversion. As such, the data might be better fitted if we relax the fault geometry and the rake angle based on a deep and comprehensive investigation of the fault. The search parameters of the fault geometry are summarized in Table 2, and the slip distribution from the source model inversion is exhibited in Figure 5. The causative fault has a strike angle of N160°E and a dip angle of 90°, which are consistent with seismological solutions and previous research (Meneses-Gutierrez et al., 2019;USGS, 2016). As shown in Figure 5, the majority of fault slip occurs at shallow depths within a slip asperity with a maximum value of~1.1 m at a depth of 5 km. The top edge of the fault is located approximately 2 km below the surface, suggesting that the earthquake did not break the surface, which is consistent with geological observations (Nishimura & Takada, 2017). The geodetic moment released was 2.4 × 10 18 Nm, equivalent to Mw 6.2. Our moment estimate is consistent with the seismic moment inverted from local strong motion data (i.e., 2.1 × 10 18 Nm) (Kubo et al., 2017) but is lower than that estimated by the USGS-NEIC (i.e., 2.8 × 10 18 Nm) (USGS, 2016).
The 3-D coseismic deformation measurements observed by five GNSS stations are used to assess the accuracy of the InSAR-derived deformation. The locations of these GNSS stations can be found in Figure 3. The InSAR results are in good agreement with the GNSS results. The method used to extract the coseismic deformation from the GNSS data can be found in Text S1 and Figure S2. Figure 6 quantitatively compares the 3-D deformation results derived from GNSS and InSAR. The RMSEs between the GNSS and InSAR results are 0.4, 1.7, and 0.4 cm on the east, north, and up components, respectively. The north component shows a lower accuracy than the east and up components, which can be ascribed to the fact that the observing geometry of InSAR measurements is less sensitive to the north deformation component (Morishita et al., 2016;Wright et al., 2004).
In addition to generating a complete 3-D deformation field, the SM-RVCE method can estimate spatially continuous SM parameters. As demonstrated in the existing literature (Grafarend & Voosoghi, 2003;Guglielmino, Nunnari, et al., 2011;Hu, Wang, Li, et al., 2014), three dimensionless invariants, that is, dilatation, differential rotation magnitude, and maximum shear stain, can be extracted from the SM parameters (equations S1-S3, Figure S7). The 3-D strain invariants of the 2016 Central Tottori earthquake (dilatation, maximum shear strain, and differential rotation magnitude) estimated from the SM parameters can be found in Figure S8, and a detailed analysis can be found in Text S2 and Figures S7-S13.

Superiority of the SM-RVCE Method
The complete 3-D surface deformation field was retrieved from spatially fragmented InSAR measurements by the SM-RVCE method. Compared with the InSAR one-dimensional deformation along the LOS direction, the 3-D deformation field allows us to distinguish vertical and horizontal motions (Fialko et al., 2001) and therefore to better constrain the fault geometry Wang et al., 2015). The SM-RVCE method can extract the 3-D deformation with considerable reliability not only in the common zone (Zone B in Figure 1a) but also in other zones (Zones A and C in Figure 1a) and behaves well in resisting the propagation of errors from InSAR observations to the 3-D deformation. Moreover, the    The WLS method, which estimates the a priori variance of InSAR measurements in a moving window with a size of 5 × 5 pixels and solves the solution on a pixel-by-pixel basis (Jung et al., 2011), is commonly used to estimate 3-D deformation fields with InSAR measurements. However, if a set of spatially fragmented InSAR measurements are used, the WLS method cannot generate reliable 3-D deformation. In this study, if we want to obtain a more reliable deformation field, the 3-D deformation can only be generated in common Zone B by exploiting all four interferograms. On the other hand, if we want to obtain a more complete deformation field, the 3-D deformation should also be reconstructed in Zones A and C by exploiting three of the four interferograms. However, there is an unexpected boundary effect between different zones (see Figures 7a-7c). Two profiles (D-D' and E-E' in Figure 3) across these three zones are selected to exhibit this discrepancy. As shown in Figure 7d, the 3-D deformation along the selected profiles derived by the WLS algorithm shows obvious jumps at the image boundary, especially on the north component along profile D-D'. However, no jump occurs in the 3-D deformation derived by the SM-RVCE method.
Furthermore, to demonstrate that the SM-RVCE method can serve as more than a filter, we apply a median filter to the InSAR measurements with a window size of 1 × 1 km and estimate the 3-D deformation from the filtered measurements using the WLS method. The results are exhibited in Figure S14, where the deformation results along the two profiles derived by the WLS and SM-RVCE methods are compared. Although the 3-D deformation results estimated from the filtered observations are quite smooth, there are steep jumps in the boundary zones.
The performance of the SM-RVCE method is also compared with that of the SM-VCE method. The 3-D deformation fields estimated by these two methods are compared in Figure S15. Some obvious discrepancies can be found locally, which can be attributed to gross errors. To quantitatively evaluate the superiority of the SM-RVCE method with respect to the SM-VCE approach, we implement a simulated experiment using synthetic measurements with the same configuration as the case study. The results show that the SM-RVCE method can achieve improvements of approximately 31%, 54%, and 2% on the east, north, and up components, respectively, indicating the superiority of the robust VCE algorithm in mitigating the gross errors of InSAR observations. Details of the simulated experiment can be found in Text S2, Figures S9-S10, and Table S1.
Because the fault that generated the 2016 Central Tottori earthquake did not rupture the surface, the SM-RVCE algorithm can perform well. However, this method might fail for events with surface rupture, such as the 2016 Kaiköura earthquake in New Zealand . In such cases, we should select homogeneous points rather than all points within the fixed window when establishing the observation functions. One of the most straightforward ways to select homogeneous points is to remove the inhomogeneous points from the other side of the fault within the fixed window. The fault can be identified from deformation maps and optical images. Another possible solution is the SM, which can be transformed to represent the spatial correlation of the InSAR deformation of adjacent points. With the assumption that the observation errors follow a normal distribution, inhomogeneous points can be semiautomatically eliminated within the fixed window.

Precision Assessment by the SM-RVCE Method
Assessing the precision is still an open issue for InSAR measurements. External observations from independent techniques such as GNSS are not always available, and comparisons are somewhat unreasonable due to the spatial-temporal discrepancy between different measurements. In addition, it is basically impossible to Figure 6. Comparison of the 3-D deformation measurements derived from GNSS and InSAR. The red, green, and blue colors denote the east, north, and up components, respectively. The horizontal error bars denote 2 × SD for the GNSS deformation, and the vertical error bars denote 2 × SD for the InSAR deformation. estimate the a priori precision for InSAR measurements due to their numerous inherent errors associated with atmospheric artifacts, spatial-temporal decorrelations, and orbital ramps. The SM-RVCE method presented in this study provides an opportunity to conduct a precision assessment for InSAR measurements by an a posteriori approach.
The SD maps of the four InSAR observations are shown in Figure 8. The SD of most areas is at the millimeter level, but a few regions have an SD reaching several centimeters. The AsL observation has the lowest precision since the area is contaminated by gross errors, so its SD values range between approximately 1 and 2 cm. In addition, higher SD values are observed in the vicinity of the epicenter for all four observations. This is reasonable since a large deformation gradient is caused by the earthquake in this region. A comparison between the SD maps (i.e., Figure 8) and the InSAR observations (i.e., Figure 2) shows that the signal exhibits similar spatial distributions in some local areas of the two results. This is also expected since the signals of the original InSAR observations in these areas are contributed by errors rather than true deformation. In addition, a similar signal pattern can also be found to the north-northwest of the epicenter for the AsR and DesL observations (see Figure 8). This could be attributed to the similar sensitivity of these two observations to the 3-D deformation. Therefore, it is difficult to distinguish these observations, especially in areas with larger deformation gradients, where inhomogeneous measurements might be involved in the SM-RVCE method. Furthermore, the precision of the estimated 3-D deformation can also be calculated by equation (17). Unlike the 3-D deformation maps, the SD maps (i.e., Figure S16) show considerable jumps at the eastern boundary transition zones. Most likely, the large SD in the east is caused by the substantial noise in the AsL observation (see Figure 8).

Slip Models of the Tottori Earthquakes
The 2000 Mw 6.8 Tottori earthquake (Semmane et al., 2005) and the 2016 event are two recent M > 6 shallow inland earthquakes that struck Tottori Prefecture, Western Japan. Both events occurred in a high P-wave velocity zone in the upper crust (Zhao et al., 2018). Comparing the characteristics of these two events, we find that the fault size, dip angle, and orientation of the 2016 event are close to those of the 2000 event. The failure of both events started at great depth and propagated upward to shallow depth without breaking the surface. Similar to the 2000 event, the source parameters of the 2016 event have been studied using both seismic and geodetic data (Kubo et al., 2017). Although different studies obtained similar source parameters and slip distribution models for the 2016 event, they have some differences. An inversion of strong motion data (Kubo et al., 2017) suggests that a large slip region is located close to the hypocenter, which is consistent with our inversion results ( Figure 5). However, our geodetic slip inversion model does not show an obvious fault slip distribution to the north-northwest of the hypocenter. Studying the relocated aftershocks, Ross  2018) found a single large asperity of~16 km 2 with a maximum slip of 5.3 m close to the hypocenter. This estimated slip value is much greater than ours and that of Kubo et al. (2017). In terms of aftershocks, we find that the majority of M > 3 aftershocks following the 2016 event (see Figure 5) was distributed surrounding the maximum slip zone, while those of the 2000 event were concentrated beneath the coseismic fault slip zone at depths between 4 and 15 km (Semmane et al., 2005). These discrepancies probably reflect the heterogeneity and complexity of the fault system.

Conclusions
By exploiting ALOS-2 ascending and descending data with both left-and right-looking operation modes, we obtain the complete and accurate 3-D coseismic deformation field of the 2016 Central Tottori earthquake. The consideration of the left-looking observation mode greatly improves the rank of the coefficient matrix of the InSAR LOS observations, making it possible to infer 3-D deformation measurements without the noisy azimuth observations provided by the POT or MAI techniques. However, the footprints of the ascending images cover different portions of the seismic event, resulting in incomplete deformation fields or obvious boundary effects by classic WLS results. The SM-RVCE method used in this study utilizes the spatial relationship of the 3-D deformation between adjacent points to establish observation functions, and thus, the 3-D deformation can be estimated even in the case of spatially fragmented observations. In addition, the SM-VCE method is improved by introducing the robust VCE algorithm, making it more flexible when confronting the gross errors in InSAR observations. Furthermore, the employment of the VCE algorithm allows us to assess the InSAR observations' precision without using any external data. It should be noted that this method might fail if there is a discontinuity, such as a surface fault rupture, in the displacement field. In such cases, we should select homogeneous points rather than all points in the fixed window when establishing the observation functions. Moreover, to extract more accurate 3-D deformation, the SM-RVCE method can also be improved in the future using an irregular zone with a self-adaptive size according to the behavior of the ground deformation.
The complete and accurate 3-D coseismic deformation field estimated by integrating left-and right-looking InSAR observations with the SM-RVCE method enables a comprehensive analysis and interpretation of the 2016 Central Tottori earthquake. This reveals that the coseismic deformation zone can be divided into four sections. In the southeast and northwest sections, a pattern of uplift and epicenter-ward horizontal deformation is observed, while the opposite deformation pattern appears in the other two sections. By inverting the complete 3-D coseismic deformation field, we find that the source model has one large slip asperity extending from 10 km near the epicenter to 2 km with a maximum strike-parallel slip of~1.1 m.