We present an original algorithm for seismic imaging, based on the depth wavefield extrapolation by the one-way wave equation. Parallel implementation of the algorithm is based on the several levels of parallelism. The input data parallelism allows processing full coverage for some area (up to one square km); thus, data are divided into several subsets and each subset is processed by a single MPI process. The mathematical approach allows dealing with each frequency independently and treating solution layer-by-layer; thus, a set of 2D cross-sections instead of the initial 3D common-offset vector gathers are processed simultaneously. This part of the algorithm is implemented suing GPU. Next, each common-offset vector image can be stacked, processed and stored independently. As a result, we designed and implemented the parallel algorithm based on the use of CPU-GPU architecture which allows computing common-offset vector images using one-way wave equation-based amplitude preserving migration. The algorithm was used to compute seismic images from real seismic land data.