Computation of Common Middle Point seismic sections and their subsequent time migration and diffraction imaging provides very important knowledge about the internal structure of 3D heterogeneous geological media and are key elements for successive geological interpretation. Full-scale numerical simulation, that computes all single shot seismograms, provides a full understanding of how the features of the image reflect the properties of the subsurface prototype. Unfortunately, this kind of simulations of 3D seismic surveys for realistic geological media needs huge computer resources, especially for simulation of seismic waves’ propagation through multiscale media like cavernous fractured reservoirs. Really, we need to combine smooth overburden with microstructure of reservoirs, which forces us to use locally refined grids. However, to resolve realistic statements with huge multi-shot/multi-offset acquisitions it is still not enough to provide reasonable needs of computing resources. Therefore, we propose to model 3D Common Middle Point seismic cubes directly, rather than shot-by-shot simulation with subsequent stacking. To do that we modify the well-known “exploding reflectors principle” for 3D heterogeneous multiscale media by use of the finite-difference technique on the base of grids locally refined in time and space. We develop scalable parallel software, which needs reasonable computational costs to simulate realistic models and acquisition. Numerical results for simulation of Common Middle Points sections and their time migration are presented and discussed.