There are a number of ways to interpret the spectral form factor (SFF), but they pretty much all follow from random matrix theory (RMT) moreso than physical first principles. The primary interpretation is that the SFF is the (Fourier transform of the) two-point correlation function of eigenvalues of the Hamiltonian, with
$$K(t) \, = \, \overline{{\rm tr} \left[\, e^{- {\rm i} \, t \, H}\, \right] {\rm tr} \left[\, e^{\, {\rm i} \, t \, H} \, \right] } \, = \, \int d\omega \, e^{\, {\rm i} \, \omega \, t} \, \rho^{\,}_2 (E,E+\omega),$$
where the overline denotes the average over an ensemble of statistically similar systems (the SFF is not so well behaved for a single realization of H, and is smoothed out considerably by introducing a small amount of disorder). The second formula appears in a reference linked in another answer, and $\rho^{\,}_2$ is the two-point density of eigenvalues of $H$ at energies $E$ and $E+\omega$. More on what this means in a bit.
The reason it is mentioned so frequently is that we basically equate "being thermal" with having the spectral properties of a random matrix ensemble (for a given Hamiltonian or time-evolution operator). The idea is that thermal / chaotic systems forget as much information as possible, except what they must remember according to symmetries; their energy levels / eigenvalues should repel one another; and microscopic details should be relatively unimportant, on average.
The spectral form factor is used as a diagnostic of RMT level statistics. Random matrices and thermal systems both show spectral rigidity in the form of level repulsion, meaning that eigenvalues of $H$ like to space themselves out (this is also true, e.g., for Hamiltonians $H$ drawn from Gaussian ensembles).
The fingerprint of this level repulsion is the scaling $K(t) = t$, which holds if $H$ is drawn from a Gaussian ensemble, for actual chaotic models, and for evolution under random unitaries (i.e. using a Floquet operator $U(T) = e^{- {\rm i} \, T \, H}$ with $U(nT) = U^n$). The latter case appears in some papers from 2018-2022 by John Chalker et al on quantum chaos and spectral statistics using random unitary circuits (e.g., this one). There are a few papers by Brian Swingle et al that could be helpful, and another by Steve Shenker and company with $\mathsf{U}(1)$ symmetry. Also, in the random unitary circuit context (but not generally), Appendix D1 of this paper derives that the spectral form factor is also equal to the average over all two point functions:
$$K(t) \, = \, \sum\limits_{a=1}^{\mathcal{D}^2} \, \left\langle \mathcal{O}^{\dagger}_a (t) \mathcal{O}^{\vphantom{\dagger}}_a (0) \right\rangle \, , ~$$
where $\mathcal{D}$ is the many-body Hilbert space dimension, the label $a$ runs over orthonormal many-body basis operators (e.g., the Pauli strings), and the average is computed at infinite temperature (since the model is Floquet, this is the only thermal state). This relation between the SFF (spectral properties) was also conjectured by Shenker, but I'm not sure that it generalizes.
The SFF can also be used to see the onset of thermalization, as in this paper. Basically, if we look at the work by Chalker et al, the expectation is that $K(t) = t$ if the Floquet operator is all-to-all, or if it acts locally on alternating bonds, e.g. However, if we only apply single-site unitaries, then instead we have $K(t) = t^N$, where $N$ is the number of sites in the lattice.
To get some intuition into thermal systems, imagine that we start from some unentangled (product) initial state and then evolve in time under some chaotic model. At early times, the interactions haven't had a chance to "kick in" so the system behaves like a free system, with $K(t) \sim t^N$. We then imagine that the interactions start to merge sites into blocks, and suppose that at time $t$ the blocks have size $\xi (t)$. Correspondingly, we have $K(t) \sim t^{N / \xi (t)}$.
We then suppose that $\xi(1)=1$ (so that in the first time step, the "thermal blocks" have size one), and then at some later "thermalization time" $\tau^{\,}_{\rm th}$, we have $\xi (\tau^{\,}_{\rm th}) = N$, meaning that the block has grown to include the full system. Since the blocks cannot continue to grow, we have $\xi (t) = N$ and therefore $K(t)=t$ for $t \geq \tau^{\,}_{\rm th}$, corresponding to the linear ramp, which is a fingerprint of chaos. In the first Chalker paper and a few others (maybe this one) there is mention of some Coulombic interaction between the energy levels that leads to this ramp, but I don't really understand this point to be frank.
Finally, for times $t > \tau^{\,}_{\rm Heis} = \mathcal{D}$, we have $K(t) = \mathcal{D}$. The time $\tau^{\,}_{\rm Heis} = \mathcal{D}$ is the "Heisenberg time", which is equal to the mean inverse level spacing $1 / \delta E \sim {\rm tr} \left[ \, \mathbb{1} \, \right]$, which is equal to the many-body Hilbert space dimension. Essentially, when $t> \mathcal{D}$, there is nothing left to probe...for $t< \mathcal{D}$, the SFF probes relations between energy levels with separation of $\mathsf{O}(1/t)$, but for $t> \mathcal{D}$, all separations have been accounted for. This is quite hand wavy, and part of the reason is that this sharp feature where $K(t) = t$ for $t<\mathcal{D}$ and $K(t) = \mathcal{D}$ for $t \geq \mathcal{D}$ does not seem to be something one can derive from RMT (at least according to the Chalker papers).
So anyway, hopefully this gives some insight into what the SFF $K(t)$ is meant to capture, why it's expected to work / the motivation for using it, and how it's derived and used. Mostly, it comes from the recognition that random matrices do a very good job of approximating thermal systems, plus the fact that the SFF has a very nice form for random matrices that is easy to benchmark against.
