A variant of the global conjugate gradient method for solving general linear systems with multiple right-hand sides is proposed. This method is called as the global conjugate gradient linear least squares (Gl-CGLS) method since it is based on the conjugate gradient least squares method(CGLS). We present how this method can be implemented for the image deblurring problems with Neumann boundary conditions. Numerical experiments are tested on some blurred images for the purpose of comparing the computational efficiencies of Gl-CGLS with CGLS and Gl-LSQR. The results show that Gl-CGLS method is numerically more efficient than others for the ill-posed problems.