In seismology, blind deconvolution aims to recover the source wavelet and the Green’s function, or parts of it (e.g., reflectivity series), from a recorded seismic trace. A multitude of algorithms exist that tackle this ill-posed problem by different approaches. Making assumptions on the phase spectra of the source wavelet and/or the statistical distribution of the reflectivity series is useful for single trace. The nature of closely spaced multichannel recordings enables a better estimation of a common source wavelet and thus increases the confidence of the results. This approach has been exploited in the past, although different types of assumptions are used for a variety of algorithms. We introduced a new method for simultaneous reconstruction of arbitrary source wavelets and local vertical reflectivity series from teleseismic earthquakes. Closely spaced receivers record vertically incident earthquake body waves and their surface-related multiples, which comprise the unknown reflectivity series. By assuming a common source wavelet for all receivers, the observation of several events resulted in a set of convolution equations relating the unknown source wavelets and unknown reflectivity series to the observed seismic trace. The overdetermined system of equations was linearized and solved by conventional inversion algorithms in the spectral domain. Synthetic tests indicated a better performance of the introduced method than conventional deconvolution in the presence of white noise, which is attributed to the constraint of a common model for all observations. Application to field data from a local deployment allowed imaging a basement reflector from teleseismic body waves, although the data were contaminated with strong coherent noise. From a practical point of view, the presented method is potentially well suited for local and regional large-scale imaging from multichannel passive seismic data.