You can look upon the field equations in different ways, i guess.
What Einstein did, is to find a generalization of Poisson's equation for the gravitational potential, written in a covariant way. The energy density is a component of the energy momentum tensor, so that's one reason that you look for a geometric tensor with 2 indices. The Riemann-tensor has 4 indices, so there has to be some contraction there. Due to the antisymmetries of the Riemann-tensor, the only sensible contraction leads to the Ricci-tensor.
If you want to go from the special to the general case, the most simple prescription is to replace partial derivatives with covariant derivatives. The energy-momentum tensor is "covariantly conserved", so the geometric tensor representing this energy distribution als should be. The Einstein tensor is the construction from the Riemanntensor which has this property:
<br />
\nabla_{\mu}G^{\mu\nu} = 0<br />
These are the contracted Bianchi identities. As far as I know, they can be seen as a consequence of the abundant degrees of freedom you have to choose your coordinates. In this sense general relativity can be looked upon as a gauge theory, with the symmetry group being general coordinate transformations.
You can also approach it via the variational principle. The most simple Lagrangian which gives non-trivial equations of motion with up to second order derivatives of the metric is the Ricci-scalar. This resulting action is called the Hilbert action, which only depends on the metric and derivatives of it. So this could give you the vacuum equations. To get the full field equations, you couple matter to gravity. This matter field is described by a Lagrangian which depends on the metric and the matter field. So we can write the total action as
<br />
S[g_{\mu\nu},\psi ] = S_{H}[g_{\mu\nu}] + S_{M}[g_{\mu\nu},\psi]<br />
I use square brackets to indice that these actions are functionals of the fields. The total variation becomes then
<br />
\delta S = \int ( \frac{\delta S_{H}}{\delta g_{\mu\nu}}\delta g_{\mu\nu} + \frac{\delta S_{M}}{\delta g_{\mu\nu}}\delta g_{\mu\nu} + \frac{\delta S_{M}}{\delta\psi}\delta\psi + \partial_{\mu}B^{\mu} )d^{4}x<br />
The B indicates the boundary terms, and the derivatives in the integrand are the functional derivatives. Because we are interested in the symmetries under general coordinate transformations, we impose the equations of motion for the matter field with putting the boundary terms B to zero. The variations are now looked upon as general coordinate transformations, which are induced by a vector field X. So we can use the Lie-derivative for this variation. The variation for a general covariant action S under a coordinate transformation can now be written as
<br />
\delta_{X}S = \int \frac{\delta S}{\delta g_{\mu\nu}}\delta g_{\mu\nu} = 0<br />
This is because the equations of motion for the matter field makes that term vanish. This variation should be zero, because the action is invariant under general coordinate transformations. If we now plug in the explicit Lie-derivative of the metric,
<br />
\delta_{X}g_{\mu\nu} = \nabla_{\mu}X_{\nu} + \nabla_{\nu}X_{\mu}<br />
and impose boundary conditions on X and use a partial integration we obtain
<br />
\nabla_{\mu}\frac{\delta S}{\delta g_{\mu\nu}} =0<br />
These aren't equations of motion ofcourse ( we used them already to derive this result ), but so called gauge identities. If S is the Hilbert action, we obtain again the contracted Bianchi-identities. These are in this way a consequence of the invariance under general coordinate transformations. If S is the matter action, we obtain that the energy momentum tensor describing that matter field is covariantly conserved. This makes the field equations consistent. Note that this freedom enables us also to add an extra metric term to the field equations, due to metric compatibility. This leads to the famous cosmological constant.
So the equations
<br />
\nabla_{\mu}T^{\mu\nu}=0<br />
can be seen as the generalization of the special relativistic case, or as a consequence of your demand of a theory of gravity being invariant under general coordinate transformations ( which can be described by diffeomorphisms )
Hope this helps a little.
