How Does Gaussian Blur Affect Image Variance
In Gaussian blur the value of each output pixel is calculated as a weighted sum of all input pixels:
$$\text{out}(x, y) = \sum^\infty_{j = -\infty} \sum^\infty_{i = -\infty} \frac{1}{{2\pi \sigma_G^2}} e^{-\frac{i^2 + j^2}{2 \sigma_G^2}}\text{in}(x+i,y+j).$$
We want to calculate the variance $\text{Var}[\text{out}(x, y)] = \sigma^2_f$ based on the variances $\text{Var}[\text{in}(x + i, y + i)] = \sigma^2_0$. We require that $\text{in}(x + i, y + i)$ are independent, which you should mention in your question. Based on the identity (variance of linear combination):
$$Var\left[\sum_i c_i X_i\right] = \sum_i c^2_i Var[X_i],$$
where $c_i$ are constants and $X_i$ are independent random variables,
$$\sigma^2_f = \sigma^2_0\sum^\infty_{j = -\infty} \sum^\infty_{i = -\infty} \left(\frac{1}{{2\pi \sigma_G^2}} e^{-\frac{i^2 + j^2}{2 \sigma_G^2}}\right)^2.$$
For a large $\sigma_G$, the squared Gaussian is rather smooth and the sum can be approximated by an integral as:
$$\sigma^2_f \approx \sigma^2_0\int^\infty_{-\infty} \int^\infty_{-\infty} \left(\frac{1}{{2\pi \sigma_G^2}} e^{-\frac{i^2 + j^2}{2 \sigma_G^2}}\right)^2\,\text{d}i\,\text{d}j
= \frac{\sigma^2_0}{4\pi \sigma_G^2}\\
\implies \sigma_f \approx \frac{\sigma_0}{2\sigma_G\sqrt{\pi}},$$
which is the same result as in Wikipedia. You can calculate the integral using Wolfram Alpha.

Note that the noise needs not be Gaussian for the above to apply.