Small angular displacements with a quaternion representation

Small angular displacements with a quaternion representation

I have the orientation of a 3D spatial object represented by a unit
quaternion:
$$ q = a_1 + a_2 i + a_3 j + a_4 k $$ $$ \|q\| = (a_1^2 + a_2^2 + a_3^2 +
a_4^2)^{1/2} = 1 $$
I'd like to perturb this orientation slightly. If this displacment is
small, I think I can get away with
$$ q'_2 = q + \Delta (u_1 + u_2 i + u_3 j + u_4 k) $$
where $\Delta << 1$ and the $u$ are uniform random numbers in the range of
$[-1,1]$. $q'_2$ is obviously not a unit quaternion anymore so I would
normalize it by
$$ q_2 = \frac{q'_2}{\|q'_2\|} $$
As $\Delta$ gets larger, I feel that this will become an increasingly
worse approximation. Is there a better way to get a random orientational
displacement?