I have read that for linear regression the Least Squares estimator is also the Maximum Likelihood Estimator under the assumption of normally distributed errors.

However is this also true for my nonlinear problem?

You are using

$$ \sqrt{(x_i - x_b)^2+(y_i - y_b)^2} = d_i +\epsilon_i$$

but it is possibly better to state it as the following instead

$$ \begin{array}{rcccl} d_i &=& \sqrt{(x_i - x_b)^2+(y_i - y_b)^2} &+& \epsilon_i \\ &=& \tilde{d}_i& + &\epsilon_i \end{array}$$

Then an observed distance $d_{i}$ is distributed according to the modeled mean distance $\tilde{d}_i$ plus some noise $\epsilon_i$.

Then the probability density (and the related likelihood) for the observations $d_i$ can be written as

$$\prod_{\forall i} \frac{1}{\sqrt{2\pi \sigma^2}} \exp \left( -\frac{(d_i-\tilde{d}_i)^2}{2\sigma^2} \right)$$

where the $\tilde{d}_i = \sqrt{(x_i - x_b)^2+(y_i - y_b)^2}$ can be modeled as functions of your unknown parameters $y_b$ and $x_b$.

Minimizing this likelihood is equivalent to minimizing the sum of squares $\sum_{\forall i} (d_i-\tilde{d}_i)^2$. It doesn't matter what how $\tilde{d}_i$ is a function of the parameters.

I have read that for linear regression the Least Squares estimator is also the Maximum Likelihood Estimator under the assumption of normally distributed errors.

However is this also true for my nonlinear problem?

You are using

$$ \sqrt{(x_i - x_b)^2+(y_i - y_b)^2} = d_i +\epsilon_i$$

but it is possibly better to state it as the following instead

$$ \begin{array}{rcccl} d_i &=& \sqrt{(x_i - x_b)^2+(y_i - y_b)^2} &+& \epsilon_i \\ &=& \tilde{d}_i& + &\epsilon_i \end{array}$$

Then an observed distance $d_{i}$ is distributed according to the modeled mean distance $\tilde{d}_i$ plus some noise $\epsilon_i$.

Then the probability density (and the related likelihood) for the observations $d_i$ can be written as

$$\prod_{\forall i} \frac{1}{\sqrt{2\pi \sigma^2}} \exp \left( -\frac{(d_i-\tilde{d}_i)^2}{2\sigma^2} \right)$$

where the $\tilde{d}_i = \sqrt{(x_i - x_b)^2+(y_i - y_b)^2}$ can be modeled as functions of your unknown parameters $y_b$ and $x_b$.

Maximizing this likelihood is equivalent to minimizing the sum of squares $\sum_{\forall i} (d_i-\tilde{d}_i)^2$. It doesn't matter in what way $\tilde{d}_i$ is a function of the parameters.

The above is just to answer your direct question. Whether your model is sound (e.g. allowing negative distances) remains an open question.