As a very general rule, whenever $X = (X_1, \ldots, X_B)$ are random variables with given covariances $\sigma_{ij}=\text{Cov}(X_i,X_j),$ then the covariance of any linear combination $\lambda \cdot X = \lambda_1 X_1 + \cdots + \lambda_B X_B$ is given by the matrix $\Sigma = (\sigma_{ij})$ via
$$\text{Cov}(\lambda X, \lambda X) = \lambda^\prime \Sigma \lambda.$$
The rest is just arithmetic.
In the present case $\sigma_{ij} = \rho\sigma^2$ when $i\ne j$ and otherwise $\sigma_{ii} = \sigma^2 = \left[\rho + (1-\rho)\right]\sigma^2$. That is to say, we may view $\Sigma$ as the sum of two simple matrices: one has $\rho$ in every entry and the other has values of $1-\rho$ on the diagonal and zeros elsewhere. This leads to an efficient calculation, because evidently
$$\Sigma = \sigma^2\left(\rho 1_B 1_B^\prime + (1-\rho)\mathbb{Id}_B \right)$$
where I have written "$1_B$" for the column vector with $B$ $1$'s in it and "$\mathbb{Id}_B$" for the $B$ by $B$ identity matrix. Whence, factoring out the scalars $\sigma^2$, $\rho$, and $1-\rho$ as appropriate, we obtain
$$\eqalign{
\text{Cov}(\lambda X, \lambda X) &= \lambda^\prime \sigma^2\left(\rho 1_B 1_B^\prime + (1-\rho)\mathbb{Id}_B \right)\lambda \\
&= \left(\lambda^\prime 1_B 1_B^\prime \lambda\right) \rho\sigma^2 + \left(\lambda^\prime \mathbb{Id}_B \lambda \right)(1-\rho)\sigma^2.
}$$
For the arithmetic mean, $\lambda = (1/B, 1/B, \ldots, 1/B)$ entailing $$\lambda^\prime 1_B 1_B^\prime \lambda = (\lambda^\prime 1_B)^2 = 1^2 = 1$$ and $$\lambda^\prime \mathbb{Id}_B \lambda = 1/B^2 + 1/B^2 + \cdots + 1/B^2 = 1/B,$$ QED.