With a little time on my hands these past couple of days, I have finally succeeded in computing this integral through contour integration; the method is very similar in spirit to my answer here. The only background needed is knowledge of Cauchy's integral theorem and the fact that $\zeta(4) = \pi^4/90$. I've made a few remarks on generalization towards the bottom.
The strategy is as follows. We integrate the principal branch of $f(z) = z\log^3{(1+e^{2i z})}$ over an appropriately chosen contour in order to prove
$$
\begin{align}
\int_{-\pi/2}^{\pi/2}x^2\log^2{(2\cos{x})}\,dx & = \int_{-\pi/2}^{\pi/2}x^4\,dx-\frac{\pi}{3}\int_0^\infty \log^3{(1-e^{-2y})}\,dy. \tag{1}
\end{align}
$$
This leaves us in good shape; due to the fact that the integrand is even, the left-hand side is twice the integral we wish to evaluate. The first integral on the right is easily computed as $\pi^5/80$. To evaluate the second integral, put $e^{-v} = 1-e^{-2y}$; mild computations give $-(e^v - 1)^{-1}\,dv = dy$, and then
$$
\begin{align}
-\frac{\pi}{3}\int_0^\infty \log^3{(1-e^{-2y})}\,dy & = \frac{\pi}{6}\int_0^\infty\frac{v^3}{e^v - 1}\,dv = \pi\zeta(4). \tag{2}
\end{align}
$$
The last equation follows from the well-known identity
$$
\Gamma(s)\zeta(s) = \int_0^\infty \frac{v^{s-1}}{e^{v}-1}\,dv
$$
which holds for $\text{Re}\,s > 1$ and can be derived easily by developing the integrand in a geometric series. Upon inserting $(2)$ into $(1)$ and simplifying we arrive at the desired result.
It remains to prove $(1)$; as mentioned, the argument is entirely analogous to one I gave in this answer to another of Cody's questions. For completeness, I include it here. Consider the region obtained from $\mathbb C$ by removing the half-lines on which $\text{Re}\,z$ is an integer multiple of $\pi/2$ and $\text{Im}\,z \leq 0$. On this region, we can define a branch of $\log{(1+e^{2iz})} = \log{(2e^{iz}\cos{z})}$; we choose the branch with imaginary part between $-\pi$ and $\pi$. Having done this, let $f(z) = z\log^3{(1+e^{2iz})}$. We wish to integrate $f(z)$ over the contour obtained by removing the corners from the rectangle determined by $-\pi/2$ and $\pi/2 + iR$, where $R > 0$. To complete the contour we replace the bottom corners with quarter-circles of radius $\delta >0$. The intention is to let $R \to \infty$ and $\delta \to 0$.
${}$
For fixed values of $\delta$ and $R$, Cauchy's theorem says that $f(z)$ integrates to zero over the contour. As $R \to \infty$, the contribution from the upper horizontal side tends to $0$ because $f(x+iR) \to 0$ uniformly for $-\pi/2 \leq x \leq \pi/2$. By writing $1+e^{2i z} = 1- e^{2i(z-\pi/2)}$ one sees that $1+e^{2i z} = O(z-\pi/2)$, and therefore that $\log{(1+e^{2iz})} = O(\log{|z-\pi/2|})$ as $z \to \pi/2$. It follows that the contribution from the left quarter-circle is $O(\delta^2\log^3{\delta})$, hence that it vanishes in the limit $\delta \to 0$. The same argument applies to the other quarter-circle, and we are left with the contributions from the vertical sides and the lower horizontal side.
After taking limits, the contribution from the vertical sides is
$$
\begin{align}
i\int_0^\infty f(iy+\pi/2)\,dy -i\int_0^\infty f(iy - \pi/2)\,dy = i\pi\int_0^\infty\log^3{(1-e^{-2y})}\,dy. \tag{3}
\end{align}
$$
Now for $x$ between $-\pi/2$ and $\pi/2$, the quantity $2\cos{x}$ is positive. This means that the unique value of $\arg{(2e^{ix}\cos{x})}$ which lies between $-\pi$ and $\pi$ is simply $x$. As we have chosen the principle branch, we obtain $\log{(2e^{ix}\cos{x})} = \log{(2\cos{x})}+ix$. Therefore the contribution from the lower horizontal side may be written
$$
\int_{-\pi/2}^{\pi/2} f(x)\,dx = \int_{-\pi/2}^{\pi/2}x\left(\log{(2\cos{x})} + ix\right)^3\,dx. \tag{4}
$$
By the preceding analysis, $(3)$ and $(4)$ sum to zero. Since $(3)$ is purely imaginary, this means in particular that the imaginary part of $(4)$ is equal to the negative of $(3)$. The last statement is equivalent to $(1)$.
Because I have spent quite some time with this question, I would like to make a few remarks in the direction of a generalization. By taking $g(z) = p(z) \log^m(1+e^{2iz})$ in place of $f(z)$ above, where $p(z)$ is a polynomial and $m \in \mathbb N$, one finds by repeating the same arguments that
$$
\begin{align}
\int_{-\pi/2}^{\pi/2} p(x)&\left(\log{(2\cos{x})} + ix\right)^m\,dx \\
&= -i \int_0^\infty \left(p(iy + \pi/2)-p(iy - \pi/2)\right)\log^m{(1-e^{-2y})}\,dy. \tag{5}
\end{align}
$$
Clever choices of $p$ and $m$ then furnish a number of integral and series identities. Moreover, the identities
$$
\int_0^\infty y^n \log{(1-e^{-2y})}\,dy = -\frac{1}{2^{n+1}}\Gamma(n+1)\zeta(n+2) \tag{6}
$$
and
$$
\int_0^\infty \log^m{(1-e^{-2y})}\,dy = \frac{(-1)^{m}}{2}\int_0^\infty \frac{y^m}{e^y - 1}\,dy = \frac{(-1)^m}{2}\Gamma(m+1)\zeta(m+1) \tag{7}
$$
provide a connection between the integrals on the left-hand side of $(5)$ and the values taken by $\zeta$ at positive integers $n\geq2$. Here $(6)$ is derived by expanding $\log{(1-e^{-2y})}$ in powers of $e^{-2y}$ while $(7)$ is derived in the same way as $(2)$.
