The GW150914 signal was observed, giving us the frequency and amplitude of the event. Because LIGO has two detectors a rough source location could be derived.
But how do these three factors allow for the mass of the black holes and their distance to be calculcated? If the wave strengths are in a square relationship to distance, then couldn't there be an infinite number of other masses and distances which would give the same signals?
The masses of the two binary objects are encoded in the frequency and frequency evolution of the gravitational waves. In the usual parameterization the two parameters that are most easily measured from the wave phase are the total mass $M=m_1+m_2$ and the "chirp mass": $$ \mathcal{M} = \frac{(m_1\, m_2)^{3/5}}{(m_1+m_2)^{1/5}} = \frac{c^3}{G}\left[\frac{5}{96}\,\pi^{-8/3}\,f^{-11/3}\,\dot{f}\right]^{3/5},$$
where $G$ and $c$ are Newton's gravitational constant and the speed of light; $f$ and $\dot{f}$ are the gravitational wave frequency and it's first derivative.
The total mass, the distance to the source, and the sources location in the sky encode in the amplitude of the waves. Once you've determined $M$ and $\mathcal{M}$ from the phase, you can use triangulation between multiple detectors to determine the sky location. Finally, with the sky location and total mass in hand you can determine the distance.
In practice all of these parameters (and several others) are fit simultaneously, and there are lots of correlations to deal with.
If you are really interested, check out LIGO document P1500218 "Properties of the binary black hole merger GW150914".
In short: because we measure both amplitude and phase.
In the amplitude $A$, distance and mass are degenerate, so you can only measure the following combination of them: $$ A \propto \frac{1}{r}\frac{(m_1m_2)^{1/2}}{(m_1+m_2)^{1/6}} $$
Meanwhile, the phase depends very sensitively on the masses of the objects, but not on distance. We can therefore constrain the masses in the expression above and break the mass-distance degeneracy to determine $r$.
The long(er) answer:
The phase and amplitude of GWs produced by compact binary mergers like GW150914 are impossible to model exactly and require numerical relativity simulations for a general solution. However, we can do a reasonably good job of approximating them in the weak-field regime, where the objects are moving sufficiently slower than light (and when we're sufficiently distant from them).
We do this by approximating their orbital dynamics with a post-Newtonian expansion in the small parameter $(v/c)^2$ (where $v$ is the orbital speed of the objects). To leading (Newtonian) order in this expansion (i.e., where $v \ll c$), the amplitude $h$ and phase $\psi$ of the gravitational waveform look like (in the frequency domain):
$$ h(f) = \frac{1}{r}\mathcal{M}^{5/6}f^{-7/6}\exp(i\psi(f)) $$ $$ \psi(f) = 2\pi f t_c - \phi_c - \frac{\pi}{4} + \frac{3}{128}(\pi\mathcal{M}f)^{-5/3} $$ where $t_c$ is the time at coalescence, $\phi_c$ is the phase at coalescence, and $\mathcal{M}$ is the chirp mass. This approximation can be improved by adding higher-order terms in $v/c$ (and there are in fact several different ways of extending the PN expansion beyond leading order).
Notice that the distance $r$ is absent from the phase $\psi$. Since the chirp mass can be determined independently (and very accurately) from the phase alone, the degeneracy in $h$ can be broken and, to the extent that the amplitude evolution can be measured (which isn't quite as well as we'd like), we can determine the distance $r$.
In practice, however, there's an additional degeneracy in amplitude with the sky location and orientation of the source binary. The strain amplitude $h$ above is approximately correct only for a face-on binary that is directly overhead a single detector; the detector's response function in fact depends on the location of the source (and its orientation), such that the amplitude for a less-than-optimally located binary will be less than this maximum.
The sky location can be determined crudely by timing triangulation, or with slightly better accuracy by incorporating phase differences between detector sites. Like Paul T says, this is most effectively done with a coherent Bayesian analysis of the detector data that fits all model parameters simultaneously (there are 15 of them).
Since the sky location is generally quite poorly measured (tens to hundreds of square degrees for typical signals), the resulting error on the distance measurement is also large: typically 10-30%.
The simple toy answer: You can extract it by maximizing the following quantity derived from Bayesian analysis:
$(s|h(\mathbf{\theta}))-(h(\mathbf{\theta})|h(\mathbf{\theta}))$, where $s$ is the signal $\theta$ are the parameters of this waveform (e.g. chirp mass, distance, etc), $(a|b)=\int \frac{a b^*}{S_n(f)} df$ with $S_n$ being the power spectral density from the LIGO detector (you can get the form from e.g. LIGO tutorial website) and $$ h(f) = \frac{1}{r}\mathcal{M}^{5/6}f^{-7/6}\exp(i\psi(f)) $$ $$ \psi(f) = 2\pi f t_c - \phi_c - \frac{\pi}{4} + \frac{3}{128}(\pi\mathcal{M}f)^{-5/3} $$ where $\mathcal{M}$ is the chirp mass, $\psi$ is the phase, $f$ frequency, $t_c$ time of coalescence, $\phi_c$ phase of coalescence To first order approximation
To answer your question on the degeracy of $M_1$ and $M_2$: It is much harder to find out the individual masses when there is noise than it is to find e.g. the chirp mass $\mathcal{M}$. This is for the exact reason that these parameters are somewhat degenerate with other parameters. Also for this reason when you look at the LIGO detection paper, these parameters have larger error bars than the chirp mass. Looking at the formula for the gravitational wave strain, amplitude and phase parameters should be possible to determine independently. Once you include both the sky localization and higher order terms to $h(f)$, you would see that there is no longer strict degeneracy between the chirp mass and individual masses.
Notes to avoid confusing
