We mainly consider the stability and reconstruction of convolution random sampling in multiply generated shift-invariant subspaces V-p(Phi) = {Sigma(k is an element of Zd) c(k)(T) Phi(center dot -k) : (c(k))(k is an element of Zd)(l(p)(Z(d)))r} of L-p(R-d), 1 < p < infinity, where Phi = (phi(1), phi(2),., phi(r))(T) with phi(i) is an element of L-p(R-d) and c = (c(1), c(2),., c(r))(T) with c(i) is an element of l(p)(Z(d)), i = 1, 2,., r. The sampling set {x(j)}(j is an element of N) is randomly chosen with a general probability distribution over a bounded cube C-K and the samples are the form of convolution {f * psi (x j)}(j is an element of N) of the signal f. Under some proper conditions for the generator Phi, convolution function psi and probability density function rho, we first approximate V-p(Phi) by a finite dimensional subspace V-N(p)(Phi) ={Sigma(r)(i=1) Sigma(vertical bar k vertical bar <= N) c(i)(k)phi(i)(. - k) : c(i) is an element of l(p) ([-N, N](d))}. Then we show that the sampling stability holds with high probability for all functions in certain compact subsets V-K(p) (Phi) = {f is an element of V-p (Phi) : integral(CK) vertical bar f * psi(x)vertical bar(p) dx >= (1 - delta) integral(Rd) vertical bar f * psi(x)vertical bar (x)vertical bar p dx} of V-p(Phi) when the sampling size is large enough. Finally, we prove that the stability is related to the properties of the random matrix generated by {phi(i) * psi}1(<= i <= r) and give a reconstruction algorithm for the convolution random sampling of functions in V-N(p)(Phi).