This is a relatively old question which lacks a formally complete answer. On finding myself in the need of the inverse of a metric and on not being able to find a proper treatment elsewhere (on casual browsing), I have decided to put a proper formal treatment here.
Following the treatment given here, one can (super-)easily derive the inverse metric to all order of perturbation theory without using ad-hoc relations. I've arranged the following in three steps.
Step - 1: Correct Statement of Problem
The metric whose inverse we intend to determine must be written in a more formal fashion:
$$ g_{\mu\nu} = \eta_{\mu\nu} + \epsilon \ ^{(1)}h_{\mu\nu} + \frac{\epsilon^2}{2!} \ ^{(2)}h_{\mu\nu} + \cdots $$
For later convenience, we move all perturbations into $H_{\mu\nu}$:
$$g_{\mu\nu} = \eta_{\mu\nu} + H_{\mu\nu} $$
This way of stating the problem is essentially different than stated by OP in the question. I hope the notation does not need any explanations.
Step-2: And the inverse is
Let us write the inverse as: b
$$ g^{\mu\nu} = (g_{\mu\nu})^{-1}$$
$$ = \eta^{\mu \alpha} \ (\delta{^\alpha_\nu} + \eta^{\alpha\beta}H_{\beta\nu})^{-1}$$
We first note, that we can contract the background metric inside the brackets: $ H{^\alpha_\nu} = \eta^{\alpha\beta}H_{\beta\nu}$. Further, to deal with the brackets, as suggested by Bob in another response, we use the binomial expansion:
$$ (1+x)^{-1} = 1 - x + x^2 -x^3 +\cdots $$
And, after a few steps of index gymnastics we reach:
$$ g^{\mu\nu} = \eta^{\mu\nu} - H^{\mu\nu} + H^{\mu\rho}H{_\rho^\nu} - H^{\mu\rho}H{_\rho^\beta}H{_\beta^\nu} + \cdots $$
Are we done?
Step-3: The expansion parameter
The beauty of this arrangement lies in the following realization:
$$ H^{\mu\nu} \xrightarrow{\text{can only give rise to terms with}} \epsilon^1, \epsilon^2, \epsilon^3 \cdots$$
$$ H^{\mu\rho}H{_\rho^\nu} \xrightarrow{\text{can only give rise to terms with}} \epsilon^2, \epsilon^3, \epsilon^4 \cdots $$
$$H^{\mu\rho}H{_\rho^\beta}H{_\beta^\nu}\xrightarrow{\text{can only give rise to terms with}} \epsilon^3, \epsilon^4, \epsilon^5 \cdots $$
Hence to get to a useful expression of the inverse, we must arrange the inverse in powers of $\epsilon$.
Doing a bit of work, we get following terms at order $\epsilon^n$:
(note overall sign comes from the last equation in step-2)
- $n=0$
$$\frac{1}{0!}(\eta^{\mu \nu}$$
- $n=1$
$$ \frac{1}{1!}(- h^{1\mu \nu}) $$
- $n=2$
$$ \frac{1}{2!}(2 h^{1}{}_{a}{}^{\nu} h^{1\mu a} - h^{2\mu \nu}) $$
- $n=3$
$$\frac{1}{3!}( -6 h^{1}{}_{a}{}^{b} h^{1}{}_{b}{}^{\nu} h^{1\mu a} + 3 h^{1\mu a} h^{2}{}_{a}{}^{\nu} + 3 h^{1}{}_{a}{}^{\nu} h^{2\mu a} - h^{3\mu \nu})$$
As should be obvious on carefully following the above treatment, the final answer neatly looks like:
$$ g^{\mu\nu} = \eta^{\mu \nu} - \epsilon h^{1\mu \nu} + \tfrac{1}{2} \epsilon^2 (2 h^{1}{}_{a}{}^{\nu} h^{1\mu a} - h^{2\mu \nu}) + \tfrac{1}{6} \epsilon^3 (-6 h^{1}{}_{a}{}^{b} h^{1}{}_{b}{}^{\nu} h^{1\mu a} + 3 h^{1\mu c} h^{2}{}_{c}{}^{\nu} + 3 h^{1}{}_{d}{}^{\nu} h^{2\mu d} - h^{3\mu \nu})$$