We introduce and study the sequence of bivariate Generalized Bernstein operators $\{\mathbf{B}_{m,s}\}_{m,s}$, $m,s\in\mathbb{N}$, \[ \mathbf{B}_{m,s}=I-(I-\mathbf{B}_m)^s,\quad \mathbf{B}_m^i=\mathbf{B}_m(\mathbf{B}_m^{i-1}), \] where $\mathbf{B}_m$ is the bivariate Bernstein operator. These operators generalize the ones introduced and studied independently in the univariate case by Mastroianni and Occorsio [Rend. Accad. Sci. Fis. Mat. Napoli extbf{44} (4) (1977), 151--169] and by Micchelli [J. Approx. Theory extbf{8} (1973), 1--18] (see also Felbecker [Manuscripta Math. extbf{29} (1979), 229--246]). As well as in the one-dimesional case, for $m$ fixed the sequence $\{\mathbf{B}_{m,s}(f)\}_s$ can be successfully employed in order to approximate ``very smooth'' functions $f$ by reusing the same data points $f\big(\frac im,\frac jm\big)$, $i=0,1,\dots,m$, $j=0,1,\dots,m$, since the rate of convergence improves as $s$ increases. A stable and convergent cubature rule on the square $[0,1]^2$, based on the polynomials $\mathbf{B}_{m,s}(f)$ is constructed. Moreover, a Nyström method based on the above mentioned cubature rule is proposed for the numerical solution of two-dimensional Fredholm integral equations on $[0,1]^2$. The method is numerically stable, convergent and the involved linear systems are well conditioned. Some algorithm details are given in order to compute the entries of the linear systems with a reduced time complexity. Moreover the procedure can be significantly simplified in the case of equations having centrosymmetric kernels. Finally, some numerical examples are provided in order to illustrate the accuracy of the cubature formula and the computational efficiency of the Nyström method.