Abstract: Depending on the matrix exponential (ME) method, a compact and highly-accurate implementation of the perfectly matched layer (PML) with the complex-frequency-shifted (CFS) scheme based upon the bilinear Z-transform (BZT) technique is developed for modeling the finite-difference time-domain (FDTD) problems on the very low frequency (VLF) subsurface sensing and the very high frequency (VHF) imaging. By using ME method, the proposal becomes a compact first-order differential matrix form so that it can completely circumvent convolution operations, formula rearrangements, or variable replacements during mathematical derivations. Furthermore, the BZT-ME-PML can not only flexibly terminate arbitrary media because of adding the DB constitutive relation, but efficiently attenuate evanescent waves and decrease late-time reflections due to incorporating the CFS scheme. Besides, the BZT-ME-PML can maintain higher absorption performances when it comes closer to truncate the FDTD domain with targets so that smaller computational regions can be adopted for modeling 3D open-domain subsurface sensing and imaging cases, resulting in palpable improvements in memory requirement and CPU time. The numerical simulations on the subsurface sensing and imaging have been carried out to illustrate the validity and accuracy of the approach.