To accurately and efficiently model the complex spatiotemporal localization of wave propagation in heterogeneous media, we present a wavelet-optimized gridding scheme for the staggered-grid finite-difference (FD) numerical method. An average-interpolating (AI) wavelet transform with its wavelet coefficients obtained from linear combinations of samples, rather than from integrals, enables the gridding scheme constructed in the physical space to determine where to refine/coarsen the grid adapted to the complex spatial localization of problems and the relatively irregular sampling of seismic signals. The AI wavelet is a perfect transform with its wavelet coefficients decaying rapidly with increasing resolution. The wavelet scale parameter is related directly to the level of meshes. We can reconstruct the solution with good accuracy using a finite number of terms of wavelet expansion under some threshold sets for refinement/coarsening in a low- or high- degree of regularity. The adaptive meshes, in response to the heterogeneities of media and the resolution scales of wave fields, lead not only to a great saving of computation time and memory but also to an enabling method for solving problems with high accuracy in some local complex zones. Numerical examples demonstrate that the proposed wavelet-optimized adaptive FD scheme is efficient and accurate for wave propagation simulation.