A new scheme is proposed for the inversion of surface waves using a continuous formulation of the inverse problem and the least squares criterion. Like some earlier schemes a Gaussian a priori covariance function controls the horizontal degree of smoothing in the inverted model, which minimizes some artifacts observed with spherical harmonic parameterizations. Unlike earlier schemes the new approach incorporates some sophisticated geometrical algorithms which dramatically increase computational efficiency and render possible the inversion of several tens of thousands of seismograms in few hours on a typical workstation. The new algorithm is also highly suited to parallelization which makes practical the inversion of data sets with more than 50,000 ray paths. The constraint on structural and anisotropic parameters is assessed using a new geometric approach based on Voronoi diagrams, polygonal cells covering the Earth's surface. The size of the Voronoi cells is used to give an indication of the length scale of the structures that can be resolved, while their shape provides information on the variation of azimuthal resolution. The efficiency of the scheme is illustrated with realistic uneven ray path configurations. A preliminary global tomographic model has been built for SV wave heterogeneities and azimuthal variations through the inversion of 24,124 fundamental and higher-mode Rayleigh waveforms. Our results suggest that the use of relatively short paths (<10,000 km) in a global inversion should minimize multipathing, or focusing/defocusing effects and provide lateral resolution of a few hundred kilometers across the globe.