Note that $x\mapsto \frac{x}{|x|^2}$ is an involution. So for any $m$ the map on functions
$$ u(x) \mapsto \frac{1}{|x|^m} u(\frac{x}{|x|^2})$$
is an involution ( see @Ivo Terek comment above). However, the Kelvin $\mathcal{K}$ has the extra property that it takes a constant function to a multiple of $\frac{1}{|x|^{n-2}}$, a harmonic function.
Note also that $\mathcal{K}$ commutes with orthogonal changes of coordinates, and behaves nicely with respect to central homotheties.
Let's see what $\mathcal{K}$ does to some known harmonic functions, say the potential of a point $u(x) = \frac{1}{|x-x_0|^{n-2}}$. When $x_0=0$ we get a constant ( $\mathcal{K}$ involution). From the above, it is enough to consider the case $x_0$ on the unit sphere, (in fact enough for $x_0 = (1,0, \ldots, 0)$). A simple calculation shows that $\mathcal K u = u$ ( one should check that). So $\mathcal{K}$ takes the potential of a point on the unit sphere to itself, and the potential of any other point to a multiple of the potential of the inverse of that point). Now the point potentials "linearly generate" the space of harmonic functions. We are done.
Note: this is more of a plea than a proof. Perhaps Kelvin realised how the transformation should work by looking at potentials, and then gave a formal proof. Calculations show that we have
$$\Delta \mathcal{K} = \frac{1}{|x|^4} \mathcal{K} \Delta $$
$\bf{Added:}$ The calculations for the point potential are simple. Let's start with an inversion map
$\mathcal{I}$ of center $O$ and modulus $\rho$. Consider two points $x$, $y$, and their transforms $x'$, $y'$ under $\mathcal{I}$. The triangles $Oxy$, and $Oy'x'$ are similar, so we have
$$\frac{x'y'}{yx} = \frac{Ox'}{Oy} = \frac{Ox'\cdot Ox}{Ox\cdot Oy}$$ that is
$$|x'-y'| = \frac{\rho}{|x|\cdot |y|} \cdot |x-y|$$
(for us $\rho =1$). This is the only calculation we need to see how $\mathcal K$ acts on point potentials.