A key stage for Kriging interpolation is the estimating of variogram model, which characterizes the spatial behavior of the variables of interest. But most traditional kriging interpolation has finite types of empirical variogram model, and sometimes, the optimal type of variogram model can not be find, which result in decreasing interpolation accuracy. In this paper, we explore the use of Multi-Gene Genetic Programming (MGGP) to automatically find an empirical variogram model that fits on an experimental variogram. Empirical variogram estimation based on MGGP, in contrast with traditional method need not select type of basic variogram model and can directly get both the functional type as well as the coefficients of the optimal variogram. The results of case study show that the proposed method can avoid the subjectivity in choosing the type of variogram models and can adaptively fit variogram according to the real data structure, which improves the interpolation accuracy of kriging significantly.