This study presents stochastic methods to simulate wave envelopes in layered random media. High-frequency Seismograms of small earthquakes are so complex due to lithospheric inhomogeneity that seismologists often analyze wave envelopes rather than wave traces to quantify the subsurface inhomogeneity. Since the statistical properties of the inhomogeneity vary regionally, it is important to develop and examine direct envelope simulation methods for non-uniform random media. As a simple example,  this study supposes plane wave propagation through two-layer random media in 2-D composed of weak and strong inhomogeneity zones. The characteristic spatial-scale of the inhomogeneity is supposed to be larger than the wavelength, where small-angle scattering around the forward direction dominates large-angle scattering. Two envelope simulation methods based on the small-angle scattering approximation are examined. One method is to solve a differential equation for the two-frequency mutual coherence function with the Markov approximation. The other is to solve the stochastic ray bending process by using the Monte Carlo method based on the Markov approximation for the mutual coherence function. The resultant wave envelopes of the two methods showed excellent coincidence both for uniform and for two-layer random media. Furthermore, we confirmed the validity of the two methods comparing with the envelopes made from the finite difference simulations of waves. The two direct envelope simulation methods presented in this study can be a mathematical base for the study of high-frequency wave propagation through randomly inhomogeneous lithosphere in seismology.