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. In order to significantly reduce the query of computer resources we propose to model these 3D seismic cubes directly rather than the shot-by-shot simulation with subsequent CMP 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. To be able to simulate realistic models and acquisition we develop scalable parallel software, which needs reasonable computational costs. Numerical results for simulation of Common Middle Points sections and their time migration are presented and discussed.