We consider a general class of symmetric or Hermitian random band matrices $H=(h_{xy})_{x,y \in \llbracket 1,N\rrbracket^d}$ in any dimension $d\ge 1$, where the entries are independent, centered random variables with variances $s_{xy}=\mathbb E|h_{xy}|^2$. We assume that $s_{xy}$ vanishes if $|x-y|$ exceeds the band width $W$, and we are interested in the mesoscopic scale with $1\ll W\ll N$. Define the {\it{generalized resolvent}} of $H$ as $G(H,Z):=(H - Z)^{-1}$, where $Z$ is a deterministic diagonal matrix with entries $Z_{xx}\in \mathbb C_+$ for all $x$. Then we establish a precise high-probability bound on certain averages of polynomials of the resolvent entries. As an application of this fluctuation averaging result, we give a self-contained proof for the delocalization of random band matrices in dimensions $d\ge 2$. More precisely, for any fixed $d\ge 2$, we prove that the bulk eigenvectors of $H$ are delocalized in certain averaged sense if $N\le W^{1+\frac{d}{2}}$. This improves the corresponding results in \cite{HeMa2018} under the assumption $N\ll W^{1+\frac{d}{d+1}}$, and in \cite{ErdKno2013,ErdKno2011} under the assumption $N\ll W^{1+\frac{d}{6}}$. For 1D random band matrices, our fluctuation averaging result was used in \cite{PartII,PartI} to prove the delocalization conjecture and bulk universality for random band matrices with $N\ll W^{4/3}$.