2. Institute of Geophysics, China Earthquake Administration, Beijing 100081, China
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 nonconvergence 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 nonlinear. For example, there is a highly nonlinear relationship between seismic surface wave dispersion curve and velocity structure of the earth. Linear inversion algorithm is not suitable for solving nonlinear problems. In recent years, with the improvement of computational methods and the realization of supercomputing, many nonlinear 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 nonlinear 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 Swave velocity and layers of the crust and upper mantle.
1 IMPROVEMENT OF SIMULATED ANNEALING ALGORITHMSSimulated 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 E_{i} 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. E_{i} is used to indicate its energy. In the new physical state j with the random perturbation, its energy is E_{j}. If the energy E_{j} < energy E_{i}, then the state j is regarded as the important state of the system. If the energy E_{j} > energy E_{i}, 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 selfadaptive 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 temperaturedependent Cauchy distribution to generate a new perturbation model, which performs largescale 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 nonlinear 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 lithosphereasthenosphere system, we improved the simulated annealing nonlinear 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 variablescale cooling function to choose largerscale cooling mode at higher temperature, and then use smallerscale 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)=T_{0}α_{i}^{k} cooling function, where T(k) is the temperature of k times, and T_{0} is initial temperature value, and α_{i} is the coefficient of cooling, generally α_{i}∈[0.7, 1].
The calculation steps are as follows:
(1) Let [v_{o}, D_{o}] be the initial model; k=1; T_{o}=T_{max} (initial maximum temperature); where v_{o} is the velocity vector of the initial model and D_{o} 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 [v_{j}, D_{o}] from the model neighborhood N([v_{i}, D_{o}]), v_{j}=v_{i}+(γT(k)((1+T(k)^{1})^{(2random(0, 1))}1)) (v_{max}v_{min}) where
Calculate the objective function ΔY_{ij}=f ([v_{j}, D_{o}])f([v_{i}, D_{o}]), where f is the dispersion operator; if ΔY_{ij}≤0, set the new model [v_{i}, D_{o}]=[v_{j}, D_{o}], otherwise, if exp(ΔY_{ij}/T(k))>random(0, 1), set the new model [v_{i}, D_{o}]=[v_{j}, D_{o}]; repeat step (2);
(3) Decrease the temperature T_{k}=T_{0}α_{1}^{k}; 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 [v_{jj}, D_{jj}] from the model neighborhood N([v_{ii}, D_{ii}]),
v_{jj}=v_{ii} +(γT(k)((1+T(k)^{1})^{(2random(0, 1))}1))(v_{max}v_{min}),
D_{jj}=D_{ii}+(γT(k)((1+T(k)^{1})^{(2random(0, 1))}1)) (D_{max}D_{mins}),
where D_{max} is the maximum thickness and D_{min} is the minimum thickness;
Calculate the objective function ΔY_{ij}=f([v_{jj}, D_{jj}])f([v_{ii}, D_{ii}]); if ΔY_{ij} < =0, set the new model [v_{ii}, D_{ii}]=[v_{jj}, D_{jj}], if exp(ΔY_{ij}/T(k))>random(0, 1), set new model [v_{ii}, D_{ii}]=[v_{jj}, D_{jj}]; repeat step (4);
(5) Decrease the temperature T_{k}=T_{k1}α_{2}^{k}; k=k+1; if the loop condition is satisfied, output the inversion model [v_{ii}, D_{ii}]; otherwise, return to step (4).
2 NUMERICAL EXPERIMENTSA 7layer 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 Pwave velocity α are obtained from the thickness h and Swave 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 Swave velocity β, Pwave 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 Swave velocity β are taken as independent variables, while Pwave velocity α and other parameters are taken as independent variables, which are derived from Swave velocity β and the depth.
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 Swave 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 nonlinear 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 nonlinear inversion algorithm can fit the dispersion curve well and get the approximate theoretical velocity model.
3 APPLICATION: SURFACE WAVE DISPERSION INVERSION IN THE WESTERN QINGHAITIBETAN PLATEAUStudying the spatial tectonic framework of the lithosphereshoal 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 QinghaiTibetan Plateau, which is located in the place of the continentcontinent 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 QinghaiTibetan 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 QinghaiTibetan Plateau and its adjacent areas and obtain the Rayleigh wave group velocity of 8150s period by surface tomography. The cells along the 86°E cross the hinterland of the QinghaiTibetan 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.
Based on the long period surface wave data of the China Seismic Network and the Global Seismic Network, after preprocessing, we applied the YanovskayaDitmar tomography (Yanovskaya T. B. et al., 1990) to inverse the surface wave group velocity at the period of 8150s in the QinghaiTibetan Plateau and adjacent areas, and the dispersion data set of the 2°×2°size cell was discretized. The Swave velocity structure obtained by linear inversion is used as the initial model (Sun Ruomei et al., 2004). The improved simulated annealing nonlinear algorithm is used to simutanously inverse the velocity and thickness for each cell. The velocity distribution of Swave across the NorthSouth section is obtained (Fig. 3). Compared with the linear inversion results, the velocity distribution obtained by the twoparameter nonlinear inversion method can relatively clearly determine the velocity offsets of Moho and lithospheric bottom. On the vertical velocity distribution map, we take the Swave velocity increasing from ~3.7km/s to ~4.3km/s as the Moho discontinuity; the depth where Swave velocity begins to reduce is considered as the bottom of the lithosphere; and the low velocity zone below the lithosphere is the asthenosphere. From Swave 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 300400km depth is wider, and the bottom of asthenosphere can only be used as a reference.
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 3540km, 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 BangongNujiang suture belt (BNS) and the Yarlung Zangbo suture belt (YZS). The depth of Moho is 6070km and a thinner lithosphere is about 160km. The thickest crust is located in the Qiangtang block (QT) in the hinterland of the QinghaiTibetan 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 crustmantle 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 150160km. 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 QinghaiTibetan 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 QinghaiTibetan 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 deepsource earthquakes (depth greater than 70km) have occured in the crust and the upper mantle lid.
4 DISCUSSION AND CONCLUSIONSUsing the improved simulated annealing nonlinear algorithm to simultaneously invert the Swave velocity structure along the 86°E section, the following deductions are obtained:
(1) Applying the improved simulated annealing algorithm, the variablescale 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 QinghaiTibetan 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 QinghaiTibetan 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.
ACKNOWLEDGEMENTThe 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.
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): 14091412. 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): 161173 (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): 19851995 (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): 402408. 
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): 678697. DOI:10.1002/cjg2.20133 
Deng Wanming, Zhong DalaiDeng Wanming, Zhong Dalai. Crustmantle transition zone and its geological significance in the lithospheric tectonic evolution[J]. Chinese Science Bulletin, 1997, 42(23): 24742482 (in Chinese). DOI:10.1360/csb199742232474 
Ding Zhifeng, He Zhengqin, Wu Jianping, Sun WeiguoDing Zhifeng, He Zhengqin, Wu Jianping, Sun Weiguo. Research on the 3D seismic velocity structures in QinghaiXizang Plateau[J]. Earthquake Research in China, 2001, 17(2): 202209 (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 M_{S}6.6 Damxung earthquake, Tibet, China[J]. Chinese Journal of Geophysics, 2010, 53(9): 20832090 (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): 39463964. 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): 403420. 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): 2121521232. 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): 1734. 
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 southnorth traverse below the Tibetan Plateau: evidence from seismology[J]. Tectonophysics, 1997, 273(1/2): 1730. 
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: 100111. 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): 233242 (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): 237330. 
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): 355364. 
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 QinghaiXizang inferred from INDEPTHⅢ teleseismic waveforms: Geological and geophysical implications[J]. Chinese Journal of Geophysics, 2006, 49(4): 10371044 (in Chinese with English abstract). 
Liu Pengcheng, Ji Chen, Hartzell S. H.Liu Pengcheng, Ji Chen, Hartzell S. H. An improved simulated annealingdownhill simplex hybrid global inverse algorithm[J]. Acta Geophysica Sinica, 1995, 38(2): 199205. 
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 prestack inversion of reservoir parameters[J]. Progress in Geophysics, 2008, 23(1): 104109. 
Lü Zhiqiang, Lei JiansheLü Zhiqiang, Lei Jianshe. 3D Swave velocity structure around the 2015 M_{S}8.1 Nepal earthquake source areas and strong earthquake mechanism[J]. Chinese Journal of Geophysics, 2016, 59(12): 45294543 (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): 317. 
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): 1365513665. 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): 10871092. DOI:10.1063/1.1699114 
Ni J., Barazangi M.Ni J., Barazangi M. Highfrequency seismic wave propagation beneath the Indian Shield, Himalayan Arc, Tibetan Plateau and surrounding regions: high uppermost mantle velocities and efficient S_{n} propagation beneath Tibet[J]. Geophysics Journal of the Royal Astronomical Society, 1983, 72(3): 665689. DOI:10.1111/j.1365246X.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): 183193. DOI:10.1016/0012821X(69)900351 
Panza G. F., Simmons G.Panza G. F., Simmons G. Synthetic Seismograms: the Rayleigh waves modal summation[J]. Journal of Geophysics, 1985, 58: 125145. 
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): 11101116 (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): 426429. DOI:10.1063/1.4823002 
Rothman D.H.Rothman D.H. Automatic estimation of large residual statics corrections[J]. Geophysics, 1986, 51(2): 332346. 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 phasevelocity spectra[J]. Geophysics, 2006, 71(4): R49R58. DOI:10.1190/1.2204964 
Schwab F., Knopoff L.Schwab F., Knopoff L. Surfacewave dispersion computations[J]. Bulletin of the Seismological Society of America, 1970, 60(2): 321344. 
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): 694707 (in Chinese with English abstract). 
Su Jinrong, Guo ZhiSu Jinrong, Guo Zhi. A study of focal mechanisms of the M_{S}6.7 earthquake of 5 October 2008 in southwestern TianshanPamier[J]. Chinese Journal of Geophysics, 2013, 56(2): 504512 (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): 193200 (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): 4353 (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): 1216 (in Chinese with English abstract). 
Sun Ruomei, Zhang Zhongjie, Teng Jiwen, Zhu Jieshou. Discussion on the lateral extrusion of Tibetan massevidence 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: 872880 (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 QinghaiXizang Plateau[J]. Acta Geophysica Sinica, 1994, 37(S1): 117130 (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: 258265 (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): 321 (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 M_{S}8.1 earthquake[J]. Chinese Journal of Geophysics, 2017, 60(1): 123141 (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): 248252 (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): 123134 (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 staggeredgrid finitedifference scheme for elastic wave modeling[J]. Chinese Journal of Geophysics, 2015, 58(7): 25082524 (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): 785789 (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): 263279. 
Wu Jianping, Ming Yuehong, Ye Tailan, Zeng RongshengWu Jianping, Ming Yuehong, Ye Tailan, Zeng Rongsheng. Upper mantle velocity structure in the QinghaiXizang Plateau from inversion of the body waveforms[J]. Acta Geophysica Sinica, 1998, 41(S1): 1525 (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. Threedimensional extrusion model of the Great Himalaya slice[J]. Geology in China, 2013, 40(3): 671680 (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): 222232 (in Chinese with English abstract). 
Yang WencaiYang Wencai. Nonlinear geophysical inversion methods: review and perspective[J]. Progress in Geophysics, 2002, 17(2): 255261 (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 M_{W}7.8 Nepal earthquake in the southern margin of QinghaiTibetan Plateau[J]. Chinese Journal of Geophysics, 2016, 59(7): 25282538 (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): 6372. DOI:10.1111/j.1365246X.1990.tb00530.x 
Yao YaoYao Yao. Improvement on nonlinear geophysical inversion simulated annealing[J]. Acta Geophysica Sinica, 1995, 38(5): 643650 (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): 16 (in Chinese with English abstract). 
Ye Xiuwei, Huang Yuanmin, Liu JipingYe Xiuwei, Huang Yuanmin, Liu Jiping. 3D Pwave velocity structure and active tectonics in the Xinfengjiang Area of Guangdong[J]. Earthquake Research in China, 2016, 32(3): 465476 (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): 780797 (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): 28332847. DOI:10.1029/2018TC004977 
Zhang Hongmei, Li Xiaolin, Liu HongZhang Hongmei, Li Xiaolin, Liu Hong. The improved method of simplexsimulated annealing residual statics[J]. Progress in Geophysics, 2004, 19(2): 341347 (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): 654660 (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 lithosphereasthenosphere system beneath the QinghaiTibetan Plateau and its adjacent areas[J]. Tectonophysics, 2014, 634: 208226. 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): 783796 (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): 557559. DOI:10.1038/366557a0 
Zhao WenjinZhao Wenjin. Geological background of the Nepal's M_{S}8.1 earthquake and its trend in the future[J]. Chinese Science Bulletin, 2015, 60(21): 19531957 (in Chinese). DOI:10.1360/N97201500506 
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): 32233236 (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 Shanxireservoir, Wenzhou[J]. Earthquake Research in China, 2011, 25(4): 468476. 
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): 38792 (in Chinese with English abstract). 