Computing surface integral on the unit sphere

I'll use the spherical-coordinates approach, and I'll assume for now that by "surface measure" you mean the Haar measure of $4\pi$ total area covering the unit sphere uniformly. If instead you are looking for an average, see further below.

The full integral would be

Integrate[Abs[Det[{{Sin[θ1] Cos[φ1], Sin[θ1] Sin[φ1], Cos[θ1]},
                   {Sin[θ2] Cos[φ2], Sin[θ2] Sin[φ2], Cos[θ2]},
                   {Sin[θ3] Cos[φ3], Sin[θ3] Sin[φ3], Cos[θ3]}}]]*
  Sin[θ1] Sin[θ2] Sin[θ3],
  {θ1, 0, π}, {φ1, 0, 2 π}, {θ2, 0, π}, {φ2, 0, 2 π}, {θ3, 0, π}, {φ3, 0, 2 π}]

but doesn't seem to evaluate.

Rotational invariance means that we can set $\theta_3=\phi_3=0$ (but keep a factor of $4\pi$ for the corresponding surface integral measure), as well as $\phi_2=0$ (but keep a factor of $2\pi$ for the corresponding line integral measure):

4π*2π*Integrate[Abs[Det[{{Sin[θ1] Cos[φ1], Sin[θ1] Sin[φ1], Cos[θ1]},
                         {Sin[θ2], 0, Cos[θ2]},
                         {0, 0, 1}}]]*
  Sin[θ1] Sin[θ2],
  {θ1, 0, π}, {φ1, 0, 2 π}, {θ2, 0, π}]

8 π^4

The numerical integral agrees with this result of $8\pi^4\approx779.273$:

NIntegrate[Abs[Det[{{Sin[θ1] Cos[φ1], Sin[θ1] Sin[φ1], Cos[θ1]},
                   {Sin[θ2] Cos[φ2], Sin[θ2] Sin[φ2], Cos[θ2]},
                   {Sin[θ3] Cos[φ3], Sin[θ3] Sin[φ3], Cos[θ3]}}]]*
  Sin[θ1] Sin[θ2] Sin[θ3],
  {θ1, 0, π}, {φ1, 0, 2 π}, {θ2, 0, π}, {φ2, 0, 2 π}, {θ3, 0, π}, {φ3, 0, 2 π}]

765 ± 17

Playing with the Method option of NIntegrate will give better precision.

If you're looking for the average instead, you need to divide by the used surface measure, $(4\pi)^3$ (three times a spherical surface of $4\pi$):

8 π^4/(4 π)^3

π/8

and if you are looking further for the average tetrahedral volume instead of the parallelepiped, then as @MikeY says you further divide by 6:

8 π^4/(4 π)^3/6

π/48

This result agrees with @MikeY's numerical answer.

Not a symbolic answer, but gives a target to shoot for.

You want the average volume of a tetrahedron with edges defined by points randomly picked from the sphere. So first the routine to randomly pick from the sphere, define a multivariate normal with mean {0,0,0} and variance = IdentityMatrix[3]. Then sample randomly from this distribution, which is evenly dispersed in all directions, and normalize so they are of length = 1.

mnd = MultinormalDistribution[{0, 0, 0}, IdentityMatrix[3]];
spherePoint := (sp = RandomVariate[mnd])/Norm[sp];

Now just compute the expected value of your function.

n = 100000; 
res = 1/n Sum[Abs@Det[{spherePoint, spherePoint, spherePoint}]/3!, {n}]

0.065389

By comparison, with n=100 the estimate is .0643, so it converges rapidly.

No obvious connection to $\pi$.