As high-resolution global coverage cannot easily be achieved by direct bathymetry, the use of gravity data is an alternative method to predict seafloor topography. Currently, the commonly used algorithms for predicting seafloor topography are mainly based on the approximate linear relationship between topography and gravity anomaly. In actual application, it is also necessary to process the corresponding data according to some empirical methods, which can cause uncertainty in predicting topography. In this paper, we established analytical observation equations between the gravity anomaly and topography, and obtained the corresponding iterative solving method based on the least square method after linearizing the equations. Furthermore, the regularization method and piecewise bilinear interpolation function are introduced into the observation equations to effectively suppress the high-frequency effect of the boundary sea region and the low-frequency effect of the far sea region. Finally, the seafloor topography beneath a sea region (117.25°-118.25° E, 13.85°-14.85° N) in the South China Sea is predicted as an actual application, where gravity anomaly data of the study area with a resolution of 1′×1′ is from the DTU17 model. Comparing the prediction results with the data of ship soundings from the National Geophysical Data Center (NGDC), the root-mean-square (RMS) error and relative error can be up to 127.4 m and approximately 3.4%, respectively.