Sparse-spike deconvolution can be viewed as an inverse problem where the locations and amplitudes of a number of spikes (reflectivity) are estimated from noisy data (seismic traces). The main objective is to find the least number of spikes that, when convolved with the available band-limited seismic wavelet estimate, fit the data within a given tolerance error (misfit). The detection of the spikes’ time lags is a highly nonlinear optimization problem that can be solved using very fast simulated annealing (SA). Amplitudes are easily estimated using linear least squares at each SA iteration. At this stage, quadratic regularization is used to stabilize the solution, to reduce its nonuniqueness, and to provide meaningful reflectivity sequences, thus avoiding the need to constrain the spikes’ time lags and/or amplitudes to force valid solutions. Impedance constraints also can be included at this stage, providing the low frequencies required to recover the acoustic impedance. One advantage of the proposed method over other sparse-spike deconvolution techniques is that the uncertainty of the obtained solutions can be estimated stochastically. Further, errors in the phase of the wavelet estimate are tolerated, for an optimum constant-phase shift is obtained to calibrate the effective wavelet that is present in the data. Results using synthetic data (including simulated data for the Marmousi2 model) and field 3D data show that physically meaningful high-resolution sparse-spike sections can be derived from band-limited noisy data, even when the available wavelet estimate is inaccurate.