Three-dimensional deformation and its uncertainty derived by integrating multiple SAR data analysis methods

Three-dimensional (3D) surface deformation data with high accuracy and resolution can help reveal the complex mechanisms and sources of subsurface deformation, both tectonic and anthropogenic. Detailed 3D deformation data are also beneficial for maintaining the position coordinates of existing ground features, which is critical for developing and advancing global positioning technologies and their applications. In seismically active regions, large earthquakes have repeatedly caused significant ground deformation and widespread damage to human society. However, the delay in updating position coordinates following deformation can hamper disaster recovery. Synthetic aperture radar (SAR) data allow high-accuracy and high-resolution 3D deformation measurements. Three analysis methods are currently available to measure 1D or 2D deformation: SAR interferometry (InSAR), split-bandwidth interferometry (SBI), and the pixel offset method. In this paper, we propose an approach to derive 3D deformation by integrating deformation data from the three methods. The theoretical uncertainty of the derived 3D deformations was also estimated using observed deformation data for each of these methods and the weighted least square (WLS) approach. Furthermore, we describe two case studies (the 2016 Kumamoto earthquake sequence and the 2016 Central Tottori earthquake in Japan) using L-band Advanced Land Observing Satellite 2 (ALOS-2) data. The case studies demonstrate that the proposed approach successfully retrieved 3D coseismic deformation with the standard error of ~ 1, ~ 4, and ~ 1 cm in the east–west, north–south, and vertical components, respectively, with sufficient InSAR data. SBI and the pixel offset method filled the gaps of the InSAR data in large deformation areas in the order of 10 cm accuracy. The derived standard errors for each pixel are also useful for subsequent applications, such as updating position coordinates and deformation source modeling. The proposed approach is also applicable to other SAR datasets. In particular, next-generation L-band SAR satellites, such as ALOS-4 and NASA-ISRO SAR (NISAR), which have a wider swath width, more frequent observation capabilities than the former L-band satellites, and exclusive main look directions (i.e., right and left) will greatly enhance the applicability of 3D deformation derivation and support the quick recovery from disasters with significant deformation.


Introduction
The ground surface can deform in three dimensions (3D) over a wide range of spatial and temporal scales and amplitudes owing to various causes ranging from large-scale natural phenomena, such as earthquakes and volcanic activity, to local landslides, surficial lateral spreading, and anthropogenic subsidence. Because it is difficult to directly observe underground conditions and motion, surface deformation data have been used to infer the complex mechanisms and sources of subsurface deformation. Although sparse point-wise deformation data have been widely used to estimate simple source models (e.g., the dislocation model in elastic halfspace (Okada 1985), the Mogi model (Mogi 1958), or a combination of the several source models), accurate and dense 3D surface deformation data should enable more complex and realistic source model derivations to better characterize subsurface Earth dynamics.
Positioning technology using the global navigation satellite system (GNSS) is already commonly used in applications such as smartphone and car navigation systems, and automated driving systems for cars and drones using positioning technology are becoming increasingly common. The more accurate the positioning system becomes, or the more the applications advance, the higher the required accuracy of the geospatial information of existing Earth surface features. When an earthquake occurs and these existing features are significantly displaced, their location information (i.e., 3D position coordinates) should be updated based on the magnitude of the displacement, for which accurate and high-resolution 3D deformation data are required. In seismically active countries like Japan, surface deformations of several meters have repeatedly been caused by large earthquakes (e.g., Himematsu and Furuya 2016;Ozawa et al. 2011). Every time a significant deformation occurs, the 3D position coordinates of control points (i.e., triangulation points and GNSS Continuously Operating Reference Stations (CORS)) and heights of benchmarks within the deformed area are updated (e.g., Hiyama et al. 2011;Nojiri et al. 2019;Ootaki et al. 2016). If the deformation is spatially smooth, the displacement at any point can be properly estimated by interpolation from the surrounding displacement data measured at the GNSS CORS without the need for field survey. In contrast, field surveys (i.e., GNSS or leveling) of all existing control points or benchmarks are required when or where complicated deformation cannot be retrieved owing to incomplete coverage by the GNSS CORS network. Field surveys require significant human resources, and the time delays involved in on-theground surveys can hamper quick recovery following a disaster. Indeed, following the 2011 Tohoku earthquake and the 2016 Kumamoto earthquake sequences, it took longer than 7 and 4 months, respectively, to update the position coordinates of the existing control points (Hiyama et al. 2011;Nojiri et al. 2019). It is even more difficult to quickly update position coordinates in countries and regions that are not covered by a dense GNSS CORS network. Moreover, to accurately measure displacement at control points in the field following an event, the control points require continuous maintenance leading up to the event, but this is not always easy. Furthermore, pointwise observations may miss unknown deformation signals between control points.
Synthetic aperture radar (SAR) data enable the highresolution measurement of 3D deformations. There are multiple methods for analyzing SAR data to derive deformation information, namely SAR interferometry (InSAR; e.g., Hanssen 2001), split-bandwidth interferometry (SBI; e.g., Bechor and Zebker 2006;Jiang et al. 2017;Mastro et al. 2020), and the pixel offset method (e.g., Michel et al. 1999). InSAR exploits the phase difference between two acquisitions and measures deformation with high accuracy in the line-of-sight (LOS) or range direction. In SBI, a double phase difference between band-split data and two acquisitions is computed and converted to deformation based on the fact that the phase change due to a location shift is proportional to the central frequency. SBI is applicable in both range and azimuth directions and is often called multiple aperture interferometry (MAI) in the azimuth. The pixel offset method measures the location shift of pixels in both range and azimuth from the cross-correlation between two SAR images. Because each method provides 1D or 2D deformation data, multiple observations from different directions (i.e., ascending/descending and right/left) need to be combined to retrieve the 3D deformation. Many studies have demonstrated 3D deformations derived from InSAR and other methods (e.g., Fialko et al. 2001;He et al. 2019aHe et al. , 2019bKobayashi et al. 2018); however, accuracy in the northsouth (NS) direction is typically much lower than in other directions (i.e., east-west (EW) and up-down (UD)) because the azimuth accuracy of MAI and the pixel offset method is much lower than that of InSAR range accuracy. While several empirical accuracy evaluations have been conducted based on comparisons with GNSS data, theoretical uncertainty has not been discussed in detail. Crucially, the uncertainty of 3D deformation derivation is important for updating and maintaining position coordinates and for subsequent deformation source modeling. Morishita et al. (2016) estimated the 3D deformation associated with a dike intrusion in Sakurajima volcano, Japan, using InSAR measurements observed from four different directions, including left-looking data. The accuracy of the 3D deformation derived from InSAR alone was higher than that with MAI or the pixel offset method, although the NS accuracy was still lower than that of the EW and UD directions because of the nearpolar orbit of the SAR satellite. However, InSAR-only approaches cannot measure large deformations (specifically in areas with large deformation gradients) because of decorrelation in interferograms. There have also been few opportunities for the further application of this approach because it requires left-looking acquisitions, which are only available in a very limited number of situations (Morishita 2019).
The L-band is more suitable for comprehensive deformation measurements (including in vegetated areas) than higher-frequency bands (i.e., C-or X-bands) because of its higher coherence (Rosen et al. 1996). Available spaceborne L-band SAR data are limited to the Japanese Earth Resources Satellite-1 (JERS-1;1992-1998, the Advanced Land Observing Satellite (ALOS;), and ALOS-2 (2014. However, new L-band satellites, such as ALOS-4 (Motohka et al. 2019(Motohka et al. , 2020 and NASA-ISRO SAR (NISAR; Jet Propulsion Laboratory 2020), which have wider swath widths and more frequent observation capabilities than the former L-band satellites, are planned to be launched and will provide abundant L-band SAR data. The main look directions of ALOS-4 and NISAR are right and left, respectively, which will enhance the opportunity for accurate 3D deformation derivation based on InSAR alone. A thorough investigation of 3D deformation analysis methods using existing L-band SAR data and their uncertainty is, therefore, important for observation planning, predicting measurement capabilities, and considering potential applications of these new satellites.
In this paper, we propose a method to derive accurate 3D deformation and its uncertainty by integrating multiple SAR analysis methods (i.e., InSAR, SBI, and the pixel offset method) using the weighted least square (WLS) approach. The approach is illustrated for two case studies of earthquakes in Japan using existing unique L-band ALOS-2 data. The uncertainty, contribution, and pros and cons of each method are discussed. The proposed method can be applied to the abundant SAR data that will be made available via next-generation satellites, enables the quick and accurate updating of geospatial information, and contributes to the rapid recovery following disaster events involving large ground deformation.

SAR interferometry (InSAR)
InSAR measures surface deformation along the range direction by computing the phase differences between two acquisitions (e.g., Hanssen 2001). To reduce speckle noise and improve coherence, multilooking is usually applied at the cost of spatial resolution. The coherence decreases as the time span between the acquisitions becomes long or when a change occurs in the surface scattering condition (e.g., by cultivation or snow cover). Decorrelation is also induced by a very large spatial gradient of deformation, which often makes it impossible to detect large deformations using InSAR. Both decorrelation effects are smaller in the L-band than in the Cand X-bands (Rosen et al. 1996;Morishita and Hanssen 2015a;He et al. 2019b). In contrast, ionospheric noise is much larger in the L-band, although the split-spectrum method (SSM) can correct most of the ionospheric phase and is almost essential for the L-band (Gomba et al. 2016;Wegmüller et al. 2018). Tropospheric noise can also be significant, but this can be mitigated to some extent using a numerical weather model or another empirical approach (Bekaert et al. 2015;Parizzi et al. 2021). Another significant error source is phase unwrapping error, the risk of which is lower in the L-band than at shorter wavelengths.

Split-bandwidth interferometry (SBI)
A SAR image has bandwidth with a central frequency both in range and azimuth, and subband images (e.g., higher and lower subbands) can be generated by a bandpass filter. Interferograms with each subband can then be produced from two acquisitions. Because the location shift induces a phase change, and its amount is proportional to the central frequency, the phase difference between two subband interferograms ( �φ H − �φ L ) divided by the central frequency difference ( f H − f L ) is proportional to the location shift (Scheiber and Moreira 2000;Jiang et al. 2017). More specifically, the displacements along the range and azimuth directions are denoted as c(�φ where c is the speed of light and v is the platform velocity (Mastro et al. 2020). These relationships are the basis of SBI (Jiang et al. 2017) and MAI in the azimuth (Bechor and Zebker 2006;Jung et al. 2009). The SBI interferogram is wrapped like the standard InSAR, but the effective wavelength is much longer, which makes phase unwrapping easy. For example, ALOS-2 Ultrafine (U) mode data with an 84-MHz bandwidth (as used in this study) split into 1/3 subbands has ~ 2.7 m/cycle in the case of SBI, whereas standard L-band InSAR has a ~ 12 cm/cycle.
The accuracy and spatial resolution of SBI are generally worse than those of InSAR, and both are vulnerable to temporal decorrelation. However, the decorrelation resulting from large deformation gradients in SBI is much milder than for InSAR because of its long effective wavelength. Therefore, SBI is more promising than InSAR for capturing large deformations. Another significant noise source in azimuth SBI is ionospheric azimuth shift, which may reach several meters (Wegmüller et al. 2006;Kim 2013). While mitigation of the azimuth shift is possible to some extent using ionospheric phases in the range estimated from the SSM, it is difficult to completely correct short-wavelength streaks (Liang and Fielding 2017;Yamashita et al. 2021).

Pixel offset
The pixel offsets (shifts) corresponding to the surface displacement in both the range and azimuth directions can be computed by incoherent cross-correlation between two SAR amplitude images (e.g., Michel et al. 1999). The accuracy of this depends on the spatial resolution of the SAR data, but is typically lower than that of InSAR data. Low-amplitude areas (i.e., smooth surfaces with few ground features) tend to have a low correlation and, therefore, sometimes return insignificant offset measurements. The median filter is generally applied to the derived offset field to reduce speckle noise or outliers. The spatial resolution is lower than that of InSAR because the cross-correlation is computed in a window that has some spatial width and must include a sufficient number of pixels. The pixel offset method is similar to SBI in some respects, and is suitable for detecting large deformations although vulnerable to ionospheric noise in the azimuth.

Theoretical error related to decorrelation
Decorrelation increases the measurement error in all three methods described in the previous sections. The theoretical error related to the decorrelation σ coh can be estimated from the absolute coherence γ and the number of independent samples (or effective looks) L , as follows: For InSAR, where is the wavelength (Rodriguez and Martin 1992;Hanssen 2001); For SBI, where B r is the ratio of the split sub-bandwidth to full bandwidth, assuming 1/3 as it brings the uncertainty closest to the Cramer-Rao bound, and p spa is the pixel spacing, which is inversely proportional to the bandwidth (Bamler and Eineder 2005;Jiang et al. 2017); For the pixel offset, (De Zan 2014). Here, note that γ is an incoherent cross-correlation. Figure 1 shows the dependency of coherence theoretical error for each method. The parameters are based on the real data used in this study, i.e., ALOS-2, U mode (84 MHz in range), beam number 7, 16 × 16 looks, L = 155, and p spa = 1.43 and 2.34 m in range and azimuth, respectively. In the pixel offset method, the theoretical error with 4 × L is also depicted by dashed lines in Fig. 1 because we used an overlapping window whose size is twice the number of looks. InSAR has an incomparably small error. The errors of SBI and the pixel offset method are relatively comparable at < 20 cm (< 0.1 pixels) if γ > 0.4, Page 5 of 19 Morishita and Kobayashi Earth, Planets and Space (2022) 74:16 and < 10 cm (< 0.05 pixels) if γ > 0.6. This is consistent with the error range reported in previous studies (Michel et al. 1999;Bechor and Zebker 2006). It is noteworthy that the theoretical error of SBI and the pixel offset in the unit of length (i.e., cm) is proportional to the spatial resolution (inversely proportional to the bandwidth), while that of InSAR is proportional to the wavelength of the SAR data. For example, in the Fine mode of ALOS-2, which has a 24 MHz bandwidth in range (i.e., 1/3 of the U mode) and is mainly used for world observations except for Japan, the range error in SBI and the pixel offset is three times that of the U mode. As this theoretical error is only associated with decorrelation, other significant error sources, such as the tropospheric and ionospheric delay, the unwrapping error in InSAR, and ionospheric azimuth shift in SBI and pixel offset should be considered separately (see Section "3D decomposition").

3D decomposition
The 3D orthogonal deformation components (i.e., EW, NS, and UD) can be retrieved from three or more independent 1D deformation datasets in different directions. Here, we follow the WLS approach proposed by Morishita et al. (2016) and Wright et al. (2004). Because the error ranges of InSAR and the other two methods are significantly different (see Section "Theoretical error related to decorrelation"), appropriate uncertainty estimation and weighting of each deformation data are essential for 3D decomposition.
The following operations were performed on a pixelby-pixel basis. First, suppose that we have N 1D deformation data D = [d 1 , . . . , d N ] T (either from InSAR, SBI, or the pixel offset in range or azimuth; note that upward in range and backward in azimuth are positive), a corresponding covariance matrix D , and EW, NS, and UD components of the unit vectors of each deformation data P = [p EW,1 , p NS,1 , p UD,1 ; . . . ; p EW,N , p NS,N , p UD,N ] , the WLS solution of the 3D deformation components x = [d EW , d NS , d UD ] T (i.e., D = Px + ε ) is as follows: and the covariance matrix of x is where the square root of the diagonal components ([σ d EW ,σ d NS ,σ d UD ]) is the standard error of the estimated 3D deformation.
Here, we assume D has only diagonal components (i.e., deformation data are independent of each other), and the variance of the deformation data can be computed by because the predominant error sources are long-wavelength atmospheric (tropospheric and ionospheric) noise and decorrelation (Akbari and Motagh 2012; Morishita and Hanssen 2015b). The σ coh for each method can be estimated from Eqs. (1, 2, 3). While it is difficult to estimate the σ atm at each pixel, approximate σ atm values representing the entire dataset can be simply calculated from the standard deviation of the spatially smoothed deformation data outside the deforming area (Morishita Page 6 of 19 Kobayashi Earth, Planets andSpace (2022) 74:16 et al. 2016). These procedures impose much more weight on accurate InSAR than the others, which enables us to jointly deal with deformation data with significantly different accuracies derived from the different methods. The decorrelated gap area of InSAR due to a large deformation gradient can be automatically filled by the joint use of SBI or the pixel offset method. The accuracy in the area without InSAR (i.e., only with SBI or the pixel offset) would be lower than in the area with InSAR, but the uncertainty of all pixels can be quantitatively obtained from Eq. (5) and exploited or considered for subsequent use (e.g., deformation source modeling or updating position coordinates). Equation (4) is overdetermined with N > 3, and the residuals of the 1D deformation data ε D can be computed as follows: Because the residual becomes large in inconsistent deformation data, it helps to identify the large noise signal included in certain deformation data (e.g., unwrapping errors, outliers, and long-wavelength noise) and correct or discard the noisy pixels or data. The residual of the deformation data with a small weight tends to increase. It may be difficult to identify truly noisy data from the residuals if N is not large enough and when weights are comparable because the residual is almost evenly distributed; a large N (redundancy) would prevent this situation.
One of the most challenging aspects of deformation measurement using SAR data is how to set a reference because the derived deformation data are always relative in a scene without an absolute reference. The most common approach is to set the reference area outside the deformation zone. Another point that needs to be carefully considered is the long-wavelength noise caused by orbit or coregistration errors, which can be approximated by a polynomial ramp (Hanssen 2001;Fattahi and Amelung 2014;Yamashita et al. 2021). This ramp is often estimated and removed by flattening the area without deformation. However, this may not work well if the non-deformation area is very small or anisotropic. There is also a risk of removing long-wavelength deformations together with long-wavelength noise. Although point-wise GNSS data, including in a deforming area, can be used for the deramp process (e.g., Kobayashi et al. 2011;Morishita et al. 2016Morishita et al. , 2018Lanari et al. 2020), this requires GNSS data and is not applicable in most areas in the world (Elliott et al. 2016).
Here, we propose an approach to tackle these two points simultaneously during 3D deformation derivation using redundant and independent deformation data. First, a common reference area was set for all deformation data. Note that while preferable, the reference area does not need to be deformation-free because the deformation in the reference area can be considered as a bias in the derived 3D deformation and, therefore, is arbitrarily retrieved after the 3D decomposition. Then, the 3D deformation and the residuals for each 1D deformation dataset are computed using Eqs. (4) and (7). Assuming that the ramp noise is random and other noise that leads to residuals (e.g., unwrapping error) is insignificant based on prior correction, the ramp noise can be estimated by a polynomial from the residuals. This includes the deforming area instead of flattening each data independently using only the non-deforming area, which reduces the risk of falsely removing any real longwavelength deformation as ramp noise. Once the ramp noise is estimated, the input 1D deformation data (and the corresponding weights) are updated by removing the ramp noise, and the 3D deformation and residuals are computed again. The ramp estimation and removal can be repeated until the residuals converge.
The derived 3D deformation may have significant noise for some pixels because of uncorrected noise or gaps in the input data. Such noisy pixels can be masked by the derived noise metrics, such as x and the root-meansquare (RMS) of ε D , using a threshold depending on the subsequent use and required accuracy of the derived 3D deformation.

Validation using real data 2016 Kumamoto earthquake sequence and 2016 Central Tottori earthquake
The Kumamoto earthquake sequence occurred in April 2016 in Kyushu, Japan (Fig. 2a). The sequence started with a foreshock with a Japan Meteorological Agency (JMA) magnitude (M j ) of 6.5 and a moment magnitude (M w ) of 6.2 on 14 April (12:26 UTC) followed by another foreshock (M j 6.4, M w 6.0) ~ 2.5 h later. These foreshocks mainly ruptured the Hinagu Fault and generated a ground deformation of > 15 cm (Kobayashi 2017). The M j 7.3 (M w 7.0) mainshock occurred on April 15 (16:25 UTC) and ruptured both the Futagawa and Hinagu Faults, generated > 2 m of deformation and numerous surface ruptures, induced liquefaction and damage to buildings and infrastructures at many sites, and caused > 270 earthquake-related deaths (Himematsu and Furuya 2016;Fujiwara et al. 2016;Wakamatsu et al. 2017a).
The M j 6.6 (M w 6.2) Central Tottori earthquake occurred on October 21, 2016, in the western part of Honshu, Japan (Fig. 2b). The associated damage was small compared with the Kumamoto earthquake sequence and, fortunately, resulted in no deaths (Amey et al. 2019). The largest displacement observed by the GNSS CORS Page 7 of 19 Morishita and Kobayashi Earth, Planets and Space (2022) 74:16 network for this event was approximately 7 cm at ~ 13 km north of the epicenter. As Fig. 2 shows, both regions are covered with dense vegetation; therefore, C-band SAR data tend to be decorrelated, while the L-band provides sufficient coherence (Jiang et al. 2017;Funning and Garcia 2018;Amey et al. 2019;Morishita 2019;He et al. 2019b). For both earthquake events, abundant pre-and post-earthquake ALOS-2 data are available including left-looking and right-looking types. Because the typical availability of L-band and left-looking data is very limited (see "Introduction" section), these two events provide an almost unique opportunity for 3D deformation studies.

Data processing
ALOS-2 has several observation modes with different bandwidths (spatial resolutions) and swath widths. In Japan, Stripmap U (84 MHz; 50 km swath) and ScanSAR Normal (14 or 28 MHz; 350 km swath) modes are mainly used for repeat basic observations. Here, we only used the U mode because the achievable accuracy of the Scan-SAR mode is much lower in comparison, although the swath width of the U mode is not sufficient to cover the area of interest by a single acquisition. In the forthcoming ALOS-4, the U mode will be increasingly applicable because its swath width will be increased by 200 km (Motohka et al. 2019(Motohka et al. , 2020. Therefore, the U mode is also appropriate for simulating the future potential of ALOS-4. We selected ALOS-2 U-mode interferometric pairs to cover the area of interest by both orbits (i.e., ascending and descending) and look directions (i.e., right and left) to achieve as high coherence as possible (Table 1). Multiple paths or beam numbers were necessary to cover the EW extent, but this will not be necessary for ALOS-4 or NISAR given their wider swath. Some of the pairs have long temporal baselines (~ 4 years at most) owing to limited data availability, especially in the left-looking direction. It should be noted that different acquisition dates allow the inclusion of different deformations other than the common coseismic one (i.e., pre-and post-seismic), especially in the Kumamoto earthquake sequence. However, pre-seismic deformation seems to be sufficiently small given that no significant deformation was observed Fig. 2 Optical images of the areas of interest in a Kumamoto and b Tottori. The orange and yellow rectangles denote footprints of the right-and left-looking SAR data, respectively. Note that the footprints do not correspond to the full observation area of the original dataset because the processed data were clipped to fit the target area. The solid and dashed white rectangles denote the enlarged areas shown in the following figures (except Fig. 6) and the excluded deformation area used for the reference, respectively. The white triangles and accompanying black arrows denote GNSS CORS and their horizontal coseismic displacement. The red crosses and lines denote the epicenter of the earthquakes and known active faults, respectively (Nakata and Imaizumi 2002) Morishita and Kobayashi Earth, Planets and Space (2022) 74:16 by the GNSS CORS network and no large seismic event occurred. Although ~ 10 cm post-seismic deformation in Kumamoto has been reported in previous studies (Fujiwara et al. 2016; Himematsu and Furuya 2020; Hashimoto 2020; Morishita 2021a), the deformed area was mostly limited where large coseismic deformation occurred and, therefore, the impact of the post-seismic deformation is relatively small. All three methods (i.e., InSAR, SBI, and pixel offset) were applied to each data pair and processed using GSISAR software (Tobita 2003;Morishita et al. 2018;Morishita 2019Morishita , 2021b. The 10 m mesh digital elevation model (DEM) prepared by the Geospatial Information Authority of Japan (GSI) was used to remove topographic phases in the InSAR. The multilook factor was 16 in both range and azimuth for all three methods. Tropospheric and ionospheric noise correction using the numerical weather model (Kobayashi 2016) and SSM (Gomba et al. 2016;Wegmüller et al. 2018) were applied, respectively, to successfully reduce the atmospheric noise in the InSAR (Additional file 1: Fig. S1). In the InSAR and SBI, a modified Goldstein filter based on coherence (Baran et al. 2003) was also applied to the wrapped interferograms followed by phase unwrapping using the statistical-cost network-flow algorithm for phase unwrapping (SNA-PHU) software version 2.0.4 (Chen and Zebker 2002).  Morishita and Kobayashi Earth, Planets and Space (2022) 74:16 Pixels with < 0.5 coherence of the filtered phase were masked before the unwrapping (Lazecký et al. 2020). In the InSAR, because very complicated phase discontinuities corresponding to surface ruptures existed northwest of Aso (Fujiwara et al. 2016) and caused many unwrapping errors that were too difficult to manually correct, this area was also masked. Other unwrapping errors in the InSAR data were identified from the residuals of the 3D decomposition (see Section "3D decomposition") and manually corrected. Unwrapping errors also existed in the SBI along the Futagawa Fault because of the large offset, but these were corrected with the aid of the pixel offset results. In the pixel offset method, the window size of the cross-correlation estimation was 32 in both range and azimuth, the data were oversampled by a factor of two during the cross-correlation estimation, pixels with < 0.2 correlation were masked, and a 2D median filter with a size of 7 (both in range and azimuth) was applied to the derived offset fields. Many azimuth deformation data from the SBI and the pixel offset method are affected by significant ionospheric noise. In particular, the ascending orbit (observed around midnight in local time) in summer often contains considerable azimuth streaks because medium-scale traveling ionospheric disturbances (MSTIDs) are highest during summer nights around Japan (Tsugawa et al. 2007;Morishita 2020). The azimuth ionospheric noise correction based on SSM (Yamashita et al. 2021) was, therefore, applied to reduce the noise to some extent, although significant noise (especially short-wavelength components) remained in some cases (Additional file 1: Fig. S2). The significantly noisy or decorrelated data (totaling 6 and 16 cases for Kumamoto and Tottori, respectively; Table 1) were excluded from the final 3D decomposition. All the data were geocoded into the common geographical grid with the pixel spacing of ~ 50 m.
The standard deviation of the smoothed deformation data based on the 2D Gaussian filter with a 1σ width of 0.5 ~ km was computed using the area excluding the expected deformation (Fig. 2) and used as σ atm (Table 1). As is mentioned above, the σ atm in the azimuth directions tend to be larger and more varied (3.5-34.7 cm) than those in the range directions (0.5-2.4 cm in InSAR, 1.1-7.9 cm in SBI and pixel offset) due to ionospheric noise. The reference area was set to the entire area, except for the expected deformation area (Fig. 2). The bilinear ramp for each dataset was estimated from the residual of the 3D decomposition and removed from the deformation data (see Section "3D decomposition"), which was repeated until the improvement in the RMS of the residual RMS of all pixels was < 0.5 mm (three times in these cases; Additional file 1: Figs. S3 and S4). Finally, some noisy pixels were masked based on the following noise metrics: [ σ d EW ,σ d NS ,σ d UD ] > [30 cm; 30 cm; 30 cm] and [1.5 cm, 6 cm, 1.5 cm] and RMS of ε D > 40 cm and 30 cm in Kumamoto and Tottori, respectively. These thresholds were determined from the frequency distribution of the noise metrics to retain most of the reliable pixels (Figs. 3  and 4). The thresholds in Kumamoto were larger than those in Tottori because they retained pixels without InSAR data (i.e., with only SBI or pixel offset) around the seismogenic faults that have relatively large x , whereas most of the pixels have InSAR data in Tottori.

Case study results
Clear 3D coseismic deformation and standard errors were obtained for both Kumamoto and Tottori (Figs. 3 and 4), and the derived 3D deformation is broadly consistent with previous studies (Himematsu and Furuya 2016;Morishita 2019;He et al. 2019b;Liu et al. 2019). The 3D deformation in the area with sufficient InSAR data was predominantly determined by InSAR with large weights, resulting in high accuracy. In Tottori, most of the areas had sufficient InSAR data and, therefore, almost uniform small standard errors (Fig. 4d-f ). In comparison, in Kumamoto, areas around the main seismogenic faults where large deformations occurred showed relatively large standard errors (Fig. 3d-f ) because the 3D deformation in these areas was derived from SBI and/or the pixel offset data. The modes of the standard errors in the EW, NS, and UD components were 1.4, 4.7, and 1.2 cm in Kumamoto (Fig. 3g-i), and 0.9, 3.8, and 0.7 cm in Tottori ( Fig. 4g-i), respectively. As expected, the NS uncertainty was higher than for the other components because the InSAR measurement directions (i.e., range) have a small NS contribution. Discontinuities may exist in the 3D deformation along the edge of each SAR data because the 1D deformation dataset differed inside and outside of the edge. However, such discrepancies were small (within the theoretical error range) along most edges (Figs. 3a-c, 4a-c, and Additional file 1: Fig. S5).

Comparison to external data
To evaluate the accuracy of the SAR-derived 3D deformation, we made a comparison with three types of external deformation data, specifically GNSS CORS, in situ offset measurements (Shirahama et al. 2016), and vertical displacement observed by a leveling survey (Ootaki et al. 2016). While GNSS CORS data are available for both Kumamoto and Tottori, the latter two data sources are only available for Kumamoto.
The coseismic displacement of the GNSS CORS was computed by subtracting the 15-day averaged pre-seismic coordinates from the 15-day averaged post-seismic coordinates (Figs. 2 and 5). The reference stations were selected outside of the deformation area. Because SAR and GNSS Page 10 of 19 Morishita and Kobayashi Earth, Planets and Space (2022) 74:16 Fig. 3 a-c 3D (i.e., EW, NS, and UD, respectively) deformation decomposed from the integration of InSAR, SBI, and pixel offset data masked based on noise metrics in Kumamoto. d-f The estimated standard errors of the 3D components. g-i Histograms of the standard errors of the 3D components. The light-gray parts denote frequency including the masked pixels Page 11 of 19 Morishita and Kobayashi Earth, Planets and Space (2022) 74:16 references were different, the average of the difference of the 3D deformation between them was considered as the reference bias. Therefore, we focused on the standard deviation of the difference. In Kumamoto, Choyo station (960701), which is located at the eastern edge of the Futagawa Fault, detected the largest deformation (~ 1 m) and showed much larger differences (~ 7 cm in the EW and UD and ~ 30 cm in NS components) compared to the other stations (< 3 cm in the EW and UD and < 7 cm in NS components). Cracks were found on the mortar basement of this station in a field survey on April 19, 2016, implying that this station may have been tilted by the earthquake sequence and, therefore, the detected displacement may be incorrect (Hiyama et al. 2016). Subsequently, this station was excluded from the comparison. The total number of GNSS stations used in Kumamoto and Tottori were 16 and
For the comparison with the in situ offset measurements (Shirahama et al. 2016), we computed the SAR-derived offset from the difference between two adjacent pixels to the in situ measurement locations. The orientation of the selected adjacent pixels depended on the strike directions (i.e., northwest and southeast adjacent pixels if the strike direction is northeast, and northeast and southwest adjacent pixels if the strike direction is northwest). The comparison shows generally good agreement considering the accuracy of the in situ measurements and the differences in the datasets (i.e., point-wise vs. pixel-wise; Fig. 6). For the masked SAR-derived 3D deformation, the RMS of the difference is 26.9 cm (n = 155) and 27.4 cm (n = 90) in the horizontal and vertical components, respectively. For the unmasked SAR-derived 3D deformation, the RMS of the difference is 35.0 cm (n = 196) and 27.1 cm (n = 101) in the horizontal and vertical components, respectively. Excluding five outliers around 32.81°N and 130.86°E that have remarkably large (> 1 m) discrepancies, the RMS for the unmasked deformation area was reduced to 29.5 cm (n = 191) in the horizontal component.
A leveling survey was conducted on the existing national benchmarks between June and August 2016 by the GSI (Fig. 7; Ootaki et al. 2016). As the previous surveys were conducted in 2003 and 2006, the observed vertical displacement includes contributions from the 10-or 13-year period before the 2016 earthquake sequence. The vertical displacements observed at GNSS CORS between 2003 and 2016 were < 2 cm; therefore, the 13-year pre-seismic vertical displacement appears to be insignificant in this area. Two benchmarks ~ 10 km west of the Hinagu Fault (around 32.73°N, 130.68°E) showed much larger subsidence (> 30 cm), and the discrepancy with SAR (> 12 cm) than the surrounding benchmarks (~ 20 cm and < 2 cm, respectively) may be because these two benchmarks were affected by liquefaction, which was common around these benchmarks (Wakamatsu et al. 2017b). Consequently, these two benchmarks were excluded from the comparison. Because the SAR and leveling references are different as well as GNSS, we focused on the standard deviation of the differences, which was 2.0 cm (n = 172) and is consistent with the expected uncertainty considering the pre-seismic vertical displacement in the leveling. All except six benchmarks were located in the area with sufficient InSAR data, and as expected, the differences for these six benchmarks were relatively high (3.7-6.3 cm). The magnitude of the difference shows no significant correlation with the magnitude of the displacement, implying that the measurement accuracy in the area with large displacements is comparable to the other areas.

Effect of weighting in the 3D decomposition
Without appropriate weight for the 3D decomposition (i.e., covariance matrix D ), the result becomes noisier because lower accuracy data (i.e., SBI and the pixel offset) are treated the same as higher-accuracy InSAR data, which degrades the final result (Additional file 1: Fig.  S7 and Table 2). The weighting enables the joint use of deformation data with differing accuracy. Moreover, the uncertainty of the 3D deformation cannot be derived without D , which is important for subsequent exploitation (e.g., masking noisy pixels, updating the position coordinates, and deformation source modeling).
Appropriate weighting also has an effect on the stability of the result even with noisy data, which would have a lower weight. Even in the case of data with significant decorrelation or azimuth ionospheric noise (Table 1, Additional file 1: Fig. S2), the results were not significantly degraded (Table 2, Additional file 1: Figs. S8 and S9). However, significantly noisy data should be excluded beforehand if possible because the RMS of the residual will be large, and masking noisy pixels using the metric is less robust (Additional file 1: Fig. S10).

Effect of iterative residual deramping
In Kumamoto and Tottori, conventional deramping (i.e., flattening the non-deformation area for each dataset independently) can be applied instead of the iterative deramping of residuals during the 3D decomposition (see Section "3D decomposition") because sufficient deformation-free areas are available in each dataset. In these case study events, the conventional deramping worked well and the results were broadly comparable (Table 2, Additional file 1: Figs. S11, and S12). However, some of the discrepancies along the data edges in the conventional deramping results were larger than those of the iterative deramping, especially in the NS component of Kumamoto (Additional file 1: Figs. S5b and S11b). This demonstrates that the proposed approach can achieve a more consistent deramp for all input data than the Page 14 of 19 Morishita and Kobayashi Earth, Planets and Space (2022)  conventional approach, which independently deramps data. Moreover, as noted in Section "3D decomposition", one of the most significant advantages of the proposed approach is that it is applicable even when a sufficient deformation-free area is unavailable because it deramps residuals that do not include deformation.

Characteristics, contributions, and limitations of each method
The 3D deformation was computed using only the data from each method (Additional file 1: Figs. S13-S18 and Table 2). In Kumamoto, the InSAR did not reveal large deformations along the seismogenic faults (Additional file 1: Fig. S13). The SBI results also tended to have low coherence and accuracy (Additional file 1: Fig. S14). The pixel offset method did detect large deformations without deterioration (Additional file 1: Fig. S15), but the true spatial resolution was lower than that of SBI because of the overlapping window and median filter (see Sections "Pixel offset" and "Theoretical error related to decorrelation"). The uncertainty from the SBI and the pixel offset method was in the order of 10 cm, which is much higher than that resulting from InSAR alone or joint use. Joint use with appropriate weighing is the best way to fill gaps and obtain the best possible accuracy for each pixel.
In Tottori, almost all significant deformations were captured by the InSAR with high accuracy (Additional file 1: Fig. S16). SBI and the pixel offset method poorly detected the ~ 10 cm coseismic deformation because of their lower levels of accuracy (Additional file 1: Figs. S17 and S18). The joint result was not significantly improved or degraded in the EW and UD components because of the large weight of the InSAR. However, the NS component was slightly improved (Table 2). Therefore, even in the case of Tottori, joint use including SBI and pixel offset provided the best result.
The U mode data which have the widest bandwidth of the L-band (84 MHz) were used in this study. As shown in Eqs. (2) and (3), σ coh is inversely proportional to the bandwidth. Therefore, the uncertainty of SBI and pixel offset would become worse than ~ 10 cm using other observation modes with a narrower bandwidth.
For the NS component, an accuracy comparable to that of the EW and UD components was not achievable using SAR data alone. Inhomogeneous accuracy between the EW and NS components will result in unexpectedly large uncertainty in the orientation angle of the horizontal deformation unless the NS deformation is significantly larger than this uncertainty, and this can lead to misinterpretation (see the area outside the deformation area in Fig. 5b). As such, quantitative standard error data are important to prevent misinterpretations.

Conclusions
We proposed a 3D decomposition method to jointly use InSAR, SBI, and pixel offset data to consistently remove polynomial ramp noise from residuals and estimate the standard error of the solution using uncertainty estimations of observed deformation data and the WLS approach. The derived 3D deformation has the standard error of ~ 1, ~ 4, and ~ 1 cm in the EW, NS, and UD components, respectively, where sufficient InSAR data are available, with no significant gaps even in high-deformation areas owing to SBI and the pixel offset method, although the accuracy is in the order of 10 cm using the highest-resolution U mode data. Two case studies (the 2016 Kumamoto earthquake sequence and 2016 Central Tottori earthquake) demonstrated that the proposed method can successfully retrieve accurate 3D coseismic deformations with theoretical standard errors comparable to the levels of uncertainty resulting from a comparison with external GNSS, in situ offset measurement, and leveling data. Importantly, the derived standard errors for each pixel are appropriate for subsequent use in 3D deformation studies, such as the updating of position coordinates or deformation source modeling.
L-band ALOS-2 U mode data with a 50-km swath width were used, and the temporal baselines were not short because of infrequent acquisitions, which led to significant temporal decorrelation in some cases (Table 1). As the data from multiple beam numbers and paths were necessary because each dataset was not sufficiently wide to cover the entire area of interest ( Fig. 1 and Table 1), some jumps occurred along the edges in the derived 3D deformation (Additional file 1: Fig. S5). However, next-generation L-band SAR satellites such as ALOS-4 and NISAR will offer a much wider swath width (~ 200 and ~ 240 km, respectively) and more frequent acquisitions than ALOS-2; these satellites will provide coseismic deformation data with a shorter temporal baseline and better coherence, and will be able to capture larger areas of interest in a single acquisition. This offers considerable potential for seamless 3D deformation derivation. Latency will also be greatly improved by a wider swath and more frequent observations. These technological developments will enable complete observations from four different directions (i.e., ascending/descending and right/left) within 14 days of a seismic event, which combined with the proposed method, will allow accurate and rapid 3D deformation analysis following major tectonic disasters.