In this work the Lippmann-Schwinger equation is used to model seismic waves in strongly scattering acoustic media. We consider the Helmholtz equation, which is the scalar wave equation in the frequency domain with constant density and variable velocity, and transform it to an integral equation of the Lippmann-Schwinger type. To directly solve the discretized problem with matrix inversion is time-consuming, therefore we use iterative methods. The Born series is a well-known scattering series which gives the solution with relatively small cost, but it has limited use as it only converges for small scattering potentials. There exist other scattering series with preconditioners that have been shown to converge for any contrast, but the methods might require many iterations for models with high contrast. Here we develop new preconditioners based on randomized matrix approximations and hierarchical matrices which can make the scattering series converge for any contrast with a low number of iterations. We describe two different preconditioners; one is best for lower frequencies and the other for higher frequencies. We use the fast Fourier transform both in the construction of the preconditioners and in the iterative solution, and this makes the methods efficient. The performance of the methods are illustrated by numerical experiments on two 2D models.