Thin plate spline (TPS) has been widely accepted as a method for smooth fitting of noisy data. However, the classical TPS always has an ill-conditioning problem when two sample points are very close. Although the modified orthogonal least squares-based TPS (TPS-M) avoids this ill-conditioning problem, it is not completely immune to over-fitting when sample points are noisy. In this paper, a regularized least squares method of thin plate spline (TPS-RLS) was developed, which adds a weight penalty term to the error criterion of orthogonal least squares (OLS). TPS-RLS combines the advantages of both regularization and OLS, which avoid the over-fitting and the ill-conditioning problems simultaneously. Numerical tests indicate that irrespective of the standard deviation of sampling errors and the number of knots, TPS-RLS is always more accurate than TPS-M for smooth fitting of noisy data, whereas TPS-M would have a serious over-fitting problem if the optimal number of knots were not determined in advance. The real-world example of fitting total station instrument data shows that among the classical interpolation methods including IDW, natural neighbor and ordinary kriging, TPS-RLS has the highest accuracy for a series of DEMs with different resolutions, especially for the coarse one. Surface modeling of DEMs with contour lines demonstrate that TPS-RLS has a better performance than the classical methods in terms of both root mean squared error and relief shaded map appearance.