I suspect you have a typo, and the denominator should be $|1-z_i \overline{z_j}|^2$. At least, I will show that the quantity is Möbius invariant with this modification. It is possible, of course, that both are invariant.
The determinant is a red herring; each individual entry in the matrix is invariant under Möbius transformations. I begin by quoting a fact from hyperbolic geometry: Define $\delta(u,v) = \frac{|u-v|^2}{(1-|u|^2)(1-|v|^2)}$ for $u$ and $v$ in $\mathbb{D}$. Then this quantity is invariant under Möbius transformations. I can't give an intuition for this fact, but its easy to give a reference.
So $1+\delta(u,v) = 1+\frac{(u-v)(\overline{u} - \overline{v})}{(1-u \overline{u}) (1-v \overline{v})} = \frac{\left( 1-u \overline{u} - v \overline{v} + u \overline{u} v \overline{v} \right) + \left( u \overline{u} - u \overline{v} - \overline{u} v + v \overline{v} \right)}{(1-u \overline{u}) (1-v \overline{v})}$ $=\frac{1- u \overline{v} - \overline{u} v + u \overline{u} v \overline{v}}{(1-u \overline{u}) (1-v \overline{v})} = \frac{(1-u \overline{v})(1-\overline{u} v)}{(1-|u|^2)(1-|v|^2)} = \frac{|1-u \overline{v}|^2}{(1-|u|^2)(1-|v|^2)}$ is invariant under Möbius tranformations.
Replacing $u$ and $v$ by $z_i$ and $z_j$, this is the reciprocal of the $(i,j)$ entry in your matrix. So every element of your matrix is Möbius invariant, as claimed.
UPDATE: I just realized a nice way to express anon's solution below. Let $\sigma$ denote Schwarz reflection in the boundary of $\mathbb{D}$. Since Schwarz reflection is a conformally invariant operation, if $\phi$ is a Möbius transformation of $\mathbb{P}^1$ preserving $\mathbb{D}$, then $\phi(\sigma(z)) = \sigma(\phi(z))$. Explicitly, $\sigma(z) = \overline{z}^{-1}$.
Then the $(i,j)$ entry in your matrix is the cross ratio $(z_i, z_j; \sigma(z_i), \sigma(z_j))$, by a computation very similar to anon's.