Can special relativity be derived from the invariance of the interval?

Yes.

Let's define that the interval needs to be the quantity constructed by adding the square of the spatial interval to the negative of the square of the temporal interval. Now, we postulate that this is frame-invariant.

To be elegant and fast, I define $\eta_{\mu\nu} = +1$ for $\mu=\nu=1,2,3$ and $\eta_{\mu\nu}=-1$ for $\mu=\nu=0$ and $\eta_{\mu\nu}=0$ for the rest of the cases. $\mu, \nu$ run from $0$ to $3$. Now, the interval is $ds^2=\eta_{\mu\nu}dx^\mu dx^\nu$ where the Einstein summation convention is in use.

Now, in a different frame, the interval will be $\eta_{\alpha \beta}dx'^{\alpha}dx'^{\beta}$. Notice that the matrix $\eta$ will remain the same because we have taken it as a definition that the quantity we call the interval needs to be constructed by adding up the square of the spatial interval to the negative of the square of the temporal interval. And now, we invoke our postulate of the invariance of this interval via demanding that $$\eta_{\mu\nu}dx^{\mu}dx^{\nu}=\eta_{\alpha\beta}dx'^{\alpha}dx'^{\beta} \tag{1}$$

Or, $$\eta_{\mu\nu}dx^{\mu}dx^{\nu}=\eta_{\alpha\beta}\dfrac{\partial x'^{\alpha}}{\partial x^\sigma}\dfrac{\partial x'^{\beta}}{\partial x^\rho}dx^\sigma dx^\rho \tag{2}$$

But since the indices $\sigma$ and $\rho$ are dummy indices, we can change them to $\mu$ and $\nu$ without any harm. So, we re-write $(2)$ as

$$\eta_{\mu\nu}dx^{\mu}dx^{\nu}=\eta_{\alpha\beta}\dfrac{\partial x'^{\alpha}}{\partial x^\mu}\dfrac{\partial x'^{\beta}}{\partial x^\nu}dx^\mu dx^\nu \tag{3}$$

The only way $(3)$ can be generically true is if the following equality holds: $$\eta_{\mu\nu}=\eta_{\alpha\beta}\dfrac{\partial x'^{\alpha}}{\partial x^\mu}\dfrac{\partial x'^{\beta}}{\partial x^\nu} \tag{4}$$ Now, as the matrix $\eta$ doesn't change with place, its differentiation with respect to any of the coordinates will be zero. Thus, $$0=\eta_{\alpha\beta}\bigg[\dfrac{\partial^2 x'^\alpha}{\partial x^{\mu}\partial {x^\kappa}}\dfrac{\partial x'^\beta}{\partial x^\nu}+\dfrac{\partial^2 x'^\beta}{\partial x^{\nu}\partial {x^\kappa}}\dfrac{\partial x'^\alpha}{\partial x^\mu}\bigg] \tag{5}$$ Since the right hand side of $(5)$ is identically zero for all the values of the independent indices, we can interchange them. So, we interchange $\mu$ and $\kappa$ and add the resulting expression to $(5)$. Similarly, we interchange $\kappa$ and $\nu$ and subtract the resulting expression from $(5)$.

The resultant is

$$0=\eta_{\alpha\beta}\bigg[\dfrac{\partial^2 x'^\alpha}{\partial x^{\mu}\partial {x^\kappa}}\dfrac{\partial x'^\beta}{\partial x^\nu}+\dfrac{\partial^2 x'^\beta}{\partial x^{\nu}\partial {x^\kappa}}\dfrac{\partial x'^\alpha}{\partial x^\mu} + \dfrac{\partial^2 x'^\alpha}{\partial x^{\kappa}\partial {x^\mu}}\dfrac{\partial x'^\beta}{\partial x^\nu}+\dfrac{\partial^2 x'^\beta}{\partial x^{\nu}\partial {x^\mu}}\dfrac{\partial x'^\alpha}{\partial x^\kappa} - \dfrac{\partial^2 x'^\alpha}{\partial x^{\mu}\partial {x^\nu}}\dfrac{\partial x'^\beta}{\partial x^\kappa}-\dfrac{\partial^2 x'^\beta}{\partial x^{\kappa}\partial {x^\nu}}\dfrac{\partial x'^\alpha}{\partial x^\mu}\bigg] \tag{6}$$

As you can see, the last term cancels the second and the fifth cancels the fourth. Notice that this argument is based on the symmetry of $\eta_{\alpha\beta}$ which is inherent in its definition. The remaining terms simplify to $$0=2\eta_{\alpha\beta} \dfrac{\partial^2x'^\alpha}{\partial x^\mu \partial x^\kappa} \dfrac{\partial x'^\beta}{\partial x^\nu} \tag{7}$$

This can be generically true only if $$ \dfrac{\partial^2x'^\alpha}{\partial x^\mu \partial x^\kappa} = 0 \tag{8} $$

The generic solution of this equation would then be $$x'^\alpha = \Lambda^\alpha _\mu x^\mu + c^\alpha \tag{9}$$ Here, $\Lambda$ and $c$ are constants in the sense they do not depend on the position in spacetime.

Therefore, $(4)$ translates to $$\eta_{\mu\nu} = \eta_{\alpha\beta}\Lambda^\alpha _\mu \Lambda^\beta _\nu \tag{10}$$ Or, in the matrix form, $$\eta = \Lambda^{T} \eta \Lambda \tag{11}$$

Thus, the allowed the set of proper homogeneous transformations is isomorphic to $SO(3, 1)$ and the set of all transformations is isomorphic to the Poincare group. All the results of special relativity can now be derived as we know the exact symmetry group of its transformations.

Edit

I gathered the arguments presented above from the wonderful GR book of Steven Weinberg, Gravitation and Cosmology.