How to change variables in a surface integral without parametrizing This is a doubt that I carry since my PDE classes. 
Some background (skippable):
In the multivariable calculus course at my university we made all sorts of standard calculations involving surface and volume integrals in $R^3$, jacobians and the generalizations of the fundamental theorem of calculus. In order to make those calculations we had to parametrize domains and calculate differentials. 
A couple of years later I took a PDE course. We worked with Evans' Partial differential equations book. This was my first experience with calculus in $\mathbb R^n$ and manipulations like 
$$\text{average}\int_{B(x,r)}f(y)\,dy= \text{average}\int_{B(0,1)}f(x+rz)\,dz.$$
This was an ordinary change of variables. $y=x+rz,\,\,dy=r^n\,dz$ and the mystery was solved. Like in that case, I was able to justify most of these formal manipulations after disentangling definitions. 
That aside, I found these quick formal calculations to be very powerful.
However,
I realized that I wasn't able to justify this:
$$\text{average} \int_{\partial B(x,r)}f(y)dS(y)= \text{average}\int_{\partial B(0,1)}f(x+rz)\,dS(z).$$
I have some vague idea of what's happening: the same substitution as before, but this time the jacobian is $r^{n-1}$ because the transformation is actually happening between regions which "lack one dimension". Also, I see some kind of pattern: a piece of arc-length in the plane is $r\,d\theta$, a piece of sphere-area is $r^2 \sin\theta \, d\phi \,d\theta$, "and so on". Maybe some measure-theoretic argument can help me: I know, roughly speaking, that for any measure $\mu$, $$\int_\Omega f\circ \phi \,d\mu=\int_{\phi(\Omega)} f \, d(\mu\circ\phi^{-1}).$$ I'd say $\phi(z)=(z-x)/r$ and $\phi^{-1}(y)=ry+x$, but I actually don't know how $dS(y)$ looks like "as a measure" (It's not a product measure or a restriction of one, but it somehow relates to Lebesgue's in $\mathbb R^n$...). Why would I conclude that $dS(y)\circ \phi^{-1}=r^{n-1}dS(z)$? I have an intuition, but either I lack the mathematical concepts and definitions to express it or I'm just too confused. Is there some theory that I could learn in order to understand? Maybe something about the measure $dS$. Is it expressible in terms of the Lebesgue measure in some way? Or set-theoretically, maybe, without having to resort to $n-1$ parameters and complicated relations?
Maybe all of this would not have been a problem if I had ever mastered n-dimensional spherical coordinates. But even so, more generally, is there a way of changing variables when I'm integrating over a subregion of "dimension$<n$" without necessarily parametrizing? 
Sorry for the vagueness, but I don't really know what to ask for exactly.
Note: I saw some of the answers to this post, but none of them were deep enough in the direction I intend.
Note II: If there are no general methods or theories, maybe restricting to linear transformations, to Lebesgue measure exclusively, or to subregions defined by simple expressions like $g(x)=C$ or $g(|x|)=C$ could get me somewhere.

Edit: I have not yet studied differential geometry, which has been mentioned in a comment. I added it to the tags.
 A: No parametrization is needed, just some careful uncovering of definitions. Maybe the following can help: 
Let $S^n(r)=\{x\in\mathbb R^{n+1}|\|x\|=r\}\subset\mathbb R^{n+1}$ (the sphere of radius $r>0$ centered at the origin). 
There is a natural volume form $\mu_{r}$ on  $S^n(r)$,  i.e. a nowhere-vanishing $n$-form, induced by the volume form $\mu=dx_1\wedge\ldots \wedge dx_{n+1}$ on  $\mathbb R^{n+1}$ and the  "outward" unit normal $N:S^n(r)\to\mathbb R^{n+1}$,  $N(x)=x/r$, and is defined as follows: 
$$\mu_r(x)(v_1,\ldots,v_n):=\mu(x)(N(x), v_1, \ldots, v_n)$$ for all $v_1, \ldots, v_n\in T_x S^n(r).$ 
(More informally, we write $\mu_r=\mu/dr$. Note that $dr(N)=1$). 
Next let  $\Phi:S^n(1)\to S^n(r)$ be defined by $\Phi(x)=rx$. Then $$\Phi^*\mu_r=r^n\mu_1. $$ This you can prove without any parametrization, just from the definitions above and generalities like  the chain rule, pull back, etc. 
It then follows, by the change of variables formula, that for any continuous function $f:S^n(r)\to\mathbb R$, $$\int_{S^n(r)}f\mu_r=\int_{S^n(1)} \Phi^*(f\mu_r)=r^n\int_{S^n(1)} (f\circ\Phi) \mu_1.$$
Makes sense? 
A: I think the measure theoretic approach works fine, note that the surface measure is n-1 dimensional Hausdorff measure, in general for the s-dimensional Hausdorff measure we have $H^{s} (rA)=r^{s} H^{s}(A)$ and this measure is translation invariant. Now use the measure theoretic change of variable formula.
A: I know this is an old question, but I thought this explanation might be helpful to some.

By definition (in $\mathbb R^3$):
$$\int_{\partial B(\boldsymbol x,r)}f(\boldsymbol y)dS(\boldsymbol y)=
\int_U f(\boldsymbol y(s,t))\left\|\frac{\partial\boldsymbol y}{\partial s}\times\frac{\partial\boldsymbol y}{\partial t}\right\|dsdt$$
Now, observe that $f(\boldsymbol y)=f(\boldsymbol x+r(\frac{\boldsymbol y-\boldsymbol x}{r}))$, and that if $\boldsymbol y(s,t)$ is a parametrization of $\partial B(\boldsymbol x,r)$ for $(s,t)\in U$, then $\frac{\boldsymbol y(s,t)-\boldsymbol x}{r}$ is a parametrization of $\partial B(\boldsymbol 0,1)$ for $(s,t)\in U$. Finally we observe that
$$\left\|\frac{\partial\boldsymbol y}{\partial s}\times\frac{\partial\boldsymbol y}{\partial t}\right\|=
r^2\left\|\frac{\partial}{\partial s} \left (\frac{\boldsymbol y-\boldsymbol x}{r} \right )\times\frac{\partial }{\partial t} \left (\frac{\boldsymbol y-\boldsymbol x}{r} \right )\right\|$$
So if we let $\boldsymbol z(s,t)=\frac{\boldsymbol y(s,t)-\boldsymbol x}{r}$, then we have
$$\int_U f(\boldsymbol y(s,t))\left\|\frac{\partial\boldsymbol y}{\partial s}\times\frac{\partial\boldsymbol y}{\partial t}\right\|dsdt=
r^2\int_U f(\boldsymbol x +r\boldsymbol z(s,t))\left\|\frac{\partial\boldsymbol z}{\partial s}\times\frac{\partial\boldsymbol z}{\partial t}\right\|dsdt\\=
r^2\int_{\partial B(\boldsymbol 0,1)}f(\boldsymbol x+r\boldsymbol z)dS(\boldsymbol z)$$

Edit by OP
As @user5753974 commented, you can generalize this if you use the fact that in $\mathbb R^n$ $$∫_{∂B(\boldsymbol x,r)}f(\boldsymbol y)dS(\boldsymbol y)=∫_{U}f(\boldsymbol y(\boldsymbol z)) \left \|\det\left (\frac{∂\boldsymbol y}{∂z_1},…,\frac{∂\boldsymbol y}{∂z_{n−1}},\boldsymbol n\right) \right \| d^{n−1}\boldsymbol z,$$ where $\boldsymbol n$ is the normal vector to the surface, and that $\boldsymbol n$ does not change when the surface is scaled and translated.
