Abstract: 3D eikonal equation is a partial differential equation for the calculation of first-arrival traveltimes and has been widely applied in many scopes such as ray tracing, source localization, reflection migration, seismic monitoring and tomographic imaging. In recent years, many advanced methods have been developed to solve the 3D eikonal equation in heterogeneous media. However, there are still challenges for the stable and accurate calculation of first-arrival traveltimes in 3D strongly inhomogeneous media. In this paper, we propose an adaptive finite-difference (AFD) method to numerically solve the 3D eikonal equation. The novel method makes full use of the advantages of different local operators characterizing different seismic wave types to calculate factors and traveltimes, and then the most accurate factor and traveltime are adaptively selected for the convergent updating based on the Fermat principle. Combined with global fast sweeping describing seismic waves propagating along eight directions in 3D media, our novel method can achieve the robust calculation of first-arrival traveltimes with high precision at grid points either near source point or far away from source point even in a velocity model with large and sharp contrasts. Several numerical examples show the good performance of the AFD method, which will be beneficial to many scientific applications.