The physical reason IS the constancy of the velocity of light... since I'm writing in a tablet the answer won't be complete, but expect to get you to the mathematical cross-road.
Constancy of velocity of light implies that
\begin{equation}
\frac{d|\vec{x}|}{dt} = c, \quad\Rightarrow\quad d|\vec{x}| = c\,dt.
\end{equation}
Since $d|\vec{x}| = \sqrt{dx^2 + dy^2 + dz^2}$, it follows that
\begin{equation}
dx^2 + dy^2 + dz^2 = c^2\,dt^2 \quad\Rightarrow\quad dx^2 + dy^2 + dz^2 - c^2\,dt^2 = 0.
\end{equation}
From here it is straightforward to see that the set (or group) of transformations preserving this quantity are those known as Lorentz transformations.
Now I leave you to analize the "generalization" to nonvanishing intervals, for massive particles. Hint: define a four-dimensional metric!
(continuation... after a few days)
The interval
As exposed previously, the physical condition of constancy of the speed of light leads to the conclusion that
All equivalent observers are connected through a transformation which
keep the quantity $$dx^2 + dy^2 + dz^2 - c^2\,dt^2 = 0.$$
This can be generalized to the preservation of the quantity $$ I = dx^2 + dy^2 + dz^2 - c^2\,dt^2, $$ called interval.
Notice that the interval can be written as
$$
I = X^t\, \eta\, X =
\begin{pmatrix}
ct & x & y & z
\end{pmatrix}
\begin{pmatrix}
-1 & 0 & 0 & 0 \\
0 & 1 & 0 & 0 \\
0 & 0 & 1 & 0 \\
0 & 0 & 0 & 1
\end{pmatrix}
\begin{pmatrix}
ct \\
x \\
y \\
z
\end{pmatrix}
$$
Invariance of the interval
In order for the interval to be invariant under a transformation $X' = M\, X$, one needs to
\begin{align}
I &= I' \notag \\
X^t\, \eta\, X &= (M\, X)^t\, \eta\, M\,X \notag \\
&= X^t\, M^t\, \eta\, M\,X \notag \\
\Rightarrow\quad \eta &= M^t\, \eta\, M. \tag{*}
\end{align}
Therefore the problem is to find a set of transformations $M$ satisfying Eq. (*).
Two-dimensional case
Finding a general 4 by 4 matrix $M$ preserving the Minkowskii metric ($\eta$) requires a lot of algebra, but one can easily find the transformation preserving the 2 by 2 restriction to the $(ct,x)$-plane.
Propose a matrix $$M = \begin{pmatrix} a & b \\ c & d \end{pmatrix},$$, and solve the equation
$$
\begin{pmatrix} -1 & 0 \\ 0 & 1 \end{pmatrix} =
\begin{pmatrix} a & c \\ b & d \end{pmatrix}
\begin{pmatrix} -1 & 0 \\ 0 & 1 \end{pmatrix}
\begin{pmatrix} a & b \\ c & d \end{pmatrix} =
\begin{pmatrix} -a^2+c^2 & -ab+cd \\ -ab+cd & -b^2+d^2 \end{pmatrix},
$$
which is simple if one uses the identities $\cosh^2\theta - \sinh^2\theta = 1$, and the condition $M(\theta\to 0) = \mathbf{1}$.
Thus, $$M = \begin{pmatrix} \cosh\theta & -\sinh\theta \\ -\sinh\theta & \cosh\theta \end{pmatrix}.$$
Relation with the velocity
In a similar fashion of Euclidean geometry, in which
$$
\frac{y}{x} = \tan\theta,
$$
one uses the transformation $M$ above to relate the $ct$ coordinate with the $x$ coordinate
$$
\frac{v}{c} \equiv \frac{x}{ct} = \mathop{\mathrm{tanh}}\theta.
$$
Now,
\begin{align}
\mathop{\mathrm{tanh}^2}\theta &= 1 - \mathop{\mathrm{sech}^2}\theta \notag \\
&= 1 - \tfrac{1}{\cosh^2\theta} \notag \\
\Rightarrow\quad \cosh\theta &= \frac{1}{\sqrt{1 - \left(\frac{v}{c}\right)^2}} \notag \\
\sinh\theta &= \frac{\frac{v}{c}}{\sqrt{1 - \left(\frac{v}{c}\right)^2}}. \notag
\end{align}
Finally, from the relation $X' = M\, X$, one obtain the usual relations
\begin{align}
x' &= -\sinh\theta\cdot t +\cosh\theta\cdot x \notag\\
&= \frac{1}{\sqrt{1 - \left(\frac{v}{c}\right)^2}}\left( x - vt \right) \notag \\
t' &= \cosh\theta \cdot t -\sinh\theta\cdot \tfrac{x}{c} \notag\\
&= \frac{1}{\sqrt{1 - \left(\frac{v}{c}\right)^2}}\left( t - \tfrac{v}{c^2}t \right). \notag
\end{align}