Here is an approach that can prove an average weak mixing of this system, ie:
$$\lim_{n \to \infty} \frac{1}{n}
\sum_{k=1}^n \mathbb{E} (\mu( T_k \dots T_1 (A) \cap B)) = \mu(A)\mu(B)$$
for all Borel sets $A$, $B$.
By a standard approximation argument, this is equivalent to:
$$
\lim_{n \to \infty}
\mathbb{E}(
\int (f \circ T_1 \circ \cdots \circ T_n)
\cdot
g
)
=
\int f \cdot g
,
$$
for all $f,g \in L^2([0,1])$.
We study
the averaging operator $P$:
$$
Pf = \mathbb{E}( f \circ T )
.
$$
The operator $P$ is self-adjoint on $L^2([0,1])$,
with operator norm 1, and
$$
P^n f = \mathbb{E}( f \circ T_1 \circ \cdots \circ T_n )
.
$$
By (an adaptation of) Von Neumann mean ergodic theorem we have, for every $f \in L^2([0,1])$,
the convergence in $L^2$:
$$
\frac{1}{n}
\sum_{k=1}^n
P^n f
\to
\pi(f)
,
$$
where $\pi$ is the orthogonal projection
on the closed subspace of $P$-invariant function.
In particular, the average weak mixing property above is equivalent to the non-existence of non-constant $P$-invariant function.
We can compute $P$ explicitly:
$$
Pf(x) =
2 \int_x^1
\int_0^t
f(u)dudt
+ x^2f(x)
.
$$
This reduces the question of average weak mixing to the 2nd order integral equation:
$$
f(x) =
2 \int_x^1
\int_0^t
f(u)dudt
+ x^2f(x)
.
$$
I found that there is no non-constant $L^2$ solution to this equation, by reducing to a differential equation (this should be verified !). This proves the average weak mixing.
I guess that the average mixing
$$\lim_{n \to \infty} \mathbb{E} (\mu( T_n \dots T_1 (A) \cap B)) = \mu(A)\mu(B)$$
could be deduced (if true) from a spectral gap of $P$ on some appropriate function space.
