Let me try a very explicit proof. Consider the problem of taking a solution modulo $p^k$ and lifting to solutions modulo $p^{k+1}$. For this, we may assume that $X \subset \mathbb{A}^m$ is smooth and affine of dimension $n$, given by polynomials $f_1, \dotsc, f_r \in \mathbb{Z}_p[x_1, \dotsc, x_m]$. Let $\mathbf{a} \in \mathbb{Z}_p^m$ satisfy $f_1(\mathbf{a}) \equiv \dotsb \equiv f_r(\mathbf{a}) \equiv 0 \pmod{p^k}$.

To find points modulo $p^{k+1}$ that lift $\mathbf{a}$, you put $\mathbf{x} = \mathbf{a} + p^k \mathbf{y}$ and look for vectors $\mathbf{y}$ modulo $p$ giving solutions to $f_i(\mathbf{x}) \equiv 0 \pmod{p^{k+1}}$. Taylor expansion gives

$(f_1(\mathbf{x}), \dotsc, f_r(\mathbf{x})) = (f_1(\mathbf{a}, \dotsc, f_r(\mathbf{a})) + p^k \mathbf{J}(\mathbf{a}) \mathbf{y} + O(p^{k+1})$

where $\mathbf{J}$ is the Jacobian matrix $(\partial(f_i)/\partial x_j)$. Dividing by $p^k$ and reducing modulo $p$ gives an inhomogeneous linear equation over $\mathbf{F}_p$:

$\mathbf{J}(\mathbf{a}) \mathbf{y} \equiv -p^{-k} (f_1(\mathbf{a}, \dotsc, f_r(\mathbf{a})) \pmod{p}$ .

That $X$ is smooth over $\mathbb{Z}_p$ at $\mathbf{a}$ implies that $\mathbf{J}(\mathbf{a})$, when reduced modulo $p$, has rank $m-n$. The space of solutions $\mathbf{y}$ is non-empty, by Hensel's Lemma, so is a linear space of dimension $n$ over $\mathbb{F}_p$. In this way you see that every solution modulo $p^k$ lifts to precisely $p^n$ solutions modulo $p^{k+1}$, and you get your formula by induction.