Stokes flows with near-touching rigid particles induce near-singular lubrication forces under relative motion, making their accurate numerical treatment challenging. With the aim of controlling the accuracy with a computationally cheap method, we present a new technique that combines the method of fundamental solutions (MFS) with the method of images. For rigid spheres, we propose to represent the flow using Stokeslet proxy sources on interior spheres, augmented by lines of image sources adapted to each near-contact to resolve lubrication. Source strengths are found by a least-squares solve at contact-adapted boundary collocation nodes. We include extensive numerical tests, and validate against reference solutions from a well-resolved boundary integral formulation. With less than 60 additional image sources per particle per contact, we show controlled uniform accuracy to three relative digits in surface velocities, and up to five digits in particle forces and torques, for all separations down to a thousandth of the radius. In the special case of flows around fixed particles, the proxy sphere alone gives controlled accuracy. A one-body preconditioning strategy allows acceleration with the fast multipole method, hence close to linear scaling in the number of particles. This is demonstrated by solving problems of up to 2000 spheres on a workstation using only 700 proxy sources per particle.