Consider the time-frequency $(t, f)$ plane interpretation. Suppose we start with $f(t)$ along time axis. Define the rotation matrix (called Elliptic group in this context) such that
$$
M_\gamma = \begin{pmatrix}\cos(\gamma) & -\sin(\gamma)\\ \sin(\gamma) & \cos(\gamma)\end{pmatrix}, \qquad \qquad \gamma = m \pi/2
$$
Thus, applying Fourier once rotates it toward frequency axes, apply it again, moves it in the opposite side of time axes (As you can verify through integration).
In general, to apply the Fourier transform $m$ times is to rotate the signal with $\gamma = \frac{m\pi}{2}$ degree (with time and frequency axes being orthogonal).
$$\mathcal{F}^m = \begin{cases}f(t)& \mbox{ if } m \equiv 0 \mod 4 \\\mathcal{F} & \mbox{ if } m \equiv 1 \mod 4\\ f(-t) & \mbox{ if } m \equiv 2 \mod 4\\ \mathcal{F}^{-1} & \mbox{ if } m \equiv 3 \mod 4\end{cases}
$$
Hence, it is not self-inverse, but rather enjoys a 4-fold periodicity property.
It should also be noted that rotation by arbitrary degrees, (i.e. $m$ not integer) is at the heart of the Fractional Fourier transform.
EDIT ($M_\gamma$ and integral transform)
Briefly, multiplication of $f(t)$ by $t$ in time domain (let that be through an operator, $Q$) corresponds to differentiation (with rotation by -i) in frequency domain (through an operator $P$). Define the (linear) operator $C_M$ such that $C_M Q = CQC^-1 = dQ - bP$ and $CPC^{-1} = -cQ + aP$ parametrized by
$$M = \begin{bmatrix}a & b\\ c& d\end{bmatrix} \in SL(2, R)$$
When $M$ is a rotation matrix, i.e. $M = M_\gamma$, and specifically $\gamma = \frac{\pi}{2}$ (so that $a, d = \cos(\pi/2) = 0, b = -c = 1$), we get the Fourier transform (applied once). For example, its action on $tf(t)$ is $C_M tf(t) = CQC^{-1}F(f) = -P F(f) = i\frac{d F}{df}$.
More importantly, composition of $C_{M_1}$ and $C_{M_2}$ corresponds to multiplication of $M_1$ and $M_2$, which explains the four-periodicity. The argument for $Q$ (with similar one for $P$) is as follows:
\begin{align}C_{M_2}Q'C_{M_2}^{-1} &= C_{M_2}C_{M_1}Q(C_{M_2}C_{M_1})^{-1} \\&= C_{M_2}(d_1Q - b_1P)C_{M_2}^{-1} \\&= c_2b_1 + d_2d_1Q - (a_2b_1 + b_2d_1)P \\&= C_{M_1M_2}QC_{M_1M_2}^{-1}\end{align}
In terms of kernel of the integral transform, $\int df F(f) C_M(t,f)$ is defined to be $C_M(t, f) = \theta_M e^{i(at^2 - 2tf + d f^2)/2b}$, which is $e^{-itf}$ for Fourier transform case, up to a constant $\theta_M$.
All the interesting properties of integral transforms are encoded through the matrix $M$, see chapter 9 in Wolf's Integral Transforms in Science and Engineering.
I'm sure the Quantum mechanics interpretation/motivation will be valuable, especially for how these came about in the first place, as hinted by the post of @Spencer here.