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 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 He and Marcozzi (Diffusion profile for random band matrices: a short proof, 2018. arXiv:1804.09446 ) that imposed the assumption $$N\ll W^{1+\frac{d}{d+1}}$$ , and the results in Erdős and Knowles (Ann Henri Poincare12(7):1227–1319, 2011; Commun Math Phys 303(2): 509–554, 2011) that imposed the assumption $$N\ll W^{1+\frac{d}{6}}$$ . For 1D random band matrices, our fluctuation averaging result was used in Bourgade et al. (J Stat Phys 174:1189–1221, 2019; Random band matrices in the delocalized phase, I: quantum unique ergodicity and universality, 2018. arXiv:1807.01559 ) to prove the delocalization conjecture and bulk universality for random band matrices with $$N\ll W^{4/3}$$ .