$$Q=\frac{\Gamma\left(\frac56\right)}{\Gamma\left(\frac13\right)}\cdot\frac{\sqrt{27}+\sqrt3\,\ln\left(16+8\,\sqrt3\right)-2\pi}{\sqrt[3]{76+44\,\sqrt3}}\cdot\sqrt\pi$$
Proof:
Looking through Table of Integrals, Series, and Products, $7^{th}$ Edition, I.S. Gradshteyn, I.M. Ryzhik, I noticed that the formula 3.255 could be taken as a starting point:
$$\int_0^1\frac{x^{\mu+\frac12}(1-x)^{\mu-\frac12}}{\left(c+2\,b\,x-a\,x^2\right)^{\mu+1}}dx=\frac{\sqrt\pi}{\left(a+\left(\sqrt{c+2\,b-a}+\sqrt{c}\right)^2\right)^{\mu+\frac12}\sqrt{c+2\,b-a}}\cdot\frac{\Gamma\left(\mu+\frac12\right)}{\Gamma(\mu+1)}$$
Fixing the parameters $a=0,\ b=1,\ c=1$, we can make an observation that the integrand in $Q$ can be represented as a derivative of the integrand in 3.255:
$$\frac{x^{5/6}}{(1-x)^{1/6}\,(1+2\,x)^{4/3}}\ln\left(\frac{1+2\,x}{x\,(1-x)}\right)=\partial_\mu\left(\frac{x^{\mu+\frac12}(1-x)^{\mu-\frac12}}{(1+2\,x)^{\mu+1}}\right)_{\mu=\frac13}$$
Now we need to calculate the derivative of the right-hand side of 3.255:
$$Q=\partial_\mu\left(\frac{\sqrt\pi}{\left(\sqrt3+1\right)^{2\mu+1}\sqrt3}\cdot\frac{\Gamma\left(\mu+\frac12\right)}{\Gamma(\mu+1)}\right)_{\mu=\frac13}\\=\sqrt{\frac\pi3}\left(\frac{\psi\left(\mu+\frac12\right)-\psi(\mu+1)-2\ln\left(\sqrt3+1\right)}{\left(\sqrt3+1\right)^{2\mu+1}}\cdot\frac{\Gamma\left(\mu+\frac12\right)}{\Gamma(\mu+1)}\right)_{\mu=\frac13}\\=\sqrt{\frac\pi3}\cdot\frac{\psi\left(\frac56\right)-\psi\left(\frac43\right)-2\ln\left(\sqrt3+1\right)}{\left(\sqrt3+1\right)^{\frac53}}\cdot\frac{\Gamma\left(\frac56\right)}{\Gamma\left(\frac43\right)}$$
Here $\psi(z)=\partial_z\ln\Gamma(z)$ is the digamma function. Using Gauss digamma theorem we can expand the values of the digamma function that appear in this formula:
$$\psi\left(\frac56\right)=\frac{\pi\,\sqrt3}2-\frac{3\ln3}2-2\ln2-\gamma$$
$$\psi\left(\frac43\right)=3-\frac{\pi\,\sqrt3}6-\frac{3\ln3}2-\gamma$$
Plugging these values back to the previous formula, we get the desired result.
Addendum (By editor): The formula $3.255$ in G&R that Vladimir quoted, can be proved by 6 substitutions, that is: $$x\to 1-t, t\to \frac1u, u\to v+1, v\to w^2, w\to y \sqrt[4]{\frac a{a+b-c}}$$
and the final one is the a V. Moll's variant of Cauchy-Schlomilch transform (i.e. dealing with integrals involving $\left(\frac{x^2}{x^4+2ax^2+1}\right)^c$), which can be find in Theorem $4.1$ of:
Amdeberhan, T. , et al. "The Cauchy-Schlomilch transformation." Mathematics (2010).
This justifies the correctness of $3.255$.
We have the following closed form.
Proposition. $$
\int_0^1\log \left(1-x\right)\log \left(-\log x\right)\:dx=\gamma-\gamma_1+\gamma_1(1,0)\tag1
$$
where $\displaystyle \gamma$ is the Euler-Mascheroni constant, where $\gamma_1$ is the Stieltjes constant,
$$\gamma_1 = \lim_{N\to+\infty}\left(\sum_{n=1}^N \frac{\log n}n-\int_1^N\frac{\log t}t\:dt\right)$$
and where $\gamma_1(a,b)$ is the poly-Stieltjes constant (see here),
$$\gamma_1(a,b) = \lim_{N\to+\infty}\left(\sum_{n=1}^N \frac{\log (n+a)}{n+b}-\int_1^N\frac{\log t}t\:dt\right)\!.$$
Proof.
One may recall the classic integral representation of the Euler gamma function
$$
\frac{\Gamma(s)}{(a+1)^s}=\int_0^\infty t^{s-1} e^{-(a+1)t}\:dt, \qquad s>0,\, a>-1. \tag2
$$ By differentiating $(2)$ with respect to $s$, putting $s=1$ and making the change of variable $x=e^{-t}$, we get
$$
\int_0^1x^a\log\left(-\log x\right)\:dx=-\frac{\gamma+\log(a+1)}{a+1},\qquad a>-1, \tag3
$$
where $\displaystyle \gamma$ is the Euler-Mascheroni constant. We are allowed to insert the standard Taylor series expansion,
$$
\log (1-x)= -\sum_{n=1}^{\infty} \frac{x^n}n, \qquad |x|<1,\tag4
$$ into the given integral, then using $(3)$ we obtain
$$
\begin{align}
\int_0^1\log \left(1-x\right)\log \left(-\log x\right)\:dx&=-\int_0^1\sum_{n=1}^{\infty}\frac{x^n}n \:\log (-\log x)\:dx\\
&=-\sum_{n=1}^{\infty} \frac1n\int_0^1 x^n\log (-\log x)\:dx\\
&=\sum_{n=1}^{\infty} \frac{\gamma+\log(n+1)}{n(n+1)}\\
&=\gamma \sum_{n=1}^{\infty}\frac1{n(n+1)}+\sum_{n=1}^{\infty} \frac{\log(n+1)}{n(n+1)}\\
&=\gamma +\sum_{n=1}^{\infty} \frac{\log(n+1)}{n(n+1)},\tag5
\end{align}
$$ and we may conclude with Theorem $2$ here to get
$$
\begin{align}
\sum_{n=1}^{\infty} \frac{\log (n+1)}{n(n+1)}=\gamma_1({1,0})-\gamma_1,\tag6
\end{align}
$$
since $\gamma_1(1,1)=\gamma_1$.
Remark. I am inclined to believe that the poly-Stieltjes constants will turn out to be a tool for many of the considered integrals.
Best Answer
Here is a different way to set up that a parameter in order to apply Feynman's trick.$$I=2\int_0^\frac{\pi}{2} \ln\left(1+\sin^2t\right) dt\overset{t=\operatorname{arccot} x}=2\int_0^\infty \frac{\ln(2+x^2)-\ln(1+x^2)}{1+x^2}dx$$ Now let us consider $$I(a)=\int_0^\infty \frac{\ln(a+x^2)-\ln(1+x^2)}{1+x^2}dx\Rightarrow I'(a)=\int_0^\infty \frac{1}{(1+x^2)(a+x^2)}dx$$ $$=\frac{1}{a-1}\left(\int_0^\infty \frac{1}{x^2+1}dx-\int_0^\infty \frac{1}{x^2+a}dx\right)=\frac{1}{a-1}\left(\frac{\pi}{2}-\frac{1}{\sqrt a}\cdot \frac{\pi}{2}\right)=\frac{\pi}{2}\frac{1}{\sqrt a(1+\sqrt a)}$$ We are looking to find $I=2I(2)$, but since $I(1)=0$ we have: $$I=2\left(I(2)-I(1)\right)=2\int_1^2 I'(a)da=\pi \int_1^2 \frac{1}{\sqrt a(1+\sqrt a)}da$$ $$\overset {\large a=x^2}=2\pi\int_1^\sqrt 2 \frac{1}{1+x}dx=2\pi \ln(1+x)\bigg|_1^\sqrt 2=2\pi \ln\frac{1+\sqrt 2}{2}$$