This is not an answer, but an amusing observation. If I didn't screw up my calculus, for the uniform distribution on the spherethe integral has a value of $\frac{\pi}{48} \approx 0.06545...$. On the other hand, since the volume is maximized when the points are orthogonal, let us try a measure concentrated at the vertices of an orthonormal basis, i.e. $\mu = \frac13 \big( \delta_{e_1}+\delta_{e_2}+ \delta_{e_3} \big)$. The expected value of the volume for this distribution will be $\frac1{27} =0.037...$, since $6$ out of $27$ choices yield a tetrahedron of volume $\frac16$, and the rest give volume $0$. We see that the uniform distribution wins.
Notice that in 2 dimensions (for the area of triangle Apex angle of a triangle as a random variable) this difference is: $\frac1{\pi} \approx 0.3183...$ vs. $\frac14= 0.25$.
I suspect that for the analogous problem in $d$ dimensions ($d$-dimensional volume of a simplex generated by random points on $S^{d-1}$) the difference would go to zero as $d\rightarrow \infty$ (due to the phenomenon of concentration of measure), but the uniform measure would be the maximizer.
1) Here is the computation of the integral. I'll work with the expected volume of parallelepiped vs tetrahedron, so that I don't have to carry $\frac16$ around. We need to compute
\begin{equation}
\int\limits_{S^{d-1}}\int\limits_{S^{d-1}}\int\limits_{S^{d-1}} | \det ( x\, y\, z) | d\sigma(x)\, d\sigma(y) \, d\sigma(z) = \int\limits_{S^{d-1}} \int\limits_{S^{d-1}} | \det ( x\, y\, z) | d\sigma(x)\, d\sigma(y),
\end{equation}
where by rotational invariance we can fix one of the variables (e.g. $z$), thus reducing to a double integral above. Since this integral is independent of the value of $z$, let us choose $z= e_3 = (0,0,1)$. Then $$| \det ( x\, y\, z) | = |x_1 y_2 - x_2 y_1 | = | \sin (\phi_x - \phi_y) | \cdot \sin \theta_y \cdot \sin \theta_y,$$ where we switched to polar coordinates. Recalling that $$ d\sigma = \frac1{4\pi} \sin \theta d\theta d\phi,$$
(we divide by $4\pi$ since $\sigma $ is a normalized measure) we reduce our integral to:
\begin{align*}
\int\limits_{S^{d-1}} \int\limits_{S^{d-1}} | \det ( x\, y\, e_3) | d\sigma(x)\, d\sigma(y)& = \frac1{16\pi^2} \bigg(\int_0^\pi \sin^2 \theta \, d\theta \bigg)^2 \cdot \bigg( \int_0^{2\pi}\int_0^{2\pi} | \sin (\phi_x - \phi_y ) | \, d\phi_x d\phi_y \bigg)\\
& = \frac1{16\pi^2} \cdot \bigg( \frac{\pi}2 \bigg)^2 \cdot 8\pi = \frac{\pi}{8}.
\end{align*}
(The integrals above are straightforward to compute.) Therefore, the expected area of the tetrahedron is $\frac{1}{6} \cdot \frac{\pi}{8} = \frac{\pi}{48} \approx 0.06545...$.
2) Here is the derivation of the recursive formula for the expected volume in $d$ dimensions. Let $E_d$ be the expected $d$-dimensional volume of the parallelepiped, spanned by $d$ independent uniformly distributed random points on the sphere $S^{d-1}$. Then the volume of the tetrahedron is simply $V_d = \frac{1}{d!} E_d$.
For any $x\in S^{d-1}$ let us write $x= (\sqrt{1-t^2} u, t)$, where $t\in [-1,1]$ is its $d^{th}$ coordinate, and $u\in \mathbb S^{d-2}$. It is a classical fact that under this change of variables $$ d\sigma (x) = \frac{\omega_{d-2}}{\omega_{d-1}} (1-t^2)^{\frac{d-3}2} dt \, d\sigma_{d-2} (u) , $$
where $\omega_{d-1} = \frac{2\pi^{d/2}}{\Gamma (d/2)}$ is the Lebesgue surface measure of the sphere $S^{d-1}$, and $\sigma_{d-2}$ is the normalized surface measure on $S^{d-2}$. Notice that this implies, in particular, that $$\int_{-1}^1 (1-t^2)^{\frac{d-3}2} dt = \frac{\omega_{d-1}}{\omega_{d-2}}. $$
We then use the same trick as above: fix $x_d = e_d = (0,0,\dots,0,1)$ to reduce the number of variables:
\begin{align*}
\mathbb E_{x_1,\dots, x_d \in S^{d-1}} | \det (x_1 \, \dots \, x_d ) | & = \mathbb E_{x_1,\dots, x_{d-1} \in S^{d-1}} | \det (x_1 \, \dots \, e_d ) |\\
& = E_{x_1,\dots, x_{d-1} \in S^{d-1}} \prod_{i=1}^{d-1} (1-t_i^2)^{1/2} \cdot \det (u_1 \, \dots \, u_{d-1} ) |\\
& = E_{t_1,\dots, t_{d-1} \in [-1,1] } \prod_{i=1}^{d-1} (1-t_i^2)^{1/2} \cdot E_{u_1,\dots, u_{d-1} \in S^{d-2}} | \det (u_1 \, \dots \, u_{d-1} ) |\\
& = \bigg( \frac{\omega_{d-2}}{\omega_{d-1}} \int_{-1}^1 (1-t^2)^{\frac{d-2}2} dt \bigg)^{d-1} \cdot E_{d-1}\\
& = \bigg( \frac{\omega_{d-2}\omega_{d}}{\omega_{d-1}^2} \bigg)^{d-1} E_{d-1} .
\end{align*}
Therefore
$$ E_d = \bigg( \frac{\omega_{d-2}\omega_{d}}{\omega_{d-1}^2} \bigg)^{d-1} E_{d-1} = \left( \frac{\big( \Gamma \big(\frac{d}{2}\big)\big)^2 }{\Gamma \big(\frac{d-1}{2}\big) \Gamma \big(\frac{d+1}{2}\big)} \right)^{d-1} E_{d-1}.$$
Let us check the case $d=3$. We have $$ E_{d-1}= E_2= \frac{1}{4\pi^2} \int_0^{2\pi}\int_0^{2\pi} |sin (\phi_1 - \phi_2) | d\phi_1 d\phi_2 = \frac{2}{\pi} ,$$ also $\omega_{d-1} = \omega_2 = 4\pi$, $\omega_{d-2} = \omega_1 = 2\pi$, $\omega_d = \omega 3 = 2\pi^2$. Therefore, $$ E_3 = \bigg( \frac{\omega_{1}\omega_{3}}{\omega_{2}^2} \bigg)^2 E_{2} = \bigg( \frac{2\pi \cdot 2\pi^2}{(4\pi)^2} \bigg)^2 \cdot \frac{2}{\pi} = \bigg( \frac{\pi}{4} \bigg)^2 \cdot \frac{2}{\pi} = \frac{\pi}{8}, $$ as we have obtained above.
