You might try HCubature with a vector-valued integrand to evaluate h(z) at a bunch of different z values simultaneously (so that it only computes G once). That is, if Z is a vector, do
H = hcubature(x -> ϕ.(x[1], Z .- x[2]) .* G(x[1], x[2]), a, b)
to get a vector H corresponding to h(z) at each z \in Z.
To decide which z values to evaluate, if your convolved function h(z) is smooth then a good choice is Chebyshev interpolation points, as implemented by ApproxFunc.jl … that is, given the values of h(z) at appropriate Chebyshev points, you can construct an exponentially accurate interpolating polynomial.
Here, however, you can’t use the built-in ApproxFun constructors, because you need to evaluate h(z) at all the Chebyshev points at once. There used to be an ApproxFun tutorial on how to do this, but it got removed (presumably it bitrotted), unfortunately. There is a FAQ that explains how to do this.