Tomostatics is a technology to estimate a source of error in seismic imageries called near-surface statics, which are spurious timing variations of seismic reflections caused by the lateral variations in topography and the thickness of the weathering layer of extremely low seismic velocities. Most tomostatics methods use first arrivals of seismic data to invert for a near-surface velocity model and estimate the static corrections based on the model. In this paper we present a 3D tomostatics method based on deformable layer tomography (DLT) using first arrival data. Because the near-surface velocity fields in many areas are composed of thickness-varying velocity layers, the DLT method focus on inverting for the geometry of layer interfaces. To account for inhomogeneous distribution of ray coverage, a multi-scale inversion is applied to invert for interface variations of different wavelengths simultaneously. Several measures are devised to establish the reference velocity models and to increase the computation efficiency. Tests with synthetic and field data demonstrated that, in areas of sufficient spatial coverage of shots and receivers, the first-arrival multi-scale tomostatics method is able to establish 3D near-surface velocity models of satisfactory levels of data fitness and model resolvability.