I would define it like this:
$$
\begin{align*}
q_{1}^{q_{2}} &= \left( q_{1} \right)^{q_{2}}\\
q_{1}^{q_{2}} &= \left( e^{\ln\left( q_{1} \right)} \right)^{q_{2}}\\
q_{1}^{q_{2}} &= e^{\ln\left( q_{1} \right) \cdot q_{2}}\\
q_{1}^{q_{2}} &= e^{\overbrace{\left(\ln\left( \left| q_{1} \right| \right) + \operatorname{sgn}\left( q_{1} - a_{1} \right) \cdot \arccos\left( \frac{a_{1}}{\left| q_{1} \right|} \right) \right) \cdot q_{2}}^{=~ q_{3}}}\\
q_{1}^{q_{2}} &= e^{q_{3}}\\
q_{1}^{q_{2}} &= e^{\Re\left( q_{3} \right)} \cdot \left( \cos\left( \left| q_{3} - \Re\left( q_{3} \right) \right| \right) + \sin\left( \left| q_{3} - \Re\left( q_{3} \right) \right| \right) \cdot \operatorname{sgn}\left( q_{3} - \Re\left( q_{3} \right) \right) \right)\\
\end{align*}
$$
where $\Re\left( q \right) = a \in \mathbb{R}$ and $\operatorname{sgn}\left( \cdot \right)$ is the Signum Function ($\operatorname{sgn}\left( q \right) = \frac{q}{\left| q \right|} \wedge q \ne 0$).
This definition would be based on a series expansion of $e^{z}$ around the point $z = 0$ ($e^{z} = \sum\limits_{k = 0}^{\infty}\left[ \frac{1}{k!} \cdot z^{k} \right]$). For a derivation of this see Exponential Function of Quaternion - Derivation.
A major problem that arises is that they are not commutative with respect to multiplication, so it makes a difference whether we left- or right-multiply $q_{2}$ in the exponent ($e^{q_{1} \cdot q_{2}} \ne e^{q_{2} \cdot q_{1}}$). Another problem is that $\left( x^{y} \right)^{z} = x^{y \cdot z}$ is not always valid, which we already know from complex numbers.
I didn't make up the formulas. They are also in the Wikipedia Article on the quaternions. I just recombined them:
$$e^{q} = e^{a} \cdot \left( \cos\left( \left| q \right| \right) + \sin\left( \left| q \right| \right) \cdot \operatorname{sgn}\left( q - a \right) \right)$$
$$\ln\left( z \right) = \ln\left( \left| z \right| \right) + \arccos\left( \frac{a}{\left| q \right|} \right) \cdot \operatorname{sgn}\left( q - a \right)$$
Properties that would emerge from the definition are (where $q_{n} = a_{n} + b_{n} \cdot i + c_{n} \cdot j + d_{n} \cdot k$):
$$
\begin{align*}
\Re\left( q' \right) = \Re\left( q_{1}^{q_{2}} \right) &= e^{a_{3}} \cdot \cos\left( \sqrt{b_{3}^{2} + c_{3}^{2} + d_{3}^{2}} \right)\\
q' - \Re\left( q' \right) &= e^{a_{3}} \cdot \sin\left( \sqrt{b_{3}^{2} + c_{3}^{2} + d_{3}^{2}} \right) \cdot \frac{b_{3} \cdot i + c_{3} \cdot j + d_{3} \cdot k}{\sqrt{b_{3}^{2} + c_{3}^{2} + d_{3}^{2}}}\\
\end{align*}
$$
But this is kind of boring, so let's take another approach:
We know that $z^{a} = \operatorname{_{1}F_{0}}\left( -a;\, \cdot;\, z - 1 \right)$ holds, where $\operatorname{_{1}F_{0}}$ is a generalized hypergeometric function given by $\operatorname{_{1}F_{0}}\left( a;\, \cdot;\, z \right) = \sum\limits_{k = 0}^{\infty}\left[ \frac{\Gamma\left( a + k \right)}{\Gamma\left( a \right) \cdot k!} \cdot z^{k} \right]$ where $\Gamma\left( \cdot \right)$ is the Complete Gamma Function.
So:
$$q_{1}^{q_{2}} = \sum\limits_{k = 0}^{\infty}\left[ \frac{\Gamma\left( -q_{2} + k \right)}{\Gamma\left( -q_{2} \right) \cdot k!} \cdot q_{1}^{k} \right]$$
I'll leave generalizing this for quaternion as a little exercise for you.
