A method is proposed for simultaneously estimating the trend and random component of soil properties at arbitrary locations using Gaussian process regression with the superposition of multiple Gaussian random fields. The proposed method is applied to the estimation of the one-dimensional spatial distributions of the trend component of three synthetic datasets. A comparison of three covariance functions, namely Gaussian, Markovian, and binary noise, indicates that Gaussian covariance is most suitable for trend estimation. The scale of fluctuation and the standard deviation of the random component of the examples are estimated using the maximum likelihood estimation method. The proposed method is also applied to the estimation of the three-dimensional spatial distribution of the trend and random components based on measured cone penetration test data. It is shown that the trend and random components at arbitrary locations can be estimated. The Whittle-Matérn covariance function is found to be more suitable than the Markovian covariance function for the estimation of the random component of the cone penetration test data based on the Akaike information criterion and the Bayesian information criterion.