<p>Due to the influence of geological conditions, acquisition environments, and economic factors, collected seismic data are often incomplete, which severely affects subsequent data processing, such as multiple wave suppression and migration imaging. Therefore, achieving high-precision reconstruction of missing data is a critical issue in seismic data processing. For this reason, this paper proposes a three-dimensional randomly missing seismic data reconstruction method based on sparse learning via maximum-likelihood principle (SLML). This method fully utilizes the harmonic structure and statistical characteristics of frequency slices of three-dimensional seismic data, thereby significantly elevating the precision of data reconstruction. First, a Fourier transform is performed on the three-dimensional seismic data along the time axis. Based on the maximum likelihood principle, the problem of reconstructing missing seismic data are transformed into a two-dimensional harmonic spectrum estimation problem for the frequency slices. The majorization-minimization (MM) algorithm is then used for iterative harmonic spectrum estimation. Finally, the inverse Fourier transform is applied to the spectrum estimation to reconstruct the missing data. In addition, in order to improve reconstruction efficiency, the Toeplitz-block-Toeplitz structure of the covariance matrix and the two-dimensional inverse fast Fourier transform (IFFT) are utilized to perform the inversion operations during the reconstruction process. Extensive experiments were conducted on both synthetic and field data, and the experimental results show that the proposed SLML method achieves higher reconstruction accuracy when compared to classical mainstream methods.</p>