Calculate $\frac{1}{H^2W^2}\int_0^H{\int_0^W{\int_0^H{\int_0^W{e^{-\alpha \sqrt{\left( x-z \right) ^2+\left( y-t \right) ^2}}}}dxdydzdt}}$ I am stucked by the following problem in part of my current research process. Could you please help me?
(You can first jump to the Background part to see what I am calculating for.)

Problem:
Given three constant real positive numbers $H,W,\alpha \in \mathbb{R}^+$, calculate the following quadruple integral:
$$
I(H,W,\alpha)=\frac{1}{H^2W^2}\int_0^H{\int_0^W{\int_0^H{\int_0^W{e^{-\alpha \sqrt{\left( x-z \right) ^2+\left( y-t \right) ^2}}}}dxdydzdt}}.
$$

Some of My Efforts:
The original problem is too difficult for me.
I just want the accurate answer and do not have to require the complete analysing process.
So I first tried to get the answer from WolframAlpha. But it seems to be failed to calculate the integral.

Background:
I was trying to calculate the average "correlation" value of a two-dimensional area $[0,H]\times[0,W]$. In my task, the "correlation" $\rho$ is modeled as:
$$
\rho(d)=e^{-\alpha d},
$$
where $d\in \mathbb{R}^+$ is the distance of two points.
For a given point $(z,t)$ satisfying $0\le z\le H$ and $0\le t \le W$, its average "correlation" is:
$$
f(z,t)=\frac{1}{HW}\int_0^H{\int_0^W{\rho \left( \sqrt{\left( x-z \right) ^2+\left( y-t \right) ^2} \right)}dxdy}.
$$
And for the whole area, the average "correlation" is:
$$
I=\frac{1}{HW}\int_0^H{\int_0^W{f\left( z,t \right)}dzdt}.
$$
Finally I get the formula in the original problem, resulting in a quadruple integral.

Some Other Aspects of Help That You May Provide:
I think that the original problem is too difficult. There may be some other operations helpful for me:

*

*I can change the modeling of $\rho$. I model it from my observed data (see the figure below). Since when $d=0$, there is $\rho=1$ (the value of each point is strongly correlated to itself), and $\rho$ is approximately decreasing but always positive in $[0,+\infty)$. I guess the formula could be $\rho(d)=e^{-\alpha d}$. Is there a better modeling?

(The horizontal axis corresponds to $d$, and the vertical axis corresponds to $\rho$. The values of $\rho$ are more reliable when $0\le d\le 300$, since the data sample number of $0\le d\le 300$ are sufficiently large and number of $d>300$ are relatively small.)



*I can change the calculation way of "average". Since the integration may be too difficult to calculate, I can use another way to measure the "average" of the area. Is there a better way? I try not to modify the original integral form, since it is easy to interpret and understand.


*Instead of directly calculating $I(H,W)$, I can also change to calculate $I(H_1,W_1)-I(H_2,W_2)$, since I care the difference of "average correlations" of two areas. However, the difference seems to be more difficult to calculate.
Since $\rho$ is decreasing in $[0,+\infty)$, I intuitively think that a small area will own a larger "average correlation", i.e., $I(H_1,W_1)<I(H_2,W_2)$ when $H_1\times W_1 > H_2\times W_2$.
As you can see, I just have a vague and qualitative cognition of my above measuring and modeling now. So many things can be adjusted to make the problem easier. (I am now looking for them.)

Could you please help me based on all the above descriptions?
(I would like to add more details if the statements are not clear enough.)
 A: Too Long For Comment:
Lemma: For any integrable function $f$, we have
$$\begin{align}
\int_{0}^{1}\mathrm{d}x\int_{0}^{1}\mathrm{d}y\,f{\left((x-y)^{2}\right)}
&=2\int_{0}^{1}\mathrm{d}x\int_{0}^{x}\mathrm{d}y\,f{\left((x-y)^{2}\right)};~~~\small{symmetry}\\
&=2\int_{0}^{1}\mathrm{d}x\int_{0}^{x}\mathrm{d}u\,f{\left(u^{2}\right)};~~~\small{\left[y=x-u\right]}\\
&=2\int_{0}^{1}\mathrm{d}u\int_{u}^{1}\mathrm{d}x\,f{\left(u^{2}\right)}\\
&=2\int_{0}^{1}\mathrm{d}u\,\left(1-u\right)f{\left(u^{2}\right)}.\\
\end{align}$$

Suppose $\left(H,W,\alpha\right)\in\mathbb{R}_{>0}^{3}$. Applying the above lemma to $\mathcal{I}$, we find
$$\begin{align}
\mathcal{I}{\left(H,W,\alpha\right)}
&=\frac{1}{H^{2}W^{2}}\int_{0}^{H}\mathrm{d}x\int_{0}^{W}\mathrm{d}y\int_{0}^{H}\mathrm{d}z\int_{0}^{W}\mathrm{d}t\,e^{-\alpha\sqrt{(x-z)^{2}+(y-t)^{2}}}\\
&=\int_{0}^{1}\mathrm{d}x\int_{0}^{1}\mathrm{d}y\int_{0}^{1}\mathrm{d}z\int_{0}^{1}\mathrm{d}t\,e^{-\alpha\sqrt{H^{2}(x-z)^{2}+W^{2}(y-t)^{2}}};~~~\small{(x,y,z,t)\mapsto(Hx,Wy,Hz,Wt)}\\
&=\int_{0}^{1}\mathrm{d}x\int_{0}^{1}\mathrm{d}z\int_{0}^{1}\mathrm{d}y\int_{0}^{1}\mathrm{d}t\,e^{-\alpha\sqrt{H^{2}(x-z)^{2}+W^{2}(y-t)^{2}}}\\
&=4\int_{0}^{1}\mathrm{d}u\int_{0}^{1}\mathrm{d}v\,\left(1-u\right)\left(1-v\right)e^{-\alpha\sqrt{H^{2}u^{2}+W^{2}v^{2}}}.\\
\end{align}$$
We can further reduce the integral  as follows:
$$\begin{align}
\mathcal{I}{\left(H,W,\alpha\right)}
&=4\int_{0}^{1}\mathrm{d}u\int_{0}^{1}\mathrm{d}v\,\left(1-u\right)\left(1-v\right)e^{-\alpha\sqrt{H^{2}u^{2}+W^{2}v^{2}}}\\
&=4\int_{0}^{1}\mathrm{d}u\int_{0}^{u}\mathrm{d}v\,\left(1-u\right)\left(1-v\right)e^{-\alpha\sqrt{H^{2}u^{2}+W^{2}v^{2}}}\\
&~~~~~+4\int_{0}^{1}\mathrm{d}u\int_{u}^{1}\mathrm{d}v\,\left(1-u\right)\left(1-v\right)e^{-\alpha\sqrt{H^{2}u^{2}+W^{2}v^{2}}}\\
&=4\int_{0}^{1}\mathrm{d}u\int_{0}^{u}\mathrm{d}v\,\left(1-u\right)\left(1-v\right)e^{-\alpha\sqrt{H^{2}u^{2}+W^{2}v^{2}}}\\
&~~~~~+4\int_{0}^{1}\mathrm{d}v\int_{v}^{1}\mathrm{d}u\,\left(1-v\right)\left(1-u\right)e^{-\alpha\sqrt{H^{2}v^{2}+W^{2}u^{2}}};~~~\small{(u,v)\mapsto(v,u)}\\
&=4\int_{0}^{1}\mathrm{d}u\int_{0}^{u}\mathrm{d}v\,\left(1-u\right)\left(1-v\right)e^{-\alpha\sqrt{H^{2}u^{2}+W^{2}v^{2}}}\\
&~~~~~+4\int_{0}^{1}\mathrm{d}u\int_{0}^{u}\mathrm{d}v\,\left(1-v\right)\left(1-u\right)e^{-\alpha\sqrt{H^{2}v^{2}+W^{2}u^{2}}}\\
&=4\int_{0}^{1}\mathrm{d}u\int_{0}^{u}\mathrm{d}v\,\left(1-u\right)\left(1-v\right)\left[e^{-\alpha\sqrt{H^{2}u^{2}+W^{2}v^{2}}}+e^{-\alpha\sqrt{W^{2}u^{2}+H^{2}v^{2}}}\right]\\
&=4\int_{0}^{1}\mathrm{d}u\int_{0}^{1}\mathrm{d}x\,u\left(1-u\right)\left(1-xu\right)\left[e^{-\alpha\sqrt{H^{2}u^{2}+W^{2}x^{2}u^{2}}}+e^{-\alpha\sqrt{W^{2}u^{2}+H^{2}x^{2}u^{2}}}\right];~~~\small{\left[v=xu\right]}\\
&=4\int_{0}^{1}\mathrm{d}x\int_{0}^{1}\mathrm{d}u\,u\left(1-u\right)\left(1-xu\right)\left[e^{-u\alpha\sqrt{H^{2}+W^{2}x^{2}}}+e^{-u\alpha\sqrt{W^{2}+H^{2}x^{2}}}\right].\\
\end{align}$$
The integral over $u$ in the last line above is elementary, so in principle we have successfully reduced $\mathcal{I}$ to a single variable integral. Unfortunately, this integral for $u$ is quite cumbersome and doesn't lead to an integral over $x$ that is any easier, so I won't bother reproducing it here. It could help with numerical procedures though.

A: You can make the computation significantly easier if you change the correlation strength to
$$\rho(d)=e^{-ad^2}$$
which shouldn't affect the results qualitatively. This choice of function comes with the added advantage that the integrals over the coordinates factorize:
$$I=\frac{1}{H^2}\left(\int_0^H dx\int_0^H dz 
e~^{-a(x-z)^2}\right)\frac{1}{W^2}\left(\int_0^W dy\int_0^W dt 
e~^{-a(y-t)^2}\right):=F(W)F(H)$$
which according to Mathematica, evaluates to:
$$F(H)=\sqrt{\frac{\pi}{aH^2}}\text{erf}(H\sqrt{a})+\frac{e^{-aH^2}-1}{aH^2}$$
