The integral of $f(0.x_1 x_2 \cdots)= 0.\sigma(x_1)\sigma(x_2)\cdots$
Use $$ \int_0^1 f(x)\,\mathrm{d}x=\int_0^1\mathcal{L}(\{f>t\})\,\mathrm{d}t $$ where $\mathcal{L}$ is the usual Lebesgue measure, and note that $\mathcal{L}(\{f>t\})$ is nonincreasing. Since $\mathcal{L}(\{f>m\cdot 10^{-n}\})=1-m\cdot 10^{-n}$ for all $0\leq m<10^n$ ($f>m\cdot 10^{-n}$ forbids exactly $m$ possible initial decimal expansion $0.x_1x_2\dots x_n$), we must have $\mathcal{L}(\{f>t\})=1-t$ for all $t\in[0,1]$.
Regard $x_n$ as function of $x$. Then
- Each $x_n$ is integrable,
- $\int_{0}^{1} \sigma(x_n) \, \mathrm{d}x = \frac{1}{10}\left(\sum_{d=0}^{9}d\right) = \int_{0}^{1} x_n \, \mathrm{d}x$, and
- $f(x) = \sum_{n=1}^{\infty} \frac{1}{10^n}\sigma(x_n)$ converges uniformly by the Weierstrass M-test.
So it follows that
$$ \int_{0}^{1} f(x) \, \mathrm{d}x = \sum_{n=1}^{\infty} \frac{1}{10^n} \int_{0}^{1} \sigma(x_n) \, \mathrm{d}x = \sum_{n=1}^{\infty} \frac{1}{10^n} \int_{0}^{1} x_n \, \mathrm{d}x = \int_{0}^{1} x \, \mathrm{d}x = \frac{1}{2}. $$
We will give an other proof that only relies on Riemann sums.
We divide $[0,1]$ into the $10^N$ intervals $[(j-1)/10^N,j/10^N],\ j=1,\ldots,10^N.$ Since $f(x)$ is given by applying a permutation of $\{0,\ldots,9\}$ to the digits of the decimal expansion of $x$, we easily observe that
$$\left\{k/10^N,k=0,\ldots,10^N\right\}=\left\{f(k/10^N),k=0,\ldots,10^N\right\}.$$
Hence,
$$\frac{1}{10^N}\sum_{k=1}^{10^N}f\left(\frac{k}{10^N}\right)=\frac{1}{10^N}\sum_{k=0}^{10^N}\frac{k}{10^N}-\frac{1}{10^N}f(0)=\frac{1}{2}-10^{-N}(f(0)+1/2).$$
The left hand side is a Riemann sum for $f$. So, taking the limit as $N\rightarrow +\infty$, we obtain that
$$\int_0^1f(x)dx=\frac{1}{2}.$$
Of course, we can, indeed, take the limits, because integrability was proved by the O.P. in the post.