<abstract><p>It is known that an efficient method for interpolation of very large scattered data sets is the method of Shepard. Unfortunately, it reproduces only the constants. In this paper, we first generalize an expansion in bivariate even order Bernoulli polynomials for real functions possessing a sufficient number of derivatives. Finally, by combining the known Shepard operator with the even order Bernoulli bivariate operator, we construct a kind of new approximated operator satisfying the higher order polynomial reproducibility. We study this combined operator and give some error bounds in terms of the modulus of continuity of high order and also with Peano's theorem. Numerical comparisons show that this new technique provides the higher degree of accuracy. Furthermore, the advantage of our method is that the algorithm is very simple and easy to implement.</p></abstract>