Well first we should note that the result you state is not true without some smoothness assumptions on $u$. Let's assume $u$ is continuous. We also assume $u$ is real-valued.
There's a standard proof of the result for surface averages that doesn't involve the things you're concerned about and that works just as well for volume
averages:
Suppose $\overline{B(0,1)}\subset\Omega$. We will show that $u$ is harmonic in $B=B(0,1)$. Let $\phi$ be the restriction of $u$ to the boundary of the unit ball, and let $$f=P[\phi],$$the Poisson integral of $\phi$. Then $f$ is harmonic in $B$ and extends continuously to $\overline B$, with values on the boundary given by $\phi=u$. We will let $f$ denote the extension of the original $f$ to $\overline B$. Now define $v:\overline B\to\mathbb R$ by $$v(x)=u(x)-f(x)\quad(x\in\overline B).$$If we can show $v=0$ in $B$ we're done, since $f$ is harmonic in $B$.
Suppose that $v$ does not vanish identically in $B$. Wlog assume that $v$ is strictly positive at some point of $B$. Define $$M=\sup_{x\in B}v(x)$$ and $$K=\{x\in B\,:\,v(x)=M\}.$$
Since $v$ is continuous in $\overline B$ and vanishes on the boundary it follows that $K$ is a nonempty compact subset of $B$. Let $p$ be a point of $K$ at minimal distance to the boundary of $B$.
Now $v$ satisfies the volume mean-value property in $B$ since both $u$ and $f$ do. So for small $r>0$ we have $$M=v(p)=\frac1{m(B(p,r))}\int_{B(p,r)}v.$$But this is impossible: $v\le M$ everywhere in $B(p,r)$ and $v<M$ in a nonempty open subset of $B(p,r)$, so $$\frac1{m(B(p,r))}\int_{B(p,r)}v<M.$$
Why is that boundary integral in the last formula zero?
What about using $q(z) =\frac12( \|z\|^2 -1 )$? Then $\nabla q = z$ and $-q(z) \ge0$ for $\|z\|\le1$, moreover $q(z) =0$ for all $\|z\|=1$. And that last inequality becomes
\begin{align*}
\varphi'(r)
&= \frac{1}{r|B_1(0)|}\int_{B_1(0)} \langle \nabla v(z),\nabla q(z)\rangle\,{\rm d}{z} \\
&= \frac{1}{r|B_1(0)|}\left( -\int_{B_1(0)} q(z)\triangle v(z)\,{\rm d}{z} + \int_{\partial B_1(0)} q(z)\frac{\partial v}{\partial \nu}(z)\,{\rm d}{S}(z)\right) \\
&= \frac{1}{r|B_1(0)|}\left( \int_{B_1(0)} \frac{1-\|z\|^2}2\triangle v(z)\,{\rm d}{z} + 0\right) \ge0
\end{align*}
Best Answer
Theorem: Let $\Omega \subset \mathbb{R}^N$, $u \in C(\Omega)$ be such that $$\frac{1}{|B(x_0,R)|}\int_{B(x_0,R)}u(y)\ dy = u(x_0) = \frac{1}{|\partial B(x_0,R)|}\int_{\partial B(x_0,R)}u\ dS$$ for every ball $\overline{B(x_0,R)} \subset \Omega$. Then $u \in C^{\infty}(\Omega)$ and it is harmonic
Proof: Consider the standard mollifier: $$\rho(x) := \begin{cases}Ce^{-\frac{1}{1 - \|x\|^2}} & \text{if $\|x\|$ < 1} \\0 & \text{otherwise.} \end{cases}$$ Here $C$ is a constant such that $\|\rho\|_{L^1} = 1.$ Let $\epsilon > 0$ and consider $$\rho_{\epsilon}(x) = \epsilon^{-N}\rho(x\epsilon^{-N}).$$ Set $\Omega_{\epsilon} = \{x \in \Omega : \text{dist}(x,\partial \Omega) > \epsilon\}$ and define for $x \in \Omega_{\epsilon}$$$u_{\epsilon}(x) = \rho_{\epsilon} * u(x) = \int_{\Omega}\rho_{\epsilon}(x - y)u(y)\ dy.$$ The following is a well know theorem in analysis, if it is new to you you can look for a proof Analysis by Lieb and Loss or anywhere else.
**Theorem:***If $u \in C(\Omega)$, then $u_{\epsilon} \to u$ uniformly on compact subsets of $\Omega$, $u_{\epsilon} \in C^{\infty}(\Omega_{\epsilon})$ and for any multindex $\alpha$ we have $$\frac{\partial^{\alpha}u_{\epsilon}}{\partial x^{\alpha}}(x) = \int_{\Omega}\frac{\partial^{\alpha}\rho_{\epsilon}}{\partial x^{\alpha}}(x - y)u(y)\ dy.$$*
Finally we can proceed with the proof!
Fix $x_0 \in \Omega_{\epsilon}$. $$u_{\epsilon}(x_0) = \int_{B(x_0,\epsilon)}\rho_{\epsilon}(x - y)u(y)\ dy = \int_{B(0,\epsilon)}\rho_{\epsilon}(z)u(x_0 - z)\ dz = $$ $$ = \int_0^{\epsilon}r^{N - 1}\int_{\partial B(0,1)}\rho_{\epsilon}(rw)u(x_0 - rw)\ dS(w)dr = $$ $$ \int_0^{\epsilon}r^{N - 1}\rho(r)\int_{\partial B(0,1)}u(x_0 - rw)\ dS(w)dr = \int_0^{\epsilon}r^{N-1}\rho_{\epsilon}(r)\frac{\alpha_N N}{|\partial B(x_0,r)|}\int_{\partial B(x_0,r)}u(y)\ dS(y)dr $$ $$ = u(x_0)\|\rho\|_{L^1} = u(x_0).$$
This proves that $u = u_{\epsilon}$ and hence $u \in C^{\infty}(\Omega_{\epsilon})$, for every $\epsilon$. We are left to prove that $u$ is harmonic. To this end consider $$f(r) = \frac{1}{|\partial B(x_0,r)|}\int_{\partial B(x_0,r)}u\ dS.$$ By assumption, $f$ is constant, hence $f' \equiv 0$. By the divergence theorem it is easy to show (and I am sure that you have already seen this result since you have proved that if $u$ is harmonic then it satisfies the mean value property) $$0 = f'(r) = \text{constant}\int_{B(x_0,r)}\Delta u(y)\ dy \to \Delta u(x_0),\ \text{if we let}\ r \to 0^+.$$ Thus $u$ is harmonic. QED