Question on covariance matrix of 2 spatial signals
If you have two signal vectors $x_1[n]$ and $x_2[n]$ each of $N$ elements, then there are two different things we can consider.


How do the quantities $\displaystyle\sum_{n=1}^N x_i[n]x_j[n],~ i, j \in \{1,2\}$ compare?  In particular, when the signals are noisy and the noises can be considered 
to be jointly stationary (or jointly wide-sense stationary), these quantities can be used to estimate the noise variances in the two signals as well as the covariance of the noises at any fixed sampling time.  This is what you get from the $2\times 2$
covariance matrix
$$R_{2\times 2} = \left[\begin{matrix} \sigma_1^2 &amp; C\\
C &amp; \sigma_2^2\end{matrix}\right].$$
The noise in $x_1[n]$ has variance 
$\sigma_1^2 = R_{1,1}$ which might be different
from $R_{2,2} = \sigma_2^2$, the variance of the noise in $x_2[n]$.  But the noises are
correlated with covariance $R_{1.2}=R_{2,1} = C$.  Now, if we plan on doing things
with just what happens at $n$, ignoring whatever might be happening at $n-1$ or $n+1$
etc., then this is all the information we need.
Unless the noise is known to be (or assumed to be) white noise so that noise
samples from different sampling instants are independent (and therefore uncorrelated)
or we simply assume uncorrelated noise samples, there is information that we 
are ignoring by not considering the correlation between $x_1[n]$ and $x_1[m]$,
samples from the same process at different times or locations, and the correlation between $x_1[n]$ and $x_2[m]$, samples from the two processes at different times
or locations.  This additional information might lead to a better estimate/solution.
Now  we have a total of $2N$ noise samples,
and therefore a $2N \times 2N$ covariance matrix to consider.  If we 
arrange matters the way the authors did, we have
$R_{\text{full}} = E[XX^T]$ where 
$$X = (x_1[1],x_1[2],\ldots,x_1[N],x_2[1],x_2[2],\ldots,x_2[N])^T 
= (\mathbf x_1,\mathbf x_2)^T$$
and so
$$R_{\text{full}} = \left[
\begin{matrix}R_{\mathbf x_1,\mathbf x_1} &amp; R_{\mathbf x_1,\mathbf x_2}\\
R_{\mathbf x_2,\mathbf x_1} &amp; R_{\mathbf x_2,\mathbf x_2}\end{matrix}
\right]$$ where $R_{\mathbf x_i,\mathbf x_j} = E[\mathbf x_i\mathbf x_j^T]$.
Note that $R_{\mathbf x_i,\mathbf x_j}$ is, in essence, the 
cross-correlation function of $(x_i[1],x_i[2],\ldots,x_i[N])$
and $(x_j[1],x_j[2],\ldots,x_j[N])$ if $i \neq j$ and the 
autocorrelation function if $i = j$.  If the noise processes are white and
uncorrelated except when $n = m$, then 
$$R_{\text{full}} \rightarrow R_{\text{simple}} = \left[ \begin{matrix}\sigma_1^2I &amp; CI\\
CI &amp; \sigma_2^2I \end{matrix}\right]$$
where $I$ is the $N\times N$ identity matrix, and 
$\sigma_1^2, \sigma_2^2$ and $C$ are as defined in Item 1 above.
How realistic this noise model might be is something for the end user
to determine. If the model is realistic, then nothing is gained
by looking at the $2N\times 2N$ matrix $R_{\text{full}}$
since all the information is there
in the $2\times 2$ matrix $R_{2\times 2}$ of Item 1 above.  Ditto if the model is
unrealistic but we don't intend to (or are unable to)
use all the information in the full 
$2N\times 2N$ matrix $R_{\text{full}}$; we will make do
with just $\sigma_1^2, \sigma_2^2$ and $C$ of Part 1
for which we don't need $R_{\text{full}}$ or $R_{\text{simple}}$,
just $R_{2\times 2}$.