Usually, taking the "average" of matrices would mean the usual arithmetic mean, $A = \frac 1n\sum_{i=1}^n A_i$. The problem with this, however, is that there is no reason to believe that this arithmetic mean would be another rotation.
One solution to this problem is to use a representation of rotations within a space that allows for linear combinations (and hence allows for arithmetic means). Jyrki's suggestion is to use the unit quaternion representation. I will suggest using a matrix logarithm.
The idea is as follows: the matrix exponential defines an onto map to $SO(3)$ from the skew-symmetric matrices (the "Lie Algebra" of $SO(3)$). If we reverse this map, average the resulting skew-symmetric "representations" of the matrix, then use the matrix exponential on this average, then in a sense we have the "average" rotation. (Aside: the exponential map over a Lie algebra is generally onto if the corresponding Lie group is connected and compact)
For a unit vector $\mathbf v = (v_x,v_y,v_z)$, let $K_\mathbf v$ denote the corresponding cross-product matrix,
$$
K_{\mathbf v} = \pmatrix{
0&-v_z & v_y\\
v_z&0&-v_x\\
-v_y & v_x & 0}.
$$
The Rodrigues rotation formula tells us that
$$
\exp(\theta K_{\mathbf v}) = I + \sin \theta K_{\mathbf v} + (1 - \cos \theta)K_{\mathbf v}^2.
$$
Reversing this formula is tricky, but here's an approach that works as long as $\theta \neq 180^\circ$. In order to reverse this formula and find values of $K$ and $\theta$ corresponding to a given rotation $R$, we can get the angle of the rotation using the of the rotation:
$$
\operatorname{Tr}(R) = 1 + 2 \cos \theta \implies \theta = \arccos\left[ \frac{\operatorname{Tr}(R) - 1}{2}\right].
$$
From there, we can use the skew-symmetric part of $R$ to find the axis:
$$
R - R^T = 2 \sin \theta K_{\mathbf v} \implies K_{\mathbf v} = \frac 1{2 \sin \theta}(R - R^T).
$$
This fails if $\theta = 0$, but in that case we can simply give the zero-matrix as the logarithm. Putting it all together, we have
$$
\theta(R) \arccos\left[ \frac{\operatorname{Tr}(R) - 1}{2}\right], \quad \log(R) = \begin{cases}
I & \theta(R) = 0,\\
\frac 1{2 \sin \theta}(R - R^T) & 0 < \theta(R) < \pi,\\
??? & \theta(R) = \pi.
\end{cases}
$$
So, given rotation matrices $R_1,\dots,R_k$, we compute the skew-symmetric matrix
$$
X = \sum_{j=1}^k \log(R_k).
$$
From there, you could compute the matrix exponential in the usual way. Alternatively, you could apply the Rodrigues formula. Note that if $X = \theta K_{\mathbf v}$, then $\theta = \sqrt{\frac 12 \text{tr}(X)}$ and $K_{\mathbf v} = \frac 1{\theta} X$.
Alternatively, I suspect that the following might lead to more consistent behavior of the average: if the rotations are all close to each other, then we can use one of the rotations as a reference point and, in essence, use the tangent space to $SO(3)$ at that point rather than the tangent space at the identity.
For a given reference rotation $R_0$ and other rotations $R_1,\dots,R_k$, this amounts to the following: compute
$$
X = \frac 1k \sum_{j=0}^k \log(R_0^T R_j)
$$
(of course we always have $R_0^TR_j = I$ for $j = 0$), and then take your "average" rotation to be $R_0 \exp(X)$.
One advantage of this method is that, as long as the rotations are close to each other, we can be sure that we won't need the logarithm in (or near) the $\theta = \pi$ case.
Another alternative to consider: you could take the usual average of the matrices, $A = \frac 1n\sum_{i=1}^n A_i$, and then get the nearest orthogonal matrix to this average. As in the solution to Wahba's problem, you can obtain this answer using the singular value decomposition of $A$. If $A$ has the SVD $A = U \Sigma V^T$, then the nearest element of $SO(3)$ is given by $R = U \tilde \Sigma V^T$, where
$$
\tilde \Sigma = \pmatrix{1&0&0\\0&1&0\\0&0&\det(U)\det(V)}.
$$
Notably, if the matrices $A_i$ are sufficiently close to each other, then the continuity of the determinant implies that we'll have $\det(U)\det(V) = \text{sgn}(\det(A)) = 1$, so that $\tilde \Sigma = I$
Also noteworthy: in the case that $\tilde \Sigma = I$, this nearest element of $SO(3)$ is also the orthogonal matrix $U$ from the polar decomposition of $A$.
