The factors $c^2$ should be no problem, so I'll focus on the question how we can tell Mathematica that all the functions are small. It boils down to introducing a bookkeeping parameter $\epsilon$ that is eventually set to $1$ after it has served its purpose of keeping track of powers of the small quantities:
replacementRule = {U :> (ϵ u[##] &),
K :> (ϵ k[##] &)};
matr = ReleaseHold[
Hold[
{{(1 - 2*U[t, x, y, r]/c^2 - 2*D[K[t, x, y, r], t]/c), -D[
K[t, x, y, r], x], -D[K[t, x, y, r], y], -D[K[t, x, y, r],
r]}, {-D[K[t, x, y, r], x], -(1 + 2*U[t, x, y, r]/c^2), 0,
0}, {-D[K[t, x, y, r], y], 0, -(1 + U[t, x, y, r]/c^2),
0}, {-D[K[t, x, y, r], r], 0, 0, -(1 + 2*U[t, x, y, r]/c^2)}}
] /. replacementRule
];
Normal[Series[Inverse[matr], {ϵ, 0, 1}]] /. ϵ ->
1 // MatrixForm
the replacementRule takes your function names and replaces them with new functions (of the same number of arguments (##) using lower case names. But these new functions are multiplied by the parameter ε.
To do the replacement, I take your initial definition and wrap it in Hold so the derivatives aren't done until the renaming has occurred, at which time ReleaseHold allows the evaluation. This is something you could of course also do by hand, but maybe not if your real application is larger.
Then the Series expansion requires only a single expansion, in terms of ε. You have to use Normal to convert the SeriesData object to a "normal" expression before finally setting ε to 1.
This is closely related to the question Multivariable Taylor expansion does not work as expectedMultivariable Taylor expansion does not work as expected, except that here we have to work with functions.
