Wave aberration of optical system describes the deviation of the exiting wavefront from perfect spherical wavefront. For field imaging, wave aberration of each field is generally different due to the existence of field-dependent aberrations induced by imperfect optical design and alignment. Wave aberration for different field is required to fully characterize the optical system. Wavefront interferometry and Shack-Hartmann wavefront sensing method are two traditional methods to measure the wave aberration. Compared to these two methods, image-based wavefront retrieval method determines wave aberration directly from image formed by the optical system with no extra hardware. This method has attracted more and more attentions for its unique advantages, and has been successfully applied in many domains such as large telescope alignment, X-ray imaging and so on. In this paper, we propose a spot image based multi-field wavefront retrieval method. As we know, wave aberration function of optical system is a function of the pupil and the field coordinates. Double Zernike expansion, which performs Zernike polynomials expansion on the pupil plane and image plane respectively, is used to fulfill field-dependent wave aberration expansion. With the coefficients of this expansion, wave aberration of every field can be obtained. The objective function and its gradient function for the multi-field wavefront retrieval are given and conjugate gradient algorithm is used to retrieve the coefficients of field-dependent wavefront expansion with the spot intensity images at different fields. Numerical simulation of this wavefront retrieval method is completed with MATLAB. The retrieval results are analyzed for both ideal and noisy conditions with different number of field points.