There are several ways to proceed. Maybe the most elegant (if you are functional-analysis minded) is to use the fact that in $\mathbb{R}^d$ $(-\Delta)^{\alpha}$, for $0 < \alpha \leq 1$ is the generator of a Markovian semigroup satisfying Sobolev inequalities with ''dimensional constant'' given by $\alpha^{-1} \, d$. This can be seen by subordination, see the Chapter II of [1].
With that in mind, assume $s_1 < s_2$ and fix
$$ \alpha = \frac{s_1}{s_2}. $$
Now, the operator $A = (-\Delta_{d_1})^{\alpha} + (-\Delta_{d_2})$ generates a Markovian semigroup over $\mathbb{R}^{d_1} \times \mathbb{R}^{d_2}$ with dimensional constant given by $d' = d_1 \alpha^{-1} + d_2$. Your anisotropic Sobolev space will be $W^{s_2,p}_{A}(\mathbb{R}^{d_1} \times \mathbb{R}^{d_2})$ and writing the Sobolev embedding for the new dimension $d'$ gives the result.
Alternatively, you could prove directly the embedding for $p=2$ using the Plancherel inequality to see that any $f \in W^{p,s_1, s_2}$ has integrable Fourier transform for the right set of parameters. Indeed
$$
\int_{\mathbb{R}^d} |\widehat{f}(\xi)| d\xi
= \int_{\mathbb{R}^d} (1 + |\xi_1|^{s_1} + |x_2|^s_2)^{-1} (1 + |\xi_1|^{s_1} + |x_2|^s_2) |\widehat{f}(\xi)| d\xi,
$$
and applying the Cauchy-Schwartz inequality for $s_1$ and $s_2$ such that $(1 + |\xi_1|^{s_1} + |x_2|^s_2)^{-1} \in L^2(\mathbb{R}^d)$ gives the desired embedding. A complex interpolation argument will then give the result for $2 < p$ and duality for the Fourier multiplier of symbol $(1 + |\xi_1|^{s_1} + |x_2|^s_2)^{-1}$ will give its $L^1 \to L^2$ boundedness.
Varopoulos, N.Th.; Saloff-Coste, L.; Coulhon, T., Analysis and geometry on groups, Cambridge Tracts in Mathematics 100. Cambridge: Cambridge University Press (ISBN 978-0-521-08801-5/pbk). xii, 156 p. (2008). ZBL1179.22009.