Just test the invariance of the action (so using integrations by part and neglecting surface terms) under "linear diffeomorphism" $h_{\mu\nu} \to h_{\mu\nu} + \partial_\mu \epsilon_\nu+\partial_\nu \epsilon_\mu$
[EDIT]
More precisions :
The linear diffeomorphism is only the linearization of the standard diffeomorphism :
$g^{'\mu\nu} = \frac{\partial x^{'\mu}}{\partial x^{\sigma}} \frac{\partial x^{'\nu}}{\partial x^{\tau}} g^{\sigma\tau}$
for a coordinate transformation $x^{'\mu} = x^\mu + \epsilon^\mu(x)$
where you plug $g^{\sigma\tau} = \eta^{\sigma\tau} + h^{\sigma\tau}$, and $\frac{\partial x^{'\mu}}{\partial x^{\sigma}} = \delta^\mu_\sigma +\partial_\sigma \epsilon^\mu$, and where you keep only terms linear in $\epsilon$.
Due to Lorentz invariance, you have only $4$ possible quadratic terms in first derivative of $h_{\mu\nu}$ , and we are searching the correct combination which leaves the action invariant under the linear diffeomorphism (which is nothing but a gauge invariance for spin $2$ field). So you may write :
$S = \int d^4x (A ~\partial_\alpha h_{\nu\beta}~\partial^\alpha h^{\nu\beta} + B ~\partial_\alpha h_\nu^\nu~\partial^\alpha h_\nu^\nu +C ~\partial_\alpha h^{\alpha\nu}~\partial^\beta h_{\beta\nu}+D ~\partial^\alpha h_\nu^\nu~\partial^\beta h_{\beta\alpha}) $
We may fix a value for $A$, here $A = \frac{-1}{32G} (\frac{1}{2})$, so we have only $3$ variables to be checked.
Now, the variation $\delta S$ of the action is going to give, after integration by parts, terms in $\epsilon \partial^3h$, that is terms linear in $\epsilon$ and $h$, with a third derivative in $h$, and there are only $3$ possible terms :
$\epsilon^\nu(\partial^2\partial^\mu h_{\mu\nu}), \epsilon^\nu (\partial_\nu \partial^2 h_\mu^\mu), \epsilon^\nu(\partial_\nu \partial^\mu \partial^\lambda h_{\mu\lambda})$
The coefficients of these terms must be zero, so you have $3$ equations for $3$ variables $B, C, D$, so you are able to calculate their values.
Integration by parts works like this, take for instance the first term, you have :
$\delta_1 S = A~ \delta(\partial_\alpha h_{\nu\beta}~\partial^\alpha h^{\nu\beta}) = 2A~( \partial_\alpha(2\partial_\nu \epsilon_\beta)\partial^\alpha h^{\nu\beta})$
Now, by applying $2$ times the integration by parts, and neglecting the surface terms (total derivatives), you get :
$\delta_1 S=4\epsilon_\beta \partial^2 \partial_\nu h^{\nu\beta}$
The result of your calculus should be : $B=-A, C=-2A, D=2A$ , and you may check that this corresponds to your original expression (developped in $h$ terms only, so by replacing $\bar h$ terms by their value)