With Domenico's clear explanation, I can actually write down more or less explicitly the DGLA describing the deformations of Einstein metrics.
First, some notation. Let $\bar{g}_{ab}$ denote a given (background) Einstein metric, with corresponding Levi-Civita connection $\bar{\nabla}_a$, Riemann tensor $\bar{R}_{abc}{}^d$, and Ricci tensor $\bar{R}_{ac} = \bar{R}_{abc}{}^b$. It is given that $\bar{\nabla}_a \bar{g}_{bc} = 0$ and $\bar{R}_{ac} = k \bar{g}_{ac}$ for a fixed constant $k$.
Let $\nabla_a$ denote an arbitrary symmetric (torsion free) affine connection, which differs from the background Levi-Civita connection as $\nabla_a v^b = \bar{\nabla}_a v^b + C^b_{ac} v^c$. Let me call $C^b_{ac} = C^b_{(ac)}$ the corresponding Christoffel tensor (or connection coefficients). The equation for an Einstein metric $g_{ab} = \bar{g}_{ab} + h_{ab}$ (with the same constant $k$) can be written in the form
\begin{align*}
\nabla_a g_{bc}
&= \bar{\nabla}_a h_{bc}
- C^d_{ab} \bar{g}_{dc} - C^d_{ac} \bar{g}_{bd}
- C^d_{ab} h_{dc} - C^d_{ac} h_{bd} ,
\\
R_{ac}[C] - k h_{ac}
&= -kh_{ac} - \bar{\nabla}_a C^b_{bc} + \bar{\nabla}_b C^b_{ac}
+ C^b_{ac} C^d_{db} - C^b_{ad} C^d_{cb} ,
\end{align*}
where $R_{ac}[C] = R_{abc}{}^b[C]$ is the usual Ricci contraction of the curvature tensor $R_{abc}{}^d[C]$ of $\nabla_a$. The first of the above equations is the $g$-compatibility condition for $\nabla_a$. Solving it for the Christoffel tensors identifies $\nabla_a$ with the Levi-Civita connection of $g_{ab}$ and plugging that solution into the second equation gives the equation $R_{abc}[g] - k g_{ac} = 0$ for an Einstein metric, due to the identity $R_{ab}[g] = \bar{R}_{ab} + R_{ab}[C]$.
Now the equations are precisely in the form that Domenico described, $f[h,C] + Q[h,C;h,C] = 0$, with $f$ linear and $Q$ quadratic
\begin{align*}
\begin{pmatrix}
f_{abc}[h,C] \\ f_{ac}[h,C]
\end{pmatrix}
&= \begin{pmatrix}
\bar{\nabla}_a h_{bc}
- C^d_{ab} \bar{g}_{dc} - C^d_{ac} \bar{g}_{bd} \\
-kh_{ac} - \bar{\nabla}_a C^b_{bc} + \bar{\nabla}_b C^b_{ac}
\end{pmatrix}
\\
\begin{pmatrix}
Q_{abc}[h,C;h',C'] \\ Q_{ac}[h,C;h',C']
\end{pmatrix}
&= \frac{1}{2} \begin{pmatrix}
- C^d_{ab} h'_{dc} - C^d_{ac} h'_{bd}
- C'^d_{ab} h_{dc} - C'^d_{ac} h_{bd} \\
+ C^b_{ac} C'^d_{db} - C^b_{ad} C'^d_{cb}
+ C'^b_{ac} C^d_{db} - C'^b_{ad} C^d_{cb}
\end{pmatrix}
\end{align*}
It remains to describe how infinitesimal symmetries (diffeomorphisms) act on the $h_{ab}$ and $C^b_{ac}$ tensor fields, as well as on $f$ and $Q$. They are generated by vector fields $u^a$. They act on tensors via the usual Lie derivative, which we find convenient to express via $\bar{\nabla}_a$, so that $\mathcal{L}_u X^a = u^b \bar{\nabla}_b X^a - X^b \bar{\nabla}_b u^a$ and $\mathcal{L}_u Y_a = u^b \bar{\nabla}_b Y_a + Y_b \bar{\nabla}_a u^b$, and they act on each other via the usual Lie bracket $[u,u'] = \mathcal{L}_u u' = -\mathcal{L}_{u'} u$. But special attention must be paid to the identities
\begin{align*}
\mathcal{L}_u \bar{\nabla}_b X^a - \bar{\nabla}_b \mathcal{L}_u X^a
& = u^c \bar{\nabla}_c \bar{\nabla}_b X^a
- (\bar{\nabla}_b X^c) \bar{\nabla}_c u^a
+ (\bar{\nabla}_c X^a) \bar{\nabla}_b u^c \\
& \quad {}
-\bar{\nabla}_b (u^c \bar{\nabla}_c X^a - X^c \bar{\nabla}_c u^a) \\
&= u^c \bar{R}_{bcd}{}^{a} X^d + X^c \bar{\nabla}_b \bar{\nabla}_c u^a
= (\bar{\nabla}_{(b} \bar{\nabla}_{c)} u^a - u^d \bar{R}_{d(bc)}{}^a) X^c , \\
% &= u^c \bar{R}_{bcd}{}^{a} X^d - X^c \bar{R}_{bcd}{}^a u^d + X^c \bar{\nabla}_c \bar{\nabla}_b u^a \\
% &= X^c \bar{R}_{cdb}{}^a u^d + X^c \bar{\nabla}_c \bar{\nabla}_b u^a , \\
\mathcal{L}_u \bar{\nabla}_b Y_a - \bar{\nabla}_b \mathcal{L}_u Y_a
& = u^c \bar{\nabla}_c \bar{\nabla}_b Y_a
+ (\bar{\nabla}_b Y_c) \bar{\nabla}_a u^c
+ (\bar{\nabla}_c Y_a) \bar{\nabla}_b u^c \\
& \quad {}
-\bar{\nabla}_b (u^c \bar{\nabla}_c Y_a + Y_c \bar{\nabla}_a u^c) \\
&= -u^c \bar{R}_{bca}{}^{d} Y_d - Y_c \bar{\nabla}_b \bar{\nabla}_a u^c
= (\bar{\nabla}_{(b} \bar{\nabla}_{a)} u^c - u^d \bar{R}_{d(ba)}{}^c) Y_c ,
\end{align*}
which identify the action of $[\mathcal{L}_u, \bar{\nabla}_b]$ as a derivation on the algebra of tensors. This means that infinitesimal diffeomorphisms generate the infinitesimal transformation $(h,C) \mapsto (h,C) + \epsilon K[u;h,C] + O(\epsilon^2)$, where
\begin{equation*}
\begin{pmatrix}
K_{ab}[u;h,C] \\ K_{ab}^c[u;h,C]
\end{pmatrix}
= \begin{pmatrix}
\bar{g}_{ac} \bar{\nabla}_b u^c + \bar{g}_{cb} \bar{\nabla}_a u^c \\
\bar{\nabla}_{(a} \bar{\nabla}_{b)} u^c - u^d \bar{R}_{d(ab)}{}^c
\end{pmatrix} .
\end{equation*}
Finally, we can put these formulas together in the definition of a DLGA $(L,d,[-,-])$. $L$ itself will break down into a sum of sections of certain tensor bundles. For simplicity, I will write $T$ to denote the space of sections of the bundle of vectors, $S^2T^*$ for symmetric covariant 2-tensors, etc. The breakdown by degree is
\begin{gather*}
\begin{array}{c|ccccc}
& 0 && 1 && 2 \\
\hline
L & T &\to& S^2 T^*\oplus S^2T^*\otimes T &\to& S^2T^* \otimes T^* \oplus S^2 T^* \\
d & & K & & f &
\end{array} , \\
\begin{array}{c|ccc}
[-,-] & 0 & 1 & 2 \\
\hline
0 & [u,u'] & K[u;h',C'] & \mathcal{L}_u \\
1 & -K[u';h,C] & 2Q[h,C;h',C'] & 0 \\
2 & -\mathcal{L}_{u'} & 0 & 0
\end{array}
\end{gather*}
I believe that this DGLA could be extended by one more degree to take the Bianchi identities into account. But I will stop here.
There are of course other ways to present the same DGLA and one can find explicit attempts in the literature of writing it down. Here's one that uses a somewhat different presentation:
Michael Reiterer, Eugene Trubowitz The graded Lie algebra of general relativity arXiv:1412.5561
