We have developed a method for the Bayesian location of seismic sequences in 3D heterogeneous velocity models. The method is based on a Bayesian algorithm for single earthquake location. Joint location of seismic sequences is performed by summing the probability density functions for individual earthquake locations. One of the main features of the method is the ability to accumulate travel times in the heterogeneous model, as a function of both the seismic stations and the grid points, simulating the seismogenic volume. The ability to compute travel times just once, separately from the location process, allows for fast 3D locations. This is a very attractive feature for real-time monitoring. Probabilistic location of seismic sequences, visualized in terms of contours of earthquake density, moment, and energy release, etc., can be obtained to give a much more seismotectonic insight than classical location algorithms. The continuous character of the output quantities is particularly indicated for seismic-network-testing purposes. It is in fact possible to compute the response of the location procedure, given the seismic network geometry, to various input earthquake distributions. Some applications of the method are also shown in this work, both to the location of synthetic earthquakes and of real sequences. Seismic sequences at Vesuvius and Campi Flegrei volcanic areas (southern Italy) have been located by the probability algorithm, allowing a better picture of the primary seismogenic structures. We also apply the method to assess the resolution of the operating monitoring networks at the two areas.