In short: because we measure both amplitude and phase.
In the amplitude $A$, distance and mass are degenerate, but aso 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 measurementdepends very sensitively on the masses of the objects, but (whichnot on distance. We can determinetherefore constrain the masses in the expression above and break the mass very precisely) breaks this-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%.