## 1 Introduction
divided the TEM diffusion field into segments based on a time sequence and applied a regularization algorithm to implement the reverse transformation of the wavefield [27]. However, this approach introduced issues with the continuity between segmented data. To address the continuity issue of the virtual wavefield, [PERSON] et al. reduced the condition number of the ill-conditioned matrix using an over-relaxation preconditioned method and solved the virtual wavefield over the full time period with a regularized conjugate gradient algorithm [28]. Nonetheless, this method exhibited low numerical accuracy and poor noise resistance. Building on this, [PERSON] et al. later employed a multi-window time-sweeping technique to suppress background noise interference in the virtual wavefield to some extent [12]. In addition, [PERSON] et al. combined regularization with singular value decomposition and introduced weighting factors to obtain a stable virtual wavefield [29]. In order to further improve the numerical accuracy of the virtual wave field, some scholars have proposed a preconditioned precise integration algorithm [23; 24], which transforms the solution of the ill-conditioned linear system into a stable process of infinite integration of the matrix exponential function, ensuring high computational efficiency while maintaining accuracy. Furthermore, by incorporating the multi-window time-sweeping technique into the preconditioned precise integration algorithm [23], the noise resistance of the virtual wavefield was significantly enhanced. However, the coherent axes of the virtual wavefield typically reflect the characteristics of shallow formations more accurately. As subsurface structures become increasingly complex, the coherent axes of the virtual wavefield in deeper regions cannot reliably represent the stratigraphic morphology, necessitating the use of migration imaging methods for spatial positioning [30]. [PERSON] et al. used the three-dimensional boundary element method combined with the Kirchhoff method to achieve the surface extension of the virtual wavefield [31]. [PERSON] et al. employed deconvolution techniques and wavelet compression to enhance the vertical resolution of migration imaging [32]. [PERSON] et al. based on a full time-domain virtual wavefield, utilized Kirchhoff integral migration, inverse synthetic aperture, and deconvolution techniques to achieve high-resolution pseudo-seismic three-dimensional imaging results [12]. [PERSON] et al., drawing on the Born approximation imaging method from the field of seismic exploration, achieved rapid imaging of the virtual wavefield [22]. By processing electromagnetic data from coal seam goat areas, they successfully identified electrical interface features. [PERSON] et al. used virtual wavefields to achieve full-waveform inversion imaging, which not only improves vertical resolution but also provides accurate parameter estimation [26]. To better address the migration imaging problem in media with significant lateral velocity variations, [PERSON] et al., starting from the time-domain wave equation, extrapolated the first-order approximation of the upgoing wave equation to obtain a time-space domain wavefield recursion formula [33]. Using the finite difference method, they solved the wavefield in subsurface spaces and obtained accurate electrical interface features. While all these imaging methods provide information on electrical interfaces within formations, they are limited by the resolution of the source field and struggle to achieve multi-resolution imaging from shallow to deep regions. In summary, this paper first transmits multiple differential pulses with varying pulse widths during the ON-Time period to obtain a multi-resolution TEM diffusion field. Then, a preconditioned precise integration wavefield reverse transformation method is employed to extract the virtual wavefield from the diffusion field.