I know how I want to answer this question. I'll write up the easy parts here, and leave the hard part for you :).

First some minor changes. It will be convenient to clear out denominators and work with $\log \left( 4 + e^{2 \pi i x} + e^{- 2 \pi i x} + e^{2 \pi i y} + e^{- 2 \pi i y} \right)$. That just changes the constant term of your Fourier series by $\log 2$. Next, it is convenient to focus on
$$ \int_0^1 \int_0^1 \log \left( 4 + e^{2 \pi i x} + e^{- 2 \pi i x} + e^{2 \pi i y} + e^{- 2 \pi i y} \right) e^{2 \pi i m x} e^{2 \pi i n y} dx dy.$$
A simple linear transformation goes between this and the cosine formulation. Let $S = \{ (z,w) : |z|=|w|=1 \}$. So we are dealing with
$$\frac{1}{(2 \pi i )^2} \int_S \log \left( 4+z+z^{-1} + w +w^{-1} \right) z^{m-1} w^{n-1} dz dw.$$
Dropping out the $4 \pi ^2$, we want to show the integrand is of the form $a \pi + b$.

**UPDATE:** Thanks to fedja for pointing out that I had oversimplified the next paragraph.

Assuming that $(m,n) \neq (0,0)$, we can integrate by parts with respect to one of the two variables, let's say $z$. Once we do that, we will have a quantity of the form
$$ (\mbox{rational number}) \cdot \int_S \frac{(z-z^{-1}) w^k z^{\ell} dw dz} {4+w+w^{-1}+z+z^{-1}}$$

So we'd like to show this quantity is of the form $a+b \pi$.

As fedja points out, we need to be careful here. Without the $z-z^{-1}$ term, the integral diverges like $\int \int ds dt/(s^2 + t^2)$ near $(-1,-1)$.

Whew! Now comes the actual hard part. Let
$$E:=\{ (z,w) \in (\mathbb{C}^*)^2 : \ 4+z+z^{-1}+w+w^{-1} =0 \}.$$
~~This is an elliptic curve with four punctures.~~ As Bjorn points out, this is a nodal cubic and can be parameterized as
$$(z,w) = \left( \frac{1-u}{u(1+u)}, \frac{u(u-1)}{1+u} \right).$$
We'll come back to this point later.

The $2$-form $dw dz/(4+z+z^{-1}+w+w^{-1})$ has a simple pole on $E$. Let $\omega$ be the $1$-form on $E$ which is the residue of that $2$-form.

I think there should be a curve $\gamma$ in $E$ such that $S$ is homotopic, in $(\mathbb{C}^*)^2 \setminus E$, to a tubular neighborhood of $\gamma$. So
$$\int_S \frac{w^k z^{\ell} (z-z^{-1}) dw dz}
{4+w+w^{-1}+z+z^{-1}} = \int_{\gamma} \omega w^k z^{\ell} (z - z^{-1}) .$$

If we substitute in the above parameterization, this will be the integral around a closed loop of some rational function in $\mathbb{Q}(u)$. In particular, we can compute this integral by residues and we will get something of the form $a+b \pi$, as desired.

Actually, it looks to me like we should just get $b \pi$. Maybe the integration by parts doesn't go as well as I hoped?

Obviously, someone should actually work this out explicitly, but I don't think it will be me.