The average here is over different realizations of random process $\eta(\tau)$. If we measured the values of noise only at discrete time instants $\tau_1, \tau_2, .., \tau_n$, we could write a joint probability of these as
$$
w(x_1, x_2, ..., x_n),
$$
which in this case should be a multivariate Gaussian distribution with diagonal covariance matrix. As we take time interval to be smaller and smaller and pass to a continuous limit, we have to average over a functional
$$
w[x(t)]
$$
and use functional calculus. A good book covering this issues is the first volume by KLyatskin. (This is unfortunately a translation from an old Russian text. If somebody can recommend an equivalent text taht is mroe readily available, I would be glad toa dd it here.)
In the equation at hand the averaging is trivial, since the correlation function is known:
$$
\int_0^t d\tau\int_0^{t'} d\tau'\langle \eta_i(\tau)\eta_j(\tau')\rangle e^{-\gamma(t+t'-\tau-\tau')} =
\int_0^t d\tau\int_0^{t'} d\tau'\delta_{i,j}\delta(\tau-\tau') e^{-\gamma(t+t'-\tau-\tau')}=\\
\delta_{i,j}\theta(t'-t)\int_0^t d\tau e^{-\gamma(t+t'-2\tau)} +
\delta_{i,j}\theta(t-t')\int_0^{t'} d\tau' e^{-\gamma(t+t'-2\tau')},
$$
where $\theta(t)$ is the Heaviside step function, which accounts for the overlap of integration ranges necessary for delta-function taking a non-zero value.
