Seismic inversion allows the prediction of subsurface properties from seismic reflection data and is a key step in reservoir modeling and characterization. Traditional seismic inversion methods usually come with a high computational cost or suffer from issues concerning the non-linearity and the strong non-uniqueness of the seismic inversion model. With the generalization of machine learning in geophysics, deep learning methods have been proposed as efficient seismic inversion methods. However, most of them lack a probabilistic approach to deal with the uncertainties inherent in the seismic inversion problems and/or rely on complete and representative training data, which is often scarcely available. To overcome these limitations, we introduce a novel seismic inversion method that explores the ability of deep convolutional neural networks to extract meaningful and complex representations from spatially structured data, combined with geostatistical stochastic simulation to efficiently invert seismicn reflection data directly for high-resolution subsurface models. Our method incorporates physics constraints, sparse direct measurements, and leverages the use of imprecise but widely distributed indirect measurements as represented by the seismic data. The geostatistical realizations provide additional information with higher spatial resolution than the original seismic data. When used as input to our inversion system, they allow the generation of multiple possible outcomes for the uncertain model. Our approach is fully unsupervised, as it does not depend on ground truth input-output pairs. In summary, the proposed method is able to: (1) provide uncertainty assessment of the predictions, (2) model the complex non-linear relationship between observed data and model, (3) extend the seismic bandwidth at both low and high ends of the frequency parameters spectrum, and (4) lessen the need for large, annotated training data. The proposed methodology is first described in the acoustic domain to invert acoustic impedance models from full-stack seismic data. Next, it is generalized for the elastic domain to invert P-wave velocity, S-wave velocity and density models from pre-stack seismic data. Finally, we show that the proposed methodology can be further extended to perform petrophysical seismic inversion in a simultaneous workflow. The method was tested on a synthetic case and successfully applied to a real three-dimensional case from a Brazilian reservoir. The inverted models are compared to those obtained from a full iterative geostatistical seismic inversion. The proposed methodology allows retrieving similar models but has the advantage of generating alternative solutions in greater numbers, providing a larger exploration of the model parameter space in less time than the geostatistical seismic inversion.