The integral can be approached by using Mellin convolution technique. Let's first start with a simple case of $a=b$. In this case we compute Mellin transform of $J_1(x)^2$ and $J_0(x) J_1(x)$: $ \mathcal{M}_s( J_1(x)^2 ) = \int_0^\infty x^{s} J_1(x)^2 \frac{\mathrm{d} x}{x} = \frac{\Gamma \left(\frac{1}{2}-\frac{s}{2}\right) \Gamma \left(\frac{s}{2}+1\right)}{2 \sqrt{\pi } \Gamma \left(1-\frac{s}{2}\right) \Gamma \left(2-\frac{s}{2}\right)} \qquad \text{for} \qquad -2 < \mathrm{Re}(s) < 1 $ and $ \mathcal{M}_s( J_1(x) J_0(x) ) = \frac{\Gamma \left(1-\frac{s}{2}\right) \Gamma \left(\frac{s+1}{2}\right)}{2 \sqrt{\pi } \Gamma \left(\frac{3}{2}-\frac{s}{2}\right)^2} \qquad \text{for} \qquad -1 < \mathrm{Re}(s) < 2 $ Now $ \begin{eqnarray} \int_0^\infty x^{-2} J_1(a x)^3 J_0(a x) \mathrm{d} x &=& \frac{a}{2} \frac{1}{2 \pi i} \int_{-1 - i \infty}^{-1 + i \infty} \frac{\Gamma \left(\frac{1}{2}-\frac{s}{2}\right) \Gamma \left(\frac{s}{2}+1\right) \Gamma \left(\frac{s}{2}+\frac{3}{2}\right) \Gamma \left(-\frac{s}{2}\right)}{2 \pi \Gamma \left(1-\frac{s}{2}\right) \Gamma \left(2-\frac{s}{2}\right) \Gamma \left(\frac{s}{2}+2\right)^2} \mathrm{d} s \\ &=& \frac{a}{2 \pi i} \int_{-\frac{1}{2} - i \infty}^{-\frac{1}{2} + i \infty} \frac{\Gamma \left(\frac{1}{2}-s\right) \Gamma \left(s+1\right) \Gamma \left(s+\frac{3}{2}\right) \Gamma \left(-s\right)}{2 \pi \Gamma \left(1-s\right) \Gamma \left(2-s\right) \Gamma \left(s+2\right)^2} \mathrm{d} s \\ &=& \frac{a}{2 \pi} G_{4,4}^{2,2}\left(1\left| \begin{array}{c} \frac{1}{2},1,2,2 \\ 1,\frac{3}{2},0,-1 \\ \end{array} \right.\right) = a \left( \frac{3 }{16} - \frac{1}{\pi^2} \right) \end{eqnarray} $
Here $G_{4,4}^{2,2}(1)$ denotes Meijer's G-function.
Now when $a \not= b$, Mellin transform of $J_1(a x) J_0(b x)$ is no longer a ratio of $\Gamma$-functions: $ \left. \mathcal{M}_s( J_1( a x) J_0(b x)) \right\vert_{-1 < \mathrm{Re}(s) <2} = \left\{ \begin{array}{cc} \frac{a 2^{s-1} b^{-s-1} \Gamma \left(\frac{s+1}{2}\right) \, _2F_1\left(\frac{s+1}{2},\frac{s+1}{2};2;\frac{a^2}{b^2}\right)}{\Gamma \left(\frac{1}{2}-\frac{s}{2}\right)} & a < b \\ \frac{2^{s-1} a^{-s} \Gamma \left(\frac{s+1}{2}\right) \, _2F_1\left(\frac{s-1}{2},\frac{s+1}{2};1;\frac{b^2}{a^2}\right)}{\Gamma \left(\frac{3}{2}-\frac{s}{2}\right)} & a > b \end{array} \right. $ Notice that Gauss hypergeometric function $ {}_2F_1$ can be represented by its defining sum, whose terms are ratios of Gamma functions themselves. Continuing this way will produce sum representation of the integral.
Integrals of this form has been discussed by W.N. Bailey in 1936, see link.