This dissertation addresses forward and inverse modeling of Rayleigh waves for near surface investigation. Results were obtained by imaging abandoned mine openings using Rayleigh waves in the laterally inhomogeneous medium. The efficient staggered grid stencil method to solve elastic wave equations using 2-D finite difference technique is presented. This numerical scheme is used to conduct a series of parametric studies on the propagation of Rayleigh waves. The first parametric study was conducted on a flat layered model of increasing and decreasing velocity with depth. A Rayleigh waves dispersion curve is found to be sensitive on a layer’s depth up to half of the minimum wavelength of Rayleigh waves. The phase velocity in the dispersion curve of Rayleigh waves is inversely and directly proportional to the frequency, depending on velocity increase or decrease with depth. The parametric study was carried out by introducing dipping layers in the model with increasing dip. The front (near the shot point) and back (at the end of receiver line) shot records are different if the subsurface contains dip. Dispersion is observed in near offset for down dip and in the far offset for up dip, computed from front and back shots respectively. Finally, a parametric study looked at subsurface anomalies with different shapes and sizes as well as their material properties. A Rayleigh wave is sensitive to very high material contrast and very low material contrast of the anomaly from it surrounding medium. The presence of a low material contrast anomaly from the surrounding medium traps the energy which causes reverberation. A Rayleigh wave is sensitive to an anomaly which is placed within the depth between one-third to half of minimum wavelength of Rayleigh wave from the surface. In order to resolve lateral heterogeneity, a new method is developed in this research which allows localization of the multichannel record in different panels. The dispersion curve of Rayleigh waves is computed in each panel using the slant stack technique. On the basis of parametric studies, an innovative inversion algorithm has been developed to minimize the error norm; ”the sum of the squares of the difference of reference and model dispersion curves” in an iterative way using a Very Fast Simulated Re-annealing (VFSR) technique.