Short Answer:
In this example you would define your differential spatial volume as:
dV = \sqrt{\mathop{det}\nolimits_3(g)}dr d\theta d\varphi
where det_3(g) is the determinant of the spatial metric when you set t=0 (defining your manifold) and dt=0 (since you want submanifold of constant t value) in matrix form:
g = \left[\begin{array}{ccc}\frac{1}{1-kr^2} & 0 & 0\\ 0 & r^2 & 0\\ 0 & 0 & r^2 \sin^2(\theta)\end{array}\right]
This depends on the fact that the spatial coordinates are orthogonal to your time coordinate and the time coordinate is basically normalized (up to unit scaling via c). Note that for space-time hyper-volumes we use \sqrt{-\mathop{det}\nolimits_4(g)} due to the signature yielding a negative determinant. This determinant is not quite a "determinant" but rather a "pseudo-determinant". It is not the determinant of an operator defined by the product of its eigen-values. The metric is a mapping from vector to co-vector which then can map a vector to a number. I'm not 100% sure how to describe this metric determinant quantity other than how it manifests below...
Very Long Answer:
The root of your problem is how to get the area of a hyper-surface in a curved space (manifold) from the metric structure. Critical to this is the idea of differential forms and the extension of the metric onto these.
Picture a volume function, called the determinant that maps, in this case four 4-vectors to a 4-volume. For a curved manifold it must map four differential vectors to a differential 4-volume d^4\mathbf{x} = det[d\mathbf{x}_1,d\mathbf{x}_2,d\mathbf{x}_3,d\mathbf{x}_4].
You also have the "dot" product defined by the metric such that ds^2 = d\mathbf{x}\bullet d\mathbf{x}. Once you have a 4-volume and a 1-volume (length) you can subtractively define a 3-volume and then a 2-volume how orthogonal (as defined by our metric) forms multiply to give higher measure forms. In this case we can define a spatial 3-volume form by how it combines with a unit, orthogonal time-like vector form to yield a space-time 4-volume form, i.e. our determinant.
dx\wedge dy\wedge dz :c dt \mapsto det[\hat{e}_1 dx, \hat{e}_2 dy, \hat{e}_3 dz,\hat{e}_4 c dt]
So, getting back to your problem. You notice that the metric for your chosen coordinate has dt coefficient independent of the other variables. This means that in this metric the coordinate t direction is orthogonal to the other coordinates, and that its differential length is a constant. This means that we can treat the t=0 spatial submanifold as a 3-manifold with the intrinsic metric obtained by setting dt=0 and t=0:
ds^2 =a^2\left[ \frac{dr^2}{1-r^2} + r^2(d\theta^2 +\sin^2(θ) d\varphi^2)\right]
Ok so far? The big question then is at this point how do you get that differential 3-volume and integrate over this 3-manifold? (similar to the question you might have asked of how to calculate the 4-volume of the space-time region from t=0 to t=42seconds or some value...)
This gets a little complicated conceptually unless you embed your manifold in some larger vector space but let's give it a go... Given you are using some coordinate system on the manifold defining a point as a continuous function of your coordinates p=p(r,\theta,\phi) there is a differential of this point valued function dp= \partial_r p dr + \partial_\theta p d\theta + \partial_\varphi p d\varphi. This differential lives in the tangent space which is where your metric is defined ds^2 = dp\bullet dp. So these "partial derivatives of position" are mappings from the point on the manifold into this tangent space.
To define differential surface and volume elements we construct differential forms from the components you can think of as "partial differentials"
d_r p = \partial_r p dr,\quad d_\theta p = \partial_\theta p d\theta, \quad d_\varphi p =\partial_\varphi p d\varphi
d\mathbf{V} = d_r p \wedge d_\theta p \wedge d_\varphi p
which is an oriented volume form with the wedge products here being totally antisymmetric tensor products we use to express partial determinants in that original determinant function. To get a simple scalar volume you need the magnitude of this thing which is where we invoke the good ole metric again.
The way the metric extends to differential forms is, if we follow the convention that the wedge of two unit vectors defines a unit area etc:
(a\wedge b \wedge c)\bullet (a\wedge b \wedge c) = <br />
(a\bullet a)(b\bullet b)(c\bullet c) + (a\bullet b)(b\bullet c)(c\bullet a) + (a\bullet c)(b\bullet a)(c\bullet b)...
- (a\bullet a)(b\bullet c)(c\bullet b)-(a\bullet c)(b\bullet b)(c\bullet a) - (a\bullet b)(b\bullet a)(c\bullet c)
Since the wedge product is totally anti-symmetric and each factor vector may have components in the direction of any or all of the other factors we must sum over all the permutations (with minus sign for odd permutations due to antisymmetry).
That's the mess for general coordinates but if you work with orthogonal coordinates and bases (but not necessarily orthonormal) i.e. with coordinates with a diagonalized metric, then only the first term will be non-zero. But in general this quantity will end up being, for your differential volume element:
d\mathbf{V}\bullet d\mathbf{V} = \mathop{det}(g) dr^2d\theta^2 d\varphi^2
where the determinant of g is the matrix determinant of the 3-metric. So we take the square root and get:
dV = \sqrt{\mathop{det}(g)}dr d\theta d\varphi
Note this is what we do in a much simplified form in Calc III, when we define the "vector" and "scalar" differental surface elements:
d\mathbf{S} = d_u \mathbf{r}\times d_v\mathbf{r} =\left(\frac{\partial \mathbf{r}}{\partial u}\times\frac{\partial \mathbf{r}}{\partial v}\right) dudv, <br />
\quad dS =<br />
|d\mathbf{S}|
