This question How to work out orthogonal polynomials for regression model and the answer https://math.stackexchange.com/a/354807/51020 explain how to build orthogonal polynomials for regression. However they only consider one dimensional functions.
How can we use (discrete) orthogonal polynomials for regression with two dimensional functions (i.e., $z = f(x, y)$)?
In three dimensions, we have the spherical harmonics, a polynomial set orthogonal over the unit sphere.
http://functions.wolfram.com/Polynomials/SphericalHarmonicY/
http://mathworld.wolfram.com/SphericalHarmonic.html
https://en.wikipedia.org/wiki/Spherical_harmonics
In two dimensions, we have the Zernike polynomials, orthogonal over the unit disk.
http://functions.wolfram.com/Polynomials/ZernikeR/
http://mathworld.wolfram.com/ZernikePolynomial.html
https://en.wikipedia.org/wiki/Zernike_polynomials
To answer your question via demonstration, construct the radial Zernike functions by applying the orthogonalizing method of Gram and Schmidt to the functions sequence $$ \left\{ r^{2k} \right\}_{k=0}^{\infty} = \left\{ 1, r^{2}, r^{4}, \dots \right\} $$ over the interval $0\le r < 1$ with monic normalization $(R(r)=1)$. The result is the set of rotationally invariant Zernike polynomials $$ \left\{ R^{0}_{2k}(r) \right\}_{k=0}^{\infty} = \left\{ 1, 2r^{2}-1, 6r^{4} - 6r +1, \dots \right\} $$ In general, $$ R^{0}_{2k}(r) = \sum_{j=0}^{k} \left(-1\right)^{j} \frac{\left( 2k-j\right)!} {j!\left(\left( k-j\right)!\right)^{2}}r^{2(k-j)} $$
The following is a cutaway diagram of the first eight terms. Notice the Fourier like quality: a constant offset term, and then functions of zero mean, with increasing frequency.
We can expand these polynomials to account for angular variation too. The full set of Zernike polynomials follows.
Start with non-negative integers $n$ for order and $m$ for angular frequency such that $n-m$ is even. The recursion relationship is simplified by defining the midpoint $\omega=\frac{1}{2}(n-m)$ and the average $\sigma=\frac{1}{2}(n+m)$. The full set of radial polynomials is: $$ R_n^m(r) = \sum_{j=0}^{\omega}{ \left(-1\right)^j \frac{(n-j)!} {j!\left(\omega-j\right)!\left(\sigma-j\right)!} r^{n-2j}} $$ The Zernike polynomials are a product of the radial functions and the exponential term $$ U_n^m(r,\theta) = R_n^{m}(r)e^{i m \theta} $$ The set complete through order $n=6$ is shown below. The order increase $n=0, 1, 2, dots$ moving down. The angular frequency increases as you move aways from the centerline.
There are a few ways to define the unit disk, $$ D_{2} = \begin{cases} \left\{z\in\mathbb{C} \colon |z| \le 1\right\} \\ \left\{\left( r, \theta \right)\in\mathbb{R}^{2} \colon 0\le r \le 1, 0 \le \theta < 2\pi\right\} \\ \left\{\left( x, y \right)\in\mathbb{R}^{2} \colon x^{2} + y^{2} \le 1\right\} \end{cases} $$ and the first few polynomials are expressed in terms of all three coordinates.
$$ \begin{array}{cclllll} % n & m & \text{Complex} && \text{Polar} && \text{Cartesian} \\\hline % 0 & 0 & 1 & Re: & 1 & Re: & 1 \\[2pt] & & & Im: & 0 & Im: & 0 \\[4pt] % 1 & 1 & z & Re: & r\cos\theta & Re: & x \\[2pt] & & & Im: & r\sin\theta & Im: & y \\[4pt] % 2 & 0 & 2|z|^2-1 & Re: & 2r^{2}-1 & Re: & 2\left(x^{2} + y^{2}\right)-1 \\[2pt] & & & Im: & 0 & Im: & 0 \\[2pt] 2 & 2 & z^{2} & Re: & r^{2}\cos\theta & Re: & (x-y) (x+y) \\[2pt] & & & Im: & r^{2}\sin\theta & Im: & 2xy \\[4pt] % 3 & 1 & z\left(3|z|^2-2\right) & Re: & \left(3r^{2}-2r\right)\cos\theta & Re: & x \left(3 \left(x^2 + y^2\right) - 2\right) \\[2pt] & & & Im: & \left(3r^{2}-2r\right)\sin\theta & Im: & y \left(3 \left(x^2 + y^2\right) - 2\right) \\[2pt] 3 & 3 & z^{3} & Re: & r^{3}\sin\theta & Re: & \phantom{-}x^{3} - 3xy^{2} \\[2pt] & & & Im: & r^{3}\sin\theta & Im: & -y^{3} + 3yx^{2} \\[4pt] \end{array} $$