Idea
I have an idea, that it is possible to generalize the winding number for surfaces of the form $f: \mathbb{R}^n \rightarrow \mathbb{R}^{n+1}$
The winding number for $n=1$ is $w_{\gamma}(x) = \oint_{\gamma} \frac{\left < L r, \dot r\right >}{|r|^2} dt$, where $L = \begin{bmatrix} 0 & -1\\ 1 & 0 \end{bmatrix}$ and $r = \gamma - x$.
I would say that $w_{\gamma}$ is similar to this calculation:
$$\hat{w}_{\gamma}(x) = \oint_{\gamma}\left < \frac{d}{dt} \frac{r}{|r|}, L \frac{r}{|r|}\right > dt$$
We can think of $\frac{r}{|r|}$ as a projection of the curve $\gamma$ to a unit circle centered at $x$. The $\frac{d}{dt} \frac{r}{|r|}$ can be thought as the tangent vector of the projection to this unit circle. The value $\left |\frac{d}{dt} \frac{r}{|r|}\right |$ is the infinitesimal angle of this projection.
So $\hat{w}_{\gamma}(x)$ can be though as the signed angle.
Questions
Does $\hat{w}_{\gamma}(x)$ has the same meaning as $w_{\gamma}(x)$?
Can we apply the same approach of $\hat{w}_{\gamma}(x)$ to generalize to surfaces $f: \mathbb{R}^n \rightarrow \mathbb{R}^{n+1}$?
What is the main approach among mathematicians to generalize the winding number?
First question attempt
I tried to change the form of $\hat{w}_{\gamma}(x)$:
$$\hat{w}_{\gamma}(x) = \oint_{\gamma}\left < \frac{d}{dt} \frac{r}{|r|}, L \frac{r}{|r|}\right > = \oint_{\gamma}\frac{1}{|r|^2}\left < \text{ort}_{r}\dot r, L \frac{r}{|r|}\right > = \oint_{\gamma}\frac{1}{|r|^4}\left (\left < \dot r, L r\right > |r|^2 - \left < r, \dot{r} \right >\left < r, Lr\right >\right )$$
$$\hat{w}_{\gamma}(x) = w_{\gamma}(x) - \oint_{\gamma}\frac{1}{|r|^4}\left< r, \dot r \right>\left< r, Lr\right>\;\;(1)$$
From (1) $\hat{w}_{\gamma}(x)$ has something to do with $w_{\gamma}(x)$.
Second question attempt
I tried to generalize to $n=2$:
Let $\gamma: U \in \mathbb{R}^2 \rightarrow \mathbb{R}^3$ be the parameterization of a closed surface in $\mathbb{R}^3$.
Then $\hat{w}_{\gamma}(x)$ will be:
$$\hat{w}_{\gamma}(x) = \iint_{\gamma} \star\hat r^{\flat} \left(\partial_{u_1} \hat r, \partial_{u_2} \hat r \right)$$
Where is $\star$ is the Hodge star operator, $\flat$ is the flat musical isomorphism, $\hat r = \frac{r}{|r|}$ and $\partial_{u_k} f = \frac{\partial f}{\partial u_k}$.
I am not a mathematician so please correct anything that is wrong, especially the differential forms parts.
Thanks in advance.
Alright, let's see if I can flesh out my comments.
Winding numbers from algebraic topology
Let $M$ be a smooth, orientable, closed (by which I mean compact) $n$-manifold. Fix an orientation on $M$. Let $f : M \to \newcommand\RR{\Bbb{R}}\RR^{n+1}$ be a smooth map. Let $x$ be a point outside the image of $M$. Without loss of generality, we can assume that $x$ is the origin. Now we have a smooth map from $\RR^{n+1}\setminus \{0\}$ to $S^n$ defined by $r\mapsto \frac{r}{\|r\|}$. Let's call this map $q$.
Taking the composition, we get a smooth map $g=q\circ f : M\to S^n$.
$g$ induces a pushforward map on homology, in particular on top level homology, giving a map $g_*:H_n(M)\to H_n(S^n)$. Since we've fixed an orientation on $M$, we have a fundamental class $[M]$, which pushes forward to a class in $H_n(S^n)$. If we fix the orientation on $S^n$ determined by outward facing normal vectors, then $H_n(S^n)\simeq \Bbb{Z}$, and $g_*[M]$ corresponds to some integer in $\Bbb{Z}$, which we can call the winding number of the map $f$ about the origin.
For example, the winding number of the inclusion of the unit sphere in $\Bbb{R}^{n+1}$ (with the outward normals orientation) is $1$, and the winding number of the inclusion of any sphere that doesn't contain the origin into $\Bbb{R}^{n+1}$ is $0$.
Computing winding numbers
You'll notice that I didn't make any use of the assumption that everything was smooth above.
We'll make use of that now.
Disclaimer: My last differential topology course was four years ago now, and it was introductory, so I can't be certain that I'm recalling things correctly. Please let me know if there are any errors in the comments, and I'll fix them.
Now we want to find which integer the class $g_*[M]$ represents in $H_n(S^n)$. Well, if $1$ is the fundamental class of $S^n$ with respect to the outwards normal orientation, and $[\omega] \in \newcommand\HdR{H_{\textrm{dR}}}\HdR^n(S^n)$ is the corresponding class of the orientation form in top-level de Rham cohomology, then $$\int_1 [\omega] = 1,$$ (with $\omega$ appropriately normalized) and $$\int_{g_*[M]}[\omega] = n,$$ where $n$ is the integer such that $n\cdot 1 = g_*[M]$.
Thus we've found an integral that allows us to compute winding numbers, unfortunately, we can't use it directly, since the whole point is that we don't know what $g_*[M]$ is.
However we can pull $\omega$ back to $M$ to get: $$\int_{g_*[M]}\omega = \int_M g^* \omega,$$ and we can compute $g^*\omega$.
With a parametrization
Now suppose we have a parametrization $\gamma : \RR^n\supseteq D \to\RR^{n+1}$, whose image is our smooth manifold $M$, realized as a submanifold of $\RR^{n+1}$, then $$\int_M g^*\omega=\int_D \operatorname{det}\left( \partial_1\frac{\gamma(\mathbf{x})}{\|\gamma(\mathbf{x})\|},\ldots, \partial_i\frac{\gamma(\mathbf{x})}{\|\gamma(\mathbf{x})\|},\ldots, \partial_n\frac{\gamma(\mathbf{x})}{\|\gamma(\mathbf{x})\|}, \frac{\gamma(\mathbf{x})}{\|\gamma(\mathbf{x})\|} \right)\, d\mathbf{x}, $$ ignoring normalization. Now $$\partial_i\frac{\gamma}{\|\gamma\|}=\frac{1}{\|\gamma\|}(\partial_i\gamma) -\frac{\gamma\cdot (\partial_i\gamma)}{\|\gamma\|^3}\gamma,$$ and the last column in our determinant is $\frac{\gamma}{\|\gamma\|}$, so we can add $\frac{\gamma\cdot(\partial_i\gamma)}{\gamma\cdot \gamma}$ times the last column to each of the first $n$ columns to get that in fact, $$\int_M g^*\omega=\int_D \operatorname{det}\left( \frac{\partial_1\gamma(\mathbf{x})}{\|\gamma(\mathbf{x})\|},\ldots, \frac{\partial_i\gamma(\mathbf{x})}{\|\gamma(\mathbf{x})\|},\ldots, \frac{\partial_n\gamma(\mathbf{x})}{\|\gamma(\mathbf{x})\|}, \frac{\gamma(\mathbf{x})}{\|\gamma(\mathbf{x})\|} \right)\, d\mathbf{x} $$ $$=\int_D \frac{1}{\|\gamma(\mathbf{x})\|^{n+1}}\operatorname{det}\left( {\partial_1\gamma(\mathbf{x})},\ldots, {\partial_i\gamma(\mathbf{x})},\ldots, {\partial_n\gamma(\mathbf{x})}, {\gamma(\mathbf{x})} \right)\, d\mathbf{x}. $$
Note that we have rederived both formulas that you've provided for computing the winding number when $n=1$.
I'm not very familiar with the Hodge star, but it appears that this formula might also agree with your 2d formula as well.
Example:
Now let's use it to do an example calculation.
We can parametrize the 2-sphere by $(\theta,\phi)\mapsto (\cos\theta\sin\phi,\sin\theta\sin\phi,\cos\phi)$ for $(\theta,\phi)\in [0,2\pi]\times [0,\pi]$.
Then, we can compute the matrix we need to take the determinant of, it's $$\begin{pmatrix}\sin\theta \sin\phi & -\cos\theta\cos\phi & \cos\theta \sin\phi \\ -\cos\theta\sin\phi & -\sin\theta\cos\phi & \sin\theta\sin\phi \\ 0 & \sin\phi & \cos\phi \end{pmatrix},$$ which has determinant $$\sin\theta\sin\phi (-\sin\theta\cos^2\phi-\sin\theta\sin^2\phi) +\cos\theta\sin\phi (-\cos\theta\cos^2\phi -\cos\theta\sin^2\phi)$$ $$=-\sin\phi(\sin^2\theta+\cos^2\theta)=-\sin\phi.$$
Then $$\int_0^{2\pi} \int_0^\pi -\sin\phi \,d\phi\,d\theta = 2\pi \left[\cos\phi \right]_0^\pi = 2\pi (-1-1)=-4\pi.$$
The negative sign should reflect that my parametrization of the sphere is backwards. Also the fact that the constant is $-4\pi$ reflects the choice not to normalize.
Note also that if I increase the domain of $\theta$ from $[0,2\pi]$ to $[0,4\pi]$, we see that $\gamma$ now parametrizes the sphere wrapped around itself twice, and the value of the integral doubles, as it should.
Note: I have to go now, so I can't finish this answer to the extent that I'd like, but I'd say it's mostly complete to my satisfaction at this point.