To overcome the difficulty in forward and inversion problem in complex geological model, including undulated topography and irregular subsurface interface, we in this paper realize a multistage triangular shortest-path method for multi-phase seismic ray tracing and combine a nonlinear inversion solver (damped minimum norm and constrained least squares problem solved by a conjugated gradient method) to simultaneously invert four elastic parameters (or four Thomsen parameters) in anisotropic media using combined direct and reflected arrival times. For trading-off between different parameter updating, we normalize the sensitivity functions (traveltimes derivatives with respect to different elastic (or Thomsen) parameter changes) to balance four updated elastic (or Thomsen) parameters in simultaneous inversion process, due to the magnitude and pattern of partial derivatives more sensitive to different phase slowness angles. In addition, in order to solve the multi-solution problems of qSV-wave inversion, a strategy of multistage inversion is proposed. For testing the efficiency and correctness of proposed inversion method, we use the picked traveltimes in synthetic seismograms as the observed traveltimes in the inversion. The results show that both inverted images by taking the picked traveltimes in synthetic seismograms as the observed data and by taking predicted traveltimes using ray tracing method as observed data have a similar imaging pictures, which verify the efficiency and correctness of the proposed inversion method. Furthermore, the inverted results by normalized kernel functions are better than that of inversion without applying the normalized kernel functions for both qP and qSV data, and the results of the multistage inversion are better than that of the simultaneous inversion without multistage inversion applied in qSV-wave inversion. In addition, both elastic and Thomsen parameter inversions for the same anisotropic model indicate that there are no big differences between the two inverted results even though the sensitivity kernels related to the Thomsen parameters have relatively large values.