Consider a scalar function of several variables $f(\mathbf{\vec{x}})$. Let us first use Cartesian coordinates $\mathbf{\vec{x}} \, = \, x_1 \, \mathbf{\hat{i}} \,+\, x_2 \, \mathbf{\hat{j}} \,+\, x_3 \, \mathbf{\hat{k}} \,$. A Taylor expansion of $f(\mathbf{\vec{x}})$ around $\mathbf{\vec{y}} \, = \, y_1 \, \mathbf{\hat{i}} \,+\, y_2 \, \mathbf{\hat{j}} \,+\, y_3 \, \mathbf{\hat{k}} \,$ looks like
\begin{equation}
f(\mathbf{\vec{x}}) \,= \, f(\mathbf{\vec{y}}) \,+ \, \sum_{i=1}^3 \Delta_i \, \partial_{i} f (\mathbf{\vec{y}}) + \frac{1}{2} \sum_{i=1}^3 \sum_{j=1}^3 \Delta_i \, \Delta_j \, \partial_{i} \partial_{j} f (\mathbf{\vec{y}}) \,+ \, \frac{1}{3!} \sum_{i=1}^3 \sum_{j=1}^3 \sum_{k=1}^3 \Delta_i \, \Delta_j \, \Delta_k \, \partial_{i} \partial_{j} \partial_{k} f (\mathbf{\vec{y}}) \,+ \, ...
\end{equation}
where
\begin{equation}
\Delta_k \,=\, x_k\,-\,y_k \, .
\end{equation}
We should focus on the case $\mathbf{\vec{y}} \, = \, \mathbf{0} \,$. Notice that in terms of the $\mathbf{\vec{x}}$ dependence, the terms in the expansion are monomials of increasing order. Using spherical coordinates as
\begin{equation}
\mathbf{\vec{x}} \, = \, r \, \sin \theta \, \cos \phi \, \mathbf{\hat{i}} \,+\, r \, \sin \theta \, \sin \phi \, \mathbf{\hat{j}} \,+\, r \, \cos \theta \, \mathbf{\hat{k}} \,
\end{equation}
we immediately realize that $\frac{\Delta_i}{r} = \frac{x_i}{r} $ is a function of only angular variables.
Notice also that the integration of a product of an odd number of these $x_i$ over the surface of the sphere vanishes. With this we mean
\begin{equation}
\left. \int_0^{2\pi}\mathrm{d}\phi \int_0^{\pi} \sin \theta \, \mathrm{d}\theta \, x_{i_1} \, x_{i_2} \cdot \cdot \cdot x_{i_n} \right|_{r\,=\,a} \, = \, 0 \, , \, \mathrm{if} \, n \,\mathrm{is} \, \mathrm{odd} \, .
\end{equation}
Even more, the only terms that lead to non-vanishing integrals are of the form
\begin{eqnarray}
&& \left. \int_0^{2\pi}\mathrm{d}\phi \int_0^{\pi} \sin \theta \, \mathrm{d}\theta \, x_1^{2 n_1} x_2^{2 n_2} x_3^{2 n_3} \right|_{r\,=\,a}
\\
&&= a^{2n_1+2n_2+2n_3}\int_0^{2\pi}\mathrm{d}\phi \int_0^{\pi} \mathrm{d}\theta \, \sin^{2 n_1 + 2 n_2 +1 } \theta \cos^{2 n_3} \theta \cos^{2 n_1} \phi \sin^{2 n_2} \phi
\\
&&= a^{2n_1+2n_2+2n_3} \frac{4\pi}{(2n_1+2n_2+2n_3+1)}\frac{(n_1+n_2+n_3)!}{n_1!n_2!n_3!}\frac{(2n_1)!(2n_2)!(2n_3)!}{(2n_1+2n_2+2n_3)!}
\end{eqnarray}
with $n_1,n_2$ and $n_3$ non-negative integers. This is one way of seeing why there is an absence of odd number of derivatives in the expansion. As you see, it does not have anything to do with the parity of the function $f(\mathbf{\vec{x}})$. From here, with some combinatorics taking care of the repeated terms being summed, one can deduce that the result for arbitrary $f(\mathbf{\vec{x}})$ is
\begin{equation}
\frac{1}{4 \pi} \left. \int_0^{2\pi}\mathrm{d}\phi \int_0^{\pi} \sin \theta \, \mathrm{d}\theta \, f(\mathbf{\vec{x}}) \, \right|_{r\,=\,a} \, = \sum_{k=0}^{\infty} \frac{a^{2k}}{(2k+1)!} \left( \nabla^2 \right)^k f(\mathbf{0}) \, .
\end{equation}
If $\left( \nabla^2 \right)^k f = 0 $ for $k>1$, the series truncates to your result. We did this for a scalar function but the result naturally extends for vector functions if you apply this process to each Cartesian component.
