The authors present the new two-grid preconditioner for the solution of realistic three-dimensional heterogeneous Helmholtz problems discretized with second-order finite difference methods with application to acoustic waveform inversion in geophysics as e.g. monitoring of pollution in groundwater, earthquake modelling or location of hydrocarbon in fractured rocks. The homogeneous Helmholtz problems can be solved by efficient multilevel solvers, e.g. by the wave-ray multigrid method or by the FETI-H non-overlapping domain decomposition method whose rate of convergence is found to be independent of the fine grid step size, the number of subdomains, and the wavenumber. This paper focuses to the approximate solution of large indefinite linear systems of equations in the frequency domain being a topic of current research. The new proposed approach is based on an iterative two-grid method acting on the Helmholtz operator where the coarse grid problem is solved inaccurately with preconditioner consisting in a multigrid method cycle applied to a complex shifted Laplacian operator. These efficient algebraic one-level preconditioners are missing in other methods. A single cycle of the new method is then used as a variable preconditioner of a flexible Krylov subspace method whose properties are analysed by means of Fourier analysis giving us also appropriate relaxation parameters in the Jacobi method that lead to acceptable smoothing factors on all the grids of a complex shifted multigrid method in three dimensions. Numerical results demonstrate the effectiveness and robustness of the algorithm on three-dimensional heterogeneous media wave propagation applications even with high frequencies (or equivalently with large wavenumbers causing steep increase in the number of iterations for other methods) on a reasonable number of cores of a distributed memory computer. The paper is written in well understandable English.