How to organize this strong coupling perturbation theory? Consider a 2d scalar field theory with quartic interaction
$$S[\phi]=\int d^2x \left((\nabla\phi)^2+m^2(\phi^2+g\phi^4)\right)$$
I want to compute the partition function
$$ Z[m,g]=\int\mathcal{D}\phi\,e^{-S[\phi]}$$
say as a function of $m,g$. I want to do this in $m^2\to\infty$ limit keeping $g$ finite.
Short statement of the question. When $m^2$ is large the saddle point methods seem appropriate. I expect that the leading term is given by the the one-loop partition function $\log\det (-\Delta+m^2)$ of the quadratic action and that the $\phi^4$ will give further $1/m^2$ corrections. However, the naive perturbation theory leads to all loop diagrams being of the same order in $m^2$. Is there a way to organize the perturbative expansion such that it gives meaningful $1/m^2$ corrections?
My attempt to do a naive perturbative expansion and why it fails.
Since there is a large parameter in the action, I try to use the saddle point expansion. The saddle point configuration is just $\phi=0$ so the action already is written for the fluctuations about the saddle point. Next, one expects quadratic term to dominate while the quartic term to produce corrections in the form of $1/m^2$ expansion. However if I try to do naive perturbation theory this turns out to be false.
Consider a simplest diagram without self-contractions which turns out to be three-loop and write it in coordinate space
$$\left<\left(gm^2\int d^2x\phi^4\right)^2\right>\simeq g^2m^4 \int d^2x'\int d^2x G^4(x-x')\simeq g^2 m^4 V \int d^2x\,\, G^4(r)$$
Here $V$ is formally the volume of space $V=\int d^2x $. If flat space it is infinite so we could put the theory in a finite box or on a closed surface but I think these details are irrelevant.
Now, naively the propagator of a heavy field should behave as $G\propto m^{-2}$ so that $G^4\propto m^{-8}$ and the whole diagram is proportional to $m^{-4}$. However, the actual propagator for the massive field in two dimensions is up to a constant
$$G(r)=K_0(mr),\qquad\qquad(-\Delta+m^2)G(r)=\delta^{(2)}(r)$$
Here $K_0(r)$ is the zeroth order modified Bessel function. It has a logarithmic singularity at $r=0$ and decays exponentially at $r\to\infty$. So actually $\int d^2x G^4(r)=\int d^2x K^4_0(mr)\propto m^{-2}$ and not $m^{-8}$. As a result the whole diagram is proportional to $m^2$. It is easy to see by similar arguments that suppressing factors $m^{-2}$ are not associated with propagators but rather with vertices in the diagrams. However as each vertex carries a factor $m^2$ coming from the action all diagrams in the perturbation theory have the same order $m^2$.
I the quartic couplig $g$ can be treated as small, then this naive perturbative expansion is sensible as expansion in powers of $g$. However if I insist on keeping $g$ of order one is there a way to reorgonanize the expansion to get $m^{-2}$ corrections described by a finite amount of diagrams?
I should perhaps note that I have very little experience with perturbative expansions of this kind. The solution may be simple, for example to use an improved propagator, or much more complex. Pointers to the literature are also very welcome.

Above I have written a simplified version of my actual problem. From discussion with the Chiral Anomaly it appears that this toy model may not be completely adequate. More specifically I'm interested in computing a partition function
$$Z[E,g]=\int \mathcal{D}_g\phi \,\,e^{-S[\phi,E,g]}$$
where
$$S[\phi,E,g]=\int d^2x\sqrt{g} \left(\nabla^\mu\phi\nabla_\mu\phi+(m^2+E(x))\frac{e^{2b\phi}-2b\phi-1}{2b^2}\right)$$
and the action has something to do with the Liouville theory and is defined on a sphere. The answer is expected to be an expansion with terms of the type $\int_x\frac{E(x)^nR(x)^{k}}{m^{2n+2k}}$ where $E(x)$ is a "variable part" of a mass and $R(x)$ is a curvature of the space. Naive perturbation theory produces this kind of terms, but all loops seem to contribute. If however one restricts to a finite power of $b$, which is an analog of $g$ in the original problem, the finite amount of diagrams will do.
 A:  Step 1: Make things well-defined 
Things work out better when we start with something well-defined. To make the problem well-defined, I'll treat the 2d space as a lattice with a finite but arbitrarily large number of sites. (Strong-coupling expansions are typically done using lattice QFT.) Then the integration variables $\phi(x)$ are ordinary real variables, one per lattice site $x$. I'll use the following abbreviations:
$$
\newcommand{\la}{\langle}
\newcommand{\ra}{\rangle}
\newcommand{\cD}{{\cal D}}
 \big(\nabla \phi(x)\big)^2 
 \equiv 
 \sum_{u} 
 \left(\frac{\phi(x+u)-\phi(x)}{\epsilon}\right)^2
\tag{1}
$$
where the $u$s are basis vectors for the lattice, both with magnitude $\epsilon$, and
$$
 \int d^2x\ L(x)\equiv \epsilon^2\sum_x L(x),
\tag{2}
$$
and
$$
 \int\cD\phi\ F[\phi]
 \equiv \int \left(\prod_x  d\phi(x)\right)\ F[\phi].
\tag{3}
$$
The path integral $\int\cD\phi$ is now an ordinary multi-variable integral over ordinary real variables $\phi(x)$, and everything is finite. In fact, all of the terms in the $1/m^2$ expansion can be evaluated in closed form. The details are shown below.
 Step 2: The large-$m^2$ expansion 
Each $\phi(x)$ is just an integration variable, and each integral in (3) is over the whole real line, so we can replace $\phi$ with $\phi/m$ to get
$$
 Z[m,g]\propto  
 \int \cD\phi\ 
 e^{-S_0[\phi]}e^{-V[\phi]}
\tag{4}
$$
with
$$
 S_0[\phi] = \int d^2x\ \phi^2
\hspace{1cm}
 V[\phi] = \int d^2x\ 
  \frac{(\nabla\phi)^2+g\phi^4}{m^2}.
\tag{5}
$$
Things work out more nicely if we work with a normalized partition function $Z'[m,g]$ whose leading term is $1$:
$$
 Z'[m,g]\equiv \frac{Z[m,g]}{Z[\infty,g]}
 =
 \frac{\int \cD\phi\ e^{-S_0[\phi]}e^{-V[\phi]}}{
  \int \cD\phi\ e^{-S_0[\phi]}}.
\tag{6}
$$
Now expand in powers of $V$, which is the same as expanding in powers of $1/m^2$:
$$
 Z'[m,g]
 =
 \sum_{n\geq 0}\frac{(-1)^n}{n!}\, V_n
\tag{7}
$$
with
$$
 V_n\equiv
 \frac{\int \cD\phi\ e^{-S_0[\phi]}\big(V[\phi]\big)^n}{
  \int \cD\phi\ e^{-S_0[\phi]}}.
\tag{8}
$$
 Step 3: Evaluating the individual terms 
The key to evaluating individual terms in the expansion is to use the factorization
$$
 e^{-S_0[\phi]}=\exp\left(-\epsilon^2\sum_x \phi^2(x)\right)
 =\prod_x \exp\left(-\epsilon \phi^2(x)\right).
\tag{8}
$$
To see how it works, use the abbreviation
$$
 V(x)\equiv \frac{\big(\nabla \phi(x)\big)^2 +g\phi^4(x)}{m^2}.
\tag{9}
$$
The quantity $V(x)$ depends on only three of the integration variables, namely $\phi(x)$ and its two neighbors $\phi(x+u)$, one for each direction $u$. Thanks to the factorization (8) and the definition (3), the first-order term $V_1$ reduces to
$$
 V_1 
 = \epsilon^2\sum_x
 \frac{\int \cD\phi\ e^{-S_0[\phi]}
 V(x)}{
  \int \cD\phi\ e^{-S_0[\phi]}}
 = \epsilon^2\sum_x
 \frac{\int \prod_{y\in V(x)}d\phi(y)\ e^{-\epsilon\phi^2(y)}
 V(x)}{
  \int \prod_{y\in V(x)}d\phi(y)\ e^{-\epsilon\phi^2(y)}}
\tag{10}
$$
where the notation $y\in V(x)$ means all of the sites involved in $V(x)$, namely  the site $x$ and its two neighbors $x+u$. After unpacking the definitions, we see that the integrals on the right-hand side of (10) can all be evaluated in closed form. The result is invariant under shifts of the reference-site $x$, so the overal sum over $x$ just gives an overall factor of $N$, the number of lattice sites. The combination $\epsilon^2 N$ is the area covered by the whole lattice.
Things become a little more interesting at second order:
$$
 V_2
 = \epsilon^2\sum_{x_1} \epsilon^2\sum_{x_2}
 \frac{\int \cD\phi\ e^{-S_0[\phi]}
 V(x_1)V(x_2)}{
  \int \cD\phi\ e^{-S_0[\phi]}}.
\tag{11}
$$
Now we get some terms in which the factors $V(x_1)$ and $V(x_2)$ don't share any integration variables (don't share any sites), and we get some terms in which they do. For the terms that don't, we can recycle the integrals that we already evaluated in $V_1$. The terms in which they do share sites are new, but again those integrals can be evaluated in closed form. This pattern continues to arbitrarily high orders in the expansion.
 Step 4: Managing the complexity 
Because of the increasing number of factors of $V(x)$ in the integrand at higher orders, the integrals become more and more complicated because of the various ways in which the factors of $V(x)$ can share sites with each other. That's where the real fun begins. This is an ancient industry, complete with its own diagrammatic notation, so you can probably find plenty of information about it by searching for the keywords strong-coupling expansion or hopping-parameter expansion. The book Quantum Fields on a Lattice by Montvay and Münster (1994) includes an introduction to the ideas. One of the early papers is Strong-coupling expansion in quantum field theory, which says this in the abstract:

We derive a simple and general diagrammatic procedure for obtaining the strong-coupling expansion of a d-dimensional quantum field theory, starting from its Euclidean path-integral representation. At intermediate stages we are required to evaluate diagrams on a lattice; the lattice spacing provides a cutoff for the theory. We formulate a simple Padé-type prescription for extrapolating to zero lattice spacing and thereby obtain a series of approximants to the true strong-coupling expansion of the theory. No infinite quantities appear at any stage of the calculation.

