This paper presents a new method Discrete Smooth Interpolation method (DSI) and its fast realization in geophysical potential data gridding. First the topological relations between the unknown points and known points are built as a linear constrain, then the DSIequation group is formulated by combining Laplace equation with the linear constrain equation. The preconditioning conjugate direction descend method is used to seek the solution and speed up the solving. To improve the efficiency of the computation, the theory of the helical coordinate system is used to convert the two dimensional data to one dimensional data and deal with them in one dimensional space. In addition, the method of Wilson Burg spectral factorization is employed to solve the linear algebraic equation derived from DSI. Through the application of theoretical model and the real raw data, it is shown that the method in this paper is suitable for the features of geophysical potential field, the effects of interpolation is satisfactory and the convergence rate is very high.