Elaborating a little bit on what I discussed in the comments: Let's assume the delays $d_1$, $d_2$ and $d_3$ are integer multiples of the sampling time interval $t_0$ so that we can treat it as an integer problem.
Then, as a first step, I would try to use the cross correlation between $y_1$ and $y_2$ to determine the delays. Ideally, it should have a peak at an offset of $d_2$ (when the shifted copies of $s_1(t)$ and $s_1(t-d_2)$ align) and another one at $d_3-d_1$ (when the shifted copies of $s_2(t-d_1)$ and $s_2(t-d_3)$ align).
Now compute $z[n] = y_1[n] - y_2[n-(d_3-d_1)]$ (using $[n]$ as a discrete time / sample index for clarity). If your estimate was correct, this cancels $s_2(t)$ and we are left with $z[n] = s_1[n] - s_1[n-\Delta]$ where $\Delta = d_3-d_1-d_2$.
Now, your signals ought to start somewhere (causality) so that we can assume $s_1[n]=0$ for $n<0$ (if you have different initial conditions, you need to work them in here). Then, we have $$z[n] = s_1[n]-s_1[n-\Delta] = s_1[n]\quad \forall n<\Delta.$$ Therefore $$s_1[n] = \begin{cases} 0 & n<0 \\ z[n] & n<\Delta \\ z[n]+s_1[n-\Delta] & \mbox{otherwise}\end{cases}.$$ This allows you to reconstruct $s_1[n]$ sample by sample.