Claim: Suppose that $f$ is a conformal map from the unit disk to the unit disk, $f(z_1) = w_1$ , $f(z_2) = w_2$, then we it can be derived that $|\frac{z_1-z_2}{1-\bar z_1 z_2}| = |\frac{w_1-w_2}{1-\bar w_1 w_2}|$. Let $z_1$ approaches $z_2$, we have $\frac{|dz|}{1-|z|^2}= \frac{|dw|}{1-|w|^2}$, so can conclude that the Riemannian metric $ds= \frac{2|dz|}{1-|z|^2}$ is invariant under the conformal mapping.
Questions:
What does $dz$ actually mean here? Is it some sort of complex form? How should we conclude that $\frac{|dz|}{1-|z|^2}= \frac{|dw|}{1-|w|^2}$ by letting $z_1$ tends to $z_2$.
Is $|dz| = (dx)^2+(dy)^2$ which is the usual Riemannian metric. I feel that to understand the claim, one needs some rigorous treatment of the complex form maybe? Is there any reference on this matter?
More details of the claim can be found in Ahlfors, 'Conformal Invariants, Topics in Geometric Function Theory' page 2.
Questions (1) and (2) have been addressed in the comments, but I am still confused about the procedure of letting '$z_1$ approach $z_2$', since it seems to be a rather vague idea, how do you really formalize this idea to arrive the conclusion such that the given metric is invariant under the conformal mapping?
The $dz$ is the canonical one-form of the complex plane: $dz=dx+idy$, if complex numbers are written as $x+iy$ for $x,y\in\mathbb R$.
This is the one-form used in complex integration. If you multiply it by a complex function $f(z)$, you get a new one-form $f(z)dz$. You can integrate this one-form over any $C^1$ curve $\gamma$ to get $\int_\gamma f(z)dz$. This will precisely match the definition of the complex integral given in any complex analysis course.
Treating $dz$ as a complex one-form is by no means necessary for complex integration. The same holds for your problem: you can formalize it using complex one-forms, but the object $dz$ is typically left undefined one relies on more informal geometric reasoning. The symmetrized tensor product of $dz=dx+idy$ and its complex conjugate $\overline{dz}=dx-idy$ is indeed $dx\otimes dx+dy\otimes dy$, the usual Riemannian metric.
More details on the invariance of the metric
Here is some elaboration on what the hyperbolic metric is in terms of Riemannian geometry and why it is preserved. You can also read from here why $|dz|$ is related to $|dw|$ as claimed. If you want more details on how to relate that claim to the one about preserved hyperbolic metric, I suggest you ask a follow-up question.
Let me denote $c(z)=\frac12(1-|z|^2)$ and let $e$ be the usual Euclidean metric (as a Riemannian metric). Then the claim is that the conformally Euclidean metric $g(z)=c^{-2}(z)e(z)$ is invariant under conformal mappings.
That is, if $f\colon D\to D$ is conformal, then $f^*g=g$. It is enough to show that for any $z\in D$ and any $v\in\mathbb R^2(=T_zD)$ we have $$ g_z(v,v) = g_{f(z)}(df_zv,df_zv). $$ This can be rephrased in terms of the Euclidean metric for which $e_z(v,v)=|v|^2$ and the conformal factor $c$: $$ c^{-2}(z)|v|^2 = c^{-2}(f(z))|df_zv|^2. $$
Now take any $C^1$ curve $\gamma$ with $\gamma(0)=x$ and $\dot\gamma(0)=v$. Then $$ |v|=\lim_{t\to0}\left|\frac{\gamma(t)-z}{t}\right| $$ and $$ |df_zv|=\lim_{t\to0}\left|\frac{f(\gamma(t))-f(z)}{t}\right|. $$ Using $z$ and $\gamma(t)$ for $z_1$ and $z_2$ in your lemma (namely $|\frac{z_1-z_2}{1-\bar z_1 z_2}| = |\frac{w_1-w_2}{1-\bar w_1 w_2}|$) gives $$ |df_zv|=\lim_{t\to0} \left|\frac{1-f(\gamma(t))\overline{f(z)}}{1-\gamma(t)\overline{z}}\right| \cdot \left|\frac{\gamma(t)-z}{t}\right|, $$ whence $$ |df_zv| = \frac{c(f(z))}{c(z)} |v|. $$ This is what we were after.
If the notation is not clear, please let me know.