Let $f:{\mathbb R}_+\mapsto {\mathbb R}$ be a smooth function with $f(0)=0.$ A problem of estimation of a functional $τ_f(Σ):= {\rm tr}(f(Σ))$ of unknown covariance operator $Σ$ in a separable Hilbert space ${\mathbb H}$ based on i.i.d. mean zero Gaussian observations $X_1,\dots, X_n$ with values in ${\mathbb H}$ and covariance operator $Σ$ is studied. Let $\hat Σ_n$ be the sample covariance operator based on observations $X_1,\dots, X_n.$ Estimators \begin{align*} T_{f,m}(X_1,\dots, X_n):= \sum_{j=1}^m C_j τ_f(\hat Σ_{n_j}) \end{align*} based on linear aggregation of several plug-in estimators $τ_f(\hat Σ_{n_j}),$ where the sample sizes $n/c\leq n_1<\dots<n_m\leq n$ and coefficients $C_1,\dots, C_n$ are chosen to reduce the bias, are considered. The complexity of the problem is characterized by the effective rank ${\bf r}(Σ):= \frac{{\rm tr}(Σ)}{\|Σ\|}$ of covariance operator $Σ.$ It is shown that, if $f\in C^{m+1}({\mathbb R}_+)$ for some $m\geq 2,$ $\|f''\|_{L_{\infty}}\lesssim 1,$ $\|f^{(m+1)}\|_{L_{\infty}}\lesssim 1,$ $\|Σ\|\lesssim 1$ and ${\bf r}(Σ)\lesssim n,$ then \begin{align*} & \|\hat T_{f,m}(X_1,\dots, X_n)-τ_f(Σ)\|_{L_2} \lesssim_m \frac{\|Σf'(Σ)\|_2}{\sqrt{n}} + \frac{{\bf r}(Σ)}{n}+ {\bf r}(Σ)\Bigl(\sqrt{\frac{{\bf r}(Σ)}{n}}\Bigr)^{m+1}. \end{align*} Similar bounds have been proved for the $L_{p}$-errors and some other Orlicz norm errors of estimator $\hat T_{f,m}(X_1,\dots, X_n).$ The optimality of these error rates, other estimators for which asymptotic efficiency is achieved and uniform bounds over classes of smooth test functions $f$ are also discussed.