This looks like a Laplace approximation. If $\sigma^2$ is small, the integral of $f(x)$ times a normal density is essentially $f(\mu)$. You could see this by expanding $f(x)$ in a power series centered at $\mu$ and keeping only the first two terms. The constant term gives the Laplace approximation and the second term integrates to zero by symmetry.
The Laplace approximation would give $$\frac{1}{2}\left\[1 + \textrm{erf}\left(\frac{\theta - \mu}{\sqrt{2q^2}}\right)\right\]$$ $$ \frac{1}{2} \left[1 + \textrm{erf}\left( \frac{\theta - \mu}{\sqrt{2q^2}} \right)\right] $$ which isn't quite the approximation in your question. But if $\sigma^2$ is sufficiently small, the two expressions are approximately equal. Perhaps the approximation in your question is a refinement of the Laplace approximation.