We develop a novel three-dimensional viscoelastic modeling method using an implicit finite difference scheme, while the properties of sparse and block matrices are utilized to make the solution of the method possible. The advantage of this method is that the introduction of implicit format makes it alleviate the numerical dispersion problem with the same accuracy and has better stability. In order to verify the validity and accuracy of the method, this paper elaborates its basic principles and analyzes the results through a series of forward simulations of simple models. In addition, we also gives an outlook on the direction of the subsequent refinement of the 3D viscoelastic implicit finite difference method as well as the potential application prospects.