Let $F$ be the cumulative distribution function of $\left|X-Y\right|$. With the notations of the opening post, we have for each $l$:
$$\tag{*} 1-F\left(\varepsilon_l\right)\leqslant \varepsilon_l.$$
Now, since the sequence $\left(\varepsilon_l\right)_{l\geqslant 1}$ is non-increasing, we get by right continuity of a cumulative distribution function that $\lim_{l\to +\infty}F\left(\varepsilon_l\right)=F\left(d\left(X,Y\right)\right)$. We conclude by taking the limit as $l$ goes to infinity in (*).
Your proof is correct. To elaborate on the parts that were causing you discomfort, I'll say the following. See that it holds to apply the triangle inequality to the integrand because $f,g,h$ all map into $\mathbb{R}$, so for all $x \in [0,1]$, let $f(x)-h(x)$ = $a_x$ and $h(x)-g(x)$ = $b_x$, $a_x, b_x \in \mathbb{R}$. Then substituting these quantities for what you have written in the beginning of your proof we have
$$\int_0^1|a_x| \, dx + \int_0^1 |b_x| \, dx$$
$$= \int_0^1 |a_x| + |b_x| \, dx.$$
Now when we apply the triangle inequality to the integrand it looks more familiar to how we usually see the triangle inequality, and now you might feel more comfortable saying
$$|a_x+b_x| \leq |a_x| + |b_x|.$$
For the second part you mentioned that seemed unjustified, consider the following lemma. It is a known theorem of calculus that if
$$h(x) \in C[a,b], h(x) \geq 0 \text{ then } \int_a^b h(x) \, dx \geq 0.$$
Since sum of two elements of $C[a,b]$ is in $C[a,b]$ the theorem implies
$$h(x)-g(x) \geq 0 \implies \int_a^b h(x) - g(x) \, dx \geq 0$$
which applying linearity to the right hand side leads us to the statement that you had reservations about, namely
$$h(x) \geq g(x) \implies \int_a^b h(x) \, dx \geq \int_a^b g(x) \, dx.$$
Applying this lemma to where we left off on the first issue we get the desired result. I hope this makes you more convinced of your proof!
Best Answer
If $|X-Y| \le \epsilon_1$ and $|Y-Z| \le \epsilon_2$, then we have $|X-Z|\le \epsilon_1 + \epsilon_2$.
Hence if $|X-Z| > \epsilon_1 + \epsilon_2$, then $|X-Y| > \epsilon_1$ or $|Y-Z|> \epsilon_2$.
$$\{ \omega:|X(\omega) - Z(\omega)| > \epsilon_1 + \epsilon_2\} \subseteq \{\omega:|X(\omega)-Y(\omega)| > \epsilon_1\} \cup \{\omega:|Y(\omega)-Z(\omega)| > \epsilon_2\}$$
Hence, by the union bound,
$$P(|X-Y| > \epsilon_1) + P(|Y-Z| > \epsilon_2) \ge P(|X-Z| > \epsilon_1 + \epsilon_2)$$