Basics of wave-mechanical model of atom
These orbitals represent the angular part of the wavefunction. The solution obtained directly from solving the Schrödinger equation produces equations containing complex numbers so cannot be drawn on normal $xyz$ axes and are hard to visualise.
The angular part of the wavefunction is given by functions called Spherical Harmonics these usually are given the symbol $Y$ and require two quantum numbers $\ell$ and $m_z$ making $Y_\ell^{m_z}$.
What is done to make the orbitals real functions, as opposed to complex ones, is to take linear combinations of the spherical harmonics, thus choosing the $z$-axis then the $\mathrm{d}_{z^2} = Y_2^0 = N(3\cos^2(\theta)-1)$) where $N$ is the normalisation constant. In cartesian coordinates ($xyz$) rather than polar ones (which are in $r,\;\theta,\;\phi$), $\displaystyle Y^0_2= N\frac{3z^2-r^2}{r^2}$ so you can see where the $z^2$ comes from and $r$ is the radial distance.
The other orbitals are also linear combinations e.g. $\mathrm{d}_{xz}=(Y_2^{1}+Y_2^{-1})/\sqrt{2}$ and in cartesian coordinates $\displaystyle Y_2^1 \sim \frac{z(x+iy)}{r^2},\; Y_2^{-1}\sim \frac{z(x-iy)}{r^2}$ so you can see where the notation $zx$ comes from and by making a linear combination the sum is no longer a complex number ($i=\sqrt{-1}$).
For completeness $\mathrm{d}_{yz}=(Y_2^{1}-Y_2^{-1})/\sqrt{2},\; \mathrm{d}_{x^2-y^2}=(Y_2^{2}+Y_2^{-2})/\sqrt{2},\; \mathrm{d}_{xy}=(Y_2^{2}-Y_2^{-2})/\sqrt{2}$  and in cartesians $\displaystyle Y_2^{\pm2}\sim \frac{(x\pm iy)^2}{r^2}$.