Via Fourier Integral:
$\def\Shi{\operatorname{Shi}}$ The hyperbolic sine integral Shi$(x)$ is invertible via complex fourier series converging on $[-L,L]=[-\Shi(l),\Shi(l)]$:
$$f(x)=\Shi(x)\implies f^{-1}(x)=\frac1{2L}\sum_{n\in\Bbb Z}e^\frac{\pi i nx}L\int_{-f^{-1}(L)}^{f^{-1}(L)} e^{-\frac{i\pi n\Shi(w)}L}\sinh(w)dw=\frac1{2\Shi(l)}\sum_{n\in\Bbb Z}e^\frac{\pi i nx}{\Shi(l)}\int_{-l}^l e^{-\frac{i\pi n\Shi(w)}{\Shi(l)}}\sinh(w)dw$$
Plotting $f^{-1}(\Shi(x))=x$ by truncating the Fourier series shows the inverse function:
$\lim\limits_{L\to\infty}$ should use the Fourier integral theorem:
$$f^{-1}(x)\mathop=^?\iint_{\Bbb R^2}e^{2\pi i(x-w)t}f^{-1}(w)dwdt=2\int_0^\infty\int_0^\infty e^{2\pi i(x-\Shi(w))t}\sinh(w)dwdt $$
Unfortunately, it is hard to compute and the integrals are not interchangeable. Trying $y=\frac wl$ and assuming one may interchange the limit:
$$f(x)\mathop=^?\sum_{n\in\Bbb Z}\int_{-1}^1\lim_{l\to\infty}e^\frac{\pi i n(x-\Shi(ly))}{\Shi(l)}\frac{l\sinh(ly)}{2\Shi(l)} dy$$
Via Sifting Property:
Using the Dirac delta $\delta(x)$ sifting property gives:
$$f^{-1}(x)=\int_a^b\delta(\Shi(w)-x)\sinh(w)dw,a\le f^{-1}(x)\le b $$
After a $\delta(y)$ integral representation:
$$f^{-1}(x)\mathop=^? \frac1{2\pi} \int_{-\infty}^\infty\int_a^b e^{it(\Shi(w)-x)}dt\sinh(w)dw\mathop=^?\frac1{2\pi}\int_a^b\int_{-\infty}^\infty e^{it(\Shi(w)-x)}\sinh(w)dwdt$$
Numerical tests are hard to process even though the integral does converge. Strong evidence for this method working is in this answer which uses the sifting property to derive an inverse Fourier transform
If both of these methods work, what is a final double integral result?