The way to proceed is to use the definition of Marginal distribution, and also Law of the unconscious statistician,
$$p_\mathbf{R}(\mathbf{r}) = \mathbb{E}_\delta\big\{p_{\mathbf{R}\mid \delta}(\mathbf{r}\mid\delta)\big\} = \int_\delta p_{\mathbf{R}\mid \delta}\left(\mathbf{r}\mid\delta\right) p_\delta(\delta) \mathrm{d} \delta = \int_{-1}^{+1} p_{\mathbf{R}\mid \delta}\left(\mathbf{r}\mid\delta\right) \frac{1}{2} \mathrm{d} \delta \tag{1}$$
Note that for a given $\delta$, the conditional $\mathbf{r}$ is $C\mathcal{N}\left(\mathbf{0}, \sigma_s^{2} \mathbf{I}_{N} + \delta^2 \sigma_n^{2} \mathbf{I}_{N}\right)$.
The remaining work is to inject the conditional pdf to the uttermost right hand side of $(1)$ and calculate the integral.
An elegant solution would be using the nice properties of Gaussian random variables: sum of (independent) Gaussian is also Gaussian. In that way, one can push the integral to the covariance matrix part and come up with the result:
$$C\mathcal{N}\left(\mathbf{0}, \sigma_s^{2} \mathbf{I}_{N} + \frac{1}{3} \sigma_n^{2} \mathbf{I}_{N}\right)$$
Bonus. Use this script to check the result of the scalar case.
sigmas = 2; sigman = 5; N = 1e5;
s = sigmas/sqrt(2)*complex(randn(N,1),randn(N,1));
n = sigman/sqrt(2)*complex(randn(N,1),randn(N,1));
delta = 1-2*rand(N,1);
r = s + delta.*n;
histogram(real(r),'Normalization','probability')
xx = linspace(-10,10,40); yy = normpdf(xx,0,sqrt((sigmas^2+sigman^2/3)/2));
hold on;
plot(xx,yy,'ro');