3D Coseismic Deformation Field and Source Parameters of the 2017 Iran-Iraq Mw7.3 Earthquake Inferred from DInSAR and MAI Measurements

The coseismic slip on the main fault related to the 2017 Iran-Iraq Mw7.3 earthquake has been investigated by previous studies using DInSAR (differential interferometric synthetic aperture radar) ground deformation measurements. However, DInSAR observation is not sensitive to the ground deformation in the along-track (AT) direction. Therefore, only the one-dimensional (1D) DInSAR coseismic deformation field measurements, derived in the LOS (line-of-sight) direction of radar, was applied in source parameters estimation. To further improve the accuracy of the fault slip inversion, the 3D (three-dimensional) coseismic deformation fields were reconstructed in the first place, by a combined use of the DInSAR and MAI (multiple aperture InSAR) measurements. Subsequently, the LOS and 3D deformation data sets were used as the constraint respectively, to perform a two-step inversion scheme to find an optimal geometry and slip distribution on the main fault. The comparative analysis indicated that the 3D coseismic deformation data sets improved the inversion-accuracy by 30%. Besides, the fault invention results revealed a deep dislocation on a NNW trending fault (the strike is 352.63◦) extending about 60 km, along the fault dips 14.76◦ to the ENE. The estimated seismic moment is 8.44 × 1019 Nm (Mw7.3), which is close to the solution provided by USGS (United States Geological Survey). The slip distributed at the depth between 12 and 18 km, and the peak slip of 6.53 m appears at the depth of 14.5 km left near the epicenter. Considering the geological structure in the earthquake region and fault source-parameters, it comes to a preliminary conclusion that the ZMFF (the Zagros Mountain Front fault) should be responsible for the earthquake. In general, this paper demonstrated the superiority of using the 3D coseismic deformation fields on source parameters estimation.


Introduction
The 2017 Mw7.3 Iran-Iraq Earthquake occurred in a region where the pattern of major plate convergence is well constrained (Figure 1).Although some studies and institutions (USGS (United States Geological Survey), GCMT (Global Centroid Moment Tensor), CENC (China Earthquake Networks Center)) published some source parameters [1][2][3][4][5][6][7], the variable slip on the main fault related with this earthquake is still uncertain.In order to accurately derive the variable slip on main fault, this present conducted detailed fault inversion using the 3D coseismic deformation fields reconstructed by a combined use of the DInSAR and Multiple Aperture InSAR (MAI) measurements, based on both the ascending and descending SAR images obtained by L-band ALOS-2 and C-band Sentinel-1A sensors respectively.
1 lists the acquisition date of each SAR image and the perpendicular baseline (B⊥), temporal baseline (BT) and orbit types of each interferometric pair.The ALOS-2 ScanSAR mode images were collected along both the ascending and descending orbit at a radar incidence angle of 39° in HH polarization and generated with a pixel spacing of 19.02 m in slant range and 25.9 m in azimuth (i.e., along track).Each image scene covers an area of about 350 km by 350 km in azimuth and range, respectively.
Meanwhile, the Senitnel-1A IW mode images were collected along both the ascending and descending orbits at radar incidence angle of 39.5° in VV polarization and generated with a pixel spacing of 5 m in slant range and 20 m in azimuth.Each image scene covers an area of about 250 km by 170 km in range and azimuth, respectively.The heading angles of the ALOS-2 and sentinel-1A ascending orbits are 347° and 344° (measured clockwise from the north), respectively; and the heading angles of the descending orbits are 195° and 193° (measured from the north), respectively.Nowadays, the geodetic data of GPS and DInSAR have been extensively used in earthquake fault inversion.Although coseismic GPS data is good at providing 3D surface displacements at individual observation locations, labor power, and material resources cost is great.Meanwhile, the spatial density of GPS sites is limited.In recent years, the DInSAR method is developing with continually launching of new SAR satellites.Because it has a few significant advantages [8], such as high resolution, shorter revisit interval, wider coverage, etc., that is widely applied in fault slip inversion at present.However, the DInSAR observation is only sensitive to displacements in the line-of-sight (LOS) direction of radar.Therefore, only one component of the surface displacement vector was brought in the fault slip inversion.As a result, uncertainties of fault source parameters arose from data resources and data incompleteness, and the accuracy of source parameters estimation is restricted.
To overcome the limitation, the Multiple Aperture InSAR (MAI) measurements was widely used to provide valuable information about azimuthal surface displacements [8][9][10][11][12].Considering the heading angle of the ascending-and descending-pass direction is varying, totally four different displacements vectors can be observed by combing the DInSAR and MAI measurements.As a result, the 3D deformation fields can be reconstructed [13][14][15][16][17][18][19][20][21] to provide accurate data for fault slip inversion and increase the inversion efficiency.
The previous study [22] measured the AT (along-track) deformation of the earthquake by the Azimuth Pixel Offset (AZO) and MAI method, respectively.The comparison results of two methods show that the MAI method is better at observing the AT (Along Track) deformation.At the same time, the 3D deformation fields of the earthquake were reconstructed by combing the DInSAR and MAI method [22].But the low frequency (e.g., L-band SAR) is more sensitive to ionosphere [23][24][25], so it was necessary to remove the relative azimuth shift caused by ionosphere to obtain better deformation fields.
In this paper, the SAR data sets of ALOS-2 and Sentinel-1a are used firstly to derive the deformation fields around the epicenter.During the procedures, we used the Precise Orbit Ephemerides for all images acquired and applied the range split-spectrum method [26,27] to estimate the differential ionospheric phase, and then mitigate the relative azimuth shift caused by ionosphere in the MAI interferogram.The accuracy of the 3D coseismic deformation fields was promoted.Then this paper focus on testing the effect of using the 3D coseismic deformation fields on fault slip inversion.In addition, the fault slip inversion will be implemented by the constrain of the LOS and 3D coseismic deformation fields, respectively.In the end, comparative analysis was performed to evaluated the results.The results demonstrate the advantage of 3D coseismic deformation fields in fault slip inversion.It should be pointed out that such method is applied for the 2017 Iraq Mw7.3 earthquake, for the first time.

Backgroud of the Study Area
The 2017 Mw7.3 Iran-Iraq earthquake stuck the Iraq-Iran border region on November 12, 2017 with an epicentral of [45.9 • E, 34.93 • N], 50 km north east of Sarpol Zahab city based on the earthquake report from the USGS (the United States Geological Survey). Figure 1 shows the tectonic in this region, the continental collision between Arabia and Eurasian plate is the driving force.The oblique convergence between the Arabia and Eurasia plates activate through a combination of strike-slip (~5 mm/y in the central part of the Talesh Mountains) and thrust faulting (~2 mm/y-6 mm/y in the northern and southern parts of the mountains), respectively [28,29].Geologically, the earthquake's epicenter was located near the Zagros Mountain Front Fault (ZMFF) where the Zagros Mountains have been uplifting since the Late Eocene by the driving force.As a result, this earthquake was triggered and characterized by thrust motion.
This earthquake event has already been studied several times by using different data sets, resulting in different fault models [1][2][3][4][5][6][7].These studies indicate that this earthquake was a shallow-slip thrusting event with a small amount of right-lateral strike slip.The major fault parameters including strike, dip, and rake are similar.But the maximum slip is quite different.Based on the deformation fields obtained by the data sets of ALOS-2, Kobayashi and Barnhart indicated that the maximum slip is approximately 5 m [1,5].Taymaz [2] indicated that the maximum slip is approximately 6.9 m based on Sentinel-1A interferometric pairs and teleseismic P-and SH-body-waveforms data.Wang indicated that the maximum slip is about 7 m while Feng indicated the maximum slip is 6 m (3,7).

Data
Figure 1 shows shaded relief map of study area and image coverage, etc.The four interferometric pairs we selected for this study, include two pairs of L-band (radar wavelength of 23.6 cm) ALOS-2 images and two pairs of C-band (radar wavelength of 5.6 cm) Sentinel-1A images (see Figure 1).Table 1 lists the acquisition date of each SAR image and the perpendicular baseline (B ⊥ ), temporal baseline (B T ) and orbit types of each interferometric pair.The ALOS-2 ScanSAR mode images were collected along both the ascending and descending orbit at a radar incidence angle of 39 • in HH polarization and generated with a pixel spacing of 19.02 m in slant range and 25.9 m in azimuth (i.e., along track).Each image scene covers an area of about 350 km by 350 km in azimuth and range, respectively.Meanwhile, the Senitnel-1A IW mode images were collected along both the ascending and descending orbits at radar incidence angle of 39.5 • in VV polarization and generated with a pixel spacing of 5 m in slant range and 20 m in azimuth.Each image scene covers an area of about 250 km by 170 km in range and azimuth, respectively.The heading angles of the ALOS-2 and sentinel-1A ascending orbits are 347 • and 344 • (measured clockwise from the north), respectively; and the heading angles of the descending orbits are 195 • and 193 • (measured from the north), respectively.
Around the epicenter, the vegetation coverage is low.At the same time, the temporal and spatial baselines are small enough (see Table 1), thus the coherence is high enough.The DInSAR and MAI methods will be processed effectively to reconstruct the 3D coseismic deformation caused by the earthquake.In ALOS-2 ScanSAR mode, SAR antenna cyclically points to several subswaths with different look angles.In each cycle, the SAR system collects a group of echoes for each subswath.The group of echoes is called a burst.By repeating this cycle over several subswaths, a wide-swath image can be acquired [30].The earthquake region was fully imaged by a burst.So, considering the efficiencies, only one burst of ALOS-2 images and whole Sentinel-1A images were used, actually.

Geodetic Data Sets
In order to provide accurate data for fault inversion, we use two result sets.Four DInSAR results provide spatially dense LOS deformation field along both the ascending and descending imaging direction.Two DInSAR results and two MAI results based on ALOS-2 interferometric pairs will be combined to reconstruct the 3D coseismic deformation-field.

Extraction of the LOS and Azimuth Coseismic Deformation Field
The previous study [22] used the conventional DInSAR method based on four images pairs of ALOS-2 and Sentinel-1A to extract the LOS deformation field.The same procedures were performed step by step.We used the Precise Orbit Ephemerides for all images acquired.A digital elevation model (DEM) derived from the Shuttle Radar Topography Mission (SRTM) with 3 arcsec resolution [31] was used to remove the topographic effects from interferograms.The DEM data of the study area is depicted as a shaded relief map in Figure 1.In order to reduce the influence of noise, the method of Goldstein [32] is used to filter interferograms.Then, for avoiding the accuracy loss caused by orbit error and correcting the orbit and phase shift, refinement and re-flattening is processed by selecting enough and high coherence reference points, which distribute uniformly and far from the deformation field in the interferograms.And phase unwrapping was processed by the minimum-cost flow (MCF) method [33].In addition, we performed the range split-spectrum method to ALOS-2 ScanSAR interferogram to estimate the differential ionospheric phase [30].
Figure 2 shows the LOS coseismic deformation fields measured by the DInSAR method based on four interferometric pairs.It has two clear deformation lobes (here refers to lobe G and H) which divided by an NNW trending bound (see black line AB in Figure 2).As the Figure 2a,c shows, the maximum value of uplift in lobe G and subsidence in lobe H are approximately 100 cm and −30 cm, respectively; and Figure 2b,d shows the maximum value of uplift in lobe G and subsidence in lobe H are approximately 50 cm and −50 cm, respectively.To extract the AT (along-track) deformation fields, we performed the conventional MAI method based on two interferometric pairs of ALOS-2 (only one burst will be used).The achieved accuracy of the AT deformation is approximately 8 cm with coherence of 0.6.Around the epicenter, the coherence is high enough.Furthermore, the coherence can be further enhanced by multilooking and filtering processes.It is suitable to use this method to measure the AT deformation of this earthquake.
Four SLC images, composed of both forward-and backward-looing SLC images from master and slave SAR data of ALOS-2, were generated by controlling the Doppler Centroids and processing bandwidth (PBW) for sub-aperture images according to the MAI processor [9][10][11][12].The conventional MAI processing for the SLC data includes azimuth common-band filtering between the master and slave SLC images and beam-splitting for the purpose of dividing forward-and backward-looking spectrum.During the MAI processing, we apply thirteen-looks and two-looks in the azimuth and range directions, respectively, to reduce the phase noise in the interferograms.The band-width of common band filtering is 482.39Hz and 731.23 Hz, two forward-looking and two backward-looking images were separated out from interferometric pairs of ALOS-2, firstly.Subsequently, both forwardand backward-looking interferograms were generated independently.Finally, MAI interferograms were obtained through dual-interferometry between forward-and backward-looking interferograms.
In particular, we used the estimated differential ionospheric phase to mitigate the relative azimuth shift caused by ionosphere in the MAI interferogram.The AT deformation measured by the MAI can be defined by the following equation [30]: To extract the AT (along-track) deformation fields, we performed the conventional MAI method based on two interferometric pairs of ALOS-2 (only one burst will be used).The achieved accuracy of the AT deformation is approximately 8 cm with coherence of 0.6.Around the epicenter, the coherence is high enough.Furthermore, the coherence can be further enhanced by multilooking and filtering processes.It is suitable to use this method to measure the AT deformation of this earthquake.
Four SLC images, composed of both forward-and backward-looing SLC images from master and slave SAR data of ALOS-2, were generated by controlling the Doppler Centroids and processing bandwidth (PBW) for sub-aperture images according to the MAI processor [9][10][11][12].The conventional MAI processing for the SLC data includes azimuth common-band filtering between the master and slave SLC images and beam-splitting for the purpose of dividing forward-and backward-looking spectrum.During the MAI processing, we apply thirteen-looks and two-looks in the azimuth and range directions, respectively, to reduce the phase noise in the interferograms.The band-width of common band filtering is 482.39Hz and 731.23 Hz, two forward-looking and two backward-looking images were separated out from interferometric pairs of ALOS-2, firstly.Subsequently, both forwardand backward-looking interferograms were generated independently.Finally, MAI interferograms were obtained through dual-interferometry between forward-and backward-looking interferograms.In particular, we used the estimated differential ionospheric phase to mitigate the relative azimuth shift caused by ionosphere in the MAI interferogram.The AT deformation measured by the MAI can be defined by the following equation [30]: where ∆x is the AT deformation, Φ MAI is the phase of the MAI interferogram, v g is the velocity of the beam footprint on the ground, K a is the azimuth frequency modulation (FM) rate, T C is the burst cycle length, and n is the difference of burst numbers.Figure 3 shows the AT deformation maps.(a) and (b) are the MAI results based on ascending and descending images, respectively.The displacements of ascending and descending are a little different, because the projected thrust motion vectors along the ascending and descending azimuth direction are different.Figure 3 indicates that the AT deformation is concentrated in the lobe G, and the measured maximum surface displacement is approximately 50 cm away from the SAR sensor.It reflects the earthquake has some thin features of strike slip.
Remote Sens. 2018, 10, x FOR PEER REVIEW 6 of 16 where Δx is the AT deformation, φMAI is the phase of the MAI interferogram, vg is the velocity of the beam footprint on the ground, Ka is the azimuth frequency modulation (FM) rate, TC is the burst cycle length, and n is the difference of burst numbers.Figure 3 shows the AT deformation maps.(a) and (b) are the MAI results based on ascending and descending images, respectively.The displacements of ascending and descending are a little different, because the projected thrust motion vectors along the ascending and descending azimuth direction are different.Figure 3 indicates that the AT deformation is concentrated in the lobe G, and the measured maximum surface displacement is approximately 50 cm away from the SAR sensor.It reflects the earthquake has some thin features of strike slip.
The deformation field covers an area of about 60 × 70 km 2 with two clear deformation lobes which divided by the trending bound AB mentioned in section 3.1.
The heading and length of the arrows indicate the direction and magnitude of the horizontal motion, respectively.The U-D component map (Figure 4a) shows somewhat similar displacement patterns with the LOS deformation maps (Figure 2) because vertical motion has dominant during the earthquake.Uplift appears in the lobe G, while subsidence appears in the lobe H.The maximum uplift and subsidence values are 100 cm and -50 cm, respectively, thus reaching up to a relative offset of 150 cm in the vertical direction.It can be confirmed that the U-D components are the most violent, and the westward motion plays a major role in the horizontal motion.The vectors reflect the characteristics of a pure thrust fault.The E-W component map (Figure 4b) clearly shows a pushing effect of plates.Both the maximum displacement in lobe G and H is approximately 50 cm toward the west.The earthquake demonstrated dominantly vertical and westward motion along a thrust fault, but it also encompassed movement in the lobe G of the N-S component map (Figure 4c), the maximum displacement in N-S direction is approximately 70 cm.The city of Sarpol-e-Zahab is located at the region that had mixed motion in the horizontal and vertical directions, so it suffered
The deformation field covers an area of about 60 × 70 km 2 with two clear deformation lobes which divided by the trending bound AB mentioned in Section 3.1.
The heading and length of the arrows indicate the direction and magnitude of the horizontal motion, respectively.The U-D component map (Figure 4a) shows somewhat similar displacement patterns with the LOS deformation maps (Figure 2) because vertical motion has dominant during the earthquake.Uplift appears in the lobe G, while subsidence appears in the lobe H.The maximum uplift and subsidence values are 100 cm and -50 cm, respectively, thus reaching up to a relative offset of 150 cm in the vertical direction.It can be confirmed that the U-D components are the most violent, and the westward motion plays a major role in the horizontal motion.The vectors reflect the characteristics of a pure thrust fault.The E-W component map (Figure 4b) clearly shows a pushing effect of plates.Both the maximum displacement in lobe G and H is approximately 50 cm toward the west.The earthquake demonstrated dominantly vertical and westward motion along a thrust fault, but it also encompassed movement in the lobe G of the N-S component map (Figure 4c), the maximum displacement in N-S direction is approximately 70 cm.The city of Sarpol-e-Zahab is located at the region that had mixed motion in the horizontal and vertical directions, so it suffered the highest intensity (the intensity is VIII) and destructions (about 12,000 buildings collapse).At the same time, the surface ruptures confirmed by Kobayashi [1] distributed in the uplift region mostly.the highest intensity (the intensity is VIII) and destructions (about 12,000 buildings collapse).At the same time, the surface ruptures confirmed by Kobayashi [1] distributed in the uplift region mostly.

Subsampling and Weighting of Data
The LOS and 3D data as shown in Figure 2 and Figure 4 consist of million data points.As the displacement field is varying smoothly, we can subsample the data with a factor of 4 and a quadtree algorithm to obtain reasonable number of data points and good spatial representation of the LOS and 3D displacements [34].The numbers of remaining data points for four LOS deformation fields are 1261, 959, 1552, 1184, respectively.Meanwhile the numbers of remaining data points for EW, NS, and U-D directions are 1012, 1124, and 869, respectively.Despite dramatic reduction in quantity, this down-sampled observations well retain the major features of the original interferograms.
We weight the data based on standard deviation of their deformation data [35].The standard deviation can be calculated by Equation (2), ( ) 1 Where u(i, j) is the observed displacement of one point, mu(i, j) is the mean value of displacements of all points, N is the number of points.After calculation, the weight factors for four subsampled LOS data sets are 0.31, 0.24, 0.2, and 0.25; the weight factors for EW, NS, and U-D directions are 0.45, 0.15, 0.4.

Results and Discussion
With the using of the LOS and 3D deformation data sets, we aim to demonstrate the effect of using the 3D coseismic deformation fields on fault slip inversion.For the 2017 Iraq Mw7.3 earthquake, based on the LOS and 3D deformation data sets, respectively, we first assumed simple uniform-slip fault and used a nonlinear optimization approach to find the optimum model parameters.Subsequently, we compared our results with the published parameters [1,7] to confirm the main fault.Finally, we took the obtained fault geometry and inverted the fault variable slip using a linear approach.

Nonlinear Optimization
Based on the assumption of homogeneous elastic half-space, the geodetic data can be associated

Subsampling and Weighting of Data
The LOS and 3D data as shown in Figures 2 and 4 consist of million data points.As the displacement field is varying smoothly, we can subsample the data with a factor of 4 and a quadtree algorithm to obtain reasonable number of data points and good spatial representation of the LOS and 3D displacements [34].The numbers of remaining data points for four LOS deformation fields are 1261, 959, 1552, 1184, respectively.Meanwhile the numbers of remaining data points for EW, NS, and U-D directions are 1012, 1124, and 869, respectively.Despite dramatic reduction in quantity, this down-sampled observations well retain the major features of the original interferograms.
We weight the data based on standard deviation of their deformation data [35].The standard deviation can be calculated by Equation (2), where u(i, j) is the observed displacement of one point, m u (i, j) is the mean value of displacements of all points, N is the number of points.After calculation, the weight factors for four subsampled LOS data sets are 0.31, 0.24, 0.2, and 0.25; the weight factors for EW, NS, and U-D directions are 0.45, 0.15, 0.4.

Results and Discussion
With the using of the LOS and 3D deformation data sets, we aim to demonstrate the effect of using the 3D coseismic deformation fields on fault slip inversion.For the 2017 Iraq Mw7.3 earthquake, based on the LOS and 3D deformation data sets, respectively, we first assumed simple uniform-slip fault and used a nonlinear optimization approach to find the optimum model parameters.Subsequently, we compared our results with the published parameters [1,7] to confirm the main fault.Finally, we took the obtained fault geometry and inverted the fault variable slip using a linear approach.

Nonlinear Optimization
Based on the assumption of homogeneous elastic half-space, the geodetic data can be associated with a shear dislocation model that contains a set of fault source parameters.The Okada model is widely applied into the inversion of the main fault source parameters [36][37][38][39][40][41][42][43].The fault source parameters can be derived by minimizing the discrepancies between the observed and simulated surface displacements.We set different objective functions to search optimum fault source parameters, including fault length, width, strike, dip, rake, slip, and depth.
When LOS deformation data sets were applied to inversion, the objective function can be expressed by Equation ( 3) where k is the number of data sets; i is the index, j is the number i data sets' index; α i is the weight factor of number I data sets; N i is the counts of observation of number i data sets; u i , f i (t) and σ i are observed deformation; simulated deformation and corresponding RMS the of number i data sets.
When the 3D deformation data sets were applied to inversion, the objective function can be defined as Equation ( 4) [20], where T denotes the Green's functions that describe how slip across rectangular dislocation causes displacement at the ground surface [36,37]; t represents a set of the fault geometric parameters; s represents a set of the slips of all the fault patches; G 3d (t)•s denotes the simulated deformation at a ground point; Σ uu is the variance-covariance matrix describing the quality of the highly spatially correlated observation data, which can be determined by using an exponential covariance function [20].
Based on the objective functions, we searched the optimum fault source parameters in intervals provided by the GCMT (Global Centroid Moment Tensor), using the Levenberg-Marquardt algorithm which added damping factor [30].We set the initial damping factor as 0.06.Table 2 shows the fault geometric parameters of our experiments and other studies [1,2]   It can be seen that our two results are quite similar with each other and closing to the published fault geometric parameters [1,2].The strike is approximately 353 • , close to NNW.The dip is approximately 15 • , thus the fault slip has thrust property.The rake is approximately 138 • , thus the fault slip has dextral property.

Variable Slip Inversion
Assuming uniform slip on fault planes is a simplification since the heterogeneity in fault displacement and rake is observed on all scales of fault [44].The fault slip is variable.Therefore, based on the results of the nonlinear inversion, we extend the fault length and width to 60 km and 30 km, then divide it into the multiple patches of 1.5 km by 1.5 km, thus having 800 fault patches.Based on the Okada model, the Green functions were calculated, which contains the predicted LOS and three direction changes at each location in the displacement data set produced by unit slip on each fault patch (linear inversion).Thus having Equation (4), where u is the observed deformation; G 1 and G 2 are the Green's function relating the surface deformation to the strike-slip and dip-slip, respectively; L 1 and L 2 are the Laplace smoothing constraint matrix of strike-slip and dip-slip, respectively; κ is the smoothing factor; s 1 and s 2 are strike-slip and dip-slip; ε is the error of observation and model building.
The Laplacian operator was used to smooth the slip solution, and its effect was controlled by the smoothing factor κ, which also controls a constraint of zero on the side and lower edges of the fault plane [44].We fixed the rake at the value derived from the uniform slip model, and use the Fast Non-Negative Least Square (FNNLS) algorithm [45] to solve Equation ( 4).The amount of slip S = (s 1 2 + s 2 2 ) 1/2 .Figure 5a,b shows the main fault slip distribution based on the 3D and LOS deformation data sets.Figure 5a shows that the slip distributes at depths between 12 and 18 km, the peak slip with value of 6.53 m appears at depth of 14.5 km; Figure 5b shows that the slip distributes at depths between 12 and 18 km too.After comparing, we can find that the slip distribution of Figure 5a is more concentrated, more reasonable and having finer details than Figure 5b, thus the constraining by the 3D deformation data sets is more successful than the constraining by the LOS deformation data sets.
Remote Sens. 2018, 10, x FOR PEER REVIEW 9 of 16 on the results of the nonlinear inversion, we extend the fault length and width to 60 km and 30 km, then divide it into the multiple patches of 1.5 km by 1.5 km, thus having 800 fault patches.Based on the Okada model, the Green functions were calculated, which contains the predicted LOS and three direction changes at each location in the displacement data set produced by unit slip on each fault patch (linear inversion).Thus having Equation (4), where u is the observed deformation; G1 and G2 are the Green's function relating the surface deformation to the strike-slip and dip-slip, respectively; L1 and L2 are the Laplace smoothing constraint matrix of strike-slip and dip-slip, respectively; κ is the smoothing factor; s1 and s2 are strikeslip and dip-slip; ε is the error of observation and model building.
The Laplacian operator was used to smooth the slip solution, and its effect was controlled by the smoothing factor κ, which also controls a constraint of zero on the side and lower edges of the fault plane [44].We fixed the rake at the value derived from the uniform slip model, and use the Fast Non-Negative Least Square (FNNLS) algorithm [45] to solve Equation ( 4).The amount of slip S = (s1 2 + s2 2 ) 1/2 .Figure 5a,b shows the main fault slip distribution based on the 3D and LOS deformation data sets.Figure 5a shows that the slip distributes at depths between 12 and 18 km, the peak slip with value of 6.53 m appears at depth of 14.5 km; Figure 5b shows that the slip distributes at depths between 12 and 18 km too.After comparing, we can find that the slip distribution of Figure 5a is more concentrated, more reasonable and having finer details than Figure 5b, thus the constraining by the 3D deformation data sets is more successful than the constraining by the LOS deformation data sets.

Comparison and Evaluation of Results
For further comparison and evaluation, we simulate the ascending and descending LOS deformation maps to calculate the residuals by our fault geometric parameters and slip based on Okada model.The 3D simulated deformation components were all projected to the LOS direction and the common area was cut for convenient comparison.
Figure 7 shows the measurements, models and residuals maps.The characteristics of simulated LOS deformation maps are well matched with the real LOS deformation maps, no matter based on the fault parameters and slip inversed from the LOS data sets or 3D data sets.The area where residuals exist of Figure 7c,f is smaller than the area of Figure 7i,l.The obvious residuals appear in the near field of Figure 7i,l.
For further comparison, we derive the profile CD on deformation maps.Figure 8 shows the displacements along the profile CD. Figure 8 shows the displacements along the profile CD.Regard smoothness and tendency of curve as a comparison-reference.It is clear that the profiles in Figure 8a,b are in agreement mostly, while the profiles in Figure 8c,d are in disagreement relatively.The differences exist due to the fault parameters inversed based on the 3D and LOS coseismic deformation data sets Figure 9 shows the residual distribution histograms.According to the Figure 9, the standard deviation (Std) of residuals in Figure 7c,f are 2.4 cm and 2.9 cm, respectively.The standard deviation of residuals in Figure 7i,l are 3.5 cm and 4.2 cm, respectively.The morphological characteristics of curves in Figure 9a,b are slimmer and higher than the curves in Figure 9c,d.
The comparative analysis indicated that the residuals of inversion based on the 3D data sets are controlled better.We regard the Std as the comparative standard.Thus, the fault inversion using the 3D data sets improved nearly 30% relatively the use of the LOS data sets.These results are within a reasonable range as defined by other studies [1,7], but our analysis is unique due to the reconstruction of constrained 3D surface movements and the derivation of the fault's optimal focal plane.These results support existing seismological interpretations and contribute to a better understanding of the main fault's behavior.

Comparison and Evaluation of Results
For further comparison and evaluation, we simulate the ascending and descending LOS deformation maps to calculate the residuals by our fault geometric parameters and slip based on Okada model.The 3D simulated deformation components were all projected to the LOS direction and the common area was cut for convenient comparison.
Figure 7 shows the measurements, models and residuals maps.The characteristics of simulated LOS deformation maps are well matched with the real LOS deformation maps, no matter based on the fault parameters and slip inversed from the LOS data sets or 3D data sets.The area where residuals exist of Figure 7c,f is smaller than the area of Figure 7i,l.The obvious residuals appear in the near field of Figure 7i,l.
For further comparison, we derive the profile CD on deformation maps.Figure 8 shows the displacements along the profile CD. Figure 8 shows the displacements along the profile CD.Regard smoothness and tendency of curve as a comparison-reference.It is clear that the profiles in Figure 8a,b are in agreement mostly, while the profiles in Figure 8c,d are in disagreement relatively.The differences exist due to the fault parameters inversed based on the 3D and LOS coseismic deformation data sets.
Figure 9 shows the residual distribution histograms.According to the Figure 9, the standard deviation (Std) of residuals in Figure 7c,f are 2.4 cm and 2.9 cm, respectively.The standard deviation of residuals in Figure 7i,l are 3.5 cm and 4.2 cm, respectively.The morphological characteristics of curves in Figure 9a,b are slimmer and higher than the curves in Figure 9c,d.

Conclusion
This paper combined DInSAR and MAI observations to reconstruct the 3D coseismic deformation associated with the 2017 Iraq Mw7.3 earthquake.From the 3D coseismic deformation filed, it can be found that the maximum coseismic displacements in U-D, N-S, and E-W direction are 100 cm, 50 cm, and 100 cm, respectively.The U-D components are the most violent, and the westward motion plays a major role in the horizontal motion.With the 3D displacements derived from the DInSAR and MAI measurements, it can be concluded that the Iraq earthquake predominantly resulted in reverse slip with small dextral slip component.Furthermore, fault inversions were performed based on the 3D and LOS coseismic deformation data sets, respectively.After comparison, we find that the 3D coseismic deformation data sets improved the inversions' accuracy by 30%.Our results reveal a deep dislocation on an NNW trending fault (the strike is 352.63°)extending about 60 km, along the fault dips 14.76° to the ENE.The estimated seismic moment is 8.44 × 10 19 Nm (Mw7.3), which is close to the solution provided by USGS (United States Geological Survey).The slip distributed at the depth between 12 and 18 km, and the peak slip of 6.53 appears at the depth of 14.5 km left near the epicenter.Further comparative analysis indicated that the inferred geodetic results derived from this study can be considered as the optimized results due the use of 3D displacements.
Considering the geological structure in the earthquake region, locations of aftershocks and fault source-parameters, it comes to a preliminary conclusion that the ZMFF (the Zagros Mountain Front fault) is responsible for the earthquake.In total, this paper demonstrates the superiority of using the 3D coseismic deformation fields on fault inversion.The comparative analysis indicated that the residuals of inversion based on the 3D data sets are controlled better.We regard the Std as the comparative standard.Thus, the fault inversion using the 3D data sets improved nearly 30% relatively the use of the LOS data sets.These results are within a reasonable range as defined by other studies [1,7], but our analysis is unique due to the reconstruction of constrained 3D surface movements and the derivation of the fault's optimal focal plane.These results support existing seismological interpretations and contribute to a better understanding of the main fault's behavior.

Conclusions
This paper combined DInSAR and MAI observations to reconstruct the 3D coseismic deformation associated with the 2017 Iraq Mw7.3 earthquake.From the 3D coseismic deformation filed, it can be found that the maximum coseismic displacements in U-D, N-S, and E-W direction are 100 cm, 50 cm, and 100 cm, respectively.The U-D components are the most violent, and the westward motion plays a major role in the horizontal motion.With the 3D displacements derived from the DInSAR and MAI measurements, it can be concluded that the Iraq earthquake predominantly resulted in reverse slip with small dextral slip component.Furthermore, fault inversions were performed based on the 3D and LOS coseismic deformation data sets, respectively.After comparison, we find that the 3D coseismic deformation data sets improved the inversions' accuracy by 30%.Our results reveal a deep dislocation on an NNW trending fault (the strike is 352.63 • ) extending about 60 km, along the fault dips 14.76 • to the ENE.The estimated seismic moment is 8.44 × 10 19 Nm (Mw7.3), which is close to the solution provided by USGS (United States Geological Survey).The slip distributed at the depth between 12 and 18 km, and the peak slip of 6.53 appears at the depth of 14.5 km left near the epicenter.Further comparative analysis indicated that the inferred geodetic results derived from this study can be considered as the optimized results due the use of 3D displacements.
Considering the geological structure in the earthquake region, locations of aftershocks and fault source-parameters, it comes to a preliminary conclusion that the ZMFF (the Zagros Mountain Front fault) is responsible for the earthquake.In total, this paper demonstrates the superiority of using the 3D coseismic deformation fields on fault inversion.

Figure 1 .
Figure 1.Shaded relief map of study area and image coverage.The yellow star indicates the epicenter of the mainshock (Mw 7.3).The red dots show the aftershocks with magnitude higher than Mw4.0 before 17 November 2017.The blue solid squares indicate the five cities with intensity reaching up to VI or higher.The black and red rectangles mark the coverage of SAR images.The dark red line represents the boundary between the Iraq and Iran.

Figure 1 .
Figure 1.Shaded relief map of study area and image coverage.The yellow star indicates the epicenter of the mainshock (Mw 7.3).The red dots show the aftershocks with magnitude higher than Mw4.0 before 17 November 2017.The blue solid squares indicate the five cities with intensity reaching up to VI or higher.The black and red rectangles mark the coverage of SAR images.The dark red line represents the boundary between the Iraq and Iran.

Figure 2 .
Figure 2. The coseismic deformation maps in the LOS directions.(a) and (b) obtained from ALOS-2 DInSAR processing.(c) and (d) obtained from Sentinel-1a DInSAR processing.The red star represents the epicenter of the mainshock; The red solid squares represent the city Halabja and Sarpol-e-Zahab.The black line AB is the deformation bound.

Figure 2 .
Figure 2. The coseismic deformation maps in the LOS directions.(a,b) obtained from ALOS-2 DInSAR processing.(c,d) obtained from Sentinel-1a DInSAR processing.The red star represents the epicenter of the mainshock; The red solid squares represent the city Halabja and Sarpol-e-Zahab.The black line AB is the deformation bound.

Figure 3 .
Figure 3.The along track deformation maps derived from MAI (multiple aperture InSAR) processing.(a) and (b) derived from ascending and descending ALOS-2 MAI processing, respectively.The red star represents the epicenter of the mainshock; The red solid squares represent the city Halabja and Sarpol-e-Zahab.The black line AB is the deformation bound.

Figure 3 .
Figure 3.The along track deformation maps derived from MAI (multiple aperture InSAR) processing.(a,b) derived from ascending and descending ALOS-2 MAI processing, respectively.The red star represents the epicenter of the mainshock; The red solid squares represent the city Halabja and Sarpol-e-Zahab.The black line AB is the deformation bound.

Figure 4 .
Figure 4.The 3D coseismic displacements maps.(a), (b), and (c) represent the deformation maps in up-down, east-west, and north-south directions.The red star represents the epicenter of the mainshock.The red solid squares represent the city Halabja and Sarpol-e-Zahab.The black line AB is the deformation bound in Figure 2. The arrows in (a) are the horizontal motion vectors including EW and N-S displacements vectors.

Figure 4 .
Figure 4.The 3D coseismic displacements maps.(a-c) represent the deformation maps in up-down, east-west, and north-south directions.The red star represents the epicenter of the mainshock.The red solid squares represent the city Halabja and Sarpol-e-Zahab.The black line AB is the deformation bound in Figure 2. The arrows in (a) are the horizontal motion vectors including EW and N-S displacements vectors.

a
Denotes the "Experiment A" and "Experiment B" represents the results based on the LOS and 3D deformation data sets, respectively.

Figure 5 .
Figure 5.The main fault slip distribution.(a) and (b) are obtained by variable slip inversion based on the 3D deformation data sets and the line-of-sight (LOS) deformation data sets, respectively.

Figure 6
Figure 6 indicates that fault has weak dextral slip motion and strongly thrust motion.Most of the slip distribute in south area of epicenter.The depth of the fault top-left corner is 10.8 km.The city Sarpol-e-Zahab locates at the fault boundary almost, that's why it suffered the highest intensity and destructions.The geologic structure background mentioned in Section 2 accounts for these characteristics of fault slip, precisely.

Figure 5 .
Figure 5.The main fault slip distribution.(a,b) are obtained by variable slip inversion based on the 3D deformation data sets and the line-of-sight (LOS) deformation data sets, respectively.

Figure 6 16 Figure 6 .
Figure 6 indicates that fault has weak dextral slip motion and strongly thrust motion.Most of the slip distribute in south area of epicenter.The depth of the fault top-left corner is 10.8 km.The city Sarpol-e-Zahab locates at the fault boundary almost, that's why it suffered the highest intensity

Figure 6 .
Figure 6.Schematic view of the source mechanism.(a,b) are source mechanism based on the 3D deformation data sets and the line-of-sight (LOS) deformation data sets, respectively.The blue star represents the epicenter of the mainshock.The red solid squares represent the city Halabja and Sarpol-e-Zahab.The arrows in (a,b) are the slip vectors.

Figure 7 .
Figure 7.The measurements, models and residuals maps.(a)-(f) derived from the 3D coseismic deformation data sets.(g)-(l) derived from the LOS coseismic deformation data sets.The black line is the profile CD.

Figure 7 .
Figure 7.The measurements, models and residuals maps.(a-f) derived from the 3D coseismic deformation data sets.(g-l) derived from the LOS coseismic deformation data sets.The black line is the profile CD.

Figure 8 .
Figure 8.The displacements along the profile CD.(a) and (b) correspond the profile CD in Figure 8(a,b) and Figure 8(d, e), respectively.(c) and (d) correspond the profile CD in Figure 8(g,h) and Figure 8(j,k), respectively.

Figure 8 .
Figure 8.The displacements along the profile CD. (a,b) correspond the profile CD in Figure 7a,b and Figure 7d,e, respectively.(c,d) correspond the profile CD in Figure 7g,h and Figure 7j,k, respectively.

Figure 9 .
Figure 9.The residual distribution histograms.(a) and (b) correspond the residual in Figure 7c,f, respectively.(c) and (d) correspond the resiudal in Figure 7i,l, respectively.

Table 1 .
The parameters of selected SAR images.

Table 2 .
Fault geometric parameters from nonlinear optimizations assuming uniform slip (this study and earlier publications).