In optical testing, such as the fringe reflection technology and Shack-Hartman wavefront sensor technology, slope of a surface is measured instead, from which the faithful surface of the test optic is obtained. Therefore, a gradient data-based wavefront reconstruction is needed. This paper shows the use of the Gram-Schmidt process for orthonormalizing the gradients of the two-dimensional Legendre polynomials. After a set of orthonormalized vector polynomials is generated in a square region, these polynomials can be used to fit the gradient data in the region. By a simple linear transformation, the fitting coefficients can be derived and transformed to the wavefront description of the two-dimension Legendre polynomials, and the wavefront and primary aberration are then obtained. Based on the zonal method, this can effectively reconstruct the high-frequency component by fitting the difference of the high-frequency error, which cannot be done by polynomial fitting. According to the computer simulation, this algorithm can primely realize the reconstruction of wavefront.