Proposition 1 :$$\color{blue}{\displaystyle \sum\limits_{k=1}^\infty \dfrac{(-1)^{k-1}}{k}\;\zeta_H(k,a)x^k \; =\; \ln\left(\dfrac{\Gamma(a)}{\Gamma(a+x)}\right)\tag1}$$
Proof :
$$\begin{align*} & \displaystyle \sum\limits_{k=0}^{\infty}(-x)^{k}\zeta_H(k+1,a)\\ & \displaystyle = \sum\limits_{k,n=0}^{\infty} \dfrac{(-x)^k}{(n+a)^{k+1}}\\ & \displaystyle = \sum\limits_{n=0}^\infty \dfrac{1}{n+a}\sum\limits_{k=0}^\infty \left(\dfrac{-x}{n+a}\right)^k \\ & \displaystyle = \sum\limits_{n=0}^\infty \dfrac{1}{n+a+x} \\ & \displaystyle = -\psi(a+x)\end{align*}$$
Now integrating,
$$\displaystyle \sum\limits_{k=1}^{\infty} \dfrac{(-1)^{k-1}}{k}\zeta_H(k,a)x^k =\ln\left(\dfrac{\Gamma(a)}{\Gamma(a+x)}\right)\\ $$
Proposition 2 :$$\color{blue}{\displaystyle \prod\limits_{n=0}^{\infty} \left(1+\dfrac{x^2}{(n+a)^2}\right)\;=\; \dfrac{\Gamma^2 (a)}{\Gamma(a+ix)\Gamma(a-ix)}\tag 2} $$
Proof :
It is sufficient to evaluate the series,
$$\begin{align*} & \displaystyle\sum\limits_{n=0}^{\infty}\ln\left(1+\dfrac{x}{n+a}\right) \\ & = \displaystyle \sum\limits_{k=1}^\infty\sum\limits_{n=0}^\infty \dfrac{(-1)^{k-1}}{k}x^k \dfrac{1}{(n+a)^k} \\ & =\displaystyle \sum\limits_{k=1}^{\infty} \dfrac{(-1)^{k-1}}{k}\zeta_H(k,a)x^k \\ & =\displaystyle \ln\left(\dfrac{\Gamma(a)}{\Gamma(a+x)}\right)\end{align*}$$
$$\displaystyle \prod\limits_{n=0}^{\infty} \left(1+\dfrac{x}{n+a}\right) = \dfrac{\Gamma(a)}{\Gamma(a+x)}$$
Therefore,
$$\displaystyle \prod\limits_{n=0}^{\infty} \left(1+\dfrac{x^2}{(n+a)^2}\right)\;=\; \dfrac{\Gamma^2 (a)}{\Gamma(a+ix)\Gamma(a-ix)}$$
Proposition 3 : If $\; \displaystyle F(s)=\int\limits_0^\infty x^{s-1}f(x)\; dx\; $ then $$\color{blue}{\displaystyle \int\limits_{-\infty}^{\infty} |F(ix)|^2 \; dx \;= \; 2\pi\int\limits_0^\infty \dfrac{|f(x)|^2}{x}\; dx\tag 3} $$
Proof :
$$\displaystyle F(it)=\int\limits_0^\infty x^{it}\dfrac{f(x)}{x}\; dx$$
Set $x=e^y$ ,
$$\displaystyle F(it)=\int\limits_0^\infty e^{ixt}f(e^x)\; dx$$
Now by properties of Fourier Transform,
$$\begin{align*} & \displaystyle f(e^t)=\int\limits_{-\infty}^\infty g(x)e^{-ixt}\; dx \\ & \displaystyle g(t)=\dfrac{1}{2\pi}\int\limits_{-\infty}^\infty f(e^x)e^{ixt}\; dx\\\end{align*}$$
$$\begin{align*}\displaystyle F(it) & =2\pi g(t) \\ \displaystyle \int\limits_{-\infty}^\infty |F(it)|^2\; dt & = 4\pi^2 \int\limits_{-\infty}^\infty |g(t)|^2\; dt \\ \displaystyle \int\limits_{-\infty}^\infty |F(it)|^2\; dt & = 2\pi \int\limits_{-\infty}^\infty g(t) \int\limits_{-\infty}^\infty e^{ixt}f(e^x)\; dx\; dt \\ \displaystyle \int\limits_{-\infty}^\infty |F(it)|^2\; dt & = 2\pi \int\limits_{-\infty}^\infty f(e^x) \int\limits_{-\infty}^\infty e^{ixt}g(t)\; dt\; dx \\ \displaystyle \int\limits_{-\infty}^\infty |F(it)|^2\; dt & = 2\pi \int\limits_{-\infty}^\infty f(e^x)\overline{f(e^x)}\; dx =2\pi\int\limits_{-\infty}^\infty |f(e^x)|^2\; dx\end{align*}$$
Now by setting $e^x=t$ we get our result,
$$\displaystyle \int\limits_{-\infty}^\infty |F(it)|^2\; dt = 2\pi\int\limits_{-\infty}^\infty \dfrac{|f(t)|^2}{t}\; dt$$
Main Problem: $$\color{blue}{\displaystyle \int\limits_{0}^{\infty}\frac {1+\dfrac {x^2}{(b+1)^2}}{1+\dfrac {x^2}{a^2}}\dfrac {1+\dfrac {x^2}{(b+2)^2}}{1+\dfrac {x^2}{(a+1)^2}}\dfrac {1+\dfrac {x^2}{(b+3)^2}}{1+\dfrac {x^2}{(a+2)^2}}\cdots \, dx=\dfrac {\sqrt{\pi}}2\dfrac {\Gamma\left(a+\frac 12\right)\Gamma\left(b+1\right)\Gamma(b-a+1)}{\Gamma(a)\Gamma\left(b+\frac 12\right)\Gamma\left(b-a+\frac 12\right)}} $$
Proof : If we denote the integral by $I$ then using $(2)$ it can be rewritten as,
$$\displaystyle I = \dfrac{\Gamma^2 (b+1)}{\Gamma^2 (a)}\dfrac{1}{2} \int\limits_{-\infty}^\infty \dfrac{|\Gamma(a+ix)|^2}{|\Gamma(b+1+ix)|^2}\; dx$$
Now by defining $\displaystyle h(x)=\dfrac{x^a (1-x)^{b-a}}{\Gamma(b-a+1)}$ for $x\in[0,1]$ and $0$ for $\forall x\notin[0,1]$ we can conclude that $\displaystyle F(s)=M[h(x)]=\dfrac{\Gamma(s+a)}{\Gamma(s+b+1)}$ and from $(3)$ it follows that,
$$\displaystyle I = \dfrac{\Gamma^2 (b+1)}{\Gamma^2 (a)}\dfrac{1}{2} \int\limits_{0}^1 \dfrac{|h(x)|^2}{x}\; dx = \dfrac {\sqrt{\pi}}2\dfrac {\Gamma\left(a+\frac 12\right)\Gamma\left(b+1\right)\Gamma(b-a+1)}{\Gamma(a)\Gamma\left(b+\frac 12\right)\Gamma\left(b-a+\frac 12\right)}$$
where last line follows from the Duplication formula , and we are done !
$$\large \color{red}{\color{blue}{\boxed{\mathfrak{PROVED}}}} $$
Best Answer
You can use residues, let $f(z)=\frac{1}{(1+z^2)(1+r^2z^2)(1+r^4z^2)...}$ this has an infite set of singularities at $z=\pm i 1/r^{n},n\in\Bbb N_0$. We can see that $\large\int_0^\infty\frac{1}{(1+x^2)(1+r^2x^2)(1+r^4x^2)...}dx=\frac{1}{2}\int_{-\infty}^\infty\frac{1}{(1+x^2)(1+r^2x^2)(1+r^4x^2)...}dx$
So we will consider the upper semi circle contour, spanning over $-\infty$ to $+\infty$, only the positive singularities are in this contour, i.e. $z=i\frac{1}{r}$.
Now $\operatorname{Res}(f(z),z=\frac{i}{r^{n}})=\large\lim_{z\to i\frac{1}{r^{n}}}(z-i\frac{1}{r^{n}})f(z)=\large\lim_{z\to i\frac{1}{r^{n}}}\frac{1}{r^n}(r^nz-i)f(z)$
$=\lim_{z\to i\frac{1}{r^{n}}}\frac{1}{r^n(r^nz+i)}\prod_{j=0,j\ne n}^\infty\frac{1}{(1+x^{2}r^{2j})}$
$=\large\frac{1}{2ir^n}\prod_{j=0,j\ne n}^\infty\frac{1}{(1-r^{2j-2n})}$
Now $\frac{1}{2}\int_{-\infty}^\infty\frac{1}{(1+x^2)(1+r^2x^2)(1+r^4x^2)...}dx=\frac{1}{2}2\pi i\sum \operatorname{Res}(f)=\frac{1}{2}\pi i\sum_{n=0}^\infty\frac{1}{ir^n}\prod_{j=0,j\ne n}^\infty\frac{1}{(1-r^{2j-2n})}$
From here I'm not sure how to reach the closed form, but hopefully this helps to show a different approach, even if you are not too well versed in integrals by residue :)