ESTIMATING THE EARTHQUAKE SOURCE PARAMETERS: SIMULATED ANNEALING VERSUS NELDER-MEAD SIMPLEX ALGORITHM

. The parameter estimation is an important process in many earth science problems in order to de(cid:133)ne the features of ground movement. It is clear that the model nonlinearity makes the estimation of parameters more di¢ cult and more challenging. In this case, metaheuristic algorithms and derivative free optimization methods are more proper than classical optimization meth-ods. In this study, Simulated Annealing (SA), a well known metaheuristic algorithm, and Nelder-Mead simplex algorithm, a derivative free optimization method, are used to estimate the earthquake source parameters. The algorithms are applied on a synthetically generated data set. The estimated parameter results show that the SA is better than Nelder-Mead simplex method.


Introduction
An earthquake is de…ned by source parameters.Modeling of the source parameters requires wide knowledge of elastic half space theory due to the elastic structure of the crust.Many past attempts to infer source geometry from deformation …elds have used elasticity theory to …nd geologically plausible models that …t the major features of the observed deformation …eld.The fault plane is used to determine surface displacements on de…nite coordinates.The correct estimation of source parameters gives good predictions for the next earthquake occurance.In this case, the estimation of the fault plane parameters becomes important as modeling.Many new methodologies have been developed in the …eld of estimating the source parameters [9,11,20,22].Derivative based searching algorithms; e.g., a Quasi-Newton method [6], derivative free optimization algorithms; e.g., Nelder-Mead simplex algorithm [5,15], heuristic optimization methods; e.g., SA algorithm [2,3,7,21] and genetic algorithms [8,18], Monte-Carlo methods [26] are used for the estimation of crust model parameters.
In this study, SA algorithm is used to estimate the earthquake source parameters using analitical solution for deformation.The SA algorithm is a random search method based on the analogy between cooling and the freezing of the metals and some liquid materials at crystal structure with minimum energy and researching the minimum of a general system.The SA is based upon that of [24], which was originally proposed as a means of …nding the equilibrium con…gurations of a collection of atoms at a given temperature.Kirkpatrick [19] was …rst proposed the algorithm basis of an optimization technique for combinatorial problems.The SA permits the e¢ cient period for the dispersion of the molecules again by decreasing the temperature gradually.Thus, the method can be considered as a minimization algorithm based on natural events and nature.Also the SA is a stochastic computational method with quick approximations to global solutions for hard optimization problems in di¤erent …elds [10,13,28,16,17,27,30,31]. The SA algorithm has advantages in dealing with the strong nonlinearities and discontinuities in the hilly structure of cost function.This paper describes the SA algorithm and explores its ability to optimize NP-hard functions in geosciences through an e¢ cient exploration of the parameter space.And also the paper gives a comparison for parameter estimation with a derivative free optimization method, called Nelder-Mead simplex method [25].The synthetically generated data set is used for the application.All the calculations are done by using Matlab.
The rest of the paper is organized as follows: In Section 2 fault plane geometry is showed and inverse problem formulation is presented.In the next section the SA algorithm is brie ‡y reviewed and the implementation of SA algorithm to fault plane parameter estimation procedure is given.This is followed by simulation results that illustrates the both proposed approaches via the Monte Carlo simulation.In the last part, numerical results of the estimations are compared and the performance of the optimization algorithms is analyzed.

Fault Plane Model and Inverse Problem Formulation
2.1.The fault plane geometry.The fault plane parameters play an important role for de…ning the characteristics of surface displacements.It is a hard work to estimate earthquake source parameters and it requires a well de…ned model.The most commonly used crustal model is the homogeneous, isotropic, linear, and elastic half-space.In spite of its limitations, the elastic half space model is widely used because of the simplicity of the expressions [7]. Figure 1 illustrates the fault plane geometry in three dimensions.In this study, a plate model, which is expressed by Aktu¼ g [4], is considered because of its simplicity and speed in computing.In the model surface relocations, caused by the slip vector, is de…ned as analitical equations for quadrangular area.The slip is consist of two components such as "lateral range"in the direction of the fault plane and the "vertical range" as vertical to lateral range.The both range functions are given below: Lateral range : Vertical range : + sin tan 1 qR I 5 sin cos where the I 1 ; I 2 ; I 3 ; I 4 ; I 5 are analitical equations and contain source parameters.These equations are consist of nonlinear combinations of the fault plane parameters.The reader is referred to the works by [4] for a detailed description of parameters.The equation systems, given above with lateral and vertical ranges, show the functional relations between the fault geometry and surface displacements.Because of the nonlinearity of these complex analitical functions, classical parameter estimation methods may fail.In this case, estimation of fault plane parameters by the SA algorithm makes certain facilities to get global solutions; e.g.providing su¢cient convergence, escaping local traps, easy implementation to real world NP hard optimization problems.

2.2.
Inverse problem formulation.One of the main aim of geophysical inversion is to identify the all models which give an acceptable loss between predicted and observed data.In terms of the loss function, a measure of the di¤erence between observed and synthetic data that varies as a function of solution parameters [23].
In this study, the inversion problem, which describes the geometry and the slip of the fault plane, is formed with analitical lateral and vertical range equations.
Let p be the nine dimensional vector whose elements p j are the geometrical source parameters wanted to be estimated denoted by g 1 g 2 g N be the N -dimensional vector whose ith element g i , i = 1; 2; :::; N is the surface relocation vector with three dimensions (x; y; z) at the ith observation point.Therefore, the observed surface displacements u can be written as a function g of the fault model parameters p where e is a vector of observational errors.The surface relocation vector u consists of two components as in the lateral direction, u l = u l x ; u l y ; u l z and in the vertical direction, u v = u v x ; u v y ; u v z .Thus, the equation 2.1 can be written as 1; 2; :::; 50 , j = x; y; z. (2. 2) The surface displacements are related nonlinearly to the fault plane parameters as seen from the lateral and vertical equations.Hence, the estimation of the fault plane parameters becomes to a nonlinear optimization problem.Estimating p from g is formulated as the minimization of the function f given where k k is the Euclidean norm.The optimal parameter values p minimize the L 2 norm of the loss function given by 2.3, which is considered as an objective function [29].Therefore, the nonlinear unconstrained optimization problem can be written as where S is the domain of parameters.
The analytical solution of the problem seems to be unavailable because of the nonlinearity of fault plane model.The size of parameter space, the existence of local minima, the continuity of objective function and the sensitivity of objective function to each of the model parameters must be considered.In order to solve the inverse problem commonly used techniques such as simple Monte Carlo methods become ine¢ cient or impractical in very large solution spaces.More recently a number of guided search techniques from the …eld of arti…cial intelligence and heuristic methods have been developed.In this study, SA is used to estimate the overall earthquake geometry using the optimization problem given with lateral and vertical equations.

Application of the SA algorithm to the estimation problem
The SA algorithm derives its name from an analogy to the cooling of metals.As a metal cools, the atoms ‡uctuate between relatively higher and lower energy levels.If the temperature is dropped slowly enough, the atoms will all reach their ground state.This cooling process is called as "annealing".However, if the temperature is dropped too quickly, the system will get trapped in a less than optimum con…guration.So the temperature is considered as the control parameter in the method.If the energy function of this physical system is instead replaced by an objective function, then the progression of this function towards the global minimum is analogous to the physical progression towards the ground state.
Convenience of application to real world problems and obtainment of good solutions makes the SA algorithm be one of the most powerful and popular heuristics to solve many optimization problems.Some factors need to be considered when designing the SA algorithm.The following elements must be provided for implementing the SA to a problem: (i) are presentation of possible solutions, (ii) a generator of random changes in solutions, (iii) a means of evaluating the probability functions, and (iv) an annealing schedule, an initial temperature and rules for lowering it as the search progresses [1].The optimization problem given in equation 2.4 will be minimized to estimate the nine fault plane parameters with the SA algorithm.The algorithm has two cycles, inner and outer.A neighboring parameter set of the current parameter set is generated, in the inner cycle of the algorithm.If the new state is better than the current state then the generated solution replaces the current solution, otherwise the solution is accepted with a criterion probability.The value of the temperature, control parameter, decreases in each iteration of the outer cycle of the algorithm.In this regard, the steps of this algorithm are brie ‡y looked into: Step 1: T 0 is chosen as the initial temperature, p 0 is chosen arbitrarily as an initial parameter set on the de…nition space, n is the dimension of fault parameter set, T min is chosen as the …nal temperature (stopping criteria of outer cycle), L is chosen as the number of neighborhood solutions generated at a certain temperature T to get equilibrium (stopping criteria of inner cycle), a and b are chosen lower and upper bound sets for parameters respectively.Firstly, the initial temperature is considered as the system temperature, T k = T 0 for k = 0 and optimal parameter value p is set to p 0 , p = p 0 .
Step 2: Under kth temperature, if the inner loop condition is met, go to Step 3 ; otherwise the new parameter set is produced at a given temperature as where, p is constrained by p 2 [a; b] .Later on, the new function value f (p j+1 ) = f j+1 and f = f j+1 f j are computed.If f < 0, the new state is accepted.Otherwise, the Metropolis criterion is followed to accept p j+1 with a probability of and Step 2 continues.
Step 3: The temperature is reduced according to a speci…ed cooling schedule, T k = T 0 c k , 0 < c < 1, where c is analogous to the Boltzmann's constant that can be used to tune the algorithm.If outer loop break condition is met computation stops and optimal parameter set is reached.Otherwise, Step 2 is repeated.
The de…nition of the starting temperature, the cooling schedule of the temperature, the iteration number at each temperature and the de…nition of the stopping criteria have great roles in the e¢ ciency of the algorithm [12].Various SA algorithms are based on the same principle explained above, e.g., Boltzmann annealing [19], fast simulated annealing [28], very fast simulated annealing [16], and adaptive simulated annealing [17].These algorithms vary in probability distribution, annealing schedule, and the generation methods of random change.In this study, geometric cooling schedule is chosen and applied to earthquake source parameter estimation scheme.

Numerical Example
In this section, a numerical example is used to investigate the performance of the SA algorithm.An operation region is de…ned for performing the simulated data.The surface relocations, taken as outputs of the estimation scheme, are generated on the experimental region for de…ned coordinates using Matlab code.It is tried to get the optimal values of fault plane parameters.In the conclusion part of the simulation study, SA and Nelder-Mead simplex algorithm are compared in terms of their convergency to optimal value of nine source parameters.

4.1.
Data.An operation region, a quadrangular area, is de…ned for performing simulated data and considered as a de…nite place of earth surface where earthquake has been occured.The simulated earthquake area is formed in ( 30; 20) ( 50; 20) coordinates with 50 geodetic points which are signed on the graph given in Figure 2.

Figure 2. De…ned 50 coordinates around the fault direction
As can be seen from the Figure 2 that the fault direct, passed along the origin with a straight line, is formed and 50 random generated coordinate values (x; y) are …xed around the fault direction.The surface relocations, denoted as [x; y; z], are generated by using Matlab code taken from the geodynamics laboratory page (http:// www.gpsg.mit.edu) of the Massachusetts Technology Institute for each coordinate.The …xed 50 coordinates and the surface relocations are given in Table 1.The fault plane parameters, used for synthetic data set and the lower and the upper bounds of the parameters are presented in Table 2.In the absence of any prior information about where the good feasible solutions might lie, it's reasonable to set each parameters midway between its lower bound and upper bound in order to start the search in the middle of the feasible region.For a sample run, initial parameter set is taken p 0 = [60 10 2.5 1.48355 4.5978 -25 -25 0 0].The initial temperature and Boltzmann's constant are considered T 0 and c = 0:9, respectively.The inner loop break condition is taken n 10, where n is the dimension of parameter set.Table 3 shows the results of SA algorithm for di¤erent stopping criteria (T min ).For each T min , the iteration number and computation time (CPU) are shown together with the objective function value (f ) and estimated values of source parameters (p ).It can be seen from the Table 3 that when the T min is decreased the loss function value gets smaller and the calculation period gets longer.The objective function takes the smallest value when the stopping criteria is choosen 10 12 .The change on temperature and the minimization of loss function value during the search process of SA are given in Figure 3.It can be seen from Figure 4 that while the temperature drops constantly from the beginning of the search, the objective function value given in equation 2.3 reaches the global minimum value.
Table 3.Estimated values of source parameters for T 0 = 100 and c = 0:9   2, while the algorithm gets the global minimum.
A single run of global converging algorithms is not su¢ cient to …nd the global solutions.The model parameters obtained by several runs may di¤er from each other because of the stochastic structure of the algorithm.Since the SA is a probabilistic search algorithm, di¤erent runs may result in dissimilar con…gurations.Therefore, several seperate runs are performed to ensure that the true global minimum has been located.The average results of 100 Monte Carlo simulations for the di¤erent initial values of fault plane parameters and average computation times (ACPU) are reported in Table 4.The initial temperature, Boltzmann's constant and stopping criteria are choosen 100, 0.9, and 10 8 respectively for all test runs.It has been observed from the simulation results that the average loss function value converges the global minimum for the arbitrary initials.It can be easily said that SA largely independent of the initial values and it can escape local minima through selective uphill moves.
Table 4.The loss function values for the di¤erent initials 4.3.Results of the Nelder-Mead simplex method.tests have been made on starting points chosen near the optimum and quite far from the optimum generated inside the parameter interval.Four scalar parameters must be speci…ed to de…ne a complete Nelder-Mead simplex method; coe¢ cients of re ‡ection , expansion , contraction , and shrinkage .These parameters are chosen to satisfy > 0, > 1, 0 < < 1, and 0 < < 1.In this study, Nelder-Mead simplex algorithm, composed by [14], is applied to the source parameter estimation scheme.The coe¢ cients of re ‡ection, expansion, contraction and shrinkage are choosen = 1, = 2, = 1 2 , and = 1 2 , respectively.The stopping criteria is considered " = 0:00001.The results of Monte Carlo simulations are given in Table 5.It can be seen from the results that the Nelder-Mead simplex method requires stringent initial estimates for the model parameters.The method is only advantageous when a good initial estimate is provided.If the initial vertices are generated far from optimum, the loss function (f ) can not reach the global minimum and CPU time gets large.Table 5.The simulation results of Nelder-Mead simplex algorithm

Conclusions
Estimation of fault plane parameters using classical optimization methods, e.g., derivative-nonderivative based optimization methods, Monte-Carlo methods, Least Squares methods, etc. is often di¢ cult because of the nonlinearity of the problem.These methods need some assumptions to implement on problems and also may fail to get global solutions of optimization problems.In this paper, SA algorithm is presented and compared with a derivative-free optimization method, Nelder-Mead simplex method, for estimating the source parameters.The SA algorithm has many advantages over other optimization methods.The method has advantages in dealing with strong nonlinearities and discontiniuties.The SA algorithm sometimes goes uphill as well as downhill unlike other direct optimization methods.This is an ability to avoid becoming trapped at local minima.
During the Monte Carlo simulation study, it's realized that assessing the performance of the Nelder-Mead simplex method is generally problematic due to the nonconverge of optimal parameter values even if it is more e¢ cient than other direct search optimization methods.The simplex method requires strict initial conditions, whereas the SA algorithm is able to explore the parameter space and focus on the most promising area without prior knowledge of its location.The simulation results show that the SA algorithm is an e¢ cient derivative-free optimization algorithm with respect to Nelder-Mead simplex method.However, the SA has some drawbacks.The most remarkable disadvantages of the SA algorithm are (i) spending much computing time to …nd the optimum solution and (ii) determining the proper cooling schedule.In order to obtain better results, the various tunable parameters used (e.g.initial temperature, …nal temperature, cooling rate etc.) needed to be chosen carefully depending on the problem variety.
Due to its versatility and independence on prior knowledge of the parameter values, the SA algorithm is particularly applicable to estimate the parameters of fault plane model that are not solvable using analytical methods.It can be said that the SA algorithm is a promising tool for estimation of parameters in geosciences.However, the utility of this procedure to estimate fault plane parameters from real data has yet to be proven.In the later studies, it's possible to improve the algorithm by using di¤erent cooling schedules.Using a hybrid optimization strategy based on SA and gradient based or direct optimization algoritm can provide more e¢ cient and accurate estimation of parameters.On the other hand as there are other well-established heuristics, such as genetic algorithms, Bee and Ant Colony, etc., a discussion of how good the estimation can be achieved by using them, may lead to a better estimation procedure.

Figure 1 .
Figure 1.The fault geometry and the relocation on the fault

Figure3.
Figure3.Cooling schedule and loss (mis…t) function value of SA for T min = 10 12

Figure 4 .
Figure 4. Estimation results of fault plane parameters

Table 1 .
The coordinates of locations and relocations Table2.The true values and the de…nition intervals of the fault plane parameters Results of the SA algorithm.The objective function is analyzed in the parameter space during the parameter optimization procedure with SA algorithm.