This is a cute problem! I toyed with it and didn't really get anywhere - I got the strong impression that it requires fields of mathematics that I am not expert in.

Indeed, given that the problem seems related to that of counting integer solutions to the equation $f(x,y) = c$, one may need to use arithmetic geometry tools (e.g. Faltings' theorem). In particular if we could reduce to the case when the genus is just 0 or 1 then presumably one could kill off the problem. (One appealing feature of this approach is that arithmetic geometry quantities such as the genus are automatically invariant (I think) with respect to invertible polynomial changes of variable such as $(x,y) \mapsto (x,y+P(x))$ or $(x,y) \mapsto (x+Q(y),y)$ and so seem to be well adapted to the problem at hand, whereas arguments based on the raw degree of the polynomial might not be.)

Of course, Faltings' theorem is ineffective, and so might not be directly usable, but perhaps some variant of it (particularly concerning the dependence on c) could be helpful. [Also, it is overkill - it controls rational solutions, and we only care here about integer ones.] This is far outside of my own area of expertise, though...

The other thing that occurred to me is that for fixed c and large x, y, one can invert the equation $f(x,y) = c$ to obtain a Puiseux series expansion for y in terms of x or vice versa (this seems related to resolution of singularities at infinity, though again I am not an expert on that topic; certainly Newton polytopes seem to be involved). In some cases (if the exponents in this series expansion are favourable) one could then use Archimedean counting arguments to show that f cannot cover all the natural numbers (this is a generalisation of the easy counting argument that shows that a 1D polynomial of degree 2 or more cannot cover a positive density set of integers), but this does not seem to work in all cases, and one may also have to use some p-adic machinery to handle the other cases. One argument against this approach though is that it does not seem to behave well with respect to invertible polynomial changes of variable, unless one works a lot with geometrical invariants.

Anyway, to summarise, it seems to me that one has to break out the arithmetic geometry and algebraic geometry tools. (Real algebraic geometry may also be needed, in order to fully exploit the positivity, though it is also possible that positivity is largely a red herring, needed to finish off the low genus case, but not necessary for high genus, except perhaps to ensure that certain key exponents are even.)

EDIT: It occurred to me that the polynomial $f(x,y)-c$ might not be irreducible, so there may be multiple components to the associated algebraic curve, each with a different genus, but presumably this is something one can deal with. Also, the geometry of this curve may degenerate for special c, but is presumably stable for "generic" c (or maybe even all but finitely many c).

It also occurs to me that one use of real algebraic geometry here is to try to express f as something like a sum of squares. If there are at least two nontrivial squares in such a representation, then f is only small when both of the square factors are small, which is a 0-dimensional set and so one may then be able to use counting arguments to conclude that one does not have enough space to cover all the natural numbers (provided that the factors are sufficiently "nonlinear"; if for instance $f(x,y)=x^2+y^2$ then the counting arguments barely fail to provide an obstruction, one has to use mod p arguments or something to finish it off...)

EDIT, FOUR YEARS LATER: OK, now I know a bit more arithmetic geometry and can add to some of my previous statements. Firstly, it's not Faltings' theorem that is the most relevant, but rather Siegel's theorem on integer points on curves - the enemy appears to be those points $(x,y)$ where $x,y$ are far larger than $f(x,y)$, and Siegel's theorem is one of the few tools available to exclude this case. The known proofs of this theorem are based on two families of results in Diophantine geometry: one is the Thue-Siegel-Roth theorem and its variants (particularly the subspace theorem), and the other is the Mordell-Weil theorem and its variants (particularly the Chevalley-Weil theorem). A big problem here is that all of these theorems have a lot of ineffectivity in them. Even for the very concrete case of Hall's conjecture on lower bounding $|x^2-y^3|$ for integers $x,y$ with $x^2 \neq y^3$, Siegel's theorem implies that this bound goes to infinity as $x,y \to \infty$, but provides no rate; as I understand it, the only known lower bounds are logarithmic and come from variants of Baker's method.

As such, a polynomial such as $f(x,y) = (x^2 - y^3 - y)^4 - y + C$ for some large constant C already looks very tough to analyse. (I've shifted $y^3$ here by $y$ to avoid the degenerate solutions to $x^2=y^3$, and to avoid some cheap way to deal with this polynomial from the abc conjecture or something.) The analogue of Hall's conjecture for $|x^2-y^3-y|$ suggests that $f(x,y)$ goes to $+\infty$ as $x,y \to \infty$ (restricting $x,y$ to be integers), but we have no known growth rate here due to all the ineffectivity. As such, we can't unconditionally rule out the possibility of an infinite number of very large pairs $(x,y)$ for which $x^2-y^3-y$ happens to be so close to $y^{1/4}$ that we manage to hit every positive integer value in $f(x,y)$ without hitting any negative ones. However, one may be able to get a conditional result assuming some sufficiently strong variant of the abc conjecture. One should also be able to exclude large classes of polynomials $f$ from working; for instance, if the curve $f(x,y)=0$ meets the line at infinity at a lot of points in a transverse manner, then it seems that the subspace theorem may be able to get polynomial bounds on solutions $(x,y)$ to $f(x,y)=c$ in terms of $c$, at which point a lot of other tools (e.g. equidistribution theory) become available.

Another minor addendum to my previous remarks: the generic irreducibility of $f(x,y)-c$ follows from Bertini's second theorem, as one may easily reduce to the case when $f$ is non-composite (not the composition of two polynomials of lower degree).

1more comment