Randomly distributed unit sources to enhance optimization in tsunami waveform inversion

In tsunami waveform inversion using the conventional Green’s function technique, an optimal solution is sometimes difficult to obtain because of various factors. This study proposes a new method to both optimize the determination of the unknown parameters and introduce a global optimization method for tsunami waveform inversion. We utilize a genetic algorithm that further enhanced by a pattern search method to find an optimal distribution of unit source locations prior to the inversion. Unlike the conventional method that characterized by equidistant unit sources, our method generates a random spatial distribution of unit sources inside the inverse region. This leads to a better approximation of the initial profile of a tsunami. The method has been tested using an artificial tsunami source with real bathymetry data. Comparison results demonstrate that the proposed method has considerably outperformed the conventional one in terms of model accuracy.


Introduction
Direct observation of sea surface deformation after the occurrence of an earthquake is still difficult to obtain; therefore, its estimation is often performed by consideration of relevant seismic information or the hydrodynamic response of the sea determined from recorded tsunami waveforms.Determination of sea surface deformation generated by earthquakes is crucial to the success of tsunami modeling.One of the most frequently used methods for determining sea surface deformation is to presume it from a fault model (Mansinha and Smylie, 1971;Okada, 1985).A more realistic approach was proposed by Satake (1987) who analyzed recorded waveforms to infer earthquake source parameters or particularly, coseismic slip, using the Green's function technique.Even though the fault model is still required, the division of a fault into smaller sub-faults allows the slip to be estimated in a heterogeneous manner, which leads to better approximation of sea surface deformation.A simpler method was actually introduced earlier by Aida (1972) for which no prior assumption of a fault model was needed.This study is in line with that of Aida because we are more interested in estimating sea surface deformation than a slip on the fault plane.The motivation behind this is that tsunami excitation can sometimes occur as a result of various factors that are independent of the associated seismic characteristics (Geist, 2002).
More recently, several studies using tsunami waveform inversion to estimate sea surface deformation without fault model assumptions have been widely developed.The basic premise is to replace the fault model by an auxiliary basis function on unit sources, which is equivalent to the sub-fault approach by Satake (1987).For instances, Baba et al. (2005) used a simplified fault model by disregarding actual earthquake parameters to produce the initial profile on each unit source, whereas Satake et al. (2005) proposed a more direct approximation using a pyramidal shape with a flat top.Other studies by Liu and Wang (2008) and Saito et al. (2010) demonstrated attempts to use Gaussian function, whereas Wu and Ho (2011) adopted a top-hat small unit source to represent the initial profile.The same approach was proposed by Tsushima et al. (2011) and Yasuda and Mase (2013) for the more practical purpose of a tsunami early warning system.However, this inversion method for tsunami waveforms possesses a limitation, in that the inverse matrix does not always exist because of the non-uniqueness of the solution.In addition to the large number of unknown parameters, which might produce many local optima on the misfit function mea-Published by Copernicus Publications on behalf of the European Geosciences Union.
sure, the search towards optimality is confined by the uniform distance of unit sources used in the regular Green's function.
Tsunami waveform inversion sometimes falls into an illposed problem, in which small errors in the observed waveforms are exceptionally amplified in the solution.Therefore, both the uniqueness and the stability of solutions are sometimes difficult to attain without appropriate treatments.The most frequently employed techniques to maintain a stable solution is to use a smoothing constraint (Gusman et al., 2010(Gusman et al., , 2013;;Saito et al., 2010).Other than that, Koike et al. (2003) suggested reducing the unknown parameters using the wavelet base to guarantee the uniqueness of the solution.However, they found later that the selection of the wavelet base was not straightforward.Another effort to overcome the issue was discussed by Voronina (2011).The study promoted a method to control numerical stability for the ill-posed problem in tsunami waveform inversion by means of singular value decomposition and r-solutions techniques.In this paper, we proposed a new approach to tackle the same problem by determining the optimal position or spatial distribution of unit sources located around the tsunami source or epicenter.A genetic algorithm (GA) as a global optimization method, combined with a pattern search (PS) method, is employed to search the mentioned positions prior to the inversion.As the selected positions are probably located in between the initial unit sources, interpolations are performed during the optimization.Therefore, the Green's function evolves dynamically at each generation of the GA and PS iteration.

Inversion method
Generally, the characteristics of tsunami propagation in deep water are linear.According to Satake (1987), even in shallow coastal areas, the first leading waves recorded at coastal tide gauges are still well simulated by the linear long wave model.Therefore, a typical linear non-dispersive shallow water equation is used in the forward modeling to compute time histories of sea surface elevation at the specified observation points: where η is the water elevation of the tsunami, V (u, v) is the depth-averaged horizontal fluid velocity vector, d is the water depth, and g is gravitational acceleration.Equation ( 1) is completed with two types of boundary conditions, as follows: where c = √ gh is the wave speed and n is the unit vector normal to the boundary.The upper and lower expression in Eq. ( 2) represent open and closed boundary respectively.The Green's function G i (x, y, t) indicating the recorded synthetic waveforms at a location (x, y) corresponding to the ith unit sources is then constructed.This tsunami Green's function is derived from 2-D linear shallow water equation, Eq. (1), solved using numerical method with specified bathymetry resolutions.The total sea surface fluctuations at (x, y, t) can be expressed as follows: where w i is the weighting factor for each unit source or the unknown parameters to be determined by the inversion and N is the number of the unit source.Equation ( 3) is developed based on the assumption of linear superposition considering the nonlinearity for tsunamis in ocean basins to be small, such that it can be neglected (Liu and Wang, 2008).In vector notation, Eq. ( 3) can be reformulated as follows: Based on the least squares method, the vector of the unknown parameters w can be acquired by solving the following inverse equation: In this study, we assume that the generation of the initial profile of the tsunami occur instantaneously.Therefore, incorporating the time variation or transient deformation of the water surface is unnecessary.

Global optimization method
The ultimate purpose of a global optimization method is to find the extreme value of a given non-convex function in a certain feasible region.Following the growth of computer science, new types of optimization based on natural processes and artificial intelligence have been developed extensively and used by scientific and engineering communities.
The reason for this is that the new optimization methods possess the interesting feature of being able to avoid local optimum solutions, which is something classical methods fail to do.The use of global optimization methods in tsunami waveform inversion is not new; relevant discussions can be found in Piatanesi and Lorito (2007) and Romano et al. (2010).They used a simulated annealing technique to solve the inverse problems.Here, we proposed a different algorithm based on a hybrid optimization of GA and PS.The hybrid technique is preferred because global optimization methods, such as GA, are capable of exploring broader search space, but not as good in fine tuning the approximation of the expected solution.Therefore, PS is employed as a local search algorithm to locate other nearby solutions that could possibly be better than the result of the previous search by GA (Payne and Eppstein, 2005;Costa et al., 2010).
The hybrid algorithm proposed in this study works by simply treating the output of the GA optimization result as the initial condition for the PS optimization.The technique is proven effective even though more fitness function evaluation is required; hence, it costs extra computational efforts.However, parallelization of either GA or PS can be easily implemented to expedite the computing time and gain substantial performance enhancement.
The formulation of an optimization problem can be expressed as follows: where x ∈ X is the vector of design parameters, f : X → R is the cost function, and X ⊂ R n is the constraint set or bounds that can be defined as follows: with −∞ ≤ l i ≤ u i ≤ ∞, for all i ∈ {1, . .., n}, where l and u are the lower and upper bounds, respectively.
In this study, we develop the proposed algorithm with two different design parameters.The first is simply to search the water elevation of each unit source without including a search of the optimum locations, thus similar to that of the ordinary least squares method.By rearranging the Eq. ( 4), the first design parameter can be expressed as an optimization problem as follows: where x is the unknown parameters w in the Eq. ( 4) and X is the constraint or bounds representing plausible values of subsidence and uplift of the water surface acted as lower and upper bounds respectively.This experiment aims to gauge the significance of the global optimization method applied to a linear system compared to the conventional method using least squares.The purpose of the second design parameters is to search the optimal location of unit sources, while the initial water elevations are calculated using the least squares inversion, where w are the calculated water elevations using least squares, Eq. ( 5), and x define a unique identification for the computational grid of the forward model located inside the inverse region X (Fig. 1).In this part, the combination of the algorithm with the least squares method is necessary to save the computational efforts.

Genetic algorithm
GA is an optimization method that searches for an optimal value of a complex function by adopting the process of natural evolution (Goldberg, 1989).It can be categorized as a type of stochastic optimization method and as a part of artificial intelligence.In GA, the model parameters or decision variables in the optimization are first transformed into a chromosome-like data structure that later evolves to form a better individual.The most common representation of design parameters in GA is their encoding into a binary string.There are three basic genetic operators in GA: selection, crossover, and mutation.As described in Wetter and Wright (2003), in GA, finite lower and upper bounds of X are discretized in the mesh M (l, 1), and all operators are such that the real value of any point is an element of M (l, 1) ∩ X.Given a non-empty set X ⊂ R n and a non-zero M ∈ N, X M is defined as the set of all sequence X with M elements, Consider B as the set containing all elements of M (l, 1) ∩ X.Similarly, B M is the set of all sequences in B with M elements, We denote k ∈ N as the generation number, M ∈ N as the population size, x k ∈ X M as the M points in the kth generation, y k ∈ B M as the binary representation of x k , and y k,i denotes the ith element of y k for i ∈ {1, . .., M}, which is a binary representation of a point in R n .The step-by-step operation of the GA is described in the following: Step (1): randomize the initial population y 0 ∈ B M of M randomly generated points.Then, evaluate all y k according to the specified cost function; thus, a fitness function : N× B M → R M computes a fitness value of y k .
Step (2): select a pair of individuals with the selection function ϑ : N × B M → B 2 .Individuals with better fitness values have a higher probability of being selected.Step (3): swap the bits (genes) between the selected individuals (chromosomes) to produce new offspring, φ : Step ( 4): the final genetic operator is termed mutation and aims to maintain genetic diversity.The mutation function ψ : N × B × [0, 1] → B alters a bit on individuals from its initial state.

Pattern Search
Similar to GA, PS is an optimization method that does not require the gradient (derivative-free) of the problem to be optimized.The method was first introduced by Hooke and Jeeves (1961).Later, Torczon (1997) conducted studies to prove the convergence of PS using the theory of positive bases.The algorithm of PS used in this study is similar to that of Wetter and Write (2003), while the design parameters are identical to the GA optimization.
For the same optimization problem as formulated in Eqs. ( 8) and ( 9), the algorithm searches a lower cost function value than f (x k ), where x k ∈ X denotes a solution at the current iterate and k ∈ N denotes the iteration number.The search takes place on the points in the set where k > 0 is the mesh size factor, s ∈ R n is a fixed parameter to scale the design parameters, and e is the unknown approximation error.The rule to select a finite number of points in X on a mesh can be defined by the following: where x 0 ∈ X is the initial iterate.
The overall procedures of PS optimization can be elaborated as follows: Step (1): initialize x 0 ∈ X and 0 > 0, for which in our case, x 0 is obtained from the GA's output.
Step (3): if f x ≥ f (x k ) for all x ∈ γ k , the search continues with x k+1 = x k and reduced mesh factor k+1 = k 2 .The search should continue until a stopping condition is satisfied; e.g., the mesh has been refined a user-specified number of time or less than a given tolerance.

Numerical experiments
We have conducted numerical experiments using an artificial tsunami source propagated on an actual bathymetry profile.An area extending from 140-145 • E and 35-41 • N is chosen as the domain of interest.The selected domain resembles that used in most studies of tsunami waveform inversion for the 2011 Tohoku tsunami.The resolution of the numerical model is 1 arcmin, which is consistent with the resolution of the bathymetry data obtained from the ETOPO1.ETOPO1 is a global relief model of the earth's surface, which includes ocean bathymetry, available from the National Geophysical Data Center of the National Oceanic and Atmospheric Administration (Amante and Eakins, 2009).There are eight artificial observation stations associated with the actual location of gauges within the study area.The tsunami source is divided into 28 unit sources that are distributed uniformly around the actual epicenter of the 2011 Tohoku tsunami (Fig. 2).

Cost function
The selection of a cost function or misfit function is essential because it directs the fate of the optimization towards the  optimal solution.Here, we use a combination of root mean square error (RMSE) and Pearson correlation coefficient (r).While the RMSE is sensitive to amplitude matching, the correlation coefficient is more sensitive to phasing between the compared series (Barnston, 1992).The RMSE serves to aggregate the individual differences of data points into a single measure of predictive power, which is defined as follows: where y is the predicted value by the model, d is the measurement data for each ith data point, or in our case, the waveform generated by the artificial tsunami source, and n is the total number of data.The Pearson correlation coefficient is defined as a division of the covariance of the two variables by the product of their standard deviations: where d and y represent the means of d and y, respectively.The correlation or relevance of the data is measured using the value R = 0.5 (r + 1) to avoid negative values in the case of a decreasing linear relationship between the compared series.
The fitness evaluation is subject to noise from various factors that might lead the optimization towards unexpected solutions.The easiest technique to overcome such problems is by means of explicit averaging over a number of samples to smooth the cost function (Jin and Branke, 2005).As the best solution or closest fit is indicated by RMSE → 0 and R → 1, the final cost function is a summation of the mean of the tth sample over the time window T of Eqs. ( 14) and ( 15), which can be written as follows:

Random initial unit source positions Interpolation
where k denotes the respective time window and N is the total number of windows.

Basis function
A Gaussian shape with 1 m amplitude is used as the basis function for each unit source (Liu and Wang, 2008).Providing a i as the amplitude of a unit source with the centroid positions of x i and y i , the basis function can be written as follows:  where z i (x, y) is the initial water surface corresponding to the ith unit source, x and y are the locations of computational grids, and σ is the spread of the blob with a length of 40 km (Fig. 3).The specified length should satisfy the long wave assumption (h/L < 0.05), i.e., the wavelength should be greater than 20 times the average water depth.

Model development
For the first design parameters, the optimization is performed merely to search the water elevation or initial amplitudes of each unit source.For this case, the Green's function is constructed based on the initial 28 unit sources, separated by a uniform distance of 60 km, which is identical to that used in the least squares inversion.Hereafter, the first model will be termed the Genetic Algorithm Pattern Search for uniform source distribution (GAPSu).The purpose of this model is simply to compare the performance of the proposed global optimization method with the traditional least squares method in the same model design and environment.
The second design parameters aim to find the optimal locations of unit sources that, at the initial state, are distributed randomly around the tsunami source.The second model will be termed the Genetic Algorithm Pattern Search for random source distribution (GAPSr).In the GAPSr model, the amplitudes are computed using the least squares inversion; therefore, the second model is actually a combination of a deterministic and stochastic optimization.However, the selected source locations might not lay precisely in the initial 28 unit sources and thus, interpolation is required to produce the sea surface fluctuations at the observation stations.Consequently, in the GAPSr model, the Green's function evolves during the optimization.Nearest neighbor-weighted interpolation is performed for estimating the phases and amplitudes of the waveforms originating from the selected locations based on the four nearest unit sources (Mulia and  Asano, 2014).An example of the interpolation results, complete with statistical evaluations in terms of RMSE and r, is shown in Fig. 4. Despite the satisfying results, based on the measure of fitness shown by the interpolation method (overall RMSE < 0.048 m and r > 0.966), the small errors might be amplified in the solution because of the ill-posed problem.Consequently, further improvements should be made to suppress the generated errors.
An artificial tsunami source is used to test our method (Fig. 6a).Instead of using a simpler profile produced by the Okada's solution, we generate a more complex shape from a superposition of 10 unit sources with random amplitudes and positions located inside the inverse region.We purposely limit the number of unit source to avoid generating shorter wavelengths than the prescribed long wave assumption.For the regular Green's function (first design parameter), using only 10 unit sources is insufficient to reconstruct the target profile.Consequently, the number of unit source should be increased.The decision of using 28 unit sources in both design parameters is for the purpose of model performance comparisons that will be further discussed in the Results and discussion section.The approximation of this initial profile is performed using unconstrained, traditional least squares inversion, GAPSu, and GAPSr.However, GAPSr is the most important part in this study; therefore, we focus our discussion on the GAPSr model.The development of the GAPSr model depicted in the flowchart (Fig. 5) can be summarized as follows: 1. Construct the initial Green's function.
2. Initialize the GAPSr model by randomly distributing the unit source locations.The search of the optimal location is bounded by the area of the inverse region.
3. Perform interpolation and update the Green's function.
4. Evaluate the fitness by performing the least squares inversion.
5. After reaching the stopping criteria, the forward numerical modeling is run again for each of the optimized unit sources to avoid errors generated from the interpolation result.Subsequently, the inversion is performed for the final time.

Results and discussion
Overall, all models can produce a relatively good estimation of the targeted sea surface deformation.This is likely because they are applied to an ideal case with artificial conditions and settings, except for bathymetry.For instance, the target source was generated from the same Gaussian shape as that used to construct the Green's function.Therefore, the task is more straightforward as fewer complexities are encountered.
In the real case, however, alternative distributions to represent the initial water height should be considered because the tsunami source does not always follow the Gaussian distribution.This may have a profound effect on the result of the first design parameter, but less significant for the second design parameter as a superposition of unit sources with random locations allows to approximate any surface profile regardless of their shape.Nevertheless, in this study, the use of the ideal case has allowed us to assess the advantage of the proposed method in a more detailed manner.
The GAPSu model yields a slightly better fit of waveforms compared with the least squares method (Table 1).This means that the global optimization method locates a better minimum value in the cost function, which is situated beyond the reach of the least squares method.However, the slight refinement by the GAPSu model over the least squares method makes it difficult to gauge the benefits of employing the method.This should not be a surprise because the model is applied to determine the coefficients in a linear system, which is relatively easy to solve using a conventional method.Moreover, the waveforms used to invert the initial sea surface deformation are generated from an artificial tsunami source instead of real measurements.Therefore, the linearity is well conserved and thus, the use of more advanced methods becomes redundant and unnecessary.In the study by Piatanesi and Lorito (2007), a global optimization method was suc- cessfully promoted for the case of tsunami waveform inversion.This was because the optimization method was applied to a nonlinear inverse problem of actual measurement data.Accordingly, the appraisal of such a method can be clearly defined.
Spurious uplifts and subsidence of the water surface profile are generated in both the least squares and GAPSu model results (see Fig. 6b and c).One may argue that the specified spatial resolution of the unit sources is too coarse to represent the complete form of the target source.A denser distribution of unit sources should improve the results; however, it might also introduce other problems.A large number of model parameters (unknown parameters), which are proportional to the degrees of freedom in the optimization, are liable to cause the solution to become easily entrapped in a local optimum.Without a smoothing constraint, the result of the tsunami waveform inversion might be bumpy and non-physical, especially for cases with high spatial resolution (Wu and Ho, 2011).In other studies on tsunami waveform inversion by Baba et al. (2005) and Wu and Ho (2011), an equality constraint was imposed to maintain the smoothness of the inverted parameters to satisfy the long wave assumption, while by Saito et al. (2010)   the constraint was used to obtain stable solutions.However, such constraint might restrict the exploration throughout the feasible search space and render the discovery of an optimum solution more difficult.Another plausible explanation for the unsatisfactory results of the least squares and GAPSu model is simply that the equidistant unit sources in the regular Green's function confines the search for optimality.This can be proven by the result of the GAPSr model, for which the same number of model parameter (28 unit sources) with random locations yields a much better estimation.
The different design parameters in the GAPSr model have considerably improved the inversion accuracy.For instance, at Gauge 1, where the best fit of the waveform is attained, the measurement of accuracy as RMSE is 0.0260, 0.0256, and 0.0094 m, and as r is 0.9973, 0.9974, 0.9996, for the least squares, GAPSu, and GAPSr models, respectively (Table 1).For further qualitative or visual assessments, comparisons of time series of the waveforms at each gauge are presented in Fig. 7.In addition, statistical evaluation results for the waveforms at all gauges are shown by scatter plots in Fig. 8. Overall, the statistical analyses on the waveforms suggest that the GAPSr model shows very good agreement with the target.These results conform to the inverted sea surface deformation resulted by the GAPSr model, which can be seen in Fig. 6d.The random location of the unit sources allows the approximation to capture the exact profile of the target source.Other than this, despite no smoothing constraint being used, the inverted sea surface deformation remains smooth and coherent.
The search for the optimal locations of the unit sources allows the least squares method to find the unique and optimal solution.Such an approach is difficult to achieve deterministically using conventional gradient methods, because there is the possibility that the constructed design parameters in the GAPSr model produce a discontinuous or non-differentiable error surface owing to the random characteristics exhibited by the artificial tsunami source (target source).The same characteristic is very likely to occur in nature.A high degree of uncertainty has been observed in tsunami sources leading to significant variations in the nearshore tsunami amplitude (Geist, 2005).Therefore, the use of the model in real case applications is encouraged to reveal the underlying dynamics in tsunami generation.However, as with typical stochastic methods, the proposed model cannot ensure a constant op-timization response time.The solution and convergence are strongly dependent on the random initial state.

Conclusions
Estimations of tsunami sea surface deformation using a global optimization method with a stochastic nature have been conducted.Our numerical experiments using the GAPSu model revealed that the use of such methods for a linear system with standard design parameters, as in ordinary tsunami waveform inversions, is redundant and promotes trivial improvements.In contrast, the different design parameters in our proposed method (GAPSr), which was applied to determine the optimum location of the unit sources prior to the inversion, demonstrated considerable improvements in the accuracy.The random location of unit sources permitted the inversion to produce a more precise approximation of the initial sea surface deformation without violating the general assumption of long wave theory.
The involvement of stochastic processes in the optimization increased the ability to reveal uncertainties in the tsunami source, which are difficult to discern using deterministic approaches.However, as the signature of typical stochastic optimizations, the optimization response time is erratic because of the strong dependency on the initial state.Using a current standard desktop computer, the required computing time varied from 5 to 10 min.Thus, a more sophisticated computer would be needed to ensure the effectiveness of the method when applied in a real-time application.Nevertheless, the results have demonstrated the efficacy of the method for post-event studies of tsunamis, because it can provide better estimations of the coseismic sea surface deformation compared with traditional tsunami waveform inversion methods.
Edited by: I. Didenkulova Reviewed by: A. Gusman and two anonymous referees

Figure 1 .
Figure 1.Example of a unique identification of the forward model computational grid used for the second design parameters.

Figure 2 .
Figure 2. Study area and bathymetry profile.Red dots indicate unit sources located throughout the inverse region.Green triangles with numbers are artificial observation stations.

Figure 4 .
Figure 4. Waveform interpolation.Left-hand figure shows selected location of a unit source indicated by a black square.Blue dots represent the four nearest unit sources used in the interpolation.Right-hand figures are comparisons of waveforms between numerical model (solid black line) and interpolation (dashed red line) at the artificial gauges.

Figure 5 .
Figure 5. Flowchart of model development procedure.

Figure 6 .
Figure 6.Sea surface deformation.Gray dots indicate the centroid of unit sources.(a) Artificial tsunami source as the target to be approximated, (b) inverted source using least squares method, (c) GAPSu model, (d) GAPSr model.

Figure 7 .
Figure 7.Comparison of waveforms at gauges.Gray bar above the time axis indicates the time range for the inversion.

Figure 8 .
Figure 8. Scatter plots of each inversion method with respect to the waveforms of the target source at all gauges.

Table 1 .
Summary of statistical evaluations based of the proposed cost function.