Why is the momentum operator imaginary?
The simplest explanation hinges on the fact that observables are represented by Hermitian operators in quantum mechanics. Once we accept this, then we can show that the momentum operator $\hat{p} = -\mathrm{i}\hbar\nabla$ is Hermitian precisely because of the factor of $\mathrm{i}$. We need to show that for any two states $f$ and $g$, $\langle f | \hat{p} | g\rangle = \langle g | \hat{p} | f\rangle^*$. For simplicity we work in one dimension, such that $\hat{p} = -\mathrm{i}\hbar(\mathrm{d}/\mathrm{d}x)$. The proof is taken from Binney and Skinner's The Physics of Quantum Mechanics (p 36), although truthfully it is quite straightforward, the main idea being integration by parts:
$$\begin{align}
\langle f | \hat{p} | g\rangle &= -\mathrm{i}\hbar \int_{-\infty}^{\infty} f^* \frac{\mathrm{d}g}{\mathrm{d}x} \,\mathrm{d}x \\[8pt]
&= -\mathrm{i}\hbar \left\{ \left[f^*g\right]_{-\infty}^{\infty} - \int_{-\infty}^{\infty} g \frac{\mathrm{d}f^*}{\mathrm{d}x} \,\mathrm{d}x \right\} \\[8pt]
&= \mathrm{i}\hbar \int_{-\infty}^{\infty} g \frac{\mathrm{d}f^*}{\mathrm{d}x} \,\mathrm{d}x \\[8pt]
&= \left( -\mathrm{i}\hbar \int_{-\infty}^{\infty} g^* \frac{\mathrm{d}f}{\mathrm{d}x} \,\mathrm{d}x \right)^* \\[8pt]
&= \langle g | \hat{p} | f \rangle ^*
\end{align}$$
(Well-behaved functions vanish as $x \to \pm \infty$, hence the term $[f^*g]_{-\infty}^{\infty}$ goes to zero.) In a similar way, it can be shown that without the the operator $\mathrm{d}/\mathrm{d}x$ itself is not Hermitian, and thus cannot be an operator corresponding to an observable quantity in QM.
The momentum of the electron can be found using the equation, $\hat p \Psi = \bar p \Psi$
This is an eigenvalue equation for the momentum, but it does not actually give you the value of the momentum. When you apply an operator to a state, you get out another state, and there's nothing that forbids this second state from being complex / imaginary. It's only the final expectation value at the end that has to be real, and this expectation value is given by $\langle \psi | \hat{p} | \psi \rangle$.
Is the momentum really imaginary?
It turns out that the factor of $\mathrm{i}$ doesn't matter, because if you calculate this correctly, then the momentum turns out to be zero. This makes sense, because all directions in the hydrogen atom are equal, and it can't be that the electron has a net momentum in some direction.
Calculating it correctly, however, is a slightly tricker topic than expected. In this case, the momentum operator $\hat{p}$ is a vector operator ($\equiv (\hat{p}_x, \hat{p}_y, \hat{p}_z)$), and consequently when you integrate it you need to be very careful to not mix up vector and scalar quantities. If you wish to do this in spherical coordinates, I first suggest going to Physics Stack Exchange and searching for questions about "momentum of hydrogen atom". You will find lots of similar questions asking about the factor of $\mathrm{i}$ and how it doesn't make sense.
Personally, I think a more direct, and slightly easier approach, is to simply integrate one (Cartesian) component of the vector at a time. Thus, we can try to evaluate $\langle \psi | \hat{p}_x | \psi \rangle$. For simplicity, I will drop the normalisation constant and use $\psi = \exp{(-r/a_0)}$. Also, note that
$$\begin{align}
\frac{\mathrm{d}r}{\mathrm{d}x} &= \frac{\mathrm{d}}{\mathrm{d}x} (x^2 + y^2 + z^2)^{1/2} \\
&= \frac{1}{2}(x^2 + y^2 + z^2)^{-1/2} \cdot (2x) \\
&= \frac{x}{r}
\end{align}$$
Thus, we have that
$$\begin{align}
\langle \psi | \hat{p}_x | \psi \rangle &= -\mathrm{i}\hbar \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \exp\left(-\frac{r}{a_0}\right) \frac{\mathrm{d}}{\mathrm{d}x}\left[\exp\left(-\frac{r}{a_0}\right)\right] \,\mathrm{d}x\,\mathrm{d}y\,\mathrm{d}z \\[8pt]
&= -\mathrm{i}\hbar \iiint \exp\left(-\frac{r}{a_0}\right) \frac{\mathrm{d}}{\mathrm{d}r}\left[\exp\left(-\frac{r}{a_0}\right)\right]\frac{\mathrm{d}r}{\mathrm{d}x} \,\mathrm{d}x\,\mathrm{d}y\,\mathrm{d}z \\[8pt]
&= -\mathrm{i}\hbar \iiint \exp\left(-\frac{r}{a_0}\right) \left(-\frac{1}{a_0}\right)\exp\left(-\frac{r}{a_0}\right)\frac{x}{r} \,\mathrm{d}x\,\mathrm{d}y\,\mathrm{d}z \\[8pt]
&= \frac{\mathrm{i}\hbar}{a_0} \iiint \exp\left(-\frac{2r}{a_0}\right) \frac{x}{r} \,\mathrm{d}x\,\mathrm{d}y\,\mathrm{d}z \\[8pt]
\end{align}$$
At this point, it's easiest to bring in symmetry arguments. Notice that $r = \sqrt{x^2 + y^2 + z^2}$ is an even function of $x$, $y$, and $z$. Thus, the integrand is odd with respect to $x$ (because of the extra factor of $x$), but even with respect to $y$ and $z$.* The entire integral therefore goes to zero.
The integrals for $\hat{p}_y$ and $\hat{p}_z$ proceed in exactly the same way.
* Although it also diverges at the origin. I'd be very glad if someone could tell me how to rigorously resolve this.
