That integral denotes that you want to integrate $1/|\nabla \epsilon(k)|$ over all wavevectors that satisfy $\epsilon(k) = \epsilon$ for any given $\epsilon$ which is taken to be a constant. Unless you find some convenient parametrization of these constant-energy surfaces, this is a hard to do analytically and one needs to numerical techniques.
Now, for the problem at hand, I do not have a full solution nor a good intuition what the result should be. But here is an outline which should work: First note $\nabla \epsilon(k) = (2 k_x,-2 k_y)$ so $|\nabla \varepsilon(k)| = 4 \sqrt{k_x^2+k_y^2}$.
We then have $$ D(\epsilon) = \int_{\epsilon(k_x,k_y)=\epsilon} d S \frac{1}{4 \sqrt{k_x^2+ k_y^2}}$$.
The form of the dispersion suggests to use a hyperbolic parametrization. Letting $k_x = a \cosh u$ and $k_y = a \sinh u$ which parametrizes the region, $k_x > 0$ and $-k_x < k_y < k_x$, we note that $\varepsilon(k) = a^2 (\cosh^2 u - \sinh^2 u) = a^2$ which implies we can set $a=\sqrt{\epsilon}$ and only integrate over $u$, where the surface element is given by $a du$ (this can be checked from the Jacobian).
You then need the integral $\int_{-\infty}^\infty 1/\sqrt{\cosh 2u} du $ which can be done with WolframAlpha, yielding $\approx 2.62$. Note that this is only one contribution to the integral, by appropriate parametrizations of $k_x$ and $k_y$ you should be able to cover the entire plane. Basically we are integrating $1/|k|$ along a level set of $k_x^2 - k_y^2$ as illustrated by the contour lines in this plot.
Perhaps it is easier to just work with polar coordinates where the integrand is trivial, but the limits of integration will not be straightforward anymore.
