Say you have a real unit $(M+1)$-vector $r_M$ which parameterizes some part of the $M$-sphere and is a function of parameters $\theta_1$ through $\theta_M$. How do we find the $M$-dimensional ``surface area" $S_M$ of vector $r_M$? For $M=2$ I know we can use \begin{equation} S_2=\int\Big|\frac{\partial r_2}{\partial \theta_1}\times\frac{\partial r_2}{\partial \theta_2}\Big|dA \end{equation}
But the cross product is only defined for $3$ (and $7$?) dimensions. What is the generalization I'm looking for here? My best guess is: \begin{equation} S_M=\int\sqrt{\sum_{i=1}^{M+1}\left(\sum_{j=1}^{M+1}\epsilon_{ij_1...j_m}\frac{\partial r_{j_1}}{\partial\theta_1}...\frac{\partial r_{j_M}}{\partial\theta_M}\right)^2}d\theta_1...d\theta_M \end{equation}
where $\epsilon$ is the Levi-Civita tensor of rank $M+1$.
The generalization of surfave integrals in higher dimensions is done using differential forms. The volume form on the euclidian sphere $\mathbb S^{N-1}\subseteq \mathbb R^{N}$ is $\omega= \sum_{i=1}^N (-1)^{i-1}x_i \text{d}x_1\wedge\dots\wedge\widehat{\text{d}x_i}\wedge\dots\wedge \text{d}x_N$. More on that in a question I asked, but since I don't know your background I'm not going to assume that knowledge and take a different approach.
In 2-dim surfaces in $\mathbb R^3$, the vector $\frac{\partial r}{\partial \theta}\times\frac{\partial r_2}{\partial \theta_2}$ represented a normal to the surface whose length measures the area of a small parallelogram. Integrating all of them gives the total area. Well, the cross product could be generalized to any dimension. The 3-7 case you are talking about is when you want the operation to be binary! We usually define the cross product of $v_1,\dots,v_{N}\in\mathbb R^{N+1}$
$$\Pi(v_1,\dots,v_{N})=\det\begin{bmatrix} e_1 & \dots & e_{N+1}\\ \cdots& v_1 & \cdots\\ \vdots & \ddots & \vdots\\ \cdots & v_{N} & \cdots \end{bmatrix}$$ With the regular abuse of notation of standard basis vectors in the first row. It's clear that the result is orthogonal to all $v_i$s since doing a dot product with any of it yields a determinant with twice the row $v_i$ which yields zero. Showing it's length rightfully represents volume depends on your definition of $N$ volume, but it's not hard to see it's linear and permutation invariant at least as it should. To calculate the length, we see the $i$th component is $$\det\begin{bmatrix} v_{1}^1& \cdots & v_{1}^{i-1} & v_{1}^{i+1}&\cdots& v_{1}^N\\ \vdots & & \vdots&\vdots & &\vdots\\ v_{N}^1& \cdots & v_{N}^{i-1} & v_{N}^{i+1}&\cdots& v_{N}^N \end{bmatrix}=\epsilon_{i,j_1,\dots,j_N} v_{1}^{j_1}\cdots v_{N}^{j_N}$$
With this in mind we can do the same proccess in 2 dimensions and get that for an $N$-dim surface $\Sigma$ parameterized with $r$ on domain $\Omega$ $$\text{vol}_N(\Sigma)=\int_\Omega \left \| \Pi\left ( \frac{\partial r}{\partial \theta_1},\dots,\frac{\partial r}{\partial \theta_N} \right ) \right \|_2 \text{d}\theta$$ Which gives an expression I reckon is close to yours but not precise, the $j$s are in the $\theta$ indices.