$$\begin{align}
&= \end{align}$$
Let's look at the diagonal of your PSF, i.e. the points where $x=y$:
$$\begin{align}
\cosh^{-2}\left(\sqrt{x^2+y^2}/w\right)&\overset{x=y}= \exp\left(-2{\sqrt{2\frac{x^2}{w^2}}}\right) + \exp\left({2\sqrt{2\frac{x^2}{w^2}}}\right)\\
&=\exp\left(2\sqrt2 \frac xw\right)+ \exp\left(-2\sqrt2 \frac xw\right)\\
&=\cosh\left(\frac{2\sqrt2}w x\right)\\&=: d_P(x)
\end{align}$$
Let us thus take the $x=y$ diagonal of your convolved image $d_g=g(x,x)$, and notice that it's the convolution of the diagonal of the original image $f(x,y)$ with $d_P(x)$.
$$\begin{align}
d_g(x) &= d_f(x) * d_P(x)\end{align}$$
We have thus reduced your two-dimensional deconvolution problem to a one-dimensional one. Nice!
Next, let's get rid of actual image content. I'm from a digital communications background, so I'd actually try to convolve with different $\cosh$ scalings, and see which one of these scalings delivers me the highest "contrast".