What is a covariance matrix?
It's the key point of array signal processing, I suppose. Say $x$ is the input vector of $[N,1]$ dimension collected from $N$ array sensors. $x(k)$ is its realization at the $k$ moment of time. By its definition covariance matrix (sometimes it's called autocorrelation matrix):


$R = E[x\cdot x^H]$ ,


where $E[]$ is expectation operator and $x^H$ is Hermitian conjugate. For the ergodic process


$R = \lim_{M\to\infty} 1/M \cdot \displaystyle\sum_{k=0}^M x(k) \cdot [x(k)]^H$ .


But in practice we can estimate $R$ with necessary precision with the snapshot of finite length. In the array processing theory snapshot is a group of vectors $x(k)$. It's the basic data block for array processing algorithms.


$R = 1/K \cdot \displaystyle\sum_{k=0}^K x(k) \cdot [x(k)]^H$ ,


where $K$ is the number of spatial vectors or snapshot size. The interesting question is what is optimal value for $K$. In models I've done, $K$ varies from 64 to 256. The estimation precision is enough for MVDR beamformer or Capon spectral estimation. There is interesting trick for adaptive beamformer design (if you have further interest). You can reduce $K$ (and accelerate adaptation process) if you use so-called diagonal loading, look:


$R = R + \sigma^2 \cdot I$,


where $\sigma$ is some variable depending on SNR (maybe someone define it more precisely?) and $I$ is identity matrix of the size $[N, N]$, same size as $R$ actually. But if you want to do spectral estimation procedure by your array (DoA estimation), you shouldn't perform diagonal loading because it will rise noise floor of the estimation and some weak signals of interest will be lost. 

Covariance matrix is the second order statistic of the random process which is measured at the array sensors. It contains information about the sources in space (number, strength, direction) and can be used for sources detection and separation. Actually the number of independent spatial signals at the input of array and the rank of $R$ is the same. Singular values decomposition (SVD) of $R$ gives the information about signal subspace which is necessary for subspace-based DoA estimation techniques like MUSIC or ESPRIT. It is worth to say algorithms based upon covariance matrix inversion/decomposition suffers from source correlation. If the sources is highly correlated (have the same waveforms or too close direction) they can't be separated. And also the system performance degrades in highly correlated scenario.

Very good reference about the topic will be


  Harry L. Van Trees - Detection, Estimation, and Modulation Theory, Optimum Array Processing


In the book more general form of covariance matrix is discussed. It's shown that it can be defined either in time domain or in frequency domain. And furthermore frequency domain interpretation is more common since it allows us to solve the problem of wideband beamformer (or DoA estimation) very easy. 

Hope this helps.