If you restrict x, say to be nonnegative, then this is a straightforward optimization problem once the integral is made into a "black box" function that only evaluates for numeric input.
f[x_?NumberQ, y_?NumberQ, r_?NumberQ] :=
NIntegrate[
1/(E^((2*x^2 - 4*x*d + 2*d^2 + 2*y^2)/(8*d + 4*r^2))*(Sqrt[
d]*(2*d + r^2))), {d, 0, Infinity}] - 1
Example:
FindMaximum[{y, f[x, y, 1/2] == 0, x >= 0}, {x, y}]
(* Out[69]= {1.26438464527, {x -> 0.340246814653, y -> 1.26438464527}} *)
