$\newcommand{\d}{\,\mathrm{d}}\newcommand{\sech}{\operatorname{sech}}\newcommand{\coth}{\operatorname{coth}}\newcommand{\csch}{\operatorname{csch}}\newcommand{\res}{\operatorname{Res}}$Surprisingly, this is actually doable.
$$\begin{align}\forall q<2,&\forall\beta>0:\\\quad\quad\int_{-\infty}^\infty\int_{-\infty}^t(t-s)^{1-q}&\tanh(\beta s)\sech^2(\beta s)\sech^2(\beta t)\d s\d t\\&=\frac{q(q^2-1)}{\pi^q\cdot\beta^{3-q}}\cdot\sec\frac{\pi q}{2}\cdot\zeta(q+1)\\&=\frac{2^{1+q}\cdot q(q^2-1)}{\beta^{3-q}}\cdot\Gamma(-q)\zeta(-q)\end{align}$$
Let's attack the inner integral first. Fix $0<\beta$ and $1<q<2$.
Fix $t\in\Bbb R$. Task:
Evaluate: $$J(t):=\int_{-\infty}^t(t-s)^{1-q}\tanh(\beta s)\sech^2(\beta s)\d s$$
Solution:
Write $\alpha:=1-q\in(-1,0)$. Define a logarithm $\log:\Bbb C^\star\to\Bbb C$ via $0\le\arg<2\pi$.
If $\Omega$ stands for $\Bbb C$ takeaway $t$ and all the zeroes of $\cosh(\beta(-))$, we define $f:\Omega\to\Bbb C$ via: $$s\mapsto\exp(\alpha\log(t-s))\tanh(\beta s)\sech^2(\beta s)$$$f$ will be holomorphic away from $(-\infty,t]$. Fixing $s\in(-\infty,t]$, we have that $f(z)\to(t-s)^\alpha\tanh(\beta s)\sech^2(\beta s)$ as $z\to s$ from below the real axis but $f(z)\to\exp(2\pi i\alpha)(t-s)^{\alpha}\tanh(\beta s)\sech^2(\beta s)$ as $z\to s$ from above the real axis.
Because of the bounds on $\alpha$, it is possible to conclude, with a little care needed to check the necessary error terms vanish, that: $$\oint_{\mathscr{H}}f(z)\d z=2\pi i\sum_{n\in\Bbb Z}\res\left(f;\frac{\pi i}{2\beta}(2n+1)\right)$$Where $\mathscr{H}$ is a Hankel contour that winds $-\infty+i0^+\to t+i0^+\to t+i0^-\to-\infty+i0^-$.
Equivalently: $$\oint_{\mathscr{H}}f(z)\d z=(\exp(2\pi i\alpha)-1)\cdot J(t)$$So that: $$J(t)=\pi\exp(-\pi i\alpha)\cdot\csc(\pi\alpha)\cdot\sum_{n\in\Bbb Z}\res\left(f;\frac{\pi i}{2\beta}(2n+1)\right)$$
Now for the residue computation. Fix $n\in\Bbb Z$. The pole here is of third order. It's easy but boring to compute that: $$\tanh\left(\frac{\pi i}{2}(2n+1)+\beta z\right)\sech^2\left(\frac{\pi i}{2}(2n+1)+\beta z\right)\\=-\coth(z)\csch^2(z)=-\beta^{-3}z^{-3}+o(1),\,z\to0$$And that, for nonnegative $n$: $$\exp\left(\alpha\log\left(t-\frac{\pi i}{2\beta}(2n+1)-z\right)\right)\\=\exp(\pi i\alpha)\cdot\exp\left(\alpha\log\left(\frac{\pi i}{2\beta}(2n+1)-t\right)\right)\\\times\left(1-\frac{\alpha z}{t-\frac{\pi i}{2\beta}(2n+1)}+\frac{1}{2}\frac{\alpha(\alpha-1)}{\left(t-\frac{\pi i}{2\beta}(2n+1)\right)^2}z^2+o(z^2)\right)$$And for negative $n$: $$\exp\left(\alpha\log\left(t-\frac{\pi i}{2\beta}(2n+1)-z\right)\right)\\=\exp(-\pi i\alpha)\cdot\exp\left(\alpha\log\left(\frac{\pi i}{2\beta}(2n+1)-t\right)\right)\\\times\left(1-\frac{\alpha z}{t-\frac{\pi i}{2\beta}(2n+1)}+\frac{1}{2}\frac{\alpha(\alpha-1)}{\left(t-\frac{\pi i}{2\beta}(2n+1)\right)^2}z^2+o(z^2)\right)$$One must take care in manipulating the logarithm as $\log(ab)=\log(a)+\log(b)$ is not always true (this is the reason the answer is different for different signs of $n$). We conclude the residue is: $$\begin{cases}-\frac{\alpha(\alpha-1)}{2\beta^3\cdot\left(t-\frac{\pi i}{2\beta}(2n+1)\right)^2}\exp(\pi i\alpha)\exp\left(\alpha\log\left(\frac{\pi i}{2\beta}(2n+1)-t\right)\right)&n\ge0\\-\frac{\alpha(\alpha-1)}{2\beta^3\cdot\left(t-\frac{\pi i}{2\beta}(2n+1)\right)^2}\exp(-\pi i\alpha)\exp\left(\alpha\log\left(\frac{\pi i}{2\beta}(2n+1)-t\right)\right)&n\le-1\end{cases}$$Concluding that: $$\small{J(t)=-\exp(-\pi i\alpha)\cdot\frac{\pi q(q-1)}{2\beta^3}\csc(\pi q)\times\sum_{n\ge0}\\\left(\exp(\pi i\alpha)\cdot\frac{\exp\left(\alpha\log\left(\frac{\pi i}{2\beta}(2n+1)-t\right)\right)}{\left(t-\frac{\pi i}{2\beta}(2n+1)\right)^2}+\exp(-\pi i\alpha)\cdot\frac{\exp\left(\alpha\log\left(\frac{\pi i}{2\beta}(2(-n-1)+1)-t\right)\right)}{\left(t-\frac{\pi i}{2\beta}(2(-n-1)+1)\right)^2}\right)}$$
Exploiting symmetry (paying careful attention to the definition of $\log$) this can be rewritten as: $$J(t)=-\frac{\pi q(q-1)}{\beta^3}\csc(\pi q)\cdot\Re\sum_{n\ge0}\frac{\exp\left(\alpha\log\left(\frac{\pi i}{2\beta}(2n+1)-t\right)\right)}{\left(t-\frac{\pi i}{2\beta}(2n+1)\right)^2}$$
As an aside, let's look at this for $t=0$. We find: $$\begin{align}J(0)&=\int_{-\infty}^0(-s)^{1-q}\tanh(\beta s)\sech^2(\beta s)\d s\\&=-\frac{\pi q(q-1)}{\beta^3}\csc(\pi q)\cdot\Re\sum_{n\ge0}\frac{\exp(\pi i\alpha/2)\left(\frac{\pi}{2\beta}\right)^{\alpha}(2n+1)^{\alpha}}{-\frac{\pi^2}{4\beta^2}(2n+1)^2}\\&=-\frac{\pi^{\alpha+1}q(q-1)}{2^{1+\alpha}\beta^{3+\alpha}}\sec(\pi q/2)(-1)\frac{4\beta^2}{\pi^2}\sum_{n\ge0}\frac{1}{(2n+1)^{1+q}}\\&=(2^{q+1}-1)\frac{q(q-1)}{2\pi^q\cdot\beta^{2-q}}\cdot\sec\left(\frac{\pi q}{2}\right)\cdot\zeta(q+1)\end{align}$$
Which agrees with the numerics, I am pleased to say! It is notable how similar the closed forms for $J(0)$ and for the overall integral are. Morally, I suppose this should be because $\sech^2(\beta t)$ rapidly decays as $t$ moves away from zero so that it is "mostly $J(0)$" which contributes to the value of the integral.
I find it highly unlikely a better closed form exists for $J(t)$ when $t\neq0$. Nevertheless we can use this formula to evaluate the required integral.
We care about: $$\begin{align}J:&=\int_{-\infty}^\infty\sech^2(\beta t)J(t)\d t\\&=-\pi q(q-1)\csc(\pi q)\cdot\beta^{q-2}\\&\times\Re\int_{-\infty}^\infty\sech^2(\beta t)\left(\sum_{n\ge0}\exp\left(-(1+q)\log\left(\frac{\pi i}{2}(2n+1)-\beta t\right)\right)\right)\d t\\&=-\frac{\pi q(q-1)}{\beta^{3-q}}\csc(\pi q)\\&\times\Re\int_{-\infty}^\infty\sech^2(t)\underset{\large\Phi(t)}{\underbrace{\sum_{n\ge0}\exp\left(-(1+q)\log\left(\frac{\pi i}{2}(2n+1)-t\right)\right)}}\d t\end{align}$$
The series - let's call it $\Phi$ - is holomorphic (viewed as a function in $t$) in the lower half plane: I justify that since the summands fall off with $n$ like $n^{-(1+q)}<n^{-2}$ so everything should locally uniformly converge, and in the lower half plane we avoid the singularities due to the logarithms. Because $1+q>1$ I believe it can be justified (the $\sech^2$ term will be extremely small except for points near the imaginary axis, and using a similar argument to the one I used at the start one can also make that issue disappear), with care, that the integral of $\Phi(\cdot)\sech^2(\cdot)$ around large semicircles in the lower half plane would vanish as the radius tends to $\infty$.
Let's denote $\int_{-\infty}^\infty\Phi(t)\sech^2(t)\d t$ by $K$. We should have: $$K=-2\pi i\cdot\sum_{m\ge0}\res\left(\Phi\cdot\sech^2;-\frac{\pi i}{2}(2m+1)\right)$$It's now time to do another residue calculation. This one is less tedious than the first due to some very lucky symmetry.
Fixing $m$: $$\sech^2\left(-\frac{\pi i}{2}(2m+1)+z\right)=-\csch^2(z)=-z^{-2}+o(1),\,z\to0$$
But what does $\Psi_m(z):=\Phi\left(-\frac{\pi i}{2}(2m+1)+z\right)$ behave like as $z\to0$? To identify the residue - noting the companion $z^{-2}$ term - we need only identify the $[z]$ coefficient in the asymptotic expansion of $\Psi_m(z)$. That is: $$\begin{align}K&=2\pi i\cdot\sum_{m\ge0}\sum_{n\ge0}[z]\exp(-(1+q)\log(\pi i(m+n+1)-z))\\&=2\pi i\cdot\sum_{n\ge1}n\cdot[z]\exp(-(1+q)\log(\pi in-z))\\&=-i\cdot\frac{2(1+q)}{\pi^{1+q}}\exp(-\pi iq/2)\cdot\zeta(q+1)\end{align}$$
And now we get to finish!
$$\begin{align}J&=-\frac{\pi q(q-1)}{\beta^{3-q}}\csc(\pi q)\cdot\Re(K)\\&=-\frac{\pi q(q-1)}{\beta^{3-q}}\csc(\pi q)\cdot\Re\left(-i\cdot\frac{2(1+q)}{\pi^{1+q}}\exp(-\pi iq/2)\cdot\zeta(q+1)\right)\\&=\frac{q(q^2-1)}{\pi^q\cdot\beta^{3-q}}\cdot\sec\frac{\pi q}{2}\cdot\zeta(q+1)\end{align}$$
It is with even greater pleasure that I can say this closed form matches the numerics.
