Is there an analogous physical interpretation for the third and higher mass moments?
While not quite an inherently "physical interpretation," the technique of inertial imaging allows one to use higher order mass moments for characterization and identification of biological molecules and molecular complexes. (See: Inertial imaging with nanomechanical systems)
The basic idea here is that sticking ("adsorbing") a molecule ("analyte") onto a nanoscale mechanical resonator ("NEMS" = nanoelectromechanical system) will shift the NEMS resonant frequencies in a way that allows us to determine the mass moments of the molecule. Different molecules (even ones with similar masses) will generally differ in their spatial size (2nd moment) and asymmetries (3rd and higher moments), so by measuring these higher moments one can improve identification of unknown analytes (by comparing to some catalogue of moments of known analytes).
More formally, consider a 1D resonant system like a NEMS doubly-clamped beam or cantilever. (We work in 1D here for simplicity, this all extends very naturally to higher dimensions.) This device will have modes $\phi_n(x)$ with resonant frequencies $\omega_n^{(0)}$. For convenience we normalize these modes so that $\int \rho_{dev}(x) |\phi_n(x)|^2 dx = M$ with $M$ the mass of the device and $\rho_{dev}(x)$ its mass density.
Now suppose some small, soft, compliant analyte adsorbs onto the device. Under these assumptions, we can consider the effect of the analyte to be some small mass perturbation $\rho(x)$ on the system. This 1D mass distribution $\rho$ is what we are going to be measuring the moments of. For things like proteins or complexes (e.g. hemoglobin) this $\rho$ is a good measure of some linear projection of the original 3D mass distribution, and carries more information about the analyte than simply measuring its total mass $\int \rho(x) dx$. This additional information can be used to characterize known proteins and help identify unknown ones.
For small masses, the effect of this perturbation will be to lower the resonant frequencies $\omega_n^{(0)}$ to some new frequencies $\omega_n$ while preserving the original vibrational modes $\phi_n$. From this, one can show that the shift in the resonant frequency $\Delta_n \equiv (\omega_n-\omega_n^{(0)})/\omega_n$ is given by $\Delta_n = -\frac{1}{2M} \int \rho(x) |\phi_n(x)|^2 dx$. These frequency shifts can be measured through a number techniques to varying degrees of precision.
The key to extracting the moments is to consider linear combinations of these frequencies. Effectively, we will be projecting the moment-generating functions $g^{(k)}(x) \sim x^k$ onto the basis of vibrational modes $\phi_n(x)$. Let $F_n \equiv \int \rho(x)|\phi_n(x)|^2 dx$ (so $\Delta_n = -\frac{1}{2M}F_n$). Then since integration is linear, a linear combination of these $F_n$ factors nicely as:
$\sum_n \alpha_n F_n = \int \rho(x) [\sum_n \alpha_n |\phi_n(r)|^2] dx \equiv \int \rho(x) g(x) dx = -\frac{1}{2M}\sum_n \alpha_n \Delta_n$.
So if we solve for a particular set $\alpha_n^{(k)}$ such that $g(x) = \sum_n \alpha_n |\phi_n(x)|^2 \sim x^k$ locally about the position of the analyte $x_0$, then we can use the frequency shifts $\Delta_n$ to determine the mass moments $\int \rho(x) x^k dx$. Given a desired region of convergence and known mode shapes $\phi_n$, it is relatively easy to compute best-fit coefficients $\alpha_n^{(k)}$. With enough (see below) modes and analytes that are small compared to the measuring device, the error from this approximation can be reduced to ~1%.
In one dimension, one generally needs to measure $k+1$ frequencies in order to extract the $k$th moment. This is sufficient for infinitesimal sized analytes; for those whose spatial size is not infinitesimal compared to the measuring device, the accuracy increases quickly as more frequencies are measured.
In the paper linked at the top of this answer, the authors use these methods to measure the mass, size (variance), and asymmetry (skewness) of some liquid droplets that agree well with optical measurements.
In statistics the third moment is used to calculate skewness. I would guess this has a physical analogy. Although I haven't thought this through, I'd guess it would be possible to take a disk and deform it asymmetrically so that the centre of mass and moment of inertia remained the same, but the third moment changed. In this case the third moment would be telling you about the asymmetry of the disk.
