Consider the inverse random source scattering problem for the two-dimensional time-harmonic elastic wave equation with an inhomogeneous, anisotropic mass density. The source is modeled as a microlocally isotropic generalized Gaussian random function whose covariance operator is a classical pseudo-differential operator. The goal is to recover the principle symbol of the covariance operator from the displacement measured in a domain away from the source. For such a distributional source, we show that the direct problem has a unique solution by introducing an equivalent Lippmann-Schwinger integral equation. For the inverse problem, we demonstrate that, with probability one, the principle symbol of the covariance operator can be uniquely determined by the amplitude of the displacement averaged over the frequency band, generated by a single realization of the random source. The analysis employs the Born approximation, asymptotic expansions of the Green tensor, and microlocal analysis of the Fourier integral operators.