Reverse time migration (RTM) has shown increasing advantages in handling seismic images of complex subsurface media, but it has not been used widely in 3D seismic data due to the large storage and computation requirements. Our prime objective was to develop an RTM strategy that was applicable to 3D vertical seismic profiling (VSP) data. The strategy consists of two aspects: storage saving and calculation acceleration. First, we determined the use of the random boundary condition (RBC) to save the storage in wavefield simulation. An absorbing boundary such as the perfect matching layer boundary is often used in RTM, but it has a high memory demand for storing the source wavefield. RBC is a nonabsorbing boundary and only stores the source wavefield at the two maximum time steps, then repropagates the source wavefield backwards at every time step, and hence, it significantly reduces the memory requirement. Second, we examined the use of the graphic processing unit (GPU) parallelization technique to accelerate the computation. RBC needs to simulate the source wavefield twice and doubles the computation. Thus, it is very necessary to realize the RTM algorithm by GPU, especially for a 3D VSP data set. GPU and central processing unit (CPU) collaborated parallel implementation can greatly reduce the computation time, where the CPU performs serial code, and the GPU performs parallel code. Because RBC does not need the same huge amount of storage as an absorbing boundary, RTM becomes practically applicable for 3D VSP imaging.