The Hilbert space of a two-qubit system is $4$-dimensional complex vector space. An arbitrary normalized vector in this space can be written as:
$$|\psi\rangle= \frac{[w_0 \, w_1 \, w_2 \, w_3]^t}{\sqrt{|w_0|^2+|w_1|^2+|w_2|^2+|w_3|^2}}$$
Since an overall phase doesn't change the state of the system, we may choose the first component $w_0$ as real, thus if we define (for $w_0\ne 0$):
$$z_i = \frac{w_i}{w_0}$$
We obtain:
$$|\psi\rangle= \frac{[1 \, z_1 \, z_2 \, z_3]^t}{\sqrt{1+|z_1|^2+|z_2|^2+|z_3|^2}}$$
The coordinate functions $z_1$, $z_2$, $z_3$ parametrize almost everywhere the complex projective space $\mathbb{C}P^3$ which is the state space of two-qubit systems; but this fact will not be needed as the Haar probability measure on the space will be derived below from scratch.
The Euclidean measure on the complex four-dimensional vector space $\mathbb{C}^4$ is given by:
$$d_{\mu_{\mathbb{C}^4}} = \prod_{k=1}^4 \frac{d\text{Re}(w_k)d\text{Im}(w_k)}{2\pi}$$
Obviously, this measure is invariant under the $4$-dimensional unitary group. The normalization condition defines a seven-dimensional sphere $S^7$, a measure on $S^7$ can be constructed as:
$$d_{\mu_{S^7 }} = \int \delta(|w_0|^2+|w_1|^2+|w_2|^2+|w_3|^2-1) \prod_{k=1}^4 \prod_{k=1}^4 \frac{d\text{Re}(w_k)d\text{Im}(w_k)}{\pi}$$
($\delta$ is the Dirac delta function which restrict the measure to a unit spherical shell). Obviously, this measure is also invariant under the $4$-dimensional unitary group, thus it is a Haar-measure.
Substituting the equations for $z_i$ and integrating over $w_0$, we obtain a Haar-measure over $\mathbb{C}P^3$, i.e., the two-qubit state space:
$$ d_{\mu_{\mathbb{C}P^3}} = \int \delta(|w_0|^2(1+|z_1|^2+|z_2|^2+|z_3|^2)-1) \prod_{k=1}^3\frac{d\text{Re}(w_k)d\text{Im}(w_k)}{\pi} = \int |w_0|^6 (1+|z_1|^2+|z_2|^2+|z_3|^2)^{-1} \delta\left(|w_0|^2- \frac{1}{1+|z_1|^2+|z_2|^2+|z_3|^2}\right) \frac{d\text{Re}(w_0d\text{Im}(w_0)}{\pi}\prod_{k=1}^3 \frac{d\text{Re}(z_k)d\text{Im}(z_k)}{\pi}$$
Performing the $w_0$ integration in polar coordinates:
We have:
$$\int |w_0|^6 \delta\left(|w_0|^2- \frac{1}{1+|z_1|^2+|z_2|^2+|z_3|^2}\right)\frac{d\text{Re}(w_0)d\text{Im}(w_0)}{\pi}\ = (1+|z_1|^2+|z_2|^2+|z_3|^2)^{-3}$$
Thus:
$$d_{\mu_{\mathbb{C}P^3}} = \frac{1}{(1+|z_1|^2+|z_2|^2+|z_3|^2)^4}\prod_{k=1}^3\frac{d\text{Re}(z_k)d\text{Im}(z_k)}{\pi}$$
In summary, we have a representation of a random two-qubit state vector together with a Haar probability measure on the state space.
Of course, repeating this construction for a single qubit, we get the usual Bloch vector and the round measure on the Bloch sphere:
$$|\psi\rangle= \frac{[1 \, z]^t}{\sqrt{1+|z|^2}}$$
$$d_{\mu_{S^2}} = \frac{1}{(1+|z|^2)^2}\frac{d\text{Re}(z)d\text{Im}(z)}{\pi}$$
The Cartesian product of two Bloch spheres is the state space of a random separable two-qubit state, the state vector is given by:
$$|\psi\rangle= \frac{[1 \, z_1]^t \otimes [1 \, z_2]^t }{\sqrt{(1+|z_1|^2)( 1+|z_2|^2)}}$$
And the corresponding measure:
$$d_{\mu_{S^2\times S^2}} = \frac{1}{(1+|z_1|^2)^2(1+|z_2|^2)^2}\prod_{k=1}^2\frac{d\text{Re}(z_k)d\text{Im}(z_k)}{\pi}$$
Practically, integrals of polynomial functions over the Haar measures $d_{\mu_{\mathbb{C}P^3}}$ $d_{\mu_{S^2\times S^2}}$ can be exactly evaluated e.g. by passing to polar coordinates.
but the jump of sampling from one to two Bloch Spheres does not seem obvious or simple to me.
