We have developed a forward modeling technique to retrieve rupture characteristics of small earthquakes (3<M<5), including rupture propagation direction, fault dimension, and rupture speed. The technique is based on an empirical Green’s function (EGF) approach, where we use data from collocated smaller events as Green’s functions to study the bigger events. We tend to choose smaller events with similar focal mechanisms for EGFs; however, we show that the events with different focal mechanisms can work equally well when corrected for radiation pattern effect. Compared to deconvolution, this forward modeling approach allows full use of both the shape and amplitude information produced by rupture propagation. Assuming a simple 1D source model, we parameterize the source time function of a target event as the convolution of two boxcars, featuring the rise time τr and the rupture time τc; we solve for τr and τc in a grid search manner by minimizing the waveform misfit between the three-component data and the synthetics constructed from the EGFs. The rupture propagation direction, fault length, and rupture speed can then be estimated by fitting the observed azimuthal pattern of τc from P and S waves. We apply the approach to the 12 largest events (Mw≥3.3) of the 2003 Big Bear sequence (excluding the mainshock) in southern California. Among them, seven events are found to exhibit robust rupture directivity. The fact that the ruptures of these events propagate in all directions reveals complicated fault geometry at depth. We compute the stress drop for each event using the resolved fault length. The results show large variations ranging from ∼1 to 90 Mpa, with no dependence on moment. However, Δσ appears inversely correlated with rupture speed Vr; in particular, events with larger Δσ tend to propagate at smaller Vr, whereas events with smaller Δσ propagate faster.