Another Approach
I'm assuming here that $A$ is roughly equal to one. If it's a lot bigger or a lot smaller, then the problem does change quite a bit.
A first estimate to getting $x(t)$ would simply to take the mean of the observations
$$x_1(t) = \frac{1}{N}\sum_{n=0}^{N-1}s_n(t)$$
All the $x_0(t)$ terms should add coherently whereas the time and phase shifted terms will be most uncorrelated.
Next we take a look at the transfer functions (I'll use capital letters for frequency domain)
$$H_n(\omega) = X(\omega)\cdot (1 + A\cdot e^{j(\theta-\omega)\tau_n})$$
That's basically a complex comb filter that has maxima at $(\theta-\omega)\tau_n = 2k\pi$ and minima at $(\theta-\omega)\tau_n = (2k+1)\pi$
Now you can use your first estimate $x_1(t)$ to normalize the individual spectra,
$$S_{n,1,norm} = \frac{S_n(\omega)}{X_1(\omega)}$$
This should roughly look like a comb filter. You use the location of the minima and maxima to estimate $\tau_{n,1}$ and the height of the maxima (or depth of the minima) to estimate $A_{n,1}$.
Now you use these estimates to refine you estimate for $x(t)$
$$x_2(t) = \frac{1}{N}\sum_{n=0}^{N-1}(s_n(t)-A_{n,1}x_1(t-\tau_{n,1})e^{j\theta \tau_{n,1}})$$
Rinse and repeat until the residual error becomes small or it stops converging.