For easier readability, I'll subdivise this answer into multiple steps.
Moreover, I'll acknowledge the fact that my solution is quite long, and that there could be a faster way to handle this question, but maybe this will inspire someone else to find such a quicker way.
Step $1$: Seeing what kind of linear function we should expect:
Before going into the proof per se, let's first start from the end: what linear function do we need to expect? This way, we know at least partially what we want to prove to progress towards the final goal.
Take $(\alpha,\beta) \in \mathbb{R}^2$, and consider the linear function $x \mapsto \alpha x + \beta$.
We then have, looking at the integrals of interest in this question:
$$\begin{split}\int_y^{y+1} (\alpha x + \beta) \mathrm{d}x &= \left[\frac{\alpha x^2}{2} + \beta x\right]_y^{y+1}\\
&= \frac{\alpha}{2}\left((y+1)^2 - y^2\right) + \beta\left(y+1 - y\right)\\
&= \alpha y + \beta + \frac{\alpha}{2}\end{split}$$
and, in a similar fashion:
$$\int_y^{y+\sqrt{2}} (\alpha x + \beta) \mathrm{d}x = \sqrt{2}\alpha y + \sqrt{2}\beta + \alpha$$
Thus, if our $f$ is indeed linear, what we need is the following:
$$\begin{cases} k = \alpha \\ l = \beta + \frac{\alpha}{2}\\ m = \sqrt{2}\alpha \\ n = \sqrt{2}\beta + \alpha\end{cases}$$
This would imply that $\alpha = k$, $\beta = l - \frac{k}{2}$, and that $m = \sqrt{2}k$, but also that:
$$\begin{split} n &= \sqrt{2}\left(\beta + \frac{k}{\sqrt{2}}\right)\\
&= \sqrt{2}\left(\beta + \frac{k}{2} - \frac{k}{2} + \frac{k}{\sqrt{2}}\right)\\
&= \sqrt{2}\left(l + k\left(\frac{1}{\sqrt{2}} - \frac{1}{2}\right)\right)\end{split}$$
We'll want to show these last two identities, so that we only have to deal with $k$s and $l$s in the final expression. However, the first one, $m = \sqrt{2}k$, is way easier to show than the one about $n$ though (or at least seems easier, I could very well have missed a shortcut!).
Step $2$: Showing that $m = \sqrt{2}k$:
Let's begin the actual proof.
We can go from $0$ to $1 + \sqrt{2}$ in $2$ ways, either by adding $1$ first and then adding $\sqrt{2}$, or by doing these actions in reverse order. This provides two expressions for the integral $\int_0^{1 + \sqrt{2}} f(x) \mathrm{d}x$:
$$\begin{split} \int_0^{1+\sqrt{2}} f(x) \mathrm{d}x &= \int_0^{1} f(x) \mathrm{d}x + \int_1^{1+\sqrt{2}} f(x) \mathrm{d}x\\
&= k \cdot 0 + l + m \cdot 1 + n = l + n + m\\
&= \int_0^{\sqrt{2}} f(x) \mathrm{d}x + \int_{\sqrt{2}}^{1+\sqrt{2}} f(x) \mathrm{d}x\\
&= m \cdot 0 + n + k \cdot \sqrt{2} + l = l + n + \sqrt{2} k\end{split}$$
We can deduce from this that, indeed: $m = \sqrt{2} k$.
Step $3$: Convexity/concavity, and thus continuity, of the function $F : y \mapsto \int_0^y f(x) \mathrm{d}x$:
Name $F$ the function $F : y \mapsto \int_0^y f(x) \mathrm{d}x$.
This is a quick step if you already know about the below lemma:
Lemma $\mathrm{A}$:
Let $h : (a,b) \to \mathbb{R}$ be an increasing (resp. decreasing) function and $c \in (a,b)$ (with $-\infty \leq a < b \leq +\infty$).
Then the function $x \in (a,b) \mapsto \int_c^x h(x) \mathrm{d}x$ is convex (resp. concave).
It then suffices to apply this lemma on $h = f$ and $(a,b,c) = (-\infty, +\infty, 0)$ to get the convexity or concavity of $F$, and thus its continuity, on the whole of $\mathbb{R}$.
For some demonstrations, you can refer to this question for example: A sufficient condition for convexity: $f$ absolutely continuous on $[a,b]$ and $f'$ increasing a.e. imply $f$ is convex
Step $4$: Density of the set $\mathcal{A} := \left\{a + b\sqrt{2} \mid (a,b) \in \mathbb{Z}\right\}$ in $\mathbb{R}$:
Name $\mathcal{A}$ the set $\mathcal{A} := \left\{a + b\sqrt{2} \mid (a,b) \in \mathbb{Z}\right\}$.
This is also a swift step if you already know about the below lemma:
Lemma $\mathrm{B}$:
Let $A$ be a subgroup of $(\mathbb{R}, +)$.
Then, either $\inf\left(A \cap \mathbb{R}_+^*\right) > 0$, in which case $A$ is of the form $a\mathbb{Z}$ for $a := \inf\left(A \cap \mathbb{R}_+^*\right)$, or $\inf\left(A \cap \mathbb{R}_+^*\right) = 0$, in which case $A$ is dense in $\mathbb{R}$.
It then suffices to apply this lemma to the set $A = \mathcal{A}$ to obtain the density of $\mathcal{A}$, as it is a subgroup of $(\mathbb{R}, +)$ that's not of the form $a\mathbb{Z}$ by the irrationality of $\sqrt{2}$.
For some demonstrations, you can refer to these questions or parse through the posts linked to/from them for example (note that this would work with any irrational, not just $\sqrt{2}$): Subgroup of $\mathbb{R}$ either dense or has a least positive element?, Proving that $m+n\sqrt{2}$ is dense in $\mathbb R$
Step $5$: Evaluating $F$ at the reals of the form $a + b\sqrt{2}$, $(a,b) \in \mathbb{Z}^2$:
Let $(a, b) \in \mathbb{Z}^2$. With the assumption we are given on $f$, just like for Step $2$, we know we can "add $1$" or "add $\sqrt{2}$, thus we can reach $a + b\sqrt{2}$ with $F$, in the sense below:
$$\begin{split} F\left(a + b\sqrt{2}\right) &= \int_0^{a + b\sqrt{2}} f(x) \mathrm{d}x\\
&= \int_0^{a} f(x) \mathrm{d}x + \int_a^{a + b\sqrt{2}} f(x) \mathrm{d}x\\
&= \sum_{i = 0}^{a - 1} \int_i^{i+1} f(x) \mathrm{d}x + \sum_{j = 0}^{b - 1} \int_{a+j\sqrt{2}}^{a + (j+1)\sqrt{2}} f(x) \mathrm{d}x\\
&= \sum_{i = 0}^{a - 1} (ki + l) + \sum_{j = 0}^{b - 1} \left(\sqrt{2}k \left(a+i\sqrt{2}\right) + n\right)\\
&= k \frac{a(a-1)}{2} + al + \sqrt{2}kab + 2k \frac{b(b-1)}{2} bn\\
&= \frac{k}{2} \left(a^2 + 2 \cdot a\sqrt{2}b + 2b^2\right) + al + bn - \frac{ka}{2} - kb\\
&= \frac{k}{2} \left(a + b\sqrt{2}\right)^2 + \left(l - \frac{k}{2}\right)a + b\sqrt{2} \left(\frac{n}{\sqrt{2}} - \frac{k}{\sqrt{2}}\right)\end{split}$$
For the sake of simplicity, what I've just written is for the case where $a$ and $b$ are positive, but you can manually check that most of the calculations are also the same for all integers $a$ and $b$.
We now introduce the quantity $r := \frac{n}{\sqrt{2}} - \left(l + k\left(\frac{1}{\sqrt{2}} - \frac{1}{2}\right)\right)$, inspired by the formula for $n$ obtained in Step $1$, so that, once introduced in the calculation above:
$$F\left(a + b\sqrt{2}\right) = \frac{k}{2}\left(a + b\sqrt{2}\right)^2 + \left(l - \frac{k}{2}\right)\left(a + b\sqrt{2}\right) + \sqrt{2}br$$
Step $6$: Obtaining $r = 0$, or equivalently $n = \sqrt{2}\left(l + k\left(\frac{1}{\sqrt{2}} - \frac{1}{2}\right)\right)$, and, thus, that $F = G$ on $\mathcal{A}$:
Name $G$ the function $G : y \in \mathbb{R} \mapsto \frac{k}{2} y^2 + \left(l - \frac{k}{2}\right) y$, and let $H := F - G$.
$G$ is polynomial thus continuous, and we have established that $F$ was continuous with Step $3$, therefore $H$ is continuous on $\mathbb{R}$.
In addition, we got from last step that:
$$\forall (a,b) \in \mathbb{Z}^2,\quad H\left(a+b\sqrt{2}\right) = \sqrt{2}br$$
In particular: $H(0) = 0$.
Let's show that $r$ is actually equal to $0$.
Let $\varepsilon > 0$.
Thanks to the continuity of $H$ at $0$, there exists $\delta > 0$ such that:
$$\forall y \in \mathbb{R},\quad |y| \leq \delta \,\Rightarrow\, |H(y)| \leq \sqrt{2}\varepsilon$$
Moreover, due to the fact that $0 = \inf\left(\mathcal{A} \cap \mathbb{R}_+^*\right)$ (see Step $4$), we can choose a sequence $\left((a_n,b_n)\right)_n \in \left(\mathbb{Z}^2\right)^\mathbb{N}$ such that the corresponding sequence $\left(a_n + b_n\sqrt{2}\right)_n$ tends to $0$ when $n$ tends to $+\infty$ with all terms strictly greater than $0$.
WLOG, we can impose that $|a_n + b_n\sqrt{2}| \leq \delta$ for all $n$.
This means that: $$\forall n \in \mathbb{N},\quad |H\left(a_n + b_n\sqrt{2}\right)| = \left|\sqrt{2}b_n r\right| \leq \sqrt{2}\varepsilon$$
Now, if all $b_n$s were equal to $0$, then we would have $a_n \to 0$ with the $a_n$s being integers, and thus the sequence $(a_n)_n$ would be eventually constant equal to $0$, which contradicts our assumption that no term of $a_n + b_n\sqrt{2}$ would be equal to $0$, therefore there exists at least one $N$ such that $b_N \neq 0$.
Yet, $b_N$ is also an integer, hence $|b_N| \geq 1$.
From this we can conclude that: $|r| \leq \frac{\varepsilon}{|b_N|} \leq \varepsilon$.
Since this is true for all $\varepsilon > 0$, we get: $|r| = 0$, thus $r = 0$. In other words: $H \equiv 0$ on $\mathcal{A}$, or also: $F = G$ on $\mathcal{A}$.
Step $7$: Final result:
$F$ and $G$ being continuous on $\mathbb{R}$ and being equal on a dense subset (it being $\mathcal{A}$), we have that $F = G$ on the whole of $\mathbb{R}$, which is to say:
$$\forall y \in \mathbb{R},\quad \int_0^y f(x) \mathrm{d}x = F(y) = \frac{k}{2} y^2 + \left(l - \frac{k}{2}\right) y = \int_0^y \left(kx + l - \frac{k}{2}\right) \mathrm{d}x$$
By a classic reasoning on simple functions, this gives: $f(x) = kx + l - \frac{k}{2}$ almost everywhere.
We now have two monotonic functions being equal almost everywhere. By this post: If two functions are both increasing and equal almost everywhere, do they have the same discontinuous points?, and after convincing ourselves that two monotonic functions equal almost everywhere have the same monotonicity in order to apply the aforementioned post's result (you can show that if an increasing function and a decreasing function are equal almost everywhere, then they are equal and constant everywhere and thus end up having the same monotonicity at the end of the day), we have that $f$ and $x \mapsto kx + l - \frac{k}{2}$ have the same discontinuity points, but the latter is a continuous function, therefore $f$ is continuous.
Finally, two continuous functions equal almost everywhere are equal everywhere, meaning that:
$$\forall x \in \mathbb{R},\quad f(x) = kx + l - \frac{k}{2}$$