In this paper, three-dimensional (3D) subsurface imaging by inversion of data obtained from the very early time electromagnetic system (VETEM) is presented. The distorted Born iterative method is used to match the internal nonlinear property of the 3D inversion problem. The frequency explored is from 0.062 to 3.052 MHz, the high end of which is beyond the induction range. No reported 3D approximation works well for this frequency range, thus the 3D stabilized bi-conjugate-gradient fast Fourier transform method is invoked as the forward solver, in which the total current formulation other than the induced current formulation is used to achieve higher stability. With the use of the reciprocity theorem, the Fréchet derivatives are computed rapidly by two calls to the forward solver. Tikhonov regularization is used so that the contribution from small singular values is suppressed to warrant a smooth solution. The selection of regularization parameter follows a heuristic rule as used in the Levenberg-Marquardt algorithm so that the iteration is stable. Processing the whole data set at one time is prohibitively time-consuming and expensive. The localized property of the geophysical problem is exploited to decompose the original large problem into a number of smaller subproblems. The nonlinear inversion is implemented subsequently on each subproblem first. The inverted conductivity profiles are then combined together to arrive at a complete conductivity image under the whole data domain. The proposed method is first tested on synthetic data and then applied to the VETEM data from the Pit 9 complex of the Department of Energy.