As Keith Conrad says, use the pigeonhole principle to find $a$ and $b$ satisfying $a^2+b^2+1 \equiv 0 \bmod p$. Set
$$I = \begin{bmatrix} 0&-1 \\ 1&0 \end{bmatrix} \qquad J = \begin{bmatrix} a&b \\ b&-a \end{bmatrix} \qquad K = \begin{bmatrix} -b&a \\ a&b \end{bmatrix}.$$
Then $I$, $J$ and $K$ obey the quanternion relations. So we have a map from $\mathbb{H}/p \mathbb{H}$ to $\text{Mat}_{2 \times 2}(\mathbb{F}_p)$ sending $i \mapsto I$, $j \mapsto J$ and $k \mapsto K$.
Assuming that $p \neq 2$, the four matrices $\text{Id}_2$, $I$, $J$ and $K$ are a basis for the vector space of $2 \times 2$ matrices, so this map is an isomorphism.
Remark: I assumed that your $\mathbb{Z}_p$ meant $\mathbb{Z}/p \mathbb{Z}$. If you meant the $p$-adics, this solution also works: You can find a solution to $a^2+b^2+1=0$ modulo $p$ by the pigeonhole principle and then lift it to the $p$-adics by Hensel's lemma.
Remark: When $p \equiv 1 \bmod 4$, you can take $a$ to be a square root of $-1$ and $b=0$, and thus recover the standard the standard description of the quaternions in terms of Pauli matrices.
$H \cong \operatorname{Mat}_2(\mathbb Z_p)$
, not H $\cong \operatorname{Mat}_2(\mathbb Z_p)$H $\cong \operatorname{Mat}_2(\mathbb Z_p)$
. I have edited accordingly. $\endgroup$