NB: Throughout what follows, we denote the standard normal density function by $\varphi(x) = (2\pi)^{-1/2} \exp(-\frac{1}{2} x^2)$ and its cumulative distribution function by $\Phi(x)$.
Method 1: Avoid calculus.
This is a consequence of a much more general result and, in fact, has nothing in particular to do with the normal distribution at all. Here is a slightly simplified version.
Let $F$ be a strictly increasing distribution function on some interval $(a,b)$ such that $F(a) = 0$ and $F(b) = 1$. We allow $a = -\infty$ and $b=\infty$, so that we can handle the case where $F$ is defined on the entire real line, as in your example.
Suppose $X$ is a random variable with distribution function $F$. Then, $Y = F(X)$ has a uniform distribution on $(0,1)$. The proof is simple. For $y \in (0,1)$,
$$
\renewcommand{\Pr}{\mathbb{P}}
\Pr(Y \leq y) = \Pr( F(X) \leq y) = \Pr(X \leq F^{-1}(y)) = F( F^{-1}(y) ) = y \> ,
$$
where the inverse $F^{-1}$ exists by the hypothesis that $F$ is strictly increasing on $(a,b)$.
Hence, $Y$ is distributed uniformly on $(0,1)$ and, as a consequence, $\newcommand{\e}{\mathbb{E}}\e Y = 1/2$.
Note that for your particular case, you start with $X \sim \mathcal{N}(a, \sigma^2)$ and so $Y = \Phi((X-a)/\sigma)$. Your problem can easily be seen to be equivalent to asking for $\e (1 - Y) = \e Y = 1/2$.
Method 2: Hammer and tongs (i.e., use calculus).
Note that your integral can be written as
$$
\int_{-\infty}^\infty \int_{-\infty}^{- (x-a)/\sigma} \varphi(y) \frac{1}{\sigma} \varphi((x-a)/\sigma) \newcommand{\rd}{\mathrm{d}} \,\rd y \, \rd x = \int_{-\infty}^\infty \int_{-\infty}^{\sigma y + a} \varphi(y) \frac{1}{\sigma} \varphi((x-a)/\sigma) \,\rd x \, \rd y \> .
$$
Now, change variables using $u = (x-a)/\sigma$, which gives the integral
$$
\int_{-\infty}^\infty \int_{-\infty}^y \varphi(u) \varphi(y) \, \rd u \,\rd y \>.
$$
Exchanging the order of the iterated integrals, we get
$$
\int_{-\infty}^\infty \int_{-\infty}^y \varphi(u) \varphi(y) \, \rd u \,\rd y = \int_{-\infty}^\infty \int_u^\infty \varphi(u) \varphi(y) \, \rd y \,\rd u \>.
$$
But, since $\varphi$ is a probability density function
$$
\int_{-\infty}^\infty \int_{-\infty}^y \varphi(u) \varphi(y) \, \rd u \,\rd y + \int_{-\infty}^\infty \int_y^\infty \varphi(u) \varphi(y) \, \rd u \,\rd y = 1 \>,
$$
and so, by symmetry, the integral must be $1/2$.
(All exchanges of order of integration are valid by Fubini's theorem.)
For every positive $a$, the change of variable $v=au$ yields $F(0)=G(b)/a$ with
$$
G(b)=\int_{-\infty}^0\varphi(v-b)\Phi(v+b)\mathrm dv.
$$
First, some easy remarks: $G(0)=\frac12\left.\Phi(v)^2\right|_{-\infty}^0=\frac18$, and, since $\Phi\leqslant1$, $G(b)\leqslant\Phi(-b)$, in particular $G(+\infty)=0$.
Now, the computation of $G$. Note that $\Phi'=\varphi$ and $\varphi'(u)=-u\varphi(u)$, hence $G'(b)=H(b)+K(b)$ with
$$
H(b)=\int_{-\infty}^0(v-b)\varphi(v-b)\Phi(v+b)\mathrm dv,
$$
and
$$
K(b)=\int_{-\infty}^0\varphi(v-b)\varphi(v+b)\mathrm dv.
$$
Integrating by parts $H(b)$ yields
$$
H(b)=\left.-\varphi(v-b)\Phi(v+b)\right|_{-\infty}^0+\int_{-\infty}^0\varphi(v-b)\varphi(v+b)\mathrm dv,
$$
that is,
$$
H(b)=-\varphi(b)\Phi(b)+K(b).
$$
Now, $\varphi(v-b)\varphi(v+b)=\varphi(b\sqrt2)\varphi(v\sqrt2)$, hence
$$
\sqrt2K(b)=\varphi(b\sqrt2)\left.\Phi(v\sqrt2)\right|_{-\infty}^0=\frac1{2}\varphi(b\sqrt2).
$$
Thus,
$$
G'(b)=-\varphi(b)\Phi(b)+2K(b)=-\varphi(b)\Phi(b)+\frac1{\sqrt2}\varphi(b\sqrt2).
$$
Together with the limit values computed above, this shows that, for every positive $a$ and for every $b$,
$$
\int_{-\infty}^0\varphi(au-b)\Phi(au+b)\mathrm du=\frac{\Phi(b\sqrt2)-\Phi(b)^2}{2a}.
$$
The same technique seems to yield, not explicit formulas for $G_x(b)$, but some funny-looking duality relations between $G_x(b)$ and $G_b(x)$, where
$$
G_x(b)=\int_{-\infty}^x\varphi(v-b)\Phi(v+b)\mathrm dv,
$$
namely, one might have (but this should be checked)
$$
G_x(b)+G_b(x)=\Phi(b\sqrt2)\Phi(x\sqrt2).
$$
Best Answer
As already explained, when $a\gt0$ the full integral is $1-\Phi\left(b/\sqrt{a^2+1}\right)$. The same approach shows that the integral considered here is $$ I(y)=P(Y\leqslant(X-b)/a,X\leqslant y), $$ where $(X,Y)$ are i.i.d. standard normal, that is, $$ I(y)=P(aY+b\leqslant X\leqslant y). $$ I see no reason to expect more explicit formulas.