Nowadays, subsalt oil and gas exploration is drawing more and more attention from the hydrocarbon industry. Hydrocarbon exploitation requires detailed geological information beneath the surface. Seismic imaging is a powerful tool employed by the hydrocarbon industry to provide subsurface characterization and monitoring information. Traditional wave-equation migration algorithms are based on the one- way-in-depth propagation using the scalar wave equation. These algorithms focus on downward continuing the upcoming waves. However, it is still really difficult for conventional seismic imaging methods, which have dip limitations, to get a correct image for the edge and shape of the salt body and the corresponding subsalt structure. The dip limitation problem in seismic imaging can be solved completely by switching to Reverse-Time Migration (RTM). Unlike old methods, which deal with the one-way wave equation, RTM propagator is two-way and, as a result, it no longer imposes dip limitations on the image. It can also handle complex waveforms, including prismatic waves. Therefore it is a powerful tool for subsalt imaging. RTM involves wave extrapolation forward and backward in time. In order to accurately and efficiently extrapolate the wavefield in heterogeneous media, I develop three novel methods for seismic wave modeling in both isotropic and tilted transversely isotropic (TTI) media. These methods overcome the space-wavenumber mixed-domain problem when solving the acoustic two-way wave equation. The first method involves cascading a Fourier Transform operator and a finite difference (FD) operator to form a chain operator: Fourier Finite Differences (FFD). The second method is lowrank finite differences (LFD), whose FD schemes are derived from the lowrank approximation of the mixed-domain operator and are represented using adapted coefficients. The third method is lowrank Fourier finite differences (LFFD), which use LFD to improve the accuracy of TTI FFD mothod. The first method, FFD, may have an advantage in efficiency, because it uses only one pair of multidimensional forward and inverse FFTs (fast Fourier transforms) per time step. The second method, LFD, as an accurate FD method, is free of FFTs and in return more suitable for massively parallel computing. It can also be applied to the FFD method to reduce the dispersion in TTI case, which results in the third method, LFFD. LFD and LFFD are based on lowrank approx- imation which is a general method to handle mixed-domain operators and can be easily applied to more complicated mixed-domain operators. I show pseudo-acoustic modeling in orthorhombic media by lowrank approximation as an example. / text