Propagation of acoustic waves in strongly heterogeneous elastic media is studied using renormalization group analysis and extensive numerical simulations. The heterogeneities are characterized by a broad distribution of the local elastic constants. We consider both Gaussian-white distributed elastic constants, as well as those with long-range correlations with a nondecaying power-law correlation function. The study is motivated in part by recent analysis of experimental data for the spatial distribution of the elastic moduli of rock at large length scales, which indicated that the distribution contains the same type of long-range correlations as what we consider in the present paper. The problem that we formulate and the results are, however, applicable to acoustic wave propagation in any disordered elastic material that contains the types of heterogeneities that we consider in the present paper. Using the Martin-Siggia-Rose method, we analyze the problem analytically and find that, depending on the type of disorder, the renormalization group (RG) flows exhibit a transition to a localized or extended regime in any dimension. We also carry out extensive numerical simulations of acoustic wave propagation in one-, two-, and three-dimensional systems. Both isotropic and anisotropic media (with anisotropy being due to stratified) are considered. The results for the isotropic media are consistent with the RG predictions. While the RG analysis, in its present form, does not make any prediction for the anisotropic media, the results of our numerical simulations indicate the possibility of the existence of a regime of superlocalization in which the waves' amplitudes decay as $\mathrm{exp}[\ensuremath{-}{(\ensuremath{\mid}\mathbf{x}\ensuremath{\mid}∕\ensuremath{\xi})}^{\ensuremath{\gamma}}]$, with $\ensuremath{\gamma}>1$, where $\ensuremath{\xi}$ is the localization length. However, further investigations may be necessary in order to establish the possible existence of such a localization regime.