The best derivation is Polyakov's, and it is found in the long string chapter of "Gauge Fields and Strings".
The key point is that the h-field in the path integral is integrated over, but it doesn't have derivative terms, so the fluctuations in the h-field just act to replace it at each point by its stationary value. The X-parts just go along for the ride when looking for stationary points of h, so you can write the action as
$$ S = \int \sqrt{h} h^{\alpha\beta} \gamma_{\alpha\beta} $$
Where $\gamma_{\alpha\beta} = \partial_\alpha X^\mu \partial_\beta X_\mu$ is the dot product of an $\alpha$ coordinate step with a $\beta$ coordinate step, i.e. it is the induced metric. The induced metric plays the role of a source term in the h path-integral (ignoring the X path integral). The stationary point condition is found by varying h (using the important determinant variation formula $\delta h = h h^{\alpha\beta} \delta h_{\alpha\beta}$ which you learn in math class as "expansion by minors" and "the inverse-minor theorem"):
$$\sqrt{h} \gamma_{\alpha\beta} + {1\over2\sqrt{h}} h h_{\alpha\beta} h^{\kappa\delta}\gamma_{\kappa\delta} $$
If you solve for h, you find that
$$ h_{\alpha\beta} = - {\gamma_{\alpha\beta}\over {1\over 2} h^{\kappa\delta}\gamma_{\kappa\delta}}$$
This might look like an incomplete solution, but the denominator on the right is a scalar, so this is saying that the tensors h and $\gamma$ are proportional
$$ h_{\alpha\beta} = A(x) \gamma_{\alpha\beta} $$
Where the proportionality constant A(x) won't make any difference (any two A choices will give solutions, and they lead to the same action).
Substitute in the extremal value for h in the action, and remember how to take an inverse matrix: $h^{\alpha\beta} = {1\over A} \gamma^{\alpha\beta}$, and you get that the action contribution for each external source $\gamma_{\alpha\beta}$ is proportional to $\sqrt{\gamma}$ no matter what $A(x)$ happens to be, which gives the Nambu-Goto action. Then you integrate the Nambu-Goto action over the remaining path-integral variables, which are the embedding coordiantes $X^\mu$.
The Nambu-Goto path integral is hard to understand in any way other than solving it classically, defining harmonic oscillators, and quantizing these by assuming they turn into standard harmonic oscillators. This is the old approach to string theory. The Polyakov action is just used to fix a gauge for h which will turn the problem into a simple sigma-model. So the equivalence between them is more of a formal thing, which relates the harmonic oscillator expansion to the vertex operators in the h formalism. It isn't necessarily a path-integral equality, because the Nambu-Goto path integral is not clearly well defined outside of turning it into Polyakov and fixing gauge for h.
There are two manifolds that are involved in string propagation.
The spacetime in which the string propagates.
The worldsheet of the string itself.
The fields $X^\mu$ are embedding coordinates of the worldsheet in the spacetime manifold. This means that for each point $(\sigma^1, \sigma^1)$ on the worldsheet, $X^\mu(\sigma^1, \sigma^2)$ gives the coordinates of that point in the spacetime manifold.
In the case you are considering, the spacetime is taken to be Minkowski, so the metric is $\eta_{\mu\nu}$. Now we could ask
"Given that the worldsheet is a two dimensional embedded submanifold of Minkowski space, is there some way that this manifold inherits its metric from the metric on the ambient spacetime?"
This question is analogous to
"Given that the sphere $S^2$ is some two-dimensional embedded submanifold of Euclidean space $\mathbb R^3$, is there some natural sense in which it inherits its metric from $\mathbb R^3$?
The answer to both of these question is yes, and the metric on the submanifold that does this is precisely the induced metric. The formula expression the induced metric for a two-dimensional submanifold of some ambient manifold with metric $g_{\mu\nu}$ (not necessarily flat) in terms of embedding coordinates is
$$
\gamma_{ab}(\sigma) = g_{\mu\nu}(X(\sigma))\partial_aX^\mu(\sigma)\partial_b X^\nu(\sigma), \qquad \sigma = (\sigma^2, \sigma^2)
$$
You are right about the derivation of the induced metric, it comes from demanding that the distance measured between points on the embedded submanifold is calculated to be the same number whether you use the ambient metric, or the induced metric. To see that the above expression for the induced metric does this, simply note that the infinitesimal distance between any two points on the embedded submanifold can be written in terms of the ambient metric and the embedding coordinates as
\begin{align}
g_{\mu\nu}(X(\sigma))d(X^\mu(\sigma))d(X^\nu(\sigma))
&= g_{\mu\nu}(X(\sigma))\partial_a X^\mu(\sigma)\partial_bX^\nu(\sigma)d\sigma^ad\sigma^b \\
&= \gamma_{a b}(\sigma)d\sigma^ad\sigma^b
\end{align}
To get some intuition for all of this, recall that expression for embedding coordinates of $S^2$ in $\mathbb R^3$ is
\begin{align}
X(\theta, \phi) &= \sin\theta\cos\phi\\
Y(\theta, \phi) &= \sin\theta\sin\phi\\
Z(\theta, \phi) &= \cos\theta
\end{align}
and using these embeddings you should be able to show that the metric on the sphere is simply
$$
\gamma_{ab}(\theta, \phi) = \mathrm{diag}(1, \sin^2\theta)
$$
Let me know if that's unclear or if you need more detail!
Best Answer
$$ \gamma_{ab} = \partial_a X^\mu \partial_b X_\mu \implies \gamma^{ab} \delta \gamma_{ab} = 2 \gamma^{ab} \partial_a X^\mu \partial_b \delta X_\mu. $$ Then, \begin{align} \delta S &= - \frac{\mu}{2} \int d^2 \zeta \sqrt{-\gamma}2 \gamma^{ab} \partial_a X^\mu \partial_b \delta X_\mu \\ &= \mu \int d^2 \zeta \partial_b [ \sqrt{-\gamma} \gamma^{ab} \partial_a X^\mu ] \delta X_\mu + \text{bdy term} \end{align} The equation of motion is then $$ \partial_b [ \sqrt{-\gamma} \gamma^{ab} \partial_a X^\mu ] = 0 \implies \nabla^2 X^\mu = 0. $$
Variations
OP asked me to clarify how the variation of $\gamma$ was derived. Think of $\gamma_{ab}$ as a function of $X$, $$ \gamma_{ab}(X) = \partial_a X^\mu \partial_b X_\mu $$ Then, the variation is the simple statement $$ \delta \gamma_{ab}(X) = \gamma_{ab} ( X + \delta X ) - \gamma_{ab} (X) $$ where on the RHS, only the linear term in $\delta X$ is kept.
In our case, we have \begin{align} \delta \gamma_{ab}(X) &= \partial_a [ X^\mu + \delta X^\mu ] \partial_b [ X_\mu + \delta X^\mu ] - \partial_a X^\mu \partial_b X_\mu \\ &= \partial_a X^\mu \partial_b \delta X_\mu + \partial_a \delta X^\mu \partial_b X_\mu + O( (\delta X)^2 ) . \end{align} It follows that $$ \gamma^{ab} \delta \gamma_{ab} = 2 \gamma^{ab} \partial_a X^\mu \partial_b \delta X_\mu $$