We have developed a new sparse-spike deconvolution (SSD) method based on Toeplitz-sparse matrix factorization (TSMF), a bilinear decomposition of a matrix into the product of a Toeplitz matrix and a sparse matrix, to address the problems of lateral continuity, effects of noise, and wavelet estimation error in SSD. Assuming the convolution model, a constant source wavelet, and the sparse reflectivity, a seismic profile can be considered as a matrix that is the product of a Toeplitz wavelet matrix and a sparse reflectivity matrix. Thus, we have developed an algorithm of TSMF to simultaneously deconvolve the seismic matrix into a wavelet matrix and a reflectivity matrix by alternatively solving two inversion subproblems related to the Toeplitz wavelet matrix and sparse reflectivity matrix, respectively. Because the seismic wavelet is usually compact and smooth, the fused Lasso was used to constrain the elements in the Toeplitz wavelet matrix. Moreover, due to the limitations of computer memory, large seismic data sets were divided into blocks, and the average of the source wavelets deconvolved from these blocks via TSMF-based SSD was used as the final estimation of the source wavelet for all blocks to deconvolve the reflectivity; thus, the lateral continuity of the seismic data can be maintained. The advantages of the proposed deconvolution method include using multiple traces to reduce the effect of random noise, tolerance to errors in the initial wavelet estimation, and the ability to preserve the complex structure of the seismic data without using any lateral constraints. Our tests on the synthetic seismic data from the Marmousi2 model and a section of field seismic data demonstrate that the proposed method can effectively derive the wavelet and reflectivity simultaneously from band-limited data with appropriate lateral coherence, even when the seismic data are contaminated by noise and the initial wavelet estimation is inaccurate
