$\def\avint{\mathop{\,\rlap{-}\!\!\int}\nolimits}$
One line in a proof I'm looking at is $$\varphi(r) = \avint_{\partial B_r(x)} u(y) d\sigma(y) = \avint_{\partial B_1(0)} u(x+rz) d\sigma(z) $$
I am not sure how they have transformed the integral from $B_r(x)$ to $B_1(0)$.
I think it might have something to do with the fact that $$\partial B_R(x) = R\partial B_1(0) + x$$ which suggest a substitution like $y=rz + x$. But then not sure how to get the area element...
Any help is appreciated.
$\def\avint{\mathop{\,\rlap{-}\!\!\int}\nolimits}$
Your initial suspicions is correct: the integral follows from the change of variables $y = x + rz$.
We can argue that the area element for the averaged integral will not change because all we've done is re-center and scale the coordinate system. Re-centering doesn't change the area element at all. Re-scaling will, but since we are averaging over the total area, this is irrelevant.
Let's support this last paragraph by a calculation. Let $|S| = |\partial B_1(0)|$ be the area ($(n-1)$-dimensional Hausdorff measure) of the unit ball. Then, $$\avint_{\partial B_r(x)} u(y) d\sigma(y) = \frac{1}{r^{n-1} |S|}\int_{\partial B_r(x)} u(y)\;d\sigma(y)$$ If we take the change of variables $y = rz + x$, then the new area element is $d\sigma(y) = r^{n-1}d\sigma(z)$ (since we are working with an $(n-1)$-dimensional area, the differential $r^{n-1}$). Thus, after doing the change of variables, we get $$\avint_{\partial B_r(x)} u(y) d\sigma(y) = \frac{r^{n-1}}{r^{n-1}|S|} \int_{\partial B_1(0)} u(rz +_ x)\;d\sigma(z) = \avint_{\partial B_1(0)} u(rz + x) \;d\sigma(z)$$