A computer implemented method of migrating seismic data comprising: (a) determining a wave field in an anisotropic subsurface of the Earth, comprising (a.1) providing a decoupled quasi-acoustic single wave mode wave equation based on spatially varied anisotropic parameters, wherein the wave equation comprises a spatial differential operator and a scalar operator and wherein the anisotropic parameters are separated out from the spatial differential operator, such that the wave equation has the form of either ∂ 2 u / ∂t 2 = S∇ 2 u or ∂ 2 u / ∂t 2 = ∇ · (S∇ u ) wherein u is the wave field, t is time, ∇ denotes the gradient, ∇ · denotes divergence, ∇ 2 denotes the Laplacian and S is the scalar operator comprising all the anisotropic parameters of the wave equation, and (a.2) numerically solving with a processor the wave equation by first defining a starting wave field that provides a boundary condition and second consecutively computing the spatial differential operator value and the scalar operator value to determine the wave field in the anisotropic subsurface for each predefined location within the subsurface at different wave propagation times; (b) providing first seismic data received at detectors and associated with the anisotropic subsurface of the Earth; and (c) performing seismic migration to process the first data using the determined wave field to produce second data associated with said subsurface. (a) determining a wave field in an anisotropic subsurface of the Earth, comprising (a.1) providing a decoupled quasi-acoustic single wave mode wave equation based on spatially varied anisotropic parameters, wherein the wave equation comprises a spatial differential operator and a scalar operator and wherein the anisotropic parameters are separated out from the spatial differential operator, such that the wave equation has the form of either ∂ 2 u / ∂t 2 = S∇ 2 u or ∂ 2 u / ∂t 2 = ∇ · (S∇ u ) wherein u is the wave field, t is time, ∇ denotes the gradient, ∇ · denotes divergence, ∇ 2 denotes the Laplacian and S is the scalar operator comprising all the anisotropic parameters of the wave equation, and (a.2) numerically solving with a processor the wave equation by first defining a starting wave field that provides a boundary condition and second consecutively computing the spatial differential operator value and the scalar operator value to determine the wave field in the anisotropic subsurface for each predefined location within the subsurface at different wave propagation times; (a.1) providing a decoupled quasi-acoustic single wave mode wave equation based on spatially varied anisotropic parameters, wherein the wave equation comprises a spatial differential operator and a scalar operator and wherein the anisotropic parameters are separated out from the spatial differential operator, such that the wave equation has the form of either ∂ 2 u / ∂t 2 = S∇ 2 u or ∂ 2 u / ∂t 2 = ∇ · (S∇ u ) wherein u is the wave field, t is time, ∇ denotes the gradient, ∇ · denotes divergence, ∇ 2 denotes the Laplacian and S is the scalar operator comprising all the anisotropic parameters of the wave equation, and (a.2) numerically solving with a processor the wave equation by first defining a starting wave field that provides a boundary condition and second consecutively computing the spatial differential operator value and the scalar operator value to determine the wave field in the anisotropic subsurface for each predefined location within the subsurface at different wave propagation times; (b) providing first seismic data received at detectors and associated with the anisotropic subsurface of the Earth; and (c) performing seismic migration to process the first data using the determined wave field to produce second data associated with said subsurface. <EOS>