Gaussian Processes are used in many applications to model spatial phenomena. Within this context, a key issue is to decide the set of locations where to take measurements so as to obtain a better approximation of the underlying function. Current state of the art techniques select such set to minimize the posterior variance of the Gaussian process. We explore the feasibility of solving this problem by proposing a novel Quadratic Unconstrained Binary Optimization (QUBO) model. In recent years this QUBO formulation has gained increasing attention since it represents the input for the specialized quantum annealer D-Wave machines. Hence, our contribution takes an important first step towards the sampling optimization of Gaussian processes in the context of quantum computation. Results of our empirical evaluation shows that the optimum of the QUBO objective function we derived represents a good solution for the above mentioned problem. In fact we are able to obtain comparable and in some cases better results than the widely used submodular technique.