The representation of spatial discrete data is crucial 1 for the data analysis of meteorology, agriculture, geological ex-2 ploration and other fields. Among various methods for scattered 3 data interpolation, the multiquadric (MQ) method is the most 4 favorable in the field of surface construction. However, the 5 classical MQ method has accuracy concern on the boundary 6 and is time consuming for large data set. This study proposes 7 an algorithm by integrating the MQ method with trend surface 8 analysis. In which, the low-order polynomial trend surface 9 equation is firstly used to model the overall trend. Then, the 10 MQ equation is applied to fit the residual surface after removal 11 the trend from the data. Our implementation can eliminate the 12 distortion in data missing areas by the classical MQ method, 13 and the modeling efficiency can be improved significantly since 14 the local MQ method divide the residual surface into a group 15 of subsurfaces. The accuracy and efficiency of the proposed 16 algorithm are validated on a synthetic model. The performance 17 of the developed algorithm is further examined on the elevation 18 data collected in Tibet and the seabed of a strait in Norway. The 19 results show that with an equivalent resolution, the developed 20 algorithm can be much more efficient than the classical MQ 21 method and well-developed Kriging method. 22 Index Terms-overall trend, residual surface, local multi-23 quadric method 24 I. INTRODUCTION 25 G EO-data collection is limited by objective factors such 26 as terrain and weather, and the collected data points 27 are always irregularly distributed on the plane [1, 2]. The 28 interpolation techniques have been developed to improve the 29 understanding of the structure and geometry of geological 30 terrain [3, 4]. These interpolation algorithms include the trend 31 surface analysis, multiquadric (MQ) method, Kriging method, 32 thin-plate spline, minimum curvature method, inverse distance 33 weighted method, and so on [5, 6, 7]. Typically, the data value 34 is approximated by the polynomial functions of the geographic 35 coordinates of the data points. Once the coefficient matrix 36