Abstract: Objectives: In the process of digital elevation model (DEM) modeling, the traditional interpolation methods are based on the assumption that the surface is smooth and continuous, only considering the spatial correlation between the sampling point and the interpolated point, while ignoring the influence of heterogeneous distribution of elevations in areas with breaklines. This causes the elevation around the breaklines to be smoothed, and distorts the constructed DEMs. Methods: Therefore, this study proposes a multivariate radial basis function interpolation methodthat takes into account spatial heterogeneity. This method couples three kinds of terrain information, including spatial distance, height difference, and normal vector, and fully considers the spatial correlation and heterogeneity between the sampling point and the interpolated point. Results: Taking 10 public datas provided by ISPRS and 1 airborne LiDAR point cloud data as examples, we compared the proposed method with the interpolation method that considers the structural tensor constraint, and three traditional interpolation methods including standard Radial Basis Function (RBF), Triangulated Irregular Network (TIN), Australian National University Digital Elevation Model(ANUDEM). Results show that the average total error of the method in this paper is the smallest, the interpolation performance is the best, and it can better maintain terrain features in the breakline area. Conclusions: In short, the quality of the treatment of breakline terrain largely affect the authenticity and accuracy of the DEM's terrain expression. The proposed method can effectively capture the spatial distribution characteristics of breakline terrain, which is conducive to achieving high-quality DEM modeling.