Asymptotic ray methods are commonly used for forward modelling in seismology and seismic exploration, but rays propagating through heterogeneous media exhibit chaotic behaviour. This means that we can observe a sensitive dependence on initial ray conditions and thus, for long times, the exponential amplification of very small perturbations. The degree of chaos in the model, which may be quantified in terms of the average Lyapunov exponent for all the rays in the model, can be approximated by the square root of the Sobolev norm of order 2 of the corresponding model parameters. In this paper, we present a new inversion technique which considers the Sobolev norm as the regularization term of the ill-conditioned seismic inverse problem. We show that in this way we can control the divergence of rays in our model during the iterations and obtain stable solutions with a negligible dependence on the chosen discretization. The robustness and limitations of the algorithm are illustrated by numerical examples produced with a real heterogeneous 3D structure. The sensitivity of the algorithm to non-Gaussian distribution of noise is also tested.