Why is the sound field intensity due to $K$ point sources given by $ I(p,\omega) = \sum_{k=1}^K \sigma_k^2(\omega) \delta(p - p_k)$?
I believe the confusion comes from the fact that equation $(*)$ is referring to a received intensity (in units of W/m^2), whereas the second equation is referring to a radiated intensity. 

The assumption behind equation $(*)$ is that the received power as a function of angle is the surface integral of the incident intensity over some aperture area $A \subset S$. Because the surface area element $dA$ has units of area, the result of the integral will be in units of power (W). 

The Dirac delta distribution $\delta(p - p_k)$ expresses that the power density is coming from exactly one direction $p_k$ and nowhere else. In this context, the Dirac delta distribution possesses the following property:
$$
\int_A \delta(p - p_k) dA' = \begin{cases} 0 &amp; p_k \notin A \\ 1 &amp; p_k \in A \end{cases},
$$
which is analogous to a "density" where all the "mass" is concentrated at one point $p_k$. 

To calculate the total received power $P$, you integrate the received intensity over the receive aperture surface $A$:
$$
P = \int_A I(p,\omega) dA'
$$

Each term of the integral will look something like this:
$$
\int_A \sigma_k^2(\omega) \delta(p - p_k) dA' = \begin{cases} 0 &amp; p_k \notin A \\ \sigma_k^2(\omega) &amp; p_k \in A \end{cases},
$$
so in other words, you are just adding up the power coming from each source that your receive aperture $A$ captures, and each source contributes power $\sigma^2_k(\omega)$. 

Implicit in this formulation is that all the sources are incoherent, i.e., the waves do not have a fixed phase relationship and thus the average power of the superposition is the sum of the contributions from the individual sources.