How can I extract a sinusoid from a very noisy signal for computing the phase shift?
You are correct. Estimating the instantaneous phase of a noisy sinusoid is NOT easy.  I suggest you design a narrow bandpass filter such that your sinusoid-of-interest is in the filter's passband.  (The better the filter the more noise that will be eliminated.) Pass your two signals through the bandpass filter to generate filtered signals $x_1[n]$ and $x_2[n]$.  Next, pass your $x_1[n]$ and $x_2[n]$ signals through a Hilbert transformer to generate $\hat{x}_1[n]$ and $\hat{x}_2[n]$.  Create two analytic (complex) signals as:

$$z_1[n] = x_1[n] + j \, \hat{x}_1[n],$$

and  $$z_2[n] = x_2[n] + j \, \hat{x}_2[n]$$

where

$$ \begin{align}
\hat{x}[n] &amp; = \mathcal{H}\{ x[n] \} \\
           &amp; = \sum\limits_{m=-\infty}^{+\infty} \frac{1 - (-1)^{m}}{\pi \, m}   x[n-m]  \\
\end{align} $$

Next, compute two instantaneous phase sequences: 

$$\phi_1[n] = \arg\{z_1[n]\}$$ 

and $$\phi_2[n] = \arg\{z_2[n]\}.$$

Finally, compare the instantaneous phase difference between the $\phi_1[n]$ and $\phi_2[n]$ sequences.