- Email: [email protected]
Deriving 3D coseismic deformation field by combining GPS and InSAR data based on the elastic dislocation model
International Journal of Applied Earth Observation and Geoinformation 57 (2017) 104–112
Contents lists available at ScienceDirect
International Journal of Applied Earth Observation and Geoinformation journal homepage: www.elsevier.com/locate/jag
Research Paper
Deriving 3D coseismic deformation field by combining GPS and InSAR data based on the elastic dislocation model Xiaogang Song a , Yu Jiang a,b,∗ , Xinjian Shan a , Chunyan Qu a a b
State Key Laboratory of Earthquake Dynamics, Institute of Geology, CEA, Beijing 100029, China China University of Petroleum (East China), Qingdao, 266580 China
a r t i c l e
i n f o
Article history: Received 26 May 2016 Received in revised form 21 November 2016 Accepted 23 December 2016 Keywords: 3-D coseismic deformation Elastic dislocation model GPS InSAR
a b s t r a c t The density of GPS measurements is usually one of the key issues in resolving 3-D coseismic deformation field from integrating GPS and Interferometric Synthetic Aperture Radar (InSAR) measurements with pure mathematic interpolation methods An approach that combines the elastic dislocation model with the Best Quadratic Unbiased Estimator (BQUE) or a robust estimation method named IGG (Institute of Geodesy and Geophysics) is proposed to reconstruct 3-D coseismic deformation field, in which only a small amount of GPS data is needed to produce a reasonable initial 3-D coseismic deformation. Then the BQUE and IGG are used to weight the InSAR and GPS measurements to avoid computational issues caused by the negative variance problem and to decrease the impact from gross errors. The Wenchuan earthquake is used to test the proposed method. We find that the developed method makes it possible to use only a few GPSs and InSAR data to recover the 3-D coseismic deformation field, which offers extensive future usage for measuring earthquake deformation, particularly in some tectonically active regions with sparse GPS measurements. © 2016 Elsevier B.V. All rights reserved.
1. Introduction The past decades have witnessed a tremendous improvement in InSAR’s ability to measure deformation of the earth’s surface. InSAR measurements provide us a wide, spatially continuous coverage to well document some geophysical phenomena such as earthquakes (Funning et al., 2005; Li et al., 2011; Elliott et al., 2013), subsidence (Dehghani et al., 2009), landslides (Catani et al., 2005), and volcanoes (Hooper et al., 2007). However, the limitation lies in that it is sensitive to surface movements only in line-of-sight direction (Wright et al., 2004), so the results cannot fully reflect the real deformation in most cases. In order to acquire the 3-D deformation, GPS and InSAR are commonly integrated, which was proposed by Ge et al. (2000) in the Double Interpolation and Double Prediction (DIDP) style. Gudmundsson et al. (2002) produced highresolution 3-D surface motion maps for Iceland by combining GPS and InSAR based on Markov random field regularization and simulated annealing. Samsonov and Tiampo (2006) and Samsonov et al. (2007), developed an analytical optimization strategy and applied
∗ Corresponding author at: China University of Petroleum (East China), Qingdao, 266580, China. E-mail addresses: [email protected] (X. Song), [email protected] (Y. Jiang), [email protected] (X. Shan), [email protected] (C. Qu). http://dx.doi.org/10.1016/j.jag.2016.12.019 0303-2434/© 2016 Elsevier B.V. All rights reserved.
it to reconstruct the 3-D deformation field of Southern California from InSAR and GPS measurements. Taking the interpolated GPS observations on the InSAR grid or a regular grid as a first guess map, these methods are strongly dependent on the availability and distribution of the GPS measurements. Due to exploitation of the densest continuous GPS networks in the world (327 GEONET sites are used), Hu et al. (2013) adopted the method of weighted least squares (WLS) and successfully retrieved the 3-D coseismic surface displacement field of the 2011 Tohoku-Oki earthquake on the onshore non-near-fault area by appropriately integrating the InSAR and GPS measurements. In the aforementioned methods, GPS data are generally interpolated by some pure mathematic interpolating algorithms without taking the mechanical property of the media into account; these methods are suitable for deforming areas with low-gradient and a sparse GPS network, but cannot deal with sudden and large ground movements caused by large earthquakes. A few researchers considered the physics of the deformation in interpolation. Based on the elastic theory, Guglielmino et al. (2011) proposed a SISTEM approach to obtain the 3-D displacement map of an active stratovolcano. However, this method is developed to study slow and continuous movement under the hypothesis of infinitesimal and homogeneous strain, such as volcanic activities and interseismic deformations, rather than coseismic deformations.
X. Song et al. / International Journal of Applied Earth Observation and Geoinformation 57 (2017) 104–112
Another key issue involved in merging heterogeneous observations is how to allocate the weight for different types of observations. The common practice is to weight the observations based on their a priori variance, but it is always difficult to know the exact variance of the observations. Hu et al. (2012) introduced the variance component estimation (VCE) approach to iteratively estimate the variances of the observations, and obtained the posterior estimation of the stochastic model of the InSAR measurements in the 3-D displacement velocity estimation from InSAR and GPS. However, the VCE approach expects enough redundant measurements and is easily end up with the negative variance problem. In addition, methods based on the least squares (LS) adjustment would also give biased estimation when gross errors exist in the observation samples. Large (moment magnitude (Mw) > 7) continental earthquakes in China often occur on the western fault systems, where the GPS density is very sparse (>50 km) although two GPS station construction projects (the Crustal Movement Observation Network of China: CMONOC-I, and the Tectonic and Environmental Observation Network of Mainland China: CMONOC-II) have been conducted in the past few decades. GPS observations are usually not dense enough to represent both the large and small scale coseismic deformation signal, which leads to significant problems in GPS interpolation (the first step in deriving 3-D surface displacement fields by integrating GPS and InSAR), especially for the deforming area with large spatial gradient in the near field of the causative fault. In this paper, we develop the fault dislocation model to invert a small number of GPS observations to recover a reasonable and reliable initial 3-D coseismic deformation field i. The Best Quadratic Unbiased Estimator (BQUE) (Sjöberg, 1983) is first used to estimate the variance components to avoid breaking the computation process due to the negative VC problem in VCE. Then GPS and InSAR are merged by a weighted least squares adjustment. Considering the contamination of atmospheric delays and other factors in the InSAR measurements, its stochastic model may not fit an exact normal distribution, so the gross error cannot be avoided. Therefore, a robust estimation method named IGG proposed by Zhou (1989) is applied to reduce the effect from gross errors. The 2008 Wenchuan earthquake is used as a case study to validate these two methods. 2. Methods Assuming a deformation field caused by an earthquake can be discretized into a m × n grid. InSAR and GPS observation is downsampled or interpolated on the grid at first. For each node of the grid, the LOS displacement observation dlos can be related to the real 3-D deformation components [ Dew model: dlos + Vlos = [ cos ˛ sin
− sin ˛ sin
T
Dud ] by the following
Dsn
− cos ] · [ Dew
Dsn
Dud ]
T
(1)
where ˛ is the azimuth angle of InSAR, and is the incidence angle. Vlos is the residuals of InSAR observation. The interpolated 3-D GPS observation [ xGPS yGPS tion components by: [ xGPS
yGPS
T
zGPS ] can be related to the 3-D deforma-
T
zGPS ] + VGPS = [ Dew
Dsn
Dud ]
T
(2)
Then, the observation equation can be formed (Eq. (3)), and to estimate the final 3-D displacement field through the least-squares (LS). L4mn×1 = D4mn×3mn X3mn×1 + V4mn×1 T
(3)
Where X = [Dew Dsn Dud ] is the unknown vector of displaceT ment on each node of the grid; L = [xgps ygps zgps dlos ] is the vector of GPS and InSAR observations; V is the vector of observation residuals; and D is the design matrix. Based on the Eq. (3), the
105
interpolated GPS and InSAR observations can then be used in a WLS model to retrieve the optimal 3-D displacement field. The key steps are (i) the interpolation of sparse GPS data appropriately, and (ii) the determination of the weighting scheme (estimating the posterior variances) for the observations. Here, we use sparse GPS data and fault dislocation model to produce a reasonable initial 3-D coseismic deformation (see Section 2.1), which will be further corrected by fusing InSAR observations. The weights for the observations are determined by using BQUE (see Section 2.2) and IGG (see Section 2.3) respectively. We will also compare the results from these two methods, and assess their accuracy by using field measurements. The overall processing strategy is also demonstrated in Fig. 1. 2.1. Recovering initial 3-D coseismic deformation field combining GPS with fault dislocation model To derive 3-D surface displacement fields by integrating GPS and InSAR, the first step is to resample the sparse GPS measurements onto the InSAR grid or another regular grid. Kriging interpolation is often used because it can preserve the theoretical spatial correlation of the measurements of adjacent data points. However, the distribution and density of the known points are always the key to maintain a good interpolating accuracy for general mathematic methods, e.g. the Kriging, bicubic spline interpolation and least squares collocation. Such methods are limited by the sparse GPS network in western China, since most faults in western China are located in remote areas where installing GPS stations is difficult and expensive. As a result, pure mathematic interpolation methods are not capable of reproducing a reliable deformation field caused by a strong earthquake. In this study, we seek for solutions that combines the physics of earthquake deformation and mathematical methods. The elastic dislocation model derived by Okada (1985, 1992) provides a good representation for surface and internal deformation generated by shear and tensile faults in a homogenous half-space or a multilayered half-space (Wang et al., 2003). It is usually used to invert the surface observations for the focal parameters and slip distribution at certain depth (Delouis et al., 2002; Feigl et al., 2002; Zhang et al., 2011). These earthquake parameters can provide us with very useful knowledge of the geophysical properties of an earthquake. Based on the inverted fault model, we can simulate the coseismic deformation in any direction. However, the deformation field from forward modeling cannot fully represent the surface deformation, regardless of the observations used in the inversion. This is because of the assumption and simplification on fault model (presumptions about the kinematics (e.g. slip and rake) and geometry (e.g. dip and strike)), inverting algorithms that we adopt, and spatial density of geodetic observations with uncertainties. A fine and accurate 3-D coseismic deformation field can not only be used to constrain earthquake source parameters, but also allows for investigations of other geophysical and geological problems such as shallow slip deficit and fault segmentation. Practically, coseismic ruptures and deformation are also important for post-disaster rescue. Observations of the earth’s surface are the most direct and effective tool to derive the surface information. In the case of insufficient observational data, a geophysical model considering the mechanical property of the media can be used to assist the reconstruction of the 3-D surface deformation conversely. This is also a new attempt applied in other fields (3-D surface displacement caused by Geothermal field (Hu et al., 2016), source fluid volumetric change of groundwater depletion (Castellazzi et al., 2016). The dislocation model describes the physical relationship between the displacement on faults and the deformation on the surface, regardless of the distribution and density of the surface geodetic observations. Here, we use the dislocation model to generate the first approximation of a 3-D coseismic deformation field,
106
X. Song et al. / International Journal of Applied Earth Observation and Geoinformation 57 (2017) 104–112
Fig. 1. The flowchart of integrating InSAR and GPS to derive the 3-D coseismic deformation field.
Fig. 2. Flow chart showing the procedure of using the dislocation model and a few GPSs to generate the first approximation of a 3-D coseismic deformation field.
with a small amount of GPS observations. The procedure is outlined in below and also in Fig. 2. Step 1 constructing the initial fault model: Fault segments can be mapped based on visual analysis of high-resolution imagery, and displacement field from InSAR or image correlation. Fault dip of each segment can be obtained from seismic data. The location, length and strike of each fault segment are adjusted iteratively in step 2. Step 2 inverting for slip distribution: Based on the initial fault model constructed in step 1, sparsely distributed GPS data are used to estimate the slip distribution on the fault plane, as well as the location, length and strike of each fault segment simultaneously. This is carried out iteratively using the SDM modeling program (Wang et al., 2009). The effect of GPS density on the inverted results will be tested in Section 3, aiming to find a stable solution with as less GPSs as possible. Step 3 synthesizing 3-D surface deformation field: based on the slip distribution and fault parameters from step 2, we could synthesize a 3-D coseismic deformation field on a regular grid by forward modeling.
vations is sufficiently high. The result also depends on the level of realistic of a priori adjustment and stochastic models, and the quality of the initial VCs (Eshagh and Sjoberg, 2008). The real cases usually do not satisfy all of these requirements. Furthermore, the strict Helmert VC estimation expects independent observational groups and enough redundant measurements, and is easy to suffer from the negative variance problem. As for coseismic deformation, InSAR measurements are usually contaminated by the atmospheric and orbital errors more or less, and the number of redundant observation is only in LOS direction, relative to three unknowns (two horizontal components and one vertical component) for each point. The Helmert VCE generally cannot well handle this situation. Thus, in this case we propose to use the Best Quadratic Unbiased Estimator from the functional model of conditional adjustment with unknown parameters (Eshagh and Sjoberg, 2008) to estimate the VC of InSAR and GPS iteratively. Given the following Gauss-Helmert model,
2.2. The BQUE of VC for InSAR and GPS
where A and B are the first and second design matrices, X is the unknown vector, ε is the residual vector, W is the misclosures vector, Q is the covariance matrix of observations, and E {• } is the statistical expectation.
Most of the VC estimation methods work well in the case where no bias exists in the observations and the redundancy of the obser-
AX + Bε = W,
with
εεT
=Q
and
E {ε} = 0
(4)
X. Song et al. / International Journal of Applied Earth Observation and Geoinformation 57 (2017) 104–112
Eq. (4) can be converted to the conditional model by adjustment multiplying both sides of the equations by I − A0 , namely
I − A0 Bε = I − A0 W,
¯ =W ¯ Bε
or
(5)
where A0 = AA− , A− stands for a specific generalized inverse of A. For the conditional adjustment (Eq. (5)), the best quadratic unbiased estimation of variance components can be written as (Eshagh and Sjoberg, 2008): ˆ = S − uT
Sij = trace RCi RCj , R=C C=
q
0
I−A
i2 Ci
0 T
C −1
Ci = BQi BT
i=1
i = 1, 2, . . ., q.
In this passage, we construct one condition adjustment equation by 4 observations, dlos , xGPS , yGPS , zGPS , on each grid point:
⎡
−1 cos ˛ sin
− sin ˛ sin
− cos
dlos
⎢ ⎢ xGPS ⎢ ⎢y ⎣ GPS
⎤
⎥ ⎥ ⎥=0 ⎥ ⎦
|v| < a
1 p (v) = ⎪ ⎩ |v| + k 0
(9)
a < |v| < b |v| ≥ b
where a and b are constant values, and depend on the observations. Empirically, we use a = 1.5, b = 2.5; k= 0.001 in the case of v = 0.
3.1. Data
(8)
YGPS 1
···
Hence, the observation equation for N points is, =0
The 12 May 2008 Mw 7.9 Wenchuan earthquake is one of the most destructive events that had ever occurred in China during the last 50 years, killed about 70,000 people (Zhang et al., 2008). This event ruptured the NE-SW striking, north-west dipping Longmenshan fault belt, which was formed by continuing eastward extrusion of the Tibet plateau to the western edge of the Sichuan Basin (Zhang et al., 2008). Many geological and geodetic studies have focused on this event and have provided abundant material to test the effectiveness of our method.
The JAXA arranged a mission of emergency and obtained a batch of L-band SAR data (ALOS/PALSAR) soon after the Wenchuan earthquake. Tracks 471–477 covered the whole surface rupture zone (Fig. 3(a)). The data were processed using GAMMA (Qu et al., 2008). Good coherence is preserved for most of the area, except for the near-fault zone where the decorrelation is caused by large surface deformation. The near-fault displacement field was therefore achieved by using the offset-tracking technique, based on which
(7)
zGPS
B• L
⎧ ⎪ ⎨1
3. Experiments
i, j = 1, 2, . . ., q,
= I−A
and
¯ T RCi RW ¯ , ui = W
absolute value of residual |v|, depending on the quadratic mean deviation (mean square error) :
(6)
where S − is a generalized inverse of S and
−1
107
where B is the design matrix
and L is the observation vector
L=
dlos 1
dlos 2
···
dlos N
XGPS 1
XGPS 1
···
XGPS N
YGPS 1
Generally,W = −B• L = / 0, because of various errors, including orbital and atmospheric errors in GPS and InSAR. We develop an iterative procedure initialing with a priori estimates of the variance components, to caulculate the variance
of the unit weight 0 2 = W T BQBT
−1
W and variance compo-
nents using Eq. (6), until0,i 2 = 0,i+1 2 . In the end, we obtain the estimation of 1 2 , 2 2 , 3 2 , 4 2 for the variance of observation dlos , xGPS , yGPS , zGPS . 2.3. The robust estimation by IGG Due to the simplicity and general good performance, the least squares (LS) adjustment has been extensively applied in geodetic data processing under the hypothesis that observation samples are subject to normal distribution. Nevertheless, gross errors often exist in observations, which would bias the result from the LS adjustment with a bigger variance. Robust estimation aims at decreasing the effect of gross errors. Krarup et al. (1980) introduced the robust estimation theory to surveying by proposing the Denmark method, based on the M estimation theory. Zhou (1989) put forward the concept of equivalent weights, made LS method formalization for M estimation by equivalent weights, and named the algorithm as IGG. In the LS adjustment (Eq. (3)), we can determine the weight p (v) of the
YGPS N
ZGPS 1
···
ZGPS 1
ZGPS N
T
we can map the full picture of the rupture trace. GPS data are from Wang et al. (2011), which comes from various sources, including CMONOC, the national key basic research and development project, and the Earthquake Administrations of Sichuan and Chongqing. Fig. 3(b) shows the location of GPS data used in the rest of the computation and accuracy control in this paper. 3.2. Calibration of InSAR data using GPS Depending on the acquisition time of the slave image, the interferometric phase contains the post-seismic deformation phase to some extend even in the far field, as well as the atmospheric effect, both biased the baseline refinement in InSAR data processing. Furthermore, unlike referenced GPS data, InSAR only gives us a relative measurement for the surface deformation. So it is necessary to make certain corrections to the spatial reference of InSAR measurements before integrating InSAR and GPS. • the first step is to project the displacement vector of each GPS station: u = Llos = s · u
ue
un
uu
T
into the InSAR LOS direction: (10)
108
X. Song et al. / International Journal of Applied Earth Observation and Geoinformation 57 (2017) 104–112
Fig. 3. (a) The study area: black vectors show the GPS velocities; black rectangles indicate the coverage of 7 PALSAR stripes; red lines show the surface rupture of the Wenchuan earthquake. (b) Blue rectangles show the fault geometry we used in our inversion. Circles stand for the location of the GPS points (the color indicates the serial number of GPS we used in calculating the deformation field when we change the number of GPS input). Black stars indicate the GPSs for checking. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
where the unit vector s = cos ˛ sin − sin ˛ sin − cos , ˛ is the azimuth angle of InSAR, and is the incidence angle. • Use a quadratic polynomial to construct the relationship between the InSAR and GPS displacements in the LOS direction, and then estimate the coefficients for each track, n = am2 + bm + c
plan
Interpolating method
Weighting method
I II III IV
Dislocation model Dislocation model Kriging Kriging
IGG BQUE IGG BQUE
(11)
Where a, b, c are the coefficients, and m, n are the LOS InSAR and GPS displacements in the LOS direction at each GPS station. • Make corrections to InSAR observations: y = ax2 + bx + c
Table 1 Combinations of different interpolating and weighting methods.
(12)
where x is the original InSAR observation, and y is the corrected value. • Mosaic 7 tracks to obtain a complete InSAR deformation field.
the focal mechanism solutions of the main shock and aftershocks (Zhang et al., 2011). • GPS interpolation. Following the procedure described in Section 2.1, GPS data is used to synthesize the 3-D surface deformation field of the Wenchuan earthquake on the regular grid by forward modeling. • Integrating InSAR and GPS. The WLS adjustment is used to integrate InSAR and GPS on each wall, and the BQUE and IGG applied to weight the observations. The results on the hanging wall and footwall are mosaicked to generate a single complete displacement field. 10 selected GPS points are used to evaluate the accuracy of the results.
3.3. Processing flowchart 4. Results and discussion The data processing scheme of the Wenchuan earthquake is similar to that in Fig. 1. The only difference is that we need to carry out data processing for hanging wall and footwall separately, and then mosaic their results in the final step. The detail description is as follow: • Downsample the InSAR data. Within the common coverage of the InSAR and GPS data, we confine our study area to the range of 30.5-32.5N and 103.2-105.2E, and grid it with a spacing of 0.2 ◦ in latitude and longitude. The corrected InSAR deformation field is downsampled onto this regular grid using bilinear interpolation. • Construct the fault model for the Wenchuan earthquake. The displacement field from offset-tracking shows the rupture trace of the Wenchuan earthquake clearly, based on which we estimate the fault geometry parameters such as the location, length and strike of each segment. The dip and rake of each segment are from
Here we combine different interpolating and weighting methods for comparison (Table 1). Following the procedure described in previous sections, we calculate the corresponding 3-D coseismic deformation fields and associate RMSEs from 10 check points. Fig. 4 shows the 3-D deformation fields from 4 combinations by using 30 (a) and 80 (b) GPS points as input. Note that the selected GPS points in each inversion and accuracy assessment should be distributed evenly over the whole deforming area, to minimize the impact of GPS distribution on the results. As shown in Fig. 4, we can see the deformation fields from plans I and II present more details, and less dependence on the number of GPSs compared with those from plans III and IV. This is because Kriging requires an appropriate theoretical semi-variogram model for each component to be interpolated. It is one of the main critical problems in geostatistics. Given the large area of the Wenchuan rupture, tens of GPS points
X. Song et al. / International Journal of Applied Earth Observation and Geoinformation 57 (2017) 104–112
are too few to describe the whole deformation area using Kriging. Instead, plans I and II, which use the fault dislocation model that accounts for the medium property and deformation mechanism, show reasonable and stable results. To assess the accuracy of our methods quantitatively, we used field observations along the rupture from Xu et al. (2008). The field measurements of horizontal and vertical fault displacements on a few points are converted into the EW, SN, and vertical components assuming NE45 striking. Figs. 5 S2 and S3 shows the comparison between the field measurements and our model simulations. The RMSEs of the Plan I (dislocation model + IGG) and II (dislocation model +BQUE) with 80 GPS points are 47.8 cm, 18.3 cm, 21.5 cm
109
and 64.8 cm, 49.8 cm, 73.7 cm for EW, SN and vertical direction respectively. Clearly, both methods improved the accuracy compared to the “Kriging + IGG” result with the accuracy of 113.9 cm, 65.3 and 202.5 cm. It indicates the near-field displacements derived from the method that integrates the dislocation model agree better with the field measurements than those from the method using Kriging interpolation. This also implies that the interpolation based on the elastic dislocation model is better than Kriging. Considering the maximum on-fault displacement is over 5m, the result with an error of less than 50 cm (Fig. 5) is satisfying. In terms of the weighting methods, the IGG works better than the BQUE based on the RMSEs. The difference lies in the fact that
Fig. 4. The 3-D coseismic deformation fields of the Wenchuan earthquake from 4 plans listed in Table 1, with (a) 30 and (b) 80 GPS points as input.
110
X. Song et al. / International Journal of Applied Earth Observation and Geoinformation 57 (2017) 104–112
Fig. 4. (Continued)
the BQUE assumes that errors are subject to the normal distribution whereas the IGG does not; as the result, the BQUE is easier to be distorted by outliers and thus requires careful elimination of outliers. The IGG shows a strong ability to decrease the weight of outliers, which is an advantage in cases where observations contain gross errors like in the Wenchuan earthquake. Our experiments demonstrate the advantages of using the fault dislocation model in interpolation that only a small amount of GPS observations is needed. It suggests that our method is very useful in many tectonically active regions where GPS stations are sparsely distributed. In order to test the sensitivity of the interpolation method to the number of GPS stations, we also vary the number
of GPS observations from 10 to 80 in steps of 10 in the slip distribution inversion. The RMSEs of the three components of surface displacement are calculated using 10 check points (the distribution is shown in Fig. 2(b)). Fig. 6 shows the corresponding RMSE in the EW, SN, and vertical directions at different number of GPS points. The RMSEs from plans I and II remain stable when using more than 20 GPS points. While the RMSEs from plans III and IV would be equal to or slightly smaller than those from plans I and II, only when the number of GPS input reaches 80. In case of the number of GPS check points increase to 30, the RMSE curves (Fig. S1) show similar statistical results. The difference is that the RMSEs from plans III and IV tend to be more stable than those from plans I
X. Song et al. / International Journal of Applied Earth Observation and Geoinformation 57 (2017) 104–112
111
Fig. 5. Field and modeled measurements of EW, SN and vertical fault displacements. Black triangles: field measurements vs. Plan I (IGG), and its RMSEs are 47.8 cm, 18.3 and 21.5 cm for EW, SN and vertical direction respectively. Green circles: field measurements vs. the Plan III (Kriging) with 80 GPS points, and its RMSEs are 113.9 cm, 65.3 and 202.5 cm for EW, SN and vertical direction respectively. The red line denotes the line ‘y=x’. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
Fig. 6. The RMSEs calculated by using 10 GPS points in (a) EW, (b) SN, (c) vertical directions with an increasing number of the input GPS points.
and II when the number of GPS input reaches 50. This indicates that we can decrease the number of GPS points significantly by incorporating the fault dislocation model and still obtain a good result, similar to that by using 80 GPS points and Kriging method. for a large deforming area (500 × 500 km) such as the Wenchuan earthquake, 20 evenly distributed GPS points are sufficient enough. This provides us information to build economical and applicable plan for those tectonically active regions with sparse GPS network. We also notice that the RMSEs (Fig. 6) calculated by using 10 GPS check points from plans III and IV have increased first and then
began to drop when the number of GPS points increased from 10 to 30. We think that it is because a small number of GPS points (e.g., ≤ 10) are not enough to present the characteristics of the widerange deformation field of the Wenchuan earthquake when using Kriging. In such case the distance between the input points and the check points plays a dominant role in calculating the RMSEs, thus RMSEs lost its statistical meaning. Fig. S1 did not show such fluctuation in the beginning, because a larger number of GPS check points (number of 30) are used to calculate the RMSEs, which reduce the dependence on the distance between the input points and the check
112
X. Song et al. / International Journal of Applied Earth Observation and Geoinformation 57 (2017) 104–112
points. In addition, there is no doubt that the accuracy for the Kriging will decrease when more GPS points are chosen as check points, which means less points are chosen for interpolation. This is the reason that the RMSEs from plans III and IV are always larger than that from plans I and II in Fig. S1. 5. Conclusions The derivation of 3-D coseismic deformation field is very important for revealing detail movement of the fault and the area affected by an earthquake. GPS observations, even for a dense GPS network, are usually not enough to recover the real 3-D coseismic deformation field (especially in the near-fault area) with large spatial coverage like that in Wenchuan earthquake example. An economical and applicable approach has been developed in this paper that can retrieval 3-D coseismic deformation from InSAR and small number of GPS data by combining the elastic dislocation model with the BQUE or the robust estimation method IGG. Compared with the pure mathematic interpolator discussed in this study, the elastic dislocation model that considers the mechanical property of the media allows us to estimate initial 3-D coseismic deformation in a more physically reasonable and accurate style. By using the BQUE and IGG, the weights of the InSAR and GPS observations have been determined optimally, in the case where large atmospheric and orbital errors exist in the InSAR observation and the redundancy of observations is not sufficiently high. An important advantage of presented approach is that the requirement for the density of GPS is significantly lowered to reconstruct a reliable 3-D coseismic deformation field. It is very useful for remotely located tectonically active regions with a sparse GPS network. The proposed method has been validated by using a series of experiments with the test case of Wenchuan earthquake. The results have clearly shown that our method is sound and can significantly improve the accuracy of the 3-D coseismic deformation reconstruction. As our method is generic, it is also possible to integrate other sources of surface movement measurements including Multiple-Aperture InSAR (MAI), offset-tracking and leveling. Overall, the method offers extensive future usage for measuring earthquake deformation. Acknowledgments This work was supported jointly by the Basic Scientific Funding of Institute of Geology, China Earthquake Administration (No.IGCEA1404), grants of State Key Laboratory of Earthquake Dynamics, China Earthquake Administration (No.LED2013A02), and the National Natural Science Foundation of China (No. 41461164002). All ALOS PALSAR data were provided and are copyrighted by JAXA. Most figures were made using the public domain Generic Mapping Tools. We thank Prof. Mehdi Eshagh in University West for instructions of the usage of the BQUE. References Castellazzi, P., Martel, R., Galloway, D.L., Longuevergne, L., Rivera, A., 2016. Assessing groundwater depletion and dynamics using GRACE and InSAR: potential and limitations. Ground Water, http://dx.doi.org/10.1111/gwat. 12453. Catani, F., Farina, P., Moretti, S., Nico, G., Strozzi, T., 2005. On the application of SAR interferometry to geomorphological studies: estimation of landform attributes and mass movements. Geomorphology 66, 119–131. Dehghani, M., Zoej, M.J.V., Entezam, I., Mansourian, A., Saatchi, S., 2009. InSAR monitoring of progressive land subsidence in Neyshabour, northeast Iran. Geophys. J. Int. 178, 47–56. Delouis, B., Giardini, D., Lundgren, P., Salichon, J., 2002. Joint inversion of InSAR, GPS, teleseismic, and strong-motion data for the spatial and temporal distribution of earthquake slip: application to the 1999 Izmit mainshock. Bull. Seismol. Soc. Am. 92, 278–299, http://dx.doi.org/10.1785/0120000806.
Elliott, J.R., Copley, A.C., Holley, R., Scharer, K., Parsons, B., 2013. The 2011 Mw 7.1 Van (Eastern Turkey) earthquake. J. Geophys. Res. 118, 1619–1637, http://dx. doi.org/10.1029/2012JB009569. Eshagh, M., Sjoberg, L., 2008. The modified best quadratic unbiased non-negative estimator (MBQUNE) of variance components. Stud. Geophys. Geod. 52, 305–320. Feigl, K., Sarti, F., Vadon, H., 2002. Estimating slip distribution for the I˙ zmit mainshock from coseismic GPS, ERS-1, RADARSAT, and SPOT measurements. Bull. Seismol. Soc. Am. 92, 138–160. Funning, G.J., Parsons, B., Wright, T.J., Jackson, J.A., Fielding, E.J., 2005. Surface displacements and source parameters of the 2003 Bam (Iran) earthquake from Envisat advanced synthetic aperture radar imagery. J. Geophys. Res. 110, B09406, http://dx.doi.org/10.1029/2004JB003338. Ge, L., Han, S., Rizos, C., 2000. The double interpretation and double prediction (DIDP) approach for InSAR and GPS integration. In: ISPRS Commission Reports, the XIXth Congress of the International Society for Photogrammetry and Remote Sensing, Amsterdam, Netherlands, pp. 205–212. Gudmundsson, S., Sigmundsson, F., Carstensen, J., 2002. Three-dimensional surface motion maps estimated from combined interferometric synthetic aperture radar and GPS data. J. Geophys. Res. 107, http://dx.doi.org/10.1029/ 2001jb000283. Guglielmino, F., Nunnari, G., Puglisi, G., 2011. Simultaneous and integrated strain tensor estimation from geodetic and satellite deformation measurements to obtain three-dimensional displacement maps. IEEE Trans. Geosci. Remote Sens. 49, 1815–1826. Hooper, A., Segall, P., Zebker, H., 2007. Persistent scatterer interferometric synthetic aperture radar for crustal displacement analysis, with application to Volcano Alcedo, Galapagos. J. Geophys. Res.-Solid Earth 112, B07407. Hu, J., Li, Z.W., Sun, Q., Zhu, J.J., Ding, X.L., 2012. Three-dimensional surface displacements from InSAR and GPS measurements with variance component estimation. IEEE Geosci. Remote Sens. Lett. 9, 754–758. Hu, J., Li, Z.W., Ding, X.L., Zhu, J.J., Sun, Q., 2013. Derivation of 3-D coseismic surface displacement fields for the 2011 Mw 9.0 Tohoku-Oki earthquake from InSAR and GPS measurements. Geophys. J. Int. 192, 573–585, http://dx.doi.org/10. 1093/gji/ggs033. Hu, J., Wang, Q., Li, Z., Zhao, R., Sun, Q., 2016. Investigating the ground deformation and source model of the yangbajing geothermal field in tibet, China with the WLS InSAR technique. Remote Sens. 8, 191, http://dx.doi.org/10.3390/ rs8030191. Krarup, T., Juhl, J., Kubik, K., 1980. Gotterdammerung over least squares adjustment. In: Proceeding of 14th Congress of the International Society of Photogrammetry, Hamburg, Germany, B 3, pp. 369–378. Li, Z., Elliott, J.R., Feng, W., Jackson, J.A., Parsons, B., Walters, R.J., 2011. The 2010 Mw 6.8 yushu (Qinghai, China) earthquake: constraints provided by InSAR and body wave seismology. J. Geophys. Res. 116, B10302, http://dx.doi.org/10. 1029/2011JB008358. Okada, Y., 1985. Surface deformation due to shear and tensile faults in a half-space. Bull. Seismol. Soc. Am. 75, 1135–1154. Okada, Y., 1992. Internal deformation due to shear and tensile faults in a half-space. Bull. Seismol. Soc. Am. 82, 1018–1040. Qu, C., Song, X., Zhang, G., 2008. Analysis on the characteristics of InSAR coseismic deformation of the Ms8.0 Wenchuan earthquake. Seismol. Geol. 30, 1076–1083. Samsonov, S., Tiampo, K., 2006. Analytical optimization of a DInSAR and GPS dataset for derivation of three-dimensional surface motion. IEEE Geosci. Remote Sens. Lett. 3, 107–111. Samsonov, S., Tiampo, K., Rundle, J., Li, Z.H., 2007. Application of DInSAR-GPS optimization for derivation of fine-scale surface motion maps of Southern California. IEEE Trans. Geosci. Remote Sens. 45, 512–521. Sjöberg, L.E., 1983. Unbiased estimation of variance-components in condition adjustment with unknowns-a MINQUE approach. Zeitschrift für Vermessungen 108, 382–387. Wang, R., Martin, F.L., Roth, F., 2003. Computation of deformation induced by earthquakes in a multi-layered elastic crust: FORTRAN programs EDGRN/EDCMP. Comp. Geosci. 29, 195–207. Wang, L., Wang, R., Roth, F., Enescu, B., Hainzl, S., Ergintav, S., 2009. Afterslip and viscoelastic relaxation following the 1999 M 7.4 Izmit earthquake from GPS measurements. Geophys. J. Int. 178, 1220–1237. Wang, Q., Qiao, X., Lan, Q., Freymueller, J., Yang, S., Xu, C., Yang, Y., You, X., Tan, K., Chen, G., 2011. Rupture of deep faults in the 2008 Wenchuan earthquake and uplift of the Longmen Shan. Nat. Geosci. 4, 634–640. Wright, T.J., Parsons, B.E., Lu, Z., 2004. Toward mapping surface deformation in three dimensions using InSAR. Geophys. Res. Lett. 31. Xu, X.W., Wen, X., Ye, J.Q., 2008. The Ms 8.0 Wenchuan earthquake surface ruptures and its seismogenic structure. Seismol. Geol. 30, 597–629. Zhang, P.Z., Xu, X.W., Wen, X.Z., Ran, Y.K., 2008. Slip rates and recurrence intervals of the Longmen Shan active fault zone and tectonic implications for the mechanism of the May 12 Wenchuan earthquake, 2008, Sichuan. China. Chin. J. Geophys. 51, 1066–1073. Zhang, G., Qu, C., Shan, X., Song, X., Zhang, G., Wang, C., Hu, J., Wang, R., 2011. Slip distribution of the 2008 Wenchuan Ms 7.9 earthquake by joint inversion from GPS and InSAR measurements: a resolution test study. Geophys. J. Int. 186, 207–220. Zhou, J.W., 1989. Classical theory of errors and robust estimation. Acta Geodetica et Cartographica Sinica 18, 115–120.
Contents lists available at ScienceDirect
International Journal of Applied Earth Observation and Geoinformation journal homepage: www.elsevier.com/locate/jag
Research Paper
Deriving 3D coseismic deformation field by combining GPS and InSAR data based on the elastic dislocation model Xiaogang Song a , Yu Jiang a,b,∗ , Xinjian Shan a , Chunyan Qu a a b
State Key Laboratory of Earthquake Dynamics, Institute of Geology, CEA, Beijing 100029, China China University of Petroleum (East China), Qingdao, 266580 China
a r t i c l e
i n f o
Article history: Received 26 May 2016 Received in revised form 21 November 2016 Accepted 23 December 2016 Keywords: 3-D coseismic deformation Elastic dislocation model GPS InSAR
a b s t r a c t The density of GPS measurements is usually one of the key issues in resolving 3-D coseismic deformation field from integrating GPS and Interferometric Synthetic Aperture Radar (InSAR) measurements with pure mathematic interpolation methods An approach that combines the elastic dislocation model with the Best Quadratic Unbiased Estimator (BQUE) or a robust estimation method named IGG (Institute of Geodesy and Geophysics) is proposed to reconstruct 3-D coseismic deformation field, in which only a small amount of GPS data is needed to produce a reasonable initial 3-D coseismic deformation. Then the BQUE and IGG are used to weight the InSAR and GPS measurements to avoid computational issues caused by the negative variance problem and to decrease the impact from gross errors. The Wenchuan earthquake is used to test the proposed method. We find that the developed method makes it possible to use only a few GPSs and InSAR data to recover the 3-D coseismic deformation field, which offers extensive future usage for measuring earthquake deformation, particularly in some tectonically active regions with sparse GPS measurements. © 2016 Elsevier B.V. All rights reserved.
1. Introduction The past decades have witnessed a tremendous improvement in InSAR’s ability to measure deformation of the earth’s surface. InSAR measurements provide us a wide, spatially continuous coverage to well document some geophysical phenomena such as earthquakes (Funning et al., 2005; Li et al., 2011; Elliott et al., 2013), subsidence (Dehghani et al., 2009), landslides (Catani et al., 2005), and volcanoes (Hooper et al., 2007). However, the limitation lies in that it is sensitive to surface movements only in line-of-sight direction (Wright et al., 2004), so the results cannot fully reflect the real deformation in most cases. In order to acquire the 3-D deformation, GPS and InSAR are commonly integrated, which was proposed by Ge et al. (2000) in the Double Interpolation and Double Prediction (DIDP) style. Gudmundsson et al. (2002) produced highresolution 3-D surface motion maps for Iceland by combining GPS and InSAR based on Markov random field regularization and simulated annealing. Samsonov and Tiampo (2006) and Samsonov et al. (2007), developed an analytical optimization strategy and applied
∗ Corresponding author at: China University of Petroleum (East China), Qingdao, 266580, China. E-mail addresses: [email protected] (X. Song), [email protected] (Y. Jiang), [email protected] (X. Shan), [email protected] (C. Qu). http://dx.doi.org/10.1016/j.jag.2016.12.019 0303-2434/© 2016 Elsevier B.V. All rights reserved.
it to reconstruct the 3-D deformation field of Southern California from InSAR and GPS measurements. Taking the interpolated GPS observations on the InSAR grid or a regular grid as a first guess map, these methods are strongly dependent on the availability and distribution of the GPS measurements. Due to exploitation of the densest continuous GPS networks in the world (327 GEONET sites are used), Hu et al. (2013) adopted the method of weighted least squares (WLS) and successfully retrieved the 3-D coseismic surface displacement field of the 2011 Tohoku-Oki earthquake on the onshore non-near-fault area by appropriately integrating the InSAR and GPS measurements. In the aforementioned methods, GPS data are generally interpolated by some pure mathematic interpolating algorithms without taking the mechanical property of the media into account; these methods are suitable for deforming areas with low-gradient and a sparse GPS network, but cannot deal with sudden and large ground movements caused by large earthquakes. A few researchers considered the physics of the deformation in interpolation. Based on the elastic theory, Guglielmino et al. (2011) proposed a SISTEM approach to obtain the 3-D displacement map of an active stratovolcano. However, this method is developed to study slow and continuous movement under the hypothesis of infinitesimal and homogeneous strain, such as volcanic activities and interseismic deformations, rather than coseismic deformations.
X. Song et al. / International Journal of Applied Earth Observation and Geoinformation 57 (2017) 104–112
Another key issue involved in merging heterogeneous observations is how to allocate the weight for different types of observations. The common practice is to weight the observations based on their a priori variance, but it is always difficult to know the exact variance of the observations. Hu et al. (2012) introduced the variance component estimation (VCE) approach to iteratively estimate the variances of the observations, and obtained the posterior estimation of the stochastic model of the InSAR measurements in the 3-D displacement velocity estimation from InSAR and GPS. However, the VCE approach expects enough redundant measurements and is easily end up with the negative variance problem. In addition, methods based on the least squares (LS) adjustment would also give biased estimation when gross errors exist in the observation samples. Large (moment magnitude (Mw) > 7) continental earthquakes in China often occur on the western fault systems, where the GPS density is very sparse (>50 km) although two GPS station construction projects (the Crustal Movement Observation Network of China: CMONOC-I, and the Tectonic and Environmental Observation Network of Mainland China: CMONOC-II) have been conducted in the past few decades. GPS observations are usually not dense enough to represent both the large and small scale coseismic deformation signal, which leads to significant problems in GPS interpolation (the first step in deriving 3-D surface displacement fields by integrating GPS and InSAR), especially for the deforming area with large spatial gradient in the near field of the causative fault. In this paper, we develop the fault dislocation model to invert a small number of GPS observations to recover a reasonable and reliable initial 3-D coseismic deformation field i. The Best Quadratic Unbiased Estimator (BQUE) (Sjöberg, 1983) is first used to estimate the variance components to avoid breaking the computation process due to the negative VC problem in VCE. Then GPS and InSAR are merged by a weighted least squares adjustment. Considering the contamination of atmospheric delays and other factors in the InSAR measurements, its stochastic model may not fit an exact normal distribution, so the gross error cannot be avoided. Therefore, a robust estimation method named IGG proposed by Zhou (1989) is applied to reduce the effect from gross errors. The 2008 Wenchuan earthquake is used as a case study to validate these two methods. 2. Methods Assuming a deformation field caused by an earthquake can be discretized into a m × n grid. InSAR and GPS observation is downsampled or interpolated on the grid at first. For each node of the grid, the LOS displacement observation dlos can be related to the real 3-D deformation components [ Dew model: dlos + Vlos = [ cos ˛ sin
− sin ˛ sin
T
Dud ] by the following
Dsn
− cos ] · [ Dew
Dsn
Dud ]
T
(1)
where ˛ is the azimuth angle of InSAR, and is the incidence angle. Vlos is the residuals of InSAR observation. The interpolated 3-D GPS observation [ xGPS yGPS tion components by: [ xGPS
yGPS
T
zGPS ] can be related to the 3-D deforma-
T
zGPS ] + VGPS = [ Dew
Dsn
Dud ]
T
(2)
Then, the observation equation can be formed (Eq. (3)), and to estimate the final 3-D displacement field through the least-squares (LS). L4mn×1 = D4mn×3mn X3mn×1 + V4mn×1 T
(3)
Where X = [Dew Dsn Dud ] is the unknown vector of displaceT ment on each node of the grid; L = [xgps ygps zgps dlos ] is the vector of GPS and InSAR observations; V is the vector of observation residuals; and D is the design matrix. Based on the Eq. (3), the
105
interpolated GPS and InSAR observations can then be used in a WLS model to retrieve the optimal 3-D displacement field. The key steps are (i) the interpolation of sparse GPS data appropriately, and (ii) the determination of the weighting scheme (estimating the posterior variances) for the observations. Here, we use sparse GPS data and fault dislocation model to produce a reasonable initial 3-D coseismic deformation (see Section 2.1), which will be further corrected by fusing InSAR observations. The weights for the observations are determined by using BQUE (see Section 2.2) and IGG (see Section 2.3) respectively. We will also compare the results from these two methods, and assess their accuracy by using field measurements. The overall processing strategy is also demonstrated in Fig. 1. 2.1. Recovering initial 3-D coseismic deformation field combining GPS with fault dislocation model To derive 3-D surface displacement fields by integrating GPS and InSAR, the first step is to resample the sparse GPS measurements onto the InSAR grid or another regular grid. Kriging interpolation is often used because it can preserve the theoretical spatial correlation of the measurements of adjacent data points. However, the distribution and density of the known points are always the key to maintain a good interpolating accuracy for general mathematic methods, e.g. the Kriging, bicubic spline interpolation and least squares collocation. Such methods are limited by the sparse GPS network in western China, since most faults in western China are located in remote areas where installing GPS stations is difficult and expensive. As a result, pure mathematic interpolation methods are not capable of reproducing a reliable deformation field caused by a strong earthquake. In this study, we seek for solutions that combines the physics of earthquake deformation and mathematical methods. The elastic dislocation model derived by Okada (1985, 1992) provides a good representation for surface and internal deformation generated by shear and tensile faults in a homogenous half-space or a multilayered half-space (Wang et al., 2003). It is usually used to invert the surface observations for the focal parameters and slip distribution at certain depth (Delouis et al., 2002; Feigl et al., 2002; Zhang et al., 2011). These earthquake parameters can provide us with very useful knowledge of the geophysical properties of an earthquake. Based on the inverted fault model, we can simulate the coseismic deformation in any direction. However, the deformation field from forward modeling cannot fully represent the surface deformation, regardless of the observations used in the inversion. This is because of the assumption and simplification on fault model (presumptions about the kinematics (e.g. slip and rake) and geometry (e.g. dip and strike)), inverting algorithms that we adopt, and spatial density of geodetic observations with uncertainties. A fine and accurate 3-D coseismic deformation field can not only be used to constrain earthquake source parameters, but also allows for investigations of other geophysical and geological problems such as shallow slip deficit and fault segmentation. Practically, coseismic ruptures and deformation are also important for post-disaster rescue. Observations of the earth’s surface are the most direct and effective tool to derive the surface information. In the case of insufficient observational data, a geophysical model considering the mechanical property of the media can be used to assist the reconstruction of the 3-D surface deformation conversely. This is also a new attempt applied in other fields (3-D surface displacement caused by Geothermal field (Hu et al., 2016), source fluid volumetric change of groundwater depletion (Castellazzi et al., 2016). The dislocation model describes the physical relationship between the displacement on faults and the deformation on the surface, regardless of the distribution and density of the surface geodetic observations. Here, we use the dislocation model to generate the first approximation of a 3-D coseismic deformation field,
106
X. Song et al. / International Journal of Applied Earth Observation and Geoinformation 57 (2017) 104–112
Fig. 1. The flowchart of integrating InSAR and GPS to derive the 3-D coseismic deformation field.
Fig. 2. Flow chart showing the procedure of using the dislocation model and a few GPSs to generate the first approximation of a 3-D coseismic deformation field.
with a small amount of GPS observations. The procedure is outlined in below and also in Fig. 2. Step 1 constructing the initial fault model: Fault segments can be mapped based on visual analysis of high-resolution imagery, and displacement field from InSAR or image correlation. Fault dip of each segment can be obtained from seismic data. The location, length and strike of each fault segment are adjusted iteratively in step 2. Step 2 inverting for slip distribution: Based on the initial fault model constructed in step 1, sparsely distributed GPS data are used to estimate the slip distribution on the fault plane, as well as the location, length and strike of each fault segment simultaneously. This is carried out iteratively using the SDM modeling program (Wang et al., 2009). The effect of GPS density on the inverted results will be tested in Section 3, aiming to find a stable solution with as less GPSs as possible. Step 3 synthesizing 3-D surface deformation field: based on the slip distribution and fault parameters from step 2, we could synthesize a 3-D coseismic deformation field on a regular grid by forward modeling.
vations is sufficiently high. The result also depends on the level of realistic of a priori adjustment and stochastic models, and the quality of the initial VCs (Eshagh and Sjoberg, 2008). The real cases usually do not satisfy all of these requirements. Furthermore, the strict Helmert VC estimation expects independent observational groups and enough redundant measurements, and is easy to suffer from the negative variance problem. As for coseismic deformation, InSAR measurements are usually contaminated by the atmospheric and orbital errors more or less, and the number of redundant observation is only in LOS direction, relative to three unknowns (two horizontal components and one vertical component) for each point. The Helmert VCE generally cannot well handle this situation. Thus, in this case we propose to use the Best Quadratic Unbiased Estimator from the functional model of conditional adjustment with unknown parameters (Eshagh and Sjoberg, 2008) to estimate the VC of InSAR and GPS iteratively. Given the following Gauss-Helmert model,
2.2. The BQUE of VC for InSAR and GPS
where A and B are the first and second design matrices, X is the unknown vector, ε is the residual vector, W is the misclosures vector, Q is the covariance matrix of observations, and E {• } is the statistical expectation.
Most of the VC estimation methods work well in the case where no bias exists in the observations and the redundancy of the obser-
AX + Bε = W,
with
εεT
=Q
and
E {ε} = 0
(4)
X. Song et al. / International Journal of Applied Earth Observation and Geoinformation 57 (2017) 104–112
Eq. (4) can be converted to the conditional model by adjustment multiplying both sides of the equations by I − A0 , namely
I − A0 Bε = I − A0 W,
¯ =W ¯ Bε
or
(5)
where A0 = AA− , A− stands for a specific generalized inverse of A. For the conditional adjustment (Eq. (5)), the best quadratic unbiased estimation of variance components can be written as (Eshagh and Sjoberg, 2008): ˆ = S − uT
Sij = trace RCi RCj , R=C C=
q
0
I−A
i2 Ci
0 T
C −1
Ci = BQi BT
i=1
i = 1, 2, . . ., q.
In this passage, we construct one condition adjustment equation by 4 observations, dlos , xGPS , yGPS , zGPS , on each grid point:
⎡
−1 cos ˛ sin
− sin ˛ sin
− cos
dlos
⎢ ⎢ xGPS ⎢ ⎢y ⎣ GPS
⎤
⎥ ⎥ ⎥=0 ⎥ ⎦
|v| < a
1 p (v) = ⎪ ⎩ |v| + k 0
(9)
a < |v| < b |v| ≥ b
where a and b are constant values, and depend on the observations. Empirically, we use a = 1.5, b = 2.5; k= 0.001 in the case of v = 0.
3.1. Data
(8)
YGPS 1
···
Hence, the observation equation for N points is, =0
The 12 May 2008 Mw 7.9 Wenchuan earthquake is one of the most destructive events that had ever occurred in China during the last 50 years, killed about 70,000 people (Zhang et al., 2008). This event ruptured the NE-SW striking, north-west dipping Longmenshan fault belt, which was formed by continuing eastward extrusion of the Tibet plateau to the western edge of the Sichuan Basin (Zhang et al., 2008). Many geological and geodetic studies have focused on this event and have provided abundant material to test the effectiveness of our method.
The JAXA arranged a mission of emergency and obtained a batch of L-band SAR data (ALOS/PALSAR) soon after the Wenchuan earthquake. Tracks 471–477 covered the whole surface rupture zone (Fig. 3(a)). The data were processed using GAMMA (Qu et al., 2008). Good coherence is preserved for most of the area, except for the near-fault zone where the decorrelation is caused by large surface deformation. The near-fault displacement field was therefore achieved by using the offset-tracking technique, based on which
(7)
zGPS
B• L
⎧ ⎪ ⎨1
3. Experiments
i, j = 1, 2, . . ., q,
= I−A
and
¯ T RCi RW ¯ , ui = W
absolute value of residual |v|, depending on the quadratic mean deviation (mean square error) :
(6)
where S − is a generalized inverse of S and
−1
107
where B is the design matrix
and L is the observation vector
L=
dlos 1
dlos 2
···
dlos N
XGPS 1
XGPS 1
···
XGPS N
YGPS 1
Generally,W = −B• L = / 0, because of various errors, including orbital and atmospheric errors in GPS and InSAR. We develop an iterative procedure initialing with a priori estimates of the variance components, to caulculate the variance
of the unit weight 0 2 = W T BQBT
−1
W and variance compo-
nents using Eq. (6), until0,i 2 = 0,i+1 2 . In the end, we obtain the estimation of 1 2 , 2 2 , 3 2 , 4 2 for the variance of observation dlos , xGPS , yGPS , zGPS . 2.3. The robust estimation by IGG Due to the simplicity and general good performance, the least squares (LS) adjustment has been extensively applied in geodetic data processing under the hypothesis that observation samples are subject to normal distribution. Nevertheless, gross errors often exist in observations, which would bias the result from the LS adjustment with a bigger variance. Robust estimation aims at decreasing the effect of gross errors. Krarup et al. (1980) introduced the robust estimation theory to surveying by proposing the Denmark method, based on the M estimation theory. Zhou (1989) put forward the concept of equivalent weights, made LS method formalization for M estimation by equivalent weights, and named the algorithm as IGG. In the LS adjustment (Eq. (3)), we can determine the weight p (v) of the
YGPS N
ZGPS 1
···
ZGPS 1
ZGPS N
T
we can map the full picture of the rupture trace. GPS data are from Wang et al. (2011), which comes from various sources, including CMONOC, the national key basic research and development project, and the Earthquake Administrations of Sichuan and Chongqing. Fig. 3(b) shows the location of GPS data used in the rest of the computation and accuracy control in this paper. 3.2. Calibration of InSAR data using GPS Depending on the acquisition time of the slave image, the interferometric phase contains the post-seismic deformation phase to some extend even in the far field, as well as the atmospheric effect, both biased the baseline refinement in InSAR data processing. Furthermore, unlike referenced GPS data, InSAR only gives us a relative measurement for the surface deformation. So it is necessary to make certain corrections to the spatial reference of InSAR measurements before integrating InSAR and GPS. • the first step is to project the displacement vector of each GPS station: u = Llos = s · u
ue
un
uu
T
into the InSAR LOS direction: (10)
108
X. Song et al. / International Journal of Applied Earth Observation and Geoinformation 57 (2017) 104–112
Fig. 3. (a) The study area: black vectors show the GPS velocities; black rectangles indicate the coverage of 7 PALSAR stripes; red lines show the surface rupture of the Wenchuan earthquake. (b) Blue rectangles show the fault geometry we used in our inversion. Circles stand for the location of the GPS points (the color indicates the serial number of GPS we used in calculating the deformation field when we change the number of GPS input). Black stars indicate the GPSs for checking. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
where the unit vector s = cos ˛ sin − sin ˛ sin − cos , ˛ is the azimuth angle of InSAR, and is the incidence angle. • Use a quadratic polynomial to construct the relationship between the InSAR and GPS displacements in the LOS direction, and then estimate the coefficients for each track, n = am2 + bm + c
plan
Interpolating method
Weighting method
I II III IV
Dislocation model Dislocation model Kriging Kriging
IGG BQUE IGG BQUE
(11)
Where a, b, c are the coefficients, and m, n are the LOS InSAR and GPS displacements in the LOS direction at each GPS station. • Make corrections to InSAR observations: y = ax2 + bx + c
Table 1 Combinations of different interpolating and weighting methods.
(12)
where x is the original InSAR observation, and y is the corrected value. • Mosaic 7 tracks to obtain a complete InSAR deformation field.
the focal mechanism solutions of the main shock and aftershocks (Zhang et al., 2011). • GPS interpolation. Following the procedure described in Section 2.1, GPS data is used to synthesize the 3-D surface deformation field of the Wenchuan earthquake on the regular grid by forward modeling. • Integrating InSAR and GPS. The WLS adjustment is used to integrate InSAR and GPS on each wall, and the BQUE and IGG applied to weight the observations. The results on the hanging wall and footwall are mosaicked to generate a single complete displacement field. 10 selected GPS points are used to evaluate the accuracy of the results.
3.3. Processing flowchart 4. Results and discussion The data processing scheme of the Wenchuan earthquake is similar to that in Fig. 1. The only difference is that we need to carry out data processing for hanging wall and footwall separately, and then mosaic their results in the final step. The detail description is as follow: • Downsample the InSAR data. Within the common coverage of the InSAR and GPS data, we confine our study area to the range of 30.5-32.5N and 103.2-105.2E, and grid it with a spacing of 0.2 ◦ in latitude and longitude. The corrected InSAR deformation field is downsampled onto this regular grid using bilinear interpolation. • Construct the fault model for the Wenchuan earthquake. The displacement field from offset-tracking shows the rupture trace of the Wenchuan earthquake clearly, based on which we estimate the fault geometry parameters such as the location, length and strike of each segment. The dip and rake of each segment are from
Here we combine different interpolating and weighting methods for comparison (Table 1). Following the procedure described in previous sections, we calculate the corresponding 3-D coseismic deformation fields and associate RMSEs from 10 check points. Fig. 4 shows the 3-D deformation fields from 4 combinations by using 30 (a) and 80 (b) GPS points as input. Note that the selected GPS points in each inversion and accuracy assessment should be distributed evenly over the whole deforming area, to minimize the impact of GPS distribution on the results. As shown in Fig. 4, we can see the deformation fields from plans I and II present more details, and less dependence on the number of GPSs compared with those from plans III and IV. This is because Kriging requires an appropriate theoretical semi-variogram model for each component to be interpolated. It is one of the main critical problems in geostatistics. Given the large area of the Wenchuan rupture, tens of GPS points
X. Song et al. / International Journal of Applied Earth Observation and Geoinformation 57 (2017) 104–112
are too few to describe the whole deformation area using Kriging. Instead, plans I and II, which use the fault dislocation model that accounts for the medium property and deformation mechanism, show reasonable and stable results. To assess the accuracy of our methods quantitatively, we used field observations along the rupture from Xu et al. (2008). The field measurements of horizontal and vertical fault displacements on a few points are converted into the EW, SN, and vertical components assuming NE45 striking. Figs. 5 S2 and S3 shows the comparison between the field measurements and our model simulations. The RMSEs of the Plan I (dislocation model + IGG) and II (dislocation model +BQUE) with 80 GPS points are 47.8 cm, 18.3 cm, 21.5 cm
109
and 64.8 cm, 49.8 cm, 73.7 cm for EW, SN and vertical direction respectively. Clearly, both methods improved the accuracy compared to the “Kriging + IGG” result with the accuracy of 113.9 cm, 65.3 and 202.5 cm. It indicates the near-field displacements derived from the method that integrates the dislocation model agree better with the field measurements than those from the method using Kriging interpolation. This also implies that the interpolation based on the elastic dislocation model is better than Kriging. Considering the maximum on-fault displacement is over 5m, the result with an error of less than 50 cm (Fig. 5) is satisfying. In terms of the weighting methods, the IGG works better than the BQUE based on the RMSEs. The difference lies in the fact that
Fig. 4. The 3-D coseismic deformation fields of the Wenchuan earthquake from 4 plans listed in Table 1, with (a) 30 and (b) 80 GPS points as input.
110
X. Song et al. / International Journal of Applied Earth Observation and Geoinformation 57 (2017) 104–112
Fig. 4. (Continued)
the BQUE assumes that errors are subject to the normal distribution whereas the IGG does not; as the result, the BQUE is easier to be distorted by outliers and thus requires careful elimination of outliers. The IGG shows a strong ability to decrease the weight of outliers, which is an advantage in cases where observations contain gross errors like in the Wenchuan earthquake. Our experiments demonstrate the advantages of using the fault dislocation model in interpolation that only a small amount of GPS observations is needed. It suggests that our method is very useful in many tectonically active regions where GPS stations are sparsely distributed. In order to test the sensitivity of the interpolation method to the number of GPS stations, we also vary the number
of GPS observations from 10 to 80 in steps of 10 in the slip distribution inversion. The RMSEs of the three components of surface displacement are calculated using 10 check points (the distribution is shown in Fig. 2(b)). Fig. 6 shows the corresponding RMSE in the EW, SN, and vertical directions at different number of GPS points. The RMSEs from plans I and II remain stable when using more than 20 GPS points. While the RMSEs from plans III and IV would be equal to or slightly smaller than those from plans I and II, only when the number of GPS input reaches 80. In case of the number of GPS check points increase to 30, the RMSE curves (Fig. S1) show similar statistical results. The difference is that the RMSEs from plans III and IV tend to be more stable than those from plans I
X. Song et al. / International Journal of Applied Earth Observation and Geoinformation 57 (2017) 104–112
111
Fig. 5. Field and modeled measurements of EW, SN and vertical fault displacements. Black triangles: field measurements vs. Plan I (IGG), and its RMSEs are 47.8 cm, 18.3 and 21.5 cm for EW, SN and vertical direction respectively. Green circles: field measurements vs. the Plan III (Kriging) with 80 GPS points, and its RMSEs are 113.9 cm, 65.3 and 202.5 cm for EW, SN and vertical direction respectively. The red line denotes the line ‘y=x’. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
Fig. 6. The RMSEs calculated by using 10 GPS points in (a) EW, (b) SN, (c) vertical directions with an increasing number of the input GPS points.
and II when the number of GPS input reaches 50. This indicates that we can decrease the number of GPS points significantly by incorporating the fault dislocation model and still obtain a good result, similar to that by using 80 GPS points and Kriging method. for a large deforming area (500 × 500 km) such as the Wenchuan earthquake, 20 evenly distributed GPS points are sufficient enough. This provides us information to build economical and applicable plan for those tectonically active regions with sparse GPS network. We also notice that the RMSEs (Fig. 6) calculated by using 10 GPS check points from plans III and IV have increased first and then
began to drop when the number of GPS points increased from 10 to 30. We think that it is because a small number of GPS points (e.g., ≤ 10) are not enough to present the characteristics of the widerange deformation field of the Wenchuan earthquake when using Kriging. In such case the distance between the input points and the check points plays a dominant role in calculating the RMSEs, thus RMSEs lost its statistical meaning. Fig. S1 did not show such fluctuation in the beginning, because a larger number of GPS check points (number of 30) are used to calculate the RMSEs, which reduce the dependence on the distance between the input points and the check
112
X. Song et al. / International Journal of Applied Earth Observation and Geoinformation 57 (2017) 104–112
points. In addition, there is no doubt that the accuracy for the Kriging will decrease when more GPS points are chosen as check points, which means less points are chosen for interpolation. This is the reason that the RMSEs from plans III and IV are always larger than that from plans I and II in Fig. S1. 5. Conclusions The derivation of 3-D coseismic deformation field is very important for revealing detail movement of the fault and the area affected by an earthquake. GPS observations, even for a dense GPS network, are usually not enough to recover the real 3-D coseismic deformation field (especially in the near-fault area) with large spatial coverage like that in Wenchuan earthquake example. An economical and applicable approach has been developed in this paper that can retrieval 3-D coseismic deformation from InSAR and small number of GPS data by combining the elastic dislocation model with the BQUE or the robust estimation method IGG. Compared with the pure mathematic interpolator discussed in this study, the elastic dislocation model that considers the mechanical property of the media allows us to estimate initial 3-D coseismic deformation in a more physically reasonable and accurate style. By using the BQUE and IGG, the weights of the InSAR and GPS observations have been determined optimally, in the case where large atmospheric and orbital errors exist in the InSAR observation and the redundancy of observations is not sufficiently high. An important advantage of presented approach is that the requirement for the density of GPS is significantly lowered to reconstruct a reliable 3-D coseismic deformation field. It is very useful for remotely located tectonically active regions with a sparse GPS network. The proposed method has been validated by using a series of experiments with the test case of Wenchuan earthquake. The results have clearly shown that our method is sound and can significantly improve the accuracy of the 3-D coseismic deformation reconstruction. As our method is generic, it is also possible to integrate other sources of surface movement measurements including Multiple-Aperture InSAR (MAI), offset-tracking and leveling. Overall, the method offers extensive future usage for measuring earthquake deformation. Acknowledgments This work was supported jointly by the Basic Scientific Funding of Institute of Geology, China Earthquake Administration (No.IGCEA1404), grants of State Key Laboratory of Earthquake Dynamics, China Earthquake Administration (No.LED2013A02), and the National Natural Science Foundation of China (No. 41461164002). All ALOS PALSAR data were provided and are copyrighted by JAXA. Most figures were made using the public domain Generic Mapping Tools. We thank Prof. Mehdi Eshagh in University West for instructions of the usage of the BQUE. References Castellazzi, P., Martel, R., Galloway, D.L., Longuevergne, L., Rivera, A., 2016. Assessing groundwater depletion and dynamics using GRACE and InSAR: potential and limitations. Ground Water, http://dx.doi.org/10.1111/gwat. 12453. Catani, F., Farina, P., Moretti, S., Nico, G., Strozzi, T., 2005. On the application of SAR interferometry to geomorphological studies: estimation of landform attributes and mass movements. Geomorphology 66, 119–131. Dehghani, M., Zoej, M.J.V., Entezam, I., Mansourian, A., Saatchi, S., 2009. InSAR monitoring of progressive land subsidence in Neyshabour, northeast Iran. Geophys. J. Int. 178, 47–56. Delouis, B., Giardini, D., Lundgren, P., Salichon, J., 2002. Joint inversion of InSAR, GPS, teleseismic, and strong-motion data for the spatial and temporal distribution of earthquake slip: application to the 1999 Izmit mainshock. Bull. Seismol. Soc. Am. 92, 278–299, http://dx.doi.org/10.1785/0120000806.
Elliott, J.R., Copley, A.C., Holley, R., Scharer, K., Parsons, B., 2013. The 2011 Mw 7.1 Van (Eastern Turkey) earthquake. J. Geophys. Res. 118, 1619–1637, http://dx. doi.org/10.1029/2012JB009569. Eshagh, M., Sjoberg, L., 2008. The modified best quadratic unbiased non-negative estimator (MBQUNE) of variance components. Stud. Geophys. Geod. 52, 305–320. Feigl, K., Sarti, F., Vadon, H., 2002. Estimating slip distribution for the I˙ zmit mainshock from coseismic GPS, ERS-1, RADARSAT, and SPOT measurements. Bull. Seismol. Soc. Am. 92, 138–160. Funning, G.J., Parsons, B., Wright, T.J., Jackson, J.A., Fielding, E.J., 2005. Surface displacements and source parameters of the 2003 Bam (Iran) earthquake from Envisat advanced synthetic aperture radar imagery. J. Geophys. Res. 110, B09406, http://dx.doi.org/10.1029/2004JB003338. Ge, L., Han, S., Rizos, C., 2000. The double interpretation and double prediction (DIDP) approach for InSAR and GPS integration. In: ISPRS Commission Reports, the XIXth Congress of the International Society for Photogrammetry and Remote Sensing, Amsterdam, Netherlands, pp. 205–212. Gudmundsson, S., Sigmundsson, F., Carstensen, J., 2002. Three-dimensional surface motion maps estimated from combined interferometric synthetic aperture radar and GPS data. J. Geophys. Res. 107, http://dx.doi.org/10.1029/ 2001jb000283. Guglielmino, F., Nunnari, G., Puglisi, G., 2011. Simultaneous and integrated strain tensor estimation from geodetic and satellite deformation measurements to obtain three-dimensional displacement maps. IEEE Trans. Geosci. Remote Sens. 49, 1815–1826. Hooper, A., Segall, P., Zebker, H., 2007. Persistent scatterer interferometric synthetic aperture radar for crustal displacement analysis, with application to Volcano Alcedo, Galapagos. J. Geophys. Res.-Solid Earth 112, B07407. Hu, J., Li, Z.W., Sun, Q., Zhu, J.J., Ding, X.L., 2012. Three-dimensional surface displacements from InSAR and GPS measurements with variance component estimation. IEEE Geosci. Remote Sens. Lett. 9, 754–758. Hu, J., Li, Z.W., Ding, X.L., Zhu, J.J., Sun, Q., 2013. Derivation of 3-D coseismic surface displacement fields for the 2011 Mw 9.0 Tohoku-Oki earthquake from InSAR and GPS measurements. Geophys. J. Int. 192, 573–585, http://dx.doi.org/10. 1093/gji/ggs033. Hu, J., Wang, Q., Li, Z., Zhao, R., Sun, Q., 2016. Investigating the ground deformation and source model of the yangbajing geothermal field in tibet, China with the WLS InSAR technique. Remote Sens. 8, 191, http://dx.doi.org/10.3390/ rs8030191. Krarup, T., Juhl, J., Kubik, K., 1980. Gotterdammerung over least squares adjustment. In: Proceeding of 14th Congress of the International Society of Photogrammetry, Hamburg, Germany, B 3, pp. 369–378. Li, Z., Elliott, J.R., Feng, W., Jackson, J.A., Parsons, B., Walters, R.J., 2011. The 2010 Mw 6.8 yushu (Qinghai, China) earthquake: constraints provided by InSAR and body wave seismology. J. Geophys. Res. 116, B10302, http://dx.doi.org/10. 1029/2011JB008358. Okada, Y., 1985. Surface deformation due to shear and tensile faults in a half-space. Bull. Seismol. Soc. Am. 75, 1135–1154. Okada, Y., 1992. Internal deformation due to shear and tensile faults in a half-space. Bull. Seismol. Soc. Am. 82, 1018–1040. Qu, C., Song, X., Zhang, G., 2008. Analysis on the characteristics of InSAR coseismic deformation of the Ms8.0 Wenchuan earthquake. Seismol. Geol. 30, 1076–1083. Samsonov, S., Tiampo, K., 2006. Analytical optimization of a DInSAR and GPS dataset for derivation of three-dimensional surface motion. IEEE Geosci. Remote Sens. Lett. 3, 107–111. Samsonov, S., Tiampo, K., Rundle, J., Li, Z.H., 2007. Application of DInSAR-GPS optimization for derivation of fine-scale surface motion maps of Southern California. IEEE Trans. Geosci. Remote Sens. 45, 512–521. Sjöberg, L.E., 1983. Unbiased estimation of variance-components in condition adjustment with unknowns-a MINQUE approach. Zeitschrift für Vermessungen 108, 382–387. Wang, R., Martin, F.L., Roth, F., 2003. Computation of deformation induced by earthquakes in a multi-layered elastic crust: FORTRAN programs EDGRN/EDCMP. Comp. Geosci. 29, 195–207. Wang, L., Wang, R., Roth, F., Enescu, B., Hainzl, S., Ergintav, S., 2009. Afterslip and viscoelastic relaxation following the 1999 M 7.4 Izmit earthquake from GPS measurements. Geophys. J. Int. 178, 1220–1237. Wang, Q., Qiao, X., Lan, Q., Freymueller, J., Yang, S., Xu, C., Yang, Y., You, X., Tan, K., Chen, G., 2011. Rupture of deep faults in the 2008 Wenchuan earthquake and uplift of the Longmen Shan. Nat. Geosci. 4, 634–640. Wright, T.J., Parsons, B.E., Lu, Z., 2004. Toward mapping surface deformation in three dimensions using InSAR. Geophys. Res. Lett. 31. Xu, X.W., Wen, X., Ye, J.Q., 2008. The Ms 8.0 Wenchuan earthquake surface ruptures and its seismogenic structure. Seismol. Geol. 30, 597–629. Zhang, P.Z., Xu, X.W., Wen, X.Z., Ran, Y.K., 2008. Slip rates and recurrence intervals of the Longmen Shan active fault zone and tectonic implications for the mechanism of the May 12 Wenchuan earthquake, 2008, Sichuan. China. Chin. J. Geophys. 51, 1066–1073. Zhang, G., Qu, C., Shan, X., Song, X., Zhang, G., Wang, C., Hu, J., Wang, R., 2011. Slip distribution of the 2008 Wenchuan Ms 7.9 earthquake by joint inversion from GPS and InSAR measurements: a resolution test study. Geophys. J. Int. 186, 207–220. Zhou, J.W., 1989. Classical theory of errors and robust estimation. Acta Geodetica et Cartographica Sinica 18, 115–120.