Estimation of north Tabriz fault parameters using neural networks and 3D tropospherically corrected surface displacement field

ABSTRACT In this paper, parameters of north Tabriz fault are studied using 3D displacement field and artificial neural networks (ANNs). We provide the 3D surface displacement along the north Tabriz fault using an integration of tropospherically corrected InSAR, GPS and precise levelling data. To perform the InSAR analysis, we use the 17 ENVISAT radar acquisitions. The line of sight (LOS) displacement field was corrected using the ERA-Interim global meteorological reanalysis models. In order to calculate the 3D displacement, we use a Simultaneous and Integrated Strain Tensor Estimation from Geodetic and Satellite Deformation Measurements approach. ANNs used to estimate the parameters. 3D displacement field and ANN algorithm yields an average slip rate of 6.1 ± 0.01 mm/year with a locking depth of 13.4 ± 0.01 km. This rate is consistent with previous geodetic estimates based on recent global positioning system measurements and InSAR analysis. In addition, the length, width and dip angle of fault are about 101.2 ± 0.01 km, 25.06 ± 0.03 km and 63 ± 0.05 deg. This study demonstrates that interseismic displacement with a sub-centimetre rate can be successfully detected by integration of multi temporal InSAR techniques, GPS and precise leveling data and shows that ANN algorithm is suitable for estimating the parameters of north Tabriz fault.


Introduction
When a moderate to strong displacement occurs, we can apply interferometric synthetic aperture radar (InSAR) technique to compute a differential interferogram. The displacement field can be more precise in three dimensions with the help of GPS and precise leveling data. If there is a fault in the region, each differential interferogram contains the information concerning the characteristics of the fault (Lundgren and Stramondo 2002). To retrieve the geometric parameters of the fault plane, we compute a huge number of synthetic interferogram, also referred to as models.
The inversion procedure stems from an artificial neural network (ANN) properly generated and trained to recognize the synthetic interferograms and to classify them as belonging to a certain subset. Furthermore, the already trained neural network is used to work with the interseismic InSAR measures. First, it singles out the synthetic interferogram best fitting the measured one, and then it is used to retrieve the fault plane parameters by this latter. More, in particular, the solution of the inverse problem means to recover the source parameters from the knowledge of InSAR surface displacement field.
To this aim, some significant results have already been achieved by means of the simulated annealing technique (Lundgren and Stramondo 2002). Inversion algorithms using geodetic data and based on the maximum likelihood criteria (Fukahata et al. 2004) or a quasi Montecarlo (Funning et al. 2005) iterative approach have been also proposed. However, due to the intrinsic ill-posedeness of the problem, some issues remain still open and more investigation is needed. ANNs have already been recognized as being a powerful tool for inversion procedure in geodesy and remote sensing applications (Del Frate et al. 2003). Other applications fields are include the urban changes (Del Frate et al. 2003), environmental studies concerning oil spill (Del Frate, Salvatori 2004), atmospheric and meteorological models (Jolivet et al. 2012, Li et al. 2009) and ozone retrieval. They are composed of many nonlinear computational elements (called neurons) operating in parallel and connected by the so-called synapses. The use of neural networks is often effective because they can simultaneously address nonlinear dependencies and complex physical behaviour.
We combine 17 ENVISAT acquisition of the area in the north-west (NW) of Iran that north Tabriz fault is located in this region. After correcting the effect of the troposphere on the interferograms using the global meteorological models, we calculate the interseismic displacement in the satellite line of sight (LOS) between the years 2003 and 2010. By integrating the displacement field and GPS observations and precise leveling data, we obtain 3D displacement field. Then, using the ANN and Levenberg-Marquardt algorithms, we estimate the parameters of the north Tabriz fault and statistical parameters. Finally, we compare the results with previous studies that have been done with other methods.

Tectonic setting and historical seismicity of the study area
The Tabriz area is part of the complex tectonic system due to the interaction between Arabia, Anatolia and Eurasia and comprising the north Anatolian fault, the east Anatolian fault and the Caucasus mountains which bounds the Zagros mountains (Jackson 1992). The north Tabriz fault is a major brittle structure which is located in the NW part of Azerbaijan. More locally, the north Tabriz fault is a clear trending strike-slip fault that runs for more than 100 km (Karakhanian et al. 2004). To the west, the north Tabriz fault joins the east-west trending right-lateral strike-slip Tasuj faults and the north-dipping Sufian reverse fault, which bounds the lake Urumieh to the north. To the east, it merges with the north Bozghush fault and the South Bozghush north-dipping reverse fault located on both sides of the Bozghush Range (Karakhanian et al. 2004). The north Tabriz fault is composed of several en-echelon right-stepping segments associated with pull-apart basins, the most important being the Tabriz basin where the city of Tabriz is located. Historical documents going back two millennia and instrumental records demonstrate that NW Iran and Eastern Turkey have been struck by numerous destructive earthquakes, the most recent largest one being the Mw 7.1, 23 October 2011 Van earthquake associated with reverse slip on a north-east-south-west trending fault (Akoglu et al. 2012). Historical documents suggest that the last two large earthquakes on the north Tabriz fault took place within 60 years in the 18th century on its adjacent segments (Ambraseys and Melville 1982). Seismotectonics map of NW Iran is shown in Figure 1.

Calculate the 3D earth surface displacement field
In order to achieve the LOS displacement field using InSAR, we must produce interferograms. Tropospheric artefacts are the main error source in interferograms (Zebker et al. 1997). One of the best ways to reduce this effect is the use of meteorological data. To minimize other errors, including Digital Elevation Model errors, orbital residual and poor coherence, a relatively new and advanced technique, i.e. multi-temporal InSAR technique has been developed. This technique uses multi temporal stacks of SAR acquisitions to generate time series of ground deformations for individual targets (Berardino et al. 2002;Ferretti et al. 2001). One of the best methods to investigate the displacement field is the Stanford method technique (StaMPS/MTI) that uses both Persistent Scatterer Interferometry (PSI) and Small-Baseline Subset (SBAS) InSAR techniques (Hooper 2008). This method combines both sets of results (PSI and SBAS) to improve phase unwrapping and the spatial sampling of the signal of interest (Sousa et al. 2011). This technique recognizes potential persistent scatterer candidates using Amplitude Dispersion Index (ADI). In addition to persistent scatterers, StaMPS also detects and uses slowly decorrelating filtered phase pixels (SDFP). The number of stable points (SP D PS C SDFP) candidates increases with increasing ADI.
Without extra constrains and/or assumptions regarding surface displacements, a single InSAR measurement generally does not allow to derive the 3D surface displacement (Li et al. 2015). In order to calculate the 3D displacement, we use an integration method (Pietrantonio and Reguzzi 2004). GPS, precise leveling and InSAR data integration is a method of calculating the 3D displacement. One of the methods for integrating GPS, precise leveling and InSAR data is based on the weighted least square (WLS) and in addition to obtaining the 3D displacement, calculate the strain tensor elements. This method is called Simultaneous and Integrated Strain Tensor Estimation from Geodetic and Satellite Deformation Measurements (SISTEM). In the following, we introduce the SISTEM approach. For this purpose, we consider the following parameters: P: an arbitrary point X 0 (x 10 , x 20 , x 30 ): position of P U i : displacements of P, i D 1, 2, 3 N: number of surrounding experimental points (EPs) X n (x 1n , x 2n , x 3n ): positions of EPs, n D 1 … N U n (u 1n , u 2n , u 3n ): displacements of EPs Under such a hypothesis, adopting a linear approach, the problem of estimating the displacement components U i of the point P from the experimental data U n can be modelled by the N equations Figure 1. Morphotectonic map of the study area within the Turkish-Iranian plateau showing active faults (black lines), focal mechanisms and GPS velocity fields relative to Eurasia with 95% confidence ellipses on shaded relief from SRTM-3 arc second data. Continuous GPS stations are indicated with black squares at vector tails. Black-dashed square shows the frame of the SAR images studied, covering major part of the Tabriz basin and the NTF (Karimzadeh et al., 2013). (Teza et al. 2008): The matrix H can be broken down into symmetric and antisymmetric parts as E is the strain tensor(symmetric part) and V is the rotation tensor (antisymmetric part) andis the tensor product. The 3n equations of type (1) give a linear matrix equation of the following form: where A: coefficient matrix x: vector of unknown parameters ([U 1 U 2 U 3 e 11 e 12 e 13 e 22 e 23 e 33 v 1 v 2 v 3 ] T ) l: observations vector ([u 11 u 21 u 31 … u 1n u 2n u 3n ] T ) A suitable method to solve the system is the WLS which gives (7) as a suitable formula to estimate the unknown parameters: P is the inverse of the data covariance matrix Cl. We use the matrix Cll which is a weighted version of the matrix Cl (Pietrantonio and Reguzzi 2004): where d (n) is the distance between the nth EP and the arbitrary point P and d 0 is a distance decaying constant defining the level of locality of the estimation; hereafter, the parameter d 0 is defined as locality. InSAR can be related to the components of the displacement vector of an arbitrary point P: where D LOS(P) : LOS displacements at the point P U i : unknown displacement vector components, i D 1, 2, 3 [S Px sp , S Py sp , S pz sp ]: unit vector pointing from the point P toward the satellite For integrating the GPS, precise leveling and InSAR data we assume the following parameters in equation (7): (13) In order to have all the data-sets (GPS, leveling and InSAR) temporally consistent, we scaled the GPS values observed over a 10-month period to a five-month one, corresponding to the leveling and InSAR ENVISAT time coverage, under the assumption of a linear evolution of the ground deformation pattern through time (Guglielmino et al. 2011).

Geophysical inverse problem and artificial neural network
Geophysical inverse problem is to obtain information about the parameters of a fault, using observations of deformation associated with the fault (Zhdanov 2002, Tarantola 2005. The inverse problem can be summarized in the following equation: where G: mathematical relationship between the observation and characterization of earth m: internal characteristics of the land d: vector of observations of the earth's surface as a function of m The purpose of the inverse problem is to find the best value for m Figure 2. In recent years, solving of geophysical inverse problem is one of the most important issues. The method of solving linear inverse problems have been well studied (Zhdanov 2002), but in many geophysical problems we are faced with complex nonlinear problems that are not directly solvable Figure 2. Geophysical inverse problem definition. (Tarantola 2005). Therefore, we use optimization methods to solve them. In this paper, we used a neural network as a method of optimization, in order to solve the nonlinear geophysical inverse problem and estimating the parameters of the fault. During the operation optimization, we use the 3D displacement field obtained as real and accurate data.
The Stuttgart Neural Network Simulator (SNNS), a very powerful piece of neural network software, was used for estimating the fault parameters. SNNS has been developed at the institute for parallel and distributed high performance systems at the university of Stuttgart, Germany, since 1989 (Zell et al. 1995). It is available from a number of different sources, or as part of various LINUX distributions. For the estimation of parameters a procedure was applied as described in the flowchart of Figure 3. The performance evaluation of SNNS learning is measured by means of statistic parameters. The more common is the sum squared error (SSE) of the learning function. It can be written as follows: where t pj is the desired output of neuron j for the pattern p, o pj is the actual output. The mean squared error (MSE) is where n is the number of observations and p the number of parameters, or weights, to be estimated. The SSE/o-units is the SSE divided by the number of outputs. In order to further study the statistical parameters of the results, we use the Levenberg-Marquardt algorithm. This algorithm is one of the best methods for statistical analysis. The Levenberg-Marquardt curve-fitting method is actually a combination of two minimization methods: the gradient descent method and the Gauss-Newton method (Lourakis 2005). The Levenberg-Marquardt method acts more like a gradient-descent method when the parameters are far from their optimal value, and act more like the Gauss-Newton method when the parameters are close to their optimal value (Lourakis 2005).  (Zell et al. 1995).
In fitting a function y^(t, p) of an independent variable t and a vector of n parameters p to a set of m data points (t i , y i ), it is customary and convenient to minimize the sum of the weighted squares of the errors (or weighted residuals) between the measured data y(t i ) and the curve-fit function y^(t i , p). This scalar-valued goodness-of-fit measure is called the chi-squared error criterion: The value w i is a measure of the error in measurement y(t i ). If the function y^ðt i ; pÞ is nonlinear in the model parameters p, then the minimization of x 2 with respect to the parameters must be carried out iteratively. The goal of each iteration is to find a perturbation h to the parameters p that reduces x 2 (Lourakis 2005).

Data-set
To calculate the displacement field in the satellite LOS we used 17 ENVISAT acquisitions between 2003 and 2010. Features acquisitions can be seen in Table 1. The temporal and perpendicular baselines for the acquisitions can be seen in Figure 4.
We used the global meteorological models to correct troposphere effect on interferograms. European Center for Medium Range Weather Forecasting is currently producing ERA-Interim, a global reanalysis of the data rich period since 1989. This reanalysis provides values of several meteorological parameters on a global »75 km day, at 0 am, 6 am, 12 pm and 6 pm UT. The vertical stratification is described on 37 pressure levels, densely spaced at low elevation (interval of 25 hPa), with the highest level around 50 km (1 hPa). For each acquisition date, we select the ERAÀI output that is the closest to the SAR acquisition time. We interpolate the temperature, water vapor and dry air partial pressure provided at each pressure level to predict the delay. A kriging interpolation in the horizontal  dimensions and a spline interpolation along altitude is then applied to produce a map of the predicted delay. In order to achieve 3D displacement field we used GPS and precise leveling data. We used 13 new survey sites installed on two profiles in the region with an average spacing on the order of 10 km. Since our sites are only survey sites, we use the values determined (Djamour et al. 2011) for the continuous GPS in NW of Iran. Finally, velocities and their 1s confidence uncertainties were estimated in ITRF2008 and then the Eurasian reference frame was defined by minimizing the horizontal velocities of IGS stations located in Europe and central Asia (Djamour et al. 2011). Despite the high accuracy in the calculation of horizontal displacement, GPS can't measure vertical displacement with high accuracy. Therefore, in order to increase the accuracy of the height of the 3D displacement field, we used precise leveling data. Precise leveling data that were used in the calculations were the results of observations at first grade stations. Distance between stations is about 0.5-5 km. The leveling surveys have been made to specifications required for first-order precise leveling. The expected standard error for the change in elevation difference between two benchmarks as measured in two surveys can be expressed as §aL 1/2 , where a D 1 mm/km 1/2 and L is the distance in kilometres along the leveling line between benchmarks. Vertical displacement rates at leveling benchmarks have been calculated from kinematic network adjustment with an accuracy better than 0.3 mm/a for 90% of the rates. About 10% of the vertical rates have been classified as outliers since their value significantly differs from estimates at adjacent benchmarks. The GPS velocities of our GPS and precise leveling sites and their uncertainties are given in Table 2. Location of GPS and precise leveling stations can be seen in Figure 5.

Results of the processing of InSAR and integration data
Considering the characteristics of acquisitions we made 30 interferograms that have at least a 0.4 correlation. Then, using meteorological data which were introduced in section 5, effect due to the troposphere was corrected on all interferograms. After doing all this processing we used (StaMPS/ MTI) technique to calculate the LOS displacement field. Figure 6 shows the mean LOS velocity field obtained from MTI time series.   As was seen in Figure 6, we used three normal profiles (two profile in the left block and one profile in the right block because the density of observations on the right side of the displacement field is less) to better illustrate the ground motion across the fault and to compare, GPS velocities projected along the satellite LOS (Figure 7).
Results agree quite well between distances of 15 and 35 km along the profile P1. Further north, the very low density of observations does not permit determination of any other additional deformation that occurs there. The middle profile (profile P2), similarly to P1, reveals a gradual LOS velocity change of about 2 mm/year between 20 and 40 km. In the area of the right profile the density of the observations is less therefore we considered the profile in area with higher density. The velocity changes 4 mm/year near the fault.
In the top of the profile (after 35 km) the rate of the displacement is about 4 mm/year. The obtained results show that the LOS velocity field is close to GPS reports.
In next step, we obtain 3D displacement field using SISTEM method, LOS displacement field, GPS and precise leveling data. The displacement along the east-west, north-south and vertical direction, obtained by integration data can be seen in Figures 8(a), (b) and (c).
As seen in Figure 8, maximum displacement along the east-west appears in eastern areas. Also maximum displacement along the north-south is seen in southern and eastern areas. The vertical displacement obtained indicates that maximum vertical displacement existed in western areas. This shows the subsidence in the area. This subsidence can be caused by a north Tabriz fault and Urmia lake in the area.

Estimating fault parameters
Fault parameters that are considered in this paper are length, width, dip angle, locking depth and slip rate of the north Tabriz fault. Using simulated data is required in order to achieve educational settings and evaluating the performance of the trained neural network. Samples are used for training and evaluation has been prepared by changing parameters in the limit. Parameter ranges and rates of change can be seen in Table 3.
We considered a fixed network for all samples which consists of 5625 points. Given the parameters of the fault and using the model, the displacement field calculated in three dimensions, and used as inputs to the ANN. Therefore, ANN has 5625 inputs and 5 outputs neurons. Therefore, 3920 samples were prepared. The 2520 sample used for training ANN and 1400 samples were used to evaluate network. MSE is used to evaluate the network performance and evaluation ends when the MSE tends to a constant value. Figure 9 shows the network used in the software.
The neural network has been trained in 600 epochs and network training ends when the estimation error becomes a constant value. Diagrams relating to the training and evaluation can be seen in Figure 10. The vertical axis represents the MSE and the horizontal axis represents the epochs.
We've entered the 3D displacement field caused by faults in the network and extracted its parameters. The results and MSE fault parameters are shown in Table 4.  Figure 9. The neural network used in this article, in order to extract the parameters of the fault. As mentioned above, we used the Levenberg-Marquardt algorithm for the statistical study. In the first step of the algorithm, the initial value for each parameter is required. We considered the initial values according to Table 5. Figure 11(a) shows the parameters in each iteration of the algorithm and table 5 shows the standard deviation of each parameter. Finally, histogram of residual displacement of the algorithm is shown in Figure 11(b). Little difference in the final results obtained from Marquardt algorithm is indicative of the high efficiency of this algorithm.
We computed the fault parameters and relevant statistical parameters using the above methods and data. Then, we compared the displacement field obtained from data with displacement field obtained from the model along the profiles that are considered. The results of this   comparison can be seen in Figure 12. Results agree quite well between distances of 18 and 30 km along the profile P1, distances of 10 and 35 km along the profile P2 and distances of 25 and 40 km along the profile P3. Table 6 compares results of this paper with the result of previous studies. In this table the names of researchers and their research results on the north Tabriz fault is expressed.
Compare the results with the results obtained in this study indicate that these values are close to each other. In conclusion we can say that the 3D displacement field obtained by integration atmospherically corrected InSAR, GPS and precise leveling data and ANN are another useful tools for the study of north Tabriz fault.

Discussion and conclusion
In this study, we used the SISTEM approach to calculate the 3D displacement field of the area in the NW of Iran along the north Tabriz fault. Input data for the SISTEM approach were displacement field obtained from InSAR, GPS and precise leveling data. To perform the InSAR analysis, we used the 17 ENVISAT radar acquisitions acquired between 2003 and 2010. We've reduced tropospheric effect on the interferograms using ERA-Interim global meteorological reanalysis models and then we used the StaMPS/MTI technique. The purpose of this study was to estimate the parameters of north Tabriz fault and evaluate the ability of ANN in this calculation. ANNs offer a number of advantages, including requiring less formal statistical training, ability to implicitly detect complex nonlinear relationships between dependent and independent variables, ability to detect all possible interactions between predictor variables, and the availability of multiple training algorithms. Due to these advantages, this geodynamic study was done based on ANN. The software that we use to run neural network was SNNS software. The goal of the SNNS is to create an efficient and flexible simulation environment for research on, and application of neural network. Processing show the fault parameters as follow: slip rate of 6.1 mm/yr, locking depth of 13.4 km, length of 101.2 km, width of 25.06 km and dip angle of 63 deg. These values are consistent with previous studies based on recent GPS and InSAR analysis. In order to perform an accurate statistical analysis we used the Levenberg-Marquardt algorithm. After using this algorithm, standard deviations were obtained as follows: standard deviation of the slip rate, locking depth, length, width and dip angle was §0.01 mm/yr, §0.01 km, §0.01 km, §0.03 km and §0.05 deg. Due to this result and comparing these results with previous studies, ANNs can be regarded as an appropriate method in order geodynamic studies. Also, the 3D displacement field obtained by the above methods was able to detect the signals of displacement in the region. Consequently, these methods can be used to study other parts of north Tabriz fault or other faults.