What would be the computationally faster way to implement this 2D numerical integration?
A very low hanging fruit, namely lifting the constant factor

0.5 * (1 / (np.sqrt(2 * np.pi) * 2))


out from the integral, reduces time by a quarter. On my system the original code took 

294.532276869


while the fixed one took

224.198880911

Another fruit, a bit above, is to abandon dblquad whatsoever. It is just a wrong tool here.

Notice that the only dependency on q is in j0(q * t), and that it does not depend on z either. From the mathematical viewpoint it means that an integral of erf(...) * exp(...) over dz can be tabulated once as a function of t, say F(t) for the lack of better name, which then can be fed into the final integration as t * j0(q*t) * F(t).

Of course you'd need a lambda to interpolate tabulated values, and to manually take care of precision, and maybe do something else that dblquad does under the hood. Nevertheless, expect a thousand-fold speedup.