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)
$$S_n(\omega) = X(\omega)\cdot (1 + A\cdot e^{j(\theta-\omega)\tau_n})$$
So transfer function between $S$ and $X$ basically a complex comb filter that has maxima at $(\theta-\omega)\tau_n = 2k\pi$. If $A< 1$ that it also has minima at $(\theta-\omega)\tau_n = (2k+1)\pi$.
Now you can use your first estimate $x_1(t)$ to estimate the transfer function,
$$H_{n,1} = \frac{S_n(\omega)}{X_1(\omega)}$$
This should roughly look like a comb filter. You can use the location of the maxima and maybe the maxima to estimate $\tau_{n,1}$ and the height of the maxima 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.
