The situation may be systematically analyzed by looking at the associated operators: for each $\epsilon > 0$, define $T_{\epsilon}$ by
$$ \forall u \in L^1(\mathbb{R}^n), \qquad T_{\epsilon}u (x) = \frac{1}{|B_{\epsilon}|} \int_{B_{\epsilon}(x)} u(z) \, \mathrm{d}z. $$
It is clear that $T_{\epsilon}$ is linear. Moreover, writing $\mathbf{1}_{A}$ for the indicator function (a.k.a. characteristic function) of the set $A$, we find that
\begin{align*}
\| T_{\epsilon} u \|_{L^1}
&\leq \int_{\mathbb{R}^n} \left( \frac{1}{|B_{\epsilon}|} \int_{B_{\epsilon}(x)} |u(z)| \, \mathrm{d}z \right) \, \mathrm{d}x \\
&= \int_{\mathbb{R}^n} \left( \frac{1}{|B_{\epsilon}|} \int_{\mathbb{R}^n} |u(z)| \mathbf{1}_{B_{\epsilon}(x)}(z) \, \mathrm{d}z \right) \, \mathrm{d}x \\
&= \int_{\mathbb{R}^n} |u(z)| \left( \frac{1}{|B_{\epsilon}|} \int_{\mathbb{R}^n} \mathbf{1}_{B_{\epsilon}(x)}(z) \, \mathrm{d}x \right) \, \mathrm{d}z \tag{Fubini-Tonelli} \\
&= \int_{\mathbb{R}^n} |u(z)| \, \mathrm{d}z
= \| u \|_{L^1}.
\end{align*}
In other words, $T_{\epsilon}$ is a continuous linear operator on $L^1(\mathbb{R}^n)$ which is a contraction. Once this is established, the standard machinery for this type of argument is to first prove the desired result on a dense subset of nice functions and then extend to all of $L^1(\mathbb{R}^n)$ by continuity of $T_{\epsilon}$'s.
Indeed, let $C_c(\mathbb{R}^n)$ denote the set of all compactly supported continuous functions on $\mathbb{R}^n$. It is well-known that $C_c(\mathbb{R}^n)$ is dense in $L^1(\mathbb{R}^n)$. Now let $u \in C_c(\mathbb{R}^n)$.
By the uniform continuity of $u$, it is easy to check that $T_{\epsilon}u \to u$ uniformly as $\epsilon \to 0^+$ on any compact subset of $\mathbb{R}^n$. Indeed, for any compact subset $K$ of $\mathbb{R}^n$,
$$ \sup_{x \in K} \left| T_{\epsilon}u(x) - u(x) \right| \leq \sup_{\substack{x \in K \\ |y - x| \leq \epsilon}} \left| u(x) - u(y) \right| $$
and this bound converges to $0$ as $\epsilon \to 0^+$ by uniform continuity of $u$.
The previous observation implies that $\| T_{\epsilon}u - u \|_{L^1} \to 0$. Indeed, let $L = \overline{\{ x \in \mathbb{R}^n : u(x) \neq 0 \}}$ be the support of $u$, which is compact. If we enlarge $L$ to get $K = \{ x \in \mathbb{R}^n : \operatorname{dist}(x, L) \leq 1 \}$, then for any $\epsilon \in (0, 1)$, $T_{\epsilon}u$ is supprted on $K$ and $T_{\epsilon}u \to u$ uniformly on $K$. So
$$ \| T_{\epsilon} u - u \|_{L^1} = \int_{K} | T_{\epsilon}u(x) - u(x)| \, \mathrm{d}x \leq \operatorname{Leb}(K) \left( \sup_{x \in K} \left| T_{\epsilon}u(x) - u(x) \right| \right), $$
where $\operatorname{Leb}(K)$ denotes the Lebesgue measure of $K$. This bound converges to $0$ as $\epsilon \to 0^+$.
So $T_{\epsilon}u \to u$ in $L^1$ when $u \in C_c(\mathbb{R}^n)$. For general $u \in L^1(\mathbb{R}^n)$, let $\varphi \in C_c(\mathbb{R}^n)$ be arbitrary and notice that
\begin{align*}
\| T_{\epsilon}u - u \|_{L^1}
&\leq \| T_{\epsilon}u - T_{\epsilon}\varphi\|_{L^1} + \| T_{\epsilon}\varphi - \varphi\|_{L^1} + \| \varphi - u \|_{L^1} \\
&\leq \| T_{\epsilon}\varphi - \varphi\|_{L^1} + 2\| \varphi - u \|_{L^1}.
\end{align*}
So, taking $\limsup$ as $\epsilon \to 0^+$, we have
$$ \limsup_{\epsilon \to 0^+} \| T_{\epsilon}u - u \|_{L^1} \leq 2\| \varphi - u \|_{L^1}. $$
But since the left-hand side is independent of the choice of $\varphi$, letting $\varphi \to u$ in $L^1(\mathbb{R}^n)$ proves the desired claim.