The multispectral and hyperspectral image fusion (MHF) technique is designed to address the challenge of integrating the spatiotemporal characteristics of multispectral images with those of hyperspectral images. In the initial stages, the image to be fused is primarily decomposed into endmembers and abundances using a matrix decomposition method. This approach, however, may disrupt the correlation of spectral data. Subsequently, tensor decomposition-based methods emerged, with the most representative being canonical polyadic decomposition and Tucker decomposition. These methods are widely applied to the MHF problem due to their good recoverability. However, they lack the ability to introduce the physical interpretation of potential factors into the framework, and it is difficult to improve the quality of the fused image by utilizing the physical properties of the endmembers. Consequently, we employ a block term tensor decomposition algorithm based on sparse regularization to estimate the optimal high spatial resolution hyperspectral image. First, the abundance information is reconstructed into a chunk matrix, and its sparsity is characterized by introducing the l2,1 norm to eliminate the scaling effect present in the model. Second, the endmembers’ matrix sparsity is facilitated by the introduction of the row sparsity of the l2,1 norm regularization, which eliminates the inverse scaling effect present in the model. Finally, the model is solved using the conjugate alternating iteration algorithm. Experiments on three standard datasets and two local datasets demonstrate that this method outperforms state-of-the-art methods.