Three-dimensional wave equation prestack depth imaging is an important tool in reconstructing images of complex subsurface structures, and it has become a technique gaining wide popularity in oil and gas industry. This is a large-scale scientific computing problem and can be considered a process of data continuation downward with the surface data or the boundary data, such as the shot-gather data. In this paper, we first discuss the decomposition of a two-way wave equation and investigate four different approaches to approximate the square-root operator. Using the known shot-gather data as input, an unconditional stable hybrid method for the wavefield extrapolation is presented. The most attractive feature of the proposed method is that it has a natural parallel characteristic and can be effectively implemented using a cluster of PCs, in which each processor performs its own shot-gather imaging independently. To demonstrate the computational efficiency and the power of the parallel hybrid algorithm, we present two case studies: one is the well-known SEG/EAEG subsalt model which has been commonly used for validation of the prestack depth imaging algorithms, and the other is the application to a 3D wavefield extrapolation problem with real data provided by the China National Petroleum Corporation. The results clearly show the capability of the proposed method, and it demonstrates that the algorithm can be effectively implemented as a practical engineering tool for 3D prestack depth imaging.