Dynamic inversion of the rupture parameters on fault system with complex geometry: A GPU parallel genetic algorithm based on BIEM

1.   Introduction

The mechanism and rupture process of an earthquake can be revealed by seismic source inversion. Inversion on earthquake source can be divided into two categories, kinematic and dynamic. Kinematic inversion on earthquake source is often used to reveal the spatio-temporal history of ruptures on earthquake faults by using recorded seismic data. And usually, a linear inversion scheme is adopted to obtain the rupture process of an earthquake shortly after it begins. Although the source kinematic inversion can rebuild the detailed rupture processes, the physical parameters describing the stress and friction property on the fault cannot be obtained. In dynamic inversion, on the contrary, the source parameters are obtained based on the dynamic simulation of the spontaneous rupture process of the fault and the inversion algorithm. Since stress field, which is responsible for the rupture process of earthquake source, is introduced, dynamic inversion is a very important tool to reveal the physical nature of an earthquake. However, researches on this field are relatively rare. The reason may lie in two aspects. One is that dynamic forward modelling itself is very time-consuming, and the other is that since a non-linear friction law is involved in the forward problem, non-linear inversion algorithm must used, which usually requires a large amount of forward simulations.

The first attempts at dynamic inversion were based on the conversion of source models obtained by kinematic inversion into dynamic models (Fukuyama and Mikumo, 1993; Bouchon, 1997; Ide and Takeo, 1997). Peyrat and Olsen (2004), Corish et al. (2007) and Di Carli et al. (2010) inverted the 2000 M_W6.7 Tottori earthquake in Japan directly for the first time. Peyrat and Olsen (2004) developed a non-linear proximity algorithm inversion method for dynamic rupture inversion of the same earthquake, based on forward modelling with the 3D finite difference method. Di Carli et al. (2010) used an elliptical patch model proposed by Vallée and Bouchon (2004) to capture the long-wavelength features of an inverted event. For these early experiences, the most common fault model in a dynamic inversion is the circular fault model that introduced by Kostrov (1964, 1966) and elliptical fault model (Peyrat and Olsen, 2004; Corish et al., 2007; Di Carli et al., 2010; Díaz-Mojica et al., 2014; Mirwald et al., 2019). The main issue in dynamic inversion is the parameterization of the source and the computation of its rupture process, which are very non-linearly related to the seismic waves emitted by the earthquake. In spite of these difficulties, with present resources, it is possible to do limited non-linear inversion for seismic sources in the single elliptical approximation. Díaz-Mojica et al. (2014) used a parallel genetic algorithm to perform dynamic source inversion on the 2011 M_W6.5 Zumpango earthquake. They used an elliptical fault model, which is described by five geometric parameters. Besides, four physical parameters are introduced in the dynamic inversion, i.e. the initial stress, the nucleation zone stress, the critical slip weakening distance, and the friction coefficient. They also analyzed some derived parameters, such as energy, magnitude, and radiation efficiency.

Physical parameters on friction, especially the critical slip weakening distance D_c, play a significant role in the dynamic process of earthquake rupture. Many studies have attempted to constrain D_c in the field for several earthquakes based on dynamic (Olsen et al., 1997; Nielsen and Olsen, 2000; Dunham and Archuleta, 2004; Peyrat and Olsen, 2004; Ma et al., 2008; Goto and Sawada, 2010) and kinematic rupture models (Ide and Takeo, 1997; Tinti et al., 2004, 2009; Galetzka et al., 2015), as well as direct estimations from near field ground motion observations (Mikumo and Yagi, 2003; Fukuyama and Mikumo, 2007; Kaneko et al., 2017). Weng and Yang (2018) used a plane fault model to perform dynamic inversion of the 2015 M_W7.8 Nepal earthquake. They applied PyLith (Aagaard et al., 2013) to conduct the dynamic rupture simulation, and combined the results of kinematic inversion to determine the friction parameters for the earthquake. Their results of inversion and the actual observations are in good agreement.

In this study, we performed a dynamic inversion to estimate the source parameters of a branched fault system. The boundary integral equation method (BIEM) with unstructured meshed is used in the forward modeling (Qian et al., 2019), and a GPU based parallel genetic algorithm (GA) is adopted in the inversion scheme. For simplicity, only two parameters, the initial stress T ⁰ and the critical slip-weakening distance D_c, are to be inverted in our present study. We performed inversion processes on different inherited parameters to explore their effects on the inversion results.

2.   Numerical method and fault model

2.1   Forward problem

In the forward problem, the boundary integral equation method (BIEM) with unstructured meshes (Qian et al., 2019) is used to simulate the rupture process. The discretized BIE can be expressed as follows

where T_mp is the stress of the element m at time p, $T_m^0$ is the initial stress, V_nq is the slip rate of the element n at time q. The integration kernel K_mp/nq is the stress Green’s function in whole space (Feng and Zhang, 2017), which represents the stress response between fault elements. The BIEM together with the slip-weakening friction law are used to obtain the slip rates and stresses on the fault elements. Given the fault geometry as well as the initial tectonic stress field and the friction parameters, the spatio-temporal evolution of the rupture process can be obtained.

In the dynamic models, we apply a slip-weakening friction law on the fault as follows,

where µ_s and µ_d are the static and dynamic friction coefficients, respectively, D is the fault slip. We set up a constant normal stress σ_n=100 MPa and constant friction coefficients µ_s=0.6 and µ_d=0.1. A reference value of D_c is defined as D₀=0.57 m.

For given parameters of medium and fault model with a certain geometry, the integral kernels are calculated, and the dynamic rupture process is obtained by combining the slip-weakening friction law in the forward modeling. The dynamic forward simulation itself is very time-consuming, so that huge resources for storage and high computation speed are required in dynamic inversion, which is impossible to perform without an efficient forward computation. In this study, we develop a set of techniques to speed up the forward simulation, making the attempt to a dynamic inversion of the parameters for dynamic ruptures feasible.

2.2   Genetic algorithm

We use the genetic algorithm (GA) (Holland, 1975; Goldberg, 1989) to achieve the dynamic inversion of the source parameters. The GA is a non-linear inversion method based on biological evolution and the basic process of GA is shown in Figure 1. To apply the GA to seismic source inversion, we first need to establish the initial parameter population, which represented in a 10-bit binary form. Then for each parameter individual in the initial parameter population, we can obtain the corresponding rupture result through forward simulation. The parametric population of each generation needs the number of forward simulations corresponding to that of its populace. Therefore, the amount of forward calculations is substantial.

Figure  1.  Flowchart of the genetic algorithm applied in this study

DownLoad: Full-Size Img PowerPoint

The rupture results, which is characterized by the final slip, obtained by forward modeling need to be compared with the results of the target parameters to select the desired parameter individuals. In the GA, the selectivity of genes is the key to the speed of its iterative convergence. How to choose the best genes will determine the efficiency and accuracy of evolution. In our research, we use the Pearson correlation coefficient (ρ) (Williams, 1996) to measure the superiority of genes,

where $R_{\rm{tar}}^{\rm{slip}}$ is the final slip of target rupture result, $R_{\rm{ga}}^{\rm{slip}}$ is the final slip of rupture result of a particular parameter individual, E is mathematical expectation and the value of ρ, calculated by MATLAB’s built-in functions, stands for the correlation of the rupture result between the parameter individual and the target parameter. In the theory of evolution, good genes have higher heritability. Therefore, we use a probability model to select genetic genes, and the selection probability of each parameter individual depends on its Pearson correlation coefficient. So there is no absolute selection range that is correlated to good genes. The greater the correlation, the greater the probability. So that the population iteration can converge to the target parameter. From our inversion experience, when the entire population converges, the correlation coefficient is generally greater than 0.99.

The inheritance of organisms usually involves two essential processes, crossover and mutation. When the rupture result of the parametric population is obtained through forward modeling, and the corresponding correlation is calculated, we can select excellent genes for inheritance. Crossover is to split gene fragments of two parameter individuals and then reassemble them into two new gene fragments. In the genetic iteration of parameter population, although a continuous selection of good genes for hybrid evolution can make the population genes better and make the parametric population converge closer to the target parameter, it may converge to a locally optimal solution. Therefore, mutations in genes are needed. Mutation refers to a probabilistic gene mutation during the genetic process, which results in new characteristics. In our algorithm, when two good genes are crossed, a smaller mutation probability will be provided. That is, some fragments in the newly generated gene chain are mutated, which makes it possible that the newly generated gene chain is significantly different from its parents. The significance of mutation lies in that if the mutated individual is not a functional gene, it will still be eliminated, and once it is a more excellent gene, its gene may dominate the generation-by-generation inheritance, making the parameter population towards a better genetic evolution.

2.3   Source model and parameters

In previous studies, the elliptic fault model, which is simple to handle, is the most frequently used one in dynamic inversion so far (Peyrat and Olsen, 2004; Corish et al., 2007; Di Carli et al., 2010; Díaz-Mojica et al., 2014; Mirwald et al., 2019). On the other hand, natural earthquakes, especially large ones, are usually occurred on faults with complex geometries, such as bend, branch, step-over, etc. In this study, we performed dynamic inversion problem for a branched fault. The branched fault model considered in this study is shown in Figure 2. The entire branched fault system consists of three planar faults: the main fault of dimension 30 km×10 km and two bifurcated faults both of dimension 15 km×10 km. Parameters of the medium and others used in forward simulations are shown in Table 1. In our forward modeling, for simplicity and to ensure the calculation efficiency, we assume that the slip is only along the x direction. To initiate rupture, a nucleation zone, called asperity, with slightly higher initial stress than the peak strength is introduced. An unbreakable barrier is set at the edge of the fault.

Table  1.  Modeling parameters used in this study

α (km/s)	β (km/s)	ρ(kg/m3)	$\dfrac{{\alpha \varDelta t}}{{\varDelta s}}$	μs	μd	$\sigma _{{\textit{zz}}}^A\left({{\rm{MPa}}} \right)$	$\tau _{{\textit{z}}x}^A\left({{\rm{MPa}}} \right)$	Δs (m)	Δt (s)	D0 (m)	ϕ (°)
5.6	3.2	3000	0.98	0.6	0.1	100	50	249	0.044	0.57	15

 | Show Table

DownLoad: CSV

Figure  2.  The branched fault model used in this study. The entire branched fault system consists of three planar faults: the main fault of dimension 30 km×10 km and two bifurcated faults both of dimension 15 km×10 km

DownLoad: Full-Size Img PowerPoint

In the GA inversion, many control parameters determine the accuracy and convergence of the algorithm. In this study, we perform a series of inversions, and the basic inversion parameters used in the simulations are listed in Table 2. Among them, N is the number of parameter population, Gen is the number of genetic iterations, l_g is the genetic chain length of the parameter individuals, and p_c and p_m are the hybridization and mutation parameters, respectively. The settings of these inversion parameters are selected by trial and error. The initial parameter population is randomly generated within the searching range.

Table  2.  Inversion model parameters used in this study

GA parameter		Inversion parameter
N	50	Target $T_{\rm{tar}}^0$ (MPa)	47
Gen	50	Target $D_{\rm{c}}^{\rm{tar}}$ (in D0)	0.8
lg	10	Searching range of T0 (MPa)	40−50
pc	0.5	Searching range of $D_{\rm{c}}^{\rm{tar}}$ (in D0)	0.6−1.1
pm	0.008

 | Show Table

DownLoad: CSV

3.   Numerical results and discussion

3.1   Dynamic inversion results

We performed the dynamic inversion of the target parameter model in Table 2 to verify the accuracy and convergence of our method. The number of parameter population and genetic iterations are both 50, which means that we will perform 2,500 forward simulations. Figure 3a shows the misfit function evolution with the genetic iterations, where the abscissa is the generation, and the ordinate is the misfit function. The blue error-bars represent the misfit function of the whole parameter population, and the red point represents the best fit model in each generation. The initial parameter population is randomly selected in the parameter space. Therefore, it can be seen that there is a massive difference between the results of the initial population and those of target. With the progress of genetic iteration, due to selection and crossover, the genes of the parameter individuals with low misfit are more prone to be retained for heredity. At the same time, those with high misfit are gradually eliminated, and the misfit function of the entire parameter population decreases gradually. Meanwhile, the best fit parameter in the population is closer to the target parameter in the iteration. After three generations of iteration, the misfit function of the parameter population and the best fit model began to converge to a specific value.

Figure  3.  Results of the dynamic inversion based on the GA. The blue error-bars represent the whole parameter population, and the red points represent the best fit model in each generation. (a) Misfit function evolution with the genetic iterations; (b) Comparison of inversion parameters T⁰ and target parameter ${T_{\rm{tar}}^0}$; (c) Comparison of inversion parameters D_c and target parameter ${D_{\rm{c}}^{\rm{tar}}}$

DownLoad: Full-Size Img PowerPoint

Figures 3b and 3c show the comparisons of the values of inverted parameters and those of target parameters during the inversion process, where the ordinate is the value of the inversion parameters and the target parameter. For the initial stress T⁰, the initial searching range of our parameter population is 40 to 50 MPa. After iterating for about 26 generations, the average value of T⁰ of the parameter population gradually converges to that of the target model. A non-dimensional critical slip-weakening distance $\overline {{D_{\rm{c}}}} = {{{D_{\rm{c}}}}}/{{{D_0}}}$ instead of D_c itself, in which D₀ is a reference value listed in Table 1, is adopted in the inversion. The initial searching range of $\overline {{D_{\rm{c}}}}$ is [0.6 1.1], and the target parameter is 0.8 D₀. After iterating for about 26 generations, the overall D_c value of the parameter population converges to that of target model. Figure 4 shows the comparison of the final slip distributions between the target model and the best fit model. The inversion results show that the accuracy and speed of convergence of our inversion are satisfactory. The parameter population tends to converge to the target ones after about 26 generations. It can be seen from Figures 3b and 3c that D_c has better convergence and accuracy than T⁰, which indicates that the effect of D_c on the rupture result is more sensitive than that of T⁰. An interesting phenomenon is that the misfit function converged fastly, but the parameters converged after 26 generations. The reason is that a simple Pearson correlation coefficient is used as the misfit function in our study, which is insensitive to small changes in parameters.

Figure  4.  Comparison of the final slip between the target model and the best fit model. The left part of the diagram are the final slip on the main plane and upper branch. In contrast, the right part represent the lower branch. (a) and (b) are the final slip distributions of the target model and the best fit model, respectively; (c) is the slip difference distribution between the target model and the best fit model

DownLoad: Full-Size Img PowerPoint

3.2   Effect of the size of population

In the dynamic inversion, the size of the initial population and the searching range play an essential role. The size of the initial population determines not only the genetic diversity but also the amount of forward computation. A large number of populations result in high stability of genetic reproduction, and easy to generate good genes and inherit it. If the population is too small, the inverted parameters are likely to fall into optimal local solutions and even fail to converge. Figure 5 shows the inversion results for the initial population of 30. It can be seen that the values of inversion parameters converge to local optimal solutions instead of those of target parameters. Figure 6 shows the comparison of the slip time functions on a specific element between the target model and the best fit model. The local optimal solution has a similar rupture result to that of the target model. Since the size of population is small, the possibility of mutations reduces, which makes the parameter samples of the entire population more concentrated.

Figure  5.  Same as Figure 3, except the size of initial parameter population is 30

DownLoad: Full-Size Img PowerPoint

Figure  6.  Comparison of the slip time functions between the target model and the best fit model. The size of initial parameter population is 30. The element is locate on x=30 km, y=5 km

DownLoad: Full-Size Img PowerPoint

Gene chain is the core of the GA. The essence of the genetic process is the selection of genes, hybridization and mutation. In our inversion algorithm, we represent the individual parameter values with a binary gene chain. Therefore, the searching range of the parameter and the length of the gene chain will jointly determine the accuracy of the parameter. Figure 7 shows the inversion results with the gene chain of length 12, while the length of gene chain is 10 in Figure 3. A longer gene chain results in that the daughter population was more similar to the that of the parents’, and reduced effects of mutations. So the inversion quickly converges to the optimal local solution in about 15 generation. Figure 8 shows the comparison of the final slip distributions between the target model and the best fit model. The rupture result of the best fit model has only little difference to the target model, which has a very high correlation with the target model. An interesting phenomenon is that the inversion parameters seem to converge to the target parameters at the beginning of the iteration, but gradually deviate and converge to the local optimal solution. This is because in the early generations, the genes of the entire population are still scattered. Although the average value is similar to the target parameter, there are still a lot of bad genes in the population.

Figure  7.  Same as Figure 3, except the length of the gene chain is 12

DownLoad: Full-Size Img PowerPoint

Figure  8.  Same as Figure 4, except the length of the gene chain is 12 here

DownLoad: Full-Size Img PowerPoint

3.3   Effect of the genetic parameters

After the forward simulation, the inversion process mainly includes three steps, namely selection, crossover and mutation (see Figure 1). In the selection process, the Pearson correlation coefficients are used to characterize the rupture results correlation between the individual parameter and the target one. At the same time, in order to highlight the excellent parameter individuals, the correlation coefficient is magnified to make the inversion converges to a more precise value. We use the correlation to assign selection probability of the parameter population and adopt the roulette model to select the parameters of the parents. The good genes with higher correlation have higher probability to be chosen as parents and to inherit their genes.

After selecting the parents of the parametric population, crossover and mutation can be used to obtain a new parameter population of the next generation. In Table 2, the crossover and mutation parameters, p_c and p_m, used in this study are listed. For the selected parents, we hybridize them on a statistical probability model, and p_c represents the probability of crossover between the parents. This allows some parents to directly retain their genes in the next generation, while others will generate new individuals by crossover. In the genetic iteration, the larger the p_c value is, the more difficult it is to keep the genes of the parents during the crossover process, and there will be more new individuals in each generation. The large p_c will make the result more likely to reach an accurate value, but at the same time, there is a risk of non-convergence. The smaller the p_c value is, the more likely to retain for the parents at the next generation, which means that the population will quickly convergence but maybe reach an optimal local solution. Figure 9 shows the inversion results of p_c=0.3. The parameter population started to converge in about ten generations, but converged to the optimal local solution, which differed significantly from the target parameter. Figure 10 shows the comparison of the rupture processes between the target model and the best fit model. The rupture slip rate difference are small on the main fault but large on bifurcation faults. Similar to Figure 7, it can still be seen here that the population seems to converge to the target parameters at the beginning of the iteration. The initial best fit model, very close to the target model, did not fully inherit its genes, thereby failing to prevent the entire population from converging towards a local optimal solution.

Figure  9.  Same as Figure 3, except the crossover parameter p_c is 0.3

DownLoad: Full-Size Img PowerPoint

Figure  10.  Comparison of the rupture processes between the target model and the best fit model. The crossover parameter p_c is 0.3. (a) The rupture slip rate of the target model. (b) The rupture slip rate difference between the target model and the best fit model

DownLoad: Full-Size Img PowerPoint

The mutation parameter p_m represents the probability of mutation in the gene chain after the crossover process. The mutation represents that a genetic variation is generated, which means some of the gene fragments of the daughter parameters are not the same as those of the parents. The occurrence of genetic variation in heredity is minimal, but it has an essential influence on the genetic process. It makes it possible for any genes to appear in the parameter population during each generation, which reduces the possibility of converging to a locally optimal solution in the inversion. The larger the p_m value is, the higher probability of gene mutation will be, and the more good genes are likely to occur in the genetic process. However, too many mutations will make the good genes of the parents less retained, which may cause non-convergence. The smaller the p_m value is, the lower the probability of a gene mutation will be, the easier it is to converge to a locally optimal solution during the genetic process. Figure 11 shows the inversion results of p_m=0.015. Due to the high probability of mutations in the entire population, its parameter range is always covered widely. And parent’s gene retention was reduced, and the entire parameter population failed to converge.

Figure  11.  Same as Figure 3, except the mutation parameter p_m is 0.015

DownLoad: Full-Size Img PowerPoint

4.   Conclusions

In this study, we introduced an inversion method to estimate the source stress and friction parameters. The inversion method is based on the GA, and the forward problem is carried out using the BIEM with the GPU parallel, which can simultaneously solve thousands of computationally expensive forward simulations. The source model is a branched fault, where a slip-weakening law governs the rupture. We inverted the stress field parameter T⁰ and friction parameter D_c in the rupture process. The inversions are based on theoretical tests by using the target parameters’ forward process result as an input to invert the target parameters. The results show that accurate results can be obtained based on our inversion method, and the parameter population often converges to the target parameter around several generations iteration. The inherited parameters have a significant impact on the inversion results choosing a suitable inherited parameter can significantly improve the efficiency and convergence of the inversion. During the inversion process, we found that the inversion result is more sensitive to D_c, which tends to converge to the target value quickly and accurately. In contrast, the stress field parameters converge slowly and prone to locally optimal solutions. This means that a small change in D_c may affect the result more than a small change in T⁰. We think that this may be caused by the research range of D_c and T⁰ in our inversion. The research range of T⁰ is large so that the slight changes will result not obvious changes.

Acknowledgments

This study is supported by the National Natural Science Foundation of China (No. 41874047) and the High-performance Computing Platform of Peking University.