The recovery of 3D volume of density from gravity field data is a key feature of geophysical and geological interpretations. Gravity field data inversion involves solving an underdetermined problem; therefore, large-scale data inversion is costly in time and memory consumption. Calculation efficiency is a primary concern for gravity field data inversion; therefore, multiple methods are considered and applied to increase the inversion efficiency. The solution for an inversion problem was formulated by incorporating constraints to obtain stable inversion results, and a new projected Barzilai-Borwein iterative algorithm was applied to accelerate convergence of the inversion method. To test the potential application of the projected Barzilai-Borwein iterative method, synthetic gravity data simulations and real data applications were carried out. Numerical performances and practical application indicate the fast convergence of projected Barzilai-Borwein iterative method and the increase in the calculated efficiency.