There are many possible ways to explain the why, but I'd like to try one.
Start with some basic definitions for the Poincaré disk model. A point of the hyperbolic plane is a point inside the unit disk. A hyperbolic line is a circle arc which is perpendicular to the unit disk. From that one can eventually deduce that the first four Eulidean axiom as well as the hyperbolic version of the axiom of parallels do hold. So you have a model of hyperbolic geometry.
What you need to study next is the set of isometries of this model. An isometry has to be a bijection between hyperbolic points, so it maps the interior of the unit disk onto itself. Furthermore, an isometry will preserve angles, and since the Poincaré models are conformal, that means Euclidean angles will be preserved.
The set of angle-preserving (conformal) transformations of the plane (or rather the Riemann sphere, but don't worry about that distinction just now) is the set of all Möbius transformations. And one of the most natural ways to represent Möbius transformations is by interpreting points in the plane as complex numbers. (Actually I consider the orojective $\mathbb{CP}^2$ view even more natural, but that's again besides the point.) This is where complex numbers enter the scene. So taken together with mapping the unit disk onto itself, this means that you are after the set of all Möbius transformations which fix the unit disk.
Now one can write down a formula for these, and then do distance comparisons. Suppose you had a hyperbolic line segment from some point $P$ to some nearby point $Q$. You could move that point $P$ into the center, using an isometry which would preserve hyperbolic distances but change Euclidean ones. Applying that transformation to $Q$ would move it close to the center, and then you could measure distances close to the center. You could move your point $Q$ infinitesimally close to $P$, and do the same thing again. What you would get is the distance metric you mentioned, up to a scale factor.
A finite distance is simply an integral over infinitely many infinitesimal distances, so you'd have to integrate from the start point ($z_1$ in your question) to the end point ($z_2$) along a geodesic path. Integrating along geodesics might be tough, though. So you can avoid that by first applying a transformation which moves one of them to the center of the unit disk.
$$
f(z)=\frac{z-z_1}{-\overline{z_1}z+1}\qquad
d\bigl(z_1,z_2\bigr)=d\bigl(0,f(z_2)\bigr)
=2\tanh^{-1}\bigl\lvert f(z_2)\bigr\rvert
=2\tanh^{-1}\bigl\lvert\frac{z_2-z_1}{1-\overline{z_1}z_2}\bigr\rvert
$$
This is the formula you quoted, except for a sign change inside the absolute value which is of no relevance. One can show that all Möbius transformations which fix the unit disk have the form
$$z\mapsto\frac{az+b}{\bar bz+\bar a}$$
so that when you know the numerator for $f$ (which should be zero for $z=z_1$) then you already know the denominator as well. You might choose $a\neq 1$, but the result would simply be the same as above followed by a rotation around the origin.
So now you'd be able to compare distances, like “that hyperbolic segment is twice as long as that other one”. The constant term in that formula is a matter of convention. It is usually chosen in such a way that the Gaussian curvature of the plane is not only constant but actually $-1$. See also this question for details on that convention.
With respect to your question about alternatives, I'd like to point out that cross ratios can be a useful tool there. Take two points in the disk and connect them by a geodesic. That geodesic will intersect the unit circle in two ideal points. Treat all four points as complex numbers and compute their cross ratio. It will not change under Möbius transformations, so it can be used to define distances. Sicne cross ratios are multiplicative and lengths are additive, you need to take the logarithm. Then choose the correct constant as above. The question I referenced there uses this cross ratio approach, so it may be useful here as well.
Note that my above post so far deliberately ignored the fact that Anti-Möbius transformations which fix the unit disk are hyperbolic isometries as well. This is however irrelevant to the question at hand. You should only keep it in mind when you build up your intuition about the set of all hyperbolic isometries.
Best Answer
The proof is relatively easy using the upper plane model. I will sketch it as it reduces to simple euclidean geometry proofs. Take a Mobius transform of the disc to the upper plane that sends $a$ to $0$ and $b$ to $\infty$; then the two circular arcs become lines through the origin $L,M$ say.
The perpendicular from a point $P \in L$ is (contained in) the Euclidean circle centered on the real axis whose tangent at the intersection $Q$ with $M$ is perpendicular on $M$. This immediately implies that the center is the origin and the distance in cause is the hyperbolic distance between $PQ$ where $P, Q$ are on the same circle with the center the origin and we need to show that is constant regardless of $P \in L$.
But it is a well-known result that if we let $A_p, B_P$ the intersection of the circle with the real axis, the hyperbolic distance from $P$ to $Q$ is $|\log \frac{|PA_P||QB_P|}{|PB_P||QA_P|}|$ (where $|PA_P|$ is just the usual Euclidean distance)
However it is obvious that $PA_P || P'A_{P'}$ for any $P, P' \in L$ as the triangles $OPA_P, OP'A_{P'}$ are isosceles at the origin ($|OP|=|OA_P|$ etc) and hence similar, so all those products ratios are constant in $P,Q$ as required (as they are just $r/r'$ or $r'/r$) and we are done!