3D Displacement Field of Wenchuan Earthquake Based on Iterative Least Squares for Virtual Observation and GPS / InSAR Observations

: The acquisition of a 3D displacement ﬁeld can help to understand the crustal deformation pattern of seismogenic faults and deepen the understanding of the earthquake nucleation. The data for 3D displacement ﬁeld extraction are usually from GPS / interferometric synthetic aperture radar (InSAR) observations, and the direct solution method is usually adopted. We proposed an iterative least squares for virtual observation (VOILS) based on the maximum a posteriori estimation criterion of Bayesian theorem to correct the errors caused by the GPS displacement interpolation process. Firstly, in the simulation examples, both uniform and non-uniform sampling schemes for GPS observation were used to extract 3D displacement. On the basis of the experimental results of the reverse fault, the normal fault with a strike-slip component, and the strike-slip fault with a reverse component, we found that the VOILS method is better than the direct solution method in both horizontal and vertical directions. When a uniform sampling scheme was adopted, the percentages of improvement for the reverse fault ranged from 3% to 9% and up to 70%, for the normal fault with a strike-slip component ranging from 4% to 8% and up to 68%, and for the strike-slip fault with a reverse component ranging from 1% to 8% and up to 22%. After this, the VOILS method was applied to extract the 3D displacement ﬁeld of the 2008 Mw 7.9 Wenchuan earthquake. In the East–West (E) direction, the maximum displacement of the hanging wall was 1.69 m and 2.15 m in the footwall. As for the North–South (N) direction, the maximum displacement of the hanging wall was 0.82 m for the southwestern, 0.95 m for the northeastern, while that of the footwall was 0.77 m. In the vertical (U) direction, the maximum uplift was 1.19 m and 0.95 m for the subsidence, which was signiﬁcantly di ﬀ erent from the direct solution method. Finally, the derived vertical displacements were also compared with the ruptures from ﬁeld investigations, indicating that the VOILS method can reduce the impact of the interpolated errors on parameter estimations to some extent. The simulation experiments and the case study of the 3D displacement ﬁeld for the 2008 Wenchuan earthquake suggest that the VOILS method proposed in this study is feasible and e ﬀ ective, and the degree of improvement in the vertical direction is particularly signiﬁcant.


Introduction
Interferometric synthetic aperture radar (InSAR) technology has become one of the methods for monitoring surface deformation due to its advantages of large spatial coverage, day-and-night

Mathematical Background
The LOS deformation obtained by InSAR technology does not represent the actual 3D surface deformation, but the projection of the East-West (E), North-South (N), and vertical (U) deformation in the LOS direction [39]. As shown in Figure 1, H is the projection of satellite flight direction on the ground, φ is the azimuth of satellite flight direction (positive clockwise from the North), and θ is the radar incidence angle at the reflection point. U e , U n , and U u are the displacements of three directions Remote Sens. 2020, 12, 977 3 of 17 of E, N, and U respectively. The relationship between the observed values L InSAR and the surface deformation in these three directions can be expressed as [40] L InSAR = − sin θ sin(φ − 3π/2) − sin θ cos(φ − 3π/2) cos θ U e U n U u T (1) The LOS deformation obtained by InSAR technology does not represent the actual 3D surface deformation, but the projection of the East-West (E), North-South (N), and vertical (U) deformation in the LOS direction [39]. As shown in Figure 1, H is the projection of satellite flight direction on the ground,  is the azimuth of satellite flight direction (positive clockwise from the North), and  is the radar incidence angle at the reflection point. e U , n U , and u U are the displacements of three directions of E, N, and U respectively. The relationship between the observed values InSAR L and the surface deformation in these three directions can be expressed as [40]    The likelihood function between the InSAR LOS displacements InSAR L and the 3D displacement X can be described as follows: where InSAR D is the absolute value of the determinant of the variance matrix InSAR D . Equation (1) is the theoretical equation, when L InSAR are the observations, and their measurement errors are included at the same time. Let S e = − sin θ sin(φ − 3π/2), S n = − sin θ cos(φ − 3π/2) and S u = cos θ, then Equation (1) can be written as , B = I n ⊗ S e S n S u , i is the ith number of the points (i = 1, 2, · · · n).
⊗ denotes the Kronecker product, and the projection vector S e S n S u denotes the average of the all points. The likelihood function between the InSAR LOS displacements L InSAR and the 3D displacement X can be described as follows: where |D InSAR | is the absolute value of the determinant of the variance matrix D InSAR . The GPS stations are interpolated into the spatial density of the LOS observations, and they are regarded as the virtual observations to constrain the 3D displacement field. The expression is showed as Thus, the prior information constrained on the 3D displacement can be presented by a probability density function (PDF) as , and |D GPS | is the absolute value of the determinant of the variance matrix D GPS .
According to the Bayesian theorem [35], the posterior PDF for the 3D displacement is calculated as follows: where the denominator is a normalizing constant independent of X, which is set as nc. Thus, substituting Equations (4) and (5) into Equation (6): where Based on the maximum a posteriori estimation criterion p(X L InSAR ) = min , which is equivalent to whereX is the estimation of the X. According to the principle of generalized least squares adjustment [42], we can obtain the expression of the 3D displacementX

The VOILS Method and its Algorithmic Flow
Although GPS points can coincide with some observation points of InSAR, the spatial density of GPS data points is still not enough, therefore it is necessary to obtain the spatial density consistent with InSAR data points by interpolation. Consequently, this paper uses the ordinary Kriging interpolation method [43] to interpolate GPS points and obtain the prior initial 3D displacement field at the observation points of InSAR, which is as shown in Equation (5). The 3D displacement field obtained by interpolation inevitably contains errors, therefore, this study designs an iterative least squares for virtual observation, which has the advantages of correcting errors in interpolated GPS displacement and reasonably fusing two types of data. The iteration process is the key step that is different from the traditional method, which is expected to realize the main objectives as described above. Figure 2 shows the flow chart of the VOILS method. It can be described as: (1) Downsampling InSAR data to obtain sparse LOS displacements; (2) Interpolating GPS data to the spatial density of InSAR downsampled data by the Kriging method as virtual observations; (3) Calculating E, N, and U component's displacement by the least squares for virtual observation (i.e., Equation (9)); (4) If the maximum absolute value of the difference between the front and back of the 3D displacement is less than the given threshold value δ, the parameterX i is the output and the iteration will be terminated, otherwise the virtual observation values E, N, and U will be updated, and the steps 2~4 will repeat.
Remote Sens. 2020, 12, 977 5 of 17 component's displacement by the least squares for virtual observation (i.e., Equation (9)); 4) If the maximum absolute value of the difference between the front and back of the 3D displacement is less than the given threshold value δ , the parameter ˆi X is the output and the iteration will be terminated, otherwise the virtual observation values E, N, and U will be updated, and the steps 2 ~ 4 will repeat.

Simulation Experiment
We adopt the following strategies (i.e., uniform and non-uniform) to select GPS observation points. In Scheme 1, we sample the GPS observation points with the uniform mode. To investigate the impact of the different density on the extracted 3D displacement field, 9, 25, 49, and 121 GPS points are selected respectively, which correspond to the spatial resolutions of 30 km, 20 km, 14 km, and 10 km of GPS observation points. In Scheme 2, we employ the non-uniform way to sample the GPS observation points in the displacement field. We use the randperm function in matrix laboratory (MATLAB) to select 9, 25, 49, and 121 points randomly according to the number of uniformly sampled GPS points.
In the comparative analysis of the results, two indicators, i.e., the root mean square error (RMSE) and the percentage of improvement of VOILS method compared with the direct solution method (Per), are used to evaluate the performances of the VOILS method. 1) RMSE is calculated by the difference d between the fitting value at each grid point and the original simulated true value. The expression is as follows: where n means the number of the grid points.
2) Per, the percentage value of improvement, is presented as follows

Simulation Experiment
We adopt the following strategies (i.e., uniform and non-uniform) to select GPS observation points. In Scheme 1, we sample the GPS observation points with the uniform mode. To investigate the impact of the different density on the extracted 3D displacement field, 9, 25, 49, and 121 GPS points are selected respectively, which correspond to the spatial resolutions of 30 km, 20 km, 14 km, and 10 km of GPS observation points. In Scheme 2, we employ the non-uniform way to sample the GPS observation points in the displacement field. We use the randperm function in matrix laboratory (MATLAB) to select 9, 25, 49, and 121 points randomly according to the number of uniformly sampled GPS points.
In the comparative analysis of the results, two indicators, i.e., the root mean square error (RMSE) and the percentage of improvement of VOILS method compared with the direct solution method (Per), are used to evaluate the performances of the VOILS method.
(1) RMSE is calculated by the difference d between the fitting value at each grid point and the original simulated true value. The expression is as follows: where n means the number of the grid points.
(2) Per, the percentage value of improvement, is presented as follows where RMSE New and RMSE Old represent the RMSE obtained by the VOILS method and the direct solution method, respectively. The uniform elastic half-space dislocation theory [44] is used to simulate the 3D surface displacement field. The deformation area is 100 km × 100 km. The model parameters of reverse fault, normal fault with a strike-slip component, and strike-slip fault with a reverse component are listed in Table 1. The simulated 3D displacement field is calculated by FORTRAN programs EDGRN/EDCMP [45]. The simulation field is divided into 51 × 51 grid cells with a 2km interval (see Figure 3, in which both horizontal and vertical coordinates range from −50 km to 50 km, and the projection coefficient [26]). The random errors adhering to normal distribution whose standard deviations are 3 mm, 5 mm, and 30 mm are added to the horizontal component, vertical component, and LOS deformation, respectively. Twenty-five selected points are displayed in the upper right corner subfigure of Figure 3.  The simulation field is divided into 51×51 grid cells with a 2km interval (see Figure 3, in which both horizontal and vertical coordinates range from -50 km to 50 km, and the projection coefficient is [Se Sn Su]= [0.340 -0.095 0.935] [26]). The random errors adhering to normal distribution whose standard deviations are 3 mm, 5 mm, and 30 mm are added to the horizontal component, vertical component, and LOS deformation, respectively. Twenty-five selected points are displayed in the upper right corner subfigure of Figure 3.

Experiment with the Reverse Fault
In Scheme 1, we carried out 200 simulation experiments. The RMSE values of the corresponding 3D displacement field are calculated with these two methods, and the results are counted to get the histogram of Figure 4. It shows the RMSE values of 3D displacement field with different GPS spatial resolutions. It can be seen that the histogram obtained by the VOILS method is closer to the origin of coordinate than that computed by the direct solution method, which implies that the VOILS method exhibits some improvement relative to the direct solution method and can retrieve the 3D displacement field better.
corresponding 3D displacement field are calculated with these two methods, and the results are counted to get the histogram of Figure 4. It shows the RMSE values of 3D displacement field with different GPS spatial resolutions. It can be seen that the histogram obtained by the VOILS method is closer to the origin of coordinate than that computed by the direct solution method, which implies that the VOILS method exhibits some improvement relative to the direct solution method and can retrieve the 3D displacement field better.  In order to quantitatively illustrate the advantages of the VOILS method and form a more intuitive comparison with the direct solution method, mean improvement percentages for Scheme 1 are listed in Table 2, which is calculated by Formula (11). It shows that the improvements in all directions are comparable under different GPS spatial resolutions. The improvement degrees of the three directions from high to low are U, E, and N, which may be related to the satellite projection coefficient. The projection coefficients of U, E, and N directions are 0.935, 0.340, and −0.095, respectively. Thus, the vertical direction contributes the most to the LOS deformation, followed by the E and N directions.
The experiments under different GPS spatial resolutions were also carried out. The improved percentages are listed in Table 1. In addition, the corresponding number of GPS points are listed in brackets. It can be seen that the percentages of improvement display a relatively stable state with the increasing number of GPS points, which indicates that the proposed method has the advantage of not being subject to GPS spatial resolution. Compared with the direct solution method, the percentages of Remote Sens. 2020, 12, 977 8 of 17 improvement for horizontal directions (i.e., E and N directions) range from 3% to 9%. Unexpectedly, the vertical direction (i.e., U direction) improvement percentage exceeds 68%. Two hundred simulation experiments were also carried out in Scheme 2, and Per averages are shown in Table 2. It can be seen that the experiment results from Scheme 2 are similar to those from Scheme 1. Moreover, we executed several random experiments with non-uniformly sampled GPS points and calculated the Per. As shown in Figure 5, the proposed method in three directions exhibits different degrees of improvement relative to the direct solution method. Evidently, the vertical direction displays obvious improvement. It can be observed that the experimental results show a relatively stable state with the variations of the random sampling GPS points.

Experiments with the Normal Fault with a Strike-Slip Component and Strike-Slip Fault with a Reverse Component
In the above experiments, the fault slip is assumed to be the reverse fault. The normal fault with a strike-slip component and strike-slip fault with a reverse component will be considered in this section further.
Firstly, we used Scheme 1 and Equation (11) to compute the Per, as shown in Table 3. It can be observed that for different types of fault movement, the improvement ratios in different directions under different GPS spatial resolutions are similar. The improvement degrees of three directions from high to low are vertical, E, and N. The results of Table 3 are similar to the ones of Table 2. In the case study of the normal fault with a strike-slip component, the improvement ratios of the horizontal direction range from 4% to 8%. Additionally, those of the vertical direction are up to 68%. As for the strike-slip fault with a reverse component, the results show that different directions in improvement ratios demonstrate significant differences. It can be seen that vertical direction has larger improvement ratios exceeding 21%, while lower improvement ratios below 2% occur in the

Experiments with the Normal Fault with a Strike-Slip Component and Strike-Slip Fault with a Reverse Component
In the above experiments, the fault slip is assumed to be the reverse fault. The normal fault with a strike-slip component and strike-slip fault with a reverse component will be considered in this section further.
Firstly, we used Scheme 1 and Equation (11) to compute the Per, as shown in Table 3. It can be observed that for different types of fault movement, the improvement ratios in different directions under different GPS spatial resolutions are similar. The improvement degrees of three directions from high to low are vertical, E, and N. The results of Table 3 are similar to the ones of Table 2. In the case study of the normal fault with a strike-slip component, the improvement ratios of the horizontal direction range from 4% to 8%. Additionally, those of the vertical direction are up to 68%. As for the strike-slip fault with a reverse component, the results show that different directions in improvement ratios demonstrate significant differences. It can be seen that vertical direction has larger improvement ratios exceeding 21%, while lower improvement ratios below 2% occur in the N direction. Despite that, the VOILS method is better than the direct solution method.  Table 3. We found that the results are similar to Table 2. Figure 6 portrays the histogram of RMSE values with a normal density curve at nine non-uniformly sampled GPS points. In addition, Figure 7 shows the improvement ratio of the VOILS method compared with the direct solution method. Obviously, RMSE values calculated by the VOILS method are smaller than those computed by the direct solution method, which indicates that the VOILS method is better than the direct solution method.  Table 3. We found that the results are similar to Table 2. Figure 6 portrays the histogram of RMSE values with a normal density curve at nine non-uniformly sampled GPS points. In addition, Figure 7 shows the improvement ratio of the VOILS method compared with the direct solution method. Obviously, RMSE values calculated by the VOILS method are smaller than those computed by the direct solution method, which indicates that the VOILS method is better than the direct solution method.  In summary, it is concluded that the improvement degree in the vertical direction is better than that in the horizontal direction. The magnitude of vertical deformation calculated by the VOILS method is closer to the true value than that computed from the direct solution method. This confirms that the VOILS method can make full use of InSAR observations. Our simulation results fully show that the VOILS method is feasible and effective. Furthermore, the proposed method displays obvious advantages in retrieving 3D displacement fields compared to the direct solution method, especially for vertical deformation.

GPS and InSAR Data
The Wenchuan earthquake occurred on May 12, 2008, on the Longmenshan fault zone in Sichuan Province (Figure 8). The earthquake caused surface rupture in the Yingxiu-Beichuan fault of the main central fault and the Guanxian-Jiangyou fault of the Qianshan fault [46]. Previous studies have illustrated that this event is a complex slip mechanism as a variable combination of a reverse and right-lateral component [46][47][48][49][50]. After the earthquake, the co-seismic deformation observation data, including GPS data and InSAR data from advanced land observation satellite (ALOS) satellite phased array type L-band synthetic aperture radar (PALSAR) images, were obtained. In this paper, we employed the downsampled 3792 LOS displacements from Xu et al. [47], and 473 GPS co-seismic displacement data calculated by Wang et al. [5]. It is noted that 297 GPS points (i.e., 284 campaigned stations and 13 continuous stations) are used in our study, with the spatial distribution of the data shown in Figure 8. In summary, it is concluded that the improvement degree in the vertical direction is better than that in the horizontal direction. The magnitude of vertical deformation calculated by the VOILS method is closer to the true value than that computed from the direct solution method. This confirms that the VOILS method can make full use of InSAR observations. Our simulation results fully show that the VOILS method is feasible and effective. Furthermore, the proposed method displays obvious advantages in retrieving 3D displacement fields compared to the direct solution method, especially for vertical deformation.

GPS and InSAR Data
The Wenchuan earthquake occurred on May 12, 2008, on the Longmenshan fault zone in Sichuan Province (Figure 8). The earthquake caused surface rupture in the Yingxiu-Beichuan fault of the main central fault and the Guanxian-Jiangyou fault of the Qianshan fault [46]. Previous studies have illustrated that this event is a complex slip mechanism as a variable combination of a reverse and right-lateral component [46][47][48][49][50]. After the earthquake, the co-seismic deformation observation data, including GPS data and InSAR data from advanced land observation satellite (ALOS) satellite phased array type L-band synthetic aperture radar (PALSAR) images, were obtained. In this paper, we employed the downsampled 3792 LOS displacements from Xu et al. [47], and 473 GPS co-seismic displacement data calculated by Wang et al. [5]. It is noted that 297 GPS points (i.e., 284 campaigned stations and 13 continuous stations) are used in our study, with the spatial distribution of the data shown in Figure 8.

Results and Comparative Analysis
The 3D displacement field of the Wenchuan earthquake obtained by the VOILS method is shown in Figures 9a-c

Results and Comparative Analysis
The 3D displacement field of the Wenchuan earthquake obtained by the VOILS method is shown in Figure 9a Figure 9a shows the E direction deformation, and it can be observed that the hanging wall moves eastward with a maximum displacement of 1.69 m, while the footwall moves westward with a maximum displacement of 2.15 m, showing a right-lateral trend. As for the N direction deformation (Figure 9b), the southwest of the hanging wall moves southward with a maximum displacement of 0.82 m and the northeast section moves northward with a maximum displacement of 0.95 m. The whole footwall moves northward with a maximum displacement of 0.77 m, and the deformation decreases in the northeast direction. The relatively large displacements of the hanging wall and footwall near the epicenter indicate that the reverse component in this area is the largest. Finally, the U direction deformation is shown in Figure 9c. The hanging wall and footwall near the fault zone appear uplift and subsidence. In addition, the maximum uplift and the maximum subsidence near the epicenter are 1.19 m and 0.95 m respectively, which means that reverse movement existed in the southwest of the fault. The 3D deformation characteristics clearly show the local movement characteristics of the fault and reflect the seismogenic fault's complexity and heterogeneity.
From Figure 9d-e, we know that the LOS displacement field is close to that of Figure 8. Meanwhile, the unit for fitting residual errors after iteration calculation is millimeters, which is at the same order magnitude of the threshold δ (2 mm) adopted in the iteration process. It is useful to compare the derived coseismic displacement field with independent GPS observations. The InSAR points with a distance of approximately 3 km around the GPS points were chosen for horizontal comparison. Figure 10 shows the horizontal displacement vectors, indicating that the agreement between the GPS observations and the derived displacements is quite good. The RMSE values between them are 6.1 cm, 2.3 cm, and 5.8 cm in the E, N, and U directions, respectively. In order to compare with the above results, the 3D displacement field (Figures 11a-c) was calculated using the direct solution method. The LOS deformation and its residual errors are shown in Figures 11d-e. Figures 11a-c show similar features as Figures 9 a-c, while the numerical values are different. It can be seen from Figure 11e that the residuals along both sides of the seismogenic fault are larger, especially in the southwest of the hanging wall. This can be ascribed to the fault being close to the epicenter, which caused larger deformation, as well as poor interpolation accuracy with less near-field GPS data. This is similar to the characteristic distribution of the LOS deformation residual map produced by Luo et al. [34]. However, the numerical values are different, which may be related to the data and the method used. In order to compare with the above results, the 3D displacement field (Figure 11a-c) was calculated using the direct solution method. The LOS deformation and its residual errors are shown in Figure 11d-e. Figure 11a-c show similar features as Figure 9a-c, while the numerical values are different. It can be seen from Figure 11e that the residuals along both sides of the seismogenic fault are larger, especially in the southwest of the hanging wall. This can be ascribed to the fault being close to the epicenter, which caused larger deformation, as well as poor interpolation accuracy with less near-field GPS data. This is similar to the characteristic distribution of the LOS deformation residual map produced by Luo et al. [34]. However, the numerical values are different, which may be related to the data and the method used.
In order to compare with the above results, the 3D displacement field (Figures 11a-c) was calculated using the direct solution method. The LOS deformation and its residual errors are shown in Figures 11d-e. Figures 11a-c show similar features as Figures 9 a-c, while the numerical values are different. It can be seen from Figure 11e that the residuals along both sides of the seismogenic fault are larger, especially in the southwest of the hanging wall. This can be ascribed to the fault being close to the epicenter, which caused larger deformation, as well as poor interpolation accuracy with less near-field GPS data. This is similar to the characteristic distribution of the LOS deformation residual map produced by Luo et al. [34]. However, the numerical values are different, which may be related to the data and the method used.  The results indicate that the deformation directions along both sides of the fault are basically opposite, and the larger displacements in the three directions are nearer to the source. Near the epicenter of the E direction, the movement of the hanging wall to the east, relative to the footwall, is dominant, and along the NE direction, the deformation of the hanging wall and footwall is opposite. In the N direction, the northward movement of the northeastern part of the hanging wall is dominant, with a small amount of reverse component. Both the hanging wall and footwall of the U direction have subsidence near the fault zone, and the uplift of the hanging wall is greater than the subsidence of the footwall. In summary, we can understand that the main fault near the epicenter is a reverse fault with a right-lateral strike-slip component, while the northeastern segment is a right-lateral strike-slip with a small amount of reverse component, which is in line with the results of other research [46][47][48][49][50].
In order to further verify the degree of improvement in the vertical direction, we compared the derived displacements to the surface rupture measurements from Xu et al. (2009) [50]. Figure 12 shows the comparison of these data. Note that the average displacements of InSAR points with a distance of about 5 km around the field points are used for comparison. It shows that the vertical displacements from the VOILS method are more consistent with the field observations, which indicates that the VOILS method can effectively reduce the impacts of the interpolated errors on parameter estimation to some extent. Meanwhile, these results indicate that the VOILS method performs better than the direct solution method in this event. Figure 11e shows that the residual magnitude of InSAR is different from that of Figure 9e. The residual magnitude of InSAR fitted by VOILS method is mm-level, which is significantly smaller than that of the direct solution method. Meanwhile, the results of this study are also different from those of the existing researches [28,[36][37][38]. Therefore, the results obtained by the VOILS method are relatively stable and the fitting is relatively good. On the other hand, the RMSE value of the direct solution method is 0.07m, representing a precision of cm-level. The above analyses show that the VOILS method can make full use of the InSAR observations, and has the ability to exert benefits from GPS and InSAR data. distance of about 5 km around the field points are used for comparison. It shows that the vertical displacements from the VOILS method are more consistent with the field observations, which indicates that the VOILS method can effectively reduce the impacts of the interpolated errors on parameter estimation to some extent. Meanwhile, these results indicate that the VOILS method performs better than the direct solution method in this event.  [50]. The black and red lines represent the displacements by making subtraction between the derived deformation on hanging wall and footwall from the VOILS method and the direct solution method, respectively.). Figure 11e shows that the residual magnitude of InSAR is different from that of Figure 9e. The residual magnitude of InSAR fitted by VOILS method is mm-level, which is significantly smaller than that of the direct solution method. Meanwhile, the results of this study are also different from those of the existing researches [28,[36][37][38]. Therefore, the results obtained by the VOILS method are relatively stable and the fitting is relatively good. On the other hand, the RMSE value of the direct solution method is 0.07m, representing a precision of cm-level. The above analyses show that the VOILS method can make full use of the InSAR observations, and has the ability to exert benefits from GPS and InSAR data.

Conclusions
In the conventional direct solution method, GPS data need to be interpolated first, which will undoubtedly introduce errors. Therefore, this paper proposes the VOILS method, which can correct the errors in GPS 3D deformation interpolation, and realize the reasonable fusion of GPS/InSAR Figure 12. Comparison of the derived vertical displacements with the surface rupture measurements (The blue bar represents the measured rupture from field investigations (Xu et al., 2009) [50]. The black and red lines represent the displacements by making subtraction between the derived deformation on hanging wall and footwall from the VOILS method and the direct solution method, respectively.).

Conclusions
In the conventional direct solution method, GPS data need to be interpolated first, which will undoubtedly introduce errors. Therefore, this paper proposes the VOILS method, which can correct the errors in GPS 3D deformation interpolation, and realize the reasonable fusion of GPS/InSAR observations in 3D displacement extraction based on the iteration process, making better use of InSAR data and calculating relatively stable 3D displacement fields.
In the simulation examples, uniform and non-uniform schemes are adopted to select GPS points. In the uniform sampling scheme, the spatial resolutions of GPS points of 30 km, 20 km, 14 km, and 10 km were considered. Accompanied with the comprehensive analyses of the reverse fault, normal fault with a strike-slip component, and strike-slip fault with a reverse component, it demonstrates that the VOILS method is better than the direct solution method in both horizontal and vertical directions. From the results of the uniform sampling scheme, we see that the percentages of improvement for the reverse fault are approximately 3~9% and 70%, for the normal fault with a strike-slip component approximately 4~8% and 68%, and for the strike-slip fault with a reverse component approximately 1~8% and 22%. Moreover, the degree of improvement in the vertical direction is more obvious. The simulation experiments show that the VOILS method can make better use of the InSAR observations, and has feasibility and validity.
The VOILS method was applied to get the 3D displacement field of the 2008 Wenchuan earthquake. The LOS residual magnitude obtained by the VOILS method is mm-level, while the magnitude of the direct solution method is cm-level. Furthermore, the results of the VOILS method are relatively stable. Based on the analyses of 3D deformation, it can be seen that the mechanism of the main shock was mainly a reverse motion with a right-lateral strike-slip component. In the northeast of the fault, the mechanism is dominated by the right-lateral strike-slip. In particular, the comparison with the rupture displacements from the field investigations shows that the VOILS method is better than the direct solution method.