Gravitational waves energy source in linearized theory By linearizing the metric in the following way (approach in most textbooks):
$g_{\mu\nu}=\eta_{\mu\nu}+h_{\mu\nu}\text{ with } |h_{\mu\nu}|\ll 1$
and choosing the transverse-traceless gauge a wave equation for the perturbation $h_{\mu\nu}^{TT}$ is obtained (with $\Box=\eta_{\mu\nu}\partial^{\mu}\partial^{\nu}$:
$\Box h_{\mu\nu}=0$
Therefore possible solutions for the perturbations $h_{\mu\nu}$ are waves with $+$- resp. $\times$-polarization.
These waves have an effect on test masses eg. deforming a ring according to the polarization of the wave. So by introducing little springs between the test masses the gravitational wave must carry energy (stored in the springs by friction). But what is the source of this energy? 
In Gravitational Waves Vol. I by Michele Maggiore the stress-energy tensor of GWs $T_{\mu\nu}^{\text{GW}}$ is derived by linearizing the metric like this:
$g_{\mu\nu}=g_{\mu\nu}^B+h_{\mu\nu}$
Where $g_{\mu\nu}^B$ is the curved background. This makes sense to me, because GWs themselves are sources of curvature. 
I'm confused by the first approach: The GWs seem to carry energy, but in the field equations $T_{\mu\nu}$ is set to zero. So where is this energy coming from? How should the waves be interpreted physically?
EDIT:
After thinking about it again: $T_{\mu\nu}$ is set to zero because we look at a region in spacetime, where there is no matter/energy density and describe how a the perturbation $h_{\mu\nu}$ behaves (like a wave). But the perturbation must be caused somewhere by some source. So the wave is coming from a region with a source (similar to the vacuum wave solution in electrodynamics). 
I'm still a bit confused what the physical interpretation of these waves is.
 A: Without exploring all the mathematical details, the essence of the energy situation is as follows.
First we suppose there is a body such as a binary star system, where the quadrupole moment is a function of time. As the stars orbit around their mutual center of mass, the spacetime around them curves accordingly, and some of this curvature propagates outwards. Meanwhile the orbits also change gradually.
When the propagating part (called by us gravitational waves) reaches other things, such as masses on springs, they produce a tidal effect, and this causes the masses to move relative to one another so the springs are compressed or extended.
If we now look at energy, then we find that the gradual change in the orbit at the source is such that the source is losing energy, the effect on the springs is that they gain energy. If there is friction then it will be converted into heat, for example. And furthermore, in the weak field limit at least, one can connect it all up and find that energy is conserved.
Add to this the fact the GR is a local theory, and you have the overall picture of energy transport from the source to the detector via the waves. And furthermore, this can all be made precise in the weak field limit, and in the few other cases that can be analyzed in full. However, this does not change the fact that energy is a subtle concept in GR. The energy associated with gravitational waves cannot be connected directly to the wave amplitude squared at a point. However, to a high (but not perfect) precision it can be connected, for low-amplitude waves, to certain spatial averages, where the average is taken over regions of space with a width of the order of the wavelength. This method only works when there is a clear division of length scale between the wavelength and the background curvature of the space on which the waves are appearing as a small perturbation.
The energy associated with gravitational waves does not appear in $T^{ab}$ (the stress energy tensor) in the full theory. Rather it is the non-linearity of the Einstein equation which, by a remarkable mathematical conjuring act, somehow accounts for the way gravitation moves energy around. In the linearized theory we start out with $T^{ab}$ zero in vacuum, but then to do the energy accounting we extend the treatment to the next order of approximation, and this can conveniently be done by including the next (i.e. 2nd order) terms in the Einstein tensor but choosing to interpret them as part of the stress-energy tensor. In this case the net effective stress-energy tensor will be non-zero in vacuum if there are gravitational waves.
A: Physical picture
The situation you are describing are gravitational waves (small perturbations to the metric) propagating on either a flat space metric $\eta_{\mu\nu}$ or some curved geometry $g^B_{\mu\nu}$. These waves do carry momentum and energy, and so there is a contribution to the full stress tensor $T_{\mu\nu}$ from these waves. However, since these waves are small, we can treat the contribution of the stress-energy tensor of the waves on the geometry (in the lingo, this is referred to as backreaction of the GWs on the metric) in a perturbative fashion. In particular, if we are working to leading order (purely linear in the metric perturbation), it is self-consistent to ignore the effect of the stress-energy tensor of the GWs, on the geometry.
Mathematical details
Recasting the exact equations
The exact Einstein's equations (in units with $c=1$) are
\begin{equation}
G_{\mu\nu} = 8\pi G_N T_{\mu\nu}
\end{equation}
where $G_N$ is Newton's constant.
Now let's consider just the left hand side. Let's write the metric as flat space $\eta_{\mu\nu}$ plus a small perturbation $h_{\mu\nu}$
\begin{equation}
g_{\mu\nu} = \eta_{\mu\nu} + \sqrt{G_N} h_{\mu\nu}
\end{equation}
The reason for putting the $\sqrt{G_N}$ will become apparent later, but keep in mind that $G_N$ is a small number and we will use it as our expansion parameter. By the way can go through all the subsequent steps with any background metric $g^B_{\mu\nu}$ instead of $\eta_{\mu\nu}$; the content is the same, but you will need to keep track of some extra details like using the covariant derivative associated with the background metric, and defining the background stress-energy tensor.
Then we can write
\begin{equation}
G_{\mu\nu} = \sqrt{G}_N \hat{\mathcal{E}}^{(1)} h_{\mu\nu} + G_N \hat{\mathcal{E}}[h]^{(2)} h_{\mu\nu} + \cdots + G_N^{N/2} \hat{\mathcal{E}}[h]^{(N)} h_{\mu\nu}
\end{equation}
where $\hat{\mathcal{E}}[h]^{(N)}$ is a differential operator containing two derivatives, and $h$ raised to the power $N-1$. Schematically, $\hat{\mathcal{E}}[h]^{(2)} h \sim h^2 \partial^2 h + h \partial h \partial h$. Note that the leading order operator, $\mathcal{E}^{(1)}$, does not depend on $h$ at all.
In Minkowski spacetime, $T_{\mu\nu}=0$. So Einstein's equations can be written as
\begin{equation}
\sqrt{G}_N \hat{\mathcal{E}}^{(1)} h_{\mu\nu} + G_N \hat{\mathcal{E}}[h]^{(2)} h_{\mu\nu} + \cdots + G_N^{N/2} \hat{\mathcal{E}}[h]^{(N)} h_{\mu\nu}= 0
\end{equation}
Note that up to this point, we actually have not made any approximations. We have written the exact Einstein equations with $T_{\mu\nu}=0$, but just in some weird choice of variables, and without explicitly computing the operators $\hat{\mathcal{E}}^{(N)}$. (Of course you can find these in the literature if you want, although my notation may not be standard).
Iterative solution
Now we develop an iterative solution for the metric perturbation $h_{\mu\nu}$. We write
\begin{equation}
h_{\mu\nu} = h^{(1)}_{\mu\nu} + \sqrt{G_N} h^{(2)}_{\mu\nu} + ... + G_N^{N/2} h^{(N)}_{\mu\nu}
\end{equation}
We then follow the following iterative procedure:

*

*Plug $h^{(1)}_{\mu\nu}$ into Einstein's equations and work to leading order in the small parameter $G_N$.

*Plug the solution $h^{(1)}_{\mu\nu}$ into Einstein's equations and solve the result for $h^{(2)}_{\mu\nu}$.

*Repeat as many times as desired.

Leading order
At leading order $\mathcal{O}(G_N^{1/2})$, we find
\begin{equation}
\hat{\mathcal{E}^{(1)}} h^{(1)}_{\mu\nu} = 0
\end{equation}
In de Donder gauge, this equation becomes
\begin{equation}
\square h^{(1)}_{\mu\nu} = 0
\end{equation}
This is a standard wave equation, and has standard wave solutions. If you are careful about keeping track of the gauge conditions and constraints, you will end up finding there are two independent solutions, corresponding to the $+$ and $\times$ polarizations.
Next to leading order
Now that we have solved this equation, we work to the next order $\mathcal{O}(G_N)$ for $h^{(2)}_{\mu\nu}$.
\begin{equation}
\hat{\mathcal{E}^{(1)}} h^{(2)}_{\mu\nu} + \hat{\mathcal{E}^{(2)}}[h^{(1)}]h^{(1)}_{\mu\nu}  = 0
\end{equation}
We can rearrange this equation as
\begin{equation}
\hat{\mathcal{E}^{(1)}} h^{(2)}_{\mu\nu} = T^{(\rm eff)}[h^{(1)}]_{\mu\nu}
\end{equation}
where $T^{\rm eff}[h^{(1)}]_{\mu\nu}$ is an effective stress-energy tensor that depends on the first order perturbation $h^{(1)}_{\mu\nu}$. Therefore, you can think of the stress-energy associated with gravitational waves as sourcing a higher order correction to the metric itself. This leads to some very interesting phenomena, such as non-linear gravitational memory.
There is an important caveat about thinking of $T^{\rm eff}[h^{(1)}]_{\mu\nu}$ as a real stress energy tensor: it is not gauge invariant. In particular, the stress energy tensor for gravitational waves depends on the choice of coordinates. Therefore it is not truly a stress-energy tensor (it is just some part of the Einstein tensor we have moved to the right hand side). Often, the physical insight you get from thinking of this term as if it were stress-energy tensor outweighs this mathematical fact, but if you aren't careful this can lead to confusion.
Having said that, as @AndrewSteane pointed out in the comments, you can obtain a gauge-invariant non-local stress-energy tensor by averaging the effective stress-energy tensor over a region consisting of several wavelengths.
What's generating the waves?
Your final question was: how can there be waves if nothing has generated the waves?
In fact this situation is not special to gravity. Mathematically the point is that solutions to the equations of motion are only unique if you specify asymptotic boundary conditions. If you specify that the metric must be asymptotically flat, then these boundary conditions will rule out plane wave solutions, and the only way to get GWs is to have a source somewhere in the spacetime (in which case $T_{\mu\nu}$ is not zero). However, mathematically it is possible to consider more general boundary conditions where plane waves can come in from infinity and go out to infinity. Physically, you can think of this as an idealization where an observer lives in some local frame and experiences a passing gravitational wave from very far away. This exact same situation can happen in any theory described by a wave equation, such as electromagnetism.
In more detail, you may notice that the curvature tensor $R_{\mu\nu\rho\sigma}$ has 20 independent components in 4 dimensions, but Einstein's equations $G_{\mu\nu}=8\pi G_{N} T_{\mu\nu}$ are only 10 equations relating components of the curvature to the stress-energy tensor. Anyway, the point is that there are components of the curvature which are not fixed by the matter distribution. These components can be written in terms of the Weyl tensor. The fact that the Weyl tensor is not fixed by the matter distribution, is what mathematically allows for gravitational waves to propagate through the spacetime even when $T_{\mu\nu}=0$.
A: Someone else can probably answer this more authoritatively and with better technical facility, but here is my rough idea of what is going on. I think you have two separate sets of issues to deal with in this interpretation.
One set of issues has to do with the fact that linearized gravity is incomplete and inconsistent. For example, linearized gravity can't describe bound systems, so it can't describe a realistic source such as two stars orbiting each other.
The other issues arise because of the way energy is represented in GR. GR (a) doesn't have a globally conserved, scalar measure of mass-energy that applies to all spacetimes, nor (b) does it have a local measure of the energy in the gravitational field. Issue a arises from the fact that energy-momentum is a vector, and you can't compare vectors that came from different places except through parallel transport, which is path-dependent. Issue b comes from the equivalence principle. Because of these issues, you should not expect to have a local energy density, such as a term in the stress-energy tensor, that represents the energy of the gravitatioal waves.
In asymptotically flat spacetimes, we do have various conserved measures of the total mass-energy of the spacetime, such as the Bondi energy and the ADM energy. The ADM energy, for example, does include energy being radiated away to null infinity by gravitational waves.
For a simple way to convince yourself that gravitational waves really do carry energy, consider that they manifest themselves in experiments as tidal forces, exactly like any other tidal forces, such as the ones that create the earth's ocean tides. Such tidal forces are certainly capable of doing work.
