How do I get the velocity of the P-wave or S-wave without using the travel time graph?
I think that it is possible with fours stations yes, not with three.

With the assumption of (roughly) homogeneous velocities we get the following forward problem. Assuming that $v_p$ and $v_s$ are constants between the Earthquake and all your receivers yields the following system:
\begin{align}
    d &amp;= \sqrt{\Delta x^2 + \Delta y^2 + \Delta z^2}=\text{3 knowns, 3 unknowns},\quad \text{for} \left\{ \begin{array}{rl} \Delta x &amp;= x_\text{earthquake}-x_\text{receiver} \\ \Delta y &amp; = y_\text{earthquake} - y_\text{receiver} \\ \Delta z &amp; = z_\text{earthquake} - z_\text{receiver}  \end{array} \begin{split} \end{split}\right.,\\
\Delta t &amp; = t_s - t_p = \text{measured},\\
p &amp;= \left(\frac{1}{v_s} - \frac{1}{v_p} \right) = \text{1 unknown}.
\end{align}
Restating the equation of your question then gives the forward problem:
$$ \Delta t = p d.$$
For example, for $v_p=2000$ m/s and $v_s=1200$ m/s, we find $p=1/3000$ and can compute e.g. with $d = 600$ that $\Delta t = 600/3000 =0.2$ s.

Now, we must set=up the inverse problem. You generally need at least as many equations as unknowns. In this case we have 4 unknowns (the velocity and the $(x,y,z)$ earthquake location). That means that we need 4 (unique) equations. I think that the simplest way of obtaining 4 unique equations is to combine 4 measurements:
$$ p\begin{pmatrix} d_1 \\ d_2 \\ d_3 \\ d_4 \end{pmatrix} = \begin{pmatrix} \Delta t_1 \\ \Delta t_2 \\ \Delta t_3 \\ \Delta t_4 \end{pmatrix}, $$
Denoting with $x_{1,2,3,4}$ and $y_{1,2,3,4}$ the latitude and longitude of the receivers, and assuming measurements at the surface, we obtain:
$$ p\begin{pmatrix} \sqrt{(x_\text{earthquake}-x_1)^2 + (y_\text{earthquake}-y_1)^2 + z_\text{earthquake}^2} \\ \sqrt{(x_\text{earthquake}-x_2)^2 + (y_\text{earthquake}-y_2)^2 + z_\text{earthquake}^2} \\ \sqrt{(x_\text{earthquake}-x_3)^2 + (y_\text{earthquake}-y_3)^2 + z_\text{earthquake}^2 } \\ \sqrt{(x_\text{earthquake}-x_4)^2 + (y_\text{earthquake}-y_4)^2 + z_\text{earthquake}^2} \end{pmatrix} = \begin{pmatrix} \Delta t_1 \\ \Delta t_2 \\ \Delta t_3 \\ \Delta t_4 \end{pmatrix}. $$

Finally, we must solve the inverse problem. Unfortunately, the system is non-linear (the unknowns are 'hidden' in the root of the square and can't be separated). You must therefore use some non-linear solver. For example, in MATLAB you can do this

% True earthquake location
x_earthquake=0;
y_earthquake=0;
z_earthquake=1000;

% True station locations
x_1 = 1000; y_1 = 1000;
x_2 =  500; y_2 = -300;
x_3 = -400; y_3 = -100;
x_4 =  -10; y_4 = 800;
X = [x_1;x_2;x_3;x_4];
Y = [y_1;y_2;y_3;y_4];

% True velocity structure
p0 = 1/3000;

% Forward equation (F(1)=x_eartquake, F(2)=y_earthquake, F(3)=z_earthquake, F(4)=p. 
t =@(F) F(4) * sqrt( (F(1)-X).^2 + (F(2)-Y).^2 + F(3).^2);

% True recordings
measured_times = t([x_earthquake,y_earthquake,z_earthquake,p0]);

% Misfit function (=0 at optimum)
G =@(F) t(F) - measured_times;

% Invert
%options = optimoptions('fsolve','FiniteDifferenceType','central'); % using these will give better results in MATLAB
%F_inv = fsolve( G, [1000,1000,1000,0.1],options)
F_inv = fsolve( G, [1000,1000,1000,0.1])
x_earthquake_inv = F_inv(1)
y_earthquake_inv = F_inv(2)
z_earthquake_inv = F_inv(3)
p_inv            = F_inv(4)


If you copy this into https://octave-online.net/, for example, it will find an earthquake at (0,0,1000) and p=1/3000, with some small errors due to the non-linear solver.