Here's a trick (which sometimes works), but it might be contained in the comments above. As stated, the Newton polyhedron (convex hull of the exponents, the latter viewed as lattice points in ${\bf R}^2$), $K$, of $P$, is contained in the double unit square in ${\bf R}^2$. First, an obvious necessary condition for nonnegativity: the coefficients of the extreme points of $K$ must be positive (as opposed to negative).

Next, for each edge of $K$, there is a corresponding polynomial (just take the exponents appearing in the edge, and use the coefficients from $P$); if $P$ is nonnegative, each of these edge polynomials must be nonnegative---but after a change of variables, the edge polynomials are really just polynomials in one variable, for which testing for nonnegativity is just the quadratic formula in your example).

This gives some necessary conditions.

Now a sufficient condition for positivity:
If $P(x,y) >$ on the interior of the positive orthant, and each of the edge polynomials is positive (which implies the coefficients at the extreme points of $K$ are all positive), then for any $Q$ with the same supporting monomials as $P$ and $Q$ having no nonnegative coefficients, there exists an integer $M$ such that the product, $Q^M P$, has no negative coefficients. (This is a generalization of Polya and Meissner's theorems; when $K$ is the double unit square, this is Meissner's theorem.)

The simplest (but probably not the best) choice for $Q$ is $\sum x^w$ where $w$ varies over the lattice points in $K$ (where $x^w = x_1^{w(1)}x_2^{w(2)}$ is the usual monomial, with $w = (w(1),w(2))$. Another choice is to replace the negative coefficients in $P$ by their absolute value. In any event, you can try testing $Q^m P$ for large (whatever that means in this context) choices of $m$.

And a variation: if you suspect your thing is only nonnegative, add $\epsilon Q = \epsilon \cdot \sum_{w \in K} x^w$ to $P$, where $0 <\epsilon$ is much smaller than all the absolute values of the nonzero coefficients of $P$, and test $Q^m \cdot (P + \epsilon Q)$ for nonnegative coefficients. If succcessful, reduce $\epsilon$. If unsuccessful, increase $\epsilon$; this might give estimates of lower bounds.

The problem, as usual, is that it can't give a negative answer, although it might suggest where to look.