Earthquake Research in China  2019, Vol. 33 Issue (3): 391-402     DOI: 10.19743/j.cnki.0891-4176.201903004
The Improved Simulated Annealing Algorithm and Application in the Surface Wave Dispersion Inversion
ZHANG Xuemei1, YANG Zhigao1, SHI Haixia1, DU Guangbao1, YANG Wen1,2, WEI Xing1, HAN Yanyan1, LI Fei1, HUANG Zhibin1, LIU Jie1
1. China Earthquake Networks Center, Beijing 100045, China;
2. Institute of Geophysics, China Earthquake Administration, Beijing 100081, China
Abstract: The velocity distribution of layers from surface wave dispersion curve is a severely non-linear program. Base on the Metropolis rule, we improved the simulated annealing algorithm to simultaneously inverse the velocities and thicknesses using the dispersion data and identified the Moho and the bottom of lithosphere. The application to the numerical examples with 5% noise shows the velocity RMS is 1.56% between the non-linear results and the original models when the condition of selecting method for temperature parameters and initial temperature are satisfied. Using the pure dispersions of Rayleigh wave, the non-linear inversion has been carried out for S-wave velocities and thicknesses of the vertical profile crossing the Indian Plate, the Qinghai-Tibetan Plateau, and the Tarim Basin. It indicated that the crustal thickness is about 70km in the Qiangtang block, while in the hinterland of the Qinghai-Tibetan Plateau the lithosphere is relatively thin (~130km) from the velocity values and their offsets.
Key words: Simulated annealing algorithm     Simultaneous inversion     Dispersion curve     Reconstruction of velocity model

INTRODUCTION

Using the observed data of seismic body wave, seismic surface wave and receiver function, especially the dispersion and attenuation characteristics of seismic surface wave, we can retrieve the velocity and thickness distribution of the deep structure of the earth and provide deep tectonic background for epicenter location, earthquake prediction, metal minerals, and oil and gas resources exploration. In depth inversion, the velocity structure of the crust and upper mantle is usually solved by linear or approximate linear inversion (Chen Guoying et al., 1985; Sun Kezhong et al., 1985; Song Zhonghe et al., 1991; Teng Jiwen et al., 1994, 1997; Hirn A. et al., 1997; Griot D. A. et al., 1998; Wu Jianping et al., 1998; Yanovskaya T. B. et al., 1990; Ding Zhifeng et al., 2001; Su Wei et al., 2002; Zhu Jieshou et al., 2004; Sun Ruomei et al., 2004; Zhong Yuyun et al., 2011; Ye Xiuwei et al., 2016). Although we can use the linear inversion method to improve the instability of problem solving and quickly solve the optimal solution of mathematical and physical equations, there are large linear errors (Dorren H. J. S., 1965). When using approximate linear algorithm to invert the variation of the velocity with depth, the calculation error of numerical derivation will lead to non-convergence and failure of iteration, and the data error will also lead to iteration divergence (Yang Wencai, 2002). When linear inversion is performed by least squares method, gradient descent method and singular value decomposition method, the partial differential or Jacobian matrix of the model is computed, and the local optimal solution is obtained through several iterations, which depends heavily on the selection of the initial model.

The relationship between model space and geophysical field space is usually non-linear. For example, there is a highly non-linear relationship between seismic surface wave dispersion curve and velocity structure of the earth. Linear inversion algorithm is not suitable for solving non-linear problems. In recent years, with the improvement of computational methods and the realization of supercomputing, many non-linear inversion methods based on Monte Carlo algorithm have been widely used in earth exploration (Yang Wencai, 1993; Yao Zhenxing et al., 1999). In the 1980s, Rothman D.H., (1986) solved the estimation of large static correction in earthquake exploration by simulated annealing method. After that, more and more geophysicists pay attention to the non-linear inversion method. The simulated annealing algorithm searches the model space without solving the partial derivative of the objective function, and chooses the appropriate initial temperature and attenuation mode to obtain the global optimal solution. This method has been widely used in seismic tomography, geophysical exploration and epicenter location etc. (Basu A. et al., 1990; Press W. H. et al., 1991; Liu Pengcheng et al., 1995; Yao Yao, 1995; Wang Shanshan et al., 1995; Ji Chen et al., 1996; Zhang Linbin et al., 1997; Martınez M. D. et al., 2000; Jing Xili et al., 2003; Zhang Hongmei et al., 2004; Ryden N. et al., 2006; Wei Chao et al., 2007; Lu Pengfei et al., 2008; Chen Xiao et al., 2010; Gao Xing et al., 2010; Su Jinrong et al., 2013; Wang Yang et al., 2015; Sun Li, 2016; Peng Guomin et al., 2018). In inversion of deep structures by surface wave dispersion, a model discretization method is usually used to determine the stratification velocity by fixing the thickness of the stratum, which results in strong uncertainty of the depth. In order to improve the recognition ability of the interface of deep structure, we use the improved simulated annealing algorithm to simultaneously invert the S-wave velocity and layers of the crust and upper mantle.

1 IMPROVEMENT OF SIMULATED ANNEALING ALGORITHMS

Simulated Annealing Algorithms (SAA) is a heuristic Monte Carlo inversion method based on statistical mechanics to simulate the annealing process of thermal equilibrium systems (Yao Yao, 1995). The solution i of the mathematical physics problem and the objective function f(i) correspond to the microscopic state and the energy Ei of the solid state system. Given the initial model of the mathematical equation, the temperature parameter T is gradually reduced, and the objective function f(i) is calculated. When the control parameter T approaches the minimum value, the global optimal solution of the mathematical equation is obtained. The solid state system undergoes a slow annealing and cooling process, so that the solid reaches a thermal equilibrium state at each temperature and finally approaches the minimum ground state of energy. According to the similarity between the optimization problem and the object annealing process, we slowly reduce the temperature T of the control parameter, so that the algorithm reaches a steady state, and the energy of the corresponding system tends to be the smallest.

Monte Carlo algorithm can be used to simulate the process of approaching thermal equilibrium, but the stochastic search algorithm has a huge amount of computation. Based on Metropolis algorithm (Metropolis N. et al., 1953), the temperature attenuation function is used to control the temperature declining process for simulated annealing of solid state systems. The initial state i of the particle's position is taken as the current state of the system. Ei is used to indicate its energy. In the new physical state j with the random perturbation, its energy is Ej. If the energy Ej < energy Ei, then the state j is regarded as the important state of the system. If the energy Ej > energy Ei, which is saffected by the thermal motion, whether the state j is an important state of the system depends on the probability of the state. Metropolis criterion is used as acceptance criterion for new states. When the temperature of a physical system is higher, we define the state with larger energy difference as an important state; when the temperature decreases, even the state with smaller energy difference can be also regarded as an important state when the temperature reaches a lower one. The way and coefficient of temperature attenuation are the key parameters, which directly affect the inversion speed and results. Too fast temperature attenuation may fall into the local region of minimum solution, and too slow temperature decreasing will affect the convergence effect. In order to avoid falling into the local minimum model space when solving the equation, we need to carefully select attenuation mode, attenuation coefficient and initial state in the geophysical inversion process.

In recent years, the simulated annealing nonlinear algorithm has made great progress. Basu A. et al., (1990) used experimental methods to determine the critical temperature T of simulated annealing when recovering the acoustic profile from the wavefield data. Wang Shanshan et al. (1995) used the Gibbs distribution and the Cauchy distribution to solve the optimal solution of the equation iteratively. Yao Yao (1995) used fuzzy prior information to determine the minimum temperature and reconstruct the objective function in order to avoid the local minimum. Jing Xili et al. (2003) proposed an self-adaptive simulated annealing algorithm by analyzing the crystal growth process, which overcomes the defects of repeated experiment selection parameters. Zhang Linbin et al. (1997) combined the generalized Gibbs distribution based on multidimensional fractal theory with the temperature-dependent Cauchy distribution to generate a new perturbation model, which performs large-scale search in high temperature conditions, and the low temperature state searches in the neighborhood, jumping out of the local minimum solution space. Wei Chao et al. (2007) introduced quantum annealing nonlinear methods based on quantum transition tunneling effect to invert the geophysical parameters. Lu Pengfei et al. (2008) calculated the candidate solutions of the approximate Cauchy distribution at the same temperature in the prestack reservoir parameter inversion, jumping out of the local minimum region and improving the simulated annealing inversion speed. Wang Jicheng et al. (2017) performed microseismic location based on simulated annealing and grid dividing. The combination of simulated annealing and other inversion methods is also applied step by step when solving the global optimal solution of the objective function (Press W. H. et al., 1991; Liu Pengcheng et al., 1995; Zhang Hongmei et al., 2004; Li Fei et al., 2017). With the improvement of simulated annealing algorithm, the convergence speed of calculation is accelerated in varying degrees, and the inversion results are more robust.

Whether linear inversion or non-linear inversion, we usually only solve one class of physical parameter. When inverting the velocity structure in each depth of the study area based on the surface wave dispersion of each cell, we usually fix the thickness of each layer and invert it to get the velocity of different layers. In order to accurately obtain the thickness distribution of the crust and lithosphere in the lithosphere-asthenosphere system, we improved the simulated annealing non-linear method to simultaneously invert the velocity and interface of deep structure base on the previous study. According to the partial derivatives of the group velocity with the respect to the velocity, depth and density, we set shear wave velocity as the main physical parameter. In the process of reducing temperature simulated annealing, we use variable-scale cooling function to choose larger-scale cooling mode at higher temperature, and then use smaller-scale cooling mode to slow down to obtain the optimal solution when reaching relatively stable state.

For the simulated annealing nonlinear inversion, we use the T(k)=T0αik cooling function, where T(k) is the temperature of k times, and T0 is initial temperature value, and αi is the coefficient of cooling, generally αi∈[0.7, 1].

The calculation steps are as follows:

(1) Let [vo, Do] be the initial model; k=1; To=Tmax (initial maximum temperature); where vo is the velocity vector of the initial model and Do is the layer thickness vector of the initial model;

(2) If the cycle stop condition is satisfied at this temperature, go to step (3); otherwise, search for [vj, Do] from the model neighborhood N([vi, Do]), vj=vi+(γT(k)((1+T(k)-1)(2random(0, 1))-1)) (vmax-vmin) where $\gamma=\left\{\begin{array}{ll}{1} & {\text { random }(0, 1)-0.5 \geqslant 0} \\ {-1} & {\text { random }(0, 1)-0.5<0}\end{array}\right.$ vmax is the maximum value of the velocity model, and vmin is the minimum value of the velocity model.

Calculate the objective function ΔYij=f ([vj, Do])-f([vi, Do]), where f is the dispersion operator; if ΔYij≤0, set the new model [vi, Do]=[vj, Do], otherwise, if exp(-ΔYij/T(k))>random(0, 1), set the new model [vi, Do]=[vj, Do]; repeat step (2);

(3) Decrease the temperature Tk=T0α1k; k=k+1; if a relatively low temperature is reached, go to step (4); otherwise, return to step (2);

(4) If the energy reaches the minimum condition at this temperature, go to step (5); otherwise, search for the new model [vjj, Djj] from the model neighborhood N([vii, Dii]),

vjj=vii +(γT(k)((1+T(k)-1)(2random(0, 1))-1))(vmax-vmin),

Djj=Dii+(γT(k)((1+T(k)-1)(2random(0, 1))-1)) (Dmax-Dmins),

where Dmax is the maximum thickness and Dmin is the minimum thickness;

Calculate the objective function ΔYij=f([vjj, Djj])-f([vii, Dii]); if ΔYij < =0, set the new model [vii, Dii]=[vjj, Djj], if exp(-ΔYij/T(k))>random(0, 1), set new model [vii, Dii]=[vjj, Djj]; repeat step (4);

(5) Decrease the temperature Tk=Tk-1α2k; k=k+1; if the loop condition is satisfied, output the inversion model [vii, Dii]; otherwise, return to step (4).

2 NUMERICAL EXPERIMENTS

A 7-layer velocity model of the crust and upper mantle with 430km thickness is set as the initial model (see Fig. 1). Based on the surface wave dispersion equation derived by Haskell's matrix (Haskell N.A., 1953), the form of the layer matrix improved by Schwab F. et al., (1970), Nur A. et al., (1969) and Panza G.F., (1985), the density ρ and P-wave velocity α are obtained from the thickness h and S-wave velocity β according to the empirical formula, then the theoretical dispersion curve of the velocity model is calculated. The group velocities in different periods have different partial derivations for the parameters such as S-wave velocity β, P-wave velocity α and density ρ (Zhu Jieshou, 1988). The ∂U/∂ρ and ∂U/∂α of the Rayleigh wave group velocity are relatively small, while the value of ∂U/∂β is relatively large. In the process of nonlinear inversion, thickness h and S-wave velocity β are taken as independent variables, while P-wave velocity α and other parameters are taken as independent variables, which are derived from S-wave velocity β and the depth.

 Fig. 1 (a) The dispersion curves for the theoretical model and the accepted model obtained with the SAA inversion and SVD inversion; (b) The accepted model obtained with the SAA inversion and SVD inversion

In Fig. 1 the velocity structure (solid line) is used as a theoretical model for forward calculation, and its dispersion curve is calculated (Fig. 1(a)). Taking the dispersion curve (solid black line) in Fig. 1(a) as measured dispersion data, the modified simulated annealing algorithm is used to simultaneously invert the distribution of the S-wave with depth up to 430km (Fig. 1(b)). There is an error in the observed dispersion curve (Zhang Xuemei et al., 2014; Chen Haopeng et al., 2014). We added 5% noise and then inverted with Singular Value Decomposition algorithm (Golub G.H. et al., 1970) and improved simulated annealing algorithm. From the improved simulated annealing nonlinear inversion, when the fitting degree of dispersion curve is 98.33%, the root mean square error between layer velocity and layer velocity of a given model is 1.56% and 1.78% for the linear inversion when the fitting degree of dispersion curve is 97.67%. The improved simulated annealing non-linear inversion results show that the plane position is close to the initial position of the theoretical model. It can be seen from the numerical simulation that the non-linear inversion algorithm can fit the dispersion curve well and get the approximate theoretical velocity model.

3 APPLICATION: SURFACE WAVE DISPERSION INVERSION IN THE WESTERN QINGHAI-TIBETAN PLATEAU

Studying the spatial tectonic framework of the lithosphere-shoal system is important for improving the theory of plate tectonics and exploring the coupling of layers and the deep dynamic processes of the Earth (Teng Jiwen et al., 1997; Teng Jiwen, 2003). The Qinghai-Tibetan Plateau, which is located in the place of the continent-continent collision between the Indian and Eurasian Plates, has a very complex deep and shallow tectonic pattern. Since the beginning of the deep exploration of the Qaidam Basin in 1958, international cooperation and geophysicists have carried out many studies of the crust and upper mantle in the Qinghai-Tibetan Plateauand achieved some important results (Chen Guoying et al., 1985; Sun Kezhong et al., 1985; Zhao Wenjin et al., 1993; Teng Jiwen et al., 1994; McNamara D. E. et al., 1994; Wittlinger G. et al., 1996; Hirn A. et al., 1997; Griot D.A. et al, 1998; Wu Jianping et al., 1998; Zeng Rongsheng et al., 2000; Ding Zhifeng et al., 2001; Su Wei et al., 2002; Sun Ruomei et al., 2004; Li Yonghua et al., 2006; Xu Zhiqin et al., 2013; Deng Qidong et al., 2014; Huang Zhouchuan et al., 2015; Zhao Wenjin, 2015; Gilligan A. et al., 2015; Zhao Junmeng et al., 2016; Zheng Chen et al., 2016; Yang Xiaoping et al., 2016; Lü Zhiqiang et al., 2016; Teng Jiwen et al., 2017; Zhang Fengxue et al., 2018). These studies pay more attention to the material properties of the deep earth, and relatively few studies on the interface of the crust, lithosphere and asthenosphere. We collect long period seismic data in the Qinghai-Tibetan Plateau and its adjacent areas and obtain the Rayleigh wave group velocity of 8-150s period by surface tomography. The cells along the 86°E cross the hinterland of the Qinghai-Tibetan Plateau, the Tarim Basin and the northern part of the Indian Plate are selected to be inverted. The profile (Fig. 2) explores the thickness of the crust, the lithospheric thickness and seismicity.

 Fig. 2 Sketch map of regional tectonics which the profile AA' is located. The blue line is the profile AA'. LS: Lhasa block; QT: Qiangtang block; BK: Bayan Har block; QDM: Qaidam Basin; TRM: Tarim Basin; MFT: Main Frontal Thrust; MBT: Main Boundary Thrust; MCT: Main Central Thrust; YZS: Yarlung Zangbo suture belt; BNS: Bangong-Nujiang suture belt; JSJS: Jinshajiang suture belt; ATF: Altyn Tagh fault

Based on the long period surface wave data of the China Seismic Network and the Global Seismic Network, after pre-processing, we applied the Yanovskaya-Ditmar tomography (Yanovskaya T. B. et al., 1990) to inverse the surface wave group velocity at the period of 8-150s in the Qinghai-Tibetan Plateau and adjacent areas, and the dispersion data set of the 2°×2°size cell was discretized. The S-wave velocity structure obtained by linear inversion is used as the initial model (Sun Ruomei et al., 2004). The improved simulated annealing non-linear algorithm is used to simutanously inverse the velocity and thickness for each cell. The velocity distribution of S-wave across the North-South section is obtained (Fig. 3). Compared with the linear inversion results, the velocity distribution obtained by the two-parameter non-linear inversion method can relatively clearly determine the velocity offsets of Moho and lithospheric bottom. On the vertical velocity distribution map, we take the S-wave velocity increasing from ~3.7km/s to ~4.3km/s as the Moho discontinuity; the depth where S-wave velocity begins to reduce is considered as the bottom of the lithosphere; and the low velocity zone below the lithosphere is the asthenosphere. From S-wave velocity and offsets along 86°E section, the Moho discontinuity and the lithospheric bottom of the section are obtained (see Fig. 4). The data used in surface wave tomography in this paper have a period of 150s and can be used to invert the velocity structure up to 400km. However, the sensitive core of velocity at 300-400km depth is wider, and the bottom of asthenosphere can only be used as a reference.

 Fig. 3 The accepted models obtained with the SAA inversion along the profile 86°E

 Fig. 4 The S-wave velocity distribution and the crustal and the lithospheric thickness along the profile 86°E. The black dash lines denote the Moho discontinuity and the lithospheric interface. The dark blue circles are the earthquakes (M≥3.0) near the profile AA'

It can be seen from the Moho along the vertical section of the 86°E and the bottom of the lithosphere (Fig. 4) that the Moho of the Ganges River drainage area of the Indian Plate (ID) is relatively shallow, between 35-40km, and the lithosphere is relatively thick, up to 185km with thin asthenosphere, and high speed upper mantle lid. The crust of Himalayan block (HM) to the north of the main boundary suture belt (MBT) is thicker than it in the ID (about 60km), the lithospheric thickness is about 155km, and the low velocity zone below the bottom is gradually thickened to the north. The Lhasa block (LS) is located between the Bangong-Nujiang suture belt (BNS) and the Yarlung Zangbo suture belt (YZS). The depth of Moho is 60-70km and a thinner lithosphere is about 160km. The thickest crust is located in the Qiangtang block (QT) in the hinterland of the Qinghai-Tibetan Plateau, with a thickness of up to 75km, while the thickness of the lithosphere is the thinnest (about 130km) of the entire section, and the asthenosphere, the low velocity layer beneath it, is well developed. In the QT, the velocity of the Moho discontinuity does not change significantly. In this area, Sn is abruptly attenuated (Ni J. et al., 1983). The Cenozoic volcanic rocks are rich in potassium, strontium and neodymium (Deng Wanming et al., 1997). We speculate that there is partial melting of the crust-mantle transition zone in the deep structure. The thickness of the crust of the Bayan Har block (BK) in the north of the QT is relatively thin (about 65km), and the thickness of the lithosphere is between 150-160km. Further northward through the Altyn Tagh fault (ATF), it enters the Tarim Basin (TRM). The thickness of the crust is thinner than the hinterland of the Qinghai-Tibetan Plateau, about 50km, while the bottom of the lithosphere is deeper, reaching 185km.

Along the 86°E section (AA'), we collected earthquakes that occurred in the vicinity of the section (±1.0°)with magnitudes greater than or equal to 3.0 since 1930. It can be seen from Fig. 4 that most of the earthquakes in this area occur in the crust above the Moho discontinuity. However, there is a depth difference of the earthquakes in the north and south of the Qinghai-Tibetan Plateau. In the QT, the Bayan Har block (BK) and the TRM, the earthquakes occurred in the crust and did not occur in the upper mantlelid. In the south of the BNS, the focal depth is relatively deep; there are earthquakes in the crust and earthquakes that pass through the Moho discontinuity and reach the upper mantle lid. Especially in the Himalayan block (HM), south of the YZS, more deep-source earthquakes (depth greater than 70km) have occured in the crust and the upper mantle lid.

4 DISCUSSION AND CONCLUSIONS

Using the improved simulated annealing non-linear algorithm to simultaneously invert the S-wave velocity structure along the 86°E section, the following deductions are obtained:

(1) Applying the improved simulated annealing algorithm, the variable-scale model changes are adopted in the process of temperature declining, and the velocity and thickness are inverted simultaneously. The results of numerical simulation are similar to those of the initial models.

(2) The improved simulated annealing algorithm still has the problem of multiple solutions in the inversion algorithm. When the shear wave velocity and its depth of each block are inverted simultaneously base on pure dispersion, there is a coupling between the depth variable and the velocity variable. In this paper, we calculate the average value of the solution in a certain error range as the optimal solution. The initial model should be constrained by geodesy, geology, geophysics and other information to gain more reasonable solutions similar to real deep structure of the earth.

(3) The Moho discontinuity and the lithospheric bottom of the 86°E section can be obtained clearly by simultaneously inversion of the thickness and velocity of each layer along the profile with the improved simulated annealing algorithm. The hinterland of the Qinghai-Tibetan Plateau has a thick crust (about 70km) and a relatively thin lithosphere (about 130km). The earthquakes distribution along this section show that there are differences between the north and the south of the Qinghai-Tibetan Plateau. In the QT, the BK and the TRM, earthquakes occur in the crust. To the south of the BNS, the focal depth is relatively deep, especially in the HM, and more earthquakes occur in the crust and the upper mantle lid.

ACKNOWLEDGEMENT

The China Earthquake Networks Center (CENC) and the International Seismological Center (ISC) provided seismic station observations for this study, and Yanovskaya provided the procedures for surface tomography, the authors would like to express their heartfelt thanks to them. The authors also thank the anonymous reviewers for their constructive suggestions.

REFERENCES
 Basu A., Frazer L.N.Basu A., Frazer L.N. Rapid determination of the critical temperature in simulated annealing inversion[J]. Science, 1990, 249(4975): 1409-1412. DOI:10.1126/science.249.4975.1409 Chen Guoying, Zeng RongshengChen Guoying, Zeng Rongsheng. The difference of lithospheric structure between Himalayan mountain and Tibetan Plateau from surface wave dispersion[J]. Acta Geophysica Sinica, 1985, 28(S1): 161-173 (in Chinese). Chen Haopeng, Zhu Liangbao, Wang Qingdong, Yang Yinghang, Zhang PanChen Haopeng, Zhu Liangbao, Wang Qingdong, Yang Yinghang, Zhang Pan. An error analysis of surface wave group velocity measurement[J]. Progress in Geophysics, 2014, 29(5): 1985-1995 (in Chinese with English abstract). Chen Xiao, Yu Peng, Zhang Luolei, Li Yang, Wang JialinChen Xiao, Yu Peng, Zhang Luolei, Li Yang, Wang Jialin. Regularized synchronous joint inversion of MT and seismic data[J]. Seismology and Geology, 2010, 32(3): 402-408. Deng Qidong, Cheng Shaoping, Ma Jin, Du PengDeng Qidong, Cheng Shaoping, Ma Jin, Du Peng. Seismic activities and earthquake potential in the Tibetan Plateau[J]. Chinese Journal of Geophysics, 2014, 57(5): 678-697. DOI:10.1002/cjg2.20133 Deng Wanming, Zhong DalaiDeng Wanming, Zhong Dalai. Crust-mantle transition zone and its geological significance in the lithospheric tectonic evolution[J]. Chinese Science Bulletin, 1997, 42(23): 2474-2482 (in Chinese). DOI:10.1360/csb1997-42-23-2474 Ding Zhifeng, He Zhengqin, Wu Jianping, Sun WeiguoDing Zhifeng, He Zhengqin, Wu Jianping, Sun Weiguo. Research on the 3-D seismic velocity structures in Qinghai-Xizang Plateau[J]. Earthquake Research in China, 2001, 17(2): 202-209 (in Chinese with English abstract). Dorren H. J. S. Estimates of Stability and Propagation of Errors in Nonlinear Inverse Problems[D]. Dutch: Utrecht University, 1965. Gao Xing, Duan Zongqi, Wang Weimin, Guo ZhiGao Xing, Duan Zongqi, Wang Weimin, Guo Zhi. Rupture process of the Oct. 6, 2008 MS6.6 Damxung earthquake, Tibet, China[J]. Chinese Journal of Geophysics, 2010, 53(9): 2083-2090 (in Chinese with English abstract). Gilligan A., Priestley K.F., Roecker S.W., Levin V., Rai S. S.Gilligan A., Priestley K.F., Roecker S.W., Levin V., Rai S. S. The crustal structure of the western Himalayas and Tibet[J]. Journal of Geophysical Research: Solid Earth, 2015, 120(5): 3946-3964. DOI:10.1002/2015JB011891 Golub G.H., Reinsch C.Golub G.H., Reinsch C. Singular value decomposition and least squares solutions[J]. Numerische Mathematik, 1970, 14(5): 403-420. DOI:10.1007/BF02163027 Griot D. A., Montagner J. P., Tapponnier P.Griot D. A., Montagner J. P., Tapponnier P. Phase velocity structure from Rayleigh and Love waves in Tibet and its neighboring regions[J]. Journal of Geophysical Research: Solid Earth, 1998, 103(B9): 21215-21232. DOI:10.1029/98JB00953 Haskell N. A.Haskell N. A. The dispersion of surface waves on multilayered media[J]. Bulletin of the Seismological Society of America, 1953, 43(1): 17-34. Hirn A., Sapin M., Lépine J. C., Diaz J., Mei JiangHirn A., Sapin M., Lépine J. C., Diaz J., Mei Jiang. Increase in melt fraction along a south-north traverse below the Tibetan Plateau: evidence from seismology[J]. Tectonophysics, 1997, 273(1/2): 17-30. Huang Zhouchuan, Wang Pan, Xu Mingjie, Wang Liangshu, Ding Zhifeng, Wu Yan, Xu Mijian, Mi Ning, Yu Dayong, Li HuaHuang Zhouchuan, Wang Pan, Xu Mingjie, Wang Liangshu, Ding Zhifeng, Wu Yan, Xu Mijian, Mi Ning, Yu Dayong, Li Hua. Mantle structure and dynamics beneath SE Tibet revealed by new seismic images[J]. Earth and Planetary Science Letters, 2015, 411: 100-111. DOI:10.1016/j.epsl.2014.11.040 Ji Chen, Yao ZhenxingJi Chen, Yao Zhenxing. The uniform design optimized method for geophysics inversion problem[J]. Acta Geophysica Sinica, 1996, 39(2): 233-242 (in Chinese with English abstract). Jing Xili, Yang Changchun, Li LiJing Xili, Yang Changchun, Li Li. Adapted simulated annealing algorithm for geophysical problems[J]. Progress in Geophysics, 2003, 18(2): 237-330. Li Fei, Zhang Xuemei, Chen Hongfeng, He ShaolinLi Fei, Zhang Xuemei, Chen Hongfeng, He Shaolin. Simiulated annealing and its application in seismic traveltime inversion of 3D velocity model[J]. Earthquake Research in China, 2017, 33(3): 355-364. Li Yonghua, Tian Xiaobo, Wu Qingju, Zeng Rongsheng Zhang RuiqingLi Yonghua, Tian Xiaobo, Wu Qingju, Zeng Rongsheng Zhang Ruiqing. The Poisson ratio and crustal structure of the central Qinghai-Xizang inferred from INDEPTH-Ⅲ teleseismic waveforms: Geological and geophysical implications[J]. Chinese Journal of Geophysics, 2006, 49(4): 1037-1044 (in Chinese with English abstract). Liu Pengcheng, Ji Chen, Hartzell S. H.Liu Pengcheng, Ji Chen, Hartzell S. H. An improved simulated annealing-downhill simplex hybrid global inverse algorithm[J]. Acta Geophysica Sinica, 1995, 38(2): 199-205. Lu Pengfei, Yang Changchun, Guo Aihua, Wang Zhenli, Liu PengLu Pengfei, Yang Changchun, Guo Aihua, Wang Zhenli, Liu Peng. Modified simulated annealing algorithm and its application in pre-stack inversion of reservoir parameters[J]. Progress in Geophysics, 2008, 23(1): 104-109. Lü Zhiqiang, Lei JiansheLü Zhiqiang, Lei Jianshe. 3-D S-wave velocity structure around the 2015 MS8.1 Nepal earthquake source areas and strong earthquake mechanism[J]. Chinese Journal of Geophysics, 2016, 59(12): 4529-4543 (in Chinese with English abstract). Martınez M.D., Lana X., Olarte J., Badal J., Canas J. A.Martınez M.D., Lana X., Olarte J., Badal J., Canas J. A. Inversion of Rayleigh wave phase and group velocities by simulated annealing[J]. Physics of the Earth and Planetary Interiors, 2000, 122(1/2): 3-17. McNamara D.E., Owens T.J., Silver P.G., Wu F.T.McNamara D.E., Owens T.J., Silver P.G., Wu F.T. Shear wave anisotropy beneath the Tibetan Plateau[J]. Journal of Geophysics Research: Solid Earth, 1994, 99(B7): 13655-13665. DOI:10.1029/93JB03406 Metropolis N., Rosenbluth A. W., Rosenbluth M.N., Teller A.H.Metropolis N., Rosenbluth A. W., Rosenbluth M.N., Teller A.H. Equation of state calculations by fast computing machines[J]. Journal of Chemical Physics, 1953, 21(6): 1087-1092. DOI:10.1063/1.1699114 Ni J., Barazangi M.Ni J., Barazangi M. High-frequency seismic wave propagation beneath the Indian Shield, Himalayan Arc, Tibetan Plateau and surrounding regions: high uppermost mantle velocities and efficient Sn propagation beneath Tibet[J]. Geophysics Journal of the Royal Astronomical Society, 1983, 72(3): 665-689. DOI:10.1111/j.1365-246X.1983.tb02826.x Nur A., Simmons G.Nur A., Simmons G. The effect of saturation on velocity in low porosity rocks[J]. Earth and Planetary Science Letters, 1969, 7(2): 183-193. DOI:10.1016/0012-821X(69)90035-1 Panza G. F., Simmons G.Panza G. F., Simmons G. Synthetic Seismograms: the Rayleigh waves modal summation[J]. Journal of Geophysics, 1985, 58: 125-145. Peng Guomin, Xu Kaijun, Du Runlin, Liu ZhanPeng Guomin, Xu Kaijun, Du Runlin, Liu Zhan. Reservoir petrophysical parameter estimation with joint inversion of MCSEM and seismic AVA data[J]. Oil Geophysical Prospecting, 2018, 53(5): 1110-1116 (in Chinese with English abstract). Press W. H., Teukolsky S. A.Press W. H., Teukolsky S. A. Simulated annealing optimization over continuous spaces[J]. Computers in Physics, 1991, 5(4): 426-429. DOI:10.1063/1.4823002 Rothman D.H.Rothman D.H. Automatic estimation of large residual statics corrections[J]. Geophysics, 1986, 51(2): 332-346. DOI:10.1190/1.1442092 Ryden N., Park C. B.Ryden N., Park C. B. Fast simulated annealing inversion of surface waves on pavement using phase-velocity spectra[J]. Geophysics, 2006, 71(4): R49-R58. DOI:10.1190/1.2204964 Schwab F., Knopoff L.Schwab F., Knopoff L. Surface-wave dispersion computations[J]. Bulletin of the Seismological Society of America, 1970, 60(2): 321-344. Song Zhonghe, An Changqiang, Chen Guoying, Chen Lihua, Zhuang Zhen, Fu Zhuwu, Lü Zhiling, Hu JiafuSong Zhonghe, An Changqiang, Chen Guoying, Chen Lihua, Zhuang Zhen, Fu Zhuwu, Lü Zhiling, Hu Jiafu. Study on 3D velocity structure and anisotropy beneath the west China from the Love wave dispersion[J]. Acta Geophysica Sinica, 1991, 34(6): 694-707 (in Chinese with English abstract). Su Jinrong, Guo ZhiSu Jinrong, Guo Zhi. A study of focal mechanisms of the MS6.7 earthquake of 5 October 2008 in southwestern Tianshan-Pamier[J]. Chinese Journal of Geophysics, 2013, 56(2): 504-512 (in Chinese with English abstract). Su Wei, Peng Yanju, Zheng Yuejun, Huang ZhongxianSu Wei, Peng Yanju, Zheng Yuejun, Huang Zhongxian. Crust and upper mantle shear velocity structure beneath the Tibetan Plateau and adjacent areas[J]. Acta Geoscientia Sinica, 2002, 23(3): 193-200 (in Chinese with English abstract). Sun Kezhong, Teng JiwenSun Kezhong, Teng Jiwen. The velocity distribution in the crust and upper mantle beneath the Xizang (Tibetan) Plateau from long period surface waves[J]. Acta Geophysica Sinica, 1985, 28(S1): 43-53 (in Chinese with English abstract). Sun LiSun Li. The improvement of the earthquake locating method of National Seismological Network[J]. Recent Developments in World Seismology, 2016(11): 12-16 (in Chinese with English abstract). Sun Ruomei, Zhang Zhongjie, Teng Jiwen, Zhu Jieshou. Discussion on the lateral extrusion of Tibetan mass-evidence of surface wave tomography. In: Zhang Zhongjie, Gao Rui, Lü Qingtian, Liu Zhenkuan (Editors). Study on the Chinese Continental Deep Structure and Geodynamics[M]. Beijing: Science Press, 2004: 872-880 (in Chinese with English abstract). Teng Jiwen, Yin Zhouxun, Liu Hongbing, Zhang Zhongjie, Hu Jiafu, Sun Kezhong, Wei JihunTeng Jiwen, Yin Zhouxun, Liu Hongbing, Zhang Zhongjie, Hu Jiafu, Sun Kezhong, Wei Jihun. The 2D and 3D lithosphere structure and continental dynamics of Qinghai-Xizang Plateau[J]. Acta Geophysica Sinica, 1994, 37(S1): 117-130 (in Chinese with English abstract). Teng Jiwen, Zhang Bingming, Hu Jiafu, Wen Yibo. The research of deep medium and structural environment for earthquake "Pregnancy". In: Chen Yuntai (Editor). Advances in Seismology in China[M]. Beijing: Seismological Press, 1997: 258-265 (in Chinese). Teng JiwenTeng Jiwen. Dynamic process of substance and energy exchanges in depths of the earth and formation of mineral resources[J]. Geotectonica et Metallogenia, 2003, 27(1): 3-21 (in Chinese with English abstract). Teng Jiwen, Ma Xueying, Zhang Xuemei, Liu Youshan, Pi JiaolongTeng Jiwen, Ma Xueying, Zhang Xuemei, Liu Youshan, Pi Jiaolong. Deep processes and dynamic responses of the generation and occurrence of the 2015 Nepal MS8.1 earthquake[J]. Chinese Journal of Geophysics, 2017, 60(1): 123-141 (in Chinese with English abstract). Wang Jicheng, Chen Zhubin, Jiang Haiyu, Lü HaoWang Jicheng, Chen Zhubin, Jiang Haiyu, Lü Hao. Based on successive very fast simulated annealing andthe grid subdivision of micro seismic localization algorithm research[J]. Science Technology and Engineering, 2017, 17(33): 248-252 (in Chinese with English abstract). Wang Shanshan, Li Canping, Li Qingren, Guan Yejun, Nie XunbiWang Shanshan, Li Canping, Li Qingren, Guan Yejun, Nie Xunbi. The seismic inversion using the fast simulated annealing algorithm[J]. Acta Geophysica Sinica, 1995, 38(S1): 123-134 (in Chinese with English abstract). Wang Yang, Liu Hong, Zhang Heng, Wang Zhiyang, Tang XiangdeWang Yang, Liu Hong, Zhang Heng, Wang Zhiyang, Tang Xiangde. A global optimized implicit staggered-grid finite-difference scheme for elastic wave modeling[J]. Chinese Journal of Geophysics, 2015, 58(7): 2508-2524 (in Chinese with English abstract). Wei Chao, Li Xiaofan, Zhang MeigenWei Chao, Li Xiaofan, Zhang Meigen. Quantum annealing optimization and geophysical inverse method[J]. Progress in Geophysics, 2007, 22(3): 785-789 (in Chinese with English abstract). Wittlinger G., Masson F., Poupinet G., Tapponnier P., Mei Jiang, Herquel G., Guilbert J., Achauer U., Xue Guanqi, Shi Danian, Lithoscope Kunlun TeamWittlinger G., Masson F., Poupinet G., Tapponnier P., Mei Jiang, Herquel G., Guilbert J., Achauer U., Xue Guanqi, Shi Danian, Lithoscope Kunlun Team. Seismic tomography of northern Tibet and Kunlun: Evidence for crustal blocks and mantle velocity contrasts[J]. Earth and Planetary Science Letters, 1996, 139(1/2): 263-279. Wu Jianping, Ming Yuehong, Ye Tailan, Zeng RongshengWu Jianping, Ming Yuehong, Ye Tailan, Zeng Rongsheng. Upper mantle velocity structure in the Qinghai-Xizang Plateau from inversion of the body waveforms[J]. Acta Geophysica Sinica, 1998, 41(S1): 15-25 (in Chinese with English abstract). Xu Zhiqin, Wang Qin, Zeng Lingsen, Liang Fenghua, Li Huaqi, Qi Xuexiang, Cai Zhihui, Li Zonghai, Cao HuiXu Zhiqin, Wang Qin, Zeng Lingsen, Liang Fenghua, Li Huaqi, Qi Xuexiang, Cai Zhihui, Li Zonghai, Cao Hui. Three-dimensional extrusion model of the Great Himalaya slice[J]. Geology in China, 2013, 40(3): 671-680 (in Chinese with English abstract). Yang WencaiYang Wencai. Nonlinear chaotic inversionof seismic traces: (I) theory and numerical experiments[J]. Acta Geophysica Sinica, 1993, 36(2): 222-232 (in Chinese with English abstract). Yang WencaiYang Wencai. Non-linear geophysical inversion methods: review and perspective[J]. Progress in Geophysics, 2002, 17(2): 255-261 (in Chinese with English abstract). Yang Xiaoping, Wu Guo, Chen Lichun, Li ChuanyouYang Xiaoping, Wu Guo, Chen Lichun, Li Chuanyou. The seismogenic structure of the April 25, 2015 MW7.8 Nepal earthquake in the southern margin of Qinghai-Tibetan Plateau[J]. Chinese Journal of Geophysics, 2016, 59(7): 2528-2538 (in Chinese with English abstract). Yanovskaya T.B., Ditmar P.G.Yanovskaya T.B., Ditmar P.G. Smoothness criteria in surface wave tomography[J]. Geophysical Journal International, 1990, 102(1): 63-72. DOI:10.1111/j.1365-246X.1990.tb00530.x Yao YaoYao Yao. Improvement on nonlinear geophysical inversion simulated annealing[J]. Acta Geophysica Sinica, 1995, 38(5): 643-650 (in Chinese with English abstract). Yao Zhenxing, Zhang LinbinYao Zhenxing, Zhang Linbin. Hybrid optimization method for acoustic impedance inversion[J]. Progress in Geophysics, 1999, 14(2): 1-6 (in Chinese with English abstract). Ye Xiuwei, Huang Yuanmin, Liu JipingYe Xiuwei, Huang Yuanmin, Liu Jiping. 3D P-wave velocity structure and active tectonics in the Xinfengjiang Area of Guangdong[J]. Earthquake Research in China, 2016, 32(3): 465-476 (in Chinese with English abstract). Zeng Rongsheng, Ding Zhifeng, Wu Qingju, Wu JianpingZeng Rongsheng, Ding Zhifeng, Wu Qingju, Wu Jianping. Seismological evidences for the multiple incomplete crustal subductions in Himalaya and southern Tibet[J]. Chinese Journal of Geophysics, 2000, 43(6): 780-797 (in Chinese with English abstract). Zhang Fengxue, Wu Qingju, Li Yonghua, Zhang Ruiqing, Sun Lian, Pan Jiatie, Ding ZhifengZhang Fengxue, Wu Qingju, Li Yonghua, Zhang Ruiqing, Sun Lian, Pan Jiatie, Ding Zhifeng. Seismic tomography ofeastern Tibet: Implications for the Tibetan Plateau growth[J]. Tectonics, 2018, 37(9): 2833-2847. DOI:10.1029/2018TC004977 Zhang Hongmei, Li Xiaolin, Liu HongZhang Hongmei, Li Xiaolin, Liu Hong. The improved method of simplex-simulated annealing residual statics[J]. Progress in Geophysics, 2004, 19(2): 341-347 (in Chinese with English abstract). Zhang Linbin, Yao Zhenxing, Ji Chen, Zhang ZhongjieZhang Linbin, Yao Zhenxing, Ji Chen, Zhang Zhongjie. Fast stimulated annealing algorithm and its application[J]. Oil Geophysical Prospecting, 1997, 32(5): 654-660 (in Chinese with English abstract). Zhang Xuemei, Teng Jiwen, Sun Ruomei, Romanelli F., Zhang Zhongjie, Panza G. F.Zhang Xuemei, Teng Jiwen, Sun Ruomei, Romanelli F., Zhang Zhongjie, Panza G. F. Structural model of the lithosphere-asthenosphere system beneath the Qinghai-Tibetan Plateau and its adjacent areas[J]. Tectonophysics, 2014, 634: 208-226. DOI:10.1016/j.tecto.2014.08.017 Zhao Junmeng, Du PinrenZhao Junmeng, Du Pinren. On the initial collision between the Indian and Eurasian continents[J]. Seismology and Geology, 2016, 38(3): 783-796 (in Chinese with English abstract). Zhao Wenjin, Nelson K.D., Che J., Quo J., Lu D., Wu C., Liu X.Zhao Wenjin, Nelson K.D., Che J., Quo J., Lu D., Wu C., Liu X. Deep seismic reflection evidence for continental underthrusting beneath southern Tibet[J]. Nature, 1993, 366(6455): 557-559. DOI:10.1038/366557a0 Zhao WenjinZhao Wenjin. Geological background of the Nepal's MS8.1 earthquake and its trend in the future[J]. Chinese Science Bulletin, 2015, 60(21): 1953-1957 (in Chinese). DOI:10.1360/N972015-00506 Zheng Chen, Ding Zhifeng, Song XiaodongZheng Chen, Ding Zhifeng, Song Xiaodong. Joint inversion of surface wave dispersion and receiver functions for crustal and uppermost mantle structure in Southeast Tibetan Plateau[J]. Chinese Journal of Geophysics, 2016, 59(9): 3223-3236 (in Chinese with English abstract). Zhong Yuyun, Zhang Zhenfeng, Kan BaoxiangZhong Yuyun, Zhang Zhenfeng, Kan Baoxiang. Simultaneous inversion of earthquake relocation and velocity structure in the Shanxi-reservoir, Wenzhou[J]. Earthquake Research in China, 2011, 25(4): 468-476. Zhu Jieshou. Calculation Method in Seismology[M]. Beijing: Seismological Press, 1988 (in Chinese). Zhu Jieshou, Cao Jiamin, Cai Xuelin, Yan ZhongqiongZhu Jieshou, Cao Jiamin, Cai Xuelin, Yan Zhongqiong. The structure of lithosphere in Eurasia and west Pacific[J]. Advance in Earth Sciences, 2004, 19(3): 387-92 (in Chinese with English abstract).