I think it can be found the following way. If $R$ is a commutative local ring then every idempotent matrix is conjugate to a diagonal idempotent $\begin{pmatrix} I_r & 0\cr 0 &0\end{pmatrix}$ . The point is that projective modules are free for local rings. This lets you write $R^n$ as a direct sum of the image and the kernel, both of which are free.
So now one has to just compute the size of the stabilizer of the standard rank $r$ diagonal idempotent under the conjugation action of $\mathrm{GL}(R)$ for $R=\mathbb{Z}_{p^{k}}$.
Added. I believe the stabilizer of the rank $r$ idempotent for a local ring is $\mathrm{GL}(R,r)\times\mathrm{GL}(R,n-r)$ like in the field case and so a formula is easily found.
Added. I compute the answer for $n\times n$ matrices over $\mathbb{Z}_{p^k}$ to be $$\sum_{r=0}^n \frac{p^{2r(n-r)(k-1)}|\mathrm{GL}(n,p)|}{|\mathrm{GL}(r,p)|\cdot|\mathrm{GL}(n-r,p)|}.$$ Here I use a matrix is invertible over a local ring iff it is over the residue field. The orders of the $\mathrm{GL}(m,p)$ are of course well known.