Start with the assumption that $f$ is Riemann-Darboux integrable and, hence, bounded.
For any $\epsilon > 0$ there exists a partition $P_\epsilon = (a=x_0,x_1, \ldots, x_{m-1},x_m=b)$ such that the upper Darboux sum satisfies
$$I \leqslant U(f,P_\epsilon) < I + \frac{\epsilon}{2}$$
Since $f$ must be bounded, there exists $M > 0$ such that $-M \leqslant f(x) \leqslant M$ and $|f(x)- f(y)| \leqslant 2M$ for all $x,y \in [a,b]$.
Let $P = (a = y_0 , y_1, \ldots , y_{r-1}, y_r = b)$ be any partition with $\|P\| < \delta = \dfrac{\epsilon}{4mM},$ and take $Q = P \cup P_\epsilon$.
Since the partition $Q$ is a refinement of $P_\epsilon$ we have $U(f,Q) \leqslant U(f,P_\epsilon)$. Furthermore, $Q$ has at most $m-1$ more partition points than $P$ since the $m+1$ points of $P_\epsilon$ have been added and the endpoints $x_0 = y_0 =a$ and $x_m = y_r = b$ coincide.
The part of the proof that requires some insight is the observation that
$$\tag{*}|U(f,P) - U(f,Q)| < 2M(m-1) \delta = 2M(m-1) \frac{\epsilon}{4mM} < \frac{\epsilon}{2},$$
which implies
$$U(f,P) < U(f,Q) + \frac{\epsilon}{2} < U(f,P_\epsilon) + \frac{\epsilon}{2}
< I + \epsilon$$
Since $U(f,P) \geqslant I$ it follows that $|U(f,P) - I| < \epsilon$. The proof that $|L(f,P) - I| < \epsilon$ is similar.
Explanation of inequality (*)
This follows because the difference between $U(f,P)$ and $U(f,Q)$ comes from the area of at most $m-1$ rectangles above the graph of $f$ with height bounded by $2M$ and width bounded by $\delta$.
For example, consider the interval $[y_j, y_{j+1}]$ of $P$ and suppose that the single point $x_k$ from $P_\epsilon$ has been added in forming $Q$ and we have $y_j < x_k < y_{j+1}$.
Let $M(\alpha,\beta) := \sup_{x \in [\alpha,\beta]}\,f(x)$
The absolute difference of upper sums has the contribution
$$|U(f,Q) - U(f,P)| = \left| \,M(y_j,x_k) (x_k - y_j)+ M(x_k,y_{j+1}) (y_{j+1} - x_k) - M(y_j,y_{j+1}) (y_{j+1} - y_j)\, \right| \\ \leqslant |M(y_j,x_k)- M(y_j,y_{j+1})| (x_k - y_j)+ |M(x_k,y_{j+1})- M(y_j,y_{j+1}) |(y_{j+1} - x_k) \\ < |M(y_j,x_k)- M(y_j,y_{j+1})|\delta + |M(x_k,y_{j+1})- M(y_j,y_{j+1})| \delta $$
Of the two terms on the RHS one must vanish where suprema coincide and in the remaining term the difference of suprema is bounded by $2M$.
Thus, $|U(f,Q) - U(f,P)| < 2M \delta$ and proceeding inductively as $m-1$ points are added we have $|U(f,Q) - U(f,P)| < 2M(m-1) \delta$.