## 1 Introduction  
In soft rock strata, the large squeezing deformation of surrounding rock often occurs during tunneling under high stress [1; 2; 3; 4; 5; 6]. Many interesting cases of tunnels, for example, the section of the F7 fault fracture zone of the Wushaoling Railway Tunnel [7], and the Saint Martin La Porte access adit of the Lyon-Turin Base Tunnel [2; 3], were recorded as representative of large squeezing deformation disasters. These disasters are mostly caused by soft rock rheology [1; 2; 3; 4]. Therefore, it is of theoretical and practical interest to investigate the large deformation mechanism of deep soft rock tunnels and to explore an efficient numerical modeling method. The accelerated creep of soft rock plays an important role in tunnel large deformation disasters under high stress. Many previous works were devoted to investigating the time-dependent deformation of surrounding rock around tunnels using viscoelastic models and viscoplastic models [5]. However, these studies only involved the attenuation creep and uniform creep of soft rocks [1; 2; 3; 4; 5; 6], and seldom involved the accelerated creep. The identification of rheological curves in the high-stress soft rock test shows that the viscosity coefficient related to the non-stationary parameter creep is often a function of the stress level [8], and the previous viscoelastic-plastic models such as the Nishihara model, cannot simulate the accelerated rheological characteristics of soft rock. Therefore, only the nonlinear rheological model with the non-stationary parameter creep can completely present the accelerated creep characteristics of soft rock under high stress. In engineering practice, with the variation of the engineering environment and tunneling process, the mechanical parameters of surrounding rock usually have uncertainty. The deformation of surrounding rock and the supporting structure is the most direct, obvious and comprehensive macroscopic manifestation on the force change of the surrounding rock support system in the time-space domain of tunneling and support. Therefore, the mechanical parameters obtained by the displacement back analysis can reflect the influence of geo-stress, the construction method, and the tunnel structure parameters as a whole. At present, the back analysis has been widely used to obtain the mechanical parameters of surrounding rocks by simulating their complicated mechanical behaviors through simple models combined with measured data [9; 10; 11; 12; 13; 14; 15]. [PERSON] and [PERSON] [16] proposed the back analysis theory for the first time, and they inverted the elastic modulus of elastic solids based on the finite element method. Since then, some traditional optimization algorithms such as Levenber-Marquard, Gauss-Newton, Bayesian and other methods were applied to the field of back analysis. [PERSON] [17] used the mixed algorithm of the variable rotation method and the simplex method to invert the mechanics of elastic-plastic rock mass. Later, some scholars proposed the parameter inverted methods for five commonly used creep models (such as Maxwell, Kelvin, generalized Kelvin and Burgers models) [18; 19]. Although the above methods have made great progress in the field of back analysis of geotechnical parameters, with the progress of machine learning theory, some scholars have proposed numerous algorithms with higher computational accuracy and computational efficiency. Compared with traditional algorithms, these algorithms can better deal with the complex nonlinear relationship in the geotechnical back analysis, so as to achieve higher accuracy. The mechanical parameters of surrounding rock may be obtained by the displacement back analysis using neural networks and machine learning, such as the genetic algorithm and particle swarm optimization [20; 21; 22; 23; 24; 25; 26; 27]. In this study, we propose a nonlinear generalized Nishihara rheological model with the non-stationary parameter creep (NGNRM), and implement it in ABAQUS.