Here's a hint to get you started:
Define $A(r) = \int_{\partial B(y,r)} u(x) \, dx$, and show that $A(r)$ is constant by differentiating with respect to $r$, and showing that the derivative is zero by using the divergence theorem to replace the integrand with a $\Delta u$.
For the subharmonic and superharmonic cases, the same technique leads to $A(r)$ increasing and decreasing, respectively.
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.$$
Best Answer
This is a re-interpretation of Ray Yang's answer, which also shows how the result can be generalized to other polygons. Introduce the function $$v(x,y)=(1-\max(|x|,|y|))^+ ,\qquad (x,y)\in\mathbb R^2$$ This is a compactly supported Lipschitz function with support $Q=[-1,1]\times [-1,1]$. Its graph is a pyramid with $Q$ as the base.
If $u$ is harmonic in a neighborhood of $Q$, then integration by parts yields $$0=\int_{\mathbb R^2} v\,\Delta u = \int_{\mathbb R^2} u\,\Delta v \tag{1}$$ By considering $u(\alpha x,\alpha y)$ with $\alpha\to 1^-$, we extend (1) to functions continuous in $Q$ and harmonic in its interior.
It remains to observe that $\Delta v$ is the distribution composed of
This follows from considering the discontinuities of the normal derivative of $v$ across the aforementioned lines; elsewhere $v$ is harmonic. One can also save the trouble of calculating the factor of $-\sqrt{2}$ by using the fact that $\int_{\mathbb R^2}\Delta v=0$.
It should be clear that there is nothing special about the square and its diagonals: any piecewise linear compactly supported function gives rise to a similar identity. For example, one can build a pyramid on top of any regular polygon $P$ and conclude that the average of a harmonic function along $\partial P$ is equal to its average over the union of segments connecting the vertices of $P$ to its center. No computations of slopes are necessary: it's clear that $(\Delta v)^+$ and $(\Delta v)^-$ are constant multiples of linear measure, and the identity $\int_{\mathbb R^2}\Delta v=0$ gives all the information we need.